Matlab系统辨识工具箱

Matlab系统辨识工具箱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-1,M]，表示信号周期长为(22logp-1)，且在1/M间隔内信号幅值不变；band缺省值为[0,1]，即生成白噪声信号。levels为1×2行向量，用来决定输入信号幅值的上下界，当信号类型为rs时，下界为高斯信号的均值减1，而上界为高斯信号的均值加1。2.数据预处理在系统辨识工具箱中，函数dtrend用来去除输入输出数据中的趋势项，函数idfilt利用Butterworth滤波器对数据进行滤波，利用idresamp函数对输入输出数据进行重新采样来插值或删减数据。这3个函数的基本调用格式分别为：zd＝dtrend(z,o,brkp)[zf,thf]＝idfilt(z,ord,Wn,hs)[zr,R_act]＝idresamp(z,R,filter_order,tol)式中：z＝[yu]为系统输入输出数据矩阵；o是趋势项的阶次，缺省值o＝0，即从数据中减去均值。o＝1，则从数据中减去线性趋势项。brkp用于包含分段线性趋势的数据。ord为滤波器的阶。Wn有几种情况：①当未指定可选参数hs，若Wn仅包含1个元素，则该参数指定低通滤波器的截止频率。若Wn＝[WlWh]，则Wn指定带通滤波器的上下限频率；②当hs＝high时，Wn指定高通滤波器的截止频率；③当hs＝stop时，Wn＝[WlWh]指定带阻滤波器的频率范围。zf为滤波后I/O数据；thf为滤波器对应Theta模型zr为重新采样数据组；R_act为实际达到的重采样因子；R为重采样因子，重新采样周期为原采样周期的R倍。R＞1时，对应的数据均应删除，R＜1时，对应插值过程。filter_order为重新采样前预采样时所用滤波器的阶数，缺省时为8；tol指定R进行有理逼近精度，缺省值为0.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＝[yu]为模型辨识的I/O数据向量或矩阵。zv＝[yrur]为模型验证的I/O数据向量或矩阵。NN为多个模型结构参数构成的矩阵，NN的每行都具有格式nn＝[nanbnk]。v的第1行为各个模型结构损失函数值，后面的各行为模型结构参数。函数ivstruc用来基于辅助变量法的ARX模型结构损失函数的计算，调用格式为：v＝ivstruc(ze,zv,NN)其中参数同arxstruc。函数selstruc用来在损失函数的基础上进行模型结构选择，调用格式为：[nn,vmod]＝selstruc(v,c)其中v由函数arxstruc或ivstruc获得的输出矩阵，为各个模型结构的损失函数。c为可选参数，用于指定模型结构选择的方式。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＝[yu]为输出输入数据矩阵。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模型nn含义同函数armax；adm、adg同函数rarx。5.模型检验与动态仿真工具箱中，用于模型验证和仿真的函数主要有compare、resid、pe、predict和idsim。其中：compare可将模型的预测输出与对象实际输出进行比较，resid用来计算和检验模型残差；pe计算预测误差；predict预测未来输出；idsim可进行模型仿真计算。经检验，如果模型不适用，则要重新进行模型结构的选择、参数估计等工作。工具箱中还提供交互式的图形界面工具，能方便实现数据的预处理、模型类型的选择和参数估计以及模型验证和比较等功能。例1：采用最小二乘法(LS)与近似最优4阶辅助变量法(IV4)进行了模型参数辨识,其Matlab程序清单如下:clear%清零工作空间U=[0,0184,0114,-0196,-0128,0193,0199,0141]’;Z=[0,01195,01513,-01272,-0142,01058,01263,0145]′;%赋系统输入输出初始值m=iddata(Z,U);%标准包含系统输入输出数据对象V=arx(m,[111])%基于最小二乘法的系统参数辨识present(V)%显示最小二乘辨识对象的信息V1=iv4(m,[111])%基于近似4阶最优辅助变量模型参数辨识present(V1)%显示辅助变量辨识。例2：考虑如下的仿真对象，其数学模型为：y(k)-1.5y(k-1)+0.705y(k-2)=0.065u(k-1)+0.084u(k-2)-0.008u(k-3)+e(k)-0.527e(k-1)+0.01e(k-2)其中，e(k)为均值为0，方差为1，服从正态分布的不相关随机噪声；输入信号u(k)采用M序列。辨识模型采用如下形式：1212312ykyk1yk2uk1uk2uk3ekek1ek2aabbbcc数据K度取N=1000。Matlab程序如下：clearall；N=1000；A=[1-1.50.705];B=[0.0650.48-0.008];C=[1-0.5270.01];M1=idpoly(A,B,C)％理想系统模型ARMAX模型step(M1,[0100]);grid;％绘制系统阶跃响应曲线％一一一一一一一一一一一一一一一一产生输入输⋯数据U=iddata([],idinput(N,'prbs'));％伪随机序列E=iddata([],idinput(N,'rgs'));％白噪声序列Y1=sim(M1,[U,E])；％输出数据Z=iddata(Y1,U);％一一一一一一一一一一一一一RMAX模型阶次的估计NN=struc(1:2,1:4,1:4);Loss_fun=arxstruc(Z,Z,NN);order=selstruc(Loss_fun,'aic');order=[order(1),order(2),1,order(3)]％一一一一一一一一ARMAX模型参数的估计(最小二乘法)Model_para=armax(Z,order);％参数估汁present(Model_para);％显示辨识结果compare(Z,Model_para);％辨识结果与实际输出比较例3：已知待辨识的对象为：y(k)-1.5y(k-1)+0.7y(k-2)=u(k-1)+0.5u(k-2)+v(k)其中v(k)为零均值白噪声，对上述对象分别进行ARX建模、相关分析，并比较2种方法得到脉冲和阶跃响应。其M程序清单为：%二阶对象进行ARX建模和相关分析、比较clf;A=[1-0.50.7];B=[0.10.5];th0=poly2th(A,B);u=idinput(300,'rbs');y=idsim([u,randn(300,1)],th0);%z=[yu];ir=cra(z);th=arx(z,[221]);%imp=[1;zeros(19,1)];irth=idsim(imp,th);%subplot(211)plot(ir,'--');holdon;plot(irth,'r')legend('相关分析','ARX建模')title('impulseresponses')%subplot(212)plot(cumsum(ir),'--');holdonplot(cumsum(irth),'r')legend('相关分析','ARX建模')title('stepresponses')例4：已知参考模型为：y(k)-1.5y(k-1)+0.7y(k-2)=u(k-1)+0.5u(k-2)+v(k)其中：系统噪声v(k)为N(0,1)分布的白噪声，控制信号u(k)采用15拍M序列，辨识模型为：y(k)+a1y(k-1)+a2y(k-2)=b1u(k-1)+b2u(k-2)+v(k)%参数辨识clearU=idinput([15160],'prbs',[01],[-11]);i=900;R0=1;S=2.^25*eye(4);Y=zeros(i,1);theta0=[1.5;-0.7;1;0.5];theta=[1;