尧图网站建设 
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 个采样点对应的真实时间是 nTT 是采样周期。所以严格写应该是 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_{i1}^{n} a_i y(n-i) \sum_{j0}^{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|_{z0.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 0g(1) -1.6667 2.6667 1g(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}{z1}$ 替换。两者都是近似方法都不是采样点精确但在工程上极其常用因为它们的变形幅度可控、稳定性好。先看后向差分。代入 $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} $$分子分母除以 $(\tauT)z$归一化$$ G(z) \frac{\frac{T}{\tauT}}{1 - \frac{\tau}{\tauT} z^{-1}} $$差分方程$$ y(n) \frac{\tau}{\tauT} y(n-1) \frac{T}{\tauT} u(n) $$注意这里跟 ZOH 有个关键差别ZOH 用的是 $u(n-1)$后向差分用的是 $u(n)$也就是当前拍输入直接进来。这会让后向差分的响应快一拍瞬时不那么平滑但好处是不需要额外记输入历史。再看双线性变换。代入 $s \frac{2}{T}\cdot\frac{z-1}{z1}$记 $k 2\tau/T$$$ G(z) \frac{1}{\tau \cdot \frac{2}{T}\frac{z-1}{z1} 1} \frac{z1}{k(z-1) (z1)} \frac{z1}{(k1)z (1-k)} $$归一化$$ G(z) \frac{1}{1k} \cdot \frac{1 z^{-1}}{1 \frac{1-k}{1k} z^{-1}} $$差分方程$$ y(n) -\frac{1-k}{1k} y(n-1) \frac{1}{1k}\left[u(n) u(n-1)\right] $$化简系数$-\frac{1-k}{1k} \frac{k-1}{k1} \frac{2\tau - T}{2\tau T}$$\frac{1}{1k} \frac{T}{2\tauT}$。所以$$ 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$三种方法的系数如下离散化方法y(n-1) 系数u(n-1) 系数u(n) 系数ZOH0.9048370.0951630后向差分0.90909100.090909双线性0.9047620.0476190.047619T 小的时候三者接近误差在千分之几的量级随便挑哪个都不影响项目。但把采样周期放大到 T 1 再看离散化方法y(n-1) 系数u(n-1) 系数u(n) 系数阶跃第一步输出ZOH0.3678790.63212100.6321后向差分0.500.50.5双线性0.3333330.3333330.3333330.6667差距一下就出来了。同一个连续对象采样周期同样是 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_{j0}^{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_0A_1A_2)e$而 $A_0A_1A_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)$。这两个口径下来覆盖了绝大多数场合。如果对象是多时间尺度比如有大惯性又有个小时间常数的执行器滞后就按最小的那个时间常数来定哪怕它很短。提示采样周期一旦定下来就必须在代码里硬编码一致。我见过最坑的一次是仿真脚本里 T0.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 变换到差分方程再到现场代码链条上每一环都可验证、可复现。我个人的体会是公式推导只占这件事三成的难度剩下七成都在采样周期选择、数值处理和实现细节上。你推得再漂亮这几处处理不好现场照样给你脸色看。 从Z变换到差分方程：位置式PID离散化与代码实现
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 个采样点对应的真实时间是 nTT 是采样周期。所以严格写应该是 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_{i1}^{n} a_i y(n-i) \sum_{j0}^{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|_{z0.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 0g(1) -1.6667 2.6667 1g(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}{z1}$ 替换。两者都是近似方法都不是采样点精确但在工程上极其常用因为它们的变形幅度可控、稳定性好。先看后向差分。代入 $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} $$分子分母除以 $(\tauT)z$归一化$$ G(z) \frac{\frac{T}{\tauT}}{1 - \frac{\tau}{\tauT} z^{-1}} $$差分方程$$ y(n) \frac{\tau}{\tauT} y(n-1) \frac{T}{\tauT} u(n) $$注意这里跟 ZOH 有个关键差别ZOH 用的是 $u(n-1)$后向差分用的是 $u(n)$也就是当前拍输入直接进来。这会让后向差分的响应快一拍瞬时不那么平滑但好处是不需要额外记输入历史。再看双线性变换。代入 $s \frac{2}{T}\cdot\frac{z-1}{z1}$记 $k 2\tau/T$$$ G(z) \frac{1}{\tau \cdot \frac{2}{T}\frac{z-1}{z1} 1} \frac{z1}{k(z-1) (z1)} \frac{z1}{(k1)z (1-k)} $$归一化$$ G(z) \frac{1}{1k} \cdot \frac{1 z^{-1}}{1 \frac{1-k}{1k} z^{-1}} $$差分方程$$ y(n) -\frac{1-k}{1k} y(n-1) \frac{1}{1k}\left[u(n) u(n-1)\right] $$化简系数$-\frac{1-k}{1k} \frac{k-1}{k1} \frac{2\tau - T}{2\tau T}$$\frac{1}{1k} \frac{T}{2\tauT}$。所以$$ 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$三种方法的系数如下离散化方法y(n-1) 系数u(n-1) 系数u(n) 系数ZOH0.9048370.0951630后向差分0.90909100.090909双线性0.9047620.0476190.047619T 小的时候三者接近误差在千分之几的量级随便挑哪个都不影响项目。但把采样周期放大到 T 1 再看离散化方法y(n-1) 系数u(n-1) 系数u(n) 系数阶跃第一步输出ZOH0.3678790.63212100.6321后向差分0.500.50.5双线性0.3333330.3333330.3333330.6667差距一下就出来了。同一个连续对象采样周期同样是 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_{j0}^{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_0A_1A_2)e$而 $A_0A_1A_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)$。这两个口径下来覆盖了绝大多数场合。如果对象是多时间尺度比如有大惯性又有个小时间常数的执行器滞后就按最小的那个时间常数来定哪怕它很短。提示采样周期一旦定下来就必须在代码里硬编码一致。我见过最坑的一次是仿真脚本里 T0.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 变换到差分方程再到现场代码链条上每一环都可验证、可复现。我个人的体会是公式推导只占这件事三成的难度剩下七成都在采样周期选择、数值处理和实现细节上。你推得再漂亮这几处处理不好现场照样给你脸色看。