多項式によるモデル予測制御器の設計法とその性質
M2012MM008 市原寛之
指導教員：大石泰章

1

はじめに

表 記 を 簡 単 に す る た め 問 題 Ox の 目 的 関 数 を(u0 , u1 , u2 , u3 , u4 ) の み を 使って 表 し ，さ ら に こ れ を平方完成して次の標準形 Px に変形をする:

モデル予測制御は産業界で最も広く使われている制御法の一つであり，制約を持つ制御対象にも適用することができる. Bemporad ら [1] は，最適な入力をあらかじめ1 T現在の状態の関数として求めておくことを提案し，これPx : minimizez Hz2によってモデル予測制御は応答の速い制御対象にも適用
subject toGz ≤ f + Sx
できるようになった．しかし最適なモデル予測制御器は，状態の区分線形な関数になり，区分の数が大きい場合に，Ox の最適化変数 (u0 , u1 , u2 , u3 , u4 ) と Px の最適化変数 z
その計算は時間がかかり，実装も難しくなる．そこで多とはアフィン変換の関係にあり，Px の最適解が得られれ項式を使って，近似をすることにより実用的なモデル予ば Ox の最適解が得られる.
測制御器を作るという方法がある [2, 3]．しかしこの方法
文献 [3] では Px の Lagrange 乗数 λ を x の多項式に限は適用例が少なく，その実用性には疑問が残る．
定して λ(x) と書き, 以下の最適化問題を考える:
以上を踏まえ，本研究ではこの方法を今までと異なる制御対象に適用し，その検証を行う．また最適なモデル予∫測制御器を求めるソフトウェアとして定評のある multiS : minimize
λ(x)T (f + Sx + GH −1 GT λ(x))dx
parametric toolbox [4] を使い，多項式近似法との比較をX行う.
subject to − GH −1 GT λ(x) ≤ f + Sx (∀x ∈ X),

2

モデル予測制御器の多項式近似法

2.1

定式化

本研究では次の線形離散時間システムを制御対象とする:

x(t + 1) = Ax(t) + Bu(t) (t = 0, 1, 2, . . . ).
ただし，

(A=)
( )
1 10,B =0 11であり，u(t) は次の制約を満たす:−1 ≤ u(t) ≤ 1.
モデル予測制御では時刻 t における入力 u(t) は以下の手順で定める．まず現在時刻における制御対象の状態 x(t)を観測し，これを以下の最適化問題 Ox のパラメータ xに代入して最適解 (u0 , u1 , u2 , u3 , u4 ) を求める．ただしこ
こでは予測ステップ数を 5 としている．そして u0 を u(t)として制御対象に入力する．時刻 t + 1 になったら再び制御対象の状態を観測し, 以下同様に繰り返す. 考える最適化問題は以下の通り:

Ox : minimize4∑

xTk xk + u2kk=0

subject tox0 = x,xk+1 = Axk + Buk (k = 0, 1, 2, 3, 4),− 1 ≤ uk ≤ 1 (k = 0, 1, 2, 3, 4).

λ(x) ≥ 0 (∀x ∈ X)
ただし X は状態空間中の制御器を求める状態の集合で，ここではX = {(x1 x2 )T | − 5 ≤ x1 ≤ 5, −5 ≤ x2 ≤ 5}とする． この問題の実行可能解 λ(x) に対して z(x) =−H −1 GT λ(x) によって多項式 z(x) を定めると，z(x) は任意の x ∈ X に対して Px の実行可能解となり，その時の S の目的関数の値は z(x) が達成する Px の目的関数値と Px の最小値との差を X 上で積分したものの上界となる．従って問題 S を解けば X 中の x に対するモデル予測制御器の多項式による近似を誤差評価つきで求めることができる．

2.2

解法

問題 S は次のようにして解ける．問題 S の目的関数は次のように変形できる．まず λ(x) は x の多項式であるから，λ(x) = Λµ(x) と書ける．ただし Λ は λ(x) の係数，µ(x) は λ(x) の単項式を並べた列ベクトルである．この表現を使うと，∫λ(x)T [f + Sx + GH −1 GT λ(x)]dx
X∫=[µ(x)T ΛT f + µ(x)T ΛT Sx+Xµ(x)T ΛT GH −1 GT Λµ(x)]dx= ξ T ΛT f + tr(S Ξ̃ΛT ) + tr(H − 2 GT ΛΞΛT GH − 2 )．
11

ここで，ξ T ,Ξ̃, Ξ はそれぞれ µ(x)T ，xµ(x)T ，µ(x)µ(x)Tを X の範囲で積分したものである. この目的関数の最小化は次のように表現することができる:

minimize ξ T ΛT f + tr(S Ξ̃ΛT ) + trV()1
V
H − 2 GT ΛΞsubject to⪰01
ΞΛT GH − 2Ξ問題 S の制約条件はすべての x ∈ X に対して満たされなくてはならないので扱うのが難しい．しかし，2 乗和多項式行列 (例えば [5]) を使うことで扱うことができる．
2 乗和多項式行列とは，変数 x = (x1 , x2 ，· · · ，xn ) に関する多項式行列 F (x) で，他の多項式行列 R(x) を用いてF (x) = R(x)T R(x) と書くことができるようなものを言う．F (x) が 2 乗和多項式ならば任意の x に対して F (x)は半正定値となる．
多項式行列 F (x) が領域 X 上で半正定値となるための十分条件を 2 乗和多項式を使って表すことを考える．領域 X は有界閉集合で多項式 g1 (x), · · · , gl (x) によってX = {x ∈ Rn |g1 (x) ≥ 0, g2 (x) ≥ 0，· · · gl (x)，≥ 0}

このように Ξ の各成分に大小の差があるために, 最適化問題が数値的に解きにくくなり，前述のような結果になったと考えられる．この点を明らかにするために本研究では Ξ の条件数を調べてみた．

3.2

条件数

状態の範囲が [−5 , 5] × [−5 , 5] の場合の Ξ の条件数は1.0870 × 108 となり， [−1 , 1] × [−1 , 1] の場合は 5.2765 ×103 となった．
この 2 つの状態による Ξ の条件数の違いは約 105 であり，非常に大きい．条件数は大きければ大きいほど数値計算上扱いづらくなる．この事が理由で領域 X が [−5 , 5] ×[−5 , 5] の場合に数値誤差が生じたと考えられる．

3.3

対策とその結果

このことを解決するため, x̃1 = x51 ，x̃2 = x52 という変数変換を問題 Px に適用した. 変数変換後の問題 S は次の形となる:と表せるとする．そして，与えられた多項式行列 F (x)S : minimize
が，2 乗和多項式行列 S0 (x), · · · , Sl (x) を使って F (x) =S0 (x) + g1 (x)S1 (x) + · · · + gl (x)Sl (x) と表せるとする．
subject to
これは F (x) ⪰ 0 (∀x ∈ X) であるための十分条件である．
さらに，次の定理も成り立つ．
定理 1 多項式 g1 (x), g2 (x), · · · , gl (x) がある緩い条件を ただし満たす時，F (x) ≻ 0 (∀x ∈ X) ならば，を満たす 2 乗和多項式行列 S0 (x), S1 (x), · · · , Sl (x) が存在する．

2 乗和多項式を使って，問題 S の制約条件をその十分条件で置き換える．特に問題 S の場合 g1 (x) = (x1 +5)(5−x1 )，g2 (x) = (x2 + 5)(5 − x2 ) とすればよい．この条件を使うことで問題 S を近似的に解くことが可能となる．また
2 乗和多項式の次数が高いほど近似精度がよくなる．

多項式近似法の適用

3.1

問題点と考察λ(x̃)T (f + 5S x̃) + GH −1 GT λ(x̃))dx̃

X̃− GH −1 GT λ(x̃) ≤ f + 5S x̃(∀x̃ ∈ X̃)λ(x̃) ≥ 0 (∀x̃ ∈ X̃)X̃ = {(x̃1 x̃2 )T | − 1 ≤ x̃1 ≤ 1, −1 ≤ x̃2 ≤ 1}.

F (x) = S0 (x) + g1 (x)S1 (x) + · · · + gl (x)Sl (x)3∫

上記の変数変換を行った結果，問題 S を解くと値は 1509となった．前述の結果と比べると絶対値は小さくなり正の値となった．領域 X を分割すればさらに目的関数を小さくすることができる．例えば X を [−5，0] × [0，5], [0, 5] ×[−5, 0]，[−5，0] × [−5，0]， [0，5] × [0，5] の 4 つの領域に分割すると目的関数値は 3.7042 となる．この場合に得られたモデル予測制御器を図 1 に示す．X の領域の分割の上にそれぞれ入力関数 u(t) が描かれている．図 1 は−1 ≤ u(t) ≤ 1 となっており妥当な結果が得られたと言える．

問題 S をそのまま解いたところ，目的関数の値は−1.146 × 1011 となった．問題 S の目的関数の値は z(x)が達成する Px の目的関数値と Px の最小値との差を X 上で積分したものの上界となるので，零以上の値とならなければならない．そこで数値誤差の関与が疑われる．この事が起きた原因としては，行列 Ξ の各成分にあると考えられる．
定義に従って計算すると，Ξ の (1, 1) 成分は∫ 5∫ 5−5−51dx1 dx2 = 100である．Ξ の (28, 28) 成分は∫ 5∫ 5−5−59x121 dx1 dx2 = 1.878 × 10である．

図 1 多項式近似法によるモデル予測制御器

4最適なモデル予測制御器

multi-parametric toolbox[4] は最適なモデル予測制御器を求めるソフトウェアとして定評があるものである．このソフトウェアを 2 章の制御対象に適用したところ図 2の制御器を得た．図 1 と似ており多項式近似法は十分良い近似ができている事がわかった．この場合の制御器は状態 X の区分 1 次関数であり，図 3 に求めた区分を示す．
図 3 の領域 1 つ 1 つの上にそれぞれ 1 次関数が定義される．図 3 を見ると，予測ステップ数が 5 と小さい場合にも分割数が 32 と多くなることが分かった．

図 4 シミュレーションにおける比較

図 2 最適なモデル予測制御器

このシミュレーションでは原点に近づくにつれ近似が悪くなる．また原点の周りでの破線の動きに違和感がある．そこで今度は分割を 9 つにし，近似を行った．
使った 分 割 は [−5, −1.7] × [−5, −1.7]，[−5, −1.7] ×[−1.7, 1.7] ，[−5, −1.7]×[1.7, 5]，[−1.7, 1.7]×[−5, −1.7],[−1.7, 1.7] × [−1.7, 1.7],[−1.7, 1.7] × [1.7, 5], [1.7, 5] ×[−5, −1.7],[1.7, 5]×[−1.7, 1.7],[1.7, 5]×[1.7, 5] である．この分割によって得られた目的関数値は 2.3228 である．目的関数値は 4 分割の場合に比べ小さくなった．得られた多項式モデル予測制御器を図 5 に示す．

図 3 状態の分割
図 5 多項式近似法によるモデル予測制御器

5

シミュレーション

5.1

4 分割の場合のシミュレーション

領域 X を [−5，0] × [0，5], [0，5] × [−5，0]，[−5，0] ×[−5, 0], [0, 5] × [0，5] のように 4 分割して，多項式近似法を適用する．得られたモデル予測制御器を使って，初期値を (1.8 1.8)T とした時のシミュレーション結果を図 4に示す．実線が multi-parametric toolbox における初期値の軌道で，破線が多項式近似法における初期値の軌道である．

そして，このモデル予測制御器を使い，再びシミュレーションを行った．実線が multi-parametric toolbox における初期値の軌道で，破線が多項式近似法における初期値の軌道である．

7

図 6 シミュレーションにおける比較 (9 分割の場合)
目的関数値は 4 分割の場合に比べ小さくなっており，先ほどより良い近似ができている．ただしまだ近似精度は十分ではない．

6

実行時間の測定おわりに

本研究ではモデル予測制御器の多項式近似法を文献とは異なる制御対象に適用した．数値誤差の対策を行い，多項式モデル予測制御器を出力した．さらに，multiparametric toolbox で求めた最適なモデル予測制御器との比較を行った．制御則の比較では十分良い近似が得られているように見えたが，シミュレーションを行ったところ，あまり良い近似となっていなかった．そのため領域の分割数を増やし，再びシミュレーションを行ったとこ
ろ，若干の改善は見られたが依然として十分な精度は得られなかった．分割数を増やすことによって目的関数値が小さくなったことより，分割数をさらに増やすと近似精度がさらに良くなると思われる．また多項式モデル予測制御器を出力するまでの実行時間の測定を行った結果，ステップ数が増えるごとに実行時間は増え，多項式近似法は実行時間が長く改善が必要なことが分かった．2 乗和多項式は一般に大きな計算量を必要とすることが知られており，この結果はこのことの反映であると思われる．

参考文献[1] A．Bemporad，M．Morari，V．Dua，and E．N．Pistikopoulos: The explicit linear quadratic regulator
for constrained systems．Automatica，vol. 38(2002)，
no. 1，pp. 3–20．

本研究では多項式モデル予測制御器の時間に関する性質を調べるために，様々な予測ステップ数に対してモデ[2] Y．Oishi: Direct design of a polynomial model preル予測制御器を出力するまでの実行時間の測定を行った．
dictive controller．Proceedings of the 7th IFAC Sym結果は図 7 の通りである．
posium on Robust Control Design，Aalborg, Denmark，June 2012，pp． 621–626．

[3] Y． Oishi: Simplified approaches to polynomial
design of model predective controllers．Proceedingof the 2013 IEEE Multi-Conference on Systems
and Control，Hyderabad，India，August 2013，pp．
960–965．
[4] M．Kvasnica，P．Grieder，andM．
Baotic:Multi-parametric toolbox (MPT)．
http://control.ee.ethz.ch/~mpt/，2004．
[5] 大石泰章: 2 乗和多項式に基づくロバスト制御系設計．
計測自動制御学会セミナーテキスト: 実践的な制御系
設計: (ポスト) ロバスト制御の最前線，大阪，2007
年 11 月，pp. 71–88．

図 7 ステップ数の変動による実行時間の変化
このグラフから，多項式近似法の実行時間は長く，しかもステップ数が増えると共に，実行時間も急速に増えることがわかる．2 乗和多項式は一般に大きな計算量を必要とすることが知られており，上の結果はこのことの反映であると思われ，改善が必要なことが分かった．