線形計画問題

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

これは目的関数

を制約条件

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

図1 線形計画問題の例

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

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

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

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

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

 Iteration   :    Best value of t so far 
 
     1                       -1.203989

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

tmin =
   -1.2040

xf =
    1.2781

yf =
    3.3444

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

 Solver for linear objective minimization under LMI constraints 

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

xs =
    1.4914

ys =
    3.7571

LPV制御…Homework

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

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

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

●そのために、は次式で表されることに着目します。

ただし

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

を、係数によって重み付けて

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

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

を、係数によって重み付けた

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

図3 LPV制御の枠組み

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

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

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

平衡状態の近傍において、リャプノフ関数とよばれる

を構成して

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

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

の漸近安定性のLMI条件は

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

と選べば

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

[4] いま、上のとを個別に安定行列としたとすると

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

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

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

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

に対して、LPV制御では

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

[5] LPVモデル

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

が成り立ちます。いま、共通解が求まっているとしますと

に注意して、LPV制御では

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

が成り立ち

が成り立ちます。これは時不変系の場合にはノルムですが、時変系の場合はゲインとよびます。

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

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

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

による閉ループ系

において、そのゲインがより小、すなわち

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

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

Minimize on subject to

この最小化問題の解 を用いて、次のようにを定めます。

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

図4 LPV制御のイメージ

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

図5 LTI制御のイメージ

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

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

演習B03…Flipped Classroom

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

%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制御

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

%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 回転体の制御

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

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

リャプノフの安定定理

[1] 　非線形システム理論に関する文献：
Hassan K.Khalil: Nonlinear Systems, 3rd ed., Prentice Hall, 2002
から次の定理（Theorem 4.1, p.114）を引用します。

すなわち、(2)と(4)を満足するリャプノフ関数の存在を示して（発見して）、非線形系(1)の漸近安定性を示すことができます。

●振り子の例について、エネルギーの変化と安定性の関係を考えてみます。その運動方程式は，摩擦力を入れた場合

です。いま、運動エネルギーと位置エネルギーの和を

と書くと、(5)を用いて

が成り立ちます。このようにエネルギーが減少することは、振り子の揺れが止まることと対応するので、エネルギーの時間変化率から安定性を判定できそうです。

リャプノフの安定定理は十分条件であり、リャプノフ関数が見つからないからといって安定性や漸近安定性が成り立たないとはかぎりません。また、リャプノフ関数の構築法について一般論はないのですが、上でみたようにエネルギーとの関連が深いです。そこで、振り子の平衡状態

の安定性を調べるために、(6)のエネルギーをリャプノフ関数の候補としてみます。まず、(6)と(7)から、はの安定性を保証するためのリャプノフ関数といえます。

つぎに、(7)から

ですが、これはがは漸近安定性を保証するためのリャプノフ関数とはならないことを示しています。しかし、実際は漸近安定なので、ほかにリャプノフ関数の候補を探せる可能性があります。

例えば、エネルギーの代わりに、つぎのリャプノフ関数の候補を考えます。

この時間微分を計算すると、(5)を用いて

となります。ここで

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

[2] 　リャプノフの安定定理を線形系

に適用して、漸近安定性を保証するリャプノフ関数を求めてみます。その候補として次を考えます。

これを微分して(17)を用いると

となります。したがって、線形行列不等式

が成り立てば、(2)と(4)が満足され、線形系(17)の漸近安定性が成り立つといえます。

●逆に、もし線形系(17)が漸近安定、すなわち

であれば、線形行列不等式(20)を満たすを次のように決めることができます。

ここで、は任意です（ある条件下ではとしてもよい）。実際

となって(20)を満足します。また、任意のに対して

ですが、とすると、が得られ、これはを意味し矛盾です。したがってでなけれならず、を得ます。以上から、線形行列不等式(20)が(17)の漸近安定性の必要十分条件になっています。

●(20)において、、とおくと、漸近安定性の等価な条件として次が得られます。

本資料では同様な多くのLMIが登場しますが、(25)のタイプをX-LMI、(26)のタイプをY-LMIと呼びます。X-LMIは特性解析のために、Y-LMIはコントローラ設計（シンセシス）のために使い分けられます。

演習B02…Flipped Classroom
　安定行列に対して、線形行列不等式(20)を満足する正定行列は、次のプログラムで計算できることを確かめよ。

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

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

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

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

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

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

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

●いま、として

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

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

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

演習B01…Flipped Classroom
図2のグラフを描け。

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

図3　SPIN0.slx

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

このSimulinkブロック線図は

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

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

[0] 3Dアニメーションの対象として、次の単振り子を考えます。

図1 単振り子

この運動方程式は

すなわち

となります。

[1] Simulink 3D Animationを使って、単振り子の3Dアニメーションを作成します。これはたとえば図2のように作成されます。

図2 単振り子アニメーションのためのSimulink図

ここで、下半分は非線形運動方程式(2)をブロック線図で表したものです。ブロックVR Sinkを開くと、図3の単振り子が運動する仮想空間（VR図）が現れます。

図3 単振り子のVR図

設定のタブをクリックすると図4の画面が現れます。

図4 単振り子のVR Sink設定

ここで、VR図を描くX3Dプログラムとして、vrpend.x3dが指定されていることがわかります。また、VR図は２つの要素（pendulum, ground）と２つの視点（view1, view2）から成ります。pendulumのところのrotationにチェックを入れていることに注意してください。これにより、VR Sinkへの入力信号（回転軸を表す３次元ベクトル、回転角を表すスカラー）の受け入れを行うことができます。回転軸は[0;0;1]となっています。

図3において座標軸の設定は図5のようになっています。左右方向がｘ軸、上下方向がy軸、画面の貫通方向がz軸です。回転軸を[0;0;1]で指定するとz軸まわりの回転を意味します。

図5 Simulink 3D Animationにおける座標軸

それでは以下に示されるvrpend.x3dについて詳しく見ていきます。

<?xml version="1.0" encoding="UTF-8"?>
<!DOCTYPE X3D PUBLIC "ISO//Web3D//DTD X3D 3.3//EN" 
"http://www.web3d.org/specifications/x3d-3.3.dtd">
<X3D profile='Immersive' version='3.3' 
xmlns:xsd='http://www.w3.org/2001/XMLSchema-instance' 
xsd:noNamespaceSchemaLocation='http://www.web3d.org/specifications/x3d-3.3.xsd'>
<head>
</head>
<Scene>

<Transform DEF='pendulum'  translation='0 0.6 0' center='0 0.5 0'>
  <Shape>
    <Appearance>
      <Material>
      </Material>
    </Appearance>
    <Box size='0.02 1 0.02'>
    </Box>
  </Shape>
  <Viewpoint 
       DEF = 'view2' 
       description = 'inside view' 
       position = '0 0 0'
       orientation = '0 1 0 1.57'> 
  </Viewpoint>
</Transform>

<Transform DEF='ground' >
  <Shape>
    <Appearance>
      <Material diffuseColor = '1.0 0.0 0.0'>
      </Material>
    </Appearance>
    <Box size='2 0.01 2'>
    </Box>
  </Shape>
</Transform>

<Viewpoint 
       DEF = 'view1' 
       description = 'outside view' 
       position = '0 0 0'
       orientation = '1 0 0 -0.5236'>      
</Viewpoint>

</Scene>
</X3D>

DEF=’pendulum’の部分について
これは振り子を描く部分で、<Box size=’0.02 1 0.02’></Box>によって、断面が２ｃｍ平方、長さが１ｍのボックス（直方体）として描いています。このボックスの重心を平行移動、または重心を通る軸まわりに回転させることができます。画面には原点があって、何も指定しなければここにボックスの重心が合わせられます。しかし、回転軸は振り子の端点とし、これを１１０ｃｍのところに固定したいので、translation=’0 0.6 0′で全体を原点から６０ｃｍ浮かせ、center=’0 0.5 0′で重心を５０ｃｍシフトさせています。x-y平面内で回転面することになります。

図6 translationとcenter

DEF=’ground’の部分について
これは床を描く部分で、<Box size=’2 0.01 2’></Box>によって、２ｍ平方、厚さ１ｃｍのボックス（直方体）として描いています。

DEF = ‘view1’の部分について
これは振り子の運動を外から観察するために視点の設定をしています。orientation=’1 0 0 -0.5236′により、回転面をx軸回りに30度傾けて見ることになります。ちなみにDEF=’pendulum’の部分においても視点が設定されています。これはちょうどブランコに乗って景色を見ていることに相当します。orientation=’0 1 0 1.57′により、視点の方向を回転面内に入れています。

[2] VRエディタvreditを用いたvrpend.x3dの作成法について述べます。

まず振り子の生成を行います。

●次図のようにメニューを辿って、ROOT の下に、Transform を生成します。

●次図のようにEdit Nameを選択して、名前を pendulum に設定します。

●次図のようにメニューを辿って、children の下に、Shape を生成します。

●次図のようにメニューを辿って、Shape の下に、Appearance を生成します。

●次図のようにメニューを辿って、Appearance の下に、Material を生成します。

●次図のようにメニューを辿って、geometry の下に、Box を生成します。

●次図のように、Box サイズを指定します。

●次図のようにメニューを辿って、内部視点 view2(inside-view) を生成します。

同様にして、床の生成と外部視点 view1(outside-view) の生成を行います。

特異値分解(Singular-Value Decomposition)

●をサイズの行列()とします。このとき、サイズの直交行列とサイズの直交行列が存在して

が成り立ちます。ここで、サイズの行列は次を満たします。

●だから、仮定より、は個の正固有値と個の零固有値をもち、互いに直交する固有ベクトルをもちます。そこで、の固有値の正の平方根を、大きい順に、のように表し、対応する固有ベクトルをを満足するようにとることができます。いま、を上のように、また

とおくと、は直交行列となり、つぎが成り立ちます。

第2式の左から、をかけて

また、とおくと、第1式からを得ます。そこで、をが直交行列となるように選ぶと

が成り立ちます。

●行列の特異値分解は

のように与えられることを確かめます。

●のサイズはですが、のサイズはであるので、を計算するととなります（サイズの行列の特異値を手計算で求めるには、とのが同じ非零固有値をもつことから、サイズの小さいほうの固有値計算を行えばよい）。これから、の固有値はで、その正の平方根が特異値で、上のの対角成分の特異値は大きい順に並べる約束ですから

のように、とが決まります。つぎに、については、視察によって

とします。ここで、とが直交しており、の制約があります。上式からが出て、と定まります。

行列のノルム

● からへの写像の「伝達特性」をどう測るかを考えます。これはスカラの場合は正比例の関係ですから、比例定数に相当する量を求める話になります。

サイズの行列の次の特異値分解を考えます。

ここで、で、とは直交行列です。

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

●次元ベクトルのノルムとして、次の３通りが知られています。

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

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

すなわち

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

すなわち行列の最大特異値（行列またはの最大固有値の正の平方根）によって表されます。

●行列のノルムについて次式が成り立ちます。

(8)より

に注意して

●一方、ベクトルの２ノルムについて次式が成り立ちます。

これらは(9)の特別な場合と考えられますが、ここでは直接導出してみます。

いま、に関する２次方程式

を考えると、この実数解は１個または０個となることから、この判別式は零または負でなければならないので

より

すなわち(12)の第２式が得られます。また、第１式は、次式から得られます。

●ちなみに、行列のフロベニウスノルムは

で定義されます。これは

と表せるので、次式が成り立ちます。

行列積のフロベニウスノルムについても

が成り立ちます。

連立１次方程式

次の線形方程式（連立１次方程式）を考えます。

ここで、はフルランク（）とします。

の場合（未知数の数が方程式の数より大きい場合）、(1)はunder-determined（劣決定）と呼ばれ、の特異値分解を代入して次のように書けます。

(2)の解候補として

を考えます。より第１項と第２項は直交することから

(3)においてとして得られるは、を最小化する最小ノルム解と呼ばれます。

の場合（未知数の数が方程式の数より小さい場合）、(1)はover-determined（過決定）と呼ばれ、の特異値分解を代入して次のように書けます。

(5)の解候補として

を考えます。において、

ここで、より第１項と第２項は直交することから

(6)はを最小化する最小２乗解と呼ばれます。

定義

●次正方行列の行列指数関数

●諸性質

とは自明でしょう。

は一般項が次式となるためです。

注意すべきは、で、指数法則は可換な行列に対してのみ成立することです。これは

が

と等しくなるためにはとが可換（）でなければならないからです。

はでとおけば出ます。

は一般項が

となることから出ます。

はの両辺を積分して

となることから出ます。

２次の行列指数関数

●の場合、の実Jordan標準形は次の３つに分類されます。

このときの行列指数関数は次式で与えられます。

は、次式において、とが可換であることから出ます。

は、次式において、とが可換であることから出ます。

３次の行列指数関数

●の場合、の実Jordan標準形は次の４つに分類されます。

このときの行列指数関数は次式で与えられます。

一般の行列指数関数

●一般の場合、の実Jordan標準形は

（）

すなわち、次の２種類のジョルダン細胞のブロック対角行列となります（, , は実数）。

このときの行列指数関数は次式で与えられます。

は、次式において、とが可換であることから出ます。

は、まず次式のようなとの和となります。

ここでクロネッカ積に関する恒等式

を用いると、とが可換、とが可換であれば、とは可換となります（の場合、ＯＫです）。また、次式が成り立つことが知られています。

高次系の漸近安定性

●行列指数関数を用いると、微分方程式

の解は次のように表されます。

ここで

と書けることに注意します（は適当な次実正方行列）。

したがって、任意のに対して、のときとなるための条件は

となります。これはのすべての固有値の実部が負を意味します。

●（）の解のグラフを見ると、の場合は、漸近安定ではないが、発散はしないので、不安定とまではいえないのではないかと思うかもしれません。したがって零の固有値を不安定とみなすのか、安定とみなすか迷うところです。しかし、において、の場合、解はとなって、の第２要素が零でない場合は発散します。したがって、一般には零の固有値は不安定とみなします。

●において、のすべての固有値の実部は負または零の場合を考えます。いま実部が零の固有値の代数的重複度をとし、次が成り立つとします。

これが任意のに対して、のときが発散しないための必要十分条件であることが知られています。

固有値と固有ベクトル

●正方行列の固有値と固有ベクトル: を満足すると
例)

●を満足するが存在するためには が必要。

●正方行列の固有値の集合:

●正方行列の固有ベクトルが１次独立ならば
例)

●実ジョルダン標準形（高橋：微分方程式入門、東京大学出版会、pp.110-112, 1988）
次正方実行列は、正則な次正方実行列により、がつぎの形式となるように変換できます。

ただし、、および

ここで、のサイズは、のサイズはです。また、は実数、指定されていない要素はすべて零です。

●の場合、は次の３つに分類されます。

●の場合、は次の４つに分類されます。

対称行列の場合

●対称行列:
例)

●対称行列の固有値は実数
例)

●対称行列の固有ベクトルは直交
例)

●対称行列が正定行列

●対称行列が準正定行列

行列の記法

●実数の集合:
●実ベクトル:
例)

●実行列:
●サイズの行列:
例1)
例2)

●サイズの零行列:
例)

●次の正方行列:
例)

●次の単位行列:
例)

●次の対角行列
例)

以下では、次の行列とを考えます。

，

●行列の和と差
例)

●行列の積
例)
●一般には、

●行列のスカラ倍
例)

●正方行列の乗
例)

●行列の転置:
例)
●

●行列のクロネッカ積
例)

●正方行列のトレース:
例)

●正方行列の行列式:
例)
●

●正則行列: を満足する正方行列
●特異行列: を満足する正方行列
●正則行列の逆行列: を満足する正則行列
例)
●

行列の階数

●
●

●ベクトルは線形独立（１次独立）:

●ベクトルは線形独立でない（１次独立でない）:

●行列の階数:
列（行）ベクトルのうち線形独立（１次独立）なベクトルの個数
例1) 横長(fat)行列が行フルランク(row full rank)

例2) 縦長(tall)行列が行フルランク(column full rank)

●線形方程式（連立１次方程式）:
例)
●が正則行列のとき、の解はとして一意に定まります。

●命題１：
この対偶をとって
●命題２：

命題１⇒は、次から明らかです。

命題１⇒の逆は、命題２⇒が次のように示されることから出ます。

●シルベスターの不等式(はの列数＝の行数)

●原、千田、佐伯、野波：ロバスト制御のためのベンチマーク問題(I)
-３慣性系に対する位置制御・速度制御
●佐伯、千田、野波、原：ロバスト制御のためのベンチマーク問題(II)
-位置制御系の設計例

[1]　次図のような３慣性系を考えます。

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

ここで、物理定数は次の値が想定されています。

物理定数	最小値	公称値	最大値
	0.0009	0.001	0.0011
	0.0009	0.001	0.0011
	0.001	0.002	0.003
	0.045	0.05	0.055
	0.0009	0.001	0.0011
	0.0014	0.007	0.035
	0.0002	0.001	0.01
	0.0009	0.001	0.0011
	828	920	1012
	72	80	88

●以下では、の３つのパラメータ誤差を考慮して、回転体３に対してインパルス外乱が加わるときのシミュレーションを行ってみます。

図1

原論文にはさまざまな問題設定がなされていますが、ここでは、回転体３に対するインパルス外乱の下で、回転体３を速やかに整定させる制御系設計を検討します。ただし、操作入力には次の振幅制限

が課され、観測出力は回転体１の変位とします。

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

ここで、、に対する次の内分式に注目します。

のとき上の状態方程式は、端点モデル

を、によって重み付けして、LPVモデル

として表されます。

[3]　状態FB

図2

[3]　出力FB

図3

演習13.3…Flipped Classroom

図1

図2

図3