応用例 | 制御系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} }$

また、評価変数として、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