下载:
-
当小型水下机器人(Autonomous Underwater Vehicle，AUV)执行任务时，能耗是影响其航行能力的主要因素，能耗主要取决于航行阻力的大小。当AUV在近水面航行时，海浪引起的纵摇运动会导致航行阻力的增加。AUV携带的能量有限，阻力的增加会导致其续航力的减少。
国际海事组织(International Maritime Organization，IMO)提出了一种新的船舶能效设计指标(Energy Efficiency Design Index，EEDI)[1]，这就要求在设计AUV的纵摇减摇控制器时，考虑AUV的节能问题。控制策略不仅要追求运动幅值最小，同时还要考虑驱动装置能量和纵摇增阻的优化[2]。基于辐射能量理论，船舶横摇的幅值将会增加航行阻力，因此在设计船舶横摇减摇鳍控制器的过程中，需要考虑由于横摇所导致的横摇增阻[3]。但是在设计AUV纵摇控制器的过程中，纵摇增阻几乎没有被涉及，因此，需提出AUV的纵摇增阻模型，并将其作为AUV纵摇控制的性能指标。
随着AUV应用得越来越广泛，关于AUV的研究也日趋深入。C-SCOUT AUV被用于海洋地理的采样工作，其动力学模型分析了其外部受力和控制翼面的水动力[4]。加拿大不列颠哥伦比亚大学研制了Dolphin MK Ⅱ半潜式AUV，通过实验方法，采用四分之一模型研究了AUV整体和控制翼面的水动力[5]。鉴于该AUV为半潜式AUV，其将在近水面工作并受到海浪力的干扰，故为保持稳定的航行姿态，需要设计纵摇控制器。
丁团结[6]针对飞机进场着陆问题进行了LQR控制器设计，解决了多控制回路的协调问题。李一波等[7]采用一种增广LQR方法设计了飞翼式无人机纵向姿态控制律，该方法将系统的输出误差和外界恒值阵风干扰信号引入到了性能指标函数中。段镇[8]对无人机的俯仰速率回路采用了鲁棒伺服LQR控制器。Bhushan等[9]对基于遗传算法的LQR控制方法的权矩阵进行了分析，认为基于遗传算法的性能指标更为稳定。H?usler等[10]在AUV的避碰路径规划中加入了能量优化的性能指标。Liang等[11]对AUV在有外干扰时的能量优化进行了分析。Petrich等[12]对AUV的纵轴简化问题进行了分析。Wang等[13]对零航速减摇的能量优化问题进行了分析，但是并未提出航行阻力的性能指标。徐建安等[14]对AUV的艏摇控制进行了能量优化，降低了艏摇速度的抖动。Chyba等[15]对AUV在两点之间的航行进行了控制方法的能量优化。Sarkar等[16]设计了一种基于能量优化的滑模控制方法，取得了较好的效果。李佳等[17]在设计潜器时，对潜器的航行阻力进行了计算。
可以看出，在以往的研究中，对于航行阻力的研究往往会在船体设计的过程中涉及，通过型线的优化来减小航行阻力，但是对纵摇运动所引起的航行阻力增加的问题研究较少。事实上，AUV的纵摇运动会导致其航行阻力的增加，控制翼面的摆动也会增加其航行阻力，航行阻力和纵摇运动及翼面摆动均呈正相关的关系，即纵摇运动越剧烈航行阻力增加越大，翼面摆动越大航行阻力增加也越大。然而，控制翼面的摆动幅值越大减纵摇的效果越明显。因此，在控制AUV的纵摇运动时，应综合考虑纵摇及控制翼面的增阻问题，使二者折衷从而使AUV的航行阻力最小，进而减小推进能量的消耗。
为实现这一目标，首先需研究航行阻力与自主式航行器姿态之间的关系，建立航行阻力与航行姿态之间的数学模型，其次，对控制翼面的增阻模型也需进行考虑。在设计控制器时，应将二者综合考虑并作为性能指标，以便使控制过程中的航行阻力最小。
本文将以Dolphin MK Ⅱ为研究对象，研究其纵摇控制策略。首先，将建立AUV的纵摇增阻模型，定量分析纵摇姿态与航行阻力之间的关系。其次，将分析翼面摆动对航行增阻的影响。最后，将纵摇增阻和翼面增阻综合考虑，并以此作为性能指标，进行AUV纵摇姿态控制器设计。设计一种LQR控制器，将纵摇增阻和翼面增阻综合起来，以达到航行姿态阻力最优。
-
AUV在近水面航行时，会受到海浪的扰动，产生纵摇运动并增加其航行阻力。图 1为AUV的结构图，本文以加拿大的Dolphin MK Ⅱ AUV为模型，建立其非线性模型，并对其受到的海浪干扰力/干扰力矩进行仿真。为控制AUV的航行姿态，首先建立AUV的六自由度非线性模型，为便于LQR控制器的设计，本文对其进行线性化。控制翼面是减纵摇的驱动装置，控制翼面采用NACA0025型翼面，并建立了控制翼面的水动力模型。
图 1 AUV的结构图
Figure 1. The structure of AUV
在对AUV进行力的分析时，可以把AUV看作是质量均匀分布的刚体，根据牛顿定律和流体力学的相关知识，从而得到AUV近水面六自由度运动方程。Dolphin MK Ⅱ AUV的参数如表 1所示。
表 1 AUV的参数
Table 1. Parameters of AUV
全长/m 质量/kg 直径/m 尾长/m 排水量/kg 最大航速/kn 8.534 4 300 1 2.35 4 600 18 由动力学知识可知，AUV的动力学模型如下：
$$ {\mathit{\boldsymbol{M}}_{{\rm{rb}}}}\frac{\partial }{{\partial t}}\mathit{\boldsymbol{v + }}{\mathit{\boldsymbol{C}}_{{\rm{rb}}}}\left( v \right)\mathit{\boldsymbol{v = }}{\mathit{\boldsymbol{\tau }}_{{\rm{rb}}}} $$ (1) 式中：Mrb为质量矩阵；Crb(v)为科氏矩阵；τrb为作用于AUV的外力/外力矩；v为AUV的速度矩阵。
一般认为τrb由3部分组成，即
$$ {\mathit{\boldsymbol{\tau }}_{{\rm{rb}}}} = {\mathit{\boldsymbol{\tau }}_{\rm{h}}} + {\mathit{\boldsymbol{\tau }}_{\rm{w}}} + {\mathit{\boldsymbol{\tau }}_{\rm{c}}} $$ (2) 式中：τh为水动力/力矩；τc为控制力/力矩；τw为海洋扰动力/力矩。
作用于AUV的水动力和水动力力矩为
$$ {\mathit{\boldsymbol{\tau }}_{\rm{h}}} = - {\mathit{\boldsymbol{M}}_{\rm{A}}}\mathit{\boldsymbol{\dot v}} - {\mathit{\boldsymbol{C}}_{\rm{A}}}\mathit{\boldsymbol{v}} - \mathit{\boldsymbol{Dv + }}{\mathit{\boldsymbol{\tau }}_{{\rm{bg}}}} $$ (3) 式中：MA为附加质量矩阵；CA为附加科氏矩阵；τbg为恢复力/力矩； $\mathit{\boldsymbol{\dot v}}$ 为加速度矩阵；D为阻尼矩阵。
AUV产生的控制力/力矩可表示为
$$ {\mathit{\boldsymbol{\tau }}_{\rm{c}}} = \mathit{\boldsymbol{Bu}} $$ (4) 式中：B为控制矩阵；u为控制向量。
因此，可以得到AUV动力学模型的矩阵形式为
$$ \left( {{\mathit{\boldsymbol{M}}_{{\rm{rb}}}} + {\mathit{\boldsymbol{M}}_{\rm{A}}}} \right)\mathit{\boldsymbol{\dot v}} + \left( {{\mathit{\boldsymbol{C}}_{{\rm{rb}}}} + {\mathit{\boldsymbol{C}}_{\rm{A}}} + \mathit{\boldsymbol{D}}} \right)\mathit{\boldsymbol{v = }}{\mathit{\boldsymbol{\tau }}_{{\rm{bg}}}} + {\mathit{\boldsymbol{\tau }}_{\rm{w}}} + \mathit{\boldsymbol{Bu}} $$ (5) -
AUV在近水面航行时，会受到海浪和海流等海洋环境的干扰，具有不确定性和随机性。其中，海浪扰动是造成AUV摇荡的主要因素。在设计AUV姿态控制器时，要有效地估计AUV近水面航行时所受的干扰力/力矩。为了确定波浪对AUV的作用力/力矩，需要对波浪力进行仿真和计算。
近水面波浪垂荡、纵摇干扰力/力矩的仿真采用线性叠加的方式。首先，选定海浪的频谱形式，本文仿真选用ITTC单参数谱。然后，采取离散化的方法，根据离散的海浪频谱求得海浪的谐波幅值，初相角为随机数，可以根据随机函数给出。确定了各个谐波和初相角之后，把各个谐波叠加起来得到仿真的长峰波海浪，根据所得的谐波函数求得海浪各方向的速度和加速度分量。最后，通过积分的方式得到海浪对AUV的垂荡、纵摇干扰力/力矩。
海浪理论研究表明，不规则的海浪可以分解成大量均匀微小的规则波。研究AUV在不规则海浪中的运动时，可以用线性叠加的方式对海浪的干扰力/力矩进行研究。由流体力学知识可知，平面进行波在地面坐标系o-xyz中的一点(x, y, 0) 所处的波面方程为
$$ {\eta _{\rm{w}}} = {\alpha _{\rm{w}}}\cos \left( {{\omega _{\rm{w}}}t - {k_x}x - {k_y}y} \right) $$ (6) 式中：αw为海浪的幅值；ωw为海浪角频率；kx=kw cos ψw，ky=kw sin ψw，其中，kw为波数，ψw为浪向角。
AUV在近水面航行时，上面所述的海浪周期与AUV航行中遭遇的波浪周期并不相等。AUV航行时遭遇的波浪周期称为遭遇周期。在AUV运动计算中，应该用遭遇角频率代替海浪角频率。遭遇角频率ωe与海浪角频率ωw之间的关系为
$$ {\omega _{\rm{e}}} = {\omega _{\rm{w}}} - \frac{{{k_{\rm{w}}}{U_{\rm{w}}}}}{{\rm{g}}}\cos \left( {{\beta _{\rm{w}}}} \right) $$ (7) 式中：Uw为海浪的合速度；g为重力加速度；βw为遭遇角。
通常假定长峰波是由很多不同波峰和波长的规则波线性叠加而成，初始相位θj是0~2π之间的随机变量，则得到不规则长峰波的数学表达式为
$$ \begin{array}{*{20}{c}} {{\eta _{\rm{w}}}\left( {x,y,t} \right) = \sum\limits_j^N {{\alpha _{{\rm{w}}j}}\cos \left( {{\omega _{{\rm{w}}j}}t - {k_{xj}}x - {k_{yj}}y + {\theta _j}} \right),} }\\ {j = 1,2, \cdots ,N} \end{array} $$ (8) 式中，αwj为第j个谐波的波幅。
采用1966年第11届国际船模试验水池会议(International Towing Tank Conference，ITTC)推荐的波谱，一般称为ITTC单参数波谱，其表达式如下：
$$ S\left( {{\omega _{\rm{w}}}} \right) = \frac{{8.1 \times {{10}^{ - 3}} \cdot {g^2}}}{{\omega _{\rm{w}}^5}}w = \exp \left( { - \frac{{3.11}}{{{h_{1/3}}\omega _{\rm{w}}^4}}} \right) $$ (9) 式中，h1/3为海浪的有义波高。
谐波的各次谐波的波幅可由下式求得：
$$ {\alpha _{\rm{w}}} = 2\sqrt {S\left( {{\omega _{\rm{w}}}} \right){\rm{d}}{\omega _{\rm{w}}}} $$ (10) 在海浪坐标系下，长峰波海浪谐波的速度分量如下：
$$ {u_{{\rm{w}}j}} = \frac{{2\pi {\alpha _j}}}{{{T_j}}}\exp \left( { - {k_j}z} \right) * \cos \left( {{\omega _{{\rm{w}}j}}t - {k_{xj}}x - {k_{yj}}y + {\theta _j}} \right) $$ (11) $$ {w_{{\rm{w}}j}} = \frac{{2\pi {\alpha _j}}}{{{T_j}}}\exp \left( { - {k_j}z} \right) * \sin \left( {{\omega _{{\rm{w}}j}}t - {k_{xj}}x - {k_{yj}}y + {\theta _j}} \right) $$ (12) 加速度分量如下：
$$ {{\dot u}_{{\rm{w}}j}} = \frac{{4{\pi ^2}{\alpha _j}}}{{T_j^2}}\exp \left( { - {k_j}z} \right) * \sin \left( {{\omega _{{\rm{w}}j}}t - {k_{xj}}x - {k_{yj}}y + {\theta _j}} \right) $$ (13) $$ {{\dot w}_{{\rm{w}}j}} = - \frac{{4{\pi ^2}{\alpha _j}}}{{T_j^2}}\exp \left( { - {k_j}z} \right) * \cos \left( {{\omega _{{\rm{w}}j}}t - {k_{xj}}x - {k_{yj}}y + {\theta _j}} \right) $$ (14) 式(11)~式(14) 中：uwj， $\dot u$ wj分别为海浪传播方向海浪谐波的速度和加速度分量；wwj， $\dot w$ wj分别为垂直于海平面方向的海浪谐波的速度和加速度分量。
根据叠加定理和坐标转换，可以得到地面坐标系o-xyz下的总速度分量和加速度分量。
总速度分量：
$$ {\left[ {{u_{\rm{w}}}} \right]_{\rm{E}}} = \left[ {\sum\limits_N {{u_{{\rm{w}}j}}} } \right]\cos \left( {{\psi _{\rm{w}}}} \right) $$ (15) $$ {\left[ {{v_{\rm{w}}}} \right]_{\rm{E}}} = \left[ {\sum\limits_N {{u_{{\rm{w}}j}}} } \right]\sin \left( {{\psi _{\rm{w}}}} \right) $$ (16) $$ {\left[ {{w_{\rm{w}}}} \right]_{\rm{E}}} = \left[ {\sum\limits_N {{w_{{\rm{w}}j}}} } \right] $$ (17) 总加速度分量：
$$ {\left[ {{{\dot u}_{\rm{w}}}} \right]_{\rm{E}}} = \left[ {\sum\limits_N {{{\dot u}_{{\rm{w}}j}}} } \right]\cos \left( {{\psi _{\rm{w}}}} \right) $$ (18) $$ {\left[ {{{\dot v}_{\rm{w}}}} \right]_{\rm{E}}} = \left[ {\sum\limits_N {{{\dot u}_{{\rm{w}}j}}} } \right]\sin \left( {{\psi _{\rm{w}}}} \right) $$ (19) $$ {\left[ {{{\dot w}_{\rm{w}}}} \right]_{\rm{E}}} = \left[ {\sum\limits_N {{{\dot w}_{{\rm{w}}j}}} } \right] $$ (20) 式(15)~式(20) 中：[uw]E，[ $\dot u$ w]E分别为x轴方向的速度和加速度分量；[vw]E，[ $\dot v$ w]E分别为y轴方向的速度和加速度分量；[ww]E, [ $\dot w$ w]E分别为z轴方向的速度和加速度分量。
由于AUV的直径相对于海浪的波长比较小，所以可以将AUV看作一细长圆柱体。而且AUV所受的波浪力是由流动在AUV周围流体的分流所引起，而非因为海浪的衍射力作用。由此可知，AUV在近水面航行时所受到的总海浪干扰力可以用Morison方程求得，该方程如下：
$$ {F_{{\rm{w}}x}} = \left[ {\frac{{{C_{\rm{d}}}}}{2}\rho AU_{\rm{w}}^2 + {C_{\rm{m}}}\nabla \rho {{\dot U}_{\rm{w}}}} \right]{\rm{d}}x $$ (21) 式中：Cd为阻力系数；Cm为附加质量系数；ρ为海水密度；A为投影面积；∇为排水量； $\dot U$ w为海浪合加速度。应用Sarpkaya的波浪力观点，考虑Dolphin MK Ⅱ半潜式AUV的操作环境，Cd和Cm分别取值为0.65和1.95。
使用Morison方程沿船长进行积分，可以得到近水面处波浪对小型AUV的垂荡干扰力Zwave、纵摇干扰力矩Mwave分别为：
$$ {Z_{{\rm{wave}}}} = \int\limits_L {\left( {{C_{\rm{d}}}\frac{{\rho {D_0}}}{2}{{\left( {{w_{\rm{w}}} - {w_j}} \right)}^2} + {C_{\rm{m}}}\frac{{\rho \pi D_0^2}}{4}\left( {{{\dot w}_{\rm{w}}} - {{\dot w}_j}} \right)} \right){\rm{d}}x} $$ (22) $$ {M_{{\rm{wave}}}} = \int\limits_L {\left( {{C_{\rm{d}}}\frac{{\rho {D_0}}}{2}{{\left( {{w_{\rm{w}}} - {w_j}} \right)}^2} + {C_{\rm{m}}}\frac{{\rho \pi D_0^2}}{4}\left( {{{\dot w}_{\rm{w}}} - {{\dot w}_j}} \right)} \right)x{\rm{d}}x} $$ (23) 式中：D0为直径; ww和 $\dot w$ w分别为船体坐标系下波浪垂荡的速度和加速度；wj， $\dot w$ j分别为AUV沿船长方向的速度和加速度分量。
可应用坐标变换将[uw]E，[vw]E和[ww]E转换到船体坐标系下：
$$ {\left[ \begin{array}{l} {u_{\rm{w}}}\\ {v_{\rm{w}}}\\ {w_{\rm{w}}} \end{array} \right]_{\rm{B}}} = {\mathit{\boldsymbol{\underline Q}} ^{ - 1}}{\left[ \begin{array}{l} {u_{\rm{w}}}\\ {v_{\rm{w}}}\\ {w_{\rm{w}}} \end{array} \right]_{\rm{E}}} $$ (24) 式中，Q为坐标变换矩阵。
式(22) 和式(23) 中的w和 $\dot w$ 可以用下列各式求得：
$$ {w_j} = w - q{x_j} $$ (25) $$ {{\dot w}_j} = \dot w - \dot q{x_j} - uq + vp $$ (26) 式中：w和 $\dot w$ 分别为AUV的垂荡速度和加速度；q和 $\dot q$ 分别为纵摇角速度和加速度；u为纵荡速度；v为横荡速度；xj为船体坐标系下的坐标。
-
在研究船舶运动控制系统时，因为船舶运动受到控制，故认为船舶受到扰动作用后在平衡位置附近做小幅度的运动，则可以在平衡位置附近对船舶运动进行线性化处理。
AUV的运动与水面运动体不同，它除了在海平面内做水平运动外，还必须在可航行的深度范围内做升沉运动。AUV的六自由度运动可分为水平面运动和垂直面运动，本文只考虑AUV的垂直面运动。垂直面运动参数取{u, w, q, θ, z}，忽略水平面的运动，即v=p=r=φ=0，忽略w, q的二次项，可以得到简化的AUV垂直面运动方程：
$$ \begin{array}{*{20}{c}} {\left[ {\begin{array}{*{20}{c}} {m - {Z_{\dot w}}}&{ - {Z_{\dot q}}}&0&0\\ { - {M_{\dot w}}}&{{I_{yy}} - {M_{\dot q}}}&0&0\\ 0&0&1&0\\ 0&0&0&1 \end{array}} \right]\left[ \begin{array}{l} {\dot w}\\ {\dot q}\\ {\dot z}\\ {\dot \theta } \end{array} \right] = }\\ {\left[ {\begin{array}{*{20}{c}} {{Z_w}}&{\left( {m - {X_{\dot u}}} \right){u_0} + {Z_q}}&0&0\\ {\left( {{X_{\dot u}} - {Z_{\dot w}}} \right){u_0} + {M_w}}&{ - {Z_{\dot q}}{u_0} + {M_q}}&0&{{z_{\rm{b}}}B}\\ 1&0&0&{ - {u_0}}\\ 0&1&0&0 \end{array}} \right]\left[ \begin{array}{l} w\\ q\\ z\\ \theta \end{array} \right] + }\\ {\left[ {\begin{array}{*{20}{c}} {{Z_{{\delta _{\rm{b}}}}}}&{{Z_{{\delta _{\rm{s}}}}}}\\ {{M_{{\delta _{\rm{b}}}}}}&{{M_{{\delta _{\rm{s}}}}}}\\ 0&0\\ 0&0 \end{array}} \right]\left[ \begin{array}{l} {\delta _{\rm{b}}}\\ {\delta _{\rm{s}}} \end{array} \right] + \left[ {\begin{array}{*{20}{c}} {{Z_{{\rm{wave}}}}}\\ {{M_{{\rm{wave}}}}}\\ 0\\ 0 \end{array}} \right]} \end{array} $$ (27) 式中：m为AUV的质量；u0为AUV的航速；Z $\dot w$ ，Z $\dot q$ ，Zw，Zq为垂荡力水动力导数；X $\dot u$ 为纵荡力水动力导数；M $\dot w$ ，M $\dot q$ ，Mw，Mq为纵摇力矩水动力导数；Iyy为AUV的转动惯量；zb为浮心的坐标；B为浮力；Zδb，Zδs为控制翼面的垂荡力水动力导数；Mδb，Mδs为控制翼面的纵摇力矩水动力导数；δb为艏部控制翼面舵角；δs为艉部控制翼面舵角； $\dot z$ 为垂向速度； $\dot \theta $ 为地面坐标系下纵摇角速度；θ为纵摇角；z为垂向位移。
-
线性二次型最优控制系统(Linear Quadratic Regulator，LQR)是指，系统的状态方程为线性方程，性能指标是状态变量和控制变量的二次型函数的系统。最优控制属于线性系统综合理论中最具重要性和最具典型性的一类优化型综合问题。优化型综合问题的特点是需要通过对指定性能指标函数取极大或极小来导出系统的控制律。
设动态系统的状态方程为
$$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{\dot x}}\left( t \right) = \mathit{\boldsymbol{A}}\left( t \right)\mathit{\boldsymbol{x}}\left( t \right) + \mathit{\boldsymbol{B}}\left( t \right)\mathit{\boldsymbol{u}}\left( t \right)}\\ {\mathit{\boldsymbol{y}}\left( t \right) = \mathit{\boldsymbol{C}}\left( t \right)\mathit{\boldsymbol{x}}\left( t \right)} \end{array} $$ (28) 式中：x(t)为n维状态向量；u(t)为m维控制向量；y(t)为l维输出向量；A(t), B(t), C(t)为时变系数矩阵。设误差向量e(t)为
$$ \mathit{\boldsymbol{e}}\left( t \right) = \mathit{\boldsymbol{z}}\left( t \right) - \mathit{\boldsymbol{y}}\left( t \right) $$ (29) 式中，z(t)为l维期望输出向量。LQR控制器的目标是寻找最优控制向量u(t)，使得任意给定的初始状态x(t0)=x0可转移到自由端状态x(tf)，并使性能指标函数具有最小值。性能指标的形式如下：
$$ \begin{array}{*{20}{c}} {J\left( u \right) = \frac{1}{2}\mathit{\boldsymbol{e}}{{\left( {{t_f}} \right)}^{\rm{T}}}\mathit{\boldsymbol{P}}\left( {{t_f}} \right)\mathit{\boldsymbol{e}}\left( {{t_f}} \right) + }\\ {\frac{1}{2}\int_0^T {\left( {\mathit{\boldsymbol{e}}{{\left( t \right)}^{\rm{T}}}\mathit{\boldsymbol{Q}}\left( t \right)\mathit{\boldsymbol{e}}\left( t \right) + \mathit{\boldsymbol{u}}{{\left( t \right)}^{\rm{T}}}\mathit{\boldsymbol{R}}\left( t \right)\mathit{\boldsymbol{u}}\left( t \right)} \right){\rm{d}}t} } \end{array} $$ (30) 式中：P(tf)为l×l对称半正定常数矩阵；Q(t)为l×l对应状态变量的加权矩阵，是对称半正定矩阵；R(t)为m×m对应控制向量的加权矩阵，是对称正定矩阵。因为R(t)是正定的，所以当u(t)≠0时，有u(t)TR(t)u(t)＞0。因为Q(t)是半正定的，所以当e(t)≠0时，有e(t)TQ(t)e(t)≥0。
对于状态调节问题，性能指标的形式如下：
$$ \begin{array}{*{20}{c}} {J\left( u \right) = \frac{1}{2}{\mathit{\boldsymbol{x}}^{\rm{T}}}\left( {{t_f}} \right)\mathit{\boldsymbol{P}}\left( {{t_f}} \right)\mathit{\boldsymbol{x}}\left( {{t_f}} \right) + }\\ {\frac{1}{2}\int_0^T {\left( {{\mathit{\boldsymbol{x}}^{\rm{T}}}\left( t \right)\mathit{\boldsymbol{Q}}\left( t \right)\mathit{\boldsymbol{x}}\left( t \right) + {\mathit{\boldsymbol{u}}^{\rm{T}}}\left( t \right)\mathit{\boldsymbol{R}}\left( t \right)\mathit{\boldsymbol{u}}\left( t \right)} \right){\rm{d}}t} } \end{array} $$ (31) 使用最小值原理求解上述问题，需要引入拉格朗日乘子，用λ(t)表示，由此可以构成哈密顿函数：
$$ \begin{array}{*{20}{c}} {H\left( {\lambda ,x,u,t} \right) = \frac{1}{2}\left[ {{\mathit{\boldsymbol{x}}^{\rm{T}}}\left( t \right)\mathit{\boldsymbol{Q}}\left( t \right)\mathit{\boldsymbol{x}}\left( t \right) + {\mathit{\boldsymbol{u}}^{\rm{T}}}\left( t \right)\mathit{\boldsymbol{R}}\left( t \right)\mathit{\boldsymbol{u}}\left( t \right)} \right] + }\\ {{\mathit{\boldsymbol{\lambda }}^{\rm{T}}}\left( t \right)\left[ {\mathit{\boldsymbol{A}}\left( t \right)\mathit{\boldsymbol{x}}\left( t \right) + \mathit{\boldsymbol{B}}\left( t \right)\mathit{\boldsymbol{u}}\left( t \right)} \right]} \end{array} $$ (32) 要实现性能指标J为最小，最优解必须满足如下条件：
1) 正则方程组：
状态方程
$$ \mathit{\boldsymbol{\dot x}}\left( t \right) = \frac{\partial }{{\partial \lambda }}H\left( {\lambda ,x,u,t} \right) = \mathit{\boldsymbol{A}}\left( t \right)\mathit{\boldsymbol{x}}\left( t \right) + \mathit{\boldsymbol{B}}\left( t \right)\mathit{\boldsymbol{u}}\left( t \right) $$ (33) 协态方程
$$ \mathit{\boldsymbol{\dot \lambda }}\left( t \right) = - \frac{\partial }{{\partial x}}H\left( {\lambda ,x,u,t} \right) = - \mathit{\boldsymbol{Q}}\left( t \right)\mathit{\boldsymbol{x}}\left( t \right) - {\mathit{\boldsymbol{A}}^{\rm{T}}}\left( t \right)\mathit{\boldsymbol{u}}\left( t \right) $$ (34) 2) 控制方程：
$$ \frac{\partial }{{\partial u}}H\left( {\lambda ,x,u,t} \right) = \mathit{\boldsymbol{R}}\left( t \right)\mathit{\boldsymbol{u}}\left( t \right) + {\mathit{\boldsymbol{B}}^{\rm{T}}}\left( t \right)\mathit{\boldsymbol{\lambda }}\left( t \right) = 0 $$ (35) 即
$$ {\mathit{\boldsymbol{u}}^ * }\left( t \right) = - {\mathit{\boldsymbol{R}}^{ - 1}}\left( t \right){\mathit{\boldsymbol{B}}^{\rm{T}}}\left( t \right)\mathit{\boldsymbol{\lambda }}\left( t \right) $$ (36) 3) 初始状态：
$$ \mathit{\boldsymbol{x}}\left( {{t_0}} \right) = {\mathit{\boldsymbol{x}}_0} $$ (37) 4) 横截条件：
$$ \mathit{\boldsymbol{\lambda }}\left( {{t_f}} \right) = \frac{\partial }{{\partial \mathit{\boldsymbol{x}}\left( {{t_f}} \right)}}\left( {\frac{1}{2}\mathit{\boldsymbol{x}}{{\left( {{t_f}} \right)}^{\rm{T}}}\mathit{\boldsymbol{Fx}}\left( {{t_f}} \right)} \right) = \mathit{\boldsymbol{Fx}}\left( {{t_f}} \right) $$ (38) 式中，F为半正定对称常数的终端加权矩阵。
为了求解上述方程，需要知道λ(t)，由于LQR问题处理的是线性问题，所以可以有如下假设：
$$ \mathit{\boldsymbol{\lambda }}\left( t \right) = \mathit{\boldsymbol{P}}\left( t \right)\mathit{\boldsymbol{x}}\left( t \right) $$ (39) 可以得到：
$$ \mathit{\boldsymbol{\dot x}}\left( t \right) = \mathit{\boldsymbol{A}}\left( t \right)\mathit{\boldsymbol{x}}\left( t \right) + \mathit{\boldsymbol{B}}\left( t \right){\mathit{\boldsymbol{R}}^{ - 1}}\left( t \right){\mathit{\boldsymbol{B}}^{\rm{T}}}\left( t \right)\mathit{\boldsymbol{P}}\left( t \right)\mathit{\boldsymbol{x}}\left( t \right) $$ (40) $$ \mathit{\boldsymbol{\dot \lambda }}\left( t \right) = \mathit{\boldsymbol{\dot P}}\left( t \right)\mathit{\boldsymbol{x}}\left( t \right) + \mathit{\boldsymbol{P}}\left( t \right)\mathit{\boldsymbol{\dot x}}\left( t \right) = - \mathit{\boldsymbol{Q}}\left( t \right)\mathit{\boldsymbol{x}}\left( t \right) - \mathit{\boldsymbol{A}}{\left( t \right)^{\rm{T}}}\mathit{\boldsymbol{P}}\left( t \right)\mathit{\boldsymbol{x}}\left( t \right) $$ (41) 则可以得到如下最优反馈控制：
$$ \mathit{\boldsymbol{u}}\left( t \right) = - {\mathit{\boldsymbol{R}}^{ - 1}}\left( t \right){\mathit{\boldsymbol{B}}^{\rm{T}}}\left( t \right)\mathit{\boldsymbol{P}}\left( t \right)\mathit{\boldsymbol{x}}\left( t \right) = - \mathit{\boldsymbol{K}}\left( t \right)\mathit{\boldsymbol{x}}\left( t \right) $$ (42) 式中：K(t)为反馈增益。
最优反馈结构图如图 2所示。
图 2 最优状态反馈控制系统的结构图
Figure 2. Structure of optimal feedback control system
黎卡提方程是一阶非线性微分方程组，一般用计算机去求数值解。P(t)可由黎卡提方程求解得到，根据有限时间调节控制规律，可以得到无限时间的控制规律如下：
黎卡提方程：
$$ \mathit{\boldsymbol{PA + }}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{P}} - \mathit{\boldsymbol{PB}}{\mathit{\boldsymbol{R}}^{ - 1}}{\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{P}} + \mathit{\boldsymbol{Q = }}0 $$ (43) 最优控制：
$$ {\mathit{\boldsymbol{u}}^ * } = - {\mathit{\boldsymbol{R}}^{ - 1}}{\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{Px}}\left( t \right) = - \mathit{\boldsymbol{Kx}}\left( t \right) $$ (44) 最优性能指标值：
$$ {\mathit{\boldsymbol{J}}^ * } = 0.5\mathit{\boldsymbol{x}}_0^{\rm{T}}\mathit{\boldsymbol{P}}{\mathit{\boldsymbol{x}}_0} $$ (45) -
建立如图 3所示的坐标系，包括地面坐标系o-xyz和船体坐标系o'-x'y'z。地面坐标系固定不动，船体坐标系随船一起运动，六自由度的运动为u, v, w, p, q, r，正方向如图中箭头所示。设初始状态为oo'重合，ox, oy, oz轴分别与ox', oy', oz'轴重合，地面坐标系中AUV在ox方向的受力为Fd，即航行阻力；船体坐标系中AUV在o'x'方向的受力为Fx′，在o'y'方向的受力为Fy′，在o'z'方向的受力为Fz′。假设AUV以航速V运动，方向与ox方向一致，且始终保持航速不变，并在x'o'z'平面内做小幅的纵摇运动。则o'x'与ox不再一致，方向上相差纵摇角θ，AUV的航速V与纵荡速度u不再相等，方向上相差纵摇角θ，纵摇运动为小幅运动，因此数值上可认为u≈V。由于是小幅的纵摇运动，依据泰勒分解，可认为，sin θ≈θ，cos θ≈1，则航行阻力
图 3 地面坐标系和船体坐标系
Figure 3. The earth-fixed coordinate and the ship-fixed coordinate
$$ {F_{\rm{d}}} = {{F'}_x} \cdot \cos \theta - {{F'}_y} \cdot \sin \theta $$ 即
$$ {F_{\rm{d}}} = {{F'}_x} - {{F'}_y} \cdot \theta $$ (46) 纵摇运动不仅会产生纵摇角，同时会对Fx′, Fz′产生影响。AUV在水中运动时，其受到的水动力可以由势流理论分析。在无限宽广的理想流体中，由于AUV的纵摇运动产生的阻力增加可以认为是由潜体辐射势产生的，则AUV所受到的水动力(即船体坐标系中的Fx′, Fy′, Fz′)为
$$ \mathit{\boldsymbol{F = }}\iint\limits_{{s_b}\left( t \right)} {\mathit{\boldsymbol{pn}}{\rm{d}}s} $$ (47) 式中：sb(t)为物体表面；n为表面的单位法线向量，指向物体内部；p为动压力，由伯努利(Bernouli)方程确定，其表达式如下：
$$ \mathit{\boldsymbol{p}} = - \rho \left( {\frac{{\partial \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}}}{{\partial t}} + \frac{1}{2}\nabla \mathit{\boldsymbol{ \boldsymbol{\varPhi} }} \cdot \nabla \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}} \right) $$ (48) 式中，Φ为潜体运动时的辐射速度势。将式(48) 代入式(47) 中，依据推理，可得
$$ \mathit{\boldsymbol{F}} = - \rho \frac{{\rm{d}}}{{{\rm{d}}t}}\iint\limits_{{s_b}\left( t \right)} {\mathit{\boldsymbol{ \boldsymbol{\varPhi} n}}{\rm{d}}s} $$ (49) 考虑AUV做六自由度非定常运动的一般情况，用u(t)表示平动速度，用ω(t)表示绕随体运动某中心点的旋转角速度，则速度势在物面上应适合下面的边界条件，即
$$ \frac{{\partial \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}}}{{\partial \mathit{\boldsymbol{n}}}} = \mathit{\boldsymbol{u}} \cdot \mathit{\boldsymbol{n + }}\left( {\mathit{\boldsymbol{\omega }} \times \mathit{\boldsymbol{r'}}} \right) \cdot \mathit{\boldsymbol{n}} $$ (50) 式中：r'是从旋转中心出发的位置向量；记u和ω在船体坐标系中的分量为u=(u1, u2, u3)和ω=(u4, u5, u6)。
对应于船体坐标系中的物理量，有u1=u, u2=v, u3=w, u4=p, u5=q, u6=r。将n和r'×n分量表示为：n=(n1, n2, n3)，r'×n=(n4, n5, n6)
则边界条件式(50) 可改写为
$$ \frac{{\partial \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}}}{{\partial \mathit{\boldsymbol{n}}}} = \sum\limits_{i = 1}^6 {{u_i}{n_i}} $$ (51) 可以得到
$$ \mathit{\boldsymbol{F}} = - \rho \sum\limits_{i = 1}^6 {\frac{{\rm{d}}}{{{\rm{d}}x}}\left[ {{u_i}\left( t \right)\iint\limits_{{s_b}\left( t \right)} {{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_i}\left( {{{x'}_1},{{x'}_2},{{x'}_3}} \right)\mathit{\boldsymbol{n}}}} \right]{\rm{d}}s} $$ (52) 式(52) 可改写为
$$ \mathit{\boldsymbol{F}} = - \rho \sum\limits_{i = 1}^6 {{{\dot u}_i}\left( t \right)\iint\limits_{{s_b}} {\mathit{\boldsymbol{n}}{\rm{d}}s}} - \rho \sum\limits_{i = 1}^6 {{u_i}\left( t \right)\mathit{\boldsymbol{\omega }} \times \iint\limits_{{s_b}} {{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_i}\mathit{\boldsymbol{n}}{\rm{d}}s}} $$ (53) AUV做纵荡运动和纵摇运动，只考虑u1, u6(ω3)变量，分离出o′x′和o′z′方向受力，可得
$$ \begin{array}{l} {{F'}_x} = - {m_{21}}{u_1}{u_6} - {m_{26}}u_6^2\\ {{F'}_z} = - {m_{11}}{u_1}{u_6} - {m_{16}}u_6^2 \end{array} $$ (54) 式中：m21，m26，m11，m16为水动力系数。
考虑到黏性阻尼力，Fx′与u12成反比，Fz′与纵摇角成正比，则Fx′和Fz′修正为
$$ \begin{array}{l} {{F''}_x} = - {m_{21}}{u_1}{u_6} - {m_{26}}u_6^2 + {m_{77}}u_1^2\\ {{F''}_z} = - {m_{11}}{u_1}{u_6} - {m_{16}}u_6^2 + {m_{88}}\theta \end{array} $$ (55) 式中：m77为u12的纵向力系数；m88为纵摇角θ的垂向力系数。
将式(55) 代入式(46)，负号表示阻力方向与ox正向相反，去掉式中的负号，可得
$$ {F_{\rm{d}}} = Au_6^2 - B\theta {u_6} + C{\theta ^2} + Du_1^2 $$ (56) 式中，u1=u≈V, u6=q，故式(56) 可改写为
$$ {F_{\rm{d}}} = A{q^2} + B\theta q + C{\theta ^2} + D{V^2} $$ (57) 上式即AUV的航行阻力模型，A, B, C, D为待定系数，取某个纵摇周期下的AUV航行阻力的CFD仿真结果，利用线性回归的方式便可求得待定系数。
AUV在水下航行时所带的能量有限，为了提高AUV的续航能力需要对AUV进行节能控制。AUV在近水面航行时由于受到海浪的干扰使姿态发生变化，由于垂直面姿态的变化导致航行阻力增加，最终导致能量消耗增加。为了维持AUV的垂直面姿态，需使AUV的前后舵偏转一定的角度，前后舵的转动也会造成阻力增加。为使AUV航行阻力最小，综合考虑这两个因素得到最优控制。
前文得到了纵摇增阻与纵摇的关系。则能量节省的指标函数可以转换成如下方程：
$$ J = \frac{1}{t}\int_0^t {\left( {{\lambda _1}{q^2} + {\lambda _2}q\theta + {\lambda _3}{\theta ^2} + {\lambda _4}u_1^2 + {\lambda _5}u_2^2} \right){\rm{d}}t} $$ (58) 式中，λi(i=1, 2, …, 5) 是使能量消耗指标最小的值。对上式做离散化处理，考虑本文中的实际需求，航行阻力及能量消耗的性能评价式可表示为
$$ J = \frac{1}{N}\sum\limits_{i = 1}^N {\left( {{\lambda _1}{q^2} + {\lambda _2}q\theta + {\lambda _3}{\theta ^2} + {\lambda _4}u_1^2 + {\lambda _5}u_2^2} \right)} $$ (59) -
首先，建立能量消耗的性能指标；然后，以此性能指标设计LQR控制器，加入海浪的干扰；最后，搭建Simulink仿真模型，对LQR控制器的控制效果进行仿真分析和检验。仿真条件如下：海况为3级，有义波高为0.88 m、平均周期为6.43 s，Cm=1.95。对AUV在水深d=3 m处以4.5 m/s的速度航行时受到遭遇角βw为45°的海浪垂荡力和纵摇力矩进行仿真，如图 4和图 5所示。干扰力和干扰力矩对比未加控制、LQR控制下AUV垂直面运动垂荡位移和纵摇角度响应的仿真曲线如图 6和图 7所示。
图 4 βw=45°, d=3 m时垂荡力的仿真结果
Figure 4. Result of heave force when βw=45°, d=3 m
图 5 βw=45°，d=3 m时纵摇力矩的仿真结果
Figure 5. Result of pitch moment when βw=45°, d=3 m
图 6 AUV垂荡位移变化曲线
Figure 6. Time histories of the heave displacement of AUV
图 7 AUV纵摇变化曲线
Figure 7. Time histories of the pitch angle of AUV
由图 6和图 7的仿真曲线可以看出，设计的AUV垂直面运动LQR控制器能够有效抑制AUV的垂直面姿态运动，具有良好的控制效果。加入LQR控制后，垂荡位移由原来的最大幅值0.981 2 m降低至0.619 6 m，降幅为36.85%。纵摇角度由原来的最大值6.83°降低至3.02°，降幅为55.77%。所设计的LQR控制器能够有效抵消海浪干扰，纵摇角和垂荡位移的最大幅值均明显降低。因此，AUV垂直面运动加入LQR控制器后运动变得平稳，航行阻力变小，降低了能量消耗。
对AUV的艏艉舵角度幅值进行15°的限位，艏艉舵转动的角度如图 8和图 9所示。
图 8 LQR控制中艏舵的摆角
Figure 8. The deflection angle of the fore control surface with LQR control
图 9 LQR控制中艉舵的摆角
Figure 9. The deflection angle of the aft control surface with LQR control
本文使用方差对控制效果进行评价，定义减纵摇效果公式如下：
$$ E = \left( {{V_{\rm{a}}} - {V_{\rm{b}}}} \right)/V_{\rm{a}}^ * 100\% $$ (60) 式中：E为减摇效果；Va为未减摇纵摇角/垂荡位移方差；Vb为减摇后纵摇角/垂荡位移方差。
用同样的方法定义减垂荡效果标准。对控制结果进行统计分析，分析结果如表 2所示。
表 2 姿态控制效果统计
Table 2. Effect statistics of attitude control
未加控制 LQR控制 垂荡最大值/m 0.981 2 0.709 7 纵摇角最大值(/°) 6.831 3 3.301 0 垂荡位移方差/m2 0.148 6 0.079 3 纵摇角方差(/°)2 6.858 5 1.534 8 减垂荡效果/% - 46.64 减纵摇效果/% - 77.62 由表 2可以看出，AUV的垂直面姿态控制效果良好，加入LQR控制后有效减小了垂荡位移和纵摇角。由定义的姿态控制评判标准(式(60))可知，LQR控制使AUV的减垂荡和减纵摇效果分别达到46.64%和77.62%。
图 10 未加控制阻力增值曲线
Figure 10. Added resistance without control
图 11 控制后阻力增值曲线
Figure 11. Added resistance with LQR control
由图 10、图 11和表 3可知，加入LQR控制后，有效减小了AUV的纵摇增阻，降低了AUV的推进器的能量消耗，实现了AUV姿态和能量的综合最优控制。
表 3 增阻控制效果统计
Table 3. Effect statistics of added resistance control
未加控制 LQR控制 阻力增值平均值/N 437.039 2 81.153 5 阻力增值标准差/N 488.747 5 82.154 8 -
本文建立了航行阻力与航行姿态之间的数学模型，将该阻力模型定义为线性二次型控制的性能指标，并以此为基础确定了LQR控制器的权矩阵，设计了LQR控制器。LQR控制方法的减垂荡和减纵摇的效果分别达到46.64%和77.62%，同时将航行阻力增值降为原来的1/6，取得了满意的控制效果。
在控制近水面AUV航行姿态时，可同时对减摇和降低航行阻力进行优化，以减小AUV的能量消耗。本文提出了将AUV航行阻力模型作为LQR控制器性能指标的方法，取得了较为满意的效果，可以为AUV减小能量消耗，降低航行阻力提供一种新的控制思路。
全文HTML
3.1. AUV动力学模型的线性化
3.2. 控制器
2017-3-111_en.pdf
|
|