C24 スイッチング関数(2)

スイッチング関数(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) }

ここでABSは、(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 { S[v_1\cdots v_{n-m}]=0 % \Rightarrow  [v_1\cdots v_{n-m}]^TS^T=0} }

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

1入力系の場合は固有ベクトル設定の自由度はありませんが、多入力系の場合は固有ベクトルの向きを設定する必要があり、これは元の状態方程式(状態変数の並び)に対して与えます。そこで(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を決定してよいことがわかります。

演習…Flipped Classroom

図1 台車駆動型倒立振子

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

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

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

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
SCILAB

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

MATLAB
%swfvpl.m
%-----
function S=swfvpl(A,B,lambda,nocomp,specpos,specent)
%lambda:n-m個の固有値
%nocomp:複素固有値の数
%specpos:固有ベクトルの非零要素の位置 (1:指定、0:指定せず)
%specent:固有ベクトルの要素の値 (0:零、1:非零指定、-1:非零任意)
 [nn,mm]=size(B);
 for i=1:nocomp
  glambda1=[lambda(2*i-1)*eye(nn)-A lambda(2*i)*eye(nn) B];
  [u,v,w]=svd(glambda1);
  nlambda1=w(1:nn,nn+1:2*nn+mm);
  plambda1=w(nn+1:2*nn,nn+1:2*nn+mm);
  alpha=[nlambda1; -plambda1]';
  beta= [plambda1;  nlambda1]';
  [u1,v1,w1]=svd(alpha);
  [u2,v2,w2]=svd(beta);
  gamma=[w1(:,nn+mm+1:2*nn)'; w2(:,nn+mm+1:2*nn)'];
  [u3,v3,w3]=svd(gamma);
  kgamma=w3(:,2*nn-2*mm+1:2*nn);
  despos=find(specpos(:,2*i-1));
  numspecr=length(despos);
  n1=[]; desent=[];
  for j=1:numspecr,
   n1(j,:)=kgamma(despos(j),:);
   desent(j)=specent(despos(j),2*i-1);
  end
  despos=find(specpos(:,2*i));
  numspecc=length(despos);
  for j=1:numspecc,
   n1(j+numspecr,:)=kgamma(despos(j)+nn,:);
   desent(j+numspecr)=specent(despos(j),2*i);
  end
  delta=n1\desent';
  vector=kgamma*delta;
  V(1:nn,2*i-1)=vector(1:nn);
  V(1:nn,2*i)=vector(nn+1:2*nn);
 end
%-----
 for i=2*nocomp+1:nn-mm,
  glambda=[lambda(i)*eye(nn)-A B];
  [u,v,w]=svd(glambda);
  nlambda=w(1:nn,nn+1:nn+mm);
  despos=find(specpos(:,i));
  numspec=length(despos);
  n2=[]; desent2=[];
  for j=1:numspec,
   n2(j,:)=nlambda(despos(j),:);
   desent2(j)=specent(despos(j),i);
  end
  delta=n2\desent2';
  V(:,i)=nlambda*delta;
 end
%-----
 [u,v,w]=svd(V);
 S=u(:,nn-mm+1:nn)';
end
%-----
%eof
SCILAB
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;
%固有ベクトルの要素の指定の有無 (1:指定、0:任意)
 specpos=rand(3,3)
%固有ベクトルの要素の値 (0:零指定、1:非零指定、-1:非零任意)       
 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;
%固有ベクトルの要素の指定の有無 (1:指定、0:任意)
 specpos=rand(3,3)
%固有ベクトルの要素の値 (0:零指定、1:非零指定、-1:非零任意)       
 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
SCILAB

C23 スイッチング関数(1)

スイッチング関数(1)

[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 \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] そこで、次の2次形式評価関数を最小にするようにMを決定する方法が提案されています。

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

を2次系式評価関数

\displaystyle{(17)\quad \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 \begin{array}{l} \hat{u}(t)=-\hat{F}x_1(t)\\ \hat{F}=M-Q_{22}^{-1}Q_{21} \end{array} }

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

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

を解いて

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

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

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

演習…Flipped Classroom

図1 台車駆動型倒立振子

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

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

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

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
SCILAB

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

MATLAB
%swflqr.m
%-----
function S=swfopt(A,B,Q)
 [nn,mm]=size(B);
 [Tr,temp]=qr(B);
 Tr=Tr';
 Tr=[Tr(mm+1:nn,:);Tr(1:mm,:)];
 Areg=Tr*A*Tr';
 Breg=Tr*B;
 A11=Areg(1:nn-mm,1:nn-mm);
 A12=Areg(1:nn-mm,nn-mm+1:nn);
 A21=Areg(nn-mm+1:nn,1:nn-mm);
 A22=Areg(nn-mm+1:nn,nn-mm+1:nn);
 B2=Breg(nn-mm+1:nn,1:mm);
%-----
 Qt=Tr*Q*Tr';
 Q11=Qt(1:nn-mm,1:nn-mm);
 Q12=Qt(1:nn-mm,nn-mm+1:nn);
 Q21=Qt(nn-mm+1:nn,1:nn-mm);
 Q22=Qt(nn-mm+1:nn,nn-mm+1:nn);
 Qhat=Q11-Q12*inv(Q22)*Q21;
 Ahat=A11-A12*inv(Q22)*Q21;
 [K,P1,E]=lqr(Ahat,A12,Qhat,Q22);
 M=inv(Q22)*(A12'*P1+Q21);
 S2=eye(mm);
 S=S2*[M eye(mm)]*Tr;
end
%-----
%eof
SCILAB
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
SCILAB

C22 SM制御則

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) \end{array} }

以下では、状態方程式とスイッチング関数は次の標準形(regular form)をとるように座標変換されていると仮定します({\rm rank}B={\rm rank}S=m)。

\displaystyle{(2)\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) + \underbrace{ \left[\begin{array}{l} f_u(t,x)\\ f_m(t,x,u) \end{array}\right] }_{f(t,x,u)} }

\displaystyle{(3)\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{(4)\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{(5)\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} \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{ \begin{array}{l} (6.1)\quad \bar{A}_{11}=A_{11}-A_{12}M\\ (6.2)\quad \bar{A}_{12}=A_{12}S_2^{-1}\\ (6.3)\quad \bar{A}_{21}=S_2(M\bar{A}_{11}+A_{21}-A_{22}M)\\ (6.4)\quad \bar{A}_{22}=S_2(M{A}_{12}+A_{22})S_2^{-1}\\ (6.5)\quad \bar{B}_{2}=S_2B_2 \end{array} }

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

このとき、SMC則

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

を、2次安定性

\displaystyle{{ \begin{array}{lll} (8.1)\quad V(s)=s^T(t)P_2s(t)&\Rightarrow \dot{V}(s)\le -s^T(t)s(t)&(t\le t_s)\\ (8.2)\quad V(x_1)=x_1^T(t)P_1x_1(t)&\Rightarrow \dot{V}(x_1)\le -x_1^T(t)Q_1x_1(t)&(t> t_s) \end{array}} }

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

[2] 可到達性の検討 (t\le t_s)

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

\displaystyle{ \begin{array}{cl} (9.1) & s(t)=0\\ \Downarrow &\\ (9.2) & \dot{s}(t)=0\\ \Downarrow &\\ (9.3) & \bar{A}_{21}x_1(t)+\bar{A}_{22}s(t)+\bar{B}_{2}u(t)=0\\ \Downarrow &\\ (9.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 &S=[0\ I]T_s, x=T_s^{-1}\bar{x}\\ (9.5) & u_{eq}(t)=-(SB)^{-1}SAx(t)} \end{array} }

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

\displaystyle{(10.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{(10.2)\quad { u_\ell(t)=-\underbrace{(SB)^{-1}(SA-\Phi S)}_{L=L_{eq}+L_{phi}}x(t)} }

のように構成します。ここで、\Phiは適当に設定された安定行列です。このとき閉ループ系は次式で与えられます。

\displaystyle{(11)\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} (12.1)\quad \dot{x}_1(t)=\bar{A}_{11}x_1(t)+\bar{A}_{12}s(t)+f_u(t,x)\\ (12.2)\quad \dot{s}(t)=\Phi s(t)+\bar{B}_{2}u_n(t)+S_2(Mf_u(t,x)+f_m(t,x,u)) \end{array} }

ここで、\Phiは安定行列なので

\displaystyle{(13)\quad \begin{array}{l} P_2\Phi+\Phi^TP_2=-I \end{array} }

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

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

すなわち

\displaystyle{(14.2)\quad { u_\ell(t)=-\underbrace{(SB)^{-1}\rho(t,x)}_{L_n}\frac{P_2s(t)}{||P_2s(t)||} }

と定めます。\rho(t,x)は、仮定

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

のもとで

\displaystyle{{ (16)\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}||} %\end{array} }}

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

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

に注意して

\displaystyle{ \begin{array}{cl} (18.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\ (17.1)\\ (18.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 &  \\ (18.3) & \ge ||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\ (17.2)  \\ (18.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\ (17.3) \\ (18.5) & \ge ||S_2|| (||M||(k_1||x(t)||+k_2)+k_3||u(t)||+\alpha(t,x))+\gamma_2}\\ \Downarrow & by\ (15.1),(15.2) \\ (18.6) & \ge ||S_2|| (||M||\,||f_u(t,x)||+||f_m(t,x,u)||)+\gamma_2} \end{array} }

を得て、次式が成り立ちます。

\displaystyle{ \begin{array}{cl} (19.1) & \dot{V}(s)=2s^T(t)P_2\dot{s}(t)\\ \Downarrow & by\ (12.2) \\ (19.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 & \\ (19.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\ (13) \\ (19.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\ (18.6) \\ (19.5) & \le -||s(t)||^2-2\gamma_2||P_2s(t)||\\ (19.6) & \le -||s(t)||^2 \end{array} }

これより(8)の第1式が成り立ちます。

●到達時刻t_sを見積もるために、(19)から

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

も成り立つことと

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

に注意して

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

よって

\displaystyle{(23)\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(P_2) V^{1/2}(s)\\ \Rightarrow \frac{d}{dt}V^{1/2}(s)< -\gamma_2\sigma_n(P_2) \end{array} }

この両辺を積分すると

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

を得て、スライディング直線に到達する時刻t_sを見積もることができます。

[2] スライディングモードの検討 (t> t_s)

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

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

を満たすP_1>0を選ぶことができます。

\displaystyle{ \begin{array}{cl} (26.1)&\dot{V}(x_1)=2x_1^T(t)P_1\dot{x}_1(t)\\ \Downarrow & by\ (5) \\ (26.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 & \\ (26.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\ (25) \\ (26.4)&=-x_1^T(t)Q_1x_1(t)+2x_1^T(t)P_1f_u(t,x)\\ \Downarrow &\\ (26.5)&\le -\sigma_n^2(Q_1)x_1^T(t)x_1(t)+2||x_1^T(t)P_1||\,||f_u(t,x)||\\ \Downarrow &\\ (26.6)&\le -\sigma_n^2(Q_1)||x_1(t)||^2+2\sqrt{x_1^T(t)P_1^2x_1(t)}\,||f_u(t,x)||\\ \Downarrow &\\ (26.7)&\le -\sigma_n^2(Q_1)||x_1(t)||^2+2\sqrt{\sigma_1^2(P_1^2)x_1^T(t)x_1(t)}\,||f_u(t,x)||\\ \Downarrow &\\ (26.8)&= -\sigma_n^2(Q_1)||x_1(t)||^2+2\sigma_1^2(P_1)||x_1(t)||\,||f_u(t,x)||\\ \Downarrow &\\ (26.9)&= -\sigma_1^2(P_1)||x_1(t)||(\frac{\sigma_n^2(Q_1)}{\sigma_1^2(P_1)}||x_1(t)||-2||f_u(t,x)||) \end{array} }

ここで

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

を仮定します。このとき

\displaystyle{(28)\quad \begin{array}{l} \dot{V}(x_1)\le -\sigma_1^2(P_1)||x_1(t)||(\frac{\sigma_n^2(Q_1)}{\sigma_1^2(P_1)}||x_1(t)||-2||f_u(t,x)||)\\ \le -\sigma_1^2(P_1)||x_1(t)||(\frac{\sigma_n^2(Q_1)}{\sigma_1^2(P_1)}||x_1(t)||)\\ =-\sigma_n^2(Q_1)||x_1(t)||^2\\ \le -x_1^T(t)Q_1x_1(t) \end{array} }

これより(8)の第2式が成り立ちます。

Note 1 標準形(2)について

列フルランクをもつ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} }

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

\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 2 2次安定性の条件(8)について
(8)は次のようにまとめられることに注意してください。

\displaystyle{(1)\quad  \begin{array}{lll} V(x')= \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right]^T }_{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] }_{x'(t)}\\ = x_1^T(t)P_1x_1(t)+s^T(t)P_2s(t)\\ \Rightarrow \dot{V}(x')\le - \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right]^T }_{x'^T(t)} \underbrace{ \left[\begin{array}{cc} Q_1 & 0\\ 0 & Q_2 \end{array}\right] }_{P} \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right] }_{x'(t)}\\ =-x_1^T(t)Q_1x_1(t)-s^T(t)Q_2s(t)\quad(Q_2=I) \end{array}} }

ここで、t>t_sではs=0ですから、(8.2)を意味します。一方、t\le t_sでは\bar{A}_{11}は安定行列なので((25)から-x_1^T(t)Q_1x_1(t)<0

\displaystyle{(2)\quad  \dot{V}(x')\le -x_1^T(t)Q_1x_1(t)-s^T(t)s(t)\le -s^T(t)s(t)\\ }

となって、(8.1)を意味します。

Note 3 条件(15)について
次式で表される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}{ccc} 0 & 1 & 0\\ 0 & 0 & \frac{K_t}{J_0}\\ 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\\ \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}{cc} \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.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) }

と表されるので、仮定(15)の第2式における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 f_u=-\xi_J\dot{\omega}=-\xi_J(\frac{J_0}{K_t}\omega_n^2\theta+2\frac{J_0}{K_t}\zeta\omega_n\dot{\theta}) }

と表されるので、仮定(15)の第1式におけるk_1k_2は次のように特定できます。

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

C21 問題設定

問題設定…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) }

いま(1)の解が、ある時刻t_sに対して

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

を満足するとき、スライディングモードが発生していると言います。

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

問題1: 動的システムの振る舞いが閉じ込められる超平面 {\cal S}=\{x:Sx=0\} を定義するスイッチング関数を決定せよ。
問題2: 有限時刻 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 { u_\ell(t)=-\underbrace{(SB)^{-1}(SA-\Phi S)}_{L=L_{eq}+L_{phi}}x(t)} }

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

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

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

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

[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 \dot{x}(t)=\underbrace{(I_n-B(SB)^{-1}S)A}_{A_{eq}}x(t)\quad(t\ge t_s) }

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

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

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

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

\displaystyle{(12)\quad { 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とおけば、21.2節の(14)となります。

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

Note 1入力系の場合…Homework

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

\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 u_\ell(t)=a_1x_1(t)+\sum_{i=2}^{n}(a_i-m_{i-1}) x_i(t) }

のように定めます。

u_nの候補として次の2つを考えます。

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

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

C12 固有ベクトル設定問題

[0] 次の2入力2次系

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

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

\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を持ちます。このように多入力系では閉ループ系の固有値を指定しても状態フィードバックは一意に定まりません。これは固有ベクトルの指定に任意性があるからです(ちなみに1入力系では一意に定まります)。そこで、上の2つの閉ループ系の固有ベクトルを求めてみると

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

の最小2乗解として\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 v=(\lambda I_n-A)^{-1}B\delta, \xi=-\delta\quad\Rightarrow\quad 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] }

の最小2乗解として\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 \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} }

演習…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;
%固有ベクトルの要素の指定の有無 (1:指定、0:任意)
 specpos=[ 1  1; 
           1  1];       
%固有ベクトルの要素の値 (0:零指定、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;
%固有ベクトルの要素の指定の有無 (1:指定、0:任意)
 specpos=[ 1  1; 
           1  1];        
%固有ベクトルの要素の値 (0:零指定、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(2,2)])
SCILAB

2^\circ 1入力系の場合は、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;
%固有ベクトルの要素の指定の有無 (1:指定、0:任意)
 specpos=rand(2,2)
%固有ベクトルの要素の値 (0:零指定、1:非零指定、-1:非零任意)       
 specent=rand(2,2) 
 F=vplace2(A,B,lambda,nocomp,specpos,specent) 
 [V,pl]=eig(A-B*F)
%-----
 lambda=[-1 1];
 nocomp=1;
%固有ベクトルの要素の指定の有無 (1:指定、0:任意)
 specpos=rand(2,2)
%固有ベクトルの要素の値 (0:零指定、1:非零指定、-1:非零任意)       
 specent=rand(2,2) 
 F=vplace2(A,B,lambda,nocomp,specpos,specent) 
 [V,pl]=eig(A-B*F) 
%-----
%eof
SCILAB

Note 固有ベクトル設定プログラム
Edwards and Spurgeonによるプログラムspvplを基に作成しています。

MATLAB
%vplace2.m
%-----
function F=vplace2(A,B,lambda,nocomp,specpos,specent)
%lambda:n-m個の固有値
%nocomp:複素固有値の数
%specpos:固有ベクトルの非零要素の位置 (1:指定、0:指定せず)
%specent:固有ベクトルの要素の値 (0:零、1:非零指定、-1:非零任意)
 [nn,mm]=size(B);
 for i=1:nocomp,
   glambda1=[lambda(2*i-1)*eye(nn)-A lambda(2*i)*eye(nn) B];
   [u,v,w]=svd(glambda1);
   nlambda1=w(1:nn,nn+1:2*nn+mm); 
   plambda1=w(nn+1:2*nn,nn+1:2*nn+mm);  
   mlambda1=w(2*nn+1:2*nn+mm,nn+1:2*nn+mm);  %%%     
   alpha=[nlambda1; -plambda1]'; 
   beta= [plambda1;  nlambda1]';
   [u1,v1,w1]=svd(alpha);
   [u2,v2,w2]=svd(beta); 
   gamma=[w1(:,nn+mm+1:2*nn)'; w2(:,nn+mm+1:2*nn)'];
   [u3,v3,w3]=svd(gamma);
   kgamma=w3(:,2*nn-2*mm+1:2*nn);
   q1=[nlambda1;-plambda1]\kgamma;  %%%  
   q2=[plambda1; nlambda1]\kgamma;  %%%    
   despos=find(specpos(:,2*i-1));
   numspecr=length(despos);
   n1=[]; desent=[];
   for j=1:numspecr,
     n1(j,:)=kgamma(despos(j),:);
     desent(j)=specent(despos(j),2*i-1);
   end
   despos=find(specpos(:,2*i));
   numspecc=length(despos);
   for j=1:numspecc,
     n1(j+numspecr,:)=kgamma(despos(j)+nn,:);
     desent(j+numspecr)=specent(despos(j),2*i);
   end
   delta=n1\desent'; 
   vector=kgamma*delta;
   V(1:nn,2*i-1)=vector(1:nn);
   V(1:nn,2*i)=vector(nn+1:2*nn);   
   W(1:mm,2*i-1)=mlambda1*q1*delta;  %%%  
   W(1:mm,2*i)=mlambda1*q2*delta;    %%%     
 end
 for i=2*nocomp+1:nn,
   glambda=[lambda(i)*eye(nn)-A B];
   [u,v,w]=svd(glambda);
   nlambda=w(1:nn,nn+1:nn+mm);
   mlambda=w(nn+1:nn+mm,nn+1:nn+mm); %%%  
   despos=find(specpos(:,i));
   numspec=length(despos);
   n2=[]; desent2=[];
   for j=1:numspec,
     n2(j,:)=nlambda(despos(j),:);
     desent2(j)=specent(despos(j),i);
   end
   delta=n2\desent2'  
   V(:,i)=nlambda*delta;
   W(:,i)=mlambda*delta;  %%% 
 end
 F=W/V;
end 
%-----
%eof
SCILAB

C11 2次安定性

リャプノフの安定定理

非線形システム理論に関する文献:
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)の漸近安定性を示すことができます。

リャプノフ方程式

●リャプノフの安定定理を線形系

\displaystyle{(5)\quad \dot{x}(t)=Ax(t)\quad(x(t)\in{\bf R}^n)}

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

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

これを微分して(5)を用いると((5)の解に沿って微分して)

\displaystyle{(7)\quad \begin{array}{l} \frac{d}{dt}V(x)=\dot{x}^T(t)Px(t)+x^T(t)P\dot{x}(t)\\ =x^T(t)A^TPx(t)+x^T(t)PAx(t)\\ =x^T(t)(A^TP+PA)x(t) \end{array} }

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

\displaystyle{(8)\quad {A^TP+PA<0\quad(P>0)}}

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

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

\displaystyle{(9)\quad x(t)=\exp(At)x(0)\rightarrow 0\quad(t\rightarrow \infty)}

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

\displaystyle{(10)\quad P=\int_0^\infty \exp(A^Tt)Q\exp(At)dt}

ここで、Q>0は任意です(ある条件下ではQ\ge0としてもよい)。実際

\displaystyle{(11)\quad \begin{array}{l} \displaystyle{A^T\int_0^\infty \exp(A^Tt)Q\exp(At)dt+\int_0^\infty \exp(A^Tt)Q\exp(At)dt\,A}\\ \displaystyle{=\int_0^\infty \frac{d}{dt}(\exp(A^Tt)Q\exp(At))\,dt}\\ \displaystyle{=\left[\exp(A^Tt)Q\exp(At)\right]_0^\infty=-Q<0} \end{array} }

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

(12)\quad \begin{array}{l} \displaystyle{v^TPv=\int_0^\infty v^T\exp(A^Tt)Q\exp(At)vdt}\\ \displaystyle{=\int_0^\infty ||Q^{1/2}\exp(At)v||dt\ge0} \end{array}

ですが、v^TPv=0とすると、Q^{1/2}\exp(At)v=0が得られ、これはv=0を意味し矛盾です。したがってv^TPv>0でなけれならず、P>0を得ます。以上から、線形行列不等式(8)が(5)の漸近安定性の必要十分条件になっています。

●一方、リャプノフ方程式と呼ばれる線形行列方程式

\displaystyle{(13)\quad PA+A^TP+Q=0\quad(Q>0)}

の解P>0を求めてリャプノフ関数を決めることができます。すなわち、線形行列方程式(13)も(5)の漸近安定性の必要十分条件となっていることに留意してください。これはスカラ系\dot{x}(t)=ax(t)の場合、次のように自明です。

\displaystyle{(14)\quad \exists p>0:2ap=-q<0\ (q>0)\quad\Leftrightarrow\quad a<0}

2次安定性

非線形系(1)に対してリャプノフ関数(6)を考えます。これを微分して(1)を用いると((2)の解に沿って微分して)

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

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

\displaystyle{(16)\quad  { %\begin{array}{l} \frac{d}{dt}V(x)=2x^T(t)Pf(x)\le -x^T(t)Qx(t)\quad(Q>0) %\end{array} }}

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

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

\displaystyle{(18)\quad {\alpha=\sigma_n(P^{-1}Q),\ \beta=\sqrt{\frac{\sigma_1(P)}{\sigma_n(P)} }}

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

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

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

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

\displaystyle{(20)\quad {\sigma_n(R)\le\frac{v^TRv}{v^Tv}\le\sigma_1(R) %,\quad\sigma_1^{-2}(R)\le\frac{v^TR^{-1}v}{v^Tv}\le\sigma_n^{-2}(R) }}

および

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

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

●そこで

\displaystyle{(22)\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{(23)\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{(24)\quad {\frac{\sigma_n(Q)}{\sigma_1(P)}=\frac{1}{\sigma_1(P)\sigma_1(Q^{-1})} \le \frac{1}{\sigma_1(PQ^{-1})}=\sigma_n(P^{-1}Q) }}

を用いて

\displaystyle{(25)\quad {\frac{d}{dt}V(x)\le -\underbrace{\sigma_n(P^{-1}Q)}_{\alpha} V(x) }}

が成り立ちます。これから

\displaystyle{(26)\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{(27)\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{(28)\quad {\sigma_n(P)x^T(t)x(t) \le \sigma_1(P) x^T(0)x(0)e^{-\alpha t} }}

すなわち

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

を得ます。これから(17)が成り立ちます。

Note (20)の導出

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

\displaystyle{(30)\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_1Rの最大固有値、特異値\sigma_nRの最小固有値であることに注意してください。このとき、R=R^{1/2}R^{1/2}を満たすRの平方根行列は

\displaystyle{(31)\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{(32)\quad \begin{array}{l} R^{-1}=U\Sigma^{-1} U^T\\ \Sigma^{-1}={\rm diag}\{\frac{1}{\sigma_1}\ge\cdots\ge\frac{1}{\sigma_n}\} \quad(\frac{1}{\sigma_n}\ge\cdots\ge\frac{1}{\sigma_1}) \end{array} }

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

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

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

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

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

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

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

\displaystyle{(36)\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{(37)\quad ||y||_2^2=||Uy'||_2^2=y'\,^TU^TUy'=y'\,^Ty'=||y'||_2^2\ \Rightarrow\ \frac{||y||_2}{||y'||_2}=1 }

すなわち

\displaystyle{(38)\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{(39)\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{(40)\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} }

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

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

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

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

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