統計数理（2005）
第 53 巻 第 2 号 361–373
c 2005 統計数理研究所


特集「計算推論―モデリング・数理・アルゴリズム―」［研究ノート］

双対逐次 2 次計画および切除平面法による状態制約最適制御問題の解法
伊藤 聡†（受付 2005 年 5 月 23 日；改訂 2005 年 8 月 8 日）

要旨

双対逐次 2 次計画と切除平面法に基づく状態不等式制約つき非線形最適制御問題の数値解法について述べる．本解法は切除平面法に基づいて緩和された有限個の汎関数状態不等式制約を持つ最適制御問題に対して双対逐次 2 次計画法を適用するものであり，これにより元の状態不等式制約つき最適制御問題は一連の有限次元 2 次計画問題に変換される．本手法は状態不等式
制約条件に対する乗数の非減少階段関数による逐次近似に基づく解法とみることもでき，特に高次の状態制約つき問題に対して有効であると考えられる．
キーワード： 最適制御問題，状態制約，逐次 2 次計画，切除平面法．
1.

はじめに—制約の緩和

状態軌道に関する不等式制約を持つ非線形最適制御問題


 T



min g0 x(T ), T +
u,T

(1.1)

subj. to



0





f0 x(t), u(t), t dtẋ(t) = f x(t), u(t), t ,g x(t), t ≤ 0,

0 ≤ t ≤ T,

x(0) = x00≤t≤Tを考える．ただし，u ∈ L2 ([0, T ], Rr )，x(t) ∈ Rn ，g0 : Rn × R → R，f0 : Rn × Rr × R → R，
f : Rn × Rr × R → R，g : Rn × R → Rm である．以下では，簡単のため，終端時刻 T は固定かつ有限また m = 1 とするが，ほぼすべての結果は簡単な変換によりそのまま一般の場合に拡張することができる（Ito and Shimizu, 1990b; Liu et al., 2001 などを参照）．ここで，関数 g0 ，f0 ，
f および g が十分滑らかであること，そして適当な制約想定などを仮定する（詳細は Ito andShimizu, 1990b; Shimizu and Ito, 1994 を参照）こととし，これらの仮定のもとで，任意の制御
入力 u ∈ L2 ([0, T ], Rr ) に対してシステム方程式を満たす状態軌道 x( · ; u) が一意に存在することを保証する．
いま，システム方程式の存在のもとで，目的汎関数および制約汎関数をそれぞれ
(1.2)
(1.3)

†





γ0 (u) := g0 x(T ; u) +

 T
0





f0 x(t; u), u(t), t dtγ(u, t) := g x(t; u), t

統計数理研究所：〒106–8569 東京都港区南麻布 4–6–7



362

統計数理 第 53 巻

第 2 号 2005と定義すると，最適制御問題（1.1）は簡単に
min γ0 (u)
u

(1.4)

γ(u, t) ≤ 0

subj. to∀t ∈ [0, T ]と書ける．ここで，[0, T ] の有限部分集合 {t1 , t2 , . . . , tl } を考え，上の状態制約式を l 個の汎関数制約式
γi (u) := γ(u, ti ) ≤ 0,

(1.5)

i = 1, 2, . . . , lに置き換えると，緩和問題
min γ0 (u)
u

(1.6)

subj. toγi (u) ≤ 0,

i = 1, 2, . . . , l

を得る．
2.

双対逐次 2 次計画

この節以降 4 節までは，緩和された最適制御問題（1.6）を解くためのアルゴリズムについて考える．逐次 2 次計画法の意味での方向発見問題は，k 番めの反復において(2.1)

1
s, B k s2γi (uk ) + ∇γi (uk ), s ≤ 0,

min ∇γ0 (uk ), s +ssubj. toi = 1, 2, . . . , lのように与えられる（ASNOP 研究会, 1991; 矢部・八巻, 1999; Gruver and Sachs, 1980 などを参照）．ここで，uk は現在の近似解，B k は Lagrange 汎関数
L(u, λ) := γ0 (u) +lλi γi (u)i=1に対して L2 ([0, T ], Rr ) 上の双線形作用素として表現される 2 階偏微分 ∇2uu L(uk , λk ) の準 Newton
法の意味での近似であり，λk ∈ Rl は補助問題（2.1）のアフィン不等式制約に対する Lagrange 乗数の現在の近似である． · , ·  は L2 ([0, T ], Rr ) における内積を意味するものとする．補助問題（2.1）の解として得られる探索方向 sk に対して，次の近似解は適当なステップ幅 αk を用いて
uk+1 := uk + αk sk(2.2)と与えられる．
方向発見補助問題（2.1）は有限個の制約条件を持つ無限次元の決定問題であるが，いわゆる線形 2 次（LQ）形式で与えられている．したがって，その双対は以下のように有限次元の 2 次計画問題となる（Shimizu et al., 1986; Ito and Shimizu, 1990b; Shimizu and Ito, 1994）．
(2.3)

T
1
1
max − λT Λk λ − β k λ − δ k
2
2

λ∈Rl

subj. toλ≥0

ただし，係数 Λk ，β k および δ k を(2.4a)

Λkij := ∇γi (uk ), (B k )−1 ∇γj (uk ),(2.4b)

βik := ∇γi (uk ), (B k )−1 ∇γ0 (uk ) − γi (uk ),i, j = 1, 2, . . . , li = 1, 2, . . . , l双対逐次 2 次計画および切除平面法による状態制約最適制御問題の解法363

δ k := ∇γ0 (uk ), (B k )−1 ∇γ0 (uk )

(2.4c)と定義した．双対補助問題（2.3）の解 λk に対して
sk = −(B k )−1 ∇γ0 (uk ) −

(2.5)

l


λki (B k )−1 ∇γi (uk )i=1が主補助問題（2.1）の解を与えることは容易に確かめられる．
3.

勾配関数

目的汎関数 γ0 は，システム方程式


(3.1)



ẋ(t; u) = f x(t; u), u(t), t ,x(0; u) = x0のもとで，式（1.2）により定義された．よく知られているように，Hamilton 関数 H0 : Rn × Rr ×
Rn × R → R を
H0 (x, u, ψ0 , t) := f0 (x, u, t) + ψ0T f (x, u, t)で定義すると，目的汎関数 γ0 の勾配は以下のように与えられる．
命題 1. 汎関数 γ0 の勾配 ∇γ0 (u) ∈ L2 ([0, T ], Rr ) は
(3.2)



∇γ0 (u)(t) = ∇u H0 x(t; u), u(t), ψ0 (t), t









= ∇u f0 x(t; u), u(t), t + ∇u f x(t; u), u(t), t ψ0 (t),0≤t≤Tで与えられる．ここで ψ0 : [0, T ] → Rn は随伴方程式

(3.3)





 ψ̇0 (t) = −∇x H0 x(t; u), u(t), ψ0 (t), t










= −∇x f0 x(t; u), u(t), t − ∇x f x(t; u), u(t), t ψ0 (t)ψ0 (T ) = ∇g0 x(T ; u)の一意解である．
他方，制約汎関数 γi ，i = 1, 2, . . . , l は，同じくシステム方程式（3.1）のもとで，
(3.4)



γi (u) = γ(u, ti ) = g x(ti ; u), tiと定義されているが，この場合はH(x, u, ψi , t) := ψiT f (x, u, t)と定義される Halmilton 関数 H : Rn × Rr × Rn × R → R を用いて以下が成立する（Shimizu andIto, 1994）．
命題 2. 汎関数 γi の勾配 ∇γi (u) ∈ L2 ([0, T ], Rr ) は
(3.5)



∇γi (u)(t) = ∇u H x(t; u), u(t), ψi (t), t= ∇u f x(t; u), u(t), t ψi (t),

0≤t≤T

364

統計数理 第 53 巻

第 2 号 2005

で与えられる．ここで ψi : [0, T ] → Rn は随伴方程式

(3.6)





ψ̇i (t) = −∇x H x(t; u), u(t), ψi (t), t







= −∇x f x(t; u), u(t), t ψi (t), 0 ≤ t ≤ tiψi (ti ) = ∇x g x(ti ; u), ti






ti < t ≤ T

ψi (t) = 0,の一意解である．
双対逐次 2 次計画の各反復における目的汎関数および制約汎関数の勾配の計算についてまと
めると以下のようになる（図 1 参照）．まず，現在の制御入力 uk ∈ L2 ([0, T ], Rr ) に対して，シ
ステム方程式（3.1）を解くことにより，これに対応する状態軌道 xk := x( · , uk ) が求まる（この段
階で汎関数値 γi (uk ), i = 0, 1, . . . , l が求まる）．次に，uk と xk に対して，l + 1 個の随伴方程式系（3.3）および（3.6）をそれぞれ解くことにより，随伴軌道 ψik ，i = 0, 1, . . . , l が求まる．uk と得られた軌道 xk および ψik ，i = 0, 1, . . . , l を式（3.2）と式（3.5）に代入することにより，勾配関数∇γi (uk )，i = 0, 1, . . . , l がそれぞれ求まる．
随伴系の右辺は Hamilton 関数の x に関する偏微分で構成されているので，関数 f および gが x に関してアフィン，すなわち
f（x，u，t）= Ax + f（u，t）
g（x，t）= cT x + g（t）

図 1.

勾配関数の計算手順．

図 2.

随伴軌道の平行移動．

双対逐次 2 次計画および切除平面法による状態制約最適制御問題の解法

365のときは，図 2 に示すように，制約汎関数 γi = γ( · , ti ) に対する随伴軌道 ψi は単に γ( · , T ) に対する随伴軌道を時間 T − ti だけ平行移動したものとなり，勾配計算は単純化される．
4.

準 Newton 更新

本節では近似作用素 B k の更新について述べる．双対逐次 2 次計画の枠組においては，式（2.4）にあるように，目的汎関数および制約汎関数の勾配は (B k )−1 ∇γi (uk )，i = 0, 1, . . . , l の形で現れる．初期近似を B 0 = I（恒等作用素）とするとき，準 Newton 法の典型的な更新規則である BFGS 公式は，逆作用素形式では以下のように表現される．
(4.1)(B k+1 )−1 = (B k )−1 +−1+(B k )−1 q k , q k pk , q k ]pk , pk [
pk , q k ]pk , (B k )−1 q k [](B k )−1 q k , pk [−kk
p , q pk , q k 

ここで
(4.2)

pk := uk+1 − uk

(4.3)

q k := ∇u L(uk+1 , λk ) − ∇u L(uk , λk )= ∇γ0 (uk+1 ) − ∇γ0 (uk ) +l

λki {∇γi (uk+1 ) − ∇γi (uk )}i=1r

とおいた．また，] · , · [ は L2 ([0, T ], R ) 上で
for each d ∈ L2 ([0, T ], Rr )]u, v[ d := v, d uと定義されるダイアディック作用素である．記述を簡潔にするため詳細は省くが，B k の正定値性を保証するためにいわゆる Powell の修正を施すものとする（Gruver and Sachs, 1980 などを参照）．
式（4.1）より，任意の d ∈ L2 ([0, T ], Rr ) に対して
(4.4)

(B k )−1 d = (B 0 )−1 d +

k−1


1+

j=0

−

(B j )−1 q j , q j pj , q j pj , d jppj , q j pj , d(B j )−1 q j , d j
p − j j (B j )−1 q jjjp , q p , q を得る．ここで，無限次元の作用素 (B k )−1 そのものを BFGS 公式（4.1）にしたがって更新することはできないが，Lasdon（1970）にあるように，関数 pj , (B j )−1 q j , j = 0, 1, . . . , k − 1 およびスカラ pj , q j , (B j )−1 q j , q j , j = 0, 1, . . . , k − 1 が過去の反復より既知であるから，式（2.4）で必要とされる 関数 (B k )−1 ∇γi (uk )，i = 0, 1, . . . , l は計算可能であることに注意されたい．以上の更新は，ここで取りあげた BFGS 公式だけでなく，Broyden 族のあらゆる更新規則に対して成立する．
1 次元探索のアルゴリズムとしては，精密ペナルティ汎関数
Q(u) := γ0 (u) + rl

max { 0, γi (u) }(r > 0)i=1に対して，Armijo の基準(4.5)

Q(uk + αsk ) ≤ Q(uk ) − ωαsk , B k sk (ω ∈ (0, 0.5), τ ∈ (0, 1))

366

統計数理 第 53 巻

第 2 号 2005

を満たすステップ幅 α ∈ (0, 1] を αk として用いることが考えられる．式（2.5）を用いると，式（4.5）における sk , B k sk  は
sk , B k sk  = ∇γ0 (uk ), (B k )−1 ∇γ0 (uk )
+2l

λki ∇γi (uk ), (B k )−1 ∇γ0 (uk )

i=1

+

l
l 


λki λkj ∇γi (uk ), (B k )−1 ∇γj (uk )

i=1 j=1と書けるが，この右辺にはまた (B k )−1 ∇γi (uk ) という形で勾配が現れる．
以上をまとめると，緩和された最適制御問題（1.6）に対する双対逐次 2 次計画のアルゴリズ
ムは以下のようになる．
アルゴリズム 1.
ステップ 1. を十分小さな正数とする．初期推定 u0 を選び，B 0 = I（恒等作用素），k = 0とする．
ステップ 2. 図 1 の手順で γi (u0 ), ∇γi (u0 ), i = 0, 1, . . . , l を求める．
ステップ 3. 式（4.4）により (B k )−1 ∇γi (uk )，i = 0, 1, . . . , l，さらに式（2.4a）,（2.4b）により係数 Λk および β k を計算する．
ステップ 4. 有限次元 2 次計画問題（2.3）を解き，その解 λki ，i = 1, 2, . . . , l に対して式（2.5）に
より探索方向 sk を求める．
ステップ 5. sk < ならば uk を問題（1.6）の最適制御として終了．
ステップ 6. uk+1 = uk + sk ，あるいは Armijo の基準（4.5）を満たすステップ幅 αk を用いて
k+1u
= uk + αk sk とする．
ステップ 7. 図 1 の手順で γi (uk+1 ), ∇γi (uk+1 ), i = 0, 1, . . . , l を，式（4.2）,（4.3）により pk ，q kを，さらに式（4.4）により (B k )−1 q k を計算する．
ステップ 8. pk , (B k )−1 q k , pk , q k , (B k )−1 q k , q k  を記憶し，k := k + 1 としてステップ 3 に
戻る．

5.

切除平面法

状態制約式を緩和する際に切除平面法の考え方を使う．最適制御問題（1.4）に対する切除平面法の基本的なアルゴリズムは以下のように書ける．
アルゴリズム 2.
ステップ 1. を十分小さな正数とする．l1 を適当な正整数とし，l1 個の要素からなる [0, T ]の部分集合
I1 = { t11 , . . . , t1l1 }を選び，k = 1 とする．
ステップ 2. 緩和問題
min γ0 (u)usubj. toγ(u，t) ≤ 0∀t ∈ Ik = { tk1 , . . . , tklk }双対逐次 2 次計画および切除平面法による状態制約最適制御問題の解法

367

をアルゴリズム 1 により解き，解 uk（および乗数 λki ，i = 1, . . . , lk ）を求める（初期推定として，
k > 1 に対しては前反復で得られた uk−1 を用いることができる）．
ステップ 3. 元の制約条件における最大侵害量
δ(uk ) := max γ(uk , t) = max g(x(t; uk ), t)t∈[0,T ]t∈[0,T ]を計算し，その大域的最適解を
tklk +1 := arg max γ(uk , t)t∈[0,T ]とする．
ステップ 4. δ(uk ) < ならば uk を問題（1.4）の最適制御として終了．
ステップ 5.
lk+1 := lk + 1, . . . , tk+1
Ik+1 := Ik ∪ {tklk +1 } = { tk+1
1
lk+1 }とし，k := k + 1 としてステップ 2 に戻る．
ステップ 1 において，可能なら l1 = 0 すなわち I1 = ∅ として構わないが，この場合は k = 1
における緩和問題は無制約最適制御問題となる．この基本アルゴリズムに対して，（i）反復が進むにつれて Ik の要素数が増大するのを防ぐため，適当な規則に基づき Ik から不要な点を削除する，（ii）ステップ 2 で緩和問題を解く際にその終了判定基準をステップ 4 のそれと連動させる
ことにより近似的に解く，などの改良を行なうことができる（Lai and Wu, 1992; Wu and Fang,
1999）．このうち（i）だけについて説明すると，ステップ 5 をステップ 5 .
Ik := { tk1 , . . . , tklk | λki > 0 }lk+1 := |Ik | + 1, . . . , tk+1
Ik+1 := Ik ∪ {tklk +1 } = { tk+1
1
lk+1 }とし，k := k + 1 としてステップ 2 に戻る．
（∀k）が成立するなら，アルゴリズムの収束性を損なわない．
に置き換えても，tklk +1 ∈ Ik+1

6.

最適性条件

適当な仮定のもとでは，最適制御問題（1.1）に対する Karush–Kuhn–Tucker（KKT）型の最適性必要条件は，システム方程式および状態不等式制約条件に加え，次のように与えられる（Girsanov, 1972; Ito and Shimizu, 1990a）．
(6.1a)(6.1b)ψ(t) = ∇g0 (x(T )) +

 T
t





{∇x f0 x(τ ), u(τ ), τ

+ ∇x f x(τ ), u(τ ), τ ψ(τ )} dτ +
 T

(6.1c)

 T

g x(t), t dλ(t) = 00(6.1d)

∇u f0 x(t), u(t), t + ∇u f x(t), u(t), t ψ(t) = 0

λ: 非減少関数

t






∇x g x(τ ), τ dλ(τ )368

統計数理 第 53 巻

図 3.

図 4.

第 2 号 2005

状態制約の接合点における乗数 λ の跳躍．

乗数 λ の非減少階段関数による近似．

ここで，ψ および λ はそれぞれシステム方程式，状態不等式制約条件に対する Lagrange 乗数である．図 3 に示すように λ が状態制約の接合点（制約境界に触れるあるいは境界から離れる点）において跳躍し得ること，したがって随伴方程式（6.1b）により ψ もこれらの点において不連続になり得ることに注意されたい（Berkovitz, 1961; Jacobson et al., 1971; 坂和, 1980; 嘉納,1987 などを参照）．なお，図 3 では乗数 λ が右連続な非減少関数として表現されているが，両端点を除けば不連続点 ti における値 λ(ti ) そのものは本質的ではなく，最適性条件（6.1）においては跳躍量 λ(ti + 0) − λ(ti − 0)（ただし両端ではそれぞれ λ(0 + 0) − λ(0)，λ(T ) − λ(T − 0)）のみが意味を持つ．
乗数 λ は有界変動関数であるが，同時に Stieltjes 測度とみることもできる．このとき，状態制約式を汎関数制約（1.5）で緩和することは，測度 λ を λ({ti }) = λi ，i = 1, 2, . . . , l なる離散測度で近似することに等しい．あるいは，関数 λ を
(6.2)

λ

l


λ i χii=1で近似することに等しいと言える．ただし χi は
χi (t) =0,0 ≤ t < ti

1,ti ≤ t ≤ T

で与えられる [ti , T ] の定義関数である．l = 2 の場合を図 4 に示す．
さて，近似（6.2）のもとでは，最適性必要条件（6.1）における乗数（随伴変数）ψ と本手法の枠

双対逐次 2 次計画および切除平面法による状態制約最適制御問題の解法

369

組における (ψ0 , ψi , λi ) は(6.3)

ψ = ψ0 +

l


λi ψi

i=1

により関連づけられる．実際，任意の u および λi , i = 1, 2, . . . , l に対して，システム方程式（3.1）
および随伴方程式系（3.3）,（3.6）により一意に定まる ψi , i = 0, 1, . . . , l を用いて，式（6.3）により関連づけられた ψ は随伴方程式（6.1b）を自動的に満たすことが容易に確かめられる．さらに，勾配（3.2）,（3.5）を考慮すると，条件（6.2）,（6.3）のもとでは，最適性必要条件（6.1）のうち残りの条件がそれぞれ
∇γ（
0 u）+
l


l


λi ∇γ（i u）= 0i=1

λi γ（i u）= 0i=1

λi ≥ 0,i = 1, 2, . . . , lと等価であることがわかるが，これは緩和された最適制御問題（1.6）の KKT 型最適性必要条件に他ならない．
7.

数値実験

例題として，次のアフィン状態制約つき LQ 問題（Jacobson and Lele, 1969）を考える．
 1

min
u

(7.1)

0

(x21 + x22 + 0.005 u2 ) dt

subj. toẋ1 = x2 ,x1 (0) = 0ẋ2 = −x2 + u,x2 (0) = −12

x1 − 8 (t − 0.5) + 0.5 ≤ 0

この状態不等式制約条件は，左辺を 2 回時間微分するとはじめて制御入力 u が陽に現れるため，2 次の状態制約と呼ばれる．この問題の最適状態軌道は時刻 t = 0.47 の近くのただ一点で制約境界に触れ，この点において乗数 λ が跳躍することが知られている．
この例題（7.1）に対してアルゴリズム 1 および 2 による数値実験を行なった．切除平面法による外側の反復（アルゴリズム 2）については，I1 = ∅ とし，u0 (t) ≡ 1 を初期推定ノミナル制御，= 10−10 すなわち
δ(uk ) = max g(x(t; uk ), t) < 10−10t∈[0,T ]を終了判定基準とし，ステップ 5 の代わりにステップ 5 を用いた．また，双対逐次 2 次計画による内側の反復（アルゴリズム 1）においては，ステップ幅は αk = 1 に固定，積分は常微分方程式系に対しては Runge–Kutta–Gill 法，内積に対しては Simpson 則を用い，終了判定基準は= 10−10 すなわち sk < 10−10 が満たされるまで反復を続けた．
実験の結果を表 1 に示す．表中には，アルゴリズム 2 の k 番めの反復に対して，Ik の要素数
lk ，Ik の各要素 tki とアルゴリズム 1 で得られた対応する λki の組，元制約の最大侵害量 δ(uk )および新たに Ik+1 に加える tklk +1 を列記した．図 5 は，最終反復（k = 9）における状態軌道 x91（破線は制約境界）と式（6.2）の意味での近似として得られた乗数 λ9 を示す．最終反復において
は t = 0.47 の近傍の 2 点で λ が跳躍しているが，切除平面法に基づく本手法は非許容点列を生370

統計数理 第 53 巻
表 1.

図 5.

第 2 号 2005

例題（7.1）に対する実験結果．

最終反復における状態軌道 x91 と乗数 λ9 ．

成する近似解法であり，δ(u9 ) < 10−10 の範囲で制約条件を僅かに侵していることに注意されたい．例題（7.1）の状態不等式制約条件において x1 を x2 に置き換えた場合のように，低次の状
態制約に対しては，最適状態軌道は一定の時間区間にわたり制約境界に張りつく形となることが多い．このような場合，乗数 λ はこの活性な時間区間において絶対連続となる．これに対して，例題（7.1）を含めた高次の状態制約つき問題においては，乗数 λ は離散的になりやすい
ため，特にアルゴリズム 2 のステップ 5 などと組み合わせたとき，実際に解くべき補助問題（2.3）の低次元化が期待できる．

双対逐次 2 次計画および切除平面法による状態制約最適制御問題の解法

8.

371

おわりに

状態制約がある場合，6 節で述べたように，制約の接合点において乗数の跳躍および随伴変数の不連続性が起こり得るため，システム方程式および随伴方程式からなる多点境界値問題を含む最適性条件を解くことは非常に難しい（このようなアプローチは間接法と呼ばれる）．本稿
で述べた手法は，切除平面法を用いて乗数 λ の逐次推定を行なうことにより，随伴変数 ψ を式（6.3）のように分解し，結果として随伴方程式（6.1b）を常微分方程式系の初期値問題すなわち（3.3）および（3.6）に変換する．
また，関数空間上で数理計画法的に解くいわゆる直接法を適用する場合，いずれかの段階で離散化（すなわち無限次元変数の有限次元近似）を行なう必要がある（最終段階でこれを行なうのが Shimizu and Ito, 1994 の方法である）．通常，当初から制御入力と制約条件（あるいは制御入力のみ）を離散化することにより非線形計画問題（あるいは半無限計画問題）に変換する手法（後者については Ito et al., 2000; Liu et al., 2001 などを参照）が使われることが多いが，結果として得られる有限次元の補助問題は潜在的に大規模なものになりがちである．これに対して，本稿で述べた手法は，特に高次の状態制約つき問題に対して，補助問題の低次元化が期待でき，有効なアプローチになり得ると考える．

参 考 文 献ASNOP 研究会
（1991）. 『パソコン FORTRAN 版非線形最適化プログラミング』，日刊工業新聞社，東京．
Berkovitz, L. D.（1961）. Variational methods in problems of control and programming, Journal of
Mathematical Analysis and Applications, 3, 145–169.
Girsanov, I. V.（1972）. Lectures on Mathematical Theory of Extremum Problems, Lecture Notes in
Economics and Mathematical Systems, 67, Springer-Verlag, Berlin.
Gruver, W. A. and Sachs, E.（1980）. Algorithmic Methods in Optimal Control, Pitman Advanced
Publishing Program, Boston, Massachusetts.
Ito, S. and Shimizu, K.（1990a）. Necessary conditions for constrained optimal control problems via
mathematical programming, Numerical Functional Analysis and Optimization, 11, 267–281.
Ito, S. and Shimizu, K.（1990b）. A solution of optimal control problems with terminal state equality
constraints by quasi-Newton method, Proceedings of the 29th SICE Annual Conference, Vol.
II, 717–720.
Ito, S., Liu, Y. and Teo, K. L.（2000）. A dual parametrization method for convex semi-infinite programming, Annals of Operations Research, 98, 189–213.
Jacobson, D. H. and Lele, M. M.（1969）. A transformation technique for optimal control problems
with a state variable inequality constraint, IEEE Transactions on Automatic Control, AC-14,
457–464.
Jacobson, D. H., Lele, M. M. and Speyer, J. L.（1971）. New necessary conditions of optimality forcontrol problems with state-variable inequality constraints, Journal of Mathematical Analysis
and Applications, 35, 255–284.

嘉納秀明（1987）. 『システムの最適理論と最適化』，コロナ社，東京．
Lai, H. C. and Wu, S. Y.（1992）. Extremal points and optimal solutions for general capacity problems,Mathematical Programming, 54, 87–113.
Lasdon, L. S.（1970）. Conjugate direction methods for optimal control, IEEE Transactions on Automatic Control, AC-15, 267–268.
Liu, Y., Ito, S., Lee, J. H. W. and Teo, K. L.（2001）. A semi-infinite programming approach tocontinuously constrained LQ optimal control problems, Journal of Optimization Theory and

372

統計数理 第 53 巻

第 2 号 2005

Applications, 108, 617–632.
坂和愛幸（1980）. 『最適化と最適制御』，森北出版，東京．
Shimizu, K. and Ito, S.（1994）. Constrained optimization in Hilbert space and a generalized dual
quasi-Newton algorithm for state-constrained optimal control problems, IEEE Transactions
on Automatic Control, AC-39, 982–986.
Shimizu, K., Fujimaki, S. and Ishizuka, Y.（1986）. Constrained optimization methods in Hilbert
spaces and their applications to optimal control problems with functional inequality constraints, IEEE Transactions on Automatic Control, AC-31, 559–564.
Wu, S.-Y. and Fang, S.-C.（1999）. Solving convex programs with infinitely many linear constraints
by a relaxed cutting plane method, Computers and Mathematics with Applications, 38, 23–33.
矢部 博，八巻直一（1999）. 『非線形計画法』，朝倉書店，東京．

Proceedings of the Institute of Statistical Mathematics Vol. 53, No. 2, 361–373 (2005)

373

A Numerical Solution of State-constrained Optimal Control Problems by
Dual Sequential Quadratic Programming with Cutting Plane Strategy
Satoshi Ito
The Institute of Statistical Mathematics

An approach based on dual sequential quadratic programming (SQP) coupled witha cutting plane strategy is proposed as a numerical method for solving optimal control
problems with state trajectory inequality constraints. This method applies the dual SQP
to a sequence of relaxed problems with a finite number of pointwise state inequality
constraints generated in a cutting plane sense, and thus converts the original optimal
control problem into a series of finite-dimensional quadratic programs. The approach
can also be viewed as a sequential estimation procedure of the multiplier for the state
trajectory constraints by nondecreasing step functions and is hence expected to perform
well on high-order state constraints.

Key words: Optimal control problem, state constraint, sequential quadratic programming, cutting
plane method.