C/GMRESによる非線形モデル予測制御1

今回は珍しく2部構成で，本記事ではシステムの立式から評価函数の設定まで，次の記事では具体的な解き方について述べます．

第1部となる本記事ではC/GMRESの核となる部分には言及しないため，そちらについては2部以降を見ていただければと思います．

---

非線形モデル予測制御(Nonlinear Model Predictive Control; NMPC)とは，その名の通り”状態方程式が非線形なシステム”に対するモデル予測制御(Model Predictive Control; MPC)です．

最近流行りの制御則で，車の自動運転やロケットの再着陸システムにも使用されているそうです．

MyEnigmaさん
モデル予測制御(Model Predictive Control:MPC)の応用例

そもそも制約有りのMPCは非線形な制御系ですが，”非線形MPC”といった場合には“制御対象が非線形”な系を指すように思います．

線形のモデル予測制御については別途，記事に起こすつもりです．

通常，線形MPCでは，システムの表現に以下のような線形な状態方程式を使うのに対して，
$$\dot{x}(t) = Ax(t)+Bu(t)$$
非線形MPCでは状態方程式の形式を問わず使用でき，以下のような形式で表現されます．
$$\dot{x}(t)=f(x(t),u(t)) = \left[\begin{array}{c} x_{2}(t) \\ ax_{1}(t)+bx_{2}(t)u_{1}(t) \end{array}\right]$$

これはセミアクティブダンパの状態方程式ですが，式中に$x_{2}(t)u_{1}(t)$という“状態変数と操作量の積”が含まれています．少し考えるとわかりますが，線形の状態方程式に書き直すことはできません．

実際，このシステムを制御する際に単純に線形MPCを適用することはできません．
このような，非線形なシステムの制御を目的としたMPCが非線形MPCです．

非線形MPCの中にも種々ありますが，ここではC/GMRESによる非線形MPCについて取り上げます．

ただ，このシステムは減衰が操作量なのでダンパとして動作させること自体はPID等でも可能です．というよりずっとONにしておけばよいのです．しかし，できるだけ少ない操作量で最大の減衰を得よう，などと考え始めるとちゃんと制御する必要があります．

---

C/GMRESは連続変形法(Continuation method)と一般化最小残差法(Generalized Minimum RESidual method)の組み合わせであることが語源です．

シー・ジーエムレスと読むそうです．

他のモデル予測制御と比較して，以下の利点が挙げられます．

1. 非線形なシステムを直接扱える
MPCを含む多くの制御系では，制御対象システムを動作点周りで線形化したり，座標変換による線形化を行う必要があります．

2. 非線形MPCの中では計算負荷が小さい
従来の制御方法と比較し，MPCはそもそも実時間の最適化を含むため計算負荷が大きい制御手法です．
非線形MPCとなるとさらに計算負荷が増えますが，C/GMRESはこれらの手法と比較して計算負荷が小さいです．

3. 制約（拘束条件）を取り扱える
MPCであれば当たり前のように思えますが，C/GMRESでもシステムの制約条件を取り扱うことができます．
これにより，使用できる限界ギリギリの”効率が最高となる稼働条件”での操業を行うことができます．

MPCでの制約の取扱いについては足立先生が翻訳された本が詳しいです．
モデル予測制御 制約のもとでの最適制御

---

C/GMRESは非線形制御の第一人者である大塚敏之先生が開発した手法です．
この先生の著作である“実時間最適化による制御の実応用”という書籍が非常にわかりやすく，本記事ではこれを実装するまでの流れについて書いていきます．
実時間最適化による制御の実応用

ここからは先にも述べましたが，当該書籍に例題として載っている”セミアクティブダンパ”を対象として説明していきます．

セミアクティブダンパは図に示すとおりですが，操作量として減衰を与えます．
減衰は速度に対して作用し，操作量の大きさによって速度の減少具合を制御します．

セミアクティブダンパを用いる利点としては，以下のようなものがあるようです．

1. アクチュエータが不要なため作りが安価である．
2. 基本的には受動動作であるため，エネルギ消費が小さい．

ここではアクチュエータそのものに関しては議論しません．

このようなシステムでは操作量が速度の減少に直結するため，最も高速に減衰させたい場合には，常に最大の操作量（減衰）を与えておけば良いことになります．

しかし，2のようなエネルギ消費量を小さくすることを目的とするのであれば，きちんと制御する必要があります．
以下では，“制御することが必要である”というモチベーションを前提とします．

---

さて，制御対象であるセミアクティブダンパが以下の式で表されるとします．
$$\left\{\begin{array}{lll} \dot{x}(t) &=& f( x(t), u(t) ) \\ C( u(t) ) &=& 0 \\ \end{array} \right.$$

実際のシステムで書き直すと，
$$\left\{\begin{array}{lll} f( x(t), u(t) ) &=& \left[\begin{array}{c} x_{2}(t) \\ ax_{1}(t)+bx_{2}(t)u_{1}(t) \end{array}\right] \\ C( u(t) ) &=& \left( u_{1}(t) - \frac{u_{max}}{2} \right)^2 + u_{2}(t)^2 - \frac{{u_{max}}^2}{4} \end{array}\right.$$
ただし，状態変数$x(t)=[Position(t), Velocity(t)]^T$です．
また，$a=-k/m$，$b=-1/m$で$m$は質量，$k$はばね定数です．

ここで，上側の式は状態方程式でシステムの特性を表します．
下側の式は制約で，システムへの操作量の制限について書かれています．

C/GMRESでは，状態方程式を状態変数と操作量で分ける必要がありません．

$\dot{x}(t) = Ax(t)+Bu(t)$でなくてもよいのです．

---

また，制約は等式制約として与えられます．

$C( u(t) ) = 0$を見てみると，これは$x^2+y^2=r^2$の式そのもので，$( u_{1}(t), u_{2}(t) )= ( u_{max}/2, 0 )$の点を中心とした，半径$u_{max}/2$の円となっています．

このとき，$u_{2}(t)$をダミー入力と呼びます．スラック変数と呼ぶこともあります．
状態方程式を見て分かる通り，システム特性を定義する式中に$u_{2}(t)$は存在しません．
つまり，実際にはシステムに影響を与えない操作量なのです．

これは，umin <= u1(t) <= umaxのような不等式制約を，等式制約として表現するための方法です．
実際にはu1(t), u2(t)の点は円の上に乗っているのですが，u2(t)を無視すると，
u1(t)は0からumaxまでを往復するような値を持っているように見えます．

最適化では$u_{1}(t), u_{2}(t)$の両方を求めますが，$u_{2}(t)$はどこにも使用されないという変数になります．

---

次に，評価函数を決めます．
モデル予測制御をされている方にとっては当たり前の函数ですが，改めて説明しておきます．

評価函数は，設計者がどのようにシステムを動かしたいかを定式化するものです．
モデル予測制御では，この評価函数を最小化するような操作量を計算します．

評価函数の多くは，以下のような2次形式で与えられます．

$$J = \frac{1}{2}x^{T}(t+T)S_{f}x(t+T)+\int_{t}^{t+T}\left\{ \frac{1}{2}\left(x^{T}(\tau)Qx(\tau)+r_{1}{u_{1}}^2(\tau)\right)-r_{2}u_{2} \right\} d\tau$$

右辺第1項は，現在の時間$t$から$T$秒だけ進んだ未来，$t+T$秒の状態変数の大きさを表します．
これを終端コストと呼びます．
$S_{f}$は設計パラメータで，たいていの場合，対角行列として与えられます．

$$\frac{1}{2}[Position(t+T), Velocity(t+T)]\left[ \begin{array}{cc} 10 & 0 \\ 0 & 1 \end{array}\right][Position(t+T), Velocity(t+T)]^T$$

このように設定すると，$Position$は$Velocity$の10倍の係数が掛かることになります．
評価函数$J$は$Position$が少しでも大きくなると，$Velocity$を大きくしたときよりも容易に大きくなります．
$J$を最小化することを考えたとき，$Velocity$よりも$Position$を小さくしたほうが効果的であることがわかるかと思います．

また，$Position$，$Velocity$は2乗されるため負の値を持ちません．
このことから，$t+T$秒後に$Position$，$Velocity$ともにゼロに近いことが最小の条件になります．

つまり右辺第1項は，$T$秒後に状態変数をゼロに近づけてね（ただし$Position$は優先的に小さくしてね）という函数なのです．

ダンパの力で振動を止めたい，と考えるときに，未来の時刻で位置と速度がゼロに近いというのは直感に反しないのではないかと思います．

---

一方，右辺第2項では，ある函数が$t$秒から$t+T$秒まで積分されています．

$$x^{T}(\tau)Qx(\tau)$$

上式は，右辺第1項と同様に状態変数の2次形式です．$Q$も設計パラメータで，同様に対角行列で与えることが多いです．
これを積分するということは，$t$秒から$T$秒の間に渡って，ずっと状態変数はゼロに近いほうがよい，という意味になります．

この式は負にならないため，これが最小になるということは$t$秒から$t+T$秒の間，全体的にずっとゼロに近いということになります．

この式がない場合には，T秒後まではどれだけ激しく振動していてもよいので，T秒の時点ではピタッと止めてね，という条件になります．これも直感的な条件かと思います．

---

次に，右辺第3項ですが，操作量$u_1$に関する条件です．

$$r_{1}{u_{1}}^2(\tau)$$

これも2次形式ですが，今回の例では操作量は1つなのでスカラとなり，単に2乗となっています．
$r_1$は設計パラメータで，これを大きくすると入力エネルギの総和を小さくする意味になります．
化学プラントや重量物の制御などではシステムの急激な変化を避けるため，このパラメータを大きく取る傾向にあるように思います．

一方で，電磁弁（ソレノイド）などのON-OFF動作のシステムでは，このパラメータを小さく取り，操作量の変化を過敏にすることでON-OFF入力での最適化操作量を擬似的に扱うこともあるようです．

大塚先生本人のsurvey
実時間最適化による非線形機械システムのモデル予測制御

また，不等式制約を扱うことで実現する方法もあるようです．

ホバークラフトのモデル予測制御
Nonlinear Receding Horizon Control of an Underactuated hovercraft

---

最後に右辺第4項ですが，以下の式で表されます．
$$r_{2}u_{2}(\tau)$$
他の項については2次形式で負の値を取らず，係数$r_2$も設計パラメータという扱いでしたが，右辺第4項については設計パラメータの項ではなく，計算を安定させるための項となっています．

$u_{1}(t)$にとって$u_{max}$は鞍点ですので，$u_{2}(t)=0$の点で安定してしまいます．

せっかくダミー入力$u_2$を導入することで不等式制約を等式制約で表現できたのに関わらず，その係数がゼロでは意味がありません．
かといって，あくまでシステムの動作に影響しない操作量ですので，大きい値を持ちたくありません．

$r_2$として計算が失敗しない程度の小さな値を与えておくとよいみたいです．

右辺の括弧内の式をまとめてステージコストと呼びます．

---

時刻$t$を起点とした評価函数$J$において，予測区間でのある時刻$\tau$を考えます．

評価函数$J$は，状態$x$を変数とした2次形式の函数であるので，どこかに$J$を最小化する状態$x$が存在します．

また，操作量$u$についての2次形式でもあるので，$J$を最小化する操作量$u$も同様に存在します．

これらの概念を図にすると，以下のような曲面になります．

この曲面上で$J$を最小化するような操作量$u$を求めることが，最適化の役割です．

実際には1変数1入力の場合には完全に下に凸な曲面（鞍点を持たない）になるような気もしますが，そのあたりはご容赦ください．

---

ここまでで，制御対象となるシステムと評価関数について定義しました．

つまり，システムが本来どのように動くかとシステムをどのように動かしたいかを定義したことになります．

評価函数の挙動について少し述べていきます．

以下のような評価関数のパラメータ設定をしてみます．

$$Q = diag( 1, 10 ) \\ S_f = diag( 1, 10 ) \\ r_1 = 1, r_2 = 0.01$$

このとき，アクティブダンパの動作は，以下のとおりです．

操作量$u_1$によって，位置$x_1$と速度$x_2$が減衰していることが見て取れます．
操作量$u_1$は，速度が大きくなるタイミングを見計らって，効果的にシステムへ入力されていることも見ることができます．

---

次に，操作量$u_1$の重み$r_1$を小さくした以下のパラメータを考えます．

$$Q = diag( 1, 10 ) \\ S_f = diag( 1, 10 ) \\ r_1 = 0.01, r_2 = 0.01$$

このとき，応答は次のように変化します．

位置$x_1$，速度$x_2$ともに先程の図よりも素早く減衰しています．
一方で，操作量$u_1$は$0.3$付近のまま，なかなか減少しません．

---

これらの評価函数（設計パラメータ）の使い分けは設計者がどのようにシステムを動かしたいかに懸かります．

減衰が遅くてもエネルギ消費の少なさを求めるなら前者でしょうし，速く減衰させたいだけであれば後者でしょう．

そもそも，評価函数の形式すら本記事で述べたとおりである必要はないため，この設定が設計者の勘所になってくるのでしょうか．

---

さて，状態と操作量の積を有する非線形系のモデル予測制御について，制御対象のモデリングから評価函数の設定まで行いました．

線形系であればここから2次計画法などのアルゴリズムで解けますが非線形系ではそうもいかないので，ハミルトン函数を利用して解くようです．

続きは第2部で，近日中に公開を目指します． (2017/09/19)
作成したSimulnkモデルはSimulink Projectという形式でGithubにアップロードしました．
ご意見いただけるとありがたいです．

Github
C/GMRES Simulink project

（readme.mdの書き方すらわからない初心者なのでリポジトリ自体に対する意見でもお願いします）