C33 応用例 | 制御系CAD

応用例…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} }$

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

$\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 { \begin{array}{l} u(t)=-Fx(t)+Gr(t)+u_\ell(t)+u_n(t)\\ u_\ell(t)=-\underbrace{(SB)^{-1}(SA_m-\Phi S)}_{L=L_{eq}+L_{\Phi}}e(t)\\ u_n(t)=-\underbrace{(SB)^{-1}\rho(t,e)}_{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制御系応答

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=swfopt(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

SCILAB

[3] 積分動作導入による追従制御

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

$\displaystyle{(20)\quad { \begin{array}{l} u(t)=u_\ell(t)+u_n(t)\\ u_\ell(t)=-\underbrace{(SB)^{-1}(SA-\Phi S)x(t)}_{L=L_{eq}+L_\Phi}\\ -\underbrace{(SB)^{-1}(\Phi S_r+S_1B_r)}_{L_r} r(t) -\underbrace{(SB)^{-1}S_r}_{L_{\dot r}} \dot{r}(t)\\ u_n(t)=-\underbrace{(SB)^{-1}\rho(t,x)}_{L_n}\frac{P_2(s(t)-S_rr(t))}{||P_2(s(t)-S_rr(t))||} \end{array}} }$

●ここでは、等価制御による閉ループ系の $A_{eq}$ が望ましい固有値・固有ベクトルをもつように設計して、次式からスイッチング関数を決める行列 $S$ を決定します。

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

$A_{eq}$ の望ましい固有値を

$\displaystyle{(22)\quad \lambda(A_{eq}):\left\{\begin{array}{cl} -5.6\pm j 4.2 & short\ period\ mode\\ -1 & pitch\ attitude\ mode\\ -0.4 & elevator\ deflection\ mode\\ -0.7 & flap\ deflextion\ mode \end{array}\right. }$

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

$\displaystyle{(23)\quad {\tt specpos}= \left[\begin{array}{ccccc} 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0\\\hline 1 & 1 & 1 & 1 & 1\\ 1 & 0 & 1 & 1 & 1\\ 0 & 1 & 1 & 0 & 1\\\hline 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0\\ \end{array}\right],\quad {\tt specent}= \left[\begin{array}{rrrrr} -1 & -1 & -1 & -1 & -1\\ -1 & -1 & -1 & -1 & -1\\\hline 0 & 0 & 1 & 0 & 1\\ 1 & -1 & 0 & 1 & 0\\ -1 & 1 & 0 & -1 & 0\\\hline -1 & -1 & -1 & -1 & -1\\ -1 & -1 & -1 & -1 & -1\\ \end{array}\right] }$

このとき次のような行列 $S$ を得ます。

$\displaystyle{(24)\quad S=\left[\begin{array}{rrrrr|rr} 0.0531 & 0.0137 & -0.1436 & -0.0260 & -0.1373 & 0.0500 & 0\\ -0.0072 & -0.0612 & 0.1635 & 0.0035 & 0.1661 & 0 & 0.0500 \end{array}\right] }$

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

$\displaystyle{(25)\quad L=\left[\begin{array}{rrrrrrr} -1.0616 & -0.2743 & 2.9379 & 0.6060 & 2.4605 & -0.4927 & -0.0900\\ 0.1440 & 1.2236 & -3.3390 & -0.2296 & -3.2769 & 0.0671 & 0.0142 \end{array}\right] }$

$\displaystyle{(26)\quad L_{r}=\left[\begin{array}{rr} -2.9499 & 0.0120\\ 0.4000 & 2.9391 \end{array}\right] }$

$\displaystyle{(27)\quad L_{\dot{r}}=\left[\begin{array}{rr} -0.1448 & 0.0013\\ 0.0196 & 0.1439 \end{array}\right] }$

図3 積分動作導によるSMC制御系応答

MATLAB

%cAITCRAFT_smci.m
%-----
 clear all, close all
 A0= [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];
 B0=[0 0;0 0;0 0;20 0;0 20]; 
 H=[1 0 1 0 0;1 0 0 0 0];
 A=[zeros(2,2) -H ;
    zeros(5,2) A0];
 B=[zeros(2,2); B0];
 Gamma=diag([-0.9,-0.7]);
%----- 
 lambda=[-5.6 4.2 -1 -0.4 -0.7];
 nocomp=1;
 specpos=[zeros(2,5);[1 1 1 1 1;1 0 1 1 1;0 1 1 0 1];zeros(2,5)]
 specent=[-ones(2,5);[0 0 1 0 1;1 -1 0 1 0;-1 1 0 -1 0];-ones(2,5)]
 S=swfvpl(A,B,lambda,nocomp,specpos,specent);
 SS=S(:,6:7)\S*0.05
 F=(SS*B)\(SS*A)
 eig(A-B*F)
%-----
 [nn,mm]=size(B);
 S1=SS(:,1:nn-mm);
 S2=SS(:,nn-mm+1:nn);
 Br=eye(nn-mm,mm);
 B2=B(nn-mm+1:nn,:); 
 Lambda=S2*B2;
 P2=diag([0.025,0.025]);
 Phi=-20/0.025*P2;
 Check=P2*Phi+Phi*P2
 L=inv(Lambda)*(SS*A-Phi*SS)
%-----
%Sr=[-0.1448 0.0013;0.0196 0.1439]
 Sr=zeros(2)
 Lr=inv(Lambda)*(Phi*Sr+S1*Br)
 Ldr=inv(Lambda)*Sr
 Ln=0.1*inv(Lambda)
%-----
 sim('AIRCRAFT_smci_2015a.mdl')
%-----
%eof

SCILAB

制御技術セミナー(応用編)目次