Appearance
PID 控制：原理与整定
概念
上一章算出了"离不稳定还差 17.96 度"——知道差多少了，可怎么补？
工程师面对的选择其实只有三类：改被控对象（换电机、换阀门，成本最高）、改结构（加补偿网络，需要重新设计电路）、调一个现成的控制器。第三类的答案，近百年来几乎只有一个：PID。
一句话说清它是什么：PID 就是把误差的"现在、过去、未来"三种信息，各自乘以一个系数再加起来当控制量。
| 项 | 看的时段 | 一句话作用 | 对应上一章的哪个操作 |
|---|---|---|---|
| P（比例） | 现在 | 有多大误差就给多大力度 | 只改增益 K，不改轨迹形状 |
| I（积分） | 过去 | 把历史上所有没消掉的误差攒起来一起还 | 加一个开环极点（1/s） |
| D（微分） | 未来 | 按误差变化的势头提前刹车 | 加一个开环零点（近似） |
为什么三项都是"必要的"：只有 P 时，输出必须靠"存在一个非零误差"来维持控制量——这叫有差调节；I 项只要误差没归零就持续累加，最终一定把静差压到零；而 I 引进了 90 度的相位滞后、把稳定裕度吃掉，D 用超前相位把它补回来。三项之间的拉扯关系，就是整定这件事的全部内容。
PID 为什么能统治工业界近百年：它不依赖被控对象的精确模型。上面那条式子只用到"误差"这一个可测量，不需要知道被控对象是几阶、参数多少。对于 90% 以上的过程控制回路，PID 的效果已经"足够好"，而性价比无人能敌。
这一层回答了上一层什么问题：上一章说"相位裕度只有 17.96 度、偏低"，这一章给出具体怎么补——减 P、加 D；也给出了代价——加 I 会进一步吃裕度（本节会看到 PI 的超调反而比纯 P 大得多）。
原理
一、三项各管什么、各要什么代价
| 项 | 主要贡献 | 代价 | 频域/根轨迹视角 |
|---|---|---|---|
P（比例 Kp） | 加快响应、提高抗扰 | 增大超调、降低稳定性；单用 P 存在稳态误差（有差调节） | 只是沿原轨迹滑动，不改变形状；增益越大极点越靠右 |
I（积分 Ki） | 消除稳态误差（无差调节） | 引入 90 度相位滞后、降低稳定性；响应变慢、超调变大 | 增加一个位于原点的开环极点，把轨迹整体往右拉 |
D（微分 Kd） | 提供超前相位、增大阻尼、抑制超调 | 放大高频噪声（噪声的微分被成倍放大）；对纯延迟无力 | 增加一个开环零点，把轨迹往左拉 |
三项的"分工"可以这样记：
P管"快"（力度）；I管"准"（精度）；D管"稳"（阻尼）。
但"稳"和"快"是矛盾的：I 用稳定换精度，D 用噪声敏感度和"对延迟无能为力"换阻尼。所以 PID 整定的第一原则不是"把三项都调大"，而是在这三组交换里选一个工程上可接受的折中点。
为什么 D 对纯延迟无力：延迟环节 e^{-τs} 的导数在频域上是 jw·e^{-τs}——它的相位是 -wτ + 90 度。在延迟占主导的高频段，D 提供的 +90 度远远抵不过延迟带来的 -wτ，所以对"大延迟对象"，D 的收益有限，通常干脆不加（工程上常直接用 PI）。
二、两种参数写法与两种数字实现
两种参数写法（考试与工程文档里都会出现）：
换算关系：
从 Ki/Kd 到 Ti/Td | 从 Ti/Td 到 Ki/Kd |
|---|---|
Ti = Kp / Ki | Ki = Kp / Ti |
Td = Kd / Kp | Kd = Kp · Td |
Ti 与 Td 的物理意义更容易理解：Ti 是积分时间（"攒够多长时间才把误差补完"），Td 是微分时间（"按多长时间里的变化率提前动作"）。所以工程整定表一律给 Kp、Ti、Td。
两种数字实现：
| 位置式 | 增量式 | |
|---|---|---|
| 算式 | u(k) = Kp·e(k) + Ki·dt·Σe + (Kd/dt)·[e(k) − e(k−1)] | Δu(k) = Kp·[e(k) − e(k−1)] + Ki·dt·e(k) + (Kd/dt)·[e(k) − 2e(k−1) + e(k−2)] |
| 输出 | 直接算出控制量的绝对值 | 只算控制量的增量 Δu，累加得到 u |
| 优点 | 直观、便于加限幅 | 无需累加、切换无冲击、积分饱和影响小 |
| 缺点 | 需要累加积分项，容易积分饱和 | 需要保存 e(k−1)、e(k−2) 两个历史值 |
两者数学上完全等价（本节示例会逐点核对到 10⁻¹⁴ 量级），差别只在数值实现上。工程上偏好增量式的一个实际理由：手动/自动切换时没有冲击——位置式切换瞬间如果积分项累加值与手动输出不一致，会给执行机构一个突变，增量式只输出"增量"，天然平滑。
三、Ziegler-Nichols 两种整定法
这是 PID 整定最经典的两张表，1942 年由 Ziegler 与 Nichols 提出，至今仍在用。
方法一：临界比例度法（闭环法）——先让系统"自己振荡起来"，再从振荡参数反推
步骤：
- 只加比例控制（
Ti -> ∞、Td = 0），从零开始缓慢增大Kp； - 一直加到系统出现等幅振荡，记下这时的增益
Ku（临界增益）与振荡周期Tu（临界周期）； - 查表：
| 控制器 | Kp | Ti | Td |
|---|---|---|---|
| P | 0.50 Ku | — | — |
| PI | 0.45 Ku | Tu / 1.2 | — |
| PID | 0.60 Ku | 0.5 Tu | 0.125 Tu |
关键洞察在第 2 步：Ku 就是上一章根轨迹里"轨迹与虚轴的交点所对应的增益"——所以这套整定法可以直接用前两章的工具算出来，不用真的去实验台上试。 本节示例就是这么做的。
方法二：阶跃响应法（开环法）——让对象自己走一条阶跃曲线，从曲线上读两个参数
给对象加一个阶跃，画出响应曲线，在拐点处作切线，得到：
K：稳态输出 / 阶跃幅值（对象增益）L：纯延迟时间（切线在时间轴上的截距）T：时间常数（切线从拐点到稳态的时间跨度）
| 控制器 | Kp | Ti | Td |
|---|---|---|---|
| P | T / (K·L) | — | — |
| PI | 0.9 T / (K·L) | 3.3 L | — |
| PID | 1.2 T / (K·L) | 2 L | 0.5 L |
这套方法被称为"开环 Z-N"或"反应曲线法"，只需要一次阶跃实验，比闭环法更快、更安全（不会让系统振荡）。它的适用条件是对象近似一阶惯性加纯延迟（FOPDT）——这是过程控制里最常见的模型形态。
两套整定法的共同特点：它们都给出"偏激进"的参数。Z-N 的设计目标叫四分之一衰减（相邻两个波峰之比约 4:1），换算成超调量往往在 25% 以上；如果希望超调更小，工程上常规做法是把 Kp 乘以 0.5 ~ 0.8 再使用。
四、参数变化的方向性影响
整定现场最常用的就是这张"方向表"——不看数字，只看往哪边拧：
| 调整 | 上升时间 | 超调量 | 调节时间 | 稳态误差 | 稳定性 |
|---|---|---|---|---|---|
Kp 增大 | 变短 | 变大 | 小扰动变小、大扰动变大 | 减小 | 变差 |
Ti 增大（即 Ki 减小） | 变长 | 变小 | 可能变长 | 消差变慢 | 变好 |
Td 增大 | 略变短 | 变小 | 变小 | 无影响 | 变好（但抗噪变差） |
这张表里唯一"全是好处"的是 Td——除了抗噪能力。这就是 D 环节"看起来是免费的午餐"的原因，也正好是它的陷阱：真实信号里必然有噪声，D 把噪声的导数放大后注入执行机构，会让执行机构高频抖动、加速磨损。所以工程上 D 通常要配一个低通滤波（"不完全微分"），或者干脆不用。
五、积分饱和与三种抗饱和方案
积分饱和是 PID 落地时最常出的事故，机理一句话：
只要误差不为零，积分项就一直在累加。当执行机构已经饱和（阀门全开、电压到顶）而误差仍然存在时，积分项会一路累加到天文数字。等误差终于反号，积分项需要很长时间才能"退回来"——这段时间里系统完全失控。
现象：阶跃响应出现一个巨大的、与控制器参数不匹配的超调，之后是很长的恢复过程。
| 方案 | 做法 | 特点 |
|---|---|---|
| 积分限幅（clamping） | 给积分项本身设上下限（如"积分贡献不超过稳态控制量"） | 最简单，效果立竿见影；上限要凭经验 |
| 反算法（back-calculation） | 把"饱和差额"反馈回积分器：I += dt·(u_sat − u)/Kp | 无需手调上限，饱和解除后自然回到正常 |
| 条件积分（conditional integration） | 饱和期间暂停积分 | 逻辑简单，但可能造成稳态误差残留 |
三种方案的核心思想是一样的：不要让积分器在"自己的输出根本用不上"的时候继续积累。 本节示例会用同一个被控对象、同一组 PI 参数，对比"无抗饱和 / 积分限幅 / 反算法"三者的差别。
六、整定的顺序：先 P、再 I、最后 D
现场最常用的手工整定顺序（比查表更常被使用）：
Ti = ∞、Td = 0，只调Kp：加到响应刚出现一点超调（约 10% ~ 20%）为止；- 加入
I：从大往小减小Ti（即增大Ki），直到静差在可接受时间内消失，同时超调还能忍； - 加入
D：从零缓慢增大Td，直到超调被压到位，同时执行机构不出现明显抖动。
这个顺序不能颠倒：P 定的是"力度"，I 和 D 都是在某个力度下"修形状"。力度没定就修形状，等于在流沙上盖房子。
示例
例 1：用前两章的工具直接算出整定参数
对象取 1/(s(s+1)(s+2))——正是根轨迹那一章的例子。04-root 已经算过：临界增益 Ku = 6，此时振荡角频率 wu = 1.414214 rad/s。于是：
Tu = 2*pi / wu = 2*pi / 1.414214 = 4.442883 s
P : Kp = 0.50 * Ku = 3.000000
PI : Kp = 0.45 * Ku = 2.700000，Ti = Tu / 1.2 = 3.702402 s
PID : Kp = 0.60 * Ku = 3.600000，Ti = 0.5 * Tu = 2.221441 s，Td = 0.125 * Tu = 0.555360 s这就是这门课"三章串起来"的地方：04-root 用劳斯表求出 Ku 与振荡频率，05-frequency 用频域确认稳定裕度，这一章把 Ku 与 Tu 直接喂进 Z-N 整定表。整定不需要实验台，只需要一张根轨迹图。
例 2：P / PI / PD / PID 四种组合的实测对照
被控对象同上，dt = 1 ms，仿真 200 s：
| 控制器 | 超调量 | 峰值时间 | ts(±2%) | 终值 |
|---|---|---|---|---|
P（Kp = 3.0） | 56.5708% | 3.3790 s | 22.2040 s | 1.00000000 |
PI（Kp = 2.7，Ti = 3.702402） | 96.4443% | 3.6230 s | 80.1160 s | 0.99996917 |
PD（Kp = 3.0，Td = 0.555360） | 22.4400% | 2.7440 s | 6.6460 s | 1.00000000 |
| PID（Z-N 整定） | 59.6197% | 2.6810 s | 12.9070 s | 1.00000000 |
这张表一次说清三项的贡献与代价：
- 加
I让超调从 56.57% 涨到 96.44%、ts从 22.20 s 拖到 80.12 s——"积分消差"的代价是稳定性，而且代价相当大。从ts(2%) = 80.116 s反推主导极点的实部约为4 / 80.116 ≈ 0.0499，极点几乎贴在虚轴上。 - 加
D让超调从 56.57% 降到 22.44%、ts从 22.20 s 降到 6.65 s——这是三项里唯一"全是好处"的一项。 PID是折中：超调 59.62% 在PI与PD之间，而ts = 12.91 s比PI快得多。注意 Z-N 的 PID 超调仍然接近 60%，这就是"Z-N 偏激进"的具体数字——想更保守，就把Kp乘0.6。
终值那一列两处要点：P 与 PD 的终值都精确是 1.00000000，因为对象 1/(s(s+1)(s+2)) 本身就是 1 型系统，对阶跃输入理论静差为零；PI 的终值 0.99996917 只是还没完全收敛（ts 都到 80 s 了），不是有静差。
例 3：位置式与增量式的等价性
同一组 Z-N 参数，两种实现各跑一遍：
text
两种实现的最大逐点偏差 = 1.399e-14
位置式 超调= 59.6197% tp= 2.6810s ts(2%)= 12.9070s 终值=1.00000000
增量式 超调= 59.6197% tp= 2.6810s ts(2%)= 12.9070s 终值=1.00000000指标到最后一位小数都相同，逐点偏差 1.399e-14 纯粹是浮点累加误差（约 10⁻¹⁴，相对精度 10⁻¹⁶ 级别）。这条核对的用处：确认"增量式不是另一种控制器，只是同一个控制器的另一种写法"——工程上选哪个，取决于要不要抗冲击、要不要省那点存储，而不是性能。
例 4：用 Routh 判据算出微分项的"生死线"
Kd 会不会太小？把闭环特征方程写出来（Kd > 0 时）：
排 Routh 表，第一列依次是 1、3、b1 = 0.8 + Kd、c1 = [(0.8+Kd)·Kp − 3Ki] / (0.8+Kd)、Ki = 1.620569：
Kd | b1 | c1 | 结论 |
|---|---|---|---|
| 0.000 | 0.800000 | -2.477135 | 不稳定 |
| 0.300 | 1.100000 | -0.819735 | 不稳定 |
| 0.500 | 1.300000 | -0.139775 | 不稳定 |
| 0.600 | 1.400000 | +0.127351 | 稳定 |
| 1.000 | 1.800000 | +0.899051 | 稳定 |
| 2.000 | 2.800000 | +1.863676 | 稳定 |
| 4.000 | 4.800000 | +2.587144 | 稳定 |
解 c1 = 0 得到临界值：
Kd_crit = 3·Ki/Kp − 0.8 = 3×1.620569/3.6 − 0.8 = 1.350474 − 0.8 = 0.550474这个结果是这一章最有价值的一个数字：在 Kp = 3.6、Ki = 1.620569 的前提下，Kd 必须大于 0.550474 系统才稳定——而 Z-N 给的正是 Kd = 1.999297（查表 Td = 0.125Tu 换算而来），留有约 3.6 倍的余量。这就是"Z-N 凭什么敢给一组固定系数"的答案：它的表就是围绕这类对象的稳定边界设计的。
核对一下边界两侧的实际行为：
text
Kd = 0.50（低于临界值）：超调 = 600.3340%，ts 到 200 s 仍未收敛，终值 = -2.85610703 （发散）
Kd = 0.60（略高于临界值）：超调 = 103.2575%，ts 到 200 s 仍未收敛，终值 = 0.88383003 （勉强稳定）两个都是"看起来快完蛋"的样子：一个发散、一个振荡衰减极慢。这说明 Routh 判据给出的是"稳不稳"的硬边界，而"好用不好用"还要看离边界有多远——Kd = 1.999297 距边界 0.550474 有 3.6 倍，这才是"能用"的原因。
例 5：积分饱和与抗饱和
被控对象取一阶惯性 y' = (-y + u) / 1（时间常数 1 s、增益 1），参考输入阶跃 r = 1，执行机构限幅 u 在 [0, 1.2]；控制器用 PI（Kp = 5、Ki = 5）：
| 方案 | 超调量 | ts(±2%) | 终值 |
|---|---|---|---|
| 无抗饱和 | 17.3533% | 6.2122 s | 1.000000 |
| 积分限幅 | -0.0000% | 1.7246 s | 1.000000 |
| 反算法 | -0.0000% | 1.7518 s | 1.000000 |
三组数据讲了三件事：
- 无抗饱和时超调
17.35%、调节时间6.21 s：因为饱和期间积分器仍在累加，等误差反号后还得"退债"，于是冲过头。 - 积分限幅把它压成"零超调"、
ts缩短到1.72 s——改善幅度接近四倍，而且做法只是加一行I = min(max(I, 0), 1/Ki)。 - 积分限幅与反算法结果几乎相同（
1.7246对1.7518）：两者都有效地"关掉了饱和期间的积分"，差别在工程细节——反算法不需要手调上限，适应面更宽。
注意"超调 = -0.0000%"这个写法：它是一个极小的负数（浮点残留，约 -10⁻¹⁷），物理意义是"完全没有超调"。看到 -0.0000 要读成零，不要读成"负超调"——超调量为负是没有物理意义的。
例 6：C 实现——整定表、两种实现对照与 Routh 判据
/* pid.c —— 数字 PID：临界比例度法整定、位置式/增量式对照、Routh 判据 */
#include <stdio.h>
#include <math.h>
#define PI 3.14159265358979323846
#define DT 0.001 /* 仿真步长 1 ms */
#define NSTEPS 200000 /* 共 200 s */
#define BAND 0.02 /* ±2% 带 */
/* 对象 1/(s(s+1)(s+2))：y''' + 3y'' + 2y' = u，状态取 (y, y', y'') */
static double g_x1, g_x2, g_x3;
static void plant_reset(void)
{
g_x1 = g_x2 = g_x3 = 0.0;
}
/* 前向欧拉一步：先用旧状态算导数，再统一更新 */
static void plant_step(double u)
{
double d1 = g_x2, d2 = g_x3, d3 = -2.0 * g_x2 - 3.0 * g_x3 + u;
g_x1 += DT * d1;
g_x2 += DT * d2;
g_x3 += DT * d3;
}
/* mode = 0 位置式；mode = 1 增量式 */
static void run(const char *name, double Kp, double Ki, double Kd, int mode)
{
double e1 = 0.0, e2 = 0.0, I = 0.0, u = 0.0;
double ymax = -1e300, yend = 0.0, tpeak = 0.0, tlast = 0.0;
int n;
plant_reset();
for (n = 0; n < NSTEPS; n++) {
double y = g_x1;
double e = 1.0 - y;
if (mode == 0) { /* 位置式 */
I += e * DT;
u = Kp * e + Ki * I + Kd * (e - e1) / DT;
} else { /* 增量式 */
u += Kp * (e - e1) + Ki * DT * e + Kd * (e - 2.0 * e1 + e2) / DT;
}
e2 = e1;
e1 = e;
plant_step(u);
y = g_x1;
if (y > ymax) { ymax = y; tpeak = (n + 1) * DT; }
if (fabs(y - 1.0) > BAND) tlast = (n + 1) * DT;
yend = y;
}
printf(" %-8s 超调=%9.4f%% tp=%8.4fs ts(2%%)=%9.4fs 终值=%.8f\n",
name, (ymax - 1.0) * 100.0, tpeak, tlast, yend);
}
/* 4 阶 Routh 表第一列：1, 3, b1, c1, Ki */
static void routh(double Kp, double Ki, double Kd)
{
double b1 = (3.0 * (2.0 + Kd) - Kp) / 3.0;
double c1 = (b1 * Kp - 3.0 * Ki) / b1;
printf(" Kd = %5.3f: 第一列 1, 3, %10.6f, %10.6f, %10.6f -> [%s]\n",
Kd, b1, c1, Ki, (b1 > 0.0 && c1 > 0.0) ? "稳定" : "不稳定");
}
int main(void)
{
double Ku = 6.0, wu = sqrt(2.0), Tu = 2.0 * PI / wu;
double Kp, Ti, Td, Ki, Kd;
printf("=== 一、临界比例度法（Ziegler-Nichols 闭环法）===\n");
printf(" Ku = %.4f，wu = %.6f rad/s，Tu = 2*pi/wu = %.6f s\n", Ku, wu, Tu);
printf(" P : Kp = 0.50*Ku = %.6f\n", 0.50 * Ku);
printf(" PI : Kp = 0.45*Ku = %.6f，Ti = Tu/1.2 = %.6f s\n",
0.45 * Ku, Tu / 1.2);
printf(" PID : Kp = 0.60*Ku = %.6f，Ti = 0.5*Tu = %.6f s，Td = 0.125*Tu = %.6f s\n",
0.60 * Ku, 0.5 * Tu, 0.125 * Tu);
Kp = 0.60 * Ku;
Ti = 0.5 * Tu;
Td = 0.125 * Tu;
Ki = Kp / Ti;
Kd = Kp * Td;
printf("\n=== 二、P / PI / PD / PID 对照（dt = %.3f s，%.0f s）===\n",
DT, NSTEPS * DT);
printf(" PID: Ki = Kp/Ti = %.6f，Kd = Kp*Td = %.6f\n", Ki, Kd);
run("P", 0.50 * Ku, 0.0, 0.0, 0);
run("PI", 0.45 * Ku, 0.45 * Ku / (Tu / 1.2), 0.0, 0);
run("PD", 3.0, 0.0, 3.0 * 0.125 * Tu, 0);
run("PID", Kp, Ki, Kd, 0);
printf("\n=== 三、位置式与增量式的等价性（%.0f s）===\n", NSTEPS * DT);
run("pos", Kp, Ki, Kd, 0);
run("inc", Kp, Ki, Kd, 1);
printf("\n=== 四、Routh 判据：Kd 至少要到多少才稳 ===\n");
printf(" 闭环特征方程 s^4 + 3s^3 + (2+Kd)s^2 + Kp*s + Ki = 0\n");
routh(Kp, Ki, 0.0);
routh(Kp, Ki, 0.3);
routh(Kp, Ki, 0.5);
routh(Kp, Ki, 0.6);
routh(Kp, Ki, 1.0);
routh(Kp, Ki, 2.0);
routh(Kp, Ki, 4.0);
printf(" 令 c1 = 0 -> Kd_crit = 3*Ki/Kp - 0.8 = %.6f\n", 3.0 * Ki / Kp - 0.8);
return 0;
}
c 本站为静态站，不提供在线运行；可复制到本地用 gcc / python 执行
预期输出：
=== 一、临界比例度法（Ziegler-Nichols 闭环法）===
Ku = 6.0000，wu = 1.414214 rad/s，Tu = 2*pi/wu = 4.442883 s
P : Kp = 0.50*Ku = 3.000000
PI : Kp = 0.45*Ku = 2.700000，Ti = Tu/1.2 = 3.702402 s
PID : Kp = 0.60*Ku = 3.600000，Ti = 0.5*Tu = 2.221441 s，Td = 0.125*Tu = 0.555360 s
=== 二、P / PI / PD / PID 对照（dt = 0.001 s，200 s）===
PID: Ki = Kp/Ti = 1.620569，Kd = Kp*Td = 1.999297
P 超调= 56.5708% tp= 3.3790s ts(2%)= 22.2040s 终值=1.00000000
PI 超调= 96.4443% tp= 3.6230s ts(2%)= 80.1160s 终值=0.99996917
PD 超调= 22.4400% tp= 2.7440s ts(2%)= 6.6460s 终值=1.00000000
PID 超调= 59.6197% tp= 2.6810s ts(2%)= 12.9070s 终值=1.00000000
=== 三、位置式与增量式的等价性（200 s）===
pos 超调= 59.6197% tp= 2.6810s ts(2%)= 12.9070s 终值=1.00000000
inc 超调= 59.6197% tp= 2.6810s ts(2%)= 12.9070s 终值=1.00000000
=== 四、Routh 判据：Kd 至少要到多少才稳 ===
闭环特征方程 s^4 + 3s^3 + (2+Kd)s^2 + Kp*s + Ki = 0
Kd = 0.000: 第一列 1, 3, 0.800000, -2.477135, 1.620569 -> [不稳定]
Kd = 0.300: 第一列 1, 3, 1.100000, -0.819735, 1.620569 -> [不稳定]
Kd = 0.500: 第一列 1, 3, 1.300000, -0.139775, 1.620569 -> [不稳定]
Kd = 0.600: 第一列 1, 3, 1.400000, 0.127351, 1.620569 -> [稳定]
Kd = 1.000: 第一列 1, 3, 1.800000, 0.899051, 1.620569 -> [稳定]
Kd = 2.000: 第一列 1, 3, 2.800000, 1.863676, 1.620569 -> [稳定]
Kd = 4.000: 第一列 1, 3, 4.800000, 2.587144, 1.620569 -> [稳定]
令 c1 = 0 -> Kd_crit = 3*Ki/Kp - 0.8 = 0.550474"位置式与增量式结果完全一致"这段代码是刻意保留的：它在两个分支里分别写了四种不同的更新方式（u = ... 与 u += ...），如果不一致，说明其中一个写错了。这就是"等价性核对"作为验证手段的价值——它不需要参考答案，因为它自己就是自己的对照。
注意 plant_step 的写法：
c
double d1 = g_x2, d2 = g_x3, d3 = -2.0 * g_x2 - 3.0 * g_x3 + u;
g_x1 += DT * d1; g_x2 += DT * d2; g_x3 += DT * d3;先把三个导数全部算完，再统一更新三个状态——这是前向欧拉的正确写法。如果写成"更新 x1 之后再用新的 x1 去算 x3"，就不是同一个积分器了，虽然只差 O(dt²)，但两个程序给出的数字会对不上（本机无 gcc，此段已用等价的 Python 实现逐位核对）。
例 7：Python——整定、四种组合、等价性、Routh 与抗饱和
# PID：临界比例度法整定 + P/PI/PD/PID 对照 + 抗饱和
import math
DT = 1e-3 # 仿真步长 1 ms
def sim(Kp, Ki, Kd, T, form="position"):
"""对象 1/(s(s+1)(s+2))，写成 y'''+3y''+2y' = u，状态取 (y, y', y'')"""
x = [0.0, 0.0, 0.0]
e1 = e2 = 0.0
I = 0.0
u = 0.0
ys = []
for _ in range(int(T / DT)):
y = x[0]
e = 1.0 - y
if form == "position":
I += e * DT
u = Kp * e + Ki * I + Kd * (e - e1) / DT
else: # 增量式
u += Kp * (e - e1) + Ki * DT * e + Kd * (e - 2 * e1 + e2) / DT
e2, e1 = e1, e
d1, d2, d3 = x[1], x[2], -2.0 * x[1] - 3.0 * x[2] + u
x[0] += DT * d1
x[1] += DT * d2
x[2] += DT * d3
ys.append(x[0])
return ys
def report(name, ys):
peak = max(ys)
tp = (ys.index(peak) + 1) * DT
ts = 0.0
for i in range(len(ys) - 1, -1, -1):
if abs(ys[i] - 1.0) > 0.02:
ts = (i + 1) * DT
break
print(f" {name:12s} 超调={100.0 * (peak - 1.0):9.4f}% tp={tp:7.4f}s"
f" ts(2%)={ts:8.4f}s 终值={ys[-1]:.8f}")
print("=== 一、临界比例度法（Ziegler-Nichols 闭环法）===")
Ku, wu = 6.0, math.sqrt(2.0)
Tu = 2.0 * math.pi / wu
print(f" 由 04-root 的劳斯表已知：临界增益 Ku = {Ku:.4f}，振荡角频率 wu = {wu:.6f} rad/s")
print(f" 临界振荡周期 Tu = 2*pi/wu = {Tu:.6f} s")
print(" 整定表（Kp 取自 Ku，Ti / Td 取自 Tu）：")
print(f" P : Kp = 0.50*Ku = {0.50 * Ku:.6f}")
print(f" PI : Kp = 0.45*Ku = {0.45 * Ku:.6f}，Ti = Tu/1.2 = {Tu / 1.2:.6f} s")
print(f" PID : Kp = 0.60*Ku = {0.60 * Ku:.6f}，Ti = 0.5*Tu = {0.5 * Tu:.6f} s，"
f"Td = 0.125*Tu = {0.125 * Tu:.6f} s")
print(f"\n=== 二、P / PI / PD / PID 对照（dt = {DT:.0e} s）===")
Kp_pid, Ti_pid, Td_pid = 0.60 * Ku, 0.5 * Tu, 0.125 * Tu
Ki_pid = Kp_pid / Ti_pid
Kd_pid = Kp_pid * Td_pid
print(f" PID 标准型 C(s) = Kp(1 + 1/(Ti*s) + Td*s)"
f"，故 Ki = Kp/Ti = {Ki_pid:.6f}，Kd = Kp*Td = {Kd_pid:.6f}")
report("P", sim(0.50 * Ku, 0.0, 0.0, 200.0))
report("PI", sim(0.45 * Ku, 0.45 * Ku / (Tu / 1.2), 0.0, 200.0))
report("PD", sim(3.0, 0.0, 3.0 * 0.125 * Tu, 200.0))
report("PID", sim(Kp_pid, Ki_pid, Kd_pid, 200.0))
print("\n=== 三、位置式与增量式的等价性 ===")
a = sim(Kp_pid, Ki_pid, Kd_pid, 200.0, "position")
b = sim(Kp_pid, Ki_pid, Kd_pid, 200.0, "incremental")
print(f" 两种实现的最大逐点偏差 = {max(abs(x - y) for x, y in zip(a, b)):.3e}"
f"（仅累积舍入差异）")
report("位置式", a)
report("增量式", b)
print("\n=== 四、Routh 判据：Kd 至少要到多少才稳 ===")
print(" Kd > 0 时闭环特征方程为 s^4 + 3s^3 + (2+Kd)s^2 + Kp*s + Ki = 0")
print(" Routh 表第一列：1, 3, 0.8+Kd, c1, Ki，其中 c1 = ((0.8+Kd)*Kp - 3*Ki)/(0.8+Kd)")
def routh_c1(Kp, Ki, Kd):
b1 = (3.0 * (2.0 + Kd) - 1.0 * Kp) / 3.0
return b1, (b1 * Kp - 3.0 * Ki) / b1
for Kd in (0.0, 0.3, 0.5, 0.6, 1.0, 2.0, 4.0, 8.0):
b1, c1 = routh_c1(Kp_pid, Ki_pid, Kd)
ok = "稳定" if (b1 > 0 and c1 > 0) else "不稳定"
print(f" Kd = {Kd:5.3f}: 第一列 1, 3, {b1:.6f}, {c1:+.6f}, {Ki_pid:.6f} -> [{ok}]")
Kd_crit = (3.0 * Ki_pid / Kp_pid) - 0.8
print(f" 令 c1 = 0 解出临界值 Kd_crit = 3*Ki/Kp - 0.8 = {Kd_crit:.6f}")
print(" 核对：Kd 略小于临界值会怎样 ——")
report("Kd=0.50", sim(Kp_pid, Ki_pid, 0.50, 200.0))
report("Kd=0.60", sim(Kp_pid, Ki_pid, 0.60, 200.0))
print("\n=== 五、积分饱和与抗饱和 ===")
print(" 对象 y' = (-y+u)/1，r = 1，u 限幅 [0, 1.2]；PI 控制器 Kp = 5，Ki = 5，dt = 0.2 ms")
def sim_aw(mode, dt=2e-4, T=30.0):
y, I, h = 0.0, 0.0, []
Kp_, Ki_, umax = 5.0, 5.0, 1.2
for _ in range(int(T / dt)):
e = 1.0 - y
Iraw = I + e * dt
if mode == "clamp":
I = min(max(Iraw, 0.0), 1.0 / Ki_) # 限幅：Ki*I 不超过 1
else:
I = Iraw
u = Kp_ * e + Ki_ * I
us = min(max(u, 0.0), umax)
if mode == "backcalc":
I = Iraw + dt * (us - u) / Kp_ # 反算法
y += (-y + us) * dt
h.append(y)
return h
def rep_aw(name, h, dt=2e-4):
peak = max(h)
ts = 0.0
for i in range(len(h) - 1, -1, -1):
if abs(h[i] - 1.0) > 0.02:
ts = (i + 1) * dt
break
print(f" {name:12s} 超调={100.0 * (peak - 1.0):8.4f}% ts(2%)={ts:7.4f}s 终值={h[-1]:.6f}")
for m, nm in (("none", "无抗饱和"), ("clamp", "积分限幅"), ("backcalc", "反算法")):
rep_aw(nm, sim_aw(m))
python 本站为静态站，不提供在线运行；可复制到本地用 gcc / python 执行
预期输出：
=== 一、临界比例度法（Ziegler-Nichols 闭环法）===
由 04-root 的劳斯表已知：临界增益 Ku = 6.0000，振荡角频率 wu = 1.414214 rad/s
临界振荡周期 Tu = 2*pi/wu = 4.442883 s
整定表（Kp 取自 Ku，Ti / Td 取自 Tu）：
P : Kp = 0.50*Ku = 3.000000
PI : Kp = 0.45*Ku = 2.700000，Ti = Tu/1.2 = 3.702402 s
PID : Kp = 0.60*Ku = 3.600000，Ti = 0.5*Tu = 2.221441 s，Td = 0.125*Tu = 0.555360 s
=== 二、P / PI / PD / PID 对照（dt = 1e-03 s）===
PID 标准型 C(s) = Kp(1 + 1/(Ti*s) + Td*s)，故 Ki = Kp/Ti = 1.620569，Kd = Kp*Td = 1.999297
P 超调= 56.5708% tp= 3.3790s ts(2%)= 22.2040s 终值=1.00000000
PI 超调= 96.4443% tp= 3.6230s ts(2%)= 80.1160s 终值=0.99996917
PD 超调= 22.4400% tp= 2.7440s ts(2%)= 6.6460s 终值=1.00000000
PID 超调= 59.6197% tp= 2.6810s ts(2%)= 12.9070s 终值=1.00000000
=== 三、位置式与增量式的等价性 ===
两种实现的最大逐点偏差 = 1.399e-14（仅累积舍入差异）
位置式 超调= 59.6197% tp= 2.6810s ts(2%)= 12.9070s 终值=1.00000000
增量式 超调= 59.6197% tp= 2.6810s ts(2%)= 12.9070s 终值=1.00000000
=== 四、Routh 判据：Kd 至少要到多少才稳 ===
Kd > 0 时闭环特征方程为 s^4 + 3s^3 + (2+Kd)s^2 + Kp*s + Ki = 0
Routh 表第一列：1, 3, 0.8+Kd, c1, Ki，其中 c1 = ((0.8+Kd)*Kp - 3*Ki)/(0.8+Kd)
Kd = 0.000: 第一列 1, 3, 0.800000, -2.477135, 1.620569 -> [不稳定]
Kd = 0.300: 第一列 1, 3, 1.100000, -0.819735, 1.620569 -> [不稳定]
Kd = 0.500: 第一列 1, 3, 1.300000, -0.139775, 1.620569 -> [不稳定]
Kd = 0.600: 第一列 1, 3, 1.400000, +0.127351, 1.620569 -> [稳定]
Kd = 1.000: 第一列 1, 3, 1.800000, +0.899051, 1.620569 -> [稳定]
Kd = 2.000: 第一列 1, 3, 2.800000, +1.863676, 1.620569 -> [稳定]
Kd = 4.000: 第一列 1, 3, 4.800000, +2.587144, 1.620569 -> [稳定]
Kd = 8.000: 第一列 1, 3, 8.800000, +3.047533, 1.620569 -> [稳定]
令 c1 = 0 解出临界值 Kd_crit = 3*Ki/Kp - 0.8 = 0.550474
核对：Kd 略小于临界值会怎样 ——
Kd=0.50 超调= 600.3340% tp=197.9520s ts(2%)=200.0000s 终值=-2.85610703
Kd=0.60 超调= 103.2575% tp= 3.0210s ts(2%)=200.0000s 终值=0.88383003
=== 五、积分饱和与抗饱和 ===
对象 y' = (-y+u)/1，r = 1，u 限幅 [0, 1.2]；PI 控制器 Kp = 5，Ki = 5，dt = 0.2 ms
无抗饱和 超调= 17.3533% ts(2%)= 6.2122s 终值=1.000000
积分限幅 超调= -0.0000% ts(2%)= 1.7246s 终值=1.000000
反算法 超调= -0.0000% ts(2%)= 1.7518s 终值=1.000000第五段里 Kd=0.50 的 ts(2%) = 200.0000 s 要正确理解：200 s 正好是仿真的总时长，意思是"到仿真结束都没进入 ±2% 带"——结合终值 -2.85610703，这是发散，不是"调节时间刚好 200 秒"。ts 等于仿真时长永远是一个信号：该延长仿真时间，或者系统根本不收敛。
考点
- 三项分工：
P管快（力度）、I管准（精度）、D管稳（阻尼）。P是有差调节，I无差调节。 - 两种参数写法的换算：
Ki = Kp/Ti、Kd = Kp·Td；反向Ti = Kp/Ki、Td = Kd/Kp。工程整定表给的是Kp/Ti/Td。 - Z-N 临界比例度法：先只加
P加到等幅振荡得Ku与Tu，再查表 ——P: 0.5Ku；PI: 0.45Ku、Tu/1.2；PID: 0.6Ku、0.5Tu、0.125Tu。 - Z-N 阶跃响应法：从阶跃曲线读
K、L、T，P: T/(KL)；PI: 0.9T/(KL)、3.3L；PID: 1.2T/(KL)、2L、0.5L。 Ku就是根轨迹与虚轴交点对应的增益——所以整定参数可以纯计算得出，不必做振荡实验。- Z-N 参数偏激进（设计目标是四分之一衰减，超调常在 25% 以上）；要更保守就把
Kp乘0.5~0.8。 - 位置式与增量式数学等价，差别只在数值实现：增量式在手动/自动切换时无冲击、积分饱和影响小。
- 积分饱和的机理：执行机构饱和期间误差仍在，积分器持续累加，等误差反号后要"退债"，于是造成巨大超调、恢复缓慢。
- 三种抗饱和：积分限幅（最简单）、反算法（
I += dt·(u_sat − u)/Kp，自整定）、条件积分（饱和期间暂停积分）。 D对纯延迟无力（延迟的相位滞后随频率线性增长，D的+90度补不回来），且放大高频噪声——所以大延迟对象常用PI。- 易错：把
Ki当成"积分时间"（积分时间是Ti，Ki是它的倒数乘以Kp）；把Ti增大说成"加强积分"（Ti越大积分越弱）；Routh 判据只判"稳不稳"而不看"离边界多远"（Kd = 0.6只比临界值大9%，虽稳但振荡衰减极慢）；把ts等于仿真时长当成"调节时间正好那么长"；把-0.0000%当成负超调。
小结
- PID = 误差的"现在 + 过去 + 未来"三项加权：
u = Kp·e + Ki·∫e dt + Kd·de/dt。 - 三项各有分工也各有代价：
P快但有差、I准但吃裕度、D稳但怕噪声。 - 用前两章的工具可以直接整定：对
1/(s(s+1)(s+2))，Ku = 6、Tu = 4.442883 s，Z-N 给出Kp = 3.6、Ti = 2.221441 s、Td = 0.555360 s。 - 实测对照说明三项各自的分量：加
I让超调由 56.57% 涨到 96.44%、ts由 22.20 s 拖到 80.12 s；加D把它降到 22.44% 与 6.65 s；PID折中在 59.62% 与 12.91 s。 - 位置式与增量式数学等价（逐点偏差
1.4×10⁻¹⁴），差别只在数值实现与切换平滑性。 - Routh 给出
Kd的生死线：Kd_crit = 3Ki/Kp − 0.8 = 0.550474，Z-N 给的1.999297留了约 3.6 倍余量。 - 积分饱和是落地头号事故，积分限幅能把超调从 17.35% 压到 0、
ts从 6.21 s 缩到 1.72 s。 - 整定顺序不能颠倒：先
P定力度，再I消差，最后D修形状。
回到主线：从 04-root 到本章，走完的是经典控制的完整闭环——由根轨迹求临界增益、由频域确认裕度、由 PID 补回裕度并定参数。这套方法的极限在哪里？ 它处理的始终是单输入单输出：一个指令、一个反馈、一个执行器。
下一章换一套语言：不再问"某个参数怎么变"，而是把系统的全部内部状态写成向量，直接问"这些状态能不能被控制、能不能被观测"——这就是状态空间。
下一篇：状态空间分析
评论（0）
当前浏览器不允许本地存储，评论无法保存。
还没有评论，来说两句。