第 2 章 モデル予測制御 11
2.5 数値シミュレーション
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 /β) 2 E T Q E + R φ = 2( F x(k)) T Q E u max /β
シミュレーションモデル
剛体衛星問題を簡略化して次の数学モデルを得る．
J p ¨ = Lu
(2.28)ここで，
p ∈ R 3
はロール，ピッチ，ヨー角度を含むベクトル，L
は入力係数行列 である．上記の式を一定時間サンプリング間隔で離散化し，次式を得る．x(k + 1) = F x(k) + Eu(k)
y(k) = Hx(k)
(2.29)ケース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とほぼ同様の結果である．PreditionStep(N)
β = 1
:kx(t f )k 2 β = 10
:kx(t f )k 2
PWM:kx(t f )k 2
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 2 β = 10
:kx(t f )k 2
PWM:kx(t f )k 2
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
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)
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
−50 0 50 100 150
Time[sec]
Torque[Nm]
β =1 β=10 PWM
Fig2.13:Responseofinputtorques(ase1)
0 20 40 60 80 100 0
0.5 1 1.5
2 x 10 4
Time[sec]
Minimum cost
β=1 β=10
Fig2.14:Timehistoryofminimumost(ase1)
0 20 40 60 80 100
0 0.1 0.2 0.3 0.4 0.5
Time[sec]
Computation time[sec]
β =1 β=10 PWM
Fig2.15:Computationtime(ase1)
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)
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
−50 0 50 100 150
Time[sec]
Torque[Nm]
β =1 β=10 PWM
Fig2.17:Responseofinputtorques(ase2)
0 20 40 60 80 100 0
2 4 6 8 10 x 10 4
Time[sec]
Minimum cost
β=1 β=10
Fig2.18:Timehistoryofminimumost(ase2)
0 20 40 60 80 100
0 0.1 0.2 0.3 0.4
Time[sec]
Computation time[sec]
β=1 β=10 PWM
Fig2.19:Computationtime(ase2)