混合整数線形計画法 (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}$

アルゴリズムの概要

'legacy' 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 [3] および Mészáros と Suhl [10] を参照してください。

線形計画法

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

f^Tx_LP ≤ f^Tx

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

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

混合整数計画法の前処理

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

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

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

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

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

カットの生成

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

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

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

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

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

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

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

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

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

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

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

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

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

目的関数の上限を取得するには、分枝限定手順で実行可能点を検出しなければなりません。分枝限定時の 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 [8]を参照してください。`'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 [5]のセクション 3.1 を参照してください。
`'none'`	`intlinprog` は実行可能点を検索しません。ソルバーは単に分枝限定検索で検出された実行可能点を取ります。

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

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

ZI round — このヒューリスティックな方法は、緩和された LP をアルゴリズムが解くときに常に実行されます。このヒューリスティックな方法では、他の制約に関する実行可能性に影響を与えることなく、各非整数変数を隣接する整数にシフトしようと試みます。詳細については、Hendel [9]を参照してください。
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[2]を参照)。
- 候補となるすべての分枝変数 (現在は非整数だが整数にしなければならない変数) を現在の擬似コストベースのスコア順に並べます。
- 最高スコアの変数から順番に、現在の分枝変数に基づいて 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 [11]および Wolsey [13]を参照してください。

HiGHS MILP アルゴリズム

HiGHS の概要
アルゴリズムのアウトライン
解決の前処理
ルートノードの評価
カットの生成
飛び込み/ダイビング
ランダム丸め、RINS、RENS
分枝限定
反復表示

HiGHS の概要

intlinprog "highs" アルゴリズムは HiGHS オープンソースソフトウェアに基づいています。intlinprog は、MATLAB^® 形式の入力とオプションを同等の HiGHS 引数に変換し、返された解も同様に標準の MATLAB 形式に変換します。

アルゴリズムのアウトライン

"highs" アルゴリズムで実行されるステップは次のとおりです。

解決の前処理
ルートノードの評価
分枝限定ノードを取得 (ない場合は停止)
停止条件まで繰り返し:
1. 停止条件まで探索:
  1. ランダム丸め、RINS、RENSを適用
  2. 飛び込み/ダイビングを実行
2. ドメインを伝播
3. ノードを枝刈りして範囲を更新し、実行不可能性が検出されるか ObjectiveCutOff オプションの値に達すると終了
4. 再開条件が満たされる場合はステップ 2 に戻る
5. 次のノードをインストール:
  1. ノードを選択して評価
  2. このノードが評価で枝刈りされる場合はステップ 5 に戻る
  3. ノードのカットを生成
  4. ドメインが実行不可能である場合はノードを切断し、開いているノードをノードキューに登録
  5. 基底を更新
  6. ステップ 4 に移動

解決の前処理

解決の前処理のステップは、冗長な変数と制約を除去し、制約行列のモデルとスパースのスケーリングを改善し、変数の範囲を強化し、モデルの主実行不可能性および双対実行不可能性を検出することを目的としています。背景については、Andersen と Andersen [3]、および Mészáros と Suhl [10]を参照してください。

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

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

整数前処理の背景については、Achterberg らの[1]を参照してください。

ルートノードの評価

ルートノードを評価するためにアルゴリズムで実行されるステップは次のとおりです。

対称性を検出して問題を単純化する。
ルート LP を評価する。
LP カットを生成して追加する (カットの生成を参照)。
ランダム丸めを適用する。
LP カットをさらに生成して追加する。
カット生成とヒューリスティックな方法をループで実行する。
再開条件をチェックし、正当であればステップ 2 から再開する。

ルートノード評価が完了すると、アルゴリズムは分枝限定に進みます。

カットの生成

カットは、intlinprog が問題に追加する別の線形不等式制約です。これらの不等式は、解が整数に近づくように LP 緩和の実行可能領域を制限しようとします。カット生成アルゴリズム (切除平面法とも呼ばれます) の背景については、Cornuéjols [6]、および Atamtürk、Nemhauser、Savelsbergh [4]を参照してください。

飛び込み/ダイビング

整数実行可能点を見つけるために、intlinprog は分枝限定ステップに類似したヒューリスティックな方法を使用します。ただし、ツリーの 1 つの分枝にのみ進み、他の分枝は作成しません。この単一分枝によってツリーフラグメントを迅速に "ダイブ" することができるため、"ダイビング" という名前が付けられています。

ダイビングヒューリスティックは通常、整数値であり、現在の解が非整数である 1 つの変数を選択します。その後、このヒューリスティックは変数が整数値になるように強制する範囲を導入し、関連する緩和された LP を再度解きます。制限する変数の選択方法がダイビングヒューリスティック間の主な相違点です。Berthold [5]のセクション 3.1 を参照してください。

ランダム丸め、RINS、RENS

新しい整数実行可能点を見つけるために、intlinprog は現在の最適な実行可能解の整数点 (可能な場合) の近傍を検索し、新しいより適切な解を見つけます。詳細については、Danna、Rothberg および Le Pape [8]を参照してください。同じく新しい整数実行可能点を見つけるために、intlinprog は緩和された問題の LP 解をノードで取り、実行可能性を維持するように整数要素を丸めます。ランダム丸めのステップにより、intlinprog で新しい実行可能点が見つかることがあります。整数実行可能点を見つけるためのもう 1 つの探索手法に RENS (Relaxation Enforced Neighborhood Search) があります。Berthold [5]を参照してください。

分枝限定

分枝限定法は MILP の解に収束しようとする一連の部分問題を作成します。部分問題は、解 f^Tx に一連の上限と下限を設定します。これらの境界のことを主境界および双対境界と呼びます。最小化問題の場合、最初の上限 (主境界) は任意の実行可能解であり、最初の下限 (双対境界) は緩和された問題の解です。最大化問題の場合は、主境界が下限になり、双対境界が上限になります。線形計画法で説明したように、最小化問題の場合、線形計画法緩和問題の解は MILP の解より低い目的関数値をもちます。また、任意の実行可能点 x_feas は以下の式を満たします。

f^Tx_feas ≥ f^Tx

(1)

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

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

アルゴリズムの分枝後に、2 つの新しいノードを探索する必要があります。intlinprog は、目的関数値を現在の大域的な範囲と比較することで一部の部分問題の解析をスキップします。

問題が実行不可能であることが証明される。
目的関数値が ObjectiveCutOff の制限に達する。
アルゴリズムが MaxTime オプションを超える。
目的関数の上限と下限の差が AbsoluteGapTolerance または RelativeGapTolerance 許容誤差未満である。
探索するノードの数が MaxNodes オプションを超える。

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

反復表示

Display オプションを既定の "iter" に設定して反復表示を選択すると、ソルバーのステップの一部が表示されます。HiGHS の反復表示は、他のソルバーの反復表示よりも広範で複雑です。さらに、HiGHS アルゴリズムでは分枝限定探索を再開できるため、その場合は反復表示も再開されることになります。

反復表示を選択するには次のようにします。

options = optimoptions("intlinprog",Display="iter");
[x,fval,exitflag,output] = intlinprog(f,intcon,A,b,Aeq,beq,...
    lb,ub,options)

プリアンブル. 反復表示では、最初に "解決の前処理" の結果が表示されます。解決の前処理のアルゴリズムは、線形制約行列の冗長な行と列を特定して削除し、問題の関連する単純化を実行して、元の問題の複雑度を軽減します。たとえば、

Presolving model
18018 rows, 26027 cols, 248579 nonzeros
15092 rows, 24343 cols, 217277 nonzeros
Objective function is integral with scale 1

ルートノードの評価. ルートノードは、整数制約を一切考慮しない線形計画法による問題の解です。最小化問題の場合、ルートノードの目的関数値が、整数制約を含む問題に対する解の目的関数値の下限になります。最小化問題で、上限 (ある場合) は、すべての制約についての実行可能点から得られます。まだ実行可能点が見つかっていない場合、上限は Inf になります。

反復表示では、解決の前処理の後に問題のサイズが表示されます。

Solving MIP model with:
   15092 rows
   24343 cols (24343 binary, 0 integer, 0 implied int., 0 continuous)
   217277 nonzeros

binary は 2 値変数の数です。
integer は整数変数の数です。
implied int は暗黙的に整数となる変数の数です。たとえば、x(1) が整数で x(1) + x(2) = 5 であれば、x(2) は暗黙的に整数になります。
continuous は連続変数の数です。

変数の合計数がモデルの列数になります。

動的制約の作成. ソフトウェアは、最初に "動的制約" を作成し、反復表示に 3 つのヘッダーを含めます。

Cuts — アクティブなカットの数
inLp — LP 行列のモデル以外の行の数
Confl. — 競合の数

ソフトウェアは、動的制約の作成を再開して制約をさらに展開できます。またこれにより、新しい状態から作成プロセスを開始して、さらに多くの制約を作成できます。

分枝限定プロセスを開始する前の最後のステップで、ソフトウェアは対称性検出の結果、および検出されたジェネレーターとオービトープの数をレポートします。たとえば、以下はプリアンブルと動的制約の作成の段階で表示される反復表示の一部を示したものです。

        Nodes      |    B&B Tree     |            Objective Bounds              |  Dynamic Constraints |       Work      
     Proc. InQueue |  Leaves   Expl. | BestBound       BestSol              Gap |   Cuts   InLp Confl. | LpIters     Time
 
         0       0         0   0.00%   0               inf                  inf        0      0      0         0     1.5s
         0       0         0   0.00%   0               inf                  inf        0      0      5      4934     3.6s
…
         0       0         0   0.00%   0               inf                  inf     4630    553    289    140296   159.6s
 
0.0% inactive integer columns, restarting
Model after restart has 15090 rows, 24338 cols (24338 bin., 0 int., 0 impl., 0 cont.), and 217075 nonzeros
 
         0       0         0   0.00%   0               inf                  inf      550      0      0    140296   160.5s
…
         0       0         0   0.00%   0               inf                  inf     5602    524    260    277323   318.0s
 
6.2% inactive integer columns, restarting
Model after restart has 14185 rows, 22816 cols (22816 bin., 0 int., 0 impl., 0 cont.), and 200423 nonzeros
 
         0       0         0   0.00%   0               inf                  inf      524      0      0    277323   318.6s
…
         0       0         0   0.00%   0               inf                  inf     4683    535    525    408393   505.8s
 
Symmetry detection completed in 3.1s
Found 215 generators and 12 full orbitope(s) acting on 730 columns

対称性検出およびオービトープの詳細については、Hojny、Pfetsch、Schmitt [7]、および Pfetsch と Rehn [12]を参照してください。ソフトウェアは、次に説明する分枝限定のステップを進めながら動的制約の追加を続けます。

分枝限定. 分枝限定法は MILP の解に収束しようとする一連の部分問題を作成します。部分問題は、解 f^Tx に一連の上限と下限を設定します。最小化問題の場合、最初の上限は任意の実行可能解であり、最初の下限は緩和された問題の解です。

分枝限定手順の実行中、intlinprog の反復表示は次のようになります。

        Nodes      |    B&B Tree     |            Objective Bounds              |  Dynamic Constraints |       Work      
     Proc. InQueue |  Leaves   Expl. | BestBound       BestSol              Gap |   Cuts   InLp Confl. | LpIters     Time
        72       0         2   1.56%   0               inf                  inf     4695    535    738    667009   686.9s
…
 T     271     107        51   1.56%   0               449              100.00%     6051    271   1767    792786   776.8s
 T     279     107        53   1.56%   0               439              100.00%     6053    271   1794    793241   777.5s
…
 L    1223     538       295   1.98%   0               434              100.00%     6984    243   7628     1689k  1580.6s
…
      1321     539       333   1.98%   0               434              100.00%     7029    243   8628     1898k  1650.6s
 
Restarting search from the root node
Model after restart has 13947 rows, 22426 cols (22426 bin., 0 int., 0 impl., 0 cont.), and 194633 nonzeros
 
      1323       0         0   0.00%   0               434              100.00%      243      0      0     1902k  1653.5s
      1323       0         0   0.00%   0               434              100.00%      243     76     10     1905k  1655.7s
…
      1694     173        52   0.00%   0               434              100.00%     9411    318   2220     3205k  2584.3s
 B    1710     167        55   0.00%   0               433              100.00%     9415    318   2263     3207k  2586.2s
      1726     224        56   0.00%   0               433              100.00%     9999    378   2307     3237k  2608.7s

左端の列に、表示のその行について、新しい実行可能点がどのようにして見つかったかを示すコードが表示されます。

- L — 主問題のヒューリスティック中に部分 MIP 問題を解いているとき
- T — ツリー探索中にノードを評価しているとき
- B — 分岐の実行中
- H — ヒューリスティックな方法による
- P — MIP を解く前の起動中
- C — 中心丸めによる
- R — ランダム丸めによる
- S — LP を解いているとき
- F — 実行可能性ポンプによる
- U — 非有界の LP から

終了. 反復表示の最後に、アルゴリズムの停止理由と結果の概要が表示されます。

Solving report
  Status            Time limit reached
  Primal bound      14
  Dual bound        0
  Gap               100% (tolerance: 0.01%)
  Solution status   feasible
                    14 (objective)
                    0 (bound viol.)
                    1.7763568394e-15 (int. viol.)
                    0 (row viol.)
  Timing            7200.02 (total)
                    3.08 (presolve)
                    0.00 (postsolve)
  Nodes             5830
  LP iterations     9963775 (total)
                    2657448 (strong br.)
                    324908 (separation)
                    2140011 (heuristics)
 
Solver stopped prematurely. Integer feasible point found.
 
Intlinprog stopped because it exceeded the time limit, options.MaxTime = 7200 . The intcon variables are integer within tolerance,
options.ConstraintTolerance = 1e-06.

Status — 反復停止の理由
Primal bound — 最小化問題の場合は目的関数値の上限、最大化問題の場合は下限
Dual bound — 最大化問題の場合は目的関数値の下限、最小化問題の場合は上限
Gap — 主境界と双対境界の間の相対ギャップ
Solution status
- 解のステータス
- 目的関数値、ラベル (objective)
- 変数範囲に対する解の最大違反、ラベル (bound viol.)
- 整数型変数の整数値からの最大違反、ラベル (int. viol.)
- 線形制約に対する解の最大違反、ラベル (row viol.)
Timing — 解の各種段階のタイミング (秒単位)
Nodes — 探索されたノードの数
LP iterations — 線形計画法の段階ごとの反復回数と総反復回数

参照

[1] Achterberg, T.,Bixby, R., Gu, Z., Rothberg, E., and Weninger, D. Presolve Reductions in Mixed Integer Programming. INFORMS J. Computing 32(2), 2019, pp. 473–506. Available at https://doi.org/10.1287/ijoc.2018.0857.

[2] 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.

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

[4] 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.

[5] Berthold, T. Primal Heuristics for Mixed Integer Programs. Technischen Universität Berlin, September 2006. Available at https://opus4.kobv.de/opus4-zib/files/1029/Berthold_Primal_Heuristics_For_Mixed_Integer_Programs.pdf.

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

[7] Hojny, C., Pfetsch, M. E., Schmitt, A. Extended formulations for column-constrained orbitopes. TU Darmstadt, Department of Mathematics, Research Group Optimization, Dolivostr. 15, 64293 Darmstadt, June 2017. Available at https://optimization-online.org/2017/06/6092/.

[8] 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.

[9] 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.

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

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

[12] Pfetsch, M. E. and Rehn, T. A computational comparison of symmetry handling methods for mixed integer programs. Mathematical Programming Computation (2019) 11:37–93. Available at https://mpc.zib.de/archive/2019/1/Pfetsch-Rehn2019_Article_AComputationalComparisonOfSymm.pdf.

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