状態方程式（物理じゃないよ）を解く

導入

　今回は状態方程式を解きます．状態方程式ってのは，熱力学で登場する$PV=nRT$ではなく
$$\dot{x}=Ax+Bu$$
です．これは制御工学で登場します．

準備

　厳密に定式化すると以下のようになります．

制御工学（特に現代制御）では，$\mathit{x}$を状態ベクトル，$\mathit{u}$を入力ベクトルと呼び，状態方程式からシステムの制御や解析を行います．たとえば，力学システムでは状態ベクトルを$x=(\text{位置},\text{速度}{)}^T$などとおきます．
　また遷移行列を定義します．

これは指数関数$e^x$の$x=0$まわりのTaylor展開：$e^x=\sum x^k/(k!)$に行列$At$を代入することで得られる形になっています．この遷移行列は，以下のような指数関数と同じ性質が成り立ちます．

導出

発見的導出

$B=O$のとき

　状態方程式は
$$\dot{\mathit{x}}=A\mathit{x}$$
となります．この微分方程式の解空間$W$は次元が1：$\dim W=1$の線形空間になります．したがって具体的に解$\mathit{\xi }(t)$を見つければ，$W$の一般解は$\mathit{x}(t)=C\mathit{\xi }(t)$となります（$C$は定数）．そこで$\mathit{x}=e^{At}\mathit{c}$（$\mathit{c}$は$\mathit{x}$と同形状の定ベクトル）としてみると，これは状態方程式を満たします．これで具体的な解：$\mathit{\xi }(t)=e^{At}\mathit{c}$が見つかったので，初期条件：$\mathit{x}(0)={\mathit{x}}_0(=C\mathit{c})$として
$$\mathit{x}(t)=e^{At}{\mathit{x}}_0$$
となります．

$B\ne O$のとき

　状態方程式は
$$\dot{\mathit{x}}=A\mathit{x}+B\mathit{u}$$
です．ここで${\displaystyle \mathit{x}={\int }_0^t{e}^{A(t-\tau )}B\mathit{u}(\tau )d\tau }$とおいてみます．すると
$$\begin{array}{rl}\dot{\mathit{x}}& =\frac{d}{dt}(e^{At}{\int }_0^t{e}^{-A\tau }B\mathit{u}(\tau )d\tau )\\ & =A{e}^{At}{\int }_0^t{e}^{-A\tau }B\mathit{u}(\tau )d\tau +e^{At}(e^{-At}B\mathit{u}(t))\\ & =A\mathit{x}+B\mathit{u}\end{array}$$
となって，状態方程式の特殊解であることがわかります．よって$\mathit{X}={\displaystyle \mathit{x}-{\int }_0^t{e}^{A(t-\tau )}B\mathit{u}(\tau )d\tau }$とおくと，$\mathit{X}$は
$$\dot{\mathit{X}}=A\mathit{X}$$
を満たします．これは前節$B=O$の場合と同じ微分方程式なので$\mathit{X}=e^{At}{\mathit{x}}_0$となります（ここで$\mathit{X}(0)=\mathit{x}(0)={\mathit{x}}_0$に注意）．よって，状態方程式の解は
$$\mathit{x}(t)=e^{At}{\mathit{x}}_0+{\int }_0^t{e}^{A(t-\tau )}B\mathit{u}(\tau )d\tau .$$

解析的導出

　状態方程式は
$$\dot{\mathit{x}}=A\mathit{x}+B\mathit{u}$$
です．両辺に$e^{-At}\in M_{n,n}(\mathbb{R})$を左から乗じると$e^{-At}\dot{\mathit{x}}=e^{-At}A\mathit{x}+e^{-At}B\mathit{u}$．定義から$e^{-At}A=A{e}^{-At}$で，定理１と積の微分法から
$$\frac{d}{dt}(e^{-At}\mathit{x})=e^{-At}B\mathit{u}.$$
両辺を積分すると
$$\begin{array}{rl}{\int }_0^t\frac{d}{d\tau }(e^{-A\tau }\mathit{x}(\tau ))d\tau & ={\int }_0^t{e}^{-A\tau }B\mathit{u}(\tau )d\tau \\ e^{-At}\mathit{x}(t)-e^{-A0}\mathit{x}(0)& ={\int }_0^t{e}^{-A\tau }B\mathit{u}(\tau )d\tau \\ e^{-At}\mathit{x}(t)& ={\mathit{x}}_0+{\int }_0^t{e}^{-A\tau }B\mathit{u}(\tau )d\tau \end{array}$$
両辺に$e^{At}\in M_{n,n}(\mathbb{R})$を左から乗じて
$$\mathit{x}(t)=e^{At}{\mathit{x}}_0+{\int }_0^t{e}^{A(t-\tau )}B\mathit{u}(\tau )d\tau .$$

具体例で確認

　導出した公式をまとめると以下のようになります．

これを具体的な例で確認してみます．

単振動系

入力なし($\mathit{u}=\mathbf{0}$)のとき

　水平で摩擦のない床に質量$m$の物体がばね(ばね定数：$k$)によって壁と接続されているとします．物体の位置をばねの伸び$X$で表すとき，状態ベクトルを$\mathit{x}=(\text{位置},\text{速度}{)}^T=(X,\dot{X}{)}^T$とおけば，状態方程式は
$$\dot{\mathit{x}}=A\mathit{x},\text{ただし，}A=\left(\begin{array}{cc}0& 1\\ -k/m& 0\end{array}\right).$$
初期条件を$\mathit{x}(0)=(X_0,v_0{)}^T(=:{\mathit{x}}_0)$とすれば，公式から$\mathit{x}=e^{At}{\mathit{x}}_0$．遷移行列$e^{At}$を計算すると，$\omega =\sqrt{\frac{k}{m}}$として
$$\begin{array}{rl}{e}^{At}& =\sum _{n=0}^{\mathrm{\infty }}\frac{t^{2n}}{(2n)!}{A}^{2n}+\sum _{n=0}^{\mathrm{\infty }}\frac{t^{2n+1}}{(2n+1)!}{A}^{2n+1}\\ & =\sum _{n=0}^{\mathrm{\infty }}\frac{t^{2n}}{(2n)!}{\left(\begin{array}{cc}-{\omega }^2& 0\\ 0& -{\omega }^2\end{array}\right)}^n+\sum _{n=0}^{\mathrm{\infty }}\frac{t^{2n+1}}{(2n+1)!}A{\left(\begin{array}{cc}-{\omega }^2& 0\\ 0& -{\omega }^2\end{array}\right)}^n\\ & =\sum _{n=0}^{\mathrm{\infty }}\frac{t^{2n}}{(2n)!}{(-{\omega }^2)}^{n}I+\sum _{n=0}^{\mathrm{\infty }}\frac{t^{2n+1}}{(2n+1)!}{(-{\omega }^2)}^{n}A\\ & =\left(\begin{array}{cc}\sum \frac{(-1{)}^n}{(2n)!}(\omega t{)}^{2n}& \frac{1}{\omega }\sum \frac{(-1{)}^n}{(2n+1)!}(\omega t{)}^{2n+1}\\ -\omega \sum \frac{(-1{)}^n}{(2n+1)!}(\omega t{)}^{2n+1}& \sum \frac{(-1{)}^n}{(2n)!}(\omega t{)}^{2n}\end{array}\right)\\ & =\left(\begin{array}{cc}\cos \omega t& \frac{1}{\omega }\sin \omega t\\ -\omega \sin \omega t& \cos \omega t\end{array}\right)\end{array}$$
よって
$$\mathit{x}=\left(\begin{array}{c}X\\ \dot{X}\end{array}\right)=e^{At}{\mathit{x}}_0=\left(\begin{array}{c}{X}_0\cos \omega t+{\displaystyle \frac{v_0}{\omega }}\sin \omega t\\ -\omega X_0\sin \omega t+v_0\cos \omega t\end{array}\right).$$

入力あり（強制振動）のとき

　先ほどの系に入力として外力$f=F\sin \mathrm{\Omega }t$を与える場合を考えます．すると運動方程式は$m\ddot{X}+kX=F\sin \mathrm{\Omega }t$ゆえ，状態ベクトル：$\mathit{x}=(X,\dot{X}{)}^T$，入力ベクトル：$\mathit{u}=f$とすれば状態方程式は
$$\dot{\mathit{x}}=A\mathit{x}+B\mathit{u},\text{ただし，}A=\left(\begin{array}{cc}0& 1\\ -k/m& 0\end{array}\right)=\left(\begin{array}{cc}0& 1\\ -{\omega }^2& 0\end{array}\right),B=\left(\begin{array}{c}0\\ 1/m\end{array}\right).$$
よって一般解は
$$\begin{array}{rl}\mathit{x}& =e^{At}{\mathit{x}}_0+{\int }_0^t{e}^{A(t-\tau )}B\mathit{u}(\tau )d\tau .\end{array}$$
第2項の積分は
$$\begin{array}{rl}{\int }_0^t{e}^{A(t-\tau )}B\mathit{u}(\tau )d\tau & ={\int }_0^t\left(\begin{array}{cc}\cos \omega (t-\tau )& \frac{1}{\omega }\sin \omega (t-\tau )\\ -\omega \sin \omega (t-\tau )& \cos \omega (t-\tau )\end{array}\right)\left(\begin{array}{c}0\\ \frac{F}{m}\sin \mathrm{\Omega }\tau \end{array}\right)d\tau \\ & =\frac{F}{m}{\int }_0^t{\left(\begin{array}{cc}\frac{1}{\omega }\sin \omega (t-\tau )\sin \mathrm{\Omega }\tau & \cos \omega (t-\tau )\sin \mathrm{\Omega }\tau \end{array}\right)}^{T}d\tau \\ & =\frac{F}{m}{\int }_0^t{\left(\frac{\cos ((\omega +\mathrm{\Omega })\tau -\omega t)-\cos ((\mathrm{\Omega }-\omega )\tau +\omega t)}{2\omega }\frac{\sin ((\omega +\mathrm{\Omega })\tau -\omega t)+\sin ((\mathrm{\Omega }-\omega )\tau +\omega t)}{2}\right)}^{T}d\tau \\ & =\frac{F}{2m}{(\frac{\sin \mathrm{\Omega }t+\sin \omega t}{\omega (\omega +\mathrm{\Omega })}-\frac{\sin \mathrm{\Omega }t-\sin \omega t}{\omega (\mathrm{\Omega }-\omega )}\frac{-\cos \mathrm{\Omega }t+\cos \omega t}{\omega +\mathrm{\Omega }}+\frac{-\cos \mathrm{\Omega }t+\cos \omega t}{\mathrm{\Omega }-\omega })}^T\end{array}$$
であるから，一般解は
$$\begin{array}{rl}\mathit{x}(t)& =\left(\begin{array}{c}{X}_0\cos \omega t+{\displaystyle \frac{v_0}{\omega }}\sin \omega t+{\displaystyle \frac{F}{2m}}({\displaystyle \frac{\sin \mathrm{\Omega }t+\sin \omega t}{\omega (\mathrm{\Omega }+\omega )}}-{\displaystyle \frac{\sin \mathrm{\Omega }t-\sin \omega t}{\omega (\mathrm{\Omega }-\omega )}})\\ -\omega X_0\sin \omega t+v_0\cos \omega t-{\displaystyle \frac{F}{2m}}({\displaystyle \frac{\cos \mathrm{\Omega }t-\cos \omega t}{\mathrm{\Omega }+\omega }}+{\displaystyle \frac{\cos \mathrm{\Omega }t-\cos \omega t}{\mathrm{\Omega }-\omega }})\end{array}\right)\\ & =\left(\begin{array}{c}{X}_0\cos \omega t+{\displaystyle \frac{v_0}{\omega }}\sin \omega t+{\displaystyle \frac{F}{m\omega ({\mathrm{\Omega }}^2-{\omega }^2)}}(\mathrm{\Omega }\sin \omega t-\omega \sin \mathrm{\Omega }t)\\ -\omega X_0\sin \omega t+v_0\cos \omega t+{\displaystyle \frac{F\mathrm{\Omega }}{m({\mathrm{\Omega }}^2-{\omega }^2)}}(\cos \omega t-\cos \mathrm{\Omega }t)\end{array}\right).\end{array}$$

ちなみに

　入力の角振動数$\mathrm{\Omega }$が固有振動数$\omega$と一致するときはどうなるか調べてみます．$\mathit{x}(t)$の第1・2項は収束が明らかですが，第3項は微妙です．そこで
$${\displaystyle \frac{F}{m\omega ({\mathrm{\Omega }}^2-{\omega }^2)}}(\mathrm{\Omega }\sin \omega t-\omega \sin \mathrm{\Omega }t)=\frac{F\mathrm{\Omega }}{m(\mathrm{\Omega }+\omega )}\times \frac{\frac{\sin \omega t}{\omega }-\frac{\sin \mathrm{\Omega }t}{\mathrm{\Omega }}}{\mathrm{\Omega }-\omega }$$
と変形し，後半$\frac{\sin \omega t/\omega -\sin \mathrm{\Omega }t/\mathrm{\Omega }}{\mathrm{\Omega }-\omega }$の収束値を調べてみます．この値の$\mathrm{\Omega }=\omega$の場合の収束値は，微分係数の定義から$g(x)=\sin (tx)/x(x>0)$とすれば$-g^{\prime }(\omega )={\displaystyle \frac{-t\cos \omega t\cdot \omega +\sin \omega t}{{\omega }^2}}$です．よって
$$\underset{\mathrm{\Omega }\to \omega }{lim}\mathit{x}(t)=X_0\cos \omega t+\frac{v_0}{\omega }\sin \omega t+\frac{F}{2m{\omega }^2}(\sin \omega t-\omega t\cos \omega t)$$
となり，最後が$t\cos \omega t$という時間発展に伴って振幅がどんどん増加する項になります．つまり入力の振動数が固有振動数と一致すると大きく振動するという結果が得られます．