连续状态方程离散化方法详解：从欧拉到ZOH的工程实践

1. 连续状态方程与离散化的底层逻辑

1.1 先搞清楚连续状态方程在干什么

接触过控制、机器人或信号处理的朋友，对状态方程应该不陌生。它的标准写法是：

[ \dot{x}(t) = A x(t) + B u(t) ] [ y(t) = C x(t) + D u(t) ]

x是系统内部状态，u是输入，y是输出，A描述状态之间的耦合关系，B描述输入如何影响状态。这个模型用微分方程描述系统的连续动态，在数学上是干净漂亮的，但在工程落地时，它面临一个非常现实的问题：几乎所有数字系统——单片机、DSP、嵌入式Linux——都是按固定周期运行代码的。你要在中断里、在实时任务里执行控制律，拿到的是某一时刻的传感器采样值，输出的是离散的指令，它根本没法真正“连续”地积分微分方程。于是，连续状态方程离散化就成了必须跨过的一道坎。

这里要澄清一个常见的误区：离散化不是把微分换个符号写成差分那么简单。它是在“用离散的采样序列近似连续动态”这个前提下，重新构造一组状态转移关系。这组关系必须保证，在采样点上的行为与连续系统的真实行为尽可能一致，同时还得考虑稳定性、计算量、延迟等一系列工程因素。

1.2 离散化到底在解决什么核心矛盾

往深了说，离散化解决的是“连续系统”与“离散控制器”之间的接口问题。你在MATLAB里设计了一个极点配置控制器，连续域仿真跑得很漂亮，相位裕度、带宽都对，但一上板子，同样一组参数就发散或者抖得厉害，十有八九问题出在离散化这一步没有处理好。

本质上，连续系统的状态方程给出的是任意时刻t的状态演化规律，而数字控制器只能在这些离散的时刻t=kTs（Ts为采样周期）采样和输出。离散化要做的事情，是把连续系统在采样点上的行为压缩成一个差分方程：

[ x[k+1] = A_d x[k] + B_d u[k] ]

A_d、B_d就是离散状态矩阵，它们必须反映这样的信息：从第k个采样时刻的状态和输入出发，经过一个采样周期Ts后，状态在第k+1个采样时刻应该落在哪里。这个“转移”过程越贴近真实微分方程的解，离散模型就越可靠。

我经常用一个生活化类比来解释这件事：连续方程相当于你对着一块匀速运动的表连续读数，你随时都知道时间；离散方程相当于你只看整点时刻的表盘，然后靠规则推断每分钟之间的变化。采样频率越高，你漏掉的过程越少；但如果你只是简单地把“一分钟”当成“一小时”来近似，那结果会错得离谱。离散化方法的选择，本质上就是在回答“我们用什么规则，从离散读数反推连续过程”。

2. 四种主流离散化方法拆解

2.1 最直观的前向欧拉法：能用但容易翻车

前向欧拉法的思路非常朴素：用当前的导数近似下一个采样时刻的状态，即

[ x[k+1] = x[k] + T_s \dot{x}[k] = x[k] + T_s(A x[k] + B u[k]) ]

整理后得到：

[ A_d = I + T_s A, \quad B_d = T_s B ]

这个方法的好处是直观、计算量小，特别适合在单片机上简单快速验证。但它的致命问题在于稳定性条件非常苛刻。连续系统稳定的条件是A的特征值都在左半平面，而前向欧拉离散后，系统稳定的条件是特征值映射到离散域后落在单位圆内，即要求每个特征值λ都满足：

[ |1 + T_s \lambda| < 1 ]

对复特征值来说，这个条件意味着采样周期必须在系统动态时间常数的量级以下。如果Ts取得太大，哪怕原连续系统是稳定的，离散化之后也会振荡甚至发散。我第一次用前向欧拉做电机转速环的时候，就吃过这个亏——连续域的PI参数明明很稳，但采样周期拉大后转速一直在震荡，后来排查下来，正是离散化引入了不稳定极点。

所以我的建议是：前向欧拉适合系统动态较慢、Ts远小于最小时间常数的情况，适合做快速原型验证，但用在正式控制器或者仿真模型里，要特别谨慎。

2.2 后向欧拉法：稳定但会压低动态

后向欧拉法的形式是：

[ x[k+1] = x[k] + T_s \dot{x}[k+1] = x[k] + T_s(A x[k+1] + B u[k]) ]

因为x[k+1]同时出现在等式两边，需要解方程，整理得：

[ A_d = (I - T_s A)^{-1}, \quad B_d = (I - T_s A)^{-1} T_s B ]

这个方法的最大优势是：无论是连续系统本身是否振荡，只要Ts为正，离散系统通常都是稳定的，因为它相当于对极点了做了向单位圆内部收缩的映射。代价是什么呢？它会引入额外的相位滞后和阻尼，让系统看起来“钝化”了——响应变慢，震荡衰减得更快。

这就带来一个需要在实操里特别注意的问题：如果你用后向欧拉离散化后再去调控制参数，你会发现连续域分析时算好的带宽和相位裕度都对不上，系统的实际响应比设计预期迟钝。遇到这种情况，不要急着怀疑硬件，先检查是不是离散化方法带来的相移太大。后向欧拉更适合数值刚性系统——比如同时存在极快和极慢两种动态的系统——这时无条件稳定比相位精度更重要。

2.3 梯形法与双线性变换（Tustin变换）：频率域与时间域的桥梁

梯形法的思路是用区间两端斜率平均来近似积分，形式上相当于把s平面通过双线性变换映射到z平面：

[ s = \frac{2}{T_s} \cdot \frac{z-1}{z+1} ]

在状态空间里，它的离散化结果可以写成：

[ A_d = (I - \frac{T_s}{2}A)^{-1}(I + \frac{T_s}{2}A) ] [ B_d = (I - \frac{T_s}{2}A)^{-1} \cdot T_s B ]

这个方法的优势在于：如果连续系统稳定，离散系统必然稳定，稳定性保持特性比前向欧拉强得多；同时它在映射时保持了连续域与离散域之间的频率对应关系（只是发生了频率压缩畸变），所以很多从频域设计滤波器、补偿器的场景非常喜欢用它。

不过，双线性变换有一个著名的“频率畸变”问题。连续系统在角频率ω处的特性，会映射到离散域的某个频率ω_d，两者之间的关系是：

[ \omega_d = \frac{2}{T_s} \tan^{-1}\left(\frac{\omega T_s}{2}\right) ]

这导致在接近奈奎斯特频率的高频段，离散化后的频率响应会被明显压缩。如果你设计的陷波滤波器正好落在系统谐振点附近，直接做双线性变换后，陷波频率可能偏掉，这时就需要对该频率做预畸变（pre-warping），让变换前后这个关键频率能够精确对齐。

2.4 零阶保持器（ZOH）离散化：最接近物理实际的方法

ZOH是工业控制里最常使用的离散化方式，也最符合真实执行机构的物理行为。它的基本假设是：在每个采样周期内，输入u保持不变——这对DAC输出、PWM输出、阀门开度指令来说就是实际的工作方式，因为它们在两个采样时刻之间就是保持恒定的。

ZOH离散化直接从连续微分方程的解出发。对状态方程两边做积分：

[ x(t+T_s) = e^{A T_s} x(t) + \int_{t}^{t+T_s} e^{A(t+T_s-\tau)} B u(\tau) d\tau ]

由于u在周期内恒等于u[k]，可以得到经典结果：

[ A_d = e^{A T_s} ] [ B_d = \int_{0}^{T_s} e^{A \tau} d\tau \cdot B ]

如果A可逆，B_d还可以写成：

[ B_d = A^{-1}(e^{A T_s} - I) B ]

其中e^{AT_s}是矩阵指数。这个解是连续微分方程在采样点上的精确解，所以ZOH离散化在“输入保持假设成立”的前提下，误差仅来自采样周期本身，而没有任何数值积分近似误差。这是它相比于三种欧拉法和双线性变换的核心优势。

但ZOH也不是没有代价。它需要计算矩阵指数，计算量比其他方法大；另外在处理多输入多输出系统时，B_d的积分式要按照输入通道逐一处理。很多嵌入式工程师看到矩阵指数就头疼，其实现在库里基本都封装好了，MATLAB里用c2d(A,B,Ts,'zoh')，Python里用scipy.signal.cont2discrete，注意选对method参数就行。

我把四种方法的特性整理成下面这张表：

3. 采样时间与矩阵指数：离散化的核心细节

3.1 采样时间怎么选才靠谱

离散化所有的误差都跟Ts相关，所以采样时间的选择是整个工程里面最关键的决策之一。理论上，采样定理告诉我们采样频率至少是系统最高频率的两倍，但对控制工程来说，两倍远远不够。实际项目中我遵循的经验是：采样频率至少取系统闭环带宽的10到20倍，或者等效地，采样周期要小于系统最小时间常数的1/10到1/5。

举个例子，假设系统开环传递函数有一个极点位于s=-100，时间常数是10ms，那采样周期至少要小于2ms，否则前向欧拉必然发散，即便换了更稳定的方法，离散化后也会明显丢失高频动态。如果你设计一个带宽为10Hz的控制器，采样频率至少到100Hz才是起步，200Hz以上才算稳妥。

还有一个实际约束是传感器与执行器的物理周期。比如IMU的更新频率是500Hz，而你PID循环跑1000Hz，那你实际能用的控制周期是2ms而不是1ms。离散化模型必须和使用周期对齐，否则你计算出来的B_d就不对。这一点在做嵌入式系统时非常容易被忽视，很多人把离散化Ts设成1ms，但实际控制任务由于中断调度、任务切换，真实周期是1.3ms甚至抖动到2ms，这样模型和现实就出现了偏差。

3.2 矩阵指数到底应该怎么算

ZOH离散化绕不开e^{ATs}。矩阵指数的计算有几种典型思路：

• 幂级数展开：e^{ATs} = I + ATs + (ATs)^2/2! + …，适合Ts比较小、矩阵范数不大的情况。但级数截断误差很难精确控制，当矩阵有快速动态时，需要很多项才能收敛。
• 缩放与平方方法（scaling and squaring）：把e^{ATs}先转化为e^{ATs/m}，通过缩放让级数快速收敛，然后再把结果连续平方m次还原。这是科学计算库中的标配算法，可靠性和精度都很好。
• 特征值分解法：如果A可以对角化，A = VΛV^{-1}，那么e^{ATs} = V e^{ΛTs} V^{-1}。这个思路非常直观，但遇到亏损矩阵（不可对角化）时就不适用。
• 增广矩阵技巧：求B_d时，为避免A求逆，可以构造分块矩阵：

[ M = \begin{bmatrix} A & B \ 0 & 0 \end{bmatrix} ]

然后计算e^{MTs}，其右上分块恰好就是B_d。这个技巧在很多数值库中实现起来特别方便，我强烈推荐，因为它既避免了A^{-1}的数值稳定性问题，也不需要担心A奇异。

上述这些算法在scipy.linalg.expm、numpy的expm、MATLAB的expm里都已经高度优化。你自己实现时最忌讳的就是在Ts很大的情况只取泰勒级数前三项，那基本算一次错一次。我第一次用C语言在嵌入式端手写ZOH离散化时，直接级数展开了六项，结果高频系统算出来的矩阵根本不对，后来查资料发现SciPy里用的是缩放平方方法来控制误差，自己重写后才算对。

3.3 离散化矩阵算完后的校验手段

很多人算完A_d、B_d就往控制器代码里塞，结果运行不对，回头根本说不清是离散化算错了还是控制器逻辑错了。我建议把校验放到计算流程里，养成习惯。

第一步，检查A_d在u=0的情况下能否复现连续系统的零输入行为。取一个初始状态x0，连续系统在t=Ts时的理论解是e^{ATs}x0，而离散系统从x0经过一步得到A_d x0，两者应该一致。

第二步，检查直流增益是否一致。连续系统从u到y的直流增益是-CA^{-1}B + D（在A可逆时），离散系统是C(I-A_d)^{-1}B_d + D，它们应该一致。如果不一致，推导基本出了问题。

第三步，结合阶跃响应对比。给输入一个单位阶跃，看连续模型与离散模型在每个采样点上的输出是否重合，如果偏差很小，说明离散化可靠；如果偏差随采样步数累积，检查Ts是否过大或者方法是否选得不合适。

4. 实操过程：一个机械系统模型的离散化全流程

4.1 从连续模型到ZOH离散化的手算演示

这里用一个最简单的惯性系统来演示。假设系统方程为：

[ \dot{x} = -2x + u ]

采样周期取Ts=0.1s。因为A=-2是一个标量，矩阵指数直接就是e^{-2×0.1}=e^{-0.2}≈0.8187。B_d的计算式是：

[ B_d = \int_{0}^{0.1} e^{-2\tau} d\tau = \frac{1-e^{-0.2}}{2} \approx 0.0906 ]

于是离散状态方程是：

[ x[k+1] = 0.8187 x[k] + 0.0906 u[k] ]

如果改用前向欧拉，得到的是x[k+1]=(1-0.2)x[k]+0.1u[k]=0.8x[k]+0.1u[k]。两者差别不大，因为系统时间常数τ=0.5s，采样周期只有它的1/5。但如果把Ts拉大到0.5s呢？ZOH得到A_d=e^{-1}≈0.3679，B_d=(1-e^{-1})/2≈0.3161；而前向欧拉会得到A_d=1-1=0，系统仍然稳定但错得离谱；如果Ts再大一点，比如0.6s，前向欧拉的A_d=1-1.2=-0.2，绝对值小于1还好，一旦Ts超过1s，A_d绝对值就大于1，系统直接变不稳定。所以你看，对于一个连续稳定的一阶系统，ZOH依然时保持稳定，而前向欧拉在Ts超过某个阈值后就会彻底失真。

4.2 二阶系统用Python完整复现

工程中更常见的是二阶系统。以质量-弹簧-阻尼系统为例：

[ m\ddot{x} + c\dot{x} + kx = F ]

取m=1kg，c=0.5N·s/m，k=10N/m，状态变量x1=x，x2=\dot{x}，状态空间矩阵为：

[ A = \begin{bmatrix} 0 & 1 \ -10 & -0.5 \end{bmatrix}, \quad B = \begin{bmatrix} 0 \ 1 \end{bmatrix} ]

取Ts=0.05s，我们用Python来做ZOH离散化：

import numpy as np from scipy.linalg import expm from scipy.signal import cont2discrete # 连续系统矩阵 A = np.array([[0, 1], [-10, -0.5]]) B = np.array([[0], [1]]) C = np.array([[1, 0]]) D = np.array([[0]]) Ts = 0.05 # 方法一：直接计算矩阵指数 Ad = expm(A * Ts) # 构造增广矩阵求Bd M = np.zeros((3, 3)) M[:2, :2] = A M[:2, 2:] = B Md = expm(M * Ts) Bd = Md[:2, 2:] print("ZOH Ad:") print(Ad) print("ZOH Bd:") print(Bd) # 方法二：用scipy自带函数验证 system = (A, B, C, D) sysd = cont2discrete(system, Ts, method="zoh") print("SciPy Ad:") print(sysd[0]) print("SciPy Bd:") print(sysd[1])

两种方法输出的矩阵数值基本一致。实际运行中这个系统在ZOH离散化后的零输入响应与连续解在采样点上的偏差极小，说明离散化精度满足需求。

做完离散化之后，可以顺便看一下前向欧拉在这个例子中的表现。A_d(I+TsA)的特征值模长如果超过1，系统就会不稳定。用代码算一下特征值，你会发现当Ts取0.05s时勉强稳定但已经开始有误差；把Ts换成0.2s，特征值模长已经大于1，离散系统发散了。这个对比是最直观的说明：采样周期选大了之后，不是“精度变差”这么简单，而是“稳定与否”的本质区别。

4.3 双线性变换在陷波滤波器的实操

除了状态方程，双线性变换还经常用在各种数字滤波器设计里。假如你要设计一个陷波频率fn=50Hz，采样频率fs=1000Hz的数字陷波器，直接在连续域设计再双线性变换，陷波频率会偏低。这时需要使用预畸变：

[ f_{pre} = \frac{f_s}{\pi} \tan\left(\frac{\pi f_n}{f_s}\right) ]

用压缩后的频率去设计连续域陷波器，再双线性变换回来，最终的数字陷波点才精确落在50Hz。我当初在控制某型电机的转速波动时，机械谐振频率非常明显，靠的就是这个预畸变技巧才把滤波器正中谐振峰，整个方案的震动噪音降了好几个dB。

5. 常见问题与排查技巧实录

5.1 离散化后系统发散

这是最常见的坑。优先检查采样周期是否过大，尤其是前向欧拉方法。判断方法很简单：求解离散系统A_d的特征值，看是否都在单位圆内。不要只看连续系统稳定就觉得离散系统稳定，它们之间没有必然的推论关系。

另外检查B_d的计算是否涉及A^{-1}，如果A本身奇异，用积分式或者增广矩阵法，不要强行求逆。有些矩阵在低维看起来没问题，但高维时条件数很大，求逆后数值出现明显误差。

5.2 阶跃响应的稳态值对不上

连续系统稳态输出与离散系统稳态输出之间差了很多，这是典型的直流增益不匹配。ZOH方法只要计算正确，稳态值是天然匹配的；但欧拉法和双线性变换都会产生不同程度的偏差。一个很实用的修正方法是，在进行控制器离散化之前，计算连续系统的直流增益，然后对B_d乘以一个缩放系数，使离散系统的稳态增益对齐。不过我更推荐直接采用ZOH，这样就不需要额外修正了。

5.3 仿真步数与控制步数不一致

很多人在Simulink或者自写仿真里，把离散化的Ts与仿真器固定步长混为一谈。仿真步长可以比Ts小很多倍，但控制器内部计算只用Ts。如果你把仿真步长当成了离散化Ts，必然导致计算结果错误。这个问题我在给研究生答疑时碰到过好几次，症状是同样的代码换个求解器结果完全不同，实际上就是步长混用造成的。

5.4 模型与实测有偏差怎么办

如果ZOH离散化做完，仿真结果和实测还是对不上，先别急着怀疑离散化。检查输入通道是否真的有“保持”性质。比如PWM输出本身是在固定周期内保持恒定，这符合ZOH假设；但如果你通过其他方式实现模拟量输出，输出可能不是理想保持型，那么ZOH模型的假设就不成立了。另一个常见来源是传感器响应延迟、执行器饱和与死区。离散化模型通常是线性的，但真实系统在这些非线性环节上有明显特征，这部分需要单独建模，不是离散化方法能解决的。

5.5 离散化误差究竟该怎么评估

经验上，评估离散化误差可以这样操作：对系统输入一个宽频激励（比如扫频信号或阶跃），同时用连续模型和离散模型仿真，计算两者输出的均方根误差。误差随Ts的变化趋势通常呈O(Ts^2)级别（ZOH和双线性变换）或O(Ts)级别（欧拉法）。如果你想精确判断自己的模型处于什么精度水平，可以算一下误差比值。对比下来，ZOH在高动态、高精度场景中的优势非常明显，这也是为什么它成为工业界事实标准的原因。

6. 实际工程中的几个重大经验

我在实际项目里反复踩过离散化的坑，最后沉淀下来的几条经验可能对大家有帮助。

第一，能用ZOH的地方尽量用ZOH。不管是控制器设计、卡尔曼滤波器的预测步，还是动力学仿真，ZOH在“输入保持”这个实际前提下是最精确的。自己实现矩阵指数并没有想象中困难，调库也行，手写级数配合缩放平方也可以，关键是不要图省事用前向欧拉替代。

第二，恒用同一个Ts。建模用1ms离散化，控制器跑2ms周期，这种混搭是最容易出问题的。建议在项目初期就确定统一的控制周期，并把所有模型的离散化都基于这个统一周期完成。如果系统多速率运行，宁可做明确的多速率离散化设计，也不要稀里糊涂地“当作相同周期”处理。

第三，把离散化结果写进单元测试。每次修改模型参数之后，自动对比连续系统与离散系统在采样点上的阶跃响应误差，保证误差保持在可接受范围内。这个做法能让你在后续算法迭代中快速发现回归问题，省去很多调试时间。

第四，高频谐振类系统务必注意频率畸变。如果你在连续域设计滤波器再转离散，一定要确认关键频率的对应关系。对匹配有强要求时，使用预畸变或直接采用基于数字域频率设计的方法。

连续状态方程离散化这件事，看起来只是数学变换里的一小步，但它连接了理论设计与实际运行世界。把这一步搞扎实了，后面无论是做仿真、做控制器还是做状态估计，都会少走很多弯路。