LPV制御 | 制御系CAD

LPV制御…Homework

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

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

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

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

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

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

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

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

ただし

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

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

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

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

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

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

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

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

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

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

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

図3 LPV制御の枠組み

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

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

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

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

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

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

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

を構成して

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

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

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

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

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

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

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

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

と選べば

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

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

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

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

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

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

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

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

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

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

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

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

に対して、LPV制御では

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

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

[5] LPVモデル

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

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

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

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

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

に注意して、LPV制御では

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

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

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

が成り立ち

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

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

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

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

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

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

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

による閉ループ系

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

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

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

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

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

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

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

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

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

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

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

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

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

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

図4 LPV制御のイメージ

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

図5 LTI制御のイメージ

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

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

演習B03…Flipped Classroom

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

MATLAB

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

図7 回転体のLPV制御

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

MATLAB

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

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

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

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