YMCPBの特性解析 | 制御系CAD

YMCPBの線形モデル

●制御対象YMCPBの運動方程式として次式を考えます。

$\displaystyle{ \begin{array}{ll} (1.1) & (M+M_x)\ddot{x}=D_b\cos\theta+(T+D_e)\cos(\theta+\theta_e)\\ &+N_L\sin\theta+N_e\sin(\theta+\theta_e)\\ (1.2)&(M+M_z)\ddot{z}=-D_b\sin\theta-(T+D_e)\sin(\theta+\theta_e)\\ &+N_L\cos\theta+N_e\cos(\theta+\theta_e)+N_B+Mg\\ &-c_{zz}\dot{z}-c_{z\theta}\dot{\theta}\\ (1.3)&(I_y+J_y)\ddot{\theta}=D_b(H_{CG}-H_D)+T(H_T-H_{De}+H_e)\\ &+D_eH_e+N_LL_L+N_eL_e+N_B(L_{CG}\cos\theta-L_B)\\ &-c_{\theta z}\dot{z}-c_{\theta\theta}\dot{\theta} \end{array} }$

すなわち

$\displaystyle{(1')\quad \begin{array}{l} \left[\begin{array}{ccc} M+M_x & 0 & 0 \\ 0 & M+M_z & 0\\ 0 & 0 & I_y+J_y \end{array}\right] \left[\begin{array}{c} \ddot{x}\\ \ddot{z}\\ \ddot{\theta} \end{array}\right]= \left[\begin{array}{ccc} 0 & 0 & 0 \\ 0 & c_{zz} & c_{z\theta}\\ 0 & c_{\theta z} & c_{\theta\theta} \end{array}\right] \left[\begin{array}{c} \dot{x}\\ \dot{z}\\ \dot{\theta} \end{array}\right]\\ +\left[\begin{array}{c} D_b\cos\theta+(T+D_e)\cos(\theta+\theta_e)\\ -D_b\sin\theta-(T+D_e)\sin(\theta+\theta_e)\\ D_b(H_{CG}-H_D)+T(H_T-H_{De}+H_e)+D_eH_e \end{array}\right]\\ +\left[\begin{array}{c} N_L\sin\theta+N_e\sin(\theta+\theta_e)\\ N_L\cos\theta+N_e\cos(\theta+\theta_e)\\ N_LL_L+N_eL_e \end{array}\right] + \left[\begin{array}{c} 0 \\ N_B+Mg\\ N_B(L_{CG}\cos\theta-L_B) \end{array}\right] \end{array} }$

ここで、 $x,z,\theta$ はそれぞれサージ、ヒーブ、ピッチを、 $T,\theta_e$ はそれぞれスラスト、船外機取付角を表しています。その他の物理パラメータの説明はここでは省略します。

●いま状態変数ベクトルと操作変数ベクトルをそれぞれ

$\displaystyle{(2)\quad \xi=\left[\begin{array}{c} x\\ z\\ \theta\\\hline \dot{x}\\ \dot{z}\\ \dot{\theta} \end{array}\right],\ \zeta=\left[\begin{array}{c} T\\ \theta_e \end{array}\right] }$

ととると、次の非線形状態方程式を得ます。

$\displaystyle{(3.1)\quad \underbrace{\frac{d}{dt}\left[\begin{array}{c} x\\ z\\ \theta\\\hline \dot{x}\\ \dot{z}\\ \dot{\theta} \end{array}\right]}_{\dot{\xi}}= \underbrace{ \left[\begin{array}{c} f_1(x,z,\theta,\dot{x},\dot{z},\dot{\theta},T,\theta_e)\\ f_2(x,z,\theta,\dot{x},\dot{z},\dot{\theta},T,\theta_e)\\ f_3(x,z,\theta,\dot{x},\dot{z},\dot{\theta},T,\theta_e)\\\hline f_4(x,z,\theta,\dot{x},\dot{z},\dot{\theta},T,\theta_e)\\ f_5(x,z,\theta,\dot{x},\dot{z},\dot{\theta},T,\theta_e)\\ f_6(x,z,\theta,\dot{x},\dot{z},\dot{\theta},T,\theta_e) \end{array}\right]}_{f(\xi,\zeta)} }$

ただし

$\displaystyle{(3.2)\quad \left\{\begin{array}{l} f_1(x,z,\theta,\dot{x},\dot{z},\dot{\theta},T,\theta_e)=\dot{x}\\ f_2(x,z,\theta,\dot{x},\dot{z},\dot{\theta},T,\theta_e)=\dot{z}\\ f_3(x,z,\theta,\dot{x},\dot{z},\dot{\theta},T,\theta_e)=\dot{\theta}\\ f_4(x,z,\theta,\dot{x},\dot{z},\dot{\theta},T,\theta_e)=\\ \quad\frac{1}{M+M_x}(D_b\cos\theta+(T+D_e)\cos(\theta+\theta_e)\\ \quad+N_L\sin\theta+N_e\sin(\theta+\theta_e))\\ f_5(x,z,\theta,\dot{x},\dot{z},\dot{\theta},T,\theta_e)=\\ \quad\frac{1}{M+M_z}(-Db\sin\theta-(T+D_e)\sin(\theta+\theta_e)\\ \quad+N_L\cos\theta+N_e\cos(\theta+\theta_e)+N_B+Mg\\ \quad-c_{zz}\dot{z}-c_{z\theta}\dot{\theta})\\ f_6(x,z,\theta,\dot{x},\dot{z},\dot{\theta},T,\theta_e)=\\ \quad\frac{1}{I_y+J_y}(D_b(H_{CG}-H_D)+T(H_T-H_{De}+H_e)\\ \quad+D_eH_e+N_LL_L+N_eL_e+N_B(L_{CG}\cos\theta-L_B)\\ \quad-c_{\theta z}\dot{z}-c_{\theta\theta}\dot{\theta}) \end{array}\right. }$

●制御対象YMCPBに対する制御目的は、あるスラスト $T^*$ と取付角 $\theta_e^*$ の下で、巡航速度 $V^*$ で走行するとき、一定の姿勢 $z^*,\theta^*$ を保つこととします。これは、次の平衡状態 $\xi^*$ を、次の平衡入力 $\zeta^*$ によって保持することとみなすことができます。

$\displaystyle{(4)\quad \xi^*=\left[\begin{array}{c} V^*t\\ z^*\\ \theta^*\\\hline V^*\\ 0\\ 0 \end{array}\right],\ \zeta^*=\left[\begin{array}{c} T^*\\ \theta_e^* \end{array}\right] }$

この平衡状態と平衡入力を定めるためには、 $T^*$ と $V^*$ を所与として、残りの $z^*,\theta^*,\theta_e^*$ を次の非線形連立方程式を解いて求めます。

$\displaystyle{(5)\quad \left\{\begin{array}{l} f_4(V^*t,z,\theta,V^*,0,0,T^*,\theta_e)=0\\ f_5(V^*t,z,\theta,V^*,0,0,T^*,\theta_e)=0\\ f_6(V^*t,z,\theta,V^*,0,0,T^*,\theta_e)=0 \end{array}\right. }$

したがって、次式が成り立つことに注意します。

$\displaystyle{(6)\quad f(\xi^*,\zeta^*)= \left[\begin{array}{c} V^*\\ 0\\ 0\\\hline 0\\ 0\\ 0 \end{array}\right]=\frac{d}{dt}\xi^* }$

●この平衡状態 $\xi^*$ と平衡入力 $\zeta^*$ の回りで、非線形状態方程式 $\dot{\xi}=f(\xi,\zeta)$ を線形近似します。

$\displaystyle{(7)\quad \dot{\xi}\simeq f(\xi^*,\zeta^*) +\underbrace{\left.\frac{\partial f(\xi,\zeta)}{\partial \xi}\right|_{\xi=\xi^*,\zeta=\zeta^*}}_{\frac{\partial f(\xi^*,\zeta^*)}{\partial \xi}}(\xi-\xi^*) +\underbrace{\left.\frac{\partial f(\xi,\zeta)}{\partial \zeta}\right|_{\xi=\xi^*,\zeta=\zeta^*}}_{\frac{\partial f(\xi^*,\zeta^*)}{\partial \zeta}}(\zeta-\zeta^*) }$

ここで(6)に注意して、線形状態方程式

$\displaystyle{(8)\quad \underbrace{\frac{d}{dt}(\xi-\xi^*)}_{\dot{x}}= \underbrace{\frac{\partial f(\xi^*,\zeta^*)}{\partial \xi}}_{A}\underbrace{(\xi-\xi^*)}_{x} +\underbrace{\frac{\partial f(\xi^*,\zeta^*)}{\partial \zeta}}_{B}\underbrace{(\zeta-\zeta^*)}_{u} }$

すなわち

$\displaystyle{(9)\quad \begin{array}{c} \underbrace{ \frac{d}{dt} \left[\begin{array}{c} x(t)-V^*t\\ z(t)-z^*\\ \theta(t)-\theta^*\\\hline \dot{x}(t)-V^*\\ \dot{z}(t)\\ \dot{\theta}(t) \end{array}\right] }_{\dot{x}(t)} = \underbrace{ \left[\begin{array}{ccc|ccc} 0 & 0 & 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 0 & 0 & 1\\\hline 0 & a_{42} & a_{43} & a_{44} & a_{45} & a_{46} \\ 0 & a_{52} & a_{53} & a_{54} & a_{55} & a_{56} \\ 0 & a_{62} & a_{63} & a_{64} & a_{65} & a_{66} \end{array}\right] }_{A} \underbrace{ \left[\begin{array}{c} x(t)-V^*t\\ z(t)-z^*\\ \theta(t)-\theta^*\\\hline \dot{x}(t)-V^*\\ \dot{z}(t)\\ \dot{\theta}(t) \end{array}\right] }_{x(t)}\\ + \underbrace{ \left[\begin{array}{cc} 0 & 0 \\ 0 & 0 \\ 0 & 0 \\\hline b_{41} & b_{42} \\ b_{51} & b_{52} \\ b_{61} & b_{62} \end{array}\right] }_{B} \underbrace{ \left[\begin{array}{c} T(t)-T^*\\ \theta_e(t)-\theta_e^*\\ \end{array}\right] }_{u(t)} \end{array}}$

を得ます。

●次の５組のデータセットが与えられています。

	艇	艇質量[kg]	長手方向重心位置[m]	備考
No1	SPTM247	2709	1.9756	INM抵抗試験レポートC2620-07CT18-RAP01による
No2	SPTM247	2813	1.925	2021/8/2重量重心測定結果より
No3	SPTM247	2845	1.937	2022/1/19重量重心測定結果より
No4	SPTM247	2945	1.937	No3重量に+100kg
No5	SPTM247	2745	1.937	No3重量に-100kg

各データセットは、133個（ $k=1,\cdots,133$ ）の平衡状態と平衡入力と $A$ 行列を含みます（ $B$ 行列と $C$ 行列は固定）。これらは、たとえば、データセットNo1に対しては図1のようにグラフ化されます。

図1 データセットNo1における平衡状態

ここで、赤い〇の平衡状態は、制御系の設計および評価のために選ばれた、次表に示す39点です。 $k=128$ の場合の平衡状態回りの線形モデルについてコントローラを設計します。

$-\theta_e^*$	run	k
-4 deg	run1	8,9,10,11,12,13,14
0 deg	run1	36,37,38,39,40,41,42,43
4 deg	run1	66,67,68,69,70,71,72,73
8 deg	run1	96,97,98,99,100,101,102,103
2 deg	run1	126,127,128,129,130,131,132,133

またデータセットNo1に対する行列 $A$ の固有値は図2のようにグラフ化されます。

図2 データセットNo1における行列 $A$ の固有値

各線形モデルに対して複素固有値が２対、実固有値が２個（うち１個は零）あることがわかります。不安定な複素固有値はピッチの振舞いに、安定な複素固有値はヒーブの振舞いに、非零実固有値は速度変動に関係しています。速度変動が発散する場合があることが気になるところです。

制御系設計用線形モデル

●制御系設計用線形モデルを得るために、まずアクチュエータについては、船外機の取付角を操作するためにレバーを引いた期間だけ一定の角速度 $\omega_e$ で回転するものとします。これを次式で表します。

$\displaystyle{(10)\quad \dot{\theta}_e(t)=\omega_e u_e(t)\quad(u_e(t)=0,\pm 1) }$

ここで、 $u_e$ は操作入力 $u(t)$ が正のとき $+1$ 、負のとき $-1$ の値をとるものとします。すなわち

$\displaystyle{(11)\quad u_e(t)={\rm sign}(u(t)) }$

したがって、制御対象の状態方程式(9)とアクチュエータ(10)を合わせて、次式を得ます。

$\displaystyle{(12)\quad \begin{array}{c} \underbrace{ \frac{d}{dt} \left[\begin{array}{c} x(t)-V^*t\\ z(t)-z^*\\ \theta(t)-\theta^*\\\hline \dot{x}(t)-V^*\\ \dot{z}(t)\\ \dot{\theta}(t)\\\hline \theta_e(t)-\theta_e^* \end{array}\right] }_{\dot{x}_a(t)} = \underbrace{ \left[\begin{array}{ccc|ccc|c} 0 & 0 & 0 & 1 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 1 & 0\\\hline 0 & a_{42} & a_{43} & a_{44} & a_{45} & a_{46} & b_{42}\\ 0 & a_{52} & a_{53} & a_{54} & a_{55} & a_{56} & b_{52}\\ 0 & a_{62} & a_{63} & a_{64} & a_{65} & a_{66} & b_{62}\\\hline 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{array}\right] }_{A_a} \underbrace{ \left[\begin{array}{c} x(t)-V^*t\\ z(t)-z^*\\ \theta(t)-\theta^*\\\hline \dot{x}(t)-V^*\\ \dot{z}(t)\\ \dot{\theta}(t)\\\hline \theta_e(t)-\theta_e^* \end{array}\right] }_{x_a(t)}\\ + \underbrace{ \left[\begin{array}{c} 0 \\ 0 \\ 0 \\\hline 0 \\ 0 \\ 0 \\\hline \omega_e \end{array}\right] }_{B_a}u_e(t) + \left[\begin{array}{c} 0 \\ 0 \\ 0 \\\hline b_{41} \\ b_{51} \\ b_{61} \\\hline 0 \end{array}\right]\underbrace{(T(t)-T^*)}_{0} \end{array}}$

また状態変数のうち、 $\theta$ 、 $\dot{\theta}$ 、 $\theta_e$ が計測できるものとします。このとき、状態方程式(12)に対する観測方程式は次式で表されます。

$\displaystyle{(13)\quad \begin{array}{l} \underbrace{ \left[\begin{array}{c} \theta(t)-\theta^*\\ \dot{\theta}(t)\\ \theta_e(t)-\theta_e^* \end{array}\right] }_{y(t)} = \underbrace{ \left[\begin{array}{ccc|ccc|c} 0 & 0 & 1 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 1 \end{array}\right] }_{C} \underbrace{ \left[\begin{array}{c} x(t)-V^*t\\ z(t)-z^*\\ \theta(t)-\theta^*\\\hline \dot{x}(t)-V^*\\ \dot{z}(t)\\ \dot{\theta}(t)\\\hline \theta_e(t)-\theta_e^* \end{array}\right] }_{x(t)} \end{array}}$

●制御対象YMCPBに対する制御目的は、あるスラスト $T^*$ と取付角 $\theta_e^*$ の下で、巡航速度 $V^*$ で走行するとき、一定の姿勢 $z^*,\theta^*$ を保つこととします。この制御目的に照らして、サージ方向の状態変数は省いてよく、また操作変数も船外機取付角に限定することができます。ヒーブ方向の状態変数を残すか省くかによって、次の２種類の線形モデルが考えられます。

５次元モデル

$\displaystyle{(14.1)\quad \begin{array}{l} \underbrace{ \frac{d}{dt} \left[\begin{array}{c} z(t)-z^*\\ \theta(t)-\theta^*\\ \dot{z}(t)\\ \dot{\theta}(t)\\ \theta_e(t)-\theta_e^* \end{array}\right] }_{\dot{x}(t)} = \underbrace{ \left[\begin{array}{ccccc} 0 & 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 1 & 0\\ a_{52} & a_{53} & a_{55} & a_{56} & b_{52}\\ a_{62} & a_{63} & a_{65} & a_{66} & b_{62}\\ 0 & 0 & 0 & 0 & 0 \end{array}\right] }_{A} \underbrace{ \left[\begin{array}{c} z(t)-z^*\\ \theta(t)-\theta^*\\ \dot{z}(t)\\ \dot{\theta}(t)\\ \theta_e(t)-\theta_e^* \end{array}\right] }_{x(t)}\\ + \underbrace{ \left[\begin{array}{c} 0 \\ 0 \\ 0 \\ 0 \\ \omega_e \end{array}\right] }_{B} u_e(t) + \left[\begin{array}{c} 0 \\ 0 \\ a_{54} \\ a_{64} \\ 0 \end{array}\right] \underbrace{(\dot{x}(t)-V^*)}_{0} + \left[\begin{array}{c} 0 \\ 0 \\ b_{51} \\ b_{61} \\ 0 \end{array}\right]\underbrace{(T(t)-T^*)}_{0} \end{array}}$

$\displaystyle{(14.2)\quad \begin{array}{l} \underbrace{ \left[\begin{array}{c} \theta(t)-\theta^*\\ \dot{\theta}(t)\\ \theta_e(t)-\theta_e^* \end{array}\right] }_{y(t)} = \underbrace{ \left[\begin{array}{ccccc} 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \end{array}\right] }_{C} \underbrace{ \left[\begin{array}{c} z(t)-z^*\\ \theta(t)-\theta^*\\ \dot{z}(t)\\ \dot{\theta}(t)\\ \theta_e(t)-\theta_e^* \end{array}\right] }_{x(t)} \end{array}}$

データセットNo1に対する行列 $A$ の固有値は図3のようにグラフ化されます。ここで、安定な固有値はヒーブの振舞いに、不安定な固有値はピッチの振舞いに関係しています。

図3 データセットNo1を用いた５次元モデルの行列 $A$ の固有値

３次元モデル

$\displaystyle{(15.1)\quad \begin{array}{l} \underbrace{ \frac{d}{dt} \left[\begin{array}{c} \theta(t)-\theta^*\\ \dot{\theta}(t)\\ \theta_e(t)-\theta_e^* \end{array}\right] }_{\dot{x}(t)} = \underbrace{ \left[\begin{array}{ccc} 0 & 1 & 0\\ a_{63} & a_{66} & b_{62}\\ 0 & 0 & 0 \end{array}\right] }_{A} \underbrace{ \left[\begin{array}{c} \theta(t)-\theta^*\\ \dot{\theta}(t)\\ \theta_e(t)-\theta_e^* \end{array}\right] }_{x(t)} + \underbrace{ \left[\begin{array}{c} 0 \\ 0 \\ \omega_e \end{array}\right] }_{B} u_e(t)\\ + \underbrace{\left[\begin{array}{cc} 0 & 0 \\ a_{62} & a_{65} \\ 0 & 0 \end{array}\right] \left[\begin{array}{c} z(t)-z^*\\ \dot{z}(t) \end{array}\right]}_{w(t)} + \left[\begin{array}{cc} 0 & 0 \\ a_{64} & b_{61} \\ 0 & 0 \end{array}\right] \underbrace{\left[\begin{array}{c} \dot{x}(t)-V^* \\ T(t)-T^* \end{array}\right]}_{0} \end{array}}$

$\displaystyle{(15.2)\quad \begin{array}{l} \underbrace{ \left[\begin{array}{c} \theta(t)-\theta^*\\ \dot{\theta}(t)\\ \theta_e(t)-\theta_e^* \end{array}\right] }_{y(t)} = \underbrace{ \left[\begin{array}{ccccc} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{array}\right] }_{C} \underbrace{ \left[\begin{array}{c} \theta(t)-\theta^*\\ \dot{\theta}(t)\\ \theta_e(t)-\theta_e^* \end{array}\right] }_{x(t)} \end{array}}$

データセットNo1に対する行列 $A$ の固有値は図4のようにグラフ化されます。これからピッチに対応する固有値は不安定であることが分かります。

図4 データセットNo1を用いた３次元モデルの行列 $A$ の固有値

●制御目的を達成する制御系を設計する立場から、２つのモデルを比較してみます。まず制御則として次式を考えます。

$\displaystyle{(16)\quad u(t)=-Fx(t)+\int_0^t(\theta_e^*-\theta_e(t))dt }$

ここで、右辺第１項は状態フィードバック、第２項は船外機取付角をある収納状態 $\theta(0)\ne0$ から所定の平衡入力値 $\theta_e^*$ に設定するための積分動作です。

５次元モデル(14.1)の場合は、ヒーブに関する状態変数を計測できないので、状態オブザーバを用いてこれらを推定して、状態フィードバックを実施することになります。

一方３次元モデル(15.1)の場合は、状態フィードバックは実施できますが、ヒーブの振舞いの影響が外乱 $w(t)$ として現れますので、これを抑制する必要があります。そのヒーブの振舞いは次式で表されます。

$\displaystyle{(17)\quad \begin{array}{l} \underbrace{ \frac{d}{dt} \left[\begin{array}{c} z(t)-z^*\\ \dot{z}(t) \end{array}\right] }_{\dot{x}'(t)} = \underbrace{ \left[\begin{array}{ccccc} 0 & 1 \\ a_{52} & a_{55} \end{array}\right] }_{A'} \underbrace{ \left[\begin{array}{c} z(t)-z^*\\ \dot{z}(t)\\ \end{array}\right] }_{x'(t)}\\ + \underbrace{ \left[\begin{array}{ccccc} 0 & 1 & 0\\ a_{53} & a_{56} & b_{52} \end{array}\right] }_{B'} \underbrace{ \left[\begin{array}{c} \theta(t)-\theta^*\\ \dot{\theta}(t)\\ \theta_e(t)-\theta_e^* \end{array}\right] }_{y(t)}\\ + \left[\begin{array}{c} 0 \\ a_{54} \end{array}\right] \underbrace{(\dot{x}(t)-V^*)}_{0} + \left[\begin{array}{c} 0 \\ b_{51} \end{array}\right]\underbrace{(T(t)-T^*)}_{0} \end{array}}$

これからヒーブの振舞いには操作入力 $u_e$ は直接影響を及ぼさず、ピッチを通して間接的に影響すること分かります。ただ、これは漸近安定（ $A'$ は安定行列）なので、ピッチの安定化を速やかにできれば（ $y(t)$ が零に収束すれば）、(17)は考慮しなくともよいことになります。しかしながら、３次元モデルに対する状態フィードバックを５次元モデルに適用して、閉ループ系が漸近安定かどうかを調べておくことが考えられます。

制御系設計に３次元モデルを用いる場合は、状態オブザーバが不要となる利点がありますが、ヒーブの振舞いの影響が抑制されていることを確かめるために、５次元モデルに対する閉ループシミュレ－ションを行うことが重要になります。

●以上の結果を、 $k=128$ の場合の平衡状態回りの線形モデルについて数値で確認しておきます。

まず元の線形モデル(12)の行列 $A$ とその固有値は次のようになります。

$\displaystyle{(18.1)\quad A= \left[\begin{array}{cccccc} 0 & 0 & 0 & 1.0000 & 0 & 0\\ 0 & 0 & 0 & 0 & 1.0000 & 0\\ 0 & 0 & 0 & 0 & 0 & 1.0000\\ 0 & 8.3458 & -6.9684 & -0.2549 & 0 & 0\\ 0 & -58.4515 & 67.9416 & 0.4686 & -0.7383 & -0.0369\\ 0 & 0.4559 & -30.3422 & -0.1922 & -0.0117 & 0.1175 \end{array}\right] }$

$\displaystyle{(18.2)\quad \left\{\begin{array}{lr} \lambda_1=& 0.0000 + 0.0000i\\ \lambda_2=& -0.2052 + 0.0000i\\ \lambda_3=& -0.4030 + 7.7086i\\ \lambda_4=& -0.4030 - 7.7086i\\ \lambda_5=& 0.0678 + 5.4080i\\ \lambda_6=& 0.0678 - 5.4080i \end{array}\right.}$

５次元モデル(14.1)における行列 $A$ とその固有値は次のようになります。

$\displaystyle{(19.1)\quad A= \left[\begin{array}{ccccc} 0 & 0 & 1.0000 & 0 & 0\\ 0 & 0 & 0 & 1.0000 & 0\\ -58.4515 & 67.9416 & -0.7383 & -0.0369 & -0.0797\\ 0.4559 & -30.3422 & -0.0117 & 0.1175 & -2.7873\\ 0 & 0 & 0 & 0 & 0 \end{array}\right] }$

$\displaystyle{(19.2)\quad \left\{\begin{array}{lrc} \lambda_1^{(5)}=& -0.3412 + 7.7012i&\approx\lambda_3\\ \lambda_2^{(5)}=& -0.3412 - 7.7012i&\approx\lambda_4\\ \lambda_3^{(5)}=& 0.0308 + 5.4151i&\approx\lambda_5\\ \lambda_4^{(5)}=& 0.0308 - 5.4151i&\approx\lambda_6\\ \lambda_5^{(5)}=& 0.0000 + 0.0000i& \end{array}\right.}$

これからサージに関わる状態変数を除いた場合、ヒーブとピッチに関わる固有値 $\lambda_3,\lambda_4,\lambda_5,\lambda_6$ がほぼ引き継がれることが分かります。

一方３次元モデル(15.1)における行列 $A$ とその固有値は次のようになります。

$\displaystyle{(20.1)\quad A= \left[\begin{array}{ccc} 0 & 1.0000 & 0\\ -30.3422 & 0.1175 & -2.7873\\ 0 & 0 & 0 \end{array}\right] }$

$\displaystyle{(20.2)\quad \left\{\begin{array}{lrc} \lambda_1^{(3)}=& 0.0587 + 5.5081i&\approx\lambda_5\\ \lambda_2^{(3)}=& 0.0587 - 5.5081i&\approx\lambda_6\\ \lambda_3^{(3)}=& 0.0000 + 0.0000i&\\ \end{array}\right.}$

これからさらにヒーブに関わる状態変数を除いた場合、ピッチに関わる固有値 $\lambda_5,\lambda_6$ がほぼ引き継がれることが分かります。

●ヒーブとピッチに関わる固有値 $\lambda_3,\lambda_4,\lambda_5,\lambda_6$ の可制御性と可観測性を、５次元モデルを用いて調べてみます。

$\displaystyle{(21)\quad \left\{\begin{array}{l} \underline{\sigma}(\left[\begin{array}{cc} B & A-\lambda_1^{(5)} I_5 \end{array}\right])=0.0048^{*}\\ \underline{\sigma}(\left[\begin{array}{cc} B & A-\lambda_2^{(5)} I_5 \end{array}\right])=0.0048^{*}\\ \underline{\sigma}(\left[\begin{array}{cc} B & A-\lambda_3^{(5)} I_5 \end{array}\right])=0.0106\\ \underline{\sigma}(\left[\begin{array}{cc} B & A-\lambda_4^{(5)} I_5 \end{array}\right])=0.0106\\ \underline{\sigma}(\left[\begin{array}{cc} B & A-\lambda_5^{(5)} I_5 \end{array}\right])=0.1140\\ \end{array}\right. }$

$\displaystyle{(22)\quad \left\{\begin{array}{l} \underline{\sigma}(\left[\begin{array}{cc} C^T & A^T-\lambda_1^{(5)} I_5 \end{array}\right])=0.0146^{*}\\ \underline{\sigma}(\left[\begin{array}{cc} C^T & A^T-\lambda_2^{(5)} I_5 \end{array}\right])=0.0146^{*}\\ \underline{\sigma}(\left[\begin{array}{cc} C^T & A^T-\lambda_3^{(5)} I_5 \end{array}\right])=0.3129\\ \underline{\sigma}(\left[\begin{array}{cc} C^T & A^T-\lambda_4^{(5)} I_5 \end{array}\right])=0.3129\\ \underline{\sigma}(\left[\begin{array}{cc} C^T & A^T-\lambda_5^{(5)} I_5 \end{array}\right])=0.9940\\ \end{array}\right. }$

(21)から、すべての固有値について行列 $\left[\begin{array}{cc} B & A-\lambda I_5 \end{array}\right]$ の最小特異値は正ですから、形式的には可制御性は成り立つと判定できます。ただ、＊印の値が示すように、ヒーブはピッチに比べて相対的に可制御性の程度が弱く、状態フィードバックによる固有値の変更が難しいと言えます。したがって、制御対象は可制御でなく、むしろ可安定と判定することが考えられます。

一方、(22)から、形式的には可観測性は成り立つと判定できます。ただ、＊印の値が示すように、ヒーブはピッチに比べて相対的に可観測性の程度が弱く、状態オブザーバによる推定が難しいと言えます。

●操作入力（取付角）により、まずピッチが制御され、それがヒーブに影響を与えるという連成の仕組みを調べるために、(15.1)と(17)を数値で確かめてみます。

まず(15.1)は

$\displaystyle{(15.1')\quad \begin{array}{l} \underbrace{ \frac{d}{dt} \left[\begin{array}{c} \theta(t)-\theta^*\\ \dot{\theta}(t)\\ \theta_e(t)-\theta_e^* \end{array}\right] }_{\dot{x}(t)} = \underbrace{ \left[\begin{array}{ccc} 0 & 1 & 0\\ %a_{63} & a_{66} & b_{62}\\ -30.3422 & 0.1175 & -2.7873\\ 0 & 0 & 0 \end{array}\right] }_{A} \underbrace{ \left[\begin{array}{c} \theta(t)-\theta^*\\ \dot{\theta}(t)\\ \theta_e(t)-\theta_e^* \end{array}\right] }_{x(t)}\\ + \underbrace{ \left[\begin{array}{c} 0 \\ 0 \\ \omega_e \end{array}\right] }_{B} u_e(t) + \underbrace{\left[\begin{array}{cc} 0 & 0 \\ %a_{62} & a_{65} \\ 0.4559 & -0.0117 \\ 0 & 0 \end{array}\right] \left[\begin{array}{c} z(t)-z^*\\ \dot{z}(t) \end{array}\right]}_{w(t)} \end{array}}$

となり、ヒーブのピッチへの影響を表す $w(t)$ に関わる係数 $a_{62}=0.4559,a_{65}=-0.017$ の値がそれほど大きくなく、３次元モデルにおいても固有値が継承されることに注意します。

一方(17)は

$\displaystyle{(17')\quad \begin{array}{l} \underbrace{ \frac{d}{dt} \left[\begin{array}{c} z(t)-z^*\\ \dot{z}(t) \end{array}\right] }_{\dot{x}'(t)} = \underbrace{ \left[\begin{array}{ccccc} 0 & 1 \\ %a_{52} & a_{55} -58.4515 & -0.7383 \end{array}\right] }_{A'} \underbrace{ \left[\begin{array}{c} z(t)-z^*\\ \dot{z}(t)\\ \end{array}\right] }_{x'(t)}\\ + \underbrace{ \left[\begin{array}{ccccc} 0 & 1 & 0\\ %a_{53} & a_{56} & b_{52} 67.9416 & -0.0369 & -0.0797 \end{array}\right] }_{B'} \underbrace{ \left[\begin{array}{c} \theta(t)-\theta^*\\ \dot{\theta}(t)\\ \theta_e(t)-\theta_e^* \end{array}\right] }_{y(t)} \end{array}}$

となります。これは漸近安定なので、ピッチの安定化が速やかに行えれば問題ないとと言えます。

●最後に、取付角からの伝達特性を調べておきます。

図5 取付角からヒーブまで（左図）とピッチまで（右図）の伝達特性

左図は２つの共振特性をもちますが、右図は（２つ目がほとんど消えていて）１つだけの共振特性をもちます。これから、ヒーブのピッチへの連成影響はほとんどないことがわかります。

●以上の数値例の結果を得るためのプログラムを以下に示します。

MATLAB

%ymcpb_ana
%-----
 clear all, close all
 dataset=[...
 "ABC_No1_20220128",...
 "ABC_No2_20220131",...
 "ABC_No3_20220128",...
 "ABC_No4_20220131",...
 "ABC_No5_20220131"];
 no=1
 load(dataset(no))
 ID=[8,9,10,11,12,13,14,...
 36,37,38,39,40,41,42,43,...
 66,67,68,69,70,71,72,73,...
 97,98,99,100,101,102,103,...
 126,127,128,129,130,131,132,133];
 id=1:length(ID);
%-----
 figure(1)
 subplot(321),plot(run_eq_Series,"o-"), hold on
 subplot(321),plot(ID,run_eq_Series(ID),"or"),title("No"+num2str(no)+": "+"run")
 subplot(323),plot(No_eq_Series,"o-"), hold on
 subplot(323),plot(ID,No_eq_Series(ID),"or"),title("No"+num2str(no)+": "+"no")
 subplot(325),plot(dx_eq_Series*3.6,"o-"), hold on
 subplot(325),plot(ID,dx_eq_Series(ID)*3.6,"or"),title("No"+num2str(no)+": "+"dxs")
 subplot(322),plot(z_eq_Series*100,"o-"), hold on
 subplot(322),plot(ID,z_eq_Series(ID)*100,"or"),title("No"+num2str(no)+": "+"zs")
 subplot(324),plot(th_eq_Series/pi*180,"o-"), hold on
 subplot(324),plot(ID,th_eq_Series(ID)/pi*180,"or"),title("No"+num2str(no)+": "+"ths")
 subplot(326),plot(the_eq_Series/pi*180,"o-"), hold on
 subplot(326),plot(ID,the_eq_Series(ID)/pi*180,"or"),title("No"+num2str(no)+": "+"thes")
%-----
 figure(2)
 for i=id, a41(i)=A0(4,1,ID(i)); end
 subplot(331),plot(id,a41,"o-",33,a41(33),"*r"),title("No"+num2str(no)+": "+"a41")
 for i=id, a42(i)=A0(4,2,ID(i)); end
 subplot(332),plot(id,a42,"o-",33,a42(33),"*r"),title("No"+num2str(no)+": "+"a42")
 for i=id, a43(i)=A0(4,3,ID(i)); end
 subplot(333),plot(id,a43,"o-",33,a43(33),"*r"),title("No"+num2str(no)+": "+"a43")
 for i=id, a51(i)=A0(5,1,ID(i)); end
 subplot(334),plot(id,a51,"o-",33,a51(33),"*r"),title("No"+num2str(no)+": "+"a51")
 for i=id, a52(i)=A0(5,2,ID(i)); end
 subplot(335),plot(id,a52,"o-",33,a52(33),"*r"),title("No"+num2str(no)+": "+"a52")
 for i=id, a53(i)=A0(5,3,ID(i)); end
 subplot(336),plot(id,a53,"o-",33,a53(33),"*r"),title("No"+num2str(no)+": "+"a53")
 for i=id, a61(i)=A0(6,1,ID(i)); end
 subplot(337),plot(id,a61,"o-",33,a61(33),"*r"),title("No"+num2str(no)+": "+"a61")
 for i=id, a62(i)=A0(6,2,ID(i)); end
 subplot(338),plot(id,a62,"o-",33,a62(33),"*r"),title("No"+num2str(no)+": "+"a62")
 for i=id, a63(i)=A0(6,3,ID(i)); end
 subplot(339),plot(id,a63,"o-",33,a63(33),"*r"),title("No"+num2str(no)+": "+"a63")
%-----
 figure(3)
 for i=id, a44(i)=A0(4,4,ID(i)); end
 subplot(331),plot(id,a44,"o-",33,a44(33),"*r"),title("No"+num2str(no)+": "+"a44")
 for i=id, a45(i)=A0(4,5,ID(i)); end
 subplot(332),plot(id,a45,"o-",33,a45(33),"*r"),title("No"+num2str(no)+": "+"a45")
 for i=id, a46(i)=A0(4,6,ID(i)); end
 subplot(333),plot(id,a46,"o-",33,a46(33),"*r"),title("No"+num2str(no)+": "+"a46")
 for i=id, a54(i)=A0(5,4,ID(i)); end
 subplot(334),plot(id,a54,"o-",33,a54(33),"*r"),title("No"+num2str(no)+": "+"a54")
 for i=id, a55(i)=A0(5,5,ID(i)); end
 subplot(335),plot(id,a55,"o-",33,a55(33),"*r"),title("No"+num2str(no)+": "+"a55")
 for i=id, a56(i)=A0(5,6,ID(i)); end, axis([0 40 -1 1])
 subplot(336),plot(id,a56,"o-",33,a56(33),"*r"),title("No"+num2str(no)+": "+"a56")
 for i=id, a64(i)=A0(6,4,ID(i)); end, axis([0 40 -1 1])
 subplot(337),plot(id,a64,"o-",33,a64(33),"*r"),title("No"+num2str(no)+": "+"a64")
 for i=id, a65(i)=A0(6,5,ID(i)); end, axis([0 40 -1 1])
 subplot(338),plot(id,a65,"o-",33,a65(33),"*r"),title("No"+num2str(no)+": "+"a65")
 for i=id, a66(i)=A0(6,6,ID(i)); end
 subplot(339),plot(id,a66,"o-",33,a66(33),"*r"),title("No"+num2str(no)+": "+"a66")
%-----
 figure(4)
 for i=id, b41(i)=B0(4,1,ID(i)); end
 subplot(321),plot(id,b41,"o-",33,b41(33),"*r"),title("No"+num2str(no)+": "+"b41")
 for i=id, b42(i)=B0(4,2,ID(i)); end
 subplot(322),plot(id,b42,"o-",33,b42(33),"*r"),title("No"+num2str(no)+": "+"b42")
 for i=id, b51(i)=B0(5,1,ID(i)); end
 subplot(323),plot(id,b51,"o-",33,b51(33),"*r"),title("No"+num2str(no)+": "+"b51")
 for i=id, b52(i)=B0(5,2,ID(i)); end
 subplot(324),plot(id,b52,"o-",33,b52(33),"*r"),title("No"+num2str(no)+": "+"b52")
 for i=id, b61(i)=B0(6,1,ID(i)); end
 subplot(325),plot(id,b61,"o-",33,b61(33),"*r"),title("No"+num2str(no)+": "+"b61")
 for i=id, b62(i)=B0(6,2,ID(i)); end
 subplot(326),plot(id,b62,"o-",33,b62(33),"*r"),title("No"+num2str(no)+": "+"b62")
%-----６次元モデル
%(1)x,(2)z,(3)th,(4)dx-V*,(5)dz,(6)dth
 sys=ss(A0(:,:,ID),B0(:,:,ID),eye(6,6),[]);
 figure(5),pzmap(sys),axis([-10,10,-10,10])
 title("No"+num2str(no)+": "+"k=ID"),hold on
 figure(5),pzmap(sys(:,:,33),"*r")
%-----５次元モデル
%(1)z,(2)th,(3)dz,(4)dth,(5)the
 k=[2,3,5,6]; i=[2,4,5]; j=[2];
 for id=ID
   A(:,:,id)=[A0(k,k,id) B0(k,j,id);zeros(1,5)];
 end 
 B=[zeros(4,1);0.0380*3];
 C=eye(5,5);
 C=C(i,:);
 sys5=ss(A(:,:,ID),B,C,[]);  
 figure(6),pzmap(sys5),axis([-10,10,-10,10])
 title("No"+num2str(no)+": "+"k=ID"),hold on
 figure(6),pzmap(sys5(:,:,33),"*r")
%-----３次元モデル
%(1)th,(2)dth,(3)the
 k=[3,6]; j=[2];
 for id=ID
   AA(:,:,id)=[A0(k,k,id) B0(k,j,id);zeros(1,3)];
 end 
 B=[zeros(2,1);0.0380*3];
 C=eye(3,3);
 sys3=ss(AA(:,:,ID),B,C,[]);  
 figure(7),pzmap(sys3),axis([-10,10,-10,10])
 title("No"+num2str(no)+": "+"k=ID"),hold on
 figure(7),pzmap(sys3(:,:,33),"*r")
%-----設計ポイント id=128,ID=33
 [Asys,Bsys]=ssdata(sys(:,:,33))
 [Asys5,Bsys5,Csys5]=ssdata(sys5(:,:,33))
 [Asys3,Bsys3,Csys3]=ssdata(sys3(:,:,33))
 pl=pole(sys(:,:,33))
 pl5=pole(sys5(:,:,33))
 pl3=pole(sys3(:,:,33))
%-----可制御性 
 w=zeros(5,1);
 for i=1:5
   S=[Bsys5 Asys5-pl5(i)*eye(5,5)];
   w(i)=min(svd(S));
 end
 C0=[pl5 w]
%-----可観測性 
 w=zeros(5,1);
 for i=1:5
   S=[Csys5; Asys5-pl5(i)*eye(5,5)];
   w(i)=min(svd(S));
 end
 OB=[pl5 w]
%-----伝達特性
%(1)x,(2)z,(3)th,(4)dx-V*,(5)dz,(6)dth
 figure(11),hold on
 sigma(sys(:,2,33),logspace(0,2))
 title("No"+num2str(no)+": "+"sys(:,2,33)")
%-----
%(2)z,(5)dz
 figure(12),hold on
 sigma(sys([2,5],2,33),logspace(0,2))
 title("No"+num2str(no)+": "+"sys([2,5],2,33)")
%-----
%(3)th,(6)dth
 figure(13),hold on
 sigma(sys([3,6],2,33),logspace(0,2))
 title("No"+num2str(no)+": "+"sys([3,6],2,33)")
%-----
%eof

参考：漸近安定性

【漸近安定性の定義とその等価な条件】
定義DA:　 $\forall x(0)\ne 0: x(t)=\exp(At)x(0)\rightarrow 0\quad(t\rightarrow\infty)$
条件A1:　 ${\rm Re}(\lambda_i(A))<0\quad(i=1,\cdots,n)$

漸近安定性とは、定義DAより、平衡状態が乱されたとき復帰できるかどうかにかかわる概念で、条件A1より、 $A$ 行列の固有値がすべて複素左半面にあることを調べて判定されます。

参考：可制御性・可安定性

【可安定性の定義とその等価な条件】
定義DS:　状態フィードバックにより安定化可能
条件S1:　 ${\rm rank}\, \left[\begin{array}{cc} B & A-\lambda I_n \end{array}\right] =n$ 　( $\lambda$ は $A$ のすべての不安定固有値)

【可制御性の定義とその等価な条件】
定義DC:　任意初期状態を，任意有限時間内に，任意状態に移動可能
条件C2:　 ${\rm rank}\, \left[\begin{array}{cccc} B & AB & \cdots & A^{n-1}B \end{array}\right]=n$
条件C3:　 $F$ を選んで， $A-BF$ の固有値を任意に設定可能
条件C4:　 $\boxed{{\rm rank}\, \left[\begin{array}{cc} B & A-\lambda I_n \end{array}\right] =n}$ 　( $\lambda$ は $A$ のすべての固有値)

可安定性は、定義DSより、状態フィードバックにより安定化可能であることを意味し、条件S1が判定条件として知られています。一方、可制御性については、条件C4が判定条件として知られています。両者の相違は、すでに左半平面にある安定な固有値について関与（再配置）するかどうかにあります。可安定性の判定では不安定な固有値のみで条件S1を判定しますが、可制御性の判定ではすべて固有値について条件C4を判定します。

それでは、 $\lambda_1,\cdots,\lambda_n$ に対して、条件S1と条件C4に出てくる行列

$\displaystyle{(3)\quad M(\lambda_i)=\left[\begin{array}{cc} B & A-\lambda_i I_n \end{array}\right] }$

の階数をどのようにして計算するかですが、これは $M$ の非零特異値の数で決定します。したがって、条件S1と条件C4は、最小特異値が正かどうかで判定します。

参考：可観測性・可検出性

【可検出性の定義とその等価な条件】
定義DD:　状態オブザーバを構成可能
条件D1:　 ${\rm rank}\, \left[\begin{array}{c} C \\ A-\lambda I_n \end{array}\right]=n$ 　( $\lambda$ は $A$ のすべての不安定固有値)

【可観測性の定義とその等価な条件】
定義DO:　任意有限時間 $[0,t]$ 上の入出力データ $u(\tau),y(\tau)\tau\in[0,t]$ から，初期状態 $x(0)$ を一意に決定できる
条件O2:　 ${\rm rank}\, \left[\begin{array}{c} C \\ CA \\ \vdots\\ CA^{n-1} \end{array}\right] =n$
条件O3:　 $H$ を選んで， $A-HC$ の固有値を任意に設定可能
条件O4:　 $\boxed{{\rm rank} \left[\begin{array}{c} C \\ A-\lambda I_n \end{array}\right]=n}$ 　( $\lambda$ は $A$ のすべての固有値)

可検出性は、定義DDより、状態オブザーバを構成可能であることを意味し、条件D1が判定条件として知られています。状態オブザーバの $A$ 行列は $A-HC$ と表され、オブザーバゲイン $H$ は、仮想システム $\dot{v}=A^Tv+C^Tw$ に対する状態ＦＢ $w=-H^Tv$ による閉ループ系 $\dot{v}=(A-HC)^Tv$ を安定化して決定します。したがって、可検出性は $\dot{v}=A^Tv+C^Tw$ の可安定性を意味します。一方、可観測性については、条件O4が判定条件として知られています。可検出性との相違は、すでに左半平面にある安定な $A$ の固有値について $A-HC$ で再配置するかどうかにあります。可検出性の判定では不安定な固有値のみで条件D1を判定しますが、可観測性の判定ではすべて固有値について条件O4を判定します。

それでは、 $\lambda_1,\cdots,\lambda_n$ に対して、条件D1と条件O4に出てくる行列

$\displaystyle{(4)\quad N(\lambda_i)=\left[\begin{array}{c} C \\ A-\lambda I_n \end{array}\right] }$

の階数をどのようにして計算するかですが、これは $N$ の非零特異値の数で決定します。したがって、条件D1と条件O4は、最小特異値が正かどうかで判定します。

参考：状態オブザーバ

●通常の状態オブザーバは

$\displaystyle{(7)\quad \dot{\hat{x}}(t)=(A-HC)\hat{x}(t)+Hy(t)+Bu(t) }$

で表され

$\displaystyle{(8)\quad \hat{x}(t)\rightarrow x(t)\quad(t\rightarrow\infty) }$

のように漸近的に元の状態を推定するものです。行列 $H$ はオブザーバゲインと呼ばれ、 $A-HC$ は安定行列とするように選ばれます。

●最小次元状態オブザーバは

$\displaystyle{(9)\quad \left\{\begin{array}{lll} \dot{\hat{z}}(t)=\hat{A}\hat{z}(t)+\hat{B}y(t)+\hat{J}u(t),\ \hat{z}(0)=0\\ \hat{x}(t)=\hat{C}\hat{z}(t)+\hat{D}y(t)\\ (\hat{z}(t)\in{\rm\bf R}^{n-p}, y(t)\in{\rm\bf R}^p, u(t)\in{\rm\bf R}^m) \end{array}\right. }$

によって表される漸近安定なシステムで、この出力が $x(t)$ に漸近するように構成するものです。そのためには、あるサイズ $n-p \times n$ の行列 $U$ に対して次式を満足させることが必要十分となります。

$\displaystyle{ \begin{array}{cl} (10.1) & \hat{A}U+\hat{B}C=UA\\ (10.2) & \hat{C}U+\hat{D}C=I_n\\ (10.3) & \hat{J}=UB \end{array} }$

参考：最小実現

一般に，不可制御かつ不可観測な状態空間表現は適当な座標変換により，つぎの正準構造をもつように変換できることが知られています。

$\displaystyle{(1)\quad \left[\begin{array}{c|c} TAT^{-1} & TB \\\hline CT^{-1} & D \end{array}\right] = \left[\begin{array}{cccc|c} A_1 & 0 & X_{13} & 0 & B_1 \\ X_{21} & A_2 & X_{23} & X_{24} & B_2 \\ 0 & 0 & A_3 & 0 & 0 \\ 0 & 0 & X_{43} & A_4 & 0 \\\hline C_1 & 0 & C_3 & 0 & 0 \end{array}\right] }$