混合整数線形計画法 (MILP) アルゴリズム
混合整数線形計画法の定義
混合整数線形計画法 (MILP) は次を含む問題です。
線形目的関数 fTx、ここで f は定数の列ベクトル、x は未知数の列ベクトルです。
範囲および線形制約、ただし非線形制約を含みません (定義については、制約の作成を参照)。
x の一部の要素を整数値とするという制限をもちます。
数式では、与えられたベクトル f、lb および ub、行列 A および Aeq、対応するベクトル b および beq、一連の指数 intcon が解となるベクトル x を求めます。
HiGHS MILP アルゴリズム
HiGHS の概要
intlinprog アルゴリズムは HiGHS オープンソース ソフトウェアに基づいています。intlinprog は、MATLAB® 形式の入力とオプションを同等の HiGHS 引数に変換し、返された解も同様に標準の MATLAB 形式に変換します。
アルゴリズムのアウトライン
アルゴリズムで実行されるステップは次のとおりです。
分枝限定ノードを取得 (ない場合は停止)
停止条件まで繰り返し:
停止条件まで探索:
ドメインを伝播
ノードを枝刈りして範囲を更新し、実行不可能性が検出されるか
ObjectiveCutOffオプションの値に達すると終了再開条件が満たされる場合はステップ 2 に戻る
次のノードをインストール:
ノードを選択して評価
このノードが評価で枝刈りされる場合はステップ 5 に戻る
ノードのカットを生成
ドメインが実行不可能である場合はノードを切断し、開いているノードをノード キューに登録
基底を更新
ステップ 4 に移動
解決の前処理
通常、問題の変数の数 (x の要素の数) を減らし、線形制約の数を減らすことができます。これらのリダクションの実行によってソルバーに時間がかかっている間に、通常解を求める全体的な時間を減らしてより大きな問題が解決可能になります。アルゴリズムは解をより数値的に安定させることができます。さらに、これらのアルゴリズムは実行不可能な問題を検出することができる場合があります。
解決の前処理のステップは、冗長な変数と制約を除去し、制約行列のモデルとスパースのスケーリングを改善し、変数の範囲を強化し、モデルの主実行不可能性および双対実行不可能性を検出することを目的としています。背景については、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 の解に収束しようとする一連の部分問題を作成します。部分問題は、解 fTx に一連の上限と下限を設定します。これらの境界のことを主境界および双対境界と呼びます。最小化問題の場合、最初の上限 (主境界) は任意の実行可能解であり、最初の下限 (双対境界) は緩和された問題の解です。最大化問題の場合は、主境界が下限になり、双対境界が上限になります。線形計画法緩和問題の解はいずれも MILP の解より低い目的関数値をもちます。また、任意の実行可能点 xfeas は以下の式を満たします。
| fTxfeas ≥ fTx | (1) |
これは fTx がすべての実行可能点で最小値であるためです。
このコンテキストでは、"ノード" は元の問題と同じ目的関数、範囲、線形制約をもつが、整数制約をもたず、線形制約または範囲の特定の変更をもつ LP です。"ルート ノード" は整数制約および線形制約または範囲の変更をもたない元の問題です。つまり、ルート ノードは初期の緩和された LP です。
範囲の開始から、分枝限定法はルート ノードからの分枝により新しい部分問題を作成します。いくつかのルールの 1 つに従って経験則的に分枝手順が取られます。各ルールは、ある変数を整数 J 以下、または J+1 以上に制限することにより問題を分割するという考え方に基づいています。これらの 2 つの部分問題は intcon で指定された整数に対応する xLP のエントリが整数ではないときに生じます。ここで、xLP は緩和された問題の解です。J を負方向の丸め (切り捨て)、J+1 を正方向の丸め (切り上げ) として取ります。生成された 2 つの問題は、さらに制限をもつため fTxLP 以上となる解をもちます。したがって、最小化問題の場合、この手順では下限が引き上げられる可能性があります。
アルゴリズムの分枝後に、2 つの新しいノードを探索する必要があります。intlinprog は、目的関数値を現在の大域的な範囲と比較することで一部の部分問題の解析をスキップします。
以下の停止条件の 1 つを満たすまで分枝限定手順を続行し、解析する部分問題を体系的に生成し、目的の上限または下限を改善しない部分問題を破棄します。
問題が実行不可能であることが証明される。
目的関数値が
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 の解に収束しようとする一連の部分問題を作成します。部分問題は、解 fTx に一連の上限と下限を設定します。最小化問題の場合、最初の上限は任意の実行可能解であり、最初の下限は緩和された問題の解です。
分枝限定手順の実行中、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://optimization-online.org/wp-content/uploads/2015/11/5209.pdf.
[13] Wolsey, L. A. Integer Programming. Wiley-Interscience, New York, 1998.
