基于问题设置线性规划
将问题转换为求解器形式
此示例说明如何使用基于问题的方法将线性问题从数学形式转换为 Optimization Toolbox™ 求解器语法。
该问题中的变量和表达式表示一家化工厂的运作模型,该模型来自 Edgar 和 Himmelblau 合著的 [1] 中的一个示例。下面两个相关联的视频对该问题进行了描述。
使用优化进行数学建模,第 1 部分以图形方式呈现该问题,说明如何生成模型描述的数学表达式。
优化建模,第 2 部分:数学模型基于问题的解说明如何将这些数学表达式转换为 Optimization Toolbox 求解器语法。此视频说明如何求解该问题,以及如何解释结果。
此示例紧跟第 2 部分的视频,重点说明如何将问题变换为求解器语法。
模型描述
第 1 部分的视频建议采用以下方法将问题转换为数学形式:
全面了解问题。
确定目标(最大化或最小化某个函数)。
标识(命名)变量。
确定约束。
确定可以控制哪些变量。
用数学表示法指定所有量。
检查模型的完整性和正确性。
有关本节中变量的含义,请参阅第 1 部分的视频。
此优化问题是在约束条件下最小化目标函数,问题和约束均以表达式列出。
目标函数是:
0.002614 HPS + 0.0239 PP + 0.009825 EP.
约束是:
2500 ≤ P1 ≤ 6250
I1 ≤ 192,000
C ≤ 62,000
I1 - HE1 ≤ 132,000
I1 = LE1 + HE1 + C
1359.8 I1 = 1267.8 HE1 + 1251.4 LE1 + 192 C + 3413 P1
3000 ≤ P2 ≤ 9000
I2 ≤ 244,000
LE2 ≤ 142,000
I2 = LE2 + HE2
1359.8 I2 = 1267.8 HE2 + 1251.4 LE2 + 3413 P2
HPS = I1 + I2 + BF1
HPS = C + MPS + LPS
LPS = LE1 + LE2 + BF2
MPS = HE1 + HE2 + BF1 - BF2
P1 + P2 + PP ≥ 24,550
EP + PP ≥ 12,000
MPS ≥ 271,536
LPS ≥ 100,623
所有变量均为正。
第一种求解方法:为每个问题变量创建优化变量
第一种求解方法涉及为每个问题变量创建一个优化变量。在创建变量时,需要包括其边界。
P1 = optimvar('P1','LowerBound',2500,'UpperBound',6250); P2 = optimvar('P2','LowerBound',3000,'UpperBound',9000); I1 = optimvar('I1','LowerBound',0,'UpperBound',192000); I2 = optimvar('I2','LowerBound',0,'UpperBound',244000); C = optimvar('C','LowerBound',0,'UpperBound',62000); LE1 = optimvar('LE1','LowerBound',0); LE2 = optimvar('LE2','LowerBound',0,'UpperBound',142000); HE1 = optimvar('HE1','LowerBound',0); HE2 = optimvar('HE2','LowerBound',0); HPS = optimvar('HPS','LowerBound',0); MPS = optimvar('MPS','LowerBound',271536); LPS = optimvar('LPS','LowerBound',100623); BF1 = optimvar('BF1','LowerBound',0); BF2 = optimvar('BF2','LowerBound',0); EP = optimvar('EP','LowerBound',0); PP = optimvar('PP','LowerBound',0);
创建问题和目标
创建一个优化问题容器。在问题中包含目标函数。
linprob = optimproblem('Objective',0.002614*HPS + 0.0239*PP + 0.009825*EP);创建并包含线性约束
问题表达式包含三个线性不等式:
I1 - HE1 ≤ 132,000EP + PP ≥ 12,000P1 + P2 + PP ≥ 24,550 | (1) |
创建以下不等式约束,并将它们包含在问题中。
linprob.Constraints.cons1 = I1 - HE1 <= 132000; linprob.Constraints.cons2 = EP + PP >= 12000; linprob.Constraints.cons3 = P1 + P2 + PP >= 24550;
该问题有八个线性等式:
I2 = LE2 + HE2LPS = LE1 + LE2 + BF2HPS = I1 + I2 + BF1HPS = C + MPS + LPSI1 = LE1 + HE1 + CMPS = HE1 + HE2 + BF1 - BF21359.8 I1 = 1267.8 HE1 + 1251.4 LE1 + 192 C + 3413 P11359.8 I2 = 1267.8 HE2 + 1251.4 LE2 + 3413 P2. | (2) |
还要包含以下约束。
linprob.Constraints.econs1 = LE2 + HE2 == I2; linprob.Constraints.econs2 = LE1 + LE2 + BF2 == LPS; linprob.Constraints.econs3 = I1 + I2 + BF1 == HPS; linprob.Constraints.econs4 = C + MPS + LPS == HPS; linprob.Constraints.econs5 = LE1 + HE1 + C == I1; linprob.Constraints.econs6 = HE1 + HE2 + BF1 == BF2 + MPS; linprob.Constraints.econs7 = 1267.8*HE1 + 1251.4*LE1 + 192*C + 3413*P1 == 1359.8*I1; linprob.Constraints.econs8 = 1267.8*HE2 + 1251.4*LE2 + 3413*P2 == 1359.8*I2;
求解问题
问题表示已完成。使用 solve 求解问题。
linsol = solve(linprob);
Optimal solution found.
检查解
计算目标函数。(您也可以通过调用 [linsol,fval] = solve(linprob) 来获取此值。)
evaluate(linprob.Objective,linsol)
ans = 1.2703e+03
如果以最低成本法运营该工厂,则成本为 1,207.30 美元。
检查解变量值。
tbl = struct2table(linsol)
tbl =
1×16 table
BF1 BF2 C EP HE1 HE2 HPS I1 I2 LE1 LE2 LPS MPS P1 P2 PP
___ ___ ______ ______ __________ __________ __________ __________ ________ ___ __________ __________ __________ ____ ______ _____
0 0 8169.7 760.71 1.2816e+05 1.4338e+05 3.8033e+05 1.3633e+05 2.44e+05 0 1.0062e+05 1.0062e+05 2.7154e+05 6250 7060.7 11239
此表太宽,其内容不易查看。堆叠变量以垂直排列它们。
vars = {'P1','P2','I1','I2','C','LE1','LE2','HE1','HE2',...
'HPS','MPS','LPS','BF1','BF2','EP','PP'};
outputvars = stack(tbl,vars,'NewDataVariableName','Amt','IndexVariableName','Var')outputvars =
16×2 table
Var Amt
___ __________
P1 6250
P2 7060.7
I1 1.3633e+05
I2 2.44e+05
C 8169.7
LE1 0
LE2 1.0062e+05
HE1 1.2816e+05
HE2 1.4338e+05
HPS 3.8033e+05
MPS 2.7154e+05
LPS 1.0062e+05
BF1 0
BF2 0
EP 760.71
PP 11239BF1、BF2和LE1为0,该值是它们的下界。I2为244,000,该值是它的上界。目标函数(成本)的非零分量是
HPS—380,328.74PP—11,239.29EP—760.71
第 2 部分的视频根据原始问题解释了这些特征。
第二种求解方法:创建一个优化变量和索引
您也可以只使用一个优化变量来求解该问题,该优化变量的索引包含问题变量的名称。通过这种方法,可一次性使所有问题变量的下界为零。
vars = {'P1','P2','I1','I2','C','LE1','LE2','HE1','HE2',...
'HPS','MPS','LPS','BF1','BF2','EP','PP'};
x = optimvar('x',vars,'LowerBound',0);设置变量边界
使用圆点表示法包括变量的边界。
x('P1').LowerBound = 2500; x('P2').LowerBound = 3000; x('MPS').LowerBound = 271536; x('LPS').LowerBound = 100623; x('P1').UpperBound = 6250; x('P2').UpperBound = 9000; x('I1').UpperBound = 192000; x('I2').UpperBound = 244000; x('C').UpperBound = 62000; x('LE2').UpperBound = 142000;
创建问题、线性约束和解
问题设置的其余部分类似于使用单独变量的设置。不同之处在于,您是通过变量的名称(如 P1)而非其索引 x('P1') 来寻址变量的。
创建问题对象,包括线性约束,并求解问题。
linprob = optimproblem('Objective',0.002614*x('HPS') + 0.0239*x('PP') + 0.009825*x('EP')); linprob.Constraints.cons1 = x('I1') - x('HE1') <= 132000; linprob.Constraints.cons2 = x('EP') + x('PP') >= 12000; linprob.Constraints.cons3 = x('P1') + x('P2') + x('PP') >= 24550; linprob.Constraints.econs1 = x('LE2') + x('HE2') == x('I2'); linprob.Constraints.econs2 = x('LE1') + x('LE2') + x('BF2') == x('LPS'); linprob.Constraints.econs3 = x('I1') + x('I2') + x('BF1') == x('HPS'); linprob.Constraints.econs4 = x('C') + x('MPS') + x('LPS') == x('HPS'); linprob.Constraints.econs5 = x('LE1') + x('HE1') + x('C') == x('I1'); linprob.Constraints.econs6 = x('HE1') + x('HE2') + x('BF1') == x('BF2') + x('MPS'); linprob.Constraints.econs7 = 1267.8*x('HE1') + 1251.4*x('LE1') + 192*x('C') + 3413*x('P1') == 1359.8*x('I1'); linprob.Constraints.econs8 = 1267.8*x('HE2') + 1251.4*x('LE2') + 3413*x('P2') == 1359.8*x('I2'); [linsol,fval] = solve(linprob);
Optimal solution found.
检查索引解
使用竖表查看解。
tbl = table(vars',linsol.x')
tbl =
16×2 table
Var1 Var2
_____ __________
'P1' 6250
'P2' 7060.7
'I1' 1.3633e+05
'I2' 2.44e+05
'C' 8169.7
'LE1' 0
'LE2' 1.0062e+05
'HE1' 1.2816e+05
'HE2' 1.4338e+05
'HPS' 3.8033e+05
'MPS' 2.7154e+05
'LPS' 1.0062e+05
'BF1' 0
'BF2' 0
'EP' 760.71
'PP' 11239参考书目
[1] Edgar, Thomas F., and David M. Himmelblau. Optimization of Chemical Processes. New York: McGraw-Hill, 1987.