简介：一套基于MATLAB的船舶航向PID控制程序，面向船舶自动化与运动控制领域，适合学习PID整定方法、船舶六自由度建模的学生和工程师。压缩包内共包含21个文件，全部为.m脚本，大小仅46KB，按功能涵盖船舶运动模型、PID控制器实现、参数调节以及避碰仿真等模块，结构清晰便于按需调用。目前已有3038人学习下载，常用于课程设计、科研验证或工程入门。通过阅读和运行源码，可掌握船舶运动方程的数值求解过程，理解比例、积分、微分三个参数对航向保持稳定性和响应速度的影响，同时体验不使用Simulink、纯M代码搭建仿真环境的思路。代码注释清晰，可灵活调整舵角限幅、控制周期等关键参数，观察航向响应曲线并总结整定规律，是理论联系实际的实用参考。

1. 从一条阶跃响应曲线说起：MATLAB 里搭船舶航向 PID 仿真到底在做什么

拿到“船运动程序 matlab_船舶；船舶航向_PID 控制”这类题目，通常不是要你在实船上做实验，而是要你在 MATLAB 里搭一套带船体模型的闭环仿真环境，先把航向 PID 控制器的参数整定明白。这件事的实用价值很清楚：实船调参成本高、风险大，一次操舵试验动辄几十分钟，而在仿真里改一个 Kp 值只要重跑一遍脚本。把船体简化成低阶模型，在纯软件环境里完成控制器选型、参数整定和极端工况测试，是船舶运动控制方向最常用的第一站。

这套方案适合两类人：一是做船舶自动驾驶、智能航行课题的学生，需要航向控制环的理论验证；二是负责舵控系统预研的工程师，需要在详细模型之前拿到一组合理的 PID 初值。本文按“船体模型 → 控制器实现 → 闭环仿真 → 参数整定与避坑 → 进阶验证”的顺序展开，每一步都给可复现的代码和参数，你照着搭完，就能看到一条从初始航向收敛到目标航向的完整曲线。

2. 先立住船体：一阶 Nomoto 模型和两个决定性的参数 K、T

做航向控制仿真，最忌讳一上来就搞六自由度全耦合模型。船舶运动确实复杂，横荡、纵荡、艏摇、横摇都耦合在一起，但航向 PID 关心的核心自由度只有一个——艏摇（yaw）。在这个前提下，业界最常用的做法是用一阶 Nomoto 模型把船体压缩成一个低阶传递函数。模型越简单，后面的控制律设计和参数整定越容易落地。

2.1 为什么选一阶 Nomoto 模型

从船舶操纵性理论推导，艏摇方向的线性化运动方程可以整理成如下传递函数：

[ \frac{\psi(s)}{\delta(s)} = \frac{K}{s(1 + T s)} ]

其中 (\psi) 是航向角，(\delta) 是舵角，K 是回转性增益（单位 1/s），T 是应舵时间常数（单位 s）。这个模型的物理含义是：给舵一个阶跃输入，航向角先以接近匀速的速率转动，受 T 影响逐渐过渡到稳定回转率 K·δ。大船 T 值大，舵响应慢；小船 T 值小，舵响应快。

选择一阶 Nomoto 的理由有两层。第一层是控制设计角度：航向 PID 只需要抓住被控对象的主导动态，K 和 T 恰好刻画了从舵角到航向响应最核心的两个特征——稳态增益和响应快慢，偏航、横荡这些因素对航向环的影响在小舵角时是二阶小量。第二层是工程习惯：船舶操纵性标准里 Nomoto 模型是公认的对比基准，实船辨识出来的 K、T 参数可以直接用。我在做预研时也常用二阶 Nomoto，但那是为了刻画艏摇加速度的初始滞后，PID 整定阶段一阶模型完全够用。

2.2 K 和 T 怎么定：实船量级与航速换算

K、T 不是随便拍的，它们和船长、航速直接相关。常见的做法是先拿到无量纲的参数 (K') 和 (T')（来自船模试验或经验公式），再按实船尺度换算：

[ K = K' \cdot \frac{V}{L}, \quad T = T' \cdot \frac{L}{V} ]

其中 V 是航速，L 是船长。实际项目中如果没有试验数据，可以参考下面的典型量级范围：

这套仿真里我习惯取 K=0.16、T=100，代表一艘中等吨位的商船。你后面自己调参时注意 K、T 会随航速变化——低速时 T 变大、K 变小，船变得“迟钝”。因此仿真阶段最好跑几个不同航速的参数包络，而不是只盯着一个工作点。

2.3 用 m 函数写出被控对象的微分方程

有了 K、T 之后，把传递函数写成状态空间形式，便于后续加非线性限幅和干扰。状态取 ([\psi, r]^T)，其中 r 是艏摇角速度，一阶 Nomoto 的微分方程为：

% ship_nomoto_model.m % 一阶 Nomoto 模型状态方程 % 输入: state = [psi; r]，舵角 delta(弧度)，参数 K,T % 输出: dstate = [dpsi; dr] function dstate = ship_nomoto_model(state, delta, K, T) psi = state(1); % 航向角，单位 rad r = state(2); % 艏摇角速度，单位 rad/s dpsi = r; dr = (K * delta - r) / T; dstate = [dpsi; dr]; end

这里 dr 的表达式就是传递函数倒回时域的结果：(T \dot{r} + r = K \delta)。写成函数而不是直接写在主脚本里，是因为后面要做 RK4 积分、加干扰项，函数形式更容易复用。

数值积分我不用 Euler，虽然 Euler 简单但稳定性差，dt 稍大就会发散。这里用标准 RK4：

% rk4_step.m % 用 RK4 推进一个仿真步长 function state_next = rk4_step(state, delta, dt, K, T) k1 = ship_nomoto_model(state, delta, K, T); k2 = ship_nomoto_model(state + 0.5*dt*k1, delta, K, T); k3 = ship_nomoto_model(state + 0.5*dt*k2, delta, K, T); k4 = ship_nomoto_model(state + dt*k3, delta, K, T); state_next = state + (dt/6) * (k1 + 2*k2 + 2*k3 + k4); end

RK4 的截断误差比 Euler 高了三个阶次，在 T=100 秒的慢系统上，dt 取 0.1 秒绰绰有余。如果后面你把舵机响应也建模进去，舵机时间常数往往只有几秒，那时 dt 要缩小到 0.01~0.05 秒，振荡频率才能被积分步长捕捉到。

3. PID 控制器：为什么用增量式，以及如何写进 MATLAB

船体模型准备好之后，接下来的核心是控制器。航向 PID 的实现方式有两种：位置式和增量式。虽然教科书上位置式写起来更直观，但实际做船舶舵控我几乎不用位置式，原因在下面展开。

3.1 位置式与增量式 PID 的关键取舍

位置式 PID 的表达式是：

[ \delta(t) = K_p e(t) + K_i \int_0^t e(\tau)d\tau + K_d \frac{de(t)}{dt} ]

输出直接就是舵角绝对值。它的第一个问题是积分饱和：大角度转向时误差大，积分项快速累积，等航向追上了，积分还在往舵角上叠加，导致明显超调，配合上舵机的 ±35° 限幅，很难收住。第二个问题是切换工况时会产生舵角跳变：如果从手动操舵切到自动舵，位置式 PID 的输出和历史舵角不一致，舵机突然跳动，这在实船上是不可接受的安全隐患。

增量式 PID 的输出是舵角增量 (\Delta \delta)，通过累加得到当前舵角：

% pid_incremental.m % 增量式 PID，输出舵角（带限幅和转舵速率限制） function [delta, e_prev, e_prev2] = pid_incremental(e, e_prev, e_prev2, ... Kp, Ki, Kd, dt, delta_prev) % e: 当前航向误差 = psi_set - psi，单位 rad % e_prev / e_prev2: 上一时刻和上上时刻的误差 % Kp / Ki / Kd: PID 增益 % dt: 控制周期，单位 s % delta_prev: 上一时刻舵角输出，单位 rad % 增量计算：比例项用误差一阶差分，微分项用二阶差分 ddelta = Kp * (e - e_prev) + Ki * e * dt + Kd * (e - 2*e_prev + e_prev2) / dt; delta = delta_prev + ddelta; % 转舵速率限制：假设舵机最大转舵速度 3°/s rate_limit = 3 * pi / 180 * dt; if delta - delta_prev > rate_limit delta = delta_prev + rate_limit; elseif delta_prev - delta > rate_limit delta = delta_prev - rate_limit; end % 舵角限幅：典型舵机最大舵角 ±35° delta_max = 35 * pi / 180; if delta > delta_max delta = delta_max; elseif delta < -delta_max delta = -delta_max; end % 更新误差历史 e_prev2 = e_prev; e_prev = e; end

注意三个增益的量纲完全不同：Kp 无量纲，Ki 的量纲是 1/s，Kd 的量纲是 s。很多人直接套位置式 PID 的数值到增量式里，结果整个系统发散，就是因为忘了 Ki 项在增量式里乘了 dt，数值被放大了几十倍。增量式中积分项本质上是 (K_i \cdot e \cdot dt)，是对误差在控制周期内的累积量做比例放大。

3.2 为什么把速率限制做进控制器

实船舵机不是理想执行器，舵叶有最大转舵速度，通常商船在 2~5°/s 之间。仿真里如果不加这个限制，PID 输出一个阶跃舵角，船模立刻响应，看起来性能很好，但实船根本做不到这么快的舵速。把速率限制加进去后，控制环里多了一个“执行器滞后”，PID 参数必须重新整定，这也是很多人仿真调得好、上船就抖的原因之一。

另外我在控制器里加了死区逻辑的话，会在主循环中做判断，而不是放在函数里：当航向误差小于 0.1° 时，不更新误差历史也不输出舵角增量，避免舵机在小误差区间来回抖。这个阈值看实际需要，一般商船的自动舵死区设在 0.1°~0.3°。

4. 联手跑闭环：主循环、航向保持与航向改变两个工况

模型和控制器都有了，现在把它们接到一个主循环里。这一步的目标是跑出一条完整的航向响应曲线，同时把数据记录下来画图。这个循环结构是后面所有工况测试和参数整定共用的底座。

4.1 主循环代码框架

% run_heading_pid.m % 船舶航向 PID 控制闭环仿真主循环 clear; clc; % 船体模型参数（中等吨位商船典型值） K = 0.16; % 回转性增益，1/s T = 100; % 应舵时间常数，s % PID 控制器参数（初始值，后面会整定） Kp = 0.05; Ki = 0.0002; Kd = 1.0; % 仿真与控制参数 dt = 0.1; % 仿真步长 = 控制周期，s simT = 1500; % 仿真总时长，s steps = round(simT / dt); % 状态初始化 psi = 5 * pi / 180; % 初始航向 5° r = 0; % 初始艏摇角速度为 0 psi_set = 5 * pi / 180; % 目标航向初始为 5° delta = 0; % 初始舵角 % 误差初始 e_prev = psi_set - psi; e_prev2 = e_prev; % 数据记录 time_log = (0:steps-1)' * dt; psi_log = zeros(steps, 1); delta_log = zeros(steps, 1); for k = 1:steps % 900s 时刻发出 25° 的航向改变指令 if time_log(k) >= 900 psi_set = 25 * pi / 180; end e = psi_set - psi; % 死区判断：误差小于 0.1° 时不操舵 if abs(e) < 0.1 * pi / 180 e_prev = e; e_prev2 = e_prev; continue; end % 调用增量式 PID [delta, e_prev, e_prev2] = pid_incremental(e, e_prev, e_prev2, ... Kp, Ki, Kd, dt, delta); % RK4 推进船体模型 state = [psi; r]; state_next = rk4_step(state, delta, dt, K, T); psi = state_next(1); r = state_next(2); % 记录 psi_log(k) = psi * 180/pi; delta_log(k) = delta * 180/pi; end % 绘图 figure; subplot(2,1,1); plot(time_log, psi_log, 'b-', 'LineWidth', 1.2); hold on; plot(time_log, ones(steps,1) * 25, 'r--'); % 目标航向 ylabel('航向 (deg)'); title('航向响应'); grid on; subplot(2,1,2); plot(time_log, delta_log, 'g-', 'LineWidth', 1); xlabel('时间 (s)'); ylabel('舵角 (deg)'); title('舵角输出'); grid on;

这个循环里有几个关键设计。第一，控制周期和仿真步长保持一致，都用 dt，避免出现“控制器 0.1 秒跑一次、模型 0.05 秒积分一次”的双速率错位，那是很多诡异振荡的来源。第二，900 秒处修改目标航向，同一个循环内既能看到 5° 小偏差的航向保持，又能看到 5° 到 25° 的大幅度航向改变。第三，死区判断直接 continue 跳过控制输出，但误差历史还是要更新，否则恢复控制时微分项会产生一个很大的跳变脉冲。

4.2 两个典型工况的验收指标

跑完仿真后，重点看两个时段的响应曲线。第一段是前面 900 秒，初始航向和目标航向一致，考察的是抗扰动能力，仿真里虽然没有外加干扰，但数值积分本身会给系统微小的激励，曲线如果稳定不动，说明闭环至少没被数值误差激发。第二段是 900 秒后的航向改变，这里要看四个指标：

商船 T=100 秒的系统，调节时间在 400~600 秒是正常的，不要拿小艇的响应速度去套大船。如果你第一次跑出来的曲线振荡很厉害，先检查 Kp，把 Kp 减半看看；如果纯发散，检查 dt 是否太大，或者 Kp 是否已经超过稳定边界。

4.3 曲线怎么“读”：超调、振荡和舵角抖动的判别

看仿真曲线不是只看“有没有收敛”，要看形态。单调缓慢爬升然后平缓收敛，说明 Kp 偏小，系统太“肉”；快速爬升但超调大、来回振荡 3 个周期以上，说明 Kp 偏大，阻尼不足；稳态阶段舵角曲线有肉眼可见的高频抖动，说明微分项放大了数值噪声。我的经验是先看舵角曲线再看航向曲线，舵角抖动往往先于航向异常出现，能提早发现问题。

5. 避坑与排查：仿真发散、微分噪声、积分饱和的 5 个现场

这套仿真框架跑起来不难，但实际调试时翻车点非常集中。以下 5 个问题是我在这个方向上处理过最多的情况，按“现象 → 原因 → 解决”逐一展开。

5.1 一加控制量仿真就发散，航向曲线直接拉飞

现象：主循环一跑，航向曲线指数式向上或向下猛冲，几秒之内就突破 ±180°，舵角顶在 35° 限幅上不回来。

原因：绝大多数情况是两个：一是 Kp 太大，闭环极点越过虚轴，系统本身不稳定；二是仿真步长 dt 过大，RK4 的稳定性也有限，尤其当舵机速率限制或限幅引入非线性后，大 dt 下数值积分误差被控制环路反复放大。

解决：先把 Ki 和 Kd 置零，Kp 从 0.02 起步，跑通稳定后再逐步加。另一步是确认 dt，对于 T=100 秒的一阶 Nomoto，dt 取 0.1 秒已经是上限，如果加了舵机二阶模型，dt 必须缩到 0.01~0.02 秒。发散射时优先检查 Kp，再看 dt，顺序不要反。

5.2 微分项一打开，舵角噪声就爆炸

现象：航向曲线看着还行，但舵角曲线出现密集的高频小幅摆动。给反馈信号加一点高斯白噪声后更明显，甚至整个舵角输出都在 ±35° 之间来回跳。

原因：微分项的本质是误差变化率放大，而测量噪声的高频成分在差分运算下被成倍放大。增量式 PID 中微分项用了 (e) 的差分，噪声被放大了 (1/dt) 倍，dt 越小噪声越猛。

解决：一是让航向反馈经过低通滤波，再做差分；二是在控制器里直接去掉微分项，用 PI 控制。对船舶航向这种大惯性对象，PD 带来的阻尼提升有限，很多实际自动舵就是纯 PI 控制加滤波。如果确实需要微分，我常用的办法是“微分先行”：只对测量航向做微分，不对设定值做微分，防止指令阶跃时微分项打出尖峰。

5.3 K、T 估计不准导致整定好的参数“水土不服”

现象：仿真里调得不错的 PID，换一组 K、T 参数后，超调量从 5% 变成 30%，甚至直接发散。

原因：K、T 不是固定值。满载和压载，高速和低速，K 可能差一倍，T 可能差三倍。整定好的 PID 只对特定参数有效，这是控制系统设计的基本边界，也是新手最容易忽略的点。

解决：把 K、T 按船舶运营范围做成参数包络，仿真时分别跑高速满载、低速满载、高速压载等组合。每个组合下 PID 都能稳定，才说明参数有鲁棒性。我在实际项目中会把 K、T 做成一个结构体数组循环仿真，避免手工改一次跑一次。

5.4 积分饱和：大角度转向后超调收不住

现象：5° 到 25° 的阶跃响应里，舵角在转向初期长时间顶在 35° 限幅，航向到达 25° 后继续冲到 28° 以上，再花很长时间拉回来。

原因：增量式 PID 虽然没有积分累积项，但舵角累加量在限幅处饱和，等效于积分饱和。误差大时控制器持续输出满舵，误差减小时输出还要“消化”之前累积的效果。

解决：两处配合使用。第一处是在控制律里做积分分离，误差超过 2° 时切断积分项；第二处是判断舵角饱和后冻结积分输出，代码里可以这样处理：

% 积分分离：误差在死区外且小于阈值时，积分项生效 if abs(e) < 2.0 * pi / 180 ddelta = Kp * (e - e_prev) + Ki * e * dt + Kd * (e - 2*e_prev + e_prev2) / dt; else ddelta = Kp * (e - e_prev) + Kd * (e - 2*e_prev + e_prev2) / dt; end

这个改动在大角度转向工况下效果立竿见影，超调能压低一半以上。

5.5 周期性海浪扰动下，航向稳态误差一直存在

现象：模型里加了正弦形式的海浪扰动力矩后，航向曲线围绕目标值做周期性摆动，平均偏差始终在 2° 以上，加 Ki 也压不下去。

原因：PID 对频率接近控制带宽的周期性扰动抑制能力有限，Ki 只能消除常值误差，对特定频率的正弦扰动增益不够。航向控制本身也不需要完全抵消海浪，实船自动舵允许 1°~3° 的航向摆动，因为频繁操舵反而耗油、磨损舵机。

解决：先确认扰动频率，再决定策略。若扰动频率远高于控制系统带宽，靠船舶本身的惯性就能滤掉；若在带宽范围内，考虑加陷波滤波器，或者对波浪频率做前馈补偿。大多数工程场景下，接受小幅周期性航向波动就是最经济的选择，不要指望 PID 把所有干扰都消干净。

6. 给 PID 参数做体检：频域裕度与随机工况验证

到这一步，你已经能跑出稳定的航向响应曲线。但稳定不等于可靠，最后一环是给这套 PID 参数做“体检”，确认它在参数变化和工况变化时依然站得住。我一般会做两个验证，都很轻量，但能提前暴露很多实船问题。

6.1 用 margin 函数检查稳定裕度

把一阶 Nomoto 模型和 PID 构成连续开环传递函数，用 MATLAB 的 margin 函数直接算增益裕度和相位裕度：

% margin_check.m % 用连续 PID 作为近似，检查闭环稳定裕度 K = 0.16; T = 100; Kp = 0.05; Ki = 0.0002; Kd = 1.0; C = pid(Kp, Ki, Kd); G = tf(K, [T 1 0]); % K / (s*(1 + T*s)) [Gm, Pm] = margin(C * G); fprintf('增益裕度: %.1f dB\n', 20*log10(Gm)); fprintf('相位裕度: %.1f deg\n', Pm);

相位裕度大于 45°、增益裕度大于 6 dB，是一般工程底线。如果算出来相位裕度只有 20°，那么前面时域仿真里看起来不错的曲线，实际上一遇到模型参数变化就会剧烈振荡。频域体检能帮你把关卡的往前提，不用等到实船出问题。

6.2 跑 20 组随机工况看鲁棒性

每个整定好的参数方案，我会跑一组随机起始航向和目标航向的批量仿真，统计超调大于 20% 的工况比例。这个习惯帮我挡掉过好几次“看起来调好了、换个角度就崩”的情况。代码不需要完整重写，只改主循环里的初始化条件和外层包一层循环：

% robustness_check.m % 随机工况鲁棒性检查（示意） rng(42); targets = -30 + 60 * rand(20,1); % 目标航向 -30° ~ 30° overshoot_count = 0; for i = 1:20 psi_set = targets(i) * pi / 180; psi = 0; % 初始航向一律从 0° 开始 % 运行和主脚本相同的仿真循环 % ...（略） % 如果最大超调超过 20%，则记录失败 if max_overshoot > 20 overshoot_count = overshoot_count + 1; end end disp(['超调>20%的工况占比: ' num2str(overshoot_count/20*100) '%']);

占比超过 10% 就说明参数整定过激，我会把 Kp 和 Kd 同时下调 20% 再重新跑。这里我还有一个习惯：把这套“批量验证”的代码长期留在脚本包里，后面往实船参数过渡时，每换一组 K、T 就跑一遍同样验证。这套流程帮我避免了不止一次把仿真参数直接带上真船的冲动。

无论是做学生课题还是工程预研，记住一点：PID 参数是模型和控制器共同作用的结果，模型不准，再好的整定方法都是空中楼阁。先花时间把 K、T 边界摸清楚，再用频域和随机工况验证兜底，比执着于某个参数的“最优值”可靠得多。希望帮到你。

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