モデル予測制御入門 3：離散時間の非線形MPC

←前章：制約付きの非線形最適化問題で学んだ最適化の手法（主双対内点法）を用いて実際非線形ダイナミクスを制御してみましょう。

離散時間の非線形ダイナミクスの例

今回は下の式で与えられる減衰振動する離散時間の一次元調和振動子に3次の補正項が入ったようなモデルを考えます。

$$z_{t+2} = 1.750 z_{t+1} - 0.902 z_t - 0.01 z_t^3 + u_t$$

初期値を$z_0=1, z_1=1$とし、制御入力を0（$u_t \equiv 0$）とすると下図のようなパターンで減衰振動をします。

この振動子をMPCで$z = 1.0$に止めてみましょう。

最適化問題の定式化

まずはダイナミクスを一階の差分方程式で表す必要があるため、新たな変数$x_t=(z_{t+1}, z_t)^T$を導入して次のように運動方程式を書き換えます。

$$\begin{align*} x^{(0)}_{t+1} & = 1.750 x^{(0)}_{t} - 0.902 x^{(1)}_{t} - 0.01 {x^{(1)}_{t}}^3 + u_t \triangleq f^{(0)}( x_{t}, u_{t})\\ x^{(1)}_{t+1} & = x^{(0)}_{t} \triangleq f^{(1)}( x_{t}, u_{t}) \end{align*}$$

また、制御入力の範囲は

$$-0.2 \leq u_t \leq 0.2$$

に限定しましょう。
制約条件は与えられたので、次にコスト関数$\mathcal{L}$を決めます。$z=z_{goal}$に静定させるとしたら、予測ホライズン長$T$に対して

$$\mathcal{L} = \alpha_z \sum_{k=0}^T (x^{(0)}_{t+k} - z_{goal})^2$$

が考えられる最も簡単なコスト関数になります。（重み係数$\alpha_z>0$）
しかし、このままでは制御入力を激しく変動させることが可能になってしまうため、$\Delta u_t = u_t - u_{t-1}$に対する正則化を追加し、またオーバーシュートを防ぐために速度$v_t = x^{(0)}_{t} - x^{(0)}_{t-1}$に対する正則化も追加して

$$\mathcal{L} = \alpha_z \sum_{k=0}^T (x^{(0)}_{t+k} - z_{goal})^2 + \alpha_{u} \sum_{k=1}^{T-1} \Delta u_t^2 + \alpha_v \sum_{k=1}^{T-1} v_t^2$$

とします。

MPCによる制御

最適化問題の定式化が出来たので、次はこれを各時刻毎に解いて実際に制御してみましょう。
実装はPythonのcasadiで行いましたが、これを使うと2章で説明した主双対内点法を用いた最適化ソルバー「IPOPT」が簡単に扱えるので非常に便利です。

実際に制御を行ってみた様子が下の動画です。

制御の際に用いた各パラメータは以下の通りです。

• ホライズン長 $T=10$
• $\alpha_z = 0.1$
• $\alpha_u = 0.5$
• $\alpha_v = 0.1$
• $z_{goal} = 1$
• $z_0=0, ~ z_1=0$

動画でもわかるようにハンチングもオーバーシュートも起きずに、また制約条件をきちんと守って制御できていることが確認出来ました。

さて、非常に簡単な問題ではありましたがMPCを用いた制御が出来るようになりました。

次の章では連続時間での非線形ダイナミクスを制御してみましょう。連続時間を直接扱うことはできないため結局は上手く離散化して近似して最適化問題を解くことになるのですが、そのテクニックの一つであるdirect collocation法を解説いたします。

次章：連続時間の非線形MPC