系统辨识实验报告xjtu

/91系统辨识实验报告高海南硕037一、实验问题单输入单输出系统：111()()()()()()AzytBzutCzet其中：11231123112()12.8512.7170.865()()10.70.22AzzzzBzzzzCzzz()~(0,1)etN1．产生必要随机数；2．设计实验，产生输出数据；3．用最小二乘辨识系统模型；4．任用一种适合有色噪声辨识算法辨识系统模型。二、实验原理1.递推最小二乘递推最小二乘辨识算法可减少运算量和数据在计算机中所占的存储量，同时也能实时辨识出系统的特性，将最小二乘转化为递推估计。最小二乘递推算法RLS的基本思想是ˆˆ1kk新的估计值（）＝老的估计值（－）＋修正项算法迭代方程为1ˆˆˆ()(1)()[()()(1)]()(1)()[()(1)()1]()(1)()()(1)TTTkkKkzkhkkKkPkhkhkPkhkPkPkKkhkPk（1）2.增广最小二乘增广最小二乘递推算法，扩充了最小二乘法的参数向量和数据向量)(kh的维数，把噪声模型的辨识同时考虑进去，因此被称为增广最小二乘法。最小二乘法的许多结论对它都是适用的，但最小二乘法只能获得模型的参数估计。如果噪声模型必须用)()(1kvzD表示时，只能用RELS算法，方可获得无偏估计，这是RLS算法所不能代替的。考虑SISO的动态系统，输入)(ku和输出)(kz是可以观测的；)z(G1是系统/92模型，用来描述系统的输入输出特性，)(ky是系统的实际输出。)z(N1是噪声模型，)(kv是均值为零的不相关随机噪声。通常)z(A)z(B)z(G1-11,)z(C)z(D)z(N111(2)式中ddccbbaann22111nn22111nn22111-nn22111zzz1)z(zzz1)z(zzz)z(zzz1)z(dddDcccCbbbBaaaA(3)若SISO系统采用平均滑动模型，即)z(A)z(C-1-1(4))()(D)()z(B)(z)z(A111kvzkuk若假定模型阶次an、bn、dn已经确定，则这类问题的辨识可用增广最小二乘法，以便获得满意的结果。令TdbaTnnnnkvkvnkukunkzkzkhdddbbbaaadba)](),1(),(,),1(),(,),1([)(],,,,,,,,,,,[212121(5)将模型(3)化为最小二乘格式)()()(zTkvkhk(6)由于)(kv是白噪声，所以利用最小二乘法即可获得参数的无偏估计。但是数据向量)(kh中包含着不可测的噪声量)(,),1(dnkvkv,它可用相应的估计值替代。置Tdbankvkvnkukunkzkzkh)](),1(),(,),1(),(,),1([)((7)式中，;0,0)(ˆkkv当k0时，)1(ˆ)()()(ˆkkhkzkvT(8)或)(ˆ)()()(ˆkkhkzkvT(9)增广最小二乘递推算法RELS/931ˆˆˆ()(1)()[()()(1)]()(1)()[()(1)()1()]()(1)()()[()(1)()]TTTTkkKkzkhkkKkPkhkhkPkhkkPkPkKkKkhkPkhkI(10)如果11，即所有采样数据都是等同加权时，增广最小二乘递推算法RELS可以写为1ˆˆˆ()(1)()[()()(1)]()(1)()[()(1)()1]()[()()](1)TTTkkKkzkhkkKkPkhkhkPkhkPkIKkhkPk(11)三、实验过程1.产生白噪声和m序列白噪声序列()et可由乘同余法自己编写函数实现，也可以由matlab自带的randn函数产生，本实验采用randn函数产生；系统的输入()ut选择m序列，其级数为4，初相为[1110],本原多项式为43f()1xxx。2.产生输出采样信号考虑噪声对输出的影响，则系统输出为2.85112.71720.86531230.710.222zkzkzkzkukukukvkvkvk3.递推最小二乘辨识给被辨识参数和P赋初值。根据（1）式计算()Kk、()k、()Pk，当参数收敛满足要求时停止迭代计算，输出辨识结果。4.增广最小二乘辨识给被辨识参数和P赋初值。根据（11）式计算()Kk、()k、()Pk，当参数收敛满足要求时停止迭代计算，输出辨识结果。四．实验结果与分析1.系统白噪声和系统输入m序列如图1所示。/94图1系统噪声及m序列2.采用递推最小二乘辨识算法实验时，给定辨识精度为E＝0.00005，m序列信号峰峰大小为10，经过157次迭代，达到给定经度。此时辨识结果为a1=-2.8599，a2=2.7339，a3=-0.8731b1=0.9875，b2=0.9681，b3=0.9705在实验时，对于m序列，其峰峰值越大，辨识的参数与真值误差越小。递推最小二乘的仿真结果如图2和图3所示。图2递推最小二乘参数结果/95图3参数误差收敛情况3.增广最小二乘辨识实验时，给定辨识精度为E＝0.00005，m序列峰峰值为2，经过40次迭代达到给定经度，此时辨识结果为a1=-2.851,a2=2.7217,a3=-0.865；b1=1，b2=1，b3=1；c1=1，c2=0.7，c3=0.22增广最小二乘的仿真结果如图4和图5所示。图4增广最小二乘参数辨识结果图5参数误差收敛情况参数的真值和估计值如下表所示：/96参数真值递推最小二乘辨识增广最小二乘辨识a1-2.851-2.8599-2.8510a22.7172.73392.7170a3-0.865-0.8731-0.8650b11.00.98751.0000b21.00.96811.0000b31.00.97051.0000c11.0-1.0000c20.7-0.7000c30.22-0.2200从上面两表对比可以看到，增广最小二乘考虑了噪声模型，与递推最小二乘算法相比，速度快、辨识结果精确，而且可以得到噪声模型参数。四．附录附录1递推最小二乘辨识matlab仿真程序clearL=200;y1=1;y2=1;y3=1;y4=0;fori=1:Lx1=xor(y3,y4);x2=y1;x3=y2;x4=y3;y(i)=y4;ify(i)0.5u(i)=-5;elseu(i)=5;endy1=x1;y2=x2;y3=x3;y4=x4;endfigure(1)i=1:20;stem(u(i));gridon;title('ÊäÈëMÐòÁеÄÇ°20¸öÖÜÆÚ')v=1*randn(L,1);z(3)=0;z(2)=0;z(1)=0;fork=4:L;z(k)=2.851*z(k-1)-2.717*z(k-2)+0.865*z(k-3)+u(k-1)+u(k-2)+u(k-3)+v(k)+0.7*v(k-1)+0.22*v(k-2);endc0=[0.0010.0010.0010.0010.0010.001]';p0=10^3*eye(6,6);E=0.000005;/97c=[c0,zeros(6,L-1)];e=zeros(6,L);fork=4:Lh1=[-z(k-1),-z(k-2),-z(k-3),u(k-1),u(k-2),u(k-3)]';k1=p0*h1*inv(h1'*p0*h1+1);c1=c0+k1*(z(k)-h1'*c0);e1=c1-c0;e2=e1./c0;e(:,k)=e2;c0=c1;c(:,k)=c1;p1(:,:,k)=p0-k1*k1'*[h1'*p0*h1+1];p0=p1(:,:,k);ifabs(e2)=EN=k;break;endN=k;enda1=c(1,1:N);a2=c(2,1:N);a3=c(3,1:N);b1=c(4,1:N);b2=c(5,1:N);b3=c(6,1:N);ea1=e(1,:);ea2=e(2,:);ea3=e(3,:);eb1=e(4,:);eb2=e(5,:);eb3=e(6,:);i=1:N;figure(2)plot(i,a1,'r',i,a2,'r:',i,a3,'r+',i,b1,'g',i,b2,'g:',i,b3,'g+');gridon;title('µÝÍÆ×îС¶þ³Ë±æʶ')figure(3)i=1:30;plot(i,ea1(1:30),'r',i,ea2(1:30),'g',i,eb1(1:30),'b',i,eb2(1:30),'r:');gridon;legend('a1','a2','a3','b1','b2','b3');title('²ÎÊýÊÕÁ²Çé¿ö')附录2增广最小二乘的matlab仿真程序clearall;clc;L=100;y1=1;y2=1;y3=1;y4=0;fori=1:Lx1=xor(y3,y4);x2=y1;x3=y2;x4=y3;y(i)=y4;ify(i)0.5u(i)=-1;elseu(i)=1;/98endy1=x1;y2=x2;y3=x3;y4=x4;endv=1*randn(1,L);z(3)=0;z(2)=0;z(1)=0;fork=4:L;z(k)=2.851*z(k-1)-2.717*z(k-2)+0.865*z(k-3)+u(k-1)+u(k-2)+u(k-3)+v(k)+0.7*v(k-1)+0.22*v(k-2);endfigure(1);subplot(2,1,2);stem(u),gridon;title('ÊäÈëMÐòÁÐ');subplot(2,1,1);plot(v),gridon;title('ÔëÉùÐòÁÐ');c0=[0.0010.0010.0010.0010.0010.0010.0010.0010.001]';p0=10^6*eye(9,9);E=5e-10;c=[c0,zeros(9,L-1)];e=zeros(9,L);N=L+1;fork=4:L;h1=[-z(k-1),-z(k-2),z(k-3),u(k-1),u(k-2),u(k-3),v(k),v(k-1),v(k-2)]';k1=p0*h1*inv(h1'*p0*h1+1);c1=c0+k1*(z(k)-h1'*c0);e1=c1-c0;e2=e1./c0;e(:,k)=e2;c0=c1;c(:,k)=c1;p1(:,:,k)=p0-k1*k1'*[h1'*p0*h1+1];p0=p1(:,:,k);ifabs(e2)=EN=k;break;endendifN~=L+1a1=c(1,1:N);a2=c(2,1:N);a3=c(3,1:N);b1=c(4,1:N);b2=c(5,1:N);b3=c(6,1:N);d1=c(7,1:N);d2=c(8,1:N);d3=c(9,1:N);ea1=e(1,1:N);ea2=e(2,1:N);ea3=e(3,1:N);eb1=e(4,1:N);eb2=e(5,1:N);eb3=e(6,1:N);ed1=e(7,1:N);ed2=e(8,1:N);ed3=e(9,1:N);figure(2);i=1:N;/99plot(i,a1,'r',i,a2,'r:',i,a3,'r+',i,b1,'b',i,b2,'b:',i,b3,'b+',i,d1,'g',i,d2,'g:',i,d3,'g+');gridon;title('Ôö¹ã×îС¶þ³ËËã·¨²ÎÊý¹À¼Æ½á¹û');figure(3);plo