はじめに
PCのフォルダを検索していたら、ずいぶん昔、苦労しながらがんばって和訳した論文を見つけた。
当時でも30年ぐらい前の論文だったが、今の二輪車の安定性解析のベースにもつながっている論文 「The Stability and Control of Motorcycles：R. S. Sharp（1971）」。
今どきネットで検索すれば原文のpdfがサイトから簡単にダウンロードできるし、それをGenAIに放り込めば、和訳どころか解説までしてくれるので、詳細はそちらにお任せするとして。。。
この和訳みつけて思い出したのが、面白かったのはもちろんだが、この論文を理解するためのベース知識がこれまた勉強にもなったなぁと、で、そのあたりの流れを含めて覚書化して残しておこうと。
（うん十年前の手書き資料も引っ張り出しながらなので、思い出すにも時間はかかりそう・・・）。
ちなみに当時和訳したのが以下（個人用なので意訳あり（+たぶん誤訳も））
[ページ送りはヘッダー＆フッターの「 ⤵」をクリック]
pdf版のダウンロードは以下より（有料：コーヒー代 かせぎ）
この論文が画期的だった点
Sharpの論文が画期的だった点は、二輪車特有の以下の３つの不安定モード
- ウォブルモード（シミー）: ステアリング軸周りの振動
- ウィーブモード : 車体全体が縫うように左右へ揺れる動き
- キャップサイズモード : 低速時に車両がゆっくりと倒れ込む挙動
に対して、「フロントボディ」と「リアボディ」を結合した4自由度モデルを採用、さらに「タイヤの緩和現象（サイドフォースの発生遅れ）」を数式に取り込み解析を行った結果、このモデルの４つの固有振動モードに上記の不安定モード３つが包括されており、解析可能であることを示したこと。
さらに、この論文はその後の二輪車シミュレーション解析手法の転換点となっている。
というのも、Sharpはこの論文において、それまで主流であった「ニュートン力学に基づき複雑な力を一つひとつベクトルで釣り合わせる」というマニュアル的な解析手法から脱却し、ラグランジュの方程式を用いた「エネルギーと座標」から機械的に運動方程式を導く手法を導入している。
この手法は「複雑な機械系を、自由度の数に応じた整然とした数式として体系化する手法」という数学的フレームワーク（マルチボディダイナミクス（MBD）の概念）が二輪車のシミュレーションに適用できることを示しており、加え、この定式化はコンピュータによる自動処理と非常に親和性が高かったことから、これを起点に二輪車の解析の世界においてMBDを使用した実用的なモデルが発展した。
つまり、現代の二輪車のシミュレーションの基礎理論（概念）の起点となる論文となっている。
マルチボディダイナミクス（MBD：Multibody Dynamics）：1950年代に宇宙工学から始まった多自由度システムの数値計算手法。複数の剛体や弾性体がジョイント（関節やヒンジ）で結合された機械システム全体の挙動をシミュレーションするための汎用的な力学解析手法
さて、詳細は別途として、本論文の挙動解析の大枠の流れをまとめてみれば以下。
1. 解析の流れ：数値モデルの作成
1. 運動方程式の算出：車両モデルの運動方程式（ラグランジェ方程式の利用）
運動エネルギー、位置エネルギー、一般化力の数式化：ラグランジェ方程式の利用の準備
４自由度モデルの各運動エネルギー（kinetic energy） \( T \)、位置エネルギー（potential energy） \( V \)、一般化力（generalized force） \( Q_i \) をそれぞれ数式化
車両モデルの運動方程式：ラグランジェ方程式の利用
ラグランジェ方程式
\( \underbrace{\dfrac{d}{dt}(\dfrac{\partial T}{\partial \dot {q_i}}) -(\dfrac{\partial T}{\partial q_i})}_{運動エネルギー項}+\underbrace{(\dfrac{\partial V}{\partial q_i})}_{位置エネルギー項}=\underbrace{Q_i} _{一般化力項} \ (i = 1,2,\cdots ,n) \)・・・①
つまり、
\( \dfrac{d}{dt}(\dfrac{\partial T}{\partial \dot {q_i}}) -(\dfrac{\partial T}{\partial q_i})+(\dfrac{\partial V}{\partial q_i})-{Q_i} = 0 \ (i = 1,2,\cdots ,n) \)・・・①’
に、数式化した各運動エネルギー（kinetic energy） \( T \)、位置エネルギー（potential energy） \( V \)、一般化力（generalized force） \( Q_i \) を代入し、これを\( y（横方向）,\psi（ヨー方向）,\phi（ロール方向）,\delta（ステアリング方向）\) の４つの自由度にて偏微分、かつその時の各一般化力 \( Q_i \)とあわせることで、以下の計4本の運動方程式を導く（各係数は和訳参照）。
<横方向>
\( (M_f + M_r)(\ddot{y}_1 + \dot{x}_1\dot{\psi}) + M_f k \ddot{\psi} + (M_f j + M_r h)\ddot{\phi} + M_f e \ddot{\delta} – Y_f – Y_r = 0 \)
<ヨー方向>
\( \begin{align}
& M_f k \ddot{y}_1 + (M_f e k + I_{fz}\cos\varepsilon)\ddot{\delta} – \dfrac{i_{fy}}{R_f}\sin\varepsilon \cdot \dot{x}_1\dot{\delta} + \{M_f j k – C_{rxz} + (I_{fz}-I_{fx})\sin\varepsilon\cos\varepsilon\}\ddot{\phi} \\
& – \left(\dfrac{i_{fy}}{R_f} + \dfrac{i_{ry}+\lambda i}{R_r}\right)\dot{x}_1\dot{\phi} + (M_f k^2 + I_{rz} + I_{fx}\sin^2\varepsilon + I_{fz}\cos^2\varepsilon)\ddot{\psi} + M_f k \dot{x}_1\dot{\psi} – l Y_f + b Y_r = 0
\end{align} \)
<ロール方向>
\(\begin{align}
& (M_f j + M_r h)\ddot{y}_1 + (M_f ej + I_{fz}\sin\varepsilon)\ddot{\delta} + \frac{i_{fy}}{R_f}\cos\varepsilon \cdot \dot{x}_1\dot{\delta} + (t Z_f – M_f eg)\delta \\
&+ (M_f j^2 + M_r h^2 + I_{rx} + I_{fx}\cos^2\varepsilon + I_{fz}\sin^2\varepsilon)\ddot{\phi} \\
&- (M_f j + M_r h)g\phi + \lbrace M_f j k – C_{rxz} + (I_{fz} – I_{fx})\sin\varepsilon\cos\varepsilon \rbrace \ddot{\psi} + \left(M_f j + M_r h + \frac{i_{fy}}{R_f} + \frac{i_{ry} + \lambda i}{R_r}\right)\dot{x}_1\dot{\psi} = 0
\end{align} \)
<ステアリング方向>
\(\begin{align} & M_f e \ddot{y}_1 + (I_{fz} + M_f e^2)\ddot{\delta} + K\dot{\delta} + (t Z_f – M_f eg)\sin\varepsilon \cdot \delta + (M_f ej + I_{fz}\sin\varepsilon)\ddot{\phi} – \frac{i_{fy}}{R_f}\cos\varepsilon \cdot \dot{x}_1\dot{\delta} \\
&+ (t Z_f – M_f eg)\phi + (M_f ek + I_{fz}\cos\varepsilon)\ddot{\psi} + \left(M_f e + \frac{i_{fy}}{R_f}\sin\varepsilon\right)\dot{x}_1\dot{\psi} + t Y_f = \tau \end{align} \)
2. 運動方程式の算出：タイヤモデルの運動方程式
（当時は緩和長の取込みが突破口だったみたいなので、ちょっとだけ詳細に記載）
フロントタイヤとリアタイヤの横力（side force）\(（Y_f, Y_r）\)
定常状態のタイヤの横力\(（Y’_f, Y’_r）\)
まず、定常状態におけるタイヤの横力は、横滑り角 \( (\alpha_f) \) とキャンバー角\( (\\phi_f) \)が微小として、コーナリング剛性\( (C_1) \)とキャンバー剛性\((C_2)\)用いた線形モデルで導く。
(「コーナリング剛性\( (N/rad) \) \(\times\) 横滑り角\( (rad) \)」と「キャンバー剛性\( (N/rad) \) \(\times\) キャンバー角\( (rad) \)」を足し合わせたもの)
つまり
- 定常状態のフロントタイヤ横力：\( \underbrace{Y’_f}_{タイヤ横力} = \underbrace{C_{f1}}_{コーナリング剛性} \cdot \underbrace{\alpha_f}_{横滑り角} + \underbrace{C_{f2} }_{キャンバー剛性}\cdot \underbrace{\phi_f}_{キャンバー角} \) ・・・②
- 定常状態のリアタイヤ横力 ：\( Y’_r = C_{r1} \cdot \alpha_r + C_{r2} \cdot \phi \) ・・・③
(リアタイヤのキャンバー角＝リアフレームのロール角 \(\phi \))
過渡特性の取り込み：タイヤの横力 (\(Y_f, Y_r\)) への緩和長 \(（\sigma）\) の取り込み
タイヤのゴムの弾性変形により、操舵による接地点のタイヤの横力は瞬間的には立ち上がらない。
この発生遅れの取込みのため、緩和長 (relaxation length:\(\sigma\))を使った1階の微分方程式を導入する。
- 過渡特性込みのフロントタイヤ横力：\( \underbrace{\dfrac{\sigma_f}{\dot{x}_1} \cdot \dot{Y}_f}_{遅れ発生による横力} + \underbrace{Y_f}_{実横力} = \underbrace{Y’_f}_{タイヤ横力} \) ・・・④
- 過渡特性込みのリアタイヤ横力 ：\( \dfrac{\sigma_r}{\dot{x}_1}\cdot \dot{Y}_r + Y_r = Y’_r \) ・・・⑤
実横力 (\(Y_f\)) に遅れ発生による横力（\( \dfrac{\sigma_f}{\dot{x}_1}\dot{Y}_f \)）が加わったものが、上述の定常横力 (\(Y’_f\)) となるタイヤモデルにすることにより、横力の発生遅れを取り込んだタイヤモデルとしている。
②③式に④⑤式を取り込めば、使用するタイヤモデルの運動方程式（フロントとリアの計２本）となる
- \( \dfrac{\sigma_f}{\dot{x}_1} \cdot \dot{Y}_f+ Y_f =C_{f1}\cdot \alpha_f + C_{f2}\cdot \phi_f \) ・・・④’
- \( \dfrac{\sigma_r}{\dot{x}_1}\cdot \dot{Y}_r + Y_r = C_{r1} \cdot \alpha_r + C_{r2} \cdot \phi \) ・・・⑤’
3. 運動方程式の型の一般化
さて、ラグランジェ方程式から導かれる運動方程式４本とタイヤの横力から導かれる運動方程式２本を、後の行列処理のために一般化しておく。
これら計６本の各方向の運動方程式 (\(y\)、\(\psi\)、\(\phi\)、\(\delta\)、\(Y_f\)、\(Y_r\)) は、（不要な要素は”０”を用いる事を前提とすれば）
運動方程式の標準的な型
\( \begin{align}
&M_{ij}\ddot{y}+M_{ij}\ddot{\psi}+M_{ij}\ddot{\phi}+M_{ij}\ddot{\delta}+M_{ij}\ddot{Y_f}+M_{ij}\ddot{Y_r}\\[6pt]
+&C_{ij}\dot{y}+C_{ij}\dot{\psi}+C_{ij}\dot{\phi}+C_{ij}\dot{\delta}+C_{ij}\dot{Y_f}+C_{ij}\dot{Y_r}\\[6pt]
+&K_{ij}{y}+K_{ij}{\psi}+K_{ij}{\phi}+K_{ij}{\delta}+K_{ij}{Y_f}+K_{ij}{Y_r}\\[6pt]
=&0 \ \ (i,j=1,2,\cdots , 6)
\end{align}\)
にて書くことが可能。
つまり、以下となる（各方向の運動方程式の型はすべて同じ形）
- \( M_{11}\ddot{y}+\cdots+M_{16}\ddot{Y_r}+C_{11}\dot{y}+\cdots+C_{16}\dot{Y_r}+K_{11}{y}+\cdots+K_{16}{Y_r}=0 \) ・・・⑥（横方向の運動方程式）
- \( M_{21}\ddot{y}+\cdots+M_{26}\ddot{Y_r}+C_{21}\dot{y}+\cdots+C_{26}\dot{Y_r}+K_{21}{y}+\cdots+K_{26}{Y_r}=0 \) ・・・⑦（ヨー方向の運動方程式）
- \( M_{31}\ddot{y}+\cdots+M_{36}\ddot{Y_r}+C_{31}\dot{y}+\cdots+C_{36}\dot{Y_r}+K_{31}{y}+\cdots+K_{36}{Y_r}=0 \) ・・・⑧（ロール方向の運動方程式）
- \( M_{41}\ddot{y}+\cdots+M_{46}\ddot{Y_r}+C_{41}\dot{y}+\cdots+C_{46}\dot{Y_r}+K_{41}{y}+\cdots+K_{46}{Y_r}=0 \) ・・・⑨（ステアリング方向の運動方程式）
- \( M_{51}\ddot{y}+\cdots+M_{56}\ddot{Y_r}+C_{51}\dot{y}+\cdots+C_{56}\dot{Y_r}+K_{51}{y}+\cdots+K_{56}{Y_r}=0 \) ・・・⑩（フロントタイヤ廻りの運動方程式）
- \( \underbrace{M_{61}\ddot{y}+\cdots+M_{66}\ddot{Y_r}}_{\mathbb{M}のベース}+ \underbrace{C_{61}\dot{y}+\cdots+C_{66}\dot{Y_r}}_{\mathbb{C}のベース}+\underbrace{K_{61}{y}+\cdots+K_{66}{Y_r}}_{\mathbb{K}のベース}=0 \) ・・・⑪（リアタイヤ廻りの運動方程式）
4. 運動方程式の行列化（数値モデル化）
ここで、\(6ｘ6 \) 行列の
\( \mathbb{M}=\left( \begin{array}{cccc}
M_{11} & M_{12} & \ldots & M_{16} \\
M_{21} & M_{22} & \ldots & M_{26} \\
\vdots & \vdots & \ddots & \vdots \\
M_{61} & M_{62} & \ldots & M_{66}
\end{array}
\right) \)、
\( \mathbb{C}=\left( \begin{array}{cccc}
C_{11} & C_{12} & \ldots & C_{16} \\
C_{21} & C_{22} & \ldots & C_{26} \\
\vdots & \vdots & \ddots & \vdots \\
C_{61} & C_{62} & \ldots & C_{66}
\end{array}
\right) \)、
\( \mathbb{K}=\left( \begin{array}{cccc}
K_{11} & K_{12} & \ldots & K_{16} \\
K_{21} & K_{22} & \ldots & K_{26} \\
\vdots & \vdots & \ddots & \vdots \\
K_{61} & K_{62} & \ldots & K_{66}
\end{array}
\right) \)、
と、状態ベクトル（\(6ｘ1 \) ベクトル）にて
\( \bf{q}=\left( \begin{array}{c}
{y} \\
{\psi} \\
{\phi} \\
{\delta}\\
{Y_f}\\
{Y_r}\\
\end{array}
\right) \)
とすれば、⑥～⑪式 をまとめて、
\( \underbrace{\left( \begin{array}{cccc}
M_{11} & M_{12} & \ldots & M_{16} \\
M_{21} & M_{22} & \ldots & M_{26} \\
\vdots & \vdots & \ddots & \vdots \\
M_{61} & M_{62} & \ldots & M_{66}
\end{array}
\right) }_{\mathbb{M}}
\underbrace{ \left( \begin{array}{c}
\ddot{y} \\
\ddot{\psi} \\
\ddot{\phi} \\
\ddot{\delta}\\
\ddot{Y_f}\\
\ddot{Y_r}\\
\end{array}
\right)} _{\bf{\ddot{q}}}
+
\underbrace{ \left( \begin{array}{cccc}
C_{11} & C_{12} & \ldots & C_{16} \\
C_{21} & C_{22} & \ldots & C_{26} \\
\vdots & \vdots & \ddots & \vdots \\
C_{61} & C_{62} & \ldots & C_{66}
\end{array}
\right) }_{\mathbb{C}}
\underbrace{ \left( \begin{array}{c}
\dot{y} \\
\dot{\psi} \\
\dot{\phi} \\
\dot{\delta}\\
\dot{Y_f}\\
\dot{Y_r}\\
\end{array}
\right)} _{\bf{\dot{q}}}
+
\underbrace{\left( \begin{array}{cccc}
K_{11} & K_{12} & \ldots & K_{16} \\
K_{21} & K_{22} & \ldots & K_{26} \\
\vdots & \vdots & \ddots & \vdots \\
K_{61} & K_{62} & \ldots & K_{66}
\end{array}
\right) }_{\mathbb{K}}
\underbrace{\left( \begin{array}{c}
{y} \\
{\psi} \\
{\phi} \\
{\delta}\\
{Y_f}\\
{Y_r}\\
\end{array}
\right)}_{\bf{q}} =0
\)
つまり、（標準化した形にしているので、当たり前ではあるが、）
\( \mathbb{M} \cdot \bf{\ddot{q}}+ \mathbb{C} \cdot \bf{\dot{q}} + \mathbb{K} \cdot \bf{q} = \rm{0} \) ・・・⑫
となる。
⑫式は、車両モデルの諸元をもとに算出したエネルギーを使用し、ラグランジェ方程式により得た運動方程式を行列化したもの（\( \mathbb{M} 、\mathbb{C} 、 \mathbb{K} \) と \(q\) を使用して行列にて一括表示したもの。）
つまり、⑫式がモーターサイクルの４自由度の「数値モデル」そのものであり、この式を使って振動モード解析（固有モード解析）を行っているのが本論文。
さて、振動モード解析へ。
標準モデルを使用していることから、一般的な振動モード解析手法もそのまま使える。
（今は論文が書かれた手計算の時代（約50年前）ではないため、以下はコンピュータ（計算ソフト）を使用を前提とした数値処理アルゴリズムによる解析手法も記載）
2. 解析の流れ：振動解析
振動解析概要
\( \mathbb{M} \cdot \bf{\ddot{q}}+ \mathbb{C} \cdot \bf{\dot{q}} + \mathbb{K} \cdot \bf{q} = \rm{0} \) （=⑫式）の \(\bf{q}\) の一般解を、
\(\bf{q}={\left( \begin{array}{c}
y_0 \cdot e^{\mu t}\\
\psi_0 \cdot e^{\mu t} \\
\phi_0 \cdot e^{\mu t} \\
\delta_0 \cdot e^{\mu t} \\
Y_{f0} \cdot e^{\mu t} \\
Y_{r0} \cdot e^{\mu t} \\
\end{array}
\right)} \) ・・・⑬
とすれば、
\( \begin{align}
&\mathbb{M} \cdot \bf{\ddot{q}}+ \mathbb{C} \cdot \bf{\dot{q}} + \mathbb{K} \cdot \bf{q} \\[6pt]
&=( \mathbb{M} \cdot \mu^2+ \mathbb{C} \cdot \mu + \mathbb{K}) \cdot \bf{q} = \rm{0} \end{align} \)
を満たす固有値 \(\mu \) （複素数）が「システムの特性根」となる。
固有値の算出について
実際の固有値算出は、車両諸元から定義した\( \mathbb{M},\mathbb{C},\mathbb{K}\)を使って、数理・数値計算ソフト（Matlab等）を使用すれば、ややこしい計算せずとも電卓レベルで算出可能。
（昔作ったMatlabのコードを見たら、（\( \mathbb{M},\mathbb{C},\mathbb{K}\)を入力し、車両速度を変えてループさせてはいるが、）実際の固有値算出は [X,egn]=polyeig(K,C,M) の一行のみ）
また今どき車両諸元も３Dモデルで設計してあれば \( \mathbb{M},\mathbb{C},\mathbb{K}\) も自動で算出できることから、シミュレーション結果まで一気通貫で実効可能。
←冒頭でも触れたMBDが ”コンピュータによる自動処理と非常に親和性が高い” の意味の所以。
固有値の解法
固有値算出については、いまどき電卓レベルで可能とはいうものの、その解法も覚書として追記。
（固有値算出のノートのかけらがあったので何処かに消えていく前に。。。）
解法１：特性方程式展開からの固有値算出概要
自分が昔使用したMatlabの [X,egn]=polyeig(K,C,M)（固有値解析コード X：固有ベクトル、egn：固有値、K,C,M：行列）は、この特性方程式から固有値の解を求めるコード 。（確か。。何の計算かが直感的にわかりやすかったので使った記憶あり）
さて、\(\bf{q}\) の一般解を、
\(\bf{q}={\left( \begin{array}{c}
y_0 \cdot e^{\mu t}\\
\psi_0 \cdot e^{\mu t} \\
\phi_0 \cdot e^{\mu t} \\
\delta_0 \cdot e^{\mu t} \\
Y_{f0} \cdot e^{\mu t} \\
Y_{r0} \cdot e^{\mu t} \\
\end{array}
\right)} \) ・・・⑬
とすれば、その一階微分 \(\bf{\dot{q}}\) は
\(\bf{\dot{q}}=\dfrac{d}{dt}{q}=\mu {\left( \begin{array}{c}
y_0 \cdot e^{\mu t}\\
\psi_0 \cdot e^{\mu t} \\
\phi_0 \cdot e^{\mu t} \\
\delta_0 \cdot e^{\mu t} \\
Y_{f0} \cdot e^{\mu t} \\
Y_{r0} \cdot e^{\mu t} \\
\end{array}
\right)}=\mu \cdot \bf{q}\) ・・・⑭
であり、二階微分 \(\bf{\ddot{q}}\) は
\(\bf{\ddot{q}}=\mu \cdot \bf{\dot{q}}=\mu^2 \cdot \bf{q}\) ・・・⑮
である。これから、
\( \mathbb{M} \cdot \bf{\ddot{q}}+ \mathbb{C} \cdot \bf{\dot{q}} + \mathbb{K} \cdot \bf{q} =( \mathbb{M} \cdot \mu^2+ \mathbb{C} \cdot \mu + \mathbb{K}) \cdot \bf{q} = \rm{0} \)
となることから、特性方程式
\(\mathbb{M} \cdot \mu^2+ \mathbb{C} \cdot \mu + \mathbb{K}= \rm{0} \) ・・・⑯
これを解けば、固有値 \(\mu\) を求める事ができる。
この解法は、モデルの自由度があがると（行列のサイズがあがると）行列式の展開計算が非常に煩雑になり（手計算では事実上不可能）非効率、数値丸めによる誤差も大きく安定性にもかける。これから、大規模モデルではあまり使われない（とのこと）。
解法2：行列展開からの固有値算出概要
行列展開による固有値算出法アルゴリズムとしては、QR法、QZ法が代表的。
概要としては、
- QR法（対象：通常の固有値問題「Standard Eigenvalue Problem」の解法）
：対象行列の型 \(AX = \mu X \) ・・・⑰ - QZ法（対象：一般化固有値問題「Generalized Eigenvalue Problem」の解法）
：対象行列の型 \((\mu B + A)X = 0 \) ・・・⑱（二つの行列 \(A\), \(B \) を同時に処理）
つまり、⑫式の \( \mathbb{M} \cdot \bf{\ddot{q}}+ \mathbb{C} \cdot \bf{\dot{q}} + \mathbb{K} \cdot \bf{q} = \bf{0} \) を
- \(AX = \mu X \) か \((\mu B + A)X = 0 \)
のどちらかの形に落とし込み、 計算ソフトの固有値解析コードに、型にあうように展開した行列\(A\)（と \(B \)） を落とし込めば、ソフトに組み込まれたアルゴリズムの利用により固有値の算出が可能。
QR法の利用
（Sharpの論文において使用されているのはQR法）
（解法１と同様に） \(\bf{q}\) の一般解を、
\(\bf{q}={\left( \begin{array}{c}
y_0 \cdot e^{\mu t}\\
\psi_0 \cdot e^{\mu t} \\
\phi_0 \cdot e^{\mu t} \\
\delta_0 \cdot e^{\mu t} \\
Y_{f0} \cdot e^{\mu t} \\
Y_{r0} \cdot e^{\mu t} \\
\end{array}
\right)} \) (=⑬式）
とすれば、\(\bf{\dot{q}}=\mu \cdot \bf{q}\)（=⑭式）、\(\bf{\ddot{q}}=\mu \cdot \bf{\dot{q}} \) （=⑮式）。
ここで、状態ベクトルとして \( X={\left( \begin{array}{c}
{\bf{q}} \\
{\bf{\dot{q}} } \\
\end{array}
\right)} \) を定義し、⑭⑮式を使えばその微分\( \dot{X}\)は、\( \dot{X}={\left( \begin{array}{c}
{\bf{\dot{q}}} \\
{\bf{\ddot{q}} } \\
\end{array}
\right)} =\mu {\left( \begin{array}{c}
{\bf{q}} \\
{\bf{\dot{q}} } \\
\end{array}
\right)} =\mu \cdot X \) となる。
さて、⑫式の \( \mathbb{M} \cdot \bf{\ddot{q}}+ \mathbb{C} \cdot \bf{\dot{q}} + \mathbb{K} \cdot \bf{q} = \rm{0} \) と \( \bf{I}\cdot \bf{\dot{q}}-\bf{I}\cdot \bf{\dot{q}}=\bf{0} \) (\(\bf{0} : \rm{0}行列、\bf{I}\) : 単位行列 (6×6) ) を使い、
再度行列化すれば
\(\begin{pmatrix} \bf{I} & \bf{0} \\ \bf{0} & \mathbb{M} \end{pmatrix} \begin{pmatrix} \bf{\dot{q}} \\ \bf{\ddot{q}} \end{pmatrix} + \begin{pmatrix} \bf{0} & -\bf{I} \\ \mathbb{K} & \mathbb{C} \end{pmatrix} \begin{pmatrix} \bf{q} \\ \bf{\dot{q}} \end{pmatrix} = \bf{0} ・・・⑲ \)
ここで、
\(\mathbb{B} = \begin{pmatrix} \bf{I} & \bf{0} \\ \bf{0} & \mathbb{M} \end{pmatrix} \)、\(\mathbb{A} = \begin{pmatrix} \bf{0} & -\bf{I} \\ \mathbb{K} & \mathbb{C} \end{pmatrix} \) とし、状態ベクトル \( X={\left( \begin{array}{c}
{\bf{q}} \\
{\bf{\dot{q}} } \\
\end{array}
\right)} \) をあわせれば、⑲式は
\(\mathbb{B}\cdot \dot{X}+\mathbb{A} \cdot X=\mathbb{B}\cdot \mu X+\mathbb{A} \cdot X=\bf{0} ・・・⑲’ \)
となる。これを変形して
\(\mathbb{A} \cdot X=- \mathbb{B}\cdot \mu X \)
\( (- \mathbb{B}^{-1}\cdot \mathbb{A} )\cdot X=\mu X \) ・・・⑳
となり、これはQR法（⑰）の標準形。つまり落とし込み完了。
ここで、
\(\mathbb{B} \cdot \mathbb{B}^{-1} = \begin{pmatrix} \bf{I} & \bf{0} \\ \bf{0} & \mathbb{M} \end{pmatrix} \begin{pmatrix} \bf{I} & \bf{0} \\ \bf{0} & \mathbb{M}^{-1} \end{pmatrix} =\bf{I} \) から \( \mathbb{B}^{-1} = \begin{pmatrix} \bf{I} & \bf{0} \\ \bf{0} & \mathbb{M}^{-1} \end{pmatrix} \)
である事を踏まえ、\(\mathbb{M} \)、\(\mathbb{C} \)、\(\mathbb{K} \) を使って⑳式を書き直せば、
\(\begin{align}
(- \mathbb{B}^{-1}\cdot \mathbb{A} )\cdot X &=- \begin{pmatrix} \bf{I} & \bf{0} \\ \bf{0} & \mathbb{M}^{-1}\end{pmatrix} \begin{pmatrix}\bf{0} & -\bf{I} \\ \mathbb{K} & \mathbb{C} \end{pmatrix} X \\[6pt]
&={\left( \begin{array}{cc}
0 & \bf{I} \\
-\mathbb{M}^{-1}\mathbb{K}&-\mathbb{M}^{-1}\mathbb{C} \\
\end{array}
\right)} \cdot X = \mu X ・・・㉑
\end{align} \)
つまり、
MATLAB等の数値計算ソフトの固有値解析コード（QR法）に、（\(AX = \mu X \)型の \(A\)として）
\( A ={\left( \begin{array}{cc}
0 & \bf{I} \\
-\mathbb{M}^{-1}\mathbb{K}&-\mathbb{M}^{-1}\mathbb{C} \\
\end{array}
\right)} \)
を使用すればよい。
QZ法の利用
⑲’式の \(\mathbb{B}\cdot \dot{X}+\mathbb{A} \cdot X=\mathbb{B}\cdot \mu X+\mathbb{A} \cdot X=\bf{0} \) を式変形し、
\((\mu \cdot \mathbb{B}+\mathbb{A}) \cdot X=\bf{0} \) ・・・㉒
とすれば、QZ法（⑱）のアルゴリズムが使用できる。
つまり、
数値計算ソフトの固有値解析コード（QZ法）に⑲’式の、（\((\mu B + A)X = 0 \)型の \(A、B\)として）
\(B = \begin{pmatrix} \bf{I} & \bf{0} \\ \bf{0} & \mathbb{M} \end{pmatrix} \)、\(A = \begin{pmatrix} \bf{0} & -\bf{I} \\ \mathbb{K} & \mathbb{C} \end{pmatrix} \)
をそのまま使用すればよい（\( A=\mathbb{A}、B=\mathbb{B} \) )。
なお、見ての通りQZ法においては、QR法では必要な逆行列計算 ( この場合 \(\mathbb{M}^{-1} \)の計算 ) が不要となるため、まるめ誤差等の計算誤差が小さく、計算精度が安定するため、この解法の方が一般的（とのこと）。
固有値・固有ベクトルが示す事
算出された固有値を観察する事により、このモーターサイクルモデルの外乱に対する挙動の減衰や発散の性質が判定できる。
1. 固有値 \(\mu \) (Eigenvalues）が示すこと
算出される固有値を、\(\mu = \sigma ± i \cdot \omega \) （共役複素数のペア）とすれば、一般解をオイラーの公式も使って式変形すれば、
\( \begin{align}
\bf{x}=&x_0 \cdot e^{\mu t}=x_0 \cdot e^{(\sigma ± i \cdot \omega) t} = x_0 \cdot e^{\sigma t} \cdot e^{±i \cdot \omega t}= x_0 \cdot e^{\sigma t} \cdot (cos \omega t ± i \cdot sin\omega t) \\[6pt]
=& x_0 \cdot e^{\sigma t} \cdot cos \omega t ± i \cdot x_0 \cdot e^{\sigma t} \cdot sin\omega t
\end{align}\)
これから、\(\mu = \sigma ± i \cdot \omega \) の
- 固有値 \(\mu \) の実部 (\(\sigma : \text{Real} 部分 \)) が、
- 正（\sigma>0）の場合：振幅( \(x_0 \cdot e^{\sigma t} \)) は時間とともに増大、つまり振動系モデルとして発散（不安定）
- 負（\sigma<0）の場合：振幅( \(x_0 \cdot e^{\sigma t} \)) は時間とともに減少、つまり振動系モデルとして収束（安定）
- 固有値 \(\mu \) の虚部 (\(\omega: \text{Imaginary} 部分\)) が、
- ゼロ (\( \omega = 0\)) の場合 ： 非振動モード（発散または減衰のみ）
「キャップサイズモード」 - 非ゼロ (\( \omega \neq 0\)) の場合 ：一定の周波数 \(\omega\) をもつ振動モード（複素共役根）
この場合に実部を考慮すれば、- 「固有値の実部が負 (\(\sigma < 0\))」と組み合わせされば「減衰しながら揺れが収まる 」
(例：収束するウィーブ・ウォブル) - 「固有値の実部が正 (\(\sigma > 0\))」と組み合わせされば「発散しながら揺れが大きくなる 」
(例：発散するウィーブ・ウォブル)
- 「固有値の実部が負 (\(\sigma < 0\))」と組み合わせされば「減衰しながら揺れが収まる 」
- ゼロ (\( \omega = 0\)) の場合 ： 非振動モード（発散または減衰のみ）
となる。
2. 固有ベクトル (Eigenvectors) が示すこと
- 各固有値 \(\mu \) に対応するモード形状（Mode Shape）が、\(\bf{q}\) に列ベクトルとして並ぶ。
- ここには、ある特定の固有値（例えばウォブルモードやウィーブモード）の周波数・減衰でシステムが振動・変動しているとき、「横すべり (\(\dot{y}_1\))」「ヨー (\(\psi\))」「ロール (\(\phi\))」「ステア (\(\delta\))」といった各自由度の変位や位相が、どのような比率（振幅・位相関係）で連動しているかが格納される。
さいごに
ざっくり言えば、４自由度モデルのモーターサイクルのモデルから得られた固有値、固有ベクトルにより、システムとしての挙動が「どのようなメカニズム（車体の傾きが主体の運動なのか、ステアの切れ角が主体の運動なのかなど）」に起因しているのかを特定・評価できる、としているのが本論文。
実際の解析手法自体は、動力学（Dynamics）の教科書であれば、例題込みで掲載されている解法なので、詳細はそちらを。
さいごのさいごに
なおこの解析にて、「発散」状況がそのモードと速度と共に結果として算出されるため、「発散」状況が着目されるが、本論文の面白いところは「収束」側にもある。
というのも、車両の外乱への収束がはやい事が「安定性が良い」と評価されるのであれば、この収束度合い（収束速度）は、その固有値の実部の負の大きさに現れる。
つまり、（固有値の実部が十分に「正」であり、振動が「発散」方向である場合はもちろんであるが）固有値の実部が０に近い値を示す場合においても、外乱の振幅に対する減衰力が低く（振幅が減衰せずに継続する）、それに伴い収束時間が長くなることになる、つまり「車両は不安定」と評価される。
逆に、固有値の実部が十分に「負」の値をとる速度域においては、発生する振動への減衰力は高く「車両は安定」と評価される。
別の見方をすれば、操舵軸に自由度をもつ二輪車が直進・直立するのは、固有値の実部が負の値を示す領域（速度域）であるが故に、車両モデルとして外乱（振動）を常時十分減衰する結果、と見ることもできる。
本論文のモデル作成のための車両諸元のベクトル計算、ラグランジェ方程式、固有値解析、等々はまた別途覚書化予定。。。