はじめに
制御工学は基本的に線形システムを対象としています。しかし、実際にはシステムの中に非線形な部分がある場合も多く、非線形な要素が扱えないのでは実用性がないので、ここでは非線形な要素の扱い方について説明します。
また、非線形な要素を平衡動作点で線形近似する方法が一般的ですが、平衡動作点を特定できない場合もあります。そのようなケースの対処方法の一例を示していきます。
非線形なシステムを状態方程式に取り込む
非線形な特性は、そのまま状態方程式に取り込んでも正しい時間応答は得られません。離散時間システムでは、離散化した動作点で非線形特性を線形近似して状態方程式の中に取り込み、線形モデルとして扱うことが必要になります。
車両の駆動力と走行抵抗のモデル
実際に車両が走行する場合、駆動力と車重、走行抵抗のバランスで加速度や走行速度が決まります。
勾配がない平坦路での駆動力と走行抵抗のモデルを考えてみます。
物理式を書き出してみると、車両の駆動力は転がり抵抗と空気抵抗、加減速時のイナーシャとバランスします。
\[F_{tr}=F_{rr}+F_{aero}+F_{w}\]
ここでは、転がり抵抗 \(F_{rr}\) は速度に依存しないので 300[N] 一定とします。
\[\text{空気抵抗 } F_{\text {aero }}=\frac{1}{2} \rho\hspace{0.1cm} C_{d} A_{f} V^{2}=K_{\text {aero }} V^{2} \hspace{0.5cm} \quad K_{\text {aero }}=0.5 \text{ とします。} \]
\[\text{加減速時のイナーシャ } F_{w}=M A=M \frac{d V}{d t} \hspace{1cm} \text{車両重量 } M \text{は 1500kg とします。} \]まとめると
\[F_{t r}=F_{r r}+K_{a e r o} V^{2}+M \frac{d V}{d t}\]
空気抵抗は車速の二乗に比例する非線形な特性です。
次にこの車両の駆動力と走行抵抗の物理式から、\(V\)と\(V^{2}\)を状態変数に選び、状態方程式を生成します。
\[\frac{d V}{d t}=-\frac{K_{\text {aero }}}{M} V^{2}+\frac{1}{M}\left(F_{t r}-F_{r r}\right)\]
\[\frac{d V^{2}}{d t}=\frac{d V^{2}}{d V} \frac{d V}{d t}=2 V \frac{d V}{d t}=-\frac{2 V K_{\text {aero }}}{M} V^{2}+\frac{2 V}{M}\left(F_{t r}-F_{r r}\right)\]
状態方程式にまとめると
\[\frac{d}{d t}\left[\begin{array}{l}V \\ V^{2}\end{array}\right]=\left[\begin{array}{cc}0 & -\frac{K_{\text {aero }}}{M} \\ 0 & -\frac{2 V K_{\text {aero }}}{M}\end{array}\right]\left[\begin{array}{l}V \\ V^{2}\end{array}\right]+\left[\begin{array}{c}\frac{1}{M} \\ \frac{2 V}{M}\end{array}\right]\left(F_{t r}-F_{r r}\right)\]
平衡動作点の近傍での作動であれば、平衡動作点で非線形性を線形近似して状態方程式に取り込めばよいです。
例えば、速度30[m/s]で走行すると、駆動力750[N]、転がり抵抗 𝐹𝑟𝑟= 300[N] 、空気抵抗 𝐹𝑎𝑒𝑟𝑜= 450[N] で平衡状態になります。速度30[m/s]の平衡動作点での挙動を確認したい場合は、速度30[m/s]での状態方程式のAとBを離散化し、得られた逐次式のPとQで時間応答を確認できます。
Excelファイルのダウンロード
非線形モデルのEXCELファイルはこちらからダウンロードできます。
EXCELの状態方程式のシートでは、状態方程式を離散時間モデルの逐次式に変換しています。
シミュレーション (動作平衡点で近似)のシートは、特定の平衡動作点で線形近似した場合のシミュレーション結果を示しています。
シミュレーション (離散時間ごとにPQ更新)のシートは、離散化したすべての動作点でPとQを再計算して線形近似した場合のシミュレーション結果を示しています。
誤差グラフのシートでは、上記の計算方法での空気抵抗の計算結果の精度比較を示しています。
平衡動作点の近傍で時間応答を算出
EXCELの状態方程式のシートの説明です。
速度30[m/s]で走行すると、駆動力750[N]で走行抵抗との平衡動作点になります。その平衡動作点での線形近似値を状態方程式のAとBのパラメータに設定し、離散化した逐次式のPとQを算出しています。
シミュレーション (動作平衡点で近似)のシートは、速度30[m/s]で走行している状態で駆動力\(F_{tr}\)を750[N]から1200[N]に上げた場合について、速度30[m/s]で線形近似したPとQのパラメータで計算したシミュレーション結果を示しています。
平衡動作点が特定できないケースの時間応答算出方法
平衡動作点が特定できない場合は、動作点ごとに非線形特性を線形近似した値を状態方程式のAとBのパラメータに設定する必要があります。時間応答の確認のためには、AとBのパラメータの更新に合わせ、離散化した逐次式のPとQのパラメータも動作点ごとに更新する必要があります。
この場合、時間応答の計算は、EXCELシート上での計算だけでは無理で、ちょっとVBAも使って離散時間（動作点）ごとにPとQを計算しています。
シミュレーション（離散時間ごとにPQ更新）のシートは、速度30[m/s]で走行している状態で駆動力\(F_{tr}\)を750[N]から1200[N]に上げた場合について、PとQのパラメータを動作点ごとに更新した場合のシミュレーション結果を示しています。
システムの非線形性に対しては、どのように対応するかはケースバイケースで考える必要があります。
VBAプログラムの説明
VBAプログラムでPとQのパラメータを動作点ごとに更新しています。
※EXCELの「開発」タグを選択し、「Visual Basic」からプログラムを確認できます。
シミュレーション結果
誤差グラフのシートでは、特定の平衡動作点で線形近似した場合と、動作点ごとに線形近似した場合について、算出される空気抵抗の誤差を比較しています。
※EXCELシート上部のスライドボリュームで駆動力 Ftr を変更できますが、変更した場合は、「非線形モデル時間応答」のボタンを押して再計算する必要があります。
動作点ごとに線形近似した場合は、動作点が動いても誤差は小さいことが分かります。
車両の駆動力と走行抵抗のモデルでは、平衡動作点が特定できないため、動作点ごとに離散化した逐次式のP,Qのパラメータを算出した方がより適切なシミュレーション結果になります。
非線形モデルをどのように状態方程式に取り込むかはシステムの特性や作動条件に応じて検討する必要があります。
最後までお読みいただき、ありがとうございました。
次は現代制御のオブザーバ（状態観測器）になります。