简介这份系统辨识教案面向自动化专业学生与相关工程技术人员围绕从输入输出数据中建立系统模型这一核心问题系统梳理辨识理论、方法与实验流程。内容涵盖实体与模型的概念区分、模型分类、Zadeh 辨识定义及数据、模型类、准则三要素并完整展开辨识的一般步骤明确目的、收集先验知识、设计试验、数据采集与预处理、选择模型类、结构辨识与参数估计、模型检验、模型转换及应用评价还以长网造纸过程定量与水份控制为例贯穿讲解。资源包内含 1 个 doc 文档约 1.42MB结构紧凑便于通读与打印。已有 148 人学习适合作为课程自学、复习备考或课程设计的参考材料帮助读者掌握最小二乘等数据拟合方法独立完成辨识实验与报告撰写。1. 从一份 40 学时教案说起系统辨识到底在解决什么问题带过几届自动化专业的课程设计后我发现一个很普遍的现象学生能背出最小二乘的公式却拿到一个数据文件不知道从哪下手。这份《系统辨识与建模》教案就是冲着这个痛点来的——它不是一本理论推导手册而是一份 40 学时的完整教学路线图从辨识三要素讲到 M 序列设计再落到数据预处理和模型检验每一步都配了可操作的流程。教案的核心逻辑很清晰先让你理解“实体与模型”的关系再教你用数据拟合的方法从输入输出记录中反推出一个能用的数学模型。它适合两类人——正在做辨识课程设计、需要一份能照着走的流程参考的在校生以及工作中需要快速建立被控对象模型、但没时间啃完 Ljung 那本厚书的工程师。考核方式也很有意思给一个数据文件编程辨识写报告。这意味着整份教案的终点是“你能独立跑通一次辨识”。2. 辨识三要素怎么落地数据、模型类与准则的工程选择2.1 数据从白噪声到 M 序列激励信号选型逻辑教案里花了大量篇幅讲随机过程、相关函数和功率谱这些不是数学装饰而是直接决定你采集到的数据能不能用。核心结论就一条要充分激励系统输入信号的频谱必须宽于系统的工作频带。白噪声的谱密度是常数理论上最理想但物理上不存在工程中常用的是 M 序列。M 序列是二位式最大长度伪随机序列用 n 级移位寄存器产生周期 N2^n-1。它的自相关函数近似 δ 函数所以以 M 序列为输入的线性系统其输入输出互相关函数序列就等于脉冲响应序列。这个性质在工程上非常值钱——你不需要做复杂的解卷积直接算互相关就能拿到系统的脉冲响应。实际设计时教案给了两套参数选择方法。第一套先确认系统过渡过程时间 T选最高工作频率 fmax(3~15)/T然后根据 1/Ts3·fmax 得出步宽 TsT/(3k)再由 4.5k2^n-1 求出阶次 n。第二套更直接选 n 后按 n·Ts(0.1~0.8)T 反算 Ts。我一般倾向第二套因为 n 通常受硬件移位寄存器位数限制先定 n 再算 Ts 更顺手。对于多输入单输出系统教案给了一个正交化技巧将同一 M 序列的不同移位序列作为各输入信号的扰动只要输出对各输入的过渡过程小于移位间隔 J就可以认为各输入信号之间是互相正交的。这招在多变量辨识里能省掉大量解耦工作。注意M 序列的幅值由安全工况确定对可中断生产的系统且试验不引起原材料浪费的可实施单独试验。连续生产线上做辨识幅值一定要先跟工艺工程师确认。2.2 模型类与准则为什么你的最小二乘结果总是不对教案把误差准则分成了三类参数误差、方程误差和输出误差。这个分类直接对应了不同的模型结构和算法选择。方程误差 J∑[e(k)]²/N其中 e(k)y(k)a·y(k-1)-b·u(k-1)-c·u(k-2)。它的特点是 J 是参数的线性函数所以能用最小二乘直接求解计算量小。但它只能作一步预测适合模型为 AYB₁U₁B₂U₂ 的形式。很多学生用最小二乘跑出结果后发现模型预测多步就发散根源就在这里——方程误差准则下的模型只保证一步预测最优。输出误差 J∑[e(k)]²/N其中 e(k)y(k)-ŷ(k)ŷ(k) 由模型递推计算。J 是参数的非线性函数需要迭代求解计算量大但当输入序列已知时可作多步预测。如果你要做仿真训练系统或预报预测必须用输出误差准则。教案还提了一个容易被忽略的点模型阶次越高J 值越小导致模型过估计。这是过拟合在系统辨识里的具体表现。解决办法不是盲目降阶而是结合模型检验步骤用交叉验证或信息准则来判断。采样间隔对模型结构的影响也值得注意采样间隔增大一倍模型阶次大致减少一半。所以如果你发现辨识出来的模型阶次高得离谱先检查采样间隔是不是设得太小。教案给了一个修正公式J∑[e(k)]ᵀW[e(k)]·T₀²/NT₀ 为采样步长。要做不同采样步长下模型的比较必须用这个修正后的误差公式。2.3 从数据文件到模型一次完整辨识的操作流程假设你拿到一个数据文件里面是输入输出序列采样间隔已知。下面是我一般会走的流程每一步都对应教案里的一个环节。第一步数据预处理。教案列了四项传感器非线性校正、滤波、重抽样、标准化。标准化最容易被跳过但它直接影响数值稳定性。import numpy as np # 假设 data 是 N×2 数组第一列输入 u第二列输出 y u data[:, 0] y data[:, 1] # 零均值化 u_mean np.mean(u) y_mean np.mean(y) U u - u_mean Y y - y_mean # 单位方差化 sigma_u np.std(U) sigma_y np.std(Y) U_std U / sigma_u Y_std Y / sigma_y # 注意标准化后的模型参数需要反演回原始量纲 # 若标准化模型为 Y_std(k) a·Y_std(k-1) B·U_std(k-1) C·U_std(k-2) # 则原始参数 b B * sigma_y / sigma_u, c C * sigma_y / sigma_u标准化这步的逻辑是降低舍入误差防止病态方程提高计算精度。教案里给了完整的反演公式标准化后的参数乘以 σy/σu 才能还原到原始量纲。我见过有人标准化后直接拿参数去物理解释结果差了三个数量级。第二步选择模型类。根据辨识目的来定。验证理论模型要求零极点、结构、参数都准确模型类同理论模型设计常规控制器要求动态响应特性、零极点、时延准确选便于分析的模型类设计数字控制器则选便于计算机运算的模型类。教案里列了六种辨识目的对应的模型类选择做课程设计时直接查表就行。第三步结构辨识与参数估计。以最常用的 ARX 模型为例import numpy as np def estimate_arx(U, Y, na, nb): ARX 模型: y(k) a1*y(k-1) ... ana*y(k-na) b1*u(k-1) ... bnb*u(k-nb) 用最小二乘估计参数 N len(Y) start max(na, nb) # 构造回归矩阵 Phi [] Y_target [] for k in range(start, N): row [] for i in range(1, na 1): row.append(-Y[k - i]) for j in range(1, nb 1): row.append(U[k - j]) Phi.append(row) Y_target.append(Y[k]) Phi np.array(Phi) Y_target np.array(Y_target) # 最小二乘: theta (Phi^T Phi)^(-1) Phi^T Y theta, residuals, rank, sv np.linalg.lstsq(Phi, Y_target, rcondNone) a_params theta[:na] b_params theta[na:] return a_params, b_params这段代码的关键在回归矩阵的构造每一行对应一个时刻 k前 na 列是过去的输出取负号后 nb 列是过去的输入。最小二乘解就是让方程误差平方和最小的参数。np.linalg.lstsq比直接求逆更稳定内部用了 SVD 分解。第四步模型检验。教案强调用输出估计误差是否具有白色性来判断辨识方法的优劣。具体做法是计算残差的自相关函数如果残差在非零延迟处接近零说明模型结构合适如果残差还有明显相关性说明模型阶次不够或结构选错了。3. 辨识实验设计从造纸过程案例拆解可复现步骤3.1 案例背景与变量关系梳理教案用了一个长网造纸过程的例子这个案例的变量关系很典型成品纸的定量 G 和水份 W 与纸浆浓度 D、纸浆流量 F、车速 V、蒸汽压力 P 有关即 Gf(D,F,V,P)Wg(D,F,V,P)。目标是建立数学模型用于计算机控制。现场调查发现的关键先验知识车速调整存在同步困难通常设为恒定流量改变到定量改变存在约 60 秒延迟响应过程约 2 秒浓度改变到定量改变存在约 120 秒延迟响应过程约 80 秒蒸汽压力改变到水份改变存在约 45 秒延迟响应过程约 60 秒。这些数字直接决定了后续实验设计的参数。流量通道的过渡过程只有 2 秒而浓度通道有 80 秒两者相差 40 倍。如果对这两个通道施加相同步宽的 M 序列流量通道的数据会严重冗余而浓度通道可能还没激励充分。3.2 分块设计与扰动信号参数计算教案的做法是分块将流量和蒸汽压力对定量、水份的影响作为一个分块试验时保持车速和纸浆浓度不变。对流量和蒸汽压力分别施加伪随机序列扰动幅度以不引起断纸为限采样频率设为 2试验时间 1000 秒。这里“采样频率为 2”指的是采样间隔 T₀0.5 秒。根据香农定理1/T₀2·fmax。流量通道过渡过程 T≈2 秒fmax≈(3~15)/21.5~7.5 Hz1/T₀2 Hz 满足要求。但浓度通道过渡过程 T≈80 秒fmax≈0.04~0.19 Hz2 Hz 的采样频率远远过采样了。教案在数据预处理里提了重抽样就是针对这种情况——对过渡过程较长的子系统用重抽样重新选择采样间隔以降低模型阶次。M 序列参数计算以流量通道为例T≈2 秒选 n5则 N31按 n·Ts(0.1~0.8)T 得 Ts(0.1~0.8)×2/50.04~0.32 秒。取 Ts0.1 秒则一个周期时长 31×0.13.1 秒满足 1.5T(2^n-1)Ts 即 33.1。试验时间 1000 秒约等于 322 个周期远超 3~5 个周期的要求。对于多输入的情况教案建议将同一 M 序列的不同移位序列作为各输入的扰动。假设流量和蒸汽压力两个通道可以用 M(k) 和 M(kJ) 作为两个扰动信号。J 的选取要满足 J·Ts1.5×max(过渡过程时间)。流量通道过渡过程 2 秒蒸汽压力通道 60 秒取 J·Ts90 秒即 J900。但 N31J 不能超过 N所以需要增大 n 或 Ts。这就是教案里说的“(2^n-1-1)·Ts1.5·∑(输出对各输入的过渡过程时间)”的约束。3.3 现场准备与数据采集的检查清单教案把现场准备单独列为一个步骤包括向企业领导申请试验时段、准备扰动信号发生器并预发信号检验、测验现场信号干扰情况、准备模数转换设备调好零迁和放大参数、现场接线。这几件事里最容易翻车的是信号零迁和放大。工业现场的传感器输出通常是 4~20mA 或 1~5V而数据采集卡的输入范围可能是 ±10V。如果不做零迁和放大有效信号可能只占 ADC 量程的 10%量化噪声会直接淹没小信号。我一般会在接线后先用万用表量一遍每个通道的静态电压确认在 ADC 量程的 30%~70% 之间再开始采集。数据采集时教案要求存盘并编写数据说明文件。数据说明文件至少包含采样间隔、通道编号与物理量对应关系、各通道量程与单位、试验工况、扰动信号参数。这份文件在后续辨识时是“后悔药”——没有它三个月后你打开数据文件根本不知道哪列是流量哪列是压力。4. 避坑与排查辨识课程设计里最常见的五个翻车点4.1 现象最小二乘跑出来的参数完全不可用模型输出发散原因最常见的是数据没有零均值化。教案在标准化一节里明确写了零均值化是标准化的第一步。如果输入输出含有较大的直流分量回归矩阵的条件数会急剧增大最小二乘解对噪声极其敏感。另一个原因是模型阶次选得过高导致过拟合参数虽然拟合了训练数据但模型极点跑到单位圆外递推几步就发散。解决先做零均值化和单位方差化再跑最小二乘。阶次选择从低到高试每次增加阶次后检查残差自相关函数直到残差不再有显著相关性为止。教案里“模型阶次越高 J 值越小导致过估计”的提醒就是针对这个问题的。4.2 现象M 序列作为输入时辨识出的脉冲响应尾部噪声很大原因M 序列的周期 N 不够大或者试验时间没有覆盖足够的周期数。教案指出 M 序列的直流分量与 N² 成正比N 不够大时直流分量会污染互相关函数的尾部。另外如果试验时间不足 3~5 个扰动信号周期互相关函数的估计方差会很大。解决增大 n 使 N 增大或者增加试验时间。教案建议持续时间为 3~5 个扰动信号周期预实验 1~2 个周期。实际操作中我一般会取 5 个周期以上尾部噪声能明显改善。4.3 现象多输入辨识时各通道的脉冲响应互相串扰原因各输入信号的扰动没有正交化。如果两个输入同时施加相同的 M 序列输出中的响应无法分离。教案给的正交化方法是同一 M 序列的不同移位序列但移位间隔 J 必须大于输出对各输入过渡过程时间的 1.5 倍。解决先确认各通道的过渡过程时间取最大的那个乘以 1.5 再除以 Ts 得到最小 J。如果 N 不够大导致 J 无法满足需要增大 n 或 Ts。另一个办法是分时施加扰动——先激励通道 1其他通道保持恒定采集完再换通道 2。这样虽然试验时间长但数据质量最好。4.4 现象数据预处理后模型精度反而下降了原因滤波器的选择或应用方式不对。教案强调“对于输入输出加相同的滤波器输入输出关系不变”。如果只对输出滤波而输入不滤波会改变系统的频率特性辨识出的模型对应的是滤波后的系统而不是原系统。重抽样时如果没有先做抗混叠滤波高频噪声会混叠到低频段。解决输入输出必须经过相同的滤波器。重抽样前先做低通滤波截止频率设为新采样频率的一半以下。标准化后的参数记得反演回原始量纲否则物理解释会出错。4.5 现象模型检验通过但实际控制效果差原因辨识目的与模型类选择不匹配。教案里列了六种辨识目的对应的模型类要求比如设计常规控制器要求动态响应特性、零极点、时延准确便于分析的模型类而监视过程参数实现故障诊断只要求参数准确能直观体现被监视过程参数的模型类。如果拿一个用于故障诊断的模型去设计控制器动态特性可能完全不对。解决回到辨识目的重新确认模型类选择。验证理论模型要求零极点、结构、参数都准确模型类同理论模型设计数字控制器要求动态响应特性、时延准确便于计算机运算的模型类。目的不同模型类的选择标准完全不同。5. 从教案到实战用互相关法估计脉冲响应的完整实现教案里有一个很实用的技巧以 M 序列为输入的线性系统其输入输出互相关函数序列等于脉冲响应序列。这个方法的工程价值在于它不需要迭代求解一次互相关计算就能拿到非参数模型。下面是我在实际项目中常用的实现。import numpy as np def estimate_impulse_response(u, y, N_period, n_shift1): 用 M 序列作为输入通过互相关估计脉冲响应 u: 输入序列M 序列取值 ±a y: 输出序列 N_period: M 序列周期长度 N 2^n - 1 n_shift: 移位步数用于多输入正交化 返回: 脉冲响应估计 g_hat a np.max(np.abs(u)) # M 序列幅值 N len(u) # 计算互相关函数 R_uy(tau) (1/N) * sum(u(k) * y(ktau)) max_lag min(N_period, N // 3) R_uy np.zeros(max_lag) for tau in range(max_lag): R_uy[tau] np.mean(u[:N - tau] * y[tau:]) # 对于 M 序列R_uu(0) a^2, R_uu(tau≠0) ≈ -a^2/N # 脉冲响应估计: g(tau) R_uy(tau) / R_uu(0) 的修正形式 # 更精确的公式: g(tau) N / ((N1)*a^2) * R_uy(tau) g_hat N / ((N 1) * a**2) * R_uy return g_hat # 使用示例 # 假设 u 是 M 序列y 是系统输出N_period31n5 # g_hat estimate_impulse_response(u, y, N_period31) # 然后对 g_hat 做阶跃响应的积分即可得到阶跃响应这段代码的核心是互相关函数的计算。M 序列的自相关函数在零延迟处为 a²在非零延迟处约为 -a²/N。所以互相关函数除以自相关函数的零延迟值再乘以 N/(N1) 的修正系数就得到脉冲响应的估计。修正系数来自 M 序列自相关函数的精确表达式当 N 充分大时趋近于 1。拿到脉冲响应后可以用它来验证参数模型的准确性——如果参数模型的脉冲响应和互相关法估计的脉冲响应在过渡过程时间段内吻合说明参数模型结构选对了。教案里“模型检验”一步的很多方法本质上都是在比较非参数模型和参数模型的一致性。还有一个细节互相关法的估计方差与试验时间成反比。教案建议 3~5 个扰动信号周期我一般会取 5~10 个周期。如果发现脉冲响应尾部噪声大先别急着换算法增加试验时间往往比换算法更有效。从那以后我每次做辨识实验都会先用互相关法快速估一版脉冲响应确认过渡过程时间和延迟跟先验知识对得上再往下走参数估计。这一步花不了多少时间但能提前发现数据采集阶段的问题避免在参数估计上反复调参。希望帮到你。本文还有配套的精品资源点击获取