「制御技術」カテゴリーアーカイブ

状態FB用LMI（H∞制御）

投稿日時: 投稿者:

状態FB用LMI（H∞制御）…Homework

[1] n次系\dot{x}=Ax+Bu,\ y=Cx+DuのH_\inftyノルムが\gammaより小となるためのLMI条件は、次の通りでした。

\displaystyle{(15)\quad \begin{array}{lll} && \sup_{\omega\in{\rm\bf R}}||\hat{G}(j\omega)||_2<\gamma \nonumber\\ &\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 \nonumber\\ &\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\nonumber \end{array} }

ただし

\displaystyle{(16)\quad \hat{G}(s)=C(sI-A)^{-1}B+D }


図1 ２ポートシステム

いま、次の２ポート表現されたn次系を考えます。

\displaystyle{(17)\quad P: \left\{\begin{array}{l} \dot{x}=Ax+B_1u_1+B_2u_2 \\ \underbrace{ \left[\begin{array}{c} y\\ u \end{array}\right] }_{y_1} = \underbrace{ \left[\begin{array}{c} C\\ 0 \end{array}\right] }_{C_1} x+ \underbrace{ \left[\begin{array}{c} 0\\ 0 \end{array}\right] }_{D_{11}}u_1+ \underbrace{ \left[\begin{array}{c} 0\\ I \end{array}\right] }_{D_{12}} u_2 \\ y_2=x \end{array}\right. }

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

\displaystyle{(18)\quad K: u_2=-Fy_2 }

による閉ループ系

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

において、そのH_\inftyノルムが\gammaより小、すなわち

\displaystyle{(19)\quad \sup_{\omega\in{\rm\bf R}} \underbrace{||\hat{P}_{CL}(j\omega)||_2}_{\sigma_1(\hat{P}_{CL}(j\omega))} <\gamma }

ただし

\displaystyle{(20)\quad \hat{P}_{CL}(s)= \left[\begin{array}{c} \hat{G}_y(s)\\ \hat{G}_u(s) \end{array}\right]= \left[\begin{array}{c} C(sI-A-B_2F)^{-1}B_1\\ -F(sI-A-B_2F)^{-1}B_1 \end{array}\right] }

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

この問題のLMI条件は、次のようになります。

\displaystyle{(21)\quad \exists Y>0:\ }
\displaystyle{ \left[\begin{array}{ccc} Y(A-B_2F)^T+(A-B_2F)Y & B_1 & Y(C_1-D_{12}F)^T \\ B_1^T & -\gamma I & D_{11}^T \\ (C_1-D_{12}F)Y & D_{11} & -\gamma I \end{array}\right]<0 }

ここで、変数変換Z=FYを行うと、次のようなYとZに関するLMIとなります。

\displaystyle{(22)\quad \exists Y>0:\ }
\displaystyle{ \left[\begin{array}{ccc} (AY-B_2Z)^T+AY-B_2Z & B_1 & (C_1Y-D_{12}Z)^T \\ B_1^T & -\gamma I & D_{11}^T \\ C_1Y-D_{12}Z & D_{11} & -\gamma I \end{array}\right]<0 }

これを解いてYとZを求め、次式によって状態フィードバックゲインを決定します。

\displaystyle{(23)\quad F=ZY^{-1} }

演習B33…Flipped Classroom
1^\circ 次のコードを参考にして、H_\infty制御（状態FB）を求める関数を作成せよ。

MATLAB 
%sf_syn_lmi5.m 
%----- 
 A=[0 1;-1 -2*0.01]; B1=[0;1]; B2=B1;
 C1=[1 0;0 0]; D11=[0;0]; D12=[0;1]; 
 [n,m]=size(B2);  
%----- 
 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],1,1);         %#2:Y   
%----- 
 lmi3=newlmi;
 lmiterm([lmi3,1,1,gam],1,1);        %#3:gam 
 lmiterm([-lmi3,1,1,0],1e3);         %#3:1000    
%----- 
 LMIs=getlmis; 
 cobj=zeros(1,decnbr(LMIs)); 
 cobj(1)=1; 
 [cost,xopt]=mincx(LMIs,cobj); 
 Y=dec2mat(LMIs,xopt,Y);  
 Z=dec2mat(LMIs,xopt,Z);  
 F=Z/Y;  
%----- 
 pl=eig(A-B2*F)
 close all,figure(1) 
 dregion(0,0,20,pi/4,20*[-1,1,-1,1]) 
 plot(real(pl),imag(pl),'*') 
 figure(2) 
 ol=ltisys(A,B1,C1,D11); 
 cl=ltisys(A-B2*F,B1,C1-D12*F,D11); 
 om=logspace(-2,2,200); 
 splot(ssub(cl,1,1),'sv',om),hold on 
 splot(ssub(ol,1,1),'sv',om),grid 
%----- 
%eof

2^\circ 上の２次振動系のH_\infty制御について、閉ループ系の固有値に制約を付けた場合の効果を確かめよ。

MATLAB 

%sf_syn_lmi6.m 
%----- 
 clear all, close all
 A=[0 1;-1 -2*0.01]; B1=[0;1]; B2=B1;
 C1=[1 0;0 0]; D11=[0;0]; D12=[0;1]; 
 [n,m]=size(B2);  
%----- 
 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 
%----- 
 alpha=0.1;
 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   
%----- 
 r=2;
 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  
%----- 
 theta=pi/4; sth=sin(theta); cth=cos(theta);
 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; 
 cobj=zeros(1,decnbr(LMIs)); 
 cobj(1)=1; 
 [cost,xopt]=mincx(LMIs,cobj); 
 gopt=dec2mat(LMIs,xopt,gam);  
 Y=dec2mat(LMIs,xopt,Y);  
 Z=dec2mat(LMIs,xopt,Z);  
 F=Z/Y;  
%----- 
 pl=eig(A-B2*F)
 close all,figure(1) 
 dregion(alpha,0,r,theta,r*[-1,1,-1,1]) 
 plot(real(pl),imag(pl),'*') 
 figure(2) 
 ol=ltisys(A,B1,C1,D11); 
 cl=ltisys(A-B2*F,B1,C1-D12*F,D11); 
 om=logspace(-2,2,200); 
 splot(ssub(cl,1,1),'sv',om),hold on 
 splot(ssub(ol,1,1),'sv',om),grid 
%----- 
%eof

状態FB用LMI（固有値制約）

投稿日時: 投稿者:

状態FB用LMI（固有値制約）…Homework

[0] n次系\dot{x}=Ax+Buに対する状態フィードバックu=-Fxによる閉ループ系\dot{x}=(A-BF)xにおいて、

\displaystyle{(1)\quad \lambda(A-BF)\subset{\cal D}_i\quad(i=1,2,3)}

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


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

[1] \lambda(A)\subset{\cal D}_1となるLMI条件は、次の通りでした。

\displaystyle{(2) \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} }

したがって、\lambda(A-BF)\subset{\cal D}_1となるLMI条件は、次のようになります。

\displaystyle{(3)\quad \exists Y>0:\ 2\alpha Y+(A-BF)Y+Y(A-BF)^T<0 }

これは、未知行列FとYの積FYをもつので、LMIではなく、lmisolverを用いて解くことができません。しかしながら、変数変換Z=FYを行うと、次のようなYとZに関するLMIとなります。

\displaystyle{(4)\quad \exists Y>0:\ 2\alpha Y+AY-BZ+(AY-BZ)^T<0 }

これを解いてYとZを求め、次式によって状態フィードバックゲインを決定します。

\displaystyle{(5)\quad F=ZY^{-1} }

[2] \lambda(A)\subset{\cal D}_2となるLMI条件は、次の通りでした。

\displaystyle{(6)\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} }

したがって、\lambda(A-BF)\subset{\cal D}_2となるLMI条件は、次のようになります。

\displaystyle{(7)\quad \exists Y>0:\ \left[\begin{array}{cc} -rY & (A-BF)Y \\ Y(A-BF)^T & -rY \end{array}\right]<0 }

ここで、変数変換Z=FYを行うと、次のようなYとZに関するLMIとなります。

\displaystyle{(8)\quad \exists Y>0:\ \left[\begin{array}{cc} -rY & AY-BZ \\ (AY-BZ)^T & -rY \end{array}\right]<0 }

これを解いてYとZを求め、次式によって状態フィードバックゲインを決定します。

\displaystyle{(9)\quad F=ZY^{-1} }

[3] \lambda(A)\subset{\cal D}_3となるLMI条件は、次の通りでした。

\displaystyle{(10) \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} }

したがって、\lambda(A-BF)\subset{\cal D}_3となるLMI条件は、次のようになります。

\displaystyle{(11)\quad \exists Y>0:\ }
\displaystyle{ \left[\begin{array}{cc} \sin\theta((A-BF)Y+Y(A-BF)^T) & \cos\theta((A-BF)Y-Y(A-BF)^T) \\ -\cos\theta((A-BF)Y-Y(A-BF)^T) & \sin\theta((A-BF)Y+Y(A-BF)^T) \end{array}\right] <0 }

ここで、変数変換Z=FYを行うと、次のようなYとZに関するLMIとなります。

\displaystyle{(12)\quad \exists Y>0:\ }
\displaystyle{ \left[\begin{array}{cc} \sin\theta(AY-BZ+(AY-BZ)^T) & \cos\theta(AY-BZ-(AY-BZ)^T) \\ -\cos\theta(AY-BZ-(AY-BZ)^T) & \sin\theta(AY-BZ+(AY-BZ)^T) \end{array}\right]<0 }

これを解いてYとZを求め、次式によって状態フィードバックゲインを決定します。

\displaystyle{(13)\quad F=ZY^{-1} }

演習B32…Flipped Classroom

1^\circ 次のコードを参考にして、\lambda(A-BF)\subset{\cal D}を達成する状態FBを求める関数を作成せよ。

MATLAB 
%sf_syn_lmi4.m
%----- 
 clear all, close all
 A=[0 1;0 0]; B=[0;1];
 [n,m]=size(B);  
%----- 
 setlmis([]); 
 Y=lmivar(1,[n 1]);  
 Z=lmivar(2,[m n]);  
%----- 
 alpha=1;
 lmi1=newlmi;  
 lmiterm([lmi1,1,1,Y],A,1,'s');     %#1:A*Y+Y*A'   
 lmiterm([lmi1,1,1,Z],-B,1,'s');    %#1:-(B*Z+Z'*B')   
 lmiterm([lmi1,1,1,Y],2*alpha,1);   %#1:2*alpha*Y   
%----- 
 r=5;
 lmi2=newlmi;  
 lmiterm([lmi2,1,1,Y],-r,1);        %#2:-r*Y  
 lmiterm([lmi2,1,2,Y],A,1);         %#2:A*Y   
 lmiterm([lmi2,1,2,Z],-B,1);        %#2:-B*Z  
 lmiterm([lmi2,2,2,Y],-r,1);        %#2:-r*Y  
%----- 
 theta=pi/6; sth=sin(theta); cth=cos(theta);
 lmi3=newlmi;  
 lmiterm([lmi3 1 1 Y],sth*A,1,'s'); %#3:sth*(A*Y+Y*A')  
 lmiterm([lmi3 1 1 Z],-sth*B,1,'s');%#3:-sth*(B*Z+Z'*B')  
 lmiterm([lmi3 1 2 Y],cth*A,1);     %#3:cth*A*Y  
 lmiterm([lmi3 1 2 Y],1,-cth*A');   %#3:-cth*Y*A'  
 lmiterm([lmi3 1 2 Z],-cth*B,1);    %#3:-cth*B*Z  
 lmiterm([lmi3 1 2 -Z],1,cth*B');   %#3:cth*Z'*B'  
 lmiterm([lmi3 2 2 Y],sth*A,1,'s'); %#3:sth*(A*Y+Y*A')  
 lmiterm([lmi3 2 2 Z],-sth*B,1,'s');%#3:-sth*(B*Z+Z'*B')  
%----- 
 lmi4=newlmi;
 lmiterm([-lmi4 1 1 Y],1,1);        %#4:Y  
%----- 
 LMIs=getlmis; 
 [tmin,xfeas]=feasp(LMIs);   
 Y=dec2mat(LMIs,xfeas,Y);  
 Z=dec2mat(LMIs,xfeas,Z);  
 F=Z/Y;  
%----- 
 pl=eig(A-B*F)
 close all,figure(1)  
 dregion(alpha,0,r,theta,r*[-1,1,-1,1]) 
 plot(real(pl),imag(pl),'*') 
%----- 
%eof

状態FB用LMI（変数変換）

投稿日時: 投稿者:

状態FB用LMI（変数変換）…Homework

[0]　n次系\dot{x}=Ax+Buに対する状態フィードバックu=-Fxによる閉ループ系\dot{x}=(A-BF)xを安定化するように状態フィードバックゲインFを求める手順を復習しておきます。

いま、A_F=A-BFのジョルダン分解V\Lambda V^{-1}を仮定して

\displaystyle{(1)\quad A-BF=\underbrace{V\Lambda V^{-1}}_{A_F} \Rightarrow AV-V\Lambda=B\underbrace{FV}_{G} }

ここで、\Lambdaが対角行列の場合、A_Fの固有値\lambda_1',\cdots,\lambda_n'に対応する固有ベクトルをv_1,\cdots,v_nとして

\displaystyle{(2)\quad A \underbrace{ \left[\begin{array}{ccc} v_1&\cdots&v_n \end{array}\right] }_{V} - \underbrace{ \left[\begin{array}{ccc} v_1&\cdots&v_n \end{array}\right] }_{V} \underbrace{ {\rm diag}\{\lambda_1',\cdots,\lambda_m'\} }_{\Lambda}}

\displaystyle{= \left[\begin{array}{ccc} (A-\lambda_1'I_n)v_1&\cdots&(A-\lambda_n'I_n)v_n \end{array}\right] =B \underbrace{ \left[\begin{array}{ccc} g_1&\cdots&g_n \end{array}\right] }_{G} }

これより、もし\lambda_1',\cdots,\lambda_n'がAの固有値と一致しないならば、固有ベクトルの表現式として次式を得ます。

\displaystyle{(3)\quad \underbrace{ \left[\begin{array}{ccc} v_1&\cdots&v_n \end{array}\right] }_{V}= \left[\begin{array}{ccc} (A-\lambda_1'I_n)^{-1}Bg_1&\cdots&(A-\lambda_n'I_n)^{-1}Bg_n \end{array}\right] }

したがって、状態フィードバックゲイン行列Fを求める次の手順が考えられます。

1) A_Fの固有値\lambda_1',\cdots,\lambda_n'をAの固有値とは異なるように指定する。

2) A_Fの固有ベクトルv_1,\cdots,v_nを決める自由度g_1,\cdots,g_nを、(3)のVが正則になるように与える。

3) 次式で Fを求める。

\displaystyle{(4)\quad \boxed{F=GV^{-1}} }

以上の方法は、複素左半平面{\rm\bf LHP}内にある望ましい固有値

\displaystyle{(5)\quad \{\lambda_1',\cdots,\lambda_n'\}\subset{\rm\bf LHP} }

をピンポイントで指定しています。これに対して、以下では、複素左半平面内に望ましい領域を指定して、状態フィードバックを設計する方法について述べます。

[1] n次系\dot{x}=Ax+Buに対する状態フィードバックu=-Fxによる閉ループ系\dot{x}=(A-BF)xにおいて、

\displaystyle{(6)\quad \lambda(A-BF)\subset{\rm\bf LHP}}

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

\lambda(A)\subset{\rm\bf LHP}となるLMI条件（Y-LMIによる）は、次の通りでした。

\displaystyle{(7) \begin{array}{lll} &&\lambda(A)\subset{\rm\bf LHP}=\{s=x+jy\in{\rm\bf C}: s+s^*<0 \}\\ &\Leftrightarrow& \exists Y>0:\ AY+YA^T<0 \end{array} }

したがって、\lambda(A-BF)\subset{\rm\bf LHP}となるLMI条件は、次のようになります。

\displaystyle{(8)\quad \exists Y>0:\ (A-BF)Y+Y(A-BF)^T<0 }

これは、未知行列FとYの積FYをもつので、LMIではなく、lmisolverを用いて解くことができません。しかしながら、変数変換Z=FYを行うと、次のようなYとZに関するLMIとなります。

\displaystyle{(9)\quad \exists Y>0:\ AY-BZ+(AY-BZ)^T<0 }

これを解いてYとZを求め、次式によって状態フィードバックゲインを決定します。

\displaystyle{(10)\quad F=ZY^{-1} }

●ちなみに、次のLMI条件（X-LMIによる）

\displaystyle{(11) \begin{array}{lll} &&\lambda(A)\subset{\rm\bf LHP}=\{s=x+jy\in{\rm\bf C}: s+s^*<0 \}\\ &\Leftrightarrow& \exists X>0:\ XA+A^TX<0 \end{array} }

を用いて、次式を得ます。

\displaystyle{(12)\quad \exists X>0:\ 2\alpha X+X(A-BF)+(A-BF)^TX<0\\ }

ここで、未知行列FとXの積はXBFの形で表され、変数変換の手法を適用することができないことに注意してください。

演習B31…Flipped Classroom
1^\circ 次のコードを参考にして、安定化状態FBを求める関数を作成せよ。

MATLAB 

%sf_syn_lmi0.m 
%----- 
 A=[0 1;0 0]; B=[0;1];
 [n,m]=size(B);  
%----- 
 setlmis([]); 
 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],-B,1,'s');    %#1:-(B*Z+Z'*B')   
%----- 
 lmi2=newlmi; 
 lmiterm([-lmi2 1 1 Y],1,1);        %#2:Y  
%----- 
 LMIs=getlmis;
 [tmin,xfeas]=feasp(LMIs);  
 Y=dec2mat(LMIs,xfeas,Y);  
 Z=dec2mat(LMIs,xfeas,Z);  
 F=Z/Y;  
%----- 
 pl=eig(A-B*F)
 close all,figure(1)  
 dregion(0,0,5,pi/2,5*[-1,1,-1,1])  
 plot(real(pl),imag(pl),'*')  
%----- 
%eof

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

投稿日時: 投稿者:

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

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


図1　直達項をもたないｎ次系

\displaystyle{ \begin{array}{l} \displaystyle{{\rm tr} (CW_cC^T)<\gamma^2}\\\\ \displaystyle{\Leftrightarrow \exists X_c>0,\ Q_c>0:\\ (1a)\quad\left[\begin{array}{cc} A^TX_c+X_cA & X_cB \\ B^TX_c & -I \end{array}\right]<0,\ \left[\begin{array}{cc} X_c & C^T \\ C & Q_c \end{array}\right]>0,\ {\rm tr}(Q_c)<\gamma^2}\\\\ \displaystyle{\Leftrightarrow \exists Y_c>0,\ Q_c>0:\\ (1b)\quad \left[\begin{array}{cc} Y_cA^T+AY_c & Y_cB \\ B^TY_c & -I \end{array}\right]<0,\ \left[\begin{array}{cc} Y_c & Y_cC^T \\ CY_c & Q_c \end{array}\right]>0,\ {\rm tr}(Q_c)<\gamma^2} \end{array}}
\displaystyle{ \begin{array}{l} \displaystyle{{\rm tr} (B^TW_oB)<\gamma^2}\\\\ \displaystyle{\Leftrightarrow \exists X_o>0,\ Q_o>0:\\ (2a)\quad\left[\begin{array}{cc} AX_o+X_oA^T & X_oC^T \\ CX_o & -I \end{array}\right]<0,\ \left[\begin{array}{cc} X_o & B \\ B^T & Q_o \end{array}\right]>0,\ {\rm tr}(Q_o)<\gamma^2}\\\\ \displaystyle{\Leftrightarrow \exists Y_o>0,\ Q_o>0:\\ (2b)\quad\left[\begin{array}{cc} Y_oA+A^TY_o & Y_oC^T \\ CY_o & -I \end{array}\right]<0,\ \left[\begin{array}{cc} Y_o & Y_oB \\ B^TY_o & Q_o \end{array}\right]>0,\ {\rm tr}(Q_o)<\gamma^2} \end{array}}

ここで、W_cとW_oはそれぞれ

(3)\quad \begin{array}{l} \displaystyle{W_c=\int_0^\infty \exp(At)BB^T\exp(A^Tt)dt}\\ \displaystyle{W_o=\int_0^\infty \exp(A^Tt)C^TC\exp(At)dt} \end{array} }

で与えられ、次式が成り立つことに注意します。

(4)\quad \begin{array}{l} {\rm tr} (CW_cC^T)={\rm tr} (B^TW_oB) \end{array} }

さて、W_cとW_oはそれぞれリャプノフ方程式

\displaystyle{(5)\quad \begin{array}{l} W_cA^T+AW_c+BB^T=0\\ W_oA+A^TW_o+C^TC=0 \end{array} }

の解であることから

\displaystyle{(6)\quad \begin{array}{l} V_cA^T+AV_c+BB^T<0\\ V_oA+A^TV_o+C^TC<0 \end{array} }

を満たすV_cとV_oに対して、X_c=V_c^{-1}、X_o=V_o^{-1}とおくと

\displaystyle{(7)\quad \begin{array}{l} A^TX_c+X_cA+X_cBB^TX_c<0 \Leftrightarrow\ \left[\begin{array}{cc} A^TX_c+X_cA & X_cB \\ B^TX_c & -I \end{array}\right]<0 \\ AX_o+X_oA^T+X_oC^TCX_o<0 \Leftrightarrow\ \left[\begin{array}{cc} AX_o+X_oA^T & X_oC^T \\ CX_o & -I \end{array}\right]<0 \end{array} }

が成り立ちます。一方、(6)と(5)から

\displaystyle{(8)\quad \begin{array}{l} V_cA^T+AV_c+BB^T=(V_c-W_c)A^T+A(V_c-W_c)<0\\ V_oA+A^TV_o+C^T=(V_o-W_o)A+A^T(V_o-W_o)<0 \end{array} }

となって、次を得ます。

\displaystyle{(9)\quad \begin{array}{l} V_c-W_c=X_c^{-1}-W_c>0\\ V_o-W_o=X_o^{-1}-W_o>0 \end{array} }

以上の準備の下で、次を示します。

\displaystyle{(10)\quad \begin{array}{l} {\rm tr} (CW_cC^T)<\gamma^2 \Leftrightarrow\ \left[\begin{array}{cc} X_c & C^T \\ C & Q_c \end{array}\right]>0,\ {\rm tr} (Q_c)<\gamma^2\\ {\rm tr} (B^TW_oB)<\gamma^2 \Leftrightarrow\ \left[\begin{array}{cc} X_o & B \\ B^T & Q_o \end{array}\right]>0,\ {\rm tr} (Q_o)<\gamma^2 \end{array} }

まず必要性については、(9)から、十分小さな\epsilon>0を選んで

\displaystyle{(11)\quad \begin{array}{l} CX_c^{-1}C^T-(CW_cC^T+\epsilon I)>0 \Leftrightarrow\ \left[\begin{array}{cc} X_c & C^T \\ C & CW_cC^T+\epsilon I \end{array}\right]>0\\ B^TX_o^{-1}B-(B^TW_oB+\epsilon I)>0 \Leftrightarrow\ \left[\begin{array}{cc} X_o & B \\ B^T & B^TW_oB+\epsilon I \end{array}\right]>0 \end{array} }

および

(12)\quad \begin{array}{l} {\rm tr} (CW_cC^T+\epsilon I)<\gamma^2\\ {\rm tr} (B^TW_oB+\epsilon I)<\gamma^2 \end{array} }

を満足させることができます。ここで、Q_c=CW_cC^T+\epsilon I、Q_o=B^TW_oB+\epsilon Iとおきます。次に十分性については、Q_c=CW_cC^T、Q_o=B^TW_oBの場合を考えれば自明です。

●以上では、次のシュール補元に関する公式を多用しました。

\displaystyle{\boxed{ \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} }} \displaystyle{\boxed{ \begin{array}{lll} && \left[\begin{array}{cc} P & M \\ M^T & Q \end{array}\right]>0\\ &\Leftrightarrow& MQ^{-1}M^T-P>0,\ Q>0\\ &\Leftrightarrow& P>0,\ M^TP^{-1}M-Q>0 \end{array} }}

演習B23…Flipped Classroom

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

MATLAB 
%ana_lmi6.m
%------
 clear all, close all
 T=1; K=1;
 A=-1/T; B=K/T; 
 C=1; D=0; 
 n=1; m=1; p=1; 
 CT=1;
 gam=1 %input('gamma =  '); 
%------
 setlmis([]);
 X=lmivar(1,[n 1]);
 Q=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],-1);      %#1:-I
%------
 lmi2=newlmi;  
 lmiterm([-lmi2 1 1 X],1,1);    %#2:X 
 lmiterm([-lmi2 1 2 0],C');     %#2:C' 
 lmiterm([-lmi2 2 2 Q],1,1);    %#2:Q
%------ 
 lmi3=newlmi;   
 lmiterm([-lmi3 1 1 X],1,1);    %#3:X>0
%------ 
 lmi4=newlmi;   
 lmiterm([-lmi4 1 1 Q],1,1);    %#4:Q>0 
%------ 
 lmi5=newlmi;   
 lmiterm([lmi5 1 1 Q],CT,1);    %#5:tr(Q)=CT*Q
 lmiterm([lmi5 1 1 0],-gam^2);  %#5:-gam^2
%------ 
 LMIs=getlmis; 
 [tmin,xfeas]=feasp(LMIs);
 X=dec2mat(LMIs,xfeas,X) 
 Q=dec2mat(LMIs,xfeas,Q) 
%------
%eof

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

投稿日時: 投稿者:

性能解析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 ２次系の周波数応答の比較

伝達特性

投稿日時: 投稿者:

伝達特性…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

投稿日時: 投稿者:

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

投稿日時: 投稿者:

セクタ制約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

投稿日時: 投稿者:

円制約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

投稿日時: 投稿者:

α制約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} }}