从数据建模：系统辨识

10.1 设计辨识实验

1. 从稳态出发：先让系统在正常工作点附近稳定下来（手动给一个固定输入），记录一段稳态数据，用来估计噪声水平。
2. 选激励信号：
   • 阶跃：最简单直观，够用来得到 FOPDT 模型和整定 PID；
   • 随机二值信号（PRBS）：在两个电平之间随机切换，包含丰富的频率成分，适合最小二乘辨识；
   • 正弦扫频：逐个频率测幅值和相位，直接得到波特图，适合找机械共振。
3. 选幅值：要足够大，让响应明显高于噪声（至少 5～10 倍）；又要足够小，保持在线性范围内、不触发保护、不影响安全。
4. 选采样频率：大约是系统带宽的 10 倍。不是越快越好：采得太快时，后面的最小二乘法会出问题（10.3 节）。
5. 同时记录输入和输出，并且时间对齐。最好多做几次，不同幅值、正反两个方向各做一次，用来检查线性和重复性。

10.2 阶跃响应法：得到 FOPDT 模型

大多数过程对象（温度、液位、压力、电机速度）的阶跃响应是「先延迟一段，再像一阶惯性环节那样上升」，可以用一阶加纯延迟模型近似：

两点法：输入阶跃 \(\Delta u\)，输出从一个稳态变到另一个稳态，变化量为 \(\Delta y\)。

1. 稳态增益 \(K=\Delta y/\Delta u\)；
2. 找出输出变化到 \(\Delta y\) 的 28.3% 和 63.2% 的时刻 \(t_{28}\)、\(t_{63}\)（从阶跃时刻算起）；
3. \(T=1.5(t_{63}-t_{28})\)，\(L=t_{63}-T\)。

为什么是这两个数：对 FOPDT 模型，\(1-e^{-1/3}=28.3\%\) 和 \(1-e^{-1}=63.2\%\) 分别出现在 \(t=L+T/3\) 和 \(t=L+T\)，两式相减就得到 \(T\)。

在组件里切换对象：

• 三阶对象 \(\frac{1}{(s+1)^3}\)：得到 \(T\approx2.1\) s、\(L\approx1.15\) s。模型和真实响应的最大偏差约 11%，看起来不算很准，但用它按 SIMC 整定的 PI（\(K_p=0.92\)，\(K_i=0.44\)），在真实对象上的超调只有 8%；
• 本来就是 FOPDT 的对象：得到的参数几乎完全正确；
• 有超调的二阶对象：两点法的结果不可信，要用二阶模型或下一节的方法；
• 加大噪声：阈值时刻会被噪声干扰，所以实际数据要先滤波，或者用多次测试的平均值。

import numpy as np

def fopdt_two_point(t, y, du):
    """t、y：阶跃之后的时间和输出（已减去阶跃前的稳态值）；du：输入阶跃的幅值"""
    dy = np.mean(y[-50:])                 # 最后一段的平均值当作新的稳态
    K = dy / du
    t28 = t[np.argmax(y >= 0.283 * dy)]
    t63 = t[np.argmax(y >= 0.632 * dy)]
    T = 1.5 * (t63 - t28)
    L = t63 - T
    return K, T, L

# 对 2/(s+1)^3 做 0.5 的阶跃测试、带测量噪声，得到 K=2.001  T=2.160  L=1.110

另一种常见做法是切线法：在响应曲线斜率最大的点作切线，切线与初始值、稳态值两条水平线的交点分别给出 \(L\) 和 \(L+T\)。手工在纸上量很方便，但对噪声更敏感。

10.3 最小二乘辨识

阶跃法只能得到很简单的模型。更一般的做法是假设一个离散模型（第8章的差分方程），用数据把系数「解」出来。最常用的是 ARX 模型：

对每个 \(k\) 写一个这样的方程，把未知系数排成向量 \(\theta=[a_1,a_2,b_1,b_2]^T\)，所有方程合起来就是 \(Y=\Phi\theta\)。方程个数远多于未知数，取使误差平方和最小的解（最小二乘解）：

import numpy as np
import control as ct

rng = np.random.default_rng(1)
G_true = ct.tf([4], [1, 1.2, 4])           # 「未知」对象，只用来产生实验数据
Ts = 0.2


def experiment(n):
    u = np.repeat(np.sign(rng.standard_normal(n // 5 + 1)), 5)[:n]    # 随机 ±1，每 5 拍变一次
    t = np.arange(n) * Ts
    y = ct.forced_response(ct.c2d(G_true, Ts, "zoh"), t, u).outputs
    return t, u, y + 0.01 * rng.standard_normal(n)                   # 加测量噪声


t, u, y = experiment(2000)
# 模型 y(k) = -a1·y(k-1) - a2·y(k-2) + b1·u(k-1) + b2·u(k-2)，写成 Y = Φθ
Phi = np.column_stack([-y[1:-1], -y[:-2], u[1:-1], u[:-2]])
a1, a2, b1, b2 = np.linalg.lstsq(Phi, y[2:], rcond=None)[0]
model = ct.tf([b1, b2], [1, a1, a2], dt=Ts)

z = np.roots([1, a1, a2]).astype(complex)  # 离散极点换算回连续域：s = ln(z)/Ts
s = np.log(z) / Ts
print(abs(s[0]), -s[0].real / abs(s[0]), (b1 + b2) / (1 + a1 + a2))
# ωn=2.003  ξ=0.303  直流增益=1.002   （真值 2、0.3、1）

# 用另一组没参与拟合的数据验证（见 10.4 节）
t2, u2, y2 = experiment(1000)
y_hat = ct.forced_response(model, t2, u2).outputs
print(100 * (1 - np.linalg.norm(y2 - y_hat) / np.linalg.norm(y2 - y2.mean())))   # 拟合度 98.9%

输出误差法：直接假设一个连续模型（比如 \(\dfrac{K\omega_n^2}{s^2+2\xi\omega_ns+\omega_n^2}\)），用它仿真出输出，调整参数让仿真输出和实测输出的差最小。计算量大一些，但对噪声和采样周期不敏感，参数也有物理意义：

from scipy.optimize import least_squares

Ts = 0.05                                  # 这次故意采得很快（接着上面的代码运行）
t, u, y = experiment(2000)

def simulate(p, t, u):
    K, wn, xi = p
    G = ct.c2d(ct.tf([K * wn**2], [1, 2 * xi * wn, wn**2]), Ts, "zoh")
    return ct.forced_response(G, t, u).outputs

res = least_squares(lambda p: simulate(p, t, u) - y, x0=[1.0, 1.0, 0.7],
                    bounds=([0, 0.01, 0.01], [10, 50, 5]))
print(res.x)
# K=1.000  ωn=2.000  ξ=0.300，在新数据上的验证拟合度 98.4%

初值 x0 不能离真值太远，可以先用阶跃响应的特征（超调量、峰值时间）粗估一个。

10.4 模型验证

用来拟合的数据上效果好，不代表模型就好（参数多了什么数据都能凑上）。一定要用另一组没参与拟合的数据检验：

• 拟合度：\(\text{fit}=100\left(1-\dfrac{\|y-\hat y\|}{\|y-\bar y\|}\right)\%\)，\(\hat y\) 是模型对同一输入的仿真输出。做控制用，一般 80% 以上就不错了；
• 看残差：\(y-\hat y\) 应该像随机噪声，如果还能看出明显的规律（比如和输入同步的波形），说明模型漏掉了某些动态；
• 换工作点、换幅值再测：如果在不同条件下参数差别很大，说明系统有明显的非线性，可能要分段建模（第12章的增益调度）。

def fit_percent(y, y_hat):
    """拟合度：100% 表示完全吻合，0% 相当于只用平均值去猜"""
    return 100 * (1 - np.linalg.norm(y - y_hat) / np.linalg.norm(y - y.mean()))

10.3 节的例子：ARX（\(T_s=0.2\) s）在新数据上的拟合度为 98.9%，采样太快时只有 58%，输出误差法为 98% 以上。

10.5 频率响应测量

机械系统（机械臂、云台、长传动轴）常有共振，阶跃测试看不清楚。这时用正弦扫频：输入一个频率逐渐变化的正弦信号，在每个频率上测输出和输入的幅值比、相位差，直接画出波特图。共振峰的位置和高度一目了然，可以据此加陷波滤波器、确定控制带宽的上限（一般要比第一个共振频率低很多）。很多伺服驱动器和飞控软件都自带扫频功能。

10.6 实际注意事项

• 先看原始数据：画出来看有没有饱和、死区、明显的非线性、传感器跳变，再决定用什么模型；
• 去掉稳态值：辨识用的是相对工作点的增量（第2章的线性化），先减去实验前的稳态值；
• 纯延迟要单独考虑：延迟会让低阶模型很难拟合，先从数据里估出延迟，再把输入平移；
• 闭环下的数据要小心：系统在反馈控制下运行时，输入和噪声相关，直接用最小二乘会有偏差。尽量在开环（手动）下做实验；
• 模型够用就行：做控制设计，FOPDT 或二阶模型往往就够了。模型的不准确由反馈和留足的稳定裕度来兜底。

本章小结

练习