离散化：位置式、增量式与采样率

采样、零阶保持和那个半拍#

数字控制器每隔 $T_s$ 读一次传感器，算出 $u_k$，然后在下一个采样点之前让它保持不变：零阶保持（ZOH）。保持器的传递函数是

代入 $s = j\omega$，把分子拆成 $e^{-j\omega T_s/2}$ 乘一个正弦：

相位是精确的 $-\omega T_s/2$，幅值是 $T_s$ 乘一个 sinc。 所以保持器在相位上就是一个半拍的纯延迟，没有近似。幅值那一项在回路的穿越频率 $\omega_c \approx 46$ rad/s、$T_s = 1$ ms 处是 $\text{sinc}(0.023) \approx 0.9999$，可以忽略。我的仿真里命令从算出来到进入电机还要多过一拍（DelayLine(1)，第 3 篇的 TRANSPORT），所以回路的有效延迟是

这是我要检验的第一个预测：用无延迟回路的 $\omega_c$ 和 $1.5T_s$ 估计相位裕度的损失，再和精确的采样模型对一对。

$z = e^{sT_s}$ 和三种近似#

连续系统的极点 $s = p$，采样后变成离散系统的极点 $z = e^{pT_s}$。这是精确的对应，但控制器里的 $1/s$（积分）和 $s$（微分）要用 $z$ 的有理函数实现，而 $s = \ln z / T_s$ 不是有理函数。所以要近似。三种常见的近似，对应 $e^{sT_s}$ 的三种展开：

• 前向欧拉： $z = e^{sT_s} \approx 1 + sT_s$，即 $s \approx \dfrac{z - 1}{T_s}$。
• 后向欧拉： $z^{-1} = e^{-sT_s} \approx 1 - sT_s$，即 $s \approx \dfrac{z - 1}{T_s z}$。
• Tustin（双线性）： $e^{sT_s} = \dfrac{e^{sT_s/2}}{e^{-sT_s/2}} \approx \dfrac{1 + sT_s/2}{1 - sT_s/2}$，即 $s \approx \dfrac{2}{T_s}\dfrac{z - 1}{z + 1}$。

用它们离散积分器 $I(s) = k_i/s$，得到三条递推：

我在第 1 篇写的 PidBasic（先 integral += e * ts 再输出）就是后向那一行。区别只在于 $I_k$ 里有没有包含当前这一拍的误差。

三种近似在频率响应上的差别。 取 $z = e^{j\theta}$，$\theta = \omega T_s$，把三个积分器的 $1/(j\omega)$ 提出来：

前向多出半拍滞后，后向多出半拍超前，Tustin 的相位是精确的 $-90°$，只有幅值有一点偏差。这半拍和 ZOH 的半拍符号相反：后向的积分器刚好抵消一部分保持器的滞后，前向则叠加上去。我猜这就是 PidBasic 那种写法在 1 kHz 下“怎么都能用”的一半原因，另一半是 1 ms 本来就很小。

稳定区域。 $z = e^{sT_s}$ 把左半平面映到单位圆内。三种近似各自的“映到单位圆内的 $s$ 区域”是：

• 前向：$|1 + sT_s| < 1$，是以 $-1/T_s$ 为圆心、$1/T_s$ 为半径的圆盘。一个稳定的连续极点落在圆外，离散化之后会变成不稳定。
• 后向：$|1 - sT_s| > 1$，是以 $+1/T_s$ 为圆心的圆的外部。一个不稳定的连续极点落在圆内，离散化之后反而稳定。
• Tustin：恰好是左半平面，稳定性一一对应。代价是频率轴被压缩，模拟频率 $\omega_a$ 和数字频率 $\Omega$ 的关系是 $\omega_a = \dfrac{2}{T_s}\tan\dfrac{\Omega T_s}{2}$，模拟的无穷远被折到 Nyquist 频率。

左图是虚轴 $s = j\omega$ 在三种映射下的像：前向是 $\text{Re}\,z = 1$ 这条竖线，永远在单位圆外；后向是以 0.5 为圆心的小圆，完全在单位圆内；Tustin 精确地落在单位圆上。右图是反过来看，哪些 $s$ 被映进了单位圆。最后这个区别后面会直接用到：Tustin 把 $s \to -\infty$ 映到 $z = -1$。

滤波微分的离散化#

纯微分 $K_d s$ 在高频增益无穷大，实际一定要接一个一阶滤波，$D(s) = K_d s/(t_f s + 1)$。分别代入后向和 Tustin。

后向： $s = \dfrac{z - 1}{T_s z}$，

分子分母同除以 $z$，交叉相乘，得到递推

Tustin： $s = \dfrac{2}{T_s}\dfrac{z - 1}{z + 1}$，分子分母同乘 $T_s(z + 1)$：

两条递推的形状一样，$d_k = a\,d_{k-1} + b\,(e_k - e_{k-1})$，区别只在极点 $a$ 和增益 $b$。连续滤波器的极点在 $s = -1/t_f$，精确的离散极点是 $e^{-T_s/t_f}$；两种近似给出的是

取 $t_f = 0.1\,T_s$（滤波器比采样还快，我选它是因为这是“只想去掉一点点毛刺”时最容易写出来的设置），精确值是 $e^{-10} = 4.5\times 10^{-5}$，后向给 $0.091$，Tustin 给 $-0.667$。后向的极点是正的，误差落在“还是一个平滑的衰减”那一侧；Tustin 的极点是负的，每一拍符号翻转。 前向欧拉在这里是 $z = 1 - T_s/t_f = -9$，在单位圆外，所以我的实现直接拒绝 FORWARD 作为微分方法。

Tustin 为什么会这样：它把 $s \to -\infty$（滤波器的高频极限）映到 $z = -1$，正好是 Nyquist 点。采样把快速旋转的车轮拍成倒转，这是同一类现象（电影里的 wagon-wheel 错觉）：映射把一个“快得没法表示”的模式折叠到 Nyquist 频率。在 $z = -1$ 处两个滤波器的增益是

比值是 $(2t_f + T_s)/(2t_f) = 6$。这是第二个预测：同样的噪声，Tustin 的 D 项应该在 Nyquist 附近放大得更多，并且符号交替。

位置式和增量式#

位置式保存积分，每拍直接算出输出：

其中 $\Delta I_k$ 是上一节表里右边的增量（前向、后向或 Tustin）。增量式保存的是上一拍的输出，每拍只算变化量：

两者在没有限幅时等价：把位置式写成 $u_k - u_{k-1}$，$k_p e_k$ 的差分是 $k_p(e_k - e_{k-1})$，$I_k - I_{k-1} = \Delta I_k$，$d_k - d_{k-1}$ 原样保留，就是增量式的右边。起点也要一致：我让 $e_{-1} = 0$、$d_{-1} = 0$、$u_{-1} = 0$，所以没有“第一拍不算微分”的特殊处理（第 5 篇的最终 PID 才会加上）。

限幅时它们就不等价了，原因是各自的状态不同。 位置式有两个量：积分 $I$ 和被限幅的输出 $u$。限幅只动 $u$，$I$ 照样累加，账本和钱包就对不上了：执行器只能给 5 A，账上却记着 1300 A。增量式只有一个量，上一拍的输出 $u_{k-1}$，限幅作用在它身上，状态本身就被夹住了。

但这个好处依赖一个条件：限幅要加在累加器上。 如果增量式只把“要发给电机的值”限幅，累加器仍然是无限制地累加，那它在数学上就是位置式，一样卷绕。我会跑这个算例来验证（sat_inc_raw）。

限幅后的增量式到底是什么？ 令 $I_k := u_k - k_p e_k - d_k$（把输出里不属于积分的部分减掉）。不限幅时代入增量式，$I_k = I_{k-1} + \Delta I_k$，就是位置式的积分。限幅时 $u_k$ 被夹在边界上，$I_k$ 被重新定义成“让输出刚好等于边界的那个值”。设位置式的未限幅输出是 $v_k = k_p e_k + I_{k-1} + \Delta I_k + d_k$，那么

这正是反算抗饱和（back-calculation），跟踪时间常数取 $T_t = T_s$，也就是一个采样周期就把积分拉回去。第 5 篇会把 $T_t$ 做成可调参数；这里我先在 report 里验证这个等价（下面“限幅”一节）。

ts 为什么必须显式出现#

连续的 $k_i$ 的量纲是 $k_p$ 除以时间，$k_d$ 是 $k_p$ 乘时间。离散化之后每拍加的是 $k_i T_s e_k$，微分是 $k_d (e_k - e_{k-1})/T_s$。据我所知，不少 RM 开源代码把 $T_s$ 直接折进增益里：ki_code = ki * ts、kd_code = kd / ts，循环里只写 integral += ki_code * e。这在固定采样率下完全正确，但这样做之后，$T_s$ 变了代码里的增益不会跟着变。比如 $k_i = 90$，从 1 ms 改到 5 ms，如果 ki_code 还是 $90 \times 0.001$，等效积分增益就是 $0.09/0.005 = 18$，只剩 1/5。我没有逐个翻那些仓库，这一条是从常见写法推的，不是审计结论。我自己的做法是把 ts 放进配置、让 update 每拍用 ts 乘，这样改采样率只改一个数。

Python：一个 update，三个开关#

模块是 pidlab/pid_discrete.py。我先把上面的三个选择做成配置里的三个独立字段，而不是十二个类：

1
class Form(Enum):
2
    POSITIONAL = 0
3
    INCREMENTAL = 1
4

5

6
class Method(Enum):
7
    FORWARD = 0
8
    BACKWARD = 1
9
    TUSTIN = 2
10

11

12
@dataclass
13
class DiscreteConfig:
14
    # ...
15
    kp: float
16
    ki: float
17
    kd: float
18
    ts: float
19
    tf: float = 0.0
20
    u_min: float = -math.inf
21
    u_max: float = math.inf
22
    form: Form = Form.POSITIONAL
23
    i_method: Method = Method.BACKWARD
24
    d_method: Method = Method.BACKWARD

三种积分、两种微分、两种形式，一共 12 个组合，每个组合是 update 里的一个分支。构造函数拒绝五种配置：ts <= 0、tf < 0、u_min > u_max，还有两条针对微分的：

1
        if cfg.d_method is Method.FORWARD:
2
            raise ValueError("derivative has no forward method")
3
        if cfg.d_method is Method.TUSTIN and cfg.tf <= 0.0:
4
            raise ValueError("Tustin derivative needs tf > 0")

前者是上面算出来的 $z = -9$，后者是 $t_f = 0$ 时 Tustin 的极点落在 $-1$，每拍翻符号永不衰减。

核心是 update。先算微分的递推，再算积分增量 di，最后按形式拼出输出：

1
    def update(self, r: float, y: float) -> float:
2
        c = self.cfg
3
        e = r - y
4
        de = e - self.e_prev
5
        d_old = self.d_state
6
        if c.d_method is Method.BACKWARD:
7
            a = c.tf / (c.tf + c.ts)
8
            self.d_state = a * d_old + c.kd / (c.tf + c.ts) * de
9
        else:
10
            a = (2.0 * c.tf - c.ts) / (2.0 * c.tf + c.ts)
11
            self.d_state = a * d_old + 2.0 * c.kd / (2.0 * c.tf + c.ts) * de
12
        if c.i_method is Method.FORWARD:
13
            di = c.ki * c.ts * self.e_prev
14
        elif c.i_method is Method.BACKWARD:
15
            di = c.ki * c.ts * e
16
        else:
17
            di = 0.5 * c.ki * c.ts * (e + self.e_prev)
18
        if c.form is Form.POSITIONAL:
19
            self.integral += di
20
            v = c.kp * e + self.integral + self.d_state
21
        else:
22
            v = self.u_prev + c.kp * de + di + (self.d_state - d_old)
23
        u = clamp(v, c.u_min, c.u_max)
24
        self.u_prev = u
25
        self.e_prev = e
26
        return u

这段和数学一节的公式是一一对应的：d_state 是 $d_k$，di 是前面表里右边的增量，两个分支就是位置式和增量式。注意两处：

• u_min、u_max 的含义随形式而变。 位置式里它们夹的是输出 v，integral 不受影响；增量式里 u_prev 存的就是夹过的输出，所以夹的是累加器。self.u_prev = u 在位置式里没有用到，写上它是为了让两种形式共用一个出口。
• 没有“第一拍不算微分”。 第 1 篇的 PidBasic 有 primed 标志，第一拍的微分为 0。这里 e_prev 从 0 开始，第一拍看到的是一个从 0 跳到 $e_0$ 的阶跃，微分会踢一脚（后面 Tustin 的算例里这一脚是 11.1 A）。我这么做是为了让位置式和增量式从第 0 拍起逐位一致；第 5 篇的 PID 会处理微分冲击。

仿真循环和精确的采样模型#

驱动 python/ch04_discretization.py 的仿真循环和第 3 篇的差别只有两处：控制器换成 PidDiscrete，以及一个只在 report 里用的 back_calc 开关，用来验证“限幅增量式等价于反算”：

1
    for k in range(len(data["t"])):
2
        r = data["r"][k]
3
        y = plant.output(x) + data["n"][k]
4
        u = pid.update(r, y)
5
        ua = min(max(u, -case.act_limit), case.act_limit)
6
        if back_calc:
7
            pid.integral += ua - u
8
        rows.append([data["t"][k], r, y, u, ua, pid.d_state])
9
        x = advance(
10
            plant.deriv, x, [delay.push(ua), data["d"][k]], ts, SUBSTEPS
11
        )

u 是控制器的输出，ua 是夹过执行器限幅后真正发给电机的命令，两列都写进 CSV。back_calc 那一行就是教科书的反算：把被夹掉的 ua - u 加回积分。

要回答“相位裕度到底掉了多少”，只用连续模型加延迟是近似。sampled_pm 把电机（连同电流环）用 scipy.linalg.expm 做精确的 ZOH 离散化，控制器用 $C(z) = k_p + k_i T_s z/(z-1)$（后向积分），再乘一拍延迟 $z^{-1}$：

1
    em = np.zeros((3, 3))
2
    em[:2, :2] = a
3
    em[:2, 2] = bv
4
    e = expm(em * ts)
5
    ad, bd = e[:2, :2], e[:2, 2]
6
    w = np.linspace(1.0, math.pi / ts * 0.999, 20000)
7
    z = np.exp(1j * w * ts)
8
    g = np.empty_like(z)
9
    for i, zi in enumerate(z):
10
        g[i] = np.linalg.solve(zi * np.eye(2) - ad, bd)[0]
11
    lz = (kp + ki * ts * z / (z - 1.0)) * g / z**TRANSPORT

expm 作用在增广矩阵 $\begin{bmatrix}A & B\\0 & 0\end{bmatrix}T_s$ 上，左上角是 $A_d = e^{AT_s}$，右上角是 $B_d$，这是 ZOH 离散化的标准技巧。然后在 $[1,\ 0.999\pi/T_s]$ 上扫频率，在 $|L|$ 穿过 1 的那个区间里线性插值出 $\omega_c$ 和相位裕度。另一个函数 sampled_radius 把同一个回路写成四维状态（角速度、电流、积分器、在途的命令）的离散状态矩阵，取特征值的最大模，大于 1 就是不稳定。

这两个函数是分析用的，不是被检验的对象：它们只用 scipy 的 expm 和 numpy 的特征值，C、C++、Rust 不需要重写。report 里的对比表、抗饱和等价性验证、D 项统计和 4 张图都在同一个驱动里。

C、C++、Rust：逐行翻译，差别在边角#

三个版本是 Python 版的逐行翻译，update 里的算术顺序也和 Python 一致，所以后面才能逐位对上。差别在这几处：

• 枚举和配置。 C 是 enum pl_form、enum pl_method 加 struct pl_pid_discrete_config，没有 typedef；C++ 是 enum class Form、enum class Method 和模板结构体 PidDiscreteConfig<T>；Rust 是 enum Form、enum Method 和 PidDiscreteConfig，另外给了一个 PidDiscreteConfig::new(kp, ki, kd, ts) 填上默认值（不限幅、后向积分、后向微分），算例里只改需要改的字段。
• 构造的失败方式。 C 的 pl_pid_discrete_init 返回 0 或 -EINVAL；C++ 的 PidDiscrete<T>::Create 返回 std::optional，失败是 std::nullopt，没有异常；Rust 的 PidDiscrete::new 返回 Result<Self, PidDiscreteError>。
• 无穷。 不限幅的 u_min、u_max 在 C 里是 -INFINITY 和 INFINITY（<math.h>），C++ 是 std::numeric_limits<T>::infinity()，Rust 是 Real::INFINITY；Python 是 math.inf。clamp 对无穷边界没有特殊处理，x < -inf 永远为假。
• 内存。 三种语言里控制器都是栈上的值：C 的 struct pl_pid_discrete 把配置拷贝进来，C++ 和 Rust 的对象也是值语义，没有堆分配。

1
  if (c->i_method == PL_METHOD_FORWARD)
2
    di = c->ki * c->ts * pid->e_prev;
3
  else if (c->i_method == PL_METHOD_BACKWARD)
4
    di = c->ki * c->ts * e;
5
  else
6
    di = 0.5 * c->ki * c->ts * (e + pid->e_prev);
7
  if (c->form == PL_FORM_POSITIONAL) {
8
    pid->integral += di;
9
    v = c->kp * e + pid->integral + pid->d_state;
10
  } else {
11
    v = pid->u_prev + c->kp * de + di + (pid->d_state - d_old);
12
  }
13
  u = pl_clamp(v, c->u_min, c->u_max);

1
  static std::optional<PidDiscrete> Create(const PidDiscreteConfig<T>& cfg) {
2
    if (cfg.ts <= T(0) || cfg.tf < T(0) || cfg.u_min > cfg.u_max) {
3
      return std::nullopt;
4
    }
5
    if (cfg.d_method == Method::kForward) return std::nullopt;
6
    if (cfg.d_method == Method::kTustin && cfg.tf <= T(0)) {
7
      return std::nullopt;
8
    }
9
    return PidDiscrete(cfg);
10
  }

1
      const T a = (T(2) * cfg_.tf - cfg_.ts) / (T(2) * cfg_.tf + cfg_.ts);
2
      d_state_ = a * d_old + T(2) * cfg_.kd / (T(2) * cfg_.tf + cfg_.ts) * de;

1
        let di = match c.i_method {
2
            Method::Forward => c.ki * c.ts * self.e_prev,
3
            Method::Backward => c.ki * c.ts * e,
4
            Method::Tustin => 0.5 * c.ki * c.ts * (e + self.e_prev),
5
        };
6
        let v = if c.form == Form::Positional {
7
            self.integral += di;
8
            c.kp * e + self.integral + self.d_state
9
        } else {
10
            self.u_prev + c.kp * de + di + (self.d_state - d_old)
11
        };
12
        let u = v.clamp(c.u_min, c.u_max);

采样率：预测对了一半#

连续回路（电流环加 $\omega_n = 30$ 的 PI，没有延迟）的穿越频率是 45.98 rad/s，相位裕度 62.62°。表里“est loss”列是前面的预测 $\omega_c \cdot 1.5T_s$，“true loss”列是 62.62° 减去精确采样模型（sampled_pm）的相位裕度，“PM cont”是用第 3 篇的 freq 模块给连续模型加 $1.5T_s$ 延迟算的：

表里 $\omega_c$ exact 的单位是 rad/s，$\rho$ 是精确采样闭环的谱半径。

1、5、10 ms 上，预测在 1.2° 以内。 估计值（3.95°、19.76°、39.52°）和真实损失（3.60°、18.59°、38.61°）每一档都略大一点，差 0.35°、1.17°、0.91°。我猜这个同号的偏差来自上面推的后向欧拉积分器的半拍超前，加上 sinc 幅值的微小下降，它们都没进 $\omega_c \cdot 1.5T_s$ 的估计；“PM exact”总比“PM cont”高 0.35 到 1.17°，方向对得上，但我没有逐项把这两个修正量分离出来，所以这只是一个符号上吻合的解释。

20 ms 上预测失效。 估计损失 79.04°，真实损失 84.34°，差了 5.3°，更重要的是相位裕度已经是 $-21.72$°，闭环特征值的谱半径 $\rho = 1.1555$，每拍增长 15.5%。估计公式用的是无延迟回路的 $\omega_c = 45.98$ rad/s，而采样模型里真实的穿越点随 $T_s$ 在往上走：46.33、47.66、49.45，到 20 ms 是 54.09 rad/s，延迟拉长又推高穿越频率，滞后就比估计的多。这一步我只是从表上看出的趋势，没有单独验证。

超调和相位裕度同向变化：23.4%、36.4%、70.8%。 最后一行的 150.42% 不是超调：电流命令碰到 ±20 A 限幅（30 拍超过 20 A，控制器输出最大 42.10 A），响应变成一个幅值固定的极限环，图里的峰值在 12 rad/s 左右。（别被那个百分比骗了，它只是 $y$ 峰值除以 5 再减 1。）

1 ms 那一行也不是完全线性的。 阶跃一来比例项就是 $k_p \times 5 = 20.83$ A，超过 20 A 限幅，最大控制器输出 21.87 A，有 7 拍被夹住（表里最后一列）。这对 23.41% 的超调贡献很小，但表面上看起来“纯线性”的数字里混进了一点限幅，我写出来是为了你复现的时候不要奇怪。

位置式、增量式：不限幅时相等#

form_pos 和 form_inc 是同一个 PID（$k_p = 4.17$、$k_i = 90$、$k_d = 0.02$、$t_f = 5$ ms），1 ms，设定值 2 rad/s，0.5 s 加 0.3 N·m 负载，控制器输出最大 15.18 A，没有触及任何限幅。输出的最大差：

1
form_pos vs form_inc (no saturation): max|dy| = 4.44e-16, max|du| = 4.00e-15, max|u| = 15.18 A

差值是浮点舍入的量级。这是两条递推在数学上相等的最直接的证据，也是阈值 $10^{-12}$ 的来源：我没有给它留余量，因为真实差值比阈值小三个数量级。

限幅：增量式只有在累加器限幅时才抗饱和#

三个算例，同一个 PI（1 ms），执行器限幅 ±5 A，设定值 0.1 s 起跳到 30 rad/s，1.0 s 回到 0，一共跑 8 s：

• sat_pos：位置式，输出夹在 ±5 A，积分不限。
• sat_inc：增量式，累加器夹在 ±5 A。
• sat_inc_raw：增量式，累加器不限（u_min = -inf、u_max = inf），只有执行器（仿真循环里的 ua）限幅。

1
        case  peak y  overshoot   max|u|  settle after drop
2
     sat_pos   38.77      29.2%      5.0             5.067s
3
     sat_inc   30.12       0.4%      5.0             0.532s
4
 sat_inc_raw   38.77      29.2%   1328.4             5.067s

上图是轮速，中图是送到电机的电流，下图是增量式累加器的值（对数刻度，$\pm 5$ A 以内是线性的）。sat_pos 和 sat_inc_raw 的轮速曲线重合：两者的最大差是 $1.146\times 10^{-13}$。 它们的区别只是名字：把增量式的累加器放开，它就变成了位置式。蓝线和绿色虚线在图上是同一条线，我在图例里特意画了一条虚线才看得出有两条。

对位置式，轮速冲到 38.77 rad/s（超调 29.2%），因为从 0.1 s 到 1.0 s 电流一直被夹在 5 A，积分器却在一直累加：sat_inc_raw 里的累加器（等于位置式未限幅的和）绝对值最大到 1328.4 A，是执行器 5 A 限幅的 266 倍。设定值回到 0 之后，积分要先花时间“倒空”，期间执行器一直饱和在反方向，轮速反向冲到约 $-33$ rad/s，再来回摆，直到设定值回零之后 5.067 s 才稳定在 $\pm 1$ rad/s 以内。

对限幅累加器的增量式，轮速到 30.12 rad/s（超调 0.4%）就停住了，下降后 0.532 s 就进入 $\pm 1$ rad/s：超调 29.2% 对 0.4%，回稳时间 5.067 s 对 0.532 s，差 9.5 倍。

限幅累加器就是反算抗饱和，$T_t = T_s$。 我在数学一节推了这个等价，这里验证：驱动里 run_case(..., back_calc=True) 让位置式 PID 每拍执行 pid.integral += ua - u（把夹掉的部分加回积分），这是教科书的反算，$T_t = T_s$。它和 sat_inc 的轮速最大差：

1
max|y_inc - y_positional_backcalc(tt = ts)| = 3.553e-15 (clamped accumulator is back-calculation)

浮点舍入的量级。所以“增量式不会卷绕”这句话，更准确的说法是：限幅的增量式是跟踪时间常数最短的反算。 这也有代价：$T_t$ 取一拍，意味着控制器一退出饱和，积分立刻被重置到“刚好让输出等于边界”，比 $T_t$ 取几个采样周期更激进，我猜在有测量噪声的回路里会让输出在边界附近抖。我没有测这一点，第 5 篇的 BACK_CALC 会把 $T_t$ 做成参数。

Tustin 微分的振铃#

d_backward 和 d_tustin：DialPosition（拨弹盘）上的 PD，$k_p = 15$、$k_d = 0.667$，$t_f = 0.1T_s = 0.1$ ms，测量噪声 $\sigma = 10^{-3}$ rad（numpy 种子 4，来自数据文件），设定值 0.01 rad。0.3 s 之后的统计：

1
       case    pole    rms u    rms d  lag-1 of d  hf share
2
 d_backward  +0.091    0.8597    0.8482   -0.4790    0.8079
3
   d_tustin  -0.667    2.8303    2.8239   -0.8261    0.9647
4

5
       case              D step response, samples 0..4  |D| at Nyquist  max|u|  rms y err
6
 d_backward    606.1     55.1      5.0      0.5      0.0          1111.1    6.29   1.04e-03
7
   d_tustin   1111.1   -740.7    493.8   -329.2    219.5          6666.7   11.10   1.04e-03

左图是 D 项对单位误差阶跃的响应（对数刻度，保留符号）：后向是 606.1、55.1、5.0、0.5，每拍乘 0.091，两拍之内就衰减了两个数量级；Tustin 是 1111.1、$-740.7$、493.8、$-329.2$、219.5，每拍乘 $-0.667$，符号翻转，衰减缓慢，第 4 拍（从 0 数起）还有 219.5。

第二行表里的 $|D|$ at Nyquist 是 1111.1 和 6666.7，比值 6.0：这就是前面推导的 $(2t_f + T_s)/(2t_f)$。同一份噪声下，Tustin 的 D 项 rms 是后向的 $2.8239/0.8482 = 3.33$ 倍，lag-1 自相关从 $-0.479$ 到 $-0.826$（纯白噪声的差分是 $-0.5$，所以后向这一项基本就是“噪声的差分”；Tustin 是它加上负极点带来的翻转），功率落在上半频段（高于 $f_s/4$）的比例从 80.8% 涨到 96.5%。中图和右图是同样的结论，频谱上 Tustin 的曲线整体高一截，尤其在 400 Hz 以上。

设定值阶跃的那一脚也不一样： 最大的控制器输出出现在阶跃那一拍，后向是 6.29 A，Tustin 是 11.10 A，对应表里的 max|u|。

但我没能从输出里看到它对轮子的影响。 两行最后一列的 rms 跟踪误差都是 $1.04\times 10^{-3}$ rad，和噪声 $\sigma = 10^{-3}$ 几乎一样，因为记录的 $y$ 是含测量噪声的读数，不是真实角度。这个指标被噪声底限住了，区分不出两者。D 项振铃的代价是在 $u$ 上：电流指令的 rms 高 3.3 倍，在每拍 $\pm$ 几安培的翻转里；M3508 的电流环带宽比 1 kHz 低得多，所以我猜电机本身不会跟着抖，但发热和 CAN 上的指令噪声是真实的代价。这是推测，我没有在真电机上测。

实验清单#

这一篇新增的是 pid_discrete 模块和 ch04 驱动；仿真循环、植物和 CSV 比对沿用前几篇。链接指向 tag tutorial-v1。