舵が切れなければ二輪車は倒れる(3) 倒立振子運動方程式の解

一階常微分方程式$\dot{x} = ax$の解が$x=C_0e^{at}$となることを利用して舵が切れない二輪車の運動方程式

$$I_{Txx}\ddot{\phi} +m_Tz_Tg\phi = T_\phi$$

を解いていきます。

運動方程式などのような二階以上の高階の常微分方程式も、一階常微分方程式に変換することで同じように解くことができるので、その方法を順を追って紹介します。

二階常微分方程式の一階化

二階常微分方程式を一階常微分方程式に変換するため、$\dot{\phi} = x$とおきます。
これを微分すると、$\ddot{\phi} = \dot{x}$であり、もとの方程式は次のようになります。

$$I_{Txx}\dot{x} +m_Tz_Tg\phi = T_\phi$$

これを$x$の定義式と連立させると、

$$\begin{cases} \dot{\phi} = x \\ \dot{x} = - \frac{m_Tz_Tg}{I_{Txx}}\phi + \frac{T_\phi}{I_{Txx}} \end{cases}$$

となりますが、これを行列形式でまとめれば次のように一階の連立常微分方程式になります。

$$\frac{d}{dt} \begin{bmatrix} \phi \\ x \end{bmatrix} = \begin{bmatrix} 0 & 1 \\ - \frac{m_Tz_Tg}{I_{Txx}} & 0 \end{bmatrix} \begin{bmatrix} \phi \\ x \end{bmatrix} + \begin{bmatrix} 0 \\ \frac{1}{I_{Txx}} \end{bmatrix} T_\phi$$

階数が下がった代わりに、未知関数は$\phi$と$x$の二つに増えてしまいましたが、線形代数で出てくる対角化という変換を行うことで、それぞれの関数について独立の常微分方程式として扱うことができ、実際に解くことができるようになります。

このような形式の常微分方程式を現代制御理論では状態方程式と呼ばれ、$[\phi\:\:x]^T = [\phi\:\:\dot{\phi}]^T$が運動の状態を表す最小情報であることから状態変数と呼ばれます。また、状態変数ベクトルに掛かる係数行列をシステム行列、外力$T_{\phi}$に掛かる係数行列を入力行列と呼んでいます。

任意関数である外力は$T_\phi = 0$として、この連立微分方程式を解くために、まずシステム行列の対角化に必要となる固有値と固有ベクトルを求めていきます。

係数行列の固有値・固有ベクトル

得られたシステム行列を$\bold{A}$とし、記述を簡潔にするために$k=-\frac{m_Tz_Tg}{I_{Txx}}$と置くことにします。

$$\bold{A} = \begin{bmatrix} 0 & 1 \\ -\frac{m_Tz_Tg}{I_{Txx}} & 0 \end{bmatrix} = \begin{bmatrix} 0 & 1 \\ k & 0 \end{bmatrix}$$

まず、行列$\bold{A}$の固有値を求めるために、固有方程式$|\bold{A}-\lambda\bold{E}|=0$を解きます。

$$\begin{vmatrix} -\lambda & 1 \\ k & -\lambda \end{vmatrix} = \lambda^2 - k = 0$$

$$\therefore \lambda = \pm\sqrt{k}$$

ここで、$z$軸は下向きが正であり、$z_T<0, \:k>0$で、固有値は実数になります。

次に、求まった固有値と対になる固有ベクトルを求めます。$\lambda_1 = \sqrt{k}, \:\lambda_2 = -\sqrt{k}$とし、それぞれに対応する固有ベクトルを$\bold{v_1}, \:\bold{v_2}$とすると、固有値・固有ベクトルの定義式より、

$\lambda_1 = \sqrt{k}$のとき、

$$(\bold{A}-\lambda_1\bold{E}) \bold{v_1} = \begin{bmatrix} -\sqrt{k} & 1 \\ k & -\sqrt{k} \end{bmatrix} \begin{bmatrix} \phi_1 \\ \dot{\phi}_1 \end{bmatrix} = \begin{bmatrix} -\sqrt{k}\phi_1 + \dot{\phi}_1 \\ k\phi_1 - \sqrt{k}\dot{\phi}_1 \end{bmatrix} = \bold{0}$$

$$\therefore \bold{v_1} = \begin{bmatrix} \phi_1 \\ \dot{\phi}_1 \end{bmatrix} = s_1 \begin{bmatrix} 1 \\ \sqrt{k} \end{bmatrix} = s_1 \begin{bmatrix} 1 \\ \lambda_1 \end{bmatrix}$$

$\lambda_2 = -\sqrt{k}$のとき、

$$(\bold{A}-\lambda_2\bold{E}) \bold{v_2} = \begin{bmatrix} \sqrt{k} & 1 \\ k & \sqrt{k} \end{bmatrix} \begin{bmatrix} \phi_2 \\ \dot{\phi}_2 \end{bmatrix} = \begin{bmatrix} \sqrt{k}\phi_2 + \dot{\phi}_2 \\ k\phi_2 + \sqrt{k}\dot{\phi}_2 \end{bmatrix} = \bold{0}$$

$$\therefore \bold{v_2} = \begin{bmatrix} \phi_2 \\ \dot{\phi}_2 \end{bmatrix} = s_2 \begin{bmatrix} 1 \\ -\sqrt{k} \end{bmatrix} = s_2 \begin{bmatrix} 1 \\ \lambda_2 \end{bmatrix}$$

※ $s_1, \:s_2$は任意定数

以上のように得られた固有ベクトルを用いて、行列の対角化を行っていきます。

係数行列の対角化

$s_1 = s_2 = 1$とし、固有ベクトル$\bold{v_1}, \:\bold{v_2}$を並べた行列$\bold{P}$を作ります。

$$\bold{P} = [ \bold{v_1} \: \bold{v_2} ] = \begin{bmatrix} 1 & 1 \\ \lambda_1 & \lambda_2 \end{bmatrix}$$

係数行列$\bold{A}$の右から$\bold{P}$を掛け、固有値・固有ベクトルの定義を用いて計算すると、

$$\begin{align*} \bold{AP} &= \bold{A}[\bold{v_1} \:\bold{v_2}] \\ &= [\bold{Av_1} \:\bold{Av_2}] \\ &= [\lambda_1\bold{v_1} \:\lambda_2\bold{v_2}] \\ &= [\bold{v_1} \:\bold{v_2}]\begin{bmatrix} \lambda_1 & 0 \\ 0 & \lambda_2 \end{bmatrix} \end{align*}$$

となり、対角項に固有値を並べた対角行列が現れます。
ここで、固有値$\lambda_1, \lambda_2$を対角項に並べた行列を、

$$\bold{\Lambda} = \begin{bmatrix} \lambda_1 & 0 \\ 0 & \lambda_2 \end{bmatrix}$$

とすると、

$$\bold{AP} = \bold{P\Lambda}$$

と表せます。

連立常微分方程式の対角化

$\bold{\Phi} = [\phi\:\:x]^T = [\phi\:\:\dot{\phi}]^T$とおき、一階化したロールの運動方程式を次のように表します。

$$\frac{d{\bold{\Phi}}}{dt} = \bold{A\Phi}$$

$\bold{AP}=\bold{P\Lambda}$より、両辺に右から$\bold{P^{-1}}$をかけると、$\bold{A}=\bold{P\Lambda P^{-1}}$が得られ、上記の運動方程式に代入すると次のようになります。

$$\frac{d{\bold{\Phi}}}{dt} = \bold{P\Lambda P^{-1}\Phi}$$

これに両辺に左から$\bold{P^{-1}}$を掛けると、

$$\frac{d}{dt}(\bold{P^{-1}\Phi}) = \bold{\Lambda P^{-1}\Phi}$$

となるので、

$$\bold{P^{-1}\Phi} = \bold{y} = \begin{bmatrix} y_1 \\ y_2 \end{bmatrix}$$

とおけば、

$$\frac{d}{dt}\begin{bmatrix} y_1 \\ y_2 \end{bmatrix} = \begin{bmatrix} \lambda_1 & 0 \\ 0 & \lambda_2 \end{bmatrix} \begin{bmatrix} y_1 \\ y_2 \end{bmatrix} = \begin{bmatrix} \lambda_1 y_1 \\ \lambda_2 y_2 \end{bmatrix}$$

$$\therefore \begin{cases} \dot{y}_1 = \lambda_1 y_1 \\ \dot{y}_2 = \lambda_2 y_2 \end{cases}$$

となります。このように対角化によって第1式は$y_1$、第2式は$y_2$についてのそれぞれ独立した常微分方程式となり、一階常微分方程式の解法で紹介した、$\dot{x} = ax$と同じ形式になっているため、解くことができます。

ロール運動方程式の解

それでは実際にロール運動方程式の解を計算していきましょう。

一階常微分方程式 $\dot{x} = ax$の解は、$x = C_0 e^{at}$ でしたので、$\bold{y}$は次のように求まります。

$$\bold{y} = \begin{bmatrix} y_1 \\ y_2 \end{bmatrix} = \begin{bmatrix} C_1 e^{\lambda_1 t} \\ C_2 e^{\lambda_2 t} \end{bmatrix}$$

また、$\bold{y} = \bold{P^{-1}\Phi}$より、

$$\bold{\Phi} = \bold{Py} = [\bold{v_1} \:\bold{v_2}] \begin{bmatrix} y_1 \\ y_2 \end{bmatrix} = C_1 e^{\lambda_1 t} \bold{v_1} + C_2 e^{\lambda_2 t} \bold{v_2}$$

となります。このように$\bold{\Phi}$は固有ベクトルを基底とし、対応する固有値で規定される指数関数を係数による線形結合で表されることがわかります。

この固有値・固有ベクトルの組が挙動の基本となる成分になることからモードと呼び、実際の物理的な挙動をモードの和で表されることをモード展開と呼んで、非常に効率よく計算できることから振動工学や制御工学など様々な分野で利用される手法となっています。^[1]
物理座標$\bold{\Phi}$に対応する$\bold{y}$はモード座標と呼ばれ、物理座標とモード座標の変換に使った固有ベクトルを並べた行列$\bold{P}$はモード行列と呼ばれています。^[2]

最後に $t=0$ における$\phi, \:\dot{\phi}$ の初期値をそれぞれ $\phi_0, \:\dot{\phi}_0$ とし、モード座標での初期値$C_1, \:C_2$を求めると、$\bold{y} = \bold{P^{-1}\Phi}$ より、

$$\begin{align*} \begin{bmatrix} C_1 \\ C_2 \end{bmatrix} &= \frac{1}{\lambda_2-\lambda_1} \begin{bmatrix} \lambda_2 & -1 \\ -\lambda_1 & 1 \end{bmatrix} \begin{bmatrix} \phi_0 \\ \dot{\phi}_0 \end{bmatrix} \\ &= -\frac{1}{2\sqrt{k}} \begin{bmatrix} -\sqrt{k} & -1 \\ -\sqrt{k} & 1 \end{bmatrix} \begin{bmatrix} \phi_0 \\ \dot{\phi}_0 \end{bmatrix} \\ &= \begin{bmatrix} \frac{1}{2}\phi_0 + \frac{1}{2\sqrt{k}}\dot{\phi}_0 \\ \frac{1}{2}\phi_0 - \frac{1}{2\sqrt{k}}\dot{\phi}_0 \end{bmatrix} \end{align*}$$

$$\therefore \begin{bmatrix} \phi \\ \dot{\phi} \end{bmatrix} = \begin{bmatrix} ( \frac{1}{2}\phi_0 + \frac{1}{2\sqrt{k}}\dot{\phi}_0 ) e^{\sqrt{k}t} + ( \frac{1}{2}\phi_0 - \frac{1}{2\sqrt{k}}\dot{\phi}_0 )e^{-\sqrt{k}t} \\ \sqrt{k} ( \frac{1}{2}\phi_0 + \frac{1}{2\sqrt{k}}\dot{\phi}_0 ) e^{\sqrt{k}t} - \sqrt{k} ( \frac{1}{2}\phi_0 - \frac{1}{2\sqrt{k}}\dot{\phi}_0 ) e^{-\sqrt{k}t} \end{bmatrix}$$

となります。

$\phi(0) = \phi_0, \:\dot{\phi}(0) = \dot{\phi}_0$ や、$\phi$を微分すると$\dot{\phi}$となること、および $\ddot{\phi}$ を計算して運動方程式に代入すると、きちんと成り立っていることが確認できます。

まとめ

以上のように、舵角$\delta=0$で固定した二輪車のロール方向の運動方程式

$$I_{Txx}\ddot{\phi} +m_Tz_Tg\phi = 0$$

の解は

$$\begin{align*} \phi(t) &= C_1 e^{\lambda_1t} + C_2 e^{\lambda_2t} \\ &= \Big(\frac{1}{2}\phi_0 + \frac{1}{2\sqrt{k}}\dot{\phi}_0\Big) e^{\sqrt{k}t} + \Big(\frac{1}{2}\phi_0 - \frac{1}{2\sqrt{k}}\dot{\phi}_0 \Big) e^{-\sqrt{k}t}, \:\: k = -\frac{m_Tz_Tg}{I_{Txx}} \end{align*}$$

として求まりました。次の記事では数値を代入してどのような結果が得られるかを見ることで、固有値・固有ベクトルと挙動との関係を見ていきたいと思います。

脚注

1. 森泰親, 「大学講義テキスト　現代制御」, コロナ社, 2020, p.28, ISBN978-4-339-03231-4

2. 長松昭男, 「モード解析入門」, コロナ社, 1993, p.95, ISBN4-339-08225-2