B02 リャプノフの安定定理

リャプノフの安定定理

[1]  非線形システム理論に関する文献:
Hassan K.Khalil: Nonlinear Systems, 3rd ed., Prentice Hall, 2002
から次の定理(Theorem 4.1, p.114)を引用します。

リャプノフの安定定理 状態方程式

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

に対して、その平衡状態 x^*=0 を含む領域 {\cal D} 考える。いま {\cal D} で連続微分可能な実数値関数 V に対して

\displaystyle{(2)\quad \left\{\begin{array}{l} V(0)=0 \\ V(x)>0\quad(x\in{\cal D},x\not=0) \end{array}\right.}

このとき

\displaystyle{(3)\quad \frac{d}{dt}V(x)\le0 \quad(x\in{\cal D})}

ならば、平衡状態 x^*=0 は安定、また

\displaystyle{(4)\quad \frac{d}{dt}V(x)<0 \quad(x\in{\cal D},x\not=0)}

ならば、平衡状態 x^*=0 は漸近安定である。

すなわち、(2)と(4)を満足するリャプノフ関数Vの存在を示して(発見して)、非線形系(1)の漸近安定性を示すことができます。

●振り子の例について、エネルギーの変化と安定性の関係を考えてみます。その運動方程式は,摩擦力を入れた場合

\displaystyle{(5)\quad  \underbrace{\frac{4}{3}m\ell^2}_J \dot{\omega}(t)=-mg\ell\sin\theta(t)-c\,\omega(t)\quad(\omega(t)=\dot{\theta}(t)) }

です。いま、運動エネルギーと位置エネルギーの和を

\displaystyle{(6)\quad  E=\frac{1}{2}J\omega^2(t)+mg\ell(1-\cos\theta(t))>0 }

と書くと、(5)を用いて

\displaystyle{(7)\quad  \dot{E}= (J\dot{\omega}(t)+mg\ell\sin\theta(t))\omega(t) =-c\,\omega^2(t) \le 0 }

が成り立つちます。このようにエネルギーが減少することは、振り子の揺れが止まることと対応するので、エネルギーの時間変化率から安定性を判定できそうです。

リャプノフの安定定理は十分条件であり、リャプノフ関数が見つからないからといって安定性や漸近安定性が成り立たないとはかぎりません。また、リャプノフ関数の構築法について一般論はないのですが、上でみたようにエネルギーとの関連が深いです。そこで、振り子の平衡状態

\displaystyle{(8)\quad  x^*= \left[\begin{array}{c} \theta^* \\ \omega^* \end{array}\right]= \left[\begin{array}{c} 0 \\ 0  \end{array}\right] }

の安定性を調べるために、(6)のエネルギーEをリャプノフ関数の候補としてみます。まず、(6)と(7)から、Ex^*=0の安定性を保証するためのリャプノフ関数といえます。

つぎに、(7)から

\displaystyle{(9)\quad  \omega(t)=0\ \Rightarrow\ \dot{E}=-c\,\omega^2(t)=0\quad (\forall \theta(t)) }

ですが、これはEx^*=0は漸近安定性を保証するためのリャプノフ関数とはならないことを示しています。しかし、実際は漸近安定なので、ほかにリャプノフ関数の候補を探せる可能性があります。

例えば、エネルギーEの代わりに、つぎのリャプノフ関数の候補を考えます。

\displaystyle{(10)\quad  V(t)=\frac{1}{2} \left[\begin{array}{c} \theta(t) \\ \omega(t) \end{array}\right]^T \left[\begin{array}{cc} \frac{1}{2}\frac{c^2}{J} & \frac{1}{2}c \\ \frac{1}{2}c & J \end{array}\right] \left[\begin{array}{c} \theta(t) \\ \omega(t) \end{array}\right] +mg\ell(1-\cos\theta(t))>0 }

この時間微分を計算すると、(5)を用いて

\displaystyle{(11)\quad  \dot{V}(t)=-\frac{1}{2}\frac{3g}{4\ell}c\,\theta(t)\sin\theta(t) -\frac{1}{2}c\,\omega^2(t) }

となります。ここで

\displaystyle{(12)\quad  \theta(t)\sin\theta(t)>0\quad (|\theta(t)|<\pi,\theta(t)\ne0) }

だから、(11)は負値をとります。したがって、(10)はリャプノフ関数となっています。

[2]  リャプノフの安定定理を線形系

\displaystyle{(17)\quad \dot{x}(t)=Ax(t)\quad(x(t)\in{\bf R}^n)}

に適用して、漸近安定性を保証するリャプノフ関数を求めてみます。その候補として次を考えます。

\displaystyle{(18)\quad V(x)=x^T(t)Px(t)\quad(P>0)}

これを微分して(17)を用いると

\displaystyle{(19)\quad \begin{array}{l} \frac{d}{dt}V(x)=\dot{x}^T(t)Px(t)+x^T(t)P\dot{x}(t)\\ =x^T(t)A^TPx(t)+x^T(t)PAx(t)\\ =x^T(t)(A^TP+PA)x(t) \end{array} }

となります。したがって、線形行列不等式

\displaystyle{(20)\quad \boxed{{A^TP+PA<0\quad(P>0)}}}

が成り立てば、(2)と(4)が満足され、線形系(17)の漸近安定性が成り立つといえます。

●逆に、もし線形系(17)が漸近安定、すなわち

\displaystyle{(21)\quad x(t)=\exp(At)x(0)\rightarrow 0\quad(t\rightarrow \infty)}

であれば、線形行列不等式(20)を満たすP>0を次のように決めることができます。

\displaystyle{(22)\quad P=\int_0^\infty \exp(A^Tt)Q\exp(At)dt}

ここで、Q>0は任意です(ある条件下ではQ\ge0としてもよい)。実際

\displaystyle{(23)\quad \begin{array}{l} \displaystyle{A^T\int_0^\infty \exp(A^Tt)Q\exp(At)dt+\int_0^\infty \exp(A^Tt)Q\exp(At)dt\,A}\\ \displaystyle{=\int_0^\infty \frac{d}{dt}(\exp(A^Tt)Q\exp(At))\,dt}\\ \displaystyle{=\left[\exp(A^Tt)Q\exp(At)\right]_0^\infty=-Q<0} \end{array} }

となって(20)を満足します。また、任意のv\ne0に対して

(24)\quad \begin{array}{l} \displaystyle{v^TPv=\int_0^\infty v^T\exp(A^Tt)Q\exp(At)vdt}\\ \displaystyle{=\int_0^\infty ||Q^{1/2}\exp(At)v||dt\ge0} \end{array}

ですが、v^TPv=0とすると、Q^{1/2}\exp(At)v=0が得られ、これはv=0を意味し矛盾です。したがってv^TPv>0でなけれならず、P>0を得ます。以上から、線形行列不等式(20)が(17)の漸近安定性の必要十分条件になっています。

●(20)において、X=PY=P^{-1}とおくと、漸近安定性の等価な条件として次が得られます。

\displaystyle{(25)\quad \exists X>0: A^TX+XA<0}
\displaystyle{(26)\quad \exists Y>0: YA^T+PY<0}

本資料では同様な多くのLMIが登場しますが、(25)のタイプをX-LMI、(26)のタイプをY-LMIと呼びます。X-LMIは特性解析のために、Y-LMIはコントローラ設計(シンセシス)のために使い分けられます。

演習B02…Flipped Classroom
1^\circ 安定行列Aに対して、線形行列不等式(20)を満足する正定行列Pは、次のプログラムで計算できることを確かめよ。

MATLAB
%ana_lmi0.m
%-----
 A=[0 1;2 -3]; n=2; 
%-----
 setlmis([]); 
 X=lmivar(1,[n 1]); 
%-----
 lmi1=newlmi; 
 lmiterm([lmi1 1 1 X],1,A,'s'); %#1:X*A+A'*X 
%-----
 lmi2=newlmi; 
 lmiterm([-lmi2 1 1 X],1,1);    %#2:X 
%-----
 LMIs=getlmis;  
 [tmin,xfeas]=feasp(LMIs); 
 X=dec2mat(LMIs,xfeas,X) 
%-----
%eof