船舶のLPV制御 | 制御系CAD

ship1

●船舶の針路制御のために、次のNomotoモデルを考えます。

$\displaystyle{ \dot{r}(t)=-\frac{1}{T(t)}r(t)+\frac{K(t)}{T(t)}\delta(t-L)+w(t)} }$

ただし、 $r(t)=\dot{\psi}(t)$ かつ

$\displaystyle{ T(t)=\frac{L}{U(t)}T',\ K(t)=\frac{U(t)}{L}K' \quad (U_1\le U(t)\le U_2) }$

ここで、舵の効果が出るまでの無駄時間 $L$ を想定しています。混同の恐れはないので船長 $L$ と同じ表記を用います。また前進速度 $U$ の変動により、時定数 $T$ と定常ゲイン $K$ が変動しますが、次式から舵は前進速度の２乗で効いてくることがわかります。

$\displaystyle{ \dot{r}(t)=-\underbrace{\left(\frac{U(t)}{L}\right)\frac{1}{T'}}_{a(t)=\frac{1}{T(t)}}r(t) +\underbrace{\left(\frac{U(t)}{L}\right)^2\frac{K'}{T'}}_{b(t)=\frac{K(t)}{T(t)}}\delta(t-L)+w(t) }$

いま、ある前進速度 $U^*$ に注目しますと、次の表現もできます。

$\displaystyle{ \dot{r}(t)=-\underbrace{\left(\frac{U(t)}{U^*}\right)\frac{1}{T^*}}_{a(t)=\frac{1}{T(t)}}r(t) +\underbrace{\left(\frac{U(t)}{U^*}\right)^2\frac{K^*}{T^*}}_{b(t)=\frac{K(t)}{T(t)}}\delta(t-L)+w(t) }$

ただし

$\displaystyle{ T^*=\frac{L}{U^*}T',\ K^*=\frac{U^*}{L}K' \quad (U_1\le U^*\le U_2) }$

●【Z試験】　舵の無駄時間および前進速度の変動は無視できるとして、いわゆるＺ試験のシミュレーションを行ってみます。

ship1a

この応答から時定数 $T$ と定常ゲイン $K$ を求めることがよく行われていますが、線形システムを特徴付けるインパルス応答またはステップ応答を求めているわけではありません。背景には、船舶は無定位系のためステップ応答を求めることができない事情があります。

ship2

しかしながら、上のような単位フィードバックを行えば、ステップ応答が求まり、その第１番目のオーバーシュートの座標から、次のように時定数と定常ゲインを求めることができます。

$\displaystyle{ (T_p, 1+p_0)\ \Rightarrow\ (\zeta',\omega_n')= \left( \sqrt{\frac{(\log p_0')^2}{(\log p_0')^2+\pi^2}}, \frac{\sqrt{(\log p_0')^2+\pi^2}}{T_p'} \right) }$

ship2

●【PID制御】いま、ある前進速度 $U^*$ に注目し、無駄時間を無視した、次のNOMOTOモデルを考えます。

$\displaystyle{ \dot{r}(t)=-\frac{1}{T^*}r(t)+\frac{K^*}{T^*}\delta(t)+w(t) }$

このとき、PID制御は次式で表されます。

$\displaystyle{ \delta(t)=K_P(\psi_c-\psi(t))+K_D\frac{d}{dt}(\psi_c-\psi(t))+K_I\int_0^t(\psi_c-\psi(\tau))\,d\tau }$

通常、Pゲイン $K_P$ 、Dゲイン $K_D$ 、Iゲイン $K_I$ は、現場でマニュアルで調整されますが、ここではNOMOTOモデルに基づいて決定してみます。

PID制御による閉ループ系は次のように表されます。

ship3

ここで、 $\dot{\psi}_c=0$ を仮定すると

$\displaystyle{ \delta(t)=-K_P\psi(t)-K_Dr(t)+K_P\psi_c+K_I\int_0^t(\psi_c-\psi(\tau))\,d\tau }$

から、PID制御による閉ループ系は次のように表されます。

ship4

いま、 $\zeta=\frac{1}{2\sqrt{T^*K^*}},\omega_n=\sqrt{\frac{K^*}{T^*}}$ とおいて、次の状態方程式を得ます。

$\displaystyle{ \underbrace{ \left[\begin{array}{c} \dot{\psi}(t) \\ \dot{r}(t) \end{array}\right] }_{\dot{x}(t)} = \underbrace{ \left[\begin{array}{cc} 0 & 1 \\ 0 & -2\zeta\omega_n \end{array}\right] }_{A} \underbrace{ \left[\begin{array}{c} \psi(t) \\ r(t) \end{array}\right] }_{x(t)} + \underbrace{ \left[\begin{array}{c} 0 \\ \omega_n^2 \end{array}\right] }_{B} \underbrace{\delta(t)}_{u(t)} + \left[\begin{array}{c} 0 \\ w \end{array}\right] }$

このとき、PID制御は、次のように、フィードフォワード項と積分動作をもつ状態フィードバックとなります。

$\displaystyle{ \underbrace{\delta(t)}_{u(t)} =- \underbrace{ \left[\begin{array}{cc} K_P & K_D \end{array}\right] }_{F} \underbrace{ \left[\begin{array}{c} \psi(t) \\ r(t) \end{array}\right] }_{x(t)} +\underbrace{K_P}_{G} \underbrace{\psi_c}_{v} +K_I \underbrace{\int_0^t(\psi_c-\psi(\tau))\,d\tau}_{x_I(t)} }$

そこで、望ましい $(T_p',1+p_0')$ を指定して、対応する $(\zeta',\omega_n')$ を求め、次式からＰＩＤゲインを決めます。

$\displaystyle{ K_P=\frac{\omega_n'^2}{\omega_n^2},\ K_D=\frac{2}{\omega_n^2}(\zeta'\omega_n'-\zeta\omega_n),\ K_I=\omega_I\frac{\omega_n'^2}{\omega_n^2} }$

ただし、 $\omega_I=\frac{w_n'}{10}$ を目安とし、外乱抑制の効果を見ながら調整します。

このとき、PID制御による閉ループ系は、次式で表されます。

$\displaystyle{ \underbrace{ \left[\begin{array}{c} \dot{\psi}(t) \\ \dot{r}(t) \\ \dot{x}_I(t) \end{array}\right] }_{\dot{x}_{CL}(t)} = \underbrace{ \left[\begin{array}{ccc} 0 & 1 & 0\\ -\omega_n'^2 & -2\zeta'\omega_n' & \omega_I\omega_n'^2 \\ -1 & 0 & 0 \end{array}\right] }_{A_{CL}} \underbrace{ \left[\begin{array}{c} \psi(t) \\ r(t) \\ x_I(t) \end{array}\right] }_{x_{CL}(t)} + \underbrace{ \left[\begin{array}{cc} 0 & 0\\ \omega_n'^2 & 1\\ 1 & 0 \end{array}\right] }_{B_{CL}} \left[\begin{array}{c} \psi_c \\ w \end{array}\right] }$

ここで、定値外乱 $\dot{w}=0$ を仮定し、両辺を微分して次式を得ます。

$\displaystyle{ \underbrace{ \left[\begin{array}{c} \ddot{\psi}(t) \\ \ddot{r}(t) \\ \ddot{x}_I(t) \end{array}\right] }_{\ddot{x}_{CL}(t)} = \underbrace{ \left[\begin{array}{ccc} 0 & 1 & 0\\ -\omega_n'^2 & -2\zeta'\omega_n' & \omega_I\omega_n'^2 \\ -1 & 0 & 0 \end{array}\right] }_{A_{CL}} \underbrace{ \left[\begin{array}{c} \dot{\psi}(t) \\ \dot{r}(t) \\ \dot{x}_I(t) \end{array}\right] }_{\dot{x}_{CL}(t)} }$

閉ループ系の漸近安定性から次式を得ます。

$\displaystyle{ \underbrace{ \left[\begin{array}{c} \psi(t) \\ r(t) \\ x_I(t) \end{array}\right] }_{x_{CL}} \rightarrow \underbrace{ \left[\begin{array}{ccc} 0 & -1 & 0\\ \omega_n'^2 & 2\zeta'\omega_n' & -\omega_I\omega_n'^2 \\ 1 & 0 & 0 \end{array}\right]^{-1} }_{-A_{CL}^{-1}} \underbrace{ \left[\begin{array}{cc} 0 \\ \omega_n'^2\psi_c +w\\ \psi_c \end{array}\right] }_{B_{CL}} = \left[\begin{array}{cc} \psi_c\\ 0\\ \frac{w}{\omega_I\omega_n'^2} \end{array}\right] }$

これは、積分動作が、未知の定値外乱を推定していることを意味しています。

それでは、まず、あるＰＤ制御による閉ループ系のシミュレーション結果を次に示します。これから、Ｉ動作がないため、定値外乱の影響を除去できていないことが分かります。

ship3

次に、あるＰＩＤ制御による閉ループ系のシミュレーション結果を示します。これは $\omega_I=2\frac{w_n'}{10}$ を使用しています。定常的に外乱は抑制されていますが、過渡的には当初の狙いから外れて少し揺れ気味になっています。

ship4

●【LQI制御】ＰＩＤゲインをLQI制御の枠組みで求めてみます。そのために、次の偏差系に注目します。

$\displaystyle{ \left[\begin{array}{c} \dot\psi(t)-\dot{\psi}_c \\ \dot r(t) \\ \dot{\delta}(t) \end{array}\right] = \left[\begin{array}{ccc} 0 & 1 & 0 \\ 0 & -2\zeta\omega_n &\omega_n^2 \\ 0 & 0 & 0 \end{array}\right] \left[\begin{array}{c} \psi(t)-\psi_c \\ r(t) \\ \delta(t) \end{array}\right] + \left[\begin{array}{c} 0 \\ 0 \\ 1 \end{array}\right] \dot{\delta} (t) }$

これに対して、次の２次形式評価関数を最小にします。

$\displaystyle{ J=\int_0^\infty(\frac{1}{M_\psi^2}(\psi(t)-\psi_c)^2+\frac{T_\psi^2}{M_\psi^2}r^2(t) +\frac{1}{M_\delta^2}\delta^2(t)+\frac{T_\delta^2}{M_\delta^2}\dot{\delta}^2(t))\,dt }$

これから、次のようにＰＩＤゲインが得られます。

$\displaystyle{ \dot{\delta}(t) =- \left[\begin{array}{ccc} f_\psi & f_r & f_\delta \end{array}\right] \left[\begin{array}{c} \psi(t)-\psi_c \\ r(t) \\ \delta(t) \end{array}\right] }$
$\displaystyle{ =- \underbrace{ \left[\begin{array}{ccc} f_\psi & f_r & f_\delta \end{array}\right] \left[\begin{array}{ccc} 0 & 1 & 0 \\ 0 & -2\zeta\omega_n &\omega_n^2 \\ 1 & 0 & 0 \end{array}\right]^{-1} }_{ \left[\begin{array}{ccc} K_P & K_D & K_I \end{array}\right] } \left[\begin{array}{c} \dot\psi(t)-\dot{\psi}_c \\ \dot{r}(t) \\ \psi(t)-\psi_c \end{array}\right] }$

$M_{th}=1, T_{th}=0.5T, M_i=0.5, T_i=0.1T$ の場合のシミュレーション結果を次に示します。PID制御と違うのは、ここではフィードフォワード項を入れていない点です。

ship5

しかしながら、無駄時間 $L=9$ を入れてみると、閉ループ系は発散してしまいます。

ship6

●【LPV制御】

ship5

いま、次のような前進速度の変動を考えます。このとき、単位フィードバックによる閉ループ系は大きく変動します。

ship7

ある前進速度 $U^*$ を基準にしたNOMOTOモデルを考えます。

$\displaystyle{ \left\{\begin{array}{l} \dot{\psi}(t)=r(t) \\ \displaystyle{\dot{r}(t)=-\left(\frac{U}{U^*}\right)\frac{1}{T^*}r(t) +\left(\frac{U}{U^*}\right)^2\frac{K^*}{T^*}\delta(t)} \end{array}\right. }$

これに、舵のダイナミクス

$\displaystyle{ \dot{\delta}(t)=-\frac{1}{T_a}\delta(t)+\frac{K_a}{T_a}u(t) }$

を接続した状態空間表現を次のように表します。

$\displaystyle{ \underbrace{ \left[\begin{array}{c} \dot{\psi}(t) \\ \dot{r}(t) \\ \dot{\delta}(t) \end{array}\right] }_{\dot{x}(t)} = \underbrace{ \left[\begin{array}{ccc} 0 & 1 & 0\\ 0 & -\left(\frac{U}{U^*}\right)\frac{1}{T^*} & \left(\frac{U}{U^*}\right)^2\frac{K^*}{T^*} \\ 0 & 0 & -\frac{1}{T_a} \end{array}\right] }_{A(U,U^2)} \underbrace{ \left[\begin{array}{c} \psi(t) \\ r(t) \\ \delta(t) \end{array}\right] }_{x(t)} + \underbrace{ \left[\begin{array}{c} 0 \\ 0 \\ \frac{K_a}{T_a} \end{array}\right] }_{B} u(t) }$
$\displaystyle{ \underbrace{ \psi(t) }_{y(t)} = \underbrace{ \left[\begin{array}{ccc} 1 & 0 & 0 \end{array}\right] }_{C} \underbrace{ \left[\begin{array}{c} \psi(t) \\ r(t) \\ \delta(t) \end{array}\right] }_{x(t)} }$

ここで、次のようなポリトピック表現ができることに注目します（本式は未発表です）。

ship6a

$\displaystyle{ A(U,U^2)=p_1(U,U^2)A_1+p_2(U,U^2)A_2+p_3(U,U^2)A_3 }$

ただし

$\displaystyle{A_1=A(U_1,U_1^2)}$
$\displaystyle{A_2=A(U_2,U_2^2)}$
$\displaystyle{A_3=A(\frac{U_1+U_2}{2},U_1U_2)}$

また、 $U_3=\frac{U_1+U_2}{2}$ を定義して

$\displaystyle{ p_1(U,U^2)=\frac{1}{p_0}\det \left[\begin{array}{cc} U-U_3 & U_2-U_3 \\ U^2-U_1U_2 & U_2^2-U_1U_2 \\ \end{array}\right]}$
$\displaystyle{ p_2(U,U^2)=\frac{1}{p_0}\det \left[\begin{array}{cc} U_1-U_3 & U-U_3 \\ U_1^2-U_1U_2 & U^2-U_1U_2 \\ \end{array}\right]}$
$\displaystyle{ p_3(U,U^2)=\frac{1}{p_0}\det \left[\begin{array}{cc} U_1-U_2 & U_2-U \\ U_1^2-U_2^2 & U_2^2-U^2 \\ \end{array}\right]}$
$\displaystyle{ p_0=\det \left[\begin{array}{cc} U_1-U_2 & U_2-U_3 \\ U_1^2-U_2^2 & U_2^2-U_1U_2 \\ \end{array}\right]}$
$\displaystyle{p_1(U,U^2)+p_2(U,U^2)+p_3(U,U^2)=1}$

これに基づいて、次のようなLPV制御の仕組みを考えます。これは端点 $(U_1,U_1^2)$ 、 $(U_2,U_2^2)$ と仮想的な端点 $(U_3,U_1U_2)$ が作る３角形領域における変動をカバーすることができます。

ship6

次の２ポートシステムを構成します。

$\displaystyle{ P: \left\{\begin{array}{l} \left[\begin{array}{c} \dot{x} \\ \dot{x}_I \end{array}\right]= \underbrace{ \left[\begin{array}{cc} A(U,U^2)& 0 \\ -C & 0 \end{array}\right] }_{\cal{A}(U,U^2)} \left[\begin{array}{c} x \\ x_I \end{array}\right] + \underbrace{ \left[\begin{array}{c} 0 \\ 1 \end{array}\right] }_{B_1} r + \underbrace{\left[\begin{array}{c} B \\ 0 \end{array}\right] }_{B_2} u\\ \underbrace{ \left[\begin{array}{c} y_{11} \\ y_{12} \end{array}\right] }_{y_1} = \underbrace{ \left[\begin{array}{cc} 0 &\omega_I\\ \omega_DCA(U,U^2) & 0 \end{array}\right] }_{C_1} \left[\begin{array}{c} x \\ x_I \end{array}\right] + \underbrace{ \left[\begin{array}{c} 0 \\ 0 \end{array}\right] }_{D_{11}} r + \underbrace{ \left[\begin{array}{c} 0 \\ \omega_DCB \end{array}\right] }_{D_{12}} u\\ \underbrace{ \left[\begin{array}{c} y \\ x_I \end{array}\right] }_{y_2} = \underbrace{ \left[\begin{array}{cc} C & 0\\ 0 & 1 \end{array}\right] }_{C_2} \left[\begin{array}{c} x \\ x_I \end{array}\right] + \underbrace{ \left[\begin{array}{c} 0 \\ 0 \end{array}\right] }_{D_{21}} r \end{array}\right.}$

$A(U,U^2)=p_1(U,U^2)A_1+p_2(U,U^2)A_2+p_3(U,U^2)A_3$

これに対して、次のコントローラを設計します。

$\displaystyle{ K_0: \left\{\begin{array}{l} \dot{x}_K=A_K(U,U^2)x_K+ \underbrace{ \left[\begin{array}{cc} B_K^{(1)}(U,U^2) & B_K^{(2)}(U,U^2) \end{array}\right] }_{B_K(U,U^2)} \left[\begin{array}{c} y \\ x_I \end{array}\right] \\ u=C_K(U,U^2)x_K + \underbrace{ \left[\begin{array}{cc} D_K^{(1)}(U,U^2) & D_K^{(2)}(U,U^2) \end{array}\right] }_{D_K(U,U^2)} \left[\begin{array}{c} y \\ x_I \end{array}\right] \end{array}\right.}$

$A_K(U,U^2)=p_1(U,U^2)A_{K1}+p_2(U,U^2)A_{K2}+p_3(U,U^2)A_{K3}$
$B_K(U,U^2)=p_1(U,U^2)B_{K1}+p_2(U,U^2)B_{K2}+p_3(U,U^2)B_{K3}$
$C_K(U,U^2)=p_1(U,U^2)C_{K1}+p_2(U,U^2)C_{K2}+p_3(U,U^2)C_{K3}$
$D_K(U,U^2)=p_1(U,U^2)D_{K1}+p_2(U,U^2)D_{K2}+p_3(U,U^2)D_{K3}$

これから、次の出力フィードバックを得ます。

$\displaystyle{ K: \left\{\begin{array}{l} \left[\begin{array}{c} \dot{x}_K \\ \dot{x}_I \end{array}\right]= \underbrace{ \left[\begin{array}{cc} A_K(U,U^2) & B_K^{(2)}(U,U^2) \\ 0 & 0 \end{array}\right] }_{A_K(U,U^2)} \left[\begin{array}{c} x_K \\ x_I \end{array}\right] + \underbrace{ \left[\begin{array}{cc} B_K^{(1)}(U,U^2) & 0\\ -1& 1 \end{array}\right] }_{B_K(U,U^2)} \left[\begin{array}{c} y \\ r \end{array}\right] \\ u= \underbrace{ \left[\begin{array}{cc} C_K(U,U^2) & D_K^{(2)}(U,U^2) \end{array}\right] }_{C_K(U,U^2)} \left[\begin{array}{c} x_K \\ x_I \end{array}\right] + \underbrace{ \left[\begin{array}{cc} D_K^{(1)}(U,U^2) & 0 \end{array}\right] }_{D_K(U,U^2)} \left[\begin{array}{c} y \\ r \end{array}\right] \end{array}\right. }$

いま、 $U=U^*$ として設計した $H_\infty$ 制御による閉ループ系のシミュレーションを次に示します。３種類の速度変動に応じてばらつきか見られます。

ship8d
ship8a

これに対して、速度変動に応じたLPV制御による閉ループ系のシミュレーションを次に示します。 $H_\infty$ 制御の場合のばらつきが抑制されていることが分かります。