整数モデリングと論理モデリング

整数制約を使用すると、次のような重要な特徴をもつモデルを作成できます。

含意 (「条件 A が成り立つならば条件 B が成り立つ。」など)
取引のコストまたは準備のコスト (「商品を 1 つも買わなければ、その商品のコストはゼロであるが、1 つ以上の商品を買えば、そのコストは取引コスト $A に商品当たり $B を加えたものである。」など)詳細は、例: 固定コストを参照してください。
論理制約 (「気密扉 A と気密扉 B を同時に開けることはできない。」など)

多くのモデリング問題は、指標変数を使用する論理モデルに対応しています。このトピックでは、指標変数および論理モデルの使い方について説明します。これらのモデルは、変数 x と定数 M が不等式 –M ≤ x ≤ M を満たすとする "Big-M" の定式化に基づいています。

最適化における制約には暗黙的な「かつ」があることを思い出してください。制約 c1、c2 および c3 が満たされるのは、3 つの制約がすべて満たされる、つまり c1 かつ c2 かつ c3 である場合となります。

Big-M の定式化

定数 M によって上限が適用される連続変数 x があるとします。

x ≤ M.

x > 0 の場合は常に z = 1 となるように、バイナリ指標変数 z を x に関連付けるとします。そのためには、制約を含めて Big-M の定式化を使用します。

x ≤ M z.

この制約によって、x > 0 の場合は必ず z = 1 となります。実関数とバイナリ指標変数を使用した論理制約の表現の説明にあるように、Big-M の定式化はさらに応用できます。

基本的な問題: 貯水槽の流れ

整数の時間にわたって貯水槽をモデル化するとします。固定範囲内における正の整数時間 t において、貯水槽は流入または放出のいずれの場合も最大 M までの連続量として、水量 xin(t) を受け入れるか、または xout(t) を放出することができます。モデルでは、必ず xin(t) > 0 ならば xout(t) = 0、および xout(t) > 0 ならば xin(t) = 0 を適用しなければなりません。この制約はどのようにモデル化するのでしょうか?1 つは次のような制約を設定する方法です。

xin(t) * xout(t) = 0.

ただし、この制約によって問題が非線形になりますが、ソルバーは通常、この種の制約を苦手としています。

この制約は、xin(t) > 0 の場合に常に 1 となるバイナリ指標変数 zin(t) と、xout(t) > 0 の場合に常に 1 となる類似のバイナリ指標変数 zout(t) を使用すると、よりきれいに実装できます。そのような変数があると仮定して制約を追加します。

zin(t) + zout(t) ≤ 1,

これにより、xin(t) と xout(t) がどちらも正になることはありません。

xin(t) > 0 の場合は常に zin(t) = 1 となるようにするには、Big-M の定式化を使用します。xin(t) は、すべての t について正の数値 M によって上限が適用されると仮定します。次の制約を含めます。

xin(t) ≤ M*zin(t).

xin(t) = 0 の場合は常に zin(t) = 0 となるようにするため、目的関数に zin(t) を含めます。このようにして目的関数を最小化すると、zin(t) が可能な限り 0 になります。

同様に、zout(t) と xout(t) を関連付けるには、次の制約を組み込みます。

xout(t) <= M*zout(t).

要約すると、xin(t) および xout(t) が同時に正にならないように制約を適用するには、各時間 t について 2 つのバイナリ変数 zin(t) と zout(t) を作成し、次の 3 つの制約を含めます。

xin(t) <= M*zin(t) % Ensures that zin(t) = 1 whenever xin(t) > 0.

xout(t) <= M*zout(t) % Ensures that zout(t) = 1 whenever xout(t) > 0.

zin(t) + zout(t) <= 1 % Ensures that zin(t) and zout(t) are not both positive.

問題ベースのアプローチでの MATLAB^® コマンドは次のとおりです。

T = 50; % Number of times
M = 40; % Maximum size of x variables
xin = optimvar('xin',T,'LowerBound',0,'UpperBound',M);
xout = optimvar('xout',T,'LowerBound',0,'UpperBound',M);
zin = optimvar('zin',T,'Type','integer','LowerBound',0,'UpperBound',1);
zout = optimvar('zout',T,'Type','integer','LowerBound',0,'UpperBound',1);
prob = optimproblem;

xinzin = xin <= M*zin;
xoutzout = xout <= M*zout;
zinzout = zin + zout <= 1;
prob.Constraints.xinzin = xinzin;
prob.Constraints.xoutzout = xoutzout;
prob.Constraints.zinzout = zinzout;

prob.Objective = sum(zin + zout);

次の点に注意してください。

最適化変数および指標変数と同様、すべての制約は長さ T をもつ制約のベクトルである。
すべての制約は、ループの中ではなく 1 つのステートメントにより定義される。これによって最高のパフォーマンスが得られる。

バイナリ変数を使用した論理制約の表現

このセクションには、論理ステートメントおよび対応する MATLAB コマンド (バイナリ変数を使用) が記載されています。ステートメントでは、変数 z、w および f がバイナリ最適化変数、すなわち各変数の型が "integer" で、下限が 0 および上限が 1 であると仮定しています。

説明	論理ステートメント	MATLAB コマンド
`z` と `w` は反対の値をもつ	`z = not w`	z = 1 - w;
`z` または `w` の少なくとも 1 つが true である	`z or w`	z + w >= 1;
`z` または `w` の多くとも 1 つが true である	`(not z) or (not w)`	z + w <= 1;
`z` が true または `w` が true の場合にのみ `f` が true となる	`f` = `z or w`	f >= z; f >= w; f <= z + w;
`z` と `w` がともに true の場合にのみ `f` が true となる	`f` = `z and w`	f <= z; f <= w; f >= z + w - 1;
`z` または `w` のいずれかが true の場合にのみ `f` が true となる	`f` = `z xor w`	f >= z - w; f >= w - z; f <= z + w; f <= 2 - (z + w);

実関数とバイナリ指標変数を使用した論理制約の表現

このセクションでは、実関数 g(x) とバイナリ変数 z を関連付けます。通常は、指標変数として z を問題に導入し、g(x) > 0 の場合の指標など、問題の何らかの側面をモデル化します。これらの制約はすべて、Big-M の定式化に基づいています。

定数 M と関数 g(x) が次の境界条件を満たすと仮定します。

–M ≤ g(x) ≤ M.

条件	制約のコード
z = 1 の場合、g(x) ≤ 0	g(x) <= M*(1 - z);
z = 1 の場合、g(x) ≥ 0	g(x) >= -M*(1 - z);
z = 1 の場合、g(x) = 0	g(x) <= M(1 - z); g(x) >= -M(1 - z);
g(x) ≤ 0 の場合、z = 1	g(x) >= -M*z;
g(x) ≥ 0 の場合、z = 1	g(x) <= M*z;
g(x) > 0 の場合、z = 1 g(x) < 0 の場合、z = 0	g(x) < 0 を示す新しいバイナリ変数 `z1` を作成します。 g(x) <= Mz; g(x) >= -Mz1; z + z1 == 1; g(x) = 0 の場合、この定式化は不確定です。

論理制約の組み合わせによる新しい式の作成

バイナリ指標変数と共に前の論理制約を使用して、新しい式を実装するコードを作成します。

条件制約のコード

条件	制約のコード
範囲制約 M を満たすスカラー関数 g(x) と h(x) について、次の条件を実装する。 g(x) ≥ 0 の場合、h(x) ≥ 0	g(x) ≥ 0 であることを示すバイナリ指標変数 `z` を導入します。次に、z = 1 ならば h(x) ≥ 0 という制約を導入します。 g(x) <= Mz; h(x) >= -M(1 - z);
g(x) = z*x (z はバイナリ変数)	この条件を次の 2 つの制約として表します。 z = 1 の場合、g(x) - x = 0 z = 0 の場合、g(x) = 0 g(x) <= x + M(1-z); g(x) >= x - M(1-z); g(x) <= Mz; g(x) >= -Mz;

範囲制約 M を満たすスカラー関数 g(x) と h(x) について、次の条件を実装する。

g(x) ≥ 0 の場合、h(x) ≥ 0

g(x) ≥ 0 であることを示すバイナリ指標変数 z を導入します。次に、z = 1 ならば h(x) ≥ 0 という制約を導入します。

g(x) <= M*z;
h(x) >= -M*(1 - z);

g(x) = z*x (z はバイナリ変数)

この条件を次の 2 つの制約として表します。

z = 1 の場合、g(x) - x = 0

z = 0 の場合、g(x) = 0

g(x) <= x + M*(1-z);
g(x) >= x - M*(1-z);
g(x) <= M*z;
g(x) >= -M*z;

例: 固定コスト

ある商品を量 x 生産する際のコストが次であるとします。

$cost = {\begin{array}{l} a + b x & if x > 0 \\ 0 & if x = 0. \end{array}$

この非線形コストは線形変数 x およびバイナリ指標変数 z を使用してモデル化できます。x > 0 の場合は常に z = 1 となるように制約を作成して、x = 0 の場合は常に z = 0 となるように目的関数に z を含めます。問題に、x ≤ M となるような範囲 M が含まれると仮定します。

x <= M*z; % Constraint, makes z = 1 when x > 0
cost = a*z + b*x;

コストを最小化すると、x = 0 の場合に z = 0 となります。

例: OR 制約

条件 A が成り立つ、または条件 B が成り立つ、または条件 C が成り立つ場合に適用される制約をモデル化する場合があります。そのためには、対応する条件 A、B、C が成り立つ場合を示すバイナリ指標変数 zA、zB および zC を作成し、次の追加制約を含めます。

zA + zB + zC >= 1;

別の例として、絶対値制約 |x| = 5 をモデル化します。これは x = 5 または x = –5 であることを意味します。x = 5 と x = –5 の場合を示す 2 つの指標変数 z1 および z2 をそれぞれ作成します。その後、次の制約を含めます。

z1 + z2 >= 1;

x = 5 の場合に z1 = 1 を設定する 1 つの方法は、次のような条件に対して 3 つの指標変数 z11、z12 および z13 を新たに導入することです。

x < 5 かつ z1 = 0 の場合 z11 = 1

x = 5 かつ z1 = 1 の場合 z12 = 1

x > 5 かつ z1 = 0 の場合 z13 = 1

続けて次の制約を導入します。

z11 + z12 + z13 = 1;

z11 を指定するには、次の 3 つの制約を使用します。

-(1 - z11) <= z1;
z1 <= (1 - z11);
x - 5 <= M(1 - z11);

z12 を指定するには、次の 4 つの制約を使用します。

-(1 - z12) <= z1 - 1;
z1 - 1 <= (1 - z12);
-M(1 - z12) <= x - 5;
x - 5 <= M(1 - z12);

z13 を指定するには、次の 3 つの制約を使用します。

-(1 - z13) <= z1;
z1 <= (1 - z13);
x - 5 >= -M(1 - z13);

モデルを完成させるには、z2 に対応する同様の制約 z21、z22、z23、および条件 x = –5 を指定します。

例: バイナリ二次問題のバイナリ線形問題への変換

次を最小化するバイナリ変数 x_i の問題について考えます。

x^TQx,

ここで、x は長さ N の最適化変数の列ベクトル、Q は指定された N 行 N 列の行列ですN が大きすぎない場合にこの問題を解くには、バイナリ変数 xij の N 行 N 列の配列を使用して、バイナリ二次問題をバイナリ線形問題に変換します。xij(i,j) = x(i)*x(j) の場合、目的関数を次のように表すことができます。

sum(Q.*xij,"all")

変数 xij が x(i)*x(j) と等しくなるようにするには、次の 3 つの線形不等式制約を使用します。このコードでは、問題ベースのアプローチを使用します。

N = 100; % Specify the number of variables
x = optimvar("x",N,Type="integer",LowerBound=0,UpperBound=1);
xij = optimvar("xij",N,N,Type="integer",LowerBound=0,UpperBound=1);
prob = optimproblem;

% Constraint xij = 1 when x(i) = 1 and x(j) = 1
prob.Constraints.f = xij >= repmat(x,1,N) + repmat(x',N,1) - 1;
% Constraint xij = 0 when x(i) = 0
prob.Constraints.g = xij <= repmat(x,1,N);
% Constraint xij = 0 when x(j) = 0
prob.Constraints.h = xij <= repmat(x',N,1);

prob.Objective = sum(Q.*xij,"all");

最適化変数および制約の問題が線形になったため、solve は intlinprog を呼び出して解を計算します。

[sol,fval] = solve(prob);

Solving problem using intlinprog.
...

intlinprog の解は、N ≤ 100 の場合はそれなりに高速ですが、N の値が大きくなると低速になるため、適度なパフォーマンスを得るには 200 程度の変数が限界です。

参考文献

最適化のモデリングに関する書籍で代表的なものは Williams [1] です。バイナリ指標変数の Big-M の定式化が数学的に完成していて、拡張できない理由の説明については、Hooker [2] を参照してください。数学的モデリングでバイナリ指標変数を使用する例については、Brown と Dell [3] および Stevens と Palocsay [4] を参照してください。

参照

[1] Williams, H. Paul. Model Building in Mathematical Programming, 5th Edition. Wiley, 2013.

[2] Hooker, John. A Principled Approach to MILP Modeling. Carnegie Mellon University, 2008. Available at https://coral.ise.lehigh.edu/mip-2008/talks/hooker.pdf.

[3] Brown, Gerald G. and Robert F. Dell. Formulating Integer Linear Programs: A Rogues' Gallery. INFORMS Transactions on Education 7 (2), 2007, pp. 153–159. Available at https://doi.org/10.1287/ited.7.2.153.

[4] Stevens, Scott P. and Susan W. Palocsay. Teaching Use of Binary Variables in Integer Linear Programs: Formulating Logical Conditions. INFORMS Transactions on Education 18 (1), 2017, pp. 28–36. Available at https://doi.org/10.1287/ited.2017.0177.