基于求解器设置线性规划

将问题转换为求解器形式

此示例说明如何使用基于求解器的方法将问题从数学形式转换为 Optimization Toolbox™ 求解器语法。虽然这里的示例是线性规划问题，但这些方法适用于所有求解器。

该问题中的变量和表达式表示一家化工厂的运作模型，该模型来自 Edgar 和 Himmelblau 合著的 [1] 中的一个示例。我们提供了两个视频来说明该问题。

• 使用优化进行数学建模，第 1 部分以图形形式展示该问题。它用图形说明如何生成模型描述的数学表达式。

• 优化建模，第 2 部分：转换为求解器形式介绍如何将这些数学表达式转换为 Optimization Toolbox 求解器语法。此视频说明如何求解该问题，以及如何解释结果。

此示例的其余部分只涉及将问题转换为求解器语法。此示例与视频优化建模，第 2 部分：转换为求解器形式。视频和示例的主要区别在于，此示例说明如何使用命名变量或索引变量，这些变量类似于哈希键。这种差异体现在将变量合并成一个向量中。

模型描述

视频使用优化进行数学建模，第 1 部分建议将问题转换为数学形式的一种方法是：

1. 全面了解问题

2. 确定目标（最大化或最小化某个函数）

3. 确定（命名）变量

4. 确定约束

5. 确定可以控制哪些变量

6. 用数学表示法指定所有量

7. 检查模型的完整性和正确性

有关本节中变量的含义，请观看视频使用优化进行数学建模，第 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
所有变量均为正。

求解方法

要求解优化问题，请执行以下步骤。

• 选择求解器

• 将变量合并成一个向量

• 编写边界约束

• 编写线性不等式约束

• 编写线性等式约束

• 编写目标

• 使用 linprog 求解问题

• 检查解

关于这些步骤，还可观看视频优化建模，第 2 部分：转换为求解器形式。

选择求解器

要找到求解此问题的合适求解器，请参考优化决策表。该表要求您根据目标函数的类型和约束的类型对问题进行分类。对于此问题，目标函数是线性的，约束也是线性的。决策表推荐使用 linprog 求解器。

如 由 Optimization Toolbox 函数处理的问题或 linprog 函数参考页中所述，linprog 求解器可求解以下形式的问题

$$\underset{x}{\min }{f}^{T}x\text{ such that }\{\begin{array}{c}A\cdot x\le b,\\ Aeq\cdot x=beq,\\ lb\le x\le ub.\end{array}$$

(1)

• f^Tx 表示由常量组成的行向量 f 乘以由变量组成的列向量 x。换言之，

  f^Tx = f(1)x(1) + f(2)x(2) + ... + f(n)x(n)，

  其中，n 是 f 的长度。

• A x ≤ b 表示线性不等式。A 是 k×n 矩阵，其中 k 是不等式的数目，n 是变量的数目（大小为 x）。b 是长度为 k 的向量。有关详细信息，请参阅线性不等式约束。

• Aeq x = beq 表示线性等式。Aeq 是一个 m×n 矩阵，其中 m 是等式的个数，n 是变量的个数（大小为 x）。beq 是长度为 m 的向量。有关详细信息，请参阅线性等式约束。

• lb ≤ x ≤ ub 表示向量 x 中的每个元素必须大于 lb 的对应元素，并且必须小于 ub 的对应元素。有关详细信息，请参阅边界约束。

如函数参考页所示，linprog 求解器的语法是

[x fval] = linprog(f,A,b,Aeq,beq,lb,ub);

linprog 求解器的输入是公式 1 中的矩阵和向量。

将变量合并成一个向量

模型描述的方程中有 16 个变量。将这些变量放入一个向量中。由变量组成的向量的名称在公式 1 中是 x。决定阶次，并基于变量构造 x 的分量。

以下代码使用变量名称组成的元胞数组来构造向量。

variables = {'I1','I2','HE1','HE2','LE1','LE2','C','BF1',...
    'BF2','HPS','MPS','LPS','P1','P2','PP','EP'};
N = length(variables); 
% create variables for indexing 
for v = 1:N 
   eval([variables{v},' = ', num2str(v),';']); 
end

执行这些命令会在工作区中创建以下命名变量：

这些命名变量表示 x 分量的索引编号。您不必创建命名变量。视频优化建模，第 2 部分：转换为求解器形式显示如何使用 x 分量的索引编号来轻松求解问题。

编写边界约束

在模型描述的方程中有四个变量包含下界，六个变量包含上界。下界：

P1 ≥ 2500
P2 ≥ 3000
MPS ≥ 271,536
LPS ≥ 100,623。

此外，所有变量均为正，这意味着其下界为零。

创建由 0 组成的下界向量 lb，然后添加其他四个下界。

lb = zeros(size(variables));
lb([P1,P2,MPS,LPS]) = ...
    [2500,3000,271536,100623];

具有上界的变量有：

P1 ≤ 6250
P2 ≤ 9000
I1 ≤ 192,000
I2 ≤ 244,000
C ≤ 62,000
LE2 ≤ 142000。

创建由 Inf 组成的上界向量，然后添加六个上界。

ub = Inf(size(variables));
ub([P1,P2,I1,I2,C,LE2]) = ...
   [6250,9000,192000,244000,62000,142000];

编写线性不等式约束

模型描述的方程中有三个线性不等式：

I1 - HE1 ≤ 132,000
EP + PP ≥ 12,000
P1 + P2 + PP ≥ 24,550。

为了使方程具有 A x≤b 形式，需要将所有变量放在不等式的左侧。所有这些方程均已采用该形式。在适当情况下，通过乘以 –1，确保每个不等式均为“小于”形式：

I1 - HE1 ≤ 132,000
-EP - PP ≤ -12,000
-P1 - P2 - PP ≤ -24,550。

在您的 MATLAB^® 工作区中，将 A 矩阵创建为 3×16 零矩阵，对应于采用 16 个变量的 3 个线性不等式。创建具有三个分量的 b 向量。

A = zeros(3,16);
A(1,I1) = 1; A(1,HE1) = -1; b(1) = 132000;
A(2,EP) = -1; A(2,PP) = -1; b(2) = -12000;
A(3,[P1,P2,PP]) = [-1,-1,-1];
b(3) = -24550;

编写线性等式约束

模型描述的方程中有八个线性方程：

I2 = LE2 + HE2
LPS = LE1 + LE2 + BF2
HPS = I1 + I2 + BF1
HPS = C + MPS + LPS
I1 = LE1 + HE1 + C
MPS = HE1 + HE2 + BF1 - BF2
1359.8 I1 = 1267.8 HE1 + 1251.4 LE1 + 192 C + 3413 P1
1359.8 I2 = 1267.8 HE2 + 1251.4 LE2 + 3413 P2.

为了使方程具有 Aeq x=beq 形式，需要将所有变量放在方程的一侧。方程变为：

LE2 + HE2 - I2 = 0
LE1 + LE2 + BF2 - LPS = 0
I1 + I2 + BF1 - HPS = 0
C + MPS + LPS - HPS = 0
LE1 + HE1 + C - I1 = 0
HE1 + HE2 + BF1 - BF2 - MPS = 0
1267.8 HE1 + 1251.4 LE1 + 192 C + 3413 P1 - 1359.8 I1 = 0
1267.8 HE2 + 1251.4 LE2 + 3413 P2 - 1359.8 I2 = 0.

现在编写对应于这些方程的 Aeq 矩阵和 beq 向量。在您的 MATLAB 工作区中，将 Aeq 矩阵创建为 8×16 零矩阵，对应于包含 16 个变量的 8 个线性方程。创建具有八个分量且值均为零的 beq 向量。

Aeq = zeros(8,16); beq = zeros(8,1);
Aeq(1,[LE2,HE2,I2]) = [1,1,-1];
Aeq(2,[LE1,LE2,BF2,LPS]) = [1,1,1,-1];
Aeq(3,[I1,I2,BF1,HPS]) = [1,1,1,-1];
Aeq(4,[C,MPS,LPS,HPS]) = [1,1,1,-1];
Aeq(5,[LE1,HE1,C,I1]) = [1,1,1,-1];
Aeq(6,[HE1,HE2,BF1,BF2,MPS]) = [1,1,1,-1,-1];
Aeq(7,[HE1,LE1,C,P1,I1]) = [1267.8,1251.4,192,3413,-1359.8];
Aeq(8,[HE2,LE2,P2,I2]) = [1267.8,1251.4,3413,-1359.8];

编写目标

目标函数是

f^Tx = 0.002614 HPS + 0.0239 PP + 0.009825 EP。

将此表达式写成由 x 向量的乘数组成的向量 f：

f = zeros(size(variables));
f([HPS PP EP]) = [0.002614 0.0239 0.009825];

使用 linprog 求解问题

您现在已经有 linprog 求解器所需的输入。调用该求解器并以设定的格式打印输出：

options = optimoptions('linprog','Algorithm','dual-simplex');
[x fval] = linprog(f,A,b,Aeq,beq,lb,ub,options);
for d = 1:N
  fprintf('%12.2f \t%s\n',x(d),variables{d}) 
end
fval

结果为：

Optimal solution found.
   136328.74 	I1
   244000.00 	I2
   128159.00 	HE1
   143377.00 	HE2
        0.00 	LE1
   100623.00 	LE2
     8169.74 	C
        0.00 	BF1
        0.00 	BF2
   380328.74 	HPS
   271536.00 	MPS
   100623.00 	LPS
     6250.00 	P1
     7060.71 	P2
    11239.29 	PP
      760.71 	EP

fval =
  1.2703e+03

检查解

fval 输出提供目标函数在任何可行点的最小值。

解向量 x 是目标函数具有最小值的点。请注意：

• BF1、BF2 和 LE1 为 0，该值是它们的下界。

• I2 为 244,000，该值是它的上界。

• f 向量的非零分量为

  • HPS — 380,328.74

  • PP — 11,239.29

  • EP — 760.71

视频优化建模，第 2 部分：转换为求解器形式根据原始问题给出这些特征的解释。

参考书目

[1] Edgar, Thomas F., and David M. Himmelblau. Optimization of Chemical Processes. McGraw-Hill, New York, 1988.

相关主题

• 基于问题设置线性规划
• 基于求解器的优化问题设置