Simulinkで始める線形システムのモデル規範適応制御(MRAC)

はじめに

本記事で対象にしている読者は

• PID制御や状態フィードバック制御の特性をある程度理解して，感覚的であれ，安定余裕や極配置などの定量指標に基づいてであれ，制御器のパラメータチューニングができる。
• ところが，制御対象の特性が(制御動作に比べて緩やかながらも)時間変化することが分かり，想定されるどんな状況でも要求仕様を満たせるロバストなパラメータを置いたり，状態に応じてパラメータを切り替えたり(ゲインスケジューリング)する必要に迫られている。
• しかし，ロバストなパラメータが極端に応答の遅い保守的なものになってしまったり，パラメータ切り替えのために何通りものパラメータ設定と切り替え条件の検討に工数を割かなければいけなかったり，パラメータ切り替え条件でチャタリング起こすんじゃないかと心配したり，で途方に暮れている。

ような制御設計者の方(というか最近の筆者)です。
そういう方向けに，今回は代表的な適応制御手法とSimulinkでの実装例を紹介します。

話の流れとしては，一般的なテキストのように理論の積み上げから実装方法を説明するのではなく，とりあえず動かせるSimulinkモデルの例をお見せして，それから適宜理論の補足をしていければと思います。

使ったツールは下記の通りです。

• Matlab
• Simulink
• Control System Toolbox

なお，世の中にある様々な制御手法の中での適応制御の位置付けに関しては，こちらの記事が分かりやすかったです。
適応制御の基礎

目標

先に述べたようなロバストパラメータ探しやゲインスケジューリングをしなくてもいいように，基本的な適応制御手法の1つであるモデル規範適応制御(MRAC)をとりあえず動かせるところを目指します。制御対象は1入力1出力の線形システムとします。また，安定性にも気を使ったパラメータ設計も行えるようにします。

ただし，ベースとなる制御はフィードフォワード制御＋状態フィードバック制御 ($u(t) = \theta_r r(t) + \theta_x^\top x(t)$) に限定します。状態$x(t)$が$[x_1(t), \dot{x}_1(t)]^\top$のように定義されていればPD制御とほぼ等価になります。
「定常偏差を消したいからI制御も含めて適応制御を使いたい！」という気持ちもあるかと思いますが，今回I制御に相当する仕組みは割愛します。本記事では触れませんが，一般的には外乱推定をしてそれを打ち消すような制御項を加えるアプローチが取られるようです。

SimulinkでMRAC

線形システムをモデル規範適応制御(MRAC)で制御するモデルをSimulink上で構成しました。

大きく，制御対象(Linear Plant)とその状態方程式の係数を変動させる(Plant Varying)ブロック，そして適応制御ブロック(Adaptive Controller)に分かれています。
簡単のために制御器が状態を直接取得していますが，現実的にはセンサやオブザーバが挟まっていると思ってください。

制御対象

今回は状態方程式

\begin{align}
\dot{x}(t) &= A x(t) + b x(t) \\
y(t) &= c x(t)
\end{align}

で表される線形システムを対象にします。

なお，適応制御器構成のためにこれらの係数の値が既知である必要はありませんが，
$b_0 = c A^{(n^*-1)}b$ ($n^{*}$は相対次数¹) の符号 $\mathrm{sign}(b_0)$ だけは知っている必要があります。

具体的な係数と初期状態は下記のように設定しました。2階の線形微分方程式で表せる制御対象であり，自由応答は不安定です。

A = [0, 1; 2, 3];
b = [0; 1];
c = [1, 0];
initial_x = [0; 0];

>> eig(A) %実部が正の固有値があれば不安定

ans =

   -0.5616
    3.5616

この係数行列は後ほど時間変化させてみます。

適応制御器

適応制御器の中身は次のようになっています。

基本構造

大きく，規範モデル，ベースとなる制御部，制御パラメータ更新部に分かれています。

規範モデルは制御目標値$r(t)$に対する理想的な応答$y_M(t)$を作るモデルです。この理想的な応答と実際の応答をぴったり一致させるのが"モデル規範"適応制御の目的です。参照モデル自体は状態方程式

\begin{align}
\dot{x}_M(t) &= A_M x_M(t) + b_M r(t) \\
y_M(t) &= c_M x_M(t)
\end{align}

やその伝達関数形式で記述します。

ベースとなる制御はシンプルであり，フィードフォワード制御と状態フィードバック制御の組み合わせ，$u(t) = \theta_r(t) r(t) + \theta_x(t)^\top x(t)$です。Simulinkモデル内では制御パラメータを$\Theta(t) = [\theta_r(t), \theta_x(t)^\top]^\top$，目標値と状態を$\omega(t)=[r(t), x(t)^\top]^\top$とまとめています。

そして，その制御パラメータ$\Theta$を適応的に調整するのがパラメータ更新部であり，まさに"適応"制御のコアとなる部分です。
パラメータの更新は規範モデルの理想的な応答$y_M(t)$と実際の応答$y(t)$の誤差$e(t) = y_M(t)-y(t)$を0に近づけるように行われます。

規範モデル設計

実際に規範モデルを作ってみます。制御仕様や要求スペックに基づいて思い思いの規範モデルを設計しましょう。

制御対象が2次なので，今回は次のような二次遅れ系を規範モデルとします。
$$\frac{1}{U(s)} = \frac{1}{\frac{s^2}{\omega_n^2} + \frac{2\delta}{\omega_n}s + 1}$$

$\delta$は減衰率であり，ステップ応答において，1未満だとオーバーシュート，1以上だと強制減衰します。

$\omega_n$は固有角周波数であり，ステップ応答において，大きいほど立ち上がりが早くなります。

今回は$\delta = 1$, $\omega_{n} = 5$ としておきます。

なお，この後適応制御のハイパーパラメータを決めるために，規範モデルを状態方程式表記したときの係数を使います。tf2ssを使って各種係数行列を出しておきましょう。

>> delta_M = 1;
omega_n_M = 5;
Ws_num = [1/omega_n_M^2, 2*delta_M/omega_n_M, 1];
[A_M, b_M, c_M] = tf2ss(1, Ws_num) % 1/W(s)

A_M =

   -10   -25
     1     0


b_M =

     1
     0


c_M =

     0    25

適応則設計

適応制御のコアの部分，パラメータの更新則は次のようにします。
$$\dot{\Theta} = \mathrm{sign}(b_0)\Lambda \omega(t) h^\top{\bf e}(t)$$

各ハイパーパラメータを適切に設定すれば制御誤差$e(t)$を0に収束させることが保証されるやり方です。

$b_0$は$b_0 = c A^{(n^*-1)}b$で表される値です ($n^{*}$は相対次数¹)。$A, b, c$が未知な場合計算できませんが，$b_0$の符号さえ分かればよいので大雑把な値が分かっていればよいと思います。

$\Lambda$ は正定行列です。今回は$\omega(t)$が3次なので，3×3行列です。機械学習の学習率に相当するものなので，値が大きいほどパラメータ更新速度が速くなります。
実はいくら大きくしてもMRACの安定論的には問題ありません。（おそらく制御器を離散時間化した場合に更新幅が大き過ぎてパラメータが暴れるという問題が起きると思います。）

$h$は誤差ベクトル${\bf e}(t) = [e(t), \dot{e}(t)]^\top$に乗じてスカラーにする係数です。今回は2×1ベクトルです。
パラメータ更新を含めた制御を安定化させるために，$h^\top (sI-A_M)^{-1}b_M$が強正実になるように決めます。
「強正実」についてはこちらで定義と必要十分条件をご確認ください。

鈴木隆. "正実性 (III) 正実性と適応制御." 計測と制御 34.11 (1995): 897-904.

必要十分条件に従えば，$h=[h_1, h_2]^\top$として，$10h_1 > h_2, h_2 > 0$ を満たせばよいことが分かります。

相対次数$n^*$を知るために今回の制御対象を伝達関数表記します。
分子分母係数は次のようになり，

>> [num, den] = ss2tf(A, b, c, 0)

num =

     0     0     1


den =

     1    -3    -2

相対次数は2と分かります。

今回適応則に関する各種パラメータを次のように設定しました。

n_relative = 2;
b0 = c*A^(n_relative-1)*b;
Lambda = 1*eye(n+1);
h = [1; 1]; % h'*(sI-A_M)^-1*b_Mが強正実となるように選ぶ

制御結果

ステップ入力に対する応答を見て，MRACの効果を確認します。
具体的には，目標値$r(t)$を10秒おきに1→0→1→0...とステップ上に変化させたときの応答を確認します。

以下の4つの例を見てみましょう。
* 制御対象変動なし，適応制御なし
* 制御対象変動あり，適応制御なし
* 制御対象変動あり，$\Lambda$小 MRAC
* 制御対象変動あり，$\Lambda$大 MRAC

制御対象変動なし，適応制御なし

まずは制御対象の変動がない場合について，適応制御を使わず，しかしパラメータを適切に場合について動作を確認します。

応答が規範モデル$A_M$に合うように極配置できるフィードバックゲイン$\theta_x(t)$を決めましょう。
placeコマンドを使ってかっこよく...

>> place(A, b, eig(A_M))
エラー: place (行 78)
"place" コマンドは、rank(B) より大きい多重度を持つ極を配置できません。

求めたいところでしたが，$A_M$の固有値が$\lambda_1=\lambda_2=-5$の二重解であるゆえに，連立方程式が作れず，フィードバックゲインが一意に決められなくて怒られたのだと思います。

少し強引ですが，ぎりぎり重解じゃなくなるように

>> place(A, b, [-5, -5+0.00001])

ans =

   26.9999   13.0000

とおいて，フィードバックゲイン$\theta_x$は$[-27, -13]^\top$ということにしましょう。

（状態フィードバック系の固有値を求める方程式
$$|\lambda I - (A+b\theta_x^\top)|=0$$に$\lambda=-5, -5+\Delta$をそれぞれ代入して$\theta_x$を求めて，$\Delta \rightarrow 0$としても良いです。（いいのか？））

フィードフォワードゲインは，状態フィードバックゲインとプラントの係数行列を用いて

 \theta_r =
\begin{bmatrix}
 -\theta_x^\top & 1 
\end{bmatrix}
\begin{bmatrix}
 A & b \\
 c & 0
\end{bmatrix}^{-1}
\begin{bmatrix}
 0 \\
 1
\end{bmatrix}

とすることで，0でない目標値$r(t)$にも追従できるようになります。

>> [27, 13, 1]*inv([A, b; c, 0])*[0; 0; 1]

ans =

    25

今回は25にすれば良さそうです。
というわけで制御パラメータは$\Theta = [\theta_r, \theta_x^\top]^\top = [25, -17, -13]^\top$と決まりました。
では動かしてみましょう。

※GIFです。
規範モデル出力にぴったり一致していることが分かります。

制御対象変動あり，適応制御なし

制御対象システムの係数$A$に以下のA_diffを加算することで，連続的に・不連続に変動させてみます。

function A_diff = gen_diff(t)

A_diff = [0, 0; -2*t/200+4*sin(2*pi*t/400), -2.5*t/200+7*sin(2*pi*t/500)]; %連続な変化
A_diff = A_diff + [0, 0; 6*floor(t/600), 7.5*floor(t/600)]; % 不連続な変化

無茶苦茶な変動の仕方ですね。モデルに基づいて決めた極配置が台無し…
果たしてどのような制御結果になるのでしょうか。

発散しないだけ偉いとはいえ，応答速度，定常誤差，オーバーシュート量，何もかもが制御対象の変動によって変わってしまいます。
適切に極配置したとはいえ，制御対象が変わればとても規範モデルに追従できる様子はありません。

制御対象変動あり，Λ小 MRAC

さて，いよいよMRACを導入してみましょう

少しはマシになったような…という応答ですね。

制御対象変動あり，Λ大 MRAC

適応速度を速くするため，思い切って$\Lambda$を10倍にしてみました。

予想通り，制御対象の変動に対して頑健な制御ができていることが分かります。
GIFでは分かりづらいですが，不連続な変動に対して一瞬出てくるオーバーシュートは，そのステップ応答1回のみで消えています。
ステップ応答数回にわたって徐々に適応していくのではなく，1回のステップ中においても適応しているのは，適応速度として凄まじいものを感じます。（強化学習の適用で，何万試行もの動作を必要とした経験があるので。）

何はともあれ，これでもうしんどいゲインスケジューリングしなくても済むんだ！
めでたしめでたし。

おまけ

適応制御，今年の10月リリースのMatlab2021bからSimulink Control DesignにModel Reference Adaptive Control(MRAC)ブロックが追加されたようで，MRACがお手軽にできるようになりました。今回紹介したFF+FB方式に加え，外乱抑制項も加えられますのでより実用的な制御が可能になりそうです。
適用例を置いておきますのでよかったら触ってみてください。

• Model Reference Adaptive Control of Aircraft Undergoing Wing Rock
  (https://jp.mathworks.com/help/slcontrol/ug/model-reference-adaptive-control-of-aircraft-undergoing-wing-rock.html)
• Model Reference Adaptive Control of Satellite Spin
  (https://jp.mathworks.com/help/slcontrol/ug/mode-rReference-adaptive-control-of-satellite-spin.html)

参考文献

• 鈴木隆. "正実性 (III) 正実性と適応制御." 計測と制御 34.11 (1995): 897-904.
• 宮里 義彦. "適応制御." コロナ社 (2018).