Predictive Functional Controlの実装

目次

 Predictive Functional Controlとは?

まずはじめに,Predictive Functional Control (PFC)とは何者か?についてです.

PFCとはモデル予測制御(Model Predictive Control; MPC)のひとつです.Richaletによって提案されました.
Richaletは著書でPFCはPID制御に置き換わる制御則に成り得るのではないかと主張しています.

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

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

 

 

モデル予測制御の分類についてはMyEnigmaさんのところが詳しいのではないかと思います.

myenigma.hatenablog.com

 

また,実装は以下の論文を参考にしています.
(というより,こちらの先生にお世話になりました)

www.jstage.jst.go.jp

www.akita-pu.ac.jp

 

 

PFCがPIDに置き換わる理由(根拠)ですが,PFCが以下のような特徴を持っているためのようです.

  • 計算負荷が極めて軽量
  • チューニングが簡単
  • 入力制約を取り扱える
  • 高次の目標値に対してオフセットフリー特性を有する

実際のところ,3つ目と4つ目の項目については通常のMPCも持っている特徴です.

もともとMPCは化学プラントなどの比較的,実行周期の遅い系に対して適用されてきました.
MPCは最適化問題を周期ごとに解く必要があるため,応答の速いメカトロニクス系への適用は困難でした.
しかし近年のCPUなどのハードウェアの進歩や,最適化問題を解くためのアルゴリズムの発展により,自動運転車やロケット燃料の自動着陸などのメカトロニクス系にも適用されるようになってきました.

myenigma.hatenablog.com

これらは複数のアクチュエータが,それぞれ異なる応答性や入出力制約を有する複雑なシステムであり,目的とする応答を達成するためには,それぞれのアクチュエータにどれだけの操作を与えるべきか高度な制御を行う必要があります.
まさしくMPCの出番というべき制御対象です.

一方で,システムの構造自体は単純なものの,MPCの持つ入力制約の取り扱いや目標値へのオフセットフリー特性などを利用したい場合があります.

アクチュエーションの分野です.


 

アクチュエーションというと広いですが,ここではモータや油圧/空圧シリンダの位置サーボ制御や,ヒータによる温度制御などの低レイヤーな制御を指します.

これらの低レイヤーな制御では,未だに殆どがPIDで制御されているのではないでしょうか.

低レイヤーな制御では特にロバスト応答速度の速さ(周波数帯域の広さ),計算負荷の小ささが要求されます.

ここで,(私たちがよく遭遇し,たびたび頭を悩ませる)次のような制御系を考えます.

  • 比較的長いむだ時間を有する
  • 入力制約がある(たとえば定格電圧が10Vである,など)
  • 目標値が高次または高速であり位相遅れが生じる

 

f:id:blockahead:20180110222154p:plain

こういった場合,従来の(PIDでの)アプローチとしては以下のようなものがあります.

  • スミス補償器でむだ時間を補償する
  • 入力制約による積分器のワインドアップを,アンチワインドアップ(補償器)によって防止する
  • 零位相誤差追従制御(ZPETC)によって,追従による位相遅れを補償する

実際のところ,これらを個別で設計することは割と面倒です.
MPCではこれらを統一的な枠組みで考えることができ便利なため,アクチュエーションの分野にも適用したくなってきます.

先ほど挙げたスミス補償器やZPETCはモデル規範(Model-based)なアプローチであり,内部に制御対象の運動モデルを組み込むことで補償を実現しています.
いちばんの面倒事である運動モデルが必要ならMPCでいいじゃないか,とは思いませんか?

そこで,MPCの良さである制約の取り扱い / むだ時間の補償 / 高次目標値への高い追従性能は残しつつ,全体的に簡略化 / 軽量化したPredictive Functional Controlはどうか,ということになります.

概念

Predictive Functional Control(PFC)では,MPCの構造を大胆に簡略化しています.
PFCの理論面での特徴としては以下のものがあります.

  • 線形時不変システムのみを対象とする
  • 自由応答強制応答に分解する
  • 制御入力を基底函数のみで表現する
  • 最適化問題を解かない
  • 入力制約は扱えるが出力制約は考慮しない

線形時不変システムのみを対象とする

MPCの中には非線形のシステムを取り扱えるものもありますが,PFCでは線形時不変(Linear Time Invariant)システムのみを取り扱います.

f:id:blockahead:20180110224100p:plain

線形システムであるとは,ある制御対象G(k)に入力U_{1}(k)U_{2}(k)を加えたときのそれぞれの出力の間に,以下が成り立つことをいいます.

\begin{align}
G(k) \{ U_{1}(k)+U_{2}(k) \} = G(k)U_{1}(k)+G(k)U_{2}(k)
\end{align}

 

システム内に非線形現象が関わってくるとPFCでの制御対象外となります.
しかしほとんどの場合においてコントローラ内部のモデルは動作点周りで線形化するのが通常ですので,よほど特殊な系でなければ問題なく適用可能かと思います.

制御対象となる実際の制御対象モデル(プラントモデル)に対して,コントローラの内部に持たせる設計用のモデルのことをノミナルモデルと呼びます.
たとえば,線形二次レギュレータ(LQR)で使用するシステムのA行列,B行列はノミナルモデルにあたります.

また,時不変システムであるとは,以下のようにシステムの特性が時間によって変化しないことを指します.

\begin{align}
G(k+i)=G(k)
\end{align}

制御工学の基礎概念ではありますが,これがかなり重要となってきます.

 

自由応答と強制応答に分解する

PFCでは,システムの応答を自由応答強制応答に分解して考えます.
まず,自由応答とはその名の通り,システムに入力を与えないときのシステムの応答を指します.

f:id:blockahead:20180110235509p:plain

たとえば運転中に車のアクセルを離したとき,車は蓄えられていた運動エネルギによってしばらく進み続けます.

しばらくすると摩擦などの減衰によって車は停止します.
これが自由応答です.

一方で,アクセルを踏み込むと車は加速していきます.
これを強制応答と呼びます.


 

実際のところ,ある時刻 kでアクセルを踏んでも,それまでの速度(運動エネルギ)は保存されているはずです.
これらを分解して考えてもよいのでしょうか.

ここでシステムの線形性を考えると,アクセルを踏んでいるときの動きは以下のように書くことができます.

\begin{align}
車 の 応 答 &= 車 の 伝 達 函 数 \times ( ア ク セ ル + 運 動 エ ネ ル ギ ) \\
&= 車 の 伝 達 函 数 \times ア ク セ ル + 車 の 伝 達 函 数 \times 運 動 エ ネ ル ギ \\
&= 強 制 応 答 + 自 由 応 答
\end{align}

つまり強制応答とは,これまでの運動とは無関係な応答と捉えることができます.

これらを踏まえるとシステム全体の応答 y(k)は,自由応答 S_{L}(k)と強制応答 S_{F}(k)を利用して以下のように書くことができます.

\begin{align}
y(k) = S_{L}(k) + S_{F}(k)
\end{align}

制御入力を基底函数のみで表現する

MPCでは制御入力(操作量)は任意の値を持ちますが,PFCでは基底函数(Basis function)の重ね合わせで表現します.
基底函数とは,ステップ(0次) / ランプ(1次) / パラボラ(2次)といった函数です.

f:id:blockahead:20180111003159p:plain

制御対象は線形であると仮定しているため,ステップ入力 U_{0}(k),ランプ入力 U_{1}(k),パラボラ入力 U_{2}(k)をそれぞれ別で与えたときの強制応答と,入力をすべて足して与えたときの強制応答 S_{F}(k)の関係は以下のとおりです.

\begin{align}
S_{F}(k) &= G(k) \{ U_{0}(k) + U_{1}(k) + U_{2}(k) \} \\
&= G(k) U_{0}(k) + G(k) U_{1}(k) + G(k) U_{2}(k)
\end{align}

ただし G(k)は伝達函数です.

ここで,操作量 MVを以下の式で表します.

\begin{align}
MV(k+i) &= U_{0}(k+i) + U_{1}(k+i) + U_{2}(k+i) \\
&= \mu_{0} + \mu_{1} \cdot i+ \mu_{2} \cdot i^{2}
\end{align}

すると現在時刻 kから i秒だけ先の強制応答は以下のように書き換えることができます.

\begin{align}
S_{F}(k+i) &= G(k+i) MV(k+i) \\
&= G(k) \{ \mu_{0} + \mu_{1} \cdot i + \mu_{2} \cdot i^{2} \} \\
&= G_{0}(k,i) \cdot \mu_{0} + G_{1}(k,i) \cdot \mu_{1} + G_{2}(k,i) \cdot \mu_{2}
\end{align}

結局,強制応答は現在の出力とは無関係になるため時刻 kに依存せず,以下のようになります.

\begin{align}
S_{F}(i) = G_{0}(i) \cdot \mu_{0} + G_{1}(i) \cdot \mu_{1} + G_{2}(i) \cdot \mu_{2}
\end{align}

ここで, G_{0}(i) G_{1}(i) G_{2}(i)はそれぞれ,ステップ入力,ランプ入力,パラボラ入力を与えたときのシステムの強制応答です.

最適化問題を解かない

通常,MPCでは線形時不変システムであるなしに関わらず最適化問題をサンプリング時刻ごとに解く必要があります.

これがMPCの計算負荷が大きいと言われる原因です.

一方でPFCではサンプリング時刻ごとに行列の掛け算をするだけです.
具体的には,最小二乗法の計算を1度行います.

PFCでは参照軌道(Reference trajectory)を利用して,ある時刻 h秒だけ先のタイミングでのみ最適化を評価します.

この評価する点のことを一致点(Coincidence point)と呼びます.

 

f:id:blockahead:20180113221251p:plain

参照軌道は図に示すように,現在の出力 CV(k)から目標値 SV(k)まで指数関数的に近づくような軌道を使用することが多いです.

\begin{align}
参照軌道 = [ SV(k) - CV(k) ] \cdot l_{h}(i)
\end{align}

 l_{h}(i)=1-e^{-ai}のような指数関数的に1に漸近する函数を使用します.

一致点の数は基底函数の次数以上必要です.

多くの場合,パラボラ状目標値まで追従できるように基底函数を3次にとるため,一致点の数も3つに取ることが多いようです.

たとえサンプリング周期が1msで評価区間が1secだとしても,参照軌道の3点しか評価しないため計算負荷が小さくなるのです.

たかだか3x3の行列の掛け算を解くことで充分なのです.

普通のMPCでは評価区間全体に渡って評価を行うため,上の条件では1000点を評価することが必要になります.

入力制約は扱えるが出力制約は考慮しない

制約の取り扱いはモデル予測制御系の制御則の大きなメリットのひとつです.
PFCでは出力制約は取り扱わず,入力制約のみを取り扱います.

MPCで計算が複雑になる原因として出力制約の影響がそれなりに大きいと思っています.

さらに,PFCでは入力制約も計算時に制約として与えるわけではなく,計算の結果得られた操作量に対してリミットを掛けることで制約として扱います.

f:id:blockahead:20180113221330p:plain

ただし,リミットをかけた後の操作量をコントローラ側にフィードバックすることで,入力制約によるワインドアップを防止しています.

制約の取り扱いについて"こんなんでいいのか"といった感じにはなるのですが,そもそも低レイヤーな制御では制約といっても入力側にアンチワインドアップ補償器が入っているだけでどうにかなってきたのですから今更かなあという思いもあります.

出力制約については,あまり出力制約が有効なアプリケーションを知らないということもあり,ここでは触れていません.
"こういうところで出力制約が有効だよ!"というのがありましたら教えていただけるとありがたいです.

実装

さて,もう少し細かい話に入っていきます.
制御対象のノミナルモデルを以下の離散時間の線形時不変システムとします.
簡単のため1入力1出力(SISO)系として考えますが,多自由度でも拡張して利用できます.

\begin{align}
x(k+1) &= Ax(k)+Bu(k) \\
y(k) &= Cx(k)
\end{align}

このとき, iサンプル先の状態は以下のようになります.

\begin{align}
x(k+1) &= Ax(k) + Bu(k) \\
x(k+2) &= Ax(k+1) + Bu(k+1) \\
&= A^{2}x(k) + ABu(k) + Bu(k+1) \\
x(k+3) &= Ax(k+2) + Bu(k+2) \\
&= A^{3}x(k) + A^{2}Bu(k) + ABu(k+1) + Bu(k+2) \\
& \vdots \\
x(k+i) &= A^{i}x(k) + [ A^{i-1}B, A^{i-2}B, \cdots, A^{2}B, AB, B] \cdot \left[ \begin{array}{c}
u(k) \\
u(k+1) \\
\vdots \\
u(k+i-3) \\
u(k+i-2) \\
u(k+i-1)
\end{array}\right]
\end{align}

同様に, iサンプル先の出力は以下のとおりです.

\begin{align}
y(k+i) &= Cx(k+i) \\
&= CA^{i}x(k) + [ CA^{i-1}B, CA^{i-2}B, \cdots, CA^{2}B, CAB, CB] \cdot \left[ \begin{array}{c}
u(k) \\
u(k+1) \\
\vdots \\
u(k+i-3) \\
u(k+i-2) \\
u(k+i-1)
\end{array}\right]
\end{align}

自由応答

自由応答 S_{L}(k+i) u(k)=\cdots=u(k+i-1)=0のときの応答であるため,以下のように得られます.

\begin{align}
S_{L}(k+i) = CA^{i}x(k)
\end{align}

一致点を h_{1} h_{2} h_{3}の3点とすると,以下のようなベクトルが得られます.

\begin{align}
S_{L}(k, h) = \left[ \begin{array}{c}
CA^{h_{1}}x(k) \\
CA^{h_{2}}x(k) \\
CA^{h_{3}}x(k)
\end{array} \right]
\end{align}

強制応答

強制応答は単位ステップ入力を考えたとき,以下の式で得られます.

\begin{align}
G_{0}(i) = [ CA^{i-1}B, CA^{i-2}B, \cdots, CA^{2}B, CAB, CB ] \cdot \left[ \begin{array}{c}
1 \\
1 \\
1 \\
\vdots
\end{array} \right]
\end{align}

単位ランプ入力の場合は以下のとおりです.

\begin{align}
G_{1}(i) = [ CA^{i-1}B, CA^{i-2}B, \cdots, CA^{2}B, CAB, CB ] \cdot \left[ \begin{array}{c}
0 \\
T_{s} \\
2T_{s} \\
\vdots
\end{array} \right]
\end{align}

パラボラは省略しますが,全体の強制応答 S_{F}(i)は以下によって得られます.

\begin{align}
S_{F}(i) = G_{0}(i) \cdot \mu_{0} + G_{1}(i) \cdot \mu_{1} + G_{2}(i) \cdot \mu_{2}
\end{align}

一致点を h_{1} h_{2} h_{3}の3点とすると,以下のような行列が得られます.

\begin{align}
S_{F}(h) = \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}

システム全体の応答

さて,システム全体の応答は以下のように得られます.

\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}}x(k) \\
CA^{h_{2}}x(k) \\
CA^{h_{3}}x(k)
\end{array} \right] + \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}

 y(k+i)を現在時刻からの増分 \Delta yを使って表すと,以下の式が得られます.

\begin{align}
\left[ \begin{array}{c}
\Delta y(k, h_{1}) \\
\Delta y(k, h_{2}) \\
\Delta y(k, h_{3})
\end{array} \right] = \left[ \begin{array}{c}
CA^{h_{1}}x(k) \\
CA^{h_{2}}x(k) \\
CA^{h_{3}}x(k)
\end{array} \right] + \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] - \left[ \begin{array}{c}
y(k) \\
y(k) \\
y(k)
\end{array} \right]
\end{align}

参照軌道

一方で参照軌道は,以下で示されます.

\begin{align}
\left[ \begin{array}{c}
Ref(k, h_{1}) \\
Ref(k, h_{1}) \\
Ref(k, h_{1})
\end{array} \right] = \left[ \begin{array}{c}
[ SV(k) - CV(k)] \cdot l_{h}(h_{1}) \\
[ SV(k) - CV(k)] \cdot l_{h}(h_{2}) \\
[ SV(k) - CV(k)] \cdot l_{h}(h_{3})
\end{array} \right]
\end{align}

システムの予測される応答 \Delta yと目標とする軌道 Refが一致していれば良いため,最終的に以下の等式が得られます.

\begin{align}
& \left[ \begin{array}{c}
[ SV(k) - CV(k)] \cdot l_{h}(h_{1}) \\
[ SV(k) - CV(k)] \cdot l_{h}(h_{2}) \\
[ SV(k) - CV(k)] \cdot l_{h}(h_{3})
\end{array} \right] \\
&= \left[ \begin{array}{c}
CA^{h_{1}}x(k) \\
CA^{h_{2}}x(k) \\
CA^{h_{3}}x(k)
\end{array} \right] + \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] - \left[ \begin{array}{c}
y(k) \\
y(k) \\
y(k)
\end{array} \right]
\end{align}

これを操作量 \muについて解くと,以下の式が得られます.

\begin{align}
\left[ \begin{array}{c}
\mu_{0} \\
\mu_{1} \\
\mu_{2}
\end{array} \right] &= \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]^{-1} \\
& \left\{ \left[ \begin{array}{c}
[ SV(k) - CV(k)] \cdot l_{h}(h_{1}) \\
[ SV(k) - CV(k)] \cdot l_{h}(h_{2}) \\
[ SV(k) - CV(k)] \cdot l_{h}(h_{3})
\end{array} \right] - \left[ \begin{array}{c}
CA^{h_{1}}x(k) \\
CA^{h_{2}}x(k) \\
CA^{h_{3}}x(k)
\end{array} \right] + \left[ \begin{array}{c}
y(k) \\
y(k) \\
y(k)
\end{array} \right] \right\}
\end{align}

長々と計算しましたが,この式がPredictive Functional Controlの全てです.

実際のところ,ステップ入力である \mu_{0}以外は開始時刻にゼロであるため,操作量 MV(k) = \mu_{0}となります.

上式の右辺のうち未知な項は,目標値 SV(k),現在値 CV(k),プラントモデルの状態 x(k),プラントモデルの出力 y(k)のみです.
それ以外についてはオフラインで計算することができます.

まず, SV CVについては問題なく得ることができます.

次にプラントモデルの状態と出力についてですが,これはサンプリング周期ごとにノミナルモデルに前回の操作量 MV(k−1)を与えることで得られます.

単純に以下の操作を行います.

\begin{align}
x(k) &= Ax(k - 1) + B \cdot MV(k - 1) \\
y(k) &= Cx(k)
\end{align}


 

いろいろ長くなりましたが,PFCではサンプリング周期ごとに以下の操作を行うだけです.

  1. 内部モデル x(k) y(k)の更新をする
     x(k) = Ax(k-1) + B \cdot MV(k-1)
  2. 現在時刻での目標値 SV(k)と現在値 CV(k)を取得する.
  3.  \muを求める式に x(k) y(k) SV(k) CV(k)を代入する
  4. 求めた \muのうち,ステップ応答の項 \mu_{0}のみをシステムに入力する

制御対象は2次遅れ系とし,ノミナルモデルと制御対象にモデル化誤差は無いものとしています.
ゲイン 1,固有周波数 1 (rad/s),減衰係数 0.5のシステムです.

\begin{align}
G(s) = \frac{1}{s^2 + s + 1 }
\end{align}

f:id:blockahead:20180113233410p:plain

いちばん左は基底函数と一致点の関係についてです.
青色で示す点が G_{0}などに相当します.

真ん中は通常のPFCでの結果です.
目標値,外乱に対してもちゃんと制御できていることがわかります.
ただ,正弦波入力では遅れが生じています.

右側は位相進み補償を入れたPFCでの結果です.
正弦波入力に対しても遅れなく追従できています.


 

この結果はSimulink project形式でGithubにアップロードしています.

github.com