前回の記事で、倒立振子の運動方程式
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*}
として求まることを示しました。
また、運動方程式を状態方程式に変換したときのシステム行列から求めた固有値・固有ベクトルが挙動の基本成分（固有モード）を表すことを紹介しました。今回は、具体的に数値を代入して計算することで、この関係を確認します。
倒立振子の挙動シミュレーション
必要な諸元はMeijaardらの論文に記載されている値を利用します。これらを用いて計算すると、以下のようになりました。
\begin{align*}
I_{Txx} &= 80.82\,\mathrm{kg\,m^2} \\
m_T &= 94.0\,\mathrm{kg},\:
z_T = -0.8612\,\mathrm{m},\:
g = 9.81\,\mathrm{m/s^2} \\
k &= 9.826\,\mathrm{s^{-2}}
\end{align*}
固有値\lambda_1,\:\lambda_2と固有ベクトル\boldsymbol{v_1},\:\boldsymbol{v_2}は次のように求まります。
\lambda_1 = 3.135, \:\boldsymbol{v_1} = \begin{bmatrix} 1 \\ 3.135 \end{bmatrix}, \;\;\;
\lambda_2 = -3.135, \:\boldsymbol{v_2} = \begin{bmatrix} 1 \\ -3.135 \end{bmatrix}
ロール角の初期値を\phi_0 = 1\,\mathrm{deg} = \pi/180\,\mathrm{rad}、ロール角速度の初期値を\dot{\phi}_0 = 0\,\mathrm{rad/s}として、
\begin{bmatrix} C_1 \\ C_2 \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} = \begin{bmatrix} \pi/360 \\ \pi/360 \end{bmatrix}
を代入し、t=0-1\,\mathrm{s}についてプロットすると次のようになります。
ロール角が指数関数的に増加し、転倒していく様子が見て取れます。これは \lambda_1 が正の実数であるためで、対応する項 C_1 e^{\lambda_1t} は時間とともに増大します。一方、\lambda_2 は負の実数であり、C_2e^{\lambda_2t} の項は時間とともに０に収束してきます。そのため、ロール角は最終的に増大する項 C_1 e^{\lambda_1t} に支配され、発散する挙動を示します。
挙動と固有値・固有ベクトルとの関係
前の記事で、物理空間でのロール角とロール角速度によるベクトル\boldsymbol{\Phi} = [\phi\:\:\dot{\phi}]^Tは固有ベクトルを基底とするベクトル空間（モード空間）で表され、
\boldsymbol{\Phi} = \begin{bmatrix} \phi \\ \dot{\phi} \end{bmatrix} = C_1 e^{\lambda_1t} \boldsymbol{v_1} + C_2 e^{\lambda_2t} \boldsymbol{v_2}
となることを紹介しました。\lambda_1が正の実数、\lambda_2が負の実数であることから、\boldsymbol{v_1}の係数は指数関数的に増大し、\boldsymbol{v_2}の係数は０へ収束していくことから、時間の経過とともに\boldsymbol{v_1}成分が支配的になります。
横軸をロール角、縦軸をロール角速度とした位相平面上に軌道をプロットすると次のようになります。時間の経過とともに軌道は直線に近づき、その傾きが固有ベクトル\boldsymbol{v_1}=[1\:\:3.135]の傾き 3.135 に漸近することが確認できます。
一方、ゼロに収束していく固有値\lambda_2と固有ベクトル\boldsymbol{v_2}の組にはどのような意味があるのでしょうか。
固有ベクトル\boldsymbol{v_2}も基底ベクトルの一つであるため、初期状態がこの方向に一致する場合には C_2 e^{\lambda_2t} の項のみが現れます。\lambda_2 は負の実数であるだめ、ロール角は時間とともに０に収束していきます。ロール角の初期値を \phi_0 = 10\,\mathrm{deg} = \pi/18\,\mathrm{rad}, ロール角速度の初期値を \dot{\phi}_0 = -31.35*\pi/180\,\mathrm{rad/s} としてプロットすると、次のようにロール角は \phi =0 の直立状態に収束していきます。
つまり、ロール角とは逆向きの初速度を適切に与えることができれば、倒れることなく直立状態に戻ることができるということを表しているのです。
| 初期状態 |
解析結果（\phi_0 = 10\,\mathrm{deg},\;\dot{\phi}_0 = -31.35\,\mathrm{deg/s}） |
 |
 |
このように、固有値が負の実数であれば挙動は時間とともに０に収束し、このモードは安定と呼ばれます。また、その絶対値は収束の速さを表し、減衰率に対応します。反対に、固有値が正の実数であれば挙動は指数関数的に増大し、このモードは不安定となります。そのため、機械・電気・制御系の設計では、すべての固有値の実部が負となるように設計することが基本となります。
不安定な固有値が存在する場合には、それに対応する固有ベクトルで表されるモードが時間とともに支配的になります。このような不安定な挙動を抑えるためには，どの固有モードが発散の原因となっているかを特定することが重要であり，固有モードの理解は設計改善の手がかりとなります。
固有値が虚数の場合
倒立振子の運動方程式に対して重力の符号を反転し、上下を逆にすると単振り子になります。この場合 k<0 となり固有値は虚数になります。そこで、k'=-kとおくと、\sqrt{k} = i\sqrt{k'} と書けるため、これを解の関数に代入すると、
\phi(t) = C_1 e^{i\sqrt{k'}t} + C_2 e^{-i\sqrt{k'}t}
となり、複素数の指数関数が現れます。この複素指数関数は三角関数に変換できることがわかっており、オイラーの公式として次の関係が成り立つことが広く知られています。
e^{i\theta} = \cos{\theta} + i\sin{\theta}
これを用いれば
\begin{align*}
\phi(t) &= C_1 e^{i\sqrt{k'}t} + C_2 e^{-i\sqrt{k'}t} \\
&= C_1 (\cos{\sqrt{k'}t} + i\sin{\sqrt{k'}t}) + C_2 (\cos{\sqrt{k'}t} - i\sin{\sqrt{k'}t}) \\
&= (C_1 + C_2)\cos{\sqrt{k'}t} + (C_1 - C_2)i\sin{\sqrt{k'}t}
\end{align*}
となります。積分定数 C_1,\:C_2 は、倒立振子の場合と同様に初期条件から求めることができます。ここで、\lambda_2 = -\lambda_1 を用いると、次のように求まります。
\begin{aligned}
\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} \\
&= \begin{bmatrix} \frac{1}{2}\left(\phi_0 + \frac{\dot{\phi}_0}{\lambda_1}\right) \\ \frac{1}{2}\left(\phi_0 - \frac{\dot{\phi}_0}{\lambda_1}\right) \end{bmatrix}
\end{aligned}
したがって、
\begin{aligned}
C_1 + C_2 &= \phi_0 \\
C_1 - C_2 &= \frac{\dot{\phi}_0}{\lambda_1} = \frac{\dot{\phi}_0}{i\sqrt{k'}}
\end{aligned}
となり、これらを代入して整理すると、
\phi(t) = \phi_0\cos{\sqrt{k'}t}+\frac{\dot{\phi}_0}{\sqrt{k'}}\sin{\sqrt{k'}t}
として三角関数で表すことができます。
更に、三角関数の合成公式を用いれば、
\phi(t) = \sqrt{\phi^2_0+\frac{\dot{\phi}^2_0}{k'}}\cos(\sqrt{k'}t - \beta), \;\;
\beta = \tan^{-1}\frac{\dot{\phi}_0}{\phi_0\sqrt{k'}}
として表すことができ、単振動となっていることがわかります。
\phi_0 = 10\,\mathrm{deg}, \:\dot{\phi}_0 = 0\,\mathrm{deg/s}としてプロットすると以下のようになります。
このように、固有値が虚数となる場合には振動が生じることを表します。このとき、固有値の虚部（虚数単位の係数）\sqrt{k'} が振動の角速度を表し、これを固有角振動数と呼びます。
このように、固有値が純虚数の場合、系は減衰も発散もせず、一定振幅の振動を続けることがわかります。
固有値が複素数の場合
上記の単振り子に対して、支点の粘性減衰（摩擦抵抗）を考慮すると、運動方程式には減衰項が加わって次のように表されます。
I_{Txx}\ddot{\phi} +c\dot{\phi} +m_Tz_Tg\phi = 0
この場合も同じように固有方程式を解いて固有値を求めると、
\lambda^2 + \frac{c}{I_{Txx}}\lambda + \frac{m_Tz_Tg}{I_{Txx}} = 0
\therefore \lambda = -\frac{c}{2I_{Txx}} \pm\sqrt{\Big(\frac{c}{2I_{Txx}}\Big)^2 - \frac{m_Tz_Tg}{I_{Txx}}}
となります。このとき、減衰係数 c が十分小さい場合には、判別式が負となり固有値は複素数になります。そこで、\lambda = -a \pm bi (a>0,\:b>0)とおいて微分方程式を解くと、解の関数は次のように求まります。
\begin{aligned}
\phi(t) &= C_1e^{(-a+bi)t} + C_2e^{(-a-bi)t} \\
&= e^{-at}\{ (C_1+C_2)\cos{bt} +(C_1-C_2)i\sin{bt} \} \\
\end{aligned}
積分定数 C_1,\:C_2 は、これまでと同様に初期条件から求められ、
\begin{aligned}
C_1 + C_2 &= \phi_0 \\
C_1 - C_2 &= -\frac{i}{b}(a\phi_0 + \dot{\phi}_0)\\
\end{aligned}
となります。これらを代入して整理すると、
\phi(t) = e^{-at} \left( \phi_0\cos{bt} + \frac{a\phi_0 + \dot{\phi}_0}{b}\sin{bt} \right)
と表されます。
c=50 \,\mathrm{Nm\,s/rad}, \:\phi_0 = 10 \,\mathrm{deg}, \:\dot{\phi}_0 = 0\,\mathrm{deg/s} としてプロットすると次のようになります。
このように固有値が複素数の場合、指数関数と三角関数の積で表される減衰振動となります。固有値の実部 -a は減衰率を、虚部 b は角振動数を表しており、これらが挙動を特徴づけています。
このように固有値が複素数の場合、時間の経過とともに振幅は指数関数的に減衰しながら振動することがわかります。
まとめ
運動方程式を変換して得られる状態方程式から、システム行列の固有値・固有ベクトルを求めることで、システムがどのようにふるまうかを知ることができます。
固有値の実部が負であればゼロに収束し安定、正であれば発散し不安定となります。不安定な固有値が一つでもあれば発散する挙動を示すため、通常はすべての固有値の実部が負となるように設計します。
また、固有値が複素数の場合には減衰振動が生じ、その実部が減衰率、虚部が角振動数を表します。
さらに、固有値に対応する固有ベクトルは、それぞれのモードにおける状態量の比（モードの形）を表しています。
次の記事から、操舵系の運動も含む二輪車の運動方程式に対して同様の解析を行い、固有値・固有ベクトルと挙動との関係をさらに詳しく見ていきます。
Discussion