ロバスト性 | 制御系CAD

線形系の２次安定性…Homework

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

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

これに出力方程式

$\displaystyle{(2)\quad y(t)=Cx(t) }$

を考慮する場合も含めると、漸近安定性の判定法は次のようにまとめられます。

【漸近安定性の定義とその等価な条件】

定義DA:　 $\forall x(0)\ne 0: x(t)=\exp(At)x(0)\rightarrow 0\quad(t\rightarrow\infty)$

条件A0:　 $\exp(At)\rightarrow 0\quad(t\rightarrow\infty)$

条件A1:　 ${\rm Re}(\lambda_i(A))<0\quad(i=1,\cdots,n)$

条件A2:　 $\exists P>0: PA+A^TP+I_n=0$

条件A3:　 $\exists P>0: PA+A^TP+C^TC=0$ 　ただし、 $(A,C)$ は可観測対

条件A4:　 $\exists P>0: PA+A^TP<0$

条件A2のリャプノフ方程式から次式を得ます。

$\displaystyle{(3)\quad x^T(t)(PA+A^TP)x(t)=-x(t)^Tx(t) }$

線形系(1)に対してリャプノフ関数

$\displaystyle{(4)\quad V(x)=x^T(t)Px(t) }$

を考えるとき、(3)は２次安定性の条件

$\displaystyle{(5)\quad \frac{d}{dt}V(x)=2x^T(t)P\dot{x}(t)=2x^T(t)PAx(t)=-x^T(t)Qx(t)\quad (Q=I) }$

を意味します。これから漸近安定な線形系(1)は２次安定であることが分かります。

●条件A2の両辺から $Q^{1/2}\ (Q>0)$ をかけて

$\displaystyle{(6)\quad \underbrace{Q^{1/2}PQ^{1/2}}_{P'}\underbrace{Q^{-1/2}AQ^{1/2}}_{A'}+\underbrace{Q^{1/2}A^TQ^{-1/2}}_{A'^T}\underbrace{Q^{1/2}PQ^{1/2}}_{P'}+\underbrace{Q^{1/2}Q^{1/2}}_{Q}=0 }$

を得ます。 $(A',Q^{1/2})$ は可観測対なので条件A3を適用して（ $'$ をはずして）

条件A2′:　 $\boxed{\exists P>0: PA+A^TP+Q=0}$ 　ただし、 $Q>0$

も漸近安定性の条件となることに注意します。

モデルの不確かさがある場合…Homework

いま漸近安定な線形系(1)において次のようなモデルの不確かさがあるとします。

$\displaystyle{(7)\quad \dot{x}(t)=(A+\Delta A)x(t)=Ax(t)+\underbrace{\Delta Ax(t)}_{\xi(t,x)}\qquad(x(0)\ne0) }$

条件A2’より、適当な $P>0$ と $Q>0$ に対して

$\displaystyle{(8)\quad PA+A^TP+Q=0 }$

が成り立ちます。このとき

$\displaystyle{(9)\quad \boxed{||\xi(t,x)||\le\frac{1}{2}\frac{\sigma_n(Q)}{\sigma_1(P)}||x(t)||} }$

ならば、リャプノフ関数

$\displaystyle{(10)\quad V(x)=x^T(t)Px(t) }$

に対して

$\displaystyle{(11)\quad \frac{d}{dt}V(x)<0 }$

となって、リャプノフ安定となります。

●実際、

$\displaystyle{(12)\quad \begin{array}{l} \frac{d}{dt}V(x)=x^T(t)P\dot{x}(t)+\dot{x}^T(t)Px(t)\\ =x^T(t)P(Ax(t)+\xi(t,x))+(x^T(t)A^T+\xi^T(t,x))Px(t)\\ =x^T(t)PAx(t)+x^T(t)A^TPx(t)+x^T(t)P\xi(t,x)+\xi^T(t,x)Px(t)\\ =x^T(t)(PA+A^TP)x(t)+2\xi^T(t,x)Px(t)\\ =-x^T(t)Qx(t)+2\xi^T(t,x)Px(t)\\ \le -x^T(t)Qx(t)+ 2||\xi(t,x)||\underbrace{||Px(t)||}_{\sqrt{x^T(t)P^2x(t)}} \end{array} }$
ここで、

$\displaystyle{(13)\quad \begin{array}{l} \sigma_n(Q)x^T(t)x(t)\le x^T(t)Qx(t) \Rightarrow \sigma_n(Q)||x(t)||^2\le x^T(t)Qx(t)\\ x^T(t)P^2x(t)\le\sigma_1(P^2)x(t)^Tx(t) \Rightarrow \sqrt{x^T(t)P^2x(t)}\le\sigma_1(P)||x(t)|| \end{array} }$

に注意して、(9)の下では

$\displaystyle{(14)\quad \begin{array}{l} \frac{d}{dt}V(x) \le -\sigma_n(Q)||x(t)||^2+ 2||\xi(t,x)||\sigma_1(P)||x(t)|| \\ =(||\xi(t,x)||-\frac{1}{2}\frac{\sigma_n(Q)}{\sigma_1(P)}||x(t)||)2\sigma_1(P)||x(t)||<0 \end{array} }$

となります。

●次の数値例を考えてみます。

$\displaystyle{(15)\quad \left[\begin{array}{c} \dot{\theta}(t) \\ \dot{\omega}(t) \end{array}\right] = \underbrace{ \left[\begin{array}{cc} 0 & 1 \\ -4 & -1 \end{array}\right] }_{A} \left[\begin{array}{c} \theta(t) \\ \omega(t) \end{array}\right] + \underbrace{\left[\begin{array}{c} 0 \\ 1 \end{array}\right]}_{D}\underbrace{\sin(2t)}_{\xi(t)} }$

$A$ は安定行列なので、リャプノフ方程式 $PA+A^TP+I_n=0$ を満足する $P>0$ は

$\displaystyle{(16)\quad P=\left[\begin{array}{cc} 2.6250 & 0.1250\\ 0.1250 & 0.6250 \end{array}\right] }$

のように求まります。

この $P$ をもつリャプノフ関数(10)に対して

$\displaystyle{(17)\quad \begin{array}{l} \frac{d}{dt}V(x)=x^T(t)P\dot{x}(t)+\dot{x}^T(t)Px(t)\\ =x^T(t)P(Ax(t)+D\xi(t))+(x^T(t)A^T+\xi^T(t)D^T)Px(t)\\ =x^T(t)PAx(t)+x^T(t)A^TPx(t)+x^T(t)PD\xi(t)+\xi(t)D^TPx(t)\\ =x^T(t)(PA+A^TP)x(t)+2\xi(t)D^TPx(t)\\ =-x^T(t)Qx(t)+2\xi(t)D^TPx(t)\\ \le -x^T(t)x(t)+ 2|\xi(t)|||D^TP||||x(t)||\\ \le -||x(t)||^2+ 2\sigma_1(PD)||x(t)|| \\ =(2\sigma_1(PD)-||x(t)||)||x(t)||<0 \end{array} }$

これから

$\displaystyle{(18)\quad ||x(t)||>2\sigma_1(PD)=1.2748 }$

のときは $\frac{d}{dt}V(x)<0$ となって原点に漸近していきますが、 $||x(t)||<1.2748$ になるとそれは保証されないことが分かります(図1参照)。

図1 モデルの不確かさがある場合の解軌道

演習C12…Flipped Classroom
$1^\circ$ 　smc17.slxを用いて、図1のシミュレーションを行え。

Note C12 振り子の安定性

振り子の運動方程式

$\displaystyle{(1)\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{(2)\quad \frac{d}{dt} \underbrace{\left[\begin{array}{c} \theta(t) \\ \omega(t) \end{array}\right]}_{x} = \underbrace{ \left[\begin{array}{c} \omega(t) \\ -\frac{mg\ell}{J}\sin\theta(t)-\frac{c}{J}\omega(t) \end{array}\right] }_{f(x)} }$

に対して、次のリャプノフ関数の候補を考えます。

$\displaystyle{(3)\quad V=\frac{1}{2} \underbrace{ \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] }_{x^T(t)Px(t)} +mg\ell(1-\cos\theta(t))>0 }$

この時間微分を計算すると

$\displaystyle{(4)\quad \begin{array}{l} \dot{V}(t)=\frac{1}{2}(x^TP\dot{x}+\dot{x}^TPx)+mg\ell\omega\sin\theta =x^TP\dot{x}+mg\ell\omega\sin\theta\\ =\left[\begin{array}{c} \theta \\ \omega \end{array}\right]^T \left[\begin{array}{cc} p_1 & p_3 \\ p_3 & p_2 \end{array}\right] \left[\begin{array}{c} \omega \\ -\frac{mg\ell}{J}\sin\theta-\frac{c}{J}\omega \end{array}\right]+mg\ell\omega\sin\theta\\ =\left[\begin{array}{cc} \theta &\omega \end{array}\right] \left[\begin{array}{cc} p_1\omega -p_3(\frac{mg\ell}{J}\sin\theta+\frac{c}{J}\omega) \\ p_3\omega -p_2(\frac{mg\ell}{J}\sin\theta+\frac{c}{J}\omega) \end{array}\right]+mg\ell\omega\sin\theta\\ =p_1\theta\omega -p_3(\frac{mg\ell}{J}\theta\sin\theta+\frac{c}{J}\theta\omega)+ p_3\omega^2 -p_2(\frac{mg\ell}{J}\omega\sin\theta+\frac{c}{J}\omega^2)+mg\ell\omega\sin\theta\\ =-\frac{mg\ell}{J}\sin\theta(p_3\theta+p_2\omega-J\omega) +(p_3-p_2\frac{c}{J})\omega^2+(p_1-p_3\frac{c}{J})\theta\omega \end{array} }$

ここで、 $p_1=\frac{1}{2}\frac{c^2}{J}$ 、 $p_2=J$ 、 $p_3=\frac{1}{2}c$ を代入して

$\displaystyle{(5)\quad \dot{V}(t)=-\frac{1}{2}\frac{mg\ell}{J}c\,\theta(t)\sin\theta(t) -\frac{1}{2}c\,\omega^2(t) }$

となります。ここで

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

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

●ちなみに(2)を次のように書き直してみます。

$\displaystyle{(7)\quad \left[\begin{array}{c} \dot{\theta}(t) \\ \dot{\omega}(t) \end{array}\right] = \underbrace{ \left[\begin{array}{cc} 0 & 1 \\ -\frac{mg\ell}{J} & -\frac{c}{J} \end{array}\right] }_{A} \left[\begin{array}{c} \theta(t) \\ \omega(t) \end{array}\right] + \underbrace{ \left[\begin{array}{c} 0 \\ \frac{mg\ell}{J} \end{array}\right] }_{D} \underbrace{ (\theta(t)-\sin\theta(t)) }_{\xi(x)} }$

我々は軸摩擦をもつ振り子の角度は十分時間がたてば零に収束することを知っています。上の議論に沿って解釈すれば、モデルの不確かさ $\xi(x)$ 自体が限りなく零に近づき（ $\sin\theta\simeq\theta$ ）、線形系となると考えられます。

Note C12 　MATLABによるリャプノフ方程式の求解

たとえば、安定行列 $\displaystyle{ A= \left[\begin{array}{cc} 0 & 1 \\ -4 & -1 \end{array}\right] }$

が与えられるとき、リャプノフ方程式 $PA+A^TP+I_n=0$ を満足する $P>0$ は

$\displaystyle{ P=\left[\begin{array}{cc} 2.6250 & 0.1250\\ 0.1250 & 0.6250 \end{array}\right] }$

のように求まります。これはMATLABで次のコマンドを用いました。

P=lyap(A’,eye(2))

ここで、P=lyap(A,Q) はリャプノフ方程式 $PA^T+AP+Q=0$ の解を求めることに注意します。この解を $(A,Q)$ に対するリャプノフ行列と呼ぶことがあります。