0 引　言

水面船舶在海上航行时，不可避免地遭遇随机风浪的干扰，由此将引起船舶六自由度运动姿态的随机运动。航向（首摇）的校正作用一般是靠自动舵控制系统承担，对于未加装减摇装置系统的船舶，其横摇运动是欠驱动系统，将会对船舶的稳定性、船上人员的舒适感、船载设备的运行性能和武备系统的命中率产生较大的影响。此外，横摇运动也将耦合影响航向的保持性能及航行的安全性^[1-3]。

船舶动力学、运动学研究结果表明^[4-5]，船舶航行时，舵叶上除了能生成首摇控制力（矩）外，还能够生成可观的横摇扶正力（矩），加之首摇运动动态特性处在较低频率域，而横摇运动相对首摇运动而言，处在较高频率域，这使得有可能寻求在不添装减摇装置系统情况下，通过合理的设计自动舵控制系统和策略，在保证航行方向稳定的同时，显著地减少船的横摇运动。

本文提出的矢量舵是在普通的后缘开襟形成一个独立的可动面—翼舵，从而与主舵形成2个相对独立的矢量控制。翼舵相当于一个可调整舵叶双侧曲率不对称度的控制面，通过翼舵控制面的转动控制，可改变舵叶双侧面曲率的不对称度，从而增大舵叶上的水动力（系数），即控制力（矩），实现提高舵效之目的。

目前，工程上应用的带有翼舵的航向控制系统，均为基于舵和翼舵是线性连动的控制面，没有实现矢量控制面，而实际上，在 $({u_r},{\alpha _r},{\beta _r})$ 三维空间中，每对应一个 ${u_r}$ 值将由多组甚至无穷多组 $({\alpha _r},{\beta _r})$ 的组合值与之对应，即对应二维 $({\alpha _r},{\beta _r})$ 平面上形成一个 ${u_r}({\alpha _r},{\beta _r})$ 空间曲面^[6-8]。因此，通过实现舵、翼舵相互独立运动的矢量控制，将能有效地提高控制效果，减小系统能耗。

本文提出了矢量舵减横摇控制系统技术，设计了 $\mu $ -鲁棒控制策略，利用改进的遗传算法构成舵角、翼舵角智能协调决策器，并进行了数字仿真试验，结果显示矢量舵减横摇控制系统能有效的减小船舶横摇运动，降低了系统的能耗。

1 船舶矢量舵减横摇控制系统建模

船舶矢量减横摇控制系统结构原理图如图1所示。

该系统主要有控制器（控制策略）、主舵角和翼舵角协调决策分配器，舵和翼舵组成的矢量舵、舵角、翼舵角矢量传动装置，舵、翼舵伺服驱动系统，横摇反馈信号传感器等构成。

船体横荡，首摇，横摇三自由度运动非线性模型为^[9]：

$ \left\{ {\begin{aligned} & \begin{array}{l} {\rm{m}}(\dot v + ur) = - {m_y}\dot y - {m_x}ur + {Y_v}v + {Y_r}r + {Y_{v\left| v \right|}}v\left| v \right| +\\ {Y_{r\left| r \right|}}r\left| r \right|+ {Y_{vvr}}{v^2}r + {Y_{vrr}}v{r^2} + {Y_r}({\alpha _r},{\beta _r}) + {Y_d}\text{，} \end{array}\\ & \begin{array}{l} {I_{\rm{z}}}\dot r = - \Delta {I_z}\dot r + {N_v}v + {N_r}r + {N_{v\left| v \right|}}v\left| v \right| + \\ {N_{r\left| r \right|}}r\left| r \right| + {N_{vvr}}{v^2}r + {N_{vrr}}v{r^2}+ {N_\varphi }\varphi + {N_{v\left| \varphi \right|}}v\left| \varphi \right| +\\ {N_{r\left| \varphi \right|}}r\left| \varphi \right| + {N_r}({\alpha _r},{\beta _r}) + {N_d}\text{，} \end{array}\\ & \begin{array}{l} {I_x}\dot p = - \Delta {I_x}\dot p - 2{N_p}p - Wh\varphi - {z_H}( - {m_y}\dot v -\\ {m_x}ur + {Y_v}v + {Y_r}r + {Y_{v\left| v \right|}}v\left| v \right|+ {Y_{r\left| r \right|}}r\left| r \right| + {Y_{vvr}}{v^2}r +\\ {Y_{vrr}}v{r^2}){\rm{ + }} {K_r}({\alpha _r},{\beta _r}) + {K_d}\text{。} \end{array} \end{aligned}} \right. $

(1)

式中： $m$ 为船舶质量； $u$ 为纵荡速度； ${I_z}$ 为船舶对Z轴的转动惯量； ${I_x}$ 为船舶对X轴的转动惯量； $v$ 为横荡速度； $r$ 为首摇角速度；P为横摇角速度； ${m_x}$ 为纵荡附加质量； ${m_y}$ 为横荡附加质量； $\Delta {I_{\rm{z}}}$ 为船舶对Z轴的附加转动惯量； $\Delta {I_x}$ 为船舶对X轴的附加转动惯量； $2{N_p}$ 为每单位横摇角速度的船舶阻尼力矩；W为船舶排水量；h横稳心高； ${z_H}$ 为船舶横向力的作用点至重心的垂向距离； $\varphi $ 为横摇角； $Y$ 、N为船舶粘性水动力系数； ${\alpha _{\rm{r}}}$ 为舵角； ${\beta _{{r}}}$ 为翼舵角； ${Y_r}({\alpha _r},{\beta _r})$ 、 ${N_r}({\alpha _r},{\beta _r})$ 、 ${K_r}({\alpha _r},{\beta _r})$ 为舵/翼舵产生的横荡力、首摇力矩和横摇力矩； ${Y_d}$ 、 ${N_d}$ 、 ${K_d}$ 为船体受到的横荡干扰力、首摇干扰力矩和横摇干扰力矩。

舵/翼舵产生的对船体的横荡力 ${Y_r}({\alpha _r},{\beta _r})$ 、首摇力矩 ${N_r}({\alpha _r},{\beta _r})$ 和横摇力矩 ${K_r}({\alpha _r},{\beta _r})$ 分布别为：

$\begin{split} & {Y_r}({\alpha _r},{\beta _r}) = (1 + {\alpha _H}){P_N}\cos {\alpha _r} \text{，}\\ & {N_r}({\alpha _r},{\beta _r}) = - (1 + {\alpha _H}){x_R}{P_N}\cos {\alpha _r} \text{，} \\ & {K_r}({\alpha _r},{\beta _r}) = (1 + {\alpha _H}){z_R}{P_N}\cos {\alpha _r} \text{。} \end{split} $

(2)

式中： ${P_N}$ 为舵/翼舵的法向力； ${x_R}$ ， ${z_R}$ 分别为舵/翼舵水动力作用点的纵向和垂向坐标。 ${\alpha _H}$ 为舵/翼舵与船体水动力影响系数，近似公式为^[10]：

${\alpha _H}{\rm{ = 1}}{\rm{.052}}{{\rm{C}}_b} + 0.125\text{。}$

(3)

考虑到舵/翼舵在水中的实际受力情况，有： ${P_N} \approx P$ ，其中 $P$ 为舵/翼舵在水流中受到的合力。

$P\cos {\alpha _r}{\rm{ = }}{P_y}$ , ${P_y}$ 为舵/翼舵的升力。近似计算公式为：

$\begin{split} & {Y_r}({\alpha _r},{\beta _r}) = (1 + {\alpha _H}){P_y} \text{，}\\ & {N_r}({\alpha _r},{\beta _r}) = - (1 + {\alpha _H}){x_R}{P_y} \text{，}\\ & {K_r}({\alpha _r},{\beta _r}) = (1 + {\alpha _H}){z_R}{P_y} \text{，}\\ & {P_y}{\rm{ = }}\frac{1}{2}\rho {\upsilon ^2}{S_p}{C_{yr}}\left( {{\alpha _r},{\beta _r}} \right) \text{。} \end{split} $

(4)

式中： $\rho $ 为海水密度； $\upsilon $ 为舵处来流速度； ${S_p}$ 为舵投影面积； ${C_{yr}}\left( {{\alpha _r},{\beta _r}} \right)$ 为舵/翼舵的升力系数。

由于受到船体伴流的影响，导致 $\upsilon $ 与正常水流的速度不一样，其中关系为： $\upsilon {\rm{ = }}\left( {1 - {\psi _c}} \right)V$ 。式中： $V$ 为船舶航速， ${\psi _c}{\rm{ = }}0.5{C_b} - 0.05$ 为船体的伴流系数。

矢量舵产生的控制力矩 ${K_r}({\alpha _r},{\beta _r})$ 选取 $m$ 阶回归模型，并且令 ${u_r}({\alpha _r},{\beta _r}) = {K_r}({\alpha _r},{\beta _r})$ ，则

$\begin{split} {u_r}({\alpha _r},{\beta _r}) = & {A_0} + {a_1}{\alpha _r} + {a_2}{\alpha _r}^2 + \cdots + {a_m}{\alpha _r}^m +\\ & {b_1}{\beta _r} + {b_2}{\beta _r}^2 +\cdots + {b_m}{\beta _r}^m +\\ & {c_1}{a_r}{\beta _r} + {c_2}{\alpha _r}^2{\beta _r} + \cdots + {c_{m - 1}}{\alpha _r}^{m - 1}{\beta _r} + \\ & {d_2}{a_r}{\beta _r}^2 + {d_3}{\alpha _r}{\beta _r}^3 + \cdots + {d_{m - 1}}{\alpha _r}{\beta _r}^{m - 1} \text{。} \end{split} $

(5)

由于舵为对称剖面，当 ${\alpha _r}{\rm{ = }}{\beta _r} = 0$ 时， ${u_r}({\alpha _r},{\beta _r}){\rm{ = }}0$ ，故 ${A_0}{\rm{ = }}0$ ，采用最小二乘法对回归系数 $( {{a_1},{a_2} \ldots ,{a_m};{b_1},}$ ${{b_2} \ldots ,{b_m};{c_1},{c_2} \ldots ,{c_{m - 1}};{d_2},{d_3} \ldots ,{d_{m - 1}}} )$ 进行参数估计，选取舵/翼舵尾弦比为0.25，剖面为NACA0021，展弦比为1.2的舵/翼舵升力系数实验曲线图谱^[11]采样数据，并进行m=1,2,3,4,5,6时回归参数的估计，图谱曲线如图2所示。根据拟合结果，当 $m \geqslant 4$ 时，其拟合误差平方 ${Q_{LS}} = {\displaystyle\sum\limits_{k = 1}^n {\left( {{u_r}(k) - {{\hat u}_r}(k)} \right)} ^2}$ 下降不明显，故选择回归模型阶数 $m{\rm{ = }}4$ ，此时，

$\begin{split} {c_r}({\alpha _r},{\beta _r}) = & 1.47{\alpha _r} + 3.04{\alpha _r}^2 - 3.82{\alpha _r}^3 - 1.46{\alpha _r}^4 + \\ & 0.29{\beta _r} + 3.86{\beta _r}^2 - 10.45{\beta _r}^3{\rm{ + }}8.37{\beta _r}^4 + \\ & {\rm{ 0}}{\rm{.50}}{a_r}{\beta _r} - 4.09{\alpha _r}^2{\beta _r} + 7.61{\alpha _r}^3{\beta _r} + \\ & 4.43{\alpha _r}{\beta _r}^2 - 4.59{\alpha _r}{\beta _r}^3 \text{。} \end{split} $

(6)

其中： ${c_r}({\alpha _r},{\beta _r})$ 为矢量舵水动力升力系数，且有：

${u_r}({\alpha _r},{\beta _r}) = \frac{1}{2}\rho {V^2}{S_r}{L_r}{c_r}\left( {{\alpha _r},{\beta _r}} \right)\text{。}$

(7)

式中： $\rho $ 为海水密度； $V$ 为舵处来流速度； ${S_r}$ 为矢量舵投影面积； ${L_r}$ 为舵力臂。

相同的方法，根据图3的扭矩系数试验曲线图谱，得到舵扭矩系数回归模型^[12]：

$\begin{split} {c_{mr}}({\alpha _r},{\beta _r}) =& 0.356\;6{\alpha _r} + 0.143\;4{\alpha _r}^2 + 1.210\;3{\alpha _r}^3 - \\ & 1.482\;8{\alpha _r}^4 + 0.381\;2{\beta _r} + 0.374\;1{\beta _r}^2 -\\ & 0.971\;2{\beta _r}^3{\rm{ + 0}}{\rm{.693\;9}}{\beta _r}^4 {\rm{ + 0}}{\rm{.111\;3}}{a_r}{\beta _r} + \\ &\! 0.720\;1{\alpha _r}^2{\beta _r}\! -\! 1.194\;8{\alpha _r}^3{\beta _r} \!+\! 2.257\;6{\alpha _r}{\beta _r}^2 -\!\!\!\!\!\!\!\!\!\! \\ & 3.103\;8{\alpha _r}{\beta _r}^3 - 0.623\;1{\alpha _r}^2{\beta _r}^2 \text{。} \end{split} $

(8)

根据图4的翼舵扭矩系数试验曲线图谱，得到翼舵扭矩系数回归模型：

$\begin{split} {c_{m3r}}({\alpha _r},{\beta _r}) =& 0.432\;4{\alpha _r} + 0.110\;0{\alpha _r}^2 - 0.511\;6{\alpha _r}^3+ \\ & 0.350\;6{\alpha _r}^4 + 0.347\;4{\beta _r} + 0.232\;1{\beta _r}^2 +\\ & 0.942\;7{\beta _r}^3 - 1.333\;5{\beta _r}^4 {\rm{ + 0}}{\rm{.241\;9}}{a_r}{\beta _r} +\\ & 1.892\;3{\alpha _r}^2{\beta _r} \!-\! 2.443\;2{\alpha _r}^3{\beta _r} \!+\! 0.424\;6{\alpha _r}{\beta _r}^2-\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \\ & 0.752\;8{\alpha _r}{\beta _r}^3 - 0.731\;0{\alpha _r}^2{\beta _r}^2 \text{。} \end{split} $

(9)

图5、图6和图7分别给出了舵/翼舵升力系数、舵扭矩系数和翼舵扭矩系数由图谱采样数据得到的曲面与拟合的回归模型计算得到的曲面之相对误差曲面（相对误差为图谱采样计算值与回归模型计算值的差值除以图谱采样计算值）。

2 $\mu $ -鲁棒控制策略设计

将式(1)中模型的非线性项和横荡影响并入到参数摄动和广义干扰项，忽略横荡的影响和方程中非线性项，得到简化的首摇/横摇线性耦合模型为^[13]：

$\left\{ \!\!\!\!\begin{array}{l} \left( {{I_x}{\rm{ + }}{a_{\varphi \varphi }}} \right)\ddot \varphi + {b_{\varphi \varphi }}\dot \varphi + {c_{\varphi \varphi }}\varphi + {a_{\varphi \psi }}\ddot \psi + {b_{\varphi \psi }}\dot \psi = {K_r}({\alpha _r},{\beta _r}) + {K_d}\text{，} \\ \left( {{I_z}{\rm{ + }}{a_{\psi \psi }}} \right)\ddot \psi + {b_{\psi \psi }}\dot \psi + {a_{\psi \varphi }}\ddot \varphi + {b_{\psi \varphi }}\dot \varphi = {N_r}({\alpha _r},{\beta _r}) + {N_d}\text{。} \\ \end{array} \right.$

(10)

记 ${x_1} = \varphi ,{x_2} = \dot \varphi ,{x_3} = \psi ,X = {\left[ {{x_1}{x_2}{x_3}} \right]^{\rm T}}$ ，则系统的状态方程为：

$\begin{split} & \dot X = \left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} 0&1&0 \end{array}} \\ { - E_1^{ - 1}{E_2}} \end{array}} \right]X + \left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} 0&0 \end{array}} \\ {E_1^{ - 1}} \end{array}} \right]\left[ \begin{array}{l} {K_r}({\alpha _r},{\beta _r}) \\ {N_r}({\alpha _r},{\beta _r}) \\ \end{array} \right] +\\ & \left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} 0&0 \end{array}} \\ {E_1^{ - 1}} \end{array}} \right]\left[ \begin{array}{l} {K_d} \\ {N_d} \\ \end{array} \right] y = CX,C = \left[ {\begin{array}{*{20}{c}} 1&0&0 \end{array}} \right] \text{。} \end{split} $

(11)

其中：

$ {E_1} = \left[\!\!\!\! {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {{I_x}{\rm{ + }}{a_{\varphi \varphi }}}&{{a_{\varphi \psi }}} \end{array}} \\ {\begin{array}{*{20}{c}} {{a_{\psi \varphi }}}&{{I_z}{\rm{ + }}{a_{\psi \psi }}} \end{array}} \end{array}} \!\!\!\!\right], {E_2} = \left[\!\!\!\! {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {{c_{\varphi \varphi }}}&{{b_{\varphi \varphi }}}&{{b_{\varphi \psi }}} \end{array}} \\ {\begin{array}{*{20}{c}} 0&{{b_{\psi \varphi }}}&{{b_{\psi \psi }}} \end{array}} \end{array}}\!\!\!\!\right]\text{。} $

(12)

将系统状态方程中由控制量舵/翼舵产生的状态 ${X_1}$ 和由干扰产生的状态 ${X_2}$ 分离，即记：

$\begin{aligned} & {{\dot X}_1} = \left[\!\!\!\! {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} 0&1&0 \end{array}} \\ { - E_1^{ - 1}{E_2}} \end{array}} \!\!\!\!\right]{X_1} + \left[\!\!\!\! {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} 0&0 \end{array}} \\ {E_1^{ - 1}} \end{array}} \!\!\!\!\right]\left[ \!\!\!\!\begin{array}{l} {K_r}({\alpha _r},{\beta _r}) \\ {N_r}({\alpha _r},{\beta _r}) \\ \end{array} \!\!\!\! \right] \text{，} \\ & {{\dot X}_2} = \left[\!\!\!\! {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} 0&1&0 \end{array}} \\ { - E_1^{ - 1}{E_2}} \end{array}} \!\!\!\!\right]{X_2} + \left[\!\!\!\! {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} 0&0 \end{array}} \\ {E_1^{ - 1}} \end{array}} \!\!\!\!\right]\left[ \!\!\!\!\begin{array}{l} {K_d} \\ {N_d} \\ \end{array} \!\!\!\!\right]\text{，} \\ & y = CX = C{X_1} + C{X_2} \text{。} \end{aligned} $

(13)

定义广义干扰 $\omega = C{X_2}$ ，由

$\left[ \!\!\!\!\begin{array}{l} {K_r}({\alpha _r},{\beta _r}) \\ {N_r}({\alpha _r},{\beta _r}) \\ \end{array} \!\!\!\! \right]{\rm{ = }}\left[\!\!\!\! \begin{array}{l} \left( {1{\rm{ + }}{a_H}} \right){z_R} \\ - \left( {1{\rm{ + }}{a_H}} \right){x_R} \\ \end{array} \!\!\!\! \right]{P_y}({\alpha _r},{\beta _r}) = {B_R}{P_y}({\alpha _r},{\beta _r})\text{，}$

(14)

得到：

$\left\{ \begin{split} & {{\dot X}_1} = \left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} 0&1&0 \end{array}} \\ { - E_1^{ - 1}{E_2}} \end{array}} \right]{X_1} + \left[ {\begin{array}{*{20}{c}} 0 \\ {E_1^{ - 1}{B_R}} \end{array}} \right]{P_y}({\alpha _r},{\beta _r})\text{，} \\ & y = C{X_1} + \omega\text{。} \\ \end{split} \right.$

(15)

这样就将系统的奇异控制问题变换为非奇异的标准控制问题。

记 ${q_1} = \left( {{c_{\varphi \varphi }},{b_{\varphi \varphi }},{b_{\varphi \psi }},{b_{\psi \varphi }},{b_{\psi \psi }}} \right)$ ， ${q_2} = \big( {{a_{\varphi \varphi }},{a_{\varphi \psi }},{a_{\psi \varphi }},}$ ${{a_{\psi \psi }}} \big),q = \left( {{q_1},{q_2}} \right)$ ，得

$M\left( q \right) = \left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {A\left( q \right)}&{B\left( q \right)} \end{array}} \\ {\begin{array}{*{20}{c}} {C\left( q \right)}&{D\left( q \right)} \end{array}} \end{array}} \right]\text{。}$

(16)

根据仿射参数摄动系统的线性分式表示问题的解法^[14]，求得：

$\begin{split} & M(g) = {F_u}({P_M},\Delta ) \text{，} \\ & \Delta = diag({c_{\varphi \varphi }},{b_{\varphi \varphi }},{b_{\varphi \psi }},{b_{\psi \varphi }},{b_{\psi \psi }},{a_{\varphi \varphi }},{a_{\varphi \psi }},{a_{\psi \varphi }},{a_{\psi \psi }}) \text{，} \\ \end{split} $

(17)

${P_M} = \left[ {\begin{array}{*{20}{c}} {{P_{M11}}}&{{P_{M12}}} \\ {{P_{M21}}}&{{P_{M22}}} \end{array}} \right]\text{，}$

(18)

假设对象参数 ${I_x}$ ， ${I_z}$ 的摄动范围为30%，得到：

$\Delta = {F_u}({P_\Delta },{\Delta _\delta })\text{，}$

(19)

${\Delta _\delta } = diag({\delta _1},{\delta _2},{\delta _3},{\delta _4},{\delta _5},{\delta _6},{\delta _7},{\delta _8},{\delta _9})\text{。} $

(20)

且 $\left| {{\delta _i}} \right| \leqslant 1$ ， ${P_\Delta } = \left[ {\begin{array}{*{20}{c}} {{P_{\Delta 11}}}&{{P_{\Delta 12}}} \\ {{P_{\Delta 21}}}&{{P_{\Delta 22}}} \end{array}} \right]$ ，利用“星积运算”，得到 $M(\delta )$ 的线性分式变换模型为：

$M(\delta ) = {F_u}(P,{\Delta _\delta })\text{。}$

(21)

其中：

$P = \left[ {\begin{split} {{F_l}({P_\Delta },{P_{M11}})}\;\;\;\;{{P_{\Delta 12}}{{(I - {P_{M11}}{P_{\Delta 22}})}^{ - 1}}{P_{M12}}} \\ {{P_{M21}}{{(I - {P_{\Delta 22}}{P_{M11}})}^{ - 1}}{P_{\Delta 12}}}\;\;\;\;\;{{F_u}({P_M},{P_{\Delta 21}})} \end{split}} \right]\text{。}$

(22)

图8为船舶舵/翼舵减摇鲁棒控制系统的原理结构图。 $\omega $ 为广义干扰， ${P_y}({\alpha _r},{\beta _r})$ 为舵/翼舵上产生的升力， ${W_p}$ 是反应系统控制目标的性能权函数，系统性能指标是对所有容许的不确定性 ${\Delta _\delta }$ 从d到Z的传递函数矩阵在 ${H_\infty }$ 范数意义下达到最小。

船减摇控制策略设计的目标为： $\left\| {\begin{array}{*{20}{c}} {{W_p}{S_r}} \\ {{W_u}{u_r}{S_r}} \end{array}} \right\| < 1$ ，其中， ${W_p}$ 为系统控制目标性能权函数， ${W_r}$ 为控制量 ${u_r}({\alpha _r},{\beta _r})$ 加权函数， $S$ 为灵敏度函数，且 ${S_r} = \displaystyle\frac{1}{{1 + G{K_{rk}}}}$ 。根据船体的横摇动态特性和首摇/横摇运动的分频特性，选择权函数为：

${W_p} = \frac{{0.39s}}{{{s^4} + 2.46{s^3} + 1.29{s^2} + 0.84s + 0.19}}\text{，} $

(23)

${W_r} = \frac{{0.03{s^2} + 0.03s + 0.01}}{{{s^2} + 100s + 0.001}}\text{。} $

(24)

利用Matlab中的D-K迭代算法求取船减摇控制器传递函数^[10]， $\mu $ 最大值为0.45，并经过降阶简化处理后，得舵减摇 $\mu $ —鲁棒控制器传递函数为：

$ {k_{rk}}\left( s \right) = \frac{{ - 8.89(s + 100)(s + 0.94)({s^2} + 0.08s + 0.003)}}{{(s + 5.10)(s + 0.63)(s + 0.05)({s^2} + 5.01s + 21.4)}}\text{。} $

(25)

3 舵角/翼舵角智能协调决策器设计

智能决策指标函数就是在矢量舵提供所需控制力矩 ${u_r}({\alpha _r},{\beta _r})$ 前提下，使系统能耗最小，即驱动能量最小^[15]。

若记 ${M_{\alpha r}}$ 为舵伺服系统驱动力矩， ${M_{\beta r}}$ 为翼舵伺服系统驱动力矩，则从k时刻到k+1时刻，系统驱动能量为：

$ \begin{split} {J_r}\left( {k + 1} \right) =& \int {_{{\alpha _r}\left( k \right)}^{{\alpha _r}\left( {k + 1} \right)}} {M_{\alpha r}}\left( {\theta ,{\beta _r}\left( k \right)} \right){\rm d}\theta + \\ & \int {_{{\alpha _r}\left( {k + 1} \right) + {\beta _r}\left( k \right)}^{{\alpha _r}\left( {k + 1} \right) + {\beta _r}\left( {k + 1} \right)}} {M_{\beta r}}\left( {{\alpha _r}\left( {k + 1} \right),\theta } \right){\rm d}\theta \text{，} \end{split} $

(26)

${M_{\alpha r}}{\rm{ + }}{M_{\alpha N}} + {M_{\alpha J}} + {M_{\alpha h}} + {M_{\alpha f}} = 0\text{，} $

(27)

${M_{\beta r}}{\rm{ + }}{M_{\beta N}} + {M_{\beta J}} + {M_{\beta h}} + {M_{\beta f}} = 0\text{。} $

(28)

其中： ${M_{\alpha N}}$ 为主舵扭矩； ${M_{\alpha J}}$ 为主舵惯性力矩； ${M_{\alpha h}}$ 为主舵恢复力矩； ${M_{\alpha f}}$ 为主舵摩擦力矩； ${M_{\beta N}}$ 为翼舵扭矩； ${M_{\beta J}}$ 为翼舵惯性力矩； ${M_{\beta h}}$ 为翼舵恢复力矩； ${M_{\alpha f}}$ 为翼舵摩擦力矩。于是有

$\begin{split} {J_r}\left( {k + 1} \right) =& \displaystyle\frac{1}{2}\rho {V^2}{S_{\alpha r}}{L_{\alpha r}}\int {_{{\alpha _r}\left( k \right)}^{{\alpha _r}\left( {k + 1} \right)}} {c_{\alpha r}}\left( {\theta ,{\beta _r}\left( {k{\rm{ + }}1} \right)} \right){\rm d}\theta + \\ & \left( {{M_{\alpha J}} + {M_{\alpha h}} + {M_{\alpha f}}} \right)\left( {{\alpha _r}\left( {k{\rm{ + }}1} \right) - {\alpha _r}\left( k \right)} \right) + \\ & \displaystyle\frac{1}{2}\rho {V^2}{S_{\beta r}}{L_{\beta r}}\int {_{{\alpha _r}\left( {k + 1} \right) + {\beta _r}\left( k \right)}^{{\alpha _r}\left( {k + 1} \right) + {\beta _r}\left( {k + 1} \right)}} {c_{\beta N}}\times \\ &\left( {{\alpha _r}\left( {k{\rm{ + }}1} \right),\theta } \right)d\theta + \left( {{M_{\beta J}} + {M_{\beta h}} + {M_{\beta f}}} \right) \times\\ &\left( {{\beta _r}\left( {k{\rm{ + }}1} \right) - {\beta _r}\left( k \right)} \right) \text{。} \end{split} $

(29)

其中： ${L_{\alpha r}}$ 和 ${L_{\beta r}}$ 分别为水动力作用点距主舵轴，翼舵轴的距离； ${c_{\alpha r}}\left( \cdot \right)$ 和 ${c_{\beta N}}\left( \cdot \right)$ 分别为主舵和翼舵的扭矩系数。

矢量舵上产生的控制力矩 ${u_r}({\alpha _r},{\beta _r})$ 是舵角和翼舵角的二元函数，当减摇调节器计算得到减横摇所需控制力矩为 ${K_r}(k + 1)$ 时，应满足：

${u_r}({\alpha _r}(k + 1),{\beta _r}(k + 1)) = {K_r}(k + 1)\text{，} $

(30)

考虑到舵角、舵角速度、翼舵角、翼舵角速度均有限制，故有：

$ \begin{split} & \left| {{\alpha _r}(k + 1)} \right| \leqslant {\alpha _{r\max }} \text{，} \\ & \left| {{{\dot \alpha }_r}(k + 1)} \right| \leqslant {{\dot \alpha }_{r\max }} \text{，} \\ & \left| {{\beta _r}(k + 1)} \right| \leqslant {\beta _{r\max }} \text{，} \\ & \left| {{{\dot \beta }_r}(k + 1)} \right| \leqslant {{\dot \beta }_{r\max }} \text{。} \end{split} $

(31)

舵角/翼舵角决策规则为在保证 ${u_r}({\alpha _r}(k + 1),{\beta _r}(k + 1))$ 同时，使系统驱动能量最小，即

$\begin{split} & {J_r}(k + 1)\min \text{，} \\ & {u_r}({\alpha _r}(k + 1),{\beta _r}(k + 1)) = {N_r}(k + 1) \text{，} \\ & \left| {{\alpha _r}(k + 1)} \right| \leqslant {\alpha _{r\max }} \text{，} \\ & \left| {{{\dot \alpha }_r}(k + 1)} \right| \leqslant {{\dot \alpha }_{r\max }} \text{，} \\ & \left| {{\beta _r}(k + 1)} \right| \leqslant {\beta _{r\max }} \text{，} \\ & \left| {{{\dot \beta }_r}(k + 1)} \right| \leqslant {{\dot \beta }_{r\max }} \text{。} \end{split} $

(32)

舵角/翼舵角最优智能决策追求的是^[16]在决策规则约束下快速精确地寻优给出 ${\alpha _{rg}}(k + 1)$ ， ${\beta _{rg}}(k + 1)$ 。本文采用改进遗传算法，通过初始种群和适应度函数选取、搜索空间范围确定、二进制编码、改进遗传算子和遗传操作等步骤，实现了舵角/翼舵角的智能协调决策分配器，其程序流程如图9所示。

4 系统仿真试验及结果

本文对某船矢量舵减横摇控制系统进行数字仿真，并给出了仿真试验结果。

仿真条件：船水动力参数见参考文献^[1]，海浪有义波高 ${H_s}$ 为3.15 m，4.2 m，浪向角 ${\psi _e}$ （迎浪时 ${\psi _e} = {0^ \circ }$ ）为30°，60°，90°，120°，150°，对象无参数摄动（ $\Delta {I_x} = 0$ ， $\Delta {I_z} = 0$ ）和参数摄动（ $\Delta {I_x} = {I_x} \times 30\% ,$ $\Delta {I_z} = {I_z} \times 30\% $ ），仿真统计结果见表1～表4。 $STD(\varphi )$ 为横摇角均方差， $\eta $ 为减摇率， $STD({\alpha _r})$ 为舵角均方根值， $STD({\beta _r})$ 为翼舵角均方根值。

仿真实验结果表明：1）船舶矢量舵减横摇控制系统能有效地减小横摇，在有利浪向下，减摇率达50%，这是非常可观的，说明了舵减摇的有效性；2）系统在横浪时减摇效果最显著；3）在尾斜浪时减摇效果最差，这是由于舵叶与水流相对速度和横摇首摇运动分频特性引起的，即舵只能减小较高频域的横摇，而对较低频域横摇无能为力；4） $\mu $ -鲁棒控制能有效抑制系统的参数摄动，而不使减摇性能受到明显影响；5）引入舵/翼舵矢量舵减摇控制系统与常规舵减摇控制系统相比，降低系统能耗约15%。

5 结　语

本文提出由舵和翼舵构成的矢量舵减横摇控制系统，设计了系统 $\mu $ -鲁棒控制器和基于改进遗传算法的舵角/翼舵角智能协调决策器，并对系统进行仿真，仿真结果表明矢量舵减横摇控制系统能有效减小横摇，降低系统能耗，且增强了抗系统参数摄动的鲁棒性。