简介这份系统辨识教案面向自动化专业学生与相关课程学习者围绕从输入输出数据中建立系统模型这一核心问题梳理了辨识理论、实验设计与工程应用之间的完整脉络。资源为单个doc文档压缩包约1.42MB内容以课程讲义形式组织涵盖参考书目、学习要求、考核办法与课堂讲授要点便于按教学进度自学或备课。文档从实体与模型的基本概念切入介绍直觉模型、物理模型与数学模型的分类并展开Zadeh关于数据、模型类和准则三要素的辨识定义进而结合长网造纸过程实例逐步讲解明确辨识目的、收集先验知识、设计试验、数据采集与预处理、模型类选择、结构辨识与参数估计、模型检验及转换等一般步骤。已有148人学习适合希望减少课堂依赖、独立掌握辨识实验设计与编程计算并借助最小二乘等方法完成建模分析的读者参考。1. 系统辨识教案从“调参玄学”到可复现的建模流程如果你在自动化专业待过大概率听过一句话“模型调得好不好全靠手感。”尤其是在现场调试时面对一个温度对象、一个电机转速环或者一个多容液位系统很多人第一反应不是去翻教材而是打开仿真软件凭经验拖几个极点、改几个增益看着响应曲线“差不多”就收工。这种做法在项目紧、工期短的时候确实能糊弄过去但一旦对象特性变了、工况漂移了或者换个人接手整套参数就成了黑匣子谁也不敢动。系统辨识教案要解决的恰恰是这个痛点。它不是让你背会最小二乘法的推导而是给你一套从实验设计、数据采集、模型结构选择到参数估计和验证的完整流程。自动化专业的学生学完能交出一份像样的课程设计一线工程师看完能直接在现场跑一遍阶跃响应把传递函数估出来而不是继续靠试凑。这篇笔记就按“理论先立住、再动手能复现”的路子把系统辨识从教案变成你工具箱里真能用的东西。2. 系统辨识的数学底子模型结构、准则函数与可辨识条件2.1 为什么先定模型结构再谈参数估计很多人一上来就急着写最小二乘代码结果数据跑出来拟合度0.99但模型一放到闭环里就发散。问题往往不在算法而在模型结构选错了。系统辨识的本质是在一个候选模型集合里找一个跟观测数据最贴合的成员。这个集合就是模型结构比如ARX、ARMAX、OE、BJ或者连续域的传递函数形式。结构选错参数估得再准也是白搭。常见做法是先根据物理机理判断对象是自平衡还是非自平衡是过阻尼还是欠阻尼有没有纯延迟。比如一个加热炉温度对象通常用一阶惯性加纯延迟来近似一个双容水箱至少得用二阶惯性。如果对象存在明显振荡二阶欠阻尼模型比一阶更合适。这一步不需要精确但方向不能错。提示模型结构不是越复杂越好。阶次每增加一阶待估参数就多一组对数据信噪比的要求也更高。现场数据往往只有几千个点高阶模型容易过拟合。2.2 最小二乘的准则函数与可辨识性条件参数估计的核心是定义一个准则函数衡量模型输出和实际输出之间的差距。最常用的是输出误差平方和J(θ) Σ (y(k) - ŷ(k|θ))²其中y(k)是实际观测输出ŷ(k|θ)是模型在参数θ下的预测输出。最小二乘就是让J(θ)最小的θ。对于ARX模型这个优化有解析解因为预测输出是参数的线性函数。但对于OE、ARMAX这些模型预测输出非线性依赖参数就得用迭代优化比如高斯-牛顿法或Levenberg-Marquardt。可辨识性条件经常被忽略。简单说输入信号必须充分激励系统的所有模态。如果你用一个恒定输入去辨识二阶系统第二个模态根本激发不出来参数估计就是病态的。常见做法是用伪随机二进制序列PRBS或者多正弦信号作为激励。PRBS实现简单在PLC里用几个移位寄存器就能生成对工业现场很友好。import numpy as np def generate_prbs(n, amplitude1.0, seed42): 生成伪随机二进制序列 n: 序列长度 amplitude: 信号幅值 seed: 随机种子保证可复现 rng np.random.default_rng(seed) # 每个时钟周期随机取 amplitude 或 -amplitude prbs rng.choice([amplitude, -amplitude], sizen) return prbs # 生成1000点PRBS幅值0.5 u generate_prbs(1000, amplitude0.5) print(u[:20])这段代码生成的是最简形式的PRBS每个采样点独立随机取正负幅值。实际工程中更常用的是基于移位寄存器的PRBS频谱更接近白噪声。参数amplitude要根据对象允许的输入范围来定太小信噪比不够太大可能激发非线性。一般取对象正常工作范围的10%到20%。2.3 连续模型与离散模型的转换关系自动化专业教材里经常混着讲连续传递函数和离散差分方程初学者容易晕。这里理一下现场采集的数据天然是离散的所以参数估计通常在离散域做。但控制器设计、频域分析又习惯用连续模型。两者之间的桥梁是零阶保持器ZOH下的离散化。假设连续传递函数G(s) K/(Ts1)在采样周期Ts下用ZOH离散得到G(z) K(1 - exp(-Ts/T)) / (z - exp(-Ts/T))反过来如果你先估出了离散模型想转回连续域可以用d2c函数但要注意采样周期不能太大否则高频信息丢失转换误差会很大。经验规则是采样频率至少是对象带宽的10倍对于一阶系统Ts取时间常数的1/10到1/5比较稳妥。import control as ct # 连续传递函数 G(s) 2 / (5s 1) G_s ct.tf([2], [5, 1]) # 采样周期 0.5 秒零阶保持器离散化 Ts 0.5 G_z ct.sample_system(G_s, Ts, methodzoh) print(连续模型, G_s) print(离散模型, G_z)control库的sample_system函数直接实现了ZOH离散化。method参数还可以选bilinear做双线性变换但ZOH更符合实际采样保持电路的物理行为。离散化后分子分母系数就是ARX模型可以直接用的参数。3. 从阶跃响应到传递函数手把手跑通一次辨识实验3.1 实验设计阶跃幅值、采样周期与数据长度怎么定阶跃响应法是现场最实用的辨识手段不需要额外设备PLC里给一个输出阶跃就行。但三个参数必须提前想清楚阶跃幅值、采样周期、数据长度。阶跃幅值一般取对象正常输入范围的5%到15%。太小了输出变化被噪声淹没太大了可能触发非线性或者安全限幅。我一般会先给一个试探性小阶跃看输出有没有明显响应再决定正式实验的幅值。采样周期按对象时间常数的1/10到1/20来选。如果你连时间常数大概多少都不知道可以先快速采一段看输出从开始变化到稳态用了多久那个时间的1/10就是粗略的采样周期。数据长度要覆盖输出从初始稳态到新稳态的整个过程通常取3到5倍时间常数。如果存在纯延迟还要额外加上延迟时间。import numpy as np import matplotlib.pyplot as plt # 模拟一个一阶惯性加纯延迟对象 # G(s) 1.5 / (8s 1) * exp(-2s) K 1.5 T 8.0 L 2.0 Ts 0.5 # 采样周期 t np.arange(0, 60, Ts) u np.zeros_like(t) u[t 5] 0.3 # 第5秒施加幅值0.3的阶跃 # 离散一阶惯性近似 y np.zeros_like(t) alpha np.exp(-Ts / T) for k in range(1, len(t)): # 纯延迟用整数个采样周期近似 delay_steps int(L / Ts) if k - delay_steps 0: y[k] alpha * y[k-1] (1 - alpha) * K * u[k - delay_steps] else: y[k] alpha * y[k-1] # 加一点测量噪声 y_meas y np.random.normal(0, 0.01, sizelen(t)) plt.plot(t, u, label输入 u) plt.plot(t, y_meas, label输出 y) plt.legend() plt.xlabel(时间 (s)) plt.grid(True) plt.show()这段代码模拟了一个典型的一阶加延迟对象。实际现场中你拿到的是u和y_meas两列数据模型参数K、T、L都是未知的。接下来的任务就是从这两列数据里把它们估出来。3.2 两点法、切线法与最小二乘拟合的实操对比教材里常讲两点法和切线法优点是手算快缺点是精度依赖选点。两点法是在阶跃响应曲线上取输出达到稳态值28.3%和63.2%的两个时刻分别对应t1和t2然后T 1.5(t2 - t1) L t2 - T这个方法对一阶系统很准但对噪声敏感。切线法是在拐点处画切线跟时间轴和稳态线的交点确定T和L手工画图误差更大。我一般直接用最小二乘拟合把连续模型离散化后做非线性优化。下面是一个完整的例子from scipy.optimize import minimize import numpy as np # 假设已有实验数据 t, u, y_meas # 这里用上面模拟的数据 def simulate_first_order(params, t, u, Ts): K, T, L params if T 0 or L 0: return np.inf * np.ones_like(t) alpha np.exp(-Ts / T) delay_steps max(0, int(L / Ts)) y_sim np.zeros_like(t) for k in range(1, len(t)): if k - delay_steps 0: y_sim[k] alpha * y_sim[k-1] (1 - alpha) * K * u[k - delay_steps] else: y_sim[k] alpha * y_sim[k-1] return y_sim def loss(params): y_sim simulate_first_order(params, t, u, Ts) return np.sum((y_meas - y_sim)**2) # 初始猜测 x0 [1.0, 5.0, 1.0] res minimize(loss, x0, methodNelder-Mead) print(估计参数 K, T, L , res.x)这里用Nelder-Mead单纯形法不需要计算梯度对初始值不敏感。参数说明K是稳态增益T是时间常数L是纯延迟。初始猜测可以随便给但T不能给0。拟合完成后一定要看残差是否白噪声如果残差还有规律说明模型结构不对。3.3 用Python的control库做ARX与OE模型估计如果数据量比较大或者对象阶次较高推荐用control库的辨识函数。它封装了ARX、OE等模型的估计流程底层是scipy的优化器。import control as ct import numpy as np # 假设 u, y_meas, Ts 已准备好 data ct.iddata(u, y_meas, Ts) # ARX模型阶次 [na, nb, nk] # na: 输出滞后阶数nb: 输入滞后阶数nk: 纯延迟 arx_model ct.arx(data, [2, 2, 1]) # OE模型阶次 [nb, nf, nk] oe_model ct.oe(data, [2, 2, 1]) print(ARX模型, arx_model) print(OE模型, oe_model) # 比较拟合度 ct.compare(data, arx_model, oe_model)arx函数的第二个参数是阶次列表[2,2,1]表示输出用前2个历史值输入用前2个历史值纯延迟1个采样周期。oe模型的参数含义类似但输出误差结构对噪声的处理更好。compare函数会画出两个模型的模拟输出和实际输出的对比拟合度用百分比表示一般到80%以上就可用。注意ARX模型假设噪声直接进入差分方程如果噪声是有色噪声ARX估计会有偏。OE模型假设噪声只影响输出测量对有色噪声更鲁棒。现场数据如果噪声明显优先试OE。4. 避坑与排查系统辨识现场翻车的五个典型场景4.1 数据里藏着非线性线性模型怎么估都不对现象拟合度死活上不去残差图呈现明显的规律性弯曲或者在不同工作点重复实验估出来的参数差异很大。原因对象存在死区、饱和、滞环等非线性特性线性模型只能在一个工作点附近近似。如果阶跃幅值跨过了非线性区数据本身就违背了线性假设。解决先做多组不同幅值的阶跃实验看稳态增益是否随幅值变化。如果变化明显要么缩小工作范围要么改用分段线性模型或Hammerstein-Wiener模型。现场最实用的做法是限制阶跃幅值保证对象工作在线性区。4.2 采样周期选太大纯延迟估计直接崩现象估计出的纯延迟L接近0或者等于采样周期但实际对象明明有延迟。原因采样周期Ts大于延迟时间L时延迟在离散数据里只体现为不到一个采样点的移位算法无法分辨。比如L0.3秒Ts0.5秒延迟信息就丢了。解决采样周期至少取延迟时间的1/5到1/10。如果PLC扫描周期限制做不到可以用过采样或者插值先把数据加密。另一个办法是先用高频采集卡单独测延迟再在模型里固定L只估K和T。4.3 输入信号激励不足参数估计病态现象最小二乘的协方差矩阵接近奇异参数估计值大得离谱或者换一组数据结果就完全变了。原因输入信号没有充分激发对象的所有模态。比如用阶跃信号辨识二阶系统第二个模态的激励很弱参数之间强相关。解决改用PRBS或者多正弦信号。PRBS的频谱在低频段能量集中适合工业对象。如果只能用阶跃至少做正反两个方向的阶跃增加信息量。另外数据长度要足够一般不少于1000个采样点。4.4 把闭环数据当开环数据用模型完全不可信现象现场设备不允许开环实验只能在闭环条件下采数据估出来的模型跟实际响应对不上。原因闭环数据里输入和输出通过控制器耦合输入信号不是独立激励最小二乘估计有偏。这是系统辨识里最经典的坑之一。解决如果必须用闭环数据需要知道控制器的结构用闭环辨识方法比如两步法或者联合输入输出辨识。更简单的做法是在控制器输出上叠加一个外部激励信号比如PRBS然后采集数据。这样输入信号里有了独立激励成分开环辨识方法就能用。4.5 模型验证只看拟合度不看残差和交叉验证现象训练数据拟合度95%换一组新数据掉到60%模型泛化能力差。原因过拟合。模型阶次选高了把噪声也拟合进去了。或者数据里存在趋势项没有去均值。解决把数据分成训练集和验证集用训练集估参数验证集算拟合度。两者差距大就是过拟合。另外画残差的自相关函数如果残差不是白噪声说明模型结构还有改进空间。我一般会要求验证集拟合度不低于训练集的80%否则就降阶或者重新选结构。5. 进阶技巧把辨识模型直接变成可用的控制器参数5.1 从辨识模型到PID参数的映射表辨识的最终目的通常是给控制器整定提供依据。如果你已经估出了一阶加延迟模型G(s) K exp(-Ls) / (Ts 1)可以直接用IMC整定公式算PID参数。下面这张表是我常用的映射关系Kp、Ti、Td分别对应比例增益、积分时间、微分时间。控制器类型KpTiTdPT / (K * L)——PIT / (K * L)T—PIDT / (K * L)TL / 2这张表假设目标是让闭环响应没有超调实际整定时可以根据需要微调。比如想要更快响应把Kp乘1.2但超调会增大。注意纯延迟L不能太小如果L接近0Kp会算得很大实际中要加限幅。5.2 用辨识模型做前馈补偿的实操步骤对于延迟大、扰动频繁的对象单靠反馈控制效果有限。用辨识出的模型做前馈可以显著改善跟踪性能。步骤很简单把辨识模型写成差分方程形式然后取逆作为前馈控制器。假设辨识得到离散模型y(k) a1 y(k-1) b1 u(k-1-d)其中d是延迟步数。前馈控制器的目标是在扰动可测时提前给出补偿输入。如果扰动w(k)到输出y(k)的通道也辨识出来了前馈增益就是两个通道稳态增益的比值。现场实施时前馈输出要加限幅并且只在扰动变化超过阈值时投入避免频繁动作。# 前馈补偿示例 # 假设扰动通道稳态增益 Kw控制通道稳态增益 Ku Kw 0.8 Ku 1.5 Kf -Kw / Ku # 前馈增益 def feedforward(w_current, w_last, Kf, threshold0.01): 简单前馈扰动变化超过阈值时输出补偿量 delta_w w_current - w_last if abs(delta_w) threshold: return Kf * delta_w return 0.0这段代码里Kf是前馈增益threshold是死区防止噪声引起频繁动作。实际投运时先把前馈增益设小一点观察效果再逐步加大。5.3 在线递推辨识让模型跟着工况慢慢变对象特性会随工况漂移比如换热器结垢后时间常数变大。离线辨识的模型用几个月就不准了。在线递推最小二乘RLS可以让模型参数缓慢更新适应这种变化。RLS的核心是每来一个新数据点就更新一次参数估计而不是重新算一遍。关键参数是遗忘因子λ通常取0.95到0.99。λ越小对新数据越敏感但参数波动也越大。我一般取0.98配合参数变化限幅防止单次异常数据把模型带偏。import numpy as np class RecursiveLeastSquares: def __init__(self, n_params, lam0.98, delta1e3): self.theta np.zeros(n_params) self.P np.eye(n_params) * delta self.lam lam def update(self, phi, y): phi: 回归向量 y: 当前输出 # 计算增益 P_phi self.P phi denom self.lam phi.T P_phi K P_phi / denom # 更新参数 self.theta self.theta K * (y - phi.T self.theta) # 更新协方差矩阵 self.P (self.P - np.outer(K, P_phi.T)) / self.lam return self.theta这个类里delta是初始协方差的对角值取大一点表示初始参数不确定性大。lam是遗忘因子越接近1越稳定。实际使用时回归向量phi由历史输入输出组成比如phi [y(k-1), y(k-2), u(k-1), u(k-2)]。每来一组新数据就调一次update参数会慢慢收敛到当前工况下的最优值。注意在线辨识一定要加参数变化率限幅。如果某次数据异常导致参数跳变可能直接让控制器发散。我一般限制每次更新参数变化不超过上一值的5%。这套流程跑下来你会发现系统辨识不再是教材里那堆公式而是一个从实验到模型再到控制器的完整链条。我自己的习惯是每接手一个新对象先花半小时做阶跃实验用最小二乘估一阶模型然后拿IMC公式算一组PID参数作为起点再根据现场响应微调。这套做法不一定最优但能保证你每次都有个靠谱的起点而不是从零开始试凑。希望帮到你。本文还有配套的精品资源点击获取