SQP 法の実行

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

• ラグランジュ関数のヘッセ行列の更新

• 二次計画問題の解法

• ライン探索とメリット関数の計算 ヘッセ行列の更新

各々の主繰り返しで、ラグランジュ行列関数のヘッセ行列の正定準 Newton 近似

H は、 をLagrange乗数の推定とする BFGS 法を使って計算されま

す。

-1 -0.5 0 0.5 1 1.5 2 2.5 3

-2 -1.5 -1 -0.5 0 0.5 1 1.5 2

o

Start Point

o

o

o o

o

oooooooooo

Constraint boundary Infeasible region

Feasible region

Solution

λ_i(i=1, ,… m)

制約付き最適化

ここで、 (3-26)

Powell [35]は、解の点で正定でない部分があったとしても、正定ヘッセ行列 を使

い続けることを推奨しています。 正定 ヘッセ行列 は、Hが正定行列で初期化され 各更新で が正である場合は、保持されます。 が正でない場合、 は、

になるように各要素毎に変更されます。この変更の一般的な目的は、

の要素に歪みを与え、正定更新にできるだけ小さく関与するようにすることで す。そのため、変更の初期ステージで、 の最大負要素を半分にする作業を 繰り返します。この手法は、 が 1e-5 より大きくなるまで続けます。この処

理後でも がまだ正になっていなければ、定数スカラ w を乗じたベクトル v

を加え、 を変更します。つまり、

(3-27) ここで、

これは、 および の場合、

その他の場合は、

wは、 が正になるまで、システマテックに増加します。

関数 fmincon, fminimax, fgoalattainおよびfseminfは、すべて、SQP 法を使 います。options パラメータの Display を 'iter' に設定すると、関数値、制約 違反の最大量のような種々の情報が与えられます。 ヘッセ行列が、正定条件を満 たしたまま、上で記述した手順の最初の方法を使って変更されるとき、情報とし て Hessian modified が表示されます。 ヘッセ行列が 2 番目の方法を使って変更 されるとき、Hessian modified twiceが表示されます。 QP サブ問題が可能領域 にない場合、infeasibleが表示されます。 このような表示は、原因を求めるもの ではなく、問題が非常に高い非線形性であるとか、通常より収束に長い時間が必 要であるということを示すものです。 時々、メッセージ no update は、 が

H_k+1 H_k qkq

k T

qk Ts

k

--- H_k^TH_k sk

TH

ks

k

---– +

=

s_k = x_k+1–x_k

q_k ∇f x( _k+1) λ_i⋅∇g_i(x_k+1)

i=1 n

∑

^{∇f x( )k λi⋅∇gi( )xk}

i=1 n

∑

+

 

 

 

 

– +

=

q_k^Ts_k q_k^Ts_k q_k

q_k^Ts_k>0 q_k

q_k.∗s·_k q_k^Ts_k

q_k^Ts_k q_k q_k = q_k+wv

v_i = ∇g_i(x_k+1)⋅g_i(x_k+1)–∇g_i( )x_k ⋅g_i( )x_k , q_k

( )_i⋅w<0 ( )q_{k i}⋅( )s_{k i}<0(i= 1,…m) v_i = 0

q_k^Ts_k

qk Ts

k

ほぼゼロであることを意味します。 これは、問題設定が間違っているか、または 非連続関数を最小化しようとしている可能性を示しています。

二次計画法での解法

SQP法の各主繰り返しで、QP問題はm行n列の行列 のi番目の行を とし、そ の型を使って解かれます。

(3-28)

Optimization Toolbox の中で使われている方法は、参考文献 [20]、[19] に記述され ている Gill et al .の方法に似ている active set strategy ( Projection 法として知られて いる ) です。 この方法は、LP 問題、QP 問題に対しても変更して使われています。

解を求める手順は、二つの部分からなっています。 最初のステップは、解が存在 するなら可能領域の中の点を計算するもので、つぎのステップは、解に収束する 可能領域の中にある点の列を形成することです。 この手法でactive set は、解を 計算する点でのactiveな制約条件(すなわち、制約の境界上での条件)の推定値 であるように維持されます。 事実上、すべての QP アルゴリズムは active set 法で

す。 構造的には非常に類似しているにもかかわらず、異なる用語で記述されてい

る多くの方法があるので、この点を強調しておきます。

は各繰り返しkで更新され、探索方向に対する基底を作成するために使われま

す。 等式制約は、常に active set の中に留まります。 変数 に対する記法は、

SQP 法の主繰り返しでの とは別なものとして使われています。探索方向 が計算され、任意のactive制約境界上に存在したまま、目的関数が最小化されま す。 に対する可能サブ空間は、基底 から作られます。この基底の列は、

active set の推定値に直交するものです (すなわち、 )。 そして の列

の任意結合に関する線形和から作られる探索方向は、active 制約の境界上にある ことが保証されます。

行列 は、行列 のQR分解の後半 列から作られます。ここで、l は

active制約の数で、l < m です。 は、つぎのように与えられます。

(3-29) ここで、

A A_i

q d( ) d∈ℜⁿ

minimize 1

2---d^THd+c^Td

=

A_id = b_i i = 1, ,… m_e A_id≤b_i i = m_e+1, ,… m

A_k

A_k

A_k dˆ_k

d_k dˆ_k

dˆ_k Z_k

A_k A_kZ_k = 0 Z_k

Z_k A_k^T m–l

Z_k Z_k = Q : l[ , + :m1 ]

制約付き最適化

いったん が見つかれば、新しい探索方向 は の最小を求めます。こ

こで、 は、active制約のヌル空間です。 つまり、 は、あるベクトル pに対し

て : の線形結合列です。

に対する代入により、２次型を p の関数として見るなら、つぎのようになり

ます。

(3-30)

p で微分すると、つぎの結果を得ます。

(3-31) は、 により定義されるサブ空間に射影される勾配なので、２次関数の射 影勾配になります。 は射影されたヘッセ行列(projected Hessian)と呼ばれ

ます。 ヘッセ行列H は、(SQPのこの方法の場合)正定であると仮定しているので、

により定義されるサブ空間での関数 q(p) の最小値は、 で生じます。

この条件は連立線形方程式の解です。

(3-32)

ステップは、つぎの式で計算されます。

(3-33) 各々の繰り返しで、目的関数の二次型の性質のために、ステップの長さ の選択 には二つの方法があります。 に沿った単位大きさのステップは、 のヌル空 間に制限された関数の最小値に向かう厳密な意味でのステップになります。 制約 条件を満たしたまま、そのようなステップをとることができると、これは、QPの

解(（3-29） 式)になります。 その他の場合、最近傍の制約に向かう に沿った

ステップが1より小さくなり、新しい制約はつぎの繰り返しで、acrive setの中に 含まれます。 任意の方向 について、制約の境界までの距離は次式で与えられま す。

(3-34) Q^TA_k^T R

0

=

Z_k dˆ_k q d( )

dˆ_k dˆ_k

Z_k dˆ_k = Z_kp dˆ_k

q p( ) 1 2---p^TZ

k THZ

kp c^TZ

kp +

=

q p( )

∇ Z

k THZ

kp Z

k Tc +

= q p( )

∇ Z_k

Z_k^THZ_k

Z_k ∇q p( ) = 0

Zk THZ

kp Z

k Tc –

=

xk+1 x

k+αdˆ_k

= where dˆ_k Z

k Tp

=

α

dˆ_k A_k

dˆ_k dˆ_k

α –(A_ix_k–b_i) A_idˆ_k

--- 

 

 

mini

= (i=1, ,… m)

これは、active set の中にない制約条件に対する定義です。そして、方向 は、

制約の境界に向いています。 すなわち、 です。 n 個の独立な 制約条件がactive setの中に含まれ、最小値の位置を含んでいないとき、Lagrange 乗数 がつぎの線形方程式の特異点でない組を満足させるように計算されます。

(3-35) すべての要素 が正ならば、 がQP(（3-29） 式)の最適解になります。 しかし、

のある要素が負ならば、そして、等式での制約に対応していないならば、対応 する要素はactive setから削除され、新しい繰り返しを探します。

初期化. アルゴリズムは、実行可能な開始点を要求します。 SQP 法において現在 の点が実行可能な領域に入っていないならば、点は、つぎの線形計画問題を解く ことにより求めることができます。

(3-36)

記述 は、行列 A の i番目の行です。 (存在するなら) （3-36） 式 での実行可能な 点は、等式制約を満足する値に x を設定することにより求めることができます。

これは、等式制約の組から作成されるoverdeterminedまたは、underdeterminedの どちらかの連立方程式を解くことにより得られます。 この問題の解が存在する場 合、余裕変数 は、この点で最大不等式の制約になります。

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

(3-37) LP 法を使って実行可能な点が求められたならば、メインの QP 部分に移行しま

す。 探索方向 は、つぎの線形方程式を解くことにより求まる を使って初期

化されます。

(3-38)

ここで、 は現在の繰り返し での目的関数の勾配です(つまり )。

dˆ_k A_idˆ_k>0 i, = 1, ,… m

λ_k A_k^Tλ_k = c

λ_k x_k

λ_k

γ ℜ∈minimize,x∈ℜⁿ γ

A_ix = b_i i = 1, ,… m_e A_ix–γ≤b_i i = m_e+1, ,… m A_i

γ

g_k

dˆ_k = –Z_kZ_k^Tg_k

dˆ

k dˆ

1

Hdˆ₁ = –g_k

g_k x_k Hx_k+c

制約付き最適化

実行可能な解がQP問題に対して求められないならば、メインのSQPルーチンに 対する探索方向 は、 を最小にするものとして決められます。

ライン探索とメリット関数

QP サブ問題の解は、ベクトル を作り、つぎの新しい繰り返しを作成するため

に使います。

(3-39) ステップ長パラメータ は、 メリット関数の中で十分な減少を生じるように決 定されます。 メリット関数は、Han [24]、Powell [35]により導入され、つぎの型を しています。

(3-40)

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

(3-41)

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

(3-42)

ここで、 は、Euclidノルムを示します。

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

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

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

dˆ_k γ

d_k

x_k+1 = x_k+αd_k α_k

Ψ( )x f x( ) r_i⋅g_i( )x

i=1 me

∑

^{ri⋅max{0 g, i( )x }}

i=me+1 m

∑

+ +

=

r_i (r_k+1)_i λ_i 1

2---(( )r_{k i}+λ_i)

 , 

 

 

i , max

= = i = 1, ,… m

r_i r_i f x∇( )

g∇ _i( )x

---=

このアルゴリズムは、制約により定義された多面体の端に沿って、ある頂点から 他の頂点まで移動し、 各ステップで、目的関数値 f^T x が減少します。 この節では、

頂点での最適解を返すオリジナルのシンプレックスアルゴリズムの改良された バージョンについて述べます。

このセクションは以下のトピックスをカバーしています。

• 3-35ページの"メインアルゴリズム"

• 3-36ページの"前処理"

• 3-36ページの"シンプレックスアルゴリズムの使用"

• 3-37ページの"基底変数と非基底変数"

• 3-38ページの"参考文献"

f^Tx

minx subject to A x⋅ ≤b

Aeq x⋅ = beq lb≤ ≤x ub

制約付き最適化

メインアルゴリズム

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

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

• 段階 2 − オリジナルの問題に対する最適解を計算します。

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

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

です。 ここで P_j(x_j) は以下のように定義されます。

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

は、

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

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

段階 2. 段階 2で、 アルゴリズムは、オリジナルの問題を解くために、段階 1から の初期点から始め、シンプレックスアルゴリズムを適用します。 各繰り返しで、ア

P P_j( )x_j

j

∑

=

P_j( )x_j

x_j–u_j 0 l_j–x_j

if x_j>u_j if l_j≤ ≤x_j u_j if l_j>x_j







=

P_j

j x

∑

min subject to A x⋅ ≤b

Aeq x⋅ = beq

ドキュメント内 optimize_print.book (ページ 138-200)