適応制御（セルフチューニングレギュレータ）-EXCELで制御工学！

はじめに

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

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

そうですよねぇ。

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

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

でも、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のコードはこんな感じです。

'/////////////////////////////////////////////////
'//
'//　逐次最小二乗推定
'//
'//   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秒の間で変動させると制御は安定せず発散してしまいます。

まとめ

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

最後までお読みいただきありがとうございました。適応制御の理解のお役に立てましたでしょうか。

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