ドキュメンテーション

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

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

線形計画法の定義

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

minxfTx

条件は以下のとおりです。 A·x ≤ b, Aeq·x = beq, l ≤ x ≤ u.

内点線形計画法

はじめに

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

メイン アルゴリズム

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

minxfTx such that {Ax=b0xu.(6-63)

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

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

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

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

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

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

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

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

xTz + sTw

"双対性ギャップ" で、(x,z,s,w) ≥ 0 のとき F の相補性部分の残差の尺度になります。

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

アルゴリズムは、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 つの方向はステップ長パラメーター α > 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 は許容誤差です。停止条件の和は 式 6-66 の最適化条件で相対的な残差の総和の尺度になります。

前処理

実際の反復アルゴリズムが開始される前にいくつかの前処理ステップが実行されます。結果の変換された問題は、次のようなものです。

  • すべての変数が 0 以下

  • すべての制約が等式制約

  • 固定された変数の除去。これらは、等しい上限と下限をもつものです。

  • 制約行列の中のすべてのゼロ点の行が除去されます。

  • 制約行列は構造的にフルランクです。

  • 制約行列の中のすべてのゼロ点の列を除去します。

  • 制約行列に単項行が多数存在する場合、関連した変数は除去される行に対して解かれます。

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

有効制約法 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],(6-67)

ここで、

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.(6-68)

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

q(p)=ZkTHZkp+ZkTc.(6-69)

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

ZkTHZkp=ZkTc.(6-70)

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

xk+1=xk+αd^k,  where d^k=ZkTp.(6-71)

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

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

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

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

A¯kTλk=c.(6-73)

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

初期化

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

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

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

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

d^k=ZkZkTgk.(6-75)

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

Hd^1=gk,(6-76)

ここで 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 variable と呼ばれる 1 変数を選択し、基底でない対応する列を追加します (定義は 基底変数と非基底変数 を参照)。

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

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

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

アルゴリズムは、段階 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 は "基底実行可能解" と呼ばれます。

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