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

前回の記事では、一階常微分方程式 $\dot{x} = ax$ の解が $x=Ce^{at}$ となることを確認しました。
ここではこの結果を利用して、舵が切れない二輪車のロール運動を表す倒立振子の運動方程式

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

を解いていきます。ただし、車体を支える力は働かないものとして $T_\phi=0$ とします。

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

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

二階常微分方程式は、その一次導関数を新しい変数として導入することで、一階の連立微分方程式に変換することができます。そこで、$\dot{\phi} = x$ とおき、これを微分した $\ddot{\phi} = \dot{x}$ をもとの方程式に代入すると、次のようになります。

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

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

$$\begin{cases} \dot{\phi} = x \\ \dot{x} = - \frac{m_Tz_Tg}{I_{Txx}}\phi \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}$$

階数が下がった代わりに、未知関数は$\phi$と$x$の二つになり、連立一階常微分方程式として表されます。

これは、線形代数で出てくる対角化という変換を行うことで、連立方程式を互いに独立な微分方程式に分解することができ、独立になればそれぞれが $\dot{x}=ax$ と同じ形になるため、前回紹介した方法で解くことができるようになります。

つまり、この連立常微分方程式を

$$\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}$$

という形式に変換できれば、それぞれ独立の一階常微分方程式となって解くことができる、ということです。

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

この連立微分方程式を解く場合、システム行列が対角化可能であれば容易に解くことができ、通常は対角化可能であることから、まずシステム行列の対角化に必要となる固有値と固有ベクトルを求めていきます。

システム行列の固有値・固有ベクトル

得られたシステム行列を$\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}$の固有値$\lambda$は、特性方程式$\det(\lambda\bold{E}-\bold{A})=0$を解くことで得られます。

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

$$\therefore \lambda^2 = k, \quad\lambda = \pm\sqrt{k}$$

ここで、座標系として $z$軸の下向きを正としているため、重心位置は $z_T<0$ であり、

$$k = -\frac{m_Tz_Tg}{I_{Txx}} > 0$$

となるため、固有値 $\lambda$ は実数になります。

次に、求まった固有値と対になる固有ベクトルを求めます。固有ベクトルを$\boldsymbol{v}$とすると、固有値・固有ベクトルの定義式より、

$$\begin{aligned} (\lambda\bold{E}-\bold{A}) \boldsymbol{v} &= \begin{bmatrix} \lambda & -1 \\ -k & \lambda \end{bmatrix} \begin{bmatrix} \phi \\ x \end{bmatrix} \\ &= \begin{bmatrix} \lambda & -1 \\ -\lambda^2 & \lambda \end{bmatrix} \begin{bmatrix} \phi \\ x \end{bmatrix} \\ &= \begin{bmatrix} \lambda\phi - x \\ -\lambda^2\phi + \lambda x \end{bmatrix} \\ &= \begin{bmatrix} \lambda\phi - x \\ -\lambda(\lambda\phi - x) \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \end{bmatrix} \end{aligned}$$

となることから、$\lambda\phi - x=0$ をみたす $[\phi\: x]^T$ が固有ベクトルとなります。そのため固有ベクトルは一意には決まらず、スカラー倍しても同じ方向を表すため無数に存在します。そこで代表として $\phi=1$とすれば

$$\boldsymbol{v} = \begin{bmatrix} 1 \\ \lambda \end{bmatrix}$$

と表すことができます。

つまり、$\lambda_1 = \sqrt{k}$ に対する固有ベクトル $\boldsymbol{v}_1$ および、 $\lambda_2=-\sqrt{k}$ に対する固有ベクトル $\boldsymbol{v}_2$ は、

$$\boldsymbol{v}_1 = \begin{bmatrix} 1 \\ \lambda_1 \end{bmatrix}, \quad \ \boldsymbol{v}_2 = \begin{bmatrix} 1 \\ \lambda_2 \end{bmatrix}$$

と表されます。このように固有ベクトルの $x = \dot{\phi}$　成分は $\phi$ 成分の固有値 $\lambda$ 倍となりますが、これは最終的に得られる解が $e^{\lambda t}$ の形式となるためであり、$\dot{\phi}=\lambda\phi$ となることに対応します。

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

システム行列の対角化

固有ベクトル$\boldsymbol{v}_1, \:\boldsymbol{v}_2$を並べた行列$\bold{P}$を作ります。

$$\bold{P} = [ \boldsymbol{v}_1 \: \boldsymbol{v}_2 ] = \begin{bmatrix} 1 & 1 \\ \lambda_1 & \lambda_2 \end{bmatrix}$$

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

$$\begin{align*} \bold{AP} &= \bold{A}[\boldsymbol{v}_1 \:\boldsymbol{v}_2] \\ &= [\bold{A}\boldsymbol{v}_1 \:\bold{A}\boldsymbol{v}_2] \\ &= [\lambda_1\boldsymbol{v}_1 \:\lambda_2\boldsymbol{v}_2] \\ &= [\boldsymbol{v}_1 \:\boldsymbol{v}_2]\begin{bmatrix} \lambda_1 & 0 \\ 0 & \lambda_2 \end{bmatrix} \\ &= \bold{P} \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 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}$ を計算して運動方程式に代入すると、きちんと成り立っていることが確認できます。

まとめ

以上のように、二階常微分方程式で表される倒立振子の運動方程式

$$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*}$$

ここで $k>0$ であるため、一つの指数項 $e^{\lambda_1t} = e^{\sqrt{k}t}$ は時間とともに増大します。

これは倒立振子が本質的に不安定な系であり，わずかな傾きでも時間とともに増大して倒れてしまうことを数式的に示しています。

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

脚注

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

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