モデル予測制御の数理 +制約の考慮

MPCはシステムモデルの入出力による予測軌道と参照軌道の最適化計算を行う．ここである時刻(ステップ)$t=k$でのシステムへの入力$\bf{u}(k)$と出力$\bf{y}(k)$が次のような入出力関係をもつとする．

$\begin{aligned}\bf{x}(k+1) &= A\bf{x}(k) + B\bf{u}(k)\\\bf{y}(k) &= C\bf{x}(k)\end{aligned}$

この時システムによっては$\bf{x}(k+1)$をそのまま$k+1$ステップ目の出力とするものも多い．この時$C=I$である．

MPCは目標とする理想軌道である参照軌道$\bf{r}(t)$に予測出力の軌道を近づけようと最適化する．ここで現在時刻である$k$ステップ目で最適化計算を実施する．$\hat{\bf{x}}$は$X$の予測値を表すとして，時刻$k$に対する$i$ステップ先の予測出力を$\hat{\bf{x}}(k+i|k)$と表記する．時刻$k$における$i$ステップ先の参照軌道を$\bf{r}(k+i|k)$，システムの予測出力を$\hat{\bf{z}}(k+i|k)$とする．MPCは有限区間$i=n$までの将来の軌道を予測し，次のコスト関数$V(k)$を最小化する入力系列$\bf{u}(t)$を求める．なお現在時刻$k$におけるシステムへの実際の入力はまだなされておらず，時刻$k$における入力はまだ予測値$\hat{\bf{u}}(k|k)$であることに注意すると，

$\begin{aligned}V(k) &= \sum_{i=0}^n||\hat{\bf{z}}(k+i|k) – \bf{r}(k+i|k)||^2 + \sum_{i=0}^{n-1}||\Delta \hat{\bf{u}}(k+i|k)||^2\end{aligned}$

である．右辺第1項は予測軌道と理想軌道の差を最小にすることを目的とする．第2項の$\Delta{\bf{u}}(k+i|k)$はある時点での入力とその1ステップ前の入力の差$u(k+i|k) -{\bf{u}}(k+i-1|k)$を表し，入力の変動も最小にすることを目的とする．

なお以降では簡単のため$\hat{\bf{z}}(k+i|k) = \hat{\bf{x}}(k+i|k)$とする．MPCの入出力モデルからステップ間の予測出力，予測入力間に、

$\begin{aligned}\hat{\bf{x}}(k+i+1|k) &= A\hat{\bf{x}}(k+i|k) + B\hat{\bf{u}}(k+i|k)\end{aligned}$

および、

$\begin{aligned}\hat{\bf{u}}(k+i|k) &= \sum_{j=0}^i\Delta\hat{\bf{u}}(k+j|k) + \bf{u}(k-1)\end{aligned}$

の関係性があることが分かる．

ここで予測ステップ数に2つの「ホライズン」という上限を設ける．1つ目の予測ホライズン(Predictive Horizon)は予測するステップ数の最大値を意味する．これを$h_p$と表記する．2つ目は制御ホライズン(Control Horizon)であり，これは予測ホライズン中の最初の何ステップだけ制御入力$\hat{\bf{u}}(k+i|k)$の値を変化させるかに相当する．制御ホライズン以降は最後の制御入力の値が一定であるとして残りの予測ホライズン中の予測出力の計算を行う．これを$h_u$と表記する．

これらを用いると前述のコスト関数は，

$\begin{aligned}V(k) &= \sum_{i=0}^{h_p}||\hat{\bf{z}}(k+i|k) – \bf{r}(k+i|k)||^2 + \sum_{i=0}^{h_u-1}||\Delta \hat{\bf{u}}(k+i|k)||^2\end{aligned}$

となる．

制御ホライズンの考えにより，制御入力は$k,k+1,…,k+h_u-1$の$h_u$個のステップでのみ変化し，以降の予測ホライズン間の$h_u \leq i\leq h_p-1$では変化しないことに注意すると，システムの予測出力$x(k+i|k)$は次のように求まる．

$\begin{aligned}\hat{\bf{x}}(k+1|k) &= A\bf{x}(k) + B\hat{\bf{u}}(k|k)\\\hat{\bf{x}}(k+2|k) &= A\hat{\bf{x}}(k+1|k) + B\hat{\bf{u}}(k+1|k)\\&= A^2\bf{x}(k) + AB\hat{\bf{u}}(k|k) + B\hat{\bf{u}}(k+1|k)\\\hat{\bf{x}}(k+h_u|k) &= A^{h_u}\bf{x}(k) + (A^{h_u-1}+…+A+I)B\Delta \hat{\bf{u}}(k|k)\\&\quad\quad+…+ (A+I)B\Delta\hat{\bf{u}}(k+h_u-1|k)\\&\quad\quad+ (A^{h_u-1}+…+A+I)B\bf{u}(k-1)\\\hat{\bf{x}}(k+h_p|k) &= A^{h_p}\bf{x}(k) + (A^{h_p-1}+…+A+I)B\Delta \hat{\bf{u}}(k|k)\\&\quad\quad+…+(A^{h_p-h_u}+…+A+I)B\Delta\hat{\bf{u}}(k+h_u-1|k)\\&\quad\quad\quad\quad+ (A^{h_p-1}+…+A+I)B\bf{u}(k-1)\\\end{aligned}$

これを適切な行列$\Psi$，$\Theta$，$ \Upsilon$，$\Delta\mathcal{U}$と$A$，$B$を用いて書き直すと，

$\begin{aligned}\left[\begin{array}{c}\hat{\bf{x}}(k+1|k)\\\vdots\\\hat{\bf{x}}(k+h_u|k)\\\hat{\bf{x}}(k+h_u1|k)\\\vdots\\\hat{\bf{x}}(k+h_p|k)\end{array}\right]&= \left[\begin{array}{c}A\\\vdots\\A^{h_u}\\A^{h_u+1}\\\vdots\\A^{h_p}\end{array}\right]\bf{x}(k)+\left[\begin{array}{c}B\\\vdots\\\sum_{i=0}^{h_u-1}A^{i}B\\\sum_{i=0}^{h_u}A^{i}B\\\vdots\\\sum_{i=0}^{h_p-1}A^{i}B\end{array}\right]\bf{u}(k-1)\\&\quad\quad+\left[\begin{array}{ccc}B & \cdots & 0 \\AB+B & \cdots & 0 \\\vdots & \ddots & \vdots \\\sum_{i=0}^{h_u-1}A^{i}B & \cdots & B \\\sum_{i=0}^{h_u}A^{i}B & \cdots & AB+B \\\vdots & \vdots & \vdots \\\sum_{i=0}^{h_p-1}A^{i}B & \cdots & \sum_{i=0}^{h_p-h_u}A^{i}B\end{array}\right]\left[\begin{array}{c}\Delta\hat{\bf{u}}(k|k)\\\vdots\\\Delta\hat{\bf{u}}(k+h_u-1|k)\end{array}\right]\\&= \Psi x(k) + \Upsilon\bf{u}(k-1) + \Theta\Delta\mathcal{U}(k)\end{aligned}$

となる．ここで$\Delta\mathcal{U}(k)$は，

$\begin{aligned}\Delta\mathcal{U}(k) = \left[\begin{array}{c}\Delta\hat{\bf{u}}(k|k)\\\vdots\\\Delta\hat{\bf{u}}(k+h_u-1|k)\end{array}\right]\end{aligned}$

であり，予測入力のステップ間の差分$\Delta \hat{\bf{u}}(k+i|k)$の関数である．さらに列ベクトル$\bf{\varepsilon}(k)$を，

$\begin{aligned}\bf{\varepsilon}(k) &=\left[\begin{array}{c}\bf{r}(k+1|k)\\\vdots\\\bf{r}(k+h_p|k)\end{array}\right] – \Psi \bf{x}(k) – \Upsilon\bf{u}(k-1)\end{aligned}$

とする．これは未来の目標軌道とシステムの応答の誤差を意味し追従誤差を表す．

また入出力への重みを与える重み行列を$Q$，$R$とする．これはそれぞれ$h_p\times h_p$行列，$h_u\times h_u$行列であり，プログラム実行中に重みを変えなければこれらは定数行列として予め求めておくことができる．

上記を用いるとコスト関数は$\Delta\mathcal{U}(k)$の関数として次のような形になる．

$\begin{aligned}V(k) &= \mathrm{const.} – \Delta\mathcal{U}^T(k)\bf{g} + \Delta\mathcal{U}^T(k)H\Delta\mathcal{U}(k)\\\bf{g} &= 2\Theta^TQ\bf{\varepsilon}(k)\\H &= \Theta^TQ\Theta + R\end{aligned}$

これより$V(k)$は凸関数であることが分かる．システムが実行中に不変であれば$2\Theta^TQ$と$H$は定数行列である．MPCではこのコスト関数を最小にするベクトル$\Delta\mathcal{U}(k)$を求め，この最初の要素である$\Delta \hat{\bf{u}}(k|k)$を取り出して$k$ステップ目での最適入力$u(k) = \Delta \hat{\bf{u}}(k|k) + \bf{u}(k-1)$としてシステムに送る．

$V(k)$は凸関数であるから，これを最小化するには$V$の微分を考えると，

$\begin{aligned}\nabla_{\Delta\mathcal{U}(k)}V = -\bf{g} + 2H\Delta\mathcal{U}(k)\end{aligned}$

を$0$にする$\Delta \mathcal{U}(k)$，すなわち凸点が最適入力である．制約が無い場合これは直ちに求まり，

$\begin{aligned}0 &= -\bf{g} + 2\bf{H}\Delta\mathcal{U}(k)\\\Delta\mathcal{U}(k) &= \frac{1}{2}H^{-1}\bf{g}\end{aligned}$

である．しかし制約がある場合制約の範囲内に凸点が入っていない場合がある．

ここで入出力に制約(constraint)がある場合を考える．入力に対する制約は$\bf{c}_1 \leq \bf{u}(k) \leq \bf{c}_2$で与えられ，出力には$\bf{c}_3 \leq \bf{x}(k) \leq \bf{c}_4$と与えられる．この制約をコスト関数$V(k)$の計算において考慮するが，そのために制約を$\Delta\mathcal{U}(k)$の関数に揃える必要がある．各入力，出力への制約は係数や定数項の関係から，行列$F$，$G$を用いて次のように書き直せる．

$\begin{aligned}F\left[\begin{array}{c}\hat{\bf{u}}(k|k)\\\vdots\\\hat{\bf{u}}(k+h_u-1|k)\\1\end{array}\right] &\leq \bf{0}\\G\left[\begin{array}{c}\hat{\bf{x}}(k+1|k)\\\vdots\\\hat{\bf{x}}(k+h_p|k)\\1\end{array}\right] &\leq \bf{0}\\\end{aligned}$

まず入力への制約を$\Delta\bf{u}(k)$の関数にする．制約の数を$q$，$\bf{u}$の次元を$m$とすると，$F$は次のような構造をもつ．

$\begin{aligned}F &= \left[\begin{array}{ccccc}F_1 & F_2 & \cdots & F_{h_u} & \bf{f}\end{array}\right]\\\end{aligned}$

$F_i$は$q\times m$行列であり，$\bf{f}$は$q\times 1$ベクトルである．これを用いて，

$\begin{aligned}F\left[\begin{array}{c}\hat{\bf{u}}(k|k)\\\vdots\\\hat{\bf{u}}(k+h_u-1|k)\\1\end{array}\right] &\leq \bf{0}\\\sum_{i=1}^{h_u}F_i\hat{\bf{u}}(k+i-1|k)+f &\leq \bf{0}\end{aligned}$

である．これに$\hat{\bf{u}}$の定義から，

$\begin{aligned}\hat{\bf{u}}(k+i-1|k) &= u(k-1) + \sum_{j=0}^{i-1}\Delta \hat{\bf{u}}(k+j|k)\end{aligned}$

の関係を用いると，

$\begin{aligned}\sum_{i=1}^{h_u}\sum_{j=i}^{h_u}F_j\Delta \hat{\bf{u}}(k+i-1|k)+\sum_{j=1}^{h_u}F_j\bf{u}(k-1) + \bf{f} \leq \bf{0}\end{aligned}$

と書き直せる．$F_i’ = \sum_{j=i}^{h_u}F_j$，$F’ = \left[\begin{array}{cccc}F_1′ & F_2′ & \cdots & F_{h_u}’\end{array}\right]$と定義すると，

$\begin{aligned}F’\left[\begin{array}{c}\Delta\hat{\bf{u}}(k|k)\\\vdots\\\Delta\hat{\bf{u}}(k+h_u-1|k)\end{array}\right] &\leq -F_1’\bf{u}(k-1) – \bf{f}\\F’\Delta\mathcal{U}(k) &\leq -F_1’\bf{u}(k-1) – \bf{f}\end{aligned}$

右辺は前ステップまでに求まる既知のベクトルである．以上より入力に対する制約を$\Delta\mathcal{U}(k)$の関数に変換することができる．

続いて出力に対する制約を$\Delta\mathcal{U}(k)$の関数に変換する．出力への制約の不等式は，予測出力ベクトルの関係式を用いて，

$\begin{aligned}G\left[\begin{array}{c}\Psi x(k) + \Upsilon\bf{u}(k-1) + \Theta\Delta\mathcal{U}(k)\\1\end{array}\right] &\leq \bf{0}\\\end{aligned}$

と書き直せる．ここで$G$の最終列をベクトル$\bf{\gamma}$，それ以外を行列$\Gamma$とみなすと，$G = \left[\begin{array}{cc}\Gamma & \bf{\gamma}\end{array}\right]$と分解できる．これを用いると，

$\begin{aligned}\left[\begin{array}{cc}\Gamma & \bf{\gamma}\end{array}\right]\left[\begin{array}{c}\Psi x(k) + \Upsilon\bf{u}(k-1) + \Theta\Delta\mathcal{U}(k)\\1\end{array}\right] &\leq \bf{0}\\\Gamma \left[\begin{array}{c}\Psi x(k) + \Upsilon\bf{u}(k-1)\end{array}\right] + \Gamma\Theta\Delta\mathcal{U}(k) + \bf{\gamma} &\leq \bf{0}\\\end{aligned}$

と書き直せる．以上より出力に対する制約も$\Delta\mathcal{U}(k)$の関数に変換することができた．

ここまで求めてきた$\Delta\mathcal{U}(k)$のコスト関数$V(k)$の最小化問題は制約付き二次計画問題(QP問題)となり，QPKWIK法などによる求値を行う．ここで$V(k)$を構成する$\bf{g}$は$h_u\times1$ベクトル，$\bf{H}$は$h_u\times h_u$行列であるから，制御ホライズンの大きさがQP問題に大きく関係すると言える．