漸近安定性

漸近安定性…Homework

[0] 次図に示すように軸支された剛体振り子は2つの静止した平衡状態(\alpha)と(\beta)を持ちます。

図1 剛体振り子

いまこの平衡状態から棒を少し傾けて手を離すと、明らかに違った振舞いをします。実際には軸回りの摩擦や空気抵抗による抗力があるので、平衡状態(\alpha)の方は元に戻りますが、平衡状態(\beta)の方は元に戻ることなく平衡状態(\alpha)に落ち着くでしょう。前者を漸近安定、後者を不安定とよびます。このような違いは平衡状態回りでの振舞いを表す線形状態方程式にどのように反映されているのか気になるところです。

この場合の運動方程式は次式となります。

\displaystyle{(1)\quad J\ddot{\theta}(t)=-mg\ell\sin\theta(t)-c\dot{\theta}(t)\quad(J=\frac{4}{3}m\ell^2) }

ここで、右辺の第2項が角速度に比例する抗力を表しています(c>0は定数)。これから平衡状態(\alpha)と(\beta)に応じて次の状態方程式を得ます。

\displaystyle{(2)\quad \underbrace{ \left[\begin{array}{l} \dot{\theta}(t) \\ \dot{\omega}(t) \end{array}\right] }_{\dot{x}(t)}= \underbrace{ \left[\begin{array}{cc} 0 & 1 \\ -\frac{3g}{4\ell}&-\frac{c}{J} \end{array}\right] }_{A=A_1} \underbrace{ \left[\begin{array}{l} \theta(t) \\ \omega(t) \end{array}\right] }_{x(t)} }

\displaystyle{(3)\quad \underbrace{ \frac{d}{dt} \left[\begin{array}{l} \theta(t)-\pi \\ \omega(t)-0 \end{array}\right] }_{\dot{x}(t)}= \underbrace{ \left[\begin{array}{cc} 0 & 1 \\ \frac{3g}{4\ell}&-\frac{c}{J} \end{array}\right] }_{A=A_2} \underbrace{ \left[\begin{array}{l} \theta(t)-\pi \\ \omega(t)-0 \end{array}\right] }_{x(t)} }

以下では、行列Aの固有値を調べて、漸近安定かどうかを判定できることを説明します。

[1] 制御対象が平衡状態にあることは線形状態方程式\dot{x}(t)=Ax(t)+Bu(t)において、x=0, u=0を意味します。そこで平衡状態が乱されてx\ne0となる時刻をt=0にとると、線形状態方程式は次式となります。

\displaystyle{(4)\quad \dot{x}(t)=Ax(t)\qquad(x(0)\ne0) }

漸近安定であることは平衡状態x=0に戻ることでしたから、次のように表すことができます。

\displaystyle{(5)\quad \boxed{x(t)\rightarrow 0\qquad(t\rightarrow\infty)} }

そのための条件を検討するため、微分方程式\dot{x}(t)=Ax(t)の解を求める必要があります。これは次式で表されます(Note A21-1参照)。

\displaystyle{(6)\quad \boxed{x(t)=\exp(At)x(0)} }

ここで、xn次元ベクトル、An次正方行列です。任意のn次正方行列Xに対して、行列指数関数\exp(X)

\displaystyle{(7)\quad \boxed{\exp(X)=I_n+X+\frac{1}{2}X^2+\cdots+\frac{1}{k!}X^k+\cdots} }

で定義されます。そこで、解の表現式(6)をn=2の場合について詳しく説明します。

[2] いま2次系に限定して、微分方程式\dot{x}(t)=Ax(t)を要素を明示して書くと

\displaystyle{(8)\quad \underbrace{ \left[\begin{array}{l} \dot{x}_1(t) \\ \dot{x}_2(t) \end{array}\right] }_{\dot{x}(t)}= \underbrace{ \left[\begin{array}{cc} a_{11} & a_{12} \\ a_{21} & a_{22} \\ \end{array}\right] }_{A} \underbrace{ \left[\begin{array}{l} x_1(t) \\ x_2(t) \end{array}\right] }_{x(t)} }

すなわち

\displaystyle{(9)\quad \left\{\begin{array}{l} \dot{x}_1(t)=a_{11}x_1(t)+a_{12}x_2(t) \\ \dot{x}_2(t)=a_{21}x_1(t)+a_{22}x_2(t) \end{array}\right. }

となります。もしAが対角行列で、a_{12}=a_{21}=0であれば

\displaystyle{(10)\quad \left\{\begin{array}{l} \dot{x}_1(t)=a_{11}x_1(t) \\ \dot{x}_2(t)=a_{22}x_2(t) \end{array}\right. }

となって、解は

\displaystyle{(11)\quad \left\{\begin{array}{l} x_1(t)=e^{a_{11}t}x_1(0) \\ x_2(t)=e^{a_{22}t}x_2(0) \\ \end{array}\right. }

すなわち

\displaystyle{(12)\quad \underbrace{ \left[\begin{array}{l} {x}_1(t) \\ {x}_2(t) \end{array}\right] }_{{x}(t)}= \underbrace{ \left[\begin{array}{cc} e^{a_{11}t} & 0 \\ 0 & e^{a_{22}t} \\ \end{array}\right] }_{\exp(At)} \underbrace{ \left[\begin{array}{l} x_1(0) \\ x_2(0) \end{array}\right] }_{x(0)} }

ですから、解の表現式(6)は納得できると思います。

[3] それではAが対角行列でない場合はどうするかですが、一対一対応の変数変換

\displaystyle{(13)\quad \left\{\begin{array}{l} {y}_1(t)=t_{11}x_1(t)+t_{12}x_2(t) \\ {y}_2(t)=t_{21}x_1(t)+t_{22}x_2(t) \end{array}\right. }

すなわち

\displaystyle{(14)\quad \underbrace{ \left[\begin{array}{l} {y}_1(t) \\ {y}_2(t) \end{array}\right] }_{y(t)}= \underbrace{ \left[\begin{array}{cc} t_{11} & t_{12} \\ t_{21} & t_{22} \\ \end{array}\right] }_{T} \underbrace{ \left[\begin{array}{l} x_1(t) \\ x_2(t) \end{array}\right] }_{x(t)} }

を考えます。またこの逆変換を次式で表しておきます。

\displaystyle{(15)\quad \left\{\begin{array}{l} {x}_1(t)=v_{11}y_1(t)+v_{12}y_2(t) \\ {x}_2(t)=v_{21}y_1(t)+v_{22}y_2(t) \end{array}\right. }

すなわち

\displaystyle{(16)\quad \underbrace{ \left[\begin{array}{l} {x}_1(t) \\ {x}_2(t) \end{array}\right] }_{x(t)}= \underbrace{ \left[\begin{array}{cc} v_{11} & v_{12} \\ v_{21} & v_{22} \\ \end{array}\right] }_{V} \underbrace{ \left[\begin{array}{l} y_1(t) \\ y_2(t) \end{array}\right] }_{y(t)} }

ここで、TV=VT=I_2が成り立つので、VTの逆行列T^{-1}V=T^{-1})となっています。逆行列をもつ行列は正則行列と呼ばれます。

すなわち、微分方程式\dot{x}(t)=Ax(t)に対して、正則行列Tによる変数変換y(t)=Tx(t)を行ないます。x(t)=Vy(t)V=T^{-1})と両辺を微分した\dot{x}(t)=V\dot{y}(t)を代入して次を得ます。

\displaystyle{(17)\quad V\dot{y}(t)=AVy(t)\Rightarrow\dot{y}(t)=\underbrace{V^{-1}AV}_{\Lambda}y(t) }

この変換された行列\Lambdaをできるだけ対角化します。その議論は、線形代数において\Lambda=V^{-1}AVを標準形とする方法として学びました。もし\exp(\Lambda t)を求めることができると

\displaystyle{(18)\quad x(t)=Vy(t)=V\exp(\Lambda t)y(0)=V\exp(\Lambda t)V^{-1}x(0) }

となって、漸近安定性(5)の判定は、

\displaystyle{(19)\quad \boxed{\exp(\Lambda t)\rightarrow 0\qquad(t\rightarrow\infty)} }

の判定に帰着されます。ここで、収束先の0n次の零行列です。

[4] さて、2次行列の標準形は次の3つに分類されます。

\displaystyle{(20)\quad \Lambda_1= \left[\begin{array}{cc} \lambda_1 & 0 \\ 0 & \lambda_2 \end{array}\right] \quad(\lambda(A)=\{\lambda_1,\lambda_2\}) }

\displaystyle{(21)\quad \Lambda_2= \left[\begin{array}{cc} \lambda_R & \lambda_I \\ -\lambda_I & \lambda_R \end{array}\right] \quad(\lambda(A)=\{\lambda_R\pm j\lambda_I\}) }

\displaystyle{(22)\quad \Lambda_3= \left[\begin{array}{cc} \lambda & 1 \\ 0 & \lambda \end{array}\right] \quad(\lambda(A)=\{\lambda,\lambda\}) }

ここで、jは虚数単位\sqrt{-1}\lambda_1,\lambda_2\lambda_R,\lambda_I\lambdaはすべて実数です。また任意のn次正方行列Xに対する行列式を{\rm det} Xで表すと

\displaystyle{(23)\quad \lambda(A)=\{\lambda:{\rm det}(\lambda I_n-A)=0\} }

すなわち、\lambda(A)は行列Aの固有値の集合です。{\rm det}(\lambda I_n-A)は特性多項式と呼ばれます。

上の3種類の\Lambdaについて\exp(\Lambda t)を具体的に計算してみます。その結果は次式となります。

\displaystyle{(24)\quad \exp(\Lambda_1 t)= \left[\begin{array}{cc} e^{\lambda_1t} & 0 \\ 0 & e^{\lambda_2 t} \end{array}\right] }

\displaystyle{(25)\quad \exp(\Lambda_2 t)=e^{\lambda_R t} \left[\begin{array}{cc} \cos(\lambda_It) & \sin(\lambda_It) \\ -\sin(\lambda_It)  & \cos(\lambda_It) \end{array}\right] }

\displaystyle{(26)\quad \exp(\Lambda_3 t)=e^{\lambda t} \left[\begin{array}{cc} 1 & t \\ 0 & 1 \end{array}\right] }

したがって、漸近安定性(5)すなわち(19)が成り立つのは、(24)のとき\lambda_1<0かつ\lambda_2<0、(25)のとき\lambda_R<0、(26)のとき\lambda<0であることがわかると思います。したがって、

漸近安定性のための条件は行列Aのすべての固有値の実部が負であること

です。対偶をとれば、一つでも実部が負の固有値があれば不安定と判定します。実部が負の固有値を安定固有値、実部が非負の固有値を不安定固有値と呼びます。

この漸近安定性の条件は、高次系についても同様に成り立ちます(Note A21-2を参照)。

上の剛体振り子の場合の漸近安定性を判定するために行列Aの固有値を調べてみます。まず、平衡状態(\alpha)の場合、次式を得ます。

(27)\quad \begin{array}{l} \displaystyle{{\rm det}(\lambda I_2-\left[\begin{array}{cc} 0 & 1 \\ -\frac{3g}{4\ell}&-\frac{c}{J} \end{array}\right]) =\lambda^2+\frac{c}{J}\lambda+\frac{3g}{4\ell}=0}\\ \displaystyle{\Rightarrow \lambda=\frac{1}{2}(-\frac{c}{J}\pm\sqrt{\frac{c^2}{J^2}-4\frac{3g}{4\ell}})} \end{array} }

これから行列Aの固有値の実部は必ず負となります(根号内が負の場合は明らか、正の場合は根号を開いても\frac{c}{J}より小)。したがって、平衡状態(\alpha)は漸近安定と判定できます。

次に、平衡状態(\beta)の場合、次式を得ます。

(28)\quad \begin{array}{l} \displaystyle{{\rm det}(\lambda I_2-\left[\begin{array}{cc} 0 & 1 \\ \frac{3g}{4\ell}&-\frac{c}{J} \end{array}\right]) =\lambda^2+\frac{c}{J}\lambda-\frac{3g}{4\ell}=0}\\ \displaystyle{\Rightarrow \lambda=\frac{1}{2}(-\frac{c}{J}\pm\sqrt{\frac{c^2}{J^2}+4\frac{3g}{4\ell}})} \end{array} }

これから行列Aの固有値は実数となり、一つは正、他は負となります。したがって、平衡状態(\beta)は漸近安定ではないと判定できます。

実は2次系の場合特性多項式の係数がすべて正であることが、漸近安定であるための条件であることがラウスフルビッツの判定法として知られています。

演習A21-1…Flipped Classroom

1^\circ 行列指数関数について次の指数法則が成り立つことを示せ。

\displaystyle{(29)\quad XY=YX\Rightarrow\exp(X+Y)=\exp(X)\exp(Y) }

2^\circ 行列指数関数の定義から、(24)を示せ。

3^\circ 行列指数関数の定義から、次式を用いて、(25)を示せ。

\displaystyle{(30)\quad \left[\begin{array}{cc} \lambda_R & \lambda_I \\ -\lambda_I & \lambda_R \end{array}\right]= \underbrace{\lambda_RI_2}_{X}+ \underbrace{\lambda_I \left[\begin{array}{cc} 0 & 1 \\ -1 & 0 \end{array}\right]}_{Y} }

4^\circ 行列指数関数の定義から、次式を用いて、(26)を示せ。

\displaystyle{(31)\quad \left[\begin{array}{cc} \lambda & 1 \\ 0 & \lambda \end{array}\right]= \underbrace{\lambda I_2}_{X}+ \underbrace{ \left[\begin{array}{cc} 0 & 1 \\ 0 & 0 \end{array}\right]}_{Y} }

演習A21-2…Flipped Classroom

与えられた行列Aに基づいて漸近安定性を調べるプログラムを次に示す。行列Pを実部が大きな順に表示されるように改良せよ。

MATLAB
%a21.m
%-----
 clear all, close all
 A=rand(10);
 p=eig(A)
%-----
 sys=ss(A,[],[],[]);
 damp(sys)
%----- 
 P=[real(p) imag(p) -cos(angle(p)) abs(p) -1./real(p)]
%-----
%eof
SCILAB
coming soon

Note A21-1 自由系の時間応答

n次自由系を表す微分方程式

\displaystyle{(1)\quad\dot{x}(t)=Ax(t)}

の解は

\displaystyle{(2)\quad x(t)=\exp(At)x(0)}

と表されます。これが解であることは元の微分方程式に代入すればすぐに確かめられ、また次のようにして導くことができます。

元の微分方程式を

\displaystyle{(3)\quad\dot{x}(t)-Ax(t)=0}

と書いて、左から積分因数と呼ばれる\exp(-At)をかけると

\displaystyle{(4)\quad\exp(-At)\dot{x}(t)-\exp(-At)Ax(t)=0}

すなわち

\displaystyle{(5)\quad\frac{d}{dt}(\exp(-At)x(t))=0}

したがって、cを定数ベクトルとして

\displaystyle{(6)\quad\exp(-At)x(t)=c\ \Rightarrow\ x(t)=\exp(At)c}

ここで、t=0とおくとcは初期値x(0)に等しいので

\displaystyle{(7)\quad x(t)=\exp(At)x(0)}

と表されます。

Note A21-2 高次系の漸近安定性

行列指数関数を用いると、微分方程式

\displaystyle{(1)\quad \dot{x}(t)=Ax(t)}

の解は次のように表されます。

\displaystyle{(2)\quad x(t)=\exp(At)x(0)=V\exp(\Lambda t)V^{-1}x(0)}

ここで

\displaystyle{(3)\quad  \begin{array}{ll} &\exp(At)=V\exp(\Lambda t)V^{-1}={\displaystyle \sum_{i=1}^p\sum_{k=1}^{m_i}t^{k-1}}e^{\lambda_it}\,R_{ik} \\ &+{\displaystyle \sum_{i=1}^q\sum_{k=1}^{\ell_i}t^{k-1}}e^{\lambda_{Ri}t}\cos\lambda_{Ii}t\,C_{ik} +{\displaystyle \sum_{i=1}^q\sum_{k=1}^{\ell_i}t^{k-1}}e^{\lambda_{Ri}t}\sin\lambda_{Ii}t\,S_{ik} \end{array}

と書けることに注意します(R_{ik},C_{ik},S_{ik}は適当なn次実正方行列)。

したがって、任意のx(0)\ne0に対して、t\rightarrow\inftyのときx(t)\rightarrow 0となるための条件は

\displaystyle{(4)\quad \exp(\Lambda t)\rightarrow 0\quad(t\rightarrow\infty)}

となります。これはAのすべての固有値の実部が負を意味します。

\dot{x}(t)=ax(t)a=1,0,0.5)の解のグラフを見ると、a=0の場合は、漸近安定ではないが、発散はしないので、不安定とまではいえないのではないかと思うかもしれません。したがって零の固有値を不安定とみなすのか、安定とみなすか迷うところです。しかし、\dot{x}(t)=Ax(t)において、A=\left[\begin{array}{cc} 0& 1\\ 0& 0 \end{array}\right]の場合、解はx(t)=\left[\begin{array}{cc} 1& t\\ 0& 1 \end{array}\right]x(0)となって、x(0)の第2要素が零でない場合は発散します。

いま、\dot{x}(t)=Ax(t)において、Aのすべての固有値の実部は負または零の場合を考え、実部が零の固有値\lambdaの重複度をq\ge2とし、次が成り立つとします。

\displaystyle{(5)\quad n-{\rm rank}(A-\lambda I_n)=q}

左辺は幾何学的重複度、左辺は代数学的重複度と呼ばれています。すべての実部が零の固有値について、幾何学的重複度と代数学的重複度が等しいことが、解が発散しないために追加的に要請される条件であることが知られています。上の例の場合、幾何学的重複度は2-1=1、代数学的重複度は2ですから、この条件は成立していないことがわかります。

Note A21-3 自由系の時間応答は固有値と固有ベクトルからどのように構成されるか

●2次自由系

\displaystyle{(1)\quad \dot{x}(t)=Ax(t)}

の解は次のように表されます。

\displaystyle{(2)\quad x(t)=\exp(At)x(0)}

2次行列Aの標準形は次の3つに分類されます。

\displaystyle{(3a)\quad \Lambda_1= \left[\begin{array}{cc} \lambda_1 & 0 \\ 0 & \lambda_2 \end{array}\right] \quad(\lambda(A)=\{\lambda_1,\lambda_2\}) }

\displaystyle{(3b)\quad \Lambda_2= \left[\begin{array}{cc} \lambda_R & \lambda_I \\ -\lambda_I & \lambda_R \end{array}\right] \quad(\lambda(A)=\{\lambda_R\pm j\lambda_I\}) }

\displaystyle{(3c)\quad \Lambda_3= \left[\begin{array}{cc} \lambda & 1 \\ 0 & \lambda \end{array}\right] \quad(\lambda(A)=\{\lambda,\lambda\}) }

ここで、jは虚数単位\sqrt{-1}\lambda_1,\lambda_2\lambda_R,\lambda_I\lambdaはすべて実数です。

このとき2次行列は、対応する固有ベクトルを用いて、次のように表されます。

\displaystyle{(4a)\quad A=\underbrace{\left[\begin{array}{cc} v_1 & v_2 \end{array}\right]}_{V} \underbrace{\left[\begin{array}{cc} \lambda_1 & 0 \\ 0 & \lambda_2 \end{array}\right]}_{\Lambda_1} \underbrace{\left[\begin{array}{cc} v_1 & v_2 \end{array}\right]^{-1}}_{V^{-1}} }

\displaystyle{(4b)\quad \begin{array}{l} A=\underbrace{\left[\begin{array}{cc} v_R+jv_I & v_R-jv_I  \end{array}\right]}_{V} \left[\begin{array}{cc} \lambda_R+j\lambda_I & 0 \\ 0 & \lambda_R-j\lambda_I \end{array}\right]\\ \times\underbrace{\left[\begin{array}{cc} v_R+jv_I & v_R-jv_I  \end{array}\right]^{-1}}_{V^{-1}}\\ =\underbrace{\left[\begin{array}{cc} v_R & v_I  \end{array}\right]}_{V'} \underbrace{\left[\begin{array}{cc} \lambda_R & \lambda_I \\ -\lambda_I & \lambda_R \end{array}\right]}_{\Lambda_2} \underbrace{\left[\begin{array}{cc} v_R & v_I   \end{array}\right]^{-1}}_{V'^{-1}} \end{array} }

\displaystyle{(4c)\quad A=\underbrace{\left[\begin{array}{cc} v & v'  \end{array}\right]}_{V} \underbrace{\left[\begin{array}{cc} \lambda & 1 \\ 0 & \lambda \end{array}\right]}_{\Lambda_3} \underbrace{\left[\begin{array}{cc} v & v' \end{array}\right]^{-1}}_{V^{-1}} }

ここで、v_1,v_2v_R,v_Ivv'はすべて2次元の実ベクトルです(v'は一般化固有ベクトル)。

これらの行列指数関数は次式となります。

\displaystyle{(5a)\quad \exp(A t)=\underbrace{\left[\begin{array}{cc} v_1 & v_2 \end{array}\right]}_{V} \underbrace{\left[\begin{array}{cc} e^{\lambda_1t} & 0 \\ 0 & e^{\lambda_2 t} \end{array}\right]}_{\exp(\Lambda_1 t)} \underbrace{\left[\begin{array}{cc} v_1 & v_2 \end{array}\right]^{-1}}_{V^{-1}} }

\displaystyle{(5b)\quad \exp(A t)=\underbrace{\left[\begin{array}{cc} v_R & v_I  \end{array}\right]}_{V'} \underbrace{e^{\lambda_R t} \left[\begin{array}{cc} \cos(\lambda_It) & \sin(\lambda_It) \\ -\sin(\lambda_It)  & \cos(\lambda_It) \end{array}\right]}_{\exp(\Lambda_2 t)} \underbrace{\left[\begin{array}{cc} v_R & v_I   \end{array}\right]^{-1}}_{V'^{-1}} }

\displaystyle{(5c)\quad \exp(A t)=\underbrace{\left[\begin{array}{cc} v & v' \end{array}\right]}_{V} \underbrace{e^{\lambda t} \left[\begin{array}{cc} 1 & t \\ 0 & 1 \end{array}\right]}_{\exp(\Lambda_3 t)} \underbrace{\left[\begin{array}{cc} v & v' \end{array}\right]^{-1}}_{V^{-1}} }

これらを(2)に代入して

\displaystyle{(6a)\quad \begin{array}{l} x(t)=\underbrace{\left[\begin{array}{cc} v_1 & v_2 \end{array}\right]}_{V} \underbrace{\left[\begin{array}{cc} e^{\lambda_1t} & 0 \\ 0 & e^{\lambda_2 t} \end{array}\right]}_{\exp(\Lambda_1 t)} \underbrace{\left[\begin{array}{cc} v_1 & v_2 \end{array}\right]^{-1}x(0)}_{\left[\begin{array}{c} c_1\\c_2 \end{array}\right]}\\ =c_1e^{\lambda_1t}v_1+c_2e^{\lambda_2t}v_2 \end{array} }

\displaystyle{(6b)\quad \begin{array}{l} x(t)=\underbrace{\left[\begin{array}{cc} v_R & v_I  \end{array}\right]}_{V'} \underbrace{e^{\lambda_R t} \left[\begin{array}{cc} \cos(\lambda_It) & \sin(\lambda_It) \\ -\sin(\lambda_It)  & \cos(\lambda_It) \end{array}\right]}_{\exp(\Lambda_2 t)} \underbrace{\left[\begin{array}{cc} v_R & v_I   \end{array}\right]^{-1}x(0)}_{\left[\begin{array}{c} c_1\\c_2 \end{array}\right]}\\ =c_1e^{\lambda_Rt}(\cos(\lambda_It)v_R-\sin(\lambda_It)v_I)+c_2e^{\lambda_Rt}(\sin(\lambda_It)v_R+\cos(\lambda_It)v_I) \end{array} }

\displaystyle{(6c)\quad \begin{array}{l} x(t)=\underbrace{\left[\begin{array}{cc} v & v' \end{array}\right]}_{V} \underbrace{e^{\lambda t} \left[\begin{array}{cc} 1 & t \\ 0 & 1 \end{array}\right]}_{\exp(\Lambda_3 t)} \underbrace{\left[\begin{array}{cc} v & v' \end{array}\right]^{-1}x(0)}_{\left[\begin{array}{c} c_1\\c_2 \end{array}\right]}\\ =c_1e^{\lambda t}v+c_2e^{\lambda t}(tv+v') \end{array} }

●一般の自由系についても、同様の式が成り立ち、Aの固有値と固有ベクトルがどのように時間応答に寄与しているかの分析ができ、モード間の連成について考察することができます。