最適制御問題の直接解法と滑空機飛行への応用に関 する研究

九州大学学術情報リポジトリ

Kyushu University Institutional Repository

最適制御問題の直接解法と滑空機飛行への応用に関 する研究

河邉, 博康

https://doi.org/10.11501/3154824

出版情報：Kyushu University, 1999, 博士（工学）, 論文博士 バージョン：

権利関係：

第4章 最適制御問題の直接解法

4.2 Modified Direct Optimization Method

4.2.1 計算結果と問題点

4.2.

ダイナミック・ソアリングの最適制御問題 を， DCNLP法とsystem constraintsを用いて計 算した. 計算 に使用した滑空機は， thermal-t仔thermal ftightの時と同じNIPPI PILATUS B4

のデ」タを使い， 初期高度をh(O)=5mとした. 部分区間数が20の時 に， DCNLP法の解を 図4.9に示し， system constraintsを用いた時の制御量を図4.10に破線で示した. DCNLP法の グラフで， 制御量である揚力係数CLとロ}ル角μの解は滑らかであるのに対し， system constraintsを使った解では尖った部分が現れている . これは， 積分公式で離散的に微分ノf

程式を解いているので， 節点での状態量zの傾きが不連続になっており， その影響が制 御量uに現れていると考えられる. よって， その影響を受けないように， 制御量を線形

補間してスムージング を行う.

第4章

最適制御問題の直接解法

2

3

4

0.2 0.4 0.6 0.8 ₁ 0.2 0.4 0.6 0.8 15

1 '0 _<L20TO五 0.8 ₁

1 4 云え!?

達

d、

『 匂3云^，

F 量、

11\ ^{/ J1}

4} '1 _，-I

^{' i}

i

J ^、 R晴

^{、Eq )}

⁵

=巴

=・501

、

od l

^t 1

1 1 -10

0.20.4t々f0.60.8 1 10 0.20.4

bTf 0.60.8 1

図4.9

DCNLP法を用いた時の状態量と制御量.

4.2.

4.2.

- ，‘、・'‘、 、、

..ι 、， 、一 ・・ ‘，，，

4・・

・

，/ ....

， 、

，，/ \_，

・ 6・

'.: �.

ヘi \. _'./\

'1 �!'， ，.

^.}�

，\

，\; ・ '.:'

、l ベJ

，、J 奇旬、ML句、聞きにと

ーいいv‘ ve古芭 』 - \ ・ 冗U

Idνトー 1且

最適制御問題の直接解法

叫にJ 第4章

0.6 0.8 .4 0.6

0.4

。も

0.2

Simpson's system constraints

• ‘ • • ， • ‘ • 、

4. 1 •

• . ‘ •

，， v 、 • 一"" • 、 • ‘ • .‘ • ‘ ‘ •.

‘‘ •

. I l

a-. • •

白，

• ， ， ， ， ， d • ， -e' ，， ， • • ， ， l l '

，

，

，

• ， • ， ， ， • 企，• ， ， ， • ， .

， ，l l

F、ぜ

量、J Rω芝とUWNh阿川口。、

1 ト /\/ 〔 八

� U

0.8 0.4

0.6

0.8 0.2 0.6 0.4

。も

0.2

Fourth-degree Gauss-Lobatto system constraints

• • • 、 ' • ι

...

• • ‘

‘ ••

、 ‘、

‘ •

、 ・.

• • • ‘ -- ‘.

• 、 • . ， F

‘He -

- Jf

‘l e 'l - ー・ 白守 ，v， •

，

• ，

，

， • ， a ， ， a， . ， ， ， ， M le'

， ，

AP !

，

， 、 • ， ん ， ， a， ，

ζJ

，、d R也」とhws~~。、

� U

-10 0.8

\40.60.8 i/ tf- '- -.� 1 - 10 02 0.40j6 " "._ _{" 'i /tf}

Fifth-degree Gauss-Lobatto system constraints

図4.10 System constraintsを用いた時の制御量.

。も 0.2

第4章 最適制御問題の直接解法 4.2.

4.2.2 制御量の線形補間法

節点での制御量の値を使って， 図4.11のように節点と節点の聞のcollocation pointsで の制御量を線形補間することにする. ³⁴⁾

Simpson's system constrainもの場合は， 簡単に次のようになる.

U，. 二(Ui+l + Ui)

"

2 ( 4.68)

Fourth-degree Gauss-Lobatto system constraintsの場合は次のような比例関係からcollocation pointsでの制御量を補間する.

(

^t，一

字方

^-t，

)

(t'+l -t，) = (U1 - Ui) (町村- ^Ui)

( 一一一

^{ムt ムt}

^夜

¹

⁾

^ω= (U1 - ui)川-Ui)

2 2

(4.69) (4.70)

Ul =川-u;)

G-法)

^{+ u， (4.71)}

同様にして

U2二(同+l-Ut)

(j

⁺

;去)

^{+ Ui (4.72)}

Fifth-degree Gauss-Lobatto system constraintsの場合は次のようになる.

十

^一

手伝

^-t，

)

^{: (t'+lーし}

( ^ヂ叫 )

^ムt=

(4.73)

(4.74)

u +

\11EE--JJ ZW 14Z4 1i一qL /flit-\

u + u

一一u _(4.75)

同様にして

U3= (U'+l - u，)

(�

⁺

�jf)

^{+ Ui}

以後， 制御量を線形補間する方法をmodi自ed direct optimization methodと呼ぶことにする.

(4.76)

第4章 最適制御問題の直接解法

Ui+j

ul/?γi !

Ui// L j i i

Ic 12(13ノ 1 i+j

図4.11制御量の線形補間法

4.2.

第4章 最適制御問題の直接解法

4.2.3 Modified Direct Optimization Method

4.2.

前節の方法で制御量を線形補間して求 め て， _NLP problemを次のように書き換える.

mllllmlze:ゆ(X_{f )}

subject to: cL三c(X)三cu

ここで， modified Simpson's system constraintsの場合， NLP problemは次のようになる.

C5(Xi ， xi+l， Ui， Ui+l) =0， .for i=O，l，2，...，ns-l

c( x) == [C 5，0， • • • ， C 5，j ，ψ0，ψj'!VO，... ，!VjJT

(4.77)

(4.78)

(4.79) ( 4.80)

ここで， modified fourth-degree Gauss-Lobatto sysもem constraintsの場合， NLP problemは次のよ うになる.

C4，l(Xi， Xi，^c， Xi+l， Ui， Ui+l) = 0， for i = 0，1，2，. ••，11s - 1 C4p2(zu忽i，c ， Xi+l ， Ui， Ui+l) = 0， for i二0，1，2，• • • ，ηs -1

(4.81)

(4.82) c(X) = ^[Cりん ，C4，1，.f，C4，2，O.... ，C4，2，1'ψ0'ψf ' tP ^0，・ ，tP.rr ， T (4.83)

ここで， modified fifth-degree Gauss-Lobatto system constraintsの場合， NLP problemは次のよう になる.

C5，I(Xj， Xj，c ， Xi+1， ^{Ui， Ui+l)} = 0_， f07' i二⁰，1，2_，• • • ，ηs - 1

C5，3(Xi， Xi，c， Xi+1， Ui， Ui+l) = 0， for i二0，1，2，• • • ，ηs - 1

( 4.84) (4.85) c(X) =

[

C5^ム0 ，.. .，C5，1，J， C5，2，0:・・. ，C5，2，J，ψ0，ψf'!VO，...，!VjJT ^(4.86) ベクトルmの拘束条件は 次のように表される.

XLくXくXu ^{( 4.87)}

ここで， CL二Cu， XL二Xuとすれば等式拘束条件となる.

CLニ[0，・ ，0，!VL， • • • ，!V Lf， Cu = [0， • • • ，0， tPu， ^{• • •} ，世U]T (4.88)

与えられた拘束条件のもとで， 目的関数を最小にするX = [x，u

f

を求 め る_.

第4章

最適制御問題の直接解法

4.2.

4.2.4

修正後の計算結果

Modified direct optimization methodを用いて， ダイナミック ・ ソアリングを解き直した 結果を図4.12の実線で示した. 実線を見ると制御量に尖った部分がなくなり， 滑らかに なってより現実的な解になっていることがわかる. ひとつ の例として， ITlodi五ed日fth-degrce Gauss-Lobatto system constraintsによって得られた結果を見ることにする. 図4.13を見ると，

アホウドリと同様に滑空機でもダイナミック・ソアリングが行われていることがわかる.

図4.14には， 状態量として， 対地速度'U， V，ωと位置x， ^y， hが，制御量として揚力係数CLと ロール角μが示しである. グラフの横軸は， 時間を終端時刻tJで割って無次元化して表し ている. 対地速度U，V，ωの初期値と終端値が， 一致しているので境界条件を満たしてい ることがわかる. 高度hについても， 初期高度と終端高度がほぼ一致して評価関数を最 小にしていることがわかる. 次に どのような制御を行っているか見てみると， ロール 角を大きくして旋回しているときに， それぞれ， 揚力係数を大きくしていることがわ かる. 特に， 高度の頂点付近では， 滑空機の速度が落ちているため， CLmaxの制限いっぱ いまで揚力係数を大きくして旋回しているのがわかる. また， 径路の終点が， 始点より も風下側になっており， このことから， 1 cycle飛行が終わるごとに少しずつ風下側へず れながら飛行を繰り返していくことがわかる. 以上より，妥当な最適解が得られている

ことがわかる.

4.2.

nhwミとMWNhqミミ

最適制御問題の直接解法

、川口)

第4章

0.6 0.8

0.6

0.4

0.4

。も

0.2

Si'mpson' s system constraints

1ゐ5

明、}

匂

CJ;)

� �

=3・501

ζ f

� u

0.8 0.6 0.4 0.2 0.8

0.4

。も 0.2

Fou吋h・degree Gauss-Lobatto system constraints

0.8

�5 ω ーで司、}

�

句足

ミー5

。

0 ，

- 10

0.4 _0.0 -0.81 -

^{J. v}

'û 0.2 0.4

_O.

t/i可 t/i可

Fifth-degree Gauss-Lobatto system constraints

. ‘ ‘ .，， ， e'

O.

0.2

modified direct optimization method system constraints

図4.12 Modified direct optimization methodを用いた時の制御量

第4章 最適制御問題の直接解法 4.2.

-』ノ例 /f t

-n

印 必 川w お 初 お 初 日 ω 5 0 印必 川切 お ぬ お 初 日 印 5 0 LI ノ 附/ rR ιμ

nU Fhu nu

A吐 nu

nu qL qd nu nu

m jリ/ nU 4g' 内G EF )

x ⑤ ιい /

A宝 nU ω

jけソ 1nU 4F 凶0 メμ1

uo ν〆

nU nu

A斗A AU

。

図4.13 Modi五ed五fth・degree Gauss-Lobatto system constraintsを用いた時の最適飛行径路.

第4章 最適制御問題の直接解法

さ 芯1

-20

-30、 夜2 0.4 0.6

3

ま2益 ²

i

^l量

1.4

0.8

r、戸、 r、 A r、，- r、c1・

4.2.

4 2

ー、E』、pH 、

-2

1 5

芝品、10

5

↓/l _{.2 0.4}

0.6

0.8 4 50

: 4EE 3 2

、。

0.2 0.4 0.6 0.8

1 3、h、、同

E q

ζ3.

二a^，

5

ミA

^ー5

-

¹

0�

ハ内

0.2 tI J 3f

O五-0.8

図4.14 Modified fifth-degree Gauss-Lobatto system constraintsを用いた時の状態量と制御量

第4章 最適制御問題の直接解法 4.2.

4.2.5 修正後の精度の比較

4.1.7節と同じ例題を， modified direcもoptimization methodで解き， 精度の比較を行った. 制 御量を線形補間しているため， 修正後は特に制御量のグラフに曲線が多い場-合， 精度 が落ちることが予想、されるが， 例題の制御量の解析解のグラフの形が線形なので1 表 4.3を見ると， どの手法も相対誤差が10-16で非常に精度がよいことがわかる. また計算 時間については， collocation pointsの制御量を線形補間した分だけ， 非線形計画法の変数 が減るので， 表4.4を見ると， 修正後は計算時間が短くなっていることがわかる.

それぞれの修正後の手法を使った結果と解析解を比較したものを図4.15から図4.17に 示した. これらを見ると どの手法も解析解と非常に良く一致しており， modified direcも optimization methodが有効であることが示された.

表4.3修正前と修正後の相対誤差の比較

System consもraints 修正前の相対誤差 修正後の相対誤差 部分区間数

Sim pson's 5.019 X 10-13 6.661 X10-16 20

Fourth 4.996 X 10-16 1.665 X10-16 20

Fifth 1.665 X 10-16 1.665 X 10-16 20

表4.4 修正前と修正後の計算時間の比較

System constraints 修正前の計算時間(s) 修正後の計算時間(s) 部分区間数

Simpson's 129.60 24.96 20

Fourth 414.45 203.4 20

Fifth 426.00 274.12 20

第4章 最適制御問題の直接解法 4.2.

〉売 =訟

ーも

t

解析解

也 ^{一一 一 一} Modified Simpson' s

図4.15 Modified Simpson 's system constraintの解と解析解の比較

第4章 最適制御問題の直接解法 ^4.2.

同 _ニ』

ーも

解析解

也 一一一一 Modified fourth-degree

図4.16 Modified fourth-degree Gauss-Lobatto system constraintsの解と解析解の比較

第4章 最適制御問題の直接解法 4.2.

同 =ー

-も

解析解

� Modified fifth-degree

図4.17 Modified fifth-degree Gauss-Lobatto system constraintsの解と解析解の比較.

第4章 最適制御問題の直接解法 4.2

4.2.6 Direct optimization methodのまとめ

本章では， direct optimization methodとして， DCNLP法と積分公式を使った方法の比較検 討を行った. 高次の積分公式を使った方法は少ない分割数と計算時間で低次の台形公ず 程度の相対誤差のオーダーの精度を得ることができた. 分割数や相対誤差のオーダー が同じであれば 積分公式を使った方法より もDCNLP法の方 が計算時間の点で優れて いることがわかった.

特に，積分公式を使った方法で複雑な最適制御問題を解くと制御量に不連続な部分が

現れた . そこで， 制御量のcollocation pointを線形補間してスムージングを行ったmodi五ed direct optirnization method を提案し， 複雑な最適制御問題を解 いたところ制御量に不連 続な部分が無くなり， 現実的な解を得ることができた. また， modified direct optimization method は， direct optimization method に比べて著し い 精度の低下もなく， 計算時間も非線 形計画法の変数が減った分だけ短くなり， 複雑な最適制御問題を解くのに非常に有効な 方法であることがわかった. さらに， 間接法では， 最適制御問題の定式化の段階で問題 を簡単化するために様々な工夫が必要であったが， 本章で示した直接法では， 複雑な最 適制御問題であっても問題の定式を簡単化することなく解くことができ， 問題 の取り 扱いが非常に容易であるという利点も示された.

第5章

ダイナミック・ ソアリングの解析結果

5.1

風速勾配の影響

風速勾配の大きさを

告

= 0.1，0.12，0.14と変化させてその影響を調べた モデルとなる 風速勾配は図5.1に示し， 風速勾配の大きさ0.1は実線で， 0.12は破線で， 0.14は一点鎖線 で表した. 対地速度uo=-8(mjs)， vo=30(mjs)， wo=O(mjs)で一定とし， 風速勾配を変化させ たときの状態量と制御量を図5.2 に示した. 本節の計算には ， DCNLP法を用い， 部分区 間数は 20とした. どの風速勾配についても速度の初期値と終端値が等しくな って境界 条件を満たしていることがわかる. また高度についても初期値と終端値が等しくなっ

て評価関数を最小にしていることがわかる. 制御量である揚力係数を見ると風速勾配 が大きいほど値の変化が少なくなっている. 風速勾配が大きいと径路の頂点付近で速度 が遅くなっても風速が大きいため 揚力係数を大きくしな くても容易に風下側へ旋回で きる. 風速勾配が小さいと径路の頂点付近で速度が遅くなると風速も小さいため， 風 下側へ旋回す るために揚力係数を揚力係数の制限まで大きくしている. 次に， 最適飛 行径路を図5.3 に示した. これを見ると， 始点で風下側へ出発した場合は， 風速勾配が 大きいほど径路全体が風下側になっている. 飛行時間は，風速勾配が0.1のとき10 .68(sec)，

0.12のとき9.7 7(sec)， 0.14のとき9.79(sec)であった.

図5.4には， 時間ごとのK.EとP.E.の値をプロットしたものと， 式(3 .7 4)を用いて計算し た空気力と風のエネルギ収支が示しである. K.E.とP.E.の関係を見ると K.E最小の時 に， P.E.最大という関係になっていることがわかる. 空気力と風のエネルギ収支のグラ

フの縦軸は高度， 横軸はエネルギを表す. 風上へ上昇中は，風から得るエネルギよりも 抵抗による仕事の方が大きいのでエネルギの値が負になっているが 風下へ降下中は

第5章 ダイナミック ・ ソアリングの解析結果 5.1

風から得るエネルギの方が抵抗による仕事よりも大きくなるのでエネルギの値が正に なっている. そして， 終点では， エネルギ収支を表す曲線が閉じて， エネルギ収支が釣 り合っていて， 風速勾配が大きいほど上昇中のエネルギ損失が少なくなっているのがわ

かる.

ιL_

10

-1・ 1 ・ 1 ・ 1 ・ 1 ・ ノ / / 1

111 1 ノ

/ /

/ /J / / / . / / / / . ， J / / / / J / / /J /

//' υy hv /γ

ダイナミック ・ ソアリングの解析結果

80

A抽Y ひてミ 屯

2 6 第5章

0.1 0.12 0.14

図5.1風速勾配のモデル(常ι= 0.1，0.12，0.14)

第5章 ダイナミック ・ ソアリングの解析結果

2 50

、 、 ー、豆

H 句-

、

、 、

、、 、 、 百d

J指2己.} 、h、、 ー 、ー、、: ^{.5 、 、}

、

、

、

、

、

江-2-0.4 O.ζ ⁿ 0.8

3叫よ司、

ノイ!;

O. --0ゥ t"l.d， ()瓜 () R 1

3k

^/-'-./

^rー!

.1叫 、\，\. ム会/

1.2トー 4

d

1υ吋ノナィ二二二_，_'二、i

0叫 \ / \、、〆/ 4

。も

^{02 04 0.6 0.8}

， ... ， ，

メ〆?�...>'"

�lf2 04 0.6 0.8

4� /-'?

官

仁

-5

0 2

0.1

0.4 0.6

0.4 O.6 0.8

一一一一一. 0.12 0.14

、 、ご

、

、

、

図5.2風速勾配を変化させた時の状態量と制御量(

守

= 0.1，0.12，0.14)

5.1.

第5章 ダイナミック・ソアリングの解析結果 ^5.1.

h(m) ��

nu qL nU A宝

χ(m)

0.1 一一一一一. _0.12

一一一一一一 _0.14

図5.3風速勾配を変化させた時の最適飛行径路(

常ι

=0.1，0.12，0.14)

第5章 ダイナミツク・ ソアリングの解析結果 _�

ミ11\

、 、

、、

、 、、

、ι、、、、ーー.-〆

0.6 08 t/iグ

ミ凶 . K

4 3

-2 o

ENERGYIω [)( 1 03] ²

0.1

一一一一一・ 0.12 0.14

図5.4風速勾配を変化させた時のK.EとP.Eおよびエネルギ収支(

常乙

⁼0.1，0.12，0.14)

第5章 ダイナミック・ソアリングの解析結果 _忌2ムー

5.2

飛行速度の影響

始点での合成対気速度をVo= 31. 0 5 ( m / s ) で一定にし，対地速度成分Uoの大きさを変え て始点での滑空機の向きを風下側， y軸方向，風上側にして， その影響を調べた 'Uo=- 8(m/s)， vo=30(m/s)，山O二O(m/s) (風下側)を実線で， uo=O(m/s)， vo=31.05(m/s)，ωo=O(m/s) (y軸

方向)を破線で， uo=8(m/s)， vo=30(m/s)，切o=O(m/s) (風上側)を一点鎖線で表し， 状態 と制御量を図5.5に，飛行径路を図5.6に，エネルギ収支を図5.7にそれぞれ示した. ただ し， y軸方向に出発する場合だけx(to)⁼ x(ij)二Oという拘束条件を加えて， ν軸上に戻っ てくるようにした. 本節の計算には， modified Simpson's system constl'aintsを用い， 部分区 間数は20 とした.

図5.5のz軸のグラフを見ると風上側へ出発すると終点も風上側になっており， ダイ ナミック・ソアリングで風上側へ飛行を継続していくことが可能であることがわかる.

風上側へ出発すると到達高度はy軸方向や風下側に出発した場合に比べて高くなってい る. 揚力係数については，風下側に出発した場合は，径路の頂点付近で揚力係数の制限 まで使っているが， その他の場合は制限までは使っていない. 図5.6の飛行径路を見ると，

始点で風下側， y軸方向， 風上側に出発すると終点が風下側， y軸上， 風上側になってい ることがわかる. 図5.7の風と空気力のエネルギ収支を見ると，風土側へ出発すると上 昇中のエネルギ損失が他の場合より少なく， 降下中に得るエネルギが他の場合よりも 大きくなっている. 風上側へ向かつて出発する方が他の場合に比べて， エネルギの点で

有利であることがわかる. 出発するときの方向や終点での位置の拘束により，始点より も風上側や風下側そしてy軸方向へ飛行を続けていくことができることがわかった. 飛 行時間は，風上側出発のとき10.93(sec)， y軸方向出発のとき11.2 9 ( sec ) ，風下側出発のとき

10 .92(sec)であった.

第5章 ダイナミック ・ ソアリングの解析結果 ιL_

..._._. _.ー- 、・、一

一/ ・一

^、

〆' 、

，" 、

， 、

" --ーー、

^、

;' /- 、 、

^司・、

一/ 今〆 一、 、、

〆 、 .、

、一

" � 、、、 、

ー・

./ I

^/" ... 司、

./ 〆/" "、 、

f 〆 / " ^、

" / " 、

〆 ， / '\. 、

ノ " / "\. 、

.� "/ " 、

10

音色しH

3

20ト " "ご子、、

ー司、

^I

，

^{I "、}

、

々、 J'/ "、 ?\ i

z_;:s 叶_ν\、ょ、γ"'-:、、\_、・、 -，'，

20

ミ10 ウ- へぎh亡、同

少ベ、\

5 4

0，1 O.

2

ヘ可、ミ与占

p、“

，、.d

B也芝句、shE

，&

。 ...:) U

t1tf

。

Uo=帽O.8(m/s) Uo=

O. o

(m/s) Uo=

O. 8

(m/s)

ルγ

図5.5UQを変化させた時の状態量と制御量(UQ₌ -8，0， 8(m/s)).

第5章 ダイナミック ・ ソアリングの解析結果

50 40

80 60

40

χ(m)

20

。

".白一一 今.・

.-・

.・ 司

白一一・' ・ '"

L" ・

，・・ ・

・ 0・ '.、h

.ム匂・、'h -"， .・一'_ ー ・ 、・-句・.

・. 、.

... ・-

.一 .

':.，

.. --，

・

.

一

-40 0

uo=-O.8(ml.ゆ

一一一一 一.

Uo= O.O(m/s)

一一一一-

Uo=

^{O. 8}

(m/s)

y(m)

図5.6

Uoを変化させた時の最適飛行径路(UQ

=

-8， O，8(m/s))

5.2.

50

;:h(m)

第5章 ダイナミック ・ ソアリングの解析結果 5.2.

[x 1 03]

詮同U時

O�8 [x 103]

を凶 . k -‘肉 、 、 、 、 、 . A崎 、 ・、

A 、 、山 1 ・

、

、

九 、 、 、

、

、 、、 i

〈べ

z hp bv

0.4--0.0

t/iグ ^。

5

O�ー←2 ENERGη'J) [x 103]

4

今3

『4

ヘ主主

Uo=・O.8(mj.砂 Uo= O.O(m/s)

Uo= O.8(mj.砂

図5.7tLoを変化させた時のK.EとP.Eおよびエネルギ収支(Uo= -8

J

⁰

J

8 _{( m} / ^s ) )

第5章 ダイナミック・ ソアリングの解析結果 L主一

次に， これまでは始点でのz軸方向の速度をωo=O(m/s)として計算してきたので， Wo を変化させた場合につい て， その影響を調べた. uo=-8(m/s)， vo=30(m/s)，ω0=5(m/s)を実 線で ， 'uo=-8(m/s)， vo=30(m/s)，山0=-5(m/s)を破線で表し ， 状態量 と制御量を図5.8に， 飛行 径路を図5.9に エネルギ収支を図5.10にそれぞれ示した. z軸は下向き正としているの でω0=5(m/s)の時は， 滑空機は下向きの飛行速度を持っていることになる.

高度hのグラフを見るとω0=5(m/s)の時は， 最初に初期高度h=5(m)よりも高度を一旦 下げてから上昇し，終点に近づいてh=5(m)よりも 下がることなく終点に到達している.

wo=-5(m/s)の時は， 最初に初期高度h=5(m)よりも高度を 下げることなく上昇し， 終点の 手前でh=5(m)よりも高度を一旦下げてから終点に到達している. ロール角を見ると始 点を出発したときはwo=5(m/s)よりもω0=-5(m/s)の方が ロール角が大きく， 終点に近づ くと wo=-5(m/s)よりもω0=5(m/s)の方が ， ロール角が大きくなっている. これは， 下向き の速度を持った状態から上昇に移るためにはロール角を小さくする必要が あるが， 上 向きの速度を持った状態から上昇に移る場合は， ロール角が大きいままでも上昇する ことができるからである. 図5.9を見るとω0=-5(m/s)の場合は， ω0=5(m/s)に比べて径路 が風上側になっていることがわかる. 図5.10のエネルギ収支を見ると， ω0=-5(m/s)の方が ω0=5(m/s)よ りも降下中にもらうエネルギが大きくなっている. 飛行時間は， wo=5(m/s) のとき10.45(sec)， ω0=-5(m/s)のとき12.11 (sec)であった.

第5章 ダイナミック・ソアリングの解析結果

ミ

_:::s

1.2ト八

。 0.2

5

;〆一一'\ 」

/

0.4 O.6 tIiグ

::^{h 『}。、

.5

、 J

、 ，

412

0.8

Wo= 5.0(m々j

一一一一一.

Wo=・5.0(m/s)

‘3 ‘ミ

、

、

、

、

、

、

図5.8ωoを変化させた時の状態量と制御量(ωo=5，-5(m/s))

5.2.

第5章 ダイナミック・ソアリングの解析結果

。

Wo= 5.0仰勾

一一一一一.

Wo=・5.0(m々ノ

図5.9 Woを変化させた時の最適飛行径路(ω0=5，・5(m/s)).

5.2.

lG _h(m)

第5章 ダイナミック ・ ソアリングの解析結果 ιL_

[xl04]

S同U叫

， ， 、

0.2 0.4 0.6 0.8 t/tf

15 y-- _、

、、、、

/

^{、、、、、}、、

、、 0.6 O.8

t/tf

5 4

、.ーa豆』、z、3

2

-1 -0-. -- 1

ENERGYI向) [x 103] 2

Wo= 5.0(m勾

一一一一一-

Wo=・5.0(m/s)

図5.10ωoを変化させた時のK.EとP.E.およびエネルギ収支(wo=5，-5(m/s))

第5章 ダイナミック ・ ソアリングの解析結果 5.3

5.3

ダイナミック・ ソアリングのまとめ

滑空機がダイナミック ・ ソアリングの飛行を継続するための条件を調べるために ， 飛 行径路の始点と終点でのエネルギが等しくなるような最適制御問題を考え， _DCNLP法 と前章で提案したmodi五ed direct optimization methodを用いて解を求めた. 空気力と風の

エネルギ収支は， 風上へ向かつて 上昇中は， 風からもらうエネルギ よりも抵抗による仕 事の方が大きいためエネルギの値は負になっているが， 風下へ降下中は， 風からもらう エネルギの方が大きいためエネルギの値は正になって最後にエネルギ収支のグラフが

閉じて， 滑空機でもエネルギ収支の釣り合ったダイナミック ・ ソアリングが可能である ことがわかった. ロール角を大きくして旋回しているときは， 揚力係数もそれぞれ大き くし， 特に飛行径路の頂点付近では， 速度が遅くなっているので揚力係数の制限一杯ま

で大きくしている. 風速勾配を大きくすると， 始点で 風下側へ出発した場合， 飛行径 路全体が風下側に移動する. 始点での滑空機の向きを風下側， y軸方向， 風上側へ向け ることで終点の位置も風下側 ， ν軸上， 風上側にすることができ， 径路全体の進んでい く方向を決めることができる. 始点で水平ではなく上向き， 下向きに出発した場合 ， 上 向きの方が下向きに比べ径路全体が風上側になっている. 滑空機が， 始点でいずれの向 きをとってもエネルギ収支の釣り合ったダイナミック・ソアリングが可能であることがわ

かった .

第6章 結論

Thermal-tcトthermalfii ghtのtotal fiight time を最小にする最適飛行径路は， まず， t'骨空機は，

出発後すぐに急降下し， 速度を速やかに増加させる. 増加した大きな 速度を維持して 飛行し， できるだけ次のth ermalに近い位置から急上昇を行って， 速度，径路角の終端条 件と一致させる. Lift co efficient rate contr olを行うことで， 滑空機がthermal内において， 最 小沈下速度時の揚力係数で飛行することができるようになった. ダイナミック ・ ソアリ ングを継続するための条件を調べるため， 径路の始点と終点でのエネルギが等しくな るような 最適制御問題を解いたところ， アホウドリと同様に滑空機でもエネルギ収支 の釣り合ったダイナミック・ ソアリングを行うことができることがわかった. また， 複雑 な最適制御問題を解くのに適している直接法として積分公式を使った解法は， 複雑な 最適制御問題を解くと制御量に不連続な部分が現れる. そこで， 制御量を線形補間す るmodified direct optimization methodを提案したところ， 不連続な部分が無くなり 精度も修

正前より も著しく落ちることもなく有効な方法であることが示された.

以上， 本研究により thermal-to-thermal ftightの最適飛行方法および滑空機によるダイナ ミック・ ソアリングが可能であることが示され， また， modi五ed optimization methodの提案 により，最適制御問題の直接解法の数値計算上の問題点が解決された.

今後は， これらの飛行方法に対して滑空機の姿勢も考慮してさらに検討を進め， 模 型滑空機を使って飛行実験を行い最適制御問題を解いた結果と比較を行う. また， 最適 化の計算手法についても， 直接法は間接法に比べ初期解の選定に関してロバストであ るが， やはりその解は，初期解の選定に大きく依存している. そこで，初期解の選定に，

大域解の探索に優れた遺伝的アルゴリズムを使用し， これによって得られた解をもとに して直接法を使い， 更に精度の良い解を求めることを考えている.

付録A

最適制御計算における評価関数の勾配

最適制御問題において， 操作量u(t)を与えると ， 式(2.15)， (2.16)よりx(t)が定まり， 式 (2.17)より評価関数 J の値が確定できる. した がって，Jはu(t)の関数であるから， J(u(t))

のように書ける. いまJに有限の値を与えるような関数u(t)の集合を考えると，最適制 御問題は この関数の集合の中からJに最小値 を与える要素を求めることであるという こと ができる.

この考えは関数空間での最小問題 という ことに導かれる. 関数空間は各u(t)をそれぞ れその中のl点として表すような空間であり， 無限次元の抽象的な空間である. このよ うな空間はヒルベルト空間 として定められ ， 勾配法， 共役勾配法などの公式がこれら

の空間にまで拡張される. ここに述べる勾配法 もこの空間での 勾配 を用いて お り， 以下 において実際の計算を示す. 7)，31)

u(t)の大きさに拘束の無い場合を考える. 探索法 の基本的なパターンは， J(x)を最小

化するのに初期値Xo と探索方向do と を適当に与え，J(xo +αdo )を最小にするαを求めてそ の値をαminと する. つぎに忽1 = Xo +αmind。を再び初期値として と り， かつ新た な探索方 向d1を何らかの方法で定めて， f(Xl +αdt)を最小にする. 以下この手続きを繰り返す.

関数の関数で あるJ(u(t)) の場合もこの基本的パターンは変わらない. まず初期値 として，何らかの時間関数uo(t)を予想する. また探索方向 を示す時間関数So(t)を定め ， J(uo(t) +α80(t))を計算する. これはl変数αだけ の関数になるのでJ(α )の最小化を行い ， αminを求め， U1 (t) = UO(t) +αmin80(t)と置いて， 以上の手続きを繰り返す.

この手続きにおいて，問題は探索方向80 をどう定めるかに あるが，勾配法においては，

J(u )のu = Uoにおける降下方向に と る. それは形式的には J のUoにおける勾配により

付録A 最適制御計算におけ る評価関数の勾配

nu-u一u

Hペ一3

ペσ一一一nu s (A.l)

と書ける. そこで， Jの勾配Juを求めること を次に考える.

f(x)の2二五における勾配fxを求めるとき，

g(t) = f(x + ty) - f(x) (A.2)

と置き ， ⁹ (t) をtで微分して ， t = 0とするとその結果は

ダ(t)It=o _{= f忽T(x) ν (A.3)}

となる. この式の右辺を見ると ， 勾配f æ(x)と五からのずれ を表すνとの内積になってい る. よってJ(u(t))の場合の勾配も以上のことと同ーの形式で考えてゆくことにする. し かし内積の定義は有限次元のものとは異なり， 2つの関数 æ(t)， y(t)の内積 を次式で定義 するものとする.

川 (収例tりt))₌

そこで， J(u(t))の勾配Juを求めるために ， u(t)が変分 η(t)によってu(t)+αη(t)に変化した ときの Jの増分

h(α) = J(u(t) +αη(t)) -J(u(t)) (A.5)

を考える. 式(A.3)との対比でゆけば

州la=O₌ (A.6)

となるようなJuを求めれば， そ れが勾配となるわけである.

h' (α) を計算する前に ， 微分方程式の拘束である式 (2.15)について考える. u(t)が与え られたときの式 (2.15)， (2.16)の解をæ(t : u)と書くことにし

αと(t)=æ(t:u+αη)一忽(t: u) (A.7)

と置くと， æ(t: u +αη)=æ(t:u)+αと(t)とu(t) +αη(t)とは式(2.15)を任意のαに対して満た さなくてはならないので

付録A 最適制御計算における評価関数の勾配

がαに対して恒等的に成り立つ. この両辺をαで微分し α=0と置くと 次の線形微分 方程式

ごい)⁼ f� (t， x (t ^{: u)， u} (t))ご(t)+fu(t，忽(t^{: u)，} U(t))Tη(t) ^(A.9) が得られる. さらにæ(t: u)， ^æ(t: u +αη)の初期値はともに式(2.16)を満たさなければな らないことから

ごい。)

=

0 (A.¹⁰)

である. ゆえに， u(t)とその変分η(t)が与えられると， x (t)の変分と(t)は式(A.9)， (A.10) から定まる. 式(A.9)は線形方程式である から， その遷移行列を4'(t，to)とすると式(A.9)，

(A.10)の解は

的)= か ^{川�(r， x(r}

^: u)， u(r))17(r)dア と書ける.

ところでJの増分に戻ると， 式(2.17)より

J(u(t) +αη(t))

=

g(tj，x(tj: ^u) +αご(tj ))

+jtf山川+必(t)州)+叫(t))dt となる. こ れをαについて微分し， α二Oと置くと

3gT(t.r，X(tJ: ^11.)) ^{r{J. \ I}

(1 r

_{åf6 (t} ， x (t : 11.)^{， 11.} ( ^t ) ) ゲ(α)1にO二 口 1. 、 と(tJ)+

I �

^。 _ご(t)

(A.ll)

(A.12)

Of6 (t ^， x ^{(t : u)， u} ( t )L _η_ _(t)トdt/J.\

1

_(A.13)

となる. 式(A.l1)を用いて式(A.13)中のと(tj )，と(t)を消去すると rtf ⁸gT(^tj

，

_忽(tf : u))�/J. J.\1T

ゲ(α)Iα=o

= l

^{n一九 1 e(tf，} t)f1u(t， x(t ^: u)， u(t))η(t)dt

+

〈fγfザ刈3川州江爪(什川ア

x

1:かp(rけ刊?パ刈tり)fι2

^川 :ι刈叫川川u叫札山川)ト川川?川刈 u叫州州州山川(収例附州tの引伽桝)リ川州)川川η州仰州川川(収例附tり付刷 )凶州d似 d計^T

1ff

^t勺1

8fl明州仰rれ刊(収川t

=

d ff川勺九 f

^{{ �げ3ωザ gT刊向Tη (}

η(収t)dt

p(収t^J， t刈tり)

fft

(いt，x (いt :u)，u(い例tり)リ) (ド例附^Tけ引)川)

t勺f ^{å灯伺州fだ爪fれ(}

^{巾川(什ア。 : 叫山川 'ρ刈u叫吋刷刷}p(r，t)dr.f�(t，忽(t: u)ρ(t) )

316 (^t，忽(t^{: u})，^u(^t))

1

8u -^{l' -， � II} r η(t) dt (A.14)

付録A 最適制御計算における評価関数の勾配

となる. 式(A.6)で示された形になったので， 勾配Ju は

r

^-'6:TI-l _L\ôg(tj，X (tj: u)) Ju = fu(t，x(t: u)ρ(t))

�

^pT(tf，t) V'y�L J

J

l δx(収t j )

+

t ρ tり))fJ)否iþ T川(什h卯川アh引，

^t

+ 10υu.(収t，忽(れt: U)， U(収tの)リ) (A.15)

となる. この式の右辺は式(2.18)を作ると き に導入した乗数ベクトルを用いるとより 簡単になる. 乗数ベクトルp(t)は式(2.20)， (2.2.1)の微分方程式を満たすが， u(t)を与え，

忽(t: u)を求めてこの方程式に代入したと きの解をp(t: u)とすれば， 式(2.20)， (2.21)はそ れぞれ

p(t: u) ⁼ -fx(t，x(t: u)，u(t))p(t: u) + fox(t，忽(t:u)，u(t)) (A.16) 3g(tf， x(tj : u))

P(tf: u) ^-

- - '1\-J�-:_/� - J

\ '

^-;; _(A.17)

ôx(tj)

となる. 式(A.16)と変分方程式(A.9)を見ると， 強制項は異なるがp(t: u)は式(A.9)の双対 の線形微分方程式を満たしていることがわかるので， 式(A.16)， 式(A.17)は式(A.11)で用 いた同じ選移行列iþ(t，to)を用いて解くことができる1 ただし，この方程式は終端の条件 式(A.17)が与えられているので， p(to : u)の値を仮定して式(A.16)を解くと

州: u ) = d M P ( t o : u ) +

J〈ケケ :)〉 ed Tη川(什h仰川 í，引，

^t

となる. ゆえに

(A.18)

附f:u) = dMP(to u)+

f

^{d(ア山(吋(ア:u)ρ(ア))dí (A則}

的o

この式を式(A.18)に代入すると， 次式となる.

時叫=ぜい.r，t)p(tj:叫-

j

^JhW判:

u)， u(í))dí

この結果と式(A.15)を見くらべると

Ju = -fu(t，x(t: u)，u(t))p(t: u) + flω(t， x(t : u)， u(t)) 式(2.18)を用いると ， 勾配は結局， 次ように表すことができる.

Ju = -Hu(t，x(t: u)，u(t)，p(t: u))

(A.21)

(A.22)

(A.23)

付録B

機体諸元と使用した数値

B.l

機体諸元

計算に使用した機体は 日本飛行機株式会社の「日飛ピラタス式B4-PCllAF型」で

ある. 35)図B.1に機体諸元， 三面図， 図B.2に滑空性能， 滑空性能曲線を示した.

その他， 使用した数値は

である .

空気密度 重力加速度 機体重量

P二0.125(kgs2/m4) g二9.8(m/s2)

mg ^{= 350(kg)}

(B.1) (B.2) (B.3)

B.1.

機体諸元と使用した数値 付録B

一冗諸

体

機 図

E。.凶-

面

機体諸元と三面図

町1

m

m

m Z

NACA _{6 43事618}

16.0 O. 4 0 10 3 0' 10 00

m

1. 0 7 m 0.425 m 0.9 37m

] 5.0

1 5. 0 1 4目l 6. 5 7 1，5 7 rtJ

長

翼 根 翼 端 笠力平均翼弦(MA C) 輿 魁

縦横比 先細比 取付角 上反角 操り下げ 高

翼 翼 巾 面 積 潟弦長

ふ宇-bエ主

全

m

nf 000匂

内Junu A

md“

図B.l 3.0 8

1. 7 1 NACA

1. 2 2 NACA 00 - 30 00

1. 3 2 翼

翼 尾巾積地

尾 全 高 面 積 翼 抱 取付角 取付角 上反角

直 水 平翼面翼

垂

付録B 機体諸元と使用した数値

滑空性能 〈最大重量 350 k9 のとき〉

hJ: l\t ftt空比 i止小沈下速度

0."

O.�

1.2

沈 1.6

下 速 度 (π〆sec )

m空速!庄(IAS)Knぺ1

8 5 7 2

作}空 比

3 5 : l 3:J 1

滑空性能曲線

60 80 100

比下速度m/scc

O. G 7 0.6 t1

滑空速度(Km/h )

120 140

図B.2 i骨空性能と滑空性能曲線

B.l.

付録B 機体諸元と使用した数値 B.2.

B.2

抵抗係数の決定

抵抗係数は， 次式のような形で与えた.

GD ⁼ CDO + Kc

i

K二 1 1T'eA

(B.4) (B.5)

CD 抵抗係数

GDO 有害抵抗係数 GL 揚力係数 A:アスペクト比 e 飛行機効率

eについてのデータがないので， 滑空機はe⁼ 0.9程度の値であると仮定してKを求め ると

J{ ⁼ 1 1

(B.6) (B.7) πεA 3.14 _{x 0.9 x 16}

= 2.21 _{X 10-2}

となった. 次にCDOを決定するために次式を用いた. 36)

(合)

^{min ，ink二}

ムム

^{。 (B.8)}

最小沈下速度の時の揚抗比の値は， j骨空性能の表より

(ま)

^{min ，ink (B.9)}

である. よって

CDO 3 _nu

-- 1『ム 司Ei

B 16 ^x 2.2_{1 x 10-2 X} 33₂ B

- 7.79 ^X 10-3

となった.

以上より， 抵抗係数は， 次式のよ5になった.

CD ⁼ 7.79 ^X 10-3 + 2.21 _{x 1o-2c}

i

^(B.12)

付録B 機体諸元と使用した数値 B.3.

B.3

Thermal圃to-thermal flightの初期値の決定

滑空機は， thermal内では， 高度を獲得するために最小沈下速度になるよう に飛行 を行う. そこで， 先ほど示した滑空性能の表から飛行速度と径路角を求めると

V(O) = 72(km/h) = 20(m/s)

γ(0) = - tan ^-1 (1/33) = -0.0303(rad)

(B.13) (B.14)

である .

また， その時の揚力係数は

mg cosγ= _L = 1/2pSV2CL (B .15) の関係から

m9 cosγ(0)

CL(0)= =0 99

pSV(0)2

、、‘，，FμU 噌BiB

である.

以上の数値を用いて計算を行った.

付録C

DCNLP(Direct Collocation with N onlinear Programming)

C.l

3次補間式の導出

最初に部分区間のみに注目して， segrnent length Sε[0，1]を考える. このsegrnent length を時間のsegrnent length Tに変換する. 時刻tとsegrnent length SとTの関係は次式のように なる. 22)

一一円bT

(C.l) (C.2) (C.3) TdS = dt

T = ^{一 -}dt dS

微分方程式とsegrnent length Sとの関係は次式のようになる.

ピ = f (C.4)

dx = f

(C.5) dt

dx dS = f

(C.6) dS dt

dx dt

(C.7) d5 = d5 f

土 = Tf (C.8)

(C.8)式の関係は， 後の補間式の導出で必要になる.

Segrnent length 5を使って， 3次の多項式を次のようにおく.

付録C DCNLP(DIRECT COLLOCATION WITH NONLINEAR PROGRAMMING)

^Cム

端点で次式のようにおく.

x(O)

_{= X1，}

x(l)

_{= X2}

生(0) こど1 生(1)= X'2

d5 ^{\ � I}

-

�.l ， d5

(C.10) (C.ll)

zの多項式をSで微分すると次式のようになる.

土(5)= C1 + 2C25 + 3C352 + 4C453

(C.12)

端点と中点でのX，土をそれぞれもとめる.

5=0のとき

x(O)

= Co = X1 土(0)⁼ C1二土l

(C.13) (C.14)

5=1のとき

x(l) =

Co + C1 + C2 + C3

土(1) =

C1+2C2+3C3

、、EE，，，、、，E，，，

vhu ρnu 唱'i 唱EA 円U 円U

，，aEE‘、，，，.‘、、

5=tのとき

11\ _ 1_ 1_ 1_ 1

川

\2 )

汁=Cn^{V '} +�C，+�Cっ+ �Cq

2

^�

4 “ 8

V • + ^�C4 ₁₆ ' 以上をまとめると 次式のような行列式になる.

(C.17)

。 。

:||:

。

1

。

X1

1

I

^C2 一

。 2

₃

J l

C3 X2

L

係数Cは次のように求まる.

「 Co

。 。 。

11 ::

C1

。 。 。

C2 C3 -11

18

-4 ⁵

-5 14 -3 l l l I I

^ZX2 ²

」

Sにcollocation pointの値を代入する.

s =

�.

_{52 =}

!.

53 =一

1 2'� 4'� - 8

(C.18)

(C.19)

(C.2

0)

付録C DCNLP(DIREC T COLLOCATION WITH NONLINEAR PROGRAMMING) ♀よー

係数Cとcollocation point Sの値を(C.9)式に代入して整理すると次式のようになる.

れ=j(町村2)一� (土1一れ)

FU Jt，‘、 nL 、、EE，，，1s・・4

ここで， segment length ^T =ムむとして(C.8)式の関係を使うと

Xc二;(zt+町村)-iムむ(fi - fi+l) (C.22)

付録C DCNLP(DIRECT COLLOCATION WITH NONLINEAR PROGRAMMING) ζL

C.2

3次補間式の微分値の導出

Segrnent length Sを使って， 3次の多項式を次のようにおく.

x(S) = Co

⁺

C1S

⁺

C2S2

⁺

C3S3

^(C.23)

この式の微分値をとると次式のようになる.

土(S)

⁼

C1 +

^{2C2S +} 3C3S2十4C4S3 (C.24) 上式の係数CとS二1/2を代入し整理すると次式のようになる.

t(i)=-;(Zl-h)一j(21+Z2)

^(C.25)

ここで， collocation pointを図4.1の表現に書きかえ， T=ムむの関係と(C.8)式を使うと， 次 式のようになる.

z;ム14=-;(Zt-h)-j伽ん1)ムti 3(Xi -Xi+I)

=一 2ムti 一一(九十点+1)

⁴

(C.26) (C.27)

付録D

非線形計画法

D.l

逐次2次計画法(successiveまたはsequential quadratic pro田 gramming method)

評価関数が2次式， 制約式がl次式である問題は2次計画問題と呼ばれる. 2次 計画問題はキューン・タツカ一条件を利用すれば， 線形の問題に変換されるため， シンプ

レックス法を拡張することによって解くことができる.

この章では非線形計画問題に対して現在最も有効な方法のーっとされている逐次2 次計画法(SQP法)を取り上げる 19)，25)これは各反復において元の問題を近似した2次 計画問題を逐次解いてゆく方法であり， 制約なし最適化問題に対する準ニュートン法の 拡張と考えることができる.

D.2

最適性の条件

次のような非線形計画問題を取り扱う.

目的関数 J(x)→最小 (D.l)

制約条件 Ci( x) = 0) i = 1， 2，・•，me

付録D _{非線形計画法} D.3

以下では， 目的関数fとすべての制約関数Ciをη次元空間Rl1で定義された2団連続微分 可能な関数と仮定し 等式制約条件と不等式制約条件の添字集合をそれぞれ

E = {1， 2，' "， me}， _J{ = {me ₊ 1， me ₊ 2，" " m}

で表す.

fを問題(D.l)の局所的最適解とする. そのとき， Ci (x*)二Oが成立している制約条件，

すなわち， x*における有効制約条件(active constraint)の勾配ベクトル マCi(x*)がl次独立 ならば， 次のキューン・タッカ一条件を満たすべクトルの組(x，U)= (x*，u*)が存在する.

マJ(x)-

L

_{Ui マ仰)=0}

i^=l (D.2)

(D.3) (D.4) Ci (x) ⁼ 0， _iεE

Ci (x)三0，同三0， Ui Ci (x) = 0， ^{ï E} J{

ただし， Uニ(Ul川2，・• ，um)Tはラグランジ、ユ乗数のベクトルである.

D.3

制約っき問題に対する準ニュートン法

この節では， まず次の等式制約のみを含む問題を考える.

目的関数 制約条件

J(x)→最小 C(x) = 0

にdD ''a，‘、

ただし， ^c( x)はベクトル関数C( x) ⁼ (C 1 ( X ) ， C 2 ( ^X ) ^， • .. ， Cm ( X ) ) Tを表すものとする. 問題(D.5)に 対するラグランジ‘ユ関数(Lagrangian)を

L(x， U) = J(^Xト

L

_UiCi(X) ^{D ハhu}

i=l

で定義し， この関数のzに関する勾配とヘッセ行列をそれぞれ

マxL(い)二 引い)-

L

^Ui \7 Ci ( X ) (D.7) i=l

\7

;

L(い) = マ2J(X)-

玄

Ui マ2 Ci (X) D 口δ

i=l

と書く. また， ベクトル関数cの点zにおけるヤコビ行列(J acobian matrix)， すなわち偏微 分係数8Cj(X)j8xiを第(i， j)成分とするようなn x m行列をマC(X)で表す.

付録D 非線形計画法 ^D.3.

さて， いま点X(k)が与えられているものとして， 次の反復点X(k+l)を求めるため， 次の ような2次計画問題を考える.

目的関数 制約条件

マバz(k))Td+

;

^{JB川→最小}

C(X{k)) +マC(X{k))T^{d = 0}

ハ同υ

D

これは， 問題(D.5)の目的関数をx = x(k)において2次関数で近似(ただし定数項れl;(k)) は無視)し， 制約条件をl次近似した問題であり， ベクトルdは現在の点X(k)からの変位 を表している. また， 問題(D.9)の目的関数の2次の係数行列B(ん)は， 後述のように， ラ グランジュ関数Lのヘッセ行列を近似した行列である. このような2次近似問題を解い て次の反復点を求めるという考え方は， 制約なし問題に対する準ニュートン法と共通 するものである.

2次計画問題(D.9)に対するキューン・タッカ一条件は

マf(X(k))+ B(k)dーマC(X(k))U= 0 C(X(k))+マC(X(k))Td二O

、、，E，，， 、、，E，，，

nU 司i 14 1i

D D

'''s‘‘、 ，，a'E、、

となるが， d = x(k+l) - X(k)， U = 'u(k+l)とおいて式(D.I0)(D.11)を変形すると

、11ノ 円4 1i

D

，，』'ー、

111111Eli--J

lk u

ぃ川 ル刊 が

「11111111』

マ

z ru Pし Jtt、

ルけ

一一

「SI--lili--」

k k z u 一 一

tE占 有ム

+ + 'K LA Z U

「11111111』 「11111111」 、、‘『，，，

ak z

d、

0

マ

T 、lk

'k rー

が Fill--till-L

マ

pu い

が得られる.

他方， 問題(D.5)は不等式制約を含まないので， そのキューン・タツカ一条件は， つぎの ようなn+m個の変数と同数の式からなる非線形連立方程式の形に表現できる.

マxL(x，^u) ^{= 0} c(x)=O

(D .13) (D ^.14)

この非線形方程式にニュートン法を直接適用すると， その反復公式は

「till--IEl--」 、、，BE，J

)

』厄

( u

、、・1''

)

、け， ι h J

(

fkz zfl

fjt、 fu ru z

マ

FIllll1111L 噌E4 「titti--tlBI』」 、、‘.，，， 、lJ 'k ra，、

z

，，，Et、 戸し nHU

マ

、、，，，，，

)

'K

パ必

T

n，‘ 、、t，，J

' )

、IJ 'k lκ rl、

( C Z

づJV ruu 〆tt、 Fしv

{ヤ

、‘，，， 、、，，， 『IllIll111J 「Ill11111L

マ

ヮ“z .可 IK 1k

，，E‘、

，，.‘、

z u

F1111111lL 一一 「l'all-SIBI--」

)

句aゐ

)

唱EA

+ + 'k lk

，，I、

rl、

z u

FEti--lElsi--』 戸hd 噌』ム

D 「illIllit--d、、I，，JEK U 刈

伏、

が

F'lliti--111」

マ

z 引

c い

一一

「Ill111111114

k k z u 一 一 + + 'κ

』κ z u

Fit-E11111Etl」「It--11111J

'k z

d、

o ' マ

が一マ。あ

川

'

で

りげ

川/

11 k

'' jK

IL

l t が い 』r- f t

k n

しハUHマ ー マ あ る

BV 斗& r--B1』1111L rli =口

付録D _{非線形計画法} D.3.

となる. 式(D.12)と式(D.15)を見比べると， 式(D.12)は式(D.15)のヘッセ行列マ;L(X(k)， U(k)) をB(k)で置き換えたものであるから， もし実際にB(k)ニマ;L(X(k)， U(k))が成り立っている ならば， ニュートン法の反復(D.15)は， 2次計画問題(D.9)の最適解とそれに対応するラ

グランジュ乗数を計算することと完全に等価であることがわかる.

以上の考察から，適当な準ニュートン法の公式を用いて， B(k)がヘツセ行列マ�L(x(ん}J U(k}) の良い近似となるように更新してゆけば， 2次計画問題(D.9)の解を逐次計算する反復 法は， ニュートン法に似た性質を持つことが期待できる.

そのような行列B(k}の更新法としては， 例えば

S(k) = X(k+l) _ X(k) 、、‘，，，，、、B-a，r

円hU 司，t 噌『よ 14ム D D ( y(k) =マxL(X(k+l)J U(k+l))ーマxL(X(k)， U(k+l)) (

によってベクトル5(k)， y(k)を定義し， これらをBFGS公式に代入することにより， B(k+l)を 定めることが考えられる. しかし， 制約なし最適化に対する準ニュートン法の考え方を 制約付きの問題に拡張する際には， 少し注意が必要である. 特に， 2次計画問題(D.9) の解を計算するという立場からすれば， 行列B(k)が常に正定値に保たれることが望ま しいが，

(y(k))T S(k) > ⁰ _{(D .18)}

が成り 立たないときには， たとえB(k)が正定値であっても， B(k+l)は正定値になるとは限 らない. 制約なし問題の場合には， ふつう式(D.18)が成立するが， 制約付き問題の場合 には， 式(D.16)， (D.17)のベクトルS(k)J y(k)が式(D.18)を満たさず， その結果BFGS公式でB(k) を更新したときB(k+l)の正定値性が失われるという事態がしばしば起こる.

式(D.18)が満たされない場合にもB(k+l)の正定値性を保証するには， 準ニュートン法の 更新公式を次のように修正すればよい. まず

I

1， (s(k))Ty(k)三O.2(S(k))TB(k)s(k)のとき

。二�

O.8(s(k))TB(k)s(k)j{(s(k))TB(k)s(k) - (s(k)fy(k)}，

| _{それ以外のとき} とおき， この8を用いてベクトル

、、.BEE'e''n同U1i D ''ea--‘‘、

。(k) = ()y(k) + (1 _ ())B(k) s(k) ，，aEI、 D nu nノ白 、‘BE，，，

を定義する. そしてBFGS公式においてy(k)をÿ(k)で置き換えた公式 B(k+l) = B(k) +。(k)(ÿ(k))T B(k) s(k) (s(k))T B(k)

(ÿ(k))TS(k) (s(k))T B(k)s(k) D 噌Eょっム