系统辨识试验报告

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

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

3、度采样间隔为2s，开始实验，加热温箱系统，并通过客户端软件实时观察温箱 温度，并记录，等待温箱温度升高到一定高度不再有明显上升变化。2. 设计辨识信号考虑到实验加热电压为60V，整个系统的加热电压为之前设置的系统加热电 压与辨识信号电压的和， 为了在较短时间内也能产生一定的变化， 辨识电压值不宜过小，选取辨识信号电压为10V。此时，加入辨识信号后的整个系统电压在 50V 和70V来回跳动。辨识序列采用 M序列，并使用庞中华系统辨识程序中的 M序列生成程序，序列长度设置为 31,结果如下图所示:30值幅列序MM序列205101520k2530O-23.当第一步的温箱电压升高到一定高度，不再有明显

4、升上变化时，便可 认为系统已经稳定。点击加入辨识序列，设置辨识信号的采样周期为30，整个辨识信号序列循环加入的次数为 2,点击确定，显示辨识信号加入进度。辨识信 号加入时间为2*30*30*2 (s)共计1小时，信号加入完成后导出.csv文件，里面 存放的是温箱温度随加入辨识信号后的控制电压变化数值，用它在MATLAB下进行后续的离线辨识。四.系统辨识对未知参数及模型的系统辨识，首先需要确定辨识方法，再判断阶次和时滞， 确定模型并辨识参数。最后验证辨识结果是否合适。通过观察温箱温度随电压变 化图形，本实验选用增广矩阵最小二乘法(RELS进行辨识参数。1.模型阶次及时滞判断本实验采用F检验准则定

5、阶，借助数理统计中的F检验法，引入统计检验准J(nJ - J(n2) N- 2n2t(n i, n2)=7J(n2)2(n2 - n”当N足够大时，统计量t渐进的服从F分布。设计中采用前20组数据用于 判断阶次和时滞参数，SISO系统中，由于选用RELS4行参数辨识，当阶次每递 增1阶，被估参数增加3 个,取置信度为0.05时，查F分布表Fo.5(3,20)=3.10。 因此可以通过增加n大小，计算t值判断阶次。时滞d的辨识与阶次辨识同步， 对任一设定的阶数n,分别假设不同的d值进行参数估计，比较残差平方总和， 使残差总和最小的d值作为时滞，时滞不宜过高，本实验检测d10。模型阶及时滞辨识框图

6、：程序见附件，输入实验数据测试结果为模型阶次n=2,阶次为2的情况下判断时滞d去不同值时的残差进行比较，发现d=1时残差值最小，因此取时滞d=1。2.参数估计参数估计选用RELS方法辨识，根据阶次判断知 na=nb=nc=2,时滞d=1。增广 矩阵法RELS是一种用于实时过程控制中系统参数估计得较好方法，可同时获得 系统参数和噪声模型的参数估计。RELS行步骤：1) 输入系统阶次n、时滞d和数据；2) 置初值B( N)和P (N),输入初始数据;3) 采样当前输入和输出，并估算三(k)；4) 计算 B( k)和 P (k)；K(k)=P(k- 1) x(k)1 + x(k)P(k- 1 )x(

7、k)G(k) = 9(k- 1) + K(k)y(k) - x(k) 9(k 1)P(k) = I - K(k)x(k) P(k - 1)5) 返回(2)直至收敛或运行完所有采样数据；增广矩阵法程序流程如下图所示：编写程序，输入实验数据进行参数估计，计算所得使用RELS4行参数辨识, 辨识结果从workspace菜单中的变量值中读取；二阶系统所得参数为ai=-0.6126, a2=-0.3880, bo=O.O39O, bi=-0.0041，b2=-0.0384, ci=-0.0062, C2=0.0064。模型验证结果如下图所示：结果显示，使用增广矩阵最好二乘法获得的辨识结果与原始数据拟合结

8、果较 好，较好的辨识出系统的特性。附件阶次及时滞辨识程序：clear all ;clc;L=20;load( tian.mat )B=tian; N1,N2=size(B);u=B(:,1);y=B(:,2);H1=zeros(L,2);for i=1:L%选择20 组数据用于判断阶次%导入数据%将tian.mat文件中的二维数组分别用 u、y表示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;j=2;while (j10)%计算阶数为 1时的残差平方和

9、H2=zeros(L,2*j);H4=;H10=; for i=1:LH1=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)-H5*estimate;V2=e*e;t1=(V1-V2)/V2*(L-2*2)/2V1=V2;for d=1:5H6=;H7=;H8=;H9=;H10=;%计算阶数为 j 时的残差平方和%验验准则%递推循环%计算迟滞 d=1:5 的残差平方和，以便从中判断时滞大小for s=1:LH6=f

10、liplr(-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/Lend%判断时滞，输岀所有的 V值，通过数值判断使 V最小的d作%为该阶次下的时滞if t13.10%判断阶次， t13.10 时， j-1 的值作为阶次，并跳岀循环N=j-1break%F分布表中，Fo.o5 (3 , 20)=3.10endj=j+1; end系统辨识程序：clear all ; clc

11、;d=1;na=2; nb=2; nc=2;L=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);P=10A6*eye (n a+nb+1+ nc);%根据阶次检测程序求得的阶次为2阶，时滞 d=1%调入数据%参数赋初值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)=theta

12、e_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)plot(1:L,thetae(1:nc,:);%参数 c1,c2 变化曲线xlabel( k ); ylabel( 参数估计 c ); legend( c_2 , c_1 );%参数 b0,b1,b2 变化曲线figure(2)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