ドキュメンテーション

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

線形計画法のアルゴリズム

線形計画法の定義

線形計画法は線形制約に従って、線形関数 fTx を最小化するベクトル x を見つける問題です。

minxfTx

次の 1 つ以上の式が成り立つ必要があります。

A·x ≤ b
Aeq·x = beq
l ≤ x ≤ u.

linprog 内点法アルゴリズム

linprog 'interior-point' アルゴリズムは、interior-point-convex quadprog アルゴリズム とよく似ています。また、linprog 'interior-point-legacy' アルゴリズムと多くの機能を共有しています。これらのアルゴリズムに共通する概要は次のとおりです。

  1. 解決の前処理。つまり、問題を単純化して標準的な形式に変換します。

  2. 初期点の生成。初期点の選択は、内点法アルゴリズムを効率的に解くために特に重要であり、このステップには時間がかかることがあります。

  3. KKT 方程式を解くための予測子修正子の反復。このステップには通常、最も時間がかかります。

解決の前処理

このアルゴリズムは、まず重複を取り除き制約を単純化することで問題を単純化しようと試みます。解決の前処理のステップで実行されるタスクは次のとおりです。

  • 範囲の上限と下限が等しい変数があるかどうかをチェックする。見つかった場合はその実行可能性をチェックし、修正してからそれらの変数を削除する。

  • 含まれている変数が 1 つだけの線形不等式制約があるかどうかチェックする。見つかった場合はその実行可能性をチェックし、線形制約を範囲に変更する。

  • 含まれている変数が 1 つだけの線形等式制約があるかどうかチェックする。見つかった場合はその実行可能性をチェックし、修正してからそれらの変数を削除する。

  • ゼロの行を含む線形制約行列があるかどうかチェックする。見つかった場合はその実行可能性をチェックし、それらの行を削除する。

  • 範囲制約と線形制約に矛盾がないかチェックする。

  • 目的関数内の線形項としてのみ含まれ、どの線形制約にも含まれない変数があるかどうかチェックする。見つかった場合はその実行可能性と有界性をチェックし、それらの変数を修正して適切な範囲制約を設定する。

  • 線形不等式制約があれば、スラック変数を追加して線形等式制約に変更する。

アルゴリズムは、実行不可能または非有界の問題を検出した場合、停止して適切な終了メッセージを表示します。

アルゴリズムが単一の実行可能点を得るような場合、これは解を表しています。

解決の前処理のステップで実行不可能または非有界の問題が検出されなかった場合、アルゴリズムは必要に応じて他のステップに進みます。最後に、アルゴリズムは元の問題を復元し、解決の前処理による変更を元に戻します。この最終ステップが解決の後処理ステップです。

簡単にするため、解決の前処理のステップで問題が解決しない場合、アルゴリズムはすべての有限の下限をゼロにシフトします。

初期点の生成

初期点 x0 を設定するために、アルゴリズムは次の処理を行います。

  1. x0ones(n,1) に初期化します。ここで、n は目的関数ベクトル f の要素数です。

  2. すべての有界の成分の下限が 0 になるように変換します。成分 i に有限の上限 u(i) がある場合は、x0(i) = u/2 になります。

  3. 要素の範囲制約が 1 つだけの場合は、必要に応じて要素に変更を加え、厳密に範囲内に収まるようにします。

  4. x0 が中央パスに近くなるようにするには、予測子修正子のステップを 1 回実行して、結果の位置およびスラック変数を変更して範囲内に十分に収まるようにします。中央パスの詳細は、Nocedal と Wright [6] (p. 397) を参照してください。

予測子修正子

fmincon内点法アルゴリズムと同様に、interior-point-convex アルゴリズムも カルーシュ・キューン・タッカー (KKT) 条件が満たされる点を見つけようとします。線形計画問題についてこれらの方程式を記述するために、前処理後の線形計画問題の標準的な形式を考えます。

minxfTx subject to {A¯x=b¯x+t=ux,t0.

  • ここでは、すべての変数に少なくとも 1 つの有限範囲があると仮定します。この仮定は、必要に応じて成分のシフトおよび符号反転を行うと、すべての x 成分の下限が 0 になることを意味します。

  • A¯ は、線形不等式および線形等式の両方を含む拡張された線形行列です。b¯ は、対応する線形等式ベクトルです。A¯ には、不等式制約を等式制約にするスラック変数 s でベクトル x を拡張するための項も含まれています。

    A¯x=(Aeq0AI)(x0s),

    ここで、x0 は元のベクトル x を意味します。

  • t は、上限を等式に変換するスラックのベクトルです。

このシステムのラグランジュ関数には、次のベクトルが含まれます。

  • y: 線形等式に関連付けられたラグランジュ乗数

  • v: 下限に関連付けられたラグランジュ乗数 (正値性制約)

  • w: 上限に関連付けられたラグランジュ乗数

ラグランジュ関数は次のようになります。

L=fTxyT(A¯xb¯)vTxwT(uxt).

そのため、このシステムの KKT 条件は次のようになります。

fA¯Tyv+w=0A¯x=b¯x+t=uvixi=0witi=0(x,v,w,t)0.

アルゴリズムは最初にニュートン・ラフソン式からステップを予測し、その後、修正子ステップを計算します。修正子は非線形相補性方程式 sizi = 0 の残差を減らそうとします。ニュートン・ラフソン ステップは次のようになります。

(0A¯T0IIA¯0000I0I00V00X000W0T)(ΔxΔyΔtΔvΔw)=(fA¯Tyv+wA¯xb¯uxtVXWT)=(rdrprubrvxrwt),(8-1)

ここで、X、V、W、および T は、ベクトル x、v、w、および t にそれぞれ対応する対角行列です。方程式の最も右辺にある残差ベクトルは以下のとおりです。

  • rd: 双対残差

  • rp: 主残差

  • rub: 上限残差

  • rvx: 下限相補性残差

  • rwt: 上限相補性残差

式 8-1 を解くには、まず対称行列形式に変換します。

(DA¯TA¯0)(ΔxΔy)=(Rrp),(8-2)

ここで、

D=X1V+T1WR=rdX1rvx+T1rwt+T1Wrub.

D および R の定義に含まれる逆行列は、行列が対角であるため、すべて簡単に計算できます。

式 8-1 から 式 8-2 を導くには、式 8-2 の 2 行目が式 8-1 の 2 行目と同じであることに注目します。式 8-2 の 1 行目を得るには、式 8-1 の最後の 2 行を Δv および Δw について解き、その後 Δt について解きます。

式 8-2 は対称ですが、項 –D があるため正定ではありません。そのため、コレスキー分解を使用して解くことができません。さらにいくつかのステップを実行して、正定であり、コレスキー分解によって効率的に解くことができる、別の方程式を導きます。

式 8-2 の 2 行目は

A¯Δx=rp

また、1 行目は

DΔx+A¯TΔy=R.

次を代入し

Δx=D1A¯TΔy+D1R

次の結果が得られます。

A¯D1A¯TΔy=A¯D1R+rp.(8-3)

通常、ニュートン ステップを求めるための最も効率的な方法は、コレスキー分解を使用して式 8-3 を Δy について解くことです。Δy の乗数である行列は明らかに対称行列であり、退化がない場合には正定になるため、コレスキー分解の使用が可能です。その後、ニュートン ステップを求めるために、結果を元の式に代入して Δx、Δt、Δv、および Δw を求めます。ただし、A¯ が密な列をもつ場合は、代わりに式 8-2 を解く方がより効率的な場合があります。linprog 内点法アルゴリズムは、列の密度に基づいて解のアルゴリズムを選択します。

アルゴリズムについての詳細は、Mehrotra [5] を参照してください。

アルゴリズムは、訂正されたニュートン ステップを計算した後さらに計算を実行して、現在のステップを長くするとともに、後続のステップをより良いものにする準備を行うことができます。これらの複数の訂正計算により、性能とロバスト性を両方とも向上させられます。詳細については、Gondzio [3] を参照してください。

停止条件

予測子修正子アルゴリズムは、実行可能 (許容誤差範囲内に収まるようにする制約を満たす) で相対ステップ サイズが小さい点に到達するまで反復します。具体的には、次のように定義します。

ρ=max(1,A¯,f,b¯).

これらのすべての条件が満たされると、アルゴリズムが停止します。

rpρTolConrdρTolFunrcTolFun,

ここで、

rc=maxi(min(|xivi|,|xi|,|vi|),min(|tiwi|,|ti|,|wi|)).

rc は実質的に、相補性の残差 xv および tw のサイズを測定します。これらの残差はそれぞれ、解でゼロ ベクトルとなるベクトルです。

レガシ内点線形計画法

はじめに

既定のレガシ内点法は LIPSOL ([52]) をベースにしています。この方法は、主双対内点法である Mehrotra 予測子修正子アルゴリズム ([47]) のバリアントです。

メイン アルゴリズム

アルゴリズムは一連の前処理ステップ (前処理を参照) を適用することから始めます。処理後、問題は次の型をしています。

minxfTx such that {Ax=b0xu.(8-4)

上限の制約は制約行列 A に暗黙的に含まれています。主スラック変数 s を加えると、式 8-4 は以下になります。

minxfTx such that {Ax=bx+s=ux0, s0.(8-5)

これは "主問題" と考えられます。x は主変数で s は主スラック変数からなります。"双対" 問題は以下のように表せます。

maxbTyuTw  such that  {ATyw+z=fz0, w0,(8-6)

ここで、y と w は双対変数で、 z は双対スラック変数です。この線形計画法に対する最適性、すなわち主式 8-5 および双対式 8-6 は以下になります。

F(x,y,z,s,w)=(Axbx+suATyw+zfxizisiwi)=0,                 x0, z0, s0, w0,(8-7)

ここで、xizi と siwi は要素単位の乗算を表します。

二次方程式 xizi = 0 および siwi = 0 は線形計画法に対する "相補性条件" と呼ばれ、他の (線形) 方程式は "実行可能条件" と呼ばれます。次の式で計算される量は、

xTz + sTw

"双対性ギャップ" で、(x,z,s,w) ≥ 0 のとき F の相補性部分の残差を計測します。

このアルゴリズムは、"主双対アルゴリズム" です。すなわち、主計画と双対計画が同時に解かれます。これは 式 8-7 の中で線形二次システム F(x,y,z,s,w) = 0 に適用されるニュートン法に似た手法と考えられます。一方、反復の x、z、w および s は正に保たれ、そのため内点法の名前をもちます(反復は 式 8-5 の中の不等式制約により表現される厳密な意味の内点領域に存在します)。

アルゴリズムは、Mehrotra が提唱する予測子修正子アルゴリズムのバリアントです。反復 v = [x;y;z;s;w] を考えてみましょう。ここでは [x;z;s;w] > 0 です。まず、いわゆる "予測" 方向を計算します。

Δvp=(FT(v))1F(v),

これはニュートン方向になります。そして "修正子" の方向は次のようになります。

Δvc=(FT(v))1F(v+Δvp)μe^,

ここで μ > 0 は "センタリング" パラメーターと呼ばれ、注意深く選択しなければなりません。 e^  は F(v) の二次方程式に対応するものに 1 をもつ、0 と 1 を要素とするベクトルです。すなわち、すべて 2 次である相補性条件のみに摂動が適用され、すべて線形である実行可能条件には適用されません。すなわち、2 次元である相補性条件のみに摂動が適用され、線形である実行可能条件には適用されません。2 つの方向はステップ長パラメーター α > 0 と結合して新たな反復 v+ を得るために v が更新されます。

v+=v+α(Δvp+Δvc),

ここで、ステップ長パラメーター α は

v+ = [x+;y+;z+;s+;w+]

[x+;z+;s+;w+] > 0.

上の予測子修正子の方向の解法の中で、アルゴリズムは A·AT のコレスキー因子を変更して、(スパースの) 直接因子分解を計算します。A が密な列をもつ場合は、代わりに Sherman-Morrison 式を使います。その解が適切でない (残差が大きすぎる) 場合、拡大系形式の LDL 分解の実行によって解を検出します。(MATLAB® ldl 関数リファレンス ページの 「例 4 — D の構造体」 を参照してください)。

そしてアルゴリズムは反復が収束するまで反復します。メイン停止条件は標準的なものです。

rbmax(1,b)+rfmax(1,f)+rumax(1,u)+|fTxbTy+uTw|max(1,|fTx|,|bTyuTw|)tol,

ここで、

rb=Axbrf=ATyw+zfru=x+su

は、主残差、双対残差、実行可能上限です。

fTxbTy+uTw

は、主目的値と双対目的値の間の差で tol は許容誤差です。停止条件の和は式 8-7 の最適化条件で相対的な残差の総和を計測します。

前処理

このアルゴリズムは、まず重複を取り除き制約を単純化することで問題を単純化しようと試みます。解決の前処理のステップで実行されるタスクは次のとおりです。

  • 範囲の上限と下限が等しい変数があるかどうかをチェックする。見つかった場合はその実行可能性をチェックし、修正してからそれらの変数を削除する。

  • 含まれている変数が 1 つだけの線形不等式制約があるかどうかチェックする。見つかった場合はその実行可能性をチェックし、線形制約を範囲に変更する。

  • 含まれている変数が 1 つだけの線形等式制約があるかどうかチェックする。見つかった場合はその実行可能性をチェックし、修正してからそれらの変数を削除する。

  • ゼロの行を含む線形制約行列があるかどうかチェックする。見つかった場合はその実行可能性をチェックし、それらの行を削除する。

  • 範囲制約と線形制約に矛盾がないかチェックする。

  • 目的関数内の線形項としてのみ含まれ、どの線形制約にも含まれない変数があるかどうかチェックする。見つかった場合はその実行可能性と有界性をチェックし、それらの変数を修正して適切な範囲制約を設定する。

  • 線形不等式制約があれば、スラック変数を追加して線形等式制約に変更する。

アルゴリズムは、実行不可能または非有界の問題を検出した場合、停止して適切な終了メッセージを表示します。

アルゴリズムが単一の実行可能点を得るような場合、これは解を表しています。

解決の前処理のステップで実行不可能または非有界の問題が検出されなかった場合、アルゴリズムは必要に応じて他のステップに進みます。最後に、アルゴリズムは元の問題を復元し、解決の前処理による変更を元に戻します。この最終ステップが解決の後処理ステップです。

簡単にするため、アルゴリズムはすべての下限をゼロにシフトします。

これらの前処理ステップがアルゴリズムの反復部分をスピードアップすることができるとき、ラグランジュ乗数が必要であれば、アルゴリズムの処理中に計算される乗数は変換された問題に使用され、元のものには使用されないため、前処理ステップは未処理でなければなりません。そのため乗数が必要でない場合、この変換の逆が計算 "されず"、計算時間が節約されるかもしれません。

linprog 有効制約法アルゴリズム

中規模の有効制約法の線形計画法は fmincon (逐次二次計画法 (SQP)) が使用する逐次二次計画法の変形です。二次計画法との違いは、二次の項がゼロに設定されることです。

SQP 法の主要な各反復で、以下の形式の QP 問題を解きます。ここで Ai は m 行 n 列の行列 A の第 i 行を示します。

mindnq(d)=cTd,Aid=bi,  i=1,...,meAidbi,  i=me+1,...,m.

Optimization Toolbox™ 関数で使われる方法は、有効制約法 (射影法として知られています) ですが、[18][17] に記述されている Gill 等による方法に似ています。この方法は修正されて、LP (線形計画法)、QP (二次計画法) のどちらにも使われています。

解を求めるステップは、2 つの部分からなっています。最初の段階は (解が存在する場合) 可能領域点を計算します。2 番目の段階は解に収束する可能領域点の列を生成します。このメソッドでは、有効制約法 A¯k が、求解点でアクティブな制約条件 (すなわち制約境界上にある制約条件) の推定値となるように維持されます。事実上、すべての QP アルゴリズムは有効制約法です。構造的には非常に類似しているにもかかわらず、異なる用語で記述されている多くの方法があるので、この点を強調しておきます。

A¯k は各反復 k で更新され、探索方向 d^k の基底を作成するために使われます。等式制約は常に有効制約法 A¯k 内にあります。変数 d^k の表記法は、ここでは SQP 法の主要反復の dk と区別するために使用されています。探索方向 d^k が計算され、任意の有効制約境界上に存在したまま目的関数が最小化されます。d^k の部分可能空間は、基底 Zk から作成されます。この基底の列は有効制約法 A¯k の推定値に直交します (すなわち、A¯kZk=0)。これにより、Zk の列結合の線形和から形成される探索方向は、有効制約法の境界上にあることが保証されます。

行列 Zk は行列 A¯kT を QR 分解した行列の最後の m – l 列から形成されます。ここで l はアクティブな制約の数であり、l < m です。つまり、Zk は次のように与えられます。

Zk=Q[:,l+1:m],(8-8)

ここで、

QTA¯kT=[R0].

Zk が見つかると、q(d) を最小にする新しい探索方向 d^k が求められます。ここで、d^k はアクティブな制約条件のヌル空間です。つまり、d^k は Zk の列の線形結合で、あるベクトル p に対して d^k=Zkp となります。

d^k を置き換え、2 次多項式を p の関数として見ると次のようになります。

q(p)=12pTZkTHZkp+cTZkp.(8-9)

p で微分すると、次の結果が得られます。

q(p)=ZkTHZkp+ZkTc.(8-10)

∇q(p) は、 Zk が定義する部分空間の射影勾配なので、二次関数の射影勾配になります。ZkTHZk は射影ヘッシアンと呼ばれます。ヘッセ行列 H は、(SQP のこの方法の場合) 正定であると仮定しているので、Zk が定義する部分空間における関数 q(p) の最小値は ∇q(p) = 0 で生じます。これは次の線形方程式系の解の場合に成り立ちます。

ZkTHZkp=ZkTc.(8-11)

ステップは次の式で計算されます。

xk+1=xk+αd^k,  where d^k=Zkp.(8-12)

各反復で目的関数の二次式の性質により、ステップの長さ α の選択には 2 つの方法があります。d^k に沿った単位長さのステップは、A¯k のヌル空間に制限された関数の最小値に向かう厳密なステップになります。制約に違反せずに、このようなステップをとることができれば、これは QP の解 (式 8-9) になります。そうでない場合、最近傍の制約に向かう d^k に沿ったステップが 1 単位より小さくなり、新しい制約は次の反復で有効制約法に含まれます。任意方向 d^k の制約境界までの距離は次式で与えられます。

α=mini{1,...,m}{(Aixkbi)Aid^k},(8-13)

これは有効制約法内にない制約に対して定義され、方向 d^k は制約の境界に向かっています。すなわち、Aid^k>0, i=1,...,m です。

n 個の独立した制約が有効制約法内に含まれ、最小値の位置を含んでいない場合、ラグランジュ乗数 λk は次の線形方程式の特異点にならないように計算されます。

A¯kTλk=c.(8-14)

λk のすべての要素が正の場合、xk は QP の最適解です (式 8-9を参照)。しかし λk に負の要素があり、等式制約に対応していない場合、その要素は有効制約法から削除され、新しい反復が行われます。

初期化

アルゴリズムは実行可能な初期点を必要とします。SQP 法における現在の点が実行可能でない場合は、次の線形計画法を解いて点を求めることができます。

minγ, xnγ  such thatAix=bi,      i=1,...,meAixγbi,  i=me+1,...,m.(8-15)

表記 Ai は行列 A の i 番目の行です。式 8-15 式での実行可能点 (存在する場合) は等式制約を満たす値に x を設定して求めることができます。これは、等式制約の組から作成される過決定または、劣決定のどちらかの線形方程式を解くことにより得られます。この問題の解が存在する場合、スラック変数 γ はこの点で最大不等式の制約になります。

LP 問題の上記 QP アルゴリズムは、各反復で探索方向を最急降下方向に設定することにより変更できます。ここで、gk は目的関数の勾配 (線形目的関数の係数に等しい) です。

d^k=ZkZkTgk.(8-16)

LP 法を使って実行可能な点が求められたならば、メインの QP 部分に移行します。探索方向 d^k は、線形連立方程式を解いて求まる探索方向 d^1 を使って初期化されます。

Hd^1=gk,(8-17)

ここで gk は現在の反復 xk での目的関数の勾配です (すなわち Hxk + c です)。

QP 問題で実行可能解が見つからない場合、メインの SQP ルーチンの探索方向 が、γ を最小化する方向として用いられます。

linprog シンプレックス アルゴリズム

1947年に George Dantzigにより発明されたシンプレックス アルゴリズムは、最も初期の有名な最適化アルゴリズムの 1 つです。このアルゴリズムは、次の線形計画法を解きます。

minxfTx such that {Axb,Aeqx=beq,lbxub.

このアルゴリズムは制約により定義された多面体の端に沿って、ある頂点から他の頂点まで移動し、各ステップで目的関数 fTx が減少します。この節では、頂点での最適解を返すオリジナルのシンプレックス アルゴリズムの改良されたバージョンについて述べます。

本節では、次のトピックを説明します。

メイン アルゴリズム

シンプレックス アルゴリズムには 2 つの段階があります。

  • 段階 1 — 初期の実行可能な基底点を計算します。

  • 段階 2 — 元の問題の最適解を計算します。

    メモ:   シンプレックス アルゴリズムでは linprog に対して初期点 x0 を提供することはできません。x0 を入力引数として渡す場合、linprogx0 を無視し、アルゴリズムの初期点を計算します。

段階 1-  このアルゴリズムは補助的な区分的線形計画法を解くことにより、初期の実行可能基底解 (定義は 基底変数と非基底変数 を参照) を見つけます。補助問題の目的関数は "線形ペナルティ関数" P=jPj(xj) です。

ここで、Pj(xj) は以下のように定義されます。

Pj(xj)={xjujif xj>uj0if ljxjujljxjif lj>xj.

P(x) は、点 x が上限と下限の条件をどれだけ超えているかを測ります。補助問題は以下になります。

minxjPj  subject to {AxbAeqx=beq.

オリジナルの問題は、補助問題が最小値 0をもつ場合に限り、実行可能な基底点をもちます。

アルゴリズムは、必要に応じてスラック変数と人為変数を追加するというヒューリスティックな方法により補助問題に対する初期点を見つけます。すると、アルゴリズムは、補助問題を解くために、この初期点をシンプレックス アルゴリズムと共に使用します。この解は、メイン アルゴリズムの段階 2 に対する初期点です。

段階 2-  アルゴリズムが元の問題を解くために、段階 1 からの初期点から始め、シンプレックス アルゴリズムを適用します。各反復で、アルゴリズムは、最適化の条件をテストし、現在の解が最適である場合、終了します。現在の解が最適でない場合、アルゴリズムは、次のことを行います。

  1. 非基底変数から "entering 変数" と呼ばれる 1 変数を選択し、基底でない対応する列を追加します (定義は 基底変数と非基底変数 を参照)。

  2. 基底変数から "leaving 変数" と呼ばれる変数を選択し、基底から対応する列を除きます。

  3. 現在の解と現在の目的関数値を更新

アルゴリズムは、解の実行可能性を保ちながら、2 つの線形システムを解くことにより、entering 変数と leaving 変数を選択します。

アルゴリズムは、段階 2 の解法プロセスで進歩がないことを検出します。そして、範囲の摂動を実行することによって続行しようとします。アルゴリズムのこの部分については、Applegate、Bixby、Chvatal、Cook の[59]を参照してください。

前処理

シンプレックス アルゴリズムでは、内点線形計画法のソルバーと同じ前処理ステップ (前処理で説明) を使用します。さらに、アルゴリズムは、他の 2 つのステップを使用します。

  1. 1 つの非零要素のみをもつ列と、それに対応する行を除きます。

  2. 各等式制約 a·x = b に対して a を Aeq の行として上限と下限が有限の場合、アルゴリズムは線形結合 a·x の上限と下限を rlb と rub として計算します。rlb または rub のいずれかが b に等しい場合、この制約は "強制制約" と呼ばれます。アルゴリズムは強制制約に依存して a·x の非ゼロの係数に対応する各変数を上限または下限に等しいものと設定します。その後、アルゴリズムはこれらの変数に対応する列を削除し、forcing constraint に対応する行を削除します。

シンプレックス アルゴリズムの使用

シンプレックス法を使用するには optimoptions を使用して [アルゴリズム] オプションを 'simplex' に設定します。

options = optimoptions(@linprog,'Algorithm','simplex');

その後、入力引数 options を指定して linprog を呼び出します。詳細は、linprog のリファレンス ページを参照してください。

linprog は前処理の手続きで実行不可能性や非有界性を検出する場合、xfval に対して空の出力引数を返します。linprog は以下の場合、現在の点を出力します。

  • 最大反復回数を超える

  • 段階 1 または 2 において、問題が実行不可能であるか、非有界であることを検出

問題が非有界である場合、linprog は非有界な方向で xfval を出力します。

基底変数と非基底変数

この節では線形計画法に対して用語として基底、非基底、および "実行可能基底解" を定義します。この定義では、問題が次の標準的な形式で与えられることを仮定します。

minxfTx such that {Ax=b,lbxub.

(A および b は元の問題で不等式を定義する行列やベクトルではないことに注意してください)。A がランク m < n の m 行 n 列の行列で、その列が {a1, a2, ..., an} であると仮定します。{ai1,ai2,...,aim} がインデックス集合 B = {i1, i2, ..., im} として A の列空間に対する基底であり、N = {1, 2, ..., n}、\B は B の補数であると仮定します。部分行列 AB"基底" と呼ばれ、相補性部分行列 AN"非基底" と呼ばれます。"基底変数" のベクトルは xB であり、"非基底変数" のベクトルは xN です。段階 2 の各反復でアルゴリズムは現在の基底の 1 列を非基底の列で置き換え、それに応じて変数 xB と xN を更新します。

x が A·x = b の解であり、xN の非基底変数すべてがそれらの下限または上限に等しいとき、x は "基底解" と呼ばれます。さらに xB の基底変数がそれらの上限と下限を満たし、したがって x が実行可能な点である場合、x は "基底実行可能解" と呼ばれます。

双対シンプレックス法アルゴリズム

上位レベルでは、linprog 'dual-simplex' アルゴリズムは実質的に、シンプレックス法アルゴリズムを "双対問題" に対して実行します。

このアルゴリズムは、まず前処理で説明されている前処理を行います。詳細については、Andersen and Andersen [1]、および Nocedal and Wright [6]の第 13 章を参照してください。この前処理によって、元の線形計画問題を式 8-4 の形に簡約します。

minxfTx such that {Ax=b0xu.

A と b は、元の制約行列を変換したバージョンです。これが主問題になります。

式 8-6 で説明されているように、双対問題とは、次を解くベクトル y と w およびスラック変数ベクトル z を求めることです。

maxbTyuTw  such that  {ATyw+z=fz0, w0.

双対問題の任意の有限解は主問題への解を与え、また主問題の任意の有限解が双対問題の解を与えることはよく知られています ([6] を参照してください)。さらに、主問題と双対問題のどちらかが非有界である場合、もう一方の問題は実行不可能です。そして、主問題または双対問題のどちらかが実行不可能な場合には、もう一方の問題は実行不可能であるか、非有界であるかのどちらかです。したがって、この 2 つの問題は、解が存在する場合には、有限解を求める点で等価です。主問題と双対問題は数学的に等価ではあるものの計算手順は異なるため、双対問題を解くことによって主問題を解く方が良い場合もあります。

退化を緩和するために (Nocedal and Wright [6] の 366 ページを参照)、双対シンプレックス法アルゴリズムは最初に目的関数に摂動を加えます。

双対シンプレックス法アルゴリズムの第 1 段階は、双対実行可能点を見つけることです。アルゴリズムは、シンプレックス法アルゴリズムの 段階 1 と同様に、補助線形計画問題を解くことによってこれを求めます。

第 2 段階では、ソルバーが entering 変数と leaving 変数を繰り返し選択します。これは、主シンプレックス アルゴリズムの 段階 2 のものと似ています。アルゴリズムは、Forrest と Goldfarb [2] により提唱された双対最急エッジ プライシング (dual steepest-edge pricing) という手法に従って leaving 変数を選択します。また、アルゴリズムは Koberstein [4] により提唱された Harris の比率判定法のバリエーションを使用して entering 変数を選択します。退化を緩和するため、アルゴリズムによって第 2 段階で追加の摂動が加えられる場合もあります。

ソルバーは、主問題の実行不可能性を低減し、かつ双対問題の実行可能性を維持できるような方法で反復を行います。反復は、低減され摂動された問題について主問題と双対問題の両方が実行可能であるような解が見つかるまで続けられます。加えられた摂動はアルゴリズムによって取り除かれます。摂動のない元の問題において、摂動を加えた問題の解では双対問題が実行不可能な場合、ソルバーは主シンプレックス法または第 1 段階のアルゴリズムを使用して双対問題の実行可能性を復元します。最後に、ソルバーが前処理の手順を元に戻し、解を元の問題に返します。

参照

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

[2] Forrest, J. J., and D. Goldfarb. Steepest-edge simplex algorithms for linear programming. Math. Programming 57, 1992, pp. 341–374.

[3] Gondzio, J. “Multiple centrality corrections in a primal dual method for linear programming.” Computational Optimization and Applications, Volume 6, Number 2, 1996, pp. 137–156. Available at http://www.maths.ed.ac.uk/~gondzio/software/correctors.ps.

[4] Koberstein, A. Progress in the dual simplex algorithm for solving large scale LP problems: techniques for a fast and stable implementation. Computational Optim. and Application 41, 2008, pp. 185–204.

[5] Mehrotra, S. “On the Implementation of a Primal-Dual Interior Point Method.” SIAM Journal on Optimization, Vol. 2, 1992, pp 575–601.

[6] Nocedal, J., and S. J. Wright. Numerical Optimization, Second Edition. Springer Series in Operations Research, Springer-Verlag, 2006.

この情報は役に立ちましたか?