0 引　言

随着科技的迅速发展，自主水下机器人（AUV）在多种领域中的应用越来越广泛，包括海洋勘探、环境监测、海底考古、水下结构物检查和军事应用等。AUV具有良好的运动控制能力，是实现各种水下任务的基础。

常见的AUV运动控制算法可以分为两大类：1）是依靠水下设备的实时误差数据进行控制的PID、ADRC控制以及基于这2种方法改进的控制器，如模糊PID、FOADRC控制器等；2）依靠机器人的动力学模型建立的相关控制器如滑模控制、自适应控制、模型预测控制等。2种方法均有一定局限性，依靠误差模型的方法随着环境信息的改变，一般控制参数也要随之进行改变，依靠水动力的方法由于水下力学信息的复杂性，难以获得精准的模型参数且计算复杂性较高。近年来，有很多改进的控制算法也被尝试于AUV的运动控制中，如Jagtap等^[1]设计了基于动态缩减模型的非线性模型预测控制器；Parijat等^[2]通过改进的立方体卡尔曼滤波算法设计了随机非线性模型预测控制方案；Steenson等^[3]设计添加执行器约束模型预测控制算法，在Delphin2 AUV并上验证了算法控制的有效性；Gomes等^[4]尽可能多地替代传统方案的在线计算负担，降低了控制系统的复杂性；Luca等^[5]设计了一种AUV控制器，控制器由点发生器、预测最优规划器、潮流估计器组成；刘昌鑫等^[6]路径参数的二阶导数作为路径虚拟控制律，将参考路径模型扩展到AUV路径跟踪预测模型中。然而，以模型预测控制为核心的控制器都面临着模型复杂不确定、随着预测步数的增加计算量庞大的问题。武建国等^[7]通过设计非线性干扰观测器（NDO）并引入滑模控制改进了AUV的三维轨迹追踪；王晓伟等^[8]使用Sigmoid函数优化了滑膜控制器，提高了AUV转弯时的跟踪精度；姚金艺等^[9]利用智能积分S面的滑模控制完成了AUV的三维路径跟踪。然而，由于AUV模型存在不确定性和参数摄动性，这些基于精确模型的滑膜方法在实际应用中的性能可能会受限。高鹏等^[10]提出一种双曲正切型航速调整策略，用以MFAC航向控制器路径点切换处的平滑过渡，但对于外部噪声的影响较为敏感。李泽宇等^[11]采用Q学习方法，对滑膜控制参数进行离线训练化，避免“维数灾”现象，其仿真结果验证了其方法的有效性；赵杰梅等^[12]提出神经网络自适应输出反馈控制器，该控制器包含动态反馈控制器、神经网络控制器、鲁棒控制器3个部分，并利用仿真分析验证了算法的有效性。但存在着神经网络对数据的依赖性较高且训练时间较长的问题。

针对上述外部环境不稳定、计算需求大、模型依赖度高的问题，本文设计以动态矩阵控制方法为核心的串级控制器用于AUV的航行控制。其目的是降低模型依赖度的同时提高控制系统的鲁棒性。在实际应用可以采用此控制方法用于各类以误差为控制参数的AUV控制系统中，提升其控制系统的鲁棒性。

1 AUV运动建模

本文研究的“云帆-07”型欠驱动AUV是在惯性坐标系$ E-{X}_{e}{Y}_{e}{Z}_{e} $和载体坐标系$ B-xyz $下进行的，惯性坐标系的原点选取为AUV发射点，载体坐标系的原点选取为AUV的浮心，AUV的重心与浮心重合（见图1）。

AUV的数学模型可表示为：

$ \left\{\begin{array}{l}\dot{\eta }={\boldsymbol{R}}(\eta )\nu \text{，}\\ {\boldsymbol{M}}\dot{\nu }+{\boldsymbol{C}}(\nu )\nu +{\boldsymbol{D}}(\nu )\nu +{\boldsymbol{g}}({\boldsymbol{\eta }})=\tau +{\tau }_{{d}}。\end{array} \right.$

(1)

式中：$ \eta ={\left[x,y,z,\theta ,\psi \right]}^{\rm T} $为广义方位参数，分别表示在惯性坐标系下AUV的纵坐标、横坐标和垂向坐标，俯仰角和偏航角；$ \nu ={\left[u,v,\omega ,q,r\right]}^{\rm T} $为广义速度参数，分别表示纵向、横向和垂向速度，纵倾和偏航角速度；$ {\boldsymbol{R}}\left({\boldsymbol{\eta}} \right) $为载体坐标系到惯性坐标系的转换矩阵。$ {\boldsymbol{M}} $为AUV的惯性矩阵，包含附加质量矩阵和刚体惯性矩阵；$ {\boldsymbol{C}}\left({\boldsymbol{\nu }}\right) $科氏力矩阵和向心力矩阵；$ {\boldsymbol{D}}\left({\boldsymbol{\nu}} \right) $是AUV在流体环境中的阻尼矩阵；$ {\boldsymbol{g}}\left({\boldsymbol{\eta}} \right) $为AUV系统恢复力和恢复力矩矩阵；$ \tau ={[{\tau }_{u},\mathrm{0,0},{\tau }_{q},{\tau }_{r}]}^{\rm T} $为控制输入量，分别表示载体坐标系下对x轴水平推力、绕y轴的扭矩、绕z轴的扭矩；$ {\tau }_{d}={[{\tau }_{1},{\tau }_{2},{\tau }_{3},{\tau }_{4},{\tau }_{5}]}^{\rm T} $为外部环境产生的外部干扰项。忽略横滚对AUV的影响，AUV的数学模型展开可表示为：

$ \left\{\begin{array}{l}\dot{x}=u\mathrm{cos}\psi \mathrm{cos}\theta -v\mathrm{sin}\psi +w\mathrm{sin}\theta \mathrm{cos}\psi ,\\ \dot{y}=u\mathrm{sin}\psi \mathrm{cos}\theta +v\mathrm{cos}\psi +w\mathrm{sin}\theta \mathrm{sin}\psi ,\\ \dot{z}=-u\mathrm{sin}\theta +w\mathrm{cos}\theta ,\\ \dot{\theta }=q,\\ \dot{\psi }=r\mathrm{sec}\theta 。\end{array} \right.$

(2)

$ \left\{\begin{array}{l}{\dot{u} }={F}_{u}u+{F}_{X}{\tau }_{\text{u}}+{\tau }_{1},\\ \dot{v}={F}_{v}v+{\tau }_{2},\\ \dot{w}={F}_{w}w+{\tau }_{3},\\ \dot{q}={F}_{q}q+{F}_{M}{\tau }_{q}+{\tau }_{4},\\ \dot{r}={F}_{r}r+{F}_{N}{\tau }_{r}+{\tau }_{5}。\end{array} \right.$

(3)

式中：$ {F_u} = \displaystyle\frac{{{X_{u\left| {\left. u \right|} \right.}}\left| u \right|}}{{m - {X_{\dot u}}}}，{F_X} = \displaystyle\frac{1}{{m - {X_{\dot u}}}}，{F_{\text{v}}} = \displaystyle\frac{{{Y_{v\left| v \right|}}\left| v \right|}}{{m - {Y_{\dot v}}}} $，$ {F_{\text{w}}} = \displaystyle\frac{{{Z_{w\left| w \right|}}\left| w \right|}}{{m - {Z_{\dot w}}}}，{F_q} = \displaystyle\frac{{{M_{q\left| q \right|}}\left| q \right|}}{{{I_y} - {M_{\dot q}}}}，{F_M} =\displaystyle\frac{1}{{{I_y} - {M_{\dot q}}}} $，$ {F_r} = \displaystyle\frac{{{N_{r\left| r \right|}}\left| r \right|}}{{{I_z} - {N_{\dot r}}}} $，$ {F_N} = \displaystyle\frac{1}{{{I_z} - {N_{\dot r}}}} $；m为AUV的质量；$ {I}_{y}、{I}_{z} $为AUV的转动惯量；$ {X}_{\dot{u}}、{Y}_{\dot{v}}、{Z}_{\dot{w}}、{M}_{\dot{q}}、{N}_{\dot{r}} $为线性阻尼系数。

2 AUV控制器设计 2.1 动态矩阵控制

DMC控制利用被控对象的单位阶跃响应序列来构建模型向量$ {\boldsymbol{a}}=[{a}_{1},\cdots {a}_{N}] $，并根据模型向量a预测未来的输出值，其中，$ {a}_{i}={\boldsymbol{a}}\left(i{T }_{s}\right) $。在此公式中，$ {a}_{i} $为单位阶跃响应的采样值；T_s为采样周期；N为建模时域。DMC控制包括模型预测、滚动优化和反馈校正3个部分^[13]。

在滚动优化时，通过性能函数在预测时域内取得最优值，计算得到未来控制量增量。在k时刻性能函数一般表示为:

$ \begin{split} {\mathrm{min}} J\left(k\right)=&\sum _{i=1}^{{p}_{i}}{{q}_{i}\left[w\right(k+i)-\tilde{y}_{M}(k+i\left|k\right)]}^{2}+\\ &\sum _{i=1}^{M}{r}_{i}\Delta u{(k+i-1)}^{2} 。\end{split}$

(4)

式中：$ w(k+i) $、$ {\stackrel{~}{y}}_{M}(k+i|k) $分别为k+i时刻的系统状态期望值和预测值；$ \Delta u(k+i-1) $为(k+i−1)时刻的控制量增量；P为预测时域；M为控制时域；$ {q}_{i} $、$ {r}_{i} $为权系数，分别表示对跟踪误差及控制量变化的抑制程度；$ \tilde{y}_{M}(k+i|k) $在k时刻为$ \tilde{y}_{M}\left(k\right) $的预测值，其计算式为：

$ \tilde{y}_{M}\left(k\right) =\tilde{y}_{P0}\left(k\right) +A\Delta {u}_{M}\left(k\right)。$

(5)

式中：$\tilde{y}_{P0}\left(k\right) $为初始预测值向量，是k时刻起控制量保持不变的系统预测值，可根据已知信息得到；A为动态矩阵，是由单位阶跃响应系数$ {a}_{i} $组成的P×M阵。

$ \boldsymbol{A}=\left[\begin{array}{ccc}{\boldsymbol{a}}_{1}& \cdots & 0\\ \vdots& \ddots & \vdots\\ {\boldsymbol{a}}_{\boldsymbol{M}}& \cdots & {\boldsymbol{a}}_{1}\\ \vdots& \ddots & \vdots\\ {\boldsymbol{a}}_{\boldsymbol{P}}& \cdots & {\boldsymbol{a}}_{\boldsymbol{P}-\boldsymbol{M}+1}\end{array}\right] 。$

(6)

将式(4)代入式(5)，整理并写成向量形式，可得新的性能函数为：

$ {\mathrm{min}}J\left(k\right) = {\left\|{w}_{p}\left(k\right) - {\tilde{y}}_{P0}\left(k\right) - A\Delta {u}_{M}\left(k\right)\right\|}_{Q}^{2} + {\left\|\Delta {u}_{M}\left(k\right)\right\|}_{R}^{2} 。$

(7)

式中：$ {w}_{p}\left(k\right) $为系统期望向量，$ {w} _{p} \left( k \right) = {\left[w\left( k + 1 ), \cdots ,w( k + P \right)\right]}^{\rm T} $;Q和R为误差权矩阵和控制权矩阵，$ {\boldsymbol{Q}}={\mathrm{diag}}({q}_{1},\cdots ,{q}_{p}) $，$ {\boldsymbol{R}}={\mathrm{diag}}({r}_{1},\cdots ,{r}_{p}) $；$ \Delta {u}_{M}\left(k\right) $为系统控制量。

在k时刻，$ {w}_{p}\left(k\right) $、$ {\stackrel{~}{y}}_{P0}\left(k\right) $均为已知，使$ J\left(k\right) $取极小的$ \Delta {u}_{M}\left(k\right) $可通过极值必要条件求得。

$ \frac{{\mathrm{d}}J\left(k\right)}{{\mathrm{d}}\Delta {u}_{M}\left(k\right)}= 0 。$

(8)

式中：$ \Delta {u}_{M}\left(k\right)={({A}^{\rm T}QA+R)}^{-1}{A}^{\rm T}Q\left[{w}_{p}(k)-{\stackrel{~}{y}}_{P0}(k)\right] $。

$ \Delta {u}_{M}\left(k\right) $为控制量向量，采用滚动控制方法，取其中实时控制增量$ \Delta {u}\left({k}\right) $构成实际控制作用于对象，$ \Delta u\left(k\right) $为$ \Delta {u}_{M}\left(k\right) $的首元素。求出控制增量$ \Delta u\left(k\right) $后，实际控制量为：

$ u\left(k\right) =u(k-1)+\Delta u\left(k\right) 。$

(9)

通过作用$ u\left(k\right) $，利用预测式(5)，可算出在其作用下未来时刻的输出预测值：

$ \tilde{y}_{N1}\left(k\right) ={\stackrel{~}{y}}_{N0}\left(k\right) +a\Delta u\left(k\right) 。$

(10)

在计算下一时刻的优化控制量前，需先检测对象的实际输出并与模型预测该时刻的预测值相比较构成输出误差：

$ e(k+1) =y(k+1)-\tilde{y}_{1}(k+1|k) ，$

(11)

$ \tilde{y}_{cor}(k+1) ={\stackrel{~}{y}}_{N1}\left(k\right) +{\boldsymbol{h}}e(k+1) ，$

(12)

$ \tilde{y}_{N0}(k+1) ={\boldsymbol{S}}{\stackrel{~}{y}}_{cor}\left(k\right)。$

(13)

式中：h为校正向量，$ {\boldsymbol{h}}=[{h}_{1},\cdots {h}_{N}] $；S为位移矩阵，通过位移矩阵将新的较正值传递给下时刻新的模型预测值，完成反馈校正操作。

2.2 模型向量获取及参数整定

在DMC控制中，采样周期、建模时域、模型向量和动态矩阵由被控对象动态特性决定^[14]。

DMC以对象的单位阶跃响应的采样值$ \left\{{a}_{i}\right\} $作为算法的出发点，其数据的光滑性对控制性能的影响非常显著，通常对阶跃响应采样值进行平滑处理以减小噪声和扰动的影响，处理后得到的模型参数$ \left\{{a}_{i}\right\} $应是光滑变化的，否则将影响控制系统的稳定性和动态性能。基于此前提，本文获取“云帆”的模型预测向量，根据其数学模型中关于航向角和深度控制的阶跃响应函数，将在湖试实验中调试效果好的PID控制器接入至AUV的数学模型中，获取航向角与深度的参考模型向量。

在AUV模型推力设定为10 N，采样频率为1 Hz，无外部干扰的情况下，获得输出为图2和图3中的阶跃响应曲线中的参考数据。

图中，PID控制器输出的控制数据，前波后波比例为4∶1，且收敛时间较快，此2组参数是实验调试中获得良好的PID参数，因为动态矩阵控制中对于模型向量平滑性的要求，所以将2组实验原数据先进行移动平均法除波动，在使用指数平滑法让数据更为光滑，最后获得变化缓慢、趋势光滑的2组模型向量。

根据阶跃响应曲线的数据，确定建模时域为N，N即系统达到稳态的时间长度；优化时域P表示为对k时刻起未来多少步内的输出值与期望值，优化时域P必须超过对象阶跃响应的时滞部分或由非最小相位特性引起的反向部分，并覆盖对象阶跃响应的主要部分；而误差权矩阵Q作为权系数则反映了对不同时刻误差的重视程度，为了使控制系统稳定，必须选择P和Q满足下式：

$ \left(\sum _{i=1}^{P}{a}_{i}{q}_{i}\right){a}_{N} > 0 。$

(14)

控制时域M为优化变量的数目，反映了控制的自由度和能力；控制权矩阵R的作用是对控制增量的剧烈变化加以适度的抑制。根据P、M、Q、R参数的含义和阶跃响应曲线特性，结合“云帆-07”控制回路特性，确定参数如表1所示。

2.3 控制器设计

动态矩阵控制器的控制结构如图4所示，主要由预测函数、反馈校正、移位更新3个部分组合而成，在处理一些无干扰的控制环境可以较好地完成控制任务。然而，在使用动态矩阵控制方法去应用解决AUV的控制问题，常伴随着这样一些问题，如不能分辨误差是由模型适配还是扰动引起的，校正参数难以选择；反馈校正时，如果扰动到输出有较长的传递时间，则控制对扰动的反应难以及时。

针对上述问题，借用串级控制的思想，把预测控制结合串级控制结构予以实现。首先，密集采样的PID控制器对于扰动有良好的抑制作用，故可在对象最易发生扰动的部分后取出相应的可测中间变量构成PID闭环结构。这一层控制采用比DMC控制高得多的频率^[15]，其目的在于快速有效地抑制突发性的扰动。这一控制回路和对象的剩余部分可视为一广义对象，在外层，再对这一广义对象实施通常的DMC控制。所设计的控制结构如图5所示。

在AUV的航行控制中，输入为期望目标函数，如期望航向追踪或深度追踪函数等，DMC为主控制器，处理更高层次的数据过程，如x、y、z位移信息，PID为从控制器主要处理AUV姿态信息的变化，如俯仰角、航向角。内环与外环的控制计算频率比为2∶1。

3 仿真与分析

为验证本研究设计的基于DMC的串级控制器的控制性能和在动态干扰下的鲁棒性，根据“云帆”AUV的数学模型进行仿真实验。实验从2个面的控制效果来验证控制器性能。在2个控制面分别设置期望航向与期望深度，观测控制器到达稳态前的调节过程与到达稳态后施加外部干扰后的抗干扰能力。定义AUV的初始位置向量为$ \eta =\left[0\ 0\ 0\ 80\ 0\right] $；定义AUV的初始速度向量为$ \dot{\eta }=\left[0\ 0\ 0\ 0\ 0\right] $；定义环境干扰函数为$ {{{y}}_{\text{1}}}(t) = 10 \times \sin (0.1\text{π} t) $，$ {{{y}}_{\text{2}}}(t) = 5 \times \sin (0.1\text{π} t) $。

3.1 水平面控制验证

设置期望航向角为90°，在系统到达稳态后，在200 s处引入外部时变干扰y₁。图6和图7为模型的航向角与横向速度控制情况。

控制器在控制航向任务时，大体趋势可以描述为，DMC-PID串级控制器相比较于PID、串级PID控制器在初始响应速度方面大抵相同、系统超调量有所减少、系统稳定的周期到达相对较快。最后，从横向速度的数据可以看出，DMC-PID的横向速度波动趋势符合航向角度控制的控制规律。

由图8航向角控制误差显示的信息可以得出，3种控制器的平均航向误差按照图序分别为：±1.8326°、±0.4577°、±0.2917°，此数据说明3种控制器在面对干扰时的不同抵抗能力。图9中对应系统的横向速度信息也反馈了其误差的趋势，其震荡的趋势与误差偏移的趋势相同。

3.2 垂直面控制验证

设置期望深度为水下1 m，在系统到达稳态后，在200 s处引入外部时变干扰y₂，图10和图11为模型深度与俯仰角控制情况。

可以看出，控制器在干扰环境中执行定深任务时，控制效果总体趋势可描述为，3种控制器在到达目标深度时的系统超量值DMC-PID控制表现良好，但其波动范围差距较小且波动数量一致，DMC-PID控制到达稳态的时间更快。最后，由系统俯仰角波动的数据可以看出，DMC-PID的俯仰角正向震荡值位于PID与串级PID之间，但其y轴的负向震荡值低于其他，且到达稳态的时间较快。

由图12定深控制误差显示的信息可以得出以下信息：3种控制器的定深波动误差范围按照图序分别约为：±0.6 m、±0.3 m、±0.2 m，图13中对应系统的俯仰角状态信息也反馈了其运动的趋势，其俯仰角波动范围分别约为0.52°±0.41°、±0.32°。

3.3 控制器性能分析

本研究采用积分时间绝对误差（ITEA）、积分时间平方误差（ISE）对控制器从开始至到达稳态的性能进行评估，数值越小代表控制器的性能越好，其中：

$ ITAE = \int {t\left| {e(t)} \right|} {\mathrm{d}}t，$

(15)

$ ITSE = \int {{\text{t}}{{[e(t)]}^2}} {\mathrm{d}}t 。$

(16)

采用积分绝对误差（IAE）与积分平均误差（ISE）来评估控制器处于干扰状态时的控制器性能，数值越小代表控制器的性能越好，其中：

$ IAE = \int {\left| {{\text{e}}(t)} \right|{\mathrm{d}}t}，$

(17)

$ ISE = {\int {\left[ {e(t)} \right]} ^2}{\mathrm{d}}t 。$

(18)

将实验数据代入式中，获得结果如图14、图15所示。

可知，在启动至稳态的控制器参数评价中，定向任务时的评价参数较为明显，按照图序，2种参数依次递减，在定深任务时，PID控制器评价参数高于其他，DMC-PID的评价参数略低于串级PID；在动态干扰下的参数评价中，定向与定深评价参数均依次递减。

通过控制器的评价参数以及前文的验证分析，得到结论，AUV系统启动至稳态的过程中，DMC-PID控制器性能较为优良，尤其于定向任务时效果更加明显；在动态干扰环境时，其抵抗干扰的能力更为优良，证明了其控制鲁棒性更高。

且通过仿真实验可以证明，此串级控制器可以满足AUV的控制需求，可以良好的完成航行任务。

4 结　语

针对欠驱动AUV的航行控制问题，基于模型预测控制策略中的动态矩阵控制方法，结合串级控制的控制思想，设计了以动态矩阵为核心的串级控制器，并且进行了定向、定深控制任务的仿真验证，实验结果表明DMC-PID控制器能够有效地提高传统PID控制方法的抗干扰能力，并且具有良好的控制性能，能够满足AUV的航行任务控制需求。