適応制御（セルフチューニングレギュレータ）

はじめに

今、適応制御について関心を持つ人がどれほどいるのか分かりませんが、、、

これまで、このブログでは、システムの振る舞いを状態方程式で記述し、その時間応答の確認や制御ロジックの適用について書いてきたものの、システムを物理式で記述すること、さらにパラメータを特定するようなモデル化は難しい という話はよくあります。

そうですよねぇ。

特にマスプロダクトでは、物バラつきや経年劣化、使われ方や使われる環境も様々で、物理式による完全なモデル化は難しいですよねぇ。

これが、様々ある制御ロジックの実用化を難しくしている理由ではないかと思います。

でも、PIDじゃ性能が足りないんだよなぁ という場合、適応制御はいかがでしょうか。

適応制御は、オンラインでシステム同定しながら、同定したモデルを使ってコントローラを構成し制御を実現します。適応制御には多くの形態がありますが、このブログではセルフチューニングレギュレータ（STR）について書いていきます。

Excelファイルのダウンロード

セルフチューニングレギュレータ（STR）のExcelファイルはこちら → リンク からダウンロードできます。適応制御(セルフチューニングレギュレータ).xlsm です。

Excelファイルの内容については順次説明していきます。

文献紹介

まず最初に、セルフチューニングレギュレータ（STR）の実用例として文献を紹介します。「適応制御理論を用いたA/Fフィードバック制御システム」です。
下記のリンクからアクセスできます。

適応制御理論を用いたA/Fフィードバック制御システム | 本田技術研究所 論文サイト (hondarandd.jp)

右側のPDFボタンを押して登録情報を入力すると論文の詳細を確認できます。

制御ロジックの概要

このブログでは、紹介した文献にあるように、セルフチューニングレギュレータ（STR）でフィードバック制御を構成してみます。

このブログで扱う仮想的な対象システム（プラント）のブロック図はこんな感じです。

プラント制御は世の中に様々なものがあると思いますが、ここではプラントに投入する燃料や反応材料の量、流量や温度、圧力などの物理量を目標に応じて制御してプラントを運転する制御システムをイメージしています。

プラント各部の制御量の名称と特性

Mobj

プラントに投入する燃料や反応材料の量、または、流量や温度、圧力などプラント制御の目標物理量

Mctr

システムの制御指示量。アクチュエータは制御指示量に応じて作動するものとします。

Mact

実際にプラントに反映できた物理量。
今回の仮想的なシステムでは、アクチュエータへの制御指示量（Mctr）と実際にプラントに反映できる物理量（Mact）の間に動特性があるものとします。また、実際にプラントに反映できた物理量（Mact）自体は検出できないものとします。

Mactの動特性について

システムの制御指示量（Mctr）と実際にプラントに反映できる物理量（Mact）の間の動特性は一次遅れのモデルとし、以下のように Mact と Mctr の時系列データのコンボリューションで算出します。 $$M_{act}(n) = a1 * M_{act}(n-1) + b1 * M_{ctr}(n)$$ また、Excelファイルでは、この動特性はスライドボリュームで時定数0.8秒から1.2秒の間で随時変更可能になっています。

Kstr

STRフィードバック補正係数です。

Ract

制御目標物理量に対し実際にプラントに反映できた物理量の比率。 $$R_{act} = M_{act}/M_{obj}$$ 制御目標物理量通りプラントに反映できれば $R_{act}=1.0$。10%不足していれば $R_{act}=0.9$ になります。また、Ractはプラントでの反応結果などからセンサで検出できるものとします。反応結果であるRactでシステムを制御する構成なのでフィードバック制御になります。

Robj

制御目標物理量と実際にプラントに反映できた物理量の比率の目標値で、常に 1.0 （Mobj = Mact）です。

システム同定アルゴリズム

適応制御では、対象となるシステムの動特性を制御サイクルごとに逐次同定します。
システム同定については、こちら → リンク で説明しています。

同定アルゴリズムも様々ありますが、このブログでは システムの同定 でも使った最も基本的なアルゴリズムである逐次最小二乗推定を使います。

逐次最小二乗推定の計算ロジックは、

ARMAモデル

$$y(n)=\sum_{k=1}^{p} a(k) y(n-k)+\sum_{k=1}^{q} b(k) u(n-k)$$ において

$$x=\left[a_{1}, a_{2}, \cdots \cdot, a_{p}, b_{1}, b_{2}, \cdots \cdot, b_{q}\right]^{T}$$ $$z(n)=[y(n-1), \cdots, y(n-p), u(n-1), \cdots, u(n-q)]^{T}$$ $$Z=\begin{bmatrix} z^{T}(n+1) \\ \vdots \\ z^{T}(N) \\  \end{bmatrix}$$ $$P_N=({Z_N}^{T} Z_N)^{-1}$$ としたときに、最小二乗法で逐次的に 𝑥 を推定するアルゴリズムは $$\hat{x}_{n+1}=\hat{x}_{n}+g(n+1)\left\{\mathrm{y}(n+1)-z^{T}(n+1) \hat{x}_{n}\right\}$$ $$g(n+1)=\frac{P_{n} z(n+1)}{1+z^{T}(n+1) P_{n} z(n+1)}$$ $$P_{n}=\left[I-g(n) z^{T}(n)\right] P_{n-1}$$ です。

詳細については 参考図書 などを参考にしていただきたいです。変更点としては、安定性に寄与する忘却係数を追加しています。

逐次最小二乗推定の計算ロジックは、ExcelファイルにVBAで実装しています。プログラムは「開発」タグを選択 →「Visual Basic」の Module1 です。

VBAのコードはこんな感じです。

逐次最小二乗推定

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356

'///////////////////////////////////////////////// '// '//　逐次最小二乗推定 '// '//   X=[a1,a2,...,ak,b1,b2,...,bk]T '//   Z(n)=[y(n-1),...,y(n-k),u(n-1),...,u(n-k)]T '//   Y(n)=∑a(k)y(n-k)+∑b(k)u(n-k) '//   Order=2k '/////////////////////////////////////////////////   Dim U(6) As Double Dim Y(6) As Double Dim X(10) As Double Dim Z(10)   As Double     ' Zt:転置行列 Dim Znext(10)  As Double Dim P(10, 10) As Double Dim Pnext(10, 10) As Double Dim G(10)    As Double Dim KPU(10, 10) As Double ' Pの更新行列 Dim PZ(10)   As Double Dim Yest As Double        ' Y推定値 Dim Ynext As Double   Dim Order As Integer      ' ARMAモデルのパラメータ数(2k:a1,a2,..,ak,b1,b2,..,bk) Dim DATASIZE As Integer   ' データサイズ Dim NSIM As Integer       ' 演算回数 Dim IdCal As Boolean   Const OrderMAX As Integer = 10       ' ARMAモデルのパラメータ最大数 Const InitialValue  As Double = 500  ' P初期値 Const Frgt As Double = 0.97          ' 忘却係数     Sub LeastSquaresEstimation()       Dim I As Integer, J As Integer, N As Integer       Order = Range("C4").Value       ' ARMAモデルのパラメータ数(2k:a1,a2,..,ak,b1,b2,..,bk)     DATASIZE = Range("C5").Value    ' テストデータ数     IdCal = Range("C6").Value       ' システム同定の可否           NSIM = DATASIZE - Order - 1     ' シミュレーション回数       '/////////////////////////////////////////////////     '  初期化     '/////////////////////////////////////////////////     For I = 1 To OrderMAX         X(I) = 0     Next I     For I = 1 To Order         For J = 1 To Order             P(I, J) = 0             If I = J Then                 '// Pnextの初期値はσI（σは十分大きな正数：InitialValue 設定）                 Pnext(I, J) = InitialValue             Else                 Pnext(I, J) = 0             End If         Next J     Next I           '/////////////////////////////////////////////////     '   パラメータ更新ステータスを初期化     '/////////////////////////////////////////////////     Application.Calculation = xlManual    ' Excelシートの更新中止     For I = 1 To DATASIZE         Cells(9 + I, 19).Value = 0        ' 初期化     Next I     Application.Calculation = xlAutomatic ' Excelシートの更新許可           '/////////////////////////////////////////////////     '   パラメータ表示セルからA1,A2,...,B1,B2,... の初期値をX()に読み込む     '/////////////////////////////////////////////////     For I = 1 To (Order / 2)         X(I) = Cells(10, 8 + I).Value       ' a1,a2,...         X(I + (Order / 2)) = Cells(10, 13 + I).Value    ' b1,b2,...     Next I       '/////////////////////////////////////////////////     ' 入出力データ(U/Y)のシートからの読み込み     '  ARMAパラメータの次数をNとすると、テストデータは(N+1)からが信頼できるデータ     '  になるため、逐次推定の初回データの設定は N+1(=(Order/2)+1)～2N(=Order)の範囲     '/////////////////////////////////////////////////     For I = 1 To (Order / 2)         U(I) = Cells(9 + (Order / 2) + I, 3).Value    ' Kstr         Y(I) = Cells(9 + (Order / 2) + I, 6).Value    ' Ract     Next I       U((Order / 2) + 1) = Cells(9 + Order + 1, 3).Value    ' Kstr     Y((Order / 2) + 1) = Cells(9 + Order + 1, 6).Value    ' Ract       Ynext = Y((Order / 2) + 1)       '//////////////////////////////////////////////////     ' 入出力データ初期値のベクトル Z への設定     '//////////////////////////////////////////////////     For I = 1 To (Order / 2)         Znext(I) = Y((Order / 2) + 1 - I)         Znext(I + (Order / 2)) = U((Order / 2) + 1 - I)     Next I       '/////////////////////////////////////////////////     ' G(1)の算出     '    G(1) = P(0)Z(1)/(1+Zt(1)P(0)Z(1))     '/////////////////////////////////////////////////     Call CalGMatrix       '//////////////////////////////////////////////////     '  ARMAパラメータの逐次更新処理     '   X(1) = X(0)+G(1)(Y(1)-Zt(1)X(0))     '//////////////////////////////////////////////////     Call RenewXMatrix       '/////////////////////////////////////////////////     '   演算結果 X[a1,a2,...,b1,b2,...] をセルに反映し、     '   EXCEL表計算によるSTR制御演算を回し表示を更新する     '/////////////////////////////////////////////////     Call UpdateExcelSheet(0)         '/////////////////////////////////////////////////     '/////////////////////////////////////////////////     '   逐次処理     '/////////////////////////////////////////////////     '/////////////////////////////////////////////////     For N = 1 To NSIM           '/////////////////////////////////////////////////         '  入出力データ(U/Y)の読み込み         '/////////////////////////////////////////////////         For I = 1 To (Order / 2)             U(I) = U(I + 1)             Y(I) = Y(I + 1)         Next I           U((Order / 2) + 1) = Cells(9 + N + Order + 1, 3).Value    ' Kstr         Y((Order / 2) + 1) = Cells(9 + N + Order + 1, 6).Value    ' Ract                   Ynext = Y((Order / 2) + 1)       ' Y(N+1)           '//////////////////////////////////////////////////         ' 入出力データのベクトル Z への設定         '//////////////////////////////////////////////////         For I = 1 To Order             Z(I) = Znext(I)              ' Z(N)             For J = 1 To Order                 P(I, J) = Pnext(I, J)    ' P(N-1)             Next J         Next I           For I = 1 To (Order / 2)             Znext(I) = Y((Order / 2) + 1 - I)             Znext(I + (Order / 2)) = U((Order / 2) + 1 - I)         Next I           '/////////////////////////////////////////////////         ' P(N)の更新         '    P(N) = [I-G(N)Zt(N)]P(N-1)         '/////////////////////////////////////////////////         Call UpdatePMatrix           '/////////////////////////////////////////////////         ' G(N+1)の算出         '    G(N+1) = P(N)Z(N+1)/(1+Zt(N+1)P(N)Z(N+1))         '/////////////////////////////////////////////////         Call CalGMatrix           '//////////////////////////////////////////////////         '  ARMAパラメータの逐次更新処理         '   X(N+1) = X(N)+G(N+1)(Y(N+1)-Zt(N+1)X(N))         '//////////////////////////////////////////////////         Call RenewXMatrix           '/////////////////////////////////////////////////         '   演算結果 X[a1,a2,...,b1,b2,...] をセルに反映し、         '   EXCEL表計算によるSTR制御演算を回し表示を更新する         '/////////////////////////////////////////////////         Call UpdateExcelSheet(N)       Next N   End Sub     '///////////////////////////////////////////////// ' ' P(N)の更新 '    P(N) = [I-G(N)Zt(N)]P(N-1) ' '/////////////////////////////////////////////////   Sub UpdatePMatrix()       Dim I As Integer, J As Integer, K As Integer     Dim Sum As Double       For I = 1 To Order         For J = 1 To Order             '/////////////////////////////////////////////////             ' KPU = I - G(k)Zt(N)             '/////////////////////////////////////////////////             If I = J Then                 KPU(I, J) = 1 - (G(I) * Z(J))        ' KPU = 1-G(N)Zt(N)             Else                 KPU(I, J) = 0 - (G(I) * Z(J))        ' KPU = 0-G(N)Zt(N)             End If         Next J     Next I           '/////////////////////////////////////////////////     '  Pnext = KPU * P(N-1)     '/////////////////////////////////////////////////     For I = 1 To Order         For J = 1 To Order             Sum = 0             For K = 1 To Order                 Sum = Sum + KPU(I, K) * P(K, J)             Next K             Pnext(I, J) = Sum / Frgt    ' Frgt:忘却係数        Next J     Next I   End Sub     '///////////////////////////////////////////////// ' ' G(N+1)の算出 '    G(N+1) = P(N)Z(N+1)/(1+Zt(N+1)P(N)Z(N+1)) ' '/////////////////////////////////////////////////   Sub CalGMatrix()       Dim I As Integer, J As Integer     Dim Sum As Double       '/////////////////////////////////////////////////     ' P(N)Z(N+1) の算出     '/////////////////////////////////////////////////     For I = 1 To Order         Sum = 0         For J = 1 To Order             Sum = Sum + Pnext(I, J) * Znext(J)    ' PZ = P(N)Z(N+1)         Next J         PZ(I) = Sum     Next I       '/////////////////////////////////////////////////     ' 1+Zt(N+1)P(N)Z(N+1) の算出     '/////////////////////////////////////////////////     Sum = Frgt     ' Frgt:忘却係数  Sum = 1+Zt(N+1)P(N)Z(N+1)           For I = 1 To Order         Sum = Sum + Znext(I) * PZ(I)    ' Zt(N+1)P(N)Z(N+1)     Next I           '/////////////////////////////////////////////////     ' G(N+1) = P(N)Z(N+1)/(1+Zt(N+1)P(N)Z(N+1))     '/////////////////////////////////////////////////     For I = 1 To Order         G(I) = PZ(I) / Sum              ' G = P(N)Z(N+1)/(1+Zt(N+1)P(N)Z(N+1))     Next I    End Sub     '////////////////////////////////////////////////// ' '  ARMAパラメータの逐次更新処理 '   X(N+1) = X(N)+G(N+1)(Y(N+1)-Zt(N+1)X(N)) ' '//////////////////////////////////////////////////   Sub RenewXMatrix()       Dim Err As Double, Sum As Double     Dim I As Integer       '//////////////////////////////////////////////////     ' 推定出力 Yest (Zt(N+1)X(N)) 算出     '//////////////////////////////////////////////////     Sum = 0     For I = 1 To Order         Sum = Sum + Znext(I) * X(I)     Next I     Yest = Sum               '//////////////////////////////////////////////////     ' 推定誤差 Err (Y(N+1)-Zt(N+1)X(N)) 算出     '//////////////////////////////////////////////////     Err = Ynext - Yest       '//////////////////////////////////////////////////     ' ARMAパラメータ更新     '//////////////////////////////////////////////////     For I = 1 To Order         X(I) = X(I) + G(I) * Err     Next I   End Sub     '///////////////////////////////////////////////// ' '   演算結果 X[a1,a2,...,b1,b2,...] をセルに反映し、 '   EXCEL表計算によるSTR制御演算の結果を更新する ' '/////////////////////////////////////////////////   Sub UpdateExcelSheet(ByVal N)       Dim I As Integer           '//  スライドボリュームによるMactの動特性変更を取得     DoEvents     DoEvents           '//  Excelシートの更新中止     Application.ScreenUpdating = False     Application.Calculation = xlManual           If IdCal = True Then         For I = 1 To (OrderMAX / 2)                   '//  ARMAパラメータ演算結果を a1,a2,...,b1,b2,... に反映             If (I <= (Order / 2)) Then                 Cells(9 + N + Order + 1, 8 + I).Value = X(I)                ' a1,a2,...                 Cells(9 + N + Order + 1, 13 + I).Value = X(I + (Order / 2)) ' b1,b2,...                           Else                 Cells(9 + N + Order + 1, 8 + I).Value = 0                 Cells(9 + N + Order + 1, 13 + I).Value = 0                           End If                       Next I     End If           '//  ARMAパラメータ更新完了     Cells(9 + N + Order + 1, 19).Value = 1           '//  Mactのabパラメータ更新     Cells(9 + N + Order + 1, 7).Value = Cells(6, 7).Value     Cells(9 + N + Order + 1, 8).Value = Cells(6, 8).Value           Cells(7, 20).Value = N  ' シミュレーション回数表示           '//  Excelシートの更新許可     Application.Calculation = xlAutomatic           Calculate  ' グラフ表示更新           Application.ScreenUpdating = True   End Sub

動特性の同定について

Mactの動特性（一次遅れ）がシステムの制御性に影響を与えますが、Mact自体は検出できないとしているので、システムの動特性を別のパラメータから特定する必要があります。

検出可能なのは、制御目標物理量に対し実際に反映できた物理量の比率 Ract です。
　$R_{act} = M_{act} / M_{obj}$
また、システムの制御指示量は $M_{ctr} = K_{str} * M_{obj}$ より、
　$K_{str} = M_{ctr} / M_{obj}$
$M_{ctr}$と$M_{act}$の動特性を、検出可能な$k_{str}$と$R_{act}$の動特性として同定します。

Excelファイルでは、制御サイクルごとに Kstr と Ract をExcelシートから読み込み、Kstr に対する Ract の応答特性を逐次最小二乗推定の計算ロジックで同定します。その計算結果は制御サイクルごとに逐次ExcelシートのARMAモデルのパラメータ（a1,a2,・・・a5, b1,b2・・・,b5）に書き出します。

STR制御ロジック

システム同定ロジックにより、Kstr に対する Ract の応答特性がARMAモデルのパラメータ（ a1,a2,・・・a5,b1,b2,・・・b5) に制御サイクルごとに逐次同定されます。

Kstr に対する Ract の応答特性は $$R_{act}(n+1) = a1*Ract(n)+a2*R_{act}(n-1)+\cdots+b1*K_{str}(n)+b2*K_{str}(n-1)+\cdots$$ です。STR制御では、制御目標物理量に対し実際に反映できた物理量の比率 Ract を制御目標 Robj =1.0 (Mobj= Mact)に制御する必要があります。次回の投入比率 Ract(n+1) が Robj=1.0 になるように $$R_{obj} = a1*R_{act}(n)+a2*R_{act}(n-1)+\cdots+b1*K_{str}(n)+b2*K_{str}(n-1)+\cdots$$ とすると、同定したARMAモデルの逆関数からSTRフィードバック補正係数 Kstr を以下のように算出できます。 $$K_{str}(n) = ( R_{obj} - (a1*R_{act}(n)+a2*R_{act}(n-1)+\cdots+b2*K_{str}(n-1)+\cdots) )/ b1$$

このように適応制御では、オンラインでシステム同定したモデルを使って制御を実現しているため、物バラつきや経年劣化などで動特性が変化する制御対象にも適用が可能です。

Excelファイルの説明

各制御パラメータは以下の計算式です。

$R_{act} = M_{act} / M_{obj}$
$M_{ctr} = K_{str} * M_{obj}$
$M_{act}(n) = a1*M_{act}(n-1)+b1*M_{ctr}(n)$ 
※$M_{act}$の動特性はスライドボリュームで時定数0.8秒から1.2秒の間で随時変更可能
$K_{str}(n)=$
  $(R_{obj} - (a1*R_{act}(n)+a2*R_{act}(n-1)+\cdots+b2*K_{str}(n-1)+\cdots) )/ b1$

Excelの計算シートでは、制御周期100ms想定で計算しています。
システム同定の処理は、VBAで実装しており、制御周期ごとに逐次最小二乗推定で計算した結果をExcel計算シートに反映します。

動かしてみる

制御なし

STR F/Bスタートボタンを押すと、 Kstr  は1.0にリセットされ下のグラフのように制御停止状態になります。 Mobj  と  Mact  は大きくずれており、 Ract  も1.0からずれています。  

STR F/B制御（通常）

その後、STR F/B制御により、 Kstr  で補正されることにより  Mobj  と  Mact  が一致するとともに、 Ract  も1.0にきれいに収束します。

STR F/B制御（Mactの動特性変動）

STR F/B制御中に、Mactの動特性をスライドボリュームで変更することもできます。
下のグラフでは、STR F/B制御稼働中にMactの動特性である時定数を 0.8秒から1.2秒の間で変動させていますが、制御的には大きな破綻は見られません。
オンラインでシステム同定している適応制御の効果と言えますが、それを確認するために、システム同定機能を止めて同じようにMactの動特性を変えてみます。
Excelシートの左上の「システム同定」の右側のセル「1」を「0」に変えるとシステム同定機能が止まります。

システム同定機能を止めた場合

システム同定機能を止めた場合の結果が下のグラフです。
Mactの動特性同定値を時定数 1.0秒固定のまま、STR F/B制御稼働中にMactの動特性である時定数を 0.8秒 から 1.2秒 の間で変動させると制御は安定せず発散してしまいます。

まとめ

適応制御は、オンラインでシステム同定したモデルを使って制御を実現しているため、物バラつきや経年劣化などで動特性が変化する制御対象にも適用が可能であることがわかります。