あらすじ
現代制御やLQRの勉強をしていると、MATLABで書かれたソースの中に
「expm」なる見慣れない関数に出会う。
上記の公式の解説を見てはみるものの、詳しいことは記載されていない。
行列指数なるものは一体何なのか？
そもそも、これが分かると何が嬉しいのか？
その辺りに関して記載していく。
行列指数とは
$A$ を$n$次正則行列として連立微分方程式を $\dot{x}(t) = Ax(t)$ で表す時、
一般解 $x(t) = e^{At} x(0)$ のうち、$e^{At}$ を行列指数と呼ぶようだ。
どうも調べてみると、「行列指数」というより「行列(の)指数関数」という
言い方が一般的なようだ。
関数expmは、この行列指数関数 $e^{At}$ を算出してくれるらしい。
つまり連立微分方程式の解が簡単に求めることができる。
ということは、状態空間方程式 $\dot{x}(t)=Ax(t) + Bu(t)$ の解も
求めることができることになる。
ホント？
どうせなので試してみる。やり方としては以下ステップを踏む。
① 連立微分方程式を手計算で解き、行列指数関数の部分を算出
② MATLABで関数expmを使って行列指数関数を算出
③ ①の結果と②の結果を突き合わせる
連立微分方程式の解き方
実際に連立微分方程式を解く前に、まずは計算方法について記述する。
以下の連立微分方程式について考える。
\left\{
\begin{array}{}
\dfrac{dx_1(t)}{dt} = ax_1(t) + bx_2(t)\\
\dfrac{dx_2(t)}{dt} = cx_1(t) + dx_2(t)
\end{array}
\right.\\
$x(t)$ と $A$ をそれぞれ
x(t) = \begin{pmatrix}
x_{1}(t) \\
x_{2}(t) \\
\end{pmatrix}
、A = \begin{pmatrix}
a & b \\
c & d \\
\end{pmatrix}
と置くと、 $\dot{x}(t) = \dfrac{dx(t)}{dt} = Ax(t)$ と表すことができる。
ここで、$A$ の固有値 $λ$ が $n$ 個存在する場合、以下の式にて対角化が可能。
（ $A$ は $2×2$ なので、この場合 $n=2$ ）
P^{-1}AP = D = \begin{pmatrix}
\lambda_1 & 0 \\
0 & \lambda_2 \\
\end{pmatrix}
なお、$P$ は $A$ の固有ベクトルから構成される正則行列とする。
対角化が可能な場合、$\dfrac{dx(t)}{dt} = Ax(t)$ を以下のように変形できる。
\begin{align}
\dfrac{dx(t)}{dt} &= Ax(t)\\
P^{-1}\dfrac{dx(t)}{dt} &= P^{-1}APP^{-1}x(t)\\
\dfrac{d}{dt}(P^{-1}x(t)) &= DP^{-1}x(t)\\
\end{align}
※ AA^{-1} = A^{-1}A = E (単位行列)\\
AE = EA = A
この時、$y(t) = P^{-1}x(t)$ と置くと、↑の式は $\dfrac{d}{dt}y(t) = \dot{y}(t) = Dy(t)$ とすることができ、更に以下のように変形する。
\dfrac{d}{dt}
\begin{pmatrix}
y_{1}(t) \\
y_{2}(t) \\
\end{pmatrix}
=
\begin{pmatrix}
\lambda_{1} & 0\\
0 & \lambda_{2} \\
\end{pmatrix}
\begin{pmatrix}
y_{1}(t) \\
y_{2}(t) \\
\end{pmatrix}\\
\\
\left\{
\begin{array}{}
\dfrac{d}{dt}y_1(t) = \lambda_1 y_1(t)\\
\dfrac{d}{dt}y_2(t) = \lambda_2 y_2(t)
\end{array}
\right.\\
算出した連立微分方程式の一般解を、それぞれ求める。
\begin{align}
\dfrac{d}{dt}y_1(t) &= \lambda_1 y_1(t)\\
dy_1(t) &= \lambda_1 y_1(t)dt\\
\dfrac{1}{y_1(t)} dy_1(t) &= \lambda_1 dt\\
\int \dfrac{1}{y_1(t)} dy_1(t) &= \int \lambda_1 dt\\
\log y_1(t) + C_1 &= \lambda_1t + C_2\\
\end{align}\\
$C_2 - C_1 = C$ と置き
\begin{align}
\log y_1(t) &= \lambda_1t + C\\
y_1(t) &= e^{\lambda_1t + C}\\
&= e^C e^{\lambda_1t}\\
\end{align}
ここで $t = 0$ とすると $y_1(0) = e^C e^{\lambda_1 * 0} = e^C$
これより $e^C = y_1(0)$ と置けるため
\therefore y_1(t) = y_1(0) e^{\lambda_1t}\\
$y_2(t)$ の式も同様に一般解を求めると、以下のようになる。
\left\{
\begin{array}{}
y_1(t) = y_1(0) e^{\lambda_1t}\\
y_2(t) = y_2(0) e^{\lambda_2t}\\
\end{array}
\right.\\
また、$y(t) = P^{-1}x(t)$ の両辺に $P$ をかけて $x(t) = Py(t)$ に変形する。
\begin{align}
x(t)&=Py(t)\\
&=P
\begin{pmatrix}
y_1(t)\\
y_2(t)\\
\end{pmatrix}\\
&=P
\begin{pmatrix}
y_1(0) e^{\lambda_1t}\\
y_2(0) e^{\lambda_2t}\\
\end{pmatrix}\\
&=P
\begin{pmatrix}
e^{\lambda_1t} & 0\\
0 & e^{\lambda_2t}\\
\end{pmatrix}
\begin{pmatrix}
y_1(0)\\
y_2(0)\\
\end{pmatrix}\\
&=P
\begin{pmatrix}
e^{\lambda_1t} & 0\\
0 & e^{\lambda_2t}\\
\end{pmatrix}y(0)
\end{align}
ここで、$y(t) = P^{-1}x(t)$ より $y(0) = P^{-1}x(0)$ を代入する。
\begin{align}
x(t) &= P
\begin{pmatrix}
e^{\lambda_1t} & 0\\
0 & e^{\lambda_2t}\\
\end{pmatrix}
P^{-1}
\begin{pmatrix}
x_1(0)\\
x_2(0)
\end{pmatrix}\\
&= P
\begin{pmatrix}
e^{\lambda_1t} & 0\\
0 & e^{\lambda_2t}\\
\end{pmatrix}
P^{-1}x(0)
\end{align}
以上の式を、連立微分方程式の一般解とする。
① 手計算で解く
以下の連立微分方程式を解く。
\left\{
\begin{array}{}
\dfrac{dx_1(t)}{dt} = -7x_1(t) + 5x_2(t)\\
\dfrac{dx_2(t)}{dt} = 2x_1(t) + 2x_2(t)
\end{array}
\right.\\
まず、最初の式より
A =
\begin{pmatrix}
-7 & 5 \\
2 & 2 \\
\end{pmatrix}
として $|A-λE|$ を計算し、固有値 $λ$ を求める。
\begin{align}
|A-λE|=\begin{vmatrix}
-7-λ & 5 \\
2 & 2-λ \\
\end{vmatrix}&=(-7-\lambda)(2-\lambda)-(5\cdot2)\\
&=(\lambda^2+5\lambda-14)-(10)\\
&=\lambda^2+5\lambda-24\\
&=(\lambda+8)(\lambda-3)\\
\\
\therefore\lambda=3, -8
\end{align}
次に、固有ベクトルを求める。
固有値と固有ベクトルの関係から $Ap_1 = λ_1p_1$
及び $Ap_2 = λ_2p_2$ が成り立つ。
まず、固有値の1つ $λ_1 = 3$ の固有ベクトルを
p_1 = \begin{pmatrix}
p_{11} \\
p_{21} \\
\end{pmatrix}
として、以下を解く。
A\vec{p_1} =
\begin{pmatrix}
-7 & 5 \\
2 & 2 \\
\end{pmatrix}
\begin{pmatrix}
p_{11} \\
p_{21} \\
\end{pmatrix}
=3\begin{pmatrix}
p_{11} \\
p_{21} \\
\end{pmatrix}
⇒
\left\{
\begin{array}{}
-7p_{11}+5p_{21}=3p_{11}\\
2p_{11}+2p_{21}=3p_{21}
\end{array}\\
\right.\\
計算すると $2p_{11}=p_{21}$ が求まるので
\therefore p_1=
\begin{pmatrix}
p_{11} \\
p_{21} \\
\end{pmatrix}
=
\begin{pmatrix}
1 \\
2 \\
\end{pmatrix}
とする。
続いて、もう一つの固有値 $λ_2=-8$ の時の固有ベクトルも同じ要領で算出する。
A\vec{p_2} =
\begin{pmatrix}
-7 & 5 \\
2 & 2 \\
\end{pmatrix}
\begin{pmatrix}
p_{12} \\
p_{22} \\
\end{pmatrix}
=-8\begin{pmatrix}
p_{12} \\
p_{22} \\
\end{pmatrix}
⇒
\left\{
\begin{array}{}
-7p_{11}+5p_{21}=-8p_{11}\\
2p_{11}+2p_{21}=-8p_{21}
\end{array}\\
\right.\\
p_{12}=-5p_{22}
\therefore p_2=
\begin{pmatrix}
p_{12} \\
p_{22} \\
\end{pmatrix}
=
\begin{pmatrix}
5 \\
-1 \\
\end{pmatrix}
以上より、固有ベクトル $p$ から構成される正則行列 $P$ が求まった。
\therefore P=
\begin{pmatrix}
p_{11} & p_{12} \\
p_{21} & p_{22} \\
\end{pmatrix}
=
\begin{pmatrix}
1 & 5 \\
2 & -1 \\
\end{pmatrix}
固有値と固有ベクトルが求まったので、 最後に $x(t)=P\begin{pmatrix}e^{\lambda_1t} & 0 \newline 0 & e^{\lambda_2t} \end{pmatrix}P^{-1}x(0)$ を解く。
\begin{align}
P\begin{pmatrix}e^{\lambda_1t} & 0 \newline 0 & e^{\lambda_2t} \end{pmatrix}P^{-1}&=
\begin{pmatrix}
1 & 5 \\
2 & -1 \\
\end{pmatrix}
\begin{pmatrix}
e^{3t} & 0 \\
0 & e^{-8t} \\
\end{pmatrix}
-
\dfrac{1}{11}
\begin{pmatrix}
-1 & 5 \\
-2 & 1 \\
\end{pmatrix}x(0)\\
&=-\dfrac{1}{11}
\begin{pmatrix}
e^{3t} & 5e^{-8t} \\
2e^{3t} & -e^{-8t} \\
\end{pmatrix}
\begin{pmatrix}
-1 & 5 \\
-2 & 1 \\
\end{pmatrix}x(0)\\
&=-\dfrac{1}{11}
\begin{pmatrix}
-e^{3t}-10e^{-8t} & -5e^{3t} + 5e^{-8t} \\
-2e^{3t}+2e^{e^{-8t}t} & -10e^{3t}-e^{-8t} \\
\end{pmatrix}x(0)
\end{align}
② MATLABで解く
$A =\begin{pmatrix}-7 & 5 \newline 2 & 2 \end{pmatrix}$ のため、これをexpmに代入するだけ。
A = [-7 5; 2 2];
expm(A)
結果は、$\begin{pmatrix}1.8263 & 9.1296 \newline 3.6519 & 18.2596 \end{pmatrix}$ となった。
③ ①と②の結果を突き合わせる
①の結果：$-\dfrac{1}{11}\begin{pmatrix}-e^{3t}-10e^{-8t} & -5e^{3t} + 5e^{-8t} \newline
-2e^{3t}+2e^{e^{-8t}t} & -10e^{3t}-e^{-8t}\end{pmatrix}$
②の結果：$\begin{pmatrix}1.8263 & 9.1296 \newline 3.6519 & 18.2596 \end{pmatrix}$
①の結果に $t = 1$ を代入して計算すると $\begin{pmatrix}1.8263 & 9.1296 \newline 3.6519 & 18.2596 \end{pmatrix}$ となり、
①と②の結果が一致することが分かった。
参照ページ
本記事の作成にあたり、以下を参考にさせていただきました。
・https://dora.bk.tsukuba.ac.jp/~takeuchi/?%E7%B7%9A%E5%BD%A2%E4%BB%A3%E6%95%B0II%2F%E9%80%A3%E7%AB%8B%E7%B7%9A%E5%BD%A2%E5%BE%AE%E5%88%86%E6%96%B9%E7%A8%8B%E5%BC%8F
・https://hb3.seikyou.ne.jp/home/E-Yama/renDE.pdf
最後に
この記事は、行列指数という言葉に出会ってから、その言葉と関連する
理論について調べ、自分なりに理解して書いたものです。
もし本記事の内容が間違っていれば、遠慮なく指摘いただけると大変助かります。