制御系CAD

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

混合感度問題

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

以下では時間関数 $x(t)$ のラプラス変換を $\displaystyle{\hat{x}(s)=\int_0^\infty x(t)e^{-st}\,dt}$ のように表します。

混合感度問題…Homework

[1] モデルの不確かさがある場合の制御系設計のために次を考えます。

図1 摂動前と摂動後の閉ループ系

上図の閉ループ系について、次の伝達関数（インパルス応答のラプラス変換、または入出力のラブラス変換の比）を定義します。

$\displaystyle{(1)\quad\hat{L}(s)=\hat{G}(s)\hat{K}(s)}$

$\displaystyle{(2)\quad\hat{R}(s)=1+\hat{L}(s)}$

$\displaystyle{(3)\quad\boxed{\hat{S}(s)=\frac{1}{\hat{R}(s)}=\frac{1}{1+\hat{L}(s)}}}$

$\displaystyle{(4)\quad\boxed{\hat{T}(s)=\frac{\hat{L}(s)}{\hat{R}(s)}=\frac{\hat{L}(s)}{1+\hat{L}(s)}}}$

ここで、 $\hat{L}(s)$ , $\hat{R}(s)$ , $\hat{S}(s)$ , $\hat{T}(s)$ はそれぞれループ伝達関数、環送差、感度関数、相補感度関数とよばれ、次の制約式が成り立ちます。

$\displaystyle{(5)\quad\boxed{\hat{S}(s)+\hat{T}(s)=1}}$

制御対象の伝達特性の不確かさを表すために次式を考えます。

$\displaystyle{(6)\quad\hat{G}_{\hat{\Delta}}(s)=(1+\hat{\Delta}(s))\hat{G}(s)}$

ここで、 $\hat{\Delta}(s)$ は乗法的摂動と呼ばれます。このとき摂動後のループ伝達関数、環送差、感度関数、相補感度関数を次のように表します。

$\displaystyle{(7)\quad\hat{L}_{{\Delta}}(s)=\hat{G}_{{\Delta}}(s)\hat{K}(s)=(1+\hat{\Delta}(s))\hat{L}(s)}$

$\displaystyle{(8)\quad\hat{R}_{\Delta}}(s)=1+\hat{L}_{{\Delta}}(s)}$

$\displaystyle{(9)\quad\hat{S}_{{\Delta}}(s)=\frac{1}{\hat{R}_{{\Delta}}(s)}}$

$\displaystyle{(10)\quad\hat{T}_{{\Delta}}(s)=\frac{\hat{L}_{{\Delta}}(s)}{\hat{R}_{{\Delta}}(s)}}$

●以下では、ラプラスの $s$ を虚軸上に限定して（ $s=j\omega$ ）議論を進めます。

まず、乗法的摂動ついては、あとで具体例を示しますが、ここでは次を仮定します。

$\displaystyle{(11)\quad\boxed{|\hat{\Delta}(j\omega)|<|W_T(j\omega)|}}$

一般に、乗法的摂動のノルムバンド $|W_T(j\omega)|$ は高周波域で大きくなります。

また、(5)より、次の制約があり、 $\hat{S}(j\omega)$ と $\hat{T}(j\omega)$ は同時には小さくできません。

$\displaystyle{(12)\quad\hat{S}(j\omega)+\hat{T}(j\omega)=1}$

そこで、低周波域で顕著な重み関数 $|\hat{W}_S|$ と高周波域で顕著な重み関数 $|\hat{W}_T|$ を用いて

$\displaystyle{(13)\quad \boxed{|\hat{W}_S(j\omega)\hat{S}(j\omega)|^2+|\hat{W}_T(j\omega)\hat{T}(j\omega)|^2}} }$

の最小化を図ります。この問題は感度関数 $\hat{S}$ と相補感度関数 $\hat{T}$ を同時に考慮しているので、混合感度問題と呼ばれています。

●これは次のようなループ伝達関数の整形をもたらします。すなわち、低周波域で支配的な $|\hat{W}_S(j\omega)|$ を用いてループゲイン $|\hat{L}(j\omega)|$ を大きくします。また、高周波域で支配的な $|\hat{W}_T(j\omega)|$ を用いてループゲイン $|\hat{L}(j\omega)|$ を小さくすることを意味します。また、安定余裕を確保するために中間周波域ではループゲイン $|\hat{L}(j\omega)|$ をなだらかに落とすことが求められます。

図2 ループ整形

[2] 混合感度問題を解くことによりモデルの不確かさを考慮したロバスト制御系を設計することができるのですが、その妥当性をナイキスト軌跡を用いて考えます。

ノミナル安定性　これはナイキストの安定判別法を表しています。

図3 NS条件

ノミナル性能　これはクリティカルポイント $-1+j0$ を半径 $|\hat{W}_S|$ の円領域に拡大し、ナイキスト軌跡がクリティカルポインから十分離れるようにしたものです。

図4 NP条件

ロバスト安定性　これはモデルの不確かさがある場合のナイキスト軌跡に幅をもたせ場合の安定判別法を表しています。

図5 RS条件

ロバスト性能　これはこれはモデルの不確かさがある場合のナイキスト軌跡に幅をもたせ場合に十分性能を持たせることを意味しています。

図6 RP条件

以上、４つの仕様のうち、ロバスト性能が最も強力で、これが満足されれば、ロバスト安定性とノミナル性能は満足されます。ロバスト性能仕様は絶対値和となり、解析的に扱いにくいので、十分条件を考慮して、混合感度問題の定式化が行われているといえます。

図7 ロバスト性能の十分条件

[3] モデルの不確かさの一例として無駄時間の場合を考えます。ここで、無駄時間 $L$ とは

$\displaystyle{(14)\quad y(t)=u(t-L)}$

のような入出力関係を言います。これをラプラス変換すると

$\displaystyle{(15)\quad \hat{y}(s)=\int_0^\infty u(t-L)e^{-st}\,dt=e^{-sL}\hat{u}(s) }$

したがって、伝達関数は

$\displaystyle{(16)\quad \frac{\hat{u}(s)}{\hat{\delta}(s)}=e^{-Ls} }$

となります。これを有理関数近似またはパデー近似(Pade approximation)することがよく行われます。

$\displaystyle{(17)\quad \begin{array}{l} 1+(-Ls)+\frac{1}{2}(-Ls)^2+\cdots\simeq\frac{1+\beta s}{1+\alpha s}\\ \Rightarrow\ 1+\alpha s+(-Ls)(1+\alpha s)+\frac{1}{2}(-Ls)^2+\cdots\simeq 1+\beta s\\ \Rightarrow\ \alpha-L=\beta,\ \frac{1}{2}L^2-\alpha L=0\\ \Rightarrow\ \alpha=\frac{1}{2}L,\beta=-\frac{1}{2}L\\ \Rightarrow\ \boxed{e^{-Ls}\simeq\frac{1-\frac{1}{2}Ls}{1+\frac{1}{2}Ls}} \end{array} }$

これは１次近似ですが、２次近似の場合は同様にして

$\displaystyle{(18)\quad \boxed{e^{-Ls}\simeq\frac{1-\frac{1}{2}Ls-\frac{1}{12}L^2s^2}{1+\frac{1}{2}Ls+\frac{1}{12}L^2s^2}}}$

となります。ステップ応答を比較してみると次のようになります。

図8 無駄時間Pade 近似の比較

それでは入力端に無駄時間がある場合の乗法的摂動を求めてみます。

$\displaystyle{(19)\quad \underbrace{\hat{G}(s)e^{-Ls}}_{\hat{G}_{\hat{\Delta}(s)}}=(1+\underbrace{e^{-Ls}-1}_{\hat{\Delta}(s)})\hat{G}(s) }$

このノルムバンドは次のように見積ることができます。

$\displaystyle{(20)\quad \begin{array}{lll} \underbrace{|e^{-Ls}-1|}_{|\hat{\Delta}(s)|} &=&|(-Ls)+\frac{1}{2}(-Ls)^2+\cdots|\\ &=&|(-Ls)||1+\frac{1}{2}(-Ls)+\cdots| \simeq\underbrace{(1+\epsilon)\left|\frac{s}{L^{-1}}\right|}_{|\hat{W}_T(s)|} \end{array} }$

次は $L=5$ の場合の例です。

図9 無駄時間のノルムバンド

演習B61…Flipped Classroom
$1^\circ$ 船舶の回頭モデル

$\displaystyle{(21)\quad \left[\begin{array}{c} \dot{\psi}(t) \\ \dot{r}(t) \end{array}\right] = \left[\begin{array}{cc} 0 & 1\\ 0 & -\frac{1}{T} \end{array}\right] \left[\begin{array}{c} \psi(t) \\ r(t) \end{array}\right] + \left[\begin{array}{c} 0 \\ \frac{K}{T} \end{array}\right] \delta(t-L) }$

に対して、適当なPIDコントローラ

$\displaystyle{(22)\quad \delta(t)=K_P(\psi_c-\psi(t))-K_Dr(t)+\int_0^t(\psi_c-\psi(\tau))d\tau }$

が与えられているとする。次のプログラムを用いて、モデル誤差のノルムバンド（図9相当）、開ループ整形（図2相当）、閉ループ整形（図7参照）を描き、ロバスト性について検討せよ。

図10 開ループ整形と閉ループ整形

MATLAB

%b61.m 
%-----
 clear all, close all
 T1=118; T2=7.8; T3=18.5; T=T1+T2-T3, K=0.185
 zeta=0.5/sqrt(T*K), wn=sqrt(K/T)
 A=[0 1;0 -2*zeta*wn]; B=[0;wn^2]; C=[1 0];
%-----   
 w=logspace(-3,1,300);
 G(:,:)=freqresp(ss(A,B,C,[]),w); G=G.'; 
 ga_G=20*log10(abs(G));
 ph_G=180/pi*atan2(imag(G),real(G));
 lag=10;
 H=exp(-lag*i*w);
 ga_H=20*log10(abs(H));
 ph_H=180/pi*atan2(imag(H),real(H));
 GH=G.*H;
 ga_GH=20*log10(abs(GH));
 ph_GH=180/pi*atan2(imag(GH),real(GH));
 figure(1)
 subplot(211),semilogx(w,ga_G,w,ga_H,'b',w,ga_GH,'r'),grid,legend('G','H','GH')
 subplot(212),semilogx(w,ph_G,w,ph_H,'b',w,ph_GH,'r'),grid,legend('G','H','GH')
 figure(2)
 plot(real(G),imag(G),'o',real(GH),imag(GH),'+',-1,0,'or'),
 axis([-2 2 -2 2]); grid
%-----  
 ga_DL=20*log10(abs(H-1));
 ga_WT=20*log10(abs(lag*i*w));
 figure(3)
 semilogx(w,ga_DL,w,ga_WT'),grid,legend('DL','WT') 
%-----   
 Kp=2.3145921; Kd=46.886428; Ki=0.0097080;
 K=Kp+Kd*(i*w)+Ki./(i*w);
 ga_K=20*log10(abs(K));
 ph_K=180/pi*atan2(imag(K),real(K));
 L=GH.*K;
 ga_L=20*log10(abs(L));
 ph_L=180/pi*atan2(imag(L),real(L));
 WS=Ki*ones(1,length(w))./(i*w);
 WT=lag*i*w;
 ga_WS=20*log10(abs(WS));
 ga_WT=20*log10(abs(WT));
 figure(4)
 %semilogx(w,ga_GH,w,ga_K,w,ga_L,w,ga_WS,w,-ga_WT),grid,legend('GH','K','L','WS','WTi') 
 semilogx(w,ga_L,w,ga_WS,w,-ga_WT),grid,legend('L','WS','WTi') 
%-----   
 S=ones(1,length(w))./(1+L);
 ga_S=20*log10(abs(S));
 ph_S=180/pi*atan2(imag(S),real(S));
 T=L./(1+L);
 ga_T=20*log10(abs(T));
 ph_T=180/pi*atan2(imag(T),real(T));
 WS=Ki*ones(1,length(w))./(i*w);
 WT=lag*i*w;
 ga_WS=20*log10(abs(WS));
 ga_WT=20*log10(abs(WT));
 mu=20*log10(abs(WS.*S)+abs(WT.*T));
 figure(5)
 %semilogx(w,ga_S,w,ga_T,w,-ga_WS,w,-ga_WT,w,mu),grid,legend('S','T','WSi','WTi','mu') 
 semilogx(w,ga_S,w,ga_T,w,-ga_WS,w,-ga_WT),grid,legend('S','T','WSi','WTi') 
%-----  
%eof

AIP

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

AIP…Homework

[1]　次図のようなアーム駆動型倒立振子の実験装置を考えます。これは某企業のエンジニアによって考案された実験装置です。以下、AIPと参照します。

図1 アーム駆動型倒立振子

外見は２重倒立振子のように見えますが、第１振子はモータ軸に直結され、第２振子が第１振子の先端に軸支されています。興味深いのはこの制御対象AIPは、次図のように複数個の平衡状態をもつので、できるだけ多くの平衡状態まわり安定化できるコントローラをどのように設計するかを考えてみます。

図2　複数個の平衡状態

まず、モータには速度駆動型のドライバが装備されているので、この制御対象AIPのモデリングを次図のように行います。

図3 速度制御型アクチュエータを持つ場合のモデリング

これから、制御対象AIPのLPVモデルを次図のように得ます。

図4 AIPのLPVモデルの導出

次図にAIPのLPV制御の概念図を示します。

図5 LPV制御の概念図

次図にAIPのLPV制御系のシミュレータを示します。

図6 LPV制御のためのSimulink図

次図にAIPの制御系の応答比較を示します。

図7 制御方式の比較

演習B54…Flipped Classroom

$1^\circ$ AIPのLPV制御系設計を行う、次のコードを説明せよ

MATLAB

%cADIP_of_gs.m
%-----
 clear all, close all
 global ell1 ell2 th10 th20
 ell1=0.13; ell2=0.15; m1=0.05; m2=0.03; g=9.8;
 Ta=0.1; Ka=1;
 th10=0; th20=0;
 th1=0; th2=70/180*pi;
 rymax=2*ell1*cos(th1); rymin=2*ell1*cos(th2); rynom=rymax;
 rxmin=2*ell1*sin(th1); rxmax=2*ell1*sin(th2); rxnom=rxmax; 
 A0=[0 1 0 0;
     3*g/(4*ell2) 0 -3*g/(4*ell2) 0;
     0 0 0 0;
     0 0 0 -1/Ta];
 A1=[0 0 0 0;
     0 0 0 0;
     0 0 0 1;
     0 0 0 0];
 Amax=A0+rymax*A1;
 Amin=A0+rymin*A1;
 Anom=A0+rynom*A1; 
 B=[0;0;0;Ka/Ta];
 CM=[1 0 0 0;
     0 0 1 0];
 DM=[0;0];
 S0=[A0 B;CM DM];
 S1=[A1 zeros(4,1);zeros(2,5)]; 
%-----
 wD=0.02; wI=0.5; C0=CM(2,:);
 AA1=[Amin zeros(4,1);-C0 0];
 AA2=[Amax zeros(4,1);-C0 0]; 
 B1=[zeros(4,1);1]; B2=[B;0];
 C1=[zeros(1,4) wI;wD*C0*Anom 0]; D11=zeros(2,1); D12=[0;wD*C0*B];
 C2=[CM zeros(2,1);zeros(1,4) 1]; D21=[0;0;0]; D22=[0;0;0];
%===== 
 alpha=0.1; r=50; th=90/180*pi;
 LMIs=of_synlmi7(AA1,AA2,B1,B2,C1,C2,D11,D12,D21,D22,alpha,r,th);
 cobj=zeros(1,decnbr(LMIs)); 
 cobj(1)=1; 
 [cost,xopt]=mincx(LMIs,cobj);  
 gopt=dec2mat(LMIs,xopt,1)
 R=dec2mat(LMIs,xopt,2); 
 S=dec2mat(LMIs,xopt,3);
%=====
 ak1=dec2mat(LMIs,xopt,4); 
 bk1=dec2mat(LMIs,xopt,5); 
 ck1=dec2mat(LMIs,xopt,6); 
 dk1=dec2mat(LMIs,xopt,7); 
 [u,sd,v]=svd(eye(size(AA1,1))-R*S); 
 sd=diag(sqrt(1./diag(sd))); 
 Ni=sd*v'; Mti=u*sd; 
 AK1=Ni*(ak1-S*(AA1-B2*dk1*C2)*R-bk1*C2*R-S*B2*ck1)*Mti; 
 BK1=Ni*(bk1-S*B2*dk1); 
 CK1=(ck1-dk1*C2*R)*Mti; 
 DK1=dk1; 
 K1=[AK1 BK1;CK1 DK1];
%-----  
 pl1=eig(AA1)
 ACL1=[AA1+B2*DK1*C2 B2*CK1;
      BK1*C2 AK1];
 plCL1=eig(ACL1)
 figure(1)
 dregion(-alpha,0,r,th,r*[-1,1,-1,1]) 
 plot(real(pl1),imag(pl1),'x',real(plCL1),imag(plCL1),'*')
%-----  
 ACL1=[AA1+B2*DK1*C2 B2*CK1;BK1*C2 AK1];
 BCL1=[B1+B2*DK1*D21;BK1*D21];
 CCL1=[C2(1,:) zeros(1,5)];
 AK1=[AK1 BK1(:,3);zeros(1,6)];
 BK1=[BK1(:,1:2) zeros(5,1); -1 0 1];
 CK1=[CK1 DK1(:,3)];
 DK1=[DK1(:,1:2) 0];
 SK1=[AK1 BK1;CK1 DK1];
%-----  
 w=logspace(-2,2,100); 
 ga_WS=20*log10(abs(wI./(i*w))); 
 ga_WT=20*log10(abs(wD*(i*w)));
 G(:,:)=freqresp(ss(Amax,B,C0,[]),w); G=G.'; ga_G=20*log10(abs(G));
 K(:,:)=freqresp(ss(AK1,BK1(:,1),CK1,DK1(:,1)),w); K=K.'; ga_K=20*log10(abs(K)); 
 T(:,:)=freqresp(ss(ACL1,BCL1(:,1),CCL1(1,:),[]),w); T=T.'; ga_T=20*log10(abs(T)); 
 ga_S=20*log10(abs(1-T)); 
 figure(2)
 subplot(121),semilogx(w,ga_G+ga_K,w,ga_WS,'b',w,-ga_WT,'r'),grid,legend('L','WS','1/WT')
 subplot(122),semilogx(w,ga_T,'r',w,-ga_WT,'r',w,ga_S,'b',w,-ga_WS,'b'),grid,legend('T','1/WT','S','1/WS')
%=====
 ak2=dec2mat(LMIs,xopt,8); 
 bk2=dec2mat(LMIs,xopt,9); 
 ck2=dec2mat(LMIs,xopt,10); 
 dk2=dec2mat(LMIs,xopt,11); 
 [u,sd,v]=svd(eye(size(AA2,1))-R*S); 
 sd=diag(sqrt(1./diag(sd))); 
 Ni=sd*v'; Mti=u*sd; 
 AK2=Ni*(ak2-S*(AA2-B2*dk2*C2)*R-bk2*C2*R-S*B2*ck2)*Mti; 
 BK2=Ni*(bk2-S*B2*dk2); 
 CK2=(ck2-dk2*C2*R)*Mti; 
 DK2=dk2; 
 K2=[AK2 BK2;CK2 DK2];
%-----  
 pl2=eig(AA2)
 ACL2=[AA2+B2*DK2*C2 B2*CK2;
       BK2*C2 AK2];
 plCL2=eig(ACL2)
 figure(3)
 dregion(-alpha,0,r,th,r*[-1,1,-1,1]) 
 plot(real(pl2),imag(pl2),'x',real(plCL2),imag(plCL2),'*') 
%-----  
 ACL2=[AA2+B2*DK2*C2 B2*CK2;BK2*C2 AK2];
 BCL2=[B1+B2*DK2*D21;BK2*D21];
 CCL2=[C2(1,:) zeros(1,5)];
 AK2=[AK2 BK2(:,3);zeros(1,6)];
 BK2=[BK2(:,1:2) zeros(5,1); -1 0 1];
 CK2=[CK2 DK2(:,3)];
 DK2=[DK2(:,1:2) 0];
 SK2=[AK2 BK2;CK2 DK2];
%-----  
 w=logspace(-2,2,100); 
 ga_WS=20*log10(abs(wI./(i*w))); 
 ga_WT=20*log10(abs(wD*(i*w)));
 G(:,:)=freqresp(ss(Amin,B,C0,[]),w); G=G.'; ga_G=20*log10(abs(G));
 K(:,:)=freqresp(ss(AK2,BK2(:,1),CK2,DK2(:,1)),w); K=K.'; ga_K=20*log10(abs(K)); 
 T(:,:)=freqresp(ss(ACL2,BCL2(:,1),CCL2(1,:),[]),w); T=T.'; ga_T=20*log10(abs(T)); 
 ga_S=20*log10(abs(1-T)); 
 figure(4)
 subplot(121),semilogx(w,ga_G+ga_K,w,ga_WS,'b',w,-ga_WT,'r'),grid,legend('L','WS','1/WT')
 subplot(122),semilogx(w,ga_T,'r',w,-ga_WT,'r',w,ga_S,'b',w,-ga_WS,'b'),grid,legend('T','1/WT','S','1/WS')
%-----
 prange=rymax-rymin; pmax=rymax; pmin=rymin; 
%sim("ADIP_of_gs")
 sim("ADIP_of_gs_2015a")
%-----
function LMIs=of_synlmi7(A1,A2,B1,B2,C1,C2,D11,D12,D21,D22,alpha,r,th)
 [n,m]=size(B2);  [p,n]=size(C2);
 setlmis([]); 
 gam=lmivar(1,[1 0]);  
 [R,xxx,Rdec]=lmivar(1,[n 1]); 
 [S,xxx,Sdec]=lmivar(1,[n 1]); 
 Ak1=lmivar(2,[n n]); 
 Bk1=lmivar(2,[n p]); 
 Ck1=lmivar(2,[m n]); 
 Dk1=lmivar(2,[m p]); 
 Ak2=lmivar(2,[n n]); 
 Bk2=lmivar(2,[n p]); 
 Ck2=lmivar(2,[m n]); 
 Dk2=lmivar(2,[m p]);  
%=====
 lmiRS1=newlmi;
 lmiterm([lmiRS1 1 1 R],A1,1,'s');       %#1:R*A'+AR 
 lmiterm([lmiRS1 1 1 Ck1],B2,1,'s');     %#1:Ck'*B2'+B2*Ck 
 lmiterm([lmiRS1 2 1 0],A1');            %#1:A' 
 lmiterm([lmiRS1 2 1 Ak1],1,1);          %#1:Ak 
 lmiterm([lmiRS1 2 1 -Dk1],C2',B2');     %#1:C2'*Dk'*B2' 
 lmiterm([lmiRS1 2 2 S],1,A1,'s');       %#1:A'*S+S*A 
 lmiterm([lmiRS1 2 2 Bk1],1,C2,'s');     %#1:C2'*Bk'+Bk*C2 
 lmiterm([lmiRS1 1 3 0],B1);             %#1:B1 
 lmiterm([lmiRS1 1 3 Dk1],B2,D21);       %#1:B2*Dk*D21 
 lmiterm([lmiRS1 2 3 S],1,B1);           %#1:S*B1 
 lmiterm([lmiRS1 2 3 Bk1],1,D21);        %#1:Bk*D21 
 lmiterm([lmiRS1 3 3 gam],-1,1);         %#1:-gam 
 lmiterm([lmiRS1 4 1 R],C1,1);           %#1:C1*R 
 lmiterm([lmiRS1 4 1 Ck1],D12,1);        %#1:D12*Ck 
 lmiterm([lmiRS1 4 2 0],C1);             %#1:C1 
 lmiterm([lmiRS1 4 2 Dk1],D12,C2);       %#1:D12*Dk*C2 
 lmiterm([lmiRS1 4 3 0],D11);            %#1:D11 
 lmiterm([lmiRS1 4 3 Dk1],D12,D21);      %#1:D12*Dk*D21 
 lmiterm([lmiRS1 4 4 gam],-1,1);         %#1:-gam 
%-----
 lmiPL11=newlmi; 
 lmiterm([lmiPL11 1 1 R],A1,1,'s');      %#2:R*A'+A*R 
 lmiterm([lmiPL11 1 1 Ck1],B2,1,'s');    %#2:Ck'*B2'+B2*Ck 
 lmiterm([lmiPL11 2 1 Ak1],1,1);         %#2:Ak 
 lmiterm([lmiPL11 1 2 0],A1);            %#2:A 
 lmiterm([lmiPL11 1 2 Dk1],B2,C2);       %#2:B2*Dk*C2 
 lmiterm([lmiPL11 2 2 S],1,A1,'s');      %#2:A'*S+S*A 
 lmiterm([lmiPL11 2 2 Bk1],1,C2,'s');    %#2:C2'*Bk'+Bk*C2 
% 
 lmiterm([lmiPL11 1 1 R],2*alpha,1);     %#2:2*alpha*R
 lmiterm([lmiPL11 2 1 0],2*alpha);       %#2:2*alpha*I 
 lmiterm([lmiPL11 2 2 S],2*alpha,1);     %#2:2*alpha*S 
%-----
 lmiPL21=newlmi; 
 lmiterm([lmiPL21 1 1 R],-r,1);          %#3:-r*R 
 lmiterm([lmiPL21 2 1 0],-r);            %#3:-r*I 
 lmiterm([lmiPL21 2 2 S],-r,1);          %#3:-r*S 
% 
 lmiterm([lmiPL21 1 3 R],A1,1);          %#3:A*R
 lmiterm([lmiPL21 1 3 Ck1],B2,1);        %#3:B2*Ck 
 lmiterm([lmiPL21 2 3 Ak1],1,1);         %#3:Ak 
 lmiterm([lmiPL21 1 4 0],A1);            %#3:A 
 lmiterm([lmiPL21 1 4 Dk1],B2,C2);       %#3:B2*Dk*C2 
 lmiterm([lmiPL21 2 4 S],1,A1);          %#3:S*A 
 lmiterm([lmiPL21 2 4 Bk1],1,C2);        %#3:Bk*C2 
% 
 lmiterm([lmiPL21 3 3 R],-r,1);          %#3:-r*R
 lmiterm([lmiPL21 4 3 0],-r);            %#3:-r*I 
 lmiterm([lmiPL21 4 4 S],-r,1);          %#3:-r*S 
%-----
 sth=sin(th); cth=cos(th);
 lmiPL31=newlmi; 
 lmiterm([lmiPL31 1 1 R],sth*A1,1,'s');  %#4:sth*(R*A'+A*R) 
 lmiterm([lmiPL31 1 1 Ck1],sth*B2,1,'s');%#4:sth*(Ck'*B2'+B2*Ck) 
 lmiterm([lmiPL31 2 1 Ak1],sth,1);       %#4:sth*(Ak) 
 lmiterm([lmiPL31 1 2 0],sth*A1);        %#4:sth*(A) 
 lmiterm([lmiPL31 1 2 Dk1],sth*B2,C2);   %#4:sth*(B2*Dk*C2) 
 lmiterm([lmiPL31 2 2 S],1,sth*A1,'s');  %#4:sth*(A'*S+S*A) 
 lmiterm([lmiPL31 2 2 Bk1],1,sth*C2,'s');%#4:sth*(C2'*Bk'+Bk*C2) 
% 
 lmiterm([lmiPL31 1 3 R],cth*A1,1);      %#1:cth*(A*R) 
 lmiterm([lmiPL31 1 3 R],1,-cth*A1');    %#1:cth*(-R*A')
 lmiterm([lmiPL31 1 3 Ck1],cth*B2,1);    %#1:cth*(B*Ck) 
 lmiterm([lmiPL31 1 3 -Ck1],-cth*B2',1); %#1:cth*(-Ck'*B') 
 lmiterm([lmiPL31 2 3 Ak1],cth,1);       %#4:cth*(Ak) 
 lmiterm([lmiPL31 1 4 -Ak1],-cth,1);     %#4:cth*(-Ak') 
 lmiterm([lmiPL31 1 4 0],A1);            %#4:cth*(A) 
 lmiterm([lmiPL31 2 3 0],-A1');          %#4:cth*(-A') 
 lmiterm([lmiPL31 1 4 Dk1],cth*B2,C2);   %#4:cth*(B2*Dk*C2) 
 lmiterm([lmiPL31 2 3 -Dk1],-cth*C2',B2');%#4:cth*(-C2'*Dk'*B2') 
 lmiterm([lmiPL31 2 4 S],1,cth*A1);      %#4:cth*(S*A) 
 lmiterm([lmiPL31 2 4 S],-cth*A1',1);    %#4:cth*(-A'*S) 
 lmiterm([lmiPL31 2 4 Bk1],1,cth*C2);    %#4:cth*(Bk*C2) 
 lmiterm([lmiPL31 2 4 -Bk1],-cth*C2',1); %#4:cth*(-C2'*Bk') 
% 
 lmiterm([lmiPL31 3 3 R],sth*A1,1,'s');  %#4:sth*(R*A'+A*R)
 lmiterm([lmiPL31 3 3 Ck1],sth*B2,1,'s');%#4:sth*(Ck'*B2'+B2*Ck) 
 lmiterm([lmiPL31 4 3 Ak1],sth,1);       %#4:sth*(Ak) 
 lmiterm([lmiPL31 3 4 0],sth*A1);        %#4:sth*(A) 
 lmiterm([lmiPL31 3 4 Dk1],sth*B2,C2);   %#4:sth*(B2*Dk*C2) 
 lmiterm([lmiPL31 4 4 S],1,sth*A1,'s');  %#4:sth*(A'*S+S*A) 
 lmiterm([lmiPL31 4 4 Bk1],1,sth*C2,'s');%#4:sth*(C2'*Bk'+Bk*C2) 
%=====
 lmiRS2=newlmi;
 lmiterm([lmiRS2 1 1 R],A2,1,'s');       %#1:R*A'+AR 
 lmiterm([lmiRS2 1 1 Ck2],B2,1,'s');     %#1:Ck'*B2'+B2*Ck 
 lmiterm([lmiRS2 2 1 0],A2');            %#1:A' 
 lmiterm([lmiRS2 2 1 Ak2],1,1);          %#1:Ak 
 lmiterm([lmiRS2 2 1 -Dk2],C2',B2');     %#1:C2'*Dk'*B2' 
 lmiterm([lmiRS2 2 2 S],1,A2,'s');       %#1:A'*S+S*A 
 lmiterm([lmiRS2 2 2 Bk2],1,C2,'s');     %#1:C2'*Bk'+Bk*C2 
 lmiterm([lmiRS2 1 3 0],B1);             %#1:B1 
 lmiterm([lmiRS2 1 3 Dk2],B2,D21);       %#1:B2*Dk*D21 
 lmiterm([lmiRS2 2 3 S],1,B1);           %#1:S*B1 
 lmiterm([lmiRS2 2 3 Bk2],1,D21);        %#1:Bk*D21 
 lmiterm([lmiRS2 3 3 gam],-1,1);         %#1:-gam 
 lmiterm([lmiRS2 4 1 R],C1,1);           %#1:C1*R 
 lmiterm([lmiRS2 4 1 Ck2],D12,1);        %#1:D12*Ck 
 lmiterm([lmiRS2 4 2 0],C1);             %#1:C1 
 lmiterm([lmiRS2 4 2 Dk2],D12,C2);       %#1:D12*Dk*C2 
 lmiterm([lmiRS2 4 3 0],D11);            %#1:D11 
 lmiterm([lmiRS2 4 3 Dk2],D12,D21);      %#1:D12*Dk*D21 
 lmiterm([lmiRS2 4 4 gam],-1,1);         %#1:-gam 
%-----
 lmiPL12=newlmi; 
 lmiterm([lmiPL12 1 1 R],A2,1,'s');      %#2:R*A'+A*R 
 lmiterm([lmiPL12 1 1 Ck2],B2,1,'s');    %#2:Ck'*B2'+B2*Ck 
 lmiterm([lmiPL12 2 1 Ak2],1,1);         %#2:Ak 
 lmiterm([lmiPL12 1 2 0],A2);            %#2:A 
 lmiterm([lmiPL12 1 2 Dk2],B2,C2);       %#2:B2*Dk*C2 
 lmiterm([lmiPL12 2 2 S],1,A2,'s');      %#2:A'*S+S*A 
 lmiterm([lmiPL12 2 2 Bk2],1,C2,'s');    %#2:C2'*Bk'+Bk*C2 
% 
 lmiterm([lmiPL12 1 1 R],2*alpha,1);     %#2:2*alpha*R
 lmiterm([lmiPL12 2 1 0],2*alpha);       %#2:2*alpha*I 
 lmiterm([lmiPL12 2 2 S],2*alpha,1);     %#2:2*alpha*S 
%-----
 lmiPL22=newlmi; 
 lmiterm([lmiPL22 1 1 R],-r,1);          %#3:-r*R 
 lmiterm([lmiPL22 2 1 0],-r);            %#3:-r*I 
 lmiterm([lmiPL22 2 2 S],-r,1);          %#3:-r*S 
% 
 lmiterm([lmiPL22 1 3 R],A2,1);          %#3:A*R
 lmiterm([lmiPL22 1 3 Ck2],B2,1);        %#3:B2*Ck 
 lmiterm([lmiPL22 2 3 Ak2],1,1);         %#3:Ak 
 lmiterm([lmiPL22 1 4 0],A2);            %#3:A 
 lmiterm([lmiPL22 1 4 Dk2],B2,C2);       %#3:B2*Dk*C2 
 lmiterm([lmiPL22 2 4 S],1,A2);          %#3:S*A 
 lmiterm([lmiPL22 2 4 Bk2],1,C2);        %#3:Bk*C2 
% 
 lmiterm([lmiPL22 3 3 R],-r,1);          %#3:-r*R
 lmiterm([lmiPL22 4 3 0],-r);            %#3:-r*I 
 lmiterm([lmiPL22 4 4 S],-r,1);          %#3:-r*S 
%-----
 sth=sin(th); cth=cos(th);
 lmiPL32=newlmi; 
 lmiterm([lmiPL32 1 1 R],sth*A2,1,'s');  %#4:sth*(R*A'+A*R) 
 lmiterm([lmiPL32 1 1 Ck2],sth*B2,1,'s');%#4:sth*(Ck'*B2'+B2*Ck) 
 lmiterm([lmiPL32 2 1 Ak2],sth,1);       %#4:sth*(Ak) 
 lmiterm([lmiPL32 1 2 0],sth*A2);        %#4:sth*(A) 
 lmiterm([lmiPL32 1 2 Dk2],sth*B2,C2);   %#4:sth*(B2*Dk*C2) 
 lmiterm([lmiPL32 2 2 S],1,sth*A2,'s');  %#4:sth*(A'*S+S*A) 
 lmiterm([lmiPL32 2 2 Bk2],1,sth*C2,'s');%#4:sth*(C2'*Bk'+Bk*C2) 
% 
 lmiterm([lmiPL32 1 3 R],cth*A2,1);      %#1:cth*(A*R) 
 lmiterm([lmiPL32 1 3 R],1,-cth*A2');    %#1:cth*(-R*A')
 lmiterm([lmiPL32 1 3 Ck2],cth*B2,1);    %#1:cth*(B*Ck) 
 lmiterm([lmiPL32 1 3 -Ck2],-cth*B2',1); %#1:cth*(-Ck'*B') 
 lmiterm([lmiPL32 2 3 Ak2],cth,1);       %#4:cth*(Ak) 
 lmiterm([lmiPL32 1 4 -Ak2],-cth,1);     %#4:cth*(-Ak') 
 lmiterm([lmiPL32 1 4 0],A2);            %#4:cth*(A) 
 lmiterm([lmiPL32 2 3 0],-A2');          %#4:cth*(-A') 
 lmiterm([lmiPL32 1 4 Dk2],cth*B2,C2);   %#4:cth*(B2*Dk*C2) 
 lmiterm([lmiPL32 2 3 -Dk2],-cth*C2',B2');%#4:cth*(-C2'*Dk'*B2') 
 lmiterm([lmiPL32 2 4 S],1,cth*A2);      %#4:cth*(S*A) 
 lmiterm([lmiPL32 2 4 S],-cth*A2',1);    %#4:cth*(-A'*S) 
 lmiterm([lmiPL32 2 4 Bk2],1,cth*C2);    %#4:cth*(Bk*C2) 
 lmiterm([lmiPL32 2 4 -Bk2],-cth*C2',1); %#4:cth*(-C2'*Bk') 
% 
 lmiterm([lmiPL32 3 3 R],sth*A2,1,'s');  %#4:sth*(R*A'+A*R)
 lmiterm([lmiPL32 3 3 Ck2],sth*B2,1,'s');%#4:sth*(Ck'*B2'+B2*Ck) 
 lmiterm([lmiPL32 4 3 Ak2],sth,1);       %#4:sth*(Ak) 
 lmiterm([lmiPL32 3 4 0],sth*A2);        %#4:sth*(A) 
 lmiterm([lmiPL32 3 4 Dk2],sth*B2,C2);   %#4:sth*(B2*Dk*C2) 
 lmiterm([lmiPL32 4 4 S],1,sth*A2,'s');  %#4:sth*(A'*S+S*A) 
 lmiterm([lmiPL32 4 4 Bk2],1,sth*C2,'s');%#4:sth*(C2'*Bk'+Bk*C2)  
%=====
 posX=-newlmi;
 lmiterm([posX 1 1 R],1,1);              %#5:R 
 lmiterm([posX 2 1 0],1);                %#5:I 
 lmiterm([posX 2 2 S],1,1);              %#5:S 
%-----
 lmiDk1=-newlmi;                        
 lmiterm([lmiDk1 1 1 0],1e2);            %#6:1e2 
 lmiterm([lmiDk1 2 2 0],1e2);            %#6:1e2 
 lmiterm([lmiDk1 2 1 Dk1],1,1);          %#6:Dk 
%-----
 lmiDk2=-newlmi;                        
 lmiterm([lmiDk2 1 1 0],1e2);            %#6:1e2 
 lmiterm([lmiDk2 2 2 0],1e2);            %#6:1e2 
 lmiterm([lmiDk2 2 1 Dk2],1,1);          %#6:Dk 
%-----
 lmig=newlmi;                          
 lmiterm([lmig,1,1,gam],1,1);            %#7:gam 
 lmiterm([lmig,1,1,0],-1e3);             %#7:1e3 
 LMIs=getlmis;
end
%-----
%eof

ACC BP

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

次の論文で提起されているロバスト制御ベンチマーク問題ACC-BPを考えます。
Bong WieDennis S. Bernsteint: Benchmark Problems for Robust Control Design

ACC BP…Homework
[1]　次図のような連結台車を考えます。

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

$\displaystyle{(1)\quad \begin{array}{l} m_1\ddot{x}_1(t)=k(x_2(t)-x_1(t))+u(t)+w_1(t)\\ m_2\ddot{x}_2(t)=k(x_1(t)-x_2(t))+w_2(t) \end{array} }$

ここで、パラメータの値は次が想定されています。

$\displaystyle{(2)\quad \begin{array}{l} m_1=m_2=1\\ 0.5\le k \le 2 \end{array} }$

●バネ定数 $k$ を乱数を100個発生させて所定の範囲内で変化させるとき、台車２に対してインパルス外乱が加わるときの $x_2-x_1$ の振舞いのシミュレーションを行ってみます。

MATLAB

%ACCBP0.m
%-----
 clear all, close all
 m1=1; m2=1; kmin=0.5; kmax=2;
 k=(kmax-kmin)/2; 
 A0=[0 0 1 0;0 0 0 1;0 0 0 0;0 0 0 0];
 A1=[0 0 0 0;0 0 0 0;-1/m1 1/m1 0 0;1/m2 -1/m2 0 0];  
 B=[0 0 0;0 0 0;1/m1 0 1/m1;0 1/m2 0];
 B1=B(:,1); B2=B(:,2); 
 C=[-1 1 0 0];
 D=0;
 S0=[A0 B1;C D];
 S1=zeros(5,5); S1(1:4,1:4)=A1;
%-----
 k=(kmax-kmin)/2;
 t0=0; t1=20; t=t0:(t1-t0)/100:t1;
 for i=1:100
   k=rand*(kmax-kmin)+kmin;
   A=A0+k*A1;
   sys=ss(A,B1,C,D);
   [y,t]=initial(sys,B2,t);
   figure(1),plot(t,y),hold on
 end
 axis([t0 t1 -2 2]),grid
%-----
%eof

図1　台車２にインパルス外乱が加わるときの $x_2-x_1$ の振舞い

原論文にはさまざまな問題設定がありますが、ここでは、台車２に対するインパルス外乱の下で、台車１の整定時間を20秒以内とする制御系を設計したいとします。ただし、操作入力は台車１に与えられ、次の振幅制限

$\displaystyle{(3)\quad \begin{array}{l} |u(t)|\le 1 \end{array} }$

が課されるとします。さらに観測出力は台車１の変位 $x_1(t)$ とします。

[2] 制御対象のLPVモデルを導出します。

$\displaystyle{(4)\quad \begin{array}{l} \underbrace{ \left[\begin{array}{cc} m_1 & 0 \\ 0 & m_2 \end{array}\right] }_{M} \underbrace{ \left[\begin{array}{c} \ddot{x}_1(t)\\ \ddot{x}_2(t) \end{array}\right] }_{\ddot{x}(t)} +\underbrace{ \left[\begin{array}{cc} k & -k \\ -k & -k \end{array}\right] }_{K(k)} \underbrace{ \left[\begin{array}{ccc} {x}_1(t)\\ {x}_2(t) \end{array}\right] }_{{\xi}(t)} =\underbrace{ \left[\begin{array}{ccc} 1\\ 0 \end{array}\right] }_{E_{2\times1}}u(t) +\underbrace{ \left[\begin{array}{ccc} w_{1}\\ w_{2} \end{array}\right] }_{w_{d}} \end{array} }$

$\displaystyle{(5)\quad \begin{array}{l} \left[\begin{array}{c} \dot{\xi}(t)\\ \ddot{\xi}(t) \end{array}\right] = \underbrace{ \left[\begin{array}{cccc} 0_{2\times2} & I_2\\ -M^{-1}K & -M^{-1}D \end{array}\right] }_{A(k)} \left[\begin{array}{c} {\xi}(t)\\ \dot{\xi}(t) \end{array}\right]\\ + \left[\begin{array}{c} 0_{2\times1}\\ M^{-1}E_{2\times1} \end{array}\right]u(t) + \left[\begin{array}{c} 0_{2\times2}\\ M^{-1} \end{array}\right]w_d \end{array} }$

ここで、 $\alpha_1\le\alpha\le\alpha_2$ に対する次の内分式に注目します。

$\displaystyle{(6)\quad \alpha=\underbrace{\frac{\alpha_2-\alpha}{\alpha_2-\alpha_1}}_{p_{1}(\alpha)}\alpha_1+ \underbrace{\frac{\alpha-\alpha_1}{\alpha_2-\alpha_1}}_{p_{2}(\alpha)}\alpha_2 }$

$\displaystyle{(7)\quad p_1(\alpha)+p_2(\alpha)=1 }$

$k_{min}\le k \le k_{max}$ のとき、上の状態方程式は、端点モデル

$\displaystyle{(8)\quad \begin{array}{l} \dot{x}(t)=\underbrace{A(k_{min})}_{A_1}x(t)+Bu(t)\\ \dot{x}(t)=\underbrace{A(k_{max})}_{A_2}x(t)+Bu(t) \end{array} }$

を、 $p_1(k),p_2(k)$ によって重み付けして、LPVモデル

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

として表されます。実際

$\displaystyle{(9)\quad k= \underbrace{\frac{k_{max}-k}{k_{max}-k_{min}}}_{p_1(k)}k_{min} + \underbrace{\frac{k-k_{min}}{k_{max}-k_{min}}}_{p_2(k)}k_{max} }$

$\displaystyle{(10)\quad A(k)= \underbrace{\frac{k_{max}-k}{k_{max}-k_{min}}}_{p_1(k)} \underbrace{A(k_{min})}_{A_1} + \underbrace{\frac{k-k_{min}}{k_{max}-k_{min}}}_{p_2(k)} \underbrace{A(k_{max})}_{A_2} }$

[3]　状態FBによる解として求める場合の２ポートシステムの例を次式に示します。

$\displaystyle{(11) P: \left\{\begin{array}{l} \left[\begin{array}{c} \dot{\xi}\\ \ddot{\xi} \end{array}\right] = \underbrace{ \left[\begin{array}{cccc} 0_{2\times2} & I_2\\ -M^{-1}K & -M^{-1}D \end{array}\right] }_{A(k)=p_1(k)A_1+p_2(k)A_2} \left[\begin{array}{c} {\xi}\\ \dot{\xi} \end{array}\right] + \underbrace{ \left[\begin{array}{c} 0_{2\times2}\\ M^{-1} \end{array}\right] }_{B_1} w_d + \underbrace{ \left[\begin{array}{c} 0_{2\times1}\\ M^{-1}E_{2\times1} \end{array}\right] }_{B_2} u\\ \underbrace{ \left[\begin{array}{c} x_2-x_1 \\ x_1\\ 2u \end{array}\right] }_{y_1} = \underbrace{ \left[\begin{array}{cc|cc} -1 & 1 & 0 & 0\\ 1 & 0 & 0 & 0\\ 0 & 0 & 0 & 0\\ \end{array}\right] }_{C_1} \left[\begin{array}{c} {\xi}\\ \dot{\xi} \end{array}\right] + \underbrace{ 0_{3\times2} }_{D_{11}} w_d + \underbrace{ \left[\begin{array}{c} 0 \\ 0 \\ 2 \end{array}\right] }_{D_{12}} u \end{array}\right. }$

バネ定数 $k$ を乱数を100個発生させて所定の範囲内で変化させるとき、台車２に対してインパルス外乱が加わるときの台車１の変位 $x_1$ と操作入力 $u$ の振舞いのシミュレーションを次に示します。

図2　状態FBによる求解結果

[4]　出力FBによる解として求める場合の２ポートシステムの例を次式に示します。

$\displaystyle{(12) P: \left\{\begin{array}{l} \left[\begin{array}{c} \dot{\xi}\\ \ddot{\xi} \end{array}\right] = \underbrace{ \left[\begin{array}{cccc} 0_{2\times2} & I_2\\ -M^{-1}K & -M^{-1}D \end{array}\right] }_{A(k)=p_1(k)A_1+p_2(k)A_2} \left[\begin{array}{c} {\xi}\\ \dot{\xi} \end{array}\right] + \underbrace{ \left[\begin{array}{c} 0_{2\times2}\\ M^{-1} \end{array}\right] }_{B_1} w_d + \underbrace{ \left[\begin{array}{c} 0_{2\times1}\\ M^{-1}E_{2\times1} \end{array}\right] }_{B_2} u\\ \underbrace{ \left[\begin{array}{c} x_2-x_1 \\ x_1\\ 2u \end{array}\right] }_{y_1} = \underbrace{ \left[\begin{array}{cc|cc} -1 & 1 & 0 & 0\\ 1 & 0 & 0 & 0\\ 0 & 0 & 0 & 0\\ \end{array}\right] }_{C_1} \left[\begin{array}{c} {\xi}\\ \dot{\xi} \end{array}\right] + \underbrace{ 0_{3\times2} }_{D_{11}} w_d + \underbrace{ \left[\begin{array}{c} 0 \\ 0 \\ 2 \end{array}\right] }_{D_{12}} u\\ \underbrace{ x_1 }_{y_2} = \underbrace{ \left[\begin{array}{cc|cc} 1 & 0 & 0 & 0 \end{array}\right] }_{C_2} \left[\begin{array}{c} {\xi}(t)\\ \dot{\xi}(t) \end{array}\right] + \underbrace{ 0_{1\times2} }_{D_{21}} w_d + \underbrace{ 0 }_{D_{22}} u \end{array}\right. }$

図3　出力FBによる求解結果

演習B53…Flipped Classroom

$1^\circ$ 図２のグラフを描く、次のコードを説明せよ

MATLAB

%cACCBP_sf_gs.m
%-----
 clear all, close all
 m1=1; m2=1; kmin=0.5; kmax=2;
 A0=[0 0 1 0;0 0 0 1;0 0 0 0;0 0 0 0];
 A1=[0 0 0 0;0 0 0 0;-1/m1 1/m1 0 0;1/m2 -1/m2 0 0];  
 C=[-1 1 0 0];
 D=0;
 a1= A0+kmin*A1; 
 a2= A0+kmax*A1;
 A1=a1; A2=a2;
 BG=[0 0 0;0 0 0;1/m1 0 1/m1;0 1/m2 0]; 
 B1=BG(:,2); %負荷外乱 
 B2=BG(:,3); %操作入力
 C1=[1 -1 0 0;1 0 0 0;0 0 0 0];  
 D11=zeros(3,1); 
 D12=[0;0;2];
%-----
 alpha=0.05; r=2; 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); 
 gopt=dec2mat(LMIs,xopt,1)
 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) 
%-----
 A0=[0 0 1 0;0 0 0 1;0 0 0 0;0 0 0 0];
 A1=[0 0 0 0;0 0 0 0;-1/m1 1/m1 0 0;1/m2 -1/m2 0 0];  
 prange=kmax-kmin; pmax=kmax; pmin=kmin; 
 t0=0; t1=20; t=t0:(t1-t0)/100:t1;
 for i=1:100
   k=rand*(kmax-kmin)+kmin; 
   A=A0+k*A1;
   p1=(pmax-k)/prange; p2=(k-pmin)/prange;
   F=p1*F1+p2*F2;
   sys=ss(A-B2*F,B1,C,D);
   [y,t]=initial(sys,B1,t);
   figure(1),plot(t,y),hold on
 end
 axis([t0 t1 -2 2]),grid
%-----
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

$2^\circ$ 図３のグラフを描く、次のコードを説明せよ。

MATLAB

%cACCBP_of_gs.m
%-----
 clear all, close all
 m1=1; m2=1; kmin=0.5; kmax=2; 
 A1=[0 0 1 0;0 0 0 1;-kmin/m1 kmin/m1 0 0;kmin/m2 -kmin/m2 0 0];
 A2=[0 0 1 0;0 0 0 1;-kmax/m1 kmax/m1 0 0;kmax/m2 -kmax/m2 0 0];
 BG=[0 0 0;0 0 0;1/m1 0 1/m1;0 1/m2 0]; 
 B1=BG(:,2); B2=BG(:,3);
 C1=[-1 1 0 0;1 0 0 0;0 0 0 0]; D11=zeros(3,1); D12=[0;0;2];
 C2=[1 0 0 0]; D21=0; D22=0;
%-----
 alpha=0.05; r=5; th=0.5*pi;
 LMIs=of_synlmi7(A1,A2,B1,B2,C1,C2,D11,D12,D21,D22,alpha,r,th);
 cobj=zeros(1,decnbr(LMIs)); 
 cobj(1)=1; 
 [cost,xopt]=mincx(LMIs,cobj);  
 gopt=dec2mat(LMIs,xopt,1)
 R=dec2mat(LMIs,xopt,2); 
 S=dec2mat(LMIs,xopt,3); 
 %-----
 ak1=dec2mat(LMIs,xopt,4); 
 bk1=dec2mat(LMIs,xopt,5); 
 ck1=dec2mat(LMIs,xopt,6); 
 dk1=dec2mat(LMIs,xopt,7); 
 [u,sd,v]=svd(eye(size(A1,1))-R*S); 
 sd=diag(sqrt(1./diag(sd))); 
 Ni=sd*v'; Mti=u*sd; 
 AK1=Ni*(ak1-S*(A1-B2*dk1*C2)*R-bk1*C2*R-S*B2*ck1)*Mti; 
 BK1=Ni*(bk1-S*B2*dk1); 
 CK1=(ck1-dk1*C2*R)*Mti; 
 DK1=dk1; 
 K1=[AK1 BK1;CK1 DK1];
 %-----
 ak2=dec2mat(LMIs,xopt,8); 
 bk2=dec2mat(LMIs,xopt,9); 
 ck2=dec2mat(LMIs,xopt,10); 
 dk2=dec2mat(LMIs,xopt,11); 
 [u,sd,v]=svd(eye(size(A2,1))-R*S); 
 sd=diag(sqrt(1./diag(sd))); 
 Ni=sd*v'; Mti=u*sd; 
 AK2=Ni*(ak2-S*(A2-B2*dk2*C2)*R-bk2*C2*R-S*B2*ck2)*Mti; 
 BK2=Ni*(bk2-S*B2*dk2); 
 CK2=(ck2-dk2*C2*R)*Mti; 
 DK2=dk2; 
 K2=[AK2 BK2;CK2 DK2];
%-----
 A0=[0 0 1 0;0 0 0 1;0 0 0 0;0 0 0 0];
 A1=[0 0 0 0;0 0 0 0;-1/m1 1/m1 0 0;1/m2 -1/m2 0 0];  
 prange=kmax-kmin; pmax=kmax; pmin=kmin; 
 t0=0; t1=20; t=t0:(t1-t0)/100:t1;
 for i=1:100
   k=rand*(kmax-kmin)+kmin; 
   A=A0+k*A1;
   p1=(pmax-k)/prange; p2=(k-pmin)/prange;
   AK=p1*AK1+p2*AK2;
   BK=p1*BK1+p2*BK2;
   CK=p1*CK1+p2*CK2;
   DK=p1*DK1+p2*DK2;
   ACL=[A+B2*DK*C2 B2*CK;BK*C2 AK];   
   BCL=[B1+B2*DK*D21;BK*D21];
   CCL=[1 0 0]*[C1+D12*DK*C2 D12*CK];
   DCL=[1 0 0]*(D11+D12*DK*D21);
   sys=ss(ACL,BCL,CCL,DCL);
   [y,t]=initial(sys,[B1;zeros(4,1)],t);
   figure(1),plot(t,y),hold on
 end
 axis([t0 t1 -2 2]),grid
%-----
function LMIs=of_synlmi7(A1,A2,B1,B2,C1,C2,D11,D12,D21,D22,alpha,r,th)
 [n,m]=size(B2);  [p,n]=size(C2);
 setlmis([]); 
 gam=lmivar(1,[1 0]);  
 [R,xxx,Rdec]=lmivar(1,[n 1]); 
 [S,xxx,Sdec]=lmivar(1,[n 1]); 
 Ak1=lmivar(2,[n n]); 
 Bk1=lmivar(2,[n p]); 
 Ck1=lmivar(2,[m n]); 
 Dk1=lmivar(2,[m p]); 
 Ak2=lmivar(2,[n n]); 
 Bk2=lmivar(2,[n p]); 
 Ck2=lmivar(2,[m n]); 
 Dk2=lmivar(2,[m p]);  
%=====
 lmiRS1=newlmi;
 lmiterm([lmiRS1 1 1 R],A1,1,'s');       %#1:R*A'+AR 
 lmiterm([lmiRS1 1 1 Ck1],B2,1,'s');     %#1:Ck'*B2'+B2*Ck 
 lmiterm([lmiRS1 2 1 0],A1');            %#1:A' 
 lmiterm([lmiRS1 2 1 Ak1],1,1);          %#1:Ak 
 lmiterm([lmiRS1 2 1 -Dk1],C2',B2');     %#1:C2'*Dk'*B2' 
 lmiterm([lmiRS1 2 2 S],1,A1,'s');       %#1:A'*S+S*A 
 lmiterm([lmiRS1 2 2 Bk1],1,C2,'s');     %#1:C2'*Bk'+Bk*C2 
 lmiterm([lmiRS1 1 3 0],B1);            %#1:B1 
 lmiterm([lmiRS1 1 3 Dk1],B2,D21);       %#1:B2*Dk*D21 
 lmiterm([lmiRS1 2 3 S],1,B1);          %#1:S*B1 
 lmiterm([lmiRS1 2 3 Bk1],1,D21);        %#1:Bk*D21 
 lmiterm([lmiRS1 3 3 gam],-1,1);        %#1:-gam 
 lmiterm([lmiRS1 4 1 R],C1,1);          %#1:C1*R 
 lmiterm([lmiRS1 4 1 Ck1],D12,1);        %#1:D12*Ck 
 lmiterm([lmiRS1 4 2 0],C1);            %#1:C1 
 lmiterm([lmiRS1 4 2 Dk1],D12,C2);       %#1:D12*Dk*C2 
 lmiterm([lmiRS1 4 3 0],D11);           %#1:D11 
 lmiterm([lmiRS1 4 3 Dk1],D12,D21);      %#1:D12*Dk*D21 
 lmiterm([lmiRS1 4 4 gam],-1,1);        %#1:-gam 
%-----
 lmiPL11=newlmi; 
 lmiterm([lmiPL11 1 1 R],A1,1,'s');      %#2:R*A'+A*R 
 lmiterm([lmiPL11 1 1 Ck1],B2,1,'s');    %#2:Ck'*B2'+B2*Ck 
 lmiterm([lmiPL11 2 1 Ak1],1,1);         %#2:Ak 
 lmiterm([lmiPL11 1 2 0],A1);            %#2:A 
 lmiterm([lmiPL11 1 2 Dk1],B2,C2);       %#2:B2*Dk*C2 
 lmiterm([lmiPL11 2 2 S],1,A1,'s');      %#2:A'*S+S*A 
 lmiterm([lmiPL11 2 2 Bk1],1,C2,'s');    %#2:C2'*Bk'+Bk*C2 
% 
 lmiterm([lmiPL11 1 1 R],2*alpha,1);    %#2:2*alpha*R
 lmiterm([lmiPL11 2 1 0],2*alpha);      %#2:2*alpha*I 
 lmiterm([lmiPL11 2 2 S],2*alpha,1);    %#2:2*alpha*S 
%-----
 lmiPL21=newlmi; 
 lmiterm([lmiPL21 1 1 R],-r,1);         %#3:-r*R 
 lmiterm([lmiPL21 2 1 0],-r);           %#3:-r*I 
 lmiterm([lmiPL21 2 2 S],-r,1);         %#3:-r*S 
% 
 lmiterm([lmiPL21 1 3 R],A1,1);          %#3:A*R
 lmiterm([lmiPL21 1 3 Ck1],B2,1);        %#3:B2*Ck 
 lmiterm([lmiPL21 2 3 Ak1],1,1);         %#3:Ak 
 lmiterm([lmiPL21 1 4 0],A1);            %#3:A 
 lmiterm([lmiPL21 1 4 Dk1],B2,C2);       %#3:B2*Dk*C2 
 lmiterm([lmiPL21 2 4 S],1,A1);          %#3:S*A 
 lmiterm([lmiPL21 2 4 Bk1],1,C2);        %#3:Bk*C2 
% 
 lmiterm([lmiPL21 3 3 R],-r,1);         %#3:-r*R
 lmiterm([lmiPL21 4 3 0],-r);           %#3:-r*I 
 lmiterm([lmiPL21 4 4 S],-r,1);         %#3:-r*S 
%-----
 sth=sin(th); cth=cos(th);
 lmiPL31=newlmi; 
 lmiterm([lmiPL31 1 1 R],sth*A1,1,'s');  %#4:sth*(R*A'+A*R) 
 lmiterm([lmiPL31 1 1 Ck1],sth*B2,1,'s');%#4:sth*(Ck'*B2'+B2*Ck) 
 lmiterm([lmiPL31 2 1 Ak1],sth,1);       %#4:sth*(Ak) 
 lmiterm([lmiPL31 1 2 0],sth*A1);        %#4:sth*(A) 
 lmiterm([lmiPL31 1 2 Dk1],sth*B2,C2);   %#4:sth*(B2*Dk*C2) 
 lmiterm([lmiPL31 2 2 S],1,sth*A1,'s');  %#4:sth*(A'*S+S*A) 
 lmiterm([lmiPL31 2 2 Bk1],1,sth*C2,'s');%#4:sth*(C2'*Bk'+Bk*C2) 
% 
 lmiterm([lmiPL31 1 3 R],cth*A1,1);      %#1:cth*(A*R) 
 lmiterm([lmiPL31 1 3 R],1,-cth*A1');    %#1:cth*(-R*A')
 lmiterm([lmiPL31 1 3 Ck1],cth*B2,1);     %#1:cth*(B*Ck) 
 lmiterm([lmiPL31 1 3 -Ck1],-cth*B2',1);  %#1:cth*(-Ck'*B') 
 lmiterm([lmiPL31 2 3 Ak1],cth,1);       %#4:cth*(Ak) 
 lmiterm([lmiPL31 1 4 -Ak1],-cth,1);     %#4:cth*(-Ak') 
 lmiterm([lmiPL31 1 4 0],A1);            %#4:cth*(A) 
 lmiterm([lmiPL31 2 3 0],-A1');          %#4:cth*(-A') 
 lmiterm([lmiPL31 1 4 Dk1],cth*B2,C2);   %#4:cth*(B2*Dk*C2) 
 lmiterm([lmiPL31 2 3 -Dk1],-cth*C2',B2');%#4:cth*(-C2'*Dk'*B2') 
 lmiterm([lmiPL31 2 4 S],1,cth*A1);      %#4:cth*(S*A) 
 lmiterm([lmiPL31 2 4 S],-cth*A1',1);    %#4:cth*(-A'*S) 
 lmiterm([lmiPL31 2 4 Bk1],1,cth*C2);    %#4:cth*(Bk*C2) 
 lmiterm([lmiPL31 2 4 -Bk1],-cth*C2',1); %#4:cth*(-C2'*Bk') 
% 
 lmiterm([lmiPL31 3 3 R],sth*A1,1,'s');  %#4:sth*(R*A'+A*R)
 lmiterm([lmiPL31 3 3 Ck1],sth*B2,1,'s');%#4:sth*(Ck'*B2'+B2*Ck) 
 lmiterm([lmiPL31 4 3 Ak1],sth,1);       %#4:sth*(Ak) 
 lmiterm([lmiPL31 3 4 0],sth*A1);        %#4:sth*(A) 
 lmiterm([lmiPL31 3 4 Dk1],sth*B2,C2);   %#4:sth*(B2*Dk*C2) 
 lmiterm([lmiPL31 4 4 S],1,sth*A1,'s');  %#4:sth*(A'*S+S*A) 
 lmiterm([lmiPL31 4 4 Bk1],1,sth*C2,'s');%#4:sth*(C2'*Bk'+Bk*C2) 
%=====
 lmiRS2=newlmi;
 lmiterm([lmiRS2 1 1 R],A2,1,'s');       %#1:R*A'+AR 
 lmiterm([lmiRS2 1 1 Ck2],B2,1,'s');     %#1:Ck'*B2'+B2*Ck 
 lmiterm([lmiRS2 2 1 0],A2');            %#1:A' 
 lmiterm([lmiRS2 2 1 Ak2],1,1);          %#1:Ak 
 lmiterm([lmiRS2 2 1 -Dk2],C2',B2');     %#1:C2'*Dk'*B2' 
 lmiterm([lmiRS2 2 2 S],1,A2,'s');       %#1:A'*S+S*A 
 lmiterm([lmiRS2 2 2 Bk2],1,C2,'s');     %#1:C2'*Bk'+Bk*C2 
 lmiterm([lmiRS2 1 3 0],B1);             %#1:B1 
 lmiterm([lmiRS2 1 3 Dk2],B2,D21);       %#1:B2*Dk*D21 
 lmiterm([lmiRS2 2 3 S],1,B1);           %#1:S*B1 
 lmiterm([lmiRS2 2 3 Bk2],1,D21);        %#1:Bk*D21 
 lmiterm([lmiRS2 3 3 gam],-1,1);         %#1:-gam 
 lmiterm([lmiRS2 4 1 R],C1,1);           %#1:C1*R 
 lmiterm([lmiRS2 4 1 Ck2],D12,1);        %#1:D12*Ck 
 lmiterm([lmiRS2 4 2 0],C1);             %#1:C1 
 lmiterm([lmiRS2 4 2 Dk2],D12,C2);       %#1:D12*Dk*C2 
 lmiterm([lmiRS2 4 3 0],D11);            %#1:D11 
 lmiterm([lmiRS2 4 3 Dk2],D12,D21);      %#1:D12*Dk*D21 
 lmiterm([lmiRS2 4 4 gam],-1,1);         %#1:-gam 
%-----
 lmiPL12=newlmi; 
 lmiterm([lmiPL12 1 1 R],A2,1,'s');      %#2:R*A'+A*R 
 lmiterm([lmiPL12 1 1 Ck2],B2,1,'s');    %#2:Ck'*B2'+B2*Ck 
 lmiterm([lmiPL12 2 1 Ak2],1,1);         %#2:Ak 
 lmiterm([lmiPL12 1 2 0],A2);            %#2:A 
 lmiterm([lmiPL12 1 2 Dk2],B2,C2);       %#2:B2*Dk*C2 
 lmiterm([lmiPL12 2 2 S],1,A2,'s');      %#2:A'*S+S*A 
 lmiterm([lmiPL12 2 2 Bk2],1,C2,'s');    %#2:C2'*Bk'+Bk*C2 
% 
 lmiterm([lmiPL12 1 1 R],2*alpha,1);     %#2:2*alpha*R
 lmiterm([lmiPL12 2 1 0],2*alpha);       %#2:2*alpha*I 
 lmiterm([lmiPL12 2 2 S],2*alpha,1);     %#2:2*alpha*S 
%-----
 lmiPL22=newlmi; 
 lmiterm([lmiPL22 1 1 R],-r,1);          %#3:-r*R 
 lmiterm([lmiPL22 2 1 0],-r);            %#3:-r*I 
 lmiterm([lmiPL22 2 2 S],-r,1);          %#3:-r*S 
% 
 lmiterm([lmiPL22 1 3 R],A2,1);          %#3:A*R
 lmiterm([lmiPL22 1 3 Ck2],B2,1);        %#3:B2*Ck 
 lmiterm([lmiPL22 2 3 Ak2],1,1);         %#3:Ak 
 lmiterm([lmiPL22 1 4 0],A2);            %#3:A 
 lmiterm([lmiPL22 1 4 Dk2],B2,C2);       %#3:B2*Dk*C2 
 lmiterm([lmiPL22 2 4 S],1,A2);          %#3:S*A 
 lmiterm([lmiPL22 2 4 Bk2],1,C2);        %#3:Bk*C2 
% 
 lmiterm([lmiPL22 3 3 R],-r,1);          %#3:-r*R
 lmiterm([lmiPL22 4 3 0],-r);            %#3:-r*I 
 lmiterm([lmiPL22 4 4 S],-r,1);          %#3:-r*S 
%-----
 sth=sin(th); cth=cos(th);
 lmiPL32=newlmi; 
 lmiterm([lmiPL32 1 1 R],sth*A2,1,'s');  %#4:sth*(R*A'+A*R) 
 lmiterm([lmiPL32 1 1 Ck2],sth*B2,1,'s');%#4:sth*(Ck'*B2'+B2*Ck) 
 lmiterm([lmiPL32 2 1 Ak2],sth,1);       %#4:sth*(Ak) 
 lmiterm([lmiPL32 1 2 0],sth*A2);        %#4:sth*(A) 
 lmiterm([lmiPL32 1 2 Dk2],sth*B2,C2);   %#4:sth*(B2*Dk*C2) 
 lmiterm([lmiPL32 2 2 S],1,sth*A2,'s');  %#4:sth*(A'*S+S*A) 
 lmiterm([lmiPL32 2 2 Bk2],1,sth*C2,'s');%#4:sth*(C2'*Bk'+Bk*C2) 
% 
 lmiterm([lmiPL32 1 3 R],cth*A2,1);      %#1:cth*(A*R) 
 lmiterm([lmiPL32 1 3 R],1,-cth*A2');    %#1:cth*(-R*A')
 lmiterm([lmiPL32 1 3 Ck2],cth*B2,1);    %#1:cth*(B*Ck) 
 lmiterm([lmiPL32 1 3 -Ck2],-cth*B2',1); %#1:cth*(-Ck'*B') 
 lmiterm([lmiPL32 2 3 Ak2],cth,1);       %#4:cth*(Ak) 
 lmiterm([lmiPL32 1 4 -Ak2],-cth,1);     %#4:cth*(-Ak') 
 lmiterm([lmiPL32 1 4 0],A2);            %#4:cth*(A) 
 lmiterm([lmiPL32 2 3 0],-A2');          %#4:cth*(-A') 
 lmiterm([lmiPL32 1 4 Dk2],cth*B2,C2);   %#4:cth*(B2*Dk*C2) 
 lmiterm([lmiPL32 2 3 -Dk2],-cth*C2',B2');%#4:cth*(-C2'*Dk'*B2') 
 lmiterm([lmiPL32 2 4 S],1,cth*A2);      %#4:cth*(S*A) 
 lmiterm([lmiPL32 2 4 S],-cth*A2',1);    %#4:cth*(-A'*S) 
 lmiterm([lmiPL32 2 4 Bk2],1,cth*C2);    %#4:cth*(Bk*C2) 
 lmiterm([lmiPL32 2 4 -Bk2],-cth*C2',1); %#4:cth*(-C2'*Bk') 
% 
 lmiterm([lmiPL32 3 3 R],sth*A2,1,'s');  %#4:sth*(R*A'+A*R)
 lmiterm([lmiPL32 3 3 Ck2],sth*B2,1,'s');%#4:sth*(Ck'*B2'+B2*Ck) 
 lmiterm([lmiPL32 4 3 Ak2],sth,1);       %#4:sth*(Ak) 
 lmiterm([lmiPL32 3 4 0],sth*A2);        %#4:sth*(A) 
 lmiterm([lmiPL32 3 4 Dk2],sth*B2,C2);   %#4:sth*(B2*Dk*C2) 
 lmiterm([lmiPL32 4 4 S],1,sth*A2,'s');  %#4:sth*(A'*S+S*A) 
 lmiterm([lmiPL32 4 4 Bk2],1,sth*C2,'s');%#4:sth*(C2'*Bk'+Bk*C2)  
%=====
 posX=-newlmi;
 lmiterm([posX 1 1 R],1,1);              %#5:R 
 lmiterm([posX 2 1 0],1);                %#5:I 
 lmiterm([posX 2 2 S],1,1);              %#5:S 
%-----
 lmiDk1=-newlmi;                        
 lmiterm([lmiDk1 1 1 0],1e2);            %#6:1e2 
 lmiterm([lmiDk1 2 2 0],1e2);            %#6:1e2 
 lmiterm([lmiDk1 2 1 Dk1],1,1);          %#6:Dk 
%-----
 lmiDk2=-newlmi;                        
 lmiterm([lmiDk2 1 1 0],1e2);            %#6:1e2 
 lmiterm([lmiDk2 2 2 0],1e2);            %#6:1e2 
 lmiterm([lmiDk2 2 1 Dk2],1,1);          %#6:Dk 
%-----
 lmig=newlmi;                          
 lmiterm([lmig,1,1,gam],1,1);            %#7:gam 
 lmiterm([lmig,1,1,0],-1e3);             %#7:1e3 
 LMIs=getlmis;
end
%-----
%eof

LPV制御（出力FB）

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

LPV制御（出力FB）…Homework

ここでは、主軸変動を伴う回転体の運動制御を例にとって、出力FB形式のGS制御について説明します。

[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} }$

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

$\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{(4)\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{(5)\quad p_1(\Omega(t))+p_2(\Omega(t))=1 }$

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

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

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

$\displaystyle{(7)\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)モデルと呼びます。

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

$\displaystyle{(8)\quad \begin{array}{l} \left\{\begin{array}{l} \dot{x}_K(t)=A_{K1}x_K(t)+B_{K1}y(t) \\ u(t)=C_{K1}x_K(t)+D_{K1}y(t) \end{array}\right.\\ \left\{\begin{array}{l} \dot{x}_K(t)=A_{K2}x_K(t)+B_{K2}y(t) \\ u(t)=C_{K2}x_K(t)+D_{K2}y(t) \end{array}\right. \end{array} }$

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

$\displaystyle{(9)\quad \begin{array}{l} \dot{x}_K(t)=\underbrace{(p_1(\Omega(t))A_{K1}+p_2(\Omega(t))A_{K2})}_{A_K(\Omega(t))}x_K(t)\\ +\underbrace{(p_1(\Omega(t))B_{K1}+p_2(\Omega(t))B_{K2})}_{B_K(\Omega(t))}y(t) \\ u=\underbrace{(p_1(\Omega(t))C_{K1}+p_2(\Omega(t))C_{K2})}_{C_K(\Omega(t))}x_K(t)\\ +\underbrace{(p_1(\Omega(t))D_{K1}+p_2(\Omega(t))D_{K2})}_{D_K(\Omega(t))}y(t) \end{array} }$

のようなLPV(Linear Parameter Varying)制御を考えます。この閉ループ系は次式で表されます。

$\displaystyle{(10)\quad \begin{array}{l} \left[\begin{array}{c} \dot{x}(t) \\ \dot{x}_K(t) \end{array}\right] =(p_1(\Omega(t)) \underbrace{ \left[\begin{array}{ccc} A_1+BD_{K1}C & BC_{K1} \\ B_{K1}C & A_{K1} \end{array}\right] }_{A_{CL1}}\\ +p_2(\Omega(t)) \underbrace{ \left[\begin{array}{ccc} A_2+BD_{K2}C & BC_{K2} \\ B_{K2}C & A_{K2} \end{array}\right] }_{A_{CL2}} )\left[\begin{array}{c} x(t) \\ x_K(t) \end{array}\right]\\ =\underbrace{(p_1(\Omega(t))A_{CL1} +p_2(\Omega(t))A_{CL2} )}_{A_{CL}(\Omega(t))} \left[\begin{array}{c} x(t) \\ x_K(t) \end{array}\right] \end{array} }$

ここで、端点コントローラによる閉ループ系

$\displaystyle{(11)\quad \begin{array}{l} \left[\begin{array}{c} \dot{x}(t) \\ \dot{x}_K(t) \end{array}\right] = \underbrace{ \left[\begin{array}{ccc} A_1+BD_{K1}C & BC_{K1} \\ B_{K1}C & A_{K1} \end{array}\right] }_{A_{CL1}} \left[\begin{array}{c} x(t) \\ x_K(t) \end{array}\right]\\ \left[\begin{array}{c} \dot{x}(t) \\ \dot{x}_K(t) \end{array}\right] = \underbrace{ \left[\begin{array}{ccc} A_2+BD_{K2}C & BC_{K2} \\ B_{K2}C & A_{K1} \end{array}\right] }_{A_{CL2}} \left[\begin{array}{c} x(t) \\ x_K(t) \end{array}\right] \end{array} }$

はそれぞれ漸近安定となるように設計しますが、閉ループ系の時変系としての安定性が保証されるかの検討が必要になります（Note B52参照）。

[2] いま、次の２ポート表現かつポリトピック表現された $n$ 次系を考えます。混乱を避けるため、 $A^{(1)}=A_1$ 、 $A^{(2)}=A_2$ と記法を改めたことに注意してください。

$\displaystyle{(12)\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)=C_2x(t) \end{array}\right. }$

これに対する出力フィードバック

$\displaystyle{(13)\quad K: \left\{\begin{array}{l} \dot{x}_K(t)=\underbrace{(p_1(\Omega(t))A_{K1}+p_2(\Omega(t))A_{K2})}_{A_K(\Omega(t))}x_K(t)\\ +\underbrace{(p_1(\Omega(t))B_{K1}+p_2(\Omega(t))B_{K2})}_{B_K(\Omega(t))}y_2(t) \\ u_2(t)=\underbrace{(p_1(\Omega(t))C_{K1}+p_2(\Omega(t))C_{K2})}_{C_K(\Omega(t))}x_K(t)\\ +\underbrace{(p_1(\Omega(t))D_{K1}+p_2(\Omega(t))D_{K2})}_{D_K(\Omega(t))}y_2(t) \end{array}\right. }$

による閉ループ系

$\displaystyle{(14)\quad P_{CL}: \left\{\begin{array}{l} \left[\begin{array}{c} \dot{x}(t) \\ \dot{x}_K(t) \end{array}\right] =\underbrace{(p_1(\Omega(t))A_{CL1} +p_2(\Omega(t))A_{CL2} )}_{A_{CL}(\Omega(t))} \left[\begin{array}{c} x(t) \\ x_K(t) \end{array}\right]\\ + \underbrace{ \left[\begin{array}{ccc} B_1 \\ 0 \end{array}\right] }_{B_{CL}} u_1(t)\\ y_1(t) = \underbrace{ \left[\begin{array}{ccc} C_1+D_{12}D_KC_2 & D_{12}C_K \end{array}\right] }_{C_{CL}} \left[\begin{array}{c} x(t) \\ x_K(t) \end{array}\right] \end{array}\right. }$

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

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

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

●そのために、各端点における $H_\infty$ 制御（出力フィードバック）を共通のリャプノフ関数をもつように解きます。その際コントローラの実装上の立場から、各端点において

$\lambda(A_{CL1}),\lambda(A_{CL2})\subset{\cal D}_1\cup{\cal D}_2\cup{\cal D}_3$

の制約を課すことがあります。この問題は次のように定式化されます。

Minimize $\gamma$ on $R,S,{\cal A}_K^{(1)},{\cal B}_K^{(1)},{\cal C}_K^{(1)},D_K^{(1)},{\cal A}_K^{(2)},{\cal B}_K^{(2)},{\cal C}_K^{(2)},D_K^{(2)}$ subject to

$\displaystyle{(16) \begin{array}{l} \left[\begin{array}{cc} R & I \\ I & S \end{array}\right] >0 \end{array} }$

$\displaystyle{(17) \left[\begin{array}{ccc} \left[\begin{array}{cc} A^{(1)}R+B_2{\cal C}_K^{(1)} & A^{(1)}+B_2D_K^{(1)}C_2\\ {\cal A}_K^{(1)} & SA^{(1)}+{\cal B}_K^{(1)}C_2 \end{array}\right]+(*)^T & \left[\begin{array}{c} B_1+B_2D_K^{(1)}D_{21} \\ SB_1+{\cal B}_K^{(1)}D_{21} \end{array}\right] & (*)^T \\ (*)^T & -\gamma I & (*)^T \\ \left[\begin{array}{cc} C_1R+D_{12}{\cal C}_K^{(1)} & C_1+D_{12}D_K^{(1)}C_2 \end{array}\right] & D_{11} & -\gamma I \end{array}\right]<0 }$

$\displaystyle{(18) \begin{array}{l} \left[\begin{array}{cc} A^{(1)}R+B_2{\cal C}_K^{(1)} & A^{(1)}+B_2D_K^{(1)}C_2\\ {\cal A}_K^{(1)} & SA^{(1)}+{\cal B}_K^{(1)}C_2 \end{array}\right] +(*)^T %\nonumber\\&& +\alpha \left[\begin{array}{cc} R & I \\ I & S \end{array}\right]<0} \\ \left[\begin{array}{cc} -r \left[\begin{array}{cc} R & I \\ I & S \end{array}\right] & \left[\begin{array}{cc} A^{(1)}R+B_2{\cal C}_K^{(1)} & A^{(1)}+B_2D_K^{(1)}C_2\\ {\cal A}_K^{(1)} & SA^{(1)}+{\cal B}_K^{(1)}C_2 \end{array}\right] \\ (*)^T & -r \left[\begin{array}{cc} R & I \\ I & S \end{array}\right] \end{array}\right] <0 \\ \left[\begin{array}{cc} \sin\theta & \cos\theta \\ -\cos\theta & \sin\theta \end{array}\right]\otimes \left[\begin{array}{cc} A^{(1)}R+B_2{\cal C}_K^{(1)} & A^{(1)}+B_2D_K^{(1)}C_2\\ {\cal A}_K^{(1)} & SA^{(1)}+{\cal B}_K^{(1)}C_2 \end{array}\right]+ (*)^T <0 \end{array} }$

$\displaystyle{(19) \left[\begin{array}{ccc} \left[\begin{array}{cc} A^{(2)}R+B_2{\cal C}_K^{(2)} & A^{(2)}+B_2D_K^{(2)}C_2\\ {\cal A}_K^{(2)} & SA^{(2)}+{\cal B}_K^{(2)}C_2 \end{array}\right]+(*)^T & \left[\begin{array}{c} B_1+B_2D_K^{(2)}D_{21} \\ SB_1+{\cal B}_K^{(2)}D_{21} \end{array}\right] & (*)^T \\ (*)^T & -\gamma I & (*)^T \\ \left[\begin{array}{cc} C_1R+D_{12}{\cal C}_K^{(2)} & C_1+D_{12}D_K^{(2)}C_2 \end{array}\right] & D_{11} & -\gamma I \end{array}\right]<0 }$

$\displaystyle{(20) \begin{array}{l} \left[\begin{array}{cc} A^{(2)}R+B_2{\cal C}_K^{(2)} & A^{(1)}+B_2D_K^{(2)}C_2\\ {\cal A}_K^{(2)} & SA^{(2)}+{\cal B}_K^{(2)}C_2 \end{array}\right] +(*)^T %\nonumber\\&& +\alpha \left[\begin{array}{cc} R & I \\ I & S \end{array}\right]<0} \\ \left[\begin{array}{cc} -r \left[\begin{array}{cc} R & I \\ I & S \end{array}\right] & \left[\begin{array}{cc} A^{(2)}R+B_2{\cal C}_K^{(2)} & A^{(2)}+B_2D_K^{(2)}C_2\\ {\cal A}_K^{(2)} & SA^{(2)}+{\cal B}_K^{(2)}C_2 \end{array}\right] \\ (*)^T & -r \left[\begin{array}{cc} R & I \\ I & S \end{array}\right] \end{array}\right] <0 \\ \left[\begin{array}{cc} \sin\theta & \cos\theta \\ -\cos\theta & \sin\theta \end{array}\right]\otimes \left[\begin{array}{cc} A^{(2)}R+B_2{\cal C}_K^{(2)} & A^{(2)}+B_2D_K^{(2)}C_2\\ {\cal A}_K^{(2)} & SA^{(2)}+{\cal B}_K^{(2)}C_2 \end{array}\right]+ (*)^T <0 \end{array} }$

この最小化問題の解 $R,S,{\cal A}_K^{(1)},{\cal B}_K^{(1)},{\cal C}_K^{(1)},D_K^{(1)},{\cal A}_K^{(2)},{\cal B}_K^{(2)},{\cal C}_K^{(2)},D_K^{(2)}$ を用いて、次式によって出力フィードバック(12)のゲインを決定します。

$\displaystyle{(21)\quad \begin{array}{lll} A_{K1}=N^{-1}({\cal A}_K^{(1)}-S(A-B_2D_K^{(1)}C_2)R-{\cal B}_K^{(1)}C_2R-SB_2{\cal C}_K^{(1)})M^{-T} \nonumber\\ B_{K1}=N^{-1}({\cal B}_K^{(1)}-SB_2D_K^{(1)}) \nonumber\\ C_{K1}=({\cal C}_K^{(1)}-D_K^{(1)}C_2R)M^{-T} \nonumber\\ (I-SR=NM^T) \end{array} }$

$\displaystyle{(22)\quad \begin{array}{lll} A_{K2}=N^{-1}({\cal A}_K^{(2)}-S(A-B_2D_K^{(2)}C_2)R-{\cal B}_K^{(2)}C_2R-SB_2{\cal C}_K^{(2)})M^{-T} \nonumber\\ B_{K2}=N^{-1}({\cal B}_K^{(2)}-SB_2D_K^{(2)}) \nonumber\\ C_{K2}=({\cal C}_K^{(2)}-D_K^{(2)}C_2R)M^{-T} \nonumber\\ (I-SR=NM^T) \end{array} }$

●この問題設定では

$\displaystyle{(23)\quad \left[\begin{array}{ccc} Y_{CL}A_{CL1}^T+A_{CL1}Y_{CL} & B_{CL} & Y_{CL}C_{CL}^T \\ B_{CL}^T & -\gamma I & D_{CL}^T \\ C_{CL}Y_{CL} & D_{CL} & -\gamma I \end{array}\right]<0 }$
$\displaystyle{(24)\quad \left[\begin{array}{ccc} Y_{CL}A_{CL2}^T+A_{CL2}Y_{CL} & B_{CL} & Y_{CL}C_{CL}^T \\ B_{CL}^T & -\gamma I & D_{CL}^T \\ C_{CL}Y_{CL} & D_{CL} & -\gamma I \end{array}\right]<0 }$

の共通解 $Y_{CL}>0$ を求めています。 $X_{CL}=Y_{CL}^{-1}$ ととれば

$\displaystyle{(25)\quad X_{CL}>0 \Leftrightarrow Y_{CL}>0 \Leftrightarrow \left[\begin{array}{cc} R & I \\ I & S \end{array}\right] >0 }$

に注意して

$\displaystyle{(26)\quad \left[\begin{array}{ccc} A_{CL1}^TX_{CL}+X_{CL}A_{CL1} & X_{CL}B_{CL} & C_{CL}^T \\ B_{CL}^TX_{CL} & -\gamma I & D_{CL}^T \\ C_{CL} & D_{CL} & -\gamma　I \end{array}\right]<0 }$
$\displaystyle{(27)\quad \left[\begin{array}{ccc} A_{CL2}^TX_{CL}+X_{CL}A_{CL2} & X_{CL}B_{CL} & C_{CL}^T \\ B_{CL}^TX_{CL} & -\gamma I & D_{CL}^T \\ C_{CL} & D_{CL} & -\gamma　I \end{array}\right]<0 }$

の共通解となっています。これにより、閉ループ系の時変系としての安定性が保証されます（Note B52参照）。

演習B52…Flipped Classroom

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

MATLAB

%cSPIN_of_gs.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")
%-----
 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)];
 C2=eye(2,2); 
 D21=zeros(2,2); 
 D22=zeros(2,2);
%----- 
 alpha=1; r=3; th=pi/4; 
 LMIs=of_synlmi7(A1,A2,B1,B2,C1,C2,D11,D12,D21,D22,alpha,r,th);
 cobj=zeros(1,decnbr(LMIs)); 
 cobj(1)=1; 
 [cost,xopt]=mincx(LMIs,cobj);  
 gopt=dec2mat(LMIs,xopt,1)
 R=dec2mat(LMIs,xopt,2); 
 S=dec2mat(LMIs,xopt,3); 
 %-----
 ak1=dec2mat(LMIs,xopt,4); 
 bk1=dec2mat(LMIs,xopt,5); 
 ck1=dec2mat(LMIs,xopt,6); 
 dk1=dec2mat(LMIs,xopt,7); 
 [u,sd,v]=svd(eye(size(A1,1))-R*S); 
 sd=diag(sqrt(1./diag(sd))); 
 Ni=sd*v'; Mti=u*sd; 
 AK1=Ni*(ak1-S*(A1-B2*dk1*C2)*R-bk1*C2*R-S*B2*ck1)*Mti; 
 BK1=Ni*(bk1-S*B2*dk1); 
 CK1=(ck1-dk1*C2*R)*Mti; 
 DK1=dk1; 
 K1=[AK1 BK1;CK1 DK1];
 ACL1=[A1+B2*DK1*C2 B2*CK1;
      BK1*C2 AK1];
 plCL1=eig(ACL1)
 %-----
 ak2=dec2mat(LMIs,xopt,8); 
 bk2=dec2mat(LMIs,xopt,9); 
 ck2=dec2mat(LMIs,xopt,10); 
 dk2=dec2mat(LMIs,xopt,11); 
 [u,sd,v]=svd(eye(size(A2,1))-R*S); 
 sd=diag(sqrt(1./diag(sd))); 
 Ni=sd*v'; Mti=u*sd; 
 AK2=Ni*(ak2-S*(A2-B2*dk2*C2)*R-bk2*C2*R-S*B2*ck2)*Mti; 
 BK2=Ni*(bk2-S*B2*dk2); 
 CK2=(ck2-dk2*C2*R)*Mti; 
 DK2=dk2; 
 K2=[AK2 BK2;CK2 DK2];
 ACL2=[A2+B2*DK2*C2 B2*CK2;
       BK2*C2 AK2];
 plCL2=eig(ACL2)
%------
 figure(1) 
 subplot(121),dregion(alpha,0,r,th,r*[-1,1,-1,1]) 
 plot(real(plCL1),imag(plCL1),'*') 
 subplot(122), dregion(alpha,0,r,th,r*[-1,1,-1,1]) 
 plot(real(plCL2),imag(plCL2),'*')  
%-----
 prange=OMmax-OMmin; pmax=OMmax; pmin=OMmin; 
%sim("SPIN_of_gs") 
 sim("SPIN_of_gs_2015a")  
%-----
function LMIs=of_synlmi7(A1,A2,B1,B2,C1,C2,D11,D12,D21,D22,alpha,r,th)
 [n,m]=size(B2);  [p,n]=size(C2);
 setlmis([]); 
 gam=lmivar(1,[1 0]);  
 [R,xxx,Rdec]=lmivar(1,[n 1]); 
 [S,xxx,Sdec]=lmivar(1,[n 1]); 
 Ak1=lmivar(2,[n n]); 
 Bk1=lmivar(2,[n p]); 
 Ck1=lmivar(2,[m n]); 
 Dk1=lmivar(2,[m p]); 
 Ak2=lmivar(2,[n n]); 
 Bk2=lmivar(2,[n p]); 
 Ck2=lmivar(2,[m n]); 
 Dk2=lmivar(2,[m p]);  
%=====
 lmiRS1=newlmi;
 lmiterm([lmiRS1 1 1 R],A1,1,'s');       %#1:R*A'+AR 
 lmiterm([lmiRS1 1 1 Ck1],B2,1,'s');     %#1:Ck'*B2'+B2*Ck 
 lmiterm([lmiRS1 2 1 0],A1');            %#1:A' 
 lmiterm([lmiRS1 2 1 Ak1],1,1);          %#1:Ak 
 lmiterm([lmiRS1 2 1 -Dk1],C2',B2');     %#1:C2'*Dk'*B2' 
 lmiterm([lmiRS1 2 2 S],1,A1,'s');       %#1:A'*S+S*A 
 lmiterm([lmiRS1 2 2 Bk1],1,C2,'s');     %#1:C2'*Bk'+Bk*C2 
 lmiterm([lmiRS1 1 3 0],B1);             %#1:B1 
 lmiterm([lmiRS1 1 3 Dk1],B2,D21);       %#1:B2*Dk*D21 
 lmiterm([lmiRS1 2 3 S],1,B1);           %#1:S*B1 
 lmiterm([lmiRS1 2 3 Bk1],1,D21);        %#1:Bk*D21 
 lmiterm([lmiRS1 3 3 gam],-1,1);         %#1:-gam 
 lmiterm([lmiRS1 4 1 R],C1,1);           %#1:C1*R 
 lmiterm([lmiRS1 4 1 Ck1],D12,1);        %#1:D12*Ck 
 lmiterm([lmiRS1 4 2 0],C1);             %#1:C1 
 lmiterm([lmiRS1 4 2 Dk1],D12,C2);       %#1:D12*Dk*C2 
 lmiterm([lmiRS1 4 3 0],D11);            %#1:D11 
 lmiterm([lmiRS1 4 3 Dk1],D12,D21);      %#1:D12*Dk*D21 
 lmiterm([lmiRS1 4 4 gam],-1,1);         %#1:-gam 
%-----
 lmiPL11=newlmi; 
 lmiterm([lmiPL11 1 1 R],A1,1,'s');      %#2:R*A'+A*R 
 lmiterm([lmiPL11 1 1 Ck1],B2,1,'s');    %#2:Ck'*B2'+B2*Ck 
 lmiterm([lmiPL11 2 1 Ak1],1,1);         %#2:Ak 
 lmiterm([lmiPL11 1 2 0],A1);            %#2:A 
 lmiterm([lmiPL11 1 2 Dk1],B2,C2);       %#2:B2*Dk*C2 
 lmiterm([lmiPL11 2 2 S],1,A1,'s');      %#2:A'*S+S*A 
 lmiterm([lmiPL11 2 2 Bk1],1,C2,'s');    %#2:C2'*Bk'+Bk*C2 
% 
 lmiterm([lmiPL11 1 1 R],2*alpha,1);     %#2:2*alpha*R
 lmiterm([lmiPL11 2 1 0],2*alpha);       %#2:2*alpha*I 
 lmiterm([lmiPL11 2 2 S],2*alpha,1);     %#2:2*alpha*S 
%-----
 lmiPL21=newlmi; 
 lmiterm([lmiPL21 1 1 R],-r,1);          %#3:-r*R 
 lmiterm([lmiPL21 2 1 0],-r);            %#3:-r*I 
 lmiterm([lmiPL21 2 2 S],-r,1);          %#3:-r*S 
% 
 lmiterm([lmiPL21 1 3 R],A1,1);          %#3:A*R
 lmiterm([lmiPL21 1 3 Ck1],B2,1);        %#3:B2*Ck 
 lmiterm([lmiPL21 2 3 Ak1],1,1);         %#3:Ak 
 lmiterm([lmiPL21 1 4 0],A1);            %#3:A 
 lmiterm([lmiPL21 1 4 Dk1],B2,C2);       %#3:B2*Dk*C2 
 lmiterm([lmiPL21 2 4 S],1,A1);          %#3:S*A 
 lmiterm([lmiPL21 2 4 Bk1],1,C2);        %#3:Bk*C2 
% 
 lmiterm([lmiPL21 3 3 R],-r,1);          %#3:-r*R
 lmiterm([lmiPL21 4 3 0],-r);            %#3:-r*I 
 lmiterm([lmiPL21 4 4 S],-r,1);          %#3:-r*S 
%-----
 sth=sin(th); cth=cos(th);
 lmiPL31=newlmi; 
 lmiterm([lmiPL31 1 1 R],sth*A1,1,'s');  %#4:sth*(R*A'+A*R) 
 lmiterm([lmiPL31 1 1 Ck1],sth*B2,1,'s');%#4:sth*(Ck'*B2'+B2*Ck) 
 lmiterm([lmiPL31 2 1 Ak1],sth,1);       %#4:sth*(Ak) 
 lmiterm([lmiPL31 1 2 0],sth*A1);        %#4:sth*(A) 
 lmiterm([lmiPL31 1 2 Dk1],sth*B2,C2);   %#4:sth*(B2*Dk*C2) 
 lmiterm([lmiPL31 2 2 S],1,sth*A1,'s');  %#4:sth*(A'*S+S*A) 
 lmiterm([lmiPL31 2 2 Bk1],1,sth*C2,'s');%#4:sth*(C2'*Bk'+Bk*C2) 
% 
 lmiterm([lmiPL31 1 3 R],cth*A1,1);      %#1:cth*(A*R) 
 lmiterm([lmiPL31 1 3 R],1,-cth*A1');    %#1:cth*(-R*A')
 lmiterm([lmiPL31 1 3 Ck1],cth*B2,1);    %#1:cth*(B*Ck) 
 lmiterm([lmiPL31 1 3 -Ck1],-cth*B2',1); %#1:cth*(-Ck'*B') 
 lmiterm([lmiPL31 2 3 Ak1],cth,1);       %#4:cth*(Ak) 
 lmiterm([lmiPL31 1 4 -Ak1],-cth,1);     %#4:cth*(-Ak') 
 lmiterm([lmiPL31 1 4 0],A1);            %#4:cth*(A) 
 lmiterm([lmiPL31 2 3 0],-A1');          %#4:cth*(-A') 
 lmiterm([lmiPL31 1 4 Dk1],cth*B2,C2);   %#4:cth*(B2*Dk*C2) 
 lmiterm([lmiPL31 2 3 -Dk1],-cth*C2',B2');%#4:cth*(-C2'*Dk'*B2') 
 lmiterm([lmiPL31 2 4 S],1,cth*A1);      %#4:cth*(S*A) 
 lmiterm([lmiPL31 2 4 S],-cth*A1',1);    %#4:cth*(-A'*S) 
 lmiterm([lmiPL31 2 4 Bk1],1,cth*C2);    %#4:cth*(Bk*C2) 
 lmiterm([lmiPL31 2 4 -Bk1],-cth*C2',1); %#4:cth*(-C2'*Bk') 
% 
 lmiterm([lmiPL31 3 3 R],sth*A1,1,'s');  %#4:sth*(R*A'+A*R)
 lmiterm([lmiPL31 3 3 Ck1],sth*B2,1,'s');%#4:sth*(Ck'*B2'+B2*Ck) 
 lmiterm([lmiPL31 4 3 Ak1],sth,1);       %#4:sth*(Ak) 
 lmiterm([lmiPL31 3 4 0],sth*A1);        %#4:sth*(A) 
 lmiterm([lmiPL31 3 4 Dk1],sth*B2,C2);   %#4:sth*(B2*Dk*C2) 
 lmiterm([lmiPL31 4 4 S],1,sth*A1,'s');  %#4:sth*(A'*S+S*A) 
 lmiterm([lmiPL31 4 4 Bk1],1,sth*C2,'s');%#4:sth*(C2'*Bk'+Bk*C2) 
%=====
 lmiRS2=newlmi;
 lmiterm([lmiRS2 1 1 R],A2,1,'s');       %#1:R*A'+AR 
 lmiterm([lmiRS2 1 1 Ck2],B2,1,'s');     %#1:Ck'*B2'+B2*Ck 
 lmiterm([lmiRS2 2 1 0],A2');            %#1:A' 
 lmiterm([lmiRS2 2 1 Ak2],1,1);          %#1:Ak 
 lmiterm([lmiRS2 2 1 -Dk2],C2',B2');     %#1:C2'*Dk'*B2' 
 lmiterm([lmiRS2 2 2 S],1,A2,'s');       %#1:A'*S+S*A 
 lmiterm([lmiRS2 2 2 Bk2],1,C2,'s');     %#1:C2'*Bk'+Bk*C2 
 lmiterm([lmiRS2 1 3 0],B1);            %#1:B1 
 lmiterm([lmiRS2 1 3 Dk2],B2,D21);       %#1:B2*Dk*D21 
 lmiterm([lmiRS2 2 3 S],1,B1);          %#1:S*B1 
 lmiterm([lmiRS2 2 3 Bk2],1,D21);        %#1:Bk*D21 
 lmiterm([lmiRS2 3 3 gam],-1,1);        %#1:-gam 
 lmiterm([lmiRS2 4 1 R],C1,1);          %#1:C1*R 
 lmiterm([lmiRS2 4 1 Ck2],D12,1);        %#1:D12*Ck 
 lmiterm([lmiRS2 4 2 0],C1);            %#1:C1 
 lmiterm([lmiRS2 4 2 Dk2],D12,C2);       %#1:D12*Dk*C2 
 lmiterm([lmiRS2 4 3 0],D11);           %#1:D11 
 lmiterm([lmiRS2 4 3 Dk2],D12,D21);      %#1:D12*Dk*D21 
 lmiterm([lmiRS2 4 4 gam],-1,1);        %#1:-gam 
%-----
 lmiPL12=newlmi; 
 lmiterm([lmiPL12 1 1 R],A2,1,'s');      %#2:R*A'+A*R 
 lmiterm([lmiPL12 1 1 Ck2],B2,1,'s');    %#2:Ck'*B2'+B2*Ck 
 lmiterm([lmiPL12 2 1 Ak2],1,1);         %#2:Ak 
 lmiterm([lmiPL12 1 2 0],A2);            %#2:A 
 lmiterm([lmiPL12 1 2 Dk2],B2,C2);       %#2:B2*Dk*C2 
 lmiterm([lmiPL12 2 2 S],1,A2,'s');      %#2:A'*S+S*A 
 lmiterm([lmiPL12 2 2 Bk2],1,C2,'s');    %#2:C2'*Bk'+Bk*C2 
% 
 lmiterm([lmiPL12 1 1 R],2*alpha,1);     %#2:2*alpha*R
 lmiterm([lmiPL12 2 1 0],2*alpha);       %#2:2*alpha*I 
 lmiterm([lmiPL12 2 2 S],2*alpha,1);     %#2:2*alpha*S 
%-----
 lmiPL22=newlmi; 
 lmiterm([lmiPL22 1 1 R],-r,1);          %#3:-r*R 
 lmiterm([lmiPL22 2 1 0],-r);            %#3:-r*I 
 lmiterm([lmiPL22 2 2 S],-r,1);          %#3:-r*S 
% 
 lmiterm([lmiPL22 1 3 R],A2,1);          %#3:A*R
 lmiterm([lmiPL22 1 3 Ck2],B2,1);        %#3:B2*Ck 
 lmiterm([lmiPL22 2 3 Ak2],1,1);         %#3:Ak 
 lmiterm([lmiPL22 1 4 0],A2);            %#3:A 
 lmiterm([lmiPL22 1 4 Dk2],B2,C2);       %#3:B2*Dk*C2 
 lmiterm([lmiPL22 2 4 S],1,A2);          %#3:S*A 
 lmiterm([lmiPL22 2 4 Bk2],1,C2);        %#3:Bk*C2 
% 
 lmiterm([lmiPL22 3 3 R],-r,1);          %#3:-r*R
 lmiterm([lmiPL22 4 3 0],-r);            %#3:-r*I 
 lmiterm([lmiPL22 4 4 S],-r,1);          %#3:-r*S 
%-----
 sth=sin(th); cth=cos(th);
 lmiPL32=newlmi; 
 lmiterm([lmiPL32 1 1 R],sth*A2,1,'s');  %#4:sth*(R*A'+A*R) 
 lmiterm([lmiPL32 1 1 Ck2],sth*B2,1,'s');%#4:sth*(Ck'*B2'+B2*Ck) 
 lmiterm([lmiPL32 2 1 Ak2],sth,1);       %#4:sth*(Ak) 
 lmiterm([lmiPL32 1 2 0],sth*A2);        %#4:sth*(A) 
 lmiterm([lmiPL32 1 2 Dk2],sth*B2,C2);   %#4:sth*(B2*Dk*C2) 
 lmiterm([lmiPL32 2 2 S],1,sth*A2,'s');  %#4:sth*(A'*S+S*A) 
 lmiterm([lmiPL32 2 2 Bk2],1,sth*C2,'s');%#4:sth*(C2'*Bk'+Bk*C2) 
% 
 lmiterm([lmiPL32 1 3 R],cth*A2,1);      %#1:cth*(A*R) 
 lmiterm([lmiPL32 1 3 R],1,-cth*A2');    %#1:cth*(-R*A')
 lmiterm([lmiPL32 1 3 Ck2],cth*B2,1);    %#1:cth*(B*Ck) 
 lmiterm([lmiPL32 1 3 -Ck2],-cth*B2',1); %#1:cth*(-Ck'*B') 
 lmiterm([lmiPL32 2 3 Ak2],cth,1);       %#4:cth*(Ak) 
 lmiterm([lmiPL32 1 4 -Ak2],-cth,1);     %#4:cth*(-Ak') 
 lmiterm([lmiPL32 1 4 0],A2);            %#4:cth*(A) 
 lmiterm([lmiPL32 2 3 0],-A2');          %#4:cth*(-A') 
 lmiterm([lmiPL32 1 4 Dk2],cth*B2,C2);   %#4:cth*(B2*Dk*C2) 
 lmiterm([lmiPL32 2 3 -Dk2],-cth*C2',B2');%#4:cth*(-C2'*Dk'*B2') 
 lmiterm([lmiPL32 2 4 S],1,cth*A2);      %#4:cth*(S*A) 
 lmiterm([lmiPL32 2 4 S],-cth*A2',1);    %#4:cth*(-A'*S) 
 lmiterm([lmiPL32 2 4 Bk2],1,cth*C2);    %#4:cth*(Bk*C2) 
 lmiterm([lmiPL32 2 4 -Bk2],-cth*C2',1); %#4:cth*(-C2'*Bk') 
% 
 lmiterm([lmiPL32 3 3 R],sth*A2,1,'s');  %#4:sth*(R*A'+A*R)
 lmiterm([lmiPL32 3 3 Ck2],sth*B2,1,'s');%#4:sth*(Ck'*B2'+B2*Ck) 
 lmiterm([lmiPL32 4 3 Ak2],sth,1);       %#4:sth*(Ak) 
 lmiterm([lmiPL32 3 4 0],sth*A2);        %#4:sth*(A) 
 lmiterm([lmiPL32 3 4 Dk2],sth*B2,C2);   %#4:sth*(B2*Dk*C2) 
 lmiterm([lmiPL32 4 4 S],1,sth*A2,'s');  %#4:sth*(A'*S+S*A) 
 lmiterm([lmiPL32 4 4 Bk2],1,sth*C2,'s');%#4:sth*(C2'*Bk'+Bk*C2)  
%=====
 posX=-newlmi;
 lmiterm([posX 1 1 R],1,1);              %#5:R 
 lmiterm([posX 2 1 0],1);                %#5:I 
 lmiterm([posX 2 2 S],1,1);              %#5:S 
%-----
 lmiDk1=-newlmi;                        
 lmiterm([lmiDk1 1 1 0],1e2);            %#6:1e2 
 lmiterm([lmiDk1 2 2 0],1e2);            %#6:1e2 
 lmiterm([lmiDk1 2 1 Dk1],1,1);          %#6:Dk 
%-----
 lmiDk2=-newlmi;                        
 lmiterm([lmiDk2 1 1 0],1e2);            %#6:1e2 
 lmiterm([lmiDk2 2 2 0],1e2);            %#6:1e2 
 lmiterm([lmiDk2 2 1 Dk2],1,1);          %#6:Dk 
%-----
 lmig=newlmi;                          
 lmiterm([lmig,1,1,gam],1,1);            %#7:gam 
 lmiterm([lmig,1,1,0],-1e3);             %#7:1e3 
 LMIs=getlmis;
end
%-----
%eof

図2 回転体のLPV制御

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

MATLAB

%cSPIN_of_hinf.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")
%-----
 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)];
 C2=eye(2,2); 
 D21=zeros(2,2); 
 D22=zeros(2,2);
 %-----
 alpha=1; r=3; th=pi/4; 
 LMIs=of_synlmi6(A,B1,B2,C1,C2,D11,D12,D21,D22,alpha,r,th);
 cobj=zeros(1,decnbr(LMIs)); 
 cobj(1)=1; 
 [cost,xopt]=mincx(LMIs,cobj);  
 gopt=dec2mat(LMIs,xopt,1)
 R=dec2mat(LMIs,xopt,2); 
 S=dec2mat(LMIs,xopt,3); 
 ak=dec2mat(LMIs,xopt,4); 
 bk=dec2mat(LMIs,xopt,5); 
 ck=dec2mat(LMIs,xopt,6); 
 dk=dec2mat(LMIs,xopt,7); 
 [u,sd,v]=svd(eye(size(A,1))-R*S); 
 sd=diag(sqrt(1./diag(sd))); 
 Ni=sd*v'; Mti=u*sd; 
 AK=Ni*(ak-S*(A-B2*dk*C2)*R-bk*C2*R-S*B2*ck)*Mti; 
 BK=Ni*(bk-S*B2*dk); 
 CK=(ck-dk*C2*R)*Mti; 
 DK=dk; 
 K=[AK BK;CK DK];
%-----
 pl=eig(A)
 ACL=[A+B2*DK*C2 B2*CK;
      BK*C2 AK];
 plCL=eig(ACL)
 figure(1)
 dregion(alpha,0,r,th,r*[-1,1,-1,1])  
 plot(real(pl),imag(pl),'x',real(plCL),imag(plCL),'*')
%sim("SPIN_of_hinf")
 sim("SPIN_of_hinf_2015a") 
%-----
function LMIs=of_synlmi6(A,B1,B2,C1,C2,D11,D12,D21,D22,alpha,r,th)
 [n,m]=size(B2);  [p,n]=size(C2);
 setlmis([]); 
 gam=lmivar(1,[1 0]);  
 [R,xxx,Rdec]=lmivar(1,[n 1]); 
 [S,xxx,Sdec]=lmivar(1,[n 1]); 
 Ak=lmivar(2,[n n]); 
 Bk=lmivar(2,[n p]); 
 Ck=lmivar(2,[m n]); 
 Dk=lmivar(2,[m p]); 
% 
 lmiRS=newlmi;
 lmiterm([lmiRS 1 1 R],A,1,'s');       %#1:R*A'+AR 
 lmiterm([lmiRS 1 1 Ck],B2,1,'s');     %#1:Ck'*B2'+B2*Ck 
 lmiterm([lmiRS 2 1 0],A');            %#1:A' 
 lmiterm([lmiRS 2 1 Ak],1,1);          %#1:Ak 
 lmiterm([lmiRS 2 1 -Dk],C2',B2');     %#1:C2'*Dk'*B2' 
 lmiterm([lmiRS 2 2 S],1,A,'s');       %#1:A'*S+S*A 
 lmiterm([lmiRS 2 2 Bk],1,C2,'s');     %#1:C2'*Bk'+Bk*C2 
 lmiterm([lmiRS 1 3 0],B1);            %#1:B1 
 lmiterm([lmiRS 1 3 Dk],B2,D21);       %#1:B2*Dk*D21 
 lmiterm([lmiRS 2 3 S],1,B1);          %#1:S*B1 
 lmiterm([lmiRS 2 3 Bk],1,D21);        %#1:Bk*D21 
 lmiterm([lmiRS 3 3 gam],-1,1);        %#1:-gam 
 lmiterm([lmiRS 4 1 R],C1,1);          %#1:C1*R 
 lmiterm([lmiRS 4 1 Ck],D12,1);        %#1:D12*Ck 
 lmiterm([lmiRS 4 2 0],C1);            %#1:C1 
 lmiterm([lmiRS 4 2 Dk],D12,C2);       %#1:D12*Dk*C2 
 lmiterm([lmiRS 4 3 0],D11);           %#1:D11 
 lmiterm([lmiRS 4 3 Dk],D12,D21);      %#1:D12*Dk*D21 
 lmiterm([lmiRS 4 4 gam],-1,1);        %#1:-gam 
%-----
 lmiPL1=newlmi; 
 lmiterm([lmiPL1 1 1 R],A,1,'s');      %#2:R*A'+A*R 
 lmiterm([lmiPL1 1 1 Ck],B2,1,'s');    %#2:Ck'*B2'+B2*Ck 
 lmiterm([lmiPL1 2 1 Ak],1,1);         %#2:Ak 
 lmiterm([lmiPL1 1 2 0],A);            %#2:A 
 lmiterm([lmiPL1 1 2 Dk],B2,C2);       %#2:B2*Dk*C2 
 lmiterm([lmiPL1 2 2 S],1,A,'s');      %#2:A'*S+S*A 
 lmiterm([lmiPL1 2 2 Bk],1,C2,'s');    %#2:C2'*Bk'+Bk*C2 
% 
 lmiterm([lmiPL1 1 1 R],2*alpha,1);    %#2:2*alpha*R
 lmiterm([lmiPL1 2 1 0],2*alpha);      %#2:2*alpha*I 
 lmiterm([lmiPL1 2 2 S],2*alpha,1);    %#2:2*alpha*S 
%-----
 lmiPL2=newlmi; 
 lmiterm([lmiPL2 1 1 R],-r,1);         %#3:-r*R 
 lmiterm([lmiPL2 2 1 0],-r);           %#3:-r*I 
 lmiterm([lmiPL2 2 2 S],-r,1);         %#3:-r*S 
% 
 lmiterm([lmiPL2 1 3 R],A,1);          %#3:A*R
 lmiterm([lmiPL2 1 3 Ck],B2,1);        %#3:B2*Ck 
 lmiterm([lmiPL2 2 3 Ak],1,1);         %#3:Ak 
 lmiterm([lmiPL2 1 4 0],A);            %#3:A 
 lmiterm([lmiPL2 1 4 Dk],B2,C2);       %#3:B2*Dk*C2 
 lmiterm([lmiPL2 2 4 S],1,A);          %#3:S*A 
 lmiterm([lmiPL2 2 4 Bk],1,C2);        %#3:Bk*C2 
% 
 lmiterm([lmiPL2 3 3 R],-r,1);         %#3:-r*R
 lmiterm([lmiPL2 4 3 0],-r);           %#3:-r*I 
 lmiterm([lmiPL2 4 4 S],-r,1);         %#3:-r*S 
%-----
 sth=sin(th); cth=cos(th);
 lmiPL3=newlmi; 
 lmiterm([lmiPL3 1 1 R],sth*A,1,'s');  %#4:sth*(R*A'+A*R) 
 lmiterm([lmiPL3 1 1 Ck],sth*B2,1,'s');%#4:sth*(Ck'*B2'+B2*Ck) 
 lmiterm([lmiPL3 2 1 Ak],sth,1);       %#4:sth*(Ak) 
 lmiterm([lmiPL3 1 2 0],sth*A);        %#4:sth*(A) 
 lmiterm([lmiPL3 1 2 Dk],sth*B2,C2);   %#4:sth*(B2*Dk*C2) 
 lmiterm([lmiPL3 2 2 S],1,sth*A,'s');  %#4:sth*(A'*S+S*A) 
 lmiterm([lmiPL3 2 2 Bk],1,sth*C2,'s');%#4:sth*(C2'*Bk'+Bk*C2) 
% 
%  lmiterm([lmiPL3 1 3 R],1,cth*A');     %#4:cth*(R*A')
%  lmiterm([lmiPL3 1 3 R],-cth*A,1);     %#4:cth*(-A*R) 
%  lmiterm([lmiPL3 1 3 -Ck],cth*B2',1);  %#4:cth*(Ck'*B2') 
%  lmiterm([lmiPL3 1 3 Ck],-cth*B2,1);   %#4:cth*(-B2*Ck) 
 lmiterm([lmiPL3 1 3 R],cth*A,1);      %#1:cth*(A*R) 
 lmiterm([lmiPL3 1 3 R],1,-cth*A');    %#1:cth*(-R*A')
 lmiterm([lmiPL3 1 3 Ck],cth*B2,1);     %#1:cth*(B*Ck) 
 lmiterm([lmiPL3 1 3 -Ck],-cth*B2',1);  %#1:cth*(-Ck'*B') 
 lmiterm([lmiPL3 2 3 Ak],cth,1);       %#4:cth*(Ak) 
 lmiterm([lmiPL3 1 4 -Ak],-cth,1);     %#4:cth*(-Ak') 
 lmiterm([lmiPL3 1 4 0],A);            %#4:cth*(A) 
 lmiterm([lmiPL3 2 3 0],-A');          %#4:cth*(-A') 
 lmiterm([lmiPL3 1 4 Dk],cth*B2,C2);   %#4:cth*(B2*Dk*C2) 
 lmiterm([lmiPL3 2 3 -Dk],-cth*C2',B2');%#4:cth*(-C2'*Dk'*B2') 
 lmiterm([lmiPL3 2 4 S],1,cth*A);      %#4:cth*(S*A) 
 lmiterm([lmiPL3 2 4 S],-cth*A',1);    %#4:cth*(-A'*S) 
 lmiterm([lmiPL3 2 4 Bk],1,cth*C2);    %#4:cth*(Bk*C2) 
 lmiterm([lmiPL3 2 4 -Bk],-cth*C2',1); %#4:cth*(-C2'*Bk') 
% 
 lmiterm([lmiPL3 3 3 R],sth*A,1,'s');  %#4:sth*(R*A'+A*R)
 lmiterm([lmiPL3 3 3 Ck],sth*B2,1,'s');%#4:sth*(Ck'*B2'+B2*Ck) 
 lmiterm([lmiPL3 4 3 Ak],sth,1);       %#4:sth*(Ak) 
 lmiterm([lmiPL3 3 4 0],sth*A);        %#4:sth*(A) 
 lmiterm([lmiPL3 3 4 Dk],sth*B2,C2);   %#4:sth*(B2*Dk*C2) 
 lmiterm([lmiPL3 4 4 S],1,sth*A,'s');  %#4:sth*(A'*S+S*A) 
 lmiterm([lmiPL3 4 4 Bk],1,sth*C2,'s');%#4:sth*(C2'*Bk'+Bk*C2) 
%-----
 posX=-newlmi;
 lmiterm([posX 1 1 R],1,1);            %#5:R 
 lmiterm([posX 2 1 0],1);              %#5:I 
 lmiterm([posX 2 2 S],1,1);            %#5:S 
%-----
 lmiDk=-newlmi;                        
 lmiterm([lmiDk 1 1 0],1e2);           %#6:1e2 
 lmiterm([lmiDk 2 2 0],1e2);           %#6:1e2 
 lmiterm([lmiDk 2 1 Dk],1,1);          %#6:Dk 
%-----
 lmig=newlmi;                          
 lmiterm([lmig,1,1,gam],1,1);          %#7:gam 
 lmiterm([lmig,1,1,0],-1e3);           %#7:1e3 
 LMIs=getlmis;
end
%-----
%eof

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

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

Note B52　閉ループ系の時変系としての安定性

●次の時変自由系の漸近安定性を考えます。

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

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

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

を構成して

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

を示すことができれば、平衡状態は漸近安定となることが知られています。

●LPVモデル

$\displaystyle{(4)\quad \underbrace{\left[\begin{array}{c} \dot{x}(t) \\ \dot{x}_K(t) \end{array}\right]}_{\dot{\bar x}(t)} =\underbrace{(p_1(\Omega(t))A_{CL1} +p_2(\Omega(t))A_{CL2} )}_{A_{CL}(\Omega(t))} \underbrace{\left[\begin{array}{c} x(t) \\ x_K(t) \end{array}\right]}_{\bar{x}(t)} }$

において、 $A_{CL1}$ と $A_{CL2}$ を個別に安定行列としたとすると

$\displaystyle{(5)\quad \exists X_1>0:X_1A_{CL1}+A_{CL1}^TX_1<0 }$
$\displaystyle{(6)\quad \exists X_2>0:X_2A_{CL2}+A_{CL2}^TX_2<0 }$

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

$\displaystyle{(7)\quad V_1(\bar{x}(t))=\bar{x}^T(t)X_1\bar{x}(t)}$
$\displaystyle{(8)\quad V_2(\bar{x}(t))=\bar{x}^T(t)X_2\bar{x}(t)}$

が得られますが、(4)のリャプノフ関数を構成したことにはなりません。一方

$\displaystyle{(9)\quad \exists X>0: \begin{array}{l} XA_{CL1}+A_{CL1}^TX<0\\ XA_{CL2}+A_{CL2}^TX<0 \end{array} }$

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

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

に対して

$\displaystyle{(11)\quad \begin{array}{l} \frac{d}{dt}V(\bar{x}(t))=\bar{x}^T(t)(A_{CL}^T(\Omega(t))X+XA_{CL}(\Omega(t)))\bar{x}(t)\\\ =\bar{x}^T(t)( ( ( p_1(\Omega(t))A_{CL1}+p_2(\Omega(t))A_{CL2} )^TX\\ +X( p_1(\Omega(t))A_{CL1}+p_2(\Omega(t))A_{CL2} ) )\bar{x}(t)\\ =p_1(\Omega(t)) \underbrace{\bar{x}^T(t)(A_{CL1}^TX+XA_{CL1})\bar{x}(t)}_{<0}\\ +p_2(\Omega(t)) \underbrace{\bar{x}^T(t)(A_{CL2}^TX+XA_{CL2})\bar{x}(t)}_{<0}<0 \end{array} }$

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

●LPVモデル

$\displaystyle{(12)\quad P_{CL}: \left\{\begin{array}{l} \underbrace{\left[\begin{array}{c} \dot{x}(t) \\ \dot{x}_K(t) \end{array}\right]}_{\dot{\bar x}(t)} =\underbrace{(p_1(\Omega(t))A_{CL1} +p_2(\Omega(t))A_{CL2} )}_{A_{CL}(\Omega(t))} \underbrace{\left[\begin{array}{c} x(t) \\ x_K(t) \end{array}\right]}_{\bar{x}(t)}\\ + \underbrace{ \left[\begin{array}{ccc} B_1 \\ 0 \end{array}\right] }_{B_{CL}} u_1(t)\\ y_1(t) = \underbrace{ \left[\begin{array}{ccc} C_1+D_{12}D_KC_2 & D_{12}C_K \end{array}\right] }_{C_{CL}} \underbrace{\left[\begin{array}{c} x(t) \\ x_K(t) \end{array}\right]}_{\bar{x}(t)} \end{array}\right. }$

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

$\displaystyle{(13)\quad \exists X_1>0:\ \left[\begin{array}{ccc} A_{CL1}^TX_1+X_1A_{CL1} & X_1B_{CL} & C_{CL}^T \\ B_{CL}^TX_1 & -\gamma I & D_{CL}^T \\ C_{CL} & D_{CL} & -\gamma　I \end{array}\right]<0 }$
$\displaystyle{(14)\quad \exists X_2>0:\ \left[\begin{array}{ccc} A_{CL2}^TX_2+X_2A_{CL2} & X_2B_{CL} & C_{CL}^T \\ B_{CL}^TX_2 & -\gamma I & D_{CL}^T \\ C_{CL} & D_{CL} & -\gamma　I \end{array}\right]<0 }$

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

$\displaystyle{(15)\quad \left[\begin{array}{ccc} A_{CL1}^TX+XA_{CL1} & XB_{CL} & C_{CL}^T \\ B_{CL}^TX & -\gamma^2 I & D_{CL}^T \\ C_{CL} & D_{CL} &-I \end{array}\right]<0 }$
$\displaystyle{ \Leftrightarrow \left[\begin{array}{cc} A_{CL1}^TX+XA_{CL1} & XB_{CL} \\ B_{CL}^TX & 0 \end{array}\right] < \left[\begin{array}{cc} C_{CL}^T & 0 \\ D_{CL}^T & I \end{array}\right] \left[\begin{array}{cc} -I & 0 \\ 0 & \gamma^2 I \end{array}\right] \left[\begin{array}{cc} C_{CL} & D_{CL} \\ 0 & I \end{array}\right] }$
$\displaystyle{(16)\quad \left[\begin{array}{ccc} A_{CL1}^TX+XA_{CL1} & XB & C_{CL}^T \\ B_{CL}^TX & -\gamma^2 I & D_{CL}^T \\ C_{CL} & D_{CL} &-I \end{array}\right]<0 }$
$\displaystyle{ \Leftrightarrow \left[\begin{array}{cc} A_{CL2}^TX+XA_{CL2} & XB_{CL} \\ B_{CL}^TX & 0 \end{array}\right] < \left[\begin{array}{cc} C_{CL}^T & 0 \\ D_{CL}^T & I \end{array}\right] \left[\begin{array}{cc} -I & 0 \\ 0 & \gamma^2 I \end{array}\right] \left[\begin{array}{cc} C_{CL} & D_{CL} \\ 0 & I \end{array}\right] }$

に注意して

$\displaystyle{(17)\quad \dot{V}(\bar{x}(t))=\frac{d}{dt}(\bar{x}^T(t)X\bar{x}(t))= }$
$\displaystyle{ \left[\begin{array}{c} \bar{x}(t) \\ u_1(t) \end{array}\right]^T \left[\begin{array}{cc} \begin{array}{c} ( p_1(\Omega(t))A_{CL1}+p_2(\Omega(t))A_{CL2} )^TX\\ +X( p_1(\Omega(t))A_{CL1}+p_2(\Omega(t))A_{CL2} ) \end{array} & XB_{CL} \\ B_{CL}^TX & 0 \end{array}\right] \left[\begin{array}{c} \bar{x}(t) \\ u_1(t) \end{array}\right] }$
$\displaystyle{ =p_1(\Omega(t)) \left[\begin{array}{c} \bar{x}(t) \\ u_1(t) \end{array}\right]^T \left[\begin{array}{cc} A_{CL1}^TX+XA_{CL1} & XB_{CL} \\ B_{CL}^TX & 0 \end{array}\right] \left[\begin{array}{c} \bar{x}(t) \\ u_1(t) \end{array}\right] }$
$\displaystyle{ +p_2(\Omega(t)) \left[\begin{array}{c} \bar{x}(t) \\ u_1(t) \end{array}\right]^T \left[\begin{array}{cc} A_{CL2}^TX+XA_{CL2} & XB_{CL} \\ B_{CL}^TX & 0 \end{array}\right] \left[\begin{array}{c} \bar{x}(t) \\ u_1(t) \end{array}\right] }$
$\displaystyle{ < \underbrace{(p_1(\Omega(t))+p_2(\Omega(t)))}_1 \left[\begin{array}{c} \bar{x}(t) \\ u_1(t) \end{array}\right]^T \left[\begin{array}{cc} C_{CL} & D_{CL} \\ 0 & I \end{array}\right]^T \left[\begin{array}{cc} - I & 0 \\ 0 & \gamma^2 I \end{array}\right] \left[\begin{array}{cc} C_{CL} & D_{CL} \\ 0 & I \end{array}\right] \left[\begin{array}{c} \bar{x}(t) \\ u_1(t) \end{array}\right] }$
$\displaystyle{ =\gamma^2 u_1^T(t)u_1(t)-y_1^T(t)y_1(t) }$

を得ます。すなわち時変系としての閉ループ系において

$\displaystyle{(18)\quad \underbrace{\frac{d}{dt}(\bar{x}^TX\bar{x})}_{\dot{V}(\bar{x})} <\underbrace{\gamma^2 u_1^Tu_1-y_1^Ty_1}_{s(u,y)} }$

が成り立ち

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

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

LPV制御（状態FB）

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

LPV制御（状態FB）…Homework

ここでは、主軸変動を伴う回転体の運動制御を例にとって、状態FB形式のGS制御について説明します。

[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} }$

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

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

ただし

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

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

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

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

$\displaystyle{(7)\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)モデルと呼びます。

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

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

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

$\displaystyle{(9)\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)制御を考えます。この閉ループ系は次式で表されます。

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

ここで、端点コントローラによる閉ループ系

$\displaystyle{(11)\quad \dot{x}(t)=(A_1-BF_1)x(t),\ \dot{x}(t)=(A_2-BF_2)x(t) }$

はそれぞれ漸近安定となるように設計しますが、閉ループ系の時変系としての安定性が保証されるかの検討が必要になります（Note B51参照）。

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

$\displaystyle{(12)\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{(13)\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{(14)\quad P_{CL}: \left\{\begin{array}{l} \dot{x}(t)=\underbrace{(p_1(\Omega(t))(A_1-B_2F_1)+p_2(\Omega(t))(A_2-B_2F_2))}_{A_F(\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{(15)\quad \sup_{u_1\in{\cal L}_2}\frac{||y_1(t)||_2}{||u_1(t)||_2}<\gamma }$

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

●そのために、各端点における $H_\infty$ 制御（状態フィードバック）を共通のリャプノフ関数をもつように解きます。その際コントローラの実装上の立場から、各端点において

$\lambda(A_1-BF_1),\lambda(A_2-BF_2)\subset{\cal D}_1\cup{\cal D}_2\cup{\cal D}_3$

の制約を課すことがあります。この問題は次のように定式化されます。

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

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

$\displaystyle{(17)\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{(18)\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{(19)\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{(20)\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{(21)\quad F_1=Z_1Y^{-1}}$
$\displaystyle{(22)\quad F_2=Z_2Y^{-1}}$

●この問題設定では

$\displaystyle{(23)\quad \left[\begin{array}{ccc} Y(A_1-B_2F)^T+(A_1-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 }$
$\displaystyle{(24)\quad \left[\begin{array}{ccc} Y(A_2-B_2F)^T+(A_2-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 }$

の共通解 $Y>0$ を求めています。 $X=Y^{-1}$ ととれば

$\displaystyle{(25)\quad \left[\begin{array}{ccc} (A_1-BF_1)^TX+X(A_1-BF_1) & XB_1 & (C_1-D_{12}F)^T \\ B_1^TX_1 & -\gamma I & D_{11}^T \\ (C_1-D_{12}F) & D_{11} & -\gamma　I \end{array}\right]<0 }$
$\displaystyle{(26)\quad \left[\begin{array}{ccc} (A_2-BF_2)^TX+X(A_2-BF_2) & XB_1 & (C_1-D_{12}F)^T \\ B_1^TX & -\gamma I & D_{11}^T \\ (C_1-D_{12}F) & D_{11} & -\gamma　I \end{array}\right]<0 }$

の共通解となっています。これにより、閉ループ系の時変系としての安定性が保証されます（Note B51参照）。

演習B51…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

図2 回転体の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

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

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

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

Note B51　閉ループ系の時変系としての安定性

●次の時変自由系の漸近安定性を考えます。

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

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

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

を構成して

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

を示すことができれば、平衡状態は漸近安定となることが知られています。

●LPVモデル

$\displaystyle{(4)\quad \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) }$

において、 $A_1-BF_1$ と $A_2-BF_2$ を個別に安定行列としたとすると

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

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

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

が得られますが、閉ループ系(4)のリャプノフ関数を構成したことにはなりません。一方

$\displaystyle{(9)\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{(10)\quad V(x(t))=x^T(t)Xx(t) }$

に対して

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

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

●LPVモデル

$\displaystyle{(12)\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=Cx(t) \end{array}\right. }$

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

$\displaystyle{(13)\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{(14)\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{(15)\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{(16)\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] }$

に注意して

$\displaystyle{(17)\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{(18)\quad \underbrace{\frac{d}{dt}(x^TXx)}_{\dot{V}(x)} <\underbrace{\gamma^2 u^Tu-y^Ty}_{s(u,y)} }$

が成り立ち

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

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

出力FB用LMI（H2制御）

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

出力FB用LMI（H2制御）…Homework

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

$\displaystyle{(1)\quad \begin{array}{l} \displaystyle{{\rm tr} (CW_cC^T)<\gamma^2\quad(W_c=\int_0^\infty \exp(At)BB^T\exp(A^Tt)dt)}\\\\ \displaystyle{\Leftrightarrow \exists X_c>0,\ Q_c>0:\\ \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:\\ \left[\begin{array}{cc} Y_cA^T+AY_c & B \\ B^T & -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}}$

図1 ２ポートシステム

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

$\displaystyle{(2)\quad P: \left\{\begin{array}{l} \dot{x}=Ax+B_1u_1+B_2u_2 \\ y_1=C_1x+D_{11}u_1+D_{12}u_2 \\ y_2=C_2x+D_{21}u_1 \end{array}\right. }$

これに対する出力フィードバック

$\displaystyle{(3)\quad K: \left\{\begin{array}{l} \dot{x}_K=A_Kx_K+B_Ky_2 \\ u_2=C_Kx_K+D_Ky_2 \end{array}\right. }$

による閉ループ系

$\displaystyle{(4) \begin{array}{l} P_{CL}: \\ \left\{\begin{array}{l} \left[\begin{array}{c} \dot{x} \\ \dot{x}_K \end{array}\right] = \underbrace{ \left[\begin{array}{ccc} A+B_2D_KC_2 & B_2C_K \\ B_KC_2 & A_K \end{array}\right] }_{A_{CL}} \left[\begin{array}{c} x \\ x_K \end{array}\right] + \underbrace{ \left[\begin{array}{ccc} B_1+B_2D_KD_{21} \\ B_KD_{21} \end{array}\right] }_{B_{CL}} u_1\\ y_1 = \underbrace{ \left[\begin{array}{ccc} C_1+D_{12}D_KC_2 & D_{12}C_K \end{array}\right] }_{C_{CL}} \left[\begin{array}{c} x \\ x_K \end{array}\right] + \underbrace{ (D_{11}+D_{12}D_KD_{21}) }_{D_{CL}} u_1 \end{array}\right. \end{array} }$

において、その $H_2$ ノルムを考えますが、直達項 $D_{CL}$ が0となるように、 $D_{11}=0$ 、 $D_K=0$ とします。すなわち閉ループ系

$\displaystyle{(5) \begin{array}{l} P_{CL}: \\ \left\{\begin{array}{l} \left[\begin{array}{c} \dot{x} \\ \dot{x}_K \end{array}\right] = \underbrace{ \left[\begin{array}{ccc} A & B_2C_K \\ B_KC_2 & A_K \end{array}\right] }_{A_{CL}} \left[\begin{array}{c} x \\ x_K \end{array}\right] + \underbrace{ \left[\begin{array}{ccc} B_1 \\ B_KD_{21} \end{array}\right] }_{B_{CL}} u_1\\ y_1 = \underbrace{ \left[\begin{array}{ccc} C_1 & D_{12}C_K \end{array}\right] }_{C_{CL}} \left[\begin{array}{c} x \\ x_K \end{array}\right] \end{array}\right. \end{array} }$

の $H_2$ ノルムが $\gamma$ より小、すなわち

$\displaystyle{(6)\quad {\rm tr} (C_{CL}W_cC_{CL}^T)<\gamma^2\quad(W_c=\int_0^\infty \exp(A_{CL}t)B_{CL}B_{CL}^T\exp(A_{CL}^Tt)dt)} }$

となるように出力フィードバックのゲイン $A_K,B_K,C_K$ を求める問題を考えます。

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

$\displaystyle{(7)\quad \begin{array}{l} \exists Y_{CL}>0,\ Q>0:\\ \left[\begin{array}{cc} Y_{CL}A_{CL}^T+A_{CL}Y_{CL} & B_{CL} \\ B_{CL}^T & -I \end{array}\right]<0,\\ \left[\begin{array}{cc} Y_{CL} & Y_{CL}C_{CL}^T \\ C_{CL}Y_{CL} & Q \end{array}\right]>0,\\ {\rm tr}(Q)<\gamma^2 \end{array} }$

出力フィードバックに関する変数変換の議論を参照して、 $Y_{CL}=\Pi_1\Pi_2^{-1}=\Pi_2^{-T}\Pi_1^T$ を代入して

$\displaystyle{(8)\quad \begin{array}{l} \exists \Pi_1\Pi_2^{-1}>0,\ Q>0:\\ \left[\begin{array}{cc} \Pi_2^{-T}\Pi_1^TA_{CL}^T+A_{CL}\Pi_1\Pi_2^{-1} & B_{CL} \\ B_{CL}^T & -I \end{array}\right]<0,\\ \left[\begin{array}{cc} \Pi_2^{-T}\Pi_1^T & \Pi_2^{-T}\Pi_1^TC_{CL}^T \\ C_{CL}\Pi_1\Pi_2^{-1} & Q \end{array}\right]>0,\\ {\rm tr}(Q)<\gamma^2 \end{array} }$

左から $\Pi_2^T$ 、右から $\Pi_2$ をかけると、次のようなLMIとなります。

$\displaystyle{(9)\quad \begin{array}{l} \exists \Pi_2^T\Pi_1>0,\ Q>0:\\ \left[\begin{array}{cc} \Pi_1^TA_{CL}^T\Pi_2+\Pi_2^TA_{CL}\Pi_1 & \Pi_2^TB_{CL} \\ B_{CL}^T\Pi_2 & -I \end{array}\right]<0,\\ \left[\begin{array}{cc} \Pi_1^T\Pi_2 & \Pi_1^TC_{CL}^T \\ C_{CL}\Pi_1 & Q \end{array}\right]>0,\\ {\rm tr}(Q)<\gamma^2 \end{array} }$

すなわち

$\displaystyle{(10)\quad \begin{array}{l} \Pi_2^TA_{CL}\Pi_1= \left[\begin{array}{cc} I & 0 \\ S & N \end{array}\right] \left[\begin{array}{ccc} A & B_2C_K \\ B_KC_2 & A_K \end{array}\right] \left[\begin{array}{cc} R & I \\ M^T & 0 \end{array}\right] \\ =\left[\begin{array}{cc} I & 0 \\ S & N \end{array}\right] \left[\begin{array}{ccc} AR+B_2C_KM^T & A \\ B_KC_2R+A_KM^T & B_KC_2 \end{array}\right] \\ =\left[\begin{array}{ccc} AR+B_2C_KM^T & A\\ SAR+SB_2C_KM^T+NB_KCR+NA_KM^T & SA+NB_KC_2 \end{array}\right] \\ =\left[\begin{array}{cc} AR+B_2{\cal C}_K & A \\ {\cal A}_K & SA+{\cal B}_KC_2 \end{array}\right] \end{array} }$

$\displaystyle{(11)\quad \begin{array}{l} \Pi_2^TB_{CL}= \left[\begin{array}{cc} I & 0 \\ S & N \end{array}\right] \left[\begin{array}{ccc} B_1 \\ B_KD_{21} \end{array}\right] \\ =\left[\begin{array}{ccc} B_1 \\ SB_1+NB_KD_{21} \end{array}\right] =\left[\begin{array}{ccc} B_1 \\ SB_1+{\cal B}_KD_{21} \end{array}\right] \end{array} }$

$\displaystyle{(12)\quad \begin{array}{l} C_{CL}\Pi_1= \left[\begin{array}{ccc} C_1 & D_{12}C_K \end{array}\right] \left[\begin{array}{cc} R & I \\ M^T & 0 \end{array}\right] \\ =\left[\begin{array}{ccc} C_1R+D_{12}C_KM^T & C_1 \end{array}\right] \\ =\left[\begin{array}{ccc} C_1R+D_{12}{\cal C}_K & C_1 \end{array}\right] \end{array} }$

ただし

$\displaystyle{(13)\quad \begin{array}{lll} {\cal A}_K=NA_KM^T+NB_KC_2R+SB_2C_KM^T+SAR \nonumber\\ {\cal B}_K=NB_K \nonumber\\ {\cal C}_K=C_KM^T \nonumber\\ (I-SR=NM^T)\nonumber \end{array} }$

したがって、次の制約条件

$\displaystyle{(14) \begin{array}{l} \left[\begin{array}{cc} R & I \\ I & S \end{array}\right] >0,\ Q>0 \end{array} }$

$\displaystyle{(15) \left[\begin{array}{ccc} \left[\begin{array}{cc} AR+B_2{\cal C}_K & A\\ {\cal A}_K & SA+{\cal B}_KC_2 \end{array}\right]+(*)^T & \left[\begin{array}{c} B_1 \\ SB_1+{\cal B}_KD_{21} \end{array}\right]\\ (*)^T & -I \end{array}\right]<0 }$

$\displaystyle{(16) \left[\begin{array}{cc} \left[\begin{array}{cc} R & I \\ I & S \end{array}\right] & (*)^T \\ \left[\begin{array}{cc} C_1R+D_{12}{\cal C}_K & C_1 \end{array}\right] & Q \end{array}\right]>0 }$

$\displaystyle{(17)\quad {\rm tr}(Q)<\gamma^2 }$

の下で、 $\gamma$ を最小化し、 $R=R^T,S=S^T,{\cal A}_K,{\cal B}_K,{\cal C}_K$ を求め、次式によって出力フィードバックのゲインを決定します。

$\displaystyle{(18)\quad \begin{array}{lll} A_K=N^{-1}({\cal A}_K-SAR-{\cal B}_KC_2R-SB_2{\cal C}_K)M^{-T} \nonumber\\ B_K=N^{-1}{\cal B}_K \nonumber\\ C_K={\cal C}_KM^{-T} \nonumber\\ (I-SR=NM^T) \end{array} }$

演習B44…Flipped Classroom
$1^\circ$ 次のコードを参考にして、 $H_2$ 制御（出力FB）を求める関数を作成せよ。

MATLAB

%of_syn_lmi7.m 
%-----
 A=[0 1;-1 -2*0.01]; B1=[0;1]; B2=B1;
 C1=[1 0;0 0]; D11=[0;0]; D12=[0;1]; 
 C2=[1 0]; D21=0; D22=0; 
 [n,m]=size(B1); p=1; 
%=====
 setlmis([]);
 gam=lmivar(1,[1 0]);  
 Q=lmivar(1,[p 1]);  
 R=lmivar(1,[n 1]); 
 S=lmivar(1,[n 1]); 
 Ak=lmivar(2,[n n]); 
 Bk=lmivar(2,[n p]); 
 Ck=lmivar(2,[m n]); 
 Dk=lmivar(2,[m p]); 
%----- (15)
 lmi1=newlmi;
 lmiterm([lmi1 1 1 R],A,1,'s');       %#1:R*A'+AR 
 lmiterm([lmi1 1 1 Ck],B2,1,'s');     %#1:Ck'*B2'+B2*Ck 
 lmiterm([lmi1 2 1 0],A');            %#1:A' 
 lmiterm([lmi1 2 1 Ak],1,1);          %#1:Ak 
 lmiterm([lmi1 2 2 S],1,A,'s');       %#1:A'*S+S*A 
 lmiterm([lmi1 2 2 Bk],1,C2,'s');     %#1:C2'*Bk'+Bk*C2 
 lmiterm([lmi1 1 3 0],B1);            %#1:B1 
 lmiterm([lmi1 2 3 S],1,B1);          %#1:S*B1 
 lmiterm([lmi1 2 3 Bk],1,D21);        %#1:Bk*D21 
 lmiterm([lmi1 3 3 0],-1);            %#1:-I
%----- (16)
 lmi2=newlmi
 lmiterm([-lmi2 1 1 R],1,1);          %#2:R 
 lmiterm([-lmi2 2 1 0],1);            %#2:I 
 lmiterm([-lmi2 2 2 S],1,1);          %#2:S 
 lmiterm([-lmi2,3,1,R],C1,1);         %#2:C1*R
 lmiterm([-lmi2,3,1,Ck],D12,1);       %#2:D12*Ck 
 lmiterm([-lmi2 3 3 Q],1,1);          %#2:Q  
%----- (17) gam-tr(Q)>0
 CT=zeros(1,p*(p+1)/2); 
 for i=1:p, j=(i+1)*i/2; CT(j)=1; end
 lmi3=newlmi;   
 lmiterm([-lmi3 1 1 gam],1,1);        %#3:gam  
 lmiterm([-lmi3 1 1 Q],-CT,1);        %#3:-tr(Q) 
%======
 LMIs=getlmis; 
 cobj=zeros(1,decnbr(LMIs)); 
 cobj(1)=1; 
 [cost,xopt]=mincx(LMIs,cobj); 
%-----
 gopt=dec2mat(LMIs,xopt,gam)
 R=dec2mat(LMIs,xopt,R); 
 S=dec2mat(LMIs,xopt,S); 
 ak=dec2mat(LMIs,xopt,Ak); 
 bk=dec2mat(LMIs,xopt,Bk); 
 ck=dec2mat(LMIs,xopt,Ck); 
 dk=dec2mat(LMIs,xopt,Dk); 
 [u,sd,v]=svd(eye(size(A,1)-S*R);  
 sd=diag(sqrt(1./diag(sd))); 
 Ni=sd*u'; Mti=v*sd; 
 AK=Ni*(ak-S*A*R-bk*C2*R-S*B2*ck)*Mti; 
 BK=Ni*bk; 
 CK=ck*Mti; 
%-----
 pl=eig(A)
 ACL=[A B2*CK;
      BK*C2 AK];
 plCL=eig(ACL)
 close all,figure(1)
 dregion(0,0,2000,pi/2,2000*[-1,1,-1,1]) 
 plot(real(pl),imag(pl),'x',real(plCL),imag(plCL),'*')
 BCL=[B1;
      BK*D21];
 CCL=[C1 D12*CK];
 DCL=zeros(2,1);
 figure(2) 
 ol=ltisys(A,B1,C1,D11); 
 cl=ltisys(ACL,BCL,CCL,DCL); 
 om=logspace(-1,1,100); 
 splot(cl,'sv',om),hold on 
 splot(ol,'sv',om),grid 
%-----
%eof

出力FB用LMI（H∞制御）

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

出力FB用LMI（H∞制御）…Homework

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

$\displaystyle{(1) \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{(2)\quad \hat{G}(s)=C(sI-A)^{-1}B+D }$

図1 ２ポートシステム

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

$\displaystyle{(3)\quad P: \left\{\begin{array}{l} \dot{x}=Ax+B_1u_1+B_2u_2 \\ y_1=C_1x+D_{11}u_1+D_{12}u_2 \\ y_2=C_2x+D_{21}u_1 \end{array}\right. }$

これに対する出力フィードバック

$\displaystyle{(4)\quad K: \left\{\begin{array}{l} \dot{x}_K=A_Kx_K+B_Ky_2 \\ u_2=C_Kx_K+D_Ky_2 \end{array}\right. }$

による閉ループ系

$\displaystyle{(5) \begin{array}{l} P_{CL}: \\ \left\{\begin{array}{l} \left[\begin{array}{c} \dot{x} \\ \dot{x}_K \end{array}\right] = \underbrace{ \left[\begin{array}{ccc} A+B_2D_KC_2 & B_2C_K \\ B_KC_2 & A_K \end{array}\right] }_{A_{CL}} \left[\begin{array}{c} x \\ x_K \end{array}\right] + \underbrace{ \left[\begin{array}{ccc} B_1+B_2D_KD_{21} \\ B_KD_{21} \end{array}\right] }_{B_{CL}} u_1\\ y_1 = \underbrace{ \left[\begin{array}{ccc} C_1+D_{12}D_KC_2 & D_{12}C_K \end{array}\right] }_{C_{CL}} \left[\begin{array}{c} x \\ x_K \end{array}\right] + \underbrace{ (D_{11}+D_{12}D_KD_{21}) }_{D_{CL}} u_1 \end{array}\right. \end{array} }$

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

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

ただし

$\displaystyle{(7)\quad \hat{P}_{CL}(s)=C_{CL}(sI-A_{CL})^{-1}B_{CL}+D_{CL} }$

となるように出力フィードバックのゲイン $A_K,B_K,C_K,D_K$ を求める問題を考えます。

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

$\displaystyle{(8)\quad \exists Y_{CL}>0:\ \left[\begin{array}{ccc} Y_{CL}A_{CL}^T+A_{CL}Y_{CL} & B_{CL} & Y_{CL}C_{CL}^T \\ B_{CL}^T & -\gamma I & D_{CL}^T \\ C_{CL}Y_{CL} & D_{CL} & -\gamma I \end{array}\right]<0 }$

出力フィードバックに関する変数変換の議論を参照して、 $Y_{CL}=\Pi_1\Pi_2^{-1}=\Pi_2^{-T}\Pi_1^T$ を代入して

$\displaystyle{(9)\quad \exists \Pi_1\Pi_2^{-1}>0:\ \left[\begin{array}{ccc} \Pi_2^{-T}\Pi_1^TA_{CL}^T+A_{CL}\Pi_1\Pi_2^{-1} & B_{CL} & \Pi_2^{-T}\Pi_1^TC_{CL}^T \\ B_{CL}^T & -\gamma I & D_{CL}^T \\ C_{CL}\Pi_1\Pi_2^{-1} & D_{CL} & -\gamma I \end{array}\right]<0 }$

左から $\Pi_2^T$ 、右から $\Pi_2$ をかけると、次のようなLMIとなります。

$\displaystyle{(10)\quad \exists \Pi_2^T\Pi_1>0:\ \left[\begin{array}{ccc} \Pi_1^TA_{CL}^T\Pi_2+\Pi_2^TA_{CL}\Pi_1 & \Pi_2^TB_{CL} & \Pi_1^TC_{CL}^T \\ B_{CL}^T\Pi_2 & -\gamma I & D_{CL}^T \\ C_{CL}\Pi_1 & D_{CL} & -\gamma I \end{array}\right]<0 }$

すなわち

$\displaystyle{(11)\quad \begin{array}{l} \Pi_2^TA_{CL}\Pi_1= \left[\begin{array}{cc} I & 0 \\ S & N \end{array}\right] \left[\begin{array}{ccc} A+B_2D_KC_2 & B_2C_K \\ B_KC_2 & A_K \end{array}\right] \left[\begin{array}{cc} R & I \\ M^T & 0 \end{array}\right] \\ =\left[\begin{array}{cc} I & 0 \\ S & N \end{array}\right] \left[\begin{array}{ccc} (A+B_2D_KC_2)R+B_2C_KM^T & A+B_2D_KC_2 \\ B_KC_2R+A_KM^T & B_KC_2 \end{array}\right] \\ =\left[\begin{array}{ccc} (A+B_2D_KC_2)R+B_2C_KM^T & \\ S(A+B_2D_KC_2)R+SB_2C_KM^T+NB_KCR+NA_KM^T & \end{array}\right. \\ \left.\begin{array}{ccc} & A+B_2D_KC_2 \\ & S(A+B_2D_KC_2)+NB_KC_2 \end{array}\right] \\ =\left[\begin{array}{cc} AR+B_2{\cal C}_K & A+B_2D_KC_2 \\ {\cal A}_K & SA+{\cal B}_KC_2 \end{array}\right] \end{array} }$

$\displaystyle{(12)\quad \begin{array}{l} \Pi_2^TB_{CL}= \left[\begin{array}{cc} I & 0 \\ S & N \end{array}\right] \left[\begin{array}{ccc} B_1+B_2D_KD_{21} \\ B_KD_{21} \end{array}\right] \\ =\left[\begin{array}{ccc} B_1+B_2D_KD_{21} \\ S(B_1+B_2D_KD_{21})+NB_KD_{21} \end{array}\right] =\left[\begin{array}{ccc} B_1+B_2D_KD_{21} \\ SB_1+{\cal B}_KD_{21} \end{array}\right] \end{array} }$

$\displaystyle{(13)\quad \begin{array}{l} C_{CL}\Pi_1= \left[\begin{array}{ccc} C_1+D_{12}D_KC_2 & D_{12}C_K \end{array}\right] \left[\begin{array}{cc} R & I \\ M^T & 0 \end{array}\right] \\ =\left[\begin{array}{ccc} (C_1+D_{12}D_KC_2)R+D_{12}C_KM^T & C_1+D_{12}D_KC_2 \end{array}\right] \\ =\left[\begin{array}{ccc} C_1R+D_{12}{\cal C}_K & C_1+D_{12}D_KC_2 \end{array}\right] \end{array} }$

ただし

$\displaystyle{(14)\quad \begin{array}{lll} {\cal A}_K=NA_KM^T+NB_KC_2R+SB_2C_KM^T+S(A+B_2D_KC_2)R \nonumber\\ {\cal B}_K=NB_K+SB_2D_K \nonumber\\ {\cal C}_K=C_KM^T+D_KC_2R \nonumber\\ (I-SR=NM^T)\nonumber \end{array} }$

したがって、次の制約条件

$\displaystyle{(15) \begin{array}{l} \left[\begin{array}{cc} R & I \\ I & S \end{array}\right] >0 \end{array} }$

$\displaystyle{(16) \left[\begin{array}{ccc} \left[\begin{array}{cc} AR+B_2{\cal C}_K & A+B_2D_KC_2\\ {\cal A}_K & SA+{\cal B}_KC_2 \end{array}\right]+(*)^T & \left[\begin{array}{c} B_1+B_2D_KD_{21} \\ SB_1+{\cal B}_KD_{21} \end{array}\right] & (*)^T \\ (*)^T & -\gamma I & (*)^T \\ \left[\begin{array}{cc} C_1R+D_{12}{\cal C}_K & C_1+D_{12}D_KC_2 \end{array}\right] & D_{11} & -\gamma I \end{array}\right]<0 }$

$\displaystyle{(17) \begin{array}{l} \left[\begin{array}{cc} AR+B_2{\cal C}_K & A+B_2D_KC_2\\ {\cal A}_K & SA+{\cal B}_KC_2 \end{array}\right] +(*)^T %\nonumber\\&& +\alpha \left[\begin{array}{cc} R & I \\ I & S \end{array}\right]<0} \\ \left[\begin{array}{cc} -r \left[\begin{array}{cc} R & I \\ I & S \end{array}\right] & \left[\begin{array}{cc} AR+B_2{\cal C}_K & A+B_2D_KC_2\\ {\cal A}_K & SA+{\cal B}_KC_2 \end{array}\right] \\ (*)^T & -r \left[\begin{array}{cc} R & I \\ I & S \end{array}\right] \end{array}\right] <0 \\ \left[\begin{array}{cc} \sin\theta & \cos\theta \\ -\cos\theta & \sin\theta \end{array}\right]\otimes \left[\begin{array}{cc} AR+B_2{\cal C}_K & A+B_2D_KC_2\\ {\cal A}_K & SA+{\cal B}_KC_2 \end{array}\right]+ (*)^T <0 \end{array} }$

の下で、 $\gamma$ を最小化し、 $R=R^T,S=S^T,{\cal A}_K,{\cal B}_K,{\cal C}_K,D_K$ を求め、次式によって出力フィードバックのゲインを決定します。

$\displaystyle{(18)\quad \begin{array}{lll} A_K=N^{-1}({\cal A}_K-S(A-B_2D_KC_2)R-{\cal B}_KC_2R-SB_2{\cal C}_K)M^{-T} \nonumber\\ B_K=N^{-1}({\cal B}_K-SB_2D_K) \nonumber\\ C_K=({\cal C}_K-D_KC_2R)M^{-T} \nonumber\\ (I-SR=NM^T) \end{array} }$

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

MATLAB

%of_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]; 
 C2=[1 0]; D21=0; D22=0; 
%-----
 setlmis([]);
 gam=lmivar(1,[1 0]);    
 R=lmivar(1,[2 1]); 
 S=lmivar(1,[2 1]); 
 Ak=lmivar(2,[2 2]); 
 Bk=lmivar(2,[2 1]); 
 Ck=lmivar(2,[1 2]); 
 Dk=lmivar(2,[1 1]); 
%-----
 lmiRS=newlmi;
 lmiterm([lmiRS 1 1 R],A,1,'s');       %#1:R*A'+AR 
 lmiterm([lmiRS 1 1 Ck],B2,1,'s');     %#1:Ck'*B2'+B2*Ck 
 lmiterm([lmiRS 2 1 0],A');            %#1:A' 
 lmiterm([lmiRS 2 1 Ak],1,1);          %#1:Ak 
 lmiterm([lmiRS 2 1 -Dk],C2',B2');     %#1:C2'*Dk'*B2' 
 lmiterm([lmiRS 2 2 S],1,A,'s');       %#1:A'*S+S*A 
 lmiterm([lmiRS 2 2 Bk],1,C2,'s');     %#1:C2'*Bk'+Bk*C2 
 lmiterm([lmiRS 1 3 0],B1);            %#1:B1 
 lmiterm([lmiRS 1 3 Dk],B2,D21);       %#1:B2*Dk*D21 
 lmiterm([lmiRS 2 3 S],1,B1);          %#1:S*B1 
 lmiterm([lmiRS 2 3 Bk],1,D21);        %#1:Bk*D21 
 lmiterm([lmiRS 3 3 gam],-1,1);        %#1:-gam 
 lmiterm([lmiRS 4 1 R],C1,1);          %#1:C1*R 
 lmiterm([lmiRS 4 1 Ck],D12,1);        %#1:D12*Ck 
 lmiterm([lmiRS 4 2 0],C1);            %#1:C1 
 lmiterm([lmiRS 4 2 Dk],D12,C2);       %#1:D12*Dk*C2 
 lmiterm([lmiRS 4 3 0],D11);           %#1:D11 
 lmiterm([lmiRS 4 3 Dk],D12,D21);      %#1:D12*Dk*D21 
 lmiterm([lmiRS 4 4 gam],-1,1);        %#1:-gam 
%-----
 posX=-newlmi;                         
 lmiterm([posX 1 1 R],1,1);            %#5:R 
 lmiterm([posX 2 1 0],1);              %#5:I 
 lmiterm([posX 2 2 S],1,1);            %#5:S 
%-----
 lmiDk=-newlmi;                        
 lmiterm([lmiDk 1 1 0],1e2);           %#6:1e2 
 lmiterm([lmiDk 2 2 0],1e2);           %#6:1e2 
 lmiterm([lmiDk 2 1 Dk],1,1);          %#6:Dk 
%-----
 lmig=newlmi;                          
 lmiterm([lmig,1,1,gam],1,1);          %#7:gam 
 lmiterm([lmig,1,1,0],-1e3);           %#7:1e3 
%-----
 LMIs=getlmis; 
 cobj=zeros(1,decnbr(LMIs)); 
 cobj(1)=1; 
 [cost,xopt]=mincx(LMIs,cobj); 
%-----
 gopt=dec2mat(LMIs,xopt,gam)
 R=dec2mat(LMIs,xopt,R); 
 S=dec2mat(LMIs,xopt,S); 
 ak=dec2mat(LMIs,xopt,Ak); 
 bk=dec2mat(LMIs,xopt,Bk); 
 ck=dec2mat(LMIs,xopt,Ck); 
 dk=dec2mat(LMIs,xopt,Dk); 
 [u,sd,v]=svd(eye(size(A,1)-S*R);  
 sd=diag(sqrt(1./diag(sd))); 
 Ni=sd*u'; Mti=v*sd; 
 AK=Ni*(ak-S*(A-B2*dk*C2)*R-bk*C2*R-S*B2*ck)*Mti; 
 BK=Ni*(bk-S*B2*dk); 
 CK=(ck-dk*C2*R)*Mti; 
 DK=dk; 
%-----
 pl=eig(A)
 ACL=[A+B2*DK*C2 B2*CK;
      BK*C2 AK];
 plCL=eig(ACL)
 close all,figure(1)
 dregion(0,0,2000,pi/2,2000*[-1,1,-1,1]) 
 plot(real(pl),imag(pl),'x',real(plCL),imag(plCL),'*')
 BCL=[B1+B2*DK*D21;
      BK*D21];
 CCL=[C1+D12*DK*C2 D12*CK];
 DCL=D11+D12*DK*D21;
 figure(2) 
 ol=ltisys(A,B1,C1,D11); 
 cl=ltisys(ACL,BCL,CCL,DCL); 
 om=logspace(-1,1,100); 
 splot(cl,'sv',om),hold on 
 splot(ol,'sv',om),grid 
%-----
%eof

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

MATLAB

%of_syn_lmi6.m 
%-----
 A=[0 1;-1 -2*0.01]; B1=[0;1]; B2=B1;
 C1=[1 0;0 0]; D11=[0;0]; D12=[0;1]; 
 C2=[1 0]; D21=0; D22=0; 
%-----
 setlmis([]);
 gam=lmivar(1,[1 0]);    
 R=lmivar(1,[2 1]); 
 S=lmivar(1,[2 1]); 
 Ak=lmivar(2,[2 2]); 
 Bk=lmivar(2,[2 1]); 
 Ck=lmivar(2,[1 2]); 
 Dk=lmivar(2,[1 1]); 
%-----
 lmiRS=newlmi;
 lmiterm([lmiRS 1 1 R],A,1,'s');       %#1:R*A'+AR 
 lmiterm([lmiRS 1 1 Ck],B2,1,'s');     %#1:Ck'*B2'+B2*Ck 
 lmiterm([lmiRS 2 1 0],A');            %#1:A' 
 lmiterm([lmiRS 2 1 Ak],1,1);          %#1:Ak 
 lmiterm([lmiRS 2 1 -Dk],C2',B2');     %#1:C2'*Dk'*B2' 
 lmiterm([lmiRS 2 2 S],1,A,'s');       %#1:A'*S+S*A 
 lmiterm([lmiRS 2 2 Bk],1,C2,'s');     %#1:C2'*Bk'+Bk*C2 
 lmiterm([lmiRS 1 3 0],B1);            %#1:B1 
 lmiterm([lmiRS 1 3 Dk],B2,D21);       %#1:B2*Dk*D21 
 lmiterm([lmiRS 2 3 S],1,B1);          %#1:S*B1 
 lmiterm([lmiRS 2 3 Bk],1,D21);        %#1:Bk*D21 
 lmiterm([lmiRS 3 3 gam],-1,1);        %#1:-gam 
 lmiterm([lmiRS 4 1 R],C1,1);          %#1:C1*R 
 lmiterm([lmiRS 4 1 Ck],D12,1);        %#1:D12*Ck 
 lmiterm([lmiRS 4 2 0],C1);            %#1:C1 
 lmiterm([lmiRS 4 2 Dk],D12,C2);       %#1:D12*Dk*C2 
 lmiterm([lmiRS 4 3 0],D11);           %#1:D11 
 lmiterm([lmiRS 4 3 Dk],D12,D21);      %#1:D12*Dk*D21 
 lmiterm([lmiRS 4 4 gam],-1,1);        %#1:-gam 
%-----
 alpha=0.1;
 lmiPL1=newlmi; 
 lmiterm([lmiPL1 1 1 R],A,1,'s');      %#2:R*A'+A*R 
 lmiterm([lmiPL1 1 1 Ck],B2,1,'s');    %#2:Ck'*B2'+B2*Ck 
 lmiterm([lmiPL1 2 1 Ak],1,1);         %#2:Ak 
 lmiterm([lmiPL1 1 2 0],A);            %#2:A 
 lmiterm([lmiPL1 1 2 Dk],B2,C2);       %#2:B2*Dk*C2 
 lmiterm([lmiPL1 2 2 S],1,A,'s');      %#2:A'*S+S*A 
 lmiterm([lmiPL1 2 2 Bk],1,C2,'s');    %#2:C2'*Bk'+Bk*C2 
% 
 lmiterm([lmiPL1 1 1 R],2*alpha,1);    %#2:2*alpha*R
 lmiterm([lmiPL1 2 1 0],2*alpha);      %#2:2*alpha*I 
 lmiterm([lmiPL1 2 2 S],2*alpha,1);    %#2:2*alpha*S 
%-----
 r=2;
 lmiPL2=newlmi; 
 lmiterm([lmiPL2 1 1 R],-r,1);         %#3:-r*R 
 lmiterm([lmiPL2 2 1 0],-r);           %#3:-r*I 
 lmiterm([lmiPL2 2 2 S],-r,1);         %#3:-r*S 
% 
 lmiterm([lmiPL2 1 3 R],A,1);          %#3:A*R
 lmiterm([lmiPL2 1 3 Ck],B2,1);        %#3:B2*Ck 
 lmiterm([lmiPL2 2 3 Ak],1,1);         %#3:Ak 
 lmiterm([lmiPL2 1 4 0],A);            %#3:A 
 lmiterm([lmiPL2 1 4 Dk],B2,C2);       %#3:B2*Dk*C2 
 lmiterm([lmiPL2 2 4 S],1,A);          %#3:S*A 
 lmiterm([lmiPL2 2 4 Bk],1,C2);        %#3:Bk*C2 
% 
 lmiterm([lmiPL2 3 3 R],-r,1);         %#3:-r*R
 lmiterm([lmiPL2 4 3 0],-r);           %#3:-r*I 
 lmiterm([lmiPL2 4 4 S],-r,1);         %#3:-r*S 
%-----
 th=0.4*pi; sth=sin(th); cth=cos(th);
 lmiPL3=newlmi; 
 lmiterm([lmiPL3 1 1 R],sth*A,1,'s');  %#4:sth*(R*A'+A*R) 
 lmiterm([lmiPL3 1 1 Ck],sth*B2,1,'s');%#4:sth*(Ck'*B2'+B2*Ck) 
 lmiterm([lmiPL3 2 1 Ak],sth,1);       %#4:sth*(Ak) 
 lmiterm([lmiPL3 1 2 0],sth*A);        %#4:sth*(A) 
 lmiterm([lmiPL3 1 2 Dk],sth*B2,C2);   %#4:sth*(B2*Dk*C2) 
 lmiterm([lmiPL3 2 2 S],1,sth*A,'s');  %#4:sth*(A'*S+S*A) 
 lmiterm([lmiPL3 2 2 Bk],1,sth*C2,'s');%#4:sth*(C2'*Bk'+Bk*C2) 
% 
%  lmiterm([lmiPL3 1 3 R],1,cth*A');   %#4:cth*(R*A')
%  lmiterm([lmiPL3 1 3 R],-cth*A,1);   %#4:cth*(-A*R) 
%  lmiterm([lmiPL3 1 3 -Ck],cth*B2',1);%#4:cth*(Ck'*B2') 
%  lmiterm([lmiPL3 1 3 Ck],-cth*B2,1); %#4:cth*(-B2*Ck) 
 lmiterm([lmiPL3 1 3 R],cth*A,1);      %#1:cth*(A*R) 
 lmiterm([lmiPL3 1 3 R],1,-cth*A');    %#1:cth*(-R*A')
 lmiterm([lmiPL3 1 3 Ck],cth*B2,1);    %#1:cth*(B*Ck) 
 lmiterm([lmiPL3 1 3 -Ck],-cth*B2',1); %#1:cth*(-Ck'*B') 
 lmiterm([lmiPL3 2 3 Ak],cth,1);       %#4:cth*(Ak) 
 lmiterm([lmiPL3 1 4 -Ak],-cth,1);     %#4:cth*(-Ak') 
 lmiterm([lmiPL3 1 4 0],A);            %#4:cth*(A) 
 lmiterm([lmiPL3 2 3 0],-A');          %#4:cth*(-A') 
 lmiterm([lmiPL3 1 4 Dk],cth*B2,C2);   %#4:cth*(B2*Dk*C2) 
 lmiterm([lmiPL3 2 3 -Dk],-cth*C2',B2');%#4:cth*(-C2'*Dk'*B2') 
 lmiterm([lmiPL3 2 4 S],1,cth*A);      %#4:cth*(S*A) 
 lmiterm([lmiPL3 2 4 S],-cth*A',1);    %#4:cth*(-A'*S) 
 lmiterm([lmiPL3 2 4 Bk],1,cth*C2);    %#4:cth*(Bk*C2) 
 lmiterm([lmiPL3 2 4 -Bk],-cth*C2',1); %#4:cth*(-C2'*Bk') 
% 
 lmiterm([lmiPL3 3 3 R],sth*A,1,'s');  %#4:sth*(R*A'+A*R)
 lmiterm([lmiPL3 3 3 Ck],sth*B2,1,'s');%#4:sth*(Ck'*B2'+B2*Ck) 
 lmiterm([lmiPL3 4 3 Ak],sth,1);       %#4:sth*(Ak) 
 lmiterm([lmiPL3 3 4 0],sth*A);        %#4:sth*(A) 
 lmiterm([lmiPL3 3 4 Dk],sth*B2,C2);   %#4:sth*(B2*Dk*C2) 
 lmiterm([lmiPL3 4 4 S],1,sth*A,'s');  %#4:sth*(A'*S+S*A) 
 lmiterm([lmiPL3 4 4 Bk],1,sth*C2,'s');%#4:sth*(C2'*Bk'+Bk*C2) 
%-----
 posX=-newlmi;
 lmiterm([posX 1 1 R],1,1);            %#5:R 
 lmiterm([posX 2 1 0],1);              %#5:I 
 lmiterm([posX 2 2 S],1,1);            %#5:S 
%-----
 lmiDk=-newlmi;                        
 lmiterm([lmiDk 1 1 0],1e2);           %#6:1e2 
 lmiterm([lmiDk 2 2 0],1e2);           %#6:1e2 
 lmiterm([lmiDk 2 1 Dk],1,1);          %#6:Dk 
%-----
 lmig=newlmi;                          
 lmiterm([lmig,1,1,gam],1,1);          %#7:gam 
 lmiterm([lmig,1,1,0],-1e3);           %#7:1e3 
%-----
 LMIs=getlmis; 
 cobj=zeros(1,decnbr(LMIs)); 
 cobj(1)=1; 
 [cost,xopt]=mincx(LMIs,cobj); 
%-----
 gopt=dec2mat(LMIs,xopt,gam)
 R=dec2mat(LMIs,xopt,R); 
 S=dec2mat(LMIs,xopt,S); 
 ak=dec2mat(LMIs,xopt,Ak); 
 bk=dec2mat(LMIs,xopt,Bk); 
 ck=dec2mat(LMIs,xopt,Ck); 
 dk=dec2mat(LMIs,xopt,Dk); 
 [u,sd,v]=svd(eye(size(A,1)-S*R);  
 sd=diag(sqrt(1./diag(sd))); 
 Ni=sd*u'; Mti=v*sd; 
 AK=Ni*(ak-S*(A-B2*dk*C2)*R-bk*C2*R-S*B2*ck)*Mti; 
 BK=Ni*(bk-S*B2*dk); 
 CK=(ck-dk*C2*R)*Mti; 
 DK=dk; 
%-----
 pl=eig(A)
 ACL=[A+B2*DK*C2 B2*CK;
      BK*C2 AK];
 plCL=eig(ACL)
 close all,figure(1)
 dregion(alpha,0,r,th,r*[-1,1,-1,1]) 
 plot(real(pl),imag(pl),'x',real(plCL),imag(plCL),'*')
 BCL=[B1+B2*DK*D21;
      BK*D21];
 CCL=[C1+D12*DK*C2 D12*CK];
 DCL=D11+D12*DK*D21;
 figure(2) 
 ol=ltisys(A,B1,C1,D11); 
 cl=ltisys(ACL,BCL,CCL,DCL); 
 om=logspace(-1,1,100); 
 splot(cl,'sv',om),hold on 
 splot(ol,'sv',om),grid 
%-----
%eof

出力FB用LMI（固有値制約）

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

出力FB用LMI（固有値制約）…Homework

[0] $n$ 次系

$\displaystyle{(1)\quad \left\{\begin{array}{l} \dot{x}=Ax+Bu\\ y=Cx \end{array}\right. }$

に対する出力フィードバック

$\displaystyle{(2)\quad \left\{\begin{array}{l} \dot{x}_K=A_Kx_K+B_Ky \\ u=C_Kx_K+D_Ky \end{array}\right. }$

による閉ループ系

$\displaystyle{(3)\quad \left[\begin{array}{c} \dot{x} \\ \dot{x}_K \end{array}\right] = \underbrace{ \left[\begin{array}{ccc} A+BD_KC & BC_K \\ B_KC & A_K \end{array}\right] }_{A_{CL}} \left[\begin{array}{c} x \\ x_K \end{array}\right] }$

において、

$\displaystyle{(4)\quad\lambda(A_{CL})\subset{\cal D}_i\quad(i=1,2,3)}$

となるように出力フィードバックのゲイン $A_K,B_K,C_K,D_K$ を求める問題を考えます。

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

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

$\displaystyle{(5) \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_{CL})\subset{\cal D}_1$ となるLMI条件は、次のようになります。

$\displaystyle{(6)\quad \exists Y_{CL}>0: 2\alpha Y_{CL}+A_{CL}Y_{CL}+Y_{CL}A_{CL}^T<0 }$

出力フィードバックに関する変数変換の議論を参照して、 $Y_{CL}=\Pi_1\Pi_2^{-1}=\Pi_2^{-T}\Pi_1^T$ を代入して

$\displaystyle{(7)\quad \exists \Pi_1\Pi_2^{-1}>0: 2\alpha \Pi_1\Pi_2^{-1}+A_{CL}\Pi_1\Pi_2^{-1}+\Pi_2^{-T}\Pi_1^TA_{CL}^T<0 }$

左から $\Pi_2^T$ 、右から $\Pi_2$ をかけると、次のようなLMIとなります。

$\displaystyle{(8)\quad \exists \Pi_2^T\Pi_1>0: 2\alpha \Pi_2^T\Pi_1+\Pi_2^TA_{CL}\Pi_1+(\Pi_2^TA_{CL}\Pi_1)^T<0 }$

すなわち

$\displaystyle{(9)\quad \exists \left[\begin{array}{cc} R & I \\ I & S \end{array}\right]>0: 2\alpha \left[\begin{array}{cc} R & I \\ I & S \end{array}\right]+ }$
$\displaystyle{ \left[\begin{array}{cc} AR+B{\cal C}_K & A+BD_KC \\ {\cal A}_K & SA+{\cal B}_KC \end{array}\right] + \left[\begin{array}{cc} AR+B{\cal C}_K & A+BD_KC \\ {\cal A}_K & SA+{\cal B}_KC \end{array}\right]^T <0 }$

これを解いて $R=R^T,S=S^T,{\cal A}_K,{\cal B}_K,{\cal C}_K,D_K$ を求め、次式によって出力フィードバックのゲインを決定します。

$\displaystyle{(10)\quad \begin{array}{lll} A_K=N^{-1}({\cal A}_K-S(A-BD_KC)R-{\cal B}_KCR-SB{\cal C}_K)M^{-T} \\ B_K=N^{-1}({\cal B}_K-SBD_K) \\ C_K=({\cal C}_K-D_KCR)M^{-T} \\ I-SR=NM^T \end{array} }$

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

$\displaystyle{(11) \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_{CL})\subset{\cal D}_1$ となるLMI条件は、次のようになります。

$\displaystyle{(12)\quad \exists Y_{CL}>0: \left[\begin{array}{cc} -rY_{CL} & A_{CL}Y_{CL} \\ Y_{CL}A_{CL}^T & -rY_{CL} \end{array}\right]<0 }$

出力フィードバックに関する変数変換の議論を参照して、 $Y_{CL}=\Pi_1\Pi_2^{-1}=\Pi_2^{-T}\Pi_1^T$ を代入して

$\displaystyle{(13)\quad \exists \Pi_1\Pi_2^{-1}>0: \left[\begin{array}{cc} -r\Pi_1\Pi_2^{-1} & A_{CL}\Pi_1\Pi_2^{-1} \\ \Pi_2^{-T}\Pi_1^TA_{CL}^T & -r\Pi_1\Pi_2^{-1} \end{array}\right]<0 }$

左から $\Pi_2^T$ 、右から $\Pi_2$ をかけると、次のようなLMIとなります。

$\displaystyle{(14)\quad \exists \Pi_2^T\Pi_1>0: \left[\begin{array}{cc} -r\Pi_2^T\Pi_1 & \Pi_2^TA_{CL}\Pi_1 \\ (\Pi_2^TA_{CL}\Pi_1)^T & -r\Pi_2^T\Pi_1 \end{array}\right]<0 }$

すなわち

$\displaystyle{(15)\quad \exists \left[\begin{array}{cc} R & I \\ I & S \end{array}\right]>0:\ }$
$\displaystyle{ \left[\begin{array}{cc} -r \left[\begin{array}{cc} R & I \\ I & S \end{array}\right] & \left[\begin{array}{cc} AR+B{\cal C}_K & A+BD_KC \\ {\cal A}_K & SA+{\cal B}_KC \end{array}\right] \\ \left[\begin{array}{cc} AR+B{\cal C}_K & A+BD_KC \\ {\cal A}_K & SA+{\cal B}_KC \end{array}\right]^T & -r \left[\begin{array}{cc} R & I \\ I & S \end{array}\right] \end{array}\right] <0 }$

これを解いて $R=R^T,S=S^T,{\cal A}_K,{\cal B}_K,{\cal C}_K,D_K$ を求め、次式によって出力フィードバックのゲインを決定します。

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

$\displaystyle{(16) \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_{CL})\subset{\cal D}_1$ となるLMI条件は、次のようになります。

$\displaystyle{(17)\quad \exists Y_{CL}>0:\ }$
$\displaystyle{ \left[\begin{array}{cc} \sin\theta(A_{CL}Y_{CL}+Y_{CL}A_{CL}^T) & \cos\theta(A_{CL}Y_{CL}-Y_{CL}A_{CL}^T) \\ -\cos\theta(A_{CL}Y_{CL}-Y_{CL}A_{CL}^T) & \sin\theta(A_{CL}Y_{CL}+Y_{CL}A_{CL}^T) \end{array}\right] <0 }$

出力フィードバックに関する変数変換の議論を参照して、 $Y_{CL}=\Pi_1\Pi_2^{-1}=\Pi_2^{-T}\Pi_1^T$ を代入して

$\displaystyle{(18)\quad \exists \Pi_1\Pi_2^{-1}>0:\ }$
$\displaystyle{ \left[\begin{array}{cc} \sin\theta(A_{CL}\Pi_1\Pi_2^{-1}+\Pi_2^{-T}\Pi_1^TA_{CL}^T) & \cos\theta(A_{CL}\Pi_1\Pi_2^{-1}-\Pi_2^{-T}\Pi_1^TA_{CL}^T) \\ -\cos\theta(A_{CL}\Pi_1\Pi_2^{-1}-\Pi_2^{-T}\Pi_1^TA_{CL}^T) & \sin\theta(A_{CL}\Pi_1\Pi_2^{-1}+\Pi_2^{-T}\Pi_1^TA_{CL}^T) \end{array}\right] <0 }$

左から $\Pi_2^T$ 、右から $\Pi_2$ をかけると、次のようなLMIとなります。

$\displaystyle{(19)\quad \exists \Pi_2^T\Pi_1>0:\ }$
$\displaystyle{ \left[\begin{array}{cc} \sin\theta(\Pi_2^TA_{CL}\Pi_1+(\Pi_2^TA_{CL}\Pi_1)^T) & \cos\theta(\Pi_2^TA_{CL}\Pi_1-(\Pi_2^TA_{CL}\Pi_1)^T) \\ -\cos\theta(\Pi_2^TA_{CL}\Pi_1-(\Pi_2^TA_{CL}\Pi_1)^T) & \sin\theta(\Pi_2^TA_{CL}\Pi_1+(\Pi_2^TA_{CL}\Pi_1)^T) \end{array}\right]<0 }$

すなわち

$\displaystyle{(20)\quad \begin{array}{l} \exists \left[\begin{array}{cc} R & I \\ I & S \end{array}\right]>0:\ \\ \left[\begin{array}{cc} \sin\theta & \cos\theta \\ -\cos\theta & \sin\theta \end{array}\right] \otimes \left[\begin{array}{cc} AR+B{\cal C}_K & A+BD_KC \\ {\cal A}_K & SA+{\cal B}_KC \end{array}\right]+ \\ (\left[\begin{array}{cc} \sin\theta & \cos\theta \\ -\cos\theta & \sin\theta \end{array}\right] \otimes \left[\begin{array}{cc} AR+B{\cal C}_K & A+BD_KC \\ {\cal A}_K & SA+{\cal B}_KC \end{array}\right])^T <0 \end{array} }$

これを解いて $R=R^T,S=S^T,{\cal A}_K,{\cal B}_K,{\cal C}_K,D_K$ を求め、次式によって出力フィードバックのゲインを決定します。

演習B42…Flipped Classroom
$1^\circ$ 次のコードを参考にして、 $\lambda(A-BF)\subset{\cal D}$ を達成する出力FBを求める関数を作成せよ。

MATLAB

%of_syn_lmi4.m 
%-----
 clear all, close all
 A=[0 1;-1 -2*0.01]; B=[0;1]; C=[1 0];
%-----
 setlmis([]);
 R=lmivar(1,[2 1]); 
 S=lmivar(1,[2 1]); 
 Ak=lmivar(2,[2 2]); 
 Bk=lmivar(2,[2 1]); 
 Ck=lmivar(2,[1 2]); 
 Dk=lmivar(2,[1 1]); 
%-----
 alpha=0.1;
 lmiPL1=newlmi; 
 lmiterm([lmiPL1 1 1 R],A,1,'s');      %#2:R*A'+A*R 
 lmiterm([lmiPL1 1 1 Ck],B,1,'s');     %#2:Ck'*B'+B*Ck 
 lmiterm([lmiPL1 2 1 Ak],1,1);         %#2:Ak 
 lmiterm([lmiPL1 1 2 0],A);            %#2:A 
 lmiterm([lmiPL1 1 2 Dk],B,C);         %#2:B*Dk*C 
 lmiterm([lmiPL1 2 2 S],1,A,'s');      %#2:A'*S+S*A 
 lmiterm([lmiPL1 2 2 Bk],1,C,'s');     %#2:C'*Bk'+Bk*C 
% 
 lmiterm([lmiPL1 1 1 R],2*alpha,1);    %#2:2*alpha*R
 lmiterm([lmiPL1 2 1 0],2*alpha);      %#2:2*alpha*I 
 lmiterm([lmiPL1 2 2 S],2*alpha,1);    %#2:2*alpha*S 
%-----
 r=5;
 lmiPL2=newlmi; 
 lmiterm([lmiPL2 1 1 R],-r,1);         %#3:-r*R 
 lmiterm([lmiPL2 2 1 0],-r);           %#3:-r*I 
 lmiterm([lmiPL2 2 2 S],-r,1);         %#3:-r*S 
% 
 lmiterm([lmiPL2 1 3 R],A,1);          %#3:A*R
 lmiterm([lmiPL2 1 3 Ck],B,1);         %#3:B*Ck 
 lmiterm([lmiPL2 2 3 Ak],1,1);         %#3:Ak 
 lmiterm([lmiPL2 1 4 0],A);            %#3:A 
 lmiterm([lmiPL2 1 4 Dk],B,C);         %#3:B*Dk*C 
 lmiterm([lmiPL2 2 4 S],1,A);          %#3:S*A 
 lmiterm([lmiPL2 2 4 Bk],1,C);         %#3:Bk*C 
% 
 lmiterm([lmiPL2 3 3 R],-r,1);         %#3:-r*R
 lmiterm([lmiPL2 4 3 0],-r);           %#3:-r*I 
 lmiterm([lmiPL2 4 4 S],-r,1);         %#3:-r*S 
%-----
 th=pi/2; sth=sin(th); cth=cos(th);
 lmiPL3=newlmi; 
 lmiterm([lmiPL3 1 1 R],sth*A,1,'s');  %#4:sth*(R*A'+A*R) 
 lmiterm([lmiPL3 1 1 Ck],sth*B,1,'s'); %#4:sth*(Ck'*B'+B*Ck) 
 lmiterm([lmiPL3 2 1 Ak],sth,1);       %#4:sth*(Ak) 
 lmiterm([lmiPL3 1 2 0],sth*A);        %#4:sth*(A) 
 lmiterm([lmiPL3 1 2 Dk],sth*B,C);     %#4:sth*(B*Dk*C) 
 lmiterm([lmiPL3 2 2 S],1,sth*A,'s');  %#4:sth*(A'*S+S*A) 
 lmiterm([lmiPL3 2 2 Bk],1,sth*C,'s'); %#4:sth*(C'*Bk'+Bk*C) 
% 
%lmiterm([lmiPL3 1 3 R],1,cth*A');     %#1:cth*(R*A')
%lmiterm([lmiPL3 1 3 R],-cth*A,1);     %#1:cth*(-A*R) 
%lmiterm([lmiPL3 1 3 -Ck],cth*B',1);   %#1:cth*(Ck'*B') 
%lmiterm([lmiPL3 1 3 Ck],-cth*B,1);    %#1:cth*(-B*Ck) 
 lmiterm([lmiPL3 1 3 R],cth*A,1);      %#1:cth*(A*R) 
 lmiterm([lmiPL3 1 3 R],1,-cth*A');    %#1:cth*(-R*A')
 lmiterm([lmiPL3 1 3 Ck],cth*B,1);     %#1:cth*(B*Ck) 
 lmiterm([lmiPL3 1 3 -Ck],-cth*B',1);  %#1:cth*(-Ck'*B')  
 lmiterm([lmiPL3 2 3 Ak],cth,1);       %#4:cth*(Ak) 
 lmiterm([lmiPL3 1 4 -Ak],-cth,1);     %#4:cth*(-Ak') 
 lmiterm([lmiPL3 1 4 0],A);            %#4:cth*(A) 
 lmiterm([lmiPL3 2 3 0],-A');          %#4:cth*(-A') 
 lmiterm([lmiPL3 1 4 Dk],cth*B,C);     %#4:cth*(B*Dk*C) 
 lmiterm([lmiPL3 2 3 -Dk],-cth*C',B'); %#4:cth*(-C'*Dk'*B') 
 lmiterm([lmiPL3 2 4 S],1,cth*A);      %#4:cth*(S*A) 
 lmiterm([lmiPL3 2 4 S],-cth*A',1);    %#4:cth*(-A'*S) 
 lmiterm([lmiPL3 2 4 Bk],1,cth*C);     %#4:cth*(Bk*C) 
 lmiterm([lmiPL3 2 4 -Bk],-cth*C',1);  %#4:cth*(-C'*Bk') 
% 
 lmiterm([lmiPL3 3 3 R],sth*A,1,'s');  %#4:sth*(R*A'+A*R)
 lmiterm([lmiPL3 3 3 Ck],sth*B,1,'s'); %#4:sth*(Ck'*B'+B*Ck) 
 lmiterm([lmiPL3 4 3 Ak],sth,1);       %#4:sth*(Ak) 
 lmiterm([lmiPL3 3 4 0],sth*A);        %#4:sth*(A) 
 lmiterm([lmiPL3 3 4 Dk],sth*B,C);     %#4:sth*(B*Dk*C) 
 lmiterm([lmiPL3 4 4 S],1,sth*A,'s');  %#4:sth*(A'*S+S*A) 
 lmiterm([lmiPL3 4 4 Bk],1,sth*C,'s'); %#4:sth*(C'*Bk'+Bk*C) 
%-----
 posX=-newlmi;
 lmiterm([posX 1 1 R],1,1);            %#5:R 
 lmiterm([posX 2 1 0],1);              %#5:I 
 lmiterm([posX 2 2 S],1,1);            %#5:S 
%-----
 LMIs=getlmis; 
 [tmin,xfeas]=feasp(LMIs);   
 R=dec2mat(LMIs,xfeas,R); 
 S=dec2mat(LMIs,xfeas,S); 
 ak=dec2mat(LMIs,xfeas,Ak); 
 bk=dec2mat(LMIs,xfeas,Bk); 
 ck=dec2mat(LMIs,xfeas,Ck); 
 dk=dec2mat(LMIs,xfeas,Dk); 
 [u,sd,v]=svd(eye(size(A,1)-S*R);  
 sd=diag(sqrt(1./diag(sd))); 
 Ni=sd*u'; Mti=v*sd; 
 AK=Ni*(ak-S*(A-B*dk*C)*R-bk*C*R-S*B*ck)*Mti; 
 BK=Ni*(bk-S*B*dk); 
 CK=(ck-dk*C*R)*Mti; 
 DK=dk; 
%-----
 pl=eig(A)
 ACL=[A+B*DK*C B*CK;
      BK*C AK];
 plCL=eig(ACL)
 figure(1)
 dregion(0,0,5,0.4*pi,5*[-1,1,-1,1])  
 plot(real(pl),imag(pl),'x',real(plCL),imag(plCL),'*')
%-----
%eof

出力FB用変数変換

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

出力FB用変数変換…Homework

[1]　 $n$ 次系

$\displaystyle{(1)\quad \left\{\begin{array}{l} \dot{x}=Ax+Bu\\ y=Cx \end{array}\right. }$

に対する出力フィードバック

$\displaystyle{(2)\quad \left\{\begin{array}{l} \dot{x}_K=A_Kx_K+B_Ky \\ u=C_Kx_K+D_Ky \end{array}\right. }$

による閉ループ系

を考えます。

出力フィードバックに関する変数変換として、以下の議論が知られています。

$\displaystyle{(4)\quad \begin{array}{l} \underbrace{ \left[\begin{array}{cc} S & N \\ N^T & S' \end{array}\right] }_{X_{CL}} \underbrace{ \left[\begin{array}{cc} R & M \\ M^T & R' \end{array}\right] }_{Y_{CL}=X^{-1}_{CL}} = \left[\begin{array}{cc} I & 0 \\ 0 & I \end{array}\right] \quad(SR+NM^T=I)\\ \underbrace{ \left[\begin{array}{cc} R & M \\ M^T & R' \end{array}\right] }_{Y_{CL}} \underbrace{ \left[\begin{array}{cc} S & N \\ N^T & S' \end{array}\right] }_{X_{CL}=Y^{-1}_{CL}} = \left[\begin{array}{cc} I & 0 \\ 0 & I \end{array}\right] \quad(RS+MN^T=I) \end{array} }$

が与えられるとき、次の命題が成立します（Note B41参照）。

$\displaystyle{(5)\quad X_{CL}>0 \Leftrightarrow Y_{CL}>0 \Leftrightarrow \left[\begin{array}{cc} R & I \\ I & S \end{array}\right] >0 }$

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

$\displaystyle{(6)\quad X_{CL}= \underbrace{ \left[\begin{array}{cc} I & S \\ 0 & N^T \end{array}\right] }_{\Pi_2} \underbrace{ \left[\begin{array}{cc} R & I \\ M^T & 0 \end{array}\right]^{-1} }_{\Pi_1^{-1}} = \underbrace{ \left[\begin{array}{cc} R & M \\ I & 0 \end{array}\right]^{-1} }_{\Pi_1^{-T}} \underbrace{ \left[\begin{array}{cc} I & 0 \\ S & N \end{array}\right] }_{\Pi_2^T} }$

$\displaystyle{(7)\quad Y_{CL}= \underbrace{ \left[\begin{array}{cc} R & I \\ M^T & 0 \end{array}\right] }_{\Pi_1} \underbrace{ \left[\begin{array}{cc} I & S \\ 0 & N^T \end{array}\right]^{-1} }_{\Pi_2^{-1}} = \underbrace{ \left[\begin{array}{cc} I & 0 \\ S & N \end{array}\right]^{-1} }_{\Pi_2^{-T}} \underbrace{ \left[\begin{array}{cc} R & M \\ I & 0 \end{array}\right] }_{\Pi_1^T} }$

$\displaystyle{(8)\quad \left[\begin{array}{cc} R & I \\ I & S \end{array}\right] = \underbrace{ \left[\begin{array}{cc} R & M \\ I & 0 \end{array}\right] }_{\Pi_1^T} \underbrace{ \left[\begin{array}{cc} I & S \\ 0 & N^T \end{array}\right] }_{\Pi_2} = \underbrace{ \left[\begin{array}{cc} I & 0 \\ S & N \end{array}\right] }_{\Pi_2^T} \underbrace{ \left[\begin{array}{cc} R & I \\ M^T & 0 \end{array}\right] }_{\Pi_1} }$

以下では、次の変換を用います。

$\displaystyle{(9)\quad \begin{array}{l} \Pi_2^TA_{CL}\Pi_1= \left[\begin{array}{cc} I & 0 \\ S & N \end{array}\right] \left[\begin{array}{ccc} A+BD_KC & BC_K \\ B_KC & A_K \end{array}\right] \left[\begin{array}{cc} R & I \\ M^T & 0 \end{array}\right] \\ =\left[\begin{array}{cc} I & 0 \\ S & N \end{array}\right] \left[\begin{array}{ccc} (A+BD_KC)R+BC_KM^T & A+BD_KC \\ B_KCR+A_KM^T & B_KC \end{array}\right] \\ =\left[\begin{array}{ccc} (A+BD_KC)R+BC_KM^T & \\ S(A+BD_KC)R+SBC_KM^T+NB_KCR+NA_KM^T & \end{array}\right. \\ \left.\begin{array}{ccc} & A+BD_KC \\ & S(A+BD_KC)+NB_KC \end{array}\right] \\ =\left[\begin{array}{cc} AR+B{\cal C}_K & A+BD_KC \\ {\cal A}_K & SA+{\cal B}_KC \end{array}\right] \end{array} }$

ただし

$\displaystyle{(10)\quad \begin{array}{lll} {\cal A}_K&=&NA_KM^T+NB_KCR+SBC_KM^T+S(A+BD_KC)R \nonumber\\ {\cal B}_K&=&NB_K+SBD_K \nonumber\\ {\cal C}_K&=&C_KM^T+D_KCR \nonumber \end{array} }$

[2] $\lambda(A)\subset{\rm\bf LHP}$ となる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\\ &\Leftrightarrow& \exists Y>0:\ AY+YA^T<0 \end{array} }$

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

$\displaystyle{(12)\quad \exists Y_{CL}>0: A_{CL}Y_{CL}+Y_{CL}A_{CL}^T<0 }$

$Y_{CL}=\Pi_1\Pi_2^{-1}=\Pi_2^{-T}\Pi_1^T$ を代入して

$\displaystyle{(13)\quad \exists \Pi_1\Pi_2^{-1}>0: A_{CL}\Pi_1\Pi_2^{-1}+\Pi_2^{-T}\Pi_1^TA_{CL}^T<0 }$

左から $\Pi_2^T$ 、右から $\Pi_2$ をかけると、次のようなLMIとなります。

$\displaystyle{(14)\quad \exists \Pi_2^T\Pi_1>0: \Pi_2^TA_{CL}\Pi_1+(\Pi_2^TA_{CL}\Pi_1)^T<0 }$

すなわち

$\displaystyle{(15)\quad \exists \left[\begin{array}{cc} R & I \\ I & S \end{array}\right]>0: }$
$\displaystyle{ \left[\begin{array}{cc} AR+B{\cal C}_K & A+BD_KC \\ {\cal A}_K & SA+{\cal B}_KC \end{array}\right] + \left[\begin{array}{cc} AR+B{\cal C}_K & A+BD_KC \\ {\cal A}_K & SA+{\cal B}_KC \end{array}\right]^T <0 }$

これを解いて $R=R^T,S=S^T,{\cal A}_K,{\cal B}_K,{\cal C}_K,D_K$ を求め、次式によって出力フィードバックのゲインを決定します。

$\displaystyle{(16)\quad \begin{array}{lll} A_K=N^{-1}({\cal A}_K-S(A-BD_KC)R-{\cal B}_KCR-SB{\cal C}_K)M^{-T} \nonumber\\ B_K=N^{-1}({\cal B}_K-SBD_K) \nonumber\\ C_K=({\cal C}_K-D_KCR)M^{-T} \nonumberr \end{array} }$

ただし、 $N^{-1}$ と $M^{-T}$ は次の分解を行って求めます。

$\displaystyle{(17)\quad I-SR=U\Sigma V^T=\underbrace{U\Sigma^{1/2}}_{N}\underbrace{\Sigma^{1/2}V^T}_{M^T} \Rightarrow N^{-1}=\Sigma^{-1/2}U^T,\ M^{-T}=V\Sigma^{-1/2} }$

または、 $I-SR=(I-RS)^T=NM^T$ に注意して

$\displaystyle{(18)\quad \begin{array}{l} I-RS=V\Sigma U^T=\underbrace{V\Sigma^{1/2}}_{M}\underbrace{\Sigma^{1/2}U^T}_{N^T}\\ \Rightarrow M^{-T}=(M^{-1})^T=(\Sigma^{-1/2}V^T)^T=V\Sigma^{-1/2},\\ N^{-1}=(N^{-T}))^T=(U\Sigma^{-1/2})^T=\Sigma^{-1/2}U^T \end{array} }$

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

MATLAB

%of_syn_lmi0.m 
%-----
 clear all, close all
 A=[0 1;0 0]; B=[0;1]; C=[1 0];
%-----
 setlmis([]); 
 [R,xxx,Rdec]=lmivar(1,[2 1]); 
 [S,xxx,Sdec]=lmivar(1,[2 1]); 
 Ak=lmivar(2,[2 2]); 
 Bk=lmivar(2,[2 1]); 
 Ck=lmivar(2,[1 2]); 
 Dk=lmivar(2,[1 1]); 
%-----
 lmiPL=newlmi; 
 lmiterm([lmiPL 1 1 R],A,1,'s');      %#1:R*A'+AR 
 lmiterm([lmiPL 1 1 Ck],B,1,'s');     %#1:Ck'*B'+B*Ck 
 lmiterm([lmiPL 2 1 Ak],1,1);         %#1:Ak 
 lmiterm([lmiPL 1 2 0],A);            %#1:A 
 lmiterm([lmiPL 1 2 Dk],B,C);         %#1:B*Dk*C 
 lmiterm([lmiPL 2 2 S],1,A,'s');      %#1:A'*S+S*A 
 lmiterm([lmiPL 2 2 Bk],1,C,'s');     %#1:C'*Bk'+Bk*C 
%-----
 posX=-newlmi;
 lmiterm([posX 1 1 R],1,1);           %#2:R 
 lmiterm([posX 2 1 0],1);             %#2:I 
 lmiterm([posX 2 2 S],1,1);           %#2:S 
%-----
 LMIs=getlmis; 
 [tmin,xfeas]=feasp(LMIs);   
 R=dec2mat(LMIs,xfeas,R); 
 S=dec2mat(LMIs,xfeas,S); 
 ak=dec2mat(LMIs,xfeas,Ak); 
 bk=dec2mat(LMIs,xfeas,Bk); 
 ck=dec2mat(LMIs,xfeas,Ck); 
 dk=dec2mat(LMIs,xfeas,Dk); 
 [u,sd,v]=svd(eye(size(A,1)-S*R);  
 sd=diag(sqrt(1./diag(sd))); 
 Ni=sd*u'; Mti=v*sd; 
 AK=Ni*(ak-S*(A-B*dk*C)*R-bk*C*R-S*B*ck)*Mti; 
 BK=Ni*(bk-S*B*dk); 
 CK=(ck-dk*C*R)*Mti; 
 DK=dk; 
%-----        
 pl=eig(A)          
 ACL=[A+B*DK*C B*CK;          
      BK*C AK];          
 plCL=eig(ACL)          
 close all,figure(1)          
 dregion(0,0,5,pi/2,5*[-1,1,-1,1])            
 plot(real(pl),imag(pl),'x',real(plCL),imag(plCL),'*')          
%-----
%eof

Note B41　補題

●　 ${\bf Projection\ Lemma}$

$\displaystyle{(1)\quad \exists \Theta:\ \Psi+P^T\Theta^TQ+Q^T\Theta P < 0 \Leftrightarrow \left\{\begin{array}{l} W_P^T\Psi W_P<0 \quad (PW_P=0)\\ W_Q^T\Psi W_Q<0 \quad (QW_Q=0) \end{array}\right. }$

実際

$\displaystyle{(2)\quad \exists U_{PQ}:\ \left[\begin{array}{c} P \\ Q \end{array}\right]U_{PQ}=0 \quad\Rightarrow\quad \left\{\begin{array}{l} \exists U_{P}:\ P\underbrace{ \left[\begin{array}{cc} U_{PQ} & U_P \end{array}\right] }_{W_P}=0\\ \exists U_{Q}:\ Q\underbrace{ \left[\begin{array}{cc} U_{PQ} & U_Q \end{array}\right] }_{W_Q}=0\\ \end{array}\right. }$

を満足する $U_{PQ},U_P,U_Q$ を定めますと、次が成り立ちます。

$\displaystyle{(3)\quad T= \left[\begin{array}{cccc} U_{PQ} & U_P & U_Q & V \end{array}\right] \quad\Rightarrow\quad \left\{\begin{array}{l} PT= \left[\begin{array}{cccc} 0 & 0 & P_3 & P_4 \end{array}\right]\\ QT= \left[\begin{array}{cccc} 0 & Q_2 & 0 & Q_4 \end{array}\right] \end{array}\right. }$

この $T$ を用いて次が成り立ちます。

$\displaystyle{(4)\quad \Psi+P^T\Theta^TQ+Q^T\Theta P < 0 \ \Leftrightarrow\ T^T\Psi T+(PT)^T\Theta^T(QT)+(QT)^T\Theta (PT) < 0 }$
$\displaystyle{ \Leftrightarrow \left[\begin{array}{cccc} \Psi_{11} & \Psi_{12} & \Psi_{13} & \Psi_{14} \\ \Psi_{12}^T & \Psi_{22} & \Psi_{23}& \Psi_{24}\\ \Psi_{13}^T & \Psi_{23}^T & \Psi_{33} & \Psi_{34}\\ \Psi_{14}^T & \Psi_{24}^T & \Psi_{34}^T & \Psi_{44} \end{array}\right] +\left[\begin{array}{cccc} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & \Theta_{11}^T & 0 & \Theta_{12}^T \\ 0 & \Theta_{21}^T & 0 & \Theta_{22}^T \end{array}\right] +\left[\begin{array}{cccc} 0 & 0 & 0 & 0 \\ 0 & 0 & \Theta_{11} & \Theta_{21} \\ 0 & 0 & 0 & 0 \\ 0 & 0 & \Theta_{12} & \Theta_{22} \end{array}\right] < 0 }$
$\displaystyle{ \Leftrightarrow \left[\begin{array}{ccc|c} \Psi_{11} & \Psi_{12} & \Psi_{13} & \Psi_{14} \\ \Psi_{12}^T & \Psi_{22} & \Psi_{23}+\Theta_{11} & \Psi_{24}+\Theta_{21} \\ \Psi_{13}^T & \Psi_{23}^T+\Theta_{11}^T & \Psi_{33} & \Psi_{34}+\Theta_{12} \\ \hline \Psi_{14}^T & \Psi_{24}^T+\Theta_{21}^T & \Psi_{34}^T+\Theta_{12}^T & \Psi_{44}+\Theta_{22}+\Theta_{22}^T \end{array}\right]<0}$
$\displaystyle{\Leftrightarrow \left\{\begin{array}{l} \left[\begin{array}{cccc} \Psi_{11} & \Psi_{12} & \Psi_{13} \\ \Psi_{12}^T & \Psi_{22} & \Psi_{23}+\Theta_{11} \\ \Psi_{13}^T & \Psi_{23}^T+\Theta_{11}^T & \Psi_{33} \end{array}\right] < 0 \\ \quad\\ \Psi_{44}+\Theta_{22}+\Theta_{22}^T-\\ \left[\begin{array}{c} \Psi_{14} \\ \Psi_{24}+\Theta_{21} \\ \Psi_{34}+\Theta_{12} \end{array}\right]^T \left[\begin{array}{cccc} \Psi_{11} & \Psi_{12} & \Psi_{13} \\ \Psi_{12}^T & \Psi_{22} & \Psi_{23}+\Theta_{11} \\ \Psi_{13}^T & \Psi_{23}^T+\Theta_{11}^T & \Psi_{33} \end{array}\right]^{-1} \left[\begin{array}{c} \Psi_{14} \\ \Psi_{24}+\Theta_{21}^T \\ \Psi_{34}+\Theta_{12}^T \end{array}\right] < 0 \end{array}\right. }$

この第２式は、与えられた $\Theta_{11},\Theta_{12},\Theta_{21}$ に対して、 $\Theta_{22}=Q_4^T\Theta P_4$ を適当に選んで、常に満足させることができます。一方、第１式は

$\displaystyle{(5)\quad \begin{array}{l} \left[\begin{array}{ccc} I & 0 & 0 \\ -\Psi_{12}^T\Psi_{11}^{-1} & I & 0 \\ -\Psi_{13}^T\Psi_{11}^{-1} & 0 & I \end{array}\right] \left[\begin{array}{cccc} \Psi_{11} & \Psi_{12} & \Psi_{13} \\ \Psi_{12}^T & \Psi_{22} & \Psi_{23}+\Theta_{11} \\ \Psi_{13}^T & \Psi_{23}^T+\Theta_{11}^T & \Psi_{33} \end{array}\right]\\ \times \left[\begin{array}{ccc} I & -\Psi_{11}^{-1}\Psi_{12} & -\Psi_{11}^{-1}\Psi_{13} \\ 0 & I & 0 \\ 0 & 0 & I \end{array}\right]<0 \end{array} }$
$\displaystyle{ \Leftrightarrow \left[\begin{array}{ccc} \Psi_{11} & 0 & 0 \\ 0 & \Psi_{22}-\Psi_{12}^T\Psi_{11}^{-1}\Psi_{12} & \Theta_{11}+\Psi_{23}-\Psi_{12}^T\Psi_{11}^{-1}\Psi_{13} \\ 0 & \Theta_{11}^T+\Psi_{23}^T-\Psi_{13}^T\Psi_{11}^{-1}\Psi_{12} & \Psi_{33}-\Psi_{13}^T\Psi_{11}^{-1}\Psi_{13} \end{array}\right] < 0 }$

ここで、 $\Theta_{11}=Q_2^T\Theta P_3$ を

$\displaystyle{(6)\quad \Theta_{11}+\Psi_{23}-\Psi_{12}^T\Psi_{11}^{-1}\Psi_{13}=0 }$

を満足するように選ぶと、次の条件を得ます。

$\displaystyle{(7)\quad \underbrace{ \left[\begin{array}{cc} \Psi_{11} & \Psi_{12} \\ \Psi_{12}^T & \Psi_{22} \end{array}\right] }_{W_P^T\Psi W_P} < 0,\ \underbrace{ \left[\begin{array}{cc} \Psi_{11} & \Psi_{13} \\ \Psi_{13}^T & \Psi_{33} \end{array}\right] }_{W_Q^T\Psi W_Q}< 0 }$

●　 ${\bf Completion\ Lemma}$

$\displaystyle{(8)\quad \begin{array}{l} \underbrace{ \left[\begin{array}{cc} S & N \\ X & S' \end{array}\right] }_{W} \underbrace{ \left[\begin{array}{cc} R & M \\ Y & R' \end{array}\right] }_{V=W^{-1}} = \left[\begin{array}{cc} I & 0 \\ 0 & I \end{array}\right] \quad(I=SR+NY=SR-NS'^{-1}XR)\\ \underbrace{ \left[\begin{array}{cc} R & M \\ Y & R' \end{array}\right] }_{V} \underbrace{ \left[\begin{array}{cc} S & N \\ X & S' \end{array}\right] }_{W=V^{-1}} = \left[\begin{array}{cc} I & 0 \\ 0 & I \end{array}\right] \quad(I=RS+MX=RS-MR'^{-1}YS) \end{array} }$

が与えられるとき、次が成り立ちます。

$\displaystyle{(9)\quad W+W^T>0 \Leftrightarrow V+V^T>0 \Leftrightarrow \left\{\begin{array}{l} R+R^T>0\\ S+S^T>0 \end{array}\right. }$

実際

$\displaystyle{(10)\quad \left\{\begin{array}{l} I=SR-NS'^{-1}XR\\ I=RS-MR'^{-1}YS \end{array}\right. \Rightarrow \left\{\begin{array}{l} S-R^{-1}=NS'^{-1}X\\ R-S^{-1}=MR'^{-1}Y \end{array}\right. }$

に注意して

$\displaystyle{(11)\quad S-R^{-1}=L_{\ell}L_r }$

のように分解し、 $N=L_{\ell}S', X=L_r$ と選ぶと

$\displaystyle{(12)\quad 0< \underbrace{ \left[\begin{array}{cc} S & L_{\ell}S' \\ L_r & S' \end{array}\right] }_W + \underbrace{ \left[\begin{array}{cc} S^T & L_r^T \\ S'\,^TL_{\ell}^T & S'\,^T \end{array}\right] }_{W^T} }$
$\displaystyle{ = \underbrace{ \left[\begin{array}{cc} S+S^T & L_r^T \\ L_r & 0 \end{array}\right] }_{\Psi} + \underbrace{ \left[\begin{array}{c} 0 \\ I \end{array}\right] }_{P^T} \underbrace{ S'\,^T }_{\Theta^T} \underbrace{ \left[\begin{array}{cc} L_\ell^T & I \end{array}\right] }_{Q} + \underbrace{ \left[\begin{array}{c} L_{\ell} \\ I \end{array}\right] }_{Q^T} \underbrace{ S' }_{\Theta} \underbrace{ \left[\begin{array}{cc} 0 & I \end{array}\right] }_{P}$
$\Leftrightarrow \left\{\begin{array}{l} \underbrace{ \left[\begin{array}{cc} I & 0 \end{array}\right] }_{W_P^T} \underbrace{ \left[\begin{array}{cc} S+S^T & L_r^T \\ L_r & 0 \end{array}\right] }_{\Psi} \underbrace{ \left[\begin{array}{c} I \\ 0 \end{array}\right] }_{W_P}<0 \quad (\underbrace{ \left[\begin{array}{cc} 0 & I \end{array}\right] }_{P} \underbrace{ \left[\begin{array}{c} I \\ 0 \end{array}\right] }_{W_P}=0)\\ \underbrace{ \left[\begin{array}{cc} I & -L_\ell \end{array}\right] }_{W_Q^T} \underbrace{ \left[\begin{array}{cc} S+S^T & L_r^T \\ L_r & 0 \end{array}\right] }_{\Psi} \underbrace{ \left[\begin{array}{c} I \\ -L_\ell^T \end{array}\right] }_{W_Q}<0 \quad (\underbrace{ \left[\begin{array}{cc} L_\ell^T & I \end{array}\right] }_{Q} \underbrace{ \left[\begin{array}{c} I \\ -L_\ell^T \end{array}\right] }_{W_Q}=0)\\ \end{array}\right. }$
$\Leftrightarrow \left\{\begin{array}{l} S+S^T<0 \\ \underbrace{S-L_\ell L_r}_{R^{-1}}+ \underbrace{S^T-L_r^TL_\ell^T}_{R^{-T}}<0 \Leftrightarrow R+R^T<0 \end{array}\right. }$

一方、

$\displaystyle{(13)\quad R-S^{-1}=L_{\ell}L_r }$

のように分解し、 $M=L_{\ell}R', Y=L_r$ と選ぶと

$\displaystyle{(14)\quad 0< \underbrace{ \left[\begin{array}{cc} R & L_{\ell}R' \\ L_r & R' \end{array}\right] }_V + \underbrace{ \left[\begin{array}{cc} R^T & L_r^T \\ R'\,^TL_{\ell}^T & R'\,^T \end{array}\right] }_{V^T} }$
$\displaystyle{ = \underbrace{ \left[\begin{array}{cc} R+R^T & L_r^T \\ L_r & 0 \end{array}\right] }_{\Psi} + \underbrace{ \left[\begin{array}{c} 0 \\ I \end{array}\right] }_{P^T} \underbrace{ R'\,^T }_{\Theta^T} \underbrace{ \left[\begin{array}{cc} L_\ell^T & I \end{array}\right] }_{Q} + \underbrace{ \left[\begin{array}{c} L_{\ell} \\ I \end{array}\right] }_{Q^T} \underbrace{ R' }_{\Theta} \underbrace{ \left[\begin{array}{cc} 0 & I \end{array}\right] }_{P}$
$\Leftrightarrow \left\{\begin{array}{l} \underbrace{ \left[\begin{array}{cc} I & 0 \end{array}\right] }_{W_P^T} \underbrace{ \left[\begin{array}{cc} R+R^T & L_r^T \\ L_r & 0 \end{array}\right] }_{\Psi} \underbrace{ \left[\begin{array}{c} I \\ 0 \end{array}\right] }_{W_P}<0 \quad (\underbrace{ \left[\begin{array}{cc} 0 & I \end{array}\right] }_{P} \underbrace{ \left[\begin{array}{c} I \\ 0 \end{array}\right] }_{W_P}=0)\\ \underbrace{ \left[\begin{array}{cc} I & -L_\ell \end{array}\right] }_{W_Q^T} \underbrace{ \left[\begin{array}{cc} R+R^T & L_r^T \\ L_r & 0 \end{array}\right] }_{\Psi} \underbrace{ \left[\begin{array}{c} I \\ -L_\ell^T \end{array}\right] }_{W_Q}<0 \quad (\underbrace{ \left[\begin{array}{cc} L_\ell^T & I \end{array}\right] }_{Q} \underbrace{ \left[\begin{array}{c} I \\ -L_\ell^T \end{array}\right] }_{W_Q}=0)\\ \end{array}\right. }$
$\Leftrightarrow \left\{\begin{array}{l} R+R^T<0 \\ \underbrace{S-L_\ell L_r}_{S^{-1}}+ \underbrace{R^T-L_r^TL_\ell^T}_{S^{-T}}<0 \Leftrightarrow S+S^T<0 \end{array}\right. }$

●特に、 $S=S^T,R=R^T,S'=S'\,^T,R'=R'\,^T$ 　かつ

$\displaystyle{(15)\quad \begin{array}{l} \underbrace{ \left[\begin{array}{cc} S & N \\ N^T & S' \end{array}\right] }_W \underbrace{ \left[\begin{array}{cc} R & M \\ M^T & R' \end{array}\right] }_{V=W^{-1}} = \left[\begin{array}{cc} I & 0 \\ 0 & I \end{array}\right] \quad(SR+NM^T=I) \\ \underbrace{ \left[\begin{array}{cc} R & M \\ M^T & R' \end{array}\right] }_{V} \underbrace{ \left[\begin{array}{cc} S & N \\ N^T & S' \end{array}\right] }_{W=V^{-1}} = \left[\begin{array}{cc} I & 0 \\ 0 & I \end{array}\right] \quad(RS+MN^T=I) \end{array} }$

のときは

$\displaystyle{(16)\quad W>0 \quad\Leftrightarrow\quad V>0 \quad\Leftrightarrow\quad R>0,\ S>0 }$

となりますが、 $R$ と $S$ は

$\displaystyle{(17)\quad \left\{\begin{array}{l} S-R^{-1}=NS'^{-1}N^T>0\\ R-S^{-1}=MR'^{-1}M^T>0 \end{array}\right. }$

を考慮して、次の命題を得ます。

$\displaystyle{(15)\quad W>0 \quad\Leftrightarrow\quad V>0 \quad\Leftrightarrow\quad \left[\begin{array}{cc} R & I \\ I & S \end{array}\right]>0 }$

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

$\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} }}$

状態FB用LMI（H2制御）

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

状態FB用LMI（H2制御）…Homework

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

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

$\displaystyle{(2)\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{(3)\quad K: u_2=-Fy_2 }$

による閉ループ系

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

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

$\displaystyle{(5)\quad {\rm tr} (C_{CL}W_cC_{CL}^T)<\gamma^2\quad(W_c=\int_0^\infty \exp(A_{CL}t)B_{CL}B_{CL}^T\exp(A_{CL}^Tt)dt)} }$

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

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

$\displaystyle{(6)\quad \begin{array}{l} \exists Y_{CL}>0,\ Q>0:\\ \left[\begin{array}{cc} Y_{CL}A_{CL}^T+A_{CL}Y_{CL} & B_{CL} \\ B_{CL}^T & -I \end{array}\right]<0,\\ \left[\begin{array}{cc} Y_{CL} & Y_{CL}C_{CL}^T \\ C_{CL}Y_{CL} & Q \end{array}\right]>0,\\ {\rm tr}(Q)<\gamma^2 \end{array} }$

すなわち

$\displaystyle{(7)\quad \begin{array}{l} \displaystyle{\exists Y>0,\ Q>0:\\ \left[\begin{array}{cc} Y(A-B_2F)^T+(A-B_2F)Y & B_1 \\ B_1^T & -I \end{array}\right]<0,\\ \left[\begin{array}{cc} Y & Y(C_1-D_{12}F)^T \\ (C_1-D_{12}F)Y & Q \end{array}\right]>0,\\ {\rm tr}(Q)<\gamma^2} \end{array}}$

ここで、変数変換 $Z_c=FY_c$ を行うと、次のような $Y_c$ と $Z_c$ に関するLMIとなります。

$\displaystyle{(8)\quad \begin{array}{l} \displaystyle{\exists Y>0,\ Q>0:\\ \left[\begin{array}{cc} (AY-B_2Z)^T+AY-B_2Z & B_1 \\ B_1^T & -I \end{array}\right]<0,\\ \left[\begin{array}{cc} Y & (C_1Y-D_{12}Z)^T \\ C_1Y-D_{12}Z & Q \end{array}\right]>0,\\ {\rm tr}(Q)<\gamma^2} \end{array}}$

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

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

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

MATLAB

%sf_syn_lmi7.m 
%----- 
%A=[0 1;-1 -2*0.01]; B1=[0;1]; B2=B1;
%C1=[1 0;0 0]; D11=[0;0]; D12=[0;1]; 
 A=1; B1=1; B2=B1; 
 C1=[1;0]; D11=[0;0]; D12=[0;1];
 [n,m]=size(B2);  p=1;
%----- 
 setlmis([]);
 gam=lmivar(1,[1 0]); 
 Q=lmivar(1,[p 1]); 
 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,0],-1);         %#1:-I 
%----- 
 lmi2=newlmi;
 lmiterm([-lmi2,1,1,Y],1,1);       %#2:Y 
 lmiterm([-lmi2,2,1,Y],C1,1);      %#2:C1*Y
 lmiterm([-lmi2,2,1,Z],-D12,1);    %#2:-D12*Z
 lmiterm([-lmi2,2,2,Q],1,1);       %#2:Q
%----- gam-tr(Q)>0
 CT=zeros(1,p*(p+1)/2); 
 for i=1:p, j=(i+1)*i/2; CT(j)=1; end
 lmi3=newlmi;   
 lmiterm([-lmi3 1 1 gam],1,1);     %#3:gam  
 lmiterm([-lmi3 1 1 Q],-CT,1);     %#3:-tr(Q)   
%----- 
 LMIs=getlmis; 
 cobj=zeros(1,decnbr(LMIs)); 
 cobj(1)=1; 
 [cost,xopt]=mincx(LMIs,cobj); 
 gam=dec2mat(LMIs,xopt,gam)
 Q=dec2mat(LMIs,xopt,Q)   
 Y=dec2mat(LMIs,xopt,Y)  
 Z=dec2mat(LMIs,xopt,Z)  
 F=Z/Y
 [Fopt,plopt]=opt(A,B1,C1(1,:),1,1)
%----- 
 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