非線形なシステムへの対応

\[F_{t r}=F_{r r}+K_{aero} V^{2}+M \frac{d V}{d t}\]

空気抵抗は車速の二乗に比例する非線形な特性です。
この物理式において、車速\(V\)を状態変数に選び、状態方程式の形に書き直すと
\[\frac{d V}{d t}= -\frac{K_{aero}}{M} V^{2}+\frac{1}{M}(F_{t r}-F_{r r}) \]

\(V^{2}\)は非線形なパラメータなので、平衡状態での線形近似を考えます。

例えば、速度30[m/s]で走行すると、駆動力 \(F_{tr}\) =750[N]、転がり抵抗 \(F_{rr}\) = 300[N] 、空気抵抗 \(F_{aero}\) = 450[N] で平衡状態になります。このような平衡動作点での速度を\(V_{e}\)、平衡動作点からの速度の変化を\(V_{d}\)、駆動力 \(F_{tr}\)から転がり抵抗\(F_{rr}\)を引いた平衡動作点での推進力を\(F_{e}\)、平衡動作点からの推進力の変化を\(F_{d}\)とすると、平衡動作点近傍での速度\(V\)と推進力\(F\)は\[V=V_{e}+V_{d} \]\[F=F_{e}+F_{d}\]

状態方程式を平衡動作点で線形近似（微分）して、平衡動作点からの偏差の式にすると\[\frac{d V_{d}}{d t}=-\frac{2 V_{e} K_{aero }}{M} V_{d}+\frac{1}{M}F_{d}\]と表せます。

平衡動作点の近傍での作動であれば、このように平衡動作点で非線形性を線形近似して状態方程式に取り込んだモデルで扱うことができます。更に、その状態方程式の A と B を離散化して得られた逐次式の P と Q で時間応答を確認することができます。

離散化についてはこちら→リンク

平衡動作点が特定できないケースの時間応答算出方法

実際には平衡動作点が特定できないケースも多いです。今回の車両の駆動力と走行抵抗のモデルでも、車両は如何なる速度でも定速走行が可能なので平衡動作点は特定できません。

更に、今回の車両の駆動力と走行抵抗のモデルでは、線形近似した状態方程式に、起点となる平衡動作点の車速 \(V_{e}\) を含んでいます。\[-\frac{2 V_{e} K_{aero }}{M} V_{d}\]

起点となる平衡動作点(例えば \(V_{e}=30[m/s]\))の近傍の車速域では、ある程度精度良く近似できますが、車速 \(V_{e}\)から 速度（動作点）がずれるとモデルの精度が低下します。

このような場合は、動作点ごとに線形近似した値を状態方程式の A と B に設定することにより精度を確保します。

時間応答の確認のためには、A と B の更新に合わせ、離散化した逐次式の P と Q も動作点ごと（離散時間ごと）に更新する必要があります。

EXCELでの時間応答の確認ですが、ちょっとVBAも使います。

Excelファイルの説明

非線形モデルのEXCELファイル（⑧非線形モデル(rev2).xls）はこちらからダウンロードできます。

EXCELの状態方程式のシートでは、状態方程式を離散時間モデルの逐次式に変換しています。

シミュレーションのシートでは、速度 V=30[m/s] で走行している状態から駆動力 \(F_{tr}\) を750[N]から1200[N]に変化させた場合について、P と Q を動作点ごと（離散時間ごと）に更新（線形近似を再計算）した場合のシミュレーション結果を示しています。

平衡動作点で線形近似

EXCELの状態方程式のシートの説明です。
速度30[m/s]で走行すると、駆動力750[N]で走行抵抗との平衡動作点になります。その平衡動作点で線形近似した状態方程式を作成。その状態方程式の A と B のパラメータを離散化した逐次式の P と Q を算出しています。

動作点ごと（離散時間ごと）に線形近似を再計算

EXCELのシミュレーション のシートは、速度30[m/s]で走行している状態から駆動力\(F_{tr}\)を750[N]から1200[N]に変化させた場合のシミュレーション結果を示しています。

シミュレーションの精度を確保するために、VBAも使って離散時間（動作点）ごとに P と Q を再計算し更新しています。

車両の駆動力と走行抵抗のモデルでシステムの非線形性の対応例を説明しましたが、非線形性を持った様々なシステムに、どのように対応するかはケースバイケースで考える必要があります。

VBAプログラムの説明

VBAプログラムでPとQのパラメータを動作点ごとに更新しています。
※EXCELの「開発」タグを選択し、「Visual Basic」からプログラムを確認できます。
VBAのコードはこんな感じです。

Sub nonlinearmodel()

    Dim CalTimes As Integer ' 演算回数
    Dim V, Vd, Fd As Double

    '/////////////////////////////////////////////////
    ' 初期化
    '/////////////////////////////////////////////////

    CalTimes = Worksheets("シミュレーション").Range("E9").Value
    V = Worksheets("シミュレーション").Cells(17, 6).Value

    For I = 1 To CalTimes
    
        ' Excelシートの更新中止
        Application.ScreenUpdating = False
        Application.Calculation = xlManual
        
        '/////////////////////////////////////////////////
        ' 動作点の速度Vを設定
        '/////////////////////////////////////////////////
    
        Worksheets("状態方程式").Range("C6").Value = V
    
        '/////////////////////////////////////////////////
        ' シミュレーションシートの逐次行列 P, Q は自動更新
        '/////////////////////////////////////////////////
    
        '/////////////////////////////////////////////////
        ' シミュレーションシートの逐次式の Fd, Vd を更新
        '/////////////////////////////////////////////////

        Vd = Worksheets("シミュレーション").Cells(16 + I, 5).Value
        Fd = Worksheets("シミュレーション").Cells(16 + I, 4).Value
    
        Worksheets("シミュレーション").Range("E12").Value = Vd
        Worksheets("シミュレーション").Range("H12").Value = Fd
    
        ' Excelシートの更新許可
        Application.Calculation = xlAutomatic
        Application.ScreenUpdating = True
    
        ' Calculate  ' グラフ表示更新
    
    
        '/////////////////////////////////////////////////
        ' 計算結果の速度変化量 Vd(N+1)を保存
        '/////////////////////////////////////////////////
    
        Vd = Worksheets("シミュレーション").Range("B12").Value
        Worksheets("シミュレーション").Cells(17 + I, 5).Value = Vd
    
        '/////////////////////////////////////////////////
        ' 動作点の速度Vを取得
        '/////////////////////////////////////////////////
    
        V = Worksheets("シミュレーション").Cells(17 + I, 6).Value

    Next I

End Sub

シミュレーション結果

EXCELのシミュレーション のシートでは、【calculate】のボタンを押すと、速度30[m/s]で走行している平衡状態から、駆動力\(F_{tr}\)を750[N]から1200[N]に変化させた場合のシミュレーション結果を表示します。P と Q を動作点ごと（離散時間ごと）に更新（線形近似を再計算）しています。

※駆動力は転がり抵抗 \(F_{rr}\)  = 300[N] を差し引いているので、450[N]から900[N]の変化になっています。

【reset】のボタンで初期状態に戻ります。

非線形モデルをどのように状態方程式に取り込むかはシステムの特性や作動条件に応じて検討する必要があります。

今回はシステムの非線形な特性の扱い方について書きました。非線形な特性を線形近似して状態方程式に取り込むことで、時間応答の確認だけでなく、この後説明していく現代制御などのロジックの適用にもつなげることができます。

最後までお読みいただき、ありがとうございました。

次は現代制御のオブザーバ（状態観測器）になります。