線形計画法の設定、ソルバーベース

ソルバー形式への問題の変換

この例では、ソルバーベースのアプローチを使用して、問題を数学的な形式から Optimization Toolbox™ ソルバーの構文に変換する方法を示します。この問題は線形計画法ですが、この手法はすべてのソルバーに適用できます。

問題に含まれる変数や式は、化学工場の運営モデルを表しており、Edgar と Himmelblau による [1] から引用しています。この問題について説明するビデオが 2 つあります。

Mathematical Modeling with Optimization, Part 1 では問題が図解され、図からモデルの説明の数式を作成する方法が説明されています。
Optimization Modeling, Part 2: Converting to Solver Form では、これらの数式を Optimization Toolbox ソルバーの構文に変換する方法が説明されています。また、問題の解法や、結果の解釈方法も示されています。

この例の残りの部分では、問題をソルバー構文に変換することのみに注目します。この例は次のビデオとほぼ同じです。Optimization Modeling, Part 2: Converting to Solver Form。ビデオと例の主な相違点は、この例ではハッシュキーと同様の名前付き変数 (つまり、インデックス変数) の使用方法が示されているという点です。これは変数を 1 つのベクトルに統合に説明されています。

モデルの説明

ビデオ Mathematical Modeling with Optimization, Part 1 では、問題を数学的な形式に変換する方法の 1 つとして、次の手順が示されています。

問題の概要の理解
ゴールの特定 (何らかの要素の最大化または最小化)
変数の識別 (命名)
制約の特定
制御できる変数の判別
すべての数量の数学的表記による指定
モデルの完全性と正確さのチェック

このセクションの変数の意味については、ビデオ Mathematical Modeling with Optimization, Part 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
すべての変数が正です。

解法

この最適化の問題を解くには、以下の手順に従います。

ソルバーの選択
変数を 1 つのベクトルに統合
範囲制約の記述
線形不等式制約の記述
線形等式制約の記述
目的の記述
linprog を使用した問題の解法
解の検証

この手順は、次のビデオでも説明されています。Optimization Modeling, Part 2: Converting to Solver Form。

ソルバーの選択

この問題に対して適切なソルバーを見つけるには、最適化の意思決定表を参考にしてください。この表では、問題を目的関数の種類と、制約の種類で分類するように求められます。この問題では、目的関数、制約ともに線形です。判定表では linprog ソルバーを使用することが推奨されています。

Optimization Toolbox の関数が扱う問題や関数 linprog のリファレンスページに説明されているように、linprog ソルバーでは次の形式の問題を解きます。

\min_{x} f^{T} x such that {\begin{matrix} A \cdot x \leq b, \\ A e q \cdot x = b e q, \\ l b \leq x \leq u b . \end{matrix}

(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 の行列とベクトルです。

変数を 1 つのベクトルに統合

モデルの説明の方程式には 16 個の変数が含まれています。これらの変数を 1 つのベクトルに統合します。変数のベクトルの名前は、式 1 の x です。次数を決定し、変数を使用して x の要素を作成します。

次のコードでは、変数の名前の cell 配列を使用してベクトルを作成します。

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

これらのコマンドを実行すると、次の名前付き変数がワークスペースに作成されます。

Listing of variables: BF1 = 8, BF2 = 9, and so forth

これらの名前付き変数は、x の要素のインデックス番号を表します。名前付き変数を作成する必要はありません。ビデオ Optimization Modeling, Part 2: Converting to Solver Form では、x の成分のインデックス番号を使用するだけで問題を解く方法を示しています。

範囲制約の記述

モデルの説明の方程式には、下限のある変数が 4 つ、上限のある変数が 6 つ含まれています。下限が指定されている変数は以下のとおりです。

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

また、すべての変数は正なので、下限はゼロになります。

下限ベクトル lb を 0 のベクトルとして作成し、4 つの下限を追加します。

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 のベクトルとして作成し、6 つの上限を追加します。

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

線形不等式制約の記述

モデルの説明の方程式には、次の 3 つの線形不等式が含まれています。

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^® ワークスペースで、16 個の変数の 3 つの線形不等式に合わせ、A 行列を 3 行 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;

線形等式制約の記述

モデルの説明の方程式には、8 つの線形等式が含まれています。

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 ワークスペースで、16 個の変数の 8 つの線形等式に合わせ、Aeq 行列を 8 行 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

ビデオ Optimization Modeling, Part 2: Converting to Solver Form では、元の問題に関して、これらの特性が説明されています。

参考文献

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