DOC《实验报告》

实验报告课程名称:系统辨识 姓名: 谈敏刚 学号:

112101123 课程号:34

一、实验目的 通过实验了解辨识方法在工程应用中的一些实际问题;

了解数据获取和数据 处理的各种方法和手段,掌握各种辨识方法的应用特点.

二、实验内容

1 数据获取 在机械制造工艺中经常采用的加热设备有多种,如热加工中所用的锻件加热炉,热处理中所用的各种热处理电炉以及金属材料高温机械性能试验用的加热炉等.这类设备的输入量是燃料供给量或电压、电流,而输出量是炉膛内腔的温度.为了采用计算机进行精确的温度控制,如按一定的规律升温或保持恒温等,必须先掌握加热设备的动态特性. 这里以高温老化试验温箱为例,以控制电压作为炉温控制系统的输入控制变量.在热稳定工况的基础上,在电压稳定值上再附加一个辨识信号,即M序列电压信号.加热炉热惯性大,升温过程较长,所以采样周期较长,M序列的周期也较长.这样施加M序列周期信号之后,记录几个周期的炉温试验数据.

2 数据处理 为了提高辨识精度,实验者必须对原始数据进行剔除坏数据、零均值化、工 频滤波等处理.

3 离线辨识 利用处理过的数据,选择某种辨识方法;

如RLS、RELS、RIV或RML等参数估计算法,以及F检验法或AIC定阶法.离线估计出模型参数和阶次,并 计算相应的模型静态增益,同时比较利用不同方法所得到的辨识结果.

三、实验步骤

1 确定实验方案,建立数据文件,编制实验程序,上机调试,记录实验结果;

实验具体步骤按指导书要求进行,这里不再详述. 设置的主要参数为:电压为45V,采样周期为2S,M序列长度为15,为80S.

2 数学模型结构辨识. 2.1根据汉格尔矩阵估计模型的阶次 设一个可观可控的SISO过程的脉冲响应序列为{个g(1),g(2),……g(L)},可以通过汉格尔(Hankel)矩阵的秩来确定系统的阶次. 令Hankel阵为: ,其中决定阵地维数,k可在1至间任意选择.则有. 如果(过程的真实阶次),那么Hankel阵的秩等于.因此可以利用Hankel阵的奇异性来确定系统的阶次. 2.2根据残差平方和估计模型的阶次 SISO过程的差分方程模型的输出残差为,数据长度L,为阶时的数据矩阵,为阶时的参数的估计量,为模型阶次估计值,为真实阶次,则残差平方和函数: 残差平方和有这样的性质:当L足够大时,随着增加先是显著地下降,当>

时,值显著下降的现象就终止.这就是损失函数法来定阶的原理. 4.3根据AIC准则估计模型的阶次 具体的定阶用法是:对不同阶次首先用极大似然法估计参数,然后计算似然函数值及值,找到使的作为.

3 实验结果 取模型为y(k)+a1*y(k-1)an*y(k-n)=b1*u(k-1-d)bn*u(k-n-d) 无时滞 a1 a2 a3 a4 a5 b1 b2 b3 b4 b5 AIC n=1 -1.0336 0.0165 -56.5880 n=2 -0.4448 -0.6208 0.0225 0.0166 -64.9677 n=3 -0.4761 -0.5468 -0.0464 0.0183 0.0210 0.0043 -65.4197 n=4 -0.4549 -0.6326 -0.1026 0.1305 0.0185 0.0210 0.0026 0.0096 63.4235 n=5 -0.4278 -0.7071 -0.1865 0.0711 0.2097 0.0210 0.0212 0.0017 0.0035 -0.0037 -61.0714 阶数为3阶 时滞d=1,结果为 a1 a2 a3 a4 a5 b1 b2 b3 b4 b5 AIC n=1 -1.0320 0.0040 -56.3187 n=2 -0.5566 -0.5059 0.0196 -0.0001 -66.8562 n=3 -0.5170 -0.5315 -0.0162 0.0209 0.0014 0.0110 -64.0782 n=4 -0.5461 -0.6187 -0.0446 0.1603 0.0186 -0.0025 0.0071 -0.0025 -60.7751 n=5 -0.5101 -0.6873 -0.2906 0.2481 0.2193 0.0220 -0.0032 0.0047 -0.0119 -0.0041 -61.7805 阶次为2阶 时滞d=2,结果为 a1 a2 a3 a4 a5 b1 b2 b3 b4 b5 AIC n=1 -1.0356 0.0002 -61.4689 n=2 -0.5779 -0.4812 0.0022 0.0081 -64.9881 n=3 -0.6503 -0.5136 0.1135 -0.0010 0.0082 -0.0053 -61.4695 n=4 -0.6713 -0.6737 -0.0339 0.3592 -0.0059 0.0047 -0.0107 -0.0065 -61.1731 n=5 -0.6395 -0.7426 -0.2224 0.3505 0.2594 -0.0117 -0.0042 -0.0180 -0.0080 -0.0124 -59.7794 阶次为2阶 时滞d=3,结果为 a1 a2 a3 a4 a5 b1 b2 b3 b4 b5 AIC n=1 -1.0342 0.0100 -62.3075 n=2 -0.5986 -0.4597 0.0088 -0.0050 -65.1569 n=3 -0.6561 -0.4957 0.1025 0.0104 -0.0076 -0.0033 -62.3620 n=4 -0.6496 -0.6522 -0.0880 0.3662 0.0028 -0.0118 -0.0040 -0.0081 -61.6653 n=5 -0.5639 -0.6280 -0.3812 0.1777 0.3989 -0.0035 -0.0229 -0.0082 -0.0108 -0.0081 -87.7509 阶次为5阶 综上所述: 辨识的结果为,系统为2阶,时滞d=1的时候,系统性能最优. y(k)= -0.5566y(k-1)-0.5059y(k-2)+0.0196u(k-2)-0.0001u(k-3) 模型验证: 多组数据交叉检验,结果为: 附程序: clc clear %清理工作间变量 L=30;

M序列的长度 y1=1;

y2=1;

y3=1;

y4=0;

%四个移位积存器的输出初始值 for i=1:L;

开始循环,长度为L x1=xor(y3,y4)第一个移位积存器的输入是第3个与第4个移位积存器的输出的 或 x2=y1;

第二个移位积存器的输入是第3个移位积存器的输出 x3=y2;

第三个移位积存器的输入是第2个移位积存器的输出 x4=y3;

第四个移位积存器的输入是第3个移位积存器的输出 y(i)=y4;

取出第四个移位积存器幅值为

0 和

1 的输出信号, if y(i)>

0.5,u(i)=-5;

%如果M序列的值为

1 时,辨识的输入信号取 -0.03 else u(i)=5;

当M序列的值为

0 时,辨识的输入信号取 0.03 end %小循环结束 y1=x1;

y2=x2;

y3=x3;

y4=x4;

为下一次的输入信号做准备 end %大循环结束,产生输入信号u z=zeros(1,30)存放加入M序列之后的温度 file=csvread('

2013-4-27-10-21-13DataDemo_Sequence_Identy.csv'

,1,3);

%导入温度数据 j=1;

for i=1:1200 if(mod(i,40)~=0)M序列的长度为15,周期为80*15,采样时间为2s,所以在一个 DA=file(i,1)-111;

状态中采样了40次,把这四十个数据减去稳态值,再取均值 z(j)=z(j)+DA;

else z(j)=z(j)/40;

j=j+1;

end end %系统阶次辨识AIC算法% bb=zeros(5,1);

AIC=zeros(4,5)存放不同阶次,不同时滞下的AIC, for d=0:3 %设时滞d从1~3,阶次1~5 H1=zeros(29-d,2);

H2=zeros(28-d,4);

H3=zeros(27-d,6);

H4=zeros(26-d,8);

H5=zeros(25-d,10);

zz1=zeros(29-d,1);

zz2=zeros(28-d,1);

zz3=zeros(27-d,1);

zz4=zeros(26-d,1);

zz5=zeros(25-d,1);

n=1;

假设为1阶for i=2+d:30 H1(i-1-d,:)=[-z(i-1) u(i-1-d)];

zz1(i-1-d,:)=[z(i)];

end aa1=inv(H1'

*H1)*H1'

*zz1 bb(1)=(zz1-H1*aa1)'

*(zz1-H1*aa1)/L;

AIC(d+1,1)=L*log(bb(1))+4*n;

n=2;

假设为2阶for i=3+d:30 H2(i-2-d,:)=[-z(i-1) -z(i-2) u(i-1-d) u(i-2-d)];

zz2(i-2-d,:)=[z(i)];

end aa2=inv(H2'

*H2)*H2'

*zz2 bb(2)=(zz2-H2*aa2)'

*(zz2-H2*aa2)/L;

AIC(d+1,2)=L*log(bb(2))+4*n;

n=3;

假设为3阶for i=4+d:30 H3(i-3-d,:)=[-z(i-1) -z(i-2) -z(i-3) u(i-1-d) u(i-2-d) u(i-3-d)];

zz3(i-3-d,:)=[z(i)];

end aa3=inv(H3'

*H3)*H3'

*zz3 bb(3)=(zz3-H3*aa3)'

*(zz3-H3*aa3)/L;

AIC(d+1,3)=L*log(bb(3))+4*n;

n=4;

假设为4阶for i=5+d:30 H4(i-4-d,:)=[-z(i-1) -z(i-2) -z(i-3) -z(i-4) u(i-1-d) u(i-2-d) u(i-3-d) u(i-4-d)];

zz4(i-4-d,:)=[z(i)];

end aa4=inv(H4'

*H4)*H4'

*zz4 bb(4)=(zz4-H4*aa4)'

*(zz4-H4*aa4)/L;

AIC(d+1,4)=L*log(bb(4))+4*n;

n=5;

假设为5阶for i=6+d:30 H5(i-5-d,:)=[-z(i-1) -z(i-2) -z(i-3) -z(i-4) -z(i-5) u(i-1-d) u(i-2-d) u(i-3-d) u(i-4-d) u(i-5-d)];

zz5(i-5-d,:)=[z(i)];

end aa5=inv(H5'

*H5)*H5'

*zz5 bb(5)=(zz5-H5*aa5)'

*(zz5-H5*aa5)/L;

AIC(d+1,5)=L*log(bb(5))+4*n;

x=min(AIC(d+1,:)) for i=1:5 if(AIC(d+1,i)==x) N=i %取同一时滞下,AIC最小的,就是所辨识出的阶次N end end end figure, plot(AIC(1,:));

grid on title('

谈敏刚,AIC定阶法,时滞d=0'

) figure, plot(AIC(2,:));

grid on title('

谈敏刚,AIC定阶法........