制御系CAD

制御を意識したモノつくりを目指して

性能解析LMI（Ｈ∞ノルム）

投稿日時: 2021年10月29日投稿者: cacsd

性能解析LMI（Ｈ∞ノルム）…Homework

[1]　直達項をもつ $n$ 次系の $H_\infty$ ノルムが $\gamma$ より小である条件は ${\rm\bf LMI}$ を用いて、次のように表されます。

図1 直達項をもつｎ次系

$\displaystyle{(1)\quad \begin{array}{l} \displaystyle{\sup_{u\in{\cal L}_2}\frac{||y(t)||_2}{||u(t)||_2} =\sup_{\omega\in{\rm\bf R}}||\hat{G}(j\omega)||_2<\gamma} \\ \displaystyle{\Leftrightarrow \exists X>0:\ \left[\begin{array}{ccc} A^TX+XA & XB & C^T \\ B^TX & -\gamma^2 I & D^T \\ C & D & -I \end{array}\right]<0} \\ \displaystyle{\Leftrightarrow \exists Y>0:\ \left[\begin{array}{ccc} YA^T+AY & B & YC^T \\ B^T & -\gamma^2 I & D^T \\ CY & D & -I \end{array}\right]<0} \end{array}}$

実際、シュール補元に関するLMIより

$\displaystyle{(2)\quad \begin{array}{l} \left[\begin{array}{ccc} A^TX+XA & XB & C^T \\ B^TX & -\gamma^2 I & D^T \\ C & D &-I \end{array}\right]<0 \\ \Leftrightarrow \left[\begin{array}{cc} A^TX+XA & XB \\ B^TX & -\gamma^2 I \end{array}\right] - \left[\begin{array}{cc} C^T\\ D^T \end{array}\right] (-I) %(-\gamma^{-1} I) \left[\begin{array}{cc} C & D \end{array}\right] <0 \\ \Leftrightarrow \left[\begin{array}{cc} A^TX+XA & XB \\ B^TX & 0 \end{array}\right] < \left[\begin{array}{cc} C^T & 0 \\ D^T & I \end{array}\right] \left[\begin{array}{cc} -I & 0 \\ 0 & \gamma^2 I \end{array}\right] \left[\begin{array}{cc} C & D \\ 0 & I \end{array}\right] \end{array} }$

が成り立つことに注意して

$\displaystyle{(3)\quad \begin{array}{l} \left[\begin{array}{ccc} A^TX+XA & XB & C^T \\ B^TX & -\gamma^2 I & D^T \\ C & D & - I \end{array}\right]<0 \\ &\Leftrightarrow \forall \left[\begin{array}{c} x \\ u \end{array}\right]\ne0: \underbrace{ \left[\begin{array}{c} x \\ u \end{array}\right]^T \left[\begin{array}{cc} A^TX+XA & XB \\ B^TX & 0 \end{array}\right] \left[\begin{array}{c} x \\ u \end{array}\right] }_{\dot{V}(x)=\frac{d}{dt}(x^TXx)} \\ < \underbrace{ \left[\begin{array}{c} x \\ u \end{array}\right]^T \left[\begin{array}{cc} C & D \\ 0 & I \end{array}\right]^T \left[\begin{array}{cc} - I & 0 \\ 0 & \gamma^2 I \end{array}\right] \left[\begin{array}{cc} C & D \\ 0 & I \end{array}\right] \left[\begin{array}{c} x \\ u \end{array}\right] }_{s(u,y)=\gamma^2 u^Tu-y^Ty} \end{array} }$

すなわち、条件

$\displaystyle{(4)\quad \underbrace{\frac{d}{dt}(x^TXx)}_{\dot{V}(x)} <\underbrace{\gamma^2 u^Tu-y^Ty}_{s(u,y)} }$

が成り立ち、これを積分して

$\displaystyle{(5)\quad \underbrace{x^T(\infty)Xx(\infty)}_{V(x(\infty))} -\underbrace{x^T(0)Xx(0)}_{V(x(0))=0} <\underbrace{\int_0^\infty (\gamma^2u^Tu-y^Ty)\,dt}_{\int_0^\infty s(u,y)dt} }$

を得ます。漸近安定性より左辺は正でなければならないので

$(6)\quad \begin{array}{l} \displaystyle{0<\gamma^2\underbrace{\int_0^\infty u^Tu\,dt}_{||u(t)||_2^2} -\underbrace{\int_0^\infty y^Ty\,dt}_{||y(t)||_2^2}}\\ \displaystyle{\Leftrightarrow ||y(t)||_2^2<\gamma^2||u(t)||_2^2}\\ \displaystyle{\Leftrightarrow \sup_{u\in{\cal L}_2}\frac{||y(t)||_2}{||u(t)||_2}<\gamma} \end{array} }$

が成り立ちます。

●入力 $\sqrt{\gamma}u(t)$ から出力 $\frac{1}{\sqrt{\gamma}}y(t)$ への入出力特性は不変であることから

$\displaystyle{(7) \begin{array}{l} \left[\begin{array}{ccc} A^TX+XA & XB & C^T \\ B^TX & -\gamma^2 I & D^T \\ C & D & -I \end{array}\right]<0 \\ \Leftrightarrow \left[\begin{array}{ccc} A^TX+XA & X(B\sqrt{\gamma}) & (\frac{1}{\sqrt{\gamma}}C)^T \\ (B\sqrt{\gamma})^TX & - \gamma^2 I & D^T \\ (\frac{1}{\sqrt{\gamma}}C) & D & -I \end{array}\right]<0 \\ \Leftrightarrow \left[\begin{array}{ccc} I & 0 & 0 \\ 0 & \sqrt{\gamma} I & 0 \\ 0 & 0 & \frac{1}{\sqrt{\gamma}} I \end{array}\right] \left[\begin{array}{ccc} A^TX+XA & XB & C^T \\ B^TX & -\gamma I & D^T \\ C & D &-\gamma I \end{array}\right] \left[\begin{array}{ccc} I & 0 & 0 \\ 0 & \sqrt{\gamma} I & 0 \\ 0 & 0 & \frac{1}{\sqrt{\gamma}} I \end{array}\right]<0 \\ \Leftrightarrow \left[\begin{array}{ccc} A^TX+XA & XB & C^T \\ B^TX & -\gamma I & D^T \\ C & D &-\gamma I \end{array}\right]<0 \end{array} }$

を得て、次が成り立ちます。

$\displaystyle{(8)\quad \begin{array}{l} \displaystyle{\sup_{u\in{\cal L}_2}\frac{||y(t)||_2}{||u(t)||_2} =\sup_{\omega\in{\rm\bf R}}||\hat{G}(j\omega)||_2<\gamma} \\ \displaystyle{\Leftrightarrow \exists X>0:\ \left[\begin{array}{ccc} A^TX+XA & XB & C^T \\ B^TX & -\gamma I & D^T \\ C & D & -\gamma I \end{array}\right]<0} \\ \displaystyle{\Leftrightarrow \exists Y>0:\ \left[\begin{array}{ccc} YA^T+AY & B & YC^T \\ B^T & -\gamma I & D^T \\ CY & D & -\gamma I \end{array}\right]<0} \end{array}}$

演習B22…Flipped Classroom

$1^\circ$ ２次振動系に対して、適当な $\gamma$ を設定し、 $H_\infty$ ノルムと$の大小関係をチェックするプログラムを作成せよ。

MATLAB

%ana_lmi5.m
%------
 clear all, close all
 zeta=0.4; Mp=0.5/zeta/sqrt(1-zeta^2)
 A=[0 1;-1 -2*zeta]; B=[0;1]; 
 C=[1 0]; D=0; n=2; 
 gam=input('gamma =  '); 
%------
 setlmis([]);
 X=lmivar(1,[n 1]);
%------
 lmi1=newlmi;   
 lmiterm([lmi1 1 1 X],1,A,'s'); %#1:X*A+A'*X 
 lmiterm([lmi1 1 2 X],1,B);     %#1:X*B 
 lmiterm([lmi1 2 2 0],-gam);    %#1:-gam 
 lmiterm([lmi1 3 1 0],C);       %#1:C 
 lmiterm([lmi1 3 2 0],D);       %#1:D 
 lmiterm([lmi1 3 3 0],-gam);    %#1:-gam
%------ 
 lmi2=newlmi;   
 lmiterm([-lmi2 1 1 X],1,1);    %#2:X 
%------ 
 LMIs=getlmis; 
 [tmin,xfeas]=feasp(LMIs);
 X=dec2mat(LMIs,xfeas,X) 
%------
%eof

図2 ２次系の周波数応答の比較

伝達特性

投稿日時: 2021年10月29日投稿者: cacsd

伝達特性…Homework

[0]　 $n$ 次系の応答の相互関係は次のようにまとめることができました。

図1 線形系の時間応答

ここで、 $n$ 次系のインパルス応答が分かれば、どのような入力を与えられても出力を計算できます。または伝達関数が分かれば同様のことが言えます。したがって $n$ 次系は漸近安定であるとして、その性能（伝達特性）の評価はインパルス応答または伝達関数を用いて行われます。伝達関数は状態空間表現のうち可制御かつ可観測な部分しか反映しないので、以下では $(A,B)$ は可制御対、 $(A,C)$ は可観測対と仮定しておきます。

[1]　 $x\in{\rm\bf R}^n$ から $y\in{\rm\bf R}^m$ への写像 $y=Ax$ の「伝達特性」をどう測るかを考えます。これはスカラの場合は正比例の関係 $y=ax$ ですから、比例定数 $a$ に相当する量を求める話になります。

図2　線形写像

サイズ $m\times n$ の行列 $A$ の次の特異値分解を考えます。

$\displaystyle{(1)\quad \boxed{\begin{array}{l} A= \underbrace{ \left[\begin{array}{cc} U_1 & U_2 \end{array}\right] }_{U(m\times m)} \underbrace{ \left[\begin{array}{cc} {\rm diag}\{\sigma_1,\cdots,\sigma_k\} & 0_{k\times (n-k)} \\ 0_{(m-k)\times k} & 0_{(m-k)\times (n-k)} \end{array}\right] }_{\Sigma= \left[\begin{array}{cc} \Sigma_1 & 0 \\ 0 & 0 \end{array}\right] \quad(\sigma_1\ge\cdots\ge\sigma_k) } \underbrace{ \left[\begin{array}{cc} V_1^T \\ V_2^T \end{array}\right] }_{V(n\times n)^T}\\ =U_1(m\times k)\Sigma_1(k\times k)V_1(n\times k)^T \end{array}} }$

ここで、 $k={\rm rank}A$ で、 $U$ と $V$ は直交行列です。

$\displaystyle{(2)\quad UU^T=U^TU=I_m,\quad VV^T=V^TV=I_n }$

したがって、次のような３つの線形写像に分解されます。

図3　特異値分解

$n$ 次元ベクトル $x\in{\rm\bf R}^n$ のノルムとして、次の３通りがよく用いられ、それぞれ１ノルム、２ノルム、 $\infty$ ノルムと呼ばれます。

$\displaystyle{(3)\quad \boxed{\left\{\begin{array}{lll} ||x||_1=|x_1|+\cdots+|x_n|&\\ ||x||_2=\sqrt{x_1^2+\cdots+x_n^2}&\Rightarrow&||x||_2^2=x^Tx\\ ||x||_\infty=\max\{|x_1|,\cdots,|x_n|\}& \end{array}\right.} }$

以下では、ベクトルのノルムとして、２番目の２ノルムを考えます。

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

$\displaystyle{(4)\quad ||x'||_2^2=||Vx||_2^2=x^TV^TVx=x^Tx=||x||_2^2\ \Rightarrow\ \frac{||x'||_2}{||x||_2}=1 }$

$\displaystyle{(5)\quad \begin{array}{l} ||y'||_2^2=||\Sigma x'||_2^2=\sigma_1^2x_1'\,^2+\cdots+\sigma_k^2x_k'\,^2+0x_{n+1}'\,^2+\cdots+0x_n'\,^2\\ \le \sigma_1^2(x_1'\,^2+\cdots+x_n'\,^2)=\sigma_1^2x'\,^Tx'=\sigma_1^2||x'||_2^2 \ \Rightarrow\ \frac{||y'||_2}{||x'||_2}\le\sigma_1 \end{array} }$

$\displaystyle{(6)\quad ||y||_2^2=||Uy'||_2^2=y'\,^TU^TUy'=y'\,^Ty'=||y'||_2^2\ \Rightarrow\ \frac{||y||_2}{||y'||_2}=1 }$

すなわち

$\displaystyle{(7)\quad \frac{||x'||_2}{||x||_2}\times\frac{||y'||_2}{||x'||_2}\times\frac{||y||_2}{||y'||_2}\le\sigma_1 \ \Rightarrow\ \frac{||y||_2}{||x||_2}=\frac{||Ax||_2}{||x||_2}\le\sigma_1 }$

したがって、線形写像 $y=Ax$ の伝達特性は、行列 $A$ の２ノルム

$\displaystyle{(8)\quad \boxed{||A||_2=\max_{x\ne 0}\frac{||Ax||_2}{||x||_2}=\sigma_1(A)} }$

すなわち行列 $A$ の最大特異値 $\sigma_1(A)$ （行列 $A^TA$ または $AA^T$ の最大固有値の正の平方根）によって測られます。

[2]　入力 $u(t)\in{\rm\bf R}^m$ から出力 $y(t)\in{\rm\bf R}^p$ への漸近安定な $n$ 次系の「伝達特性」をどう測るかを考えます。

図4　ｎ次系

$n$ 次系の入出力応答は次式で表されます。ただし、 $x(0)=0$ とします。

$\displaystyle{(9)\quad y(t)=\int_0^t G(t-\tau)u(\tau)d\tau\quad(G(t)=C\exp(At)B) }$

これをフーリエ変換しますと次式を得ます。

$\displaystyle{(10)\quad \hat{y}(j\omega)=\hat{G}(j\omega)\hat{u}(j\omega)\quad(\hat{G}(j\omega)=C(j\omega I_n-A)^{-1}B) }$

いま時間関数 $x(t)\in{\rm\bf R}^n$ のノルムを次のように定義し、このノルムが有界な時間関数の集合を ${\cal L}_2$ と表記します。

$\displaystyle{(11)\quad ||x(t)||_2=\sqrt{\int_0^\infty x^T(t)x(t)dt}<\infty }$

このとき、次のパーセバルの定理

$\displaystyle{(12)\quad \int_0^\infty x^T(t)x(t)\,dt=\frac{1}{2\pi}\int_{-\infty}^\infty {\hat x}^{*T}(j\omega){\hat x}(j\omega)\,d\omega \quad ({\hat x}(j\omega)=\int_{-\infty}^\infty x(t)e^{-j\omega t}dt) }$

および周波数伝達関数行列 $\hat{G}(j\omega)$ の２ノルムに関する不等式

$\displaystyle{(13)\quad \frac{||\hat{G}(j\omega)\hat{u}(j\omega)||_2}{||\hat{u}(j\omega)||_2} \le||\hat{G}(j\omega)||_2=\sigma_1(\hat{G}(j\omega)) }$

に注意して、次式が成り立ちます。

$(14)\quad \begin{array}{l} \displaystyle{\underbrace{\int_0^\infty y^T(t)y(t)\,dt}_{||y(t)||_2^2} =\frac{1}{2\pi}\int_{-\infty}^\infty \hat{y}^{*T}(j\omega)\hat{y}(j\omega)\,d\omega} \\ \displaystyle{=\frac{1}{2\pi}\int_{-\infty}^\infty \hat{u}^{*T}(j\omega)\hat{G}^{*T}(j\omega)\hat{G}(j\omega)\hat{u}(j\omega)\,d\omega} \\ \displaystyle{\le\frac{1}{2\pi}\int_{-\infty}^\infty ||\hat{G}(j\omega)||_2^2\hat{u}(j\omega)^{*T}\hat{u}(j\omega)\,d\omega} \\ \displaystyle{\le\left(\sup_{\omega\in{\rm\bf R}}||\hat{G}(j\omega)||_2\right)^2\frac{1}{2\pi}\int_{-\infty}^\infty \hat{u}(j\omega)^{*T}\hat{u}(j\omega)\,d\omega} \\ =\displaystyle{\left(\sup_{\omega\in{\rm\bf R}}||\hat{G}(j\omega)||_2\right)^2\underbrace{\int_0^\infty u^T(t)u(t)\,dt}_{||u(t)||_2^2}} \end{array} }$

すなわち

$\displaystyle{(15)\quad \frac{||y(t)||_2}{||u(t)||_2}\le\sup_{\omega\in{\rm\bf R}}||\hat{G}(j\omega)||_2 }$

これに基づいて、 $n$ 次系の ${\bf H}_\infty$ ノルムを次式で定義します。

$\displaystyle{(16)\quad \boxed{\sup_{u\in{\cal L}_2}\frac{||y(t)||_2}{||u(t)||_2} =\sup_{\omega\in{\rm\bf R}}||\hat{G}(j\omega)||_2 =\sup_{\omega\in{\rm\bf R}}\sigma_1(\hat{G}(j\omega))} }$

●いま、任意の複素数 $c$ の特異値は $\sigma_1(c)=|c|$ であることに注意すると、スカラ系の $H_\infty$ ノルムは、ゲイン曲線の最大値であることがわかります。多変数系に対して $\sigma_1(\hat{G}(j\omega))$ をプロットした曲線は $\sigma$ プロットと呼ばれています。

●次の直達項をもつｎ次系を考えます。

図5 直達項をもつｎ次系

このとき、インパルス応答と周波数伝達関数は、それぞれ

$\displaystyle{(17)\quad G(t)=C\exp(At)B+D }$

$\displaystyle{(18)\quad \hat{G}(j\omega)=C(j\omega I_n-A)^{-1}B+D }$

となりますが、上の議論はそのまま成り立ち、 $H_\infty$ ノルムを定義できます。

[3]　漸近安定な $n$ 次系の「伝達特性」を測るもう一つの方法は、すべてのインパルス応答の２乗面積の総和によるものです。 $m$ 入力 $p$ 出力 $n$ 次系のインパルス応答は

$\displaystyle{(19)\quad g_{ij}(t)=C_i\exp(At)B_j\quad(i=1,\cdots,p; j=1,\cdots,m) }$

の $pm$ 個あります。いま任意の行列 $X=[x_{ij}]$ のすべての要素の２乗和 ${\rm tr}XX^T$ はまたは ${\rm tr}X^TX$ で表されることに注意します。すべてのインパルス応答の２乗面積の総和は

$\displaystyle{(20)\quad {\rm tr}(\int_0^\infty G(t)G^T(t)dt)={\rm tr} (C\underbrace{\int_0^\infty \exp(At)BB^T\exp(A^Tt)dt}_{W_c}C^T) }$

または

$\displaystyle{(21)\quad %\begin{array}{l} {\rm tr}(\int_0^\infty G^T(t)G(t)dt)={\rm tr} (B^T\underbrace{\int_0^\infty \exp(A^Tt)C^TC\exp(At)dt}_{W_o}B) %\end{array} }$

のように表されます。ここで、 $W_c$ と $W_o$ はそれぞれリャプノフ方程式

$\displaystyle{(23)\quad W_cA^T+AW_c+BB^T=0 }$

$\displaystyle{(24)\quad W_oA+A^TW_o+C^TC=0 }$

の解であることに注意しておきます。 ${\rm tr}XX^T={\rm tr}X^TX$ に注意して

$\displaystyle{(25)\quad {\rm tr} (CW_cC^T)={\rm tr} (B^TW_oB) }$

の正の平方根を $n$ 次系の ${\bf H}_2$ ノルムと呼びます。

●直達項をもつｎ次系の場合は、インパルス応答

$\displaystyle{(26)\quad g_{ij}(t)=C_i\exp(At)B_j+D_{ij}\quad(i=1,\cdots,p;\ j=1,\cdots,m) }$

の２乗面積の総和は発散してしまうので、 $H_2$ ノルムはできないことに注意します。

Note B21 １次系のパーセバルの定理　

１次系

$\displaystyle{(1)\quad \left\{\begin{array}{l} \dot{x}(t)=-\frac{1}{T}x(t)+\frac{K}{T}u(t) \\ y(t)=x(t) \end{array}\right. }$

のインパルス応答と周波数伝達関数の絶対値は

$\displaystyle{(2)\quad g(t)=\frac{K}{T}e^{-\frac{1}{T}t},\ |\hat{g}(j\omega)|=\frac{K}{\sqrt{1+\omega^2 T^2}} }$

となります。まずパーセバルの定理

$\displaystyle{(3)\quad \int_0^\infty g^2(t)\,dt =\frac{1}{2\pi}\int_{-\infty}^\infty |\hat{g}(j\omega)|^2\,d\omega }$

が成り立つことは次のように確かめられます。

$\displaystyle{(4)\quad \int_0^\infty g^2(t)\,dt =\int_0^\infty \frac{K^2}{T^2}e^{-\frac{2}{T}t}\,dt} =\frac{K^2}{T^2}\left[-\frac{T}{2}e^{-\frac{2}{T}t}\right]_0^\infty =\frac{K^2}{2T} }$

$\displaystyle{(5)\quad \frac{1}{2\pi}\int_{-\infty}^\infty |\hat{g}(j\omega)|^2\,d\omega =\frac{1}{2\pi}\int_{-\infty}^\infty \frac{K^2}{1+\omega^2 T^2}\,d\omega =\frac{1}{2\pi}{K^2}\left[\frac{1}{T}\tan^{-1}\omega T\right]_{-\infty}^\infty =\frac{1}{2\pi}\frac{K^2}{T}\pi }$

また１次系のノルムは、複素数 $\hat{g}(j\omega)$ の特異値がその絶対値であることに注意して

$\displaystyle{(6)\quad \sup_{u\in{\cal L}_2}\frac{||y(t)||_2}{||u(t)||_2} =\sup_{\omega\in{\rm\bf R}}|\hat{g}(j\omega)| =\sup_{\omega\in{\rm\bf R}}\frac{K}{\sqrt{1+\omega^2 T^2}} =K }$

D制約LMI

投稿日時: 2021年10月29日投稿者: cacsd

D制約LMI

[1] 次の命題が成り立ちます。

$\displaystyle{ \begin{array}{lll} &&\lambda(A)\subset {\cal D}=\{s\in{\rm\bf C}: L+sM+s^*M^T<0\ (L=L^T)\}\\ &\Leftrightarrow& \exists X>0:\ L\otimes X+M\otimes(XA)+M^T\otimes(A^TX)<0\\ &\Leftrightarrow& \exists Y>0:\ L\otimes Y+M\otimes(AY)+M^T\otimes(YA^T)<0 \end{array} }$

●十分性（ $\Leftarrow$ ）は、クロネッカ積の公式 $(A\otimes B)(C\otimes D)=(AC\otimes BD)$ を用いて、 $Av=\lambda v$ （ $v^HA^T=\lambda^*v^H$ ）のとき

$\displaystyle{(1)\quad \begin{array}{l} (I_n\otimes v^H)(L\otimes X+M\otimes(XA)+M^T\otimes(A^TX))(I_n\otimes v)\\ =(I_n\otimes v^H)(L\otimes Xv+M\otimes(XAv)+M^T\otimes(A^TXv))\\ =L\otimes v^HXv+M\otimes(v^HXAv)+M^T\otimes(v^HA^TXv)\\ =L\otimes v^HXv+M\otimes(v^HX\lambda v)+M^T\otimes(\lambda^*v^HXv)\\ =\underbrace{v^HXv}_{>0}\underbrace{(L +\lambda M+\lambda^*M^T)}_{<0}<0 \end{array} }$

が成り立つことから出ます。

●必要性（ $\Rightarrow$ ）を示すために、まず $A$ は対角化可能とします。このとき

$\displaystyle{(2)\quad L+\lambda_i M+\lambda_i^*M^T<0\quad (i=1,\cdots,n) }$

ならば、明らかにこれらをブロック対角にもたせた次式が成り立ちます。

$\displaystyle{(3)\quad L\otimes I_n+M\otimes {\rm diag}\{\lambda_1,\cdots,\lambda_n\}+M^T\otimes {\rm diag}\{\lambda_1^*,\cdots,\lambda_n^*\}<0 }$

一般に、 $A$ が対角化可能ではないときは、そのジョルダン分解を

$\displaystyle{(4)\quad \begin{array}{l} A=TJT^{-1}\\ J={\rm diag}\{J_1,\cdots,J_r\}\\ J_i={\rm diag}\{J(\lambda_i,s_{i1}),\cdots,J(\lambda_i,s_{ir_i})\}\in{\rm\bf R}^{m_i\times m_i}\\ (m_1+\cdots+m_r=n;\ s_{i1}+\cdots+s_{ir_i}=m_i) \end{array} }$

ただし、ジョルダン細胞を

$\displaystyle{(5)\quad J(\lambda,k)=\left[\begin{array}{cccc} \lambda & 1 & & 0\\ 0 & \lambda & \ddots & \\ \vdots & \ddots & \ddots & 1 \\ 0 & \cdots & 0 & \lambda \end{array}\right]\in{\rm\bf R}^{k\times k} }$

とします。いま

$\displaystyle{(6)\quad T_\epsilon=\left[\begin{array}{cccc} 1 & 0 & \cdots & 0\\ 0 & \frac{1}{\epsilon} & \ddots & \vdots\\ \vdots & \ddots & \ddots & 0\\ 0 & \cdots & 0 & \frac{1}{\epsilon^{k-1}} \end{array}\right],\ T_\epsilon^{-1}=\left[\begin{array}{cccc} 1 & 0 & \cdots & 0\\ 0 & {\epsilon} & \ddots & \vdots\\ \vdots & \ddots & \ddots & 0\\ 0 & \cdots & 0 & {\epsilon^{k-1}} \end{array}\right] }$

を用いて

$\displaystyle{(7)\quad \begin{array}{l} T_\epsilon J(\lambda,k)T_\epsilon^{-1}\\ =T_\epsilon \left[\begin{array}{cccc} \lambda & 1 & & 0\\ 0 & \lambda & \ddots & \\ \vdots & \ddots & \ddots & 1 \\ 0 & \cdots & 0 & \lambda \end{array}\right] \left[\begin{array}{cccc} 1 & 0 & \cdots & 0\\ 0 & {\epsilon} & \ddots & \vdots\\ \vdots & \ddots & \ddots & 0\\ 0 & \cdots & 0 & {\epsilon^{k-1}} \end{array}\right]\\ =\left[\begin{array}{cccc} 1 & 0 & \cdots & 0\\ 0 & \frac{1}{\epsilon} & \ddots & \vdots\\ \vdots & \ddots & \ddots & 0\\ 0 & \cdots & 0 & \frac{1}{\epsilon^{k-1}} \end{array}\right] \left[\begin{array}{cccc} \lambda & \epsilon & & 0\\ 0 & \lambda\epsilon & \ddots & \\ \vdots & \ddots & \ddots & \epsilon^{k-1} \\ 0 & \cdots & 0 & \lambda\epsilon^{k-1} \end{array}\right]\\ =\left[\begin{array}{cccc} \lambda & \epsilon & & 0\\ 0 & \lambda & \ddots & \\ \vdots & \ddots & \ddots & \epsilon \\ 0 & \cdots & 0 & \lambda \end{array}\right] \end{array} }$

を得ます。したがって、 $A$ の各ジョルダン細胞に対する $T_\epsilon$ をブロック対角に持たせた $\hat{T}_\epsilon$ を用いて

$\displaystyle{(8)\quad \Lambda_\epsilon = \underbrace{\hat{T}_\epsilon T}_{S}A\underbrace{T^{-1}\hat{T}_\epsilon^{-1}}_{S^{-1}} }$

は十分小さな $\epsilon$ に対して対角行列となります。したがって、(3)と同様に次式を得ます。

$\displaystyle{(9)\quad \left.L\otimes I_n+M\otimes \Lambda_\epsilon+M^T\otimes \Lambda_\epsilon^*\right|_{\epsilon\rightarrow0}<0 }$

これから、 $S=\left.\hat{T}_\epsilon T|_{\epsilon\rightarrow0}$ とおき、 $\Lambda_\epsilon=SAS^{-1}$ （ $\Lambda_\epsilon^*=S^{-H}A^TS^H$ ）を用いて

$\displaystyle{(10)\quad \begin{array}{l} (I_n\otimes S^H)(L\otimes I_n+M\otimes(SAS^{-1})+M^T\otimes(S^{-T}A^TS^T)(I_n\otimes S)\\ =(I_n\otimes S^H)(L\otimes S+M\otimes(SA)+M^T\otimes(S^{-H}A^TS^HS))\\ =L\otimes \underbrace{S^HS}_{X}+M\otimes(\underbrace{S^HS}_{X}A)+M^T\otimes(A^T\underbrace{S^HS}_{X})<0 \end{array} }$

が成り立ち、必要性が示されたことになります。

セクタ制約LMI

投稿日時: 2021年10月29日投稿者: cacsd

セクタ制約LMI…Homework

[1] 次の領域 ${\cal D}_3$ を考えます。

図1 領域 ${\cal D}_3$

このとき、次の命題が成り立ちます。

$\displaystyle{(12)\quad \begin{array}{lll} &&\lambda(A)\subset {\cal D}_3=\{s=x+jy\in{\rm\bf C}:\\ && \left[\begin{array}{cc} \sin\theta & \cos\theta \\ -\cos\theta & \sin\theta \end{array}\right]s + \left[\begin{array}{cc} \sin\theta & \cos\theta \\ -\cos\theta & \sin\theta \end{array}\right]^Ts^* <0 \}\nonumber\\ &&\Leftrightarrow \exists X>0:\ \left[\begin{array}{cc} \sin\theta(XA+A^TX) & \cos\theta(XA-A^TX) \\ -\cos\theta(XA-A^TX) & \sin\theta(XA+A^TX) \end{array}\right] <0 \nonumber\\ &&\Leftrightarrow \exists Y>0:\ \left[\begin{array}{cc} \sin\theta(AY+YA^T) & \cos\theta(AY-YA^T) \\ -\cos\theta(AY-YA^T) & \sin\theta(AY+YA^T) \end{array}\right] <0\nonumber \end{array} }$

実際、まず領域 ${\cal D}_3$ が上図を表すことは、次のように確かめられます。

$\displaystyle{(13)\quad \begin{array}{lll} && \left[\begin{array}{cc} \sin\theta & \cos\theta \\ -\cos\theta & \sin\theta \end{array}\right]s + \left[\begin{array}{cc} \sin\theta & \cos\theta \\ -\cos\theta & \sin\theta \end{array}\right]^Ts^* <0 \nonumber\\ &\Leftrightarrow& \left[\begin{array}{cc} (s+s^*)\sin\theta & (s-s^*)\cos\theta \\ -(s-s^*)\cos\theta & (s+s^*)\sin\theta \end{array}\right] <0 \nonumber\\ &\Leftrightarrow& \left[\begin{array}{cc} 2x\sin\theta & 2jy\cos\theta \\ -2jy\cos\theta & 2x\sin\theta \end{array}\right] <0 \nonumber\\ &\Leftrightarrow& x\sin\theta-jy\cos\theta\frac{1}{x\sin\theta}(-j)y\cos\theta<0,\ x\sin\theta<0 \nonumber\\ &\Leftrightarrow& x^2\sin^2\theta-y^2\cos^2\theta>0,\ x<0 \nonumber\\ &\Leftrightarrow& \tan\theta>\frac{|y|}{-x},\ x<0 \nonumber \end{array} }$

次に、十分性は、 $Av=\lambda v$ とすると、次のように確かめられます。

$\displaystyle{(14) \begin{array}{lll} &&\left[\begin{array}{cc} v^{*T} & 0^{T}\\ 0^{T} & v^{*T} \end{array}\right] \left[\begin{array}{cc} \sin\theta(XA+A^TX) & \cos\theta(XA-A^TX) \\ -\cos\theta(XA-A^TX) & \sin\theta(XA+A^TX) \end{array}\right] \left[\begin{array}{cc} v & 0\\ 0 & v \end{array}\right] \nonumber\\ &=& \left[\begin{array}{cc} \sin\theta(v^{*T}X(\lambda v)+(\lambda^*v^{*T})Xv) & \cos\theta(v^{*T}X(\lambda v)-(\lambda^*v^{*T})Xv) \\ -\cos\theta(v^{*T}X(\lambda v)-(\lambda^*v^{*T})Xv) & \sin\theta(v^{*T}X(\lambda v)+(\lambda^*v^{*T})Xv) \end{array}\right] \nonumber\\ &=& \underbrace{ \left[\begin{array}{cc} (\lambda+\lambda^*)\sin\theta & (\lambda-\lambda^*)\cos\theta \\ -(\lambda-\lambda^*)\cos\theta & (\lambda+\lambda^*)\sin\theta \end{array}\right] }_{<0} \underbrace{v^{*T}Xv}_{>0}<0 \nonumber \end{array} }$

さらに、 $y_1=Xx_1,y_2=Xx_2$ とおくと、次が成り立ちます。

$\displaystyle{(15) \begin{array}{lll} && \left[\begin{array}{cc} x_1^T & x_2^T \end{array}\right] \left[\begin{array}{cc} \sin\theta(XA+A^TX) & \cos\theta(XA-A^TX) \\ -\cos\theta(XA-A^TX) & \sin\theta(XA+A^TX) \end{array}\right] \left[\begin{array}{c} x_1\\ x_2 \end{array}\right] <0 \nonumber\\ &&\quad (\forall x_1,x_2\ne0)\nonumber\\ &\Leftrightarrow& \left[\begin{array}{cc} x_1^TX & x_2^TX \end{array}\right]\\ && \left[\begin{array}{cc} \sin\theta(AX^{-1}+X^{-1}A^T) & \cos\theta(AX^{-1}-X^{-1}A^T) \\ -\cos\theta(AX^{-1}-X^{-1}A^T) & \sin\theta(AX^{-1}+X^{-1}A^T) \end{array}\right] \left[\begin{array}{c} Xx_1\\ Xx_2 \end{array}\right] <0 \nonumber\\ &\Leftrightarrow& \left[\begin{array}{cc} y_1^T & y_2^T \end{array}\right] \left[\begin{array}{cc} \sin\theta(AY+YA^T) & \cos\theta(AY-YA^T) \\ -\cos\theta(AY-YA^T) & \sin\theta(AY+YA^T) \end{array}\right] \left[\begin{array}{c} y_1\\ y_2 \end{array}\right] <0 \nonumber\\ &&\quad (\forall y_1,y_2\ne0)\nonumber \end{array} }$

必要性については、あとで述べます。

演習B13-1…Flipped Classroom
$1^\circ$ セクタ制約を調べる次のコードを説明せよ。

MATLAB

%ana_lmi3.m
%-----
 clear all, close all
 A=[0 1;-2 -2]; n=2;
%-----
 setlmis([]);  
 X=lmivar(1,[n 1]);  
%-----
 th=pi/6; sth=sin(th); cth=cos(th); 
 lmi1=newlmi;  
 lmiterm([lmi1 1 1 X],1,sth*A,'s'); %#1:sth*(X*A+A'*X) 
 lmiterm([lmi1 1 2 X],1,cth*A);     %#1:cth*X*A 
 lmiterm([lmi1 1 2 X],-cth*A',1);   %#1:-cth*A'*X 
 lmiterm([lmi1 2 2 X],1,sth*A,'s'); %#1:sth*(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

[2] 次の領域 ${\cal D}={\cal D}_1\cup{\cal D}_2\cup{\cal D}_3$ を考えます。

図2 領域 ${\cal D}={\cal D}_1\cup{\cal D}_2\cup{\cal D}_3$

このとき、次の命題が成り立ちます。

$\displaystyle{(16)\quad \begin{array}{lll} && \lambda(A)\subset {\cal D}={\cal D}_1\cap{\cal D}_2\cap{\cal D}_3\nonumber\\ &&\Leftrightarrow \exists X>0:\nonumber\\ & & \left\{\begin{array}{l} 2\alpha X+XA+A^TX<0 \\ \left[\begin{array}{cc} -rX & XA \\ A^TX & -rX \end{array}\right]<0 \\ \left[\begin{array}{cc} \sin\theta(XA+A^TX) & \cos\theta(XA-A^TX) \\ -\cos\theta(XA-A^TX) & \sin\theta(XA+A^TX) \end{array}\right] <0 \end{array}\right.\nonumber\\ &&\Leftrightarrow \exists Y>0:\nonumber\\ & & \left\{\begin{array}{l} 2\alpha Y+AY+YA^T<0 \\ \left[\begin{array}{cc} -rY & AY \\ YA^T & -rY \end{array}\right]<0 \\ \left[\begin{array}{cc} \sin\theta(AY+YA^T) & \cos\theta(AY-YA^T) \\ -\cos\theta(AY-YA^T) & \sin\theta(AY+YA^T) \end{array}\right] <0 \end{array}\right.\nonumber \end{array} }$

演習B13-2…Flipped Classroom
$1^\circ$ 次を参考にして、 ${\cal D}$ 制約を満たす状態FBを求める関数を作成せよ。

MATLAB

%ana_lmi4.m 
%-----
 clear all, close all
 A=[0 1;-2 -3]; n=2;
%-----
 setlmis([]);  
 X=lmivar(1,[n 1]);  
%-----
 alpha=0.5;
 lmi1=newlmi; 
 lmiterm([lmi1 1 1 X],1,A,'s');     %#1:XA+A'*X   
 lmiterm([lmi1,1,1,X],2*alpha,1);   %#1:2*alpha*X
%
 q=0; r=3;
 lmi2=newlmi; 
 lmiterm([lmi2 1 1 X],-r,1);        %#2:   
 lmiterm([lmi2 1 2 X],1,A);         %#2:2*alpha*X  
 lmiterm([lmi2 2 2 X],-r,1);        %#2: 
%
 th=pi/4; sth=sin(th); cth=cos(th); 
 lmi3=newlmi;  
 lmiterm([lmi3 1 1 X],1,sth*A,'s'); %#3:sth*(X*A+A'*X) 
 lmiterm([lmi3 1 2 X],1,cth*A);     %#3:cth*X*A 
 lmiterm([lmi3 1 2 X],-cth*A',1);   %#3:-cth*A'*X 
 lmiterm([lmi3 2 2 X],1,sth*A,'s'); %#3:sth*(X*A+A'*X) 
%
 lmi4=newlmi;  
 lmiterm([-lmi4 1 1 X],1,1);        %#4:X  
%----- 
 LMIs=getlmis;  
 [tmin,xfeas]=feasp(LMIs);  
 X=dec2mat(LMIs,xfeas,X)
%-----
%eof

円制約LMI

投稿日時: 2021年10月29日投稿者: cacsd

円制約LMI…Homework

[1] 次の領域 ${\cal D}_2$ を考えます。

図1 領域 ${\cal D}_2$

このとき、次の命題が成り立ちます。

$\displaystyle{(8)\quad \begin{array}{lll} &&\lambda(A)\subset {\cal D}_2=\{s=x+jy\in{\rm\bf C}: \left[\begin{array}{cc} -r & s \\ s^* & -r \end{array}\right] <0 \}\nonumber\\ &\Leftrightarrow& \exists X>0:\ \left[\begin{array}{cc} -rX & XA \\ A^TX & -rX \end{array}\right]<0 \nonumber\\ &\Leftrightarrow& \exists Y>0:\ \left[\begin{array}{cc} -rY & AY \\ YA^T & -rY \end{array}\right]<0\nonumber \end{array} }$

●実際、まず領域 ${\cal D}_2$ が上図を表すことは、次のように確かめられます。

$\displaystyle{(9)\quad \begin{array}{lll} && \left[\begin{array}{cc} -r & s \\ s^* & -r \end{array}\right]<0\\ &\Leftrightarrow& (-r)-s\frac{1}{(-r)}s^*<0,\ -r<0 \nonumber\\ &\Leftrightarrow& r^2-(x+jy)(x-jy)>0 \nonumber\\ &\Leftrightarrow& x^2+y^2<r^2 \nonumber\end{array}}$

●次に、十分性は、 $Av=\lambda v$ とすると、次のように確かめられます。

$\displaystyle{(10)\quad \begin{array}{lll} &&\left[\begin{array}{cc}v^{*T}&0^{T}\\0^{T}&v^{*T}\end{array}\right] \left[\begin{array}{cc}-rX&XA\\A^TX&-rX\end{array}\right] \left[\begin{array}{cc}v&0\\0&v\end{array}\right]\nonumber\\ &=& \left[\begin{array}{cc} -rv^{*T}Xv&v^{*T}X(\lambda v)\\ (\lambda^*v^{*T})Xv&-rv^{*T}Xv& \end{array}\right]\nonumber\\ &=& \underbrace{ \left[\begin{array}{cc} -r&\lambda\\\lambda^*&-r \end{array}\right] }_{<0} \underbrace{v^{*T}xv}_{>0}<0 \nonumber \end{array} }$

さらに、 $y_1=Xx_1,y_2=Xx_2$ とおくと、次が成り立ちます。

$\displaystyle{(11)\quad \begin{array}{lll} && \left[\begin{array}{cc} x_1^T & x_2^T \end{array}\right] \left[\begin{array}{cc} -rX & XA \\ A^TX & -rX \end{array}\right] \left[\begin{array}{c} x_1\\ x_2 \end{array}\right] <0 \quad (\forall x_1,x_2\ne0)\nonumber\\ &\Leftrightarrow& \left[\begin{array}{cc} x_1^TX & x_2^TX \end{array}\right] \left[\begin{array}{cc} -rX^{-1} & AX^{-1} \\ X^{-1}A^T & -rX^{-1} \end{array}\right] \left[\begin{array}{c} Xx_1\\ Xx_2 \end{array}\right] <0 \nonumber\\ &\Leftrightarrow& \left[\begin{array}{cc} y_1^T & y_2^T \end{array}\right] \left[\begin{array}{cc} -rY & AY \\ YA^T & -rY \end{array}\right] \left[\begin{array}{c} y_1\\ y_2 \end{array}\right] <0 \quad (\forall y_1,y_2\ne0)\nonumber \end{array} }$

必要性については、あとで述べます。

演習B12…Flipped Classroom
$1^\circ$ 円制約を調べる次のコードを説明せよ。

MATLAB

%ana_lmi2.m
%-----
 A=[0 1;-2 -3]; n=2;
%-----
 setlmis([]);  
 X=lmivar(1,[n 1]);  
%-----
 q=0; r=1.5;
 lmi1=newlmi; 
 lmiterm([lmi1 1 1 X],-r,1);     %#1:   
 lmiterm([lmi1 1 2 X],1,A);      %#1:2*alpha*X  
 lmiterm([lmi1 2 2 X],-r,1);     %#1: 
%
 lmi2=newlmi;  
 lmiterm([-lmi2 1 1 X],1,1);     %#2:X  
%-----
 LMIs=getlmis;  
 [tmin,xfeas]=feasp(LMIs);  
 X=dec2mat(LMIs,xfeas,X)
%-----
%eof

$2^\circ$ 次の領域を表す ${\rm\bf LMI}$ を導出し、これをチェックするプログラムを作成せよ。

$\displaystyle{(17)\quad \lambda(A)\subset{\cal D}_4=\{s=x+jy\in{\rm\bf C}: \frac{x^2}{a^2}+\frac{y^2}{b^2}<1\} }$

$\displaystyle{(18)\quad \begin{array}{lll} &&\frac{x^2}{a^2}+\frac{y^2}{b^2}<1\Leftrightarrow\frac{(\frac{s+s^*}{2})^2}{a^2}+\frac{(\frac{s-s^*}{2j})^2}{b^2}<1 \\ &&\Leftrightarrow\ b^2(s+s^*)^2-a^2(s-s^*)^2<4a^2b^2\\ &&\Leftrightarrow\ (b(s+s^*)+a(s-s^*))(b(s+s^*)-a(s-s^*))<4a^2b^2 \\ &&\Leftrightarrow\ -4a^2b^2+((b+a)s+(b-a)s^*))((b-a)s+(b+a)s^*))<0 \\ &&\Leftrightarrow\ \left[\begin{array}{cc} -4a^2b^2 & (b+a)s+(b-a)s^* \\ (b-a)s+(b+a)s^* & -1 \end{array}\right]<0 \end{array} }$

MATLAB

%ana_lmi2a.m
%-----
 A=[0 1;-1 -2]; n=2;
%-----
 setlmis([]);  
 X=lmivar(1,[n 1]);  
%-----
 a=2; b=1;
%  LMI1=-[-4*a^2*b^2*X (b+a)*X*A+(b-a)*A'*X;
%         (b-a)*X*A+(b+a)*A'*X -X];   
 lmi1=newlmi; 
 lmiterm([lmi1 1 1 X],-4*a^2*b^2,1);  %#1:-4*a^2*b^2*X  
 lmiterm([lmi1 1 2 X],b+a,A);         %#1:(b+a)*X*A
 lmiterm([lmi1 1 2 X],(b-a)*A',1);    %#1:(b-a)*A'*X
 lmiterm([lmi1 2 2 X],-1,1);          %#1:-X
%
 lmi2=newlmi;  
 lmiterm([-lmi2 1 1 X],1,1);     %#2:X  
%-----
 LMIs=getlmis;  
 [tmin,xfeas]=feasp(LMIs);  
 X=dec2mat(LMIs,xfeas,X)
%-----
%eof

$3^\circ$ 　次の領域を表す ${\rm\bf LMI}$ を導出し、これをチェックするプログラムを作成せよ。

$\displaystyle{(19)\quad \lambda(A)\subset{\cal D}_5=\{s=x+jy\in{\rm\bf C}: \frac{x^2}{a^2}-\frac{y^2}{b^2}>1, x<-a\} }$

$\displaystyle{(20)\quad \begin{array}{lll} &&\frac{x^2}{a^2}-\frac{y^2}{b^2}>1,\ x<-a\ \Leftrightarrow\ \frac{(\frac{s+s^*}{2})^2}{a^2}-\frac{(\frac{s-s^*}{2j})^2}{b^2}>1,\ \frac{s+s^*}{2}<-a}\\ &&\Leftrightarrow\ b^2(s+s^*)^2+a^2(s-s^*)^2>4a^2b^2,\ s+s^*+2a<0\\ &&\Leftrightarrow\ b^2(s+s^*+2a)(s+s^*-2a)+a^2(s-s^*)^2>0,\ s+s^*+2a<0\\ &&\Leftrightarrow\ b(s+s^*-2a)-\frac{-a^2(s-s^*)^2}{b(s+s^*+2a)}<0,b(s+s^*+2a)<0\\ &&\Leftrightarrow\ \left[\begin{array}{cc} b(s+s^*+2a) & a(s-s^*) \\ -a(s-s^*) & b(s+s^*-2a) \end{array}\right]<0 \end{array} }$

MATLAB

%ana_lmi2b.m
%-----
 A=[0 1;-1 -2]; n=2;
%-----
 setlmis([]);  
 X=lmivar(1,[n 1]);  
%-----
 a=2; b=sqrt(2); f=sqrt(a^2+b^2);
%  LMI1=-[b*(X*A+A'*X+2a*X) a*(X*A-A'*X);
%        -a*(X*A-A'*X) b*(X*A+A'*X-2a*X)];
 lmi1=newlmi; 
 lmiterm([lmi1 1 1 X],b,A,'s');       %#1:b*(X*A+A'*X)
 lmiterm([lmi1 1 1 X],2*a*b,1);       %#1:2*a*b*X
 lmiterm([lmi1 1 2 X],a,A);           %#1:a*X*A
 lmiterm([lmi1 1 2 X],-a*A',1);       %#1:-a*A'*X
 lmiterm([lmi1 2 2 X],b,A,'s');       %#1:b*(X*A+A'*X)
 lmiterm([lmi1 2 2 X],-2*a*b,1);      %#1:-2*a*b*X
%
 lmi2=newlmi;  
 lmiterm([-lmi2 1 1 X],1,1);     %#2:X  
%-----
 LMIs=getlmis;  
 [tmin,xfeas]=feasp(LMIs);  
 X=dec2mat(LMIs,xfeas,X)
%-----
%eof

α制約LMI

投稿日時: 2021年10月29日投稿者: cacsd

α制約LMI…Homework

[0] 自由系 $\dot{x}=Ax$ が漸近安定であるための必要十分条件は、行列 $A$ の固有値がすべて複素左平面 ${\rm\bf LHP}$ に存在することでした。でも平衡状態 $x=0$ に復帰するその過程はいつも好ましいとは言えません。たとえば、平衡状態 $x=0$ に復帰する速さが適切であること、その過程が振動的でないこと、そして適切なサンプリング周期をもつ（絶対値が極端に大きな固有値をもたない）ことが必要です。そこで行列 $A$ の固有値が含まれる望ましい領域 $\cal{D}$

$\lambda(A)\subset{\cal D}\subset{\rm\bf LHP}$

を考え、そのための条件を線形行列不等式 ${\rm\bf LMI}$ (Linear Matrix Inequality)で表します。

●準備として、次の命題をチェックします。

$\displaystyle{(1)\quad \begin{array}{lll} && \left[\begin{array}{cc} P & M \\ M^T & Q \end{array}\right]<0\\ &\Leftrightarrow& P-MQ^{-1}M^T<0,\ Q<0\\ &\Leftrightarrow& P<0,\ Q-M^TP^{-1}M<0 \end{array} }$

ここに出てくる $P-MQ^{-1}M^T$ または $Q-M^TP^{-1}M$ はシュール補元(Shure Complement)と呼ばれています。

実際、次のように示されます。

$\displaystyle{(2)\quad \begin{array}{lll} & & \underbrace{ \left[\begin{array}{cc} P & M \\ M^T & Q \end{array}\right] }_{R} \\ &=& \left[\begin{array}{cc} I & MQ^{-1} \\ 0 & I \end{array}\right] \left[\begin{array}{cc} P-MQ^{-1}M^T & 0 \\ 0 & Q \end{array}\right] \left[\begin{array}{cc} I & 0 \\ Q^{-1}M^T & I \end{array}\right] \\ &=& \underbrace{ \left[\begin{array}{cc} I & 0 \\ M^TP^{-1} & I \end{array}\right] }_{T^T} \underbrace{ \left[\begin{array}{cc} P & 0 \\ 0 & Q-M^TP^{-1}M \end{array}\right] }_{S} \underbrace{ \left[\begin{array}{cc} I & P^{-1}M \\ 0 & I \end{array}\right] }_{T} \end{array} }$

これから

$\displaystyle{(3)\quad x^TRx=\underbrace{x^TT^T}_{y^T}S\underbrace{Tx}_{y}=y^TSy}$

を得て、次が成り立ちます。

$\displaystyle{(4)\quad \forall x\ne0: x^TRx<0\Leftrightarrow \forall y=Tx\ne0: y^TSy<0}$

[1] 次の領域 ${\cal D}_1$ を考えます。

図1 領域 ${\cal D}_1$

このとき、次の命題が成り立ちます。

$\displaystyle{ \begin{array}{lll} &&\lambda(A)\subset {\cal D}_1=\{s=x+jy\in{\rm\bf C}: 2\alpha+s+s^*<0 \}\\ &\Leftrightarrow& \exists X>0:\ 2\alpha X+XA+A^TX<0\\ &\Leftrightarrow& \exists Y>0:\ 2\alpha Y+AY+YA^T<0 \end{array} }$

●実際、まず領域 ${\cal D}_1$ が上図を表すことは、次のように確かめられます。

$\displaystyle{(5)\quad \begin{array}{lll} &&2\alpha+s+s^*<0\nonumber\\ &\Leftrightarrow& 2\alpha+(x+jy)+(x-jy)<0\nonumber\\ &\Leftrightarrow& x<-\alpha\nonumber \end{array} }$

●次に、十分性は、 $Av=\lambda v$ とすると、次のように確かめられます（*は複素共役）。

$\displaystyle{(6)\quad \begin{array}{lll} &&v^{*T}(2\alpha X+XA+A^TX)v\nonumber\\ &=&2\alpha v^{*T}Xv+v^{*T}X(\lambda v)+(\lambda^*v^{*T})Xv \nonumber\\ &=&(2\alpha+\lambda+\lambda^*)v^{*T}Xv \nonumber\\ &=&2\underbrace{(\alpha+{\rm Re}(\lambda))}_{<0}\underbrace{v^{*T}Xv}_{>0}<0 \nonumber \end{array} }$

さらに、 $y=Xx$ とおくと、次が成り立ちます。

$\displaystyle{(7)\quad \begin{array}{lll} && x^T(2\alpha X+XA+A^TX)x<0\quad (\forall x\ne0) \nonumber\\ &\Leftrightarrow& \underbrace{x^TX}_{y^T} (2\alpha \underbrace{X^{-1}}_{Y}+A\underbrace{X^{-1}}_{Y}+\underbrace{X^{-1}}_{Y}A^T) \underbrace{Xx}_y<0\quad (\forall y\ne0) \nonumber \end{array} }$

必要性については、あとで述べます。

演習B11…Flipped Classroom
$1^\circ$ α制約を調べる次のコードを説明せよ。

MATLAB

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

Note B11 Inversion Lemma
$\displaystyle{(1)\quad \left[\begin{array}{cc} M_{11} & M_{12} \\ M_{21} & M_{22} \end{array}\right] = \left[\begin{array}{cc} I & M_{12}M_{22}^{-1} \\ 0 & I \end{array}\right] \left[\begin{array}{cc} M_{11}-M_{12}M_{22}^{-1}M_{21} & 0 \\ 0 & M_{22} \end{array}\right] \left[\begin{array}{cc} I & 0 \\ M_{22}^{-1}M_{21} & I \end{array}\right] }$

$\displaystyle{(2)\quad \left[\begin{array}{cc} M_{11} & M_{12} \\ M_{21} & M_{22} \end{array}\right]^{-1} = \left[\begin{array}{cc} \Phi^{-1} & -\Phi^{-1}M_{12}M_{22}^{-1} \\ -M_{22}^{-1}M_{21}\Phi^{-1} & M_{22}^{-1}+M_{22}^{-1}M_{21}\Phi^{-1}M_{12}M_{22}^{-1} \end{array}\right] }$

$\displaystyle{(3)\quad \Phi\stackrel{\rm def}{=}M_{11}-M_{12}M_{22}^{-1}M_{21}\ (M_{22}\ {\rm nonsingular}) }$

および

$\displaystyle{(4)\quad \left[\begin{array}{cc} M_{11} & M_{12} \\ M_{21} & M_{22} \end{array}\right] = \left[\begin{array}{cc} I & 0 \\ M_{21}M_{11}^{-1} & I \end{array}\right] \left[\begin{array}{cc} M_{11} & 0 \\ 0 & M_{22}-M_{21}M_{11}^{-1}M_{12} \end{array}\right] \left[\begin{array}{cc} I & M_{11}^{-1}M_{12} \\ 0 & I \end{array}\right] }$

$\displaystyle{(5)\quad \left[\begin{array}{cc} M_{11} & M_{12} \\ M_{21} & M_{22} \end{array}\right]^{-1} = \left[\begin{array}{cc} M_{11}^{-1}+M_{11}^{-1}M_{12}\Psi^{-1}M_{21}M_{11}^{-1} & -M_{11}^{-1}M_{12}\Psi^{-1} \\ -\Psi^{-1}M_{21}M_{11}^{-1} & \Psi^{-1} \end{array}\right] }$

$\displaystyle{(6)\quad \Psi\stackrel{\rm def}{=}M_{22}-M_{21}M_{11}^{-1}M_{12}\ (M_{11}\ {\rm nonsingular}) }$

が成り立つことから、次を得ます。

$\displaystyle{(7)\quad \Phi^{-1}=(M_{11}-M_{12}M_{22}^{-1}M_{21})^{-1}= M_{11}^{-1}+M_{11}^{-1}M_{12}\Psi^{-1}M_{21}M_{11}^{-1} }$

$\displaystyle{(8)\quad \Psi^{-1}=(M_{22}-M_{21}M_{11}^{-1}M_{12})^{-1}= M_{22}^{-1}+M_{22}^{-1}M_{21}\Phi^{-1}M_{12}M_{22}^{-1} }$

特に、次が成り立ちます。

$\displaystyle{(9)\quad\boxed{ (A-LBR)^{-1}=A^{-1}+A^{-1}L(B^{-1}-RA^{-1}L)^{-1}RA^{-1} }}$

線形計画問題

投稿日時: 2021年10月18日投稿者: cacsd

線形計画問題

[1] 　次の線形計画問題を考えます。

$\displaystyle{\boxed{ \begin{array}{l} {\rm\bf maxmize}\ f(x,y)=x+y\\ {\rm\bf subject\ to}\ \left\{\begin{array}{l} x\ge0\\ y\ge0\\ 5x+6y\le30\\ 3x+2y\le12 \end{array}\right. \end{array}} }$

これは目的関数

$\displaystyle{(1)\quad f(x,y)=x+y}$

を制約条件

$\displaystyle{(2)\quad \left\{\begin{array}{l} x\ge0\\ y\ge0\\ 5x+6y\le30\\ 3x+2y\le12 \end{array}\right. }$

の下で最小化する決定変数 $x,y$ を求める問題として定式化されています。制約条件を満たす $x,y$ の集合を実行可能領域と呼びます。この問題は図1のように図示され、実行可能領域は凸多角形となり、これを通過する直線 $y=-x+z$ のうち最大の $y$ 切片を与える端点として求められます。

図1 線形計画問題の例

●この問題を制約条件付最小化問題

$\displaystyle{\boxed{ \begin{array}{l} {\rm\bf minimize}\ f(x,y)=-(x+y)\\ {\rm\bf subject\ to}\ \left\{\begin{array}{l} x\ge0\\ y\ge0\\ 5x+6y\le30\\ 3x+2y\le12 \end{array}\right. \end{array}} }$

として、MATLABで解く場合のコードを次に示します。

MATLAB

%lp.m
%-----
 clear all close all
 x=-1:10; y=-1:10; 
 plot(x,-5/6*x+5,'b',x,-3/2*x+6,'b',x,-x+3,'r')
 axis([0 6 0 6]),grid, hold on
%------ LMIの設定
 setlmis([]); 
 x=lmivar(1,[1 1]); y=lmivar(2,[1 1]); 
 lmi1=newlmi;                   %#1: x>=0
 lmiterm([-lmi1 1 1 x],1,1);    %x in RHS 
 lmi2=newlmi;                   %#2: y>=0 
 lmiterm([-lmi2 1 1 y],1,1);    %y in RHS 
 lmi3=newlmi;                   %#3: 5x+6y-30<=0
 lmiterm([lmi3 1 1 x],5,1);     %5x in LHS
 lmiterm([lmi3 1 1 y],6,1);     %6y in LHS
 lmiterm([lmi3 1 1 0],-30);     %-30 in LHS
 lmi4=newlmi;                   %#3: 3x+2y-12<=0
 lmiterm([lmi4 1 1 x],3,1);     %3x in LHS
 lmiterm([lmi4 1 1 y],2,1);     %2y in LHS
 lmiterm([lmi4 1 1 0],-12);     %-12 in LHS
 LMIs=getlmis;  
%----- 実行可能解の存在
 [tmin,xfeas]=feasp(LMIs); tmin
 xf=dec2mat(LMIs,xfeas,x), yf=dec2mat(LMIs,xfeas,y)
 plot(xf,yf,'*')
%----- 目的関数の設定と最小解求解
 cobj=-[1 1]; 
 [cost,xopt]=mincx(LMIs,cobj); cost
 xs=dec2mat(LMIs,xopt,x), ys=dec2mat(LMIs,xopt,y)
 plot(xs,ys,'o')
%-----
%eof

実行可能解の存在をチェックする関数feapの実行結果は次のようになり、 $t<0$ の値が得られているので、一つの実行可能解が示されています。


 Solver for LMI feasibility problems L(x) < R(x)
    This solver minimizes  t  subject to  L(x) < R(x) + t*I
    The best value of t should be negative for feasibility

 Iteration   :    Best value of t so far 
 
     1                       -1.203989

 Result:  best value of t:    -1.203989
          f-radius saturation:  0.000% of R =  1.00e+09
 

tmin =
   -1.2040

xf =
    1.2781

yf =
    3.3444

制約条件付最小化問題を解く関数mincxの実行結果は次のようになり、最小解が示されています。


 Solver for linear objective minimization under LMI constraints 

 Iterations   :    Best objective value so far 
     1                  -5.075328
***                 new lower bound:    -5.611584
     2                  -5.211637
***                 new lower bound:    -5.312330
     3                  -5.248462
***                 new lower bound:    -5.298322
 Result:  feasible solution of required accuracy
          best objective value:    -5.248462
          guaranteed relative accuracy:  9.50e-03
          f-radius saturation:  0.000% of R =  1.00e+09 
 
cost =
   -5.2485

xs =
    1.4914

ys =
    3.7571

LPV制御

投稿日時: 2021年10月18日投稿者: cacsd

LPV制御…Homework

[1] 回転体は次の運動方程式で表されます。

$\displaystyle{(1)\quad \left\{\begin{array}{l} J_1\dot{\omega}_1(t)-\omega_2(t)\Omega(t)(J_2-J_3)=\tau_1(t) \\ J_2\dot{\omega}_2(t)-\omega_1(t)\Omega(t)(J_3-J_1)=\tau_2(t) \end{array}\right. }$

ここで、次のパラメータ変動を想定します。

$\displaystyle{(2)\quad \Omega_{min}<\Omega(t)<\Omega_{max} }$

次の状態方程式を得ます。

$\displaystyle{(3) \underbrace{ \left[\begin{array}{c} \dot{\omega}_1(t) \\ \dot{\omega}_2(t) \end{array}\right] }_{\dot{x}(t)} = \underbrace{ \left[\begin{array}{cc} 0 & \Omega(t)\frac{J_2-J_3}{J_1} \\ \Omega(t)\frac{J_3-J_1}{J_2} & 0 \end{array}\right] }_{A(\Omega(t))} \underbrace{ \left[\begin{array}{c} \omega_1(t) \\ \omega_2(t) \end{array}\right] }_{x(t)} + \underbrace{ \left[\begin{array}{cc} \frac{1}{J_1} & 0 \\ 0 & \frac{1}{J_2} \end{array}\right] }_{B} \underbrace{ \left[\begin{array}{c} \tau_1(t) \\ \tau_2(t) \end{array}\right] }_{u(t)} }$

●そのために、 $A(\Omega(t))$ は次式で表されることに着目します。

$\displaystyle{(5)\quad A(\Omega(t))= \underbrace{\frac{\Omega_{max}-\Omega(t)}{\Omega_{max}-\Omega_{min}}}_{p_1(\Omega(t))} \underbrace{A(\Omega_{min})}_{A_1} + \underbrace{\frac{\Omega(t)-\Omega_{min}}{\Omega_{max}-\Omega_{min}}}_{p_2(\Omega(t))} \underbrace{A(\Omega_{max})}_{A_2} }$

ただし

$\displaystyle{(6)\quad p_1(\Omega(t))+p_2(\Omega(t))=1 }$

このとき上の状態方程式は、端点(vertex)モデルとよばれる

$\displaystyle{(7)\quad \dot{x}(t)=A_1x(t)+Bu(t),\ \dot{x}(t)=A_2x(t)+Bu(t) }$

を、係数 $p_1(\Omega),p_2(\Omega)$ によって重み付けて

$\displaystyle{(8)\quad \dot{x}(t)=\underbrace{(p_1(\Omega(t))A_1+p_2(\Omega(t))A_2)}_{A(\Omega(t))}x(t)+Bu(t) }$

のように得られることが分かります。これをLPV(Linear Parameter Varying)モデルと呼びます。

[2] 状態フィートバックによるパラメータ変動抑制について考えます。そのために端点モデルに対応した端点コントローラとよばれる

$\displaystyle{(9)\quad u(t)=-F_1x(t),\ u(t)=-F_2x(t) }$

を、係数 $p_1(\Omega(t)),p_2(\Omega(t))$ によって重み付けた

$\displaystyle{(10)\quad u(t)=-\underbrace{(p_1(\Omega(t))F_1+p_2(\Omega(t))F_2)}_{F(\Omega(t))}x(t) }$

のようなLPV(Linear Parameter Varying)制御を考えます。これは次のように図示されます。

図3 LPV制御の枠組み

この閉ループ系は次式で表されます。

$\displaystyle{(11)\quad \dot{x}(t)=\underbrace{(p_1(\Omega(t))(A_1-BF_1)+p_2(\Omega(t))(A_2-BF_2))}_{A_F(\Omega(t))}x(t) }$

ここで、 $A_1-BF_1$ と $A_2-BF_2$ は安定行列となるように設計しますが、閉ループ系の時変系としての安定性が保証されるかの検討が必要になります。

[3] そのための準備として、次の時変自由系の漸近安定性を考えます。

$\displaystyle{(12)\quad \dot{x}(t)=A(t)x(t) }$

平衡状態 $x=0$ の近傍 ${\cal X}$ において、リャプノフ関数とよばれる

$\displaystyle{(13)\quad \left\{\begin{array}{ll} V(x)>0 & (x\in{\cal X}, x\ne0) \\ V(x)=0 & (x=0) \end{array}\right. }$

を構成して

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

を示すことができれば、平衡状態は漸近安定となることが知られています。リャプノフ関数の存在は十分条件であること、またリャプノフ関数も時変でなく時不変としていることに注意してください。

●ちなみに、時不変自由系

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

の漸近安定性 $\lambda(A)\subset{\bf LHP}$ のLMI条件は

$\displaystyle{(16)\quad \exists X>0:XA+A^TX<0\Leftrightarrow\exists Y>0:AY+YA^T<0 }$

でした（前者がX-LMI、後者がY-LMI）。このときは

$\displaystyle{(17)\quad V(x(t))=x^T(t)Xx(t)}$

と選べば

$\displaystyle{(18)\quad \frac{d}{dt}V(x(t))=\dot{x}^T(t)Xx(t)+x^T(t)X\dot{x}(t)=x^T(t)(A^TX+XA)x(t)<0 }$

となっていることが分かります。

[4] いま、上の $A_1-BF_1$ と $A_2-BF_2$ を個別に安定行列としたとすると

$\displaystyle{(19)\quad \exists X_1>0:X_1(A_1-BF_1)+(A_1-BF_1)^TX_1<0 }$
$\displaystyle{(20)\quad \exists X_2>0:X_2(A_2-BF_2)+(A_2-BF_2)^TX_2<0 }$

から、個別にリャプノフ関数の候補

$\displaystyle{(21)\quad V_1(x(t))=x^T(t)X_1x(t)}$
$\displaystyle{(22)\quad V_2(x(t))=x^T(t)X_2x(t)}$

が得られますが、閉ループ系

$\displaystyle{(23)\quad\dot{x}(t)=A_F(t)x(t)}$

のリャプノフ関数を構成したことにはなりません。一方

$\displaystyle{(24)\quad \exists X>0: \begin{array}{l} X(A_1-BF_1)+(A_1-BF_1)^TX<0\\ X(A_2-BF_2)+(A_2-BF_2)^TX<0 \end{array} }$

となる共通の $X>0$ を見つけることができれば

$\displaystyle{(25)\quad V(x(t))=x^T(t)Xx(t) }$

に対して、LPV制御では

$\displaystyle{(26)\quad \begin{array}{l} \frac{d}{dt}V(x(t))=x^T(t)(A_F^T(\Omega(t))X+XA_F(\Omega(t)))x(t)\\\ =x^T(t)( ( ( p_1(\Omega(t))(A_1-BF_1)+p_2(\Omega(t))(A_2-BF_2) )^TX\\ +X( p_1(\Omega(t))(A_1-BF_1)+p_2(\Omega(t))(A_2-BF_2) ) )x(t)\\ =p_1(\Omega(t)) \underbrace{x^T(t)((A_1-BF_1)^TX+X(A_1-BF_1))x(t)}_{<0}\\ +p_2(\Omega(t)) \underbrace{x^T(t)((A_2-BF_2)^TX+X(A_2-BF_2))x(t)}_{<0}<0 \end{array} }$

となって、リャプノフ関数を構成したことになります。

[5] LPVモデル

$\displaystyle{(27)\quad \left\{\begin{array}{l} \dot{x}(t)=\underbrace{(p_1(\Omega(t))(A_1-BF_1)+p_2(\Omega(t))(A_2-BF_2))}_{A_F(t)}x(t)+Bu(t) \\ y(t)=Cx(t)+Du(t) \end{array}\right. }$

に対して、各端点モデルの閉ループ系において、 $H_\infty$ ノルムは $\gamma$ より小であるとすると

$\displaystyle{(28)\quad \exists X_1>0:\ \left[\begin{array}{ccc} (A_1-BF_1)^TX_1+X_1(A_1-BF_1) & X_1B & C^T \\ B^TX_1 & -\gamma I & D^T \\ C & D & -\gamma　I \end{array}\right]<0 }$
$\displaystyle{(29)\quad \exists X_2>0:\ \left[\begin{array}{ccc} (A_2-BF_2)^TX_2+X_2(A_2-BF_2) & X_2B & C^T \\ B^TX_2 & -\gamma I & D^T \\ C & D & -\gamma　I \end{array}\right]<0 }$

が成り立ちます。いま、共通解 $X_1=X_2=X>0$ が求まっているとしますと

$\displaystyle{(30)\quad \left[\begin{array}{ccc} (A_1-BF_1)^TX+X(A_1-BF_1) & XB & C^T \\ B^TX & -\gamma^2 I & D^T \\ C & D &-I \end{array}\right]<0 }$
$\displaystyle{ \Leftrightarrow \left[\begin{array}{cc} (A_1-BF_1)^TX+X(A_1-BF_1) & XB \\ B^TX & 0 \end{array}\right] < \left[\begin{array}{cc} C^T & 0 \\ D^T & I \end{array}\right] \left[\begin{array}{cc} -I & 0 \\ 0 & \gamma^2 I \end{array}\right] \left[\begin{array}{cc} C & D \\ 0 & I \end{array}\right] }$
$\displaystyle{(31)\quad \left[\begin{array}{ccc} (A_1-BF_1)^TX+X(A_1-BF_1) & XB & C^T \\ B^TX & -\gamma^2 I & D^T \\ C & D &-I \end{array}\right]<0 }$
$\displaystyle{ \Leftrightarrow \left[\begin{array}{cc} (A_2-BF_2)^TX+X(A_2-BF_2) & XB \\ B^TX & 0 \end{array}\right] < \left[\begin{array}{cc} C^T & 0 \\ D^T & I \end{array}\right] \left[\begin{array}{cc} -I & 0 \\ 0 & \gamma^2 I \end{array}\right] \left[\begin{array}{cc} C & D \\ 0 & I \end{array}\right] }$

に注意して、LPV制御では

$\displaystyle{(32)\quad \dot{V}(x(t))=\frac{d}{dt}(x^T(t)Xx(t))= }$
$\displaystyle{ \left[\begin{array}{c} x(t) \\ u(t) \end{array}\right]^T \left[\begin{array}{cc} \begin{array}{c} ( p_1(\Omega(t))(A_1-BF_1)+p_2(\Omega(t))(A_2-BF_2) )^TX\\ +X( p_1(\Omega(t))(A_1-BF_1)+p_2(\Omega(t))(A_2-BF_2) ) \end{array} & XB \\ B^TX & 0 \end{array}\right] \left[\begin{array}{c} x(t) \\ u(t) \end{array}\right] }$
$\displaystyle{ =p_1(\Omega(t)) \left[\begin{array}{c} x(t) \\ u(t) \end{array}\right]^T \left[\begin{array}{cc} (A_1-BF_1)^TX+X(A_1-BF_1) & XB \\ B^TX & 0 \end{array}\right] \left[\begin{array}{c} x(t) \\ u(t) \end{array}\right] }$
$\displaystyle{ +p_2(\Omega(t)) \left[\begin{array}{c} x(t) \\ u(t) \end{array}\right]^T \left[\begin{array}{cc} (A_2-BF_2)^TX+X(A_2-BF_2) & XB \\ B^TX & 0 \end{array}\right] \left[\begin{array}{c} x(t) \\ u(t) \end{array}\right] }$
$\displaystyle{ < \underbrace{(p_1(\Omega(t))+p_2(\Omega(t)))}_1 \left[\begin{array}{c} x(t) \\ u(t) \end{array}\right]^T \left[\begin{array}{cc} C & D \\ 0 & I \end{array}\right]^T \left[\begin{array}{cc} - I & 0 \\ 0 & \gamma^2 I \end{array}\right] \left[\begin{array}{cc} C & D \\ 0 & I \end{array}\right] \left[\begin{array}{c} x(t) \\ u(t) \end{array}\right] }$
$\displaystyle{ =\gamma^2 u^T(t)u(t)-y^T(t)y(t) }$

を得ます。すなわち時変系としての閉ループ系において、消散性の条件

$\displaystyle{(33)\quad \underbrace{\frac{d}{dt}(x^TXx)}_{\dot{V}(x)} <\underbrace{\gamma^2 u^Tu-y^Ty}_{s(u,y)} }$

が成り立ち

$\displaystyle{(34)\quad \sup_{u\in{\cal L}_2}\frac{||y(t)||_2}{||u(t)||_2}<\gamma }$

が成り立ちます。これは時不変系の場合には $H_\infty$ ノルムですが、時変系の場合は ${\cal L}_2$ ゲインとよびます。

(28)と(29)の共通解 $X_1=X_2=X>0$ が求まれば、(24)は自動的に満足されるので、閉ループ系の時変系としての安定性が保証されることになります。

[6] いま、次の２ポート表現かつポリトピック表現された $n$ 次系を考えます。

$\displaystyle{(35)\quad P: \left\{\begin{array}{l} \dot{x}(t)=\underbrace{(p_1(\Omega(t))A_1+p_2(\Omega(t))A_2)}_{A(\Omega(t))}x(t)+B_1u_1(t)+B_2u_2(t) \\ \underbrace{ \left[\begin{array}{c} y(t)\\ u(t) \end{array}\right] }_{y_1(t)} = \underbrace{ \left[\begin{array}{c} C\\ 0 \end{array}\right] }_{C_1} x(t)+ \underbrace{ \left[\begin{array}{c} 0\\ 0 \end{array}\right] }_{D_{11}}u_1(t)+ \underbrace{ \left[\begin{array}{c} 0\\ I \end{array}\right] }_{D_{12}} u_2(t) \\ y_2(t)=x(t) \end{array}\right. }$

これに対する状態フィードバック

$\displaystyle{(36)\quad K: u_2(t)=-\underbrace{(p_1(\Omega(t))F_1+p_2(\Omega(t))F_2)}_{F(\Omega(t))}y_2(t) }$

による閉ループ系

$\displaystyle{(37)\quad P_{CL}: \left\{\begin{array}{l} \dot{x}(t)=(A(\Omega(t))-B_2F(\Omega(t)))x(t)+B_1u_1(t) \\ y_1(t)= \underbrace{ \left[\begin{array}{c} C\\ -F \end{array}\right] }_{C_1-D_{12}F} x(t) \end{array}\right. }$

において、その ${\cal L}_2$ ゲインが $\gamma$ より小、すなわち

$\displaystyle{(38)\quad \sup_{u_1\in{\cal L}_2}\frac{||y_1(t)||_2}{||u_1(t)||_2}<\gamma }$

となるように状態フィードバックゲイン $F_1,F_2$ を求める問題を考えます。

●詳しい説明はあとで行いますが、この問題は次のように、線形行列不等式を制約条件にもつ線形計画問題として定式化されます。

Minimize $\gamma$ on $Y,Z_1,Z_2$ subject to

$\displaystyle{(39)\quad Y>0 }$

$\displaystyle{(40)\quad \left[\begin{array}{ccc} (A_1Y-B_2Z_1)^T+A_1Y-B_2Z_1 & B_1 & (C_1Y-D_{12}Z_1)^T \\ B_1^T & -\gamma I & D_{11}^T \\ C_1Y-D_{12}Z_1 & D_{11} & -\gamma I \end{array}\right]<0 }$

$\displaystyle{(41)\quad \begin{array}{l} 2\alpha Y+A_1Y-B_2Z_1+(A_1Y-B_2Z_1)^T<0 \\ \left[\begin{array}{cc} -rY & A_1Y-B_2Z_1 \\ (A_1Y-B_2Z_1)^T & -rY \end{array}\right]<0 \\ \left[\begin{array}{cc} \sin\theta(A_1Y-B_2Z_1+(A_1Y-B_2Z_1)^T) & \\ -\cos\theta(A_1Y-B_2Z_1-(A_1Y-B_2Z_1)^T) & \end{array}\right. \\ \left.\begin{array}{cc} & \cos\theta(A_1Y-B_2Z_1-(A_1Y-B_2Z_1)^T) \\ & \sin\theta(A_1Y-B_2Z_1+(A_1Y-B_2Z_1)^T) \end{array}\right]<0 \end{array} }$

$\displaystyle{(42)\quad \left[\begin{array}{ccc} (A_2Y-B_2Z_2)^T+A_2Y-B_2Z_2 & B_1 & (C_1Y-D_{12}Z_2)^T \\ B_1^T & -\gamma I & D_{11}^T \\ C_1Y-D_{12}Z_2 & D_{11} & -\gamma I \end{array}\right]<0 \\ }$

$\displaystyle{(43)\quad \begin{array}{l} 2\alpha Y+A_1Y-B_2Z_2+(A_1Y-B_2Z_2)^T<0 \\ \left[\begin{array}{cc} -rY & A_1Y-B_2Z_2 \\ (A_1Y-B_2Z_2)^T & -rY \end{array}\right]<0 \\ \left[\begin{array}{cc} \sin\theta(A_2Y-B_2Z_2+(A_2Y-B_2Z_2)^T) & \\ -\cos\theta(A_2Y-B_2Z_2-(A_2Y-B_2Z_2)^T) & \end{array}\right. \\ \left.\begin{array}{cc} & \cos\theta(A_2Y-B_2Z_2-(A_2Y-B_2Z_2)^T) \\ & \sin\theta(A_2Y-B_2Z_2+(A_2Y-B_2Z_2)^T) \end{array}\right]<0 \end{array} }$

この最小化問題の解 $Y,Z_1,Z_2$ を用いて、次のように $F_1,F_2$ を定めます。

$\displaystyle{(38)\quad F_1=Z_1Y^{-1}}$
$\displaystyle{(39)\quad F_2=Z_2Y^{-1}}$

[6] LPV制御のイメージを次に示します。LPVモデルとして表すことができる制御対象に対して、端点コントローラを求め、これを変動パラメータでスケジューリングする仕組みを表しています。

図4 LPV制御のイメージ

これに対して従来型のPID制御、LQI制御などのLTI制御は、次のように表すことができます。このようにピンポイントの制御が可能なので、一般には優れた性能を示しますが、パラメータ変動を考慮していないのでロバスト性は劣ると考えられます。

図5 LTI制御のイメージ

また、LPV制御の副産物として、スケジューリングを止めるとLTI制御方式が得られ、実装上の観点からメリットがあります。

図6 LPV制御においてスケジューリングを止めた場合

演習B03…Flipped Classroom

$1^\circ$ 　回転体のパラメータ変動を状態FBによるLPV制御で抑制したシミュレーションプログラムを以下に示す。これを実行し効果を確かめよ。

MATLAB

%cSPIN_sf_gs.m
%-----
 clear all close all
 J1=1; J2=1; J3=0.5;
 OMnom=2*pi; OMmin=0*OMnom; OMmax=2*OMnom;
 A1=[0 (J2-J3)/J1;(J3-J1)/J2 0];
 B=diag([1/J1 1/J2]); C=eye(2); D=zeros(2,2);
 S0=[zeros(2,2) B;C D];
 S1=zeros(4,4); S1(1:2,1:2)=A1;
%-----
 J1=1; J2=1; J3=0.5; 
 OMnom=2*pi; OMmin=0*OMnom; OMmax=2*OMnom;
 A1= OMmin*[0 (J2-J3)/J1;(J3-J1)/J2 0]; 
 A2= OMmax*[0 (J2-J3)/J1;(J3-J1)/J2 0];
 B=diag([1/J1 1/J2]); 
 B1=B; B2=B;
 C1=[eye(2,2);zeros(2,2)]; 
 D11=zeros(4,2); 
 D12=[zeros(2,2);eye(2,2)];
 alpha=1; r=3; th=pi/4; 
 LMIs=sf_synlmi7(A1,A2,B1,B2,C1,D11,D12,alpha,r,th);
 cobj=zeros(1,decnbr(LMIs)); 
 cobj(1)=1; 
 [cost,xopt]=mincx(LMIs,cobj); 
 Y=dec2mat(LMIs,xopt,2);  
 Z1=dec2mat(LMIs,xopt,3);
 Z2=dec2mat(LMIs,xopt,4); 
 F1=Z1/Y,pl1=eig(A1-B2*F1)
 F2=Z2/Y,pl2=eig(A2-B2*F2) 
%------
 figure(1) 
 subplot(121),dregion(alpha,0,r,th,7*[-1,1,-1,1]) 
 plot(real(pl1),imag(pl1),'*') 
 subplot(122), dregion(alpha,0,r,th,7*[-1,1,-1,1]) 
 plot(real(pl2),imag(pl2),'*')  
%------ 
 prange=OMmax-OMmin; pmax=OMmax; pmin=OMmin; 
 sim("SPIN_sf_gs") 
%-----
function LMIs=sf_synlmi7(A1,A2,B1,B2,C1,D11,D12,alpha,r,th)
 [n,m]=size(B2);  
 sth=sin(th); cth=cos(th);
 setlmis([]);
 gam=lmivar(1,[1 0]); 
 Y=lmivar(1,[n 1]); 
 Z1=lmivar(2,[m n]); 
 Z2=lmivar(2,[m n]); 
%
 lmi11=newlmi;
 lmiterm([lmi11,1,1,Y],A1,1,'s');      %#1:A*Y+Y*A' 
 lmiterm([lmi11,1,1,Z1],-B2,1,'s');    %#1:-(B2*Z+Z*B2') 
 lmiterm([lmi11,1,2,0],B1);            %#1:B1 
 lmiterm([lmi11,2,2,gam],-1,1);        %#1:-gam 
 lmiterm([lmi11,3,1,Y],C1,1);          %#1:C1*Y 
 lmiterm([lmi11,3,1,Z1],-D12,1);       %#1:D12*Z 
 lmiterm([lmi11,3,2,0],D11);           %#1:D11 
 lmiterm([lmi11,3,3,gam],-1,1);        %#1:-gam 
 lmi21=newlmi;  
 lmiterm([lmi21,1,1,Y],A1,1,'s');      %#2:A*Y+Y*A'   
 lmiterm([lmi21,1,1,Z1],-B2,1,'s');    %#2:-(B2*Z+Z'*B2')   
 lmiterm([lmi21,1,1,Y],2*alpha,1);     %#2:2*alpha*Y   
 lmi31=newlmi;  
 lmiterm([lmi31,1,1,Y],-r,1);          %#3:-r*Y  
 lmiterm([lmi31,1,2,Y],A1,1);          %#3:A*Y   
 lmiterm([lmi31,1,2,Z1],-B2,1);        %#3:-B2*Z  
 lmiterm([lmi31,2,2,Y],-r,1);          %#3:-r*Y  
 lmi41=newlmi;  
 lmiterm([lmi41,1,1,Y],sth*A1,1,'s');  %#4:sth*(A*Y+Y*A')  
 lmiterm([lmi41,1,1,Z1],-sth*B2,1,'s');%#4:-sth*(B2*Z+Z'*B2')  
 lmiterm([lmi41,1,2,Y],cth*A1,1);      %#4:cth*A*Y  
 lmiterm([lmi41,1,2,Y],1,-cth*A1');    %#4:-cth*Y*A'  
 lmiterm([lmi41,1,2,Z1],-cth*B2,1);    %#4:-cth*B2*Z  
 lmiterm([lmi41,1,2,-Z1],1,cth*B2');   %#4:cth*Z'*B2'  
 lmiterm([lmi41,2,2,Y],sth*A1,1,'s');  %#4:sth*(A*Y+Y*A')  
 lmiterm([lmi41,2,2,Z1],-sth*B2,1,'s');%#4:-sth*(B2*Z+Z'*B2')  
%
 lmi12=newlmi;
 lmiterm([lmi12,1,1,Y],A2,1,'s');      %#1:A*Y+Y*A' 
 lmiterm([lmi12,1,1,Z2],-B2,1,'s');    %#1:-(B2*Z+Z*B2') 
 lmiterm([lmi12,1,2,0],B1);            %#1:B1 
 lmiterm([lmi12,2,2,gam],-1,1);        %#1:-gam 
 lmiterm([lmi12,3,1,Y],C1,1);          %#1:C1*Y 
 lmiterm([lmi12,3,1,Z2],-D12,1);       %#1:D12*Z 
 lmiterm([lmi12,3,2,0],D11);           %#1:D11 
 lmiterm([lmi12,3,3,gam],-1,1);        %#1:-gam 
 lmi22=newlmi;  
 lmiterm([lmi22,1,1,Y],A2,1,'s');      %#2:A*Y+Y*A'   
 lmiterm([lmi22,1,1,Z2],-B2,1,'s');    %#2:-(B2*Z+Z'*B2')   
 lmiterm([lmi22,1,1,Y],2*alpha,1);     %#2:2*alpha*Y   
 lmi32=newlmi;  
 lmiterm([lmi32,1,1,Y],-r,1);          %#3:-r*Y  
 lmiterm([lmi32,1,2,Y],A2,1);          %#3:A*Y   
 lmiterm([lmi32,1,2,Z2],-B2,1);        %#3:-B2*Z  
 lmiterm([lmi32,2,2,Y],-r,1);          %#3:-r*Y  
 lmi42=newlmi;  
 lmiterm([lmi42,1,1,Y],sth*A2,1,'s');  %#4:sth*(A*Y+Y*A')  
 lmiterm([lmi42,1,1,Z2],-sth*B2,1,'s');%#4:-sth*(B2*Z+Z'*B2')  
 lmiterm([lmi42,1,2,Y],cth*A2,1);      %#4:cth*A*Y  
 lmiterm([lmi42,1,2,Y],1,-cth*A2');    %#4:-cth*Y*A'  
 lmiterm([lmi42,1,2,Z2],-cth*B2,1);    %#4:-cth*B2*Z  
 lmiterm([lmi42,1,2,-Z2],1,cth*B2');   %#4:cth*Z'*B2'  
 lmiterm([lmi42,2,2,Y],sth*A2,1,'s');  %#4:sth*(A*Y+Y*A')  
 lmiterm([lmi42,2,2,Z2],-sth*B2,1,'s');%#4:-sth*(B2*Z+Z'*B2')  
%
 lmi5=newlmi;
 lmiterm([-lmi5,1,1,Y],1,1);           %#5:Y   
 lmi6=newlmi;
 lmiterm([lmi6,1,1,gam],1,1);          %#6:gam 
 lmiterm([-lmi6,1,1,0],1e3);           %#6:1000    
 LMIs=getlmis; 
end
%-----
%eof

図7 回転体のLPV制御

$2^\circ$ 　回転体のパラメータ変動を状態FBによる $H^\infty$ 制御で抑制したシミュレーションプログラムを以下に示す。これを実行し $1^\circ$ と比較せよ。

MATLAB

%cSPIN_sf_hinf.m
%-----
 clear all close all
 J1=1; J2=1; J3=0.5;
 OMnom=2*pi; OMmin=0*OMnom; OMmax=2*OMnom;
 A1=[0 (J2-J3)/J1;(J3-J1)/J2 0];
 B=diag([1/J1 1/J2]); C=eye(2); D=zeros(2,2);
 S0=[zeros(2,2) B;C D];
 S1=zeros(4,4); S1(1:2,1:2)=A1;
%-----
 A=OMnom*[0 (J2-J3)/J1;(J3-J1)/J2 0];
 B=diag([1/J1 1/J2]); 
 B1=B; B2=B;
 C1=[eye(2,2);zeros(2,2)]; 
 D11=zeros(4,2); 
 D12=[zeros(2,2);eye(2,2)];
 alpha=1; r=3; th=pi/4; 
 LMIs=sf_synlmi6(A,B1,B2,C1,D11,D12,alpha,r,th);
 cobj=zeros(1,decnbr(LMIs)); 
 cobj(1)=1; 
 [cost,xopt]=mincx(LMIs,cobj); 
 Y=dec2mat(LMIs,xopt,2);  
 Z=dec2mat(LMIs,xopt,3);  
 F=Z/Y;  
 pl=eig(A-B2*F)
 sim("SPIN_sf_hinf")
 close all,figure(1) 
 dregion(alpha,0,r,th,5*[-1,1,-1,1]) 
 plot(real(pl),imag(pl),'*') 
%-----
function LMIs=sf_synlmi6(A,B1,B2,C1,D11,D12,alpha,r,th)
 [n,m]=size(B2);  
 sth=sin(th); cth=cos(th);
 setlmis([]);
 gam=lmivar(1,[1 0]); 
 Y=lmivar(1,[n 1]); 
 Z=lmivar(2,[m n]); 
 lmi1=newlmi;
 lmiterm([lmi1,1,1,Y],A,1,'s');      %#1:A*Y+Y*A' 
 lmiterm([lmi1,1,1,Z],-B2,1,'s');    %#1:-(B2*Z+Z*B2') 
 lmiterm([lmi1,1,2,0],B1);           %#1:B1 
 lmiterm([lmi1,2,2,gam],-1,1);       %#1:-gam 
 lmiterm([lmi1,3,1,Y],C1,1);         %#1:C1*Y 
 lmiterm([lmi1,3,1,Z],-D12,1);       %#1:D12*Z 
 lmiterm([lmi1,3,2,0],D11);          %#1:D11 
 lmiterm([lmi1,3,3,gam],-1,1);       %#1:-gam 
 lmi2=newlmi;  
 lmiterm([lmi2,1,1,Y],A,1,'s');      %#2:A*Y+Y*A'   
 lmiterm([lmi2,1,1,Z],-B2,1,'s');    %#2:-(B2*Z+Z'*B2')   
 lmiterm([lmi2,1,1,Y],2*alpha,1);    %#2:2*alpha*Y   
 lmi3=newlmi;  
 lmiterm([lmi3,1,1,Y],-r,1);         %#3:-r*Y  
 lmiterm([lmi3,1,2,Y],A,1);          %#3:A*Y   
 lmiterm([lmi3,1,2,Z],-B2,1);        %#3:-B2*Z  
 lmiterm([lmi3,2,2,Y],-r,1);         %#3:-r*Y  
 lmi4=newlmi;  
 lmiterm([lmi4,1,1,Y],sth*A,1,'s');  %#4:sth*(A*Y+Y*A')  
 lmiterm([lmi4,1,1,Z],-sth*B2,1,'s');%#4:-sth*(B2*Z+Z'*B2')  
 lmiterm([lmi4,1,2,Y],cth*A,1);      %#4:cth*A*Y  
 lmiterm([lmi4,1,2,Y],1,-cth*A');    %#4:-cth*Y*A'  
 lmiterm([lmi4,1,2,Z],-cth*B2,1);    %#4:-cth*B2*Z  
 lmiterm([lmi4,1,2,-Z],1,cth*B2');   %#4:cth*Z'*B2'  
 lmiterm([lmi4,2,2,Y],sth*A,1,'s');  %#4:sth*(A*Y+Y*A')  
 lmiterm([lmi4,2,2,Z],-sth*B2,1,'s');%#4:-sth*(B2*Z+Z'*B2')  
 lmi5=newlmi;
 lmiterm([-lmi5,1,1,Y],1,1);         %#5:Y   
 lmi6=newlmi;
 lmiterm([lmi6,1,1,gam],1,1);        %#6:gam 
 lmiterm([-lmi6,1,1,0],1e3);         %#6:1000    
 LMIs=getlmis; 
end
%-----
%eof

図8 回転体の $H_\infty$ 制御

$3^\circ$ 　 $1^\circ$ のLPV制御において、次のようにスケジューリングを止めた場合、これを実行し $1^\circ$ と $2^\circ$ と比較せよ。

図9 回転体のLPV制御においてスケジューリングを止めた場合

リャプノフの安定定理

投稿日時: 2021年10月18日投稿者: cacsd

リャプノフの安定定理

[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)から、 $E$ は $x^*=0$ の安定性を保証するためのリャプノフ関数といえます。

つぎに、(7)から

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

ですが、これは $E$ が $x^*=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=P$ 、 $Y=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

パラメータ変動システム

投稿日時: 2021年10月18日投稿者: cacsd

パラメータ変動システム…Homework

[1] 　次図のような回転体の運動を考えます。

図1　パラメータ変動システムの例

これは次の運動方程式で表されます。

$\displaystyle{(1)\quad \left\{\begin{array}{l} J_1\dot{\omega}_1(t)-\omega_2(t)\Omega(t)(J_2-J_3)=\tau_1(t) \\ J_2\dot{\omega}_2(t)-\omega_1(t)\Omega(t)(J_3-J_1)=\tau_2(t) \end{array}\right. }$

ここで、次のパラメータ変動を想定します。

$\displaystyle{(2)\quad \Omega_{min}<\Omega(t)<\Omega_{max} }$

次の状態方程式を得ます。

●いま、 $\Omega_{min}=0,\Omega_{nom}=2\pi,\Omega_{max}=4\pi$ として

$\displaystyle{(4)\quad \Omega(t)=\Omega_{nom}(1+\sin(2\pi t)) }$

のようなパラメータ変動下での零入力応答をシミュレーションしてみます。

図2　パラメータ変動システムの応答シミュレーション例

もしパラメータ変動がない場合はきれいな正弦波となりますから、パラメータ変動がある場合はかなりの動特性の変動が表れています。

演習B01…Flipped Classroom
$1^\circ$ 図2のグラフを描け。

MATLAB

%lspin.m
%-----
 clear all close all
 J1=1; J2=1; J3=0.5;
 OMnom=2*pi; OMmin=0*OMnom; OMmax=2*OMnom;
 sim("SPIN0")
%-----
 A1=[0 (J2-J3)/J1;(J3-J1)/J2 0];
 B=diag([1/J1 1/J2]); C=eye(2); D=zeros(2,2);
 S0=[zeros(2,2) B;C D];
 S1=zeros(4,4); S1(1:2,1:2)=A1;
 sim("SPIN")
%-----
%eof

図3　SPIN0.slx

Note B01　パラメータ変動システムの応答シミュレーション

このSimulinkブロック線図は

$\displaystyle{(1)\quad \left\{\begin{array}{ll} \dot{x}(t)=-(1+\sin\omega t)x(t)+u(t)\\ y(t)=x(t) \end{array}\right. }$

を次式のように書き換えて、Productを利用したものです。

$\displaystyle{(2) \left[\begin{array}{c} \dot{x}(t) \\ y(t) \end{array}\right] = ((1+\sin\omega t)\left[\begin{array}{cc} -1 & 0\\ 0 & 0 \end{array}\right] + \left[\begin{array}{cc} 0 & 1\\ 1 & 0 \end{array}\right]) \left[\begin{array}{c} x(t) \\ u(t) \end{array}\right] }$

それでは演習B01のSPIN.slxはどのように構成すればよいでしょうか？