ARTICLE · INTELLIGENCE

战地情报 · 详情页

来自尧图项目组的一线实战观察与深度解析

全驱动船舶轨迹跟踪的DSC+滑模+扰动观测器MATLAB仿真

全驱动船舶轨迹跟踪的DSC+滑模+扰动观测器MATLAB仿真

简介:面向无人船与水下机器人控制方向的学习者,这套MATLAB仿真资料聚焦全驱动船舶在模型不确定与外部扰动下的轨迹跟踪问题,将自适应动态面控制、滑模控制与扰动观测器结合,覆盖圆轨迹与直线轨迹两类仿真场景,适合具备一定非线性控制基础、希望复现并理解算法实现的研究生或工程人员。压缩包共20个文件,约56KB,以15个.m脚本为主体,另有2个txt文本、2个mdl模型及1个r2012a文件,脚本用于搭建船舶数学模型、设计控制器与观测器并输出仿真结果,模型与文本文件便于对照理解系统结构。已有1313人学习,读者可从中获得控制器设计思路、扰动估计与补偿的实现方式,以及轨迹跟踪仿真的参数调试与结果分析参考,适合作为课程作业、论文复现或课题预研的入门素材。

1. 从一条抖起来的跟踪曲线说起:这套方案要解决什么

调全驱动船舶轨迹跟踪的仿真,很多人卡在同一个地方:海况一加,位置误差曲线开始周期性鼓包;把滑模增益往上推,误差压下去了,舵角却开始高频跳变,控制力矩像噪声一样。问题不在增益大小,而在于把「模型不确定」「外部扰动」「执行器抖振」三件事挤在同一个滑模面上去消化。带扰动观测器的自适应动态面滑模控制(DSC + SMC + DOB)做的就是拆解:动态面负责把反步法的微分爆炸按住,扰动观测器把风浪流这些慢变项单独估出来前馈掉,自适应律只负责兜住残余上界,滑模面因此可以取小增益、低抖振。这套 MATLAB 方案面向做毕设、课设或交付验证的读者,从三自由度模型、控制器结构一路写到可运行的脚本和验收指标,拿到手能改参数、能出图、能答辩。

2. 被控对象与问题建模:三自由度船舶模型与扰动假设怎么定

建模这一层如果偷懒,后面所有整定都是玄学。全驱动船舶的「全驱动」不是形容词,而是指纵荡、横荡、艏摇三个自由度都有独立执行器(主推、侧推、舵桨组合),控制输入维度等于自由度维度,惯性矩阵可逆,控制律可以直接反解出来。欠驱动船只有纵向推力和艏摇力矩,横荡靠耦合,控制律要绕一大圈。这个差别决定了本方案能写成简洁的反馈线性化加鲁棒补偿结构。

2.1 三自由度运动学与动力学方程

北东坐标系下的位置与艏向记作 η = [x, y, ψ]ᵀ,船体坐标系下的纵荡、横荡、艏摇速度记作 ν = [u, v, r]ᵀ,两者之间靠旋转矩阵连接:

η̇ = J(ψ) ν M ν̇ + C(ν) ν + D(ν) ν = τ + d(t)

J(ψ) 是标准正交旋转阵,第三行第三列为 1,其余耦合项由 cosψ、sinψ 构成,正交性带来一个很实用的性质:J⁻¹ = Jᵀ,后面算期望速度时不用求逆。M 为含附加质量的惯性矩阵,对称正定;C(ν) 是科氏向心力矩阵,必须构造成反对称,即 C = -Cᵀ,这个性质在稳定性证明里用来消掉 νᵀCν 项;D(ν) 为阻尼矩阵,低速段通常取线性近似,高速段再补二次项。τ 是控制输入,d(t) 是风、浪、流合成的环境扰动。

参数取值直接决定仿真跑不跑得动。一套常用于验证的量级如下,注意这些是验证用参数,不是某条实船的数据:

符号含义典型取值设置要点
M₁₁纵向惯性(含附加质量)200附加质量约占刚体质量 10%~30%
M₂₂ / M₃₃横向惯性 / 转动惯量250 / 80横荡通常大于纵荡
M₂₃ = M₃₂横荡-艏摇耦合惯性30非对角项不要随手填 0
D₁₁ / D₂₂ / D₃₃线性阻尼系数70 / 100 / 50由阻尼比反推更稳
D₂₃ = D₃₂耦合阻尼20影响转弯时的侧滑

科氏矩阵不是随便填的,它由 M 的元素和 ν 共同决定,写错符号是后面观测器发散的头号原因:

% 由惯性矩阵元素构造反对称科氏矩阵,保证 nu'*C*nu = 0 function C = coriolis(nu, M) C = [ 0, 0, -M(2,2)*nu(2) - M(2,3)*nu(3); 0, 0, M(1,1)*nu(1); M(2,2)*nu(2) + M(2,3)*nu(3), -M(1,1)*nu(1), 0 ]; end

这段的三处非零元素互为镜像,取负号的位置必须成对出现。写完之后建议加一行断言assert(norm(nu'*C*nu) < 1e-10)自检,反对称性不成立就说明符号抄错了。别小看这一行,它能省掉半天的排查时间。

2.2 环境扰动怎么描述才算合理

扰动模型决定了观测器能估到什么程度。工程上常用三段叠加:常值分量对应定常风与均匀流,低频正弦对应波浪的诱导漂移,慢变随机项对应阵风。写成时域表达式更直观:

d₁(t) = 1.2 + 0.8·sin(0.30t) d₂(t) = -0.9 + 0.6·sin(0.22t + 0.5) d₃(t) = 0.4 + 0.3·sin(0.18t + 1.1)

频率都压在 0.3 rad/s 以下,对应周期 20 秒以上,这是有意的——真实海况的能量集中在低频段,观测器带宽只要高于这个频段就能跟得住。如果你把扰动频率设到 5 rad/s 以上,再用带宽 8 的观测器去估,估不出来不是算法的问题,是假设不成立。

设计控制器时需要三条前提:d(t) 有界;ḋ(t) 有界(慢变假设);执行器不饱和或饱和时间足够短。前两条是观测器误差一致最终有界的充分条件,第三条是闭环稳定性的隐含条件。答辩时被问到「你的方法对扰动有什么要求」,把这三条讲清楚,比背公式管用。

2.3 误差坐标与控制目标

位置误差取 z₁ = η - η_d,艏向分量必须折算到 (-π, π],否则船从 179° 转到 -179° 时误差会跳到 358°,控制器会命令船原地转一整圈。速度层不直接对 ν - ν_d 做反馈,而是对 ν - α_1d,其中 α_1d 是虚拟控制律 α_1 经过一阶低通滤波后的结果,这个下标 d 是 dynamic surface 的来源。

跟踪目标定成分层的形式更好验收:瞬态阶段要求误差单调收敛、超调小于初始误差的 20%;稳态阶段要求位置误差 RMS 小于 1 m,艏向误差小于 0.02 rad。这两个数字不是拍脑袋,1 m 相对 60 m 量级的回转半径约 1.7%,已经优于多数港内作业的定位精度要求;0.02 rad 约 1.15°,对应舵角分辨率的合理量级。

参考轨迹用匀速圆周最省事,也最容易暴露问题:半径 60 m、角速度 0.05 rad/s,艏向沿切线方向,期望速度可由 Jᵀ 直接反解。相比直线轨迹,圆周轨迹在整个仿真周期内都保持恒定的向心加速度需求,一旦控制器增益不够或观测器滞后,误差曲线会立刻鼓起来,问题藏不住。

3. 控制器结构:动态面、滑模与扰动观测器的三层分工

三个模块各管一段,边界划清楚才不会互相打架。动态面处理的是「虚拟控制律求导太难」,滑模处理的是「模型不精确时误差怎么压住」,扰动观测器处理的是「慢变外扰不要用高增益硬扛」。把它们串起来的顺序是:运动学层给出虚拟速度指令,动态面滤波得到光滑的导数估计,动力学层在滤波后的指令上做滑模跟踪,观测器把扰动估计前馈进控制律。

3.1 反步法为什么需要动态面

标准反步法的推导路径是:针对 z₁ 设计 Lyapunov 函数,取虚拟控制 α₁ = Jᵀ(-c₁z₁) + ν_d,再对 z₂ = ν - α₁ 设计第二步。问题出在对 α₁ 求导——它含有 J(ψ) 的导数、η_d 的二阶导、以及 ν_d 的导数,展开后项数会翻好几倍,而且需要知道参考轨迹的加速度,实际系统里这个量要么测不到,要么噪声极大。

动态面的做法是把解析求导换成物理滤波:让 α₁ 通过一个时间常数为 τ_f 的一阶低通环节,得到 α_1d,用 (α₁ - α_1d)/τ_f 近似 α̇_1d。代价是引入滤波误差 y₂ = α_1d - α₁,需要在 Lyapunov 函数里加一项 0.5y₂ᵀy₂ 来补偿。这笔交易非常划算:求导爆炸没了,代价只是一项可以用 τ_f 调小来压缩的误差。

滤波器离散实现只有两行,但有个坑:必须用前向欧拉且步长与仿真步长一致,不能用连续系统的解析解代替,否则滤波器的实际带宽和你以为的不一样:

% 一阶低通滤波:alpha 为虚拟控制律,alpha_f 为滤波后的指令 alpha_dot = (alpha - alpha_f) / ctrl.tf; % 用滤波前后差值近似导数 alpha_f = alpha_f + dt * alpha_dot; % 前向欧拉,dt 与仿真步长严格一致

τ_f 的选法有讲究。取太小,滤波器带宽过高,等价于退回反步法,求导噪声全部透传;取太大,α_1d 跟不上 α₁,跟踪出现明显滞后,稳态误差压不下去。经验区间是闭环带宽的 1/10 到 1/5,本方案闭环带宽约 2 rad/s,取 τ_f = 0.05 s 对应带宽 20 rad/s,实测既不抖也不拖。

3.2 非线性扰动观测器:把慢变外扰单独拎出来

观测器的思路是构造一个辅助状态,让估计误差自己收敛,而不是靠控制器的鲁棒项去硬吃。常用形式是:

d̂ = ζ + K_o M ν ζ̇ = -K_o ζ - K_o ( -C(ν)ν - D(ν)ν + τ + K_o M ν )

把两式对时间求导整理,可以得到估计误差的动态方程 d̃̇ = ḋ - K_o d̃。取 Lyapunov 函数 V = 0.5d̃ᵀd̃,求导得 V̇ ≤ -λ_min(K_o)‖d̃‖² + ‖d̃‖·‖ḋ‖。当 ‖d̃‖ > ‖ḋ‖_max / λ_min(K_o) 时 V̇ < 0,误差被拉进一个球内,球半径与扰动变化率成正比、与观测器增益成反比。这个结论说明两件事:K_o 越大稳态估计越准;但 K_o 同时也放大测量噪声,因为观测器输入是 ν 和 τ,任何一个带噪都会被乘以 K_o。

整定上,K_o 取闭环带宽的 3 到 5 倍是稳妥区间。本方案闭环带宽约 2 rad/s,取 K_o = 8 对应 8 rad/s,既能跟住 0.3 rad/s 的扰动频段,又不会把数值噪声放大到控制量里。如果发现 d̂ 曲线毛刺很重,先别降 K_o,检查一下 ν 是不是直接用了带噪信号,加个一阶滤波再进观测器往往更有效。

上一节提到的tau_prev也是这里的细节:观测器需要当前的 τ,而 τ 又要用 d̂,形成代数环。工程上的做法是用上一采样周期的控制量,延迟一个步长对 0.01 s 的步长而言完全可以接受,稳定性分析里把它归入有界扰动项处理。

3.3 自适应增益与滑模面的耦合关系

滑模面取 s = ν - α_1d,趋近律采用等效控制加切换项的结构。如果把 τ 展开写成:

τ = C(ν)ν + D(ν)ν + M α̇_1d - d̂ - M K₁ s - (M K₂ + k̂)·sat(s/ε)

前三项是对消已知非线性,-d̂ 是对消估计出的外扰,-M K₁ s 是线性反馈负责指数收敛,最后一项是切换项负责兜住残余不确定。k̂ 是自适应切换增益,按 k̂̇ = Γ|s| - σk̂ 更新,Γ 为自适应速率,σ 为泄漏系数。

泄漏项不是可有可无的装饰。|s| 在稳态时接近零但不会严格等于零,如果没有 -σk̂ 这一项,k̂ 会缓慢单调增长,跑几千秒后切换增益大到离谱,控制量随之饱和。加上泄漏项之后,k̂ 的稳态值收敛到 |s|·Γ/σ 附近的平衡点,既有界又能自适应调整。

饱和函数 sat(s/ε) = s/(|s| + ε) 替换符号函数,是把抖振从源头掐掉的关键。ε 取得太小等于没换,取得太大滑模面附近的等效增益下降、稳态误差上升。经验上 ε 取 0.01 到 0.05 之间,与采样步长同量级偏大一点,本方案取 0.02。

把上述所有项装进一个 Lyapunov 函数做整体分析:V = 0.5z₁ᵀz₁ + 0.5sᵀMs + 0.5d̃ᵀd̃ + 0.5tr(k̃ᵀΓ⁻¹k̃) + 0.5y₂ᵀy₂,逐项求导后,z₁ 项由 -c₁‖z₁‖² 压制,s 项由 -λ_min(K₁)‖s‖² 压制,d̃ 项由观测器增益压制,k̃ 项由自适应律压制,y₂ 项需要处理 α̇₁ 的有界性——这是全场唯一需要假设的地方,通常假设 α̇₁ 存在上界 M₂,用 -y₂²/τ_f + y₂M₂ 配成完全平方放缩。结论是闭环所有信号一致最终有界,界的大小与 M₂、‖ḋ‖_max、ε 成正比。答辩时这条链条能一口气讲下来,基本就稳了。

4. MATLAB 落地实现:从模型脚本到仿真主循环

代码结构按「一个主脚本 + 四个函数文件」组织,主脚本负责时间推进和记录,函数各管一摊。这样拆的好处是调参只改主脚本顶部的参数区,不用翻进函数里;换被控对象只改 ship_dyn,控制器一行不动。

4.1 文件组织与主仿真脚本

目录建议这样分:main_sim.m主脚本,ship_dyn.m船舶动力学,rk4_step.m积分器,controller.m控制律与观测器,plot_results.m出图与指标统计。全部放在同一目录即可,不依赖任何工具箱——艏向折算自己写 wrap 函数,别用 wrapToPi,那个函数属于 Robotics System Toolbox,换台机器就跑不起来。

主脚本顶部的参数区是所有整定的入口:

%% 主仿真脚本:全驱动船舶轨迹跟踪(DSC + 滑模 + 扰动观测器) clear; clc; close all; dt = 0.01; T = 120; t = 0:dt:T; N = numel(t); % ---- 被控对象参数 ---- M = [200 0 0; 0 250 30; 0 30 80]; % 对称正定,含附加质量 D = [ 70 0 0; 0 100 20; 0 20 50]; % 线性阻尼 param.M = M; param.D = D; param.Minv = inv(M); % ---- 控制器参数 ---- ctrl.c1 = diag([0.40, 0.40, 0.60]); % 位置误差反馈增益 ctrl.tf = 0.05; % 动态面滤波时间常数 ctrl.K1 = diag([30, 30, 20]); % 滑模线性增益 ctrl.K2 = diag([ 2, 2, 1]); % 切换增益基值 ctrl.eps = 0.02; % 饱和函数边界层 ctrl.gamma = diag([0.8, 0.8, 0.6]); % 自适应速率 ctrl.sigma = 0.05; % 泄漏系数,防增益无界 ctrl.Ko = diag([8, 8, 8]); % 观测器增益 % ---- 初始状态 ---- eta = [55; 10; 0]; nu = [0; 0; 0]; alpha_f = zeros(3,1); zeta = zeros(3,1); dhat = zeros(3,1); khat = [0;0;0]; tau_prev = zeros(3,1);

初始位置刻意放在离参考轨迹起点 55 m 的地方,这不是随便定的——初始误差足够大,才能看出瞬态阶段有没有超调、观测器要多久收敛、自适应增益会不会冲高。如果一上来就贴着轨迹跑,所有曲线都是平的,什么也验证不出来。

主循环里记录量和状态推进是分开的,顺序不能乱:先算控制,再记录,最后积分推进。如果先积分再算控制,相当于控制器晚了一步,dt 越小影响越小,但曲线上的相位滞后是真实存在的:

for k = 1:N tk = t(k); [eta_d, nu_d] = ref_traj(tk); % 参考轨迹与期望速度 d = disturbance(tk); % 扰动真值,仅用于记录与对比 [tau, alpha_f, dhat, zeta, khat, tau_prev] = controller( ... eta, nu, eta_d, nu_d, alpha_f, dhat, zeta, khat, ... tau_prev, ctrl, param, dt); % 记录 log.eta(:,k)=eta; log.eta_d(:,k)=eta_d; log.nu(:,k)=nu; log.tau(:,k)=tau; log.d(:,k)=d; log.dhat(:,k)=dhat; log.khat(:,k)=khat; log.s(:,k)=nu-alpha_f; % 对象推进 xnew = rk4_step(@(tt,xx) ship_dyn(tt, xx, tau, d, param), tk, [eta;nu], dt); eta = xnew(1:3); nu = xnew(4:6); end plot_results(t, log, ctrl);

参考轨迹函数把圆周和艏向一次算清,注意 J 用 Jᵀ 求逆那一行:

function [eta_d, nu_d] = ref_traj(tk) R = 60; w = 0.05; psi_d = w*tk + pi/2; % 艏向沿切线 eta_d = [R*cos(w*tk); R*sin(w*tk); psi_d]; eta_d_dot = [-R*w*sin(w*tk); R*w*cos(w*tk); w]; Jd = [cos(psi_d), -sin(psi_d), 0; sin(psi_d), cos(psi_d), 0; 0, 0, 1]; nu_d = Jd' * eta_d_dot; % 正交阵逆等于转置 end

4.2 被控对象与积分器封装

船舶动力学函数只做一件事:给定当前状态和控制量,算导数。科氏矩阵调用前面写好的 coriolis,别在主脚本里内联展开,否则换参数时要改两处:

function dx = ship_dyn(~, x, tau, d, p) eta = x(1:3); nu = x(4:6); psi = eta(3); J = [cos(psi), -sin(psi), 0; sin(psi), cos(psi), 0; 0, 0, 1]; C = coriolis(nu, p.M); dnu = p.Minv * (tau + d - C*nu - p.D*nu); % 解出加速度 dx = [J*nu; dnu]; end

积分器用四阶龙格库塔,不用 ode45 的原因是:ode45 是变步长,控制器的离散滤波依赖固定 dt,变步长会让滤波器的实际时间常数在每个周期都不一致,跑出来的曲线不可重复。固定步长 + 自己写 RK4,好处是结果逐位可复现,答辩演示时不会出现「昨天跑的图今天对不上」的尴尬:

function x_new = rk4_step(f, tk, x, dt) k1 = f(tk, x); k2 = f(tk + dt/2, x + dt/2*k1); k3 = f(tk + dt/2, x + dt/2*k2); k4 = f(tk + dt, x + dt*k3); x_new = x + dt/6*(k1 + 2*k2 + 2*k3 + k4); end

步长选择有一条硬约束:dt 必须明显小于动态面滤波时间常数 τ_f,通常要求 dt ≤ τ_f/20。τ_f = 0.05 对应 dt ≤ 0.0025,本方案取 dt = 0.01 其实是偏大的,但因为滤波器用前向欧拉、被控对象用 RK4,实测仍然稳定。保守起见,如果要把 τ_f 调到 0.02,dt 得同步降到 0.001。

4.3 控制器函数:一次把虚拟控制、观测、滑模算完

控制器函数是整个方案的核心,按顺序六步走:算误差、算虚拟控制、滤波、算滑模面、更新观测器、算控制量。顺序错了结果就是错的,比如先算滑模面再滤波,s 里用的就是旧指令:

function [tau, alpha_f, dhat, zeta, khat, tau_prev] = controller( ... eta, nu, eta_d, nu_d, alpha_f, dhat, zeta, khat, ... tau_prev, ctrl, param, dt) % 1) 位置误差,艏向折算到 (-pi, pi] z1 = eta - eta_d; z1(3) = mod(z1(3) + pi, 2*pi) - pi; % 2) 运动学层虚拟控制律 psi = eta(3); J = [cos(psi), -sin(psi), 0; sin(psi), cos(psi), 0; 0, 0, 1]; alpha = J' * (-ctrl.c1 * z1) + nu_d; % 3) 动态面一阶滤波 alpha_dot = (alpha - alpha_f) / ctrl.tf; alpha_f = alpha_f + dt * alpha_dot; % 4) 滑模面 s = nu - alpha_f; % 5) 扰动观测器(用上一拍控制量打断代数环) C = coriolis(nu, param.M); dhat = zeta + ctrl.Ko * param.M * nu; zeta_dot = -ctrl.Ko*zeta ... - ctrl.Ko*( -C*nu - param.D*nu + tau_prev + ctrl.Ko*param.M*nu ); zeta = zeta + dt * zeta_dot; % 6) 滑模控制律 + 自适应切换增益 sat_s = s ./ (abs(s) + ctrl.eps); % 饱和函数替代 sign k2v = diag(ctrl.K2); tau = C*nu + param.D*nu + param.M*alpha_dot - dhat ... - param.M*ctrl.K1*s - (k2v + khat).*sat_s; khat = khat + dt*( ctrl.gamma*abs(s) - ctrl.sigma*khat ); khat = max(khat, 0); % 增益不允许为负 tau_prev = tau; % 传给下一拍观测器 end

几个参数的作用值得单独说清。ctrl.c1决定位置误差收敛快慢,三个分量分别对应纵、横、艏向,艏向取 0.6 略大是因为角度误差对推力需求更敏感。ctrl.K1是滑模面的线性反馈增益,取值受执行器饱和限制,30 对应时间常数约 0.03 s,已经接近采样周期的三倍,再大就会出现数值抖振。ctrl.gamma与ctrl.sigma的比值决定自适应增益的稳态值,比值越大 k̂ 稳态越高、鲁棒性越强但抖振越明显,0.8/0.05 = 16 是个折中。

最后一行khat = max(khat, 0)是保护,不是装饰。当 s 长期处于零附近且有噪声时,泄漏项可能把 k̂ 拉到负值,负的切换增益等于给系统注入正反馈,一跑就发散。这行代码拦住的就是这种情况。

4.4 出图与指标统计

四张图是这套方案的标准交付物:轨迹对比、三轴误差、扰动估计对比、控制输入。缺任何一张,验收时都会被打回来。轨迹图必须axis equal,否则圆会画成椭圆,视觉上判断不了跟踪形状。

function plot_results(t, log, ctrl) e = log.eta - log.eta_d; e(3,:) = mod(e(3,:) + pi, 2*pi) - pi; % 艏向误差统一折算 figure('Name','轨迹跟踪'); plot(log.eta_d(1,:), log.eta_d(2,:), 'k--'); hold on; plot(log.eta(1,:), log.eta(2,:), 'b', 'LineWidth', 1.2); axis equal; grid on; legend('期望','实际','Location','best'); figure('Name','跟踪误差'); subplot(3,1,1); plot(t, e(1,:)); ylabel('e_x / m'); grid on; subplot(3,1,2); plot(t, e(2,:)); ylabel('e_y / m'); grid on; subplot(3,1,3); plot(t, e(3,:)); ylabel('e_\psi / rad'); xlabel('t / s'); grid on; figure('Name','扰动估计'); subplot(3,1,1); plot(t, log.d(1,:), 'k', t, log.dhat(1,:), 'r--'); ylabel('d_1'); grid on; % 稳态段指标:取后 30% 时间窗 idx = round(0.7*numel(t)) : numel(t); rms_e = sqrt(mean(e(:,idx).^2, 2)); fprintf('稳态RMS: ex=%.3f m, ey=%.3f m, epsi=%.4f rad\n', rms_e); fprintf('峰值力矩: %.1f / %.1f / %.1f\n', max(abs(log.tau(:,idx)),[],2)); end

指标统计只取后 30% 时间窗,是因为前段包含初始误差收敛过程,把瞬态也算进 RMS 会让数字虚高,掩盖稳态性能的真实水平。这个细节在答辩时会被问,主动说明比被动解释好。

5. 避坑与排查:五类最容易翻车的现场

调试阶段踩的坑基本集中在五个地方,按出现频率排序,每条都给出「现象 → 原因 → 解决」的路径,照着对号入座能省掉大量时间。

5.1 滤波器时间常数与控制增益的整定顺序搞反

现象:位置误差能收敛,但控制量在稳态段有持续的周期性振荡,周期与参考轨迹的圆周周期一致。原因:先调了 c₁ 和 K₁,最后才随手给 τ_f 一个值,导致滤波器带宽低于滑模环带宽,α_1d 跟不上 α₁,滑模面 s 里混进了滤波滞后引起的周期性成分,控制器把这个滞后当成误差去追。解决办法是把整定顺序倒过来:先定 τ_f,再定 c₁,最后定 K₁。具体做法是先让 K₁ 和 K₂ 置零、只留运动学层,调 c₁ 直到位置误差以期望速率收敛;然后固定 c₁ 打开滤波,从 τ_f = 0.1 往下试,找到误差曲线刚开始出现高频纹路的上一个值;最后打开滑模层调 K₁。

5.2 符号函数直接怼进舵角,抖振压不住

现象:控制力矩曲线呈方波状高频跳变,步长减半后跳变更严重。原因:用了 sign(s) 而不是饱和函数,s 在零附近来回穿越时控制量在两个极值之间瞬间切换,离散仿真里表现为逐拍翻转。这个坑的血泪经验是:只要执行器是物理舵机或推进器,符号函数就不该出现在最终版代码里。解决分三步:把 sign 换成 sat(s/ε);ε 从 0.01 起试,观察稳态误差和抖振的折中;如果仍有残差,把 s 本身过一道一阶滤波再进饱和函数,代价是牺牲一点响应速度。

5.3 观测器坐标系与对象不一致

现象:扰动估计 d̂ 不收敛到真值,而是收敛到一个明显偏移的量,或者直接发散。原因:观测器公式里的 ν 必须取船体坐标系速度,如果误用了北东坐标系下的 η̇,或者科氏矩阵的符号写反(C 应该是 -Cᵀ 那个方向),观测器内部的状态方程与被控对象就不匹配,误差动态 d̃̇ = ḋ - K_o d̃ 的推导前提直接崩掉。解决:先把控制量置零、扰动置常值,单独跑观测器,看 d̂ 能不能收敛到常值;这一步过了再接控制器。另外务必检查 coriolis 函数的反对称性,用nu'*C*nu是否为零来判定,比肉眼核对快得多。

5.4 积分步长与滤波器带宽失配

现象:把 dt 从 0.01 改到 0.005 想让曲线更光滑,结果反而发散了。原因:前向欧拉离散化一阶滤波器有一阶稳定条件 dt < 2τ_f,但实际工程中还要留裕度,真正可用的是 dt ≤ τ_f/20。dt 减半的同时如果忘了同步调整 τ_f,或者反过来把 τ_f 调得很小却没减 dt,滤波器就会进入数值不稳定区。解决办法是固定一条规则:任何一次修改 τ_f 或 dt,都用 dt/τ_f ≤ 0.05 这个比值校验一遍。仿真前打印一行比值,养成习惯就不会忘。

5.5 参考轨迹本身不可跟踪

现象:位置误差始终收敛不到零,稳态有一个固定偏差,增大任何增益都无效。原因:参考轨迹的加速度需求超出了执行器能力,或者艏向指令与轨迹切线方向不匹配,导致横荡方向需要持续输出侧向力,而全驱动船的侧推能力通常远小于主推。解决:先做可达性检查,把参考轨迹代入模型反解所需的 τ,看是否超出饱和限;超了就降低轨迹速度或曲率。另外艏向指令最好直接取轨迹切线方向 atan2(ẏ_d, ẋ_d),人为设定一个与轨迹无关的 ψ_d 会让船「横着走」,误差自然收不住。

6. 验收与答辩:四张曲线、三个指标怎么摆出来

验收这件事,核心不是证明曲线好看,而是证明「这套控制律在什么条件下、以什么代价、达到什么精度」。把这三件事说清楚,评审基本不会再追问。建议按下面的表格准备,每一项都要能当场从仿真结果里调出来:

验收项定义方式合格线对应图
稳态位置精度后 30% 时间窗 RMS≤ 1.0 m误差曲线
稳态艏向精度后 30% 时间窗 RMS≤ 0.02 rad误差曲线
扰动估计残差‖d - d̂‖ 的稳态均值≤ 真值幅值的 15%扰动估计图
控制量平滑度τ 的差分绝对值最大值不出现逐拍翻转控制输入图
瞬态超调误差峰值与初值之比≤ 20%误差曲线前段

指标统计代码可以直接在主脚本末尾加一段,注意逻辑索引和维度对齐:

e = log.eta - log.eta_d; e(3,:) = mod(e(3,:) + pi, 2*pi) - pi; idx = round(0.7*N) : N; % 稳态窗口 rms_pos = sqrt(mean(sum(e(1:2,idx).^2, 1))); % 位置合成误差 rms_psi = sqrt(mean(e(3,idx).^2)); res_d = mean(sqrt(sum((log.d(:,idx)-log.dhat(:,idx)).^2, 1))); smooth = max(max(abs(diff(log.tau(:,idx), 1, 2)))); fprintf('位置RMS %.3f m | 艏向RMS %.4f rad | 扰动残差 %.3f | 力矩跳变 %.2f\n', ... rms_pos, rms_psi, res_d, smooth);

答辩环节最常被追问的三个问题,提前准备好答案。第一个是「为什么不用纯反步法」,回答要点是虚拟控制律求导需要参考轨迹二阶导,实船上这个量不可靠,动态面用滤波替代解析求导,代价是引入了需要补偿的滤波误差。第二个是「观测器需要知道扰动导数上界吗」,回答是不需要知道具体值,只需要它存在且有界,估计误差的收敛球半径与这个上界成正比,工程上通过加大 K_o 来压缩。第三个是「自适应增益会不会无限增长」,回答是有泄漏项时 k̂ 的稳态值被限制在 Γ|s|/σ 附近,加上非负截断,不会发散。

还有一个容易被忽略的演示技巧:把扰动估计图和真实扰动画在同一张子图上,用实线画真值、虚线画估计,两条线贴合的程度比任何文字说明都有说服力。如果残差在某个频段偏大,主动指出这是观测器带宽的固有局限,并给出「提高 K_o 的代价是噪声放大」的权衡说明,比等评审问出来强。

从那以后我每次交这类方案,都会强制走一遍「先关扰动、再开扰动、最后加噪声」的三段式验证:第一段确认标称模型下能收敛,第二段确认观测器能跟上慢变扰动,第三段确认测量噪声进来后控制量不失控。三关都过,再去调参优化指标,顺序

本文还有配套的精品资源,点击获取

RELATED READING

延伸阅读

更多一线实战笔记与深度复盘,助您持续精进