モデル予測制御を用いた目標値追従制御
M2011MM070 竹本徹
指導教員：高見勲

1

はじめに化学プラントなどのプロセス制御で既に広く使われているモデル予測制御 (Model Predictive Control) は, リア
ルタイムで最適制御問題を解くことで得られる制御入力を用いた制御法である. 有限時間の最適制御問題を解くためには相応の計算時間が必要となる. そのため従来は化学プラントなどの比較的遅いシステムに採用されていた.
予測制御は制約を考慮した最適化という側面と調整が容易で直観的であるという側面も持っている. 計算機の高速化と同時に, 最適化アルゴリズムの向上のおかげで高速に凸最適化問題の解を求めることが可能になった. このような背景から発展を遂げ, 現在では比較的遅いシステムから速いシステムへの適用と研究が活発化しており, 発展が期待されている [1].
モデル予測制御は有限時間内で N ステップ区間ごとに制御入力の最適値を算出し, さらに区間をずらしそこを初期状態とした制御入力の最適値を算出する. このような手順を繰り返すため, 外乱などの影響を受けてもそこからの目標軌道への最適な制御入力の値を再度算出しなおすことができる. このためモデル予測制御は経路計画や目標軌道への追従などのより高度な制御を実現することが可能である [2].
以上を踏まえ, 本研究ではサーボ系の目標値追従問題に対する, モデル予測制御系設計方法を検討し, その成果を２次元ボールバランサに適用する事を目的としている. ２次元ボールバランサはプレートの上においたボールをギヤ操作を用いて所望の位置に移動させたり, 意図した軌道を描かせることが出来る. これらのことにモデル予測制御は特に有効であり, またモデル予測制御の有用性を検証するのに,2 次元ボールバランサは十分の実験装置であると考えている.

2

図 1 制御対象モデリング

3

２次元ボールバランサは x 軸方向,y 軸方向で独立したモータを搭載している. このため x 軸方向,y 軸方向を別々に制御することが出来る. 本研究では x 軸方向のモデリングを行い, 制御性能を検証する.
ボールの重心位置を x[m], プレートの一辺の長さをLp [m], ロードギアの仰角を θl,x [rad], レバーアームのオフセットを rarm [m] として２次元ボールバランサの x 軸方向の概略図を図 2 に, ギヤトレインの概略図を図 3 示す.

xLoadGear

２次元ボールバランサは, 独立した２つの DC モータを駆動させそれぞれに備えられたギヤを操作し, 上部に設置された正方形のプレートを操り, ボールを所望の位置に移動させたり軌道を描くことができる実験装置である. この２次元ボールバランサの特徴として, プレートの上部にCCD カメラが搭載されており, ボールの位置情報を検知しその情報をフィードバックし制御を実現するビジュア
ルフィードバックシステムを備えた実験装置である.
２次元ボールバランサを図 1 に示す.

rb
θ l,x

制御対象

Lprarmαx

図 2 ２次元ボールバランサの x 軸方向の概略図

3.1

制御系

状態空間表現を式 (1) とする.
{żx (t) = Azx (t) + Bux (t)yx (t) = Czx (t)(1)

制御は有限時間区間の最適制御問題を, 時刻が進むにつれて評価時間区間を先にずらして進めていくということを行う制御であり Receding Horizon Control とも呼ばれる.
今日 モデル予測制御は比較的応答の遅い化学プラントの定温・定圧制御などのプロセス制御において, もっとも応用された高度な制御方式となっている. また, 近年では, 高度なソフトウェアの完備によって機械システムという速い動特性をもつシステムへ適用させることが可能となってきている. しかも, オープンループの最適問題を入力や状態の制約を加味しながら解くことも可能となった [4].
図 3 ギヤトレインの概略図
ラグランジュの運動方程式より導出した A, B, C 行列を式 (2) とする.

0
 0
A= 00100000016rarm g5Lp00
B=t

B

状態量を zx (t) =(2)0θl,x (t) θ̇l,x (t)

制御量を yx (t) = x(t), 操作量をモータ電圧 ux (t) =Vm,x (t) とする.
なお, y 軸に関する状態空間表現もプレートが正方形でありギヤなども同じものを使用しているため,x 軸と同様に表すことができるので省略する [3].
ここで, それぞれの変数を Rm [m] はモータ抵抗, g[m/s2 ]は重力加速度, Am [Nm/A] はアクチュエータゲイン,
Beq [kg · m2 /s] は等価粘性摩擦係数, Jeq [kg · m2 ] は総合慣性モーメントとする.

モデル予測制御

4.1

基本概念

tTTt

T

モデル予測制御の基本概念を図 4 を用いて示す. 現在
時刻 t から有限時間先の T までにおいて最適な制御入力列とその予測値を求める. そして得られた最適な入力列の中で最初の制御入力 u(t) のみを現在時刻 t において入力する (図中の黒塗りの部分). 次の時刻では, 再び現在時刻の状態を初期値として現在時刻 t から有限時間先の T までにおいて最適な制御入力列とその予測値を求め, その中で最初の制御入力 u(t) のみを入力する. さらに次の時
刻では, 同様な手順を繰り返していく. 各時刻においては,開ループの最適制御問題を解いているが, 各時刻において初期値をフィードバックして考えていることから, 開ループ系から閉ループ系にすることができる [5]. それを摸擬したものが図 5,6 である.

]T4u
入力列



[x(t) ẋ(t)u

入力列y参照軌道

図 4 基本概念]0ry

[1 0ur参照軌道

入力列

AmRm Jeq

C=y参照軌道eq− Jeq000r

モデル予測制御とは, 各時刻において有限時間区間先の未来までの最適な入力軌道をオンラインで計算し,1 ステップ分, すなわち最初のサンプル時間分だけ最適な入力をシステムに加え, 遷移した状態を初期状態とし, 再び同じ有限時間区間先までの最適な入力軌道を計算し, 1 ステップ分最適な入力をシステムに加えるということを繰り返し行う. 最適な入力というのは, 制御対象に対して, ある評価関数を定めて, その評価関数の値 (コスト) を最小化するというオープンループ最適制御問題を解くことで得られる. この操作をサンプル時間ごと繰り返し行うことで,フィードバック制御則を生成する. このようにモデル予測

r Optimal ux
Plant
Controlr Modelxu
Predictive Plant
Control

図 5 最適制御

図 6 モデル予測制御

モデル予測制御を適用した２次元ボールバランサにおける簡単な考え方を以下の図 7 に示す．
制御対象である２次元ボールバランサのプレートの上にボールが乗っている様子である. 黒く塗りつぶしているのが制御されるボールである. また, 破線がシミュレーションにおける目標軌道であり, 実線が実際にボールが動いた軌道である.
ここで, 図 7 における左側の図は原点にボールが停止しており, その状態から円運動における目標軌道を予測している状態である. この状態から入力である電圧が加わることでモータのギア操作によりボールが動き出す. ボールが動き出した様子を表しているのが右側の図である. ボール

この時最適性の必要条件は, 随伴変数を λ = ∂J∂x , ハミルT
トニアンを H = L + λ f と定義すると第一変分の停留条件から

λk (t) = λk+1 (t) + HxT [xk (t), λk+1 (t), uk (t)]∆τ (7)図 7 モデル予測制御
は予測された目標軌道にあう入力が印加され円運動を開始するが, 何らかの影響により目標軌道を外れてしまうと仮定する. ここで, あくまでも円運動を続けるために, 前もって予測された軌道に戻ろうとする. しかし, ボールは速やかに目標軌道には戻れない. したがって, この目標軌道からずれた位置の状況から再び円運動に戻るたに算出された経路から修正されたコースをたどる. 元の目標軌道に戻り, 再び円運動を開始する. このように右側の図は「制御量の将来挙動の予測―予測が理想とする動きに近くなるような操作量修正―挙動の予測」という手続きの繰り返しを行い制御を実現する. これらがモデル予測制御のアルゴリズムの基本である．

4.2ẋ = f (x(t), u(t))x(0) = x0







Hu [xk (t), λk+1 (t), uk (t)] = 0(9)F [U (t), x(t), t] :=HuT [x0 (t), λ1 (t), u0 (t)]HuT [x1 (t), λ2 (t), u1 (t)]..
.
HuT [xN −1 (t), λN (t), uN −1 (t)]


=0
(10)

この連立方程式の解 U (t) が評価関数 (4) 式の意味の最適性の必要条件を満たす入力列となる. よってこの連立方
程式の解の U (t) のうち, 初めの値の U (t)1 をシステムを入力する事で制御を行う.

xmin ≤ x(t + 1) ≤ xmax

(11)

ymin ≤ y(t + 1) ≤ ymax

(12)umin ≤ u(t) ≤ umax

(13)

さらに状態 x, 入力 u, 出力 y にそれぞれ制約条件を付加する事により制御性能を向上する事が出来る [8].

4.3

2 次元ボールバランサでの制御系設計

２次元ボールバランサでの評価関数を (14),(15) 式とする
∫ t+T1J(u, t) =l(x(τ ), u(τ ))dt(14)2tl(x(τ ), u(τ )) = ((x(τ ) − xref (τ ))T Q(x(τ ) − xref (τ ))+ u(τ )T Ru(τ ))(3)

この時, 各時刻 t において T 秒だけ未来までの評価関数を(4) 式とする.
∫ t+TJ = φ(t + T ) +L(x(τ ), u(τ ))∆τ(4)(8)のように得られる. よって, 制御入力の系列 U (t) と初期状態 x(t) を与えれば,(5),(7),(8) 式より状態変数と随伴変数Nの系列 xi (t)Ni=0 , λi (t)i=0 が決まる. これより (9) 式を U (t)に関する方程式とみなし, 次のように表す.

制約条件付きのモデル予測制御

評価区間を長く取ることで, 無限時間の最適制御問題がもつ安定特性やロバスト性など望ましい特性を得られる
が, 計算量が多くなり, 機械系などの動特性が速いシステ
ムに適用できない [6]. モデル予測制御は有限時間区間の最適制御問題とホライズンが無限大である最適制御問題の流れから発展した考え方である. モデル予測制御の評価関数は考えているシステムによって, 自由に設定することができるが, 評価関数を最小とする入力を計算していくということからも, 評価関数の設定は意義のあるものでなければいけない [7].
まず, 状態方程式が次の形であらわされるシステムを考える.

λN (t) = φTx [xN (t)](15)

また, この時のハミルトン関数を (16) 式とする.

H=1((x(τ ) − xref (τ ))T Q(x(τ ) − xref (τ ))2+u(τ )T Ru(τ )) + λT (Ax + Bu)(16)t(4) 式を最小にするような制御入力を求める. しかしこのままでは解析的に入力を計算する事が困難である. したがって次のように離散化する.
xk+1 (t) = xk (y) + f (xf (t), uk (t))∆τJ = φ(xN (t)) +N−1∑

L(xk (t), uk (t))∆τ(5)(6)

k=0

ただし時間刻みは ∆τ = T /N であり,xk (t) とは t 時刻のk ステップ先の予測値. J を最小にする制御入力の系列
U (t) = {u0 (t), u1 (t), ‥‥, uN −1 (t)} を決定する.

ここで Q ≥ 0, R > 0 とする. Q は状態への重みであり,対角要素に値を持つ行列とする. また,R は入力への重みで, 制御が強い場合には数値を大きくし, 逆に制御が弱い場合には数値を小さくする.
xref (τ ) は目標値である. 状態と目標値の偏差をとり評価する事により目標値へ追従させる.
これらのパラメータを試行錯誤により決定させ, シミュレーションに用いる.

Q = diag(10, 100, 1, 1)(17)

R=1(18)

状態に対する制約条件をプレートの大きさとし, ボールが移動できる距離を制限する.
入力電圧にも制約条件をかけ, 過度な入力電圧の発生を抑える.

シミュレーション

5

Matlab/Simulink に搭載されている toolbox である
Multi-Parametric Toolbox[9] を用いてシミュレーションを行った.
Multi-Parametric Toolbox を 用 い る こ と で m-ﬁle やSimulink 上で自由に制約条件などの設計が出来る [10].
. x 軸,y 軸共に同じ制御系で動作するため,x 軸方向のみのシミュレーションを行った. そのシミュレーション結果を示す.
原点から x 軸正方向に 0.12m 移動させたシミュレーショ
ン結果を図 8 に, またこの時の入力 u を図 9 に示す.
目標値に速やかに追従し, オーバーシュートのない安定したシミュレーション結果を得た. 入力は制約条件の中におさまり制御を実現している. これをもとに実験装置においても適用させる.
0.2

voltage[v]position[m]0.10.080.06

x−position

0.040.02
0012time[s]3u−0.204

図 8 x 軸方向のデータ

1u-constraint-0.25≦u≦0.25
234time[s]図 9 入力 u のデータ

実験結果

6

実験でも同様に Matlab/Simulink に搭載されているtoolbox である Multi-Parametric Toolbox[9] を用いて行った.
シミュレーション時と同様に原点から x 軸正方向に 0.12m移動させた実験結果を図 10 に, またこの時の入力 u を図11 に示す. ２次元ボールバランサの特徴として急激な入力電圧の印加によって予期せぬ挙動を起こす事がある. したがって入力に対する制約条件において実験装置を動かせる最低限度程度の条件を課した. その数値が 0.25[V] である.

0.25
0.2

0−0.05x−posotion−0.1

x-constraint-0.13≦x≦0.13

−0.15
−0.2
00.5

time[s]11.5

図 10 x 軸方向の実験結果

voltage[v]position[m]

0.1
0.05

おわりに

本研究の成果を以下に示す.
サーボ系の目標値追従問題に対し, モデル予測制御を用いて制御系を設計する事が出来た.
また, 実験において直線運動の制御性能を検証し, サーボ系においてモデル予測制御が有効である事を立証した.

参考文献[1] Jan M.Maciejowski(著)，足立修一 (訳)，管野政明(訳)：『モデル予測制御-制約のもとでの最適制御(Predictive Control with Constraints)』. 東京電機
大学出版局，2005

0−0.1

x-constraint-0.13≦x≦0.13

7[2] 大嶋正裕， 小河守正：『モデル予測制御-I: 基礎編: 発展の歴史と現状』．システム/制御/情報 :システム制
御情報学会誌 46(5), 286-293, 2002-05-15.

0.120.1

シミュレーション時よりも速応性が高く目標値に追従している.
データ上ではオーバーシュートしているように見てとれるがこれは制御対象であるボールが完全な球体でないため目標値に到達後プレート上で多少ぶれた為に観測されていると考えられる. 実験における印加電圧が制約条件内の 0.25[V] におさまっており, 制御を実現しているの事が図 11 より確認できる.

0.1
0.05
00u-constraint-0.25≦u≦0.25
0.5
1time[s][4] 加納 学，大嶋 正裕：『モデル予測制御-II : 線形モデ
ル予測制御』．システム/制御/情報 : システム制御情報学会誌 46(7), 418-424, 2002-07-15.
[5] 大嶋正裕， 加納学：
『モデル予測制御-IV: 制約条件付きのモデル予測制御』．システム/制御/情報 : システム制御情報学会誌 46(11), 702-708, 2002-11-15.
[6] M.Lazar，W.P.M.H.Heemels，S.Weiland，A.Bemporad： 『Stabilizing Model Predictive
Control of Hybrid System』 ．IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 51,
NO. 11, NOVEMBER 2006
[7] 五十嵐裕司：『情報構造を持つ離散線形システムに対する分散 Receding Horizon 制御』東京工業大学卒業
論文,2006[8] Jaroslav Pekar，Vladimir Havlena：『DESIGN AND
ANALYSIS OF MODEL PREDICTIVE CONTROL USING MPT TOOLBOX』Czech[9] M.Kvasnica,P.Grieder,M.Baotic,F.Christophersen：『Multi-Parametric Toolbox(MPT)』Czech,March30,2006

u0.15

−0.05

[3] 西脇; 『最適レギュレータを用いた２自由度ボールバランサの位置決め制御』南山大学卒業論文, 2011.

1.5

図 11 入力 u の実験結果[10] 向井正和， 畑中健志， 藤田政之：『ハイブリッドシステムの制御-IV：モデル予測制御』．システム/制御/
情報 : システム制御情報学会誌 51(7), 148-155, 2007.