-
摘要:目的
针对存在船舶模型参数未知、外界扰动未知以及舵机约束等问题，提出一种基于数据驱动的在线辨识船舶参数，迭代解析计算最优控制量的航迹跟踪控制方法。
方法构建双螺旋桨船的三自由度动力学方程，通过采集船舶的运动数据，设计扩张状态观测器−多新息递推最小二乘交互式算法。再将辨识得到的船舶运动模型在采样周期内近似为定常线性模型，将船舶航迹跟踪问题转变成带约束和干扰的线性二次型优化控制问题。通过引入加权矩阵与罚函数，构建包含轨迹误差、外界干扰量和控制量约束不等式的二次型性能指标，并运用精细积分法获得矩阵黎卡提微分方程的解析解，得到有限时间状态调节器的迭代计算式。
结果实现了在线辨识船舶运动模型参数和估计未知扰动，并设计出一种“启动后不用管”的航迹跟踪控制算法，降低了参数辨识和控制算法对实验设计的严格要求。
结论通过Matlab数值仿真分析权重矩阵 \boldsymbol{Q}，\boldsymbol{R} 和S对航迹跟踪精度的影响，验证了参数辨识和控制算法的有效性。
Abstract:ObjectiveAiming at the problems of unknown ship model parameters and external disturbance and servo constraints, this paper proposes a method for the data-driven online identification of ship parameters and iterative analytical calculation of the optimal control quantity of track tracking control.
MethodA three degrees of freedom dynamics equation of a double propeller ship is constructed, and the extended state observer-multiple innovation recursive least squares interactive algorithm is designed by collecting the motion data of the ship. By approximating the identified ship motion model to a time-invariant linear model in the sampling period, the ship trajectory tracking problem can be transformed into a linear quadratic optimization control problem with constraints and disturbances. The weighted matrix and penalty function are introduced to construct the quadratic performance index including trajectory error, external disturbance, and control constraint inequality. The precise integration method is then used to obtain the analytical solution of the matrix Riccati differential equation and the iterative calculation formula of the finite time state regulator.
ResultsTh online identification of the ship motion model parameters and estimation of unknown disturbances are achieved, and a trajectory tracking control algorithm with "no need to worry after startup" is designed, reducing the strict requirements of parameter identification and control algorithms for experimental design.
ConclusionUsing MATLAB to carry out numerical simulation and analyze the influence of weight matrix \boldsymbol{Q},\boldsymbol{R} and S on trajectory tracking accuracy, the results verify the effectiveness of the parameter identification and control algorithm.
-
0. 引 言
航迹跟踪控制是船舶运动控制中一个较为典型的问题。当前的研究从不同角度出发，提出多种解决方法，推动了船舶智能航行理论的发展及实际应用技术的进步。然而，从实际应用的角度来看，部分现有的航迹跟踪控制方法仍存在以下问题：首先，这些方法以船舶运动数学模型为基础设计控制算法，然而实际船舶通常无法提供船舶运动模型，即使能够提供，也无法保证模型的参数和结构在运行过程中保持不变；其次，外部扰动的随机性以及执行器的非线性延迟特性未得到客观地表达，通常的处理方式是人为假定扰动的随机分布形式，并假设其在整个控制过程保持不变，同时，执行器特性建模未充分考虑其工作机理；最后，所推导的控制算法大多基于理想条件得出，导致算法中含有大量经验参数，而这些参数在实际应用中缺乏标准化的确定流程或逻辑计算方法，通常是针对特定船型和特定环境在仿真实验中反复尝试获得，难以在实船上直接应用，即控制参数缺乏自适应性。为解决上述问题，本文拟设计一种“下水后不用管”的采集、辨识、控制一体化自适应方法体系，仅依赖船舶运动的实时状态数据驱动，无需人为假设或难以获取的先验信息，从而推动船舶运动控制算法向实用化方向发展。
在船舶运动控制中，实现航迹跟踪的常用控制器方法包括鲁棒控制器、滑模控制器、比例−积分−微分（PID）控制器、模型预测控制器（MPC）、强化学习控制和线性二次调节器（LQR）等。Weng等[1]提出基于数据驱动级联控制的复合鲁棒跟踪控制方法，克服了由不确定性、扰动和未建模动力学引起的集总扰动。郭杰等[2]设计一种基于多模态快速非奇异终端滑模的自抗扰控制器，通过设计基于多模态的自抗扰控制器（ADRC）控制律，将航迹控制问题转化为易于实现的航向控制问题。宋利飞等[3]引入深度强化学习理论，设计智能体环境状态、动作和奖励函数以在线调整PID参数，增强智能体训练后期的稳态性能，提高了路径跟踪精度。张子昌等[4]针对无人水下航行器（AUV）的水平面和垂直面控制，设计了基于状态空间形式的模型预测控制算法。Wang等[5]提出一种数据驱动的性能规定强化学习控制方案，可以避免模型不准确、扰动难以估计的问题，同时追求最优控制和规定的跟踪精度。Nguyen等[6]提出一种船舶运动的最优控制器，建立了基于线性二次正则算法的船舶运动控制系统。
然而，这些控制器都是以船舶模型参数已知且不变为前提。在实际工程中，船舶参数大多未知，且由于船舶长时间运行，舵机性能、船体阻力、操纵性参数等会发生变化，水流、船舶吃水的变化也会实时改变船舶操纵性参数[7]。因此有必要实现船舶操纵模型参数的在线辨识。
目前常用的系统辨识方法包括最小二乘法[8]、卡尔曼滤波[9]、支持向量机[10]、人工神经网络[11]和多新息系统辨识[12]。这些辨识方法在限定场景和船型下均取得了良好效果，但在实际场景中，由于船舶所受环境扰动的不确定性，会严重影响辨识参数的精度。
本文将从数据出发，完全基于船舶的运动数据，如纵向速度、横向速度、艏摇角速度、艏摇角、螺旋桨转速等，在不确定环境干扰下，设计一种结合扩张状态观测器和多新息递推最小二乘辨识理论的在线交互式辨识算法，同时进行扰动估计和参数辨识。该算法无需完全遵循固定的航迹要求，大大降低辨识算法对实验设计的严格要求。从控制角度出发，将辨识后的模型转变成线性二次型，提出含有约束不等式和干扰的目标函数，采用极小值原理和连续动态规划设计含有约束和干扰的有限时间状态调节器。
1. 船舶运动模型及问题描述
1.1 船舶运动模型的建立
船舶在海面上的运动状态如图1所示，主要包括横摇、纵荡和横荡以及艏摇，一般忽略船舶的横摇（横倾角速度p = 0）和纵摇（纵倾角速度q = 0）。
因此，根据《船舶操纵运动数学模型》[13]可描述如下：
\left\{\begin{aligned} & \dot{x}=u\mathrm{cos}\psi -\upsilon \mathrm{sin}\psi \\& \dot{y}=u\mathrm{sin}\psi +\upsilon \mathrm{cos}\psi \\& \dot{\psi }=r \end{aligned}\right. (1) \left\{\begin{aligned} &\left(m-{X}_{\dot{u}}\right)\dot{u}-\left(m-{Y}_{\dot{v}}\right)vr={X}_{\mathrm{H}}+{X}_{\mathrm{R}}+{X}_{\mathrm{P}}\\& \left(m-{Y}_{\dot{v}}\right)\dot{v}-\left(m-{X}_{\dot{u}}\right)ur={Y}_{\mathrm{H}}+{Y}_{\mathrm{R}}+{Y}_{\mathrm{P}}\\& \left({I}_{Z}-{N}_{\dot{r}}\right)\dot{r}={N}_{\mathrm{H}}+{N}_{\mathrm{R}}+{N}_{\mathrm{P}}\end{aligned}\right. (2) 其中，
\left\{\begin{aligned} &{X}_{\mathrm{H}}={X}_{vv}{v}^{2}+{X}_{rr}{r}^{2}+{X}_{vr}vr+{X}_{\mathrm{H}0}\left({F}_{r}\right)\\& {Y}_{\mathrm{H}}={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}_{v\left|r\right|}v\left|r\right|\\& {N}_{\mathrm{H}}= {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}_{\left|v\right|r}\left|v\right|r\end{aligned}\right. (3) 式中： x，y，\psi 分别为船体坐标系下相对风浪流的纵向 、横向速度及艏向角速度； {X}_{\mathrm{H}}，{Y}_{\mathrm{H}}，{N}_{\mathrm{H}} 分别为船体产生的纵向力、侧向力、偏航力矩； {X}_{\mathrm{R}}， {Y}_{\mathrm{R}}，{N}_{\mathrm{R}} 分别为舵产生的纵向力、侧向力、偏航力矩； {X}_{\mathrm{P}}，{Y}_{\mathrm{P}}，{N}_{\mathrm{P}} 分别为螺旋桨产生的纵向力、侧向力、偏航力矩； {X}_{vv}，{X}_{rr}，{X}_{vr}，{Y}_{v}，{Y}_{r}，{Y}_{v\left|v\right|}，{Y}_{r\left|r\right|}， {Y}_{v\left|r\right|}，{N}_{v}，{N}_{r}，{N}_{v\left|v\right|}，{N}_{r\left|r\right|}，{N}_{\left|v\right|r} 为流体动力导数，在模型中为常数； {X}_{\mathrm{H}0}\left({F}_{r}\right) 为直航阻力，可以近似[14]为
{X}_{\mathrm{H}0}\left({F}_{r}\right)={X}_{uu}{u}^{2} (4) 一般认为螺旋桨推力在y轴向分量 {Y}_{\mathrm{P}}=0 ，且本研究采用双桨推进无人艇进行船舶仿真实验，不存在舵，无人艇只依靠左右螺旋桨前进和转向，因此 {X}_{\mathrm{R}}=0,{Y}_{\mathrm{R}}=0,{N}_{\mathrm{R}}=0 。
在给出正车情况下，螺旋桨力及力矩的计算模型和舵机特性模型[15]有
\left\{\begin{aligned} &{X}_{\mathrm{P}}=\left(1-{t}_{\mathrm{P}}\right)\left({T}_{\mathrm{P}}+{T}_{\mathrm{S}}\right)\\& {Y}_{\mathrm{P}}=0\\& {N}_{\mathrm{P}}=\left(1+{d}_{\mathrm{N}\mathrm{P}}\right)\left({T}_{\mathrm{P}}-{T}_{\mathrm{S}}\right){\boldsymbol{y}}_{\mathrm{P}}\end{aligned}\right. (5) T=\rho {n}^{2}{D}_{\mathrm{P}}^{4}{k}_{\mathrm{T}}\left({J}_{\mathrm{P}}\right) (6) {k}_{\mathrm{T}}\left({J}_{\mathrm{P}}\right)={a}_{0}+{a}_{1}{J}_{\mathrm{P}}+{a}_{2}{J}_{\mathrm{P}}^{2} (7) {J}_{\mathrm{P}}=\frac{\left(1-{w}_{\mathrm{P}}\right)u}{n{D}_{\mathrm{P}}} (8) 式中： {t}_{\mathrm{P}} 为螺旋桨的推力减额系数； {T}_{\mathrm{S}} 和 {T}_{\mathrm{P}} 分别为左、右螺旋桨的推力； {d}_{\mathrm{N}\mathrm{P}} 为螺旋桨对转向力矩的影响系数； {\boldsymbol{y}}_{\mathrm{P}} 为1/2倍螺旋桨间距； \rho 为水的密度； {D}_{\mathrm{P}} 为螺旋桨的直径； {k}_{\mathrm{T}}\left({J}_{\mathrm{P}}\right) 为推力系数； {J}_{\mathrm{P}} 为进速系数； {w}_{\mathrm{P}} 为螺旋桨的伴流系数；n为转速。
对式（5）～式（8）进一步简化可得螺旋桨水动力模型：
\left\{\begin{aligned} &{X}_{\mathrm{P}}={\check{X}}_{uu}{u}^{2}+{\check{X}}_{u\mathrm{P}}u{n}_{\mathrm{P}}+{\check{X}}_{u\mathrm{S}}{un}_{\mathrm{S}}+{\check{X}}_{\mathrm{P}\mathrm{P}}{n}_{\mathrm{P}}^{2}+{{\check{X}}_{\mathrm{S}\mathrm{S}}n}_{\mathrm{S}}^{2}\\& {Y}_{\mathrm{P}}=0\\& {N}_{\mathrm{P}}={N}_{uu}{u}^{2}+{\check{N}}_{u\mathrm{P}}u{n}_{\mathrm{P}}+{\check{N}}_{u\mathrm{S}}{un}_{\mathrm{S}}+{\check{N}}_{\mathrm{P}\mathrm{P}}{n}_{\mathrm{P}}^{2}+{{\check{N}}_{\mathrm{S}\mathrm{S}}n}_{\mathrm{S}}^{2}\end{aligned} \right. (9) 式中：{\check{X}}_{uu} ， {\check{X}}_{u\mathrm{P}} ， {\check{X}}_{u\mathrm{S}} ， {\check{X}}_{\mathrm{P}\mathrm{P}} ， {\check{X}}_{\mathrm{S}\mathrm{S}} ， {\check{N}}_{uu} ， {\check{N}}_{u\mathrm{P}} ， {\check{N}}_{u\mathrm{S}} ， {\check{N}}_{\mathrm{P}\mathrm{P}} ， {\check{N}}_{\mathrm{S}\mathrm{S}} 为简化后的综合系数； {n}_{\mathrm{S}}\mathrm{和}{n}_{\mathrm{P}} 分别为左、右螺旋桨转速。
1.2 参数辨识问题
考虑到环境干扰以及测量噪声干扰等因素的存在，通过将状态量与辨识参数分开，可将式（2）中每一维转换为关于辨识参数的线性模型
{\boldsymbol{y}}={\bf\textit{φ}}^{\mathrm{T}}{\boldsymbol{\theta}} +{\boldsymbol{w}} (10) 式中：y为每一维的输出量，即 \dot{u} ， \dot{v} 或 \dot{r} ； \boldsymbol{\theta}\in\boldsymbol{\mathbf{\mathrm{\mathbf{R}}}}^{l\times1} ，为包含船舶水动力参数及螺旋桨参数等模型参数的向量，其中l为参数的维数；{\bf\textit{φ}}\in {\mathbf{R}}^{l\times 1} ，为信息向量，是由 u，v，r，{n}_{\mathrm{S}}，{n}_{\mathrm{P}} 构成的函数；w为船舶受环境影响导致的扰动量，包含了风浪流等因素。
若有t组测量数据，且数据量足够多，每一组数据都可得到输出量 {\boldsymbol{y}}\left(t\right) 、信息量{\bf\textit{φ}}\left(t\right) 和扰动量 {\boldsymbol{w}}\left(t\right) 。将测量数据代入式（10）可得矩阵方程
{\boldsymbol{Y}}_{t}={\boldsymbol{X}}_{t}{\boldsymbol{\theta}} +{\boldsymbol{w}}_{t} (11) 式中 ：{\boldsymbol{Y}}_{t}\in {\mathbf{R}}^{t\times 1} ， {\boldsymbol{X}}_{t}\in {\mathbf{R}}^{t\times l} ， {\boldsymbol{w}}_{t}\in {\mathbf{R}}^{t\times 1} 。若 \hat{\boldsymbol{\theta }} 为 \boldsymbol{\theta } 的参数辨识估计值， \hat{y}\left(t\right) 为 \boldsymbol{y}\left(t\right) 的拟合值，且满足 \hat{y}\left(t\right)={\bf\textit{φ}}^{\mathrm{T}}\left(t\right)\hat{\boldsymbol{\theta }} ，则拟合值向量 {\hat{\boldsymbol{Y}}}_{t} 有
{\hat{\boldsymbol{Y}}}_{t}=\left[\begin{matrix}\hat{y}\left(1\right)\\ \hat{y}\left(2\right)\\ \vdots \\ \hat{y}\left(t\right)\end{matrix}\right]={\boldsymbol{X}}_{t}\hat{\boldsymbol{\theta }}=\left[\begin{matrix}{\bf\textit{φ}}^{\mathrm{T}}\left(1\right)\hat{\boldsymbol{\theta }}\\ {\bf\textit{φ}}^{\mathrm{T}}\left(2\right)\hat{\boldsymbol{\theta }}\\ \vdots \\ {\bf\textit{φ}}^{\mathrm{T}}\left(t\right)\hat{\boldsymbol{\theta }}\end{matrix}\right] (12) 参数识别是在使用输入和输出数据建立数学模型的基础上，找到符合方程（12）的最佳参数 \hat{\boldsymbol{\theta }} 的过程。最小二乘法是一种常用的识别方法，其利用方差最小化原理，通过最小化剩余的标准函数来求解参数 \hat{\boldsymbol{\theta }} 。对于辨识模型，取准则函数可以得到
\begin{split} & J(\hat{\boldsymbol{\theta }})={\|{\boldsymbol{Y}}_{t}-{\hat{\boldsymbol{Y}}}_{t}\|}^{2}=\sum _{i=1}^{t}{\|[y(i)-{\bf\textit{φ}}^{\mathrm{T}}(i)\hat{\boldsymbol{\theta }}]\|}^{2}=\\& \qquad\qquad {({\boldsymbol{Y}}_{t}-{\boldsymbol{X}}_{t}\hat{\boldsymbol{\theta }})}^{\mathrm{T}}({\boldsymbol{Y}}_{t}-{\boldsymbol{X}}_{t}\boldsymbol{\theta }) \end{split} (13) 为了最小化J\left(\hat{\boldsymbol{\theta }}\right) ，设J的偏导数为0，可以得到
\frac{\partial J\left(\hat{\boldsymbol{\theta }}\right)}{\partial \hat{\boldsymbol{\theta }}}=-{\boldsymbol{X}}_{t}^{\mathrm{T}}\left({\boldsymbol{Y}}_{t}-{\boldsymbol{X}}_{t}\hat{\boldsymbol{\theta }}\right)=0 (14) 通过式（14），可以得到
\hat{\boldsymbol{\theta }}={\left({\boldsymbol{X}}_{t}^{\mathrm{T}}{\boldsymbol{X}}_{t}\right)}^{-1}{\boldsymbol{X}}_{t}^{\mathrm{T}}{\boldsymbol{Y}}_{t} (15) 将式（11）代入式（15），并对式（15）的两边进行数学期望，得到
E\left\{\hat{\boldsymbol{\theta }}\right\}=E\left\{\boldsymbol{\theta }\right\}+E\left\{{\left({\boldsymbol{X}}_{t}^{\mathrm{T}}{\boldsymbol{X}}_{t}\right)}^{-1}{\boldsymbol{X}}_{\mathrm{t}}^{\mathrm{T}}\right\}E\left\{{\boldsymbol{w}}_{t}\right\} (16) 式（16）给出了扰动 {\boldsymbol{w}}_{t} 对参数辨识的影响。因此，如何处理扰动影响是决定辨识算法鲁棒性的关键。
1.3 航迹跟踪问题
航迹跟踪是在给定期望航迹为大地坐标系中由一系列航路点组成的连续曲线，船舶从任意出发点到达期望航迹并以给定期望速度以尽量小的跟踪误差沿期望航迹航行[16]。为便于表示航迹跟踪误差，现以坐标原点为初始时刻航迹参考点 {P}_{0} ，建立平面直角坐标系，如图2所示。
将期望航迹曲线表示为参数形式，定义航迹上参考点坐标 {P}_{k}\left({x}_{\mathrm{d}}\right(t)，{\boldsymbol{y}}_{\mathrm{d}}(t\left)\right) ，其中k表示航迹上第k个点，U为参考点速度，若 {\alpha }_{k} 为参数路径上任一点处和上一路径点连线与y轴的夹角，规定顺时针为正，表达式为
{\alpha }_{k}=\arctan\frac{{{x}_{\mathrm{d}}\left(t\right)}'}{{{\boldsymbol{y}}_{\mathrm{d}}\left(t\right)}'},\;\;\; {\alpha }_{k}\in \left(-{\text{π}},{\text{π}}\right) (17) 在直角坐标系中，船舶当前位置为 P({x}_{0},{\boldsymbol{y}}_{0}) ，根据几何关系可得航迹纵向跟踪误差为
e\left(t\right)=-\mathrm{sin}{\alpha }_{k}(x-{x}_{\mathrm{d}}\left(t\right))+\mathrm{cos}{\alpha }_{k}\left(y-{\boldsymbol{y}}_{\mathrm{d}}\left(t\right)\right) (18) 若期望实现船舶的航迹跟踪，其最终控制目标公式表示为
\underset{t\to \mathrm{\infty }}{\mathrm{lim}}e\left(t\right)=0 (19) 2. 模型参数辨识
2.1 在线辨识模型
以双螺旋桨船为实验对象，通过系统输入一定变化的左右螺旋桨转速 {n}_{\mathrm{S}}，{n}_{\mathrm{P}} 来获取一系列船舶的纵荡速度u、横荡速度v、艏摇角速度r和艏摇角 \psi 等数据。若设船舶采样周期为h，t为测量次数，通过使用前向差分法可近似加速度量 \dot{u} ， \dot{v} ， \dot{r} ，考虑存在外界干扰的影响，设存在对加速度量的扰动分别为 {w}_{u}\left(t\right)，{w}_{v}\left(t\right)，{w}_{r}\left(t\right) ，则可以将式（2）的船舶动力学模型重新写成
\left\{\begin{aligned} &\frac{{{u}}\left(t+1\right)-{{u}}\left(t\right)}{h}={\mathcal{R}}_{u}\left(t\right){\boldsymbol{\theta }}_{u}+{w}_{u}\left(t\right)\\& \frac{v\left(t+1\right)-v\left(t\right)}{h}={\mathcal{R}}_{v}\left(t\right){\boldsymbol{\theta }}_{v}+{w}_{v}\left(t\right)\\& \frac{r\left(t+1\right)-r\left(t\right)}{h}={\mathcal{R}}_{r}{\left(t\right)\boldsymbol{\theta }}_{r}+{w}_{r}\left(t\right)\end{aligned}\right. (20) \left\{\begin{aligned} &{\mathcal{R}}_{u}=\left[{vr,u}^{2},{v}^{2},{{r}^{2},u{n}_{\mathrm{P}},{un}_{\mathrm{S}},{n}_{\mathrm{P}}^{2},n}_{\mathrm{S}}^{2}\right]\\& {\mathcal{R}}_{v}=\left[ur,v,r,v\left|v\right|,r\left|r\right|,v\left|r\right|\right]\\& {\mathcal{R}}_{r}=\left[v,r,v\left|v\right|,r\left|r\right|,r\left|v\right|,{u}^{2}, {u{n}_{\mathrm{P}},{un}_{\mathrm{S}},{n}_{\mathrm{P}}^{2},n}_{\mathrm{S}}^{2}\right]\end{aligned}\right. (21) \left\{\begin{aligned} &{\boldsymbol{\theta }}_{u}={\left[\frac{\left(m-{Y}_{\dot{v}}\right)+{X}_{vr}}{m-{X}_{\dot{u}}},\frac{{{X}_{uu}+\check{X}}_{uu}}{m-{X}_{\dot{u}}},\frac{{X}_{vv}}{m-{X}_{\dot{u}}},\frac{{X}_{rr}}{m-{X}_{\dot{u}}},\frac{{\check{X}}_{up}}{m-{X}_{\dot{u}}},\frac{{\check{X}}_{us}}{m-{X}_{\dot{u}}},\frac{{\check{X}}_{pp}}{m-{X}_{\dot{u}}},\frac{{\check{X}}_{ss}}{m-{X}_{\dot{u}}}\right]}^{\mathrm{T}}\\& {\boldsymbol{\theta }}_{v}={\left[\frac{m-{X}_{\dot{u}}}{m-{Y}_{\dot{v}}},\frac{{Y}_{v}}{m-{Y}_{\dot{v}}},\frac{{Y}_{r}}{m-{Y}_{\dot{v}}},\frac{{Y}_{v\left|v\right|}}{m-{Y}_{\dot{v}}},\frac{{Y}_{r\left|r\right|}}{m-{Y}_{\dot{v}}},\frac{{Y}_{v\left|r\right|}}{m-{Y}_{\dot{v}}}\right]}^{\mathrm{T}}\\& {\boldsymbol{\theta }}_{r}={\left[\frac{{N}_{v}}{{I}_{{\textit{z}}{\textit{z}}}-{N}_{\dot{r}}},\frac{{N}_{r}}{{I}_{{\textit{z}}{\textit{z}}}-{N}_{\dot{r}}},\frac{{N}_{v\left|v\right|}}{{I}_{{\textit{z}}{\textit{z}}}-{N}_{\dot{r}}},\frac{{N}_{r\left|r\right|}}{{I}_{{\textit{z}}{\textit{z}}}-{N}_{\dot{r}}},\frac{{N}_{v\left|r\right|}}{{I}_{{\textit{z}}{\textit{z}}}-{N}_{\dot{r}}},\frac{{\check{N}}_{up}}{{I}_{{\textit{z}}{\textit{z}}}-{N}_{\dot{r}}},\frac{{\check{N}}_{us}}{{I}_{{\textit{z}}{\textit{z}}}-{N}_{\dot{r}}},\frac{{\check{N}}_{up}}{{I}_{{\textit{z}}{\textit{z}}}-{N}_{\dot{r}}},\frac{{\check{N}}_{pp}}{{I}_{{\textit{z}}{\textit{z}}}-{N}_{\dot{r}}},\frac{{\check{N}}_{ss}}{{I}_{{\textit{z}}{\textit{z}}}-{N}_{\dot{r}}}\right]}^{\mathrm{T}}\end{aligned}\right. (22) 式中： {\mathcal{R}}_{u}\in {\mathbf{R}}^{8\times 1} ， {\mathcal{R}}_{v}\in {\mathbf{R}}^{6\times 1} ， {\mathcal{R}}_{r}\in {\mathbf{R}}^{10\times 1} ，为信息向量； {\boldsymbol{\theta }}_{u}\in {\mathbf{R}}^{8\times 1} ， {\boldsymbol{\theta }}_{v}\in {\mathbf{R}}^{6\times 1} ， {\boldsymbol{\theta }}_{r}\in {\mathbf{R}}^{10\times 1}， 为参数向量。为后期编写方便，将式（22）用简单符号表示，则有
\left\{\begin{aligned} &{\boldsymbol{\theta }}_{u}={\left[{\theta }_{11},{\theta }_{12},{\theta }_{13},{\theta }_{14},{\theta }_{15},{\theta }_{16},{\theta }_{17},{\theta }_{18}\right]}^{\mathrm{T}}\\& {\boldsymbol{\theta }}_{v}={\left[{\theta }_{21},{\theta }_{22},{\theta }_{23},{\theta }_{24},{\theta }_{25},{\theta }_{26}\right]}^{\mathrm{T}}\\& {\boldsymbol{\theta }}_{r}={\left[{\theta }_{31},{\theta }_{32},{\theta }_{33},{\theta }_{34},{\theta }_{35}, {\theta }_{36},{\theta }_{37},{\theta }_{38},{\theta }_{39},{\theta }_{30}\right]}^{\mathrm{T}}\end{aligned}\right. (23) 2.2 扩张状态观测器
扩张状态观测器是一种非线性滤波器，通过递归地生成状态估计值和协方差估计值来实现它的目标，其主要思想是将非线性动态系统离散化为一个线性动态系统，并利用一个当前时刻系统状态估计值和相关一阶导数组成的扩张状态向量来描述该线性系统，是估计模型内部不确定性和外界扰动的有效方式。
若记 \boldsymbol{\nu }={\left[u,v,r\right]}^{\mathrm{T}} ， \boldsymbol{w}=[{w}_{u}\left(t\right),{w}_{v}\left(t\right) , {{w}_{r}\left(t\right)]}^{\mathrm{T}} ， {\boldsymbol{\mathcal{L}}}={\left[{\mathcal{R}}_{u}{\boldsymbol{\theta }}_{u},{\mathcal{R}}_{v}{\boldsymbol{\theta }}_{v},{\mathcal{R}}_{r}{\boldsymbol{\theta }}_{r}\right]}^{\mathrm{T}}， 系统引入观测变量 {\boldsymbol{y}}_{\mathrm{g}} ，则可将船舶动力学模型重新改写为
\left\{\begin{aligned} &\dot{\boldsymbol{\nu }}={\cal{L}}+{\boldsymbol{w}}\\& {\boldsymbol{y}}_{\mathrm{g}}={\boldsymbol{\nu}} \end{aligned}\right. (24) 针对式（24）设计一阶线性扩张状态观测器，则有
\left\{\begin{aligned} &{\boldsymbol{e}}_{\mathrm{o}}=\hat{\boldsymbol{\nu }}-{\boldsymbol{y}}_{\mathrm{g}}\\& \dot{\hat{\boldsymbol{\nu }}}=-{\alpha }_{1}{\boldsymbol{e}}_{\mathrm{o}}+{\cal{L}}+\hat{\boldsymbol{w}}\\& \dot{\hat{\boldsymbol{w}}}=-{\alpha }_{2}{\boldsymbol{e}}_{\mathrm{o}}\end{aligned}\right. (25) \left\{\begin{aligned} &{\alpha }_{1}=2\sigma \\& {\alpha }_{2}={\sigma }^{2}\end{aligned}\right. (26) 式中： \hat{\boldsymbol{\nu }} 和 \hat{\boldsymbol{w}} 分别为对状态量 \boldsymbol{\nu } 和扰动w的估计值； {\boldsymbol{e}}_{\mathrm{o}} 为状态估计误差； {\alpha }_{1} 和 {\alpha }_{2} 为调节参数； \sigma 为在线整定且 \sigma 取值范围很广。
设采样时间间隔h，可将式（25）离散化得到扩张状态观测器扰动估计算法，则有
\left\{\begin{aligned} &{\boldsymbol{e}}_{\mathrm{o}}\left(t\right)=\hat{\boldsymbol{\nu }}\left(t\right)-{\boldsymbol{y}}_{\mathrm{g}}\left(t\right)\\& \hat{\boldsymbol{\nu }}\left(t+1\right)=\hat{\boldsymbol{\nu }}\left(t\right)+h(-{\alpha }_{1}{\boldsymbol{e}}_{\mathrm{o}}\left(t\right)+{\cal{L}}\left(t\right)+\hat{\boldsymbol{w}}\left(t\right))\\& \hat{\boldsymbol{w}}\left(t+1\right)=\hat{\boldsymbol{w}}\left(t\right)+h\times \left(-{\alpha }_{2}{\boldsymbol{e}}_{\mathrm{o}}\left(t\right)\right)\end{aligned}\right. (27) 该算法在实际运行过程中计算 \mathcal{L}\left(t\right) 用到辨识参数 \boldsymbol{\theta } ，可以使用辨识估计值 \hat{\boldsymbol{\theta }} 替代。
2.3 多新息递推最小二乘算法
多新息辨识是基于传统时域辨识方法，使用以当前时刻为基准，具有一定时间长度的历史信息数据来进行递推辨识，以达到参数快速收敛至较高精度的目的[17]。
根据式（20）对每一维进行单新息递推，则统一形式为
{\boldsymbol{y}}_{*}\left(t\right)={\mathcal{R}}_{\mathcal{*}}\left(t\right){\boldsymbol{\theta }}_{\mathrm{*}}+{\boldsymbol{{w}}}_{\mathrm{*}}\left(t\right) (28) 式中：*代表下标 u,v,r 的集合； {\boldsymbol{y}}_{*}\left(t\right) 为模型输出矩阵； {\boldsymbol{\theta }}_{\mathrm{*}} 为模型参数矩阵； {\mathcal{R}}_{\mathcal{*}}\left(t\right) 为模型输入输出数据形成的信息矩阵； {\boldsymbol{w}}_{\mathrm{*}}\left(t\right) 为扰动值矩阵。
为实现在线参数辨识，单新息递推最小二乘辨识的递推公式[18]为
\left\{\begin{aligned} &{\boldsymbol{e}}_{\rm{c}}\left(t\right)={\boldsymbol{y}}_{*}\left(t\right)-{\mathcal{R}}_{\mathcal{*}}\left(t\right)\hat{\boldsymbol{\theta }}\left(t-1\right)\\& {\boldsymbol{P}}^{-1}\left(t\right)={\boldsymbol{P}}^{-1}\left(t-1\right)+{{\mathcal{R}}_{\mathcal{*}}}^{\mathrm{T}}\left(t\right){\mathcal{R}}_{\mathcal{*}}\left(t\right)\\& \hat{\boldsymbol{\theta }}\left(t\right)=\hat{\boldsymbol{\theta }}\left(t-1\right)+{\boldsymbol{P}}\left(t\right){\mathcal{R}}_{\mathcal{*}}\left(t\right)e\left(t\right)\end{aligned}\right. (29) 式中： {\boldsymbol{e}}_{\rm{c}}\left(t\right) 为测量结果所产生的新息； \boldsymbol{P}\left(t\right) 为递推修正依赖于每次更新的协方差矩阵； \hat{\boldsymbol{\theta }}\left(t\right)为向量 {\boldsymbol{\theta }}_{\mathrm{*}} 的估计值。
为充分地利用之前的信息并提高参数估计收敛速度，记 p\in {\mathbf{N}}^{+} 为新息数量，对其中的输出矩阵 {\boldsymbol{y}}_{*}\left(t\right) 、信息矩阵 {\mathcal{R}}_{\mathcal{*}}\left(t\right) 进行扩展为
\left\{\begin{aligned} &Y\left(p,t\right)={\left[{\boldsymbol{y}}_{*}\left(t\right),{\boldsymbol{y}}_{*}\left(t-1\right),\cdots ,{\boldsymbol{y}}_{*}\left(t-p+1\right)\right]}^{\mathrm{T}}\\& {\boldsymbol{X}}\left(p,t\right)=\left[{\mathcal{R}}_{\mathcal{*}}\left(t\right),{\mathcal{R}}_{\mathcal{*}}\left(t-1\right),\cdots ,{\mathcal{R}}_{\mathcal{*}}\left(t-p+1\right)\right]\end{aligned}\right. (30) 则新息向量可以为
{\boldsymbol{E}}\left(p,t\right)=\left[\begin{matrix}{\boldsymbol{e}}_{\rm{c}}\left(t\right)\\ {\boldsymbol{e}}_{\rm{c}}\left(t-1\right)\\ \vdots \\ {\boldsymbol{e}}_{\rm{c}}\left(t-p+1\right)\end{matrix}\right]={\boldsymbol{Y}}\left(p,t\right)-{\boldsymbol{X}}^{\mathrm{T}}\left(p,t\right)\hat{\boldsymbol{\theta }}\left(t-1\right) (31) 最终可得到多新息递推最小二乘算法为
\left\{\begin{aligned} &{\boldsymbol{E}}\left(p,t\right)={\boldsymbol{Y}}\left(p,t\right)-{\boldsymbol{ X}}^{\mathrm{T}}\left(p,t\right)\hat{\boldsymbol{\theta }}\left(t-1\right)\\& {\boldsymbol{P}}^{-1}\left(t\right)={\boldsymbol{P}}^{-1}\left(t-1\right)+{\boldsymbol{X}}\left(p,t\right){\boldsymbol{X}}^{\mathrm{T}}\left(p,t\right)\\& \hat{\boldsymbol{\theta }}\left(t\right)=\hat{\boldsymbol{\theta }}\left(t-1\right)+{\boldsymbol{P}}\left(t\right){\boldsymbol{X}}\left(p,t\right){\boldsymbol{E}}\left(p,t\right)\end{aligned}\right. (32) 式（32）所需的扰动估计值可通过设计的扩张状态观测器获取，进而可以对辨识模型（20）应用算法式（32），辨识 {\boldsymbol{\theta }}_{u} ， {\boldsymbol{\theta }}_{v} ， {\boldsymbol{\theta }}_{r} 。因此扩张状态观测器对模型中的未知扰动进行估计，应同时与多新息递推最小二乘算法辨识得到的模型参数不断交互，从而提高无人艇运动状态受到外部不确定扰动情况下的辨识精度，其算法结构图如图3所示。
3. 最优状态调节器设计
3.1 线性二次型问题
当研究对象为线性系统时，所采用的性能指标为状态变量和控制变量的二次型函数积分。求解使性能指标最优化的控制量的问题，被称为线性二次型问题。线性二次型调节器（LQR）是解决该问题的一种方法。LQR能够设计出状态反馈控制器K，从而使二次型目标函数J取最小值。线性二次型最优控制的解可以写成统一的表达形式，其求解过程较为规范，且可构成一个闭环线性状态反馈控制，具有结构简单、便于操作以及控制精度高的特点，因而在实际工程问题中得到了广泛应用[19]。
然而，船舶模型为非线性时变系统，无法直接应用于线性二次型问题。此外，LQR控制器通常不对控制量进行约束，也不考虑外界干扰的影响。因此，提出一种递推计算的状态调节器。在采样间隔时间内，船舶的运动状态可近似为常值，从而可将船舶运动模型视为线性时变系统。通过引入加权矩阵与罚函数，将船舶轨迹误差、外界干扰量、控制量以及约束不等式构成复杂的二次型性能指标。借助连续系统的极小值原理和精细积分法，推导出最优控制量与船舶运动参数之间的解析关系。
3.2 状态调节器设计
若记 \boldsymbol{\eta }={\left[x,y,\psi \right]}^{\mathrm{T}} ， \boldsymbol{\nu }={\left[u,v,r\right]}^{\mathrm{T}} ， \boldsymbol{w}=[{\hat{w}}_{u}\left(t\right),\hat{w}_{v}\left(t\right), {\hat{w}}_{r}\left(t\right)]^{\mathrm{T}} ，通过参数辨识后的船舶三自由度运动模型可简化为
\left\{\begin{aligned} &\dot{\boldsymbol{\eta }}={\boldsymbol{R}}\left(\psi \right){\boldsymbol{\nu}} \\& \dot{\boldsymbol{\nu }}+{\boldsymbol{O}}\left(\boldsymbol{\nu }\right){\boldsymbol{\nu}} +{\boldsymbol{F}}\left(\boldsymbol{\nu }\right){\boldsymbol{\nu}} ={\boldsymbol{\tau}} +{\boldsymbol{w}}\end{aligned}\right. (33) {\boldsymbol{R}}\left(\psi \right)=\left[\begin{matrix} \mathrm{cos}\psi & -\mathrm{sin}\psi & 0\\ \mathrm{sin}\psi & \mathrm{cos}\psi & 0\\ 0& 0& 1 \end{matrix}\right] (34) {\boldsymbol{O}}\left(\boldsymbol{\nu }\right)=-\left[\begin{matrix}{\theta }_{12}u& {\theta }_{11}r+{\theta }_{12}v& {\theta }_{14}r\\ {\theta }_{21}r& {\theta }_{22}& {\theta }_{23}\\ {\theta }_{36}u& {\theta }_{31}& {\theta }_{32}\end{matrix}\right] (35) {\boldsymbol{F}}\left(\boldsymbol{\nu }\right)=-\left[\begin{matrix} 0& 0& 0\\ 0& {\theta }_{24}\left|v\right|+{\theta }_{26}\left|r\right|& {\theta }_{25}\left|r\right|\\ 0& {\theta }_{33}\left|v\right|& {\theta }_{25}\left|r\right|+{\theta }_{35}\left|v\right| \end{matrix}\right] (36) {\boldsymbol{\tau}} =\left[\begin{matrix} {\tau }_{\mathrm{c}u}\\ {\tau }_{\mathrm{c}v}\\ {\tau }_{\mathrm{c}r} \end{matrix}\right]=\left[\begin{matrix} {\theta }_{15}{un}_{\mathrm{P}}+{\theta }_{16}{un}_{\mathrm{S}}+{\theta }_{17}{n}_{\mathrm{P}}^{2}+{\theta }_{18}{n}_{\mathrm{S}}^{2}\\ 0\\ {\theta }_{37}{u}^{2}+{\theta }_{37}{un}_{\mathrm{P}}+{\theta }_{38}{un}_{\mathrm{S}}+{\theta }_{39}{n}_{\mathrm{P}}^{2}+{\theta }_{30}{n}_{\mathrm{S}}^{2} \end{matrix}\right] (37) 式中： \boldsymbol{R}\left(\psi \right)\in {\mathbf{R}}^{3\times 3} ， \boldsymbol{O}\left(\boldsymbol{\nu }\right)\in {\mathbf{R}}^{3\times 3} ， \boldsymbol{F}\left(\boldsymbol{\nu }\right)\in {\mathbf{R}}^{3\times 3} ， \boldsymbol{\tau }\in {\mathbf{R}}^{3\times 1} ，为控制力和力矩。若取系统状态变量 \boldsymbol{x}\left(t\right)=[x,y, \psi ,u,v,r]^{\mathrm{T}} ，控制变量 \boldsymbol{u}={[{\tau }_{\mathrm{c}u},{\tau }_{\mathrm{c}v},{\tau }_{\mathrm{c}r}]}^{\mathrm{T}} ，则可将式（33）改写为
\left\{\begin{aligned} &{\dot{\boldsymbol{x}}\left(t\right)}=f\left(\boldsymbol{x}\left(t\right),\boldsymbol{u}\left(t\right),t\right)={\boldsymbol{A}}\left(t\right){\boldsymbol{x}}\left(t\right)+{\boldsymbol{B}}\hat{\boldsymbol{u}}\left(t\right)\\& {\boldsymbol{y}}\left(t\right)={\boldsymbol{Cx}}\left(t\right)\end{aligned}\right. (38) {\boldsymbol{A}}\left(t\right)=\left[\begin{matrix} {{0}_{3\times 3}} & {\boldsymbol{R}\left(\psi \right)}\\ {{0}_{3\times 3}} & {-[\boldsymbol{O}\left(\boldsymbol{\nu }\right)+\boldsymbol{F}\left(\boldsymbol{\nu }\right)]} \end{matrix}\right] (39) {\boldsymbol{B}}={\left[{0}_{3\times 3}\;\;\; {\boldsymbol{{I}}}_{3\times 3}\right]}^{\mathrm{T}} (40) \hat{\boldsymbol{u}}\left(t\right)={\boldsymbol{u}}\left(t\right)+{\boldsymbol{w}}\left(t\right) (41) 式中 ：\boldsymbol{y}\left(t\right) 为l维输出向量；C为定常矩阵。
在直接航迹控制方法中，船舶的运动状态被作为系统输出，若 {\boldsymbol{y}}_{l}\left(t\right) 表示l维航迹上的参考点状态，则 \boldsymbol{e}\left(t\right)={\boldsymbol{y}}_{l}\left(t\right)-\boldsymbol{y}\left(t\right) 为误差向量。此时，需确定最优控制 \boldsymbol{u}\left(t\right) ，以使下列二次型性能指标达到极小值
\begin{split} & J={\bf\textit{φ}}(\boldsymbol{x}\left({t}_{f}\right),\boldsymbol{u}\left({t}_{f}\right)+\int _{{t}_{0}}^{{t}_{f}}\boldsymbol{L}\left(x,u,t\right)\mathrm{d}t=\\& \frac{1}{2}{{[\boldsymbol{y}}_{l}\left({t}_{f}\right)-\boldsymbol{y}\left({t}_{f}\right)]}^{\mathrm{T}}\boldsymbol{S}{[\boldsymbol{y}}_{l}\left({t}_{f}\right)+-\boldsymbol{y}({t}_{f}\left)\right]+ \\& \frac{1}{2}\int _{0}^{{t}_{f}}{[\boldsymbol{e}}^{\mathrm{T}}\left(t\right)\boldsymbol{Q}\left(t\right)\boldsymbol{e}\left(t\right)+{\hat{\boldsymbol{u}}}^{\mathrm{T}}\left(t\right)R\left(t\right)\hat{\boldsymbol{u}}\left(t\right)]{\rm{d}}t \end{split} (42) 式中：S为 l\times l 维对称非负定常矩阵； \boldsymbol{Q}\left(t\right) 为 l\times l 维对称非负定时变常矩阵； \boldsymbol{R}\left(t\right) 为 3\times 3 维对称非负定时变常矩阵；末端时刻 {t}_{f} 为固定且为有限值。在线性二次型问题中， 由于采用二次型目标函数进行最优控制，船舶的转速和舵角存在限制， \boldsymbol{u}\left(t\right) 则存在不等式约束 {\boldsymbol{u}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\le \boldsymbol{u}\le {\boldsymbol{u}}_{\mathrm{m}\mathrm{a}\mathrm{x}} ，并可写成一般形式 \boldsymbol{h}\left(u\right)=\boldsymbol{a}u+{\boldsymbol{b}}\le 0 ，若i为不等式约束的个数，则a为 i\times 3 维矩阵，b为 i\times 1 维矩阵。设存在加权矩阵 \boldsymbol{G}=\mathrm{diag}({g}_{1},{g}_{2}，\dots ，{g}_{i}) 且满足
{g}_{i}=\left\{\begin{aligned} &\zeta > 0,&& h\left({u}_{i}\right) > 0\\& 0, && h\left({u}_{i}\right)\le 0\end{aligned}\right. (43) 则可将目标函数扩展为
\begin{gathered} J=\frac{1}{2}{{[\boldsymbol{y}}_{l}\left({t}_{f}\right)-\boldsymbol{y}\left({t}_{f}\right)]}^{\mathrm{T}}{\boldsymbol{S}}{[\boldsymbol{y}}_{l}\left({t}_{f}\right)-{\boldsymbol{y}}\left({t}_{f}\right)] +\\ \frac{1}{2}\int _{0}^{{t}_{f}}{\boldsymbol{e}}^{\mathrm{T}}\left(t\right)\boldsymbol{Q}\left(t\right)\boldsymbol{e}\left(t\right){+\hat{\boldsymbol{u}}}^{\mathrm{T}}\left(t\right)\boldsymbol{R}\left(t\right)\hat{\boldsymbol{u}}\left(t\right)+ {\boldsymbol{h}\left(u\right)}^{\mathrm{T}}{\boldsymbol{Gh}}\left(u\right){\rm{d}}t \end{gathered} (44) 若 {\boldsymbol{u}}^{*}\left(t\right) 和 {t}_{f}^{*} 是使性能指标取最小值的最优解， {\boldsymbol{x}}^{*}\left(t\right) 为相应的最优轨线，则存在n维向量函数 \boldsymbol{\lambda }\left(t\right) ，可利用最优控制理论中Hamilton函数
\begin{gathered} H=\boldsymbol{L}\left(\boldsymbol{x}\left(t\right),\boldsymbol{u}\left(t\right),t\right)+{\boldsymbol{\lambda }}^{\mathrm{T}}\left(t\right)f\left(\boldsymbol{x}\left(t\right),\boldsymbol{u}\left(t\right),t\right)=\\ \frac{1}{2}{\left[{\boldsymbol{y}}_{l}\left(t\right) - \boldsymbol{C}\boldsymbol{x}\left(t\right)\right]}^{\mathrm{T}}\boldsymbol{Q}\left(t\right)\left[{\boldsymbol{y}}_{l}\left(t\right)-\boldsymbol{C}\boldsymbol{x}\left(t\right)\right] + \frac{1}{2}{\hat{\boldsymbol{u}}}^{\mathrm{T}}\left(t\right)\boldsymbol{R}\left(t\right)\hat{\boldsymbol{u}}\left(t\right)+\\ \frac{1}{2}{\boldsymbol{h}}^{\mathrm{T}}\left(u\right)Gh\left(u\right)+{{\boldsymbol{x}}^{\mathrm{T}}\left(t\right)\boldsymbol{A}}^{\mathrm{T}}\left(t\right){\boldsymbol{\lambda }}\left(t\right)+{\hat{\boldsymbol{u}}}^{\mathrm{T}}\left(t\right){\boldsymbol{B}}^{\mathrm{T}}{\boldsymbol{\lambda }}\left(t\right) \end{gathered} (45) 由H的极值条件有 \partial H/\partial \hat{\boldsymbol{u}}=0 ，且考虑到 \dot{\boldsymbol{h}}\left(u\right)=\boldsymbol{a} ，则可得最优控制
{\boldsymbol{u}}^{*}\left(t\right)=-{\tilde {{\boldsymbol{R}}}}^{-1}\left(t\right)\left[{\boldsymbol{B}}^{\mathrm{T}}\boldsymbol{\lambda }\left(t\right)+\tilde {\boldsymbol{b}}\right] (46) 式中： \tilde {\boldsymbol{R}}\left(t\right)=\boldsymbol{R}\left(t\right)+{\boldsymbol{a}}^{\mathrm{T}}\boldsymbol{G}\boldsymbol{a} ； \tilde {\boldsymbol{b}}={\boldsymbol{a}}^{\mathrm{T}}\boldsymbol{G}\boldsymbol{b} 。
在每一时刻，辨识算法需先于执行控制算法，两者之间存在时间延迟，所以外界造成的扰动量可以作为下一个时刻的估计值，则控制变量为
{\boldsymbol{u}}\left(t\right)={\boldsymbol{u}}^{*}\left(t\right)-\hat{\boldsymbol{w}}\left(t+1\right) (47) 由正则方程可得
{\dot{\boldsymbol{x}}\left(t\right)}=\frac{\partial H}{\partial \lambda }={\boldsymbol{A}}\left(t\right)x\left(t\right)-{\boldsymbol{B}}{\tilde {R}}^{-1}\left(t\right){\boldsymbol{B}}^{\mathrm{T}}{\boldsymbol{\lambda }}\left(t\right) -{\boldsymbol{B}}{\tilde {R}}^{-1}\left(t\right)\tilde {\boldsymbol{b}} (48) {\dot{\boldsymbol{\lambda }}\left(t\right)}=-\frac{\partial H}{\partial \boldsymbol{x}}={\boldsymbol{C}}^{\mathrm{T}}{\boldsymbol{Q}}\left(t\right)\left[{\boldsymbol{y}}_{l}\left(t\right)-\boldsymbol{C}\boldsymbol{x}\left(t\right)\right] -{\boldsymbol{A}}^{\mathrm{T}}\left(t\right){\boldsymbol{\lambda }}\left(t\right) (49) 由横截条件，有
{\boldsymbol{\lambda }}\left({t}_{f}\right)=\frac{\partial }{\partial \boldsymbol{x}\left({t}_{f}\right)}\left[\frac{1}{2}{\boldsymbol{e}}^{\mathrm{T}}\left({t}_{f}\right){\boldsymbol{S}}e\left({t}_{f}\right)\right] = -{\boldsymbol{C}}^{\mathrm{T}}{\boldsymbol{S}}\left[{\boldsymbol{y}}_{l}\left({t}_{f}\right)-\boldsymbol{C}\boldsymbol{x}\left({t}_{f}\right)\right] (50) 根据正则方程为线性方程，则可假设有
{\boldsymbol{\lambda }}\left(t\right)={\boldsymbol{P}}\left(t\right)x\left(t\right)-{\boldsymbol{g}}\left(t\right) (51) 式中： \boldsymbol{P}\left(t\right) 为 l\times l 维对称非负定实矩阵； \boldsymbol{g}\left(t\right) 为 l\times 1 维伴随向量。
对式（51）进行求导，可得
\begin{split} & {\dot{\boldsymbol{\lambda }}\left(t\right)}={\dot{\boldsymbol{P}}\left(t\right)}\boldsymbol{x}\left(t\right)+\boldsymbol{P}\left(t\right){\dot{\boldsymbol{x}}\left(t\right)}-{\dot{\boldsymbol{g}}\left(t\right)}= {\dot{\boldsymbol{P}}\left(t\right)}\boldsymbol{x}\left(t\right)+\\&\quad \boldsymbol{P}\left(t\right)[\boldsymbol{A}\left(t\right)\boldsymbol{x}\left(t\right)-\boldsymbol{B}\left(t\right){\tilde {R}}^{-1}\left(t\right){\boldsymbol{B}}^{\mathrm{T}}\left(t\right){\boldsymbol{\lambda }}\left(t\right)-\\& \qquad\qquad {\boldsymbol{B}}\left(t\right){\tilde {R}}^{-1}\left(t\right)\tilde {\boldsymbol{b}}]-\dot{\boldsymbol{g}\left(t\right)} \end{split} (52) 将式（48）代入式（52），并令式（49）与式（52）两端相等，则有
\begin{split} & -{\dot{\boldsymbol{P}}\left(t\right)}={\boldsymbol{P}}\left(t\right){\boldsymbol{A}}\left(t\right)+{\boldsymbol{A}}^{\mathrm{T}}\left(t\right){\boldsymbol{P}}\left(t\right)+{\boldsymbol{C}}^{\mathrm{T}}{\boldsymbol{Q}}\left(t\right){\boldsymbol{C}}-\\& \qquad\qquad {\boldsymbol{P}}\left(t\right){\boldsymbol{B}}\left(t\right){\tilde {R}}^{-1}\left(t\right){\boldsymbol{B}}^{\mathrm{T}}\left(t\right){\boldsymbol{P}}\left(t\right) \end{split} (53) \begin{split} & -{\dot{\boldsymbol{g}}\left(t\right)}={\boldsymbol{P}}\left(t\right){\boldsymbol{B}}\left(t\right){\tilde {R}}^{-1}\left(t\right)\tilde {\boldsymbol{b}}+{\boldsymbol{C}}^{\mathrm{T}}{\boldsymbol{Q}}\left(t\right){\boldsymbol{y}}_{l}\left(t\right)+\\& \left[{\boldsymbol{A}}^{\mathrm{T}}\left(t\right)-\boldsymbol{P}\left(t\right)\boldsymbol{B}\left(t\right){\tilde {R}}^{-1}\left(t\right){\boldsymbol{B}}^{\mathrm{T}}\left(t\right)\boldsymbol{P}\left(t\right)\right]{\boldsymbol{g}}\left(t\right) \end{split} (54) 式 （53）为矩阵黎卡提微分方程，式 （54） 为伴随矩阵。
当式（51）中 t=t\mathit{\mathit{\mathit{\mathrm{_{\mathit{f}}}}}} 并与式（50）两端相等，可得
{\boldsymbol{P}}\left({t}_{f}\right)={\boldsymbol{C}}^{\mathrm{T}}{\boldsymbol{SC}},{\boldsymbol{g}}\left({t}_{f}\right)={\boldsymbol{C}}^{\mathrm{T}}{\boldsymbol{S}}{\boldsymbol{y}}_{l}\left({t}_{f}\right) (55) 为获取最优控制，需求解矩阵黎卡提微分方程及伴随矩阵方程。借助MATLAB的龙格−库塔法求解伴随矩阵方程，可得伴随向量g(t)。因此，求解有限时间的矩阵黎卡提微分方程是获得最优解的关键环节。
3.3 精细积分法
精细积分法[20]摒弃了传统求解动力方程常用的差分格式，是一种基于 {2}^{N} 类算法的精细积分方法，主要思想是将原本很小的步长细分成 {2}^{N} 小段，在每一小段内利用泰勒级数展开对指数矩阵进行精细求解，由于细分后的步长很小，截断误差可以忽略不计[21]。此方法适用于解决复杂定积分问题，可用于二点边值问题以及矩阵微分方程组的求解。在该研究中，由于采用数值迭代求解每一时刻矩阵黎卡提方程的解，故在每一时刻可将 \boldsymbol{A}\left(t\right)，\boldsymbol{B}\left(t\right) 矩阵视为定常矩阵， \boldsymbol{Q}\left(t\right)，\boldsymbol{R}\left(t\right) 等加权矩阵亦可视为定常矩阵。因此式（53）可写成为
\left\{\begin{aligned} &-\dot{\boldsymbol{P}\left(t\right)}={\boldsymbol{PA}}+{\boldsymbol{A}}^{\mathrm{T}}{\boldsymbol{P}}+{\boldsymbol{C}}^{\mathrm{T}}{\boldsymbol{QC}}-{\boldsymbol{PB}}\tilde {\boldsymbol{R}}^{-1}{\boldsymbol{B}}^{\mathrm{T}}{\boldsymbol{P}}\\& {\boldsymbol{P}}\left({t}_{f}\right)={\boldsymbol{C}}^{\mathrm{T}}{\boldsymbol{SC}}\end{aligned}\right. (56) 设时程积分的步长为 \eta ，对 t\mathrm{_{\mathit{f}}} 时段进行分段，有 t_0=0,t_1=\eta,t_2=2\eta,\dots,t\mathit{_{\mathrm{\mathit{f}}}}=k\mathrm{_{\mathit{f}}}\eta ，精细积分法基于步长 \eta 实施，则取 \tau =\eta /{2}^{N} (N取 20 )进行泰勒展开：
\left\{\begin{aligned} &{\boldsymbol{E}}\left(\tau \right)\approx {e}_{1}\tau +{e}_{2}{\tau }^{2}+{e}_{3}{\tau }^{3}+{e}_{4}{\tau }^{4}\\& {\boldsymbol{G}}\left(\tau \right)\approx {g}_{1}\tau +{g}_{2}{\tau }^{2}+{g}_{3}{\tau }^{3}+{g}_{4}{\tau }^{4}\\& F\left(\tau \right)\approx I+{\boldsymbol{F}}'\left(\tau \right)\\& {\boldsymbol{F}}'\left(\tau \right)={f}_{1}\tau +{f}_{2}{\tau }^{2}+{f}_{3}{\tau }^{3}+{f}_{4}{\tau }^{4}\end{aligned}\right. (57) 其中：
{e}_{1}={\boldsymbol{C}}^{\mathrm{T}}\boldsymbol{Q}\boldsymbol{C} {e}_{2}={(\boldsymbol{A}}^{\mathrm{T}}{e}_{1}+{e}_{1}\boldsymbol{A})/2 {e}_{3}={({f}_{2}}^{\mathrm{T}}{e}_{1}+{e}_{1}{f}_{2}+{\boldsymbol{A}}^{\mathrm{T}}{e}_{1}\boldsymbol{A})/3 {e}_{4}={({f}_{3}}^{\mathrm{T}}{e}_{1}+{e}_{1}{f}_{3}+{\boldsymbol{A}}^{\mathrm{T}}{e}_{1}{f}_{2}+{({f}_{3}}^{\mathrm{T}}{e}_{1}\boldsymbol{A}\left)\right)/4 {g}_{1}=\boldsymbol{B}{{\tilde {\boldsymbol{R}}}^{-1}}{\boldsymbol{B}}^{\mathrm{T}} {g}_{2}={({g}_{1}\boldsymbol{A}}^{\mathrm{T}}+\boldsymbol{A}{g}_{1})/2 {g}_{3}={({{g}_{1}f}_{2}}^{\mathrm{T}}+{f}_{2}{g}_{1}+{\boldsymbol{A}{g}_{1}\boldsymbol{A}}^{\mathrm{T}})/3 {g}_{4}={({{g}_{1}f}_{3}}^{\mathrm{T}}+{f}_{3}{g}_{1}+\boldsymbol{A}{g}_{1}{{f}_{2}}^{\mathrm{T}}+{f}_{2}{g}_{1}{\boldsymbol{A}}^{\mathrm{T}})/4 {f}_{1}=\boldsymbol{A} {f}_{2}={(\boldsymbol{A}}^{2}-{g}_{1}{e}_{1})/2 {f}_{3}=\frac{\boldsymbol{A}{f}_{2}-{g}_{2}{e}_{1}-{g}_{1}{e}_{1}\boldsymbol{A}}{3} {f}_{4}={(\boldsymbol{A}f}_{3}-{g}_{3}{e}_{1}+{g}_{2}{e}_{1}\boldsymbol{A}-{g}_{1}{e}_{1}{f}_{2})/4 细分的每一区段需要首尾相连进行区段合并消元，区段合并消元公式[22]为：
\left\{\begin{aligned} &{\boldsymbol{E}}_{\mathrm{c}}={\boldsymbol{E}}_{1}+{\boldsymbol{F}}_{1}^{\mathrm{T}}{\left({\boldsymbol{E}}_{2}^{-1}+{\boldsymbol{G}}^{-1}\right)}^{-1}{\boldsymbol{F}}_{1}\\& {\boldsymbol{F}}_{\mathrm{c}}={\boldsymbol{F}}_{2}{\left(I+{\boldsymbol{G}}_{1}{\boldsymbol{E}}_{2}\right)}^{-1}{\boldsymbol{F}}_{1}\\& {\boldsymbol{G}}_{\mathrm{c}}={\boldsymbol{G}}_{2}+{\boldsymbol{F}}_{2}{\left({\boldsymbol{G}}_{1}^{-1}+{\boldsymbol{E}}_{2}\right)}^{-1}{\boldsymbol{F}}_{2}^{\mathrm{T}}\end{aligned}\right. (58) 在区段合并消元中 {\boldsymbol{F}}'\left(\tau \right) 不应与I直接相加，否则会极大降低精度，因此式（58）中 {\boldsymbol{F}}_{\mathrm{c}} 应改为
\begin{split} & \;\;{\boldsymbol{F}}_{\mathrm{c}}'={\boldsymbol{F}}_{1}'+{\boldsymbol{F}}_{2}'+{\boldsymbol{F}}_{2}'{\boldsymbol{F}}_{1}' -\left(\boldsymbol{I}+{\boldsymbol{F}}_{2}'\right)\left(\boldsymbol{I}+{\boldsymbol{F}}_{1}'\right)\cdot\\& \left[{\left(\boldsymbol{I}+{\boldsymbol{G}}_{1}{\boldsymbol{E}}_{2}\right)}^{-1}{\boldsymbol{G}}_{1}{\boldsymbol{E}}_{2}+{\boldsymbol{G}}_{1}{\boldsymbol{E}}_{2}{\left(\boldsymbol{I}+{\boldsymbol{G}}_{1}{\boldsymbol{E}}_{2}\right)}^{-1}\right]/2 \end{split} (59) 时段 \eta 的 {2}^{N} 算法步骤如下：
1） 输入参数A ,\boldsymbol{B},\boldsymbol{C},\boldsymbol{Q},\boldsymbol{R},\boldsymbol{S},\eta ,N ，按照式（57）计算得 \boldsymbol{E}\left(\tau \right),\boldsymbol{G}\left(\tau \right),{\boldsymbol{F}}'\left(\tau \right) ，并将其分别送存入 {\boldsymbol{E}}_{\mathrm{c}},{\boldsymbol{G}}_{\mathrm{c}},{\boldsymbol{F}}_{\mathrm{c}}' 。
2） \mathrm{f}\mathrm{o}\mathrm{r}\left(i=1;i\le N;i++\right)\{{\boldsymbol{E}}_{1}={\boldsymbol{E}}_{\rm{c}};{\boldsymbol{E}}_{2}={\boldsymbol{E}}_{\rm{c}}；{\boldsymbol{G}}_{1}= {\boldsymbol{G}}_{\mathrm{c}} ；{\boldsymbol{G}}_{2}={\boldsymbol{G}}_{\mathrm{c}};{\boldsymbol{F}}_{1}'={\boldsymbol{F}}_{\mathrm{c}}';{\boldsymbol{F}}_{1}'={\boldsymbol{F}}_{\mathrm{c}}' ；按式（58）和式（59）重新计算 {\boldsymbol{E}}_{\mathrm{c}},{\boldsymbol{G}}_{\mathrm{c}},{\boldsymbol{F}}_{\mathrm{c}}';\} 。
3） {\boldsymbol{F}}_{\mathrm{c}}=\boldsymbol{I}+{\boldsymbol{F}}_{\mathrm{c}}' ，此时计算得到的 {\boldsymbol{E}}_{\mathrm{c}},{\boldsymbol{G}}_{\mathrm{c}},{\boldsymbol{F}}_{\mathrm{c}} 分别为 \boldsymbol{E}\left(\eta \right),\boldsymbol{G}\left(\eta \right),\boldsymbol{F}\left(\eta \right) 。
通过上述步骤，可计算得到 \boldsymbol{E}\left(\eta \right),\boldsymbol{G}\left(\eta \right),\boldsymbol{F}\left(\eta \right) ，进而可通过递推求解黎卡提方程的解。设长度 \eta 为区段1，长度 k\eta 为区段2，依据区段合并消元公式，可计算出 (k+1)\eta 对应的区段。随后，将 (k+1)\eta 重新设为区段2，区段1保持不变，从而能够计算出 (k+2)\eta 的区段。依此类推，可计算得到 \boldsymbol{E}\left(k\mathit{\mathit{\mathrm{_{\mathrm{\mathit{\mathrm{\mathit{f}}}}}}}}\eta\right), \boldsymbol{G}\left(k\mathrm{\mathrm{_{\mathrm{\mathit{\mathrm{\mathit{f}}}}}}}\eta\right),\boldsymbol{F}\left(k\mathit{\mathit{\mathrm{_{\mathrm{\mathit{\mathrm{\mathit{f}}}}}}}}\eta\right) 。然而，由于存在初值条件 \boldsymbol{P}\left(t_{\mathrm{\mathit{f}}}\right)= \boldsymbol{C}^{\mathrm{T}}\boldsymbol{S}\boldsymbol{C} ，根据文献[12]，还存在一小区段满足 {\boldsymbol{E}}_{2}= {\boldsymbol{C}}^{\mathrm{T}}\boldsymbol{S}\boldsymbol{C},{\boldsymbol{F}}_{2}=\boldsymbol{I},{\boldsymbol{G}}_{2}=0 ，将该小区段设为1，任一长度 {\eta ～k}_{f}\eta 设为区段2，从而可得到P，每一段合并公式
{\boldsymbol{P}}\left(t\right)={\boldsymbol{E}}+{\boldsymbol{F}}^{\mathrm{T}}{\left(\boldsymbol{I}+{\boldsymbol{C}}^{\mathrm{T}}\boldsymbol{S}\boldsymbol{C}\boldsymbol{G}\right)}^{-1}{\boldsymbol{C}}^{\mathrm{T}}{\boldsymbol{SCF}} (60) 精细积分算法的具体流程如图4所示。
通过上述算法，可以计算出在时间t时，若使目标函数J最小，矩阵黎卡提方程的解为 \boldsymbol{P}\left(t\right) ，代入式（54）可以计算得到 \boldsymbol{g}\left(t\right) ，则最优控制为
{\boldsymbol{u}}^{\mathrm{*}}\left(t\right)=-{\tilde {R}}^{-1}\left(t\right)\left\{{\boldsymbol{B}}^{\mathrm{T}}\left(t\right)[\boldsymbol{P}\left(t\right)\boldsymbol{x}\left(t\right)\right.-{\boldsymbol{g}}\left(t\right)]+\tilde {\boldsymbol{b}}\} (61) 从式（61）可知，船舶轨迹跟踪最优控制与系数矩阵 \boldsymbol{B},\boldsymbol{R},\boldsymbol{P},\boldsymbol{g},\boldsymbol{b},\boldsymbol{G} 等有关，在每一周期的运算中，若船舶的状态保持不变且约束已提前设定，则轨迹跟踪的实现取决于矩阵 \boldsymbol{R},\boldsymbol{P},\boldsymbol{G} 。最优控制解存在的条件是矩阵黎卡提方程有解。根据该条件可得， \left\{\boldsymbol{A}\right.,\left.\boldsymbol{B}\right\} 需可控且 [\boldsymbol{B},\boldsymbol{A}\boldsymbol{B}] 应为满秩。根据精细积分法，矩阵黎卡提方程的可解性与系数矩阵 \boldsymbol{Q},\boldsymbol{R},\boldsymbol{S},\boldsymbol{G} 有关。由式（44）可知，矩阵Q为航迹误差在目标函数中的权重，矩阵S为最终状态在目标函数中的权重，R为控制量在目标函数中的权重，G为约束在目标函数中的权重，且G的取值可与R相同。因此，调节 \boldsymbol{Q},\boldsymbol{R},\boldsymbol{S} 的值是决定航迹跟踪性能的关键。
为实现船舶的在线辨识和控制，将参数辨识算法与控制算法关联，其具体步骤如下：
1） 将获得的船舶运动数据输入设计的扩张状态观测器−多新息递推最小二乘算法，可以获得t时刻的辨识参数 \hat{\boldsymbol{\theta }}\left(t\right) 和扰动估计 \hat{\boldsymbol{w}}\left(t\right) 。
2） 将辨识的参数和扰动估计代入状态方程中，将航迹跟踪问题改写为线性二次型优化控制问题。
3） 输入参数 \boldsymbol{C},\boldsymbol{Q},\boldsymbol{R},\boldsymbol{S},\eta ,N ，利用精细积分法可求解矩阵 \boldsymbol{P}\left(t\right) ，进而通过式（54）和式（61）计算最优控制 {u}^{\mathrm{*}}\left(t\right) 。
4） 设计的扩张状态观测器−多新息递推最小二乘算法可获得 t+1 时刻的辨识参数 \hat{\boldsymbol{\theta }}\left(t+1\right) 和扰动估计 \hat{\boldsymbol{w}}(t+1) 。
5） 将 {u}^{\mathrm{*}}\left(t\right) 代入式（47），可得控制变量 {\boldsymbol{u}}\left(t\right) 。将 {\boldsymbol{u}}\left(t\right) 代入状态方程，可得船舶轨迹 x\left(t\right) 。最终，将 x\left(t\right) 更新至运动数据中。
6） 按照上述循环过程，参数辨识算法为控制算法提供辨识的参数和估计的扰动量，而控制算法则为辨识算法提供船舶运动的数据，以此实现在线辨识和控制的目标。
参数辨识算法与控制算法的联合策略流程图如图5所示。
![]()
图 5 参数辨识与控制算法的流程图Figure 5. Flow chart of parameter identification and control algorithm4. 仿真实验
4.1 参数辨识仿真
仿真实验利用海洋系统模拟器（MSS）工具箱中Otter号无人艇作为参数辨识仿真实验对象。为提高参数辨识的准确性，将无人艇的预定运动轨迹为连续的型实验，无人艇的指令信号为左、右螺旋桨转速 ({n}_{\mathrm{S}},{n}_{\mathrm{P}}) 随时间（采样周期h=0.2 s，采样周期数为0~3 000）的变化情况，如图6所示。
由于该无人艇由双螺旋桨驱动，艇前进运动与转向运动依靠于左右螺旋桨的转速大小，假设仿真中存在如下扰动量：
\left\{\begin{aligned} &{w}_{u}\left(t\right)=0.1\mathrm{sin}\left(0.01t+3{g}_{1}\left(t\right)\right)\\& {w}_{v}\left(t\right)=0.1\mathrm{cos}\left(0.01t+3{g}_{2}\left(t\right)\right)\\& {w}_{r}\left(t\right)=0.1\mathrm{cos}\left(0.01t+2{g}_{3}\left(t\right)\right)\end{aligned}\right. (62) 式中： {g}_{1} ， {g}_{2} ， {g}_{3} 为 \left[\mathrm{0,1}\right] 的随机变量且相互独立。
在式（62）所描述的扰动条件下，无人艇按照预定的螺旋桨转速运行。随后，将船舶仿真运动数据逐步输入至扩张状态观测－多新息递推最小二乘交互式算法（EMS）中。通过这一过程，可不断更新干扰估计值与参数辨识值。最终，将更新后的值与船舶仿真运动数据进行对比，得到图7～图9。
从图7～图9可以看出，参数辨识算法所估算的纵荡、横荡和艏摇角速度与实际数据的差异较小。这表明该参数辨识算法估计出的状态量 \hat{\boldsymbol{\nu }}\left(t\right) 和扰动量 \hat{\boldsymbol{w}}\left(t\right) 具有一定的准确性。当迭代周期次数约大于1 500时，该算法辨识出的参数趋于稳定，如表1所示。
表 1 辨识参数Table 1. Identification parameters序号 参数 数值 序号 参数 数值 1 {\boldsymbol{\theta }}_{11} 0.265 138 669 258 724 13 {\boldsymbol{\theta }}_{25} 0.305 785 017 250 565 2 {\boldsymbol{\theta }}_{12} −0.593 139 836 433 817 14 {\boldsymbol{\theta }}_{26} 0.817 640 200 880 232 3 {\boldsymbol{\theta }}_{13} 0.055 061 829 481 853 15 {\boldsymbol{\theta }}_{31} 0.061 005 725 705 921 4 {\boldsymbol{\theta }}_{14} 0.193 328 379 010 969 16 {\boldsymbol{\theta }}_{32} −0.247 108 873 113 790 5 {\boldsymbol{\theta }}_{15} 0.004 535 423 216 801 17 {\boldsymbol{\theta }}_{33} −0.071 435 501 716 979 6 {\boldsymbol{\theta }}_{16} 0.001 500 027 456 607 18 {\boldsymbol{\theta }}_{34} −0.067 124 874 659 099 7 {\boldsymbol{\theta }}_{17} −0.00 001 031 479 303 19 {\boldsymbol{\theta }}_{35} 0.150 310 762 377 552 8 {\boldsymbol{\theta }}_{18} −2.101 476 849 154 E-09 20 {\boldsymbol{\theta }}_{36} 0.156 863 665 033 070 9 {\boldsymbol{\theta }}_{21} −0.057 189 322 465 896 21 {\boldsymbol{\theta }}_{37} −0.000 528 418 571 230 10 {\boldsymbol{\theta }}_{22} −0.107 906 938 569 009 22 {\boldsymbol{\theta }}_{38} −0.001 535 487 343 785 11 {\boldsymbol{\theta }}_{23} −0.021 188 406 958 892 23 {\boldsymbol{\theta }}_{39} 7.379 676 517 343 E-06 12 {\boldsymbol{\theta }}_{24} −0.343 953 417 948 503 24 {\boldsymbol{\theta }}_{30} −6.976 693 955 860 E-07 4.2 控制算法仿真
在辨识的起始阶段，辨识出的参数波动幅度较大，这可能会对控制算法仿真分析产生较大干扰。为避免这种干扰，在控制算法仿真中，船舶运动参数选用表1中的数据，外界扰动量w选用辨识算法周期数为
1500 ～2250 的扰动估计值。该控制算法的关键在于选取合适的加权矩阵 \boldsymbol{Q},\boldsymbol{R},\boldsymbol{S} ，这直接关系到最优控制解的精度。对于航迹跟踪，目标是在控制量不超过船舶自身限制的前提下，使消耗的能量最小，同时尽可能最小化船舶跟踪误差。因此，采用式（18）计算每一周期的航迹误差并求其误差平方总和。
由于采用直接航迹误差，系统的输出仅为船舶的实际坐标点。因此，可选取矩阵 \boldsymbol{C}=\mathrm{diag}\{\mathrm{1,1}, \mathrm{0,0},\mathrm{0,0}\} ， l=6 。由于采用迭代求解最优控制，在每次迭代中将船舶状态量视为常量，因此每次迭代的时间应尽可能小，同时为了降低算法的复杂度，选取每次迭代时间为 h=0.2\;\mathrm{s} ，即式（44）中tf取0 ～ 0.2\; \mathrm{s} 。采用精细积分法时，选取时段 \eta =0.025 ，则有kf取0 ～ 8 。
设t为采样时间，( {x}_{\mathrm{d}}\left(t\right), {\boldsymbol{y}}_{\mathrm{d}}\left(t\right)) 为期望参考点，且参考点是以速度U进行变化，在 t=0 时，船舶初始状态 \boldsymbol{x}\left(0\right)={[\mathrm{0,0},\mathrm{0,0.01,0},0]}^{\mathrm{T}} ，则可设期望输出向量 {\boldsymbol{y}}_{l}\left(t\right)={\left[{x}_{\mathrm{d}}\left(t\right), {\boldsymbol{y}}_{\mathrm{d}}\left(t\right),\mathrm{0,0},\mathrm{0,0}\right]}^{\mathrm{T}} ，为使船舶跟踪更加平缓，参考轨迹的速度U引入双曲正切函数，设 {u}_{\mathrm{d}} 为轨迹期望速度，则有
U={u}_{\mathrm{d}}\times \mathrm{tanh}\left(\frac{t}{100}\right) (63) 由于左右螺旋桨的转速限制为 0～120\;\mathrm{r}\mathrm{a}\mathrm{d}/\mathrm{s} ，根据表1数据和式（37）可得：
\left\{\begin{aligned} &0{\le \tau }_{\mathrm{c}u}\le 0.5\\& {-0.25\le \tau }_{\mathrm{c}r}\le 0.15\end{aligned}\right. (64) 从而将式（64）转换为如下不等式：
\left\{\begin{aligned} &{-\tau }_{\mathrm{c}u}\le 0\\& {\tau }_{\mathrm{c}u}-0.5\le 0\\& {-\tau }_{\mathrm{c}r}-0.25\le 0\\& {\tau }_{\mathrm{c}r}-0.15\le 0\end{aligned}\right. (65) 为将式（65）转化为形如 \boldsymbol{h}\left(u\right)=\boldsymbol{a}u+\boldsymbol{b}\le 0 的矩阵不等式，其中 \boldsymbol{u}={[{\tau }_{\mathrm{c}u},{\tau }_{\mathrm{c}v},{\tau }_{\mathrm{c}r}]}^{\mathrm{T}} ，则不等式约束个数 i=4 ，矩阵a， b如下：
\boldsymbol{a}=\left[\begin{matrix} -1& 0& 0\\ 1& 0& 0\\ 0& 0& -1\\ 0& 0& 1\end{matrix}\right],\;\;\boldsymbol{b}=\left[\begin{matrix} 0\\ -0.5\\ -0.25\\ -0.15\end{matrix}\right] 为得到使轨迹误差最小的矩阵 \boldsymbol{Q},\boldsymbol{R},\boldsymbol{S} ，选用 {x}_{\mathrm{d}}\left(t\right)=Ut\mathrm{c}\mathrm{o}\mathrm{s}\left(\dfrac{{\text{π}}}{4}\right)\mathrm{和}{\boldsymbol{y}}_{\mathrm{d}}\left(t\right)=Ut\mathrm{s}\mathrm{i}\mathrm{n}({\text{π}}/4) 为参考轨迹，参考轨迹点的期望速度 u_{\mathrm{d}}=1\ \mathrm{m}/\mathrm{s} ，参数 \zeta =6\; 000 。采用控制变量法过调节矩阵 \boldsymbol{Q},\boldsymbol{R},\boldsymbol{S} 值，并进行多组实验，每一组实验迭代 1\;500 次。根据矩阵 \boldsymbol{Q},\boldsymbol{R},\boldsymbol{S} 实际意义和经验选取下列基础数据：
\left\{\begin{aligned} &{\boldsymbol{Q}}_{0}=500\times {\mathrm{diag}}\left\{\mathrm{1,1},\mathrm{1,1},\mathrm{1,1}\right\}\\& {\boldsymbol{R}}_{0}= 6\; 000\times {\mathrm{diag}}\left\{\mathrm{0.4,5},5\right\}\\& {\boldsymbol{S}}_{0}=500\times {\mathrm{diag}}\left\{\mathrm{1,1},\mathrm{1,1},\mathrm{1,1}\right\}\end{aligned}\right. (66) 1） 固定 \boldsymbol{Q}={\boldsymbol{Q}}_{0},{\boldsymbol{R}=\boldsymbol{R}}_{0} ，S分别为 {\boldsymbol{S}}_{0},\;2{\boldsymbol{S}}_{0},3{\boldsymbol{S}}_{0}, 4{\boldsymbol{S}}_{0} 做4组实验，可得到每一组实验的船舶轨迹、轨迹纵向误差和船舶状态，并将其进行对比，得到图10和图11。
![]()
图 10 不同S值船舶轨迹与纵向误差对比图Figure 10. Comparison of ship trajectory and longitudinal error with different S values![]()
图 11 不同S 值艏向角与速度变化对比图Figure 11. Comparison of heading angle and velocity changes with different S values为能够更清晰分析不同S值对航迹跟踪效果的影响，在不同S值下，将峰值纵向误差、纵向误差平方总和、峰值艏向角、平稳艏向角时间、平稳速度时间整理成表2。其中平稳时间的衡量标准是误差与最终稳定数值的比值在 1\text{%} 内的时间，且后面数据也保持在该范围内。
表 2 不同S数据对比Table 2. Data comparison table of different S参数 峰值纵向
误差/m误差平方
总和/m2峰值
艏向角/rad平稳艏向角
时间/s平稳速度
时间/s{\boldsymbol{S}}_{0} −2.80 591 1.03 67.2 122.7 {2\boldsymbol{S}}_{0} −2.12 246 1.09 70.9 71.4 {3\boldsymbol{S}}_{0} −1.78 151 1.13 73.4 52.4 4{\boldsymbol{S}}_{0} −1.57 114 1.15 75.7 42.0 从表2和图10、图11中可以看出，随着S值的增大，误差峰值和误差平方总和均减小，平稳速度所需时间缩短，但峰值艏向角越来越大，平稳艏向角时间延长。因此，增加S可以加快响应时间，但随着其不断增大，会导致航向稳定时间变长。
2） 固定 \boldsymbol{S}={\boldsymbol{S}}_{0},{\boldsymbol{R}=\boldsymbol{R}}_{0} ，Q分别为 {\boldsymbol{Q}}_{0}, 3{\boldsymbol{Q}}_{0}, 6{\boldsymbol{Q}}_{0}, 9{\boldsymbol{Q}}_{0} ，进行4组实验。每组实验得到的船舶轨迹、轨迹纵向误差和船舶状态，均进行了对比，结果如图12和图13所示。
![]()
图 12 不同Q值船舶轨迹与纵向误差对比图Figure 12. Comparison of ship trajectory and longitudinal error with different Q values![]()
图 13 不同Q值艏向角与速度变化对比图Figure 13. Comparison of heading angle and velocity changes with different Q values在不同Q值下得到的峰值纵向误差、纵向误差总和、峰值艏向角、平稳艏向角时间、平稳速度时间如表3所示。
表 3 不同Q值数据对比表Table 3. Data comparison table of different Q参数 峰值纵向
误差/m误差平方
总和/m2峰值
艏向角/rad平稳艏向角
时间/s平稳速度
时间/s{\boldsymbol{Q}}_{0} −2.80 591 1.032 59.1 122.7 {3\boldsymbol{Q}}_{0} −2.51 440 1.038 62.8 115.1 {6\boldsymbol{Q}}_{0} −2.20 305 1.046 65.3 109.2 9{\boldsymbol{Q}}_{0} −1.97 224 1.052 67.2 107.1 从表3和图12、图13可以看出，随着Q值增大，误差峰值和误差平方总和减小，平稳速度所需时间缩短，但峰值艏向角增大，平稳艏向角时间延长。因此，增加Q值可缓慢加快响应时间，但同时会增大航向稳定时间。
3） 固定 \boldsymbol{Q}={\boldsymbol{Q}}_{0},{\boldsymbol{S}=\boldsymbol{S}}_{0} ，R分别为 {\boldsymbol{R}}_{0},{\boldsymbol{R}}_{1}=2\;000\times \mathrm{diag}\left\{\mathrm{0.4,5},5\right\},{\boldsymbol{R}}_{2}=4\;000\times \mathrm{diag}\left\{\mathrm{0.4,5},5\right\} ，{\boldsymbol{R}}_{3}=8\;000\times \mathrm{diag}\left\{\mathrm{0.4,5},5\right\} ，进行4组实验。可得到每组实验的船舶轨迹、轨迹纵向误差和船舶状态，并进行对比，结果如图14和图15所示。
![]()
图 14 不同R值船舶轨迹与纵向误差对比图Figure 14. Comparison of ship trajectory and longitudinal error with different R values![]()
图 15 不同R值艏向角与速度变化对比图Figure 15. Comparison of heading angle and velocity changes with different R values在不同R值下得到的峰值纵向误差、纵向误差总和、峰值艏向角、平稳艏向角时间、平稳速度时间如表4所示。
表 4 不同R值数据对比表Table 4. Data comparison table of different R参数 峰值纵向
误差/m误差平方
总和/m2峰值
艏向角/rad平稳艏向角
时间/s平稳速度
时间/s{\boldsymbol{R}}_{1} −1.69 132 1.13 71.6 53.6 {\boldsymbol{R}}_{2} −2.34 334 1.07 63.8 88.6 {\boldsymbol{R}}_{0} −2.78 591 1.03 67.2 122.7 {\boldsymbol{R}}_{3} −3.11 887 1.00 77.3 150.0 从表4和图14、图15中可以看出，随着R值的增大，误差峰值和误差平方总和均增大，平稳速度所需要时间延长，而峰值艏向角减小，在 {\boldsymbol{R}}_{2} 后平稳艏向角时间延长，因此通过减小\boldsymbol{R} 的值会减少响应时间和航向稳定时间。
综上所述可总结为：通过增大权重矩阵 \boldsymbol{Q},\boldsymbol{S} 可加快系统响应速度，但若增大过多，会导致稳定时间延长；通过减小权重矩阵R，可加快系统响应速度，但减小过多，也会导致稳定时间延长。
为验证算法对其他参考轨迹的跟踪性能，设定参考轨迹为{x}_{\mathrm{d}}\left(t\right)=Ut ，{\boldsymbol{y}}_{\mathrm{d}}\left(t\right)=50\mathrm{sin}\left(Ut/100\right) ，且参考轨迹点的期望速度 {u}_{\mathrm{d}}=1\;\mathrm{m}/\mathrm{s} ，参数 \zeta =6\; 000 。为实现更快的系统响应，选取权重矩阵 \boldsymbol{Q},\boldsymbol{R},\boldsymbol{S} 分别为：
\left\{\begin{aligned} & {\boldsymbol{Q}}=3 \;000\times {\mathrm{diag}}\left\{\mathrm{1,1},\mathrm{1,1},\mathrm{1,1}\right\}\\& {\boldsymbol{R}}= 4\; 000\times {\mathrm{diag}}\left\{\mathrm{0.4,5},5\right\}\\& {\boldsymbol{S}}=1\; 000\times {\mathrm{diag}}\left\{\mathrm{1,1},\mathrm{1,1},\mathrm{1,1}\right\}\end{aligned} \right. (67) 1500 次迭代后，得到如图16所示的船舶航行轨迹（黑线是参考轨迹，红线为船舶实际运行轨迹），以及船舶与参考轨迹点的纵向误差图17。从图16和图17中可以看出，该算法的曲线跟踪效果较好，曲线最高点处的误差约为0.7 m，且平稳跟踪的纵向误差不超过0.7 m，处于十分精确的范围内。
为验证控制算法与辨识算法的联合性能，选取参数辨识实验周期数为
1500 ～2250 时的航迹为参考轨迹点，外界扰动量w同样选取周期数为1500 ～2250 的扰动估计值。权重矩阵 \boldsymbol{Q},\boldsymbol{R},\boldsymbol{S} 选取式（67）所示的形式。为跟踪参考轨迹点的每一时刻状态，矩阵 \boldsymbol{C}=\mathrm{diag}\left\{\mathrm{1,1},\mathrm{1,1},\mathrm{0,0}\right\} ，采用参数辨识算法和控制算法联合策略进行迭代。运行 150\;\mathrm{s} 后，可得到船舶航行轨迹如图18所示，其中黑线为参考轨迹，红线为船舶实际运行轨迹。从图18中可以看出，参考轨迹与实际运行轨迹在拐点处相差较大，其余部分相差较小。总体而言，该控制算法能实现轨迹跟踪，辨识算法与控制算法的联合策略具有可行性。
因此，本研究提出的船舶航迹跟踪控制算法无需使用者具备专业的船舶运动建模与控制背景知识，仅需设置加权矩阵 \boldsymbol{Q},\boldsymbol{R},\boldsymbol{S} 即可运行。此外，该方法具有广泛的环境适应性，与扰动的具体概率分布形式基本无关，且算法中的参数无需外部人为调整，能够达到“启动后不用管”的航迹跟踪控制状态。
5. 结 语
本文基于船舶仿真运动的数据，提出一种含扩张状态观测器的新息递推最小二乘算法。该辨识算法有效解决了船舶复杂动力系统导致的船舶模型参数未知问题，并能够实时估计外界干扰量。针对参数辨识稳定后的系统状态空间模型，设计了船舶航迹跟踪的有限时间状态调节器。该控制算法有效解决了船舶运动方程的非线性导致的难于优化控制难题，并获得了最优控制与船舶运动参数的解析关系。
为验证所提算法的有效性，搭建了双螺旋桨驱动的仿真无人艇，并进行了连续“Z实验”。结果表明，参数辨识算法估计出的状态量 \hat{\boldsymbol{\nu }}\left(t\right) 和扰动量 \hat{\boldsymbol{w}}\left(t\right) 具有一定的准确性。通过控制变量法对比分析不同权重矩阵 \boldsymbol{Q},\boldsymbol{R},\boldsymbol{S} 对系统性能的影响，结果表明，增大权重矩阵 \boldsymbol{Q}\mathrm{和}\boldsymbol{S} 或减小权重矩阵R可以加快系统响应速度，但过度增加 \boldsymbol{Q},\boldsymbol{S} 或减小R会导致稳定时间延长。
最后，通过曲线轨迹和仿真航迹跟踪实验，证明该系统响应较快，航迹跟踪较稳定。因此，本文所提出的方法能够在未知先验状态下较好地实现航迹跟踪控制。然而，该控制方法的极限性能尚待进一步研究，指标参数设计对控制效果的影响也仍需进一步探讨。
-
![]()
图 5 参数辨识与控制算法的流程图
Figure 5. Flow chart of parameter identification and control algorithm
![]()
图 10 不同S值船舶轨迹与纵向误差对比图
Figure 10. Comparison of ship trajectory and longitudinal error with different S values
![]()
图 11 不同S 值艏向角与速度变化对比图
Figure 11. Comparison of heading angle and velocity changes with different S values
![]()
图 12 不同Q值船舶轨迹与纵向误差对比图
Figure 12. Comparison of ship trajectory and longitudinal error with different Q values
![]()
图 13 不同Q值艏向角与速度变化对比图
Figure 13. Comparison of heading angle and velocity changes with different Q values
![]()
图 14 不同R值船舶轨迹与纵向误差对比图
Figure 14. Comparison of ship trajectory and longitudinal error with different R values
![]()
图 15 不同R值艏向角与速度变化对比图
Figure 15. Comparison of heading angle and velocity changes with different R values
表 1 辨识参数
Table 1 Identification parameters
序号 参数 数值 序号 参数 数值 1 {\boldsymbol{\theta }}_{11} 0.265 138 669 258 724 13 {\boldsymbol{\theta }}_{25} 0.305 785 017 250 565 2 {\boldsymbol{\theta }}_{12} −0.593 139 836 433 817 14 {\boldsymbol{\theta }}_{26} 0.817 640 200 880 232 3 {\boldsymbol{\theta }}_{13} 0.055 061 829 481 853 15 {\boldsymbol{\theta }}_{31} 0.061 005 725 705 921 4 {\boldsymbol{\theta }}_{14} 0.193 328 379 010 969 16 {\boldsymbol{\theta }}_{32} −0.247 108 873 113 790 5 {\boldsymbol{\theta }}_{15} 0.004 535 423 216 801 17 {\boldsymbol{\theta }}_{33} −0.071 435 501 716 979 6 {\boldsymbol{\theta }}_{16} 0.001 500 027 456 607 18 {\boldsymbol{\theta }}_{34} −0.067 124 874 659 099 7 {\boldsymbol{\theta }}_{17} −0.00 001 031 479 303 19 {\boldsymbol{\theta }}_{35} 0.150 310 762 377 552 8 {\boldsymbol{\theta }}_{18} −2.101 476 849 154 E-09 20 {\boldsymbol{\theta }}_{36} 0.156 863 665 033 070 9 {\boldsymbol{\theta }}_{21} −0.057 189 322 465 896 21 {\boldsymbol{\theta }}_{37} −0.000 528 418 571 230 10 {\boldsymbol{\theta }}_{22} −0.107 906 938 569 009 22 {\boldsymbol{\theta }}_{38} −0.001 535 487 343 785 11 {\boldsymbol{\theta }}_{23} −0.021 188 406 958 892 23 {\boldsymbol{\theta }}_{39} 7.379 676 517 343 E-06 12 {\boldsymbol{\theta }}_{24} −0.343 953 417 948 503 24 {\boldsymbol{\theta }}_{30} −6.976 693 955 860 E-07 
下载: 导出CSV
表 2 不同S数据对比
Table 2 Data comparison table of different S
参数 峰值纵向
误差/m误差平方
总和/m2峰值
艏向角/rad平稳艏向角
时间/s平稳速度
时间/s{\boldsymbol{S}}_{0} −2.80 591 1.03 67.2 122.7 {2\boldsymbol{S}}_{0} −2.12 246 1.09 70.9 71.4 {3\boldsymbol{S}}_{0} −1.78 151 1.13 73.4 52.4 4{\boldsymbol{S}}_{0} −1.57 114 1.15 75.7 42.0
下载: 导出CSV
表 3 不同Q值数据对比表
Table 3 Data comparison table of different Q
参数 峰值纵向
误差/m误差平方
总和/m2峰值
艏向角/rad平稳艏向角
时间/s平稳速度
时间/s{\boldsymbol{Q}}_{0} −2.80 591 1.032 59.1 122.7 {3\boldsymbol{Q}}_{0} −2.51 440 1.038 62.8 115.1 {6\boldsymbol{Q}}_{0} −2.20 305 1.046 65.3 109.2 9{\boldsymbol{Q}}_{0} −1.97 224 1.052 67.2 107.1
下载: 导出CSV
表 4 不同R值数据对比表
Table 4 Data comparison table of different R
参数 峰值纵向
误差/m误差平方
总和/m2峰值
艏向角/rad平稳艏向角
时间/s平稳速度
时间/s{\boldsymbol{R}}_{1} −1.69 132 1.13 71.6 53.6 {\boldsymbol{R}}_{2} −2.34 334 1.07 63.8 88.6 {\boldsymbol{R}}_{0} −2.78 591 1.03 67.2 122.7 {\boldsymbol{R}}_{3} −3.11 887 1.00 77.3 150.0
下载: 导出CSV
-
[1] WENG Y P, NAN D, WANG N, et al. Compound robust tracking control of disturbed quadrotor unmanned aerial vehicles: A data-driven cascade control approach[J]. Transactions of the Institute of Measurement & Control, 2022, 44(4): 941–951. doi: 10.1177/01423312211043675
[2] 郭杰, 刘轶华, 马利华. 基于多模态快速非奇异终端滑模的船舶航迹跟踪自抗扰控制[J]. 中国航海, 2020, 43(2): 7–13. doi: 10.3969/j.issn.1000-4653.2020.02.002 GUO J, LIU Y H, MA L H. Active disturbance rejection control for ship trajectory tracking with multimodal fast nonsingular terminal sliding mode strategy[J]. Navigation of China, 2020, 43(2): 7–13 (in Chinese). doi: 10.3969/j.issn.1000-4653.2020.02.002
[3] 宋利飞, 许传毅, 郝乐, 等. 基于改进DDPG算法的无人艇自适应控制[J/OL]. 中国舰船研究: 1−9 (2022-11-21) [2023-11-01]. https://dx.doi.org/10.19693/j.issn.1673-3185.03122. SONG L F, XU C Y, HAO L, et al. Adaptive control of unmanned surface vehicle based on improved DDPG algorithm[J/OL]. Chinese Journal of Ship Research: 1−9 (2022-11-21) [2023-11-01] (in Chinese). https://dx.doi.org/10.19693/j.issn.1673-3185.03122.
[4] 张子昌, 徐雪峰, 侯成刚. 基于MPC算法的AUV空间航迹跟踪控制[J]. 舰船科学技术, 2020, 42(12): 86–91. doi: 10.3404/j.issn.1672-7649.2020.12.017 ZHANG Z C, XU X F, HOU C G. AUV space trackway-pointfollowing control based on MPC[J]. Ship Science and Technology, 2020, 42(12): 86–91 (in Chinese). doi: 10.3404/j.issn.1672-7649.2020.12.017
[5] WANG N, GAO Y, ZHANG X F. Data-driven performance-prescribed reinforcement learning control of an unmanned surface vehicle[J]. IEEE Transactions on Neural Networks and Learning Systems, 2021, 32(12): 5456–5467. doi: 10.1109/TNNLS.2021.3056444
[6] NGUYEN M C, VUONG D P. The optimal control system of the ship based on the linear quadratic regular algorithm[J]. International Journal of Electrical and Computer Engineering, 2020, 10(5): 4562–4568. doi: 10.11591/ijece.v10i5.pp4562-4568
[7] 乐美龙. 船舶操纵性预报与港航操纵运动仿真[M]. 上海: 上海交通大学出版社, 2004. LE M L. Ship Manoeuvrability Prediction and simulation of port and waterway manoeuvring motion[M]. Shanghai: Shanghai Jiaotong University Press, 2004 (in Chinese).
[8] 朱曼, 文元桥, 孙吴强, 等. 一种基于扩展状态观测器的智能船舶Nomoto模型参数辨识方法[J]. 中国舰船研究, 2023, 18(3): 75–85. doi: 10.19693/j.issn.1673-3185.02552 ZHU M, WEN Y Q, SUN W Q, et al. Extended state observer-based parameter identification of Nomoto model for autonomous vessels[J]. Chinese Journal of Ship Research, 2023, 18(3): 75–85 (in Chinese). doi: 10.19693/j.issn.1673-3185.02552
[9] 赵大明, 施朝健, 彭静. 应用扩展卡尔曼滤波算法的船舶运动模型参数辨识[J]. 上海海事大学学报, 2008, 29(3): 5–9. doi: 10.3969/j.issn.1672-9498.2008.03.002 ZHAO D M, SHI C J, PENG J. Parameter identification to motion model of ship by extended Kalman filter[J]. Journal of Shanghai Maritime University, 2008, 29(3): 5–9 (in Chinese). doi: 10.3969/j.issn.1672-9498.2008.03.002
[10] 徐锋, 邹早建, 徐小卡, 等. 基于支持向量机的船舶操纵运动黑箱建模[J]. 北京航空航天大学学报, 2013, 39(11): 1553–1557. doi: 10.13700/j.bh.1001-5965.2013.11.026 XU F, ZOU Z J, XU X K, et al. Black-box modeling of ship manoeuvring motion based on support vector machines[J]. Journal of Beijing University of Aeronautics and Astronautics, 2013, 39(11): 1553–1557 (in Chinese). doi: 10.13700/j.bh.1001-5965.2013.11.026
[11] DAI S L, WANG C, LUO F. Identification and learning control of ocean surface ship using neural networks[J]. IEEE Transactions on Industrial Informatics, 2012, 8(4): 801–810. doi: 10.1109/TII.2012.2205584
[12] 陈中伟, 沈艳霞. 基于多新息最小二乘的感应电机参数辨识策略[J]. 江南大学学报(自然科学版), 2010, 9(5): 531–535. doi: 10.3969/j.issn.1671-7147.2010.05.006 CHEN Z W, SHEN Y X. Parameters identification of the induction motor by using multiple innovation least squares[J]. Journal of Jiangnan University (Natural Science Edition), 2010, 9(5): 531–535 (in Chinese). doi: 10.3969/j.issn.1671-7147.2010.05.006
[13] CB/Z 809-2016. 船舶操纵运动数学模型[S]. 中华人民共和国工业和信息化部, 2016.07. 11. CB/Z 809-2016. Mathematical model of ship manoeuvring motion [S]. Ministry of Industry and Information Technology of the People's Republic of China, July 11, 2016 (in Chinese).
[14] 贾欣乐, 杨盐生. 船舶运动数学模型: 机理建模与辨识建模[M]. 大连: 大连海事大学出版社, 1999. JIA X L, YANG Y S. Modeling of ship motion mathematical model mechanism and identification[M]. Dalian: Dalian Maritime University Press, 1999 (in Chinese).
[15] 杨炜帆. 水面无人艇的建模与运动特性仿真[D]. 大连: 大连海事大学, 2013. YANG W F. Unmanned surface vehicle modeling and motion simulation[D]. Dalian: Dalian Maritime University, 2013 (in Chinese).
[16] 董蛟, 刘忠, 张建强, 等. 基于干扰观测的欠驱动无人艇自适应航迹跟踪控制算法[J]. 系统工程与电子技术, 2019, 41(7): 1606–1616. doi: 10.3969/j.issn.1001 DONG J, LIU Z, ZHANG J Q, et al. Adaptive course tracking control of underactuated USV based on disturbance observation[J]. Systems Engineering and Electronics, 2019, 41(7): 1606–1616 (in Chinese). doi: 10.3969/j.issn.1001
[17] 谢朔. 基于多新息辨识理论的USV运动模型参数辨识[D]. 武汉: 武汉理工大学, 2017. XIE S. Motion model parameters identification of USV based on multi-innovation theory[D]. Wuhan: Wuhan University of Technology, 2017 (in Chinese).
[18] 丁峰. 系统辨识(6): 多新息辨识理论与方法[J]. 南京信息工程大学学报: 自然科学版, 2012, 4(1): 1–28. doi: 10.3969/j.issn.1674-7070.2012.01.001 DING F. System identification. Part F: Multi-innovation identification theory and methods[J]. Journal of Nanjing University of Information Science and Technology: Natural Science Edition, 2012, 4(1): 1–28 (in Chinese). doi: 10.3969/j.issn.1674-7070.2012.01.001
[19] 刘明亮, 朱江淼. 数字信号处理对电子测量与仪器的影响研究[J]. 电子测量与仪器学报, 2014, 28(10): 1041–1046. doi: 10.13382/j.jemi.2014.10.001 LIU M L, ZHU J M. Investigation in influence of digital signal processing to electronic measurement and instrument[J]. Journal of Electronic Measurement and Instrumentation, 2014, 28(10): 1041–1046 (in Chinese). doi: 10.13382/j.jemi.2014.10.001
[20] 钟万勰. 矩阵黎卡提方程的精细积分法[J]. 计算结构力学及其应用, 1994, 11(2): 113–119. ZHONG W X. The precise integration for matrix Riccati equation[J]. Computational Structural Mechanics and Applications, 1994, 11(2): 113–119 (in Chinese).
[21] 马丽华. 基于精细积分法的电力系统仿真研究[D]. 哈尔滨: 哈尔滨工业大学, 2007. MA L H. Research on simulation of power system based on precise integration algorithm[D]. Harbin: Harbin Institute of Technology, 2007 (in Chinese).
[22] 钟万勰. 线性二次最优控制的精细积分法[J]. 自动化学报, 2001, 27(2): 166–173. doi: 10.16383/j.aas.2001.02.004 ZHONG W X. The precise integration of LQ control problems[J]. Acta Automatica Sinica, 2001, 27(2): 166–173 (in Chinese). doi: 10.16383/j.aas.2001.02.004
计量
- 文章访问数: 358
- HTML全文浏览量: 20
- PDF下载量: 61
