企业官网建设流程全解析

1. 为什么搞控制的最后都要回到差分方程

Z变换这门东西，在学校里是被当纯数学工具教的：定义、收敛域、性质、反变换，考完试就还给老师了。可一旦你真的去写控制代码，就会发现根本躲不开它。你在 Matlab 或者纸上手推出一个 G(z)，看着挺漂亮，但烧进单片机的那段程序里，压根没有 z 这个变量，也没有 z^{-1} 这个算子，只有一个 while(1) 主循环和几个 float 变量。Z变换方程和差分方程之间那一次"翻译"，就是纸面设计和现场代码之间最容易出岔子的地方。我见过太多人 Z 传递函数推得一点没错，代码写出来响应就是不对，波形要么发散要么迟钝，追根溯源，全卡在这次翻译上。

这篇东西就是把这个翻译过程掰开揉碎讲一遍。适合谁看？写过 PID 但没系统推过离散化的、被"Z变换转差分方程"卡过一晚上的、以及想搞清楚"位置式 PID 用差分方程到底怎么落地"的人。核心就一件事：拿到一个 Z 传递函数，怎么一步一步变成可以在代码里跑的递推式，中间有哪几个参数必须算准，哪几个坑必须绕开。不堆公式，但我保证每一步都能自己推出来。

1.1 z^{-1} 的本质是一个采样周期的记忆

整个转换的地基只有一条性质，就是移序定理。对于采样序列 x(n)，在零初始条件下：

$$ Z{x(n-1)} = z^{-1} X(z) $$

$$ Z{x(n-2)} = z^{-2} X(z) $$

翻译成人话：z^{-1} 不是某种抽象算子，它的物理含义就是"延迟一个采样周期"。你在程序里想用得上一拍的值，就是 z^{-1}；想用得上前两拍的值，就是 z^{-2}。仅此而已。这条性质记住，后面所有的推导都只是这条定理的反复使用。

这里有个前提很多人会漏掉：移序定理成立需要初始条件为零，也就是 x(-1) = x(-2) = ... = 0。实际项目里，系统上电那一刻所有历史状态本来就该清零，所以这个前提天然满足。但如果你是从某个已知的中间状态开始仿真，比如做在线切换控制器，那就必须把这部分初值补回去，否则递推出来的前几拍会跟真实值差一截。这一点我在做控制器在线切换的时候吃过亏，切换瞬间输出跳一下，查了半天才发现是历史状态没对齐。

采样序列 x(n) 里的 n 代表第 n 个采样点，对应的真实时间是 nT，T 是采样周期。所以严格写应该是 x(nT)，但工程习惯上省掉 T，写成 x(n)。你看到 x(n-1)，脑子里要自动翻译成"上一个采样时刻的值"。

1.2 从 s 域到 z 域，中间隔着一个采样开关

在讲转换之前，得先把位置摆正。你手上的 G(z) 是从哪来的？它不是凭空冒出来的，是从连续被控对象 G(s) 经过采样和保持之后得到的。这条链路是：连续对象 G(s) → 采样开关（周期 T）→ 零阶保持器 ZOH → 离散化 → G(z)。

很多人推差分方程推错，根本原因不是数学差，而是搞混了自己手上的 G(z) 到底是"控制器的 Z 传递函数"还是"被控对象离散化之后的 Z 传递函数"。这两者的物理含义完全不同，差分方程里的系数含义也完全不同。控制器是你自己设计的，被控对象是客观存在的，你只能去逼近它。

零阶保持器这一步的物理意义值得说清楚。数字控制器算出一拍输出之后，这个值在整个采样周期内保持不变，直到下一拍更新，这就是零阶保持。它的传递函数是：

$$ G_{ZOH}(s) = \frac{1 - e^{-sT}}{s} $$

所以被控对象从 s 域到 z 域的精确映射是：

$$ G(z) = (1 - z^{-1}) \cdot Z\left{\frac{G(s)}{s}\right} $$

注意这里 Z{·} 里面除了一个 s，这个 1/s 就是保持器带来的积分效应。忘了除这个 s，是新手推 ZOH 离散化最常见的错误。记住这个公式，后面一阶惯性环节那一节我会完整走一遍。

2. 拿到一个 Z 传递函数，三条路把它拆成差分方程

先约定一个通用形式。工程上遇到的 Z 传递函数，绝大多数都能写成下面这种分子分母都是 z^{-1} 多项式的形式：

$$ G(z) = \frac{Y(z)}{U(z)} = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2} + \cdots + b_m z^{-m}}{1 + a_1 z^{-1} + a_2 z^{-2} + \cdots + a_n z^{-n}} $$

这里我把分母的首项强制归一化成 1，这是个好习惯。如果原始的 G(z) 分母首项不是 1，先除一下，不然后面交叉相乘会多出一堆没必要的除法，代码里也容易写错整除。分子最高阶留了 m 次，分母留了 n 次，工程上一般 n ≥ m，否则会出现输入超前于输出的非因果系统，物理上实现不了。

从上面这个形式出发，有三条路可以走到差分方程。三条路得到的最终结果一定一致，但适用场景和计算复杂度不同。

2.1 直接交叉相乘法：最省事也最不容易错

最直接的办法，是把 Y(z)/U(z) = 分子/分母 这个等式交叉相乘，把分母乘到左边：

$$ (1 + a_1 z^{-1} + a_2 z^{-2} + \cdots) Y(z) = (b_0 + b_1 z^{-1} + b_2 z^{-2} + \cdots) U(z) $$

然后用移序定理，把每一个 z^{-k} 对应的项翻译成时域的 x(n-k)，直接得到：

$$ y(n) = -\sum_{i=1}^{n} a_i y(n-i) + \sum_{j=0}^{m} b_j u(n-j) $$

这就是差分方程最原始的形态，也是我推荐第一优先使用的办法。原因有三点：第一，不需要记忆任何反变换表；第二，整个过程是机械的，几乎不可能出错；第三，得到的系数直接就是代码里要用的系数，一一对应，不用二次换算。

举个具体的例子。设：

$$ G(z) = \frac{z}{(z-0.5)(z-0.8)} $$

先把分母乘出来：$(z-0.5)(z-0.8) = z^2 - 1.3z + 0.4$。于是：

$$ G(z) = \frac{z}{z^2 - 1.3z + 0.4} $$

分子分母同时除以 z^2，把首项归一化成 1：

$$ G(z) = \frac{z^{-1}}{1 - 1.3z^{-1} + 0.4z^{-2}} $$

交叉相乘：

$$ (1 - 1.3z^{-1} + 0.4z^{-2}) Y(z) = z^{-1} U(z) $$

移序翻译：

$$ y(n) = 1.3 y(n-1) - 0.4 y(n-2) + u(n-1) $$

到这一步就完成了。你可以立刻拿它去写代码，三个系数分别是 1.3、-0.4、1，一目了然。

提示：交叉相乘之后，所有 y 相关项移到左边，u 相关项留在右边，然后用移序定理。别把两边搞反，搞反了就是"求输入"而不是"求输出"了。

2.2 部分分式展开加查表：适合手推和验证

当你想知道系统的解析行为、想看极点分布、想在纸上快速估个响应形状的时候，部分分式法比直接法更直观。

做法是先把 G(z)/z 展开成部分分式（注意，是 G(z) 除以 z，不是 G(z) 本身），再乘回来，然后查常见 Z 变换对。为什么是除以 z？因为标准的 Z 变换表里，$(z/(z-p))$ 对应的时域序列是 $p^n$，先除 z 是为了让每一项都能配成这种标准形式。

还是用上面那个例子：

$$ \frac{G(z)}{z} = \frac{1}{(z-0.5)(z-0.8)} $$

设展开为：

$$ \frac{1}{(z-0.5)(z-0.8)} = \frac{A}{z-0.5} + \frac{B}{z-0.8} $$

用留数法求系数：

$$ A = \left.(z-0.5)\frac{1}{(z-0.5)(z-0.8)}\right|_{z=0.5} = \frac{1}{0.5-0.8} = -\frac{1}{0.3} \approx -3.3333 $$

$$ B = \frac{1}{0.8-0.5} = \frac{1}{0.3} \approx 3.3333 $$

所以：

$$ G(z) = -\frac{3.3333z}{z-0.5} + \frac{3.3333z}{z-0.8} $$

查表反变换，得到脉冲响应序列：

$$ g(n) = -3.3333 \times 0.5^n + 3.3333 \times 0.8^n,\quad n \geq 0 $$

这个结果的含义很清楚：系统的极点就在 0.5 和 0.8 两个位置，脉冲响应是两个指数衰减分量的叠加。数值验证一下：g(0) = -3.3333 + 3.3333 = 0，g(1) = -1.6667 + 2.6667 = 1，g(2) = -0.8333 + 2.1333 = 1.3。跟前面直接法递推出来的 y(1)=1、y(2)=1.3 完全对上。两条路殊途同归，这也是我用来交叉验证的手段。

部分分式法的价值在于，它把系统行为拆成了一个个一阶模态。你能一眼看出哪个极点主导响应、哪个衰减得快。0.5 那个极点衰减得比 0.8 快，所以稳态附近是 0.8 那个模态说了算。这种洞察在调参数的时候很有用。

2.3 长除法取幂级数：直接给你一段脉冲响应

第三条路是把 G(z) 按 z^{-1} 的幂次展开成无穷级数：

$$ G(z) = g_0 + g_1 z^{-1} + g_2 z^{-2} + g_3 z^{-3} + \cdots $$

展开之后，每个系数 g_k 就是脉冲响应的第 k 个采样值。原因是输入为单位脉冲 U(z) = 1 时，Y(z) = G(z)，而 Y(z) 的系数按定义就是输出的采样序列。

还是那个例子：

$$ G(z) = \frac{z^{-1}}{1 - 1.3z^{-1} + 0.4z^{-2}} $$

记 x = z^{-1}，做长除：

$$ \frac{x}{1 - 1.3x + 0.4x^2} = x + 1.3x^2 + 1.29x^3 + \cdots $$

系数依次是 1、1.3、1.29，对应 g(1)=1、g(2)=1.3、g(3)=1.29，跟前面两个方法得到的完全一致。

长除法适合干两件事：一是快速验证你推导的差分方程对不对，跑几拍看看数值是不是这个序列；二是当你只想关心有限拍响应、或者要在 FPGA 里做纯流水线实现时，直接展开成有限项级数往往比递归更省资源。

注意：长除法得到的幂级数是无穷级数，截断就代表近似。截断长度不够，会有明显的稳态误差；截断太长，又浪费算力。一般截到幅值降到峰值的千分之一以下就够用。

三条路的分工我总结一下：日常写代码用直接交叉相乘法，纸上分析用部分分式法，快速验证用长除法。三者交叉使用，基本不会出错。

3. 一个一阶惯性环节，三种离散化方法的差分方程对比

讲理论容易飘，拿一个现场最最常见的对象走一遍：一阶惯性环节。

$$ G(s) = \frac{1}{\tau s + 1} $$

这个模型覆盖了电机转速、温度、液位、电池充放电等一大票场景。它的三种常见离散化方式——零阶保持器法、后向差分法、双线性变换法——得到的差分方程系数不一样，实际响应也不一样。搞清楚这三者的差异，是避免"仿真对、实机错"的关键。

3.1 零阶保持器法：采保电路的精确等效

零阶保持器法（ZOH）是在采样点上精确成立的方法。它的物理前提是：真实的 D/A 输出加上被控对象，确实构成了一个零阶保持链路。这是数字控制里最符合实际物理过程的一种离散化方式。

套用公式：

$$ G(z) = (1 - z^{-1}) \cdot Z\left{\frac{G(s)}{s}\right} $$

先算 G(s)/s：

$$ \frac{G(s)}{s} = \frac{1}{s(\tau s + 1)} = \frac{1}{s} - \frac{\tau}{\tau s + 1} = \frac{1}{s} - \frac{1}{s + 1/\tau} $$

查表做 Z 变换：

$$ Z\left{\frac{1}{s}\right} = \frac{z}{z-1},\quad Z\left{\frac{1}{s + 1/\tau}\right} = \frac{z}{z - e^{-T/\tau}} $$

代入：

$$ Z\left{\frac{G(s)}{s}\right} = \frac{z}{z-1} - \frac{z}{z - e^{-T/\tau}} $$

再乘保持器因子 $(1 - z^{-1}) = (z-1)/z$：

$$ G(z) = \frac{z-1}{z}\left(\frac{z}{z-1} - \frac{z}{z-a}\right) = 1 - \frac{z-1}{z-a} $$

其中 $a = e^{-T/\tau}$。通分整理：

$$ G(z) = \frac{z-a - z + 1}{z-a} = \frac{1-a}{z-a} = \frac{(1-a)z^{-1}}{1 - a z^{-1}} $$

交叉相乘得到差分方程：

$$ y(n) = a \cdot y(n-1) + (1-a) \cdot u(n-1),\quad a = e^{-T/\tau} $$

这个形式非常干净：只有一个系数 $a$，另一个系数是 $1-a$，两者之和恒等于 1，保证直流增益为 1。参数 $\tau$ 越大（惯性越大），$a$ 越接近 1，系统响应越慢，符合直觉。

3.2 后向差分与双线性变换：现场最常用的两种近似

后向差分法，是把 s 用 $(z-1)/(Tz)$ 替换；双线性变换（也叫 Tustin 变换），是把 s 用 $\frac{2}{T}\cdot\frac{z-1}{z+1}$ 替换。两者都是近似方法，都不是采样点精确，但在工程上极其常用，因为它们的变形幅度可控、稳定性好。

先看后向差分。代入 $s = (z-1)/(Tz)$：

$$ G(z) = \frac{1}{\tau \frac{z-1}{Tz} + 1} = \frac{Tz}{\tau(z-1) + Tz} = \frac{Tz}{(\tau + T)z - \tau} $$

分子分母除以 $(\tau+T)z$，归一化：

$$ G(z) = \frac{\frac{T}{\tau+T}}{1 - \frac{\tau}{\tau+T} z^{-1}} $$

差分方程：

$$ y(n) = \frac{\tau}{\tau+T} y(n-1) + \frac{T}{\tau+T} u(n) $$

注意这里跟 ZOH 有个关键差别：ZOH 用的是 $u(n-1)$，后向差分用的是 $u(n)$，也就是当前拍输入直接进来。这会让后向差分的响应"快一拍"，瞬时不那么平滑，但好处是不需要额外记输入历史。

再看双线性变换。代入 $s = \frac{2}{T}\cdot\frac{z-1}{z+1}$，记 $k = 2\tau/T$：

$$ G(z) = \frac{1}{\tau \cdot \frac{2}{T}\frac{z-1}{z+1} + 1} = \frac{z+1}{k(z-1) + (z+1)} = \frac{z+1}{(k+1)z + (1-k)} $$

归一化：

$$ G(z) = \frac{1}{1+k} \cdot \frac{1 + z^{-1}}{1 + \frac{1-k}{1+k} z^{-1}} $$

差分方程：

$$ y(n) = -\frac{1-k}{1+k} y(n-1) + \frac{1}{1+k}\left[u(n) + u(n-1)\right] $$

化简系数，$-\frac{1-k}{1+k} = \frac{k-1}{k+1} = \frac{2\tau - T}{2\tau + T}$，$\frac{1}{1+k} = \frac{T}{2\tau+T}$。所以：

$$ y(n) = \frac{2\tau - T}{2\tau + T} y(n-1) + \frac{T}{2\tau + T}\left[u(n) + u(n-1)\right] $$

双线性变换的特点是分子里同时出现了 $u(n)$ 和 $u(n-1)$，它相当于对输入做了一个梯形积分近似，所以又叫梯形法。它对高频段的频率畸变比较严重，但低频段保真度高，稳定性好，是被控对象离散化里非常受欢迎的方法。

3.3 三组系数代进相同数据，结果差多少

光看系数不够直观，代数值算一遍。取 $\tau = 1$、$T = 0.1$，三种方法的系数如下：

T 小的时候三者接近，误差在千分之几的量级，随便挑哪个都不影响项目。但把采样周期放大到 T = 1 再看：

差距一下就出来了。同一个连续对象，采样周期同样是 1 秒，ZOH 认为第一拍就该冲到 0.63，后向差分只冲到 0.5，双线性冲到 0.667。这三个都是"对"的，因为它们各自的近似准则不同。你在仿真里用哪套系数，实机上就必须用哪套，混着用必然出问题。我遇到过有人 Matlab 里用 ZOH 离散化算增益，C 代码里照着后向差分写递推，结果闭环增益对不上，排查了两天才发现是这里。

还有个更隐蔽的坑：双线性变换法在 T 很大时，分子里的系数会出现负值。比如 $T > 2\tau$ 时，$2\tau - T$ 变成负数，极点跑到负半轴，阶跃响应会出现明显的振荡（虽然是衰减的）。这不是错，是这个方法本身的特性，但不理解的人会以为程序写坏了。

4. 位置式 PID 的差分方程到底是怎么来的

前面讲的都是被控对象。现在换到控制器这边，题目里热词点到的"位置式 PID 用差分方程"是个高频问题，我把完整推导走一遍。

位置式 PID 的连续形式是：

$$ u(t) = K_p e(t) + K_i \int_0^t e(\tau),d\tau + K_d \frac{de(t)}{dt} $$

其中 $K_i = K_p / T_i$，$K_d = K_p T_d$，$T_i$ 是积分时间，$T_d$ 是微分时间。这三个参数在整定时是直接给的。

4.1 先把连续 PID 写成 Z 传递函数

离散化有三个算子要选：积分用矩形近似还是梯形近似，微分用后向差分还是前向差分。现场最常见的组合是积分用矩形法、微分用后向差分，因为计算量小、代码短。

积分项离散化（矩形法，累加）：

$$ \int_0^t e(\tau),d\tau \approx T \sum_{j=0}^{n} e(j) $$

它的 Z 变换是 $T \cdot E(z)/(1 - z^{-1})$，因为累加这个操作在 z 域里就是一个 $1/(1-z^{-1})$ 环节。

微分项离散化（后向差分）：

$$ \frac{de(t)}{dt} \approx \frac{e(n) - e(n-1)}{T} $$

它的 Z 变换是 $\frac{1 - z^{-1}}{T} E(z)$。

把这三项加起来，写成 U(z)/E(z)：

$$ G_c(z) = K_p + \frac{K_i T}{1 - z^{-1}} + \frac{K_d}{T}(1 - z^{-1}) $$

这就是位置式 PID 的 Z 传递函数。注意这里 G_c(z) 是控制器，输入是误差 e，输出是控制量 u。

4.2 交叉相乘得到位置式递推式

现在用第 2 节讲的直接交叉相乘法。先通分：

$$ G_c(z) = \frac{K_p(1-z^{-1}) + K_i T + \frac{K_d}{T}(1-z^{-1})^2}{1 - z^{-1}} $$

把分子展开，$(1-z^{-1})^2 = 1 - 2z^{-1} + z^{-2}$：

$$ \text{分子} = K_p - K_p z^{-1} + K_i T + \frac{K_d}{T} - \frac{2K_d}{T}z^{-1} + \frac{K_d}{T}z^{-2} $$

按 z 的幂次合并：

$$ \text{分子} = \left(K_p + K_i T + \frac{K_d}{T}\right) + \left(-K_p - \frac{2K_d}{T}\right)z^{-1} + \left(\frac{K_d}{T}\right)z^{-2} $$

记三个系数：

$$ A_0 = K_p + K_i T + \frac{K_d}{T} $$

$$ A_1 = -K_p - \frac{2K_d}{T} $$

$$ A_2 = \frac{K_d}{T} $$

于是：

$$ (1 - z^{-1})U(z) = (A_0 + A_1 z^{-1} + A_2 z^{-2})E(z) $$

移序翻译成差分方程：

$$ u(n) - u(n-1) = A_0 e(n) + A_1 e(n-1) + A_2 e(n-2) $$

整理成位置式：

$$ u(n) = u(n-1) + A_0 e(n) + A_1 e(n-1) + A_2 e(n-2) $$

到这一步，位置式 PID 的递推形式就彻底出来了。它看起来是个"增量式"，但算的其实是绝对位置 u(n)，所以叫位置式递推。你仔细看，它是用"上一拍输出 + 增量"的方式实现的，本质上还是位置式。

如果换成大家更熟悉的系数写法，把 $K_i = K_p/T_i$、$K_d = K_p T_d$ 代回去：

$$ A_0 = K_p\left(1 + \frac{T}{T_i} + \frac{T_d}{T}\right) $$

$$ A_1 = -K_p\left(1 + \frac{2T_d}{T}\right) $$

$$ A_2 = K_p \frac{T_d}{T} $$

这三个系数直接就是代码里的三个常量，整定完之后算一次，运行时不用重复计算。验证一下稳态特性：当误差恒定时 $e(n)=e(n-1)=e(n-2)=e$，增量等于 $(A_0+A_1+A_2)e$，而 $A_0+A_1+A_2 = K_p T/T_i = K_i T$，正好是每拍积分累积量，逻辑自洽。

4.3 位置式与增量式的换算关系和各自适用场景

很多人搞不清"位置式 PID 的差分方程"和"增量式 PID"到底是不是一回事。答案：它们是同一个控制器在不同状态变量下的两种实现。增量式 PID 的输出是 $\Delta u(n) = u(n) - u(n-1)$，同一条差分方程，换个写法而已。

增量式更常用在带积分饱和风险的场合，因为它的输出是增量，天然可以通过限幅增量来防止饱和。位置式更适合输出直接对应执行器绝对位置的场合，比如阀门开度、PWM 占空比。两者的转换只需要一个累加：

$$ u(n) = u(n-1) + \Delta u(n) $$

工程上还有个细节：位置式递推里 $u(n-1)$ 这一项在浮点运算下会不断累积误差，长时间运行可能出现极缓慢的漂移。我的做法是在执行器侧加一个软限幅，并且每隔一段时间用一次绝对式重算（把累加器清一次、用当前误差重建积分量），能有效遏制这种漂移。这个技巧教科书上一般不讲，但在连续运行几个月都不停机的设备上很重要。

5. 写进代码时才会遇到的坑

推完公式只是完成了一半，真正让人头疼的都在代码落地环节。下面这几条，全是我自己踩出来的。

5.1 采样周期 T 的选取：不是越小越好

新手有个直觉：T 越小越接近连续，肯定越好。实际完全不是这么回事。

T 太小，会带来两个麻烦。第一，微分项 $\frac{K_d}{T}$ 系数被急剧放大，而实际采样值是被量化和噪声污染的，微分项会把噪声放大成剧烈的控制抖动，执行器跟着抖，机械件很快就磨损。第二，浮点运算的字长有限，T 很小时相邻两拍误差几乎相等，相减之后有效位数大量丢失，微分项算出来全是舍入噪声。

T 太大，则会丢掉高频动态，闭环带宽上不去，响应迟钝甚至振荡。

我的经验规则是：根据被控对象的主要时间常数 $\tau$ 来定，取 $T \approx \tau / 10$ 到 $\tau / 5$；或者根据期望闭环带宽 $f_c$，取 $T \le 1/(10 f_c)$。这两个口径下来，覆盖了绝大多数场合。如果对象是多时间尺度（比如有大惯性又有个小时间常数的执行器滞后），就按最小的那个时间常数来定，哪怕它很短。

提示：采样周期一旦定下来就必须在代码里硬编码一致。我见过最坑的一次是仿真脚本里 T=0.01，嵌入式代码里定时器配成了 0.02，两边系数没改，结果闭环增益凭空变了，找了一下午。

5.2 积分饱和与微分噪声：递推式的两个死穴

积分饱和的成因很直接：当输出已经到达执行器上限，误差还在持续累积，积分项就会一直往上加，等到误差反向时，控制量要花很长时间才能从饱和区退出来，表现为超调巨大、恢复缓慢。

解决思路是抗积分饱和：当输出被限幅时，停止积分累加，或者把积分项反向修正回来。工程上最省事的做法是"钳位法"——输出限幅之后，用限幅后的输出反算应该保留的积分量，把累加器按差值回退。这样不需要额外的状态判断，代码量小，效果好。

微分噪声的抑制有两招。第一招是加低通滤波，也叫"不完全微分"，把纯微分 $\frac{K_d s}{1}$ 换成 $\frac{K_d s}{1 + T_f s}$，$T_f = T_d / N$，N 一般取 5 到 10。离散化（仍用后向差分）之后：

$$ \text{微分项}(z) = \frac{K_d (1-z^{-1})}{T + T_f - T_f z^{-1}} $$

它不再是一个纯差分，而是一阶滤波后的差分，高频噪声被压下去，同时低频微分作用几乎不变。第二招是微分先行，把微分作用从误差上移到反馈量上，这样设定值突变时微分项不会产生冲击。

这两招我在温度控制和电机电流环上都用过，效果立竿见影。尤其是电流环，采样噪声大，不加不完全微分的话 PWM 输出会一直滋滋响。

5.3 浮点递推的累积误差与状态初值

浮点递推的误差累积问题，前面提过一次，这里展开说。像 $y(n) = a y(n-1) + (1-a)u(n-1)$ 这种形式，如果 $a$ 非常接近 1（大惯性对象就是），递推过程中的相对误差会被不断放大。单精度 float 只有约 7 位有效十进制数字，跑几万拍之后，末位误差就可能累积到肉眼可见。

解决办法有两个。一是能用双精度就用双精度，在工控和桌面级应用上，double 的开销完全可以接受。二是在算法层面做归一化，比如定期用绝对式重算，或者采用结构上更数值友好的实现形式（下面 5.4 会讲）。

状态初值的问题更隐蔽。差分方程里一般默认初始状态为零，但真实系统往往不是。比如一个温度控制器，启动时炉温已经是 200 度，如果你把历史状态全清零，控制器会以为当前温度是从零开始上升，积分项会往一个荒谬的方向跑。所以上电初始化时，一定要用实际测量值去填充历史状态，而不是简单地清零。这是从"能跑"到"跑得对"之间的一道坎。

5.4 实现结构：直接 I 型、直接 II 型与环形缓冲

同一个差分方程，有不同的程序实现结构。

直接 I 型就是最朴素的写法，每一个 $u(n-k)$ 和 $y(n-k)$ 都用一个独立变量存，代码读起来跟公式一一对应，调试时最容易对照。缺点是需要的历史变量最多，内存占用大。

直接 II 型（也叫典范型）是把公共延迟单元合并，用更少的存储单元实现同样的传递函数。它把中间变量统一成一条延迟链，非常适合内存紧张的嵌入式场景。代价是代码逻辑不那么直观，一旦系数符号出错很难定位。

实际操作上，当阶数很低（一阶、二阶，PID 也就是二阶），我一般就用直接 I 型，变量命名清楚点，可读性优先。当阶数上到四五阶（比如多阶滤波器），才考虑直接 II 型。

环形缓冲是另一个常用技巧。当差分方程里要用的历史项很多时，用一个固定长度的数组加一个写指针，每次采样把新值写进去、指针加一取模，就不用为每一拍单独建变量了。这样代码可以写得非常通用——一个函数处理任意长度的差分方程，系数从配置里读，不用为每个对象改代码。

/* 通用直接I型差分方程实现，环形缓冲版 */ #define MAX_ORDER 8 typedef struct { float b[MAX_ORDER + 1]; float a[MAX_ORDER + 1]; float u_hist[MAX_ORDER + 1]; float y_hist[MAX_ORDER + 1]; int order; } DiffEq; float diffeq_step(DiffEq *de, float u_in) { float y = 0.0f; int i; /* 移位历史输入 */ for (i = de->order; i > 0; i--) { de->u_hist[i] = de->u_hist[i - 1]; } de->u_hist[0] = u_in; /* 分子部分：b0*u(n) + b1*u(n-1) + ... */ for (i = 0; i <= de->order; i++) { y += de->b[i] * de->u_hist[i]; } /* 分母部分：-a1*y(n-1) - a2*y(n-2) - ... （a0 归一化为1） */ for (i = 1; i <= de->order; i++) { y -= de->a[i] * de->y_hist[i - 1]; } /* 移位历史输出 */ for (i = de->order; i > 0; i--) { de->y_hist[i] = de->y_hist[i - 1]; } de->y_hist[0] = y; return y; }

这段代码把"Z变换转差分方程"的成果直接变成了可复用的运行时。系数填进去就能跑，不用改代码。注意这里的a[0]假定已经是 1，前面归一化那一步没白做。

6. 一个可复现的验证流程：Python 先跑，C 再落地

公式推完、代码写完，最后一步是验证。我习惯的流程是先在 Python 里跑通，确认数值行为符合预期，再往嵌入式平台移植。这条流程能挡掉绝大部分低级错误。

6.1 用 Python 对照 ZOH 与差分递推的结果

拿第 3 节的一阶惯性环节做例子，把 ZOH 的解析结果和递推实现对照：

import numpy as np T = 0.1 # 采样周期 tau = 1.0 # 时间常数 N = 100 # ZOH 离散化系数 a = np.exp(-T / tau) b = 1.0 - a u = np.ones(N) # 单位阶跃输入 y = np.zeros(N) for n in range(1, N): y[n] = a * y[n-1] + b * u[n-1] # 对比：直接看几拍 print("前5拍输出:", y[:5]) print("稳态输出:", y[-1]) # 理论稳态值应为 1.0

跑出来稳态值应该收敛到 1.0，前几拍是 0、0.0952、0.1813、0.2592……如果第一步不是 0，说明输入延迟放错了位置——这里 u(n-1) 意味着第一拍输出是 0，因为 u(-1)=0。这个"第一拍为 0"的特征是区分 ZOH 和后向差分的最快判据，后者第一拍就是非零的。

再拿一阶惯性环节做个后向差分的对照，就能直观看到两者的差异：

y_back = np.zeros(N) c1 = tau / (tau + T) c2 = T / (tau + T) for n in range(N): y_back[n] = c1 * (y_back[n-1] if n > 0 else 0.0) + c2 * u[n] print("后向差分前5拍:", y_back[:5])

对比两组数就能确认自己没把两套系数搞混。

6.2 移植到 C 时的定点/浮点注意事项

从 Python 到 C，跨的不只是语言，还有数值表示。如果目标平台是成本敏感的单片机，用的是定点运算（Q15、Q31 那种），事情就复杂一层。

定点的核心问题是量程和精度。差分方程的系数范围可能跨好几个数量级，比如 $K_p$ 是几十，而 $\frac{K_d}{T}$ 可能是几百。如果统一用 Q15 表示，小数点位置固定，很容易溢出或者精度不够。工程上的常规做法是给每一项单独定标，或者干脆把整个递推式做一次缩放，把所有系数乘一个公共因子，让它落在定点能表示的范围里，最后输出时再除回来。

浮点平台就简单多了，但要注意编译器的浮点运算选项。有些平台默认用单精度，你得显式写 double 或者开 double 支持，不然精度不够。还有一点：浮点运算在中断里执行时间可能很长，如果控制周期很短，要评估一下最坏情况下的运算时间，必要时把部分计算放到主循环里做。

6.3 实测中我踩过的几个具体坑

最后说几个具体的教训，都是"课本上不会写、踩过才知道"的类型。

第一个是历史变量没做边界保护。刚开始写通用差分方程的时候，指针或者数组下标忘了取模，跑了几十万拍之后数组越界，现象是每隔一段时间控制器输出突然抽一下。这种错误在仿真里时间短看不出来，只有长期运行才暴露。后来我养成了习惯：访问历史数组的地方一律先做边界检查，或者用取模操作。

第二个是限幅位置放错。抗积分饱和的钳位逻辑，必须用"限幅之后的输出"去反算积分量。我早期写成了用限幅之前的输出反算，结果饱和判定永远不成立，抗饱和完全失效。这个坑很隐蔽，因为逻辑上看起来没错，但状态变量的更新时机错了。

第三个是从 Z 传递函数直接抄系数的符号错误。差分方程里 y 的历史项系数是取负号的，$y(n) = -\sum a_i y(n-i) + \cdots$，而 Z 传递函数分母里 $a_i$ 是正号。这两个符号关系一旦搞混，写进代码里就是全反，系统直接发散。我的对策是永远在代码注释里把原始的 Z 传递函数写上去，和系数一一对应，review 的时候能一眼看出符号对不对。

第四个是离散化方法和采样周期不匹配。前面 3.3 节讲过，T 大的时候三种离散化差异巨大。有一次我用 ZOH 算的系数，但实机上采样周期被临时改成了原来的一半（因为要提带宽），系数忘了重算，等效于把对象时间常数错认了一半，闭环勉强能跑但响应明显不对，后来重新算系数才恢复正常。采样周期和差分方程系数是死绑定的，改一个必须改另一个。

第五个是初次部署时不加软启动。位置式 PID 上电瞬间，积分项从零开始，但实际误差可能很大（比如设定值已经给了，反馈还没起来），控制量会瞬间冲到上限。加个简单的斜坡软启动，让设定值在几秒内从当前测量值爬升到目标值，能显著降低上电冲击，对各种执行器都友好。

这套流程走下来，从 Z 变换到差分方程再到现场代码，链条上每一环都可验证、可复现。我个人的体会是，公式推导只占这件事三成的难度，剩下七成都在采样周期选择、数值处理和实现细节上。你推得再漂亮，这几处处理不好，现场照样给你脸色看。