简介面向船舶工程与控制工程学习者的动力定位仿真项目围绕MATLAB与Simulink搭建船舶动力定位模拟器可用于课程实验、控制策略验证与船员操作认知训练。内容从船舶六自由度运动模型出发涵盖纵向、横向、垂向及摇荡运动的状态空间描述并结合控制器、海浪干扰与白噪声设置借助示波器观察位移、速度和控制反馈帮助理解定位精度、响应时间与稳定性之间的关系。压缩包共13个文件以mat数据、m脚本和mdl、slx Simulink模型为主另有少量工程缓存与配置文件整体约58KB便于在MATLAB环境中直接打开和拆解学习。已有1305人学习适合希望把动力定位理论落到仿真模型上的读者可据此复现船舶运动与被控响应测试PID等控制方案在不同海况扰动下的表现并进一步优化参数。1. 船舶动力定位模拟器一条不抛锚的船怎么在 Simulink 里稳住把一条船钉在海面上某个坐标点艏向还要卡在 ±0.5° 以内这件事靠的不是锚是推力。船舶动力定位DP做的就是这件事位置和艏向给定期望值控制器算出需要的合力与合力矩推力分配器再把它摊到几个推进器上风、浪、流一边在推你的船系统一边在往回拽。MATLAB/Simulink 版本的模拟器把这条闭环完整搬进了桌面不用等海上试验窗口也不用碰水池。这套东西最适合三类人做船舶运动控制课题的学生、做 DP 系统联调和上位机逻辑验证的工程师以及要给操控人员做演示培训的人。它不能替代实船试验但它能在你上船之前把控制器参数、推力饱和、模式切换这些容易翻车的地方先跑一遍。下面按拆包看骨架、跑通第一次闭环、复现三个典型工况、排错的顺序来。2. 拆开包先看骨架三自由度运动方程和 Simulink 分层的对应关系这类仿真包的目录结构大多来自同一个习惯物理对象、环境扰动、控制器、分配器各占一个子系统顶层只做连线和信号路由。先弄清这一层后面改参数才知道该动哪个文件否则拿到包直接双击 slx面对一堆 From/Goto 标签很容易半小时都找不到船体模型在哪。2.1 目录里每个文件负责什么解压后别急着打开模型先看层级。常见的组织方式是这样# 解压后先看层级不要直接双击 slx unzip -q dp_simulator.zip -d dp_simulator find dp_simulator -maxdepth 2 -type f | sort典型输出会长成这样文件名各家不同分层逻辑基本一致dp_simulator/ ├── dp_main.slx % 顶层闭环参考 - 控制器 - 分配 - 船体 - 环境 ├── init/ │ ├── dp_init.m % 参数唯一入口返回一个 struct │ └── vessel_data.m % 船型主尺度与推进器布置 ├── sub/ │ ├── vessel3dof.slx % 三自由度船舶运动模型 │ ├── env_load.slx % 风、流、二阶波浪漂移 │ ├── controller.slx % DP 控制器PID 或 LQR │ └── allocation.slx % 推力分配QP ├── scripts/ │ ├── run_capability.m % 能力图批量扫描 │ └── plot_result.m % 结果后处理 └── data/看清楚三件事就够了参数从哪来init 目录、船体模型在哪sub 目录、顶层模型的采样率是多少。第三点最容易被忽略顶层如果统一用连续求解器控制器和船体就被迫用同一个步长后面想做 0.1 s 的控制器采样会很别扭。提示先把整个包复制一份再改尤其是 slx。Simulink 模型的自动保存会覆盖原文件改崩了没有后悔药。2.2 把 M、D、J 写成 MATLAB Function为什么先做船体再做控制器动力定位工作在零速到低速区间横荡、纵荡、艏摇三个自由度足够描述位置保持能力横摇纵摇垂荡在这个场景下不参与控制回路。所以模型是标准的三自由度形式运动学η̇ R(ψ)·ν动力学M·ν̇ D·ν τ τ_env其中η [x; y; ψ]是北东坐标系下的位置和艏向ν [u; v; r]是船体坐标系下的速度τ是推进器合力τ_env是环境载荷。M是刚体质量矩阵加附加质量D是阻尼矩阵。我一般的顺序是先把船体单独跑开环给一个恒定的纵向力看速度是不是收敛到一个稳定值再给一个恒定的艏摇力矩看角速度是不是也被阻尼吃掉。这一步过了才敢接控制器。反过来先调控制器一旦发散你根本分不清是增益错了还是船体矩阵符号错了。function [eta_dot, nu_dot] vessel3dof(eta, nu, tau) %VESSEL3DOF 三自由度纵荡-横荡-艏摇船舶模型 % 输入: % eta [x; y; psi] 北东坐标系位置与艏向 [m; m; rad] % nu [u; v; r] 船体坐标系速度 [m/s; m/s; rad/s] % tau [X; Y; N] 合力与合力矩 [N; N; N*m] % 输出: % eta_dot, nu_dot 状态导数 % 惯性矩阵 M M_RB M_A已含附加质量 M [5.3122e6, 0, 0; 0, 8.2831e6, 0; 0, 0, 3.7454e9]; % 线性阻尼矩阵 D低速区线性化含艏摇-横荡耦合项 D [5.0242e4, 0, 0; 0, 2.7229e5, -4.3933e6; 0, -4.3933e6, 4.1894e8]; psi eta(3); R [cos(psi), -sin(psi), 0; sin(psi), cos(psi), 0; 0, 0, 1]; eta_dot R * nu; % 船体系速度 - 北东系位置变化率 nu_dot M \ (tau - D * nu); % 左除不要写 inv(M)*() end几个参数的含义和取舍要说清楚。M(1,1)对应纵荡方向的总惯性量级是排水量的 1.1 到 1.3 倍M(3,3)是艏摇惯量通常按m·L²/12估算再加附加量。D里D(2,3)和D(3,2)这一对非对角项是横荡-艏摇耦合去掉它船在受横向力时会表现得“太干净”和真实船舶的低速响应差得比较远做能力图扫描时结果会偏乐观。代码里用了M \ (...)而不是inv(M)*(...)。这不是洁癖M的三个对角元差了三个数量级显式求逆在长时间仿真里会放大舍入误差跑到几十分钟后位置开始缓慢漂移排查起来非常费劲。如果你手上的包保留了科里奥利项C(ν)·ν也不用删低速下它的贡献很小但要确认符号C(ν)是反对称的符号写反会往系统里注入能量表现为速度越跑越大。这就是典型的玄学问题看曲线完全看不出原因。2.3 风、浪、流三块载荷建模口径不同力不能直接相加三种环境力的物理性质完全不一样混在一起写是常见错误。风的本质是相对风速产生的气动力必须用船速去修正风速船的航速会改变相对风向角相对风向角又决定气动力系数。流不是“加一个力”它改变的是阻尼项里的相对速度正确做法是把D·ν换成D·(ν - ν_c)。波浪分两部分一阶波浪力是周期性振荡幅值不小但平均值为零控制器不该去追它测量里要滤掉二阶波浪漂移力是慢变的平均值不为零是真正需要抵抗的扰动。function tau_env env_load(eta, nu, p) %ENV_LOAD 风 流产生的船体系合力二阶波浪漂移单独用成形滤波器给 % p 由 dp_init 生成 psi eta(3); R [cos(psi), -sin(psi), 0; sin(psi), cos(psi), 0; 0, 0, 1]; % ---------- 1. 风先算相对风速矢量 ---------- Vw_n [p.wind.speed*cos(p.wind.dir); p.wind.speed*sin(p.wind.dir)]; % 北东系 Vs_n R(1:2,1:2) * nu(1:2); % 船速北东系 Vrw_n Vw_n - Vs_n; % 相对风北东系 Vrw_b R(1:2,1:2) * Vrw_n; % 转到船体系 U_rw norm(Vrw_b); % 相对风速 [m/s] gamma atan2(Vrw_b(2), Vrw_b(1)); % 相对风向角 [rad] % 气动力系数按相对风向角查表0~2pi 周期 Cx interp1(p.wind.tab_ang, p.wind.tab_Cx, mod(gamma, 2*pi), linear, extrap); Cy interp1(p.wind.tab_ang, p.wind.tab_Cy, mod(gamma, 2*pi), linear, extrap); Cn interp1(p.wind.tab_ang, p.wind.tab_Cn, mod(gamma, 2*pi), linear, extrap); q 0.5 * p.rho_air * U_rw^2; % 动压 [Pa] tau_wind q * [Cx * p.vessel.AF; Cy * p.vessel.AL; Cn * p.vessel.AL * p.vessel.L]; % ---------- 2. 流不产生“外力”改的是相对速度 ---------- Vc_n [p.curr.speed*cos(p.curr.dir); p.curr.speed*sin(p.curr.dir)]; Vc_b R(1:2,1:2) * Vc_n; nu_c [Vc_b(1); Vc_b(2); 0]; % 流速在船体系的分量 tau_curr -p.D * (-nu_c); % 等价于阻尼项里用 (nu - nu_c) tau_env tau_wind tau_curr; end这段里有几个参数必须在dp_init里给全AF是水线以上正投影面积AL是侧投影面积L是垂线间长。气动力系数表tab_Cx/Cy/Cn一般按 0 到 360° 每 10° 一个点注意它是以相对风向角为自变量的周期函数查表时一定要做mod(gamma, 2*pi)否则负角度会走extrap外插出来的力可能大得离谱这条我踩过。二阶波浪漂移我没放进这个函数因为它需要状态。常见做法是用一个二阶低通滤波器驱动白噪声滤波器的中心频率取波浪谱的峰值频率阻尼比取 0.1 左右输出近似为慢变漂移力。把它单独做成一个带状态的模块好处是你可以直接调节海况等级而不用改风力模型。3. 让闭环第一次转起来初始化脚本、增益和推力分配的三层参数模型能单独跑之后接下来是把它接成闭环。这一章是真正决定仿真能不能出结果的环节绝大多数“跑不通”都卡在这三层中的某一层。3.1 初始化脚本的书写顺序把 base workspace 的隐式依赖去掉很多包的习惯是写一个脚本跑完之后变量全在 base workspace 里Simulink 模型直接引用变量名。这种写法对外人极不友好你永远不知道某个变量是脚本产生的、还是上次调试残留的。我自己改这类包第一件事是把脚本改成函数返回一个 struct模型里只认这一个变量名。function p dp_init() %DP_INIT 全参数唯一入口。Simulink 模型里只使用 p 这一个变量 p struct(); % ---- 仿真时间 ---- p.sim.Ts_phys 0.01; % 被控对象积分步长 [s] p.sim.Ts_ctrl 0.10; % 控制器采样周期 [s] p.sim.Tend 1800; % 单次仿真时长 [s] % ---- 船型主尺度 ---- p.vessel.L 76.2; % 垂线间长 [m] p.vessel.B 18.8; % 型宽 [m] p.vessel.T 6.25; % 吃水 [m] p.vessel.Cb 0.65; % 方形系数 p.vessel.AF 320; % 水线以上正投影面积 [m^2] p.vessel.AL 780; % 水线以上侧投影面积 [m^2] % ---- 船体矩阵 ---- p.M [5.3122e6, 0, 0; 0, 8.2831e6, 0; 0, 0, 3.7454e9]; p.D [5.0242e4, 0, 0; 0, 2.7229e5, -4.3933e6; 0, -4.3933e6, 4.1894e8]; % ---- 推进器布置[x位置, y位置, 最大推力, 类型] ---- % 类型 1 全回转本算例固定指向0 隧道式侧推 p.thr [-35, -6, 1.0e6, 1; -35, 6, 1.0e6, 1; 60, 0, 4.5e5, 0]; % ---- 环境缺省值 ---- p.rho_air 1.225; p.wind.speed 8; % [m/s] p.wind.dir deg2rad(45);% 来向北东系 p.curr.speed 0.6; % [m/s] p.curr.dir deg2rad(200); end顺序上有个讲究先给仿真时间再给船型再给矩阵最后给环境和控制器。因为船型主尺度会喂给风载荷计算而矩阵是独立量两者顺序错了不影响结果但会影响你排查时的阅读顺序。所有带单位的量都在注释里标出单位这条看着傻等你隔两个月回来改参数时就知道值不值。统一用 struct 的另一个好处是批量仿真。后面要扫几十个风浪方向只需要复制 struct 改一个字段不用在脚本里反复assignin(base, ...)。3.2 从 PID 到 LQR增益怎么定为什么先跑 PIDDP 控制器最朴素的形态是三个独立 PID纵向位置误差给纵向力横向位置误差给横向力艏向误差给力矩。坑在于误差是在北东系算的力却要在船体系输出中间必须过一个R(ψ)这个旋转一旦漏掉船在艏向 90° 附近会表现出“横向通道失灵”。先用 PID 跑通再考虑 LQR理由很实际PID 你能从曲线上一眼看出是比例大了还是积分大了LQR 出问题你面对的是六个状态权重得反复试。function K lqr_gain(p) %LQR_GAIN 在零速、零艏向工作点线性化得到 3x6 状态反馈增益 % 状态 x [eta; nu]输入 u tau A [zeros(3), eye(3); zeros(3), -p.M \ p.D]; B [zeros(3); p.M \ eye(3)]; % 关键Q/R 用物理量程归一化不要全填 1 sigma_pos 5; % 允许位置误差 5 m sigma_psi 3*pi/180; % 允许艏向误差 3 deg sigma_F 3e5; % 允许力误差 300 kN sigma_N 3e7; % 允许力矩误差 3e7 N*m Q diag([1/sigma_pos^2, 1/sigma_pos^2, 1/sigma_psi^2, ... 1e-6, 1e-6, 1e-6]); % 速度项只需很小权重即可稳定 R diag([1/sigma_F^2, 1/sigma_F^2, 1/sigma_N^2]); K lqr(A, B, Q, R); endQ和R的填法是这段代码的核心。用默认的eye(6)和eye(3)去算位置项量级 1 m和力矩项量级 1e7 N·m在同一套权重下完全失衡出来的增益会让推进器一直打满。按“我能容忍多大误差”折算成1/σ²是唯一靠得住的做法。速度项的权重给个 1e-6 就够它只负责抑制超调给大了会让系统变得迟钝。另外A里我把psi直接按 0 处理了。DP 工作点附近艏向变化不大线性化误差可以接受如果你要验证大角度回转工况这个增益就得跟着艏向调度。3.3 QP 推力分配等式约束、不等式约束和权重矩阵控制器给的是合力与合力矩推进器给的是各自的推力中间隔着一个几何矩阵B。三个推进器只有三个未知量时可以直接解析求解但一旦加上推力上限、变化率限制就必须用二次规划。function [u, exitflag] alloc_qp(tau_d, p) %ALLOC_QP 三推进器推力分配固定指向简化版 % 决策变量 u [T1; T2; T3] 单位 N % T1/T2 为艉部两台全回转本算例指向纵向T3 为艏部侧推指向横向 B [ 1, 1, 0; % 纵向合力 0, 0, 1; % 横向合力 6, -6, 60]; % 艏摇力矩-y_i*Fx_i 或 x_i*Fy_i Tmax p.thr(:,3); % 能耗项按最大推力平方归一化让三台推进器可比 W diag(0.01 ./ Tmax.^2); % 指令偏差项同样按量程归一化300 kN / 3e7 N*m 的容忍度 Q diag([1/3e5^2, 1/3e5^2, 1/3e7^2]); H 2 * (B * Q * B W); f -2 * B * Q * tau_d; % 不等式约束推力变化率限幅避免推力锯齿 du_max 0.02 * Tmax; % 每周期最大变化 2% 额定推力 Aineq [ eye(3); -eye(3) ]; bineq [ du_max; du_max ]; opts optimoptions(quadprog, Display, off, Algorithm, interior-point-convex); [u, ~, exitflag] quadprog(H, f, Aineq, bineq, [], [], -Tmax, Tmax, [], opts); endB矩阵第三行的数值是几何量艉部推进器在y ±6 m产生纵向力时对艏摇的贡献是-y·Fx所以一个是 6、一个是 -6艏部侧推在x 60 m产生横向力时贡献x·Fy 60。这个 6 和 60 差了一个数量级正是权重必须归一化的原因如果Q全填 1QP 会把绝大部分优化预算花在满足 2.7e7 N·m 量级的力矩上横向力被严重牺牲表现出来就是定点保持时横向误差收敛不了。du_max那组约束是让推力曲线好看的关键。仿真里推力高频抖动的根源往往不是控制器而是分配器每个周期都在两个近似最优解之间跳加上变化率约束就平滑了。需要说明的边界这个版本把全回转推进器的指向固定了。真实的全回转推进器角度也是决策变量B会随角度变化得用序列二次规划或者“先解推力再解角度”的两级结构。做能力图那种大范围扫描时固定指向的结果会偏乐观尤其是横向来风工况通常要打个折扣看。3.4 第一次仿真结束必须看的四条曲线仿真跑完别只看船的轨迹图。第一条是位置误差|η(1:2) - η_d(1:2)|的时间历程看它是不是单调收敛、有没有稳态偏差。第二条是艏向误差量级应该和你在Q里设的sigma_psi差不多。第三条是三台推进器的推力使用率|T_i|/Tmax_i如果长期贴在 1.0说明这个风浪条件下根本没有可行解模型的“定点保持”是假的。第四条是指令合力与分配后实际合力的差值B·u - τ_d正常情况下应该接近零如果常年有几十 kN 的残差回去检查Q的量纲。% 后处理四条曲线一次出 figure(Position, [100 100 900 700]); subplot(2,2,1); plot(out.t, out.pos_err); grid on; ylabel(位置误差 [m]); subplot(2,2,2); plot(out.t, rad2deg(out.psi_err)); grid on; ylabel(艏向误差 [deg]); subplot(2,2,3); plot(out.t, abs(out.u)./Tmax); grid on; ylabel(推力使用率); subplot(2,2,4); plot(out.t, out.tau_res); grid on; ylabel(分配残差 [N]);这四张图一起看能区分三种情况误差不收敛且推力没饱和是控制器问题误差不收敛且推力顶格是能力不足误差收敛但推力抖是分配器或采样率问题。4. 三个典型工况复现定点保持、低速循迹、DP 能力图参数对了之后模拟器的价值才开始体现——它能批量跑工况。这一章给三个最常用的复现脚本每个都对应一个实际验收指标。4.1 定点保持斜坡指令和三条验收标准最基础也最常被做错的是定点保持测试。直接给一个阶跃的位置指令推进器瞬间打满积分项迅速累积船冲过头之后要花很长时间才拉回来得到的“超调量”完全是测试方法造成的不代表控制器性能。正确做法是斜坡指令目标点分 60 到 120 s 缓慢走过去。p dp_init(); p.ctrl.mode DP; p.ref.mode ramp; % 斜坡参考避免阶跃打满推进器 p.ref.t0 60; % 60 s 开始机动 p.ref.T 120; % 斜坡时长 p.ref.eta_d [50; 30; deg2rad(15)]; out sim(dp_main.slx); % 指标提取稳态误差取最后 10% 时间窗 idx out.t 0.9*p.sim.Tend; ess_pos rms(vecnorm(out.eta(idx,1:2) - p.ref.eta_d(1:2), 2, 2)); ess_psi rms(abs(out.eta(idx,3) - p.ref.eta_d(3)));验收上我一般卡三条稳态位置误差小于 2 m艏向误差小于 0.5°全过程推力使用率不超过 80%。第三条最容易被忽略但最有价值——长期在 90% 以上运行的定点保持遇到一阵突风就没有余量了实船上这是要报警的。注意斜坡时长p.ref.T不是越长越好。太长会把机动过程拉得很慢看起来“很稳”实际上掩盖了控制器响应不足的问题。按 0.3 到 0.5 m/s 的等效移动速度折算比较合适。4.2 低速循迹LOS 制导和模式切换的衔接动力定位不只有定点保持低速循迹沿一条测线以 0.5 m/s 左右的速度走也是常见工况比如水下设备作业、海底管线巡检。这时参考点不是固定的而是一个沿航线移动的目标点。制导律用视线法LOS就够了算当前位置到目标航路点的方向加上一个前视距离做平滑。function psi_d los_heading(eta, wp, Delta) %LOS_HEADING 视线制导输出期望艏向 % Delta 为前视距离典型取 1.5 到 3 倍船长 dx wp(1) - eta(1); dy wp(2) - eta(2); if hypot(dx, dy) Delta psi_d eta(3); % 进入到达圆内保持艏向避免原地打转 else psi_d atan2(dy, dx); end end真正麻烦的是从定点模式切到循迹模式的那一瞬间。两个模式的参考点如果不连续位置环的积分器会突然看到一个很大的误差输出一个巨大的力矩船艏直接甩一下。处理办法是切换时做两件事把当前实际位置直接赋给循迹模式的初始参考点实现参考连续同时把控制器的积分项清零或者做无扰切换。这两步在 Simulink 里可以用一个Reset端口配合Enabled Subsystem实现比在脚本里改初值可靠。4.3 能力图扫描一次跑几十个工况出表比看图快DP 能力图回答的是“这个风浪条件下我还能不能保持位置”。做法很直接固定一组海况扫风向和风速看每个组合下稳态位置误差是否超限。wind_dirs 0:15:345; % 24 个风向 wind_speeds 6:2:24; % 10 档风速 cap nan(numel(wind_dirs), numel(wind_speeds)); for i 1:numel(wind_dirs) for j 1:numel(wind_speeds) p dp_init(); p.wind.dir deg2rad(wind_dirs(i)); p.wind.speed wind_speeds(j); p.sim.Tend 1200; % 稳态工况不用跑太久 out sim(dp_main.slx); idx out.t 0.8*p.sim.Tend; % 取后 20% 平均值作为稳态 err mean(vecnorm(out.eta(idx,1:2) - p.ref.eta_d(1:2), 2, 2)); cap(i,j) err; end end % 判定误差小于 5 m 视为保持成功 pass cap 5;240 个工况串行跑会很慢。我的做法是用parfor把外层风向循环并行内层风速串行改动量最小。也可以用SimulationInput数组配合parsim更适合把参数打包成结构体的情况。结果我倾向于先出一张表而不是直接画极坐标图。表里pass为 false 的格子一眼就能看出集中在哪个风向区间比看图快本文还有配套的精品资源点击获取