系统辨识举例.docx

1、系统辨识举例例3.3 考虑仿真对象 (3.41)其中，是服从正态分布的白噪声N。输入信号采用4阶M序列，幅度为1。选择如下形式的辨识模型 (3.42)设输入信号的取值是从k =1到k =16的M序列，则待辨识参数为=。其中，被辨识参数、观测矩阵z L、H L的表达式为 , ， (3.43)程序框图如图3.2所示。Matlab6.0仿真程序如下：%二阶系统的最小二乘一次完成算法辨识程序，在光盘中的文件名：FLch3LSeg1.mu=-1,1,-1,1,1,1,1,-1,-1,-1,1,-1,-1,1,1; %系统辨识的输入信号为一个周期的M序列z=zeros(1,16); %定义输出观测值的长度

2、for k=3:16 z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2); %用理想输出值作为观测值endsubplot(3,1,1) %画三行一列图形窗口中的第一个图形stem(u) %画输入信号u的径线图形subplot(3,1,2) %画三行一列图形窗口中的第二个图形i=1:1:16; %横坐标范围是1到16，步长为1plot(i,z) %图形的横坐标是采样时刻i, 纵坐标是输出观测值z, 图形格式为连续曲线subplot(3,1,3) %画三行一列图形窗口中的第三个图形stem(z),grid on %画出输出观测值z的径线图形，并显示坐标网格u,

3、z %显示输入信号和输出观测信号%L=14 %数据长度HL=-z(2) -z(1) u(2) u(1);-z(3) -z(2) u(3) u(2);-z(4) -z(3) u(4) u(3);-z(5) -z(4) u(5) u(4);-z(6) -z(5) u(6) u(5);-z(7) -z(6) u(7) u(6);-z(8) -z(7) u(8) u(7);-z(9) -z(8) u(9) u(8);-z(10) -z(9) u(10) u(9);-z(11) -z(10) u(11) u(10);-z(12) -z(11) u(12) u(11);-z(13) -z(12) u(13

4、) u(12);-z(14) -z(13) u(14) u(13);-z(15) -z(14) u(15) u(14) %给样本矩阵HL赋值ZL=z(3);z(4);z(5);z(6);z(7);z(8);z(9);z(10);z(11);z(12);z(13);z(14);z(15);z(16) % 给样本矩阵z L赋值%Calculating Parameters c1=HL*HL; c2=inv(c1); c3=HL*ZL; c=c2*c3 %计算并显示%Display Parameters a1=c(1), a2=c(2), b1=c(3),b2=c(4) %从中分离出并显示a1 、a

5、2、 b1、 b2%End程序运行结果： u = -1，1，-1，1，1，1，1，-1，-1，-1，1，-1，-1，1，1z = 0，0，0.5000，0.2500，0.5250，2.1125， 4.3012，6.4731，6.1988，3.2670，-0.9386， -3.1949，-4.6352，6.2165，-5.5800，-2.5185HL = 0 0 1.0000 -1.0000 -0.5000 0 -1.0000 1.0000 -0.2500 -0.5000 1.0000 -1.0000 -0.5250 -0.2500 1.0000 1.0000 -2.1125 -0.5250 1

6、.0000 1.0000 -4.3012 -2.1125 1.0000 1.0000 -6.4731 -4.3012 -1.0000 1.0000 -6.1988 -6.4731 -1.0000 -1.0000 -3.2670 -6.1988 -1.0000 -1.0000 0.9386 -3.2670 1.0000 -1.0000 3.1949 0.9386 -1.0000 1.0000 4.6352 3.1949 -1.0000 -1.0000 6.2165 4.6352 1.0000 -1.0000 5.5800 6.2165 1.0000 1.0000ZL = 0.5000，0.250

7、0，0.5250，2.1125，4.3012，6.4731，6.1988，3.2670，-0.9386，-3.1949， -4.6352，-6.2165，-5.5800，-2.5185Tc = -1.5000，0.7000，1.0000，0.5000Ta1 = -1.5000a2 = 0.7000b1 = 1.0000b2 =0.5000 从仿真结果表3.1可以看出，由于所用的输出观测值没有任何噪声成分，所以辨识结果也无任何误差。例3.5 考虑图3.6所示的仿真对象，图中， 是服从N分布的不相关随机噪声。且 ， ，选择图3.6所示的辨识模型。仿真对象选择如下的模型结构 (3.68)其中，是服从

8、正态分布的白噪声N。输入信号采用4位移位寄存器产生的M序列，幅度为0.03。按式 (3.69)构造h (k)；加权阵取单位阵；利用式(3.61)计算K(k)、和P(k)，计算各次参数辨识的相对误差，精度满足要求式(3.67)后停机。最小二乘递推算法辨识的Malab6.0程序流程如图3.7所示。下面给出具体程序。%最小二乘递推算法辨识程序, 在光盘中的文件名：FLch3RLSeg3.m%FLch3RLSeg3clear %清理工作间变量L=15; % M序列的周期y1=1;y2=1;y3=1;y4=0; %四个移位寄存器的输出初始值for i=1:L;%开始循环，长度为L x1=xor(y3,y

9、4); %第一个移位寄存器的输入是第三个与第四个移位寄存器的输出的“或” x2=y1; %第二个移位寄存器的输入是第一个移位寄存器的输出 x3=y2; %第三个移位寄存器的输入是第二个移位寄存器的输出 x4=y3; %第四个移位寄存器的输入是第三个移位寄存器的输出 y(i)=y4; %取出第四个移位寄存器的幅值为0和1的输出信号，即M序列 if y(i)0.5,u(i)=-0.03; %如果M序列的值为1, 辨识的输入信号取“-0.03” else u(i)=0.03; %如果M序列的值为0, 辨识的输入信号取“0.03” end %小循环结束 y1=x1;y2=x2;y3=x3;y4=x4;

10、 %为下一次的输入信号做准备end %大循环结束，产生输入信号u figure(1); %第一个图形stem(u),grid on %显示出输入信号径线图并给图形加上网格z(2)=0;z(1)=0; %设z的前两个初始值为零for k=3:15; %循环变量从3到15 z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2); %输出采样信号 end %RLS递推最小二乘辨识c0=0.001 0.001 0.001 0.001; %直接给出被辨识参数的初始值,即一个充分小的实向量p0=106*eye(4,4); %直接给出初始状态P0，即一个充分大的实数单位矩阵

11、E=0.000000005; %取相对误差E=0.000000005c=c0,zeros(4,14); %被辨识参数矩阵的初始值及大小e=zeros(4,15); %相对误差的初始值及大小for k=3:15; %开始求K h1=-z(k-1),-z(k-2),u(k-1),u(k-2); x=h1*p0*h1+1; x1=inv(x); %开始求K(k) k1=p0*h1*x1;%求出K的值 d1=z(k)-h1*c0; c1=c0+k1*d1; %求被辨识参数c e1=c1-c0; %求参数当前值与上一次的值的差值 e2=e1./c0; %求参数的相对变化 e(:,k)=e2; %把当前相

12、对变化的列向量加入误差矩阵的最后一列 c0=c1; %新获得的参数作为下一次递推的旧参数 c(:,k)=c1; %把辨识参数c 列向量加入辨识参数矩阵的最后一列 p1=p0-k1*k1*h1*p0*h1+1; %求出 p(k)的值 p0=p1; %给下次用 if e2c = 0.0010 0 0.0010 -0.4984 -1.2328 -1.4951 -1.4962 -1.4991 -1.4998 -1.4999 0.0010 0 0.0010 0.0010 -0.2350 0.6913 0.6941 0.6990 0.6998 0.6999 0.0010 0 0.2509 1.2497 1

13、.0665 1.0017 1.0020 1.0002 0.9999 0.9998 0.0010 0 -0.2489 0.7500 0.5668 0.5020 0.5016 0.5008 0.5002 0.5002图3.8 最小二乘递推算法的参数辨识仿真-1.5000 -1.5000 -1.5000 -1.4999 -1.49990.7000 0.7000 0.7000 0.7000 -0.70000.9999 0.9999 0.9999 0.9999 0.99990.5000 0.5000 0.5000 0.5000 0.5000e = 0 0 0 -499.4200 1.4734 0.212

14、8 0.0007 0.0020 0.0004 0.0000 0 0 0 0 -235.9916 -3.9416 0.0042 0.0070 0.0012 0.0001 0 0 249.8612 3.9816 -0.1466 -0.0607 0.0003 -0.0018 -0.0003 -0.0001 0 0 -249.8612 -4.0136 -0.2443 -0.1143 -0.0007 -0.0016 -0.0012 -0.00010.0001 0.0000 -0.0000 -0.0000 0.00000.0001 -0.0000 0.0000 0.0000 0.00000.0001 0.

15、0000 0.0000 0.0000 -0.0000-0.0004 0.0000 -0.0000 0.0000 -0.0000表3.2 最小二乘递推算法的辨识结果参 数 a1 a2 b1 b2真 值 -1.5 0.7 1.0 0.5 估计值 -1.4999 0.7 0.9999 0.5000仿真结果表明，大约递推到第十步时，参数辨识的结果基本达到稳定状态，即a1=-1.4999, a2= 0.7000, b1=0.9999, b2=0.5000。此时，参数的相对变化量E0.000000005。从整个辨识过程来看，精度的要求直接影响辨识的速度。虽然最终的精度可以达到很小，但开始阶段的相对误差会很

16、大，从图3.8(3)图形中可看出，参数的最大相对误差会达到三为数第4章第4章例：考虑图4.2表示的线性时不变SISO系统。根据卷积定理，系统的输出y(k)与输入序列的关系可表示成 (4.40)其中，组成系统的脉冲响应。系统在伪随机码输入作用下的输出响应如表4.1所示。利用这些数据辨识系统的脉冲响应，当N取3时，有 (4.41)令 (4.42)若系统的真实脉冲响应记作g0, 则有 (4.43)根据式(4.16)，脉冲响应估计值的递推算式为 (4.44)其中，脉冲响应估计值的初始值取，权矩阵取如下形式 (4.45a) (4.45b)确定性问题梯度校正递推算法辨识的流程如图4.3所示。下面给出Mal

17、ab6.0程序。%梯度校正参数辨识程序，在光盘中的文件名：FLch4GAeg1.m clearu=-1,-1,-1,-1,1,1,1,-1,1,1,-1,-1,1,-1,1,-1,-1,-1,-1,1;y=0,-2,-6,-7,-7,-3,5,7,3,-1,5,3,-5,-3,1,-1,1,-5,-7,-7;%画出u和y图形figure(1), subplot(2,1,1), stem(u), subplot(2,1,2), stem(y), hold onk=1:20；plot(k,y)%给出初始值h1=-1,0,0; h2=-1,-1,0; g=0,0,0; I=1,0,0;0,1/2,0

18、;0,0,1/4;h=h1,h2,zeros(3,16);%计算样本数据h(k)for k=3:18 h(:,k)=u(k),u(k-1),u(k-2);end%计算权矩阵R(k)和g的估计值for k=1:18a=h(1,k)2+(h(2,k)2)/2+(h(3,k)2)/4; a1=1/a; R=a1*I; %按照式(4.45a)和(4.45b)计算权矩阵g(:,k+1)=g(:,k)+R*h(:,k)*(y(k+1)-h(:,k)*g(:,k); %按照式(4.44)计算脉冲响应的估计值end%绘图g1=g(1,:); g2=g(2,:); g3=g(3,:);figure(2); k=

19、1:19; subplot(121); plot(k,g1,r,k,g2,g,k,g3,b), grid on%计算模型输出ym及系统输出与模型输出之间的误差Eyfor k=1:18ym(k)=h(:,k)*g(:,k); Ey(k)=y(k+1)-ym(k);endk=1:18; subplot(122); plot(k,Ey), grid ong, ym, Ey %显示脉冲响应估计值、模型输出及系统输出与模型输出之间的误差figure(3); x=0:1:3; y=0,g(1,18),g(2,18),g(3,18); xi=linspace(0,3); yi=interp1(x,y,xi,

20、cubic);plot(x,y,o,xi,yi,m), grid on %画出脉冲响应估计值及其三次插值曲线 程序执行结果：g =0 2.0000 4.6667 5.2381 5.2381 1.5374 2.1438 2.2393 1.9658 1.95590 0 1.3333 1.6190 1.6190 3.4694 3.7726 3.8204 3.9571 3.96210 0 0 0.1429 0.1429 1.0680 0.9164 0.9403 1.0087 1.00622.0063 1.9918 1.9980 1.9953 1.9995 1.9995 1.9995 2.0001 2.

21、00323.9873 3.9946 3.9977 3.9990 3.9969 3.9969 3.9969 3.9972 3.99880.9936 0.9972 0.9957 0.9963 0.9974 0.9974 0.9974 0.9972 0.9980ym = 0,-2.0000,-6.0000,-7.0000,3.4762,3.9388,6.8328,2.5213,-0.9826,4.9118,2.9746,-4.9891,-2.9953, 1.0073, -1.0000,1.0000,-4.9990,-6.9945Ey =-2.0000,-4.0000,-1.0000,0.0000,-

22、6.4762,1.0612,0.1672,0.4787,-0.0174,0.0882,0.0254,0.0109, -0.0047,-0.0073,-0.0000,0,-0.0010,-0.0055 图4.4 确定性问题梯度校正参数辨识仿真(被辨识参数的个数为3) 图4.5 确定性问题梯度校正参数辨识仿真(被辨识参数的个数为5)图4.4表明，递推到第10步时，被辨识参数基本上达到了稳定状态，即, ; 此时系统的输出与模型的输出误差也基本达到稳定状态，即。由于被辨识参数的个数较少，递推校正算法的收敛性比较好，也就是说，输入输出的观测数据量已足够。但从图4.4 Figure 3中可看出，由于被辨识

23、参数的个数较少，它还不足以充分显示全部的系统脉冲响应。为了充分地显示出系统的脉冲响应，可以把被辨识参数的个数增加到5个。图4.5给出了被辨识参数的个数为5时的辨识结果。图4.5基本显示出了系统的全部脉冲响应过程。但在图4.5的辨识中，还是利用上面给出的20对输入输出数据。通过图4.4与图4.5的比较可以看出，在观测数据量一定的情况下，随着被辨识参数的个数的增加，梯度校正辨识算法的收敛性变差，这是由于观测数据不足，递推步数有限造成的。所以，随着被辨识参数的个数的增加，若要保持梯度校正辨识算法的收敛效果，需要同步增加观测数据量。第七章的用改进的神经网络MBP算法辨识 例7.1 对具有随机噪声的二阶

24、系统的模型辨识对具有随机噪声的二阶系统的模型辨识，进行标幺化以后系统的参考模型差分方程为 (7.90)式中，为随机噪声。由于神经网络的输出最大为1，所以，被辨识的系统应先标幺化，这里标幺化系数为5。利用图7.5正向建模(并联辨识)结构，神经网络选用3-9-9-1型，即输入层i，隐层j包括2级，输出层k的节点个数分别为3、9、9、1个；由于神经网络的最大输出为1，因此在辨识前应对原系统参考模型标么化处理，辨识结束后再乘以标么化系数才是被辨识系统的辨识结果。1)编程如下:%w10ij表示第一隐层权值，w11ij表示；w120ij表示第二隐层权值%， w121ij表示；w20j表示输出层权值，w21

25、j表示%；q表示隐层阈值；p表示输出层阈值；置标幺化系数f1=5等w10ij=.01 .01 .02; .1 .11 .02; .01 0 .1; .11 .01 .02;.1 .1 .02; .11 .1 .1;.1 .1 .1;0 .1 .1;.1 0 .1; w11ij=.1 .2 .11; .02 .13 .04; .09 .08 .08; .09 .1 .06; .1 .11 .02; .06 0 .1;.1 .1 .1;0 .1 0;.1 .1 .1; w20j=.01;.02;.1;.2;.1;.1;.1;.1;.1; w21j=0; 0.1; .1; .02;0;.1;.1;.

26、1;.1;q0j=.9 .8 .7 .6 .1 .2 .1 .1 .1; q120j=q0j;q11j=.5 .2 .3 .4 .1 .2 .1 .1 .1;q12j=q11j;w121ij=w20j*q0j;w120ij=w20j*q11j; f1=5;q2j=0; % threshold valuep0=.2;k1=1;p1=.3;w=0;xj=1 1 1; % inputs error=0.0001;a1=1 1 1 1;n=1;e1=0;e0=0;e2=0;e3=0;e4=0;yo=0;ya=0;yb=0;y0=0;y1=0;y2=0;y3=0;u=0;u1=0;u2=0.68;u3=

27、.780;u4=u3-u2;k1=1;kn=28;e3=.055; z1=0;z12=0; q123j=0; t2j=0; o12j=0;r=0;r1=0; s=0.1;d2j=0;%+% calculating output of the hidden layer v1=randn(40,1);for m=1:40 s1=0.1*v1(m) yn=.3366*y2+.6634*u1+s*s1;y1=y2;y2=yn; yp=yn;u0=u1;u1=u2; ya(m)=yn; for k=1:100% calculating output of the hidden layer(1) for i=1:9 x1=w11ij(i,1)*xj(:,1)+w11ij(i,2)*xj(:,2)+w11ij(i,3)*xj(:,3); x=x1+q11j(:,i); o=1/1+exp(-x); o11j(i)=o; end% calculating output of the hidden layer(2) for i=1:9 for j=1:9 z1=z1+w121ij(i,j)*o11j(:,j); end z=z1+q12j(:,i); o=