系统辨识实验报告

1、系统的数学建模与辨识实验报告编号： 8 姓名： 学号： 1 电加热炉动态特性辨识实验报告一 实验目的通过实验了解辨识方法在工程应用中的一些实际问题；了解数据获取和数据处理的各种方法和手段，掌握各种辨识方法的应用特点。二 实验内容1. 数据获取 按照电加热计算机控制系统使用说明，加热温箱并实时监测温度。当温度升高到一定高度基本不再发生变化时，加入辨识信号。本设计辨识信号采用伪随机二维序列中的M序列，由于设计的M序列过长，从中截取部分作为辨识信号。设定辨识信号的采样周期和整个辨识信号循环次数，加入辨识信号记录温度随电压变化的数据结果。2. 数据处理 观察实验数据，可发现扰动较严重，需要对数据进行预

2、处理以便进行后续辨识。对采集数据进行滤波等处理，实验输入为M序列长度为30，由于实际采集周期为2秒，所得数据远大于需要用到的数据，同时为了提高精度，减少扰动因素影响，数据处理过程中，当输入一个M序列元素后，采集到的数据中取最后5个数据求其平均作为该输入的输出结果。数据处理最终结果为30组输入和输出，将其生成为.mat文件以供后续离线辨识使用。3. 离线辨识利用处理过的数据选择某种辨识方法进行参数估计，并判断阶次及迟滞。辨识完成后，进行模型验证。三 实验操作 1设定相关实验参数，启动试验设备。本实验设置的加热电压为60V，温度采样间隔为2s，开始实验，加热温箱系统，并通过客户端软件实时观察温箱温

3、度，并记录，等待温箱温度升高到一定高度不再有明显上升变化。2. 设计辨识信号考虑到实验加热电压为60V，整个系统的加热电压为之前设置的系统加热电压与辨识信号电压的和，为了在较短时间内也能产生一定的变化，辨识电压值不宜过小，选取辨识信号电压为10V。此时，加入辨识信号后的整个系统电压在50V和70V来回跳动。辨识序列采用M序列，并使用庞中华系统辨识程序中的M序列生成程序，序列长度设置为31，结果如下图所示：3. 当第一步的温箱电压升高到一定高度，不再有明显升上变化时，便可认为系统已经稳定。点击加入辨识序列，设置辨识信号的采样周期为30，整个辨识信号序列循环加入的次数为2，点击确定，显示辨识信号加

4、入进度。辨识信号加入时间为2*30*30*2（s）共计1小时，信号加入完成后导出.csv文件，里面存放的是温箱温度随加入辨识信号后的控制电压变化数值，用它在MATLAB下进行后续的离线辨识。四 系统辨识对未知参数及模型的系统辨识，首先需要确定辨识方法，再判断阶次和时滞，确定模型并辨识参数。最后验证辨识结果是否合适。通过观察温箱温度随电压变化图形，本实验选用增广矩阵最小二乘法（RELS）进行辨识参数。1. 模型阶次及时滞判断本实验采用F检验准则定阶，借助数理统计中的F检验法，引入统计检验准则:tn1,n2=Jn1-Jn2Jn2N-2n22(n2-n1)当N足够大时，统计量t渐进的服从F分布。设计

5、中采用前20组数据用于判断阶次和时滞参数，SISO系统中，由于选用RELS进行参数辨识，当阶次每递增1阶，被估参数增加3个，取置信度为0.05时，查F分布表F0.05（3,20）=3.10。因此可以通过增加n大小，计算t值判断阶次。时滞d的辨识与阶次辨识同步，对任一设定的阶数n，分别假设不同的d值进行参数估计，比较残差平方总和，使残差总和最小的d值作为时滞，时滞不宜过高，本实验检测d10。模型阶及时滞辨识框图：程序见附件，输入实验数据测试结果为模型阶次n=2，阶次为2的情况下判断时滞d去不同值时的残差进行比较，发现d=1时残差值最小，因此取时滞d=1。2. 参数估计参数估计选用RELS方法辨识

6、，根据阶次判断知na=nb=nc=2,时滞d=1。增广矩阵法RELS是一种用于实时过程控制中系统参数估计得较好方法，可同时获得系统参数和噪声模型的参数估计。RELS执行步骤：1) 输入系统阶次n、时滞d和数据；2) 置初值（N）和P（N），输入初始数据；3) 采样当前输入和输出，并估算（k）；4) 计算（k）和P（k）；Kk=Pk-1x(k)1+xkPk-1x(k)k=k-1+Kkyk-x(k)(k-1)Pk=I-KkxkP(k-1)5) 返回（2）直至收敛或运行完所有采样数据；增广矩阵法程序流程如下图所示：编写程序，输入实验数据进行参数估计，计算所得使用RELS进行参数辨识，辨识结果从wor

7、kspace菜单中的变量值中读取；二阶系统所得参数为a1=-0.6126，a2=-0.3880，b0=0.0390，b1=-0.0041，b2=-0.0384，c1=-0.0062，c2=0.0064。模型验证结果如下图所示：结果显示，使用增广矩阵最好二乘法获得的辨识结果与原始数据拟合结果较好，较好的辨识出系统的特性。附件阶次及时滞辨识程序：clear all; clc; L=20; %选择20组数据用于判断阶次load(tian.mat) %导入数据B=tian;N1,N2=size(B);u=B(:,1);y=B(:,2); %将tian.mat文件中的二维数组分别用u、y表示H1=zer

8、os(L,2);for i=1:L H1(i,1)=-y(i); H1(i,2)=u(i);endestimate=inv(H1*H1)*H1*(y(2:L+1);e=(y(2:L+1)-H1*estimate;V1=e*e; %计算阶数为1时的残差平方和j=2; while(j10)H2=zeros(L,2*j);H4=;H10=;for i=1:L H1=fliplr(-y(i:i+j-1); H3=fliplr(u(i:i+j-1); H2=H1;H3; H5=H4;H2; H4=H5;endestimate=inv(H5*H5)*H5*(y(j+1:L+j);e=(y(j+1:L+j)

9、-H5*estimate;V2=e*e; %计算阶数为j时的残差平方和t1=(V1-V2)/V2*(L-2*2)/2 %F检验准则V1=V2; %递推循环for d=1:5 %计算迟滞d=1:5的残差平方和，以便从中判断时滞大小 H6=;H7=;H8=;H9=;H10=;for s=1:L H6=fliplr(-y(s:s+j-1); H7=fliplr(u(s+d:s+j-1+d); H8=H6;H7; H9=H10;H8; H10=H9;endestimate1=inv(H9*H9)*H9*(y(j+1:L+j);e=(y(j+1:L+j)-H9*estimate1;V5(d)=e*e/L

10、 %判断时滞，输出所有的V5值，通过数值判断使V5最小的d作end %为该阶次下的时滞 if t13.10 %判断阶次，t13.10时，j-1的值作为阶次，并跳出循环 N=j-1 %F分布表中，F0.05(3，20)=3.10 break endj=j+1;end系统辨识程序：clear all; clc;d=1; na=2; nb=2; nc=2; %根据阶次检测程序求得的阶次为2阶，时滞d=1L=29;load(tian.mat) %调入数据B=tian;N1,N2=size(B);u=B(:,1);y=B(:,2);xik=y;thetae_1=zeros(na+nb+1+nc,1);

11、%参数赋初值P=106*eye(na+nb+1+nc);for k=na+1:L %RELS计算，计算完所有数据 phie=fliplr(xik(k-na:k-1);u(k+d-nb:k+d);(-y(k-na:k-1); K=P*phie/(1+phie*P*phie); thetae(:,k)=thetae_1+K*(y(k)-phie*thetae_1); P=(eye(na+nb+1+nc)-K*phie)*P; xik(k)=y(k)-phie*thetae(:,k); thetae_1=thetae(:,k); endfigure(1) %参数c1,c2变化曲线plot(1:L,t

12、hetae(1:nc,:);xlabel(k); ylabel(参数估计c);legend(c_2,c_1); figure(2) %参数b0,b1,b2变化曲线plot(1:L,thetae(nc+1:nc+nb+1,:);xlabel(k); ylabel(参数估计b);legend(b_2,b_1,b_0); figure(3) %参数a1,a2变化曲线plot(1:L,thetae(nc+nb+2:nc+nb+na+1,:);xlabel(k); ylabel(参数估计a);legend(a_2,a_1); d(1)=y(1); %根据辨识结果绘制辨识后曲线d(2)=y(2);m(1)=0;m(2)=0;for k=na+1:L d(k)=0.6126*d(k-1)+0.3880*d(k-2)