混合整数線形計画法 (MILP) アルゴリズム

混合整数線形計画法の定義

混合整数線形計画法 (MILP) は次を含む問題です。

線形目的関数 f^Tx、ここで f は定数の列ベクトル、x は未知数の列ベクトルです。
範囲および線形制約、ただし非線形制約を含みません (定義については、制約の作成を参照)。
x の一部の要素を整数値とするという制限をもちます。

数式では、与えられたベクトル f、lb および ub、行列 A および Aeq、対応するベクトル b および beq、一連の指数 intcon が解となるベクトル x を求めます。

$\min_{x} f^{T} x subject to {\begin{array}{l} x (intcon) are integers \\ A \cdot x \leq b \\ A e q \cdot x = b e q \\ l b \leq x \leq u b . \end{array}$

アルゴリズムの概要

intlinprog は、この基本的な方法を使用して混合整数線形計画法を解きます。intlinprog は任意の段階で問題を解くことができます。ある段階で問題が解かれた場合、intlinprog は後の段階を実行しません。

線形計画法の前処理を使用して問題のサイズを縮小します。
線形計画法を使用して初期の緩和された問題 (非線形) を解きます。
混合整数計画法の前処理を実行し、混合整数問題の LP 緩和を厳しくします。
カットの生成を実行し、混合整数問題の LP 緩和をさらに厳しくしてみます。
ヒューリスティックな方法を使用して整数実行可能解を求めます。
分枝限定法アルゴリズムを使用して最適解を体系的に求めます。このアルゴリズムは、整数変数の制限された範囲の取り得る値によって LP 緩和を解きます。最適な目的関数値の更新された範囲のシーケンスを生成しようとします。

線形計画法の前処理

混合整数線形計画法の定義により、一連の線形不等式と線形等式をエンコードする行列 A および Aeq と対応するベクトル b および beq があります。

$\begin{matrix} A \cdot x \leq b \\ A e q \cdot x = b e q . \end{matrix}$

これらの線形制約は解 x を制限します。

通常、問題の変数の数 (x の要素の数) を減らし、線形制約の数を減らすことができます。これらのリダクションの実行によってソルバーに時間がかかっている間に、通常解を求める全体的な時間を減らしてより大きな問題が解決可能になります。アルゴリズムは解をより数値的に安定させることができます。さらに、これらのアルゴリズムは実行不可能な問題を検出することができる場合があります。

前処理手順は、冗長な変数と制約を除去し、制約行列のモデルとスパースのスケーリングを改善し、変数の範囲を強化し、モデルの主実行不可能性および双対実行不可能性を検出することを目的としています。

詳細については、Andersen と Andersen [2] および Mészáros と Suhl [8] を参照してください。

線形計画法

初期の "緩和された" 問題は、混合整数線形計画法の定義と同じ目的関数と制約をもつが、整数制約をもたない線形計画問題です。x_LP (緩和された問題の解) と x (整数制約をもつ元の問題の解) を呼び出します。以下になることは明らかです。

f^Tx_LP ≤ f^Tx

これは、x_LP がより少ない制限で同じ関数を最小化するためです。

この初期の緩和された LP (ルートノード LP) と分枝限定アルゴリズムを使用中に生成されたすべての LP 緩和は線形計画法の手法を使用して解決されます。

混合整数計画法の前処理

混合整数計画法の前処理中に、intlinprog は線形不等式 A*x ≤ b を整合性制限と共に解析して以下を判定します。

問題が実行不可能。
一部の範囲を厳しくできる。
一部の不等式は冗長であるため、無視または削除できる。
一部の不等式を強化できる。
一部の整数変数を固定できる。

IntegerPreprocess オプションでは、intlinprog がいくつかの手順を取るか、すべての手順を取るか、ほとんど何の手順も取らないかを選択できます。x0 引数を含めると、intlinprog は前処理でその値を使用します。

混合整数計画法の前処理の主なゴールは、その後の分枝限定法の計算を簡素化することです。前処理では、分枝限定法が解析する可能性がある一部の無意味な部分問題の候補をすばやく事前確認して除去します。

整数前処理についての詳細は、[10]を参照してください。

カットの生成

カットは、intlinprog が問題に追加する別の線形不等式制約です。これらの不等式は、解が整数に近づくように LP 緩和の実行可能領域を制限しようとします。CutGeneration オプションにより、intlinprog が使用するカットのタイプを制御できます。

'basic' カットは次のとおりです。

混合整数丸めカット
ゴモリーカット
クリークカット
カバーカット
フローカバーカット

さらに、問題が純粋に整数 (すべての変数が整数値) である場合は、intlinprog は以下のカットも使用します。

厳密なフバタルゴモリーカット
ゼロハーフカット

'intermediate' カットには、すべての 'basic' カットに加え、次のカットが含まれます。

単純なリフト射影カット
単純なピボット削減カット
縮小分割カット

'advanced' カットには、縮小分割カットを除くすべての 'intermediate' カットに加え、次のカットが含まれます。

厳密なフバタルゴモリーカット
ゼロハーフカット

純粋に整数の問題の場合、'intermediate' では最も多くのカットタイプが使用されますが (縮小分割カットを使用するため)、'advanced' では使用されません。

もう 1 つのオプション CutMaxIterations は、intlinprog がカットの生成を反復するときの上限を指定します。

カット生成アルゴリズム (切除平面法とも呼ばれます) の詳細については、Cornuéjols[5]を参照してください。クリークカットの詳細については、Atamtürk、Nemhauser および Savelsbergh[3]を参照してください。

実行可能解を求めるためのヒューリスティックな方法

目的関数の上限を取得するには、分枝限定手順で実行可能点を検出しなければなりません。分枝限定時の LP 緩和の解は整数で実行可能であり、改善された上限を元の MILP に提供できます。特定の手法では、分枝限定前や分枝限定時に、よりすばやく実行可能点を検出できます。intlinprog はこうした手法をルートノードで使用し、分枝限定反復中は使用しません。これらの手法はヒューリスティックなものです。つまり、これらの手法は成功する可能性も、失敗する可能性もあるアルゴリズムです。

ヒューリスティックな方法には、"開始" ヒューリスティックがあります。この方法は、ソルバーが初期または新規の整数実行可能解を見つけるのに役立ちます。ヒューリスティックな方法には、"改善" ヒューリスティックもあります。この方法では、整数実行可能点から始まり、より良い整数実行可能点 (つまり、目的関数値がより低い点) を検出しようと試みます。intlinprog の改善ヒューリスティックは、'rins'、'rss'、1-opt、2-opt、およびガイド付きダイビングです。

'Heuristics' オプションを使用して intlinprog のヒューリスティックな方法を設定します。以下のオプションがあります。

オプション	説明
`'basic'` (既定の設定)	ソルバーは、丸めヒューリスティックを異なるパラメーターで 2 回実行し、ダイビングヒューリスティックを異なるパラメーターで 2 回実行した後、`'rss'` を実行します。先のヒューリスティックな方法で十分に適した整数実行可能解が導かれている場合、ソルバーは、後のヒューリスティックな方法を実行しません。
`'intermediate'`	ソルバーは、丸めヒューリスティックを異なるパラメーターで 2 回実行した後、ダイビングヒューリスティックを異なるパラメーターで 2 回実行します。整数実行可能解が存在する場合、ソルバーは `'rins'` の後に `'rss'` を実行します。`'rss'` で新しい解が見つかった場合、ソルバーは `'rins'` を再び実行します。先のヒューリスティックな方法で十分に適した整数実行可能解が導かれている場合、ソルバーは、後のヒューリスティックな方法を実行しません。
`'advanced'`	ソルバーは、丸めヒューリスティックを異なるパラメーターで 2 回実行した後、ダイビングヒューリスティックを異なるパラメーターで 2 回実行します。整数実行可能解が存在する場合、ソルバーは `'rins'` の後に `'rss'` を実行します。`'rss'` で新しい解が見つかった場合、ソルバーは `'rins'` を再び実行します。先のヒューリスティックな方法で十分に適した整数実行可能解が導かれている場合、ソルバーは、後のヒューリスティックな方法を実行しません。
`'rins'` または同等の `'rins-diving'`	`intlinprog` は現在の最適な実行可能解の整数点 (可能な場合) の近傍を検索し、新しいより適切な解を見つけます。詳細は、Danna、Rothberg および Le Pape [6]を参照してください。`'rins'` を選択すると、ソルバーは丸めヒューリスティックを異なるパラメーターで 2 回実行し、ダイビングヒューリスティックを異なるパラメーターで 2 回実行した後、`'rins'` を実行します。
`'rss'` または同等の `'rss-diving'`	`intlinprog` は、`'rins'` およびローカル分枝からアイデアを組み合わせるハイブリッド手順を適用して整数実行可能解を検索します。`'rss'` を選択すると、ソルバーは丸めヒューリスティックを異なるパラメーターで 2 回実行し、ダイビングヒューリスティックを異なるパラメーターで 2 回実行した後、`'rss'` を実行します。先のヒューリスティックな方法で十分に適した整数実行可能解が導かれている場合、ソルバーは、後のヒューリスティックな方法を実行しません。これらの設定は、同じヒューリスティックな方法を `'basic'` として実行します。
`'round'`	`intlinprog` は緩和された問題の LP 解をノードで取り、実行可能性を維持するように整数要素を丸めます。`'round'` を選択すると、ソルバーはルートノードで丸めヒューリスティックを異なるパラメーターで 2 回実行し、ダイビングヒューリスティックを異なるパラメーターで 2 回実行します。その後、ソルバーはいくつかの分枝限定ノードで丸めヒューリスティックのみを実行します。
`'round-diving'`	ソルバーは `'round'` と同様に機能しますが、ルートノードだけではなくいくつかの分枝限定ノードで (丸めヒューリスティックに加えて) ダイビングヒューリスティックも実行します。
`'diving'`	`intlinprog` は分枝限定ステップに類似したヒューリスティックな方法を使用します。ただし、ツリーの 1 つの分枝にのみ進み、他の分枝は作成しません。この単一分枝によってツリーフラグメントを迅速に "ダイブ" することができるため、"ダイビング" という名前が付けられています。現在、`intlinprog` は、6 つのダイビングヒューリスティックを次の順序で使用します。ベクトル長ダイビング係数ダイビング非整数ダイビング擬似コストダイビング直線探索ダイビングガイド付きダイビング (ソルバーが少なくとも 1 つの整数実行可能点を既に見つけている場合に適用) ダイビングヒューリスティックは通常、整数値であり、現在の解が非整数である 1 つの変数を選択します。その後、このヒューリスティックは変数が整数値になるように強制する範囲を導入し、関連する緩和された LP を再度解きます。制限する変数の選択方法がダイビングヒューリスティック間の主な相違点です。Berthold [4]のセクション 3.1 を参照してください。
`'none'`	`intlinprog` は実行可能点を検索しません。ソルバーは単に分枝限定検索で検出された実行可能点を取ります。

'intermediate' と 'advanced' の主な相違点は、'advanced' は分枝限定反復中により頻繁にヒューリスティックな方法を実行することです。

上記の表に加えて、Heuristics オプションが 'basic'、'intermediate'、または 'advanced' のとき、以下のヒューリスティックな方法が実行されます。

ZI round — このヒューリスティックな方法は、緩和された LP をアルゴリズムが解くときに常に実行されます。このヒューリスティックな方法では、他の制約に関する実行可能性に影響を与えることなく、各非整数変数を隣接する整数にシフトしようと試みます。詳細は、Hendel[7]を参照してください。
1-opt — このヒューリスティックな方法は、新しい整数実行可能解をアルゴリズムが見つけたときに常に実行されます。このヒューリスティックな方法では、目的関数値を低減しながら、他の制約に関する実行可能性に影響を与えることなく、各整数変数を隣接する整数にシフトしようと試みます。
2-opt — このヒューリスティックな方法は、新しい整数実行可能解をアルゴリズムが見つけたときに常に実行されます。2-opt は、同じ制約に影響するすべての整数変数のペアを検出します。これは、A または Aeq 制約行列の同じ行に非ゼロのエントリがあることを意味します。ペアごとに、2-opt は、1 つの整数実行可能解をとり、4 つすべての移動手段 (上上、上下、下上、下下) を使用して変数ペアの値を上下に移動し、より適切な目的関数値をもつ実行可能隣接解を探します。アルゴリズムは、ペア内の変数ごとに制約を満たし、目的関数値も改善するシフトの最大サイズ (同じ大きさ) を計算することによって、各整数変数ペアをテストします。

ヒューリスティックな方法の最初の段階では、Heuristics が 'none' でない限り、または x0 引数に初期の整数実行可能点を提供しない限り、intlinprog は、"自明" ヒューリスティックを実行します。自明ヒューリスティックでは、実行可能性に関して以下の点を確認します。

すべての零点
上限
下限 (非ゼロの場合)
"ロック" 点

"ロック" 点は、すべての変数の上限と下限が有限である問題に関してのみ定義されます。各変数の "ロック" 点は、次のように選択される上限または下限となります。各変数 j について、線形制約行列 A(:,j) に含まれる、対応する正の要素の数をカウントし、対応する負の要素の数で減算します。結果が正であれば、その変数の下限 lb(j) を使用します。そうでない場合は、その変数の上限 ub(j) を使用します。"ロック" 点が各変数の線形不等式制約の最大数を満たすよう試みられますが、それが見つからないこともあります。

各ヒューリスティックな方法が実行可能解を得て完了した後、intlinprog は出力関数とプロット関数を呼び出します。詳細は、intlinprog の出力関数とプロット関数の構文を参照してください。

x0 引数を含めると、intlinprog でより良い整数実行可能点が見つかるまで、'rins' およびガイド付きダイビングヒューリスティックでその値が使用されます。そのため、x0 を指定する場合には、'Heuristics' オプションを 'rins-diving' に設定するか、'rins' を使用する別の設定にすることによって、良好な結果を得ることができます。

分枝限定法

分枝限定法は MILP の解に収束しようとする一連の部分問題を作成します。部分問題は、解 f^Tx に一連の上限と下限を設定します。最初の上限は任意の実行可能解であり、最初の下限は緩和された問題の解です。上限の詳細は、実行可能解を求めるためのヒューリスティックな方法を参照してください。

線形計画法で説明したように、線形計画法緩和問題の解は MILP の解より低い目的関数値をもちます。また、任意の実行可能点 x_feas は以下の式を満たします。

f^Tx_feas ≥ f^Tx

これは f^Tx がすべての実行可能点で最小値であるためです。

このコンテキストでは、"ノード" は元の問題と同じ目的関数、範囲、線形制約をもつが、整数制約をもたず、線形制約または範囲の特定の変更をもつ LP です。"ルートノード" は整数制約および線形制約または範囲の変更をもたない元の問題です。つまり、ルートノードは初期の緩和された LP です。

範囲の開始から、分枝限定法はルートノードからの分枝により新しい部分問題を作成します。いくつかのルールの 1 つに従って経験則的に分枝手順が取られます。各ルールは、ある変数を整数 J 以下、または J+1 以上に制限することにより問題を分割するという考え方に基づいています。これらの 2 つの部分問題は intcon で指定された整数に対応する x_LP のエントリが整数ではないときに生じます。ここで、x_LP は緩和された問題の解です。J を負方向の丸め (切り捨て)、J+1 を正方向の丸め (切り上げ) として取ります。生成された 2 つの問題は、さらに制限をもつため f^Tx_LP 以上となる解をもちます。したがって、この手順では下限が引き上げられる可能性があります。

分枝限定法のパフォーマンスは分割する変数を選択するルール (分枝ルール) によって異なります。アルゴリズムは、以下の BranchRule オプションで設定できるこのようなルールを使用します。

'maxpscost' — 最大 "擬似コスト" をもつ非整数変数を選択します。
擬似コスト
変数 i の擬似コストは、i が現在点 x の i 要素の小数部と組み合わせた分枝変数として選択された場合の下限の変更の経験的推定に基づきます。小数部 p は、次の下部と上部の 2 つの部分で構成されます。
p_i^– = x(i) – ⌊x(i)⌋
p_i⁺ = 1 – p_i^–.
x_i^– は x(i) ≤ ⌊x(i)⌋ となるように制限された線形計画法の解とし、目的関数の変更は次のように定義されるものとします。
Δ_i^– = f^Tx_i^– – f^Tx
同様に、Δ_i⁺ は、問題が x(i) ≥ ⌈x(i)⌉ となるように制限された場合の目的関数の変更です。
変数 x_i における目的の単位変更当たりのゲインは次のようになります。
$d_{i}^{-} = \frac{Δ_{i}^{-}}{p_{i}^{-}} or d_{i}^{+} = \frac{Δ_{i}^{+}}{p_{i}^{+}} .$
s_i^– と s_i⁺ は、この時点までの分枝限定アルゴリズムを使用中の d_i^– と d_i⁺ の経験的平均値とします。経験値は、観測値が得られる前の条件として目的係数 f(i) の絶対値に初期化されます。次に、'maxpscost' ルールが i ノードに分枝され、何らかの正の重み w⁺ および w^– を付けるために次の値に最大化されます。
w^– * p_i^– * s_i^– + w⁺ * p_i⁺ * s_i⁺.
おおまかに言うと、このルールは下限を最大限に増加させる可能性がある係数を選択します。
'strongpscost' — 'maxpscost' と同様ですが、変数ごとに 1 に初期化される擬似コストの代わりに、ソルバーは、擬似コストがより信頼できる推定値となった後のみ変数で分枝を試みます。より信頼できる推定値を取得するために、ソルバーは以下を実行します (Achterberg、Koch、および Martin[1]を参照)。
- 候補となるすべての分枝変数 (現在は非整数だが整数にしなければならない変数) を現在の擬似コストベースのスコア順に並べます。
- 最高スコアの変数から順番に、現在の分枝変数に基づいて 2 つの緩和された線形計画法を実行します (変数がまだ分枝計算に使用されていない場合)。ソルバーはこれらの 2 つの解を使用して、現在の分枝変数の擬似コストを更新します。ソルバーは、早期にこのプロセスを停止して分枝の選択における時間を節減できます。
- k 個の連続する変数に対して、現在の擬似コストベースの最高スコアが変わらなくなるまで、リスト内の変数の選択を続行します。ここで、k は内部で選択された値で、通常は 5 ～ 10 の範囲です。
- 擬似コストベースの最高スコア変数で分枝を行います。ソルバーは、早期の擬似コスト推定手順の間に、既にこの変数に基づいて緩和された線形計画法を計算している場合があります。
追加の線形計画法の解が原因で、'strongpscost' 分枝の各反復では既定の 'maxpscost' より長い時間がかかります。ただし、分枝限定法による反復の数は通常減少するため、'strongpscost' 手法では全体の時間を短縮できます。
'reliability' — 'strongpscost' と同様ですが、初期化されていない擬似コスト分枝に対してのみ緩和された線形計画法を実行する代わりに、'reliability' は変数ごとに最大 k2 回まで計画法を実行します。ここで、k2 は 4 や 8 などの小さい整数です。したがって、'reliability' は 'strongpscost' に比べて、分枝に時間がかかりますが、分枝限定法による反復数は減る可能性があります。
'mostfractional' — 1/2 に最も近い小数部をもつ変数を選択します。
'maxfun' — 目的ベクトルf の絶対値に対応する最大値をもつ変数を選択します。

アルゴリズムの分枝後に、2 つの新しいノードを探索する必要があります。アルゴリズムは、以下のルールのいずれかを使用して、利用可能なすべてのノードから探索するノードを選択します。

'minobj' — 最小目的関数値をもつノードを選択します。
'mininfeas' — 実行不可能な整数の最小合計値をもつノードを選択します。これは、ノードの実行不可能なすべての整数要素 x(i) に対して、以下の p_i^– と p_i⁺ のどちらか小さい方を加算することを意味します。
p_i^– = x(i) – ⌊x(i)⌋
p_i⁺ = 1 – p_i^–.
'simplebestproj' — "最良の投影" をもつノードを選択します。
最良の投影
x_B は、これまでに検出された最良の整数実行可能点を表し、x_R は、ルートノードでの緩和された LP の解を表し、x は検証するノードを表すものとします。in(x) は、ノード x での実行不可能な整数の合計を表します ('mininfeas' を参照)。最良の投影ルールは以下を最小化します。
$f^{T} x + \frac{f^{T} x_{B} - f^{T} x_{R}}{i n (x_{R})} i n (x) .$
これまでに整数実行可能点が検出されていない場合、f^Tx_B = 0 に設定します。

intlinprog は、目的関数の最大公約数 (GCD) など、元の問題からの情報を考慮することにより、一部の部分問題の解析を省略します。

以下の停止条件の 1 つを満たすまで分枝限定手順を続行し、解析する部分問題を体系的に生成し、目的の上限または下限を改善しない部分問題を破棄します。

アルゴリズムが MaxTime オプションを超える。
目的関数の上限と下限の差が AbsoluteGapTolerance または RelativeGapTolerance 許容誤差未満である。
探索するノードの数が MaxNodes オプションを超える。
整数実行可能点の数が MaxFeasiblePoints オプションを超える。

分枝限定手順についての詳細は、Nemhauser と Wolsey [9] および Wolsey [11] を参照してください。

参照

[1] Achterberg, T., T. Koch and A. Martin. Branching rules revisited. Operations Research Letters 33, 2005, pp. 42–54. Available at https://opus4.kobv.de/opus4-zib/files/788/ZR-04-13.pdf.

[2] Andersen, E. D., and Andersen, K. D. Presolving in linear programming. Mathematical Programming 71, pp. 221–245, 1995.

[3] Atamtürk, A., G. L. Nemhauser, M. W. P. Savelsbergh. Conflict graphs in solving integer programming problems. European Journal of Operational Research 121, 2000, pp. 40–55.

[4] Berthold, T. Primal Heuristics for Mixed Integer Programs. Technischen Universität Berlin, September 2006. Available at https://www.zib.de/groetschel/students/Diplom-Berthold.pdf.

[5] Cornuéjols, G. Valid inequalities for mixed integer linear programs. Mathematical Programming B, Vol. 112, pp. 3–44, 2008.

[6] Danna, E., Rothberg, E., Le Pape, C. Exploring relaxation induced neighborhoods to improve MIP solutions. Mathematical Programming, Vol. 102, issue 1, pp. 71–90, 2005.

[7] Hendel, G. New Rounding and Propagation Heuristics for Mixed Integer Programming. Bachelor's thesis at Technische Universität Berlin, 2011. PDF available at https://opus4.kobv.de/opus4-zib/files/1332/bachelor_thesis_main.pdf.

[8] Mészáros C., and Suhl, U. H. Advanced preprocessing techniques for linear and quadratic programming. OR Spectrum, 25(4), pp. 575–595, 2003.

[9] Nemhauser, G. L. and Wolsey, L. A. Integer and Combinatorial Optimization. Wiley-Interscience, New York, 1999.

[10] Savelsbergh, M. W. P. Preprocessing and Probing Techniques for Mixed Integer Programming Problems. ORSA J. Computing, Vol. 6, No. 4, pp. 445–454, 1994.

[11] Wolsey, L. A. Integer Programming. Wiley-Interscience, New York, 1998.