制約付き非線形最適化アルゴリズム

制約付き最適化の定義

制約付きの最小化は、スカラー関数 f(x) の局所的最小値のベクトル x を見つける問題です。この x には制約があります。

$$\underset{x}{\min }f(x)$$

次の 1 つ以上の式が成り立つ必要があります。c(x) ≤ 0, ceq(x) = 0, A·x ≤ b, Aeq·x = beq, l ≤ x ≤ u。半無限計画法で使用される制約もあります。fseminf の問題の定式化とアルゴリズムを参照してください。

fmincon の信頼領域 Reflective 法アルゴリズム

非線形最小化に対する信頼領域法

Optimization Toolbox™ のソルバーで使用される多くのメソッドは、"信頼領域法" を基にしています。信頼領域法はシンプルなものですが最適化では重要な概念です。

最適化の信頼領域法のアプローチを理解するために制約なし最小化問題を考え、f(x) を最小化します。ここで関数はベクトル引数を取り、スカラーを出力します。n 空間の点 x を想定し、より小さい関数値へ移動して最小化を行う場合を考えてみましょう。基本的な概念は、シンプルな関数 q で f を近似することです。この関数は、点 x の近傍 N で関数 f の挙動をよく表すものです。この近傍が信頼領域です。テスト ステップ s が N における最小化 (または近似最小化) によって計算されます。信頼領域の部分問題を次に示します。

$$\underset{s}{\min }\left\{q(s),s\in N\right\}.$$

(1)

f(x + s) < f(x) の場合、現在の点が x + s に更新されます。そうでない場合は現在の点は変更されず、信頼領域 N は縮小され、テスト ステップの計算が繰り返されます。

f(x) を最小化するための信頼領域を決める上で重要な問題は、近似 q (現在の点 x で定義) の選択および計算方法、信頼領域 N の選択および変更方法、信頼領域の部分問題を解く精度です。この節では制約なしの問題に焦点をあてます。変数上に制約があるために複雑度が増す問題については、後節で取り扱います。

標準的な信頼領域法 ([48]) では二次近似 q は、x における F のテイラー近似のはじめの 2 項によって決められます。近傍 N は球形または楕円体です。数学的に、信頼領域の部分問題は、次のように表現できます。

$$\min \left\{\frac{1}{2}{s}^{T}Hs+s^{T}g\text{ such that }\Vert Ds\Vert \le \Delta \right\},$$

(2)

ここで g は現在の点 x における f の勾配です。H はヘッセ行列 (2 次導関数の対称行列) です。D は対角スケーリング行列であり、Δ は正のスカラーです。‖ . ‖ は 2 ノルムです。式 2を解くために適したアルゴリズムがあります ([48] を参照)。このようなアルゴリズムは、H のすべての固有値と以下の永年方程式に適用されるニュートン過程の計算を一般的に含みます。

$$\frac{1}{\Delta }-\frac{1}{\Vert s\Vert }=0.$$

このようなアルゴリズムにより 式 2 の正確な解が与えられます。しかし、H の因子分解に比例して時間が必要になります。そのために、大規模問題に対しては別のアプローチが必要となります。式 2 に基づく近似とヒューリスティックな方法は文献 ([42] と [50]) で提案されています。Optimization Toolbox のソルバーがサポートする近似アプローチとして、信頼領域法の部分問題を 2 次元の部分空間 S ([39] と [42]) に制限します。部分空間 S が計算されると、(部分空間では問題は 2 次元になるため) 完全な次元の固有値 / 固有ベクトルの情報が必要な場合でも 式 2 を解く作業は簡単になります。ここでの主な仕事は、部分空間を決定することになります。

2 次元の部分空間 S は、以下に示す前処理付き共役勾配法過程を用いて決められます。ソルバーは s₁ と s₂ で決まる線形空間として S を定義します。ここで、s₁ は勾配 g の方向にあります。s₂ は近似ニュートン方向、すなわち以下の解

$$H\cdot s_2=-g,$$

(3)

または負の曲率の方向です。

$$s_2^T\cdot H\cdot s_2<0.$$

(4)

このように S を選択する背景の考え方は、最急降下方向または負の曲率方向にグローバルな収束を進めることと、ニュートン ステップが存在する場合は、これを介して迅速にローカルな収束を達成することです。

信頼領域法の概念を使用した制約なしの最小化の概要は以下になります。

1. 2 次元の信頼領域の部分問題の定式化

2. テスト ステップ s を決めるため、式 2 を解きます。

3. f(x + s) < f(x) の場合、x = x + s とします。

4. Δ を調節します。

これら 4 つのステップは、収束するまで繰り返されます。信頼領域の大きさ Δ は標準的な規則に基づいて調整されます。特に、使用するステップが適用されない場合、すなわち f(x + s) ≥ f(x) の場合、内点集合は小さくなります。詳細は [46] と [49] を参照してください。

Optimization Toolbox のソルバーは特定の関数 f の重要かつ特別なケースをいくつか扱います。非線形最小二乗法、二次関数、線形最小二乗法を考えてみましょう。しかし、根底に存在するアルゴリズムは、一般的な場合と同じです。これらの特別な場合は後続の節で説明します。

前処理付き共役勾配法

線形方程式系 Hp = –g の大きな対称正定値システムを解く一般的な方法は、前処理付き共役勾配法 (PCG) です。この反復法は、形式 H·v の行列ベクトル積を計算する機能を必要とします。ここで v は任意のベクトルです。対称正定値行列 M は H の "前提条件子" です。すなわち、M = C² です。ここで、C^–1HC^–1 は、条件数の良い行列または分類分けされた固有値をもつ行列です。

最小化の過程で、ヘッセ行列 H が対称と仮定します。しかし、H は、強い最小化子の近傍の中でのみ正定値であることが保証されます。PCG アルゴリズムは、負または 0 の曲率方向が検出された場合に終了します (d^THd ≤ 0)。PCG 出力方向 p は、負の曲率方向またはニュートン システム Hp = –g への近似解のどちらかです。いずれにせよ、p は信頼領域法のアプローチで使用される 2 次元部分空間を定義するために役立ちます (非線形最小化に対する信頼領域法を参照)。

線形等式制約

線形制約は、制約なしの最小化問題に対して記述された状況を複雑にします。しかし、前述の基本的な考え方は、明確かつ効率的な形で踏襲されます。Optimization Toolbox のソルバーの信頼領域法は、厳密に実行可能な反復処理を行います。

一般的な線形等号制約付き最小化問題は、次のように表現されます。

$$\min \left\{f(x)\text{ such that }Ax=b\right\},$$

(5)

ここで、A は m 行 n 列 (m ≤ n) の行列です。いくつかの Optimization Toolbox のソルバーは A を前処理し、A^T の LU 分解をベースにした手法 (参考文献 [46]) を使って、厳密な意味で線形従属の部分を除去します。ここで A はランク m とします。

式 5 を解くためのメソッドは制約なしアプローチと 2 つの点で異なります。まず初期実行可能点 x₀ は、スパースの最小二乗ステップを使って Ax₀ = b となるように計算されます。第 2 に、近似退化ニュートン ステップ (または A のヌル空間における負の曲率方向) を計算するために、PCG アルゴリズムは RPCG (退化前処理付き共役勾配法) に置き換えられます ([46] を参照)。キーになる線形代数ステップは、次の型のシステムを解くことです。

$$\left[\begin{array}{cc}C& {\tilde{A}}^T\\ \tilde{A}& 0\end{array}\right]\left[\begin{array}{c}s\\ t\end{array}\right]=\left[\begin{array}{c}r\\ 0\end{array}\right],$$

(6)

ここで、$$\tilde{A}$$ は A を近似 (ランクを失わないという条件で A の小さな非ゼロがゼロに設定される) し、C は H へのスパース対称正定値近似、すなわち C = H です。詳細は [46] を参照してください。

ボックス制約

ボックス制約問題は、次のように表現されます。

$$\min \left\{f(x)\text{ such that }l\le x\le u\right\},$$

(7)

ここで、l は下限を表すベクトル、u は上限を表すベクトルです。l の要素のいくつか (またはすべて) は –∞ にすることができ、u の要素のいくつか (またはすべて) は ∞ にすることができます。この方法は厳密な意味で実行可能点の列を生成します。ロバストな収束挙動を達成しながら、可能領域を維持するために 2 つの手法が使われます。1 つ目の手法ではスケーリングされた変更ニュートン ステップが (2 次元の部分空間 S を定義するために) 制約のないニュートン ステップと置き換わります。2 つ目の手法では反射がステップ サイズを増加させるために使われます。

スケーリングされた変更ニュートン ステップは 式 7 式のキューン・タッカーの必要条件をチェックして生成されます。

$${\left(D(x)\right)}^{-2}g=0,$$

(8)

ここで、

$$D(x)=\text{diag}\left({\left|v_k\right|}^{-1/2}\right),$$

ベクトル v(x) は 1 ≤ i ≤ n の範囲で次のように定義されます。

• g_i < 0 および u_i < ∞ の場合、v_i = x_i – u_i とします。

• g_i ≥ 0 および l_i > –∞ の場合、v_i = x_i – l_i とします。

• g_i < 0 および u_i = ∞ の場合、v_i = –1 とします。

• g_i ≥ 0 および l_i = –∞ の場合、v_i = 1 とします。

非線形システム 式 8 には微分不可能な点があります。v_i = 0 のとき微分不可能になります。厳密に実行可能性を維持することにより、このような点を避けることができます。すなわち、l < x < u の制約を適用します。

式 8 により与えられる非線形方程式系に対する、スケーリングされた変更ニュートン ステップ s_k は、

$$\widehat{M}D{s}^N=-\widehat{g}$$

(9)

の第 k 反復における線形システムの解として定義されています。ここで

$$\widehat{g}=D^{-1}g=\text{diag}\left({\left|v\right|}^{1/2}\right)g,$$

(10)

および以下となります。

$$\widehat{M}=D^{-1}H{D}^{-1}+\text{diag}(g)J^v.$$

(11)

ここで J^v は |v| のヤコビアンの役割をします。対角行列 J^v の個々の対角要素は 0、–1、1 のいずれかに等しくなります。l と u のすべての要素が有限ならば、J^v = diag(sign(g)) になります。g_i = 0 の点では、v_i は微分可能にならない可能性があります。そのような点では、$$J_{ii}^v=0$$ が定義されます。このような微分不可能性は、v_i がどの値を取るかは重要でないため、問題とはなりません。さらに、|v_i| はこの点でもまだ不連続となりますが、関数 |v_i|·g_i は連続です。

2 つ目の手法では反射がステップ サイズを増加させるために使われます。(単一) 反射ステップは、次のように定義されます。範囲制約に交差するステップ p を与え、p に交差する最初の範囲制約を与えます。すなわち、これを i 番目の範囲制約 (i 番目の上限または i 番目の下限) であると仮定します。i 番目の要素を除いて、反射ステップは p^R = p になります。ここで p^R_i = –p_i です。

fmincon アクティブ セット アルゴリズム

はじめに

制約付き最適化における一般的な目的は、問題をより解き易い部分問題に変換し、部分問題を繰り返しの過程の基準として使用することです。初期に開発された多くの手法の特徴は、制約付き問題を、制約境界に近接または超過している制約にペナルティ関数を使って、基本的な非制約問題に変換することです。このように制約付き問題は、制約された問題に収束する（シーケンスの）範囲で制約されていないパラメータ化された最適化のシーケンスを使用して解決されます。現在、これらの方法は比較的効率が悪いと考えられており、カルーシュ・キューン・タッカー (KKT) 方程式の解に注目した方法に代行されています。KKT 方程式は制約付き最適化問題に対し、最適性に対する必要条件になっています。問題がいわゆる凸計画問題ならば、f(x) と G_i(x), i = 1,...,m は凸関数になり、KKT 方程式はグローバルな解に対して必要十分になります。

一般的な問題 GP (式 1) に関して、キューン・タッカー方程式は

$$\begin{array}{c}\nabla f\left(x^{*}\right)+{\displaystyle \sum _{i=1}^m{\lambda }_i\cdot \nabla G_i\left(x^{*}\right)}=0\\ {\lambda }_i\cdot G_i\left(x^{*}\right)=0,i=1,\mathrm{...},m_e\\ {\lambda }_i\ge 0,i=m_e+1,\mathrm{...},m,\end{array}$$

(12)

式 1 の元の制約とは別に次のように表すことができます。

最初の方程式は、解の点で目的関数とアクティブな制約条件との間で勾配がキャンセルされることを示しています。勾配をキャンセルするためには、目的関数と制約勾配との大きさの違いを平衡させるために、ラグランジュ乗数 (λ_i, i = 1,...,m) を必要とします。アクティブな制約のみがこのキャンセル操作に含まれるので、非アクティブな制約はこの操作に含まれることはないはずであり、与えられたラグランジュ乗数はゼロに等しくなります。これはキューン・タッカー式の最後の 2 つの方程式によって陰的に述べられています。

KKT の解は多くの非線形計画法アルゴリズムの基礎となっています。これらのアルゴリズムは、ラグランジュ乗数を直接計算しようとしています。制約付き準ニュートン法は、準ニュートン更新手法を使用する KKT 方程式に関して 2 次の情報を集めることにより、超線形収束を保証します。QP 部分問題がそれぞれの主な反復で解かれるので、これらの方法は一般に逐次二次計画法 (SQP: Sequential Quadratic Programming) と呼ばれています (Iterative Quadratic Programming 法、Recursive Quadratic Programming 法、Constrained Variable Metric 法としても知られています)。

'active-set' アルゴリズムは大規模なアルゴリズムではありません。大規模アルゴリズムと中規模アルゴリズム を参照してください。

逐次二次計画法 (SQP)

逐次二次計画法は、非線形計画法の中で最先端の技法です。たとえば、Schittkowski [36] は効率性、精度、解の求まる割合を改良して、多くのテスト問題について他のすべての手法と共にテストしました。

Biggs [1]、Han [22] および Powell ([32]、[33]) の成果を基に、制約付き最適化のニュートン法に非常に似た手法が制約なしの最適化に対して適用されます。各々の主反復で、準ニュートン更新法を使ってラグランジュ関数のヘッシアンで近似がなされます。そして、この手法は直線探索手法に対して、探索方向を求めるために使われる QP 部分問題を作成するのに使われます。SQP の概要は Fletcher [13]、Gill 等 [19]、Powell [35] および Schittkowski [23] に記述されています。ここでは、一般的な方法について記述しています。

GP (式 1) で問題の記述を与えると、基本的な考え方はラグランジュ関数の二次近似をもとにした QP 部分問題の定式化です。

$$L(x,\lambda )=f(x)+{\displaystyle \sum _{i=1}^m{\lambda }_i\cdot g_i(x)}.$$

(13)

ここで 式 1 は範囲付きの制約が不等式として表現されることを仮定しています。非線形制約を線形化することで、QP 部分問題を得ることができます。

QP 部分問題

$$\begin{array}{l}\underset{d\in {\Re }^n}{\min }\frac{1}{2}{d}^T{H}_{k}d+\nabla f{\left(x_k\right)}^{T}d\\ \nabla g_i{\left(x_k\right)}^{T}d+g_i\left(x_k\right)=0,i=1,\mathrm{...},m_e\\ \nabla g_i{\left(x_k\right)}^{T}d+g_i\left(x_k\right)\le 0,i=m_e+1,\mathrm{...},m.\end{array}$$

(14)

この部分問題は、どの QP アルゴリズムでも解くことができます (例は、二次計画法での解法 を参照)。解は次の新しい反復を作るために使われます。

x_{k + 1} = x_k + α_kd_k.

ステップ長パラメーター α_k はメリット関数内で十分な減少が得られるように適当な直線探索手法により決められます (ヘッセ行列の更新 を参照)。行列 H_k はラグランジュ関数 (式 13) のヘッセ行列の正定値近似になります。BFGS 法が最もポピュラーですが、準ニュートン法のいずれかを使っても H_k は更新できます (ヘッセ行列の更新を参照)。

非線形の制約付き問題は、SQP 法を使った制約なし問題よりも少ない反復回数で解けることもあります。この理由の一つは、可能領域上に限定されているためで、オプティマイザーが探索方向とステップ長に関して情報に基づく決定を行っていることが一因しています。

Rosenbrock 関数に付加的な非線形不等式制約 g(x) を適用することを考慮します。

$$x_1^2+x_2^2-1.5\le 0.$$

(15)

制約なしの場合は 140 回の反復であるのに対して、SQP 法を使って解くと 96 回の反復になります。図 5-3、非線形制約付き Rosenbrock 関数における SQP 法は x = [–1.9,2.0] から始め、x = [0.9072,0.8228] の解の点までの経路を示しています。

図 5-3、非線形制約付き Rosenbrock 関数における SQP 法

SQP 法の実装

SQP 法は、次の小節で簡単に論議されている 3 つの部分から作られています。

• ヘッセ行列の更新

• 二次計画法での解法

• 初期化

• 直線探索とメリット関数

ヘッセ行列の更新.  主要な各反復でラグランジュ関数のヘッシアンの正定値準ニュートン近似 H は、λ_i, i = 1,...,m をラグランジュ乗数の推定とする BFGS 法を使って計算されます。

$$H_{k+1}=H_k+\frac{q_k{q}_k^T}{q_k^T{s}_k}-\frac{H_k{s}_k{s}_k^T{H}_k^T}{s_k^T{H}_k{s}_k},$$

(16)

ここで、

$$\begin{array}{c}{s}_k=x_{k+1}-x_k\\ q_k=\left(\nabla f\left(x_{k+1}\right)+{\displaystyle \sum _{i=1}^m{\lambda }_i\cdot \nabla g_i\left(x_{k+1}\right)}\right)-\left(\nabla f\left(x_k\right)+{\displaystyle \sum _{i=1}^m{\lambda }_i\cdot \nabla g_i\left(x_k\right)}\right).\end{array}$$

Powell [33] は解の点で正定値でない部分があったとしても正定値ヘッシアンを使い続けることを推奨しています。正定値ヘッシアンは、$$q_k^T{s}_k$$ が各更新で正であり、H が正定値行列で初期化された場合に保持されます。$$q_k^T{s}_k$$ が正でない場合、q_k は $$q_k^T{s}_k>0$$ になるように各要素ごとに変更されます。この変更の一般的な目的は、q_k の要素に歪みを与え、正定値更新にできるだけ小さく関与するようにすることです。そのため、変更の初期ステージで q_k*s_k の最小の負要素を半分にする作業を反復します。この手続きは $$q_k^T{s}_k$$ が小さな負の許容値以上になるまで続けられます。この処理後でも $$q_k^T{s}_k$$ がまだ正になっていなければ、定数スカラー w を乗じたベクトル v を加え、q_k を変更します。

$$q_k=q_k+wv,$$

(17)

ここで、

$$\begin{array}{l}{v}_i=\nabla g_i\left(x_{k+1}\right)\cdot g_i\left(x_{k+1}\right)-\nabla g_i\left(x_k\right)\cdot g_i\left(x_k\right)\\ \text{ if }{\left(q_k\right)}_i\cdot w<0\text{ and }{\left(q_k\right)}_i\cdot {\left(s_k\right)}_i<0,i=1,\mathrm{...},m\\ v_i=0\text{ otherwise,}\end{array}$$

w は $$q_k^T{s}_k$$ が正になるまでシステマチックに増加します。

関数 fmincon、fminimax、fgoalattain、fseminf は、すべて SQP 法を使います。options の Display を 'iter' に設定すると、関数値、制約違反の最大量のような種々の情報が与えられます。ヘッシアンの正定値を維持するために、先の手順の最初の方法を使って変更する必要がある場合、Hessian modified が表示されます。前のアプローチの 2 番目の方法を使って再度ヘッシアンを変更する必要がある場合は、Hessian modified twice が表示されます。QP 部分問題が実行不可能な場合、infeasible が表示されます。このような表示は、原因を求めるものではなく、問題が非常に高い非線形性であるとか、通常より収束に長い時間が必要であるということを示すものです。場合によってはメッセージ no update が表示され、$$q_k^T{s}_k$$ がほぼゼロであることが示されます。これは、問題設定が間違っているか、または非連続関数を最小化しようとしている可能性を示しています。

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

$$\begin{array}{c}\underset{d\in {\Re }^n}{\min }q(d)=\frac{1}{2}{d}^{T}Hd+c^{T}d,\\ A_{i}d=b_i,i=1,\mathrm{...},m_e\\ A_{i}d\le b_i,i=m_e+1,\mathrm{...},m.\end{array}$$

(18)

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

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

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

行列 Z_k は行列 $${\overline{A}}_k^T$$ を QR 分解した行列の最後の m – l 列から形成されます。ここで l はアクティブな制約の数であり、l < m です。つまり、Z_k は次のように与えられます。

$$Z_k=Q\left[:,l+1:n\right],$$

(19)

ここで、

$$Q^T{\overline{A}}_k^T=\left[\begin{array}{c}R\\ 0\end{array}\right].$$

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

$${\widehat{d}}_k$$ を置き換え、二次式を p の関数として見ると次のようになります。

$$q(p)=\frac{1}{2}{p}^T{Z}_k^{T}H{Z}_{k}p+c^T{Z}_{k}p.$$

(20)

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

$$\nabla q(p)=Z_k^{T}H{Z}_{k}p+Z_k^{T}c.$$

(21)

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

$$Z_k^{T}H{Z}_{k}p=-Z_k^{T}c.$$

(22)

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

$$x_{k+1}=x_k+\alpha {\widehat{d}}_k,\text{ where }{\widehat{d}}_k=Z_{k}p.$$

(23)

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

$$\alpha =\underset{i\in \left\{1,\mathrm{...},m\right\}}{\min }\left\{\frac{-\left(A_i{x}_k-b_i\right)}{A_i{\widehat{d}}_k}\right\},$$

(24)

これはアクティブ セットにない制約に対して定義され、方向 $${\widehat{d}}_k$$ は制約の境界に向かっています。すなわち、$$A_i{\widehat{d}}_k>0,i=1,\mathrm{...},m$$ です。

n 個の独立した制約がアクティブ セットに含まれ、最小値の位置を含んでいない場合、ラグランジュ乗数 λ_k は次の線形方程式の特異点にならないように計算されます。

$${\overline{A}}_k^T{\lambda }_k=c+H{x}_k.$$

(25)

λ_k のすべての要素が正の場合、x_k は QP の最適解です (式 18を参照)。しかし λ_k に負の要素があり、等式制約に対応していない場合、その要素はアクティブ セットから削除され、新しい反復が行われます｡

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

$$\begin{array}{l}\underset{\gamma \in \Re ,x\in {\Re }^n}{\min }\gamma \text{ such that}\\ A_{i}x=b_i,i=1,\mathrm{...},m_e\\ A_{i}x-\gamma \le b_i,i=m_e+1,\mathrm{...},m.\end{array}$$

(26)

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

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

$${\widehat{d}}_k=-Z_k{Z}_k^T{g}_k.$$

(27)

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

$$H{\widehat{d}}_1=-g_k,$$

(28)

ここで g_k は現在の反復 x_k での目的関数の勾配です (すなわち Hx_k + c です)。

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

直線探索とメリット関数.  QP 部分問題の解はベクトル d_k を生成し、次の新しい反復を作成するために使用されます。

$$x_{k+1}=x_k+\alpha d_k.$$

(29)

ステップ長パラメーター α_k はメリット関数を十分小さくするように決定されます。この実装では、Han [22] および Powell [33] により導入された、次の型のメリット関数が使用されます。

$$\Psi (x)=f(x)+{\displaystyle \sum _{i=1}^{m_e}{r}_i\cdot g_i(x)}+{\displaystyle \sum _{i=m_e+1}^m{r}_i\cdot \max [0,g_i(x)]}.$$

(30)

Powell は、ペナルティ パラメーターを設定することを推奨しています。

$$r_i={\left(r_{k+1}\right)}_i=\underset{i}{\max }\left\{{\lambda }_i,\frac{{\left(r_k\right)}_i+{\lambda }_i}{2}\right\},i=1,\mathrm{...},m.$$

(31)

これは、QP 解において前の反復では active であったが現在 inactive となっている制約からの正の寄与を可能にするものです。これを実装するにあたり、ペナルティ パラメーター r_i は最初に次のように設定されています。

$$r_i=\frac{\Vert \nabla f(x)\Vert }{\Vert \nabla g_i(x)\Vert },$$

(32)

ここで $$\Vert \Vert$$ はユークリッド ノルムを示します。

これは、解の点において制約が active になる場合に、ペナルティ パラメーターへの小さな勾配をもつ制約からの大きな寄与を補償します。

fmincon SQP アルゴリズム

sqp アルゴリズム (およびほぼ同一の sqp-legacy アルゴリズム) は active-set アルゴリズムと同様です (詳細についてはfmincon アクティブ セット アルゴリズムを参照)。基本的な sqp アルゴリズムは、Nocedal と Wright の [31] の第 18 章に説明されています。

sqp アルゴリズムは実質的に sqp-legacy アルゴリズムと同一ですが、異なる実装をもちます。通常、sqp は sqp-legacy より実行時間が速く、メモリ使用量は少なくなります。

sqp と active-set アルゴリズムの最も重要な差異は、以下のとおりです。

範囲に関する厳密な実行可能性

sqp アルゴリズムは、範囲によって制約されている領域内ですべての反復ステップを取ります。さらに、有限差分ステップも範囲に従います。範囲は厳密ではありません。ステップは境界上に正確に置くことができます。この厳密な実行可能性は、目的関数または非線形制約関数が未定義または範囲によって制約される領域外で複雑であるときに便利です。

非倍精度結果に対するロバスト性

反復時に、sqp アルゴリズムは失敗するステップを取るよう試みることができます。こうすることにより、与えた目的関数または非線形制約関数は NaN、Inf、または複雑な値を返します。この場合、アルゴリズムはいまよりも小さなステップを取るよう試みます。

リファクタリングされた線形代数ルーチン

sqp アルゴリズムは異なるセットの線形代数ルーチンを使用して、二次計画法部分問題 式 14 を解きます。これらのルーチンは、active-set のルーチンよりも、メモリ使用量およびスピードの両面においてより効率的です。

再定式化された実行可能性ルーチン

sqp アルゴリズムは、制約が満たされないときの 式 14 の解に対して、2 つの新しい方法をもっています。

• sqp アルゴリズムは目的関数と制約関数を 1 つのメリット関数に結合します。このアルゴリズムは、緩和された制約に従ってメリット関数を最小化しようと試みます。変更されたこの問題により実行可能解が得られる場合があります。ただし、この方法は元の問題よりも多くの変数をもつので、式 14における問題のサイズは増加します。サイズが増加されたことにより、部分問題の解を得るためにより多くの時間を必要としています。これらのルーチンは、Spellucci の [60] と Tone の [61] の記事に基づいています。sqp アルゴリズムでは、メリット関数 式 30 のペナルティ パラメーターが [41] の提案に従って設定されます。

• 非線形制約が満たされず、試みたステップにより制約違反が増大するという状況を考えてみます。sqp アルゴリズムは、制約に対して二次近似を使用して実行可能性を得るように試みます。この二次近似の手法により実行可能解が得られる場合があります。ただし、この手法は非線形制約関数の評価をさらに多く必要とするので、解を得るためにより多くの時間を必要としています。

fmincon の内点法アルゴリズム

バリア関数

制約付き最小化の内点法アプローチとは、一連の近似最小化問題を解くことです。元の問題は次のようになります。

$$\underset{x}{\min }f(x),\text{ subject to }h(x)=0\text{ and }g(x)\le 0.$$

(33)

0 より大きい各 μ に対し、近似問題は次のようになります。

$$\begin{array}{l}\underset{x,s}{\min }{f}_{\mu }(x,s)=\underset{x,s}{\min }f(x)-\mu {\displaystyle \sum _i\ln \left(s_i\right)},\text{ subject to }s\ge 0,h(x)=0,\text{ and }g(x)+s=0.\\ \end{array}$$

(34)

不等式制約 g と同じだけのスラック変数 s_i があります。反復解を実行可能領域の内部にとどめるために、s_i は正に制限されます。μ がゼロに減少すると、f_μ は f の最小値に到達します。追加された対数項はバリア関数と呼ばれます。このメソッドは [40]、[41]、[51] に記述されています。

近似問題 式 34 は等式制約付き問題になります。これは元の非線形制約付き問題 式 33 を解くより簡単です。

近似問題を解くために、各反復でアルゴリズムは 2 種類のステップ タイプの中からいずれかを使用します。

• (x、s) の直接ステップ。このステップは線形近似の近似問題に対して KKT 式の 式 2 と 式 3 を解きます。これはニュートン ステップとも呼ばれます。

• 信頼領域法を使用する "CG" (共役勾配法) ステップ。

既定の設定では、このアルゴリズムはまず直接ステップを試します。直接ステップが使えない場合、CG ステップを試します。直接ステップが使えない場合として、近似問題が現在の反復で局所的に凸でない場合があります。

反復の都度、このアルゴリズムは "メリット関数" を以下のように減少させます。

$$f_{\mu }(x,s)+\nu \Vert \left(h(x),g(x)+s\right)\Vert .$$

(35)

パラメーター $$\nu$$ は解が可能領域に入るようにするため、反復数と共に増加します。実行したステップがメリット関数を減少させない場合、アルゴリズムは実行したステップを棄却して新しいステップを試します。

反復 x_j において、目的関数または非線形制約関数のどちらかが複素数値、NaN、Inf、またはエラーを返すと、アルゴリズムは x_j を棄却します。棄却は、メリット関数が十分に減少しなかった場合と同じ効果があります。次に、アルゴリズムは別のより短いステップを試みます。try-catch でエラーを発生する可能性のあるコードをすべてラップします。

function val = userFcn(x)
try
    val = ... % code that can error
catch
    val = NaN;
end

目的関数と制約は初期点で適切な値 (double) を生成しなければなりません。

直接ステップ

以下の変数は直接ステップの定義に使用されます。

• H は f_μ のラグランジュ関数のヘッシアンを示します。

  $$H={\nabla }^{2}f(x)+{\displaystyle \sum _i{\lambda }_i{\nabla }^2{g}_i(x)}+{\displaystyle \sum _j{y}_j{\nabla }^2{h}_j(x)}.$$

  (36)

• J_g は制約関数 g のヤコビアンを示します。

• J_h は制約関数 h のヤコビアンを示します。

• S = diag(s).

• λ は次の制約に関連するラグランジュ乗数ベクトルを示します。 g

• Λ = diag(λ).

• y は h に関連するラグランジュ乗数ベクトルを示します。

• e は g と同じサイズの 1 のベクトルを示します。

式 38は直接ステップ (Δx, Δs) を定義します。

$$\left[\begin{array}{cccc}H& 0& J_h^T& J_g^T\\ 0& \Lambda & 0& S\\ J_h& 0& 0& 0\\ J_g& I& 0& 0\end{array}\right]\left[\begin{array}{c}\Delta x\\ \Delta s\\ \Delta y\\ \Delta \lambda \end{array}\right]=-\left[\begin{array}{c}\nabla f+J_h^{T}y+J_g^T\lambda \\ S\lambda -\mu e\\ h\\ g+s\end{array}\right].$$

(37)

この方程式は線形化されたラグランジュ関数を使用して式 2と式 3を解こうとすることで直接得られます。

2 つ目の変数 Δs を S^-1 で事前に乗算しておくことで、方程式を対称化できます。

$$\left[\begin{array}{cccc}H& 0& J_h^T& J_g^T\\ 0& S\Lambda & 0& S\\ J_h& 0& 0& 0\\ J_g& S& 0& 0\end{array}\right]\left[\begin{array}{c}\Delta x\\ S^{-1}\Delta s\\ \Delta y\\ \Delta \lambda \end{array}\right]=-\left[\begin{array}{c}\nabla f(x)+J_h^{T}y+J_g^T\lambda \\ S\lambda -\mu e\\ h\\ g(x)+s\end{array}\right].$$

(38)

(Δx, Δs) に対してこの式を解くために、このアルゴリズムは行列の LDL 分解を行います。(MATLAB^® ldl 関数リファレンス ページの 例 3 — D の構造体 を参照。) これは最も計算に時間がかかるステップです。この因子分解の 1 つの結果は、予測されたヘッシアンが正定値であるかどうかの判定基準になります。正定値でない場合、このアルゴリズムは共役勾配ステップに記述されている共役勾配ステップを使用します。

範囲パラメーターの更新

近似問題式 34を元の問題に接近させる場合、反復の続行に伴って範囲パラメーター μ が 0 に向かって減少する必要があります。このアルゴリズムには 2 つの範囲パラメーター更新オプションがあります。これらは BarrierParamUpdate オプションの 'monotone' (既定) または 'predictor-corrector' を使用して指定します。

'monotone' オプションは、近似問題が直前の反復で十分な精度で解かれると、パラメーター μ を 1/100 または 1/5 に減少させます。1 回または 2 回の反復だけで十分な精度を達成した場合は 1/100、それ以外の場合は 1/5 が使用されます。精度の測定は次のテスト (式 38の右側のすべての項のサイズが μ 未満であるかどうかを決定する) によります。

$$\max \left(\Vert \nabla f(x)+J_h^{T}y+J_g^T\lambda \Vert ,\Vert S\lambda -\mu e\Vert ,\Vert h\Vert ,\Vert g(x)+s\Vert \right)<\mu .$$

メモ

次のいずれかの場合、fmincon は BarrierParamUpdate 設定を 'monotone' にオーバーライドします。

• 問題に等式制約 (範囲制約を含む) がない場合。

• SubproblemAlgorithm オプションが 'cg' の場合。

範囲パラメーター μ を更新するための 'predictor-corrector' アルゴリズムは、線形計画法の予測子修正子アルゴリズムに似ています。

予測子修正子ステップは、ニュートン ステップにおける線形化エラーを補正することにより、既存のフィアッコ・マコーミック (モノトーン) アプローチを加速できます。予測子修正子アルゴリズムには 2 つの効果があります。それは、ステップの方向を頻繁に改善すると同時に、範囲パラメーターをセンタリング パラメーター σ と適応的に更新することで、中央パスに沿った反復を促すというものです。なぜ中央パスによってステップ サイズが大きくなり、その結果として収束が速くなるのかについては、線形計画法用の予測子修正子ステップについて論じた Nocedal と Wright [31] を参照してください。

予測子ステップは、μ =0 (つまりバリア関数なし) で線形化されたステップを使用します。

$$\left[\begin{array}{cccc}H& 0& J_h^T& J_g^T\\ 0& \Lambda & 0& S\\ J_h& 0& 0& 0\\ J_g& I& 0& 0\end{array}\right]\left[\begin{array}{c}\Delta x\\ \Delta s\\ \Delta y\\ \Delta \lambda \end{array}\right]=-\left[\begin{array}{c}\nabla f+J_h^{T}y+J_g^T\lambda \\ S\lambda \\ h\\ g+s\end{array}\right].$$

ɑ_s と ɑ_λ が、非負制約に違反しない最も大きなステップ サイズとなるように定義します。

$$\begin{array}{l}{\alpha }_s=\max \left(\alpha \in (0,1]:s+\alpha \Delta s\ge 0\right)\\ {\alpha }_{\lambda }=\max \left(\alpha \in (0,1]:\lambda +\alpha \Delta \lambda \ge 0\right).\end{array}$$

次に、予測子ステップから相補性を計算します。

$${\mu }_P=\frac{\left(s+{\alpha }_s\Delta s\right)\left(\lambda +{\alpha }_{\lambda }\Delta \lambda \right)}{m},$$

(39)

ここで、m は制約数です。

最初の修正子ステップは、ニュートンの根を求めるための線形化では無視された 2 次項を補正します。

$$\left(s+\Delta s\right)\left(\lambda +\text{}\Delta \lambda \right)=\underset{\text{Linear term set to 0}}{\underbrace{s\lambda +\text{}s\Delta \lambda +\lambda \Delta s}}+\underset{\text{Quadratic}}{\underbrace{\Delta s\Delta \lambda }}.$$

2 次エラーを修正するには、修正子ステップ方向用の線形システムを解きます。

$$\left[\begin{array}{cccc}H& 0& J_h^T& J_g^T\\ 0& \Lambda & 0& S\\ J_h& 0& 0& 0\\ J_g& I& 0& 0\end{array}\right]\left[\begin{array}{c}\Delta x_{\text{cor}}\\ \Delta s_{\text{cor}}\\ \Delta y_{\text{cor}}\\ \Delta {\lambda }_{\text{cor}}\end{array}\right]=-\left[\begin{array}{c}0\\ \Delta s\Delta \lambda \\ 0\\ 0\end{array}\right].$$

2 つ目の修正子ステップは、中心化ステップです。中心化補正は、方程式の右側にある変数 σ に基づきます。

$$\left[\begin{array}{cccc}H& 0& J_h^T& J_g^T\\ 0& \Lambda & 0& S\\ J_h& 0& 0& 0\\ J_g& I& 0& 0\end{array}\right]\left[\begin{array}{c}\Delta x_{\text{cen}}\\ \Delta s_{\text{cen}}\\ \Delta y_{\text{cen}}\\ \Delta {\lambda }_{\text{cen}}\end{array}\right]=-\left[\begin{array}{c}\nabla f+J_h^{T}y+J_g^T\lambda \\ S\lambda -\mu e\sigma \\ h\\ g+s\end{array}\right].$$

ここで、σ は次のように定義されます。

$$\sigma ={\left(\frac{{\mu }_P}{\mu }\right)}^3,$$

ここで、μ_P は方程式式 39で定義され、$$\mu =s^T\lambda /m$$ です。

範囲パラメーターの減少速度が極端に速くなるとアルゴリズムが不安定になる可能性があるため、アルゴリズムのセンタリング パラメーター σ は 1/100 を下回らないようにします。このアクションにより、範囲パラメーター μ の減少を 1/100 以下に抑えられます。

アルゴリズム上、最初の修正子ステップと中心化ステップは互いに独立しているため、これらのステップは一緒に計算されます。また、予測子ステップについても、2 つの修正子ステップについても、左側の行列の内容は同じです。したがって、アルゴリズム上、行列は一度因子分解されているため、その因子分解はこれらのステップすべてに使用されていることになります。

アルゴリズムは、推奨されている予測子修正子ステップがメリット関数の値 式 35 を増加 (相補性を 2 倍以上増加) させる場合、または計算された慣性が間違っている (非凸問題に見える) 場合、このステップを棄却することがあります。このような場合、アルゴリズムは違うステップまたは共役勾配ステップを取るように試みます。

共役勾配ステップ

近似問題 式 34 を解く共役勾配アプローチは他の共役勾配計算と同様です。この場合、このアルゴリズムはx と s を調整し、スラック s を正に保持します。このアプローチは線形化された制約の下で信頼領域法の近似問題の二次近似を最小化します。

特に R は信頼領域の半径を示すようにし、他の変数を 直接ステップ で定義されているようにします。このアルゴリズムは λ を正の状態にしたまま最小二乗的に以下の KKT 式を近似して解き、ラグランジュ乗数を得ます。

$${\nabla }_{x}L={\nabla }_{x}f(x)+{\displaystyle \sum _i{\lambda }_i\nabla g_i(x)}+{\displaystyle \sum _j{y}_j\nabla h_j(x)}=0,$$

そしてステップ (Δx, Δs) をとり、以下を近似して解きます。

$$\underset{\Delta x,\Delta s}{\min }\nabla f^T\Delta x+\frac{1}{2}\Delta x^T{\nabla }_{xx}^{2}L\Delta x+\mu e^T{S}^{-1}\Delta s+\frac{1}{2}\Delta s^T{S}^{-1}\Lambda \Delta s,$$

(40)

これは、次の線形制約に従います。

$$g(x)+J_g\Delta x+\Delta s=0,h(x)+J_h\Delta x=0.$$

(41)

式 41を解くために、アルゴリズムは R でスケーリングされた半径の範囲内で線形化された制約のノルムを最小化しようとします。そして制約を使って 式 40 を解きます。この制約は半径 R の信頼領域内で s を厳密に正の状態にし、式 41 を解いた残差を適合させます。アルゴリズムと微分の詳細は [40]、[41]、[51] を参照してください。共役勾配の他の説明は前処理付き共役勾配法を参照してください。

実行可能性モード

EnableFeasibilityMode オプションが true のとき、反復しても実行不可能性が十分に早く減少しない場合、アルゴリズムは実行可能性モードに切り替わります。この切り替えは、ノーマル モードでは実行不可能性が減少せず、共役勾配モードへの切り替えにも失敗した場合に発生します。したがって、ソルバーが実行可能性モードを使わずに実行可能解を見つけられない場合に最良のパフォーマンスを得るには、SubproblemAlgorithm を 'cg' に設定して実行可能モードを使用します。これを行うことで、ノーマル モードで無駄な検索が行われなくなります。

実行可能性モードのアルゴリズムは、Nocedal、Öztoprak、および Waltz [1] に基づいています。このアルゴリズムは目的関数を無視しますが、代わりに、不等式制約関数の正の部分および等式制約関数の絶対値の総計として定義された実行不可能性を最小化しようと試みます。それぞれが不等式、等式の正の部分、等式の負の部分に対応する緩和変数 $$r=\left[r_I,r_e^{+},r_e^{-}\right]$$ に関して、この問題は次のようになります。

$$\underset{x,r}{\min }{1}^{T}r=\underset{x,r}{\min }{\displaystyle \sum r}$$

以下の制約に従います。

$$\begin{array}{c}{r}_I\ge c(x)\\ r_e^{+}-r_e^{-}=ceq(x)\\ r\ge 0.\end{array}$$

この緩和された問題を解く場合、ソフトウェアは対数のバリア関数による内点法定式化とスラック $$s=\left[s_I,s_R,s_e^{+},s_e^{-}\right]$$ を使用して最小化を行います。

$$\underset{x,r,s}{\min }{1}^{T}r-\mu {\displaystyle \sum \left(\log \left(s_{I,i}\right)+\log \left(s_{R,i}\right)\right)}-{\displaystyle \sum \left(\log \left(s_{e,i}^{+}\right)+\log \left(s_{e,i}^{-}\right)\right)}$$

以下の制約に従います。

$$\begin{array}{c}ceq(x)=r_e^{+}-r_e^{-}\\ c(x)=r_I-s_I\\ r_e^{+}=s_e^{+}\\ r_e^{-}=s_e^{-}\\ r_I=s_R\\ s\ge 0.\end{array}$$

緩和された問題の解法プロセスは、現在の範囲パラメーター値に初期化された μ から開始します。スラック変数 s_I は、メイン モードから継承された、現在の不等式スラック値に初期化されます。変数 r は次のように初期化されます。

$$\begin{array}{c}{r}_I=\max \left(s_I+c(x),\mu \right)\\ r_e^{+}=\max \left(ceq(x),\mu \right)\\ r_e^{-}=-\min \left(ceq(x),-\mu \right).\end{array}$$

残りのスラック変数は次のように初期化されます。

$$\begin{array}{l}{s}_R=r_I\\ s_e^{+}=r_e^{+}\\ s_e^{-}=r_e^{-}.\end{array}$$

実行可能性モードのアルゴリズムは、この初期点から始まりますが、通常の内点法アルゴリズムのコードを再利用します。変数 r は線形であり、したがってそれらに関連する 2 階微分はゼロになるため、このプロセスでは特別なステップ計算が必要です。つまり、実行可能性問題の目的関数のヘッシアンでは、ランク落ちします。したがって、このアルゴリズムはニュートン ステップを取ることができません。代わりに、最急降下方向ステップを取ります。アルゴリズムはこれらの変数について目的関数の勾配を最初に求め、制約のヤコビアンのヌル空間に勾配を投影し、生成されたベクトルのステップ長が適切になるようにベクトルを再スケーリングします。このステップは、実行不可能性の減少に役立つ可能性があります。

実行可能性モードのアルゴリズムは、実行不可能性が 1/10 に減少すると終了します。実行可能性モードが終了すると、アルゴリズムは変数 x と s_I をメイン アルゴリズムに渡し、他のスラック変数と緩和変数 r を破棄します。

参照

[1] Nocedal, Jorge, Figen Öztoprak, and Richard A. Waltz. An Interior Point Method for Nonlinear Programming with Infeasibility Detection Capabilities. Optimization Methods & Software 29(4), July 2014, pp. 837–854.

内点法のアルゴリズムのオプション

ここでは、内点法アルゴリズムのいくつかのオプションの効果と意味を説明します。

• HonorBounds — true に設定すると、各反復は設定した範囲制約を満たします。false に設定すると、アルゴリズムは中間の反復中に範囲違反します。

• HessianApproximation — 以下に設定された場合

  • 'bfgs': fmincon が密な準ニュートン近似によってヘッシアンを計算します。

  • 'lbfgs': fmincon が制限されたメモリの大規模の準ニュートン近似によってヘッシアンを計算します。

  • 'fin-diff-grads': fmincon が勾配の有限差分によってヘッシアンとベクトルの積を計算します。他のオプションが適切に設定される必要があります。

• HessianFcn — fmincon が HessianFcn で指定された関数ハンドルを使用して、ヘッシアンを計算します。詳細は、ヘッシアンを含めるを参照してください。

• HessianMultiplyFcn — ヘッシアンとベクトルの積を評価する別の関数を指定します。詳細については、ヘッシアンを含めるとヘッセ乗算関数を参照してください。

• SubproblemAlgorithm — 直接ニュートン ステップを使うかどうかを判定します。既定の設定 'factorization' はこのタイプのステップを使用します。'cg' を設定すると共役勾配ステップのみを使用します。

オプションの完全な一覧については、fmincon の options、「内点法アルゴリズム」を参照してください。

fminbnd アルゴリズム

fminbnd は MATLAB をインストールすると使用できるソルバーです。このソルバーは範囲内の 1 次元で局所的最小値を求めるソルバーです。これは導関数をベースにしません。その代わり、黄金分割法と放物線内挿法を使用します。

fseminf の問題の定式化とアルゴリズム

fseminf の問題の定式化

fseminf は fmincon とは異なるタイプの制約を使って最適化問題を解きます。fmincon の式は以下になります。

$$\underset{x}{\min }f(x)$$

c(x) ≤ 0, ceq(x) = 0, A·x ≤ b, Aeq·x = beq, and l ≤ x ≤ u が条件になります。

fseminf は上記の他に以下の半無限制約を追加します。1 次元または 2 次元の範囲空間または四角形 I_j 内の w_j と、連続関数 K(x, w) のベクトルに対して制約は以下になります。

K_j(x, w_j) ≤ 0 (すべての w_j∈I_j について)

fseminf 問題の "次元" は制約集合 I の最大次元を意味します。すべての I_j が区間ならば 1 であり、少なくとも 1 つの I_j が四角形ならば 2 になります。K のベクトル サイズはこの概念の次元を示しません。

半無限計画法と呼ばれる理由は変数 (x と w_j) の数は有限ですが、制約数が無限であるためです。これは x の制約が連続区間または四角形 I_j の集合上にあるためです。これは点数が無限にあるため、制約が無限にあります。無限個の点 w_j に対して K_j(x, w_j) ≤ 0 です。

制約が無限にある問題を解くのは不可能だと思うかもしれません。fseminf は問題を 2 つの部分からなる同様な問題に再定式することによってこのような問題を扱います。2 つの部分とは最大化と最小化です。半無限制約は以下のように再定式します。

$$\underset{w_j\in I_j}{\max }{K}_j(x,w_j)\le 0\text{ for all }j=1,\mathrm{...},\left|K\right|,$$

(42)

ここで |K| はベクトル K の要素数です。すなわち、半無限制約関数の数です。x を固定するために、これは範囲区間または四角形上の一般最大化になります。

fseminf はさらに、ソルバーが使用する x ごとに、区分的な二次または三次近似 κ_j(x, w_j) を関数 K_j(x, w_j) にすることで、問題を簡略化します。fseminf は、式 42にあるように、K_j(x, w_j) の代わりに、内挿関数 κ_j(x, w_j) の最大のみを考慮します。これは元の問題 (半無限制約関数の最小化) を有限数の制約付き問題に変形します。

サンプリング点.  半無限制約関数はサンプリング点 (二次近似または三次近似に使用する点) を与える必要があります。これを行うには以下を行わなければなりません。

• サンプリング点 w 間の初期空間 s

• s からサンプリング点 w を生成する方法

初期空間 s は |K| 行 2 列の行列です。s の j 行目は制約関数 K_j の隣接したサンプリング点間隔を示します。K_j が 1 次元の w_j に依存する場合、s(j,2) = 0 に設定します。fseminf は後続の反復で行列 s を更新します。

fseminf は行列 s を使用して、サンプリング点 w を生成します。これは近似 κ_j(x, w_j) を作成するために使用されます。s から w を生成する手順は、同じ間隔または四角形 I_j を最適化は維持する必要があります。

サンプリング点の作成例.  2 つの半無限制約 K₁ と K₂ をもつ問題を考えてみましょう。K₁ は 1 次元の w₁ をもち、K₂ は 2 次元の w₂ をもちます。以下のコードは w₁ = 2 から 12 までのサンプリング集合を生成します。

% Initial sampling interval
if isnan(s(1,1))
    s(1,1) = .2;
    s(1,2) = 0;
end

% Sampling set
w1 = 2:s(1,1):12;

fseminf は制約関数をはじめて呼び出すときに s を NaN として指定します。これを確認して、初期サンプリング間隔を設定することができます。

以下のコードは、各要素を 1 から 100 までにして正方形内の w₂ からサンプリング集合を生成します。2 番目の要素より 1 番目の要素の方が初期サンプルが多くなります。

% Initial sampling interval
if isnan(s(1,1))
   s(2,1) = 0.2;
   s(2,2) = 0.5;
end
 
% Sampling set
w2x = 1:s(2,1):100;
w2y = 1:s(2,2):100;
[wx,wy] = meshgrid(w2x,w2y);

前のコードの一部は以下のように簡略化できます。

% Initial sampling interval
if isnan(s(1,1))
    s = [0.2 0;0.2 0.5];
end

% Sampling set
w1 = 2:s(1,1):12;
w2x = 1:s(2,1):100;
w2y = 1:s(2,2):100;
[wx,wy] = meshgrid(w2x,w2y);

fseminf のアルゴリズム

fseminf は基本的に半無限計画問題を fmincon が扱う問題に変形します。fseminf は半無限計画問題を解くために以下のステップをとります。

1. x の現在値での fseminf は、内挿 κ_j(x, w_j,i) が局所的最大値となる条件において、すべての w_j,i を識別します (最大値は固定された x に対する w の変動を参照します)。

2. fseminf は以下の fmincon 問題の解において 1 反復ステップをとります。

   $$\underset{x}{\min }f(x)$$

   条件は c(x) ≤ 0, ceq(x) = 0, A·x ≤ b, Aeq·x = beq, and l ≤ x ≤ u です。ここで、c(x) は、すべての w_j∈I_j に対して取られた κ_j(x, w_j) のすべての最大値で拡張され、これは κ_j(x, w_j,i) の j と i に対する最大値と等価です。

3. fseminf は、新しい点 x が停止判定に見合うか (反復を停止するか) を確認します。停止判定に合わない場合は手順 4 に進みます。

4. fseminf は半無限制約の離散化を更新しなければならないかを確認し、近似的にサンプリング点を更新します。これにより、更新された近似 κ_j(x, w_j) が提供されます。次に手順 1 を続行します。