兵工自动化 软件开发与应用 O. I. Automation 2006 年第 25 卷第 10 期 Exploitation and Application of Software 2006, Vol. 25, No. 10 文章编号：1006-1576（2006）10-0088-03 基于 Matlab 系统辨识工具箱的系统建模 齐晓慧，田庆民，董海瑞 （军械工程学院 光学与电子工程系，河北 石家庄 050003） 摘要：Matlab 系统辨识工具箱的功能主要包括参数/非参数化模型辨识工具、模型验证工具、递推参数估计、各 种模型类的建立和转换函数及集成多种功能的图形用户界面等。故基于该工具箱进行了系统建模，主要包括获取观 测数据、数据预处理、选择模型结构、参数估计、模型检验与动态仿真等步骤。 关键词：系统建模；系统辨识工具箱；Matlab 中图分类号：TP273 文献标识码：A System Modeling Based on System Identification Toolbox in Matlab (Dept. of Optics & Electronics Engineering, Ordnance Engineering College, Shijiazhuang 050003, China) QI Xiao-hui, TIAN Qing-min, DONG Hai-rui Abstract: The function in system identification toolbox of Matlab mainly includes parametric or nonparametric model identification, model validation, regressive parameter estimation, model class definition and conversion, and graphic interface. So, based on the toolbox, the system was established models. The main steps include data collection, data preprocess, model structure selection, parameters estimation, model validation and dynamic simulation, etc. Keywords: System modeling; System identification toolbox; Matlab 0 引言 为定量分析和研究某个对象，需用数学模型描 述将其模型化。当对象较简单时，可应用有关知识 建立对象的机理模型，并达到较高精度。当对象较 复杂时，建立其机理模型困难。可通过系统辨识获 得对象的 I/O 数据，并依此建立对象的数学模型。 利用该模型，涉及到观测数据的获取、数据检验、 模型结构选择、参数估计、模型检验和修改等。由 于实际建模中存在不确定性，获取对象数据时受条 件和环境的制约，从数据获取到模型的建立经反复 探索，计算量大，用手工难以完成。故借助 Matlab 系统辨识工具箱进行建模，可简化计算过程、提高 效率，其图形功能使建模过程直观、形象。 1 系统辨识工具箱基本内容 该 工 具 箱 的 功 能 包 括 参 数 化 和 非 参 数 化 模 型 辨识工具和验证工具、递推参数估计、各模型类的 建立和转换函数，及集成多种功能的用户界面。 1.1 观测数据的获取 观测数据含输入、输出、噪声等。而工具箱提 供系统辨识的输入信号函数 idinput，调用格式为： u＝idinput (N, type, band, levels) N 为生成输入信号的数据长度。type 为输入信 号类型，包含：高斯随机信号 type＝rs，二值随机信 号 type＝rbs，二值伪随机信号 type＝prbs，正弦扫 描信号 type＝sine。band 为 1×2 行向量，即信号带 宽。当信号类型为 rs、rbs 和 sine 时，band＝[低频, 高频]，其中低频和高频为 Nyquist 标准频率，其值 在 0～1 之间；当信号类型为 prbs 时，band＝[2logp, M]，表示信号周期长为 (22log p-1)，且在 1/M 间隔 内信号幅值不变；band 缺省值为[0,1]，即生成白噪 声信号。levels 为 1×2 行向量，用来决定输入信号 幅值的上下界，当信号类型为 rs 时，下界为高斯信 号的均值减 1，而上界为高斯信号的均值加 1。 1.2 数据预处理 系统建模时，要求输入输出数据的统计特性与 统计时间起点无关，且均值为 0。而实际测量直接得 到的数据是随机时间序列，须进行数据预处理，主 要包括消除数据的趋势项、对测量数据进行滤波和 重新采样。在系统辨识工具箱中，函数 dtrend 用来 去 除 输 入 输 出 数 据 中 的 趋 势 项 ， 函 数 idfilt 利 用 Butterworth 滤波器对数据进行滤波，利用 idresamp 函数对输入输出数据进行重新采样来插值或删减数 据。这 3 个函数的基本调用格式分别为： zd＝dtrend (z, o, brkp) [zf,thf]＝idfil t(z, ord, Wn, hs) [zr, R_act]＝idresamp (z, R ,filter_order, tol) 式中：z＝[y u] 为系统输入输出数据矩阵；o 是 趋势项的阶次，缺省值 o＝0，即从数据中减去均值。 o＝1，则从数据中减去线性趋势项。brkp 用于包含 收稿日期：2006-06-08；修回日期：2006-06-30 作者简介：齐晓慧（1962-），女，辽宁人，南京理工大学在读博士，从事武器控制系统的分析、仿真与设计研究。 ·88· 

兵工自动化 软件开发与应用 O. I. Automation 2006 年第 25 卷第 10 期 Exploitation and Application of Software 2006, Vol. 25, No. 10 分段线性趋势的数据。ord 为滤波器的阶。Wn 有几 种情况：① 当未指定可选参数 hs，若 Wn 仅包含 1 个元素，则该参数指定低通滤波器的截止频率。若 Wn＝[Wl Wh]，则 Wn 指定带通滤波器的上下限频 率；② 当 hs＝high 时，Wn 指定高通滤波器的截止 频率；③ 当 hs＝stop 时，Wn＝[Wl Wh] 指定带阻 滤波器的频率范围。zf 为滤波后 I/O 数据；thf 为滤 波器对应 Theta 模型；zr 为重新采样数据组；R_act 为实际达到的重采样因子；R 为重采样因子，重新 采样周期为原采样周期的 R 倍。R＞1 时，对应的数 据均应删除，R＜1 时，对应插值过程。filter_order 为重新采样前预采样时所用滤波器的阶数，缺省时 为 8；tol 指定 R 进行有理逼近精度，缺省值为 0.1。 1.3 模型结构的选择 系统辨识工具箱提供对多种模型类的支持，包 括非参数和参数等模型类。非参数模型类包括脉冲 响 应 和 频 域 描 述 等 模 型 。 参 数 模 型 类 有 ARX 、 ARMAX、BJ、输出误差和状态空间等模型。 系统辨识工具箱提供通用参数模型结构和参数 存储与表示形式，即 Theta 模型格式。该格式以矩 阵形式存储各种参数模型结构和参数，并能方便地 与各种参数模型进行相互转换。 系 统 辨 识 工 具 箱 提 供 的 模 型 结 构 选 择 函 数 有 struc、arxstruc、ivstruc 和 selstruc。 函数 struc 生成 ARX 结构参数，调用格式为： NN＝struc (NA, NB, NK) 其中，NA、NB 分别为 ARX 模型多项式 A(q)、 B(q) 的阶次范围；NK 为 ARX 模型纯时延的大小范 围；NN 为模型结构参数集构成的矩阵。 函数 arxstruc 用来计算 ARX 模型结构的损失函 数，即归一化的输出预测误差平方和，调用格式为： v＝arxstruc (ze, zv, NN) 其中，ze＝[y u] 为模型辨识的 I/O 数据向量或 矩阵。zv＝[yr ur] 为模型验证的 I/O 数据向量或矩 阵。NN 为多个模型结构参数构成的矩阵，NN 的每 行都具有格式 nn＝[na nb nk]。v 的第 1 行为各个模 型结构损失函数值，后面的各行为模型结构参数。 函数 ivstruc 用来基于辅助变量法的 ARX 模型结 构损失函数的计算，调用格式为： v＝ivstruc (ze, zv, NN) 其中参数同 arxstruc。函数 selstruc 用来在损失 函数的基础上进行模型结构选择，调用格式为： [nn, vmod]＝selstruc (v, c) 其中 v 由函数 arxstruc 或 ivstruc 获得的输出矩 阵，为各个模型结构的损失函数。c 为可选参数，用 于指定模型结构选择的方式。 1.4 参数估计 系统辨识工具箱中，支持的参数模型包括 AR、 ARX、ARMAX、BJ、状态空间和输出误差等模型， 含一次完成和递推辨识等算法。一次完成算法的参 数模型辨识函数有 ar、arx、armax、ivx 等。用递推 算法进行参数模型辨识的函数有 rarx、rarmax 等。 函数 ar 用于 AR 模型辨识，调用格式为： th＝ar (y, n , approach) y 为对象在白噪声作用下的输出；n 为 AR 模型 的阶次；approach 指定参数估计的最小二乘类方法， 取值包括：① fb 为前向-后向方法（缺省时默认）； ② ls 为标准的最小二乘法；③ yw 为 Yule-Walker 方法；④ burg 为 Burg 方法；⑤ gl 为几何网络法。 函数 arx、ivx 用于 ARX 模型的辨识，前者为基 于最小二乘法，后者为用辅助变量法进行系统参数 的估计，调用格式分别为： th＝arx (z, nn) th＝ivx (z, nn) 其中，z＝[y u]为输出输入数据矩阵。nn＝[na, nb, nk] 分别给出 ARX 模型中多项式 A(q-1)、B(q-1) 的阶次和系统时延。 函数 armax 对 ARMAX 模型进行辨识，格式为： th＝armax (z, nn, ’trace’) 其中，z 的含义同前。nn＝[na, nb, nc, nk] 给出 ARMAX 模型中多项式 A(q-1)、B(q-1)、C(q-1)阶次和 系统时延。’trace’ 指定显示迭代优化过程信息。 rarx 为基于递推最小二乘算法的 AR 或 ARX 模 型辨识函数，调用格式为： thm＝rarx (z, nn, adm, adg) 其中 z 的含义同前。当 z＝y 时，可对 AR 模型 进行辨识；nn 为 AR 或 ARX 模型阶次，对 AR 模型 nn 为 A(q-1) 的阶次，对 ARX 模型，与函数 arx 中对 应的含义一样；adm、adg 指定采用递推最小二乘算 法的类型，分别有：① adm＝ff、adg=lam，对应具 有遗忘因子的递推最小二乘法；② adm＝ug、adg ＝gam，对应未归一化的梯度算法；③ adm＝ng、 adg＝gam，对应归一化的梯度算法；④ adm＝kf、 adg＝R1，采用基于 Kalman 滤波器的递推算法。 rarmax 为采用递推算法进行 ARMA 或 ARMAX 模型的辨识函数，用于单 I/O 系统，调用格式为： thm＝rarmax (z, nn, adm, adg) 其中，z 的含义同前。当 z＝y 时，可对 ARMA 模型进行辨识；nn 为 ARMA 或 ARMAX 模型的阶 次，对 ARMA 模型 nn＝[na, nc]，对 ARMAX 模型 ·89· 

兵工自动化 软件开发与应用 O. I. Automation 2006 年第 25 卷第 10 期 Exploitation and Application of Software 2006, Vol. 25, No. 10 nn 含义同函数 armax；adm、adg 同函数 rarx。 1.5 模型检验与动态仿真 plot(cumsum(irth),'r') legend('相关分析','ARX 建模') title('step responses') 由给定的观测数据建立系统数学模型后，还需 要进行检验，看模型是否适用，如果不适用，则要 修改模型结构，重新进行参数估计等。 工具箱中，用于模型验证和仿真的函数主要有 compare、resid、pe、predict 和 idsim。其中：compare 可 将 模 型 的 预 测 输 出 与 对 象 实 际 输 出 进 行 比 较 ， resid 用来计算和检验模型残差；pe 计算预测误差； predict 预测未来输出；idsim 可进行模型仿真计算。 经检验，如果模型不适用，则要重新进行模型结 构的选择、参数估计等工作。 工具箱中还提供交互式的图形界面工具，能方 便实现数据的预处理、模型类型的选择和参数估计 以及模型验证和比较等功能。 2 建模举例及结论 利用系统辨识工具箱提供的辨识函数，可很方 便地完成系统建模工作。 1 0 Impulse responses 相关分析 ARX 建模 -1 2 1 0 5 10 15 20 25 Step responses 相关分析 ARX 建模 0 0 5 10 15 20 25 图 1 ARX 建模和相关分析法得到的响应曲线 [算例 1] 已知待辨识的对象为： )k(v)2k(u5.0)1k(u − 其中 v(k) 为零均值白噪声，对上述对象分别进 行 ARX 建模、相关分析，并比较 2 种方法得到脉 冲和阶跃响应。仿真曲线如图 1，其 M 程序清单为： )2k(y7.0)1k(y5.1)k(y +− +− + − − = % 二阶对象进行 ARX 建模和相关分析、比较 clf; A=[1 -0.5 0.7]; B=[0 1 0.5]; th0=poly2th(A,B) ; u=idinput (300,'rbs'); y=idsim ([u, randn(300,1)],th0); % z=[y u]; ir= cra(z); th=arx (z,[2 2 1]) ; % imp=[1;zeros(19,1)]; irth=idsim(imp,th); % subplot(211) plot(ir,'--') ;hold on;plot(irth,'r') legend('相关分析','ARX 建模') title('impulse responses') % subplot(212) plot(cumsum(ir),'--') ; hold on ·90· k(y7.0)1k(y5.1)k(y [算例 2] 已知参考模型为： − 其中：系统噪声 v(k) 为 N(0,1) 分布的白噪声， )k(v)2k(u5.0)1k(u +− +− )2 + − = − 控制信号 u(k) 采用 15 拍 M 序列，辨识模型为： )k(y +− k(y2a )1k(y1a + 辨识的参数曲线如图 2，M 程序清单为： k(u2b)1k(u1b +− )k(v )2 )2 + − = − % 参数辨识 clear U=idinput([15 1 60],'prbs',[0 1],[-1 1]); % i=900; R0=1; S=2.^25*eye(4); Y=zeros(i,1); theta0=[1.5;-0.7;1;0.5]; % theta=[1;1;1;1]; V=idinput(i,'rgs',[0 1],[0 1]); % for m=3:i X=[Y(m-1) Y(m-2) U(m-1) U(m-2)]; Y(m)=X*theta0+V(m); F=S'*X'; Beta=R0+F'*F; Arfa=1/(Beta+(R0*Beta)^0.5); K=S*F/Beta; theta(:,(m-1))=theta(:,(m-2))+K*(Y(m)-X*theta(:,(m-2))); S=(eye(4)-Arfa*Beta*K*X)*S/R0^0.5; end % result=[(-theta(1:2,(i-1)));theta(3:4,(i-1))]; % figure; plot(1:(i-1),-theta(1,:),'--',1:(i-1),-theta(2,:),'-.',1:(i-1), theta(3,:),':',1:(i-1),theta(4,:),'-'); A=plot(1:(i-1),-theta(1,:),'--',1:(i-1),-theta(2,:),'-.',1:(i-1), theta(3,:),':',1:(i-1),theta(4,:),'-'); legend(A,'a1','a2','b1','b2'); title('参数辨识曲线'); hold on 可见，借助 Matlab 系统辨识工具箱可以使建模 变得简单、方便。 1 0 -1 -2 b1 a2 b2 a1 0 200 400 600 800 图 2 辨识的参数曲线 参考文献： [1] 李言俊, 张科. 系统辨识理论及应用[M]. 北京: 国防工 业出版社, 2003. [2] 王江涛. Matlab 在自适应控制中的应用[J]. 电工技术与 自动化, 2003, (4): 51-53. [3] 徐昕, 李涛, 等. Matlab 工具箱应用指南[M]. 北京: 电 子工业出版社, 2000.