0 引　言

隐蔽性是潜艇的重要特性之一，随着水下探测技术的发展，提高潜艇声隐身能力已成为各国潜艇技术发展的中心任务之一。螺旋桨噪声、机械噪声和水动力噪声是潜艇在水下航行时的主要噪声，各噪声源具有非线性、时变性和不确定性，以及相互影响和耦合作用^[1]。随着对潜艇绕流场和水动力噪声研究不断深入，潜艇水下运动时绕流场变化和操纵翼面对尾流场作用，进而对螺旋桨噪声和水动力噪声产生影响的机理已不断得到试验验证，使得潜艇低噪声航行工况成为操纵控制的一种必备模式^[2]。

Thomas F.Brooks等^[3]以NACA0012型翼面为研究对象，开展了0°～25.2°攻角下0.5×10⁶～4.6×10⁶雷诺数情况下翼面自噪声风动测试试验，提出了一套完整的平稳流中翼面自噪声半经验预报方法。Javier García^[4]和Diana Ovalle等^[5]设计了一种基于6自由度非线性耦合水动力方程的梯度下降寻优非线性深度和航向操纵控制算法，应用于操纵翼面优化、液压舵机功耗减少及降低操舵水动力噪声、减少复杂场景下潜艇过量运动等试验验证。试验结果揭示了潜艇运动时，操纵翼面涡流分离导致水动力噪声的发生机理。

国内相关研究人员借鉴国外操纵控制对流场扰动和水动力噪声的影响研究成果，从操纵控制对流场影响分析、低噪声操纵策略等方面开展了相关研究。李士强^[6]仿真分析了单独舵、直航打舵以及斜航工况下潜艇伴流对舵水动力性能的影响。赵桥生等^[7]开展了针对水下航行体机动航行时，流场扰动及水动力噪声的技术研究。王京齐等^[8]对低噪声安静操纵控制规律进行了探讨，提出了尽量减少对螺旋桨上流尾流场扰动等低噪声操纵控制的特殊要求。上述研究成果指明了“慢操舵、操小舵、缓变化”的低噪声操舵控制原则。

颜俐等^[9]为消除围壳舵产生的噪声，进而增加声呐探测距离，提出了低噪声航行工况下单尾舵模糊控制器，实现潜艇在垂直面深度和纵倾控制的操舵控制方案。黄利华等^[10]以降低舵机噪声、提高控制精度为目标，验证了潜艇低噪声航行工况下对转舵速率和舵角幅值约束的有效性。林超等^[11]在潜艇垂直面线性最优控制中，引入前置滤波器对指令信号进行柔化处理，提出了一种二自由度操纵控制策略，该操纵策略可降低操舵速度、减小操舵角和操舵次数。综合分析可知，上述控制策略和算法设计均需离线设计，仍不能较好适应多种复杂场景下潜艇低噪声操纵控制需求。

本文针对潜艇低噪声操纵控制需求，基于多约束条件模型预测控制理论，提出了一种根据操纵任务特点，实时设置深度变化率、纵倾角、舵角和舵速等约束条件，在线设计控制律参数的潜艇垂直面运动控制算法。为有效解决多种复杂场景下潜艇低噪声操纵控制问题提供一种思路方法。

1 潜艇垂直面运动的离散化差分状态方程

由潜艇操纵性^[12]理论可知，潜艇在水中空间运动，弱机动时可分解成2个平面运动，即潜艇水平面运动和垂直面运动。潜艇垂直面运动的非线性方程非常复杂，但低速运动潜艇，其垂直面运动可简化为如下线性状态方程和输出方程：

$ \begin{split} \underbrace {\left[ {\begin{array}{*{20}{c}} {\dot w} \\ {\dot q} \\ {\dot \theta } \\ {\dot \eta } \end{array}} \right]}_{\dot x} =& \underbrace {\left[ {\begin{array}{*{20}{c}} {{a_{11}}(V)}&{{a_{12}}(V)}&{{a_{13}}}&0 \\ {{a_{21}}(V)}&{{a_{22}}(V)}&{{a_{23}}}&0 \\ 0&1&0&0 \\ 1&0&{{a_{43}}(V)}&0 \end{array}} \right]}_{{A_c}}\underbrace {\left[ {\begin{array}{*{20}{c}} w \\ q \\ \theta \\ \eta \end{array}} \right]}_{{x_c}} +\\ & \underbrace {\left[ {\begin{array}{*{20}{c}} {{b_{11}}(V)}&{{b_{12}}(V)} \\ {{b_{21}}(V)}&{{b_{22}}(V)} \\ 0&0 \\ 0&0 \end{array}} \right]}_{{B_c}}\underbrace {\left[ {\begin{array}{*{20}{c}} {{\delta _b}} \\ {{\delta _s}} \end{array}} \right]}_u，\end{split} $

(1)

$ \underbrace {\left[ {\begin{array}{*{20}{c}} \theta \\ \eta \end{array}} \right]}_{\text{y}} = \underbrace {\left[ {\begin{array}{*{20}{c}} 0&0&1&0 \\ 0&0&0&1 \end{array}} \right]}_{{C_c}}\left[ \begin{gathered} \begin{array}{*{20}{c}} w \\ q \\ \theta \end{array} \\ \eta \\ \end{gathered} \right] 。$

(2)

式中： $ w $ 为垂速； $ q $ 为纵倾速度； $ \theta $ 为纵倾角； $ \eta $ 为下潜深度； $ {\delta _b} $ 为围壳舵（或首舵）舵角； $ {\delta _s} $ 为为升降舵舵角； $ V $ 为潜艇航速， $ {a_{ij}}(V) $ 和 $ {b_{ij}}(V) $ 是与航速相关的参数（以后为叙述方便，省略 $ (V) $ ）。 ${x_m} = {[w,q,\theta ,\eta ]^{\rm{T}}}$ ， $u = {[{\delta _b},{\delta _s}]^{\rm{T}}}$ ， $y = {\left[ {\theta ,\eta } \right]^{\rm{T}}}$ ，具体可参照文献[13]。

式（1）和式（2）简记为：

$ \begin{gathered} \dot x = {A_c}{x_c} + {B_c}u ，\\ y = {C_c}{x_c}。\\ \end{gathered} $

(3)

上式离散化后，可得：

$ \begin{split} & {x_m}(k + 1) = {A_m}{x_m}(k) + {B_m}u(k) ，\\ & y(k) = {C_m}{x_m}(k) 。\\ \end{split} $

(4)

式中： ${A_m} \approx I + {A_c} \cdot {T_s}$ ， ${B_m} \approx {B_c} \cdot {T_s}$ ， $ {C_m} = {C_c} $ ， $ I $ 为单位矩阵， $ {T_s} $ 为采样时间。 ${x_m}(k) = {[w(k),q(k),\theta (k),\eta (k)]^{\rm{T}}}$ ， $u(k) = {[{\delta _b}(k),{\delta _s}(k)]^{\rm{T}}}$ ， $y(k) = {\left[ {\theta (k),\eta (k)} \right]^{\rm{T}}}$ 。

将式(4)两端差分运算后，引入

$ \begin{split} &\Delta {x_m}\left( {k + 1} \right) = {x_m}\left( {k + 1} \right) - {x_m}\left( k \right) ，\\ & \Delta {x_m}\left( k \right) = {x_m}\left( k \right) - {x_m}\left( {k - 1} \right)，\\ &{\Delta u\left( k \right) = u\left( k \right) - u\left( {k - 1} \right)} 。\end{split} $

(5)

则可得如下增量式状态方程和状态输出方程：

$ \begin{split} & {\underbrace {\left[ {\begin{array}{*{20}{c}} {\Delta {x_m}\left( {k + 1} \right)} \\ {y\left( {k + 1} \right){\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} } \end{array}} \right]}_{{x_o}(k + 1)} = \underbrace {\left[ {\begin{array}{*{20}{c}} {{A_m}}&0 \\ {{C_m}{A_m}}&1 \end{array}} \right]}_{{A_o}}\underbrace {\left| {\begin{array}{*{20}{c}} {\Delta w(k)} \\ {\Delta q(k)} \\ {\Delta \theta (k)} \\ {\Delta \eta (k)} \\ {\theta (k)} \\ {\eta (k)} \end{array}} \right|}_{{x_o}(k)} + \underbrace {\left[ {\begin{array}{*{20}{c}} {{B_m}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} } \\ {{C_m}{B_m}} \end{array}} \right]}_{{B_o}}\Delta u\left( k \right)}，\\ & {\underbrace {y\left( k \right)}_{{y_o}(k)} = \underbrace {\left[ {\begin{array}{*{20}{c}} 0&0&1&0&0&0 \\ 0&0&0&1&0&0 \\ 0&0&0&0&1&0 \\ 0&0&0&0&0&1 \end{array}} \right]}_{{C_o}}\left| {\begin{array}{*{20}{c}} {\Delta w(k)} \\ {\Delta p(k)} \\ {\Delta \theta (k)} \\ {\Delta \eta (k)} \\ {\theta (k)} \\ {\eta (k)} \end{array}} \right|} 。\end{split} $

(6)

上式简记为：

$ \begin{split} &{{x_o}\left( {k + 1} \right) = {A_o}{x_o}\left( k \right) + {B_o}\Delta u\left( k \right)}，\\ &{{y_o}\left( k \right) = {C_o}{x_o}\left( k \right)} 。\end{split}$

(7)

式（7）可用于后续状态观测器的设计。

2 多约束条件下垂直面运动控制算法设计 2.1 深度控制器设计

潜艇一般用围壳舵控制深度，用尾升降舵控制纵倾^[14]。因此，在潜艇垂直面运动控制律设计中，多数采用分离设计思想，即围壳舵—深度控制律和尾升降舵—纵倾控制律设计方法。本文根据上述原则进行垂直面运动控制律设计。

从式（7）中抽取相应状态和参数，可获得下述深度控制律设计所需状态方程：

$ \begin{split} \underbrace {\left[ \begin{gathered} \Delta w(k + 1) \\ \Delta q(k + 1) \\ \Delta \theta (k + 1) \\ \Delta \xi (k + 1) \\ \eta (k + 1) \\ \end{gathered} \right]}_{x(k + 1)} =& \underbrace {\left[ {\begin{array}{*{20}{c}} {1 + {a_{11}}{T_s}}&{{a_{12}}{T_s}}&{{a_{13}}{T_s}}&0&0 \\ {{a_{21}}{T_s}}&{1 + {a_{22}}{T_s}}&{{a_{23}}{T_s}}&0&0 \\ 0&{{T_s}}&1&0&0 \\ {{T_s}}&0&{ - V{T_s}}&1&0 \\ {{T_s}}&0&{ - V{T_s}}&1&1 \end{array}} \right]}_A，\\ & \underbrace {\left[ \begin{gathered} \Delta w(k) \\ \Delta q(k) \\ \Delta \theta (k) \\ \Delta \eta (k) \\ \eta (k) \\ \end{gathered} \right]}_{x(k)} + \underbrace {\left[ \begin{gathered} {b_{11}}{T_s} \\ {b_{21}}{T_s} \\ 0 \\ 0 \\ 0 \\ \end{gathered} \right]}_B\underbrace {\Delta {\delta _b}(k)}_{\Delta {u_b}(k)}。\\[-65pt] \end{split} $

(8)

为了对深度变化率施加约束，则可设置输出方程为：

$ z（k）=\underset{C}{\underbrace{\left[\begin{array}{ccccc}0& 0& 0& 1& 0\end{array}\right]}}\left[\begin{array}{l}\Delta w(k)\\ \Delta q(k)\\ \Delta \theta (k)\\ \Delta \eta (k)\\ \eta (k)\end{array}\right]。$

(9)

式（8）和式（9）可简记为：

$ \begin{array}{l}x(k+1)=Ax(k)+B\Delta {u}_{b}(k)，\\ z(k)=Cx(k)。\end{array} $

(10)

考虑指令深度为 $ r({k_i}) $ ，构造指令向量 ${x_r}({k_i}) = {[0,0,0,0,r({k_i})]^{\rm{T}}}$ 。则可将式（10）中的状态变量变换为：

$ x(k) = {[\Delta w(k),\Delta q(k),\Delta \theta (k),\Delta \eta (k),\eta (k) - r({k_i})]^{\rm{T}}}。$

(11)

式（10）中的输出相应变换为深度偏差。后续控制律代价函数中将采用式（11）所示的状态向量。

假设在 $ k $ 时刻的 $ m $ 个采样间隔后解算的预测指令舵角增量为：

$ \Delta {u_b}\left( {k + m} \right) = {L_b}{\left( m \right)^{\rm{T}}}{\Theta _b}。$

(12)

式中： ${L_b}{\left( m \right)^{\rm{T}}}$ 为拟合深度MPC控制器的离散化Laguerre函数向量， $ {\Theta _b} $ 为Laguerre参数^[15-16]。

假设在 $ k $ 时刻，式(12)的状态为 $ x\left( k \right) $ ，则在此时刻后的 $ m $ 个采样间隔的预测状态为：

$ \begin{gathered} \begin{array}{*{20}{c}} {x\left( {k + {m \mathord{\left/ {\vphantom {m k}} \right. } k}} \right) = {A^m}x\left( k \right) + \sum\limits_{j = 0}^{m - 1} {{A^{m - j - 1}}} B{L_b}{{\left( j \right)}^{\rm{T}}}{\Theta _b}，{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} } \\ = {A^m}x\left( k \right) + \phi {{\left( m \right)}^{\rm{T}}}{\Theta _b} 。\end{array} \\ \phi {\left( m \right)^{\rm{T}}} = \sum\limits_{j = 0}^{m - 1} {{A^{m - j - 1}}} B{L_b}{\left( j \right)^{\rm{T}}} 。\\ \end{gathered} $

(13)

为设计深度预测测控制律，设代价函数为：

$ J(k) = \sum\limits_{m = 1}^{{N_p}} {(x{{({k_i} + m|{k_i})}^{\rm{T}}}Q} x({k_i} + m|{k_i}) + {\Theta ^{\rm{T}}}R\Theta 。$

(14)

式中：权矩阵 $ Q $ 为 $ 5 \times 5 $ 维半正定矩阵， $ R $ 为正定矩阵。

将式（13）所述预测状态代入上述代价函数，则利用Laguerre函数的正交特性，可得：

$ \begin{gathered} J(k) = {\Theta _b}^{\rm{T}}(\sum\limits_{m = 1}^{{N_p}} {\phi (m)} Q\phi {(m)^{\rm{T}}} + R){\Theta _b} + \\ 2{\Theta _b}^{\rm{T}}(\sum\limits_{m = 1}^{{N_p}} {\phi (m)} Q{A^m})x(k) + \sum\limits_{m = 1}^{{N_p}} {x{{(k)}^{\rm{T}}}} {({A^{\rm{T}}})^m}Q{A^m}x(k) = \\ {\Theta _b}^{\rm{T}}\Omega {\Theta _b}{\kern 1pt} + {\kern 1pt} {\kern 1pt} 2{\Theta _b}^{\rm{T}}\varPhi x(k) + \sum\limits_{m = 1}^{{N_p}} {x{{(k)}^{\rm{T}}}} {({A^{\rm{T}}})^m}Q{A^m}x(k)。\\ \end{gathered} $

(15)

式中：

$ \begin{gathered} \varOmega = \left(\sum\limits_{m = 1}^{{N_p}} {\phi (m)} Q\phi {(m)^{\rm{T}}} + R \right)，\\ \varPhi = \left(\sum\limits_{m = 1}^{{N_p}} {\phi (m)} Q{A^m}\right) 。\end{gathered} $

令上式对 $ {\Theta _b} $ 求导，可得无约束条件下，线性模型预测控制律为：

$ \Delta {u_b}(k + 1) = {{\boldsymbol{L}}_b}{\left( 0 \right)^{\text{T}}}{\varOmega ^{ - 1}}\varPhi x(k)。$

(16)

若在深度控制过程中，为降低操舵对绕流场影响，则代价函数求解需考虑以下约束：

1）舵角限制： $ {u_{\min }} \leqslant {\delta _s}(k) \leqslant {u_{^{\max }}} $

舵角约束不等式为：

$ {u_{\min }} \leqslant \sum\limits_{i = 0}^{m - 1} {{L_1}{{(i)}^{\rm{T}}}} {\Theta _b} + u({k_i} - 1) \leqslant {u_{\max }}。$

(17)

将 $ m = 1,2, \ldots ,{n_c} $ 对应的舵角约束整理成矩阵不等式，则可得：

$ \begin{split} & \underbrace {\left[ \begin{gathered} L{(0)^{\rm{T}}} \\ \sum\limits_{i = 0}^1 {L{{(i)}^{\rm{T}}}} \\ \vdots \\ \sum\limits_{i = 0}^{m - 1} {L{{(i)}^{\rm{T}}}} \\ \end{gathered} \right]}_{{M_u}}{\Theta _b} \leqslant \underbrace {\left[ \begin{gathered} {u_{\max }} - u({k_i} - 1) \\ {u_{\max }} - u({k_i} - 1) \\ \vdots \\ {u_{\max }} - u({k_i} - 1) \\ \end{gathered} \right]}_{{U_{\max }} - \bar u({k_i} - 1)} \text{，} \\ &- \underbrace {\left[ \begin{gathered} L{(0)^{\rm{T}}} \\ \sum\limits_{i = 0}^1 {L{{(i)}^{\rm{T}}}} \\ \vdots \\ \sum\limits_{i = 0}^{m - 1} {L{{(i)}^{\rm{T}}}} \\ \end{gathered} \right]}_{{M_u}}，\eta \leqslant \underbrace {\left[ \begin{gathered} {u_{\min }} + u({k_i} - 1) \\ {u_{\min }} + u({k_i} - 1) \\ \vdots \\ {u_{\min }} + u({k_i} - 1) \\ \end{gathered} \right]}_{ - {U_{\min }} + \bar u({k_i} - 1)} 。\end{split} $

简记为：

$ \begin{gathered} {M_u}{\Theta _b} \leqslant {U_{\max }} - \bar u({k_i} - 1) ，\\ - {M_u}{\Theta _b} \leqslant - {U_{\min }} + \bar u({k_i} - 1) 。\\ \end{gathered} $

(18)

2）舵速限制： $ \Delta {u_{\min }} \leqslant \Delta u \leqslant \Delta {u_{\max }} $

舵速约束不等式为：

$ \Delta {u_{\min }} \leqslant L{(m)^{\rm{T}}}{\Theta _b} \leqslant \Delta {u_{\max }}。$

(19)

将 $ m = 1,2, \ldots ,{n_c} $ 对应的舵速约束整理为矩阵不等式形式，则有：

$ \underbrace {\left[ \begin{gathered} L{(1)^{\rm{T}}} \\ L{(2)^{\rm{T}}} \\ \vdots \\ L{({n_c})^{\rm{T}}} \\ \end{gathered} \right]}_{{M_\Delta }}{\Theta _b} \leqslant \underbrace {\left[ \begin{gathered} \Delta {u_{\max }} \\ \Delta {u_{\max }} \\ \vdots \\ \Delta {u_{\max }} \\ \end{gathered} \right]}_{\Delta {U_{\max }}} \text{，} \underbrace {\left[ \begin{gathered} - L{(1)^{\rm{T}}} \\ - L{(2)^{\rm{T}}} \\ \vdots \\ - L{({n_c})^{\rm{T}}} \\ \end{gathered} \right]}_{ - {M_\Delta }}{\Theta _b} \leqslant \underbrace {\left[ \begin{gathered} - \Delta {u_{\min }} \\ - \Delta {u_{\min }} \\ \vdots \\ \Delta {u_{\min }} \\ \end{gathered} \right]}_{ - \Delta {U_{\min }}}。$

简记为：

$ \begin{gathered} {M_\Delta }{\Theta _b} \leqslant \Delta {U_{\max }}，\\ - {M_\Delta }{\Theta _b} \leqslant - \Delta {U_{\min }}。\\ \end{gathered} $

(20)

将式（18）和上式联立后简记为：

$ {M_H}{\Theta _b} \leqslant {\gamma _H} 。$

(21)

3）深度变化率限制： $ \Delta {\eta _{\min }} \leqslant \Delta \eta \leqslant \Delta {\eta _{\max }} $

在变深过程中，为防止深度变化率过大，可对深度变化率施加约束。变深过程中某些时间段内，可适当突破该约束，故该约束为软约束。

由式（10）和预测状态（13）可得：

$ \begin{gathered} x\left( {k + {m \mathord{\left/ {\vphantom {m k}} \right. } k}} \right) = {A^m}x\left( k \right) + \phi {\left( m \right)^{\rm{T}}}{\Theta _b} ，\\ z\left( {k + {m \mathord{\left/ {\vphantom {m k}} \right. } k}} \right) = C{A^m}x\left( k \right) + C\phi {\left( m \right)^{\rm{T}}}{\Theta _b}。\\ \end{gathered} $

(22)

由此可得：

$ \begin{gathered} \Delta {\eta _{\min }} \leqslant z({k_i} + m|{k_i}) \leqslant \Delta {\eta _{\max }} ，\\ - C\phi {\left( m \right)^{\rm{T}}}{\Theta _b} \leqslant - \Delta {\eta _{\min }} + C{A^m}x({k_i}) ，\\ C\phi {\left( m \right)^{\rm{T}}}{\Theta _b} \leqslant \Delta {\eta _{\max }} - C{A^m}x({k_i})。\\ \end{gathered} $

则有：

$ \begin{gathered} \underbrace { - \left[ \begin{gathered} C\phi {\left( 1 \right)^{\rm{T}}} \\ C\phi {\left( 2 \right)^{\rm{T}}} \\ \vdots \\ C\phi {\left( {Np} \right)^{\rm{T}}} \\ \end{gathered} \right]}_{ - {M_{\Delta \eta }}}{\Theta _b} \leqslant \underbrace { - {{\left[ \begin{gathered} 1 \\ 1 \\ \vdots \\ 1 \\ \end{gathered} \right]}_{Np}}\Delta {\eta _{\min }} + \left[ \begin{gathered} CA \\ C{A^2} \\ \vdots \\ C{A^{Np}} \\ \end{gathered} \right]x({k_i})}_{ - {\gamma _{\Delta \eta \_\min }}} ，\\ \underbrace {\left[ \begin{gathered} C\phi {\left( 1 \right)^{\rm{T}}} \\ C\phi {\left( 2 \right)^{\rm{T}}} \\ \vdots \\ C\phi {\left( {Np} \right)^{\rm{T}}} \\ \end{gathered} \right]}_{{M_{\Delta \eta }}}{\Theta _b} \leqslant \underbrace {{{\left[ \begin{gathered} 1 \\ 1 \\ \vdots \\ 1 \\ \end{gathered} \right]}_{Np}}\Delta {\eta _{\max }} - \left[ \begin{gathered} CA \\ C{A^2} \\ \vdots \\ C{A^{Np}} \\ \end{gathered} \right]x({k_i})}_{{\gamma _{\Delta \eta \_\max }}} 。\\ \end{gathered} $

联立为矩阵不等式，简记为：

$ {M_S}{\Theta _b} \leqslant {\gamma _S} 。$

(23)

联立式（21）和式（23），有

$ \underbrace {\left[ {\begin{array}{*{20}{c}} {{M_H}} \\ {{M_S}} \end{array}} \right]}_M{\Theta _b} \leqslant \underbrace {\left[ {\begin{array}{*{20}{c}} {{\gamma _H}} \\ {{\gamma _S}} \end{array}} \right]}_\gamma 。$

(24)

至此，多约束条件下，深度MPC控制器的求解问题被转换为二次规划问题：

$ \left\{ \begin{gathered} {J_1}(k) = \frac{1}{2}{\Theta _b}^{\text{T}}{\boldsymbol{\varOmega }}{\Theta _b}{\kern 1pt} + {\kern 1pt} {\kern 1pt} {\Theta _b}^{\text{T}}{\boldsymbol{\varPhi x}}(k)，\\ {\boldsymbol{M}}{\varTheta _b} \leqslant {\boldsymbol{\gamma }}。\\ \end{gathered} \right. $

(25)

该二次规划问题采用活动集法求解，可利用Hildreth二次规划方法求解^[15-17]。故深度MPC控制器可设计为：

$ \Delta {u_b}(k + 1) = {{\boldsymbol{L}}_b}{\left( 0 \right)^{\text{T}}}{\Theta _b} 。$

(26)

2.2 纵倾控制器设计

从式（7）中抽取与纵倾相关的和参数，可获得下述纵倾控制律设计所需状态方程：

$ \begin{split} \left[ \begin{gathered} \Delta w(k + 1) \\ \Delta q(k + 1) \\ \Delta \theta (k + 1) \\ \theta (k + 1) \\ \end{gathered} \right] =& \left[ {\begin{array}{*{20}{c}} {1 + {a_{11}}T}&{{a_{12}}{T_s}}&{{a_{13}}{T_s}}&0 \\ {{a_{21}}{T_s}}&{1 + {a_{22}}{T_s}}&{{a_{23}}T}&0 \\ 0&{{T_s}}&1&0 \\ 0&{{T_s}}&1&1 \end{array}} \right]，\\ & \left[ \begin{gathered} \Delta w(k) \\ \Delta q(k) \\ \Delta \theta (k) \\ \theta (k) \\ \end{gathered} \right] + \left[ \begin{gathered} {b_{12}}{T_s} \\ {b_{22}}{T_s} \\ 0 \\ 0 \\ \end{gathered} \right]\Delta {\delta _s}(k) 。\end{split} $

(27)

考虑到变深过程中纵倾角对潜艇流场，尾升降舵对螺旋桨尾部紊流场有较大影响，故在纵倾控制律设计中，需将纵倾和尾升降舵作为约束进行相应的理。

为导出纵倾约束条件，设置如下输出方程：

$ {y_\theta }(k) = \left[ {0,0,0,1} \right]\left[ {\begin{array}{*{20}{c}} {\Delta w(k)} \\ {\Delta q(k)} \\ {\Delta \theta (k)} \\ {\theta (k)} \end{array}} \right]。$

(28)

与深度控制器设计方法类似，将纵倾、尾升降舵舵速（舵增量）、尾升降舵角作为约束条件，可设计尾升降舵——纵倾控制器如下：

$ \Delta {u_s}(k + 1) = {{\boldsymbol{L}}_\theta }{\left( 0 \right)^{\text{T}}}{\Theta _\theta }。$

(29)

3 垂直面运动状态观测器设计

由式（8）和式（22）可知，在深度控制器和纵倾控制器中，需要用到 $ \Delta w(k) $ 和 $ \Delta q(k) $ 信息，而在潜艇操纵控制系统中，没有垂速 $ w $ 和纵倾变化率 $ q $ 的量测装置，故不能用二者量测信号直接求取 $ \Delta w(k) $ 和 $ \Delta q(k) $ 。因此，需设计相应状态观测器来提取 $ \Delta w(k) $ 和 $ \Delta q(k) $ 的观测信号。

由式（7）和式（8）可知， $ \Delta w(k) $ 和 $ \Delta q(k) $ 具有较强耦合，且围壳舵和尾升降舵对二者均有影响。因此，需采用耦合设计思想设计相应的状态观测器。

与最优控制理论类似，在模型预测控制中，观测器和控制器可分别设计。采用基于对偶原理观测器设计^[18]方法，构造式（6）和式（7）描述系统的对偶系统方程如下：

$ \begin{gathered} z(k + 1) = {A_o}^ * z(k) + {C_o}^ * \upsilon (k) ，\\ \vartheta (k) = {B_o}^ * z(k)。\\ \end{gathered} $

(30)

针对该对偶系统，运用类似前述控制器设计方法，构造多维Laguerre函数，求解观测器增益。设计如下全维状态观测器：

$ \begin{gathered} {{\hat x}_o}\left( {k + 1} \right) = {A_o}{{\hat x}_o}\left( k \right) + {B_o}\Delta u\left( k \right) + \\ {K_o}({y_o}(k) - {C_o}{{\hat x}_o}\left( k \right))。\\ \end{gathered} $

(31)

式中， $ {K_o} $ 为对偶系统控制器增益的共轭转置。

$ \begin{split} & {K_o}^ * = {L_o}{(0)^{\rm{T}}}\left(\sum\limits_{m = 1}^{{N_{po}}} {{\phi _o}(m)} {Q_o}{\phi _o}{(m)^{\rm{T}}} + {R_o}\right)\times\\ & \left(\sum\limits_{m = 1}^{{N_{po}}} {{\phi _o}(m)} {Q_o}{A_o}^ * \right), \\ & {\phi _o}{\left( m \right)^{\rm{T}}} = \sum\limits_{j = 0}^{m - 1} {A{{_o^ * }^{m - j - 1}}} C_o^ * {L_o}{\left( j \right)^{\rm{T}}} 。\end{split} $

(32)

4 仿真结果与分析

为验证本文提出控制算法的有效性，利用文献[4]潜艇6自由度方程作为潜艇运动数学模型进行有约束和无约束条件下的仿真验证实验。该潜艇的部分参数如表1所示。

控制器和观测器的Laguerre矩阵维数均设置为7，控制器的Laguerre网络极点均设置为 $ \alpha = 0.55 $ ；观测器的Laguerre网络极点设置为 $ \alpha = [0.55,0.55,0.55,0.55] $ ，控制器预测步设置为 $ {N_p} = 80 $ ，预测控制步长为 $ {N_c} = 0.2{N_p} $ 。观测器预测步长与预测控制步长一致，控制器采样时间为1s。

离散时间MPC控制器为在线实时设计，故约束条件均可由操纵人员依据战术任务而实时设置。

1）6 kn航速下潜60 m操纵控制仿真实验

仿真实验中，设置初始深度为60 m，指令深度为120 m，围壳舵舵角约束为 $ {{ - 2}}{{\text{0}}^ \circ } \leqslant {\delta _b} \leqslant {\text{2}}{{\text{0}}^ \circ } $ ，尾升降舵舵角约束为 $ {{ - }}{{\text{5}}^ \circ } \leqslant {\delta _s} \leqslant {{\text{5}}^ \circ } $ ，两舵舵角增量（即最大舵速）约束均设置为 $ {{ - }}{{\text{2}}^ \circ } \leqslant \Delta \delta \leqslant {{\text{2}}^ \circ } $ 。在变深过程中，设置纵倾约束为 $ {{ - }}{{\text{2}}^ \circ } \leqslant \theta \leqslant {{\text{2}}^ \circ } $ ，深度增量约束设置为 $ - 0.5\;{\rm{m}} \leqslant \Delta \eta \leqslant 0.5\;{\rm{m}} $ （即深度变化率不超过0.5 $ {\rm{m/s}} $ ），仿真时间设置为600 s。操纵控制仿真实验结果如图1所示。

2）12 kn航速上浮60 m仿真实验

仿真实验中，设置初始深度为120 m，指令深度为60 m，围壳舵舵角约束为 $ {{ - 2}}{{\text{0}}^ \circ } \leqslant {\delta _b} \leqslant {\text{2}}{{\text{0}}^ \circ } $ ，尾升降舵舵角约束为 $ {{ - }}{{\text{5}}^ \circ } \leqslant {\delta _s} \leqslant {{\text{5}}^ \circ } $ ，两舵舵角增量（即最大舵速）约束均设置为 $ {{ - }}{{\text{2}}^ \circ } \leqslant \Delta \delta \leqslant {{\text{2}}^ \circ } $ ，变深过程中，设置纵倾约束为 $ {{ - }}{{\text{2}}^ \circ } \leqslant \theta \leqslant {{\text{2}}^ \circ } $ 。深度增量约束为 $ - 1\;{\rm{m}} \leqslant \Delta \eta \leqslant 1\;{\rm{m}} $ （即深度变化率不超过1.0 m/s），仿真时间设置为600 s。操纵控制仿真实验结果如图2所示。

由图1（a）和图2（a）可知，潜艇在6 kn和12 kn航速下的变深控制过程中运动平稳、无超调、无振荡。结合相应纵倾变化曲线和深度变化率曲线可知，在变深操纵控制过程中，潜艇纵倾和深度变化率均在约束范围内。由此表明：潜艇变深过程中，纵倾和深度增量均没有超越各自约束值，深度控制和纵倾控制算法能够高效合理解算出围壳舵角和尾升降舵角指令。对比图1（a）和图2（a）中无约束变深曲线和多约束条件变深曲线可知，在变深动态过程中，多约束条件下过渡过程时间较无约束条件有所增加。结合纵倾变化曲线可知，过渡过程时间的增加与纵倾动态变化密切相关。如图1(b)中，无约束条件下，最大纵倾值为−3.2°，而多约束条件下最大纵倾值仅为−2.0°；由图1(c)和图2(c)可知，无约束条件下，深度变化率均大于多约束条件，表明深度变化率约束条件对深度控制律的求解有影响，纵倾约束条件对纵倾控制律的求解有影响。而由图1（d）和图2（d）可知，在变深过程中，围壳舵角和尾升降舵舵角均在各自约束范围内，表明所设计的多约束控制器合理有效。结合图1(b)和图1（d）也可看到，在62 s～64 s时间段内，在有约束条件下，纵倾在62 s时刻达到了约束下限。此时，为了将纵倾角保持在约束范围内，尾升降舵在随后2个采样步长内，从3.6°正舵角变化到了−3.0°负舵角，最大舵速为4.5°，超出最大舵速约束值。如文献[17]所述，其原因是采用的Hildreth二次规划算法属于逐个元素搜索方法，当有效集个数大于决策变量中元素个数或有效约束线性相关，则迭代会在计数其最大值时终止，但算法不会出错，只是解算结果会得到一个违背约束条件的近优解。由此也可表明，本文的Hildreth二次规划算法合理。从图1和图2还可看出，本文提出的在线设计状态观测器可实时、准确跟踪实际信号，表明基于对偶原理和模型预测控制理论，实时设计的状态观测器具有良好的收敛性和稳定性。

5 结　语

本文针对潜艇低噪声操纵控制问题，提出了一种在深度变化率、纵倾、舵角和舵速等多种约束条件下，垂直面运动模型预测控制算法。探索出一种多约束条件下深度控制器、纵倾控制器和垂直面运动状态观测器的在线设计方法，实现了潜艇垂直面低噪稳定运动控制。仿真结果表明：1）基于Laguerre函数的模型控制设计方法可在线设计相应深度控制器、纵倾控制器，控制效果良好，易于实现；2）基于对偶原理，结合Laguerre函数的模型控制设计方法可在线设计垂直面运动状态观测器，所设计的状态观测器能够实时精确地跟踪深度与纵倾，并可对深度增量和纵倾增量进行实时精确提取，具有较好的工程实用价值。