電気通信大学大学院情報理工学研究科博士

(1)

入力状態安定性を保証する量子化入力による宇宙機のモデル予測制御

淺川岳也

電気通信大学大学院情報理工学研究科博士

(

工学

)

の学位申請論文

2015

年

3

月

(2)

入力状態安定性を保証する量子化入力による宇宙機のモデル予測制御

論文審査委員

主査木田隆教授委員田中一男教授委員新誠一教授委員明愛国教授委員樋口幸治准教授

(3)

著作権所有者淺川岳也

2015

(4)

Thisdotoralthesisstudiesamodelpreditiveontrol(MPC)oflinearsystemsdriven

by quantized ontrol inputs. We formulate the problem as an integer quadrati pro-

grammingproblem,andnewlyderiveonditionsfortheterminalostwhihisobtained

by solving linear matrix inequalities (LMI) so as to ensure the losed-loop stability.

Some numerialandexperimentalstudiesare performedtoverifythe proposedMPC.

InChapter1,werstgivestudyobjetivesandproblemdesriptions. Chapter2derive

aLMIonditionforahievingthelosed-loopstabilityoflineartimeinvariantsystems

(LTI) both for ontinuous ontrol inputs by reation wheels (RW) and on-off ontrol

inputsbyreationontrolsystems(RCS).ByusingthesolutionofLMIastheterminal

ostmatrixinoptimizationproblems,itisshownthatthelosed-loopsystemisasymp-

totiallystable(AS)forontinuousvalueinputsandtheonditionisageneralizationof

the well-knownonditionsderivedbyaRiatiequation. In theaseofon-offinputs,

it is also shown that the losed-loopsystem results in the Input-State-Stability (ISS)

ratherthanAS.Additionallywealsoshowthatitanbegeneralizedtothemulti-stage

quantizedinputase.Asanumerialstudy,weapplyMPCtospaeraftattitudeontrol

problemandomparetheontrolperformanefor ontrolinputsofontinuous,on-off

andtheirombinationandshowtheeffetivenessoftheombineduseofRCSandRW

frombothaspetsofagilityandauray. Wealsoshowtheeffetivenessofproposed

MPC underthestateonstraintbyapplyingtothespaeraftformationightproblem

withollisionavoidaneability.InChapter3,weextendtheresultintwoways. Firstis

theextensiontolinearparametervaryingspaeraft(LPV)andderivetheLMIterminal

ost onditions in order to ensure the losed-loopISS. Other is the robust stabilizing

MPC against dynami model errors. We then perform some experimental results for

LTI, LPVand robust MPCusing asingle-axisair stageand examinetheimplementa-

tionoftheon-linealgorithms.Chapter5isforonlusionsandperspetives. Thebasi

fatsaresummarizedintheAppendies.

(5)

概要

本研究は，宇宙機の大角度姿勢変更によく用いられるスラスタ・クラスタ（RCS）による最適制御を目的としている．RCSは，角運動量の交換装置であり連続値制御入力を発生するリアクションホイール（RW）に比べて発生できるトルクが大きいという長所をもつ反面，オンオフの2値となるトルクのみを生成する非線形アクチュエータであることが最適制御設計を困難としてきた．そこで，ここでは整数二次計画問題の解法によるモデル予測制御の適用を考え，入力状態安定性（ISS）を達成するアルゴリズムを構築する．そして，数値シミュレーションおよび実験でその有効性を示す．各章の構成は以下となっている．

まず1章で研究目的・問題設定の記述を行い，2章で線形時不変（LTI）の宇宙機を対象にして，RWによる連続値入力，RCSによるオンオフ入力の双方について閉ループ系を安定とする終端拘束行列が満たすべき条件を線形行列不等式（LMI）として導出する．このうち連続値入力の条件は従来の漸近安定となるためのリカッチ方程式条件を包含する一般化された条件であること，オンオフ入力の場合は漸近安定ではなくISS条件であることを示す．そして，RWとRCSの組み合わせ使用によって，宇宙機の高速かつ漸近安定な大角度姿勢変更が可能となることを数値例を用いて検証する．さらに得られた結果は多段の量子化入力に拡張できることを示す．また衝突回避機能をもつ軌道上フォーメーションのような宇宙システムへの適用結果を述べる．

3章では，これらの結果を宇宙機がパラメータ変動を持つ線形パラメータ変動

（LPV）系である場合に拡張して，やはり端点モデルに対する連立LMIから求めた共通解を終端拘束に用いることでISSを保証できる事を示す．更に動的なモデル誤差に対してロバスト安定なモデル予測制御則を提案する．これらの結果を用いることによって，宇宙機本体に対して回転する大型柔軟な太陽電池パネルなどをもつ現実的な宇宙機の量子化入力によるモデル予測制御が可能となる．大型人工

衛星ETS-VIIIの数値例を用いてその有効性を示す．

4章では，静空気圧で浮上させた1軸回転自由度を持つエアステージの実験装置を用いた上述のオンオフ入力によるLTI，LPVおよびロバストモデル予測制御の実験結果を報告する．実験の目的は，モデル予測制御のオンライン最適化計算を市販の計算機に実装可能であるかを検証することであり，いずれの場合にも妥当なサンプリング周期で実現できることを明らかにする．

5章は，結論と課題であり，宇宙システムのオンオフ制御に対する展望，宇宙開発以外の分野への応用可能性について述べる．最後に，宇宙機の運動方程式，入力状態安定性など基礎的な事実を付録に一括してまとめる．

(6)

第

1

章序論

1.1

研究背景

近年では，民間や大学で開発された小型衛星をピギーバック

1

で打ち上げるなど，宇宙機開発はより身近なものになっている．一方で大型宇宙機によるミッションはますます大規模なものになっている，例えば，既にJAXAが開発運用中の通信衛星ETS-VIII(Fig.1.1)は3[ton℄の大型衛星であり，また現在計画中の天文観測衛星ASTRO-H(Fig.1.2)は質量2.4[ton℄，望遠鏡進展後，全長14mにも及ぶ大型のX

線天文衛星となる予定である．これは日本の有する天文衛星では最大クラスである．またNASAが開発中のジェイムズ・ウェッブ宇宙望遠鏡(JWST，Fig.1.3)は，

質量6.2[ton℄になる見込みで，ハッブル宇宙宇宙望遠鏡が11[ton℄に対して小型化

したことになるが，そのミッション要求からJWSTはより遠い軌道に投入され，その基線距離を考えるとこれもまた大型の宇宙機と言える．このような天文衛星は観測のために指向方向を自在に変更することが求められ，指向方向への制御精度は観測の性能に直接影響し，目標方向の指向の速さも同時に必要である．そのためこれらの大型衛星の三軸姿勢制御の設計は，ミッションの達成に非常に大きく影響するものである．

こうした宇宙機の大型化や構造の複雑化に伴って，三軸姿勢制御にはますます大きい課題が与えられている．まず，大型化に伴う構造振動制御が重要である．制御性能に大きな影響を及ぼす低周波数帯域に，宇宙機の構造から生じる弾性振動の周波数が近づきスピルオーバー不安定と呼ばれる，制御器と構造振動の干渉による不安定さが発生することがよく知られている[1,2℄．これらを回避する制御系設計手法には，線形二次レギュレータ(LQR)に周波数遮断フィルターを用いて，弾性振動周波数が励起されにくくすることで安定性を確保する方法[3,4℄や，力学系の二階行列微分方程式の特性に着目して，制御器を設計するDiretVeloityDiret

Feedbak(DVDFB)を用いる方法[5,6℄，または周波数特性を考慮することができる

H ∞

制御による制御系設計[7℄などが検討されてきた．次に，上記の制御系設計では直接対応できないが大型宇宙機では容易に発生しうる入力飽和やアクチュエータの特性の考慮，すなわち制御入力量に対する制約も重要な要素である．

宇宙機のアクチュエータは主に内力アクチュエータと外力アクチュエータに区分される．内力アクチュエータは例えば，リアクションホイール(RW)やコント

(9)

ロールモーメンタムジャイロ(CMG)などの回転体を有するアクチュエータがそれにあたる．ここでリアクションホイールは，この回転体と宇宙機の角運動量交換によってトルクを生成するもので，姿勢安定化や姿勢制御を行うことができるが，

問題点として衛星の大きさに対して比較的小さいレンジのトルク(数十[mNm℄や

数[Nm℄)しか生成できないことや，回転体の有する角運動量が飽和した場合には

外力アクチュエータにより角運動量を宇宙機の外に放出するアンローディングが必要となることがあげられる．また外力アクチュエータにおいて主に用いられるスラスタは，高圧のガスやイオンを高速で噴射し推力を得るアクチュエータである．多くの場合スラスタはクラスター状に集められてリアクションコントロールシ

ステム(RCS)として衛星に搭載される．スラスタはアンローディングで使用され

ることや軌道変更にも用いられていることから分かるように，内力アクチュエータに比べて大きなトルク(数[Nm℄や数十[Nm℄クラス)を生成することが可能である．一方で，推薬弁の開閉が直接的な操作量となるためオンオフ信号によって所望の制御入力を実現しなくてはならない．そのためこのオンオフ特性から擬似的に連続的なトルクを得るために，PulseWidthModulation(PWM)を行って線形制御則と組み合わせて実装する方法が広く用いられている[1℄．しかしこの方法では，

制御性能と安定性を厳密に保証することは困難となるといった問題がある．

このように，近年求められている姿勢制御系の設計には様々な要素を複合的に考慮しなくてはならず，柔軟に設計仕様を組み込むことができる制御系設計が求められている．そこで本論文では，上記のさまざまな問題を近年注目されているモデル予測制御[8,9℄を用いて解決することを目的とする．モデル予測制御は，理論的には最適制御問題[10℄をもとにしていて，制御系設計がしやすいという利点がある．また近年では，線形システムのみに限らず非線形システムやハイブリットシステムを対象とする理論研究[11,12℄[13℄が進んでいる．それらの結果を受けて，衛星の姿勢制御[14, 15,16℄やフォーメーションフライト[17, 18℄などの宇宙システムを対象とした問題についても盛んに研究されている．

オンオフ入力を必要とするシステムに対するモデル予測制御は前述の[14,15,16℄

や，宇宙機以外の分野では自動車の電空クラッチアクチュエータを対象とした文献

[19℄，同様に自動車の暖房を対象とした研究[20℄がある．文献[14℄では，モデル予測制御の一つであるExpliitModelPreditiveControl[21℄を用いて制御系を構築し，剛体衛星の姿勢制御問題に適用している．この文献では，実際のRCSへの信号は不感帯を有するBang-Bangモデュレーションによる近似を行っており，PWM

を用いて実装する手法と同様に安定性や制御性能に劣化が生じる可能性がある．一

方，文献[15℄の手法はExpliitMPCを用いて入力変動分に対する制約を与えるこ

とで，

{−1, 0, 1}

の3値のオンオフ入力を直接の制御入力として，RCSを用いた剛体衛星の姿勢制御問題を取り扱っている．これらの方法は実装する上で，計算負荷が小さくできる点で優れている．しかしながら，安定性についてはまだ十分に議論がなされていない．そこで本研究では，安定性についての議論も含めて3値のオンオフ入力を直接の制御入力とするオンオフモデル予測制御を考える．

(10)

(11)

Fig1.3:ImageofJamesWebbSpaeTelesope(C)NASA

(12)

1.2

研究目的

本研究の目的を次のように設定する．

1. リアクションコントロールシステムを主たるアクチュエータとして高速姿勢マヌーバを実現する．

2. そのモデル予測制御系の構築法と安定性を明らかにする．

3. LPVモデルや柔軟構造物を有する宇宙機モデルについても安定性を保証した制御系を構築する．

4. 地上実験を用いて，オンライン計算や制御器の実装が有効であるかを確認する．

1.3

問題設定

本論文では，次の3つのモデルに対する制御系設計を通じて，宇宙機の姿勢制御問題に取り組む．まず第2章では剛体衛星の姿勢制御問題を考える，これは3

値の離散入力によるオンオフモデル予測制御をもっとも簡単なモデルに適用して，

その有効性や安定性について議論するためである．

次に，第3章で線形パラメータ変動システムとなる柔軟衛星の姿勢制御問題に拡張する．実際の宇宙機は軌道上でソーラーパネルを太陽方向に指向させることが多く，その結果としてモデルのダイナミクスの変動が生じる．更に柔軟構造物の振動を考慮して，モデル低次元化とそれに伴う残余モードに対するロバスト安定化問題を併せて検討する．また第4章では第3章で得られた結果をもとに地上実験を行い，提案手法の有効性とオンライン計算アルゴリズムの実装可能性を確認する．

1.3.1

剛体衛星

軌道上の剛体衛星の運動方程式として平衡点近傍で線形化した次式を考える．

J p ¨ = Lu

^(1.1)

ここで，

p ∈ R ³

はロール，ピッチ，ヨーを含む角度ベクトルであるものとし，

J

は慣性モーメント行列，

L

は制御入力

u ∈ R ³

の係数行列である．本研究のモデル予測制御は離散時間システムを扱うため，上記の式を零次ホールドを仮定して一定時間サンプリング間隔で離散化し，次式を得る．

x(k + 1) = F x(k) + Eu(k)

(13)

この時，状態量は

x ∈ R ⁶

であり，制御入力

u ∈ R ³

としてリアクションホイール

(RW)とリアクションコントロールシステム(RCS)の2種類のアクチュエータを想定する．したがって，制御入力

u(k)

は

u = u rcs + u rw

^(1.3)

となり，それぞれのアクチュエータの合力で表現され，それらの制約は，

u rcs : {−1, 0, 1} ∈ Z ³

(1.4)

u rw : u _{rw,min} ≤ u rw ≤ u _{rw,max} ∈ R ³

(1.5)

で与えられる．ここで，

u {rw,min}

と

u {rw,max}

はリアクションホイールの最大最小出力を表す．

1.3.2

パラメータ変動を有する柔軟衛星

宇宙機が持つ太陽電池パネルは，常に太陽方向を指向するように宇宙機の軌道周回に同期して回転する．その結果，システムは次の非拘束モードモデル(付録B.1)

で記述される線形パラメータ変動モデルとなる．

ξ ¨ + Λ(θ) ˙ ξ + Σ(θ)ξ = Γ(θ)u

^(1.6)

ここで，

ξ ∈ R ⁿ

は剛体および弾性振動のモード座標，

Λ(θ) ≥ 0

はモード減衰項であり，

Σ(θ) ≥ 0

はモード剛性項，そして

Γ(θ) ∈ R ^n×m

は入力係数行列である．

上記のダイナミクスを離散型LPVシステムに変換した次のシステムについて考える．

x ^∗ (k + 1) = F (θ k )x ^∗ (k) + E(θ k )u(k)

^(1.7) ここで，

x ^∗ ∈ R ²ⁿ

は状態変数，

u ∈ R ^m

は制御入力，

θ k ∈ θ

は

k

ステップにおける既知の変動パラメータである．そして，

(F (θ _k ), E(θ _k ))

はすべての

θ _k

において可制御であり，

F (θ k )

^,

E(θ k )

はいずれも有界で，端点数

p

のポリトープ表現を用いて，

つぎの表記が可能であると仮定する．

F (θ k ) =

p

X

i=1

ρ i (θ k )F i , E(θ k ) =

p

X

i=1

ρ i (θ k )E i

^(1.8)

ここで，スカラ関数

ρ i (θ k )

は

p

X

i=1

ρ i (θ k ) = 1, ρ i (θ k ) ≥ 0 ∀i

^(1.9)

を満たす．またこの問題ではRCSのみの使用を考えることとして，制御入力

u(k)

の入力制約集合

Υ

は

u ∈ Υ = {−1, 0, 1} ^m ⊂ Z ^m

(1.10)

であるとする．

(14)

1.3.3

柔軟衛星

柔軟構造物を有する宇宙機は，多くの場合数値計算をもとに非拘束モードモデルが得られる．

ξ ¨ + Λ ˙ ξ + Σξ = Γu

^(1.11)

ここで，

ξ ∈ R ⁿ ^c

は剛体および弾性振動のモード座標．

Σ ≥ 0

はモード剛性項，

Λ ≥ 0

はモード減衰項であり，そして

Γ ∈ R ^n×m

通常，モード次数

n c

は非常に大きいので制御系設計には，低次元化モデルを用いる．これを離散型システムに変換した状態方程式

x(k + 1) = F x(k) + Eu(k) + w

y(k) = Hx(k)

^(1.12)

について考えるものとする．この時，状態量

x ∈ R ²ⁿ ^c

であり(

n _c < n

⁾，

w ∈ R ⁿ ^c

は低次元化した際のモデル誤差を主とする外乱である．RCSにより生じる制御入力

u(k) ∈ R ^m

は，次の入力制約集合

Υ

u ∈ Υ = {−1, 0, 1} ^m ⊂ Z ^m

(1.13)

を持つものとする．

(15)

1.4

の最大特異値，最小特異値．

(16)

第

2

章モデル予測制御

2.1

はじめに

本章では，離散LTIシステムを対象とした3値量子化入力(オンオフ)によって制御する新しいモデル予測制御手法を提案する．概要では，これまでの連続値入力によるMPCにも共通する最適化問題の定式化を行う．その後，提案手法が入力状

態安定(ISS)を保証することを示す．数値シミュレーションでは，連続入力(RW)

とオンオフ入力(RCS)を併用した場合の性能について，また，

n

値量子化入力の最適化問題に拡張できることを示す．最後に，状態制約が容易に扱えることを利用したフォーメーションフライト問題への適用結果を示す．

2.2

概要

以下の線形離散時間時不変システムについて考えるものとする．

x(k + 1) = f(x(k), u(k)) = F x(k) + Eu(k)

^(2.1)

ここで

x(k) ∈ R ⁿ

と

u(k) ∈ R ^m

はそれぞれ時刻

k

における状態量ベクトル，制御入力ベクトルであり，

(F, E)

は可制御であるとする．モデル予測制御の概要は一般に以下のようになる．

まず，現在の時刻

k

における状態量

x(k)

をもとに

x(k + 1)

から

x(k + N )

への

N

ステップ先までの予測状態軌道

X (k) = [x(k + 1) ^T , . . . , x(k + N ) ^T ] ^T

を構成する．ここで，

N

は予測ステップ数である．同様に対応する制御入力列に関しては時刻

k

から

k + N − 1

まで考え，U

(k) = [u(k) ^T , . . . , u(k + N − 1) ^T ] ^T

と記述する.

そして，X

(k)

と

U (k)

に対してある評価関数

V ( X (k), U (k))

の各サンプリング時間ごとの最適化を行い，得られた入力列

U ô (k) = [u ô (k) ^T , . . . , u ô (k + N − 1) ^T ] ^T

のうち，ステップ

k

時点の入力列をシステムに印可する(すなわち，

u(k) ^T = u ^o (k) ^T

とする)．これをステップ毎に繰り返す．

予測状態軌道

X (k)

は，システムが離散線形である場合には，次の方法で簡単に求めることができる．まず，式(2.1)のシステムに対して，現在の状態量

x(k)

から予測状態軌道

X (k)

は，

X (k) = F x(k) + EU (k)

(17)

となる．ここで，それぞれの係数行列

F，E

は，

F =





 F ¹ F ²

.

F ^N





 , E =







F ⁰ E 0 . . . 0 F ¹ E F ⁰ E 0

^.^.^.

.

. .

.

0 F ^N ⁻¹ E . . . F ¹ E F ⁰ E







(2.3)

である．次に評価関数は二次形式評価関数として

X (k)

と

U (k)

から，

V _N (x(k), U (k)) = X (k) ^T Q X (k) + U (k) ^T R U (k)

=

k+N−1

X

i=k+1

x(i) ^T Qx(i) + u(i − 1) ^T Ru(i − 1) +x(k + N ) ^T P x(k + N )

(2.4)

= ( F x(k) + EU (k)) ^T Q( F x(k) + EU (k)) + U (k) ^T R U (k)

= ǫ ^T Qǫ + 2ǫ ^T Q EU (k) + U (k) ^T ( E ^T Q E + R) U (k)

を得る．ここで，予測状態軌道のうち定数となる項

F x(k)

は

ǫ

とおいた．以上より，

この評価関数の制御入力

U

による最適化は

2ǫ ^T Q EU (k) + U (k) ^T ( E ^T Q E + R) U (k)

についてのみ考えればよいことがわかる．

よって，本論文で用いる最適化問題は，

min U { U (k) ^T Φ U (k) + φ ^T U (k)}

s.t.

x : x(k) ∈ Ξ, ∀k u : u(k) ∈ Υ, ∀k

U : [u(k) ^T , . . . , u(k + N − 1) ^T ] ^T

(2.5)

(18)

である．ただし，

Φ = E ^T Q E + R φ = 2( F x(k)) ^T Q E

とした．この最適化問題の評価関数が凸であるための条件は

Φ = Φ ^T ≥ 0

のときである．

2.3

モデル予測制御の安定性

ここで，オンオフ入力で制御されるLTIシステム(2.1)の入力状態安定(ISS)条件を導出する．入力状態安定の概念は，E.D.Sontagによって導入された安定条件

[23, 24,25℄で，リヤプノフ安定性の枠組みをもとに，外乱入力による非線形性を許す安定条件となっている．そのため非線形制御問題の安定解析に対する有用なツールとしてよく用いられている．更に，JiangとWangによって離散時間システムに拡張された[26℄ことによって，非線形モデル予測制御やロバストモデル予測制御の安定性解析にも利用されている[27,28℄．

モデル予測制御を用いる際の準備として，まず終端状態量

x(k + N )

と終端制御入力

u(k + N − 1)

をそれぞれ

x N

，

u N

として，以下の仮定をおく[29,30℄．仮定1

[A1℄

0 ∈ X f ⊂ Ξ

かつ

X f

は閉集合

[A2℄

u N (x N ) ∈ Υ, ∀x N ∈ X f

[A3℄

f(x N , u N ) ∈ X f , ∀x N ∈ X f

ここで，

X _f

は終端制約集合を，

f (x _N , u _N )

はシステムのダイナミクス(2.1)を表す．仮定1は，終端状態量

x N

が

X f

の中にあるとき，

x N

は状態制約集合

Ξ

を満足し，それに対応する終端制御入力

u N (x N )

は入力状態集合

Υ

を満足すること，そして，

X _f

は不変集合であることを意味する．

この仮定のもとで，つぎの補題が得られる．

補題1 モデル予測制御の評価関数(2.4)をISSリアプノフ関数の候補とする．仮定1のもとで，クラス

K _∞

関数

α ₃

およびクラス

K

関数

σ

が存在して，次の不等式

F ^∗ (x N ) + L(x N , u N ) ≤ −α 3 (|x N |) + σ(|µ|), ∀x N ∈ X f

^(2.6)

を満足する時，閉ループシステムの原点は入力状態安定となる．ここで，終端評価の差分

F ^∗ (x N )

と状態量と入力の増分

L(x N , u N )

は，それぞれ

F ^∗ (x N ) = |x N+1 | ² _P − |x N | ² _P

L(x N , u N ) = |x N | ² _Q + |u N | ² _R

^(2.7)

(19)

証明1 最適状態軌道

X ^o (k)

とその最適入力列

U ^o (k)

が

X ô (k) = {x ô (k + 1) ^T , x ô (k + 2) ^T , . . . , x ô (k + N ) ^T } ^T U ô (k) = {u ô (k) ^T , u ô (k + 1) ^T , . . . , u ô (k + N − 1) ^T } ^T

と与えられるとすると，1ステップ先の軌道は

X (k +1) = {x ô (k +2) ^T , . . . , x ô (k +N ) ^T , x(k +N +1) ^T } ^T U (k +1) = {u ô (k +1) ^T , . . . , u ô (k +N −1) ^T , u(k+N ) ^T } ^T

となり，その評価関数は現在の評価関数の値

V ( X ^o (k), U ^o (k))

を用いて，

V ( X (k +1), U (k +1)) = V ( X ^o (k), U ^o (k))

− |x ô (k +1)| ² _Q −|u ô (k)| ² _R −|x ô (k +N)| ² _P

^(2.8)

+ |x ^o (k +N )| ² _Q +|u ^o (k +N −1)| ² _R +|x(k +N +1)| ² _P

と記述することができる．ここで，リアプノフの安定定理と同様に，1ステップにおける評価関数の差分を

∆V

として

∆V = V ( X (k + 1), U (k + 1)) − V ( X ^o (k), U ^o (k))

^(2.9)

とおく．式(2.7)を用いて式(2.8)は

∆V + L(x ô (k +1),u ô (k)) = F ^∗ (x ô (k +N))+ L(x ô (k+N ),u ô (k +N −1))

^(2.10)

と書ける．式(2.10)の状態変数

x ^o (k + N )

と入力

u ^o (k + N − 1)

を

x N

と

u N

とおくと，式(2.10)中の

L(x ^o (k + 1), u ^o (k))

が非負であることから

∆V ≤ F ^∗ (x N ) + L(x N , u N ), ∀x N ∈ X f

^(2.11)

が得られる．したがって，もしクラス

K _∞

関数

α ₃

およびクラス

K

関数

σ

が存在して，式(2.6)を満足するならば，

V N (x(k), u(k))

はISS-リアプノフ関数となり，ISS

に関する定理3(付録C)より閉ループシステムの原点は入力状態安定となる．

この補題を用いて，次の定理を得る．

定理1 状態量と入力の重み行列をそれぞれ

Q > 0

，

R ≥ 0

とし，終端制御入力を定数ゲイン入力

Kx N

と補償入力

ν

から次式のように構成されているとする．

u N (x N ) = Kx N + ν

^(2.12)

仮定1のもとで，あるスカラー値

γ > 0

が存在して，終端評価行列

P

が条件

S ^T P S − P + Q + K ^T RK ≤ −γI n , ∀x N ∈ X f , u N (x N ) ∈ Υ

^(2.13)

(20)

x ^N

-1 1 0 u

X

x

sat

N

K v

Fig2.1:On-offinput

x ^N

0 u

X

x

sat

N

K

v

1 -1

Fig2.2:

n

quantizationinput

を満たすとき，最適化問題(2.5)を解くオンオフモデル予測制御からなる閉ループシステムはISSとなる．

ここで

X f ⊆ {x| |K i x| ≤ 1 ∀i}

であり，

S

は，

S = F + EK, |λ max (S)| < 1

^(2.14)

である．

証明2 まずモデル予測制御の評価関数(2.4)がISSリアプノフ関数の候補であると仮定し，補題1を用いて安定性を示す．式(2.6)の左辺は，

F ^∗ (x _N ) + L(x _N , u _N ) = x ^T _N ₊₁ P x _N ₊₁ −x ^T _N P x _N + x ^T _N Qx _N +u ^T _N Ru _N

． (2.15)

ここで終端状態量

x N

が終端状態制約集合

X f

を満たし，終端制御入力

u N

が入力制約集合

Υ

を満たすとき，補償入力

ν

のすべての要素

ν i

を式(2.12)より以下の範囲で選ぶことができる(Fig.2.1)．

ν i ∈ {ν i | |ν i | ≤ 1}

^(2.16)

したがって，式(2.15)は，式(3.4)，(2.12)，(2.14)より，

(2.15) = x ^T _N (S ^T P S − P + Q + K ^T RK )x N + F

^(2.17)

を満たす．ただし

F

は，

F = 2ν ^T (E ^T P S + RK)x N + ν ^T (E ^T P E +R)ν

^(2.18)

(21)

である．終端評価行列

P

がLMI(2.13)を満たすとき，以下の不等式が成り立つ．

(2.17) ≤ −γx ^T _N x _N + F

≤ −γ|x N | ² +(kEk ² kP k+kRk)|ν| ² + 2(kEkkP kkSk + kRkkKk)|ν||x N |

^(2.19) ここで式(2.19)の右辺に非負の

(|x _N |/z −z(kEkkP kkSk + kRkkK k)|ν|) ² ≥ 0

を加えると，

(2.19) ≤ −γ|x _N | ² + |x _N | ² /z ²

+ z ² (kEkkP kkSk + kRkkK k) ² |ν| ² + (kEk ² kP k + kRk)|ν| ²

^(2.20) が得られる．したがって式(2.15)は

F ^∗ (x N ) + L(x N , u N )≤ −α 3 (|x N |) + σ(|ν|)

^(2.21) を満たす．ここで，

α 3 ( |x N | )=(γ − 1/z ² )|x N | ²

σ( |ν| )=( z ² (kEkkP kkSk + kRkkKk) ² +kEk ² kP k + kRk )|ν| ²

である．明らかに

γ > 1/z ² > 0

を満たす任意のスカラ変数

z

に対して，

α 3 (|x N |)

はクラス

K _∞

関数であり，

σ(|ν|)

はクラス

K

関数である．したがって評価関数

V _N

はISSリアプノフ関数であり，補題1より閉ループシステムはISSである．また，

終端制御則に

u N (x N ) = Kx N + ν

を用いた評価関数

V N (x(k), U (k))

と最適オンオフモデル予測制御の評価関数

V N (x(k), U ^o (k))

を比較すると，最適性より

V N (x(k), U (k)) ≥ V N (x(k), U ^o (k)), ∀k ∈ Z ₊

(2.22)

となることは明らかである．したがって最適化問題 (2.5)を解く最適オンオフモデ

ル予測制御もまたISSを保証する．

これより連続値入力の場合の次の安定条件が得られる．

系1 連続値入力が実現できるときは，

ν = 0

とすることができ，それに伴って

式(2.6)中の

σ = 0

u N (x N ) = Kx N

とする，ここで

K = (E ^T P E + R) ⁻¹ E ^T P F

．

式(2.13)へ代入すると，次式を得る．

x ^T _N (F ^T P F −P +Q−F ^T P E(E ^T P E+R) ⁻¹ E ^T P F )x N = 0, ∀x N ∈ X ¯ sat ⊆ X f .

^(2.24) ここで，

X ¯ sat = {x ∈ R ⁿ |u min ≤ u N (x) ≤ u max }

は入力制約に違反しない状態制約

集合である．

このDARE(2.23)の解

P

は証明から明らかなように，LMI(2.13)を満たす解に包含される．したがって，LMI条件(2.13)が，MPCの終端評価の設計においてより柔軟性があるといえる．

系3 式(2.12)について，

n

値量子化入力の時(Fig.2.2)，3値量子化(オンオフ)入力と比べて

ν

は明らかに同じか，より小さくすることができる．

2.4

モデル予測制御のアルゴリズム

提案するモデル予測制御のアルゴリズムをまとめて示す．

オフライン設計:

R 1.

評価関数中の任意の設計パラメータ行列

Q, R > 0

を設定する．

R 2. F + EK

のスペクトル半径を1以下とするゲイン

K

を設定する

1

．

R 3.

^LMI^(2.13)を解くことで終端評価行列

P

を得る．

R 4.

予測ステップ数に対応する拡大設計パラメータ行列

Q, R

を構成する．

オンライン最適化:

S 1.

予測ステップ

N

先までの予測状態軌道

X

(2.3)を計算する．

S 2.

現在の状態量

x(k) ˆ

を取得し，最適化問題(2.5)を解く．

S 3.

最適化し得られた入力列

U ^o

のうち最初のステップをシステムの入力

u(k)

として加える．

S 4. S 1

に戻る．

1

ゲイン

K

をできる限り小さく設定することで終端状態制約集合

X ¯ f

を大きくすることができ，

(23)

2.5

数値シミュレーション

2.5.1

剛体宇宙機の姿勢制御

(RW+RCS)

剛体衛星の姿勢運動は次の運動方程式で表せる．

J p ¨ = Lu

^(2.25)

ここでは慣性モーメント行列

J ∈ R ^3×3

のみを考え，減衰項

D

や剛性項

K

は存在しないものとする．

p ∈ R ³

はロール，ピッチ，ヨー角度を含むベクトル，

L

このモデルに対し，Table2.1をシミュレーションパラメータとして，3ケースに分けて数値シミュレーションを行いオンオフモデル予測制御の有効性を検証する．

数値シミュレーションはYALMIP[32℄とMPTtoolbox[33℄を用いてMatlab環境上で実行した．

ケース1(終端評価行列

P

の比較):

まず，RW のみを用いた連続値入力の MPC 問題を考え，終端評価行列

P

を式(2.13)のLMIより求めた場合と，従来手法である式 (2.23)のDAREの解を用いた場合を比較する．Fig.2.3はロール角の1[deg℄のステップ応答とそれに伴う入力応答を示す．他の軸については目標角度は全て零を設定している．またこの時の姿勢誤差のノルムをFig.2.4に示す．両方のMPCにおいて目標角度に収束し漸近安定となっていることが確認できる．また他軸との力学的カップリングについては無視できるものであることもわかる．しかしながら，DAREを用いたMPCの方の応答は振動的であり入力飽和が度々生じていることが確認できる，また結果的には両者の整定時間は同様になっていることも分かる．これは，DAREの一意の解が入力制約に対して十分大きくないことに起因し，終端評価が小さく見積もられているためである．これに対しLMI(2.13)より得られた終端評価はゲイン行列

K

を小さく見積もることで，終端評価行列

P

をより大きく設定することができる．

結果的にLMIをもとにしたMPCのほうがDAREをもとにしたMPCよりも設計に対して柔軟性を持ち優れた結果となることを示した．

ケース2(アクチュエータの組合せ，ステップ応答):

ケース1ではRWのみでの制御に関して2つの異なる終端評価を用いてモデル予測制御を行った．ケース2では，ロール角に対する1[deg℄へのステップ応答について，RCSのみを用いた場合，RCSとRWを併用した場合について，そしてRW

のみを用いた前回の結果Fig.2.3(LMI)を比較し，それぞれのアクチュエータの組み合わせに着目し制御性能の比較を行う．Fig.2.5の姿勢角及び入力応答結果から，

(24)

RCSを用いた姿勢マヌーバのみが入力状態安定(ISS)となっていること，またRW

単体での姿勢マヌーバは漸近安定とはなるがRCSに比べて，応答が遅く望ましい性能が得られないことが確認できる．ISS安定となるRCSのロール角姿勢誤差は

約0.05[deg℄であった．RCSとRWを併用した場合のみ両者のアクチュエータ特性

を活かし，応答が速く，かつ漸近安定となっていることが確認できた．Fig.2.6には制御性能の評価として姿勢誤差のノルムを示す．Fig.2.7にはそれぞれの入力量について図示した．RCSは最初に10[se℄ほど作動し，状態量を安定領域近傍に移動させ，それ以降はRWのみでシステムを漸近安定としていることが確認でき

る.Fig.2.8にはMatlab環境下ではあるが各ステップでの最適化時間を示す．明ら

かにサンプリング周期1[se℄以内で最適化計算が収まっていることがわかる．

ケース3(アクチュエータの組合せ，高速スイッチングマヌーバ):

最後に今まで提案したMPC(RW,RCS,RW+RCS)を用いて，ロール角度を目標角度0[deg℄と10[deg℄とし，100秒周期毎に反復するスイッチングマヌーバ問題を考える．ただし，ピッチ角，ヨー角については初期0[deg℄を維持することを制御目標とする．明らかに，RWのみで制御しようとする場合，繰り返し指令に対して応答が遅いために追従できていないこと，RCSのみでの制御の場合は追従はするが偏差が残ってしまうことがFig.2.9から確認できる．このことは角度誤差に関するFig.2.10からも確認することができる．RCSとRWを併用した際のそれぞれの入力はFig.2.11に示す．これらの結果から全ての要求を満たしているのはRCSと

RWを併用する場合のみである．

(25)

Parameter Parametervalue

Inertia

J [kgm ² ]







6.0701 −0.0163 −0.0974

−0.0163 6.9443 0.0119

−0.0974 0.0119 8.9646





 × 10 ⁴

Preditionsteps

N

¹⁰

Samplingtime

[sec]

¹

RCSoutput

[Nm]

¹⁰⁰

RWoutput

[Nm]

^0.04

Stageostmatrix

Q

^=10E-1

I n

^,

R

^=10E-3

[10,0,0℄(100[se℄)

→

[0,0,0℄(200[se℄)

→

[10,0,0℄(300[se℄)

(26)

0 500 1000 1500 2000 2500 3000

−1

−0.5 0 0.5 1 1.5 2

Time[sec]

A ngl es [de g]

0 1000 2000 3000

−0.05 0 0.05

Time[sec]

Input s[N m ]

Fig2.3:Responsesofattitude(top)andontrolinput(bottom): DARE(),LMI()

(27)

0 500 1000 1500 2000 2500 3000 0

0.2 0.4 0.6 0.8 1

Time[sec]

N orm of e rror

Fig2.4:Normsofattitudeerror: DARE(),LMI()

(28)

0 50 100 150 200 250 300

−0.2 0 0.2

RW

0 50 100 150 200 250 300

0 0.5 1

RCS

0 50 100 150 200 250 300

0 0.5 1

Time[sec]

RW + RCS

0 100 200 300

−0.05 0 0.05

RW

0 100 200 300

−100 0 100

RCS

0 100 200 300

−100 0 100

Time[sec]

RW + RCS

Fig2.5:Responsesofattitude(top)andontrolinput(bottom):x(),y(---),z(

·

⁾

(29)

0 50 100 150 200 250 300 0

0.2 0.4 0.6 0.8 1 1.2

Time[sec]

N orm of e rror

Fig2.6:Normsofattitudeerror: RW(),RCS(---),RW+RCS(

·

⁾

(30)

0 100 200 300

−100 0 100

RCS

0 100 200 300

−0.05 0 0.05

RW

Time[sec]

Fig2.7:ControlinputsofRCS(top)andRW(bottom):x(),y(---),z(

·

⁾

0 50 100 150 200 250 300

0 0.1 0.2 0.3 0.4 0.5

Time[sec]

Com put at ion t im e[s ec ]

· · ·

(31)

0 50 100 150 200 250 300

−0.2 0 0.2

RW

0 50 100 150 200 250 300

0 5 10

RCS

0 50 100 150 200 250 300

0 5 10

Time[sec]

RW + RCS

0 100 200 300

−0.05 0 0.05

RW

0 100 200 300

−100 0 100

RCS

0 100 200 300

−100 0 100

Time[sec]

RW + RCS

Fig2.9:Attituderesponses(top)andontrolinputs(bottom):x(),y(---),z(

·

⁾

(32)

0 50 100 150 200 250 300 0

0.2 0.4 0.6 0.8 1 1.2

Time[sec]

N orm of e rror

Fig2.10:Normsofattitudeerror:RW(),RCS(---),RW+RCS(

·

⁾

(33)

0 100 200 300

−100 0 100

RCS

0 100 200 300

−0.05 0 0.05

RW

Time[sec]

Fig2.11:ControlinputsofRCS(top)andRW(bottom):x(),y(---),z(

·

⁾

(34)

2.5.2

剛体宇宙機に対する

3

値量子化

(

オンオフ

)

と n ^{値量子化によ} る制御

前節で求めた系3を剛体宇宙機の姿勢制御問題へ適用し，3値量子化(オンオフ)

入力と

n

値量子化入力を比較して効果を確認する．

量子化入力

ここでは，制御アクチュエータが正負方向の

n

値量子化(Fig.2.2)出力が可能なスラスタの組として考える．すなわち，

u _i = α

β u _max , α = 0, ±1, . . . , ±β

^(2.26) ここで，

α/β

はスラスタの噴射のデューティ比を表す．

β = 1

のとき

u

は不感帯を持つオンオフ制御となる．それ以外の時は，

u max

に関して

β

段階の入力に分割される階段値入力となる．ただし

β

は正の整数である．これを入力制約集合

Υ = {u _i |u _i = ^α _β u _max , α = 0, ±1, . . . , ±β}

とすると，最適化問題は次の整数二次計画問題に書き直すことができる．

min U { U (k) ^T Φ U (k) + φ ^T U (k)}

s.t.

x : x(k) ∈ Ξ, ∀k

u : u(k) ∈ {−β, −β + 1, . . . , 0, . . . , β − 1, β} ^m ⊆ Υ, ∀k U : [u(k) ^T , . . . , u(k + N − 1) ^T ] ^T

(2.27)

ここで，

Φ = (u max /β) ² E ^T Q E + R φ = 2( F x(k)) ^T Q E u _max /β

シミュレーションモデル

剛体衛星問題を簡略化して次の数学モデルを得る．

J p ¨ = Lu

^(2.28)

ここで，

p ∈ R ³

はロール，ピッチ，ヨー角度を含むベクトル，

L

は入力係数行列である．上記の式を一定時間サンプリング間隔で離散化し，次式を得る．

x(k + 1) = F x(k) + Eu(k)

y(k) = Hx(k)

^(2.29)

(35)

ケース1(状態フィードバック)

最適化問題(2.27)に基づく，状態フィードバックによるモデル予測制御のシミュレーション結果を示す．アクチュエータの特性は，式(2.26)に従うものとし，3値量子化(オンオフ)制御

(β = 1)

，

n

値量子化制御

(β = 10)

の場合を考える．また，比較のために決定変数に連続値をとり，

(β = 10)

の分解能で整数丸めを行うモデル予測

制御(PWM)の結果も併せて示す．全ての制御則で制御サンプリング周波数は1Hz

とし，予測ステップ数はN=10とする．評価関数の重み行列は，

Q = I 6 , R = 0.1×I 3

とした．

Fig.2.12，Fig.2.13にロール角を10[deg℄姿勢変更した際の応答を示す．それぞれ上から順に，ロール，ピッチ，ヨーの角度と制御トルクである．またTable.2.2，

Table.2.3にはそれぞれ，サンプリング周波数を1Hzに固定し予測ステップ数を変

化させた場合と，予測ステップ数を固定しサンプリング周波数を変化させた場合の

t f = 100[sec]

の

x(t f )

のノルムを示す．漸近安定ではないが原点を含む領域に留まっていることがわかる．各ステップにおける評価関数の最小値の時間応答を

Fig.2.14に，オンライン最適化に要する計算時間をFig.2.15に示す．

n

値量子化制御がオンオフ制御より高性能であること，Matlab使用での結果であるが，オンライン計算時間はサンプリング周期1[se℄以内に収まっていることがわかる．また，

(β = 10)

とPWMを角度精度について比較すると，本手法のオンオフ値を用いて

直接最適化するモデル予測制御の方が制御性能が高いことが確認できる．

ケース2(出力フィードバック)

次に角度のみが観測できる場合を仮定し，出力フィードバックによるモデル予測制御のシミュレーション結果を示す．出力フィードバックを実現するために，

ˆ

x(k + 1) = F x(k) + ˆ Eu + G(y − H x(k)) ˆ

^(2.30)

で表される状態オブザーバを用いる．ここで

G

はオブザーバゲインである．ケース1と同条件での結果をFig.2.16からFig.2.19に示す．オンライン計算時間が増加していること以外は，ケース1とほぼ同様の結果である．

(36)

PreditionStep(N)

β = 1

^:

kx(t f )k ² β = 10

^:

kx(t f )k ²

^PWM:

kx(t f )k ²

1

10.1E+1 6.46E+0 6.55E+0

3

14.9E−2 1.84E−2 4.12E−2

5

14.9E−2 1.84E−2 4.50E−2

10

14.9E−2 1.84E−2 4.50E−2

15

14.9E−2 1.84E−2 3.68E−2

Table2.3:DiffereneofsamplingfrequenywithN=10

Samplingfrequeny[Hz℄

β = 1

^:

kx(t f )k ² β = 10

^:

kx(t f )k ²

^PWM:

kx(t f )k ²

1

14.9E−2 1.84E−2 4.50E−2

2

3.13E−2 9.76E−3 8.77E−3

5

3.83E−2 5.85E−3 9.01E−3

10

2.54E−2 8.73E−4 5.73E−3

(37)

0 20 40 60 80 100

−5 0 5 10

Time[sec]

Roll angle[deg]

β=1 β =10 PWM

0 20 40 60 80 100

−0.025

−0.02

−0.015

−0.01

−0.005 0 0.005 0.01

Time[sec]

Pitch angle[deg]

β =1 β =10 PWM

0 20 40 60 80 100

−0.02

−0.01 0 0.01 0.02 0.03 0.04 0.05

Time[sec]

Yaw angle[deg]

β =1 β=10 PWM

Fig2.12:Responseofattitudeangles(ase1)

(38)

0 20 40 60 80 100

−150

−100

−50 0 50 100 150

Time[sec]

Torque[Nm]

β=1 β =10 PWM

0 20 40 60 80 100

−150

−100

−50 0 50 100 150

Time[sec]

Torque[Nm]

β =1 β =10 PWM

0 20 40 60 80 100

−150

0 20 40 60 80 100

−2 0 2 4 6 8 10

Time[sec]

Roll angle[deg]

β=1 β =10 PWM

0 20 40 60 80 100

−0.025

−0.02

−0.015

−0.01

−0.005 0 0.005 0.01

Time[sec]

Pitch angle[deg]

β =1 β =10 PWM

0 20 40 60 80 100

−0.02

−0.01 0 0.01 0.02 0.03 0.04 0.05

Time[sec]

Yaw angle[deg]

β =1 β=10 PWM

Fig2.16:Responseofattitudeangles(ase2)

(41)

0 20 40 60 80 100

−150

−100

−50 0 50 100 150

Time[sec]

Torque[Nm]

β=1 β =10 PWM

0 20 40 60 80 100

−150

−100

−50 0 50 100 150

Time[sec]

Torque[Nm]

β =1 β =10 PWM

0 20 40 60 80 100

−150

−100

(43)

2.5.3

フォーメーションフライト

オンオフ入力を用いたMPCのもう一つの応用例として，複数機の人工衛星のフォーメーションフライト問題を考える．ここでは衝突回避則を状態制約として与えてシミュレーションを行う．

複数機の運動方程式

まず，円軌道の1機の人工衛星の運動はCWH方程式(付録B.2)により次の状態方程式で記述できる．



 

 

˙ x =

"

0 3 I 3

−K −G

# x +

"

0 3

I ₃

# u y = x

(2.31)

ここで，

x = [p ^T , p ˙ ^T ] ^T

であり，

p ∈ R ³

は相対位置ベクトル，

p x

は軌道面に対する接線方向，

p _y

は法線方向，

p _z

は地心方向である．また

G, K

は

G =







0 0 −2n v

0 0 0

2n v 0 0





 , K =







0 0 0

0 n ² _v 0 0 0 −3n ² _v





 .

n v

は軌道角速度である．状態方程式(2.31)は次式のように複数機のフォーメーションフライトのモデルに拡張することができる．

(

˙

x m = I s ⊗

"

0 3 I 3

−K −G

#

x m + 1 _s ⊗

"

0 3

I ₃

#

u m

^(2.32)

ここで，

⊗

はクロネッカー積を示し，添え字

s

は宇宙機の個数である，

x m ∈ R ^6s

と

u _m ∈ R ^3s

はそれぞれ拡張した状態ベクトル及び入力ベクトルである．モデル予測制御ではこの状態方程式(2.32)を一定サンプリング周期の離散時間モデルに変換し，離散状態方程式(2.1)に変換して使用する．

衝突回避則

次に，衝突回避則の導入のために，まず2機の衛星に対する立方体状の安全領域を考える．

p ^l

と

p ^lm

はそれぞれ2機の位置ベクトルを表し，

or

条件を用いて次の条件を得る．

|p ^l _x − p ^lm _x | ≥ σ x or |p ^l _y − p ^lm _y | ≥ σ y or |p ^l _z − p ^lm _z | ≥ σ z

^(2.33)

(44)

バイナリー変数

δ i

と十分大きい正の実数

M

を導入して，

or

条件を

and

条件へと書き直して次の条件が成り立つ[34℄．

p ^l _x − p ^lm _x ≥ σ x − Mδ 1

p ^lm _x − p ^l _x ≥ σ x − Mδ 2

p ^l _y − p ^lm _y ≥ σ _y − Mδ ₃ p ^lm _y − p ^l _y ≥ σ y − Mδ 4

p ^l _z − p ^lm _z ≥ σ z − Mδ 5

p ^lm _z − p ^l _z ≥ σ _z − Mδ ₆

6 X

i=1

δ i ≤ 5

(2.34)

この

and

条件を行列形式の不等式としてまとめると

Lx ≥ s − MI 6 d , 1 ^T ₆ d ≤ 5

^(2.35)

となる．ここで

L

は位置ベクトル

p ^l

と

p ^lm

の組み合わせを与える行列で，

s

は

σ {x,y,z}

に関するベクトルである．また

d

はバイナリ変数の組をなすベクトルである.従って式(1.4)と式(2.35)は最適化変数

u(k)

に対する拘束となる. これらの結果から，

式(2.34)を予測ステップ数

N

分拡大することによって最適化問題は次の形でまと

められる．

min U , D { U (k) ^T Φ U (k) + φ ^T U (k)}

s.t.

x : x(k + 1) = F x(k) + Eu(k) U : u(k), u(k + 1), . . . , u(k + N − 1) u i : u i ∈ {−1, 0, 1}

D : L { F x(k) + EU (k)} ≥ S − MD ( 1 _N ⊗ I n p ⊗ 1 ₆ ^T ) D ≤ 1 _n _p _N ⊗ 5

(2.36)

ここで，

Φ = E ^T Q E + R

，

φ ^T = 2( F x(k) − X _r ) ^T Q E

とおいた．また，

n p

は組み合わせ数．係数行列

L , M

と拡大ベクトル

S

はそれぞれ，

L =





 L 1

L ₂

.

L _N







, M = diag(MI 6n p , . . . , MI 6n p ) ∈ R ⁶ⁿ ^p ^N×6n ^p ^N , S =





 s ₁ s ₂

.

s _N







(2.37)

となり，

L 1 , . . . , L N

および

s ₁ , . . . , s _N

は各ステップ毎の制約を表し，Dはバイナリー変数組である．

(45)

−20 0

20 −20 0 20

−1

−0.5 0 0.5 1

X−axis Z−axis

Y−axis

S1 S2

Fig2.20:3-Dplotofresponses

2機の宇宙機のフォーメーションフライト

数値シミュレーションでは，式(2.32)に従う2機(S1とS2)の宇宙機のフォーメーションフライトを考える．

数値モデルのパラメータとして

µ = 3.986 × 10 ¹⁴ [m ³ /s ² ]

，円軌道半径を

r = 7000

[km℄，入力の最大値を

0.02

^[N/kg℄として，システムは零次ホールドサンプリングとして

1

^[se℄周期で離散化したモデルを用いる．

評価関数のパラメータは

N = 10

，

Q= I n

，

R = 0.01I m

に設定した．衝突回避の安全領域

σ

は3方向全てに対して1[m℄とした．それぞれの衛星の初期位置は衛星 1(S1)，衛星2(S2)に対してそれぞれ

[−10, 0, −10]

と

[−10, 0, 10]

^[m℄を設定し，

目標位置はそれぞれ

[10, 0, 10]

と

[10, 0, −10]

^[m℄とする制御を行う．この問題設定ではS1とS2の位置は軌道平面上で交わることになるため，衝突回避を要する問題である．

Figs.2.20-2.21はそれぞれの衛星の軌道を表す．Fig.2.23はそれぞれの衛星の目標位置との差をノルムを用いて評価する．Fig.2.20とFig.2.23の結果から，衝突回避を達成し目標位置近傍に収束させることができたことがわかる．Fig. 2.22では全ての衛星のスラスターが目標位置近傍においても噴射を続けていることが確認できるが，このフォーメーションシステムは入力制約と安定性を満たしている．

電気通信大学大学院 情報理工学研究科 博士

入力状態安定性を保証する量子化入力による 宇宙機のモデル予測制御

淺川 岳也