【数学建模】基于问题的线性规划和混合整数规划求解
基于问题的线性规划和混合整数规划求解(problem_based linear programming)。
在MatLab中，线性规划类问题的求解基本上有两种解决方案，最简单的是直接调用求解器(solver)求解，这叫做solver-based linear programming，求解的命令是linprog和intlinprog。这种方案简单，但需要我们手动列出所有系数矩阵、向量(Ax<=b;Aeq.x<=beq;and so on)。当约束增多，这个工作几乎是不可行的。
MatLab提供了基于问题的求解方案(problem_based linear programming)。这种方案更加直观，缺点是需要自己一步步实现，它实际上上也是调用了求解器，使用单纯形法、内点法等方法求解(可以指定)。
基于问题求解步骤如下：
1. 定义优化问题
prob = optimproblem('ObjectiveSense','maximize');
%参数指定优化问题为最大化目标函数
2. 定义决策变量
x1 = optimvar('name'); %可以这样一个个定义决策变量，也可以一次定义多个
X = optimvar('name',['x1','x2','x3']); %一次定义三个变量
x = optimvar('name',3,4); %当然也可以定义多维的
x = optimvar('name','LowerBound','0','UpperBound','3'); %设置上下界
x = optimivar('name','Type','integer'); %当然也可以设置整数变量
3. 确定目标函数
prob.Objective = x1^2+x2+100；%这就是使用上面的决策变量
%也可以使用下面的方法
obj = x1^2+x2+100;
prob.Objective = obj;
4. 定义约束并添加进本优化问题
const1 = x1 + x2 == 100; %注意使用逻辑运算
const2 = x1 + 2*x2 <= 100;
prob.Constraints.const1 = const1; %添加进本问题的约束，.const1是任意名
prob.Constraints.const2 = const2;
%也可以定义const时先不指定<=100，引进prob中时再指定。
5. solve求解
sol = solve(prob); %直接solve，返回sol，sol.x1是决策变量x1
[sol,fval] = solve(prob);
6. 求最优值
val = evaluate(obj,sol); %参数是目标函数和solve求解结果，返回目标函数最优值
7. 查看构建的约束问题
showproblem(prob);
下面针对一个具体的例子求解。
【Question】混合整数-线性规划求解。
购买下面材料获得25吨的钢材，使得钢材中含有5%的碳和5%的钼的化学成分,当然，下面三种种材料都是有特定成本的。
其中：有四种锭刚才可供购买，但之多买一种；三种合金任意买，可以是小数；还有废料也可购买。
这个问题是从卡尔-亨利 Westerberg, 克韦斯 Bjorklund, Eskil Hultman, "混合整数规划在瑞典钢厂的应用。接口1977年2月卷 7, 2 页 39–43, 其摘要是在http://interfaces.journal.informs.org/content/7/2/39.abstract.
| 锭 | 重量吨 | %Carbon | %Molybdenum | 成本/吨 |
|---|---|---|---|---|
| 1 | 5 |
5 |
3 |
350元 |
| 2 | 3 |
4 |
3 |
330元 |
| 3 | 4 |
5 |
4 |
310元 |
| 4 | 6 |
3 |
4 |
280元 |
| 合金 | %Carbon | %Molybdenum | 成本/吨 |
|---|---|---|---|
| 1 | 8 |
6 |
500元 |
| 2 | 7 |
7 |
450元 |
| 3 | 6 |
8 |
400元 |
| 废料 | 3 |
9 |
100元 |
基于problem的解决方案如下：
1 steelprob = optimproblem; 2 ingots = optimvar('ingots',4,'Type','integer','LowerBound',0,'UpperBound',1); 3 alloys = optimvar('alloys',3,'LowerBound',0); 4 scrap = optimvar('scrap','LowerBound',0); 5 %定义目标函数 6 %先对相应变量的成本创建表达式 7 weightIngots = [5,3,4,6]; 8 costIngots = weightIngots.*[350,330,310,280]; 9 costAlloys = [500,450,400]; 10 costScrap = 100; 11 %目标函数 12 cost = costIngots * ingots + costAlloys*alloys + costScrap*scrap; 13 %包含进问题 14 steelprob.Objective = cost; 15 %创建约束 16 %第一个约束是总重量为25吨 17 const1 = weightIngots*ingots + sum(alloys) + scrap >=25; 18 %第二个约束是Carbon的重量达标 19 carbonIngots = [5 3 5 3]/100; 20 carbonAlloys = [8 7 6]/100; 21 carbonScrap = 3/100; 22 const2 = (weightIngots.*carbonIngots)*ingots + (carbonAlloys*alloys) + carbonScrap*scrap >= 1.25; 23 %第三个约束是Molybdenum达标 24 muIngots = [3 3 4 4]/100; 25 muAlloys = [6 7 8]/100; 26 muScrap = 9/100; 27 const3 = (weightIngots.*muIngots)*ingots + (muAlloys*alloys) + muScrap*scrap >= 1.25; 28 %约束包含进问题 29 steelprob.Constraints.const1 = const1; 30 steelprob.Constraints.const2 = const2; 31 steelprob.Constraints.const3 = const3; 32 [sol,fval] = solve(steelprob)
