多入力多出力システムのPredictive Functional Control

はじめに

以前，Predictive Functional Control（PFC）についての記事を書きましたが，
あくまで単入力単出力（Single input single output; SISO）システムについてのみでした．

これを多入力多出力（Multi input multi output; MIMO）システムに拡張することが今回の目的です．

hamachannel.hatenablog.com

PFCの基礎に関する文献は考案者であるJ. Richaletの文献が詳しいです．

翻訳書もあったようなのですが，現在では絶版しています．

Predictive Functional Control: Principles and Industrial Applications (Advances in Industrial Control)

作者: Jacques Richalet,Donal O'Donovan,Karl E.Åstroem
出版社/メーカー: Springer
発売日: 2009/05/25
メディア: ハードカバー
この商品を含むブログを見る

基本的な概念は通常のモデル予測制御（Model predictive control; MPC）と同じですので，MPCを知っているとMIMOシステムへの拡張は容易です．

MIMOシステムのPFC

SISOシステムでのPFC

以前の記事で述べたとおり，PFCではシステムの応答 $y$ を自由応答と強制応答に分解します．

右辺第1項が自由応答，右辺第2項が強制応答です．

\begin{align}
\left[ \begin{array}{c}
y(k+h_{1}) \\
y(k+h_{2}) \\
y(k+h_{3})
\end{array} \right] = \left[ \begin{array}{c}
CA^{h_{1}} \\
CA^{h_{2}} \\
CA^{h_{3}}
\end{array} \right] x(k) + \left[ \begin{array}{ccc}
G_{0}(h_{1}) & G_{1}(h_{1}) & G_{2}(h_{1}) \\
G_{0}(h_{2}) & G_{1}(h_{2}) & G_{2}(h_{2}) \\
G_{0}(h_{3}) & G_{1}(h_{3}) & G_{2}(h_{3})
\end{array} \right] \left[ \begin{array}{c}
\mu_{0} \\
\mu_{1} \\
\mu_{2}
\end{array} \right]
\end{align}

ここで， $A$ , $C$ はシステムの状態方程式表現での係数行列， $h_{n}$ は予測ホライズン上の一致点（参照点）， $\mu_{k}$ は基底函数（Basis function）で操作量を表します．

Basis functionにおいて， $\mu_{0}$ ， $\mu_{1}$ ， $\mu_{2}$ はそれぞれステップ，ランプ，パラボラ入力を表します．

また $G_{k}(h_{n})$ は，操作量 $\mu_{k}$ を与えたときの一致点 $h_{n}$ 時点でのプラント出力を表します．

たとえば， $h = [ 10, 20, 30 ]^{T}$ とすると， $G_{0}(h_{2})$ は「ゼロ状態からステップ入力を与えた際の20サンプル先のシステム出力」を表します．

PFCの文献の多くはシステムを伝達函数で表現していますが，このような状態方程式ベースで書いておくと拡張は容易です．

MIMOシステムへの拡張

MIMOシステムへの拡張は簡単で，2入力2出力システムの場合は自由応答と強制応答を以下のように（文字通り）拡張するだけです．

\begin{align}
\left[ \begin{array}{c}
y_{1}(k+h_{1}) \\
y_{2}(k+h_{1}) \\
y_{1}(k+h_{2}) \\
y_{2}(k+h_{2}) \\
y_{1}(k+h_{3}) \\
y_{2}(k+h_{3})
\end{array} \right] &= \left[ \begin{array}{c}
CA^{h_{1}} \\
CA^{h_{2}} \\
CA^{h_{3}}
\end{array} \right] x(k) \\
& + \left[ \begin{array}{cccccc}
G_{110}(h_{1}) & G_{120}(h_{1}) & G_{111}(h_{1}) & G_{121}(h_{1}) & G_{112}(h_{1}) & G_{122}(h_{1}) \\
G_{210}(h_{1}) & G_{220}(h_{1}) & G_{211}(h_{1}) & G_{221}(h_{1}) & G_{212}(h_{1}) & G_{222}(h_{1}) \\
G_{110}(h_{2}) & G_{120}(h_{2}) & G_{111}(h_{2}) & G_{121}(h_{2}) & G_{112}(h_{2}) & G_{122}(h_{2}) \\
G_{210}(h_{2}) & G_{220}(h_{2}) & G_{211}(h_{2}) & G_{221}(h_{2}) & G_{212}(h_{2}) & G_{222}(h_{2}) \\
G_{110}(h_{3}) & G_{120}(h_{3}) & G_{111}(h_{3}) & G_{121}(h_{3}) & G_{112}(h_{3}) & G_{122}(h_{3}) \\
G_{210}(h_{3}) & G_{220}(h_{3}) & G_{211}(h_{3}) & G_{221}(h_{3}) & G_{212}(h_{3}) & G_{222}(h_{3}) \\
\end{array} \right] \left[ \begin{array}{c}
\mu_{10} \\
\mu_{20} \\
\mu_{11} \\
\mu_{21} \\
\mu_{12} \\
\mu_{22}
\end{array} \right]
\end{align}

ここで $y_{i}$ はシステムの $i$ 番目の出力， $\mu_{jk}$ はシステムの $j$ 番目の操作量を表します．

たとえば先で述べた一致点 $h$ を利用すると， $G_{120}(h_{2})$ は「ゼロ状態からシステムの2番めの入力としてステップ入力を与えた際の，20サンプル先でのシステムの1番目の出力」となります．（いわゆる干渉項です）

より単純な例として，一致点 $h$ を1点，basis functionをステップ $\mu_{j0}$ のみとした場合には以下のような式が得られます．
左上の $2 \times 2$ 要素のみが残るようなイメージです．

\begin{align}
\left[ \begin{array}{c}
y_{1}(k+h) \\
y_{2}(k+h) \\
\end{array} \right] = \left[ \begin{array}{c}
CA^{h} \\
\end{array} \right] x(k) + \left[ \begin{array}{cc}
G_{110}(h) & G_{120}(h) \\
G_{210}(h) & G_{220}(h)
\end{array} \right] \left[ \begin{array}{c}
\mu_{10} \\
\mu_{20} \\
\end{array} \right]
\end{align}

Basis functionの次数を1にしたときの2入力2出力システムの強制応答 $G_{ijk}(h_{n})$ は下図のようになります．

サンプリング周期 $Ts = 1$ s，一致点は $h = [ 10, 15, 30 ]^{T}$ の3点を想定しています．

f:id:blockahead:20180728192443p:plain

実際には一致点は立ち上がりの範囲で取ることが多く，この例だと最大でも $h = [ 2, 3, 6 ]^{T}$ くらいかと思います．

これらの強制応答は，シミュレーションによって得ることもできますし，以下のように畳み込みで解析的に求めることもできます．

\begin{align}
G_{*20}(h_{2}) =\left[ \begin{array}{cccccc}
CA^{h_{2}-1}B & CA^{h_{2}-2}B & \cdots & CA^{2}B & CAB & CB
\end{array} \right] \left[ \begin{array}{c}
\left[ \begin{array}{c}
0 \\
1
\end{array} \right] \\
\left[ \begin{array}{c}
0 \\
1
\end{array} \right] \\
\vdots \\
\left[ \begin{array}{c}
0 \\
1
\end{array} \right]
\end{array} \right]
\end{align}

通常のMPCでは右辺の横長の行列をさらに縦に並べたToeplitz行列を保持しますが，PFCでは畳み込まれたあとの結果だけでよいため，メモリの使用量と計算負荷が小さいという特徴があります．

シミュレーション（MATLAB/Simulink）

対象とするシステム

ここで対象とするシステムは，以下のような単純な2自由度のSpring-mass-damperシステムです．

両質量の位置が測定できるものとしています．

f:id:blockahead:20180729101832p:plain

状態方程式

\begin{align}
\frac{d}{dt} \left[ \begin{array}{c}
x_{1} \\
x_{2} \\
\dot{x_{1}} \\
\dot{x_{2}}
\end{array} \right] = \left[ \begin{array}{cccc}
0 & 0 & 1 & 0 \\
0 & 0 & 0 & 1 \\
-\frac{K_{1}+K_{2}}{M_{1}} & \frac{K_{2}}{M_{1}} & -\frac{D_{1}+D_{2}}{M_{1}} & \frac{D_{2}}{M_{1}} \\
\frac{K_{2}}{M_{2}} & -\frac{K_{2}}{M_{2}} & \frac{D_{2}}{M_{2}} & -\frac{D_{2}}{M_{2}}
\end{array} \right] \left[ \begin{array}{c}
x_{1} \\
x_{2} \\
\dot{x_{1}} \\
\dot{x_{2}}
\end{array} \right] + \left[ \begin{array}{cc}
0 & 0 \\
0 & 0 \\
\frac{1}{M_{1}} & 0 \\
0 & \frac{1}{M_{2}}
\end{array} \right] \left[ \begin{array}{c}
f_{1} \\
f_{2}
\end{array} \right]
\end{align}

観測方程式

\begin{align}
\left[ \begin{array}{c}
x_{1} \\
x_{2}
\end{array} \right] = \left[ \begin{array}{cccc}
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0
\end{array} \right] \left[ \begin{array}{c}
x_{1} \\
x_{2} \\
\dot{x_{1}} \\
\dot{x_{2}}
\end{array} \right]
\end{align}

ここでは， $M$ ， $D$ ， $K$ などすべての値は $1$ に設定しています．