企业官网建设流程全解析

简介：面向船舶工程与控制工程学习者的动力定位仿真项目，围绕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：增益怎么定，为什么先跑 PID

DP 控制器最朴素的形态是三个独立 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); end

Q和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); end

B矩阵第三行的数值是几何量：艉部推进器在y = ±6 m，产生纵向力时对艏摇的贡献是-y·Fx，所以一个是 +6、一个是 -6；艏部侧推在x = 60 m，产生横向力时贡献x·Fy = 60。这个 6 和 60 差了一个数量级，正是权重必须归一化的原因：如果Q全填 1，QP 会把绝大部分优化预算花在满足 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 的格子一眼就能看出集中在哪个风向区间，比看图快

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