第7章

用 Python 仿真控制系统

学完这一章你应该能

用 python-control 建立传递函数、组合框图、求极点和零点；
计算阶跃响应和性能指标、波特图、稳定裕度、根轨迹、奈奎斯特图；
自己写一个「离散对象 + 数字控制器」的仿真循环（app.py 就是这样工作的）；
用代码验证前面各章的手算结果，知道常见的坑。

前置知识：第2～6章的概念。Python 只需要会写函数、循环和用 numpy 数组。

为什么要会用代码？手算能让你理解原理，但实际工作中的系统往往是五阶、十阶，还带延迟和饱和，没人手算。工程师的日常是：先建模，在电脑上仿真、调参，确认没问题了再上真机。这一章的代码都在 python-control 0.10.2 上运行过，注释里给出的输出就是实际结果，你可以对照着跑。

7.1 准备环境

需要 python-control（控制理论工具箱，接口和 MATLAB 的 Control System Toolbox 很像）、numpy、matplotlib。三种方式任选一种：

不改动这个仓库的依赖，临时用：在项目目录执行 uv run --with control 你的脚本.py，uv 会临时装上 python-control 再运行；
把它加进项目依赖：uv add control（会修改 pyproject.toml 和 uv.lock）；
其他环境：pip install control matplotlib numpy scipy。

7.2 建立模型

import numpy as np
import control as ct

s = ct.tf('s')                           # 先定义 s，后面可以像写公式一样写传递函数
G = 1 / (s * (s + 1) * (s + 3))          # 1/(s(s+1)(s+3))
G_same = ct.tf([1], [1, 4, 3, 0])        # 同一个模型，用系数写（降幂排列）
L = 4 * G                                # 开环传递函数，K = 4
T = ct.feedback(L, 1)                    # 单位负反馈的闭环
print(T)                                 # 4 / (s^3 + 4 s^2 + 3 s + 4)
print(ct.poles(T))                       # [-3.4675, -0.2662+1.0405j, -0.2662-1.0405j]
print(ct.dcgain(T))                      # 1.0

常用操作：

操作	写法
串联、并联	`G1 * G2`、`G1 + G2`
负反馈（反馈支路为 H）	`ct.feedback(G, H)`；正反馈加 `sign=1`
零极点形式	`ct.zpk([-2], [0, -1, -5], 10)` 即 \(\frac{10(s+2)}{s(s+1)(s+5)}\)
极点、零点、直流增益	`ct.poles(G)`、`ct.zeros(G)`、`ct.dcgain(G)`
约去公共零极点	`ct.minreal(G)`
转成状态空间	`ct.ss(G)`（第11章）
离散化	`ct.c2d(G, Ts, method="zoh")`（第8章）

7.3 时域响应

import matplotlib.pyplot as plt

t = np.linspace(0, 30, 3001)
resp = ct.step_response(T, t)            # 阶跃响应
plt.plot(resp.time, resp.outputs)
plt.grid(); plt.show()

info = ct.step_info(T)                   # 性能指标，调节时间默认按 2% 误差带
print(info['Overshoot'], info['PeakTime'], info['SettlingTime'])
# 42.5  3.30  13.45   和第4章、第5章算出来的一致

# 任意输入：比如斜坡
r = t.copy()
resp = ct.forced_response(T, t, r)
print(r[-1] - resp.outputs[-1])          # 斜坡输入的稳态误差 ≈ 1/Kv

比较几个增益下的响应，就是第3章、第4章做的事：

for K in [1, 4, 8]:
    r = ct.step_response(ct.feedback(K * G, 1), t)
    plt.plot(r.time, r.outputs, label=f"K={K}")
plt.legend(); plt.grid(); plt.show()

7.4 频域：波特图、稳定裕度、奈奎斯特图

ct.bode_plot(L, dB=True, display_margins=True)   # 波特图，并标出两个裕度
plt.show()

gm, pm, wcg, wcp = ct.margin(L)
print(20 * np.log10(gm), pm, wcg, wcp)
# 9.54 (dB)  29.8 (°)  1.732 (相位穿越频率)  0.932 (截止频率)

ct.nyquist_plot(L)                        # 奈奎斯特图
plt.show()

mag, phase, omega = ct.frequency_response(L, [0.5, 1, 2])
print(mag, np.degrees(phase))
# [2.3527 0.8944 0.2481]  [-126.03 -153.43 172.87]

注意相位的折算：上面 ω=2 时真实相位是 −187.13°，但函数返回的是 172.87°（折算到了 ±180° 之间）。读数时要心里有数，需要连续的相位曲线时用 np.unwrap(phase)。另外 margin 返回的幅值裕度 gm 是倍数，不是 dB。

7.5 根轨迹

ct.root_locus_plot(G)                     # 画根轨迹，可以在图上点选看增益
plt.show()

for K in [1, 4, 12, 20]:                  # 数值版：直接算几个增益下的闭环极点
    print(K, np.round(ct.poles(ct.feedback(K * G, 1)), 3))
# 12 [-4.+0.j  -0.+1.732j  -0.-1.732j]   ← 临界增益，极点在虚轴上
# 20 [-4.363+0.j  0.182+2.133j  0.182-2.133j]   ← 进入右半平面

7.6 自己写仿真循环

python-control 的函数适合分析线性系统。可是实际控制器是离散的（每隔 \(T_s\) 算一次），还有限幅、抗饱和这些非线性，这时最直接的办法是自己写一个循环。仓库里的 app.py 和 plant.py 就是这样做的：把连续对象离散化成状态方程，每一步做一次矩阵乘法。

import numpy as np
import control as ct

# 被控对象 1/(s+1)^3，离散化成每 Ts 秒更新一次的状态方程（零阶保持）
Ts = 0.01
plant = ct.c2d(ct.ss(ct.tf([1], [1, 3, 3, 1])), Ts, method="zoh")
A, B, C, D = plant.A, plant.B, plant.C, plant.D


class PID:
    def __init__(self, kp, ki, kd, out_limit=None):
        self.kp, self.ki, self.kd = kp, ki, kd
        self.out_limit = out_limit
        self.integral = 0.0
        self.prev_measurement = None

    def update(self, setpoint, measurement, dt):
        error = setpoint - measurement
        if self.prev_measurement is None:
            d = 0.0
        else:
            d = -(measurement - self.prev_measurement) / dt
        self.prev_measurement = measurement
        u = self.kp * error + self.ki * self.integral + self.kd * d
        saturated = self.out_limit is not None and abs(u) >= self.out_limit
        if not saturated or error * u < 0:
            self.integral += error * dt
        return u


def simulate(pid, T=40.0, umax=10.0):
    n = int(T / Ts)
    x = np.zeros((A.shape[0], 1))
    t = np.arange(n) * Ts
    y = np.zeros(n)
    u = np.zeros(n)
    for k in range(n):
        y[k] = (C @ x).item()
        r = 1.0 if t[k] >= 1 else 0.0                          # t = 1 s 时设定值阶跃
        u[k] = np.clip(pid.update(r, y[k], Ts), -umax, umax)   # 执行器限幅
        x = A @ x + B * u[k]                                   # 对象前进一步
    return t, y, u


t, y, u = simulate(PID(4.8, 2.65, 2.18))
print(f"超调 {100 * (y.max() - 1):.1f}%")                     # 超调 53.0%
t, y, u = simulate(PID(4.8, 2.65, 2.18, out_limit=1.5), umax=1.5)
print(f"限幅 1.5 + 抗饱和：超调 {100 * (y.max() - 1):.1f}%")   # 2.0%

这个循环的结构（读测量值 → 算控制量 → 限幅 → 对象前进一步）和单片机里的控制循环一模一样，只是对象换成了真实的硬件。把这段代码改一改，就能仿真你自己的项目：换对象模型、加噪声、加延迟、加扰动都很容易。第6章和第9章的 PID 实验台就是用 JavaScript 写的同一个循环。

7.7 用代码验证前面的例子

下面每一条都可以自己写几行代码验证，括号里是应该得到的结果：

第3章按指标设计：\(G=\dfrac{3.63}{s(0.276s+1)}\) 的闭环阶跃响应（超调 16.3%，峰值时间 1.0 s）；
第3章测速反馈：\(\dfrac{4}{s^2+2.8s+4}\)（超调 4.6%）与 \(\dfrac{4}{s^2+1.2s+4}\)（超调 37.2%）对比；
第4章临界增益：\(\dfrac{K}{s(s+1)(s+5)}\) 在 \(K=30\) 时的闭环极点（\(\pm j2.24\)）；
第5章：\(\dfrac{10}{s(0.1s+1)(0.5s+1)}\) 的裕度（截止频率 4.08 rad/s，相位裕度约 3.9°）；
第6章超前校正：\(\dfrac{4(1.56s+1)}{s(s+1)(s+3)(0.267s+1)}\) 的相位裕度（约 50°）；
第6章滞后校正：\(\dfrac{6(19.4s+1)}{s(s+1)(s+3)(66s+1)}\) 的相位裕度（约 49°）。

常见的坑

系数要按降幂写，缺的项补 0：\(s^3+3s\) 写成 [1, 0, 3, 0]；
margin 的幅值裕度是倍数；相位会被折算到 ±180°；
仿真时间太短，step_info 算出的调节时间不对；时间步长太大，曲线不准；
连续系统和离散系统不能直接相乘或反馈（离散系统带 dt）；
自己写循环时，先读测量、再算控制、最后让对象前进，顺序搞错相当于多了或少了一拍延迟；
脚本里画完图要 plt.show()，否则窗口不出来。

本章小结

建模：ct.tf、s = ct.tf('s')、ct.zpk、ct.feedback；分析：poles、step_response、step_info、bode_plot、margin、root_locus_plot、nyquist_plot。
线性分析用库函数；有离散控制器、限幅、非线性时自己写仿真循环。
任何手算结果，都应该习惯用几行代码验证一下。

练习

1. 用 python-control 求 \(\dfrac{K}{s(s+2)(s+5)}\) 的临界增益（提示：对 K 循环求闭环极点，或者用 margin）。

G = ct.tf([1], [1, 7, 10, 0])
gm, pm, wcg, wcp = ct.margin(G)
print(gm, wcg)    # 70.0  3.162：K=1 时幅值裕度 70 倍，所以临界增益 70，振荡频率 √10

手算验证：特征方程 \(s^3+7s^2+10s+K=0\)，劳斯表 \(s^1\) 行 \(\frac{70-K}{7}\)，\(K_u=70\)；辅助方程 \(7s^2+70=0\)，\(\omega=\sqrt{10}\)。

2. 把 7.6 节的仿真改成一阶对象 \(\frac{1}{s+1}\)、只用 P 控制，验证稳态误差 \(\frac{1}{1+K_p}\)。

把 ct.tf([1], [1, 3, 3, 1]) 换成 ct.tf([1], [1, 1])，用 PID(4, 0, 0)，末尾输出约 0.8，误差 0.2 = 1/(1+4)。

3. 在 7.6 节的循环里加入测量噪声和一个负载扰动，应该改哪几行？

噪声：把 pid.update(r, y[k], Ts) 里的 y[k] 换成 y[k] + np.random.normal(0, 0.01)（控制器只能看到带噪声的测量）。扰动：把 x = A @ x + B * u[k] 换成 x = A @ x + B * (u[k] + d)，其中 d 在某个时刻之后取 −0.5 之类的值（负载扰动加在对象输入端）。