OpenNH

日常のひとこま(自分用のメモとかあれこれ)

セミアクティブサスペンションの最適制御とMATLABサンプルコード

Butsuen氏、TetsuroN氏の”THE DESIGN OF SEMI-ACTIVE SUSPENSION S FOR AUTOMOTIVE VEHICLES”(1989)を参考にしてセミアクティブサスペンションの最適制理論を記載しています。
DSpace@MIT: The design of semi-active suspensions for automotive vehicles

セミアクティブサスペンションの単輪モデル

f:id:FounderLeis:20181001151321p:plain
Semi-active suspension quarter model

運動方程式
 m_2\ddot{z}_2 = -c(\dot{z}_2-\dot{z}_1)-k_2(z_2-z_1)+f_c  \tag{1}
 m_1\ddot{z}_1 =  c(\dot{z}_2-\dot{z}_1)+k_2(z_2-z_1)-c(\dot{z}_1-\dot{z}_0)+k_1(z_1-z_0)-f_c  \tag{2}

状態方程式
\begin{align}
\dot{x} &= Ax + Bf_c+ W\dot{z}_0 \\
&= Ax + B(x)c_{semi} + W\dot{z}_0 \tag{3}
\end{align}

状態変数
\[ x =
\left[ \begin{array}{c}
z_2-z_1 \\
\dot{z}_2 \\
z_1-z_0 \\
\dot{z}_1
\end{array}\right]
= \left[ \begin{array}{c}
{x_1 \\ x_2 \\ x_3 \\ x_4 }
\end{array}
\right] \tag{4}
\]
{\displaystyle x_1 } : サスペンションストローク
{\displaystyle x_2 } : ばね上上下速度
{\displaystyle x_3 } : タイヤたわみ
{\displaystyle x_4 } : ばね下上下速度


\[
A = \left[ \begin{array}{cccc}
0 & 1 & 0 & -1 \\
-\frac{k_2}{m_2} & -\frac{c_2}{m_2} & 0 & \frac{c_2}{m_2} \\
0 & 0 & 0 & 1 \\
\frac{k_2}{m_1} & \frac{c_2}{m_1} & -\frac{k_1}{m_1} & -\frac{c_2}{m_1} \\
\end{array}\right] ~~~~
\]
\[
B(x) = \left[ \begin{array}{c}
0 \\ -\frac{1}{m_2} \\ 0 \\ \frac{1}{m_1}
\end{array}\right]
(x_2 - x_4)
=-{B}_0 (x_2 - x_4),~~~~
W = \left[ \begin{array}{c}
0 \\ 0 \\ -1\\ 0
\end{array}\right]
\]

最適制御理論

評価関数を以下のように設計します。
\begin{equation}
J = \lim_{T \to \infty} \frac{1}{T}E \left[ \int_0^T \left( \ddot{z}^2_2 + \rho_1 (z_2-z_1)^2 + \rho_2 \dot{z}^2_2 + \rho_3 (z_1-z_0)^2 + \rho_4 \dot{z}^2_1 \right) dt \right] \tag{5}
\end{equation}
{\displaystyle ~~~~~~~~~~\rho_1, \rho_2, \rho_3, \rho_4} : 重み要素
ばね上加速度、サスペンションストローク、ばね上速度、タイヤたわみ、ばね下速度を評価値としている。基本的に、乗り心地にはばね上加速度、接地性にはタイヤたわみ、制御制約にサスペンションストロークの要素を用います。

上記の評価関数 {\displaystyle  J } は、状態量表記にするために以下のように書き換えられます。
\begin{equation}
J = \lim_{T \to \infty} \frac{1}{T}E \left[ \int_0^T {x^T(Q_0+Q_{semi})x} dt \right] \tag{6}
\end{equation}

ここで {\displaystyle  Q_{semi} } は可変ダンパー係数 {\displaystyle  c_{semi}} を含む項を表しています。

路面入力(外乱入力)に対して、状態量{\displaystyle  x }はランダムプロセスとなりますが、初期値は {\displaystyle  \dot{z}_r=0 }{\displaystyle  x(0) \neq \bf{0}}を想定します。

最適減衰係数{\displaystyle  c^*_{semi}}は,制約条件[式(8), (9)]、初期条件[式 (10)]のもとで評価関数{\displaystyle  J}[式(7)]を最小化するようにして得られます。
\[ J = \int_0^\infty~x^T(Q_0+Q_{semi})x~ dt \tag{7} \]
\[ \dot{x} = Ax+B(x)c_{semi} \tag{8} \]
\[ 0 \leq c_{semi} \leq c_{max} \tag{9} \]
\[ x(0) = x_0 \tag{10} \]

ここで、
\[ Q_0 = \left[ \begin{array}{cccc}
\frac{k^2_2}{m^2_2} + \rho_1 & \frac{c_2k_2}{m^2_2} & 0 & -\frac{c_2k_2}{m^2_2} \\
\frac{c_2k_2}{m^2_2} & \frac{c^2_2}{m^2_2} + \rho_2 & 0 & -\frac{c^2_2}{m^2_2} \\
0 & 0 & \rho_3 & 0 \\
-\frac{c_2k_2}{m^2_2} & -\frac{c^2_2}{m^2_2} & 0 & \frac{c^2_2}{m^2_2} + \rho_4 \\
\end{array}\right]
\tag{11} \]
\[ Q_{semi} = \left[ \begin{array}{cccc}
0 & \frac{c_{semi}k_2}{m^2_2} & 0 & -\frac{c_{semi}k_2}{m^2_2} \\
\frac{c_{semi}k_2}{m^2_2} & \frac{(2c^2_2+c_{semi})c_{semi}}{m^2_2} & 0 & -\frac{(2c^2_2+c_{semi})c_{semi}}{m^2_2} \\
0 & 0 & 0 & 0 \\
-\frac{c_{semi}k_2}{m^2_2} & -\frac{(2c^2_2+c_{semi})c_{semi}}{m^2_2} & 0 & \frac{(2c^2_2+c_{semi})c_{semi}}{m^2_2} \\
\end{array}\right]
\tag{12} \]

次に、式(7)を以下のように書き換えます。
\[ J = \int_0^\infty~(x^TQ_0x + 2c_{semi}^TS(x)x + c_{semi}^TR(x)c_{semi})~ dt \tag{13} \]
ここで、
\[ R(x) = \frac{1}{m_2^2} (x_2-x_4)^2 = R_0(x_2-x_4)^2 \tag{14} \]
\[ S(x) = -\left[ -\frac{k_2}{m_2^2} ~~~~-\frac{c_2}{m_2^2} ~~~~0 ~~~~\frac{c_2}{m_2^2} \right](x_2-x_4) = -S_0(x_2-x_4) \tag{15} \]

この時、
・可変減衰係数{\displaystyle  c_{semi}}に対する拘束条件を考えない場合 (i.e. {\displaystyle  -\infty \le c_{semi} \le \infty})、最適減衰係数{\displaystyle  c^*_{semi}}は以下の関係式で得ることが出来ます。
\begin{align}
c^*_{semi} &= -R(x)^{-1}[B^T(x)P+S(x)]x ~~~~~~~~~~ &if&~~~~~ x_2 \neq x_4 \tag{16}\\
c^*_{semi} &= 0 ~~~~~~~~~~&if&~~~~~ x_2 = x_4 \tag{17}
\end{align}

ここで、{\displaystyle P}は式(18)に示すRiccati方程式の解です。
\[ P\bar{A} + \bar{A}^T P + \bar{Q} - PB_0R_0^-1B_0^TP = 0 \tag{18} \]
ただし、
\begin{align}
\bar{A} &= A - B_0 R_0^{-1} S_0 \\
\bar{Q} &= Q_0 - S_0^T R_0^{-1} S_0
\end{align}
次いで、
・可変減衰係数{\displaystyle  c_{semi}}に対する拘束条件を考える場合 (i.e. {\displaystyle  0 \le c_{semi} \le c_{max}})、最適減衰係数{\displaystyle  c^*_{semi}}は以下の関係式で得ることが出来ます。
この式がセミアクティブサスペンション最適制御での制御値となります。
\begin{align}
\rm(i)~~~\it c^*_{semi} &= 0 ~~~~~~~~~~&if&~~~~~ (B^T(x)P+S(x))x \ge 0 \tag{19} \\
\rm(ii)~~~\it c^*_{semi} &= -R(x)^{-1}[B^T(x)P+S(x)]x ~~~~~~~~~~ &if&~~~~~ -R(x)c_{max} < (B^T(x)P+S(x))x < 0 \tag{20}\\
\rm(iii)~~~\it c^*_{semi} &= c_{max}~~~~~~~~~~&if&~~~~~ (B^T(x)P+S(x))x \le -R(x)c_{max} \tag{21}
\end{align}

f:id:FounderLeis:20181002163128p:plain

非線形減衰力特性

サスペンションのダンパーはストローク速度によって発生力が変化しますが、それは線形ではなくダンパーのヒステリシスによって非線形に変化していきます。

そこで、Gordon氏の論文[*1]に記載されていた非線形減衰力特性マップを利用しようと思います。論文の中では減衰力切り替え指令値として{\displaystyle  n}という無次元変数が使われていましたが、数式的に意味を持たせるために減衰比{\displaystyle  \zeta}を導入しました。

ある減衰比{\displaystyle  \zeta}とサスペンションストローク速度{\displaystyle  v_s}が与えられた場合のダンパー発生力{\displaystyle  F_{damp}}は以下のようにして計算されます。
{  F_{damp} = (200+4800\zeta)arctan(2v_s) \tag{22} }
\[ \zeta = \frac{c_{semi}}{2\sqrt{m_2k_2}} \]

f:id:FounderLeis:20181002173659p:plain
Non-linear damper force map

 
 

MATLAB/Simulinkによる制御シミュレーション

GitHubMATLAB/Simulinkのファイルを置いてあるのでダウンロードして使用してください。
github.com

上記のプログラムを実行すると以下のような結果が表示されるはずです。
<線形ダンパーを用いた結果>
f:id:FounderLeis:20181004135754p:plain

非線形ダンパーを用いた結果>
f:id:FounderLeis:20181004135803p:plain