0 引　言

随着人类对海洋探索的推进，深海定位技术变得越发重要，其中动力定位系统（Dynamic Positioning System，DP）就是深海定位技术的重要组成部分，且因其能在复杂的海况及超深水领域发挥作用，重要性日益提升。动力定位系统是指通过自身推进器之间的相互配合，来抵消外界风浪流干扰的系统，主要由位置测量系统，控制系统和推力系统等子系统组成^[1]，其中，控制系统和推力系统是整个动力定位系统的核心。Fossen^[2]对动力定位的众多算法进行了详细介绍，并给出了众多算例。近年来随着智能化发展，众多学者致力于先进控制算法的研究，如非线性反步法^[3]、动态面控制方法^[4]、模糊自适应控制^[5]、神经网络与鲁棒控制相结合的方法^[6]、基于深度强化学习的方法^[7]等，都展现了各自的控制特点。动力定位船舶在控制过程中有众多物理限制，如推进器角度的变化、推力的变化，这些变化都是有界的，且非瞬时的。模型预测控制（Model Predictive Control，MPC）因其在提前处理约束方面的优越性能，以及其在线滚动优化的特点^[8]，越来越多的学者研究将其应用到动力定位船舶上。

最初康斯伯格公司提出将MPC技术用于动力定位系统，主要因为其对目标函数的设计有利于降低能源的消耗，其研究结果表明该策略可使每年的燃料使用量降低近20%^[9]。随后又因MPC可处理约束系统及在线滚动优化等优点，与动力定位技术所面临的问题相适应，MPC开始成为动力定位领域研究的一大热点。Fannemel^[10]为解决线性化模型难以准确描述动力定位船舶，提出将非线性模型预测控制（The Nonlinear Model Predictive Control，NMPC）应用到动力定位上，并对其应用于动力定位系统的实际价值进行评估。Sotnikova和Veremey^[11]则对变化缓慢的非线性模型进行了研究，在可预测时域内进行线性化，并假设线性模型和非线性模型的误差是定值，以此将复杂的动力定位船数学模型进行了线性化处理。

但目前大部分学者专注于对非线性模型预测的研究，却很少有学者关注模型预测控制的参数自整定问题，大部分还依赖于专家经验和以往数据里进行调参。而MPC的参数整定，是MPC控制策略至关重要的问题，通常在不同工况下，对参数的要求不同，在环境干扰多变的工况下，MPC的参数能够随船舶状态自动调整，对提高定位性能具有重要意义。本文提出一种策略，使用模糊算法对模型预测控制中的参数进行实时参数调整，在DP船状态发生改变时，模型预测参数也会发生相应的变化，以此来提高模型预测控制的精度。在此基础上设计一款非线性干扰观测器（Nonlinear disturbance observer，NDO）用于估计未知的环境力，实现对环境力的补偿，进一步提升系统的鲁棒性。

1 DP船数学模型 1.1 动力学模型和运动学模型

通常对于水面动力定位船舶只需考虑水平面三自由度运动，即横荡，纵荡和首摇。在大地坐标系中船舶位置和欧拉角描述为 $ \eta = {\left[ {x,{\text{ }}y,{\text{ }}\psi } \right]^{\rm{T}}} $ ，随体坐标系中船体的运动描述为 $ v = {\left[ {u,{\text{ }}v,{\text{ }}r} \right]^{\rm{T}}} $ 。运动学模型和动力学模型可描述为：

$ \dot \eta = {\boldsymbol{R}}(\psi )\nu ，$

(1)

$ {\boldsymbol{M}}\dot \nu + {\boldsymbol{C}}(\nu )\nu + {\boldsymbol{D}}(\nu )v = {\tau _{env}} + {\tau _c}。$

(2)

式中： $ {\boldsymbol{M}} $ 为质量矩阵，包含有附加质量； $ {\boldsymbol{C}}(\nu ) $ 为科氏力和向心力作用矩阵； $ {\boldsymbol{D}}(\nu ) $ 为船舶阻尼系数矩阵，当船速较低时，可简化成常量 $ {\boldsymbol{D}} $ ， $ {\tau _{env}} $ 和 $ {\tau _c} $ 分别为外干扰和船舶控制力。转换矩阵 ${\boldsymbol{ R}}(\psi ) $ 表示如下：

$ {\boldsymbol{R}}(\psi ) = \left[ {\begin{array}{*{20}{c}} {\cos \psi }&{ - \sin \psi }&0 \\ {\sin \psi }&{\cos \psi }&0 \\ 0&0&1 \end{array}} \right]。$

(3)

1.2 非线性干扰观测器设计

动力定位船舶会受到风浪流等多种环境力干扰，考虑这些环境干扰将导致动力定位产生较大误差，而面对时变未知环境干扰时，控制器很难实时获得精确环境干扰力。目前常用的解决方法是为控制系统设计干扰观测器，本文设计一款非线性干扰观测器，并对其稳定性进行证明。观测器设计如下：

$ \left\{ \begin{gathered} {{\hat \tau }_{env}} = Z + K{\boldsymbol{M}}v ，\\ \dot Z = - KZ - K( - {\boldsymbol{C}}(\nu )\nu - {\boldsymbol{D}}(\nu )v + {\tau _c} + K{\boldsymbol{M}}v。\end{gathered} \right. $

(4)

其中，Z为辅助函数， $ {\hat \tau _{env}} $ 为对干扰估计量，K为正定对角矩阵。对式（2）进行变换：

$ {\boldsymbol{M}}\dot v = - {\boldsymbol{C}}(\nu )\nu - {\boldsymbol{D}}(\nu )v + {\tau _c} + {\tau _{env}}，$

(5)

于是有：

$ \begin{split} {{\dot {\hat \tau} }_{env}} =& - KZ - K( - {\boldsymbol{C}}(v)v - {\boldsymbol{D}}(v )v + {\tau _c} + KMv) + K{\boldsymbol{M}}\dot v = \\ &- K(Z + K{\boldsymbol{M}}v) + K{\tau _{env}} = K({\tau _{env}} - {{\hat \tau }_{env}}) = K{{\tilde \tau }_{env}} ，\end{split} $

(6)

因此：

$ {\dot {\hat {\tau} }_{env}} = {\dot \tau_{env}} - {{\hat \tau }_{env}} = {\dot \tau _{env}} - K{{\tilde \tau }_{env}}。$

(7)

设计李雅普诺夫函数如下：

$ {V_0} = \frac{1}{2}{\tilde \tau ^{\rm{T}}}_{env}{\tilde \tau _{env}} ，$

(8)

则

$ \begin{split} {{\dot V}_0} =& {{\tilde \tau }^{\rm{T}}}_{env}{{\dot {\tilde \tau }}_{env}} = {{\tilde \tau }^{\rm{T}}}_{env}{{\dot \tau }_{env}} - {{\tilde \tau }^{\rm{T}}}_{env}K{{\tilde \tau }_{env}} \leqslant\\ &- {\lambda _{\min }}(K){{\tilde \tau }^{\rm{T}}}_{env}{{\tilde \tau }_{env}} + \varepsilon {{\tilde \tau }^{\rm{T}}}_{env}{{\tilde \tau }_{env}} +\\ &\frac{1}{{4\varepsilon }}{{\dot \tau }^{\rm{T}}}_{env}{{\dot \tau }_{env}} \leqslant - {k_0}{V_0} + \frac{1}{{4\varepsilon }}k_1^2 。\end{split} $

(9)

其中， $ {k_0} = 2({\lambda _{\min }}(K) - \varepsilon ) $ ， $ {\lambda _{\min }}(K) $ 为正定对角矩阵K的最小特征值， $ \varepsilon $ 为任意可设计小量正值，满足 $ {\lambda _{\min }}(K) > \varepsilon $ 。假设波浪变化有界，即对任意 $ |{\dot {\tilde \tau }_{env}}| \leqslant {k_1} $ 。

解不等式可得：

$ 0 \leqslant {V_0} \leqslant \frac{{k_1^2}}{{4\varepsilon {k_0}}} + \left({V_0}(0) - \frac{{k_1^2}}{{4\varepsilon {k_0}}}\right){e^{ - {k_0}t}}。$

(10)

即有：

$ ||{\tilde \tau _{env}}|| \leqslant \sqrt {\frac{{k_1^2}}{{2\varepsilon {k_0}}} + 2*\left({V_0}(0) - \frac{{k_1^2}}{{4\varepsilon {k_0}}}\right){e^{ - {k_0}t}}}。$

(11)

可知， $ {\tilde \tau _{env}} $ 最终收敛到 $ R = \sqrt {\dfrac{{k_1^2}}{{2\varepsilon {k_0}}}} $ 的零点邻域中，可通过选取适合的 $ {k_0} $ 来加快收敛速度，以及减小收敛误差。

1.3 非线性模型预测模型

为将模糊算法应用于模型预测控制，选择如下非线性集中控制模型：

$ \begin{split} & \dot X(t) = f(X(t),U(t),\\ &{\tau _{env}}(t)) = A(t)X(t) + B(t)U(t) + G(t){\tau _{env}}(t)，\\ & Y(t) = C(t)X(t) 。\end{split} $

(12)

其中： $ X(t) = {[\eta {(t)^{\rm{T}}},\nu {(t)^{\rm{T}}}]^{\rm{T}}} $ 为状态向量；A(t)、B(t)、G(t)、C(t)为状态空间方程的时系数。

$ A(t) = \left[ \begin{array}{cc} {0_{3 \times 3}}& R(\psi (t)) \\ {0_{3 \times 3}}& {M^{ - 1}}( - C(\nu (t)) + D) \end{array} \right]，$

(13)

$ B(t) = \left[ \begin{gathered} {0_{3 \times 3}} \\ {M^{ - 1}} \\ \end{gathered} \right]T(\alpha (t)) ，$

(14)

$ G(t) = \left[ \begin{gathered} {0_{3 \times 3}} \\ {M^{ - 1}} \\ \end{gathered} \right]，$

(15)

$ C = \left[ {{I_{3 \times 3}}\quad{{\text{0}}_{3 \times 3}}} \right] 。$

(16)

其中， $ T(\alpha (t)) = [{T_1}({\alpha _1}(t)),{T_2}({\alpha _2}(t)),......,{T_n}({\alpha _n}(t))] $ ，为不同推进器之间状态（推进器角度和推力）与推进器控制合力之间的关系。 $ {T_i}({\alpha _i}(t)) \in {R^{3 \times 1}} $ ， $ {T_i}({\alpha _i}(t)) $ 可表示如下：

$ {T_i}({\alpha _i}(t)) = \left[ \begin{gathered} {\text{ }}\cos {\alpha _i}(t) \\ {\text{ }}\sin {\alpha _i}(t) \\ l{x_i}\sin {\alpha _i}(t) - l{y_i}\cos {\alpha _i}(t) \\ \end{gathered} \right] 。$

(17)

基于采样时间 ${T_s}$ ，对式(6)离散化后有如下形式：

$ \begin{split} & X(k + 1) = X(k) + \int\nolimits_{k{T_s}}^{(k + 1){T_s}} {f(X(t),U(t),{\tau _{env}}(t)){\rm{d}}t} ，\\ & Y(k) = CX(k) 。\end{split} $

(18)

其中， $ k $ 表示当前时间的步数。此时可设计目标函数如下：

$ \begin{split} \min J(U) = & \sum\limits_{j = 1}^{{N_P} - 1} {(||Y(k + j|k) - r(k + j)||_{{Q_{err}}}^2) + } \\ &||Y(k + {N_P}|k) - r(k + {N_P})||_{{q_{ter}}{Q_{err}}}^2 \times\\ &\sum\limits_{j = 0}^{{N_C} - 1} {||U(k + j|k)||_{{Q_{in}}}^2 + } \sum\limits_{j = 0}^{{N_C} - 1} {||\Delta U(k + j|k)||_{{Q_{din}}}^2} 。\end{split} $

(19)

$ \begin{split} & {U_{\min }} \leqslant U(k + j|k) \leqslant {U_{\max }} ，\\ &\Delta {U_{\min }} \leqslant \Delta U(k + j|k)\leqslant\Delta {U_{\max}} 。\end{split} $

(20)

其中： $ {N_P} $ 、 $ {N_C} $ 分别为MPC算法中的预测步数和控制步数，且 $ {N_C} \leqslant {N_P} $ ； $ {Q_{err}}、{Q_{in}}、{Q_{din}} $ ，均为权重矩阵。终端系数 $ {q_{ter}} $ 的设置，是为了最优化问题在预测时域终点有可行解。 $ r(k + j) $ 为参考位置，根据已知路径进行选择，对动力定位问题 $ r(k + j) $ 在时域内设为常数，此处为简化问题设置为[0,0,0°]。考虑到在动力定位问题中，将推力和推力角度加入到目标函数中，主要是为提前考虑推进器的物理限制，以及使推进器的运动更加平缓，在此处可对权重矩阵 $ {{\boldsymbol{Q}}_{in}}、{{\boldsymbol{Q}}_{din}} $ 进行预先设置，即在整个动力定位过程中并不进行变化，而位置权重 $ {Q_{err}} $ 是动力定位过程中最重要参数。在不同位置偏差和不同速度时，需进行一定修正变化。如位置偏差大、速度变化大时，权重 $ {Q_{err}} $ 需适当增大，当DP船受环境干扰较小，较为平稳时，则需注重降低DP船的输出功率。此时 $ {Q_{err}} $ 适当减小，相对的 ${Q_{in}}、{Q_{din}}$ 就会增大，有利于DP船的平稳和节约能源。为应对DP船水平面三自由度，假定 $ {Q_{err}} = {\rm{diag}}[{q_1} {q_2}{\text{ }}{q_3}] $ ，以 $ {q_3} $ 为例进行模糊化规则设计，选择 $ {q_3} $ 是因为对首摇的控制相对横荡和纵荡更加困难。

1.4 模糊化规则设计

本文采用的模糊化方法是将模糊化规则和控制表对应，正常来说，需对输入输出给出模糊集和论域，并在此基础上确定隶属函数。参考Zhao^[12]对角度偏差 $ {\psi _e} $ ，角度偏差变化速度即角速度 $ r $ ，以及输出参数 $ {q_3} $ 的模糊集以及隶属函数分别规定如下：

$ {\psi _e} $ 的模糊集为{NB,NM,NS,ZO,PS,PM,PB}。其中，P为正，N为负，ZO为0，B为大，M为中，S为小，其隶属函数如图1所示。

$ r $ 的模糊集为{N,P}，其中对角度偏差速度进行了最简化处理，只选择了正负，这是因为在动力定位中DP船通常变化缓慢，只需进行正负区分来和 $ {\psi _e} $ 进行匹配，隶属函数如下：

$ {u_{rP}} = \left\{ \begin{aligned} &{\text{0}}，r < 0 ，\\ &{\text{1}} r \geqslant 0。\end{aligned} \right. $

(21)

$ {u_{rN}} = \left\{ \begin{aligned} & 1，r < 0 ，\\ & 0，r \geqslant 0 。\end{aligned} \right. $

(22)

$ {q_3} $ 的模糊集选择为{S,B},即可变参数有大或者小2种模糊情况。为了能设计 $ {q_3} $ 的隶属函数，假设 $ {q_3} $ 的变化存在一个范围，范围为 $ \left[ {{q_{3\min }}{\text{ }}{q_{3\max }}} \right] $ ，此处的范围通常由实验获得，引入 $ {q_{3x}} \in [0,1] $ 。

$ {q_{3x}} = \frac{{{q_3} - {q_{3\min }}}}{{{q_{3\max }} - {q_{3\min }}}}。$

(23)

此时，对 $ {q_3} $ 的处理转化成了对 $ {q_{3x}} $ 的处理，有模糊集{S,B}，与之隶属函数如图2所示。图中 $ u $ 与 $ x = {q_{3x}} $ 函数表达关系如下：

$ \left\{ \begin{aligned} & {u_{Small}} = - \frac{1}{4}\ln x，\\ & {u_{Big}} = - \frac{1}{4}\ln (1 - x) 。\end{aligned} \right. $

(24)

需注意， $ u $ 与 $ x = {q_{3x}} $ 的范围均属于[0,1]。

模糊规则本质上是一系列从专家和实际经验中获得的if-then语句，例如当 $ {\psi _e} $ 为PB即正大， $ r $ 为P时，此时误差大，且误差速度与误差同方向，误差正不断变大，此时需较大的 $ {q_3} $ 参数来提升角度误差的权重，来抑制角度误差持续变大，可表述为if PB and P then B，如此即有如表1所示模糊控制表，共有14条规则，表中每一条规则都可如此解释。

本文采用的是单值模糊器和乘积推理机，因此表1的14条规则所对应的具体值都可如下式表示：

$ {u_i} = {u_{\psi ei}}*{u_{ri}}。$

(25)

同时：

$ \sum\limits_{i = 1}^{14} {{u_i} = 1} 。$

(26)

其中， $ {u_{\psi ei}} $ 根据相应角度偏差由图1得出， $ {u_{ri}} $ 则根据相应偏差变化值，此处即角速度值，由式（13）和式（14）得出。而第i条规则 $ {q_{3xi}} $ 的值则可根据 $ {u_i} $ 从图2得出，由此可求出 $ {q_{3x}} $ 如下式：

$ {q_{3x}} = \sum\limits_{i = 1}^{14} {{q_{3xi}}*{u_i}} 。$

(27)

当得出 $ {q_{3x}} $ 后，根据式（15）可得：

$ {q_3} = ({q_{\max }} - {q_{\min }})*{q_{3x}} + {q_{\min }}。$

(28)

需要说明的是，该方法同样适用于横荡和纵荡。

2 仿真DP船参数和环境载荷

DP模型船主要参数如表2所示^[13]，螺旋桨布置如图3所示。

环境干扰参照Fossen^[14]选择如下：

$ {\tau _{env}} = {R^{\rm{T}}}(\psi )b 。$

(29)

基于一阶马尔科夫模型构建低频环境力：

$ \dot b = - {\boldsymbol{T}}_{env}^{ - 1}b + \varGamma \varpi 。$

(30)

其中， ${T_{env}}、\varGamma$ 为可以指定的对角时不变矩阵； ${T_{env}}$ 为时间偏置正定矩阵。此模型可用来描述变化缓慢的环境外载荷，设置 ${{\boldsymbol{T}}_{env}} = {\rm{diag}}[100,100,100],{\boldsymbol{\varGamma}} = {\rm{diag}}[0.3, 0.2,0.1]$ ， $\varpi \in {R^{3 \times 1}}$ 为零均值高斯白噪声。

3 模拟仿真与讨论 3.1 非线性干扰观测器模拟

对设计的非线性观测器进行仿真模拟，设计参数 $ K = {\rm{diag}}[2.78{\text{ 2}}{\text{.78 2}}{\text{.78]}} $ ，由图4可看出，观测器对时变外干扰进行了有效跟踪，观测值和实际值基本一致，初步验证了本文所设计观测器的有效性。

为了进一步验证观测器的性能，图5给出了 $ {\tilde \tau _{env}} $ 的曲线。可看出观测误差，与图4进行比对，发现观测误差最大值发生在环境干扰变化增大的区间内，但观测器仍能把误差限制在可接受范围内。因此本文所设计干扰观测器满足使用精度要求，在后续使用时，当环境干扰变化增大时观测误差需重点关注。

3.2 定位能力比较

将对参数固定的非线性模型预测策略（Nonlinear Model Predictive Control，NMPC）和模糊自适应模型预测控制策略（Fuzzy MPC）进行比较，这2种策略的区别仅在于Fuzzy策略对NMPC策略中的部分固定参数，进行了模糊化处理，在控制器中增加模糊化模块，其他各个参数设计保持一致，故可选取各参数如下： $ {N_P} = 6，{N_C} = 3$ ， $ {Q_{in}} = {\rm{diag}}[100\;\;\;100\;\;100\;\; 100\;\;0\;\;0\;\;0\;\;0] $ ， $ {Q_{din}} = {\rm{diag}}[1\;\,1\,\;1\,\;1\,\;75\,\;75\,\;75\,\;75]，{q_{ter}} = 10 $ ，NMPC策略中选取 ${Q_{err}}= {\rm{diag}}[750\;\;750\;\;750]$ 。

从图6可明显看出，Fuzzy策略的定位精度更高，波动幅度也小很多，在时间步长150步之前，2种策略控制精度趋势较为一致，Fuzzy策略精度稍高，而150步之后，开始出现较为不同的趋势。从环境干扰可看出，在150步左右时 $ {\tau _{env1}} $ 和 $ {\tau _{env2}} $ 发生较大变化，而 $ {\tau _{env3}} $ 保持相对平稳，表现在定位精度上，则是2种策略出现较大差距。这是因为Fuzzy策略会随着环境干扰而变化，将适应调节与位置精度相关的参数，因此当DP船出现误差增大情况，控制策略中的模糊模块将自适应增大与位置精度相关的参数，参数的增大将有效抑制定位误差持续增大。表3给出了Fuzzy策略和NMPC策略的最大误差和误差标准差，X方向Fuzzy策略误差减小19.5%，Y方向Fuzzy策略误差减小19%，首摇角Fuzzy策略减小16.7%。从表3误差的标准差也可发现，Fuzzy策略的标准差也远小于NMPC策略，X方向，Y方向和首摇角Fuzzy策略分别减小了30.1%、46.4%、11.8%，这也说明了Fuzzy策略不仅在最大误差上小很多，在控制的平稳性上也更为优越。

从图（7）可看出2种策略的速度变化，NMPC策略的速度标准差，分别为0.0052 m/s、0.0046 m/s和0.2747°/s，Fuzzy策略的速度标准差为0.0053 m/s、0.0040 m/s，0.3889°/s，这2种策略在速度的平稳性上是近似的，而Fuzzy策略却能获得更好的定位精度，说明对参数的模糊化调控能更及时应对外界干扰。

为更进一步衡量动力定位的精度，引入参数，积分绝对误差（Integral Absolute Error，IAE）^[15]。

$ IAE(t) = \sum\limits_{k = 0}^{k = t} {\sqrt {{e^T}(k)e(k)} }。$

(31)

式中： $ e(k) = [{x_e}(k)/5,{y_e}(k)/5,{\psi _e}(k)/25] $ 为k时刻的标准化位置误差； $ {x_e}、{y_e}、{\psi _e} $ 分别代表在X方向，Y方向，首摇这3个维度的位置误差。相同模拟时域，IAE越小代表控制性能越好。从图8可看出2种策略的变化，在最开始时2种策略IAE基本一致，模拟结束时，NMPC策略的IEA值为46.47，Fuzzy策略的IEA值为34.12，精度提升为26.6%，说明了Fuzzy策略的有效性。

从图9可看出，2种策略在推力及推力变化上是相似的，由此可推断，在获得更好推力效果的同时，螺旋桨的功率变化并不大。为了能更进一步看出功率消耗的情况，现假定这里使用的螺旋桨均为定螺距桨，则功率和推力有如下关系：

$ P = \frac{{2\text{π} {K_Q}}}{{K_T^{3/2}{\rho ^{1/2}}D}}{T^{\frac{3}{2}}} = \zeta {T^{\frac{3}{2}}} 。$

(32)

其中：K_T为推进器的推力系数；K_Q为螺旋桨的扭矩系数；η为推进器的敞水效率；ρ为流体的密度；D为螺旋桨的直径； $ \zeta $ 为各种螺旋桨参数的简化，对2种策略分别进行计算；NMPC策略消耗总功率为214.65 $ \zeta $ ；Fuzzy策略消耗总功率为223.78 $ \zeta $ ，相对功率消耗提高4.25%，相比精度的提升，功率的消耗是值得的。

4 结　语

本文针对MPC应用于DP船时的参数整定困难问题，设计了模糊自适应模型预测控制策略，引入模糊算法来实现模型预测控制的参数自整定。在此基础上为进一步提升控制系统鲁棒性，设计了一种非线性干扰观测器用于估计未知时变环境力，从而实现对环境力的补偿。通过理论和模拟仿真验证，在相同的环境干扰和模拟时域内，Fuzzy策略精度提升显著。通过IEA值可看出，在模拟结束时，精度提升26.6%。同时反映定位平稳性的标准差也小于NMPC策略，说明新的控制策略控制的更加平稳，定位精度更高，所消耗功率仅提升4.25%，验证了该策略的有效性，为解决模型预测控制的参数整定问题提供了一定参考性。