制御系CAD

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

応用例

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

応用例…Homework

[1]　ある航空機の線形状態方程式として、次を考えます。

$\displaystyle{(1)\quad \begin{array}{l} \dot{x}(t)=Ax(t)+Bu(t)\\ x(t)= \left[\begin{array}{ll} \theta & pitch\ angle\\ q & pitch\ rate\\ \alpha & angle\ of\ atack\\ \eta & elevator\ deflection\\ \delta & flap\ deflection \end{array}\right],\quad u(t)= \left[\begin{array}{ll} \eta_c & elevator\ command\\ \delta_c& flap\ comannd \end{array}\right]\\ A=\left[\begin{array}{rrrrr} 0 & 1 & 0 & 0 & 0\\ 0 &-1.99 & -13.41 & -18.95 & -3.60\\ 0 & 1 & -1.74 & -0.08 & -0.59\\ 0 & 0 & 0 & -20 & 0\\ 0 & 0 & 0 & 0 & -20 \end{array}\right],\quad B=\left[\begin{array}{rr} 0 & 0 \\ 0 & 0 \\ 0 & 0 \\ 20 & 0\\ 0 & 20 \end{array}\right] \end{array} }$

$A$ 行列の固有値は次のように求められます。

$\displaystyle{(2)\quad \lambda(A):\left\{\begin{array}{cl} -1.864\pm j 3.660 & short\ period\ mode\\ 0 & pitch\ attitude\ mode\\ -20 & elevator\ actuator\ mode\\ -20 & flap\ actuator\ mode \end{array}\right. }$

状態変数として、pitch angle $\theta$ の代わりに、fligt path angle $\gamma=\theta-\alpha$ を用いると、次式となります。

$\displaystyle{(3)\quad \begin{array}{l} \dot{x}(t)=Ax(t)+Bu(t)\\ x(t)= \left[\begin{array}{ll} \gamma=\theta-\alpha & fligt\ path\ angle\\ q & pitch\ rate\\ \alpha & angle\ of\ atack\\ \eta & elevator\ deflection\\ \delta & flap\ deflection \end{array}\right],\quad u(t)= \left[\begin{array}{ll} \eta_c & elevator\ command\\ \delta_c& flap\ comannd \end{array}\right]\\ A=\left[\begin{array}{rrrrr} 0 & 0 & 1.74 & 0.08 & 0.59\\ 0 &-1.99 & -13.41 & -18.95 & -3.60\\ 0 & 1 & -1.74 & -0.08 & -0.59\\ 0 & 0 & 0 & -20 & 0\\ 0 & 0 & 0 & 0 & -20 \end{array}\right],\quad B=\left[\begin{array}{rr} 0 & 0 \\ 0 & 0 \\ 0 & 0 \\ 20 & 0\\ 0 & 20 \end{array}\right] \end{array} }$

また、評価変数として、pitch angle $\theta$ とfligt path angle $\gamma$ をとると、その出力方程式は次式のように与えられます。

$\displaystyle{(4)\quad \begin{array}{l} y(t)=Hx(t)\\ H= \left[\begin{array}{ccccc} 1 & 0 & 1 & 0 & 0\\ 1 & 0 & 0 & 0 & 0 \end{array}\right] \end{array} }$

制御目的は、次図のように、pitch angle $\theta$ とfligt path angle $\gamma$ の間の非干渉化を達成し、独立して設定値に合せることです。

図1 航空機の操縦法

[2] モデル規範による追従制御

この場合の制御則は次のように構成されます。

$\displaystyle{(5)\quad { \boxed{\begin{array}{l} u(t)=-Fx(t)+Gr(t)+\underbrace{u_\ell(t)}_{linear\ control}+\underbrace{u_n(t)}_{switching\ component}\\ u_\ell(t)=-\underbrace{(SB)^{-1}(SA_m-\Phi S)}_{L=L_{eq}+L_{phi}}x(t) \quad(\Phi:stable\ matrix)\\ u_n(t)=-\underbrace{(SB)^{-1}\rho(t,x)}_{L_n}\frac{P_2s(t)}{||P_2s(t)||}\quad(P_2\Phi+\Phi^TP_2=-I) \end{array}}} }$

●追従すべき規範モデルは

$\displaystyle{(6)\quad \begin{array}{l} \dot{x}(t)=A_mx(t)+B_mu(t)\\ y(t)=Hx(t) \end{array}} }$

ただし

$\displaystyle{(7)\quad \begin{array}{l} A_m=A-BF\\ B_m=BG\\ G=-(H(A-BF)^{-1}B)^{-1} \end{array} }$

いまこの規範モデルにもたせるべき望ましい固有値を

$\displaystyle{(8)\quad \lambda(A_m):\left\{\begin{array}{cl} -5.6\pm j 4.2 & short\ period\ mode\\ -1 & pitch\ attitude\ mode\\ -20 & elevator\ actuator\ mode\\ -20 & flap\ actuator\ mode \end{array}\right. }$

とし、また望ましいモード分布を次により指定します。

$\displaystyle{(9)\quad {\tt specpos}= \left[\begin{array}{ccc|cc} 1 & 1 & 1 & 0 & 0\\ 1 & 0 & 1 & 0 & 0\\ 0 & 1 & 1 & 0 & 0\\\hline 0 & 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 0 & 1\\ \end{array}\right],\quad {\tt specent}= \left[\begin{array}{rrr|rr} 0 & 0 & 1 & -1 & -1\\ 1 & -1 & 0 & -1 & -1\\ -1 & 1 & 0 & -1 & -1\\\hline -1 & -1 & -1 & 1 & -1\\ -1 & -1 & -1 & -1 & 1\\ \end{array}\right] }$

ここで、 ${\tt specpos}$ は望ましい固有ベクトルの要素を指定する場合１、指定しない場合０で表した行列、 ${\tt specエent}$ は望ましい固有ベクトルの要素が非零指定の場合１、零の場合０、非零任意の場合－１で表した行列です。

このとき以下に示すプログラムを使って、次が得られます。

$\displaystyle{(10)\quad \begin{array}{l} F=WV^{-1} =\left[\begin{array}{rrrrr} -1.9338 & -0.5744 & -2.1188 & 0.4750 & 0.1066\\ 1.9571 & 0.2253 & 3.1273 & -0.0694 & -0.0515 \end{array}\right]\\ W=\left[\begin{array}{rrrrr} 0.8445 & 1.9935 & -0.2742 & 0 & 0\\ -2.7055 & 1.2780 & 0.1767 & 0 & 0 \end{array}\right]\\ V=\left[\begin{array}{rrrrr} 0 & -0.0000 & 0.6667 & 0.0030 & -0.0313\\ 1.0000 & -9.5000 & -0.3333 & 0.9964 & 0.1493\\ -0.9286 & 1.0000 & -0.3333 & -0.0528 & 0.0238\\ -1.8252 & -2.2364 & 0.2886 & 1.0000 & -0.0649\\ 2.9860 & -2.6459 & -0.1860 & -0.0825 & 1.0000 \end{array}\right] \end{array} }$

●ここでは、スイッチング関数

$\displaystyle{(11)\quad s(t)= \underbrace{ \left[\begin{array}{cc} M & I \end{array}\right] }_{S} x(t) }$

を構成するために

$\displaystyle{(12)\quad \begin{array}{l} \dot{x}_1(t)=A_{11}x_1(t)+A_{12}\hat{u}(t) \end{array} }$

を２次形式評価関数

$\displaystyle{(13)\quad J=\int_0^\infty(x_1^T(t)Q_{11}x_1(t)+\hat{u}^T(t)Q_{22}\hat{u}(t))dt }$

を最小にする

$\displaystyle{(14)\quad \hat{u}(t)=-Mx_1(t) }$

を求めます。

$\displaystyle{(15)\quad Q_{11}={\rm diag}\{10,5,5\},\ Q_{22}={\rm diag}\{20,20\} }$

と選ぶとき、スイッチング関数を決める行列 $S$ が、次のように得られます。

$\displaystyle{(16)\quad S=\left[\begin{array}{rrr|rr} 0.7025 & 0.4209 & 0.1894 & -1 & 0\\ -0.0810 & 0.0635 & 0.0267 & 0 & -1\\ \end{array}\right] }$

● $P_2\Phi+\Phi^TP_2=-I$ を満たす $\Phi$ と $P_2$ は、たとえば次のように定めます。

$\displaystyle{(17)\quad \underbrace{ \left[\begin{array}{rr} 0.025 & 0\\ 0 & 0.025 \end{array}\right] }_{P_2} \underbrace{ \left[\begin{array}{rr} 20 & 0\\ 0 & 20 \end{array}\right] }_{\Phi} +(*)^T=-I }$

このときSMC則における行列 $L$ は次のように計算されます。

$\displaystyle{(18)\quad L=\left[\begin{array}{rrrrr} -1.2285 & -0.1858 & -2.1624 & 0.0782 & 0.0460\\ 1.8631 & 0.2827 & 3.0804 & -0.1299 & -0.0060 \end{array}\right] }$

図2 モデル規範法によるSMC制御系応答

$\displaystyle{{ \boxed{\begin{array}{l} u(t)=-Fx(t)+Gr(t)+\underbrace{u_\ell(t)}_{linear\ control}+\underbrace{u_n(t)}_{switching\ component}\\ u_\ell(t)=-\underbrace{(SB)^{-1}(SA_m-\Phi S)}_{L=L_{eq}+L_{phi}}x(t) \quad(\Phi:stable\ matrix)\\ u_n(t)=-\underbrace{(SB)^{-1}\rho(t,x)}_{L_n}\frac{P_2s(t)}{||P_2s(t)||}\quad(P_2\Phi+\Phi^TP_2=-I) \end{array}}} }$

MATLAB

%cAIRCRAFT_smcm.m
%-----
 clear all, close all
 A= [0 0 1.74 0.08 0.59;
     0 -1.99 -13.41 -18.95 -3.6;
     0 1 -1.74 -0.08 -0.59;
     0 0 0 -20 0;
     0 0 0 0 -20];
 B=[0 0;0 0;0 0;20 0;0 20];
%----- 
 lambda=[-5.6 4.2 -1 -20 -20];
 nocomp=1;
%固有ベクトルの非零要素の位置 (１：指定、０：指定せず)
 specpos=[ 1  1  1  0  0; 
           1  0  1  0  0; 
           0  1  1  0  0;
           0  0  0  1  0;
           0  0  0  0  1]
%固有ベクトルの要素の値 (０：零、１：非零指定、－１：非零任意)       
 specent=[ 0  0  1 -1 -1;
           1 -1  0 -1 -1;
          -1  1  0 -1 -1;
          -1 -1 -1  1 -1;
          -1 -1 -1 -1  1]
 F=vplace2(A,B,lambda,nocomp,specpos,specent)       
 pl=eig(A-B*F)
%-----
 H=[1 0 1 0 0;1 0 0 0 0];
 Am=A-B*F;
 G=-inv(H*(Am\B));
 Bm=B*G;
%-----
 Q=diag([10,5,5,20,20]);
 S=swflqr(A,B,Q)
%-----
 P2=diag([0.025,0.025]);
 Phi=-20/0.025*P2;
 Check=P2*Phi+Phi*P2
 L=inv(S*B)*(S*Am-Phi*S)
 Ln=inv(S*B)
%-----
 sim('AIRCRAFT_smcm_2015a.mdl')
%-----
%eof

モデル規範による追従SM制御

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

モデル規範による追従SM制御…Homework

[1]　次の状態方程式で表される制御対象を考えます。

$\displaystyle{(1)\quad \dot{x}(t)=Ax(t)+Bu(t) }$

この状態が理想的なモデル

$\displaystyle{(2)\quad \dot{x}_m(t)=A_mx_m(t)+B_mr(t) }$

の状態を追従するように、すなわち

$\displaystyle{(3)\quad e(t)=x(t)-x_m(t)\rightarrow 0 \quad(t\rightarrow\infty) }$

となるように制御則を決定したいとします。

●追従誤差のダイナミックスは

$\displaystyle{(4)\quad \begin{array}{l} %\bullet \dot{e}(t)=Ax(t)+Bu(t)-(A_mx_m(t)+B_mr(t))\\ %=Ax(t)-Ax_m(t)+Ax_m(t)-A_mx_m(t)+Bu(t)-B_mr(t)\\ %=Ae(t)+(A-A_m)x_m(t)+Bu(t)-B_mr(t)\\ %\bullet \dot{e}(t)=Ax(t)+Bu(t)-(A_mx_m(t)+B_mr(t))\\ =Ax(t)-A_mx(t)+A_mx(t)-A_mx_m(t)+Bu(t)-B_mr(t)\\ =A_me(t)+(A-A_m)x(t)+Bu(t)-B_mr(t) \end{array} }$

となります。いま条件

$\displaystyle{(5)\quad \begin{array}{l} (BB^\dag-I)(A-A_m)=0\Leftrightarrow B\underbrace{B^\dag(A-A_m)}_{F}=A-A_m\\ (BB^\dag-I)B_m=0\Leftrightarrow B\underbrace{B^\dag B_m}_{G}=B_m \end{array} }$

を仮定すると（ $B^\dag$ は $BB^\dag B=B$ を満たす疑似逆行列）、制御則

$\displaystyle{(6)\quad u(t)=-Ke(t)-Fx(t)+Gr(t) }$

の下では、誤差方程式は次式となります。

$\displaystyle{(7)\quad \dot{e}(t)=(A_m-BK)e(t) }$

実際、

$\displaystyle{(8)\quad \begin{array}{l} \dot{e}(t)=A_me(t)+(A-A_m)x(t)+Bu-B_mr(t)\\ =A_me(t)+(A-A_m)x(t)+B(-Ke(t)-Fx(t)+Gr(t))-B_mr(t)\\ =A_me(t)+(A-A_m)x(t)-BKe(t)-BB^\dag(A-A_m)x(t)\\ +BB^\dag B_mr(t)-B_mr(t)\\ =(A_m-BK)e(t)-\underbrace{(BB^\dag-I)(A-A_m)}_{0}x(t)+\underbrace{(BB^\dag -I)B_m}_{0}r(t)\\ \end{array} }$

●以上から、制御則(6)の設計の指針として、望ましいモデル(2)を(1)に対する状態フィードバックと入力変換で得るものとします。すなわち

$\displaystyle{(9)\quad \boxed{\begin{array}{l} A_m=A-BF\\ B_m=BG \end{array}} }$

ここで $G$ は、 $r$ から出力

$\displaystyle{(10)\quad y(t)=Hx(t) }$

までの定常ゲイン

$\displaystyle{(11)\quad %\begin{array}{l} \left\{\begin{array}{l} 0=(A-BF)x(\infty)+BGr(\infty)\\ y(\infty)=Hx(\infty) \end{array}\right. \Rightarrow y(\infty)=\underbrace{-H(A-BF)^{-1}BG}_{I}r(\infty) %\end{array} }$

が単位行列となるように、次式で決めます。

$\displaystyle{(12)\quad \boxed{G=-(H(A-BF)^{-1}B)^{-1}} }$

ちなみに $K$ については、(7)が漸近安定であることが前提となり、望ましいモデル(2)より、追従誤差は速いダイナミックスをもつことが必要であることに注意します。

[2]　誤差方程式(4)すなわち

$\displaystyle{(13)\quad \dot{e}(t)=A_me(t)+(A-A_m)x(t)+Bu(t)-B_mr(t) }$

に対して、スライディングモード制御を適用することを考えます。スイッチング関数を

$\displaystyle{(14)\quad s(t)=Se(t) }$

とするとき、スライディングモード時は

$\displaystyle{(15)\quad Se(t)=0\Rightarrow S\dot{e}(t)=0 }$

すなわち

$\displaystyle{(16)\quad S\dot{e}(t)=S(A_me(t)+(A-A_m)x(t)+Bu(t)-B_mr(t))=0 }$

が成り立ち、等価制御は次式となります。

$\displaystyle{(17)\quad u_{eq}(t)=-(SB)^{-1}S(A_me(t)+(A-A_m)x(t)-B_mr(t)) }$

このとき追従誤差のダイナミックスは次式となります。

$\displaystyle{(18)\quad \begin{array}{l} \dot{e}(t)=A_me(t)+(A-A_m)x(t)+Bu(t)-B_mr(t)\\ =A_me(t)+(A-A_m)x(t)\\ +B(-(SB)^{-1}S(A_me(t)+(A-A_m)x(t)-B_mr(t)))-B_mr(t)\\ =(I-B(SB)^{-1}S)(A_me(t)+(A-A_m)x(t)-B_mr(t))\\ =(I-B(SB)^{-1}S)(A_me(t)+BFx(t)-BGr(t))\\ =(I-B(SB)^{-1}S)A_me(t)+\underbrace{(I-B(SB)^{-1}S)B}_{0}(Fx(t)-Gr(t))\\ =\underbrace{(I-B(SB)^{-1}S)A_m}_{A_{eq}}e(t) \end{array} }$

したがって、スライディングモード制御を適用する場合、全体の制御則は、次式で表されます。

$\displaystyle{(19)\quad { \boxed{\begin{array}{l} u(t)=-Fx(t)+Gr(t)+\underbrace{u_\ell(t)}_{linear\ control}+\underbrace{u_n(t)}_{switching\ component}\\ u_\ell(t)=-\underbrace{(SB)^{-1}(SA_m-\Phi S)}_{L=L_{eq}+L_{phi}}x(t) \quad(\Phi:stable\ matrix)\\ \displaystyle{u_n(t)=-\underbrace{(SB)^{-1}\rho(t,x)}_{L_n}\frac{P_2s(t)}{||P_2s(t)||}\quad(P_2\Phi+\Phi^TP_2=-I)} \end{array}}} }$

演習C51…Flipped Classroom

MATLAB

%cCIP_smcm.m
%-----
 clear all, close all
 global mc m ell g th0
 mc=1; m=0.1; ell=0.2; g=9.8;
 ths=0; %input('ths   = <0,180> ')/180*pi;
 th0=3/180*pi; %input('th(0) = <0,180> ')/180*pi;
%-----
 A=[zeros(2,2) eye(2);zeros(2,4)];
 A(3,1)=0; 
 A(3,2)=-3*m*g/(m+4*mc); 
 A(4,1)=0; 
 A(4,2)=3*(m+mc)*g/((m+4*mc)*ell);  
 B=zeros(4,1);
 B(3)=4/(m+4*mc);
 B(4)=-3/((m+4*mc)*ell); 
 H=[1 0 0 0];
%-----
 lambda=[-0.8*1.5 0.5 -1.5*1.5 1];
 nocomp=2;
 specpos=rand(4,4);
 specent=rand(4,4);
 F=vplace2(A,B,lambda,nocomp,specpos,specent)       
 pl=eig(A-B*F)
%-----
 H=[1 0 0 0];
 Am=A-B*F;
 G=-inv(H*(Am\B));
 Bm=B*G;
%-----
 Mr=0.5; Mth=3/180*pi; 
 Tc=1; Tpen=0.25*2*pi*sqrt(2*ell/g);
 Q=diag([1/Mr^2,1/Mth^2,Tc^2/Mr^2,Tpen^2/Mth^2]);
 S=swflqr(A,B,Q) 
%-----
 Phi=-0.1;
 P2=0.5*inv(-Phi);
 Check=P2*Phi+Phi*P2
 P2S=P2*S;
 Leq=inv(S*B)*S*A;
 LPhi=-inv(S*B)*Phi*S;
 L=Leq+LPhi;
 rho=1;
 Ln=inv(S*B)*rho;
%-----
 x0=[0;th0*0;0;0];
 sim('CIP_smcm_2015a.mdl')
%-----
%eof

固有ベクトル設定問題

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

[0] 次の２入力２次系

$\displaystyle{(1)\quad \dot{x}= \underbrace{ \left[\begin{array}{cc} 0 & 0\\ 0 & -1 \end{array}\right] }_{A} x+ \underbrace{ \left[\begin{array}{cc} 1 & 1\\ 1 & -1 \end{array}\right] }_{B} u }$

に対して２つの異なる状態フィードバック

$\displaystyle{(2)\quad %\begin{array}{l} u=- \underbrace{ \left[\begin{array}{cc} 1 & 1\\ 1 & -1 \end{array}\right] }_{F_1} x,\quad u=- \underbrace{ \left[\begin{array}{cc} 1 & 0\\ 1 & 2 \end{array}\right] }_{F_2} x %\end{array} }$

による閉ループ系は、それぞれ

$\displaystyle{(3)\quad \dot{x}= \underbrace{ \left[\begin{array}{cc} -2 & 0\\ 0 & -3 \end{array}\right] }_{A-BF_1} x,\quad \dot{x}= \underbrace{ \left[\begin{array}{cc} -2 & 2\\ 0 & -3 \end{array}\right] }_{A-BF_2} x }$

となって、同じ固有値-2,-3を持ちます。このように多入力系では閉ループ系の固有値を指定しても状態フィードバックは一意に定まりません。これは固有ベクトルの指定に任意性があるからです（ちなみに１入力系では一意に定まります）。そこで、上の２つの閉ループ系の固有ベクトルを求めてみると

$\displaystyle{(4)\quad V_1=\left[\begin{array}{cc} 1 & 0\\ 0 & 1 \end{array}\right],\quad V_2=\left[\begin{array}{cc} 1 & 1\\ 0 & 1 \end{array}\right] }$

また応答は次式となります。

$\displaystyle{(5)\quad \dot{x}= \underbrace{ \left[\begin{array}{cc} e^{-2t} & 0\\ 0 & e^{-3t} \end{array}\right] }_{\exp((A-BF_1)t)=V_1\exp(\Lambda t)V_1^{-1}} x(0) }$
$\displaystyle{(6)\quad \dot{x}= \underbrace{ \left[\begin{array}{cc} 1 & 1\\ 0 & 1 \end{array}\right] \left[\begin{array}{cc} e^{-2t} & 0\\ 0 & e^{-3t} \end{array}\right] \left[\begin{array}{cc} 1 & -1\\ 0 & 1 \end{array}\right] }_{\exp((A-BF_2)t)=V_2\exp(\Lambda t)V_2^{-1}} x(0) }$

明らかに(5)ではモード $e^{-2t}$ と $e^{-3t}$ が分離していますが、(6)ではそれらがミックスされています。このように、多入力系の場合の状態フィードバックの設計は、モードをどのように分布させるかがポイントとなります。

[1] 実固有値の場合

$\displaystyle{(7)\quad (A-BF)v=\lambda v \Rightarrow \underbrace{ \left[\begin{array}{cc} \lambda I-A & B \end{array}\right] }_{G(\lambda)} \left[\begin{array}{c} v\\ \xi \end{array}\right] =0 \quad(\xi=Fv) }$

より、 $v$ と $\xi$ は、サイズ $n\times(n+m)$ の行列 $G(\lambda)$ の零化空間の基底行列に対して、適当な $m$ 次元ベクトル $\delta$ を用いて

$\displaystyle{(8)\quad G(\lambda) \left[\begin{array}{c} N(\lambda)\\ M(\lambda)\\ \end{array}\right] =0 \Rightarrow \left[\begin{array}{c} v\\ \xi \end{array}\right]= \left[\begin{array}{c} N(\lambda)\\ M(\lambda)\\ \end{array}\right]\delta }$

のように表されます。そこで望ましい固有ベクトル $v$ を与えて

$\displaystyle{(9)\quad N(\lambda)\delta=v }$

を解いて $\delta$ を求めます。ただし、望ましい固有ベクトル $v$ のすべての成分を指定できるわけではないので

$\displaystyle{(10)\quad \tilde{R}N(\lambda)\delta=\tilde{R}v =\left[\begin{array}{c} \tilde{v}^s\\ \tilde{v}^u\\ \end{array}\right] \Rightarrow \left[\begin{array}{cc} I_q & 0 \end{array}\right] \tilde{R}N(\lambda)\times\delta=\tilde{v}^s }$

の最小２乗解として $\delta=\delta^*$ を定めます。ここで、セレクタ行列 $\tilde{R}$ により、望ましい固有ベクトル $v$ の指定された $q$ 個の成分だけを集めた $\tilde{v}^s$ が得られるとしています。

このとき、割当可能な $v^*$ と $\xi^*$ は次式で与えられます。

$\displaystyle{(11)\quad \left[\begin{array}{c} v^*\\ \xi^* \end{array}\right]= \left[\begin{array}{c} N(\lambda)\\ M(\lambda)\\ \end{array}\right]\delta^* }$

●もし $A-\lambda I_n$ が正則ならば

$\displaystyle{(12)\quad (\lambda I_n-A)\underbrace{(\lambda I_n-A)^{-1}B}_{N(\lambda)}+B\underbrace{(-I_m)}_{M(\lambda)}=0 }$

より、 $N(\lambda)$ と $M(\lambda)$ が特定できます。(8)より

$\displaystyle{(13)\quad \begin{array}{l} v=(\lambda I_n-A)^{-1}B\delta\\ \xi=-\delta \end{array} \quad\Rightarrow\quad \boxed{v=(A-\lambda I_n)^{-1}B\xi} }$

となって、指定した固有ベクトル $v$ を近似する $\xi$ を直接求めることができます。

●具体的に、(1)に対して、固有値-2,-3を持たせる状態フィードバックを、対応する固有ベクトルが(4)の $V_1$ となるように求めてみます。

$\displaystyle{(14)\quad \begin{array}{l} \underbrace{ \left[\begin{array}{cc} 0-(-2) & 0\\ 0 & -1-(-2) \end{array}\right]^{-1} \left[\begin{array}{cc} 1 & 1\\ 1 & -1 \end{array}\right] }_{(A-\lambda I_n)^{-1}B} \xi_1 =\underbrace{ \left[\begin{array}{c} 1\\ 0 \end{array}\right]}_{V_1(:,1)}\\ \Rightarrow \xi_1= \left[\begin{array}{cc} 1/2 & 1/2\\ 1 & -1 \end{array}\right]^{-1} \left[\begin{array}{c} 1\\ 0 \end{array}\right] = \left[\begin{array}{cc} 1 & 1/2\\ 1 & -1/2 \end{array}\right] \left[\begin{array}{c} 1\\ 0 \end{array}\right] = \left[\begin{array}{cc} 1\\ 1 \end{array}\right] \end{array} }$

$\displaystyle{(15)\quad \begin{array}{l} \underbrace{ \left[\begin{array}{cc} 0-(-3) & 0\\ 0 & -1-(-3) \end{array}\right]^{-1} \left[\begin{array}{cc} 1 & 1\\ 1 & -1 \end{array}\right] }_{(A-\lambda I_n)^{-1}B} \xi_2 =\underbrace{ \left[\begin{array}{c} 0\\ 1 \end{array}\right]}_{V_1(:,2)}\\ \Rightarrow \xi_2= \left[\begin{array}{cc} 1/3 & 1/3\\ 1/2 & -1/2 \end{array}\right]^{-1} \left[\begin{array}{c} 0\\ 1 \end{array}\right] = \left[\begin{array}{cc} 3/2 & 1\\ 3/2 & -1 \end{array}\right] \left[\begin{array}{c} 0\\ 1 \end{array}\right] = \left[\begin{array}{cc} 1\\ -1 \end{array}\right] \end{array} }$

状態フィードバックゲイン $F$ は

$\displaystyle{(16)\quad \left[\begin{array}{cc} \xi_1 & \xi_2 \end{array}\right] =FV_1 \quad\Rightarrow\quad F= \left[\begin{array}{cc} \xi_1 & \xi_2 \end{array}\right]V_1^{-1} = \left[\begin{array}{cc} 1 & 1\\ 1 & -1 \end{array}\right] }$

[2] 複素共役固有値の場合

$\displaystyle{(17)\quad \begin{array}{l} (A-BF)(v_R\pm jv_I)=(\lambda_R\pm j\lambda_I)(v_R\pm jv_I)\\ \Rightarrow (A-BF)v_R=\lambda_Rv_R-\lambda_Iv_I,\ (A-BF)v_I =\lambda_Rv_I+\lambda_Iv_R \\ \Rightarrow \underbrace{ \left[\begin{array}{ccc} \lambda_R I-A & \lambda_I I & B \end{array}\right] }_{G(\lambda)} \left[\begin{array}{cc} v_R & v_I\\ -v_I & v_R\\ \xi_R & \xi_I \end{array}\right] =0 \quad(\xi_R=Fv_R,\xi_I=Fv_I) \end{array} }$

より、 $v_R,v_I$ と $\xi_R,\xi_I$ は、サイズ $n\times(2n+m)$ の行列 $G(\lambda)$ の零化空間の基底行列に対して、適当な $m$ 次元ベクトル $\delta_R,\delta_I$ を用いて

$\displaystyle{(18)\quad G(\lambda) \left[\begin{array}{c} N(\lambda)\\ P(\lambda)\\ M(\lambda) \end{array}\right] =0 \Rightarrow \left[\begin{array}{cc} v_R & v_I\\ -v_I & v_R\\ \xi_R & \xi_I \end{array}\right]= \left[\begin{array}{c} N(\lambda)\\ P(\lambda)\\ M(\lambda) \end{array}\right] \left[\begin{array}{cc} \delta_R & \delta_I \end{array}\right] }$

のように表されます。そこで望ましい $v_R,v_I$ を与えて

$\displaystyle{(19)\quad \left[\begin{array}{c} v_R \\ v_I \end{array}\right]= \underbrace{ \left[\begin{array}{cc} N(\lambda) \\ -P(\lambda) \end{array}\right] }_{\alpha} \delta_R= \underbrace{ \left[\begin{array}{cc} P(\lambda)\\ N(\lambda) \end{array}\right] }_{\beta} \delta_I }$

の解として $\delta_R,\delta_I$ を求めますが、 $\alpha$ と $\beta$ が独立ではないので、その求解には工夫が必要です。ここでは、Edwards and Spurgeon に従って、その手順を述べます。

まず、行列 $\alpha$ に対して、 $\alpha^T$ の零化空間の基底行列を $K_\alpha$ で表すと、 $\alpha^TK_\alpha=0$ より $K_\alpha^T\alpha=0$ を得て、 $K_\alpha^T$ の零化空間の基底行列は $\alpha$ 、すなわち ${\cal N}(K_\alpha^T)={\cal R}(\alpha)$ となることに注意します。同様に行列 $\beta$ に対して、 $\beta^T$ の零化空間の基底行列を $K_\beta$ で表すと、 $K_\beta^T$ の零化空間の基底行列は $\beta$ 、すなわち ${\cal N}(K_\beta^T)={\cal R}(\beta)$ となります。そこで、

$\displaystyle{(20)\quad \underbrace{ \left[\begin{array}{c} K_\alpha^T\\ K_\beta^T \end{array}\right] }_{\gamma^T} K_\gamma=0 \Rightarrow {\cal R}(K_\gamma)\subset {\cal N}(K_\alpha^T)\cap {\cal N}(K_\beta^T)={\cal R}(\alpha)\cap {\cal R}(\beta) }$

に注意して、 $\gamma=[K_\alpha,K_\beta]$ の零化空間の基底行列 $K_\gamma$ を用いて、望ましい $v_R,v_I$ を与えて

$\displaystyle{(21)\quad \left[\begin{array}{c} v_R \\ v_I \end{array}\right]= K_\gamma \delta }$

を解いて $\delta$ を求めます。ただし、望ましい $v_R,v_I$ のすべての成分を指定できるわけではないので

$\displaystyle{(22)\quad \tilde{R}K_\gamma\delta=\tilde{R}v =\left[\begin{array}{c} \tilde{v}_R^s\\ \tilde{v}_I^s\\ \tilde{v}_R^u\\ \tilde{v}_I^u\\ \end{array}\right] \Rightarrow \left[\begin{array}{cc} I_{2q} & 0 \end{array}\right] \tilde{R}K_\gamma\times\delta= \left[\begin{array}{c} \tilde{v}_R^s\\ \tilde{v}_I^s\\ \end{array}\right] }$

の最小２乗解として $\delta=\delta^*$ を定めます。ここで、セレクター行列 $\tilde{R}$ により、望ましい $v_R,v_I$ の指定された $2q$ 個の成分だけを集めた $\tilde{v}_R^s,\tilde{v}_I^s$ が得られるとしています。

このとき、割当可能な $v_R^*,v_I^*$ は次式で与えられます。

$\displaystyle{(23)\quad \left[\begin{array}{c} v_R^*\\ v_I^* \end{array}\right]= K_\gamma\delta^* }$

これは(14)の $\delta_R^*,\delta_I^*$ を、それぞれ $\alpha^{\dag}K_\gamma\delta^*$ と $\beta^{\dag}K_\gamma\delta^*$ で定めれば（ $\dag$ は疑似逆行列を表す）、対応する $\xi_R^*,\xi_I^*$ は次式で与えられます。

$\displaystyle{(24)\quad \left[\begin{array}{c} \xi_R^*\\ \xi_I^* \end{array}\right]= \left[\begin{array}{c} M(\lambda)\alpha^{\dag}K_\gamma\delta^*\\ M(\lambda)\beta^{\dag}K_\gamma\delta^*\\ \end{array}\right] }$

[3] 複素共役固有値の場合（別のアプローチ）

(17)は次式のようにも表すことができます。

$\displaystyle{(25)\quad \begin{array}{l} (A-BF)(v_R\pm jv_I)=(\lambda_R\pm j\lambda_I)(v_R\pm jv_I)\\ \Rightarrow (A-BF)v_R=\lambda_Rv_R-\lambda_Iv_I,\ (A-BF)v_I =\lambda_Rv_I+\lambda_Iv_R \\ \Rightarrow %\underbrace{ \left[\begin{array}{cccc} \lambda_R I-A & -\lambda_I I & B & 0\\ \lambda_I I & \lambda_R I-A & 0 & B\\ \end{array}\right] %}_{G(\lambda)} \left[\begin{array}{cc} v_R \\ v_I \\ \xi_R \\ \xi_I \end{array}\right] =0 \quad(\xi_R=Fv_R,\xi_I=Fv_I) \end{array} }$

したがって、[2]において

$\displaystyle{(26)\quad \begin{array}{l} G(\lambda)\leftarrow \left[\begin{array}{cccc} \lambda_R I-A & -\lambda_I I & B & 0\\ \lambda_I I & \lambda_R I-A & 0 & B \end{array}\right]\\ v\leftarrow\left[\begin{array}{c} v_R\\ v_I \end{array}\right],\ \xi\leftarrow\left[\begin{array}{c} \xi_R\\ \xi_I \end{array}\right],\ \delta\leftarrow\left[\begin{array}{c} \delta_R\\ \delta_I \end{array}\right] \end{array} }$

と置き換えて同様の議論ができそうです（[2]の議論と等価であるかは検討中です）。

●もし $A-\lambda_R I_n$ が正則ならば、公式

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

を参照すると

$\displaystyle{(28)\quad \left[\begin{array}{cc} \lambda_R I-A & -\lambda_I I\\ \lambda_I I & \lambda_R I-A \end{array}\right]^{-1} }$

を計算できます。そこで

$\displaystyle{(29)\quad \begin{array}{l} \left[\begin{array}{cc} \lambda_R I-A & -\lambda_I I\\ \lambda_I I & \lambda_R I-A \end{array}\right] \underbrace{ \left[\begin{array}{cc} \lambda_R I-A & -\lambda_I I\\ \lambda_I I & \lambda_R I-A \end{array}\right]^{-1} \left[\begin{array}{cc} B & 0\\ 0 & B \end{array}\right] }_{N(\lambda)}\\ + \left[\begin{array}{cc} B & 0\\ 0 & B \end{array}\right] \underbrace{(-I_{2m})}_{M(\lambda)}=0 \end{array} }$

より、 $N(\lambda)$ と $M(\lambda)$ が特定できます。(8)に相当する式より

$\displaystyle{(30)\quad \boxed{ \left[\begin{array}{c} v_R\\ v_I \end{array}\right] = \left[\begin{array}{cc} A-\lambda_R I & \lambda_I I\\ -\lambda_I I & A-\lambda_R I \end{array}\right]^{-1} \left[\begin{array}{cc} B & 0\\ 0 & B \end{array}\right] \left[\begin{array}{c} \xi_R\\ \xi_I \end{array}\right]} }$

となって、指定した固有ベクトルの実部 $v_R$ と虚部 $v_I$ を近似する $\xi_R$ と $\xi_I$ を直接求めることができます。

●具体的に、(1)に対して、固有値 $-1\pm j$ を持たせる状態フィードバックを、対応する固有ベクトルが

$\displaystyle{(31)\quad \underbrace{ \left[\begin{array}{c} 1\\ 0 \end{array}\right]}_{v_R} \pm j \underbrace{ \left[\begin{array}{c} 0\\ 1 \end{array}\right]}_{v_I} }$

となるように求めてみます。

$\displaystyle{(32)\quad \begin{array}{l} \underbrace{ \left[\begin{array}{cccc} 0-(-1) & 0 & 1 & 0\\ 0 & -1-(-1) & 0 & 1\\ -1 & 0 & 0-(-1) & 0\\ 0 & -1 & 0 & -1-(-1) \end{array}\right]^{-1} }_{\left[\begin{array}{cc} A-\lambda_R I & \lambda_I I\\ -\lambda_I I & A-\lambda_R I \end{array}\right]^{-1}}\\ \times\underbrace{ \left[\begin{array}{cccc} 1 & 1 & 0 & 0\\ 1 & -1& 0 & 0\\ 0 & 0 & 1 & 1\\ 0 & 0 & 1 & -1 \end{array}\right] }_{\left[\begin{array}{cc} B & 0\\ 0 & B \end{array}\right]} \left[\begin{array}{c} \xi_R\\ \xi_I \end{array}\right] =\underbrace{ \left[\begin{array}{c} 1\\ 0\\ 0\\ 1 \end{array}\right]}_{\left[\begin{array}{c} v_R\\ v_I \end{array}\right]}\\ \Rightarrow \left[\begin{array}{c} \xi_R\\ \xi_I \end{array}\right]= \left[\begin{array}{cccc} 1/2 & 1/2 & -1/2 & -1/2\\ 0 & 0 & -1 & 1\\ 1/2 & 1/2 & 1/2 & 1/2\\ 1 & -1& 0 & 0 \end{array}\right]^{-1} \left[\begin{array}{c} 1\\ 0\\ 0\\ 1 \end{array}\right]= \left[\begin{array}{c} 1\\ 0\\ -1/2\\ -1/2 \end{array}\right] \end{array} }$

状態フィードバックゲイン $F$ は

$\displaystyle{(33)\quad \begin{array}{l} \left[\begin{array}{cc} \xi_R & \xi_I \end{array}\right] =F \left[\begin{array}{cc} v_R & v_I \end{array}\right]\\ \Rightarrow F= \left[\begin{array}{cc} \xi_R & \xi_I \end{array}\right] \left[\begin{array}{cc} v_R & v_I \end{array}\right]^{-1} = \left[\begin{array}{cc} 1 & -1/2\\ 0 & -1/2 \end{array}\right] \end{array} }$

演習C42…Flipped Classroom
$1^\circ$ 　(16)と(33)の計算を行う次のプログラムを実行せよ。

MATLAB

%test_vplace2.m
%-----
 clear all, close all
 A=[0 0;
    0 -1];
 B=[1 1;
    1 -1];
%-----
 lambda=[-2 -3];
 nocomp=0;
%固有ベクトルの要素の指定の有無 (１：指定、０：任意)
 specpos=[ 1  1; 
           1  1];       
%固有ベクトルの要素の値 (０：零指定、１：非零指定、－１：非零任意)       
 specent=[ 1  0;
           0  1];
 F=vplace2(A,B,lambda,nocomp,specpos,specent)
 [V,pl]=eig(A-B*F)
 VV=V*diag([1/V(2,1),1/V(1,2)])
 disp("---")
%-----
 lambda=[-1 1];
 nocomp=1;
%固有ベクトルの要素の指定の有無 (１：指定、０：任意)
 specpos=[ 1  1; 
           1  1];        
%固有ベクトルの要素の値 (０：零指定、１：非零指定、－１：非零任意)       
 specent=[ 1  0;
           0  1];        
 F=vplace2(A,B,lambda,nocomp,specpos,specent) 
 [V,pl]=eig(A-B*F) 
 VV=V*diag([1/V(1,1),1/V(1,2)])

$2^\circ$ 　１入力系の場合は、specposとspecentは乱数で与えてよい、なぜか？

MATLAB

%test_vplace2a.m
%-----
 clear all, close all
 A=[0 1;
    0 0];
 B=[0;
    1];
%-----
%lambda=[-1 -1];
 lambda=[-2 -3];  
 nocomp=0;
%固有ベクトルの要素の指定の有無 (１：指定、０：任意)
 specpos=rand(2,2)
%固有ベクトルの要素の値 (０：零指定、１：非零指定、－１：非零任意)       
 specent=rand(2,2) 
 F=vplace2(A,B,lambda,nocomp,specpos,specent) 
 [V,pl]=eig(A-B*F)
%-----
 lambda=[-1 1];
 nocomp=1;
%固有ベクトルの要素の指定の有無 (１：指定、０：任意)
 specpos=rand(2,2)
%固有ベクトルの要素の値 (０：零指定、１：非零指定、－１：非零任意)       
 specent=rand(2,2) 
 F=vplace2(A,B,lambda,nocomp,specpos,specent) 
 [V,pl]=eig(A-B*F) 
%-----
%eof

Note C42 固有ベクトル設定プログラム
これはEdwards and Spurgeonによるプログラムspvplを基に作成できます。

スイッチング関数(2)

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

スイッチング関数(2)

[1] 状態方程式とスイッチング関数は次の標準形をとるように座標変換されていると仮定します。

$\displaystyle{(1)\quad \underbrace{ \left[\begin{array}{c} \dot x_1(t)\\ \dot x_2(t) \end{array}\right] }_{\dot x(t)} = \underbrace{ \left[\begin{array}{cc} A_{11} & A_{12} \\ A_{21} & A_{22} \\ \end{array}\right] }_{A} \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x(t)} + \underbrace{ \left[\begin{array}{c} 0\\ B_2 \end{array}\right] }_{B} u(t) }$

$\displaystyle{(2)\quad s(t)= \underbrace{ \left[\begin{array}{cc} S_1 & S_2 \\ \end{array}\right] }_{S} \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x(t)} = \underbrace{S_2 \left[\begin{array}{cc} M & I \\ \end{array}\right] }_{S} \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x(t)} \ (M=S_2^{-1}S_1) }$

このとき、等価制御

$\displaystyle{(3)\quad u_{eq}(t)=-(SB)^{-1}SAx(t) }$

による閉ループ系は次式で表されました。

$\displaystyle{(4)\quad \dot{x}(t)=\underbrace{(I_n-B(SB)^{-1}S)A}_{A_{eq}}x(t)\quad(t\ge t_s) }$

ここで $A$ 、 $B$ 、 $S$ は、(1)と(2)で与えられるとすると

$\displaystyle{(5)\quad \begin{array}{l} A_{eq}= \underbrace{ \left[\begin{array}{cc} I & 0 \\ -M & 0 \\ \end{array}\right] }_{I_n-B(SB)^{-1}S} \left[\begin{array}{cc} A_{11} & A_{12} \\ A_{21} & A_{22} \\ \end{array}\right]= \left[\begin{array}{cc} A_{11} & A_{12} \\ -MA_{11} & -MA_{12} \\ \end{array}\right]\\ = \left[\begin{array}{cc} I & 0 \\ -M & I \\ \end{array}\right] \left[\begin{array}{cc} \bar{A}_{11} & A_{12} \\ 0 & 0 \\ \end{array}\right] \left[\begin{array}{cc} I & 0 \\ -M & I \\ \end{array}\right]^{-1} \end{array} }$

となって、 $\bar{A}_{11}$ が安定行列なので、 $A_{eq}$ は $m$ 個の零固有値を持ちます。一方、 $A_{eq}$ の非零固有値 $\lambda_i\ (i=1,\cdots,n-m)$ に対応する固有ベクトルを $v_i$ とすると、

$\displaystyle{(6)\quad SA_{eq}=0\Rightarrow SA_{eq}v_i=0\Rightarrow \lambda_iSv_i=0 \Rightarrow Sv_i=0 }$

すなわち

$\displaystyle{(7)\quad { \boxed{S[v_1\cdots v_{n-m}]=0 %　\Rightarrow [v_1\cdots v_{n-m}]^TS^T=0}} }$

が成り立ちます。したがって、まず等価制御による閉ループ系の $A_{eq}$ が望ましい固有値・固有ベクトルをもつように設計して、(7)から $S$ を決定することが考えられます。

● １入力系の場合は固有ベクトル設定の自由度はありませんが、多入力系の場合は固有ベクトルの向きを設定する必要があり、これは元の状態方程式(状態変数の並び)に対して与えます。そこで(7)が元の状態方程式ではどう変わるかを調べておきます。

元の状態方程式

$\displaystyle{(8)\quad \begin{array}{l} \dot{x}(t)=Ax(t)+Bu(t)+f(t,x,u)\\ (x(t)\in{\rm\bf R}^n, u(t)\in{\rm\bf R}^m, f(t,x,u)\in{\rm\bf R}^n) \end{array} }$

に対する標準形(1)は、適当な直交行列 $T_r$ 座標による座標変換

$\displaystyle{(9)\quad x_{r}(t)=T_rx(t) \Leftrightarrow x(t)=T_r^Tx_r(t) }$

を行って（ $T_r^{-1}=T_r^T$ に注意）

$\displaystyle{(10)\quad \begin{array}{l} \underbrace{ \left[\begin{array}{c} \dot{x}_{r1}(t)\\ \dot{x}_{r2}(t) \end{array}\right] }_{\dot{x}_r(t)} = \underbrace{ \left[\begin{array}{cc} A_{r11} & A_{r12} \\ A_{r21} & A_{r22} \\ \end{array}\right] }_{A_r=T_rAT_r^T} \underbrace{ \left[\begin{array}{c} x_{r1}(t)\\ x_{r2}(t) \end{array}\right] }_{x_r(t)} + \underbrace{ \left[\begin{array}{c} 0\\ \Sigma_1V^T \end{array}\right] }_{B_r=T_rB} u(t)\\ + \underbrace{ \left[\begin{array}{l} f_{ru}(t,x)\\ f_{rm}(t,x,u) \end{array}\right] }_{f_r(t,x,u)=T_rf(t,x,u)} \end{array} }$

のように得られたものでした。また(2)は次式に

$\displaystyle{(11)\quad s(t)= \underbrace{ \left[\begin{array}{cc} S_{r1} & S_{r2} \\ \end{array}\right] }_{S_r=ST_r^T} \underbrace{ \left[\begin{array}{c} x_{r1}(t)\\ x_{r2}(t) \end{array}\right] }_{x_r(t)=T_rx(t)} }$

に、(7)は次式に相当します。

$\displaystyle{(12)\quad S_r[v_{r1}\cdots v_{r,n-m}]=0 }$

いま、元の状態空間での固有ベクトル $v_1\cdots v_{n-m}$ は、座標変換後には

$\displaystyle{(13)\quad Av_i=\lambda_iv_i\Rightarrow T_rAT_r^TT_rv_i=\lambda_iT_rv_i\Rightarrow v_{ri}=T_rv_i }$

に変わることに注意します。(12)から

$\displaystyle{(14)\quad S_rT_rT_r^T[v_{r1}\cdots v_{r,n-m}]=S[v_1\cdots v_{n-m}]=0 }$

を得ます。これから元の状態方程式に対するスイッチング関数と固有ベクトルについても、(12)と同様の関係が成り立つと言えます。すなわち固有ベクトル $v_1\cdots v_{n-m}$ は元の状態空間で与えておいて、(14)から $S$ を決定してよいことがわかります。

演習C41…Flipped Classroom

$1^\circ$ スイッチング関数決定プログラムswfvplを用いて、図1の台車駆動型倒立振子を安定化する次のSMC則を求めよ。

$\displaystyle{{ \boxed{\begin{array}{l} u(t)=\underbrace{u_\ell(t)}_{linear\ control}+\underbrace{u_n(t)}_{switching\ component}\\ u_\ell(t)=-\underbrace{(SB)^{-1}(SA-\Phi S)}_{L=L_{eq}+L_{phi}}x(t) \quad(\Phi:stable\ matrix)\\ u_n(t)=-\underbrace{(SB)^{-1}\rho(t,x)}_{L_n}\frac{P_2s(t)}{||P_2s(t)||}\quad(P_2\Phi+\Phi^TP_2=-I) \end{array}}} }$

図1 台車駆動型倒立振子

MATLAB

%cCIP_smc2.m
%-----
 clear all, close all
 global mc m ell g th0
 mc=1; m=0.1; ell=0.2; g=9.8;
 ths=0; %input('ths   = <0,180> ')/180*pi;
 th0=3/180*pi; %input('th(0) = <0,180> ')/180*pi;
%-----
 A=[zeros(2,2) eye(2);zeros(2,4)];
 A(3,1)=0; 
 A(3,2)=-3*m*g/(m+4*mc); 
 A(4,1)=0; 
 A(4,2)=3*(m+mc)*g/((m+4*mc)*ell);  
 B=zeros(4,1);
 B(3)=4/(m+4*mc);
 B(4)=-3/((m+4*mc)*ell); 
 C=diag([100 180/pi])*eye(2,4); 
%-----
 lambda=[-0.8 0.5 -3];
 nocomp=1;
 specpos=rand(4,3)
 specent=rand(4,3)
 S=swfvpl(A,B,lambda,nocomp,specpos,specent)
%-----
 Phi=-0.1;
 P2=0.5*inv(-Phi);
 Check=P2*Phi+Phi*P2
 P2S=P2*S;
 Leq=inv(S*B)*S*A;
 LPhi=-inv(S*B)*Phi*S;
 L=Leq+LPhi;
 rho=0.6;
 Ln=inv(S*B)*rho;
%-----
 x0=[0;th0;0;0];
 sim('CIP_smc_2015a.mdl')
%-----
%eof

Note C41 スイッチング関数決定プログラム２
Edwards and Spurgeon によって、(26)から、スイッチング関数の行列 $S$ を求めるプログラムswfvplが開発されています。

MATLAB

%test_swfvpl.m
%-----
 clear all, close all
 A=[0 1 0;
    0 0 1;
    0 0 0];
 B=[0;
    0;
    1];
%-----
 lambda=[-2 -3];
 nocomp=0;
%固有ベクトルの要素の指定の有無 (１：指定、０：任意)
 specpos=rand(3,3)
%固有ベクトルの要素の値 (０：零指定、１：非零指定、－１：非零任意)       
 specent=rand(3,3) 
 S=swfvpl(A,B,lambda,nocomp,specpos,specent) 
 Leq=inv(S*B)*S*A;
 [V,pl]=eig(A-B*Leq) 
%-----
 lambda=[-1 1];
 nocomp=1;
%固有ベクトルの要素の指定の有無 (１：指定、０：任意)
 specpos=rand(3,3)
%固有ベクトルの要素の値 (０：零指定、１：非零指定、－１：非零任意)       
 specent=rand(3,3) 
 S=swfvpl(A,B,lambda,nocomp,specpos,specent) 
 Leq=inv(S*B)*S*A;
 [V,pl]=eig(A-B*Leq) 
%-----
%eof

補遺：LQ制御

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

LQ制御…Homework

[1]　可制御な制御対象

$\displaystyle{(1)\quad %\left\{\begin{array}{ll} \dot{x}(t)=Ax(t)+Bu(t)\quad (x(t)\in{\bf R}^n,u(t)\in{\bf R}^m)\\ %z(t)=Cx(t)&(z(t)\in{\bf R}^q) %\end{array}\right. }$

を安定化する状態フィードバック

$\displaystyle{(2)\quad u(t)=-Fx(t) }$

の決定法を考えます。一つの方法は，閉ループ系

$\displaystyle{(3)\quad %\left\{\begin{array}{l} \dot{x}(t)=(A-BF)x(t) \\ %z(t)=Cx(t) %\end{array}\right. }$

の時間応答に関する評価規範として，２次形式評価関数

$(4)\quad \boxed{\begin{array}{l} \displaystyle{J=\int_0^\infty (z^T(t)Qz(t)+u^T(t)Ru(t))\,dt\quad (Q>0,R>0)}\\ \displaystyle{z(t)=Cx(t)\quad (z(t)\in{\bf R}^q)} \end{array}}$

を設定し，これを最小化する問題を解くことです。ただし、 $(A,C)$ は可観測対とします。これによる状態フィードバックのゲイン行列 $F$ は，リッカチ方程式

$\displaystyle{(5)\quad {\boxed{\Pi A+A^T\Pi-\Pi BR^{-1}B^T\Pi+C^TQC=0}} }$

の解 $\Pi>0$ を用いて，次式で与えられます。

$\displaystyle{(6)\quad {\boxed{F=R^{-1}B^T\Pi}} }$

この証明はNoteに示しています。

Note A52-1 行列による微分　

いま、任意の行列 $X$ の $(i,j)$ 要素を $x_{ij}=[X]_{ij}$ で表すとき、スカラ関数 $f$ を行列変数 $X$ の各要素で微分して得られる行列を $\displaystyle{[\frac{\partial f}{\partial X}]_{ij}=\frac{\partial f}{\partial x_{ij}}}$ で定義します。このとき、行列のトレースについて、次が成り立ちます。

$\displaystyle{(1)\ {\rm tr}AB={\rm tr}BA}$
$\displaystyle{(2)\ \frac{\partial}{\partial X}{\rm tr}AXB=A^TB^T}$
$\displaystyle{(3)\ \frac{\partial}{\partial X}{\rm tr}AX^TB=BA}$
$\displaystyle{(4)\ \frac{\partial}{\partial X}{\rm tr}AXBX^T=A^TXB^T+AXB}$

実際、

$\displaystyle{(1)\ \sum_{i}[AB]_{ii}=\sum_{i}\sum_{j}a_{ij}b_{ji} =\sum_{j}\sum_{i}b_{ji}a_{ij}=\sum_{j}[BA]_{jj}}$

$\displaystyle{(2)\ [\frac{\partial}{\partial X}{\rm tr}AXB]_{ij} =\frac{\partial}{\partial x_{ij}}\sum_{k}[AXB]_{kk} =\frac{\partial}{\partial x_{ij}}\sum_{k}\sum_{i,j}a_{ki}x_{ij}b_{jk}}$
$\displaystyle{=\sum_{k}b_{jk}a_{ki}=[BA]_{ji}=[A^TB^T]_{ij}}$

$\displaystyle{(3)\ [\frac{\partial}{\partial X}{\rm tr}AX^TB]_{ij} =\frac{\partial}{\partial x_{ij}}\sum_{k}[AX^TB]_{kk} =\frac{\partial}{\partial x_{ij}}\sum_{k}\sum_{i,j}a_{ki}x_{ji}b_{jk}}$
$\displaystyle{=\sum_{k}b_{ik}a_{kj}=[BA]_{ij}}$

$\displaystyle{(4)\ \frac{\partial}{\partial X}{\rm tr}AXBX^T =\frac{\partial}{\partial X}{\rm tr}(AXB)X^T +\frac{\partial}{\partial X}{\rm tr}AX(BX^T)}$
$\displaystyle{=AXB+A^TXB^T}$

Note A52-2 LQ制御問題の解法　

可制御かつ可観測な $n$ 次系

$\displaystyle{(1)\quad \dot{x}(t)=Ax(t)+Bu(t),\ y(t)=Cx(t) }$

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

$\displaystyle{(2)\quad u(t)=-Fx(t) }$

による閉ループ系

$\displaystyle{(3)\quad \begin{array}{l} \dot{x}(t)=(A-BF)x(t),\ z(t)=Cx(t)\\ A_F=A-BF:\ stable\ matrix \end{array} }$

に対して、評価関数

$\displaystyle{(4)\quad J=\int_0^\infty(z^T(t)Qz(t)+u^T(t)Ru(t))\,dt }$

を最小化するように $F$ を決める問題を考えます。

閉ループ系における状態の振る舞いは次式で与えられます。

$\displaystyle{(5)\quad x(t)=\exp({A_Ft})x(0) }$

ここで、１次系の場合は初期状態は $x(0)\ne0$ であればよかったのですが、一般の場合はインパルス応答となるように $B$ の列ベクトル $B^{(i)}$ を考えます。各インパルス応答

$\displaystyle{(6)\quad x^{(i)}(t)=\exp({A_Ft})B^{(i)} }$

に対する評価関数 $J$ の総和は

$(7)\quad \begin{array}{l} \displaystyle{\sum\int_0^\infty x^{(i)}\,^T(t)(C^TQC+F^TRF)x^{(i)}(t)\,dt}\\ \displaystyle{=\sum B^{(i)}\,^T\underbrace{\int_0^\infty \exp(A_F^Tt)(C^TQC+F^TRF)\exp(A_Ft)\,dt}_{\Pi}B^{(i)}}\\ \displaystyle{=\sum B^{(i)}\,^T\Pi B^{(i)}={\rm tr}\ \Pi BB^T} \end{array} }$

と書けます。いま $\Gamma$ をラグランジュの未定定数として

$\displaystyle{(8)\quad J'={\rm tr}\ \Pi x(0)x^T(0)+{\rm tr}\ ((\Pi A_F+A_F^T\Pi+C^TQC+F^TRF)\Gamma) }$

を最小化する問題を考えます。ここで、制約条件は、リャプノフ方程式と呼ばれる

$\displaystyle{(9)\quad \Pi A_F+A_F^T\Pi+C^TQC+F^TRF=0 }$

ですが、 $\Pi>0$ は $A_F$ が安定行列を意味することに注意します。

そこで、必要条件として次を得ます。

$\displaystyle{(11)\quad \frac{\partial J'}{\partial\Pi}&=&BB^T+\Gamma A_F^T+A_F\Gamma=0\Rightarrow \Gamma>0 }$
$\displaystyle{(12)\quad \frac{\partial J'}{\partial F}&=&2(RF-B^T\Pi)\Gamma=0\Rightarrow F=R^{-1}B^T\Pi }$
$\displaystyle{(13)\quad \frac{\partial J'}{\partial\Gamma}&=&\Pi A_F+A_F^T\Pi+C^TQC+F^TRF=0 }$

ここで、第２式から得られる $F=R^{-1}B^T\Pi$ を第３式に代入して

$\displaystyle{(14)\quad \begin{array}{l} \Pi (A-BR^{-1}B^T\Pi)+(A-BR^{-1}B^T\Pi)^T\Pi+C^TQ\nonumber\\ +(R^{-1}B^T\Pi)^TRR^{-1}B^T\Pi=0 \end{array} }$

すなわち、リッカチ方程式と呼ばれる $\Pi$ の行列方程式

$\displaystyle{(15)\quad \Pi A+A^T\Pi-\Pi BR^{-1}B^T\Pi+C^TQC=0 }$

を得ます。これから $\Pi>0$ を求めて、 $F$ は

$\displaystyle{(16)\quad F=R^{-1}B^T\Pi }$

のように得られます。このような制御方式をLQ制御と呼びます。

一方、十分性の議論は次のように行われます。まず、被積分項は次のように表すことができます。

$\displaystyle{(17)\quad \begin{array}{lll} &&y^TQy+u^TRu=(u+R^{-1}B^T\Pi x)^TR(u+R^{-1}B^T\Pi x)\nonumber\\ &&-\underbrace{\frac{d}{dt}x^T\Pi x}_{2x^T\Pi\dot{x}} \end{array} }$

実際、右辺に $\dot{x}=Ax+Bu$ を代入し、リッカチ方程式を用いると

$\displaystyle{(18)\quad \begin{array}{lll} &&{\rm RHS}=u^TRu+2x^T\Pi Bu+x^T\Pi BR^{-1}B^T\Pi x-2x^T\Pi(Ax+Bu)\nonumber\\ &&=u^TRu+x^T(\Pi A+A^T\Pi+C^TQC)x-2x^T\Pi Ax={\rm LHS} \end{array} }$

したがって、上記の両辺を積分して

$\displaystyle{(19)\quad J=\int_0^\infty (u+R^{-1}B^T\Pi x)^TR(u+R^{-1}B^T\Pi x)\,dt-\left[x^T\Pi x\right]_0^\infty }$

を得ます。ここで、 $x(t)=\exp(A_Ft)x(0)\rightarrow 0\ (t\rightarrow\infty)$ を前提とするので

$\displaystyle{(20)\quad J=\int_0^\infty (u+R^{-1}B^T\Pi x)^TR(u+R^{-1}B^T\Pi x)\,dt+x^T(0)\Pi x(0) }$

を得ます。これから $u=-R^{-1}B^T\Pi x$ が評価関数を最小化することが分かります。

スイッチング関数(1)

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

スイッチング関数(1)

[1] 状態方程式とスイッチング関数は次の標準形をとるように座標変換されていると仮定します。

これに対して、座標変換

$\displaystyle{(3)\quad \begin{array}{l} \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right] }_{\bar{x}(t)} = \underbrace{ \left[\begin{array}{cc} I & 0 \\ S_1 & S_2 \\ \end{array}\right] }_{T_s} \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x(t)}\\ \Leftrightarrow \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x(t)} = \underbrace{ \left[\begin{array}{cc} I & 0 \\ -S_2^{-1}S_1 & S_2^{-1} \\ \end{array}\right] }_{T_s^{-1}} \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right] }_{\bar{x}(t)} \end{array} }$

を行って、次式を得ます。

$\displaystyle{(4)\quad \underbrace{ \left[\begin{array}{c} \dot x_1(t)\\ \dot s(t) \end{array}\right] }_{\dot{\bar{x}}(t)} = \underbrace{ \left[\begin{array}{cc} \bar{A}_{11} & \bar{A}_{12} \\ \bar{A}_{21} & \bar{A}_{22}\\ \end{array}\right] }_{\bar{A}=T_s A T_s^{-1}} \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right] }_{\bar{x}(t)} + \underbrace{ \left[\begin{array}{c} 0\\ \bar{B}_2 \end{array}\right] }_{\bar{B}=T_sB} u(t) }$

$\displaystyle{(5)\quad \begin{array}{l} \bar{A}_{11}=A_{11}-A_{12}M\\ \bar{A}_{12}=A_{12}S_2^{-1}\\ \bar{A}_{21}=S_2(M\bar{A}_{11}+A_{21}-A_{22}M)\\ \bar{A}_{22}=S_2(M{A}_{12}+A_{22})S_2^{-1}\\ \bar{B}_{2}=S_2B_2 \end{array} }$

ここで、 $s(t)=0$ 、すなわち、状態が超平面内に拘束されているとすると

$\displaystyle{(6)\quad s(t)=0 \Rightarrow Sx(t)=S_1x_1(t)+S_2x_2(t)=0\Rightarrow x_2(t)=-\underbrace{S_2^{-1}S_1}_Mx_1(t)}}$

となって、 $x_2$ の振る舞いは $x_1$ によって決まります。一方、 $x_1$ の振る舞いは

$\displaystyle{(7)\quad \dot{x}_1(t)=\underbrace{(A_{11}-A_{12}M)}_{\bar{A}_{11}}x_1(t)}$

によって決まります。これは

$\displaystyle{(8)\quad \dot{x}_1(t)=A_{11}x_1(t)+A_{12}x_2(t)}$

に対する安定化状態フィードバック

$\displaystyle{(9)\quad x_2(t)=-Mx_1(t)}$

による閉ループ系とみなすことができます。

したがって、 $\bar{A}_{11}$ が安定行列となるようにスイッチング関数を選ぶためには、(8)に対する状態フィードバックによる安定化問題を解けばよいことが分かります。

[2]　そこで、次の２次形式評価関数を最小にするように $M$ を決定する方法が提案されています。

$\displaystyle{(10)\quad \boxed{J=\int_0^\infty \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right]^T }_{x^T(t)} \underbrace{ \left[\begin{array}{cc} Q_{11} & Q_{12} \\ Q_{21} & Q_{22} \end{array}\right] }_{Q} \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x(t)} dt} }$

これは、恒等式

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

を用いて

$\displaystyle{(12)\quad \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right]^T \left[\begin{array}{cc} Q_{11} & Q_{12} \\ Q_{21} & Q_{22} \end{array}\right] \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] =\left[\begin{array}{c} x_1(t)\\ Q_{12}Q_{22}^{-1}x_1(t)+x_2(t) \end{array}\right]^T }$
$\displaystyle{ \times \left[\begin{array}{cc} Q_{11}-Q_{12}Q_{22}^{-1}Q_{21} & 0 \\ 0 & Q_{22} \end{array}\right] \left[\begin{array}{c} x_1(t)\\ Q_{12}Q_{22}^{-1}x_1(t)+x_2(t) \end{array}\right] }$

より

$\displaystyle{(13)\quad J=\int_0^\infty(x_1^T(t)\underbrace{(Q_{11}-Q_{12}Q_{22}^{-1}Q_{21})}_{\hat{Q}}x_1(t)}$
$\displaystyle{+\underbrace{(Q_{12}Q_{22}^{-1}x_1(t)+x_2(t))^T}_{\hat{u}^T}\underbrace{Q_{22}}_{\hat{R}}\underbrace{(Q_{12}Q_{22}^{-1}x_1(t)+x_2(t))}_{\hat{u}})dt}$

と書けます。ここで、(8)と(6)をそれぞれ

$\displaystyle{(14)\quad \dot{x}_1(t)=\underbrace{(A_{11}-A_{12}Q_{22}^{-1}Q_{21})}_{\hat{A}}x_1(t)+\underbrace{A_{12}}_{\hat{B}}\underbrace{(Q_{22}^{-1}Q_{21}x_1(t) +x_2(t))}_{\hat{u}(t)}}$

$\displaystyle{(15)\quad \underbrace{Q_{22}^{-1}Q_{21}x_1(t) +x_2(t)}_{\hat{u}(t)}=-\underbrace{(M-Q_{22}^{-1}Q_{21})}_{\hat{F}}x_1(t)}$

と書き直してみますと、次のような最適化制御問題の定式化ができていることがわかります。すなわち

$\displaystyle{(16)\quad \boxed{\begin{array}{l} \dot{x}_1(t)=\hat{A}x_1(t)+\hat{B}\hat{u}(t)\\ \hat{A}=A_{11}-A_{12}Q_{22}^{-1}Q_{21}\\ \hat{B}=A_{12} \end{array}} }$

を２次系式評価関数

$\displaystyle{(17)\quad \boxed{\begin{array}{l} J=\int_0^\infty(x_1^T(t)\hat{Q}x_1(t)+\hat{u}^T(t)\hat{R}\hat{u}(t))dt\\ \hat{Q}=Q_{11}-Q_{12}Q_{22}^{-1}Q_{21}\\ \hat{R}=Q_{22} \end{array}} }$

を最小にするように

$\displaystyle{(18)\quad \boxed{\begin{array}{l} \hat{u}(t)=-\hat{F}x_1(t)\\ \hat{F}=M-Q_{22}^{-1}Q_{21} \end{array}} }$

を求める問題です。この解は、リッカチ方程式

$\displaystyle{(19)\quad \boxed{\Pi\hat{A}+\hat{A}^T\Pi-\Pi\hat{B}\hat{R}^{-1}\hat{B}^T\Pi+\hat{Q}=0}}$

を解いて

$\displaystyle{(20)\quad \boxed{\hat{F}=\hat{R}^{-1}\hat{B}^T\Pi}}$

を求め、次式から $M$ を決めればよいことがわかります。

$\displaystyle{(21)\quad \boxed{M=\hat{F}+Q_{22}^{-1}Q_{21}}}$

演習C31…Flipped Classroom

$1^\circ$ スイッチング関数決定プログラムswflqrを用いて、図1の台車駆動型倒立振子を安定化する次のSMC則を求めよ。

図1 台車駆動型倒立振子

MATLAB

%cCIP_smc1.m
%-----
 clear all, close all
 global mc m ell g th0
 mc=1; m=0.1; ell=0.2; g=9.8;
 ths=0; %input('ths   = <0,180> ')/180*pi;
 th0=3/180*pi; %input('th(0) = <0,180> ')/180*pi;
%-----
 A=[zeros(2,2) eye(2);zeros(2,4)];
 A(3,1)=0; 
 A(3,2)=-3*m*g/(m+4*mc); 
 A(4,1)=0; 
 A(4,2)=3*(m+mc)*g/((m+4*mc)*ell);  
 B=zeros(4,1);
 B(3)=4/(m+4*mc);
 B(4)=-3/((m+4*mc)*ell); 
 C=diag([100 180/pi])*eye(2,4); 
%-----
 Mr=0.5; Mth=3/180*pi; 
 Tc=1; Tpen=0.25*2*pi*sqrt(2*ell/g);
 Q=diag([1/Mr^2,1/Mth^2,Tc^2/Mr^2,Tpen^2/Mth^2]);
 S=swflqr(A,B,Q)
%-----
 Phi=-0.1;
 P2=0.5*inv(-Phi);
 Check=P2*Phi+Phi*P2
 P2S=P2*S;
 Leq=inv(S*B)*S*A;
 LPhi=-inv(S*B)*Phi*S;
 L=Leq+LPhi;
 rho=1.5;
 Ln=inv(S*B)*rho;
%-----
 x0=[0;th0;0;0];
 sim('CIP_smc_2015a.mdl')
%-----
%eof

Note C31 スイッチング関数決定プログラム１
Edwards and Spurgeon によって、(19)、(20)、(21)から、スイッチング関数の行列 $S$ を求めるプログラムswflqrが開発されています。

MATLAB

%test_swflqr.m
%-----
 clear all, close all
 A=[0 1 0;
    0 0 1;
    0 0 0];
 B=[0;
    0;
    1];
%-----
 Q=eye(3,3);
 S=swflqr(A,B,Q)
%-----
%eof

SM制御則

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

SM制御則…Homework

[1]　制御対象の状態方程式として次式を考えます。

$\displaystyle{(1)\quad \begin{array}{l} \dot{x}(t)=Ax(t)+Bu(t)+f(t,x,u)\\ (x(t)\in{\rm\bf R}^n, u(t)\in{\rm\bf R}^m, f(t,x,u)\in{\rm\bf R}^n, 1\le m< n, {\rm rank}B=m) \end{array} }$

ここで、 $f(t,x,u)$ はモデル誤差、非線形要素、外乱などの影響を表し、有界かつ未知とします。これに対し、次のスイッチング関数を定義します。

$\displaystyle{(2)\quad s(t)=Sx(t)\quad(s(t)\in{\rm\bf R}^m, {\rm rank}S=m) }$

以下では、状態方程式とスイッチング関数は次のSM標準形(regular form)をとるように座標変換されているとします(Note C22-1参照)。

$\displaystyle{(3)\quad \boxed{\underbrace{ \left[\begin{array}{c} \dot x_1(t)\\ \dot x_2(t) \end{array}\right] }_{\dot{x}(t)} = \underbrace{ \left[\begin{array}{cc} A_{11} & A_{12} \\ A_{21} & A_{22} \\ \end{array}\right] }_{A} \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x(t)} + \underbrace{ \left[\begin{array}{c} 0\\ B_2 \end{array}\right] }_{B} u(t) + \underbrace{ \left[\begin{array}{l} f_u(t,x)\\ f_m(t,x,u) \end{array}\right] }_{f(t,x,u)}} }$

$\displaystyle{(4)\quad \boxed{s(t)= \underbrace{ \left[\begin{array}{cc} S_1 & S_2 \\ \end{array}\right] }_{S} \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x(t)} = \underbrace{S_2 \left[\begin{array}{cc} M & I_m \\ \end{array}\right] }_{S} \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x(t)} \ (M=S_2^{-1}S_1)} }$

ここで、 $B_2$ と $S_2$ はともにサイズ $m\times m$ の正則行列とします。したがって、 $f_m$ はマッチング条件を満たしますが、 $f_u$ はマッチング条件を満たさないことに留意します。

これに対して、座標変換

$\displaystyle{(5)\quad \begin{array}{l} \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right] }_{\bar{x}(t)} = \underbrace{ \left[\begin{array}{cc} I & 0 \\ S_1 & S_2 \\ \end{array}\right] }_{T_s} \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x(t)}\\ \Leftrightarrow \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x(t)} = \underbrace{ \left[\begin{array}{cc} I & 0 \\ -S_2^{-1}S_1 & S_2^{-1} \\ \end{array}\right] }_{T_s^{-1}} \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right] }_{\bar{x}(t)} \end{array} }$

を行うために、まず(5)を(3)に代入して

$\displaystyle{(6)\quad \begin{array}{l} \underbrace{ \left[\begin{array}{cc} I & 0 \\ -M & S_2^{-1} \\ \end{array}\right] }_{T_s^{-1}} \underbrace{ \left[\begin{array}{c} \dot x_1(t)\\ \dot s(t) \end{array}\right] }_{\dot{\bar{x}}(t)}\\ = \underbrace{ \left[\begin{array}{cc} A_{11} & A_{12} \\ A_{21} & A_{22} \\ \end{array}\right] }_{A} \underbrace{ \left[\begin{array}{cc} I & 0 \\ -M & S_2^{-1} \\ \end{array}\right] }_{T_s^{-1}} \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right] }_{\bar{x}(t)} + \underbrace{ \left[\begin{array}{c} 0\\ B_2 \end{array}\right] }_{B} u(t)\\ + \underbrace{ \left[\begin{array}{l} f_u(t,x)\\ f_m(t,x,u) \end{array}\right] }_{f(t,x,u)} \end{array} }$

左から $T_s$ をかけて

$\displaystyle{(7)\quad \begin{array}{l} \underbrace{ \left[\begin{array}{c} \dot x_1(t)\\ \dot s(t) \end{array}\right] }_{\dot{\bar{x}}(t)} = \underbrace{ \left[\begin{array}{cc} I & 0 \\ S_1 & S_2 \\ \end{array}\right] }_{T_s} \underbrace{ \left[\begin{array}{cc} A_{11}-A_{12}M & A_{12}S_2^{-1} \\ A_{21}-A_{22}M & A_{22}S_2^{-1} \\ \end{array}\right] }_{AT_s^{-1}} \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right] }_{\bar{x}(t)}\\ + \underbrace{ \left[\begin{array}{cc} I & 0 \\ S_1 & S_2 \\ \end{array}\right] }_{T_s} \underbrace{ \left[\begin{array}{c} 0\\ B_2 \end{array}\right] }_{B} u(t) + \underbrace{ \left[\begin{array}{cc} I & 0 \\ S_1 & S_2 \\ \end{array}\right] }_{T_s} \underbrace{ \left[\begin{array}{l} f_u(t,x)\\ f_m(t,x,u) \end{array}\right] }_{f(t,x,u)} \end{array} }$

すなわち、次式を得ます。

$\displaystyle{(8)\quad \boxed{ \begin{array}{l} \underbrace{ \left[\begin{array}{c} \dot x_1(t)\\ \dot s(t) \end{array}\right] }_{\dot{\bar{x}}(t)} = \underbrace{ \left[\begin{array}{cc} \bar{A}_{11} & \bar{A}_{12} \\ \bar{A}_{21} & \bar{A}_{22}\\ \end{array}\right] }_{\bar{A}=T_s A T_s^{-1}} \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right] }_{\bar{x}(t)} + \underbrace{ \left[\begin{array}{c} 0\\ \bar{B}_2 \end{array}\right] }_{\bar{B}=T_sB} u(t)\\ + \underbrace{ \left[\begin{array}{c} f_u(t,x)\\ S_1f_u(t,x)+S_2f_m(t,x,u) \end{array}\right] }_{T_sf(t,x,u)} \end{array}} }$

ただし

$\displaystyle{(8')\quad \left\{\begin{array}{l} \bar{A}_{11}=A_{11}-A_{12}M\quad(M=S_2^{-1}S_1)\\ \bar{A}_{12}=A_{12}S_2^{-1}\\ \bar{A}_{21}=S_2(M\bar{A}_{11}+A_{21}-A_{22}M)\\ \bar{A}_{22}=S_2(M{A}_{12}+A_{22})S_2^{-1}\\ \bar{B}_{2}=S_2B_2 \end{array}\right. }$

以下では、 $\bar{A}_{11}$ がサイズ $n-m\times n-m$ の安定行列となるようにスイッチング関数が選ばれていると仮定します。

●このとき、スライディングモード制御則（SM制御則、SMC則）

$\displaystyle{(9)\quad u(t)=\underbrace{u_\ell(t)}_{linear\ control}+\underbrace{u_n(t)}_{switching\ component} }$

を、２次安定性

$\displaystyle{(10)\quad \boxed{ \begin{array}{lll} V(\bar{x})= \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right]^T }_{\bar{x}^T(t)} \underbrace{ \left[\begin{array}{cc} P_1 & 0\\ 0 & P_2 \end{array}\right] }_{P} \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right] }_{\bar{x}(t)}\\ \Rightarrow \dot{V}(\bar{x})\le - \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right]^T }_{\bar{x}^T(t)} \underbrace{ \left[\begin{array}{cc} Q_1 & 0\\ 0 & I \end{array}\right] }_{Q} \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right] }_{\bar{x}(t)} \end{array}}} }$

すなわち

$\displaystyle{(11)\quad \begin{array}{lll} V(\bar{x})= \underbrace{x_1^T(t)P_1x_1(t)}_{V(x_1)}+\underbrace{s^T(t)P_2s(t)}_{V(s)}\\ \Rightarrow \dot{V}(\bar{x})\le -x_1^T(t)Q_1x_1(t)-s^T(t)s(t) \end{array}} }$

が成り立つように決定します（ $P_1>0$ , $P_2>0$ , $Q_1>0$ ）。

[2] 可到達性の検討

ここでは、スライディングモード制御則(9)の具体的な表現を求め、スライディングモードに到達する時刻 $t_s$ の見積もりを行います。

●等価制御は、(7)において $f_u=0, f_m=0$ として

$\displaystyle{ \begin{array}{cl} (12.1) & s(t)=0\\ \Downarrow &\\ (12.2) & \dot{s}(t)=0\\ \Downarrow &\\ (12.3) & \bar{A}_{21}x_1(t)+\bar{A}_{22}s(t)+\bar{B}_{2}u(t)=0\\ \Downarrow &\\ (12.4) & u_{eq}(t)=-\underbrace{\bar{B}_{2}^{-1}}_{([0\ I]T_sB)^{-1}} \underbrace{\left[\begin{array}{cc} \bar{A}_{21} & \bar{A}_{22} \\ \end{array}\right]}_{[0\ I]T_sAT_s^{-1}} \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right] }_{\bar{x}(t)}\\ \Downarrow & by\ S=[0\ I]T_s, x=T_s^{-1}\bar{x}\\ (12.5) & u_{eq}(t)=-(SB)^{-1}SAx(t)} \end{array} }$

のように得られます。(9)の線形制御 $u_\ell$ は、この等価制御をベースして

$\displaystyle{(13.1)\quad u_\ell(t)=-\bar{B}_{2}^{-1} (\left[\begin{array}{cc} \bar{A}_{21} & \bar{A}_{22} \\ \end{array}\right]\bar{x}(t)-\Phi s(t)) }$

すなわち

$\displaystyle{(13.2)\quad \boxed{u_\ell(t)=-\underbrace{(SB)^{-1}(SA-\Phi S)}_{L=L_{eq}+L_{phi}}x(t)}} }$

のように構成します。ここで、 $\Phi$ は適当に設定されたサイズ $m\times m$ の安定行列です。

●このとき閉ループ系は次式で与えられます。

$\displaystyle{(14)\quad \begin{array}{l} \left[\begin{array}{c} \dot{x}_1(t)\\ \dot{s}(t) \end{array}\right] = \left[\begin{array}{cc} \bar{A}_{11} & \bar{A}_{12} \\ 0 & \Phi \\ \end{array}\right] \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right]\\ + \left[\begin{array}{c} f_u(t,x)\\ \bar{B}_{2}u_n(t)+S_2(Mf_u(t,x)+f_m(t,x,u)) \end{array}\right] \end{array} }$

すなわち

$\displaystyle{ \begin{array}{l} (15.1)\quad \dot{x}_1(t)=\bar{A}_{11}x_1(t)+\bar{A}_{12}s(t)+f_u(t,x)\\ (15.2)\quad \dot{s}(t)=\Phi s(t)+\bar{B}_{2}u_n(t)+\underline{S_2(Mf_u(t,x)+f_m(t,x,u))} \end{array} }$

閉ループ系の状態は $x_1$ と $s$ から成りますが、(15.2)は $s$ だけに関係しているので、まずその下線部に抗してスライディングモード $s=0$ を達成する $u_n$ を明らかにします。

さて、 $\Phi$ は安定行列なので

$\displaystyle{(16)\quad \begin{array}{l} \boxed{P_2\Phi+\Phi^TP_2=-I_m} \end{array} }$

を満たす $P_2>0$ を選ぶことができます。これを用いて、(7)のスイッチング項 $u_n$ を

$\displaystyle{(17.1)\quad u_n(t)=-\underbrace{\bar{B}_{2}^{-1}}_{(SB)^{-1}}\rho(t,x)\frac{P_2s(t)}{||P_2s(t)||} }$

すなわち

$\displaystyle{(17.2)\quad \boxed{u_n(t)=-\underbrace{(SB)^{-1}\rho(t,x)}_{L_n}\frac{P_2s(t)}{||P_2s(t)||}} }$

と定めます。

● $\rho(t,x)$ は、仮定

${\displaystyle{ \begin{array}{l} (18.1)\quad ||f_u(t,x)|| \le k_1||x(t)||+k_2\\ (18.2)\quad ||f_m(t,x,u)|| \le k_3||u(t)||+\alpha(t,x)\\ (18.3)\quad k_3 ||\bar{B}_{2}||\,||\bar{B}_{2}^{-1}||\,||B_2^{-1}|| < 1 \end{array} }$

のもとで

$\displaystyle{(19)\quad \boxed{ %\begin{array}{l} \rho(t,x)\ge \frac{||S_2|| (||M||(k_1||x(t)||+k_2)+k_3||u_\ell(t)||+\alpha(t,x))+\gamma_2}{1-k_3 ||\bar{B}_{2}||\,||\bar{B}_{2}^{-1}||\,||B_2^{-1}||} %\end{array} }}$

のように選びます。このとき

$\displaystyle{ \begin{array}{l} (20.1)\quad S_2=S_2B_2 B_2^{-1} \Rightarrow ||S_2||<||\bar{B}_{2}||\,||B_2^{-1}||\\ (20.2)\quad ||u_n(t)||\le\rho(t,x)||\bar{B}_{2}^{-1}||\\ (20.3)\quad ||u(t)||\le ||u_\ell(t)||+||u_n(t)|| \end{array} }$

に注意して

$\displaystyle{ \begin{array}{cl} (21.1) & \rho(t,x)\ge k_3 \rho(t,x)||\bar{B}_{2}||\,||\bar{B}_{2}^{-1}||\,||B_2^{-1}||\\ & +||S_2|| (||M||(k_1||x(t)||+k_2)+k_3||u_\ell(t)||+\alpha(t,x))+\gamma_2}\\ \Downarrow & by\ (20.1)\\ (21.2) & \ge k_3 \rho(t,x)||S_2||\,||\bar{B}_{2}^{-1}||\\ & + ||S_2|| (||M||(k_1||x(t)||+k_2)+k_3||u_\ell(t)||+\alpha(t,x))+\gamma_2}\\ \Downarrow & \\ (21.3) & = ||S_2|| \{||M||(k_1||x(t)||+k_2)\\ & +k_3(||u_\ell(t)||+ \rho(t,x)||\bar{B}_{2}^{-1}||)+\alpha(t,x)\}+\gamma_2}\\ \Downarrow & by\ (20.2) \\ (21.4) & \ge ||S_2|| (||M||(k_1||x(t)||+k_2) +k_3(||u_\ell(t)||+||u_n(t)||)+\alpha(t,x))+\gamma_2}\\ \Downarrow & by\ (20.3) \\ (21.5) & \ge ||S_2|| (||M||(k_1||x(t)||+k_2)+k_3||u(t)||+\alpha(t,x))+\gamma_2}\\ \Downarrow & by\ (18.1),(18.2) \\ (21.6) & \ge ||S_2|| (||M||\,||f_u(t,x)||+||f_m(t,x,u)||)+\gamma_2} \end{array} }$

を得て、次式が成り立ちます。(15.2)は $s$ だけに関係しているので、リャプノフ関数 $V(s)=s^T(t)P_2s(t)$ の微分について調べます。

$\displaystyle{ \begin{array}{cl} (22.1) & \dot{V}(s)=2s^T(t)P_2\dot{s}(t)\\ \Downarrow & by\ (15.2) \\ (22.2) & =2s^T(t)P_2 (\Phi s(t)-\rho(t,x)\frac{P_2s(t)} {||P_2s(t)||}+S_2(Mf_u(t,x)+f_m(t,x,u)))\\ \Downarrow & \\ (22.3) & =s^T(t)(P_2\Phi+\Phi^TP_2)s(t)\\ & +2s^T(t)P_2(-\rho(t,x)\frac{P_2s(t)} {||P_2s(t)||}+S_2(Mf_u(t,x)+f_m(t,x,u)))\\ \Downarrow & by\ (16) \\ (22.4) & \le-||s(t)||^2\\ & -2||P_2s(t)||(\rho(t,x)-||S_2||(||M||||f_u(t,x)||+||f_m(t,x,u)||)\\ \Downarrow & by\ (21.6) \\ (22.5) & \le -||s(t)||^2-2\gamma_2||P_2s(t)||\\ (22.6) & \le -||s(t)||^2 \end{array} }$

これより(15.2)における２次安定性が成り立つことが分かります。２次安定であれば状態 $s$ は0に漸近します。しかし興味深いのは、ある時刻 $t_s$ 以降は $s=0$ を保持するような振舞いをすることです。

●その到達時刻 $t_s$ を見積もるために、(22.5)から

$\displaystyle{(23)\quad \dot{V}(s)　\le -2\gamma_2||P_2s(t)|| }$

も成り立つことと

$\displaystyle{(24)\quad ||P_2s(t)||^2=(P_2^{1/2}s(t))^TP_2(P_2^{1/2}s(t)) \ge \sigma_n(P_2)||P_2^{1/2}s(t)||^2=\sigma_n(P_2)V(s) }$

に注意して

$\displaystyle{(25)\quad \dot{V}(s)\le -2\gamma_2\sigma_n^{1/2}(P_2)V^{1/2}(s) }$

よって

$\displaystyle{(26)\quad \begin{array}{l} \dot{V}(s)=\frac{d}{dt}(V^{1/2}(s))^2=2V^{1/2}(s)\frac{d}{dt}V^{1/2}(s)< -2\gamma_2\sigma_n^{1/2}(P_2) V^{1/2}(s)\\ \Rightarrow \frac{d}{dt}V^{1/2}(s)< -\gamma_2\sigma_n^{1/2}(P_2) \end{array} }$

この両辺を積分すると

$\displaystyle{(27)\quad \begin{array}{l} \underbrace{V^{1/2}(s(t_s))}_0-V^{1/2}(s(0))< -\gamma_2\sigma_n^{1/2}(P_2) t_s\\ \Leftrightarrow \boxed{t_s\le \frac{V^{1/2}(s(0))}{\gamma_2\sigma_n^{1/2}(P_2)}} \end{array} }$

を得て、 $s=0$ に到達する時刻 $t_s$ を見積もることができます。

[2] スライディングモードの検討

次に $s=0$ のもとで、(15.1)から $x_1$ の振舞いについて調べます。

$\bar{A}_{11}$ は安定行列なので

$\displaystyle{(28)\quad \begin{array}{l} \boxed{P_1\bar{A}_{11}+\bar{A}_{11}^TP_1=-Q_1<0} \end{array} }$

を満たす $P_1>0$ を選ぶことができます。そこで、(15.1)に対するリャプノフ関数 $V(s)=x_1^T(t)P_1x_1(t)$ の微分について調べます。

$\displaystyle{ \begin{array}{cl} (29.1)&\dot{V}(x_1)=2x_1^T(t)P_1\dot{x}_1(t)\\ \Downarrow & by\ (7) \\ (29.2)&=2x_1^T(t)P_1(\bar{A}_{11}x_1(t)+\bar{A}_{12}\underbrace{s(t)}_{0}+f_u(t,x))\\ \Downarrow & \\ (29.3)&=x_1^T(t)(P_1\bar{A}_{11}+\bar{A}_{11}^TP_1)x_1(t)+2x_1^T(t)P_1f_u(t,x)\\ \Downarrow & by\ (28) \\ (29.4)&=-x_1^T(t)Q_1x_1(t)+2x_1^T(t)P_1f_u(t,x)\\ \Downarrow &by\ \sigma_n(Q_1)x_1^Tx_1\le x_1^TQ_1x_1\\ (29.5)&\le -\sigma_n(Q_1)x_1^T(t)x_1(t)+2||x_1^T(t)P_1||\,||f_u(t,x)||\\ \Downarrow &\\ (29.6)&\le -\sigma_n(Q_1)||x_1(t)||^2+2\sqrt{x_1^T(t)P_1^2x_1(t)}\,||f_u(t,x)||\\ \Downarrow &by\ x_1^TP_1^2x_1\le\sigma_1(P_1^2)x_1^Tx_1\\ (29.7)&\le -\sigma_n(Q_1)||x_1(t)||^2+2\sqrt{\sigma_1(P_1^2)x_1^T(t)x_1(t)}\,||f_u(t,x)||\\ \Downarrow &by\ \sigma_1(P_1^2)=\sigma_1^2(P_1)\\ (29.8)&= -\sigma_n(Q_1)||x_1(t)||^2+2\sigma_1(P_1)||x_1(t)||\,||f_u(t,x)||\\ \Downarrow &\\ (29.9)&= -\sigma_1(P_1)||x_1(t)||(\frac{\sigma_n(Q_1)}{\sigma_1(P_1)}||x_1(t)||-2||f_u(t,x)||) \end{array} }$

ここで

$\displaystyle{(30)\quad \boxed{ ||f_u(t,x)||<\frac{1}{2}\frac{\sigma_n(Q_1)}{\sigma_1(P_1)}||x_1(t)|| }}$

のとき

$\displaystyle{(31)\quad \dot{V}(x_1)\le -\sigma_1(P_1)||x_1(t)||(\frac{\sigma_n(Q_1)}{\sigma_1(P_1)}||x_1(t)||-2||f_u(t,x)||)<0 }$

となって、リヤプノフ安定性が成り立ちます。したがって、ロバスト性での議論と同様に、 $x_1$ が発散することはないことが分かります。

Note C22-1 SM標準形について

列フルランクをもつ $B$ は、サイズ $n\times n$ の直交行列 $U$ とサイズ $m\times m$ の直交行列 $V$ を用いて

$\displaystyle{(1)\quad B=U\Sigma V^T}$

のように特異値分解できます。ただし

$\displaystyle{(2)\quad \Sigma= \left[\begin{array}{cc} \Sigma_1 \\ 0_{n-m\times m} \end{array}\right] }$

$\displaystyle{(3)\quad \Sigma_1={\rm diag}\{\sigma_1,\cdots,\sigma_m\}\quad(\sigma_1\ge\cdots\ge\sigma_m>0) }$

いま

$\displaystyle{(4)\quad U= \left[\begin{array}{cc} U_1(n\times m) & U_2(n\times n-m) \end{array}\right] }$

と分割するとき、本論(1)に対する座標変換

$\displaystyle{(5)\quad x_{r}(t)= \underbrace{ \left[\begin{array}{cc} U_2 & U_1 \end{array}\right] }_{T_r} x(t) \Leftrightarrow x(t)=T_r^Tx_r(t) }$

を行うと（ $T_r^{-1}=T_r^T$ に注意）

$\displaystyle{(6)\quad \begin{array}{l} \underbrace{ \left[\begin{array}{c} \dot{x}_{r1}(t)\\ \dot{x}_{r2}(t) \end{array}\right] }_{\dot{x}_r(t)} = \underbrace{ \left[\begin{array}{cc} A_{r11} & A_{r12} \\ A_{r21} & A_{r22} \\ \end{array}\right] }_{A_r=T_rAT_r^T} \underbrace{ \left[\begin{array}{c} x_{r1}(t)\\ x_{r2}(t) \end{array}\right] }_{x_r(t)} + \underbrace{ \left[\begin{array}{c} 0\\ \Sigma_1V^T \end{array}\right] }_{B_r=T_rB} u(t)\\ + \underbrace{ \left[\begin{array}{l} f_{ru}(t,x)\\ f_{rm}(t,x,u) \end{array}\right] }_{f_r(t,x,u)=T_rf(t,x,u)} \end{array} }$

のように本論(3)が得られます。ちなみに本論(4)は次式に相当します。

$\displaystyle{(7)\quad s(t)= \underbrace{ \left[\begin{array}{cc} S_{r1} & S_{r2} \\ \end{array}\right] }_{S_r=ST_r^T} \underbrace{ \left[\begin{array}{c} x_{r1}(t)\\ x_{r2}(t) \end{array}\right] }_{x_r(t)} }$

これから、標準形に対して求めたスイッチング関数を元の状態方程式に関するものに戻すには

$\displaystyle{(8)\quad S=S_rT_r }$

とすればよいことが分かります。

Note C22-2 条件(18)について
次式で表されるDCモータを考えます。

$\displaystyle{(1)\quad \left\{\begin{array}{l} J_0\ddot{\theta}=K_ti_a\\ L_0\dot{i}_a+Ri_a=v_a-K_e\dot{\theta} \end{array}\right. }$

これより、次の状態方程式を得ます。

$\displaystyle{(2)\quad \underbrace{ \left[\begin{array}{c} \dot{\theta}\\ \dot{\omega}\\ \dot{i}_a \end{array}\right] }_{\dot{x}} = \underbrace{ \left[\begin{array}{cc|c} 0 & 1 & 0\\ 0 & 0 & \frac{K_t}{J_0}\\\hline 0 & -\frac{K_e}{L_0} & -\frac{R}{L_0} \end{array}\right] }_{A} \underbrace{ \left[\begin{array}{c} \theta\\ \omega\\ i_a \end{array}\right] }_{x} + \underbrace{ \left[\begin{array}{c} 0\\ 0\\\hline \frac{1}{L_0} \end{array}\right] }_{B} \underbrace{ v_a }_{u} }$

またスイッチング関数として次式を考えます。

$\displaystyle{(3)\quad s= \underbrace{ \left[\begin{array}{ccc} \frac{J_0}{K_t}\omega_n^2 & 2\frac{J_0}{K_t}\zeta\omega_n & 1 \end{array}\right] }_{S}x \quad(M= \left[\begin{array}{ccc} \frac{J_0}{K_t}\omega_n^2 & 2\frac{J_0}{K_t}\zeta\omega_n \end{array}\right]) }$

さらに次式のようなパラメータの不確かさを考えます。

$\displaystyle{(4)\quad \left\{\begin{array}{l} \underbrace{(1+\xi_J)J_0}_{J}\ddot{\theta}=K_ti_a\\ \underbrace{(1+\xi_L)L_0}_{L}\dot{i}_a+Ri_a=v_a-K_e\dot{\theta} \end{array}\right. \Rightarrow \left\{\begin{array}{l} \ddot{\theta}=\frac{K_t}{J_0}i_a\underbrace{-\xi_J\ddot{\theta}}_{f_{u}}\\ \dot{i}_a=-\frac{R}{L_0}i_a-\frac{K_e}{L_0}\dot{\theta}+\frac{1}{L_0}v_a\underbrace{-\xi_L\dot{i}_a}_{f_m} \end{array}\right. }$

ただし、次を仮定します。

$\displaystyle{(5)\quad |\xi_J|<0.5,\ |\xi_L|<0.1 }$

これにより状態方程式は次式となります。

$\displaystyle{(6)\quad \dot{x}=Ax+Bu+ \left[\begin{array}{c} 0\\ f_{u}\\ f_{m} \end{array}\right]%} }$

まず $f_m$ は

$\displaystyle{(7)\quad f_m=-\xi_L\dot{i}_a=-\xi_L(-\frac{R}{L_0}i_a-\frac{K_e}{L_0}\dot{\theta}+\frac{1}{L_0}v_a) }$

と表されるので、仮定(18)の第２式における $k_3$ と $\alpha$ は次のように特定できます。

$\displaystyle{(8)\quad |f_m|<\frac{\xi_L}{L_0}|v_a|+ \frac{\xi_L}{L_0}(R|i_a|+K_e|\dot{\theta}|) <\underbrace{\frac{1}{10L_0}}_{k_3}|v_a|+ \underbrace{\frac{1}{10L_0}(R|i_a|+K_e|\dot{\theta}|)}_{\alpha} }$

また $f_u$ については、閉ループ系では

$\displaystyle{(9)\quad \begin{array}{l} f_u=-\xi_J\dot{\omega} =-\xi_J\frac{K_t}{J_0}(i_a+\frac{J_0}{K_t}\omega_n^2\theta+2\frac{J_0}{K_t}\zeta\omega_n\dot{\theta}) \end{array} }$

と考えられるので、仮定(18)の第１式における $k_1$ と $k_2$ は次のように特定できます。

$\displaystyle{(10)\quad |f_u|<\xi_J||M||||x|| <\underbrace{0.5||M||}_{k_1}||x||+\underbrace{0}_{k_2} }$

さて、(19)を見積もってみます。

$\displaystyle{(11)\quad \begin{array}{l} \rho(t,x)\ge \frac{||S_2|| (||M||(k_1||x(t)||+k_2)+k_3||u_\ell(t)||+\alpha(t,x))+\gamma_2}{1-k_3 ||\bar{B}_{2}||\,||\bar{B}_{2}^{-1}||\,||B_2^{-1}||}\\ =\frac{0.5||M||^2||x(t)||+\frac{1}{10L_0}||u_\ell(t)||+\frac{1}{10L_0}(R|i_a|+K_e|\dot{\theta}|)+\gamma_2}{1-\frac{1}{10L_0} L_0}\\ =\frac{5L_0||M||^2||x(t)||+||u_\ell(t)||+R|i_a|+K_e|\dot{\theta}|+10L_0\gamma_2}{9L_0} \end{array} }$

以上に基づいて、DCモータに対するSMCのシミュレーションを行ってみます。

MATLAB

%cDCM_smc.m
%-----
 clear all, close all
 R=1.2; L0=0.05; Ke=0.6; Kt=0.6; J=0.1352;
 A=[0 1 0;
    0 0 Kt/J;
    0 -Ke/L0 -R/L0];
 B=[0;0;1/L0];
 C=eye(3); D=zeros(3,1);
 A11=A(1:2,1:2); A12=A(1:2,3);
 A21=A(3,1:2);   A22=A(3,3);
%
 L=0.046;
 Ad=[0 1 0;
     0 0 6;
     0 -Ke/L -R/L];
 Bd=[0;0;1/L];
%-----
 wn=2; zeta=0.95;
 M=[wn^2 2*zeta*wn]/A12(2,1);
 S=[M 1];
 Leq=inv(S*B)*S*A
%
 Phi=-2
 LPhi=-inv(S*B)*Phi*S
% 
 P2=-1/(Phi*2);
 P2S=P2*S;
 gam2=0.01
 rho=(1+Ke*1+R*1+5*L0*norm(M)^2*1+10*L0*gam2)/(9*L0)
 Ln=inv(S*B)*rho
%-----
 x0=[1;0;0];
 sim('DCM_smc')
%-----
%eof

図1 DCモータに対するSMCのシミュレーション例

Note C22-4 スライディングモード制御系の安定性について

スライディングモード制御系の振舞いは、スライディングモードへの突入動作と、そこでの滑り動作から成ります。前者は、 $s=0$ を達成することですから、２次安定性

$\displaystyle{(1)\quad V(s)=s^T(t)P_2s(t)&\quad\Rightarrow \dot{V}(s)\le -s^T(t)s(t)&}\quad(t\le t_s) }$

を検討します。これが可到達条件となります。ただ、これではまだ原点 $x=0$ への収束は保証されないので、２次安定性

$\displaystyle{(2)\quad V(x_1)=x_1^T(t)P_1x_1(t)&\quad\Rightarrow \dot{V}(x_1)\le -x_1^T(t)Q_1x_1(t)&}\quad(t> t_s) }$

を検討します。

Note C22-5 安定行列 $\bar{A}_{11}=A_{11}-A_{12}M\ (M=S_2^{-1}S_1)$ について

●「 $(A,B)$ が可制御対 $\Leftrightarrow$ $(A_{11},A_{12})$ も可制御対」

実際、
$\displaystyle{(1)\quad \begin{array}{l} {\rm rank}\left[\begin{array}{cc} B & \lambda I_n-A \end{array}\right]\\ = {\rm rank}\left[\begin{array}{ccc} 0 & A_{11}-\lambda I_m & A_{12}\\ B_{2} & A_{21} & A_{22}-\lambda I_m \end{array}\right]\\ = {\rm rank}\left[\begin{array}{ccc} A_{11}-\lambda I_m & A_{12} \end{array}\right]+m \end{array} }$

これより

$\displaystyle{(2)\quad {\rm rank}\left[\begin{array}{cc} B & \lambda I_n-A \end{array}\right]=n \Leftrightarrow {\rm rank}\left[\begin{array}{ccc} A_{12} & A_{11}-\lambda I_m \end{array}\right]=n-m }$

●「 $\bar{A}_{11}$ の固有値はスイッチング関数を出力方程式としたときの不変零点に等しい」

まず、線形系

$\displaystyle{(2)\quad \begin{array}{l} \dot{x}(t)=Ax(t)+Bu(t)\\ y(t)=Cx(t)\\ (x(t)\in{\rm\bf R}^n, u(t),y(t)\in{\rm\bf R}^m) \end{array} }$

の不変零点とは、行列 $\left[\begin{array}{cc} A & B\\ C & 0 \end{array}\right]$ の一般化固有値すなわち

$\displaystyle{(3)\quad \left[\begin{array}{cc} A & B\\ C & 0 \end{array}\right] \left[\begin{array}{cc} v_\lambda\\ u_\lambda \end{array}\right] =\lambda \left[\begin{array}{cc} I_n & 0\\ 0 & 0 \end{array}\right] \left[\begin{array}{cc} v_\lambda\\ u_\lambda \end{array}\right]\quad (\left[\begin{array}{cc} v_\lambda\\ u_\lambda \end{array}\right]\ne0) }$

を満足する $\lambda$ として定義されます。これから次式が成り立ちます。

$\displaystyle{(4)\quad \left\{\begin{array}{l} v_\lambda=(\lambda I_n-A)^{-1}Bu_\lambda\\ Cv_\lambda=0 \end{array}\right. \Rightarrow C(\lambda I_n-A)^{-1}Bu_\lambda=0 }$

詳しい説明は省きますが、これは指数関数入力 $e^{\lambda t}u_\lambda$ に対する応答が零になることを意味しています。

(3)が非自明な解をもつためには次式が必要十分となります。

$\displaystyle{(5)\quad {\rm det}\left[\begin{array}{cc} \lambda I_n-A & -B\\ -C & 0 \end{array}\right]=0 }$

さて、次の線形系を考えます。

$\displaystyle{(6)\quad \begin{array}{l} \dot{x}(t)=Ax(t)+Bu(t)\\ s(t)=Sx(t)\\ (x(t)\in{\rm\bf R}^n, u(t),s(t)\in{\rm\bf R}^m) \end{array} }$

この不変零点を求めます。

$\displaystyle{(7)\quad \begin{array}{l} {\rm det}\left[\begin{array}{cc} \lambda I_n-A & -B\\ -S & 0 \end{array}\right]\\ = {\rm det}\left[\begin{array}{ccc} \lambda I_{n-m}-A_{11} & -A_{12} & 0\\ -A_{21} & \lambda I_m-A_{22} & -B_{2}\\ -S_1 & -S_2 & 0 \end{array}\right]=0\\ \Leftrightarrow {\rm det}\left[\begin{array}{ccc} \lambda I_{n-m}-A_{11} & -A_{12} \\ -S_1 & -S_2 \end{array}\right]=0\\ \Leftrightarrow{\rm det} \left[\begin{array}{cc} I_{n-m} & -A_{12}S_2^{-1} \\ 0 & I_m \end{array}\right]\left[\begin{array}{cc} \lambda I_{n-m}-\bar{A}_{11} & 0 \\ 0 & -S_2 \end{array}\right] \left[\begin{array}{cc} I_{n-m} & 0 \\ M & I_m \end{array}\right]\\ ={\rm det}(\lambda I_{n-m}-\bar{A}_{11}){\rm det}(-S_2)=0\\ \Leftrightarrow {\rm det}(\lambda I_{n-m}-\bar{A}_{11})=0 \end{array} }$

すなわち、スイッチング関数を出力方程式とするシステム(6)の不変零点は、行列 $\bar{A}_{11}$ の固有値と一致します。この行列 $\bar{A}_{11}$ はスライディングモードの振舞いを決めるので、 $\bar{A}_{11}$ は安定行列とすることが前提となります。したがって、どのようなスイッチング関数が望ましいかは、どのような超平面上をどのようにスライディングさせるかに関わっていると言えます。

問題設定

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

問題設定…Homework

[1] 制御対象の状態方程式として次式を考えます。

$\displaystyle{(1)\quad \begin{array}{l} \dot{x}(t)=Ax(t)+Bu(t)+f(t,x,u)\\ (x(t)\in{\rm\bf R}^n, u(t)\in{\rm\bf R}^m, f(t,x,u)\in{\rm\bf R}^n, 1\le m< n, {\rm rank}B=m) \end{array} }$

ここで、 $f(t,x,u)$ はモデル誤差、非線形要素、外乱などの影響を表し、有界かつ未知とします。これに対し、次のスイッチング関数を定義します。

$\displaystyle{(2)\quad \boxed{s(t)=Sx(t)}\quad(s(t)\in{\rm\bf R}^m, {\rm rank}S=m) }$

以下では、(1)の解が、ある時刻 $t_s$ に対して

$\displaystyle{(3)\quad s(t)=0\quad(t\ge t_s) }$

を満足するような状況、すなわちスライディングモードを考えます。

以上の準備の下で、スライディングモード制御問題は次のように記述されます。

問題１: 動的システムの振舞いを閉じ込める超平面 ${\cal S}=\{x:Sx=0\}$ を定義するスイッチング関数を決定せよ。
問題２: 有限時刻 $t_s$ において状態を超平面内に拘束し、引き続き留まらせる( $\forall t>t_s: s(t)=0$ )ことのできるスライディングモード制御則（SMC則）を設計せよ。

この問題の解、SMC則は、次のように表されます。

$\displaystyle{(4)\quad { u(t)=\underbrace{u_\ell(t)}_{linear\ control}+\underbrace{u_n(t)}_{switching\ component}} }$

$\displaystyle{(5)\quad { \boxed{u_\ell(t)=-\underbrace{(SB)^{-1}(SA-\Phi S)}_{L=L_{eq}+L_{phi}}x(t)}} }$

$\displaystyle{(6)\quad { \boxed{u_n(t)=-\underbrace{(SB)^{-1}\rho(t,x)}_{L_n}\frac{P_2s(t)}{||P_2s(t)||}}} }$

ここで、 $SB$ は正則行列であることを仮定しています。また設計パラメータは、安定行列 $\Phi$ 、スカラー $\rho>0$ です。 $P_2$ は次式から定まる正定行列です。

$\displaystyle{(7)\quad P_2\Phi+\Phi^TP_2=-I_m }$

以下ではこの解の妥当性を２次安定化の観点から検討します。

[2]　いま(1)において $f(t,x,u)=0$ 、すなわち

$\displaystyle{(8)\quad \dot{x}(t)=Ax(t)+Bu(t) }$

とし、 $SB$ の正則性を仮定します。 $t\ge t_s$ において、スライディングモード時の制御則は等価制御と呼ばれ、次式のように求められます。

$\displaystyle{(9)\quad \begin{array}{l} s(t)=0\\ \Rightarrow\dot{s}(t)=0\\ \Rightarrow S\dot{x}(t)=SAx(t)+SBu(t)=0\\ \Rightarrow u_{eq}(t)=-\underbrace{(SB)^{-1}SA}_{L_{eq}}x(t) \end{array} }$

これによる閉ループ系は次式となります。

$\displaystyle{(10)\quad \boxed{\dot{x}(t)=\underbrace{(I_n-B(SB)^{-1}S)A}_{A_{eq}}x(t)}\quad(t\ge t_s) }$

ちなみに、(4)の第１項 $u_\ell$ は、この等価制御をベースして

$\displaystyle{(11)\quad u_\ell(t)=-\underbrace{(SB)^{-1}SA}_{L_{eq}}x(t)-\underbrace{(SB)^{-1}(-\Phi) S}_{L_{phi}}x(t) }$

のように構成されています。この $\Phi$ の役割についてはあとで述べます。

[3]　いま(1)において $f(t,x,u)\ne 0$ の場合、適合サイズをもつ行列 $R$ を用いて

$\displaystyle{(12)\quad　{ \boxed{f(t,x,u)=BR\xi(t,x)}} }$

のように表されると仮定します。このとき、等価制御は

$\displaystyle{(13)\quad \begin{array}{l} s(t)=0\\ \Rightarrow\dot{s}(t)=0\\ \Rightarrow S\dot{x}(t)=SAx(t)+SBu(t)+SBR\xi(t,x)=0\\ \Rightarrow u(t)=-(SB)^{-1}(SAx(t)+SBR\xi(t,x)) \end{array} }$

となります。これによる閉ループ系は次式となります。

$\displaystyle{(14)\quad \begin{array}{l} \dot{x}(t)=Ax(t)-B(SB)^{-1}(SAx(t)+SBR\xi(t,x))+BR\xi(t,x)\\ =\underbrace{(I_n-B(SB)^{-1}S)A}_{A_{eq}}x(t)+ \underbrace{(I_n-B(SB)^{-1}S)B}_0R\xi(t,x)\\ =\underbrace{(I_n-B(SB)^{-1}S)A}_{A_{eq}}x(t) \end{array} }$

すなわちマッチング条件とよばれる(12)が成り立つとき、スライディングモード時は $f(t,x,u)$ の影響を受けないことが分かります。

演習…Flipped Classroom
$1^\circ$ 次の剛体振子の状態方程式

$\displaystyle{ \left[\begin{array}{c} \dot{x}_1(t)\\ \dot{x}_2(t) \end{array}\right]= \underbrace{ \left[\begin{array}{cc} 0 & 1\\ 0 & 0 \end{array}\right] }_{A} \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right]+ \underbrace{ \left[\begin{array}{c} 0\\ 1 \end{array}\right] }_{B} u(t)+ \underbrace{ \left[\begin{array}{c} 0\\ -a\sin x_1(t) \end{array}\right] }_{f} }$

に対して、スイッチング関数を

$\displaystyle{ s(t)= \underbrace{ \left[\begin{array}{cc} m & 1 \end{array}\right] }_{S} \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }$

と選ぶとき、SMC則の線形制御の部分は次式となることを示せ。

$\displaystyle{ u_\ell(t)=- \left[\begin{array}{cc} -m\Phi & m-\Phi \end{array}\right] \left[\begin{array}{cc} x_1(t)\\ x_2(t) \end{array}\right] }$

※ $-\Phi$ を改めて $\Phi$ とおけば、可到達条件のところで述べた次式(14)となります。

$\displaystyle{(14)\quad u(t)=- \left[\begin{array}{cc} m\Phi & m+\Phi \end{array}\right] \left[\begin{array}{cc} x_1(t)\\ x_2(t) \end{array}\right] -\rho\frac{s(t)}{|s(t)|} }$

$2^\circ$ 上の剛体振子の状態方程式についてマッチング条件(12)を調べよ。

Note C21１入力系の場合

●ここでは次の１入力系を考えます。

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

$\displaystyle{(2)\quad A= \left[\begin{array}{cccc|c} 0 & 1 & 0 & \cdots & 0 \\ 0 & 0 & 1 & 0 & \vdots \\ \vdots & & \ddots & \ddots & 0 \\ 0 & \cdots & \cdots & 0 & 1 \\\hline -a_1 & -a_2 & \cdots & \cdots & -a_n \end{array}\right],\quad B=\left[\begin{array}{c} 0 \\ 0 \\ \vdots \\ 0 \\\hline 1 \end{array}\right] }$

$\displaystyle{(3)\quad \begin{array}{r} \Delta A(t)= \left[\begin{array}{cccc|c} 0 & 0 & 0 & \cdots & 0 \\ 0 & 0 & 0 & 0 & \vdots \\ \vdots & & \ddots & \ddots & 0 \\ 0 & \cdots & \cdots & 0 & 0 \\\hline -\delta_1(t) & -\delta_2(t) & \cdots & \cdots & -\delta_n(t) \end{array}\right]\\ k_i^-\le\delta_i(t)\le k_i^+\quad(i=1,\cdots,n) \end{array} }$

また次のスイッチング関数を考えます。

$\displaystyle{(4)\quad s(t)=Sx(t) }$

$\displaystyle{(5)\quad S=\left[\begin{array}{ccccc} m_1 & m_2 & \cdots & m_{n-1} & 1 \end{array}\right] }$

このとき、SMC則

$\displaystyle{(6)\quad u(t)=\underbrace{u_\ell(t)}_{linear\ control}+\underbrace{u_n(t)}_{switching\ component} }$

を可到達条件を満たすように構成することを考えます。

●線形制御 $u_\ell$ は、

$(7)\quad \begin{array}{l} \displaystyle{\dot{s}(t)=\sum_{i=1}^{n-1}m_i\dot{x}_i(t)+\dot{x}_n(t)}\\ \displaystyle{=\sum_{i=1}^{n-1}m_i x_{i+1}(t)-\sum_{i=1}^{n}(a_i+\delta_i(t))x_i(t)+u(t)}\\ \displaystyle{=-a_1x_1(t)+\sum_{i=2}^{n}(m_{i-1}-a_i) x_i(t)- \sum_{i=1}^{n}\delta_i(t)x_i(t)+u(t)} \end{array}$

において、 $\delta_i=0$ とおいて

$\displaystyle{(8)\quad \boxed{u_\ell(t)=a_1x_1(t)+\sum_{i=2}^{n}(a_i-m_{i-1}) x_i(t)} }$

のように定めます。

● $u_n$ の候補として次の２つを考えます。

$\displaystyle{(9)\quad \boxed{u_n(t)=-\rho(t,x)\frac{s(t)}{|s(t)|}} }$

$\displaystyle{(10)\quad \boxed{u_n(t)=\sum_{i=1}^{n}k_ix_i(t)-\eta\frac{s(t)}{|s(t)|}} }$

以下では、これらが適当な条件の下で可到達条件 $s\dot{s}<0$ を満足しSMC則になることを示します。

(9)の場合：

$\displaystyle{(11)\quad \begin{array}{l} \rho(t,x)=\sum_{i=1}^{n}\max\{|k_i^+|,|k_i^-|\}|x_i(t)|+\eta\\ \ge |\sum_{i=1}^{n}\delta_i(t)x_i(t)|+\eta \end{array} }$

を満足する $\rho(t,x)$ を選ぶと

$\displaystyle{(12)\quad \begin{array}{l} s(t)\dot{s}(t)=s(t)(-\sum_{i=1}^{n}\delta_i(t)x_i(t)+u_n(t))\\ =s(t)(-\sum_{i=1}^{n}\delta_i(t)x_i(t)-\rho(t,x)\frac{s(t)}{|s(t)|})\\ =-s(t)\sum_{i=1}^{n}\delta_i(t)x_i(t)-|s(t)|\rho(t,x)\\ \le |s(t)|(|\sum_{i=1}^{n}\delta_i(t)x_i(t)|-\rho(t,x))\\ < -\eta|s(t)| \end{array} }$

(10)の場合：

$\displaystyle{(13)\quad k_i= \left\{\begin{array}{ll} k_i^- & (sx_i>0)\\ k_i^+ & (sx_i<0) \end{array}\right. }$

のように $k_i$ を選ぶと

$\displaystyle{(14)\quad \begin{array}{l} s(t)\dot{s}(t)=s(t)(-\sum_{i=1}^{n}\delta_i(t)x_i(t)+u_n(t))\\ =s(t)(-\sum_{i=1}^{n}\delta_i(t)x_i(t)+\sum_{i=1}^{n}k_ix_i(t)-\eta\frac{s(t)}{|s(t)|})\\ =\sum_{i=1}^{n}\underbrace{(k_i-\delta_i(t))s(t)x_i(t)}_{<0}-\eta|s(t)|\\ < -\eta|s(t)| \end{array} }$

ロバスト性

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

線形系の２次安定性…Homework

対象が平衡状態にあることは線形状態方程式 $\dot{x}(t)=Ax(t)+Bu(t)$ において、 $x=0, u=0$ を意味します。そこで平衡状態が乱されて $x\ne0$ となる時刻を $t=0$ にとると、線形状態方程式は次式となります。

$\displaystyle{(1)\quad \dot{x}(t)=Ax(t)\qquad(x(0)\ne0) }$

これに出力方程式

$\displaystyle{(2)\quad y(t)=Cx(t) }$

を考慮する場合も含めると、漸近安定性の判定法は次のようにまとめられます。

【漸近安定性の定義とその等価な条件】

定義DA:　 $\forall x(0)\ne 0: x(t)=\exp(At)x(0)\rightarrow 0\quad(t\rightarrow\infty)$

条件A0:　 $\exp(At)\rightarrow 0\quad(t\rightarrow\infty)$

条件A1:　 ${\rm Re}(\lambda_i(A))<0\quad(i=1,\cdots,n)$

条件A2:　 $\exists P>0: PA+A^TP+I_n=0$

条件A3:　 $\exists P>0: PA+A^TP+C^TC=0$ 　ただし、 $(A,C)$ は可観測対

条件A4:　 $\exists P>0: PA+A^TP<0$

条件A2のリャプノフ方程式から次式を得ます。

$\displaystyle{(3)\quad x^T(t)(PA+A^TP)x(t)=-x(t)^Tx(t) }$

線形系(1)に対してリャプノフ関数

$\displaystyle{(4)\quad V(x)=x^T(t)Px(t) }$

を考えるとき、(3)は２次安定性の条件

$\displaystyle{(5)\quad \frac{d}{dt}V(x)=2x^T(t)P\dot{x}(t)=2x^T(t)PAx(t)=-x^T(t)Qx(t)\quad (Q=I) }$

を意味します。これから漸近安定な線形系(1)は２次安定であることが分かります。

●条件A2の両辺から $Q^{1/2}\ (Q>0)$ をかけて

$\displaystyle{(6)\quad \underbrace{Q^{1/2}PQ^{1/2}}_{P'}\underbrace{Q^{-1/2}AQ^{1/2}}_{A'}+\underbrace{Q^{1/2}A^TQ^{-1/2}}_{A'^T}\underbrace{Q^{1/2}PQ^{1/2}}_{P'}+\underbrace{Q^{1/2}Q^{1/2}}_{Q}=0 }$

を得ます。 $(A',Q^{1/2})$ は可観測対なので条件A3を適用して（ $'$ をはずして）

条件A2′:　 $\boxed{\exists P>0: PA+A^TP+Q=0}$ 　ただし、 $Q>0$

も漸近安定性の条件となることに注意します。

モデルの不確かさがある場合…Homework

いま漸近安定な線形系(1)において次のようなモデルの不確かさがあるとします。

$\displaystyle{(7)\quad \dot{x}(t)=(A+\Delta A)x(t)=Ax(t)+\underbrace{\Delta Ax(t)}_{\xi(t,x)}\qquad(x(0)\ne0) }$

条件A2’より、適当な $P>0$ と $Q>0$ に対して

$\displaystyle{(8)\quad PA+A^TP+Q=0 }$

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

$\displaystyle{(9)\quad \boxed{||\xi(t,x)||\le\frac{1}{2}\frac{\sigma_n(Q)}{\sigma_1(P)}||x(t)||} }$

ならば、リャプノフ関数

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

に対して

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

となって、リャプノフ安定となります。

●実際、

$\displaystyle{(12)\quad \begin{array}{l} \frac{d}{dt}V(x)=x^T(t)P\dot{x}(t)+\dot{x}^T(t)Px(t)\\ =x^T(t)P(Ax(t)+\xi(t,x))+(x^T(t)A^T+\xi^T(t,x))Px(t)\\ =x^T(t)PAx(t)+x^T(t)A^TPx(t)+x^T(t)P\xi(t,x)+\xi^T(t,x)Px(t)\\ =x^T(t)(PA+A^TP)x(t)+2\xi^T(t,x)Px(t)\\ =-x^T(t)Qx(t)+2\xi^T(t,x)Px(t)\\ \le -x^T(t)Qx(t)+ 2||\xi(t,x)||\underbrace{||Px(t)||}_{\sqrt{x^T(t)P^2x(t)}} \end{array} }$
ここで、

$\displaystyle{(13)\quad \begin{array}{l} \sigma_n(Q)x^T(t)x(t)\le x^T(t)Qx(t) \Rightarrow \sigma_n(Q)||x(t)||^2\le x^T(t)Qx(t)\\ x^T(t)P^2x(t)\le\sigma_1(P^2)x(t)^Tx(t) \Rightarrow \sqrt{x^T(t)P^2x(t)}\le\sigma_1(P)||x(t)|| \end{array} }$

に注意して、(9)の下では

$\displaystyle{(14)\quad \begin{array}{l} \frac{d}{dt}V(x) \le -\sigma_n(Q)||x(t)||^2+ 2||\xi(t,x)||\sigma_1(P)||x(t)|| \\ =(||\xi(t,x)||-\frac{1}{2}\frac{\sigma_n(Q)}{\sigma_1(P)}||x(t)||)2\sigma_1(P)||x(t)||<0 \end{array} }$

となります。

●次の数値例を考えてみます。

$\displaystyle{(15)\quad \left[\begin{array}{c} \dot{\theta}(t) \\ \dot{\omega}(t) \end{array}\right] = \underbrace{ \left[\begin{array}{cc} 0 & 1 \\ -4 & -1 \end{array}\right] }_{A} \left[\begin{array}{c} \theta(t) \\ \omega(t) \end{array}\right] + \underbrace{\left[\begin{array}{c} 0 \\ 1 \end{array}\right]}_{D}\underbrace{\sin(2t)}_{\xi(t)} }$

$A$ は安定行列なので、リャプノフ方程式 $PA+A^TP+I_n=0$ を満足する $P>0$ は

$\displaystyle{(16)\quad P=\left[\begin{array}{cc} 2.6250 & 0.1250\\ 0.1250 & 0.6250 \end{array}\right] }$

のように求まります。

この $P$ をもつリャプノフ関数(10)に対して

$\displaystyle{(17)\quad \begin{array}{l} \frac{d}{dt}V(x)=x^T(t)P\dot{x}(t)+\dot{x}^T(t)Px(t)\\ =x^T(t)P(Ax(t)+D\xi(t))+(x^T(t)A^T+\xi^T(t)D^T)Px(t)\\ =x^T(t)PAx(t)+x^T(t)A^TPx(t)+x^T(t)PD\xi(t)+\xi(t)D^TPx(t)\\ =x^T(t)(PA+A^TP)x(t)+2\xi(t)D^TPx(t)\\ =-x^T(t)Qx(t)+2\xi(t)D^TPx(t)\\ \le -x^T(t)x(t)+ 2|\xi(t)|||D^TP||||x(t)||\\ \le -||x(t)||^2+ 2\sigma_1(PD)||x(t)|| \\ =(2\sigma_1(PD)-||x(t)||)||x(t)||<0 \end{array} }$

これから

$\displaystyle{(18)\quad ||x(t)||>2\sigma_1(PD)=1.2748 }$

のときは $\frac{d}{dt}V(x)<0$ となって原点に漸近していきますが、 $||x(t)||<1.2748$ になるとそれは保証されないことが分かります(図1参照)。

図1 モデルの不確かさがある場合の解軌道

演習C12…Flipped Classroom
$1^\circ$ 　smc17.slxを用いて、図1のシミュレーションを行え。

Note C12 振り子の安定性

振り子の運動方程式

$\displaystyle{(1)\quad \underbrace{\frac{4}{3}m\ell^2}_J \dot{\omega}(t)=-mg\ell\sin\theta(t)-c\,\omega(t)\quad(\omega(t)=\dot{\theta}(t)) }$

から得られる非線形状態方程式

$\displaystyle{(2)\quad \frac{d}{dt} \underbrace{\left[\begin{array}{c} \theta(t) \\ \omega(t) \end{array}\right]}_{x} = \underbrace{ \left[\begin{array}{c} \omega(t) \\ -\frac{mg\ell}{J}\sin\theta(t)-\frac{c}{J}\omega(t) \end{array}\right] }_{f(x)} }$

に対して、次のリャプノフ関数の候補を考えます。

$\displaystyle{(3)\quad V=\frac{1}{2} \underbrace{ \left[\begin{array}{c} \theta(t) \\ \omega(t) \end{array}\right]^T \left[\begin{array}{cc} \frac{1}{2}\frac{c^2}{J} & \frac{1}{2}c \\ \frac{1}{2}c & J \end{array}\right] \left[\begin{array}{c} \theta(t) \\ \omega(t) \end{array}\right] }_{x^T(t)Px(t)} +mg\ell(1-\cos\theta(t))>0 }$

この時間微分を計算すると

$\displaystyle{(4)\quad \begin{array}{l} \dot{V}(t)=\frac{1}{2}(x^TP\dot{x}+\dot{x}^TPx)+mg\ell\omega\sin\theta =x^TP\dot{x}+mg\ell\omega\sin\theta\\ =\left[\begin{array}{c} \theta \\ \omega \end{array}\right]^T \left[\begin{array}{cc} p_1 & p_3 \\ p_3 & p_2 \end{array}\right] \left[\begin{array}{c} \omega \\ -\frac{mg\ell}{J}\sin\theta-\frac{c}{J}\omega \end{array}\right]+mg\ell\omega\sin\theta\\ =\left[\begin{array}{cc} \theta &\omega \end{array}\right] \left[\begin{array}{cc} p_1\omega -p_3(\frac{mg\ell}{J}\sin\theta+\frac{c}{J}\omega) \\ p_3\omega -p_2(\frac{mg\ell}{J}\sin\theta+\frac{c}{J}\omega) \end{array}\right]+mg\ell\omega\sin\theta\\ =p_1\theta\omega -p_3(\frac{mg\ell}{J}\theta\sin\theta+\frac{c}{J}\theta\omega)+ p_3\omega^2 -p_2(\frac{mg\ell}{J}\omega\sin\theta+\frac{c}{J}\omega^2)+mg\ell\omega\sin\theta\\ =-\frac{mg\ell}{J}\sin\theta(p_3\theta+p_2\omega-J\omega) +(p_3-p_2\frac{c}{J})\omega^2+(p_1-p_3\frac{c}{J})\theta\omega \end{array} }$

ここで、 $p_1=\frac{1}{2}\frac{c^2}{J}$ 、 $p_2=J$ 、 $p_3=\frac{1}{2}c$ を代入して

$\displaystyle{(5)\quad \dot{V}(t)=-\frac{1}{2}\frac{mg\ell}{J}c\,\theta(t)\sin\theta(t) -\frac{1}{2}c\,\omega^2(t) }$

となります。ここで

$\displaystyle{(6)\quad \theta(t)\sin\theta(t)>0\quad (|\theta(t)|<\pi,\theta(t)\ne0) }$

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

●ちなみに(2)を次のように書き直してみます。

$\displaystyle{(7)\quad \left[\begin{array}{c} \dot{\theta}(t) \\ \dot{\omega}(t) \end{array}\right] = \underbrace{ \left[\begin{array}{cc} 0 & 1 \\ -\frac{mg\ell}{J} & -\frac{c}{J} \end{array}\right] }_{A} \left[\begin{array}{c} \theta(t) \\ \omega(t) \end{array}\right] + \underbrace{ \left[\begin{array}{c} 0 \\ \frac{mg\ell}{J} \end{array}\right] }_{D} \underbrace{ (\theta(t)-\sin\theta(t)) }_{\xi(x)} }$

我々は軸摩擦をもつ振り子の角度は十分時間がたてば零に収束することを知っています。上の議論に沿って解釈すれば、モデルの不確かさ $\xi(x)$ 自体が限りなく零に近づき（ $\sin\theta\simeq\theta$ ）、線形系となると考えられます。

Note C12 　MATLABによるリャプノフ方程式の求解

たとえば、安定行列 $\displaystyle{ A= \left[\begin{array}{cc} 0 & 1 \\ -4 & -1 \end{array}\right] }$

が与えられるとき、リャプノフ方程式 $PA+A^TP+I_n=0$ を満足する $P>0$ は

$\displaystyle{ P=\left[\begin{array}{cc} 2.6250 & 0.1250\\ 0.1250 & 0.6250 \end{array}\right] }$

のように求まります。これはMATLABで次のコマンドを用いました。

P=lyap(A’,eye(2))

ここで、P=lyap(A,Q) はリャプノフ方程式 $PA^T+AP+Q=0$ の解を求めることに注意します。この解を $(A,Q)$ に対するリャプノフ行列と呼ぶことがあります。

２次安定性

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

リャプノフの安定定理

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

リャプノフの安定定理　状態方程式

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

に対して、その平衡状態 $x^*=0$ を含む領域 ${\cal D}$ を考える。いま ${\cal D}$ で連続微分可能な実数値関数 $V$ に対して

$\displaystyle{(2)\quad \left\{\begin{array}{l} V(0)=0 \\ V(x)>0\quad(x\in{\cal D},x\not=0) \end{array}\right.}$

このとき

$\displaystyle{(3)\quad \frac{d}{dt}V(x)\le0 \quad(x\in{\cal D})}$

ならば、平衡状態 $x^*=0$ は安定、また

$\displaystyle{(4)\quad \frac{d}{dt}V(x)<0 \quad(x\in{\cal D},x\not=0)}$

ならば、平衡状態 $x^*=0$ は漸近安定である。

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

２次安定性

非線形系(1)に対してリャプノフ関数

$\displaystyle{(5)\quad V(x)=x^T(t)Px(t)\quad(P>0)}$

を考えます。これを微分して(1)を用いると（(1)の解に沿って微分して）

$\displaystyle{(6)\quad %\begin{array}{l} \frac{d}{dt}V(x)=\dot{x}^T(t)Px(t)+x^T(t)P\dot{x}(t)=2x^T(t)P\dot{x}(t)=2x^T(t)Pf(x) %\end{array} }$

となります。このとき２次安定とは

$\displaystyle{(7)\quad { \boxed{\frac{d}{dt}V(x)=2x^T(t)Pf(x)\le -x^T(t)Qx(t)\quad(Q>0)} }}$

が成り立つことを言います。これは、次式を示すことができ、漸近安定性を意味するからです。

$\displaystyle{(8)\quad \boxed{||x(t)||\le \beta||x(0)||e^{-\frac{1}{2}\alpha t}} }$

ただし

$\displaystyle{(9)\quad \boxed{\alpha=\frac{\sigma_n(Q)}{\sigma_1(P)},\ \beta=\sqrt{\frac{\sigma_1(P)}{\sigma_n(P)}} %\alpha=\frac{\sigma_n(Q)}{\sigma_1(P)}\le\sigma_n(P^{-1}Q),\ }$

以下では、これを導出します。

●予備知識として、任意の対称行列 $R>0$ に対して、特異値（固有値に等しい）を

$\displaystyle{(10)\quad \sigma_1(R)\ge\cdots\ge\sigma_n(R) }$

と表記し、また任意の $v\ne0$ に対して

$\displaystyle{(11)\quad \sigma_n(R)\le\frac{v^TRv}{v^Tv}\le\sigma_1(R) }$

および

$\displaystyle{(12)\quad \sigma_1(R^{-1})=\frac{1}{\sigma_n(R)} }$

が成り立つこと使います。

●そこで

$\displaystyle{(13)\quad \begin{array}{l} \sigma_n(Q)x^T(t)x(t)\le x^T(t)Qx(t)\\ x^T(t)Px(t)\le\sigma_1(P)x(t)^Tx(t) \end{array} }$

に注意して、次式を得ます。

$\displaystyle{(14)\quad \frac{d}{dt}V(x)\le -x^T(t)Qx(t) \le -\sigma_n(Q)x^T(t)x(t) \le -\sigma_n(Q) \frac{x^T(t)Px(t)}{\sigma_1(P)} \le -\alpha V(x) }$

これから

$\displaystyle{(15)\quad V(x)\le V(0)e^{-\alpha t} \Leftrightarrow x^T(t)Px(t)\le x^T(0)Px(0)e^{-\alpha t} }$

ここで

$\displaystyle{(16)\quad {\begin{array}{l} %\sigma_n^2(P)x^Tx\le x^TPx\le\sigma_1^2(P)x^Tx\\ %\sigma_n^2(P)x^T(0)x(0)\le x^T(0)Px(0)\le\sigma_1^2(P)x(0)^Tx(0) \sigma_n(P)x^T(t)x(t)\le x^T(t)Px(t)\\ x^T(0)Px(0)\le\sigma_1(P)x(0)^Tx(0) \end{array} }}$

に注意して

$\displaystyle{(17)\quad {\sigma_n(P)x^T(t)x(t) \le \sigma_1(P) x^T(0)x(0)e^{-\alpha t} }}$

すなわち

$\displaystyle{(18)\quad {\underbrace{x^T(t)x(t)}_{||x(t)||^2} \le \underbrace{\frac{\sigma_1(P)}{\sigma_n(P)}}_{\beta^2} \underbrace{x(0)^Tx(0)}_{||x(0)||^2} e^{-\alpha t} }}$

を得ます。これから(8)が成り立ちます。なお一般には

$\displaystyle{(19)\quad \underbrace{\frac{\sigma_n(Q)}{\sigma_1(P)}}_{\alpha} =\frac{1}{\sigma_1(Q^{-1})\sigma_1(P)} \le \frac{1}{\sigma_1(Q^{-1}P)} =\underbrace{\sigma_n(P^{-1}Q)}_{\alpha'} }$

ですが、 $P$ と $Q$ を対角行列に選ぶときは等号が成り立ち、(8)において $\alpha$ ではなく $\alpha'$ を用いることができます。もう少し詳しくは、対称行列 $X=P$ と $Y=Q^{-1}$ が可換（ $XY=YX$ ）で、同時対角化可能（ $TXT^{-1}$ と $TYT^{-1}$ を対角行列とする $T$ が存在）の場合には等号が成り立ち、(8)において $\alpha$ ではなく $\alpha'$ を用いることができます。

Note C11 (11),(12)の導出

●任意の正定行列 $R>0$ （サイズ $n\times n$ ）に対して、次の特異値分解を考えます。

$\displaystyle{(1)\quad \begin{array}{l} R=U\Sigma U^T\\ \Sigma={\rm diag}\{\sigma_1,\cdots,\sigma_n\}\quad(\sigma_1\ge\cdots\ge\sigma_n)\\ UU^T=U^TU=I_n \end{array} }$

ここで、 $R$ は対称行列なので固有値は実数、しかも正定なのでそれらはすべて正であること、したがって特異値 $\sigma_1$ は $R$ の最大固有値、特異値 $\sigma_n$ は $R$ の最小固有値であることに注意してください。このとき、 $R=R^{1/2}R^{1/2}$ を満たす $R$ の平方根行列は

$\displaystyle{(2)\quad \begin{array}{l} R^{1/2}=U\Sigma^{1/2} U^T\\ \Sigma^{1/2}={\rm diag}\{\sqrt{\sigma_1},\cdots,\sqrt{\sigma_n}\} \end{array} }$

で与えられます。また、 $R$ の逆行列は

$\displaystyle{(3)\quad \begin{array}{l} R^{-1}=U\Sigma^{-1} U^T\\ \Sigma^{-1}={\rm diag}\{\frac{1}{\sigma_1},\cdots,\frac{1}{\sigma_n}\} \quad(\frac{1}{\sigma_n}\ge\cdots\ge\frac{1}{\sigma_1}) \end{array} }$

のように表されます。これから(12)は自明です。

いま、 $R^{1/2}$ を用いて線形写像 $y=R^{1/2}x$ を考えると、これは３つの線形写像に分解されます。

$\displaystyle{(4)\quad \begin{array}{l} x'=U^Tx\\ y'=\Sigma^{1/2} x'\\ y=Uy' \end{array} }$

任意の $n$ 次元ベクトル $x\in{\rm\bf R}^n$ のノルムとして、次の２ノルムを考えます。

$\displaystyle{(5)\quad ||x||_2=\sqrt{x_1^2+\cdots+x_n^2}&\Rightarrow&||x||_2^2=x^Tx\\ }$

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

$\displaystyle{(6)\quad ||x'||_2^2=||U^Tx||_2^2=x^TUU^Tx=x^Tx=||x||_2^2\ \Rightarrow\ \frac{||x'||_2}{||x||_2}=1 }$

$\displaystyle{(7)\quad \begin{array}{l} ||y'||_2^2=||\Sigma^{1/2} x'||_2^2=\sigma_1x_1'\,^2+\cdots+\sigma_nx_n'\,^2\\ \le \sigma_1(x_1'\,^2+\cdots+x_n'\,^2)=\sigma_1x'\,^Tx'=\sigma_1||x'||_2^2 \ \Rightarrow\ \frac{||y'||_2}{||x'||_2}\le\sqrt{\sigma_1} \end{array} }$

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

すなわち

$\displaystyle{(9)\quad \frac{||x'||_2}{||x||_2}\times\frac{||y'||_2}{||x'||_2}\times\frac{||y||_2}{||y'||_2}\le\sqrt{\sigma_1} \ \Rightarrow\ \frac{||y||_2}{||x||_2}=\frac{||R^{1/2}x||_2}{||x||_2}\le\sqrt{\sigma_1} }$

一方

$\displaystyle{(10)\quad \begin{array}{l} ||y'||_2^2=||\Sigma^{1/2} x'||_2^2=\sigma_1x_1'\,^2+\cdots+\sigma_nx_n'\,^2\\ \ge \sigma_n(x_1'\,^2+\cdots+x_n'\,^2)=\sigma_nx'\,^Tx'=\sigma_n||x'||_2^2 \ \Rightarrow\ \frac{||y'||_2}{||x'||_2}\ge\sqrt{\sigma_n} \end{array} }$

より、次式が成り立ちます。

$\displaystyle{(11)\quad \frac{||x'||_2}{||x||_2}\times\frac{||y'||_2}{||x'||_2}\times\frac{||y||_2}{||y'||_2}\ge\sqrt{\sigma_n} \ \Rightarrow\ \frac{||y||_2}{||x||_2}=\frac{||R^{1/2}x||_2}{||x||_2}\ge\sqrt{\sigma_n} }$

したがって、(7)と合わせて、次式を得ます。

$\displaystyle{(12)\quad \sqrt{\sigma_n} \le \frac{||R^{1/2}x||_2}{||x||_2}=\frac{\sqrt{x^TR^{1/2}R^{1/2}x}}{\sqrt{x^Tx}} \le\sqrt{\sigma_1} }$

これを２乗して次が得られます。

$\displaystyle{(13)\quad \sigma_n \le \frac{x^TRx}{x^Tx} \le\sigma_1 }$

すなわち(11)が導出されました。