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

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

Kyushu University Institutional Repository

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

河邉, 博康

https://doi.org/10.11501/3154824

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

権利関係：

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

河遺 博康

目次

1 序論

1.1 研究の背景と目的 1.2 本論文の構成

1

2

2.9 運動方程式の再定式化 2.10最適制御問題の再定式化 . _• 2.11 Li仇Coe伍cient Rateの決定

2.12拘束条件がある場合とない場合の比較 2.13 (dCL/dx)mαzの影響

2.14 VTの影響 .

2.15 Thermal-to-thermal flightのまとめ

4 4 6

8

8 0 1 4 9 9 0 3 4 9 3 7

ーム 1i 可i 1i 1i nL 円ノU η4 ηL 円J qu 2 Thermal-to・thermal flight

2.1 Thermal-to-thermal flightの研究目的

2.2 Thermal-to-thermal flightの運動方程式 2.3 評価関数の定式化

2.4 最適制御問題の定式化

2.5 共役勾配法による最適制御計算法 2.6 Thermal-to-thermal flightの最適制御問題

2.7 計算結果と考察

2.8 Li乱Coe伍cient Rate Control

3 ダイナミック・ ゾアリング

3.1 ダイナミック・ソアリングとは

00 00 つiu qd

目次

3.3 運動方程式の無次元化の必要性 3.4 無次元運動方程式の導出 .

3.5 ダイナミック ・ ソアリングの最適制御問題 3.6 ダイナミック ・ ソアリングのエネルギ収支

ハhυ ハhU 口問U H凶d A吐 Aせ A斗A A斗・A

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

4.1.3 DCNLP (Direct Collocation with N onlinear Programming) 4.1.4 Quadrature Rules

4.1.5 System Constraints 4.1.6 Quadrature Rule Scheme

4.1.7 それぞれの手法の精度の比較 4.1.8 精度を比較した結果

'L 守i 1i 円，&

円4

5 5 5 5 5 4.1 Direct Optirnization Method

4.1.1 一般的な最適制御問題

4.1.2 終端時間自由問題から終端時間固定問題への変換

5.5 55

4.2 _Modi五edDirect Optimization Method

4.2.1 計算結果と問題点

4.2.2 制御量の線形補間法

nu -- 4A nu nu qJ ハO ρO PO 円i 月i 司t

4.2.3 Modified Direct Optimization Method 4.2.4 修正後の計算結果

4.2.5 修正後の精度の比較. _•

4.2.6 Direct optimization me出odのまとめ

75 円O nU Aせ 可i nHU nHυ

5 ダイナミック・ ソアリングの解析結果 5.1 風速勾 配の影響

5.2 飛行速度の影響

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

VD Ea -- ny oo nRu

nuuod

6 結論 100

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

B機体諸元と使用した数値 B.1機体諸元

105 105

目次

B.2抵抗係数の決定 ^{• • • • • • • • • • • • • • • • •} • • • • • • • • • • 10 B.3 Thermal-to-thermal ftightの初期値の決定 ^{• • • • • • • • • • • • • • • • 109}

C DCNLP(Direct Collocation with No凶near Prograluming) 110 C.l 3次補間式の導出 . ^{• • • • • • • . ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ー 一 110} C.2 3次補間式の微分値の導出 . ^{• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 113}

D非線形計画法 114

D.1逐次2次計画法(successiveまたはsequential quadratic programmi時method) ^{. .} .. ¹¹⁴ D.2最適性の条件. ^{• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 114} D.3制約っき問題に対する準ニュートン法 .^{• • • • • • • • • • • • • • • • • • • • • • • 11.5}

E Gauss-Lobatto Quadrature Rules 119

E.1 Gauss-Lobatto Quadrature Rule ^{. . ・ ・ ・ ・ 白 . • . . • . • . • . . • . • • . . • . • • . . . • 119} E.2 Fourth-degree Gauss-Lobatto Ruleの補間式の導出 . ^{• • • • • • • • • • • • • • • • • • • 12:3} E.3 Fourth-degree Gauss-Lobatto System Constraintsの導出 ^{• • • • • . • 125} E.4 Fifth-degree Gauss-Lobatto Ruleの補間式の導出 ^{• •• • • • • • • • • 127} E.5 Fifth-degree Gauss-Lobatto System Constra.intsの導出 . ^{. 目 . . .} . . . 130

第1章 序論

1.1

研究の背景と目的

滑空機が行うcross country soaringとは， 滑空機が周り の大気からエネルギーを取り出 し， 長い距離に渡って飛行を継続することである. その方法には， ridge soaring， wave soar­

ing， thermal-to-thermal flight， dynamic soaring (ダイナミック・ソアリング)等の方法がある.

Ridge soaringとは， 強い水平風が， 山の尾根に沿って上向きに曲げられ， その上昇風を利 用して， 山脈と平行に飛行する方法である. vVave soaring は， 山の風下側にできるlee wave と呼ばれる風下波を利用する方法である. Thermal-to-thermal flight は， therma.lと呼ばれる 熱上昇気流を探し出し， 高度を獲得するためにthermal内で旋回飛行に移り， 螺旋を描き ながら上昇を行う. そして， 雲などの状況から判断して次のthermalを発見したら， 次の thermalまで滑空飛行を行い， 滑空飛行によって失った高度を次のthermal内で再び螺旋を 描きながら上昇することで回復させる方法である. ダイナミック ・ ソアリングは アホ

ウドリの様な海鳥が海面付近で風速勾配を持った風に対して， 風上に向かつて上昇し上 空で低速旋回した後， 風下に向かつて急降下し海面すれすれで高速旋回を行うという ものである. こ れらの方法の中で， 従来からよく行われてきたthermal-to-ther・mal flightに ついて， 更に最適な飛行方法を見いだしたり， ダイナミック・ソアリングのように， これ まで行われることの無かった飛行方法を滑空機に応用することができるならば， 航空ス ポーツ界に対して多大な貢献をすることができ， 生物の飛行の原理を解明する上でも 大変意義のあるものとなる. そこ で本研究では， 最適制御理論を使ってもhermal-to-therrnal flihgt の最適飛行方法を求めるとともに， 滑空機でダイナミック ・ソアリングを行う可能 性について検討を行う ものである.

第l章 序論 1.2.

もうひとつ， ここで必要となってくる重要な研究課題は， これら滑空機の最適制御問 題を解くための解法の研究である. 最適制御問題を解くための解法には大きく分けて 2つの解法がある. 変分法から導出される多点境界値問題を解く間接法と， 最適制御問 題を非線形計画問題に変換して解く直接法がある . 間接法で状態量不等式拘束条件が あるような多点境界値問題を数値的に解くには 特別な工夫1) ，2)が必要で困難を伴う場

合が多い. また， 初期解の選定に関して敏感で， 問題に対する十分な考察に基づいて 初期解を与える必要がある. 一方， 直接法は， 時間的に連続な系を離散化し， 状態量や 制御量を非線形計画問題の最適化の変数に直接置いて， 状態方程式を補間式や積分公 式に変換して非線形計画問題の等式拘束条件に取り込むという方法である . 直接法の 利点は， 拘束条件の取り扱いが非常に容易で， 初期解の選定に関しても間接法に比べ てロバストな事である. 近年， 複雑な最適制御問題の解法として注目されている直接 法は， 本研究のような滑空機の最適制御問題の解法として適していると考えられる.

そこで， 直接法の中でも， よく使われているDCNLP法やquadrature rule schemeについて 比較検討を行いその性質を明らかにし， さらに数値計算上の問題点を解決したmodlfied direct optim ization methodの提案を行なう. また， thermal-to-thermal flightに対しては， 変分 法による解法で解くために， どのような定式上の工夫がされているかを示し， さらに 複雑で変分法による定式が困難なダイナミック・ソアリングに対して直接法を適用する ことがいかに有利であるかを示し， modified direct optimization methodの有効性を示す.

1.2

本論文の構成

本論文の構成は以下の通りである.

まず， 第l章は序論で， 本研究の背景及び目的について述べ， 各章の内容について概 説する.

第2章では， 滑空機のthermal-tかthermal flightの最短時間飛行問題を最適制御問題とし て定式化を行い， 変分法に基づく間接法によって解きその結果を解析する. 間接法で複 雑な最適制御問題を解く際に， 問題を簡単化するために様々な工夫が必要であること を示す.

第3章では， アホウドリが行うダイナミック・ ソアリングと呼ばれる飛行方法を滑空 機で行うために運動方程式を導出し， 飛行を継続するための最低条件を求めるために 最適制御問題として定式化を行う. ダイナミック・ソアリングのエネルギ収支の解析に

第1章 序論 1.2

第4章では， 直接法としてDCNLP法とquadrature rule schemeについて， 解析解のわかっ ている例題を解いて， 精度， 計算時間などを比較し手法の性質を明らかにする . さら に， ダイナミック・ ソアリングのような複雑な最適制御問題に直接法を適用したときに 生じる数値計算上の問題点を解決するために， modified direct optimization methodを提案

し， その有効性を示す.

第5章では， 改善された直接法を用いてダイナミック ・ ソアリングの最適制御問題を 解いて滑空機によるダイナミック ・ ソアリングの可能性について解析を行う.

第6章では， 今回の研究の総括を行い， さらに， 今後の課題について述べる .

第2章

Thermal-to-thermal flight

2.1

Thermal-to-thermal flightの研究目的

図2.1にthermal-もo-thermalflightの飛行径路を示した. 本章の研究の目的としては， 距離 Xfで離れた2つのthermal内のA地点からB地点まで飛行速度Vで滑空飛行を行う時間 と， 滑空飛行で失った高度YをB地点からC地点まで上昇速度竹で上昇することによって

回復させるのに必要な時間とを合計したもotal flight timeを最小にする滑空機の最適飛行 径路を求めることである.

ここで上昇速度 竹とは， thermalの上向きの速度と滑空機の沈下速度との差である.

そのとき， 滑空機は高度獲得のためにthermal内では最小沈下速度で飛行を行っている とする. また， 図2.1では， thermalに幅があるように描かれているが， thermalは線とし て考え， 滑空機は上昇速度時でB地点からC地点まで螺旋を描かず， まっすぐに上昇を 行うとする. また， 2つのもhermal聞の距離は， 十分離れていると考える.

第2章 THERMAL-TO-THERMAL FLIGHT

A

Y

Xf

図2.1 Thermal-to・thermalflightの飛行径路.

V : 飛行速度

Xf : 2つのもhermal聞の距離 y 高度

C

B

斗

竹 (thermalの速度)ー(滑空機の沈下速度)

2.1.

VT

第2章 THERMAL-TO-THERMAL FLIGHT 2.2.

2.2

Thermal-to-thermal flightの運動方程式

運動方程式を導くにあったって， いくつかの仮定を行う. まず， 重力加速度9， 空気密 度pは一定とし， 滑空機は質点として考え， X-Y平面内で運動すると考える. そして，

thermal内での運動は， 一定の垂直方向の速度怜で上昇を行うとする. Thermal-to-thermal

fiightの座標系を図2.2に示す. 図中で， Lは揚力， Dは抗力，Vは飛行速度，γは径路角，Y は高度， Xは飛行距離， 9は重力加速度， mは機体質量を表している.

運動方程式は次のように書ける. ^3)�6)

dV dt -D (V， CL)/m -9 sinγ dγ

一 [L(V， CL)/m - 9 cosγ]/V dt

dY dt ^一 Vsin γ dX

dt V cosγ

(2.1) (2.2)

(2.3) (2.4)

終端時間固定の最適制御問題の定式で解けるように 時間tを飛行距離Xに置き換える.

そのため式(2.4)で式(2.1)， (2.2)， (2.3)をそれぞれ割ると， 運動方程式は

のようになる.

， dV

V' = ^:� dX = -[D(V， CL)/m + 9 siw-y]/V cosγ γ'-生 = [L(V， CL)/(mcosì) - g]/v2 dX

y'=一一dY 二 一一γ

dX ^{一一 ，}

(2.5) (2.6) (2.7)

第2章 THERMAL-TO-THERMAL FLIGHT _2.2.

L y

V

X D

mg

図2.2 Thermal-to-thermal ftightの座標系

第2章 THERMAL-TO-THERMAL FLIGHT 2.3.

2.3

評価関数の定式化

式(2.8)の評価関数においてt 右辺第1項目は2つのもherma.l聞の滑空時間を表し， 右 辺第2項目はthermal内でもとの高度に到達するまでの上昇時間を表す. そして， それら を合計したtotal flight timeを最小にする.

J =

Afhf)d川+川市[げY

^(O)

一

印川川川fバ刈)リν]ν/

=

J _一一一 dX

^l

f: f ) . 二f品-UXftmγdX

=

[' [じ:J�1ldX

(2.8) (2.9)

(2.10)

、‘，.，，，噌li144

fa'E、 q，b

終端での境界条件は，

V(

X

^{j) = V(O)}

γ(Xj) =γ(0)

^{(2.13) (2.12)}

であることが望ましいが， 厳密に満足することは要求せずに評価関数の中に終端値 評 価として次式のように表現する.

X ，α

匂I一 FJ n 一γ

二 S一ω ル v 一h 匂lFIll--L m丘昨 一X 作fl。

+ + x

nu x v 一一7，d v

(2.14)

2.4

最適制御問題の定式化

最適制御問題を変分法を用いた方法で解くと， 乗数ベクトルについての微分方程 式が導入されて， 状態方程式と連立した2点境界値問題となる. この2点境界値問題 を解けば， 最適軌道や最適操作量が求められるが， 多くの場合にそれは困難な問題と なる. したがって， 種々の数値 解法が工夫されている.

そこで最も よく用いられている方法は勾配法によるものであり， これは変分法を無 次元空間 での最小探索法と見て， 勾配法， 共役勾配法をヒルベルト空間などの抽象空 間に拡張し， この結果を最適制御問題に適用するものである. この方法ではまず最適

第2章 THERMAL-TO-THERMAL FLIGHT 2.4.

乗数ベクトルの微分方程式を終端から逆時間について解いて， これらの結果より評価 関数の勾配を計算して， 操作量の予想値を修正するものである. この方法は， 主として 操作量に拘束のない問題について適用される.

本章で考える最適制御問題をまず定式化しておく. 7)，8)ここでは， 状態量には拘束条 件がつかない次の問題を考える.

微分方程式

æ (t) = f(t， æ(t)， u (t)) (2.15)

を満たし， 初期条件

x (to)ニ忽O (2.16)

を満たす忽(t) = [Xl(t)，X2(t)， スn(t)]T，u(t) ⁼ [Ul(t)，U2(t)，... ，U，.(t)]Tのうちから

J二9MMH

j : ^:

^{山(t)， u (t)) dt (2.17)}

を最小にする最適軌道互い)，u(t) を求める. このために最適性の必要条件から， 乗数ベク トルp (t)= ^{[Pl (t)}，p2(t)，..' ，Pn(t)]Tを導入してハミルトン関数を次のように置く.

H = -fo(t， x， u) ₊ pT (t)f(t， x， u) ，，z'E、、 nJム 可ti 、、‘‘『，，口δ オイラーの方程式， 横断性条件から

P(t)= - BE åH

二 一fXT (t，忽 (t)，u(t))p(t) + fω(t， x (t)， u (t)) (tf)= -

J

_8x(tJ) L=Pf

rJ

となり， またuについての拘束条件がないときは 8H

åu

(2.19) (2.20) (2.21)

(2.22) 拘束条件のあるときは， Hをuが最大化するという最大原理の条件より， uを求めて

u(t) = h (t，x(t)，p(t)) (2.23)

第2章 THERMAL-TO-THERMAL FLIGHT 2.5.

のように解ける.

式(2.23)のu(t)を式(2.15)， 式(2.20)の右辺に代入すると， 式(2.15)， (2.20) はx(t) とp(t)につ いての連立微分方程式で， 式(2.16)で初期値が， 式(2.21)で終端値が与えられた2点境界 値問題 となる.

以上の問題に対して共役勾配法による計算アルゴリズムを考える.

2.5

共役勾配法による最適制御計算法

共役勾配法のlつFletcher- Reeves法を最適制御の計算に応用したアルゴリズムを次に 示す. (ここで用いられる評価関数の勾配は付録Aを参照)

• Fletcher-Reeves法による最適制御 計算アルゴリズム

(1)任意の時間関数uo(t)，t E [to， tJ]を与える.

(2)式(2.15)， (2.16)を用いて ， æo(t) = æ(t : uo)を得る. uo{t)，xo(t)を用いて， 式

(2.20)， (2.21)をt=りからt= toまで逆時間に解き， Po(t) = p(t : uo)を得る. これ より勾配Juoを計算 し， ⁸⁰ = -Juoとする. ⁱ = 0 と置く.

(3) i = 0のときは(4)へ， z三1のとき， ßi = (JUi， JUi)/(Jui-l， JUiー1)を求めて，

8i = -JUi + ßi8i-lと置く.

(4) J(Ui+α8i)を最小にするαを求め， これをαzとし， Ui+l二Ui+αiSiと置く.

(5)町+l(t)=æ(t:叫+l)，Pi+l(t)= p(t : Ui+l)を(2)におけると同じように計算 して ， 勾配JUi+1を計算 する.

(6) JUi+1が0に十分近ければ停止. そうでなければi=i+1として ， ステップ(3)

第2章 THERMAL-TO-THERMAL FLIGHT 2.6.

2.6

Thermal-to圃thermal flightの最適制御問題

微分方程式

V'二一[D(V， CL)/m

+ 9

sinγ]/V cos

'1 '1'

= [L(V， CL)/(mcosì ) - g]/V2 y'二tanγ

(2.24) (2.25) (2.26)

を満たし， 初期条件

V=V(O)，ì=γ(0)， y二y(0) _(2.27) を満たす次の

状態ベクトル:x = [V，γ ，Yf 制御変数:u = CL

(2.28) (2.29)

のうちから， 次の評価関数

X ，α

u z x nU

Z14

X fitl 。

+

，v，d x z x nuu 一一Td

(2.30₎

を最小にする最適軌道æ， Uを求める.

ここで， 式(3.12)のそれぞれの項は

9

= (V- V(0))21x=x，

+

(γ一γ(0))2Ix=x，

ん=

[

^{竹一V目的}lよrV cos γ

l

(2.31) (2.32)

である.

乗数ベクトルpを導入して， ハミルトン関数を次式のように置く.

H ⁼ - fo(X， x， u)

+

pT (X)f(X，夙叫 (2.33)

第2章 THERMAL-TO-THERMAL FLIGHT

ここで それぞれの項は

である.

ん=

[斤-

_lケVCOSγ ^Vsinγ

l

pT = [P1，p2，P3]

f=[V'，，'，Y'f

オイラーの方程式より， 乗数ベクトルの微分方程式は 次式のようになる.

i>(X) = - f æ ^T (X， æ， U)p(X) + fox (X，夙u) ここで， それぞれの項とその要素は

位θV二 色θVL 位δv二

[

^{F F}

H |

f æ ^T =

I

F21 F22 F23

I

⁼

I等努許

F31 F32 F33

笠θY二 むδYL 笠θYL

九1=-

bぷs

^{γ(GDo+ KGl)一}

^会

^tan-y

1

F1? = 2g 12一戸-

F13 = 0

I _pSV 1/'1 <，/'1')' Slll， . 9 1

21 = 一|_L τ_LTn 一(GDO+ I{Gl)三平+ Jァ|

- COSゐγ v cOS'"γJ

F22 =

[ ^凸

2m _cos2， ^sm

l

F23= 4-

cos'"γ

F31 = 0 F32二O F33 = 0

p = [Pl， P2 ， P3

f

2.6.

(2.34) (2.35) (2.36)

(2.37)

(2.38)

(2.39)

(2.40) (2.41) (2.42)

(2.43)

(2.44) (2.45) (2.46) (2.47) (2.48)

第2章 THERMAL-TO-THERMAL FLIGHT

fox

=

である.

V2 COSγ

l

�旦ユ一_{V COS2γ}

一

^{一一」一一­}_{VT COS2γ}

。

横断性の条件より， 次式のようになる.

p(Xf) 8g

- 8x(Xf)

r

^{2(V -} V(O))lx=x，

= 一I

2(γ- ì(O)) Ix=x，

1 I

評価関数の勾配は， 次式のようになる.

ム ー a互

'"

-

8u

= -f�(X，忽，

^{u) p (X) + f Ou (X ， x， 'U)}

ここで， それぞれの項とその要素は

である.

f�二[Fll， F12， F日

= イ[ ^{� 笠竺豆と，}

8CL ' 8CL ' 8CL

^8ι2乙L

^δ8Y' - mcosγ

Fl ? = �

1 2 - 2m cosγ ^-

F13 ⁼ 0

fou

=

⁰

p = [P1，P2，P3 f

以上の問題に対して共役勾配法を用いて計算を行った.

2.6.

(2.49)

(2.50)

(2.51)

(2.52) (2.53)

(2.54) (2.55) (2.56) (2.57) (2.58) (2.59) (2.60)

第2章 THERMAL-TO-THERMAL FLIGHT 2.7.

2.7

計算結果と考察

計算に使用した滑空機は， 日本飛行機株式会社のNIPPI PILATUS B4である. (付録 B参照) 計算例として ， 図2.3， 図2.4， 図2.5 に， XJニ500， 1000， 2000(m)， 1ケ=2， 4， 6(m/s) の場 合を示している. (線種については表2.1を参照)これは高度損失， 揚力係数， 飛行速度，

径路角のグラフで， それぞれのグラフの横軸は飛行距離Xをthermal聞の距離X.rで害IJって 無次元化して 表して いる.

図2.4の実線を例にとると， まず， 高度損失のグラフを見ると， 最初に定めた初期条件 の状態から出発した後， すぐに機体を急降下させて いるのがわかる. そのとき， 径路角 が下向きに大きくとられ， 飛行速度も増加して いる. その後， 揚力係数? 径路角， 飛行 速度ほぽ一定の飛行を行い， 次のthermalに近づくと 飛行速度と径路角の終端条件を 初期条件と等しくするために， 揚力係数を増加させて 急上昇を行って 増加した飛行速 度を減少させ， 終端条件と初期条件が合うようにして いる. そして， 次のもherrnalに到着

した後は， 時という速度で， もとの高度に到達するまで垂直方向に上昇を行って いる.

どの図を見て も， ì令の値が大きいほど飛行径路を大きく降下させているのがわかる.

1ケの値が大きいほど径路角をより 大きく下向きにとって ， 機体を降下させ， 飛行速度を 大きくとって いるのがわかる. つまり， VTの値が大きいほど， 次のもhermal内での高度の 回復が早くなるので， 多少高度損失が大きくなって も， 滑空機を大きく降下させて 飛行 速度を大きくした方が， total flight timeを短くすることができる. 表2.2には， それぞれ のグラフのtotal flight tirneを示して いる. これを見ると， \今の値が大きいほどtotal flight timeが短くなって いるのがわかる.

揚力係数のグラフを見ると， 滑空機はthermal内において は最小沈下速度で飛行する ということから， その時の揚力係数の値CL(O) ⁼ 0.99を初期値として 与えたにもかかわ らず， CL(O) ⁼ 0.3 からグラフが始まって いる. これは， CLの値を0.99から0.3へ急激に減少 させて 機体を急降下させて いるために， 揚力係数が2つのtherma.l内で最小沈下速度の 時の値に一致しなくなった.

そこで， 次節では揚力係数の値が0.99から始まり， その値の変化が急激な変化をせず，

2つのthermal内で最小沈下速度の時の揚力係数の値になるように揚力係数の変化率に 拘束条件を加えることにした.

第2章 THERMAL-TO-THERMAL FLIGHT 2.7.

表2.1線種

i

^VT(m/s)

I

^{線 種}

n4一Aヨ一円。

表2.2Total flight time

Fig. No. lケ(m/s) XJ(m) total flight time (sec)

I

2 24.25

2.3 4 500 20.30

6 18.83

2 46.07

2.4 4 1000 37.15

6 33.47

2 90.23

2.5 4 2000 70.27

6 62.21

L一一一

2.7

l. 5 O. 0

-内HUnu l tu --T

2. 5 2. 0 5. 5 S. 0 L 5 L 0 3. 5 3. 0

{凶\E)

〉

10. 0 _{_ 1} 寧10 •

THERMAL-TO-THERMAL FLIGHT

8. 0 6. 0

4. 0 2. 0

第2章

� ^-0.5

1/) ω O

�

-l.

0

11.1 勺.... P

ミーl.

5

- 2. 0 O. 0 事lO _t

O. 0

10. 0 車10 .

10. 0 _{_ 1} 事10 . 8.0

8. 0 6. 0

6. 0

X/Xf

グ/Jh / /・/γ 〆//ノ/;

， /ノd人 //

、、‘ ~、\

4.0

2.0 。

図2.3 Xf=500(m)，昨=2，4，6(m/s).

2. 0

O. 6

O. 2

ト‘

-0. 2 -0. ( O. (

0.0

(官同h)

l O. 0 _{_ 1} 事10 . 8. 0

6. 0 X/Xf 4， 0 2. 0

l. 2

O. 0 1. 0

O. 4 O. 2

-0. 2 O. 6 O. 8

ド斗o

2.7.

I nHU --E‘ AHU --v hhu THERMAL-TO-THERMAL FLIGHT

第2章

事lO _J O. 0

5. 5

L 5 3. ⁵

2. 5 (. 0

3. ⁰ 5. 0 (凶\E}

〉

) ーO.E ₅

凶M

J ーしoO

，。(lJ ョ

ミーし5

2. 0 l. ⁵

O. 0 L O. 0

事L 0 • 8. ⁰

6. 0 .. 0

2 .0 国 2 . 0

O. 0 10. 0

牢10 .

，_---.:，...'

/' "...

，/ /'

，/，/

" ，/ /

， • 人 /

いd1/ ・

\ノ

6.0 8.0 .. 0

2. 0

O. 6 O. 4 O. 2

O. 0 炉、

-0. 2 (勺同-H) l. 2

l. 0

O. 2 O. 8 O. 6

O. ₄

O. 0

』斗υ

10. 0 寧10 . 8. ⁰

6.0 X/Xf 4. 0 2. 0

図2.4 X _f =1000(m)， \ケ=2，4，6(m/s)

-0. 4

-0. 6 O. 0 10. 0

宇L 0 • 8. ⁰

6. ⁰ X/Xf 4. 0 2 . 0

-0. 2

2.7.

t

a--‘ 《Hu AHU --v EU

5， 5 THERMAL-TO-THERMAL FLIGHT

...ーー 司、-- 、、

，/ ーーーーーーーー ^\

/ --- ----_ーー- ^\

I>__' ...勺

j/ 、

i

5， 0 4. 5 L 0 3， 5 3， 0 2， 5

{凶\日}

〉

第2章

事10 ₁ O. 0

'''，，，・

一

/ //'一 J

~

~

~\ ~\ 』『\\ \\ \ \ \

\

\

\\ \、\\、\\・、

、

‘

、、.ぃ、、

、、仏、、.、，、

) ー O.日 5

ω ν、

� - l， O 0

'ロω ロ

ミー1.

⁵

2. 0 1. 5

0 ， 0 1 0. 0 _{_ 1}

事10 . 8，0

6. 0 4. 0

2. 0 - 2， 0

0， 0 1 0， 0 _.

事lO 1 .

10， 0 事lO • 8.0

6. 0 8.0 6. 0

X/Xf 4. 0

4. 0 2.0

図2.5 Xj=2000(m)，時=2，4，6(m/s)

2. 0

O. 6 O. 4

O. 0

-0. ⁶ O. 0 (句mNH)

、』

- 0 . 2 O. 2

-0. (

1 0 . 0 事l O • 8. 0

6. 0 X/Xf 4. 0 2， 0

1. 2 l. 0

O. 0 O. 6 O. 4 O. 2

- 0. 2 O. 8

...:l u

第2章 THERMAL-TO-THERMAL FLIGHT 2.8

2.8 Lift Coefficient Rate Control

前節までの計算結果を見ると， 2つのthermal内において揚力係数の値は， 初めに 設定した 最小沈下速度の時の値になってい ない. そ こで ， 揚力係数の値が急激な変化を 起こさず，もhermal内では最小沈下速度の時の揚力係数の値になるようにしたい.

したがって， 揚力係数の変化率が ， ある値の範囲内におさまるように次のような不等 式拘束条件を加えた.

m \11111/

く一 /Illt-\ ぬ一白

色白 、‘，』，， 'lよ 円ハU 円ノ白 ，，aEE、、

2.9

運動方程式の再定式化

運動方程式の再定式化を行う際に， 不等式拘束条件を等式拘束条件に変換し ， 今 まで述べてきた不等式拘束条件無しの最適制御問題の定式で解けるようにする.

そのため， 揚力係数に対して， 新しくダミー変数Udを導入し不等式拘束条件を等式 拘束条件に変換する. そうすることによって今まで制御変数であった揚力係数は状態変 数になり， ダミー変数が制御変数となって， 状態変数がlつ増えた不等式拘束条件なし の最適制御問題となる.

以上のことから運動方程式の再定式化を行うと

dV _dt

一一

-D(V， CL)/m - 9 sin γ (2.62) dγ [L(V， CL)!m

-

^{9 cosγ]/V (2.63)}

一

dt

dY V sin γ (2.64)

dt

dX _dt

一一

V cosγ (2.65)

dCL ^dt

^一 (す)山

^{mαz (2.66)}

のようになる.

先ほどと同様， 終端時間固定の最適制御問題の定式で解ける ように， 時間tを飛行距 離Xに置き換える. そのため式(2.65)で式(2.62)， (2.63)， (2.64)， (2.66)をそれぞれ割ると

，

運動方程式は

V'=一dX dV = 一[D(V， CL)/m + 9 sinγ]/V cosγ

，'=

�

^{ヱ ニ [L(V，}CL)!(m cosγ)

_

^g)/v2

dX

(2.67) (2.68)

第2章 THERMAL-TO-THERMAL FLIGHT 2.10.

y'

=

^{一 dY} _dX ^一

/'1' dCL

一 ^{一 ­}

VL

-

^dX

^{Ju 、、EEJ'1ょ}

u

nL n pD z a qL m

\Ell-/

γ 仇一ぽ n

，α一G

U

/{i\

(2.69) (2.70) のようになる.

2.10 最適制御問題の再定式化

微分方程 式

V' = -[D(V， CL)/m ⁺ g sinγ]/V cosγ イ=[L(V，CL)/(mcosγ)

_

^g]/V2

y':二もanì

C� =

( ^舎 にz

^(2山一1)

(2.71) (2.72) (2.73) (2.74) を満たし， 初期条件

V=V(O)，ì=γ(0)， Y = Y(O)， CL二CL(O) (2.75) を満たす次の

状態ベクトル: x = [V，γ，Y，CLY 制御変数: U = Ud

(2.76) (2.77)

のうちから， 次の評価関数

rx，

J = g(X， x(Xf)) ⁺

I

fo(X，夙u)dX

を最小にする最適軌道x， 五を求める.

ここで， 式(2.78)のそれぞれの項は

(2.78)

g = (V- V(0))21x=x， ⁺ (γ一γ(O))2Ix=x， ⁺ (CL - CL(0))2Ix=x， ^(2.79)

第2章 THERMAL-TO-THERMAL FLIGHT 2.10.

である.

式(2.83) のような乗数ベクトルを導入して， ハミルトン関数を次式のように置く.

H ₌ -fo(X， ^{x， u}) + pT (^X)j(X^{， X， u}) 〆lI、 nJH nRu 、、11'11ム ここで， それぞれの項は

o = [V�卜吟ト 一 Vsin Lよ今rV cos γ

pT = [Pl，P2，P3，P4]

f= [V'，γ/3Y/?cilT

(2.82) (2.83) (2.84)

である.

オイラーの方程式より， 乗数ベクトルの微分方程式は 次式のようになる.

þ(X)=-fæT(X，夙u)p(X)+ foæ (X， æ， ^u) (2.85 ) ここで， それぞれの項とその要素は

θθVV '

むθV: δδYV I 包θv^L 九1 F12 F13ん1 �θV _FY むθγ: δθYγ '

竺θγ^L 九九2 F23 九

|

fæT =

I

F31 F32 F33 F34 I I白θθVY' ， 色δY; θθYY ' 竺θY^L 九1 F42 F43 F44 J

I

δδCYL I 笠δC LL θθCVL ' 2θC2L L

F11=-lps山+KCI) -J2tan-y _2mcosì

1

F12 = V3 2g

F13 = ⁰ F14 = ⁰

F21=- (CDo+KCI)osっzγ+V cosっz γ

(2.86)

(2.87)

(2.88) (2.89) (2.90) (2.91)

第2章 THERMAL-TO-THERMAL FLIGHT 2.10.

F22二

[ヂ

^{cos γ si}

^t

^γ

l

F23= -_cosゐγ

L

^-

F24 = 0

(2.92)

F32 = ⁰

(2.93) (2.94) (2.95)

(2.96) (2.97) (2.98) (2.99) ん1 = ⁰

F33

二

^O

L一C

/一C一喝 一一vhJ一O-d11 FL-nb 一~l F一md一∞

ρ一一η

O 一 一2

0 一一一一一一一一 4 1 2 3 九 九

九 九

ん4 = 0

(2.100) (2.101) (2.102) (2.103)

p = [Pl， P2 ， P3， P4

Y

1

V2 cosγ

51旦ユ一一一一」一一 V C052 -y VT C052γ

foæ = ^(2.104)

。

。

である .

横断性の条件より， 次式のようになる.

8g (2.105)

p(Xf) = 8x(Xf)

2(V - V(O))lx=x，

2(γ-γ(O))lx=x，

(2.106)

。

第2章 THERMAL-TO-THERMAL FLIGHT

評価関数の勾配は， 次式のようになる.

ム一一θH... - åu

= -

f� (X， x， u)p(X) +ん(X，忽， u)

ここで， それぞれの項とその要素は

f�

= [Fll， F12， F13， F14]

=

[ôV'位。y'位l

ÔUd ^I ÔUd ^I ÔUd ^I ÔUd Fl1 ⁼ 0

F12 ⁼ 0 F13 = 0

(ぬ)

F14 ⁼

�

_{�X (、4}sin同cos Ud) max

101.1. ⁼ 0

P = [Pl ，P2，P3，P4

f

である.

以上の問題に対して共役勾配法を用いて計算を行った.

2.11

Lift Coefficient Rateの決定

2.11

(2.107) (2.108)

(2.109) (2.110) (2.111) (2.112) (2.113) (2.114) (2.115) (2.116)

(dCL/dX)max は， 飛行距離Xに対する揚力係数CLの変化率を表す. 計算を行うl緊に，

まず， この(dCL/dX)maxの値を与えなければならない.

計算結果のグラフの横軸は， 飛行距離Xを2つのthermal間の距離X.fで割って無次元化 して表しているので， ここでも， 計算の便宜上， XをXJで割って無次元化して表してお

無次元量をzと置くと，

であり， Xは

Z二一一X

X.f

X二xXf

(2.117)

(2.118)

第2章 THERMAL-TO-THERMAL FLIGHT

となる.

よって(dCL/dX )max は，

となる.

(会t.x

₌

(ポ�)t.x

=

( _会 t

_.

)

_f

2，12，

(2，119) (2，120)

このように表すと，例えば，(dCL/dx)max = 4の時，揚力係数CLのグラフの横軸X

=

_X/Xj

の変化量1に対して縦軸CLの変化量が最大で4 になることが容易に判断できる.

以降， (dCL/dx)maxの値や， 2つのthermal聞の距離Xfの値， thermal内での上昇速度l今 の値を変化させ， 先ほどと同じ初期値を用いて計算を行う.

2.12

拘束条件がある場合とない場合の比較

計算例として， まず，拘束条件のある場合とない場合の比較を行った. thermal聞の距 離XJ =500， 1000， 2000(m)， thermal内での上昇速度時= 4 (m/s) で， 拘束条件(dCL/dx)mα.τ=4 がある場合とない場合について比較したものを図2.6から図2，8に示した. (線種は表2，:3 を参照)これらのグラフには， 新たにダミー変数Udのグラフが加わっている.

CLは，初期条件CL(O)ニ0.99 から始まり，拘束条件の下限値に近い傾きで減少しながら，

機体は降下を始める. 例えば， 図2.8のCLのグラフを見ると， 横軸X

=

X/XfがO から0.1 まで変化する聞に縦軸のCLの値は， 約0.4 変化しており， (dCL/dx)mαzニ 4の拘束条件を満 たしていることがわかる. 拘束がある場合は，拘束がない場合に比べて径路が緩やか に降下しており， そのため飛行速度も緩やかに増加を始めている. CLの値がある程度 小さくなって， CLの傾きも小さくなってくると， 径路角を下向きにさらに大きくとり機 体を大きく降下させ， 飛行速度をさらに増加させている. 次の七hermalに近ずくと ， 拘束 条件があるために急上昇ができないので， 拘束がない場合に比べて早めに揚力係数が 増加し始め径路角も大きくなって上昇を開始していることがわかる. ある程度CLの値 が大きくなり CLの傾きが拘束の上限値に近ずいてくると径路角は小さくなり， 拘束が ない場合に比べて径路が緩やかに上昇している. 飛行速度， 径路角も次第に減少して ゆき， thermalに到達したときには，飛行速度，径路角， 揚力係数が初期条件にほぼ一致 し， CLについてはthermal内において， 最小沈下速度の時のCLの値になるという条件を 満たしている. また， 揚力係数の値が振動しているのは， 拘束条件をS2n関数で与えて

第2章 THERMAL-TO-THERMAL FLIGHT 2.12.

それぞれの図について， 表2.4にtotal flight timeを示している. これを見ると， 拘束条 件がある場合はない場合に比べて， どれもtotal flight timeは長くなっている. 例えば， 図 2.8 の速度のグラフを見ると， 拘束がなければ， 出発してすぐに急降下して速度を大き くとり， 次のthermalのすぐ手前で急上昇して大きな速度で飛行する区間が長くなるが，

拘束があると急降下， 急上昇ができないため， 加速が終わるとすぐに減速に転じてし まい， 大きな速度で飛行する区間が短くなるので， total flight timeは拘束がない場合に 比べて長くなることがわかる.

Fig. No.

2.6 2.7 2.8

表2.3線種

拘束条件 線 種 (dCL/dx)max = 4

拘束なし ^{t一一一一一一}

表2.4 Total flight time

lケ(m/_s) _{Xf(m) 拘束あり(sec) 拘束なし(日c)}

I

500 22.51 20.30

4 1000 41._{10 37.15}

2000 78.05 70.27

2.12.

THERMAL-TO-THERMAL FLIGHT

第2章

10. 0 _{_ .} 事10 . 8.0

t AHU AHU

---E・

phuw

令-v

1. 5

O. ⁰

5. 5 5. 0 4. 5 4. 0 3. 5 3. 0 2.5 2. 0

(ω\巨)

〉 宇lO _f

O. 0

� ^-0.5

ω vl O

，....l _ ^{l. 0}

cu

"0

-2

ミー1.

⁵

内HUAHM nHU 内ft・ •

4. 0 6.0 2. 0

1

O. ⁰ _{_.}

事10 . 8. 0

6.0 4. 0

2.0

O. ⁶

O. 4

& 6

O. 2

0.0

ト・・

- 0. ² - 0. ⁴ {古川h)

/一、 � � /_-、"

/ 、、ー ---- \

/ 一 一 一一一- _\

l. 2 1. 0 O. ⁸

O. ⁴ O. 6

O. 2 O. 0

炉斗 υ

10. 0 _{_.} 事lO •

8.0 6. 0

(. 0 L O. 0 事L0 _1 • 2.0

8. 0 6. 0

4.0 2. 0

-0. 2 -0. 4

O. 0

X/Xf 2. 0

l. ⁵

1. 0

O. 5

O. 0

-0. 5

(勺国同)

"0 D

l O. 0 .

事lO .

図2.6Xj=500(m)，ìケ=4(m/s) ，(dCL/ dx )max=4

8. 0 6. 0 X/Xf

4. 0 2.0

-l. ⁰ O. 0

七園七七七七 ^七・

2.12.

THERMAL-TO-THERMAL FLIGHT 第2章

‘、、‘‘、、、、、、、

\ \

\ \ / ~

/ /

，， /

，，， ，f

，，，

事L0 _I 6. 0

3. 5 3. 0

，‘1u

-aa---

5. 5 5. 0

L 0 (ω\巴)

〉

，f ，，

，，，

J〆/ / / / 寧10 ₁

O. 0

.ß ^-0. 5

- 1. 0

'ロ4J ロ

ミー1.

⁵

ωωoJ

2. 5

8. 0 6. 0

2. 0 4.0 - 2. 0

O. 0 L O. 0

事L0

• 8.0

6. 0 c 0

2. 0 10. 0

寧L0

•

0.6

O. 2

ト・.

-0. 2 O. 0 O. 4

{勺伺H}

1.

²

L. 0 O. 8

O. 0 O. 6 O. 4 O. 2 J ιJ

-0. 4 -0. 2

-0. 6 O. 0 L O. 0

牢10 • 8.0

6. 0 4. 0

2.0 10. 0

本L0

• 8.0

6.0

X/Xf 4. 0 2. 0

2. 0

1. 5

1. 0

O. 5 O. 0

(勺司h)

可ロ

-0. 5

10. -0 宇L0 •

図2.7 XJ=lOOO(m)，�ケ=4(m/s) ，(dCL/ dx )max=4

8.0 6. 0

X/Xf L 0 2. 0

ー1.0 O. 0

2.12.

THERMAL-TO-THERMAL FLIGHT

第2章

寧10 ₂ O. 0

寧10 I

6. 0 5. 5 5. 0

3. 5 3. 0 4. 5 L 0 (凶\巨)

、、、、、 \ \ \ \ \ \ ‘\ 、 〉

_g ^-0. 5

ω ul O

� -l， 0

勺CLI

.... ロ

ミ-l，

5

2. 5 2. ₀

10. 0 寧10 _• 8. 0

6. 0 L 0

-2. 0 2.0

O. 0 10. 0

本10 .

2. 0

6. 0 8.0

6.0 4. 0

4. 0 2.0

O. 6

O. 2 O. 0

-0. ₆ 0.0 {匂悶旬}

ト、

-0. 2 -0. ( O. •

10. 0

* 10 . 8. 0

6. 0 .. 0

2. 0 1. 2

1. 0 O. 8 O. ₆

- O. 2 O. _• O. 2 O. 0

� u

10. 0 ‘

牢10 . 8. 0

X/Xf 2. 0

1. 5 1. 0 O. 5 O. 0

(勺市川H)

"c 2

1 O. 0 寧10 .

図2.8 Xj=2000(m) ，Vr=4(m/s) ，(dCL/dx)max=4

6. 0 8.0 X/Xf

•. 0 2. 0

-0. 5 -1.0

O. 0

第2章 THERMAL-TO-THERMAL FLIGHT 主辺ー

2.13

(dCL/dx)mαzの影響

XJ• �今の値を一定にし• (dCL/dx)maxの値のみを変化させたものを図2.9から図2.11に示 し， その影響を調べた. (線種は 表2.5を参照)

図2.10を例にとると• (dCL/dx)maxの値が大きいほどCLの傾きが大きくなるので，揚力

係数が大きく変化する. そのため， (dCL/dx)maxの値が大きいほど機体をすみやかに降 下させて速度を大きくとることができ， 次のもhermalにより近いところから上昇するこ とができる. よって， (dCL/dx)maxの値が大きいほど， より大きな速度を得ることができ，

その大きな速度で飛行できる区間も長くなるので， total ftight timeを短くすることがで きる. 表2.6を見ると， (dCL/dx)maxの値が大きいほどtotal ftight timeが短くなっているのが

わかる.

Fig. No.

2.9

2.10

2.11

表2.5線種

|

^(dω

2

4 ^{I一一一一一一} 8 Iーーーーーーーー

表2.6 Total flight time

時(m/s) ^XJ(_m) (dCL/ dX)max total flight time(sec)

2 26.68

2 500 4 26.04

8 25.77

2 51.27

2 1000 4 49.43

8 47.79

2 100.57

2 2000 4 94.97

8 92.66

THERMAL-TO-THERMAL FLIGHT 主旦ι

第2章

t AMU aHU

--

-- phu

5. 5 5. 0 (. 5 (. 0 3. 5 3. 0 2. 5 コ\巨)

〉

� ^-0.5

凶凶 O

，....l -1. 0

'0 乱b

;:l

ミー1.

⁵

事L0 _I O. 0

10. 0 寧10 . 8.0

4.0 6.0 2.0

O. 6 O. • O. 2

O. 0

ト‘

-0. 2 - 0. ( (唱団口}

10. 0 _{_ .} 事1 0 •

2. 0

8. 0

8.0 6. 0

6. 0 (. 0

L 0 2. 0

L 2

1.

⁰

O. 8 O. 6

O. 2 O. 0 -0. 2 O. 4 ーA 仁J

10. 0 寧10 . B.O

6.0 X/Xf 4. 0 1 0. 0 _{_ .} 2.0

亭lO •

2. 0

1. 5

1.

⁰

O. 5 O. 0 - 0 . 5

(司伺H)

て33

ー1. 0

O. 0 10. 0 .

寧lO •

図2.9 Xj=500(m)】竹=2(m/s)，(dCL/ dx )max=2，4，8.

B. 0 6. 0

x/xr .. 0 2. 0