系统辨识及其matlab仿真(一些噪声和辨识算法)

1】随机序列产生程序2】白噪声产生程序【3】M序列产生程序4】二阶系统一次性完成最小二乘辨识程序5】实际压力系统的最小二乘辨识程序6】递推的最小二乘辨识程序7】增广的最小二乘辨识程序8】梯度校正的最小二乘辨识程序9】递推的极大似然辨识程序10】Bayes辨识程序11】改进的神经网络MBP算法对噪声系统辨识程序12】多维非线性函数辨识程序的Matlab程序【13】模糊神经网络解耦Matlab程序【14】F-检验法部分程序【1】随机序列产生程序A=6;x0=1;M=255;fork=1:100x2=A*x0;x1=mod(x2,M);v1=x1/256;v(:,k)=v1;x0=x1;v0=v1;endv2=vk1=k;%grapherk=1:k1;plot(k,v,k,v,'r');xlabel('k'),ylabel('v');title('(O,l)均匀分布的随机序列')【2】白噪声产生程序A=6;x0=1;M=255;f=2;N=100;fork=1:Nx2=A*x0;x1=mod(x2,M);v1=x1/256;v(:,k)=(v1-0.5)*f;x0=x1;v0=v1;endv2=vk1=k;%grapherk=1:k1;plot(k,v,k,v,'r');xlabel('k'),ylabel('v');title('(-l,+l)均匀分布的白噪声')【3】M序列产生程序Xl=l;X2=0;X3=l;X4=0;%移位寄存器输入Xi初T态(0101),Yi为移位寄存器各级输出m=60;%置M序列总长度fori=1:m%1#Y4=X4;Y3=X3;Y2=X2;Y1=X1;X4=Y3;X3=Y2;X2=Y1;X1=xor(Y3,Y4);%异或运算ifY4=0U(i)=-1;elseU(i)=Y4;endendM=U%绘图i1=ik=1:1:i1;plot(k,U,k,U,'rx')xlabel('k')ylabel('M序列')title('移位寄存器产生的M序列')【4】二阶系统一次性完成最小二乘辨识程序%FLch3LSeg1u=-1,1,-1,1,1,1,1,-1,-1,-1,1,-1,-1,1,1;%系统辨识的输入信号为一个周期的M序列z=zeros(1,16);%定义输出观测值的长度fork=3:16z(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),gridon%画出输出观测值z的经线图形，并显示坐标网格u,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)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)%给样本矩阵zL赋值%calculatingparameters%计算参数c1=HL'*HL;c2=inv(c1);c3=HL'*ZL;c=c2*c3%计算并显示%DISPLAYPARAMETERSa1=c(1),a2=c(2),bl=c(3),b2=c(4)%从中分离出并显示al、a2、bl、b2%End【5】实际压力系统的最小二乘辨识程序%FLch3LSeg2clear%工作间清零V=54.3,61.8,72.4,88.7,118.6,194.0'%赋初值V,并显示P=61.2,49.5,37.6,28.4,19.2,10.1'%赋初值P,并显示%logP=-alpha*logV+logbeita=-logV,1alpha,log(beita)'=HL*sita%注释P、V之间的关系fori=1:6;%循环变量的取值为从1到6Z(i)=log(P(i);%赋系统的输出采样值end%循环结束ZL=Z'%给zL赋值HL=-log(V(1),1;-log(V(2),1;-log(V(3),1;-log(V(4),1;-log(V(5),1;-log(V(6),1%给HL赋值%calculatingparameters%计算参数c1=HL'*HL;c2=inv(c1);c3=HL'*ZL;c4=c2*c3%计算%SeparationofParameters%分离变量alpha=c4(1)%为c4的第1个元素beita=exp(c4(2)%为以自然数为底的c4的第2个元素的指数【6】递推的最小二乘辨识程序%FLch3RLSeg3clear%清理工作间变量L=15;%M序列的周期y1=1;y2=1;y3=1;y4=0;%四个移位积存器的输出初始值fori=1"%开始循环，长度为Lx1=xor(y3,y4);%第一个移位积存器的输入是第3个与第4个移位积存器的输出的或x2=y1;%第二个移位积存器的输入是第3个移位积存器的输出x3=y2;%第三个移位积存器的输入是第2个移位积存器的输出x4=y3;%第四个移位积存器的输入是第3个移位积存器的输出y(i)=y4;%取出第四个移位积存器幅值为"0"和"1"的输出信号，ify(i)>0.5,u(i)=-0.03;%如果M序列的值为"1"时，辨识的输入信号取-0.03elseu(i)=0.03;%当M序列的值为"0"时，辨识的输入信号取0.03end%小循环结束y1=x1;y2=x2;y3=x3;y4=x4;%为下一次的输入信号做准备end%大循环结束，产生输入信号ufigure(1);%第1个图形stem(u),gridon%以径的形式显示出输入信号并给图形加上网格z(2)=0;z(1)=0;%取z的前两个初始值为零fork=3:15;%循环变量从3到15z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2);%给出理想的辨识输出采样信号end%RLS递推最小二乘辨识c0=0.0010.0010.0010.001'%直接给出被辨识参数的初始值,即一个充分小的实向量p0=10人6*eye(4,4);%直接给出初始状态P0，即一个充分大的实数单位矩阵E=0.000000005;%相对误差E=0.000000005c=c0,zeros(4,14);%被辨识参数矩阵的初始值及大小e=zeros(4,15);%相对误差的初始值及大小fork=3:15;%开始求Kh1=-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;%求被辨识参数ce1=c1-c0;%求参数当前值与上一次的值的差值e2=e1./c0;%求参数的相对变化e(:,k)=e2;%把当前相对变化的列向量加入误差矩阵的最后一列c0=c1;%新获得的参数作为下一次递推的旧参数c(:,k)=c1;%把辨识参数c列向量加入辨识参数矩阵的最后一列p1=p0-k1*k1'*h1'*p0*h1+1;%求出p(k)的值p0=p1;%给下次用ife2<=Ebreak;%若参数收敛满足要求，终止计算end%小循环结束end%大循环结束c%显示被辨识参数e%显示辨识结果的收敛情况%分离参数a1=c(1,;a2=c(2,;b1=c(3,;b2=c(4,;ea1=e(1,;ea2=e(2,;eb1=e(3,;eb2=e(4,;figure(2);%第2个图形i=1:15;%横坐标从1到15plot(i,a1,'r',i,a2,':',i,b1,'g',i,b2,':')%画出a1，a2，b1，b2的各次辨识结果title('ParameterIdentificationwithRecursiveLeastSquaresMethod')%图形标题figure(3);%第3个图形i=1:15;%横坐标从1到15plot(i,ea1,'r',i,ea2,'g',i,eb1,'b',i,eb2,'r:')%画出a1，a2，b1，b2的各次辨识结果的收敛情况title('IdentificationPrecision')%图形标题【7】增广的最小二乘辨识程序%FLch3ELSeg4clearL=60;%四位移位积存器产生的M序列的周期yl=l;y2=l;y3=l;y4=0;%四个移位积存器的输出初始值fori=l7;xl=xor(y3,y4);%第一个移位积存器的输入信号x2=y1;%第二个移位积存器的输入信号x3=y2;%第三个移位积存器的输入信号x4=y3;%第四个移位积存器的输入信号y(i)=y4;%第四个移位积存器的输出信号，幅值"0"和"1"ify(i)>0.5,u(i)=-1;%M序列的值为"1"时，辨识的输入信号取-1elseu(i)=1;%M序列的值为"0"时，辨识的输入信号取1endy1=x1;y2=x2;y3=x3;y4=x4;%为下一次的输入信号作准备endfigure(1);%画第一个图形subplot(2,1,1);%画第一个图形的第一个子图stem(u),gridon%画出M序列输入信号v=randn(1,60);%产生一组60个正态分布的随机噪声subplot(2,1,2);%画第一个图形的第二个子图plot(v),gridon;%画出随机噪声信号R=corrcoef(u,v);%计算输入信号与随机噪声信号的相关系数r=R(1,2);%取出互相关系数u%显示输入型号V%显示噪声型号z=zeros(7,60);zs=zeros(7,60);zm=zeros(7,60);zmd=zeros(7,60);%输出采样、不考虑噪声时系统输出、不考虑噪声时模型输出、模型输出矩阵的大小z(2)=0;z(1)=0;zs(2)=0;zs(1)=0;zm(2)=0;zm(1)=0;zmd(2)=0;zmd(1)=0;%给输出采样、不考虑噪声时系统输出、不考虑噪声时模型输出、模型输出赋初值%增广递推最小二乘辨识c0=0.0010.0010.0010.0010.0010.0010.001'%直接给出被辨识参数的初始值,即一个充分小的实向量p0=10人6*eye(7,7);%直接给出初始状态P0，即一个充分大的实数单位矩阵E=5.0e-15;%取相对误差Ec=c0,zeros(7,59);%被辨识参数矩阵的初始值及大小e=zeros(7,60);%相对误差的初始值及大小fork=3:60;%开始求Kz(k)=1.5*z(k-l)-0.7*z(k-2)+u(k-l)+0.5*u(k-2)+v(k)-v(k-l)+0.2*v(k-2);%系统在M序列输入下的输出采样信号hl=-z(k-l),-z(k-2),u(k-l),u(k-2),v(k),v(k-l),v(k-2)'%为求K(k)作准备x=h1'*p0*h1+1;x1=inv(x);k1=p0*h1*x1;%Kd1=z(k)-h1'*c0;cl=cO+kl*dl;%辨识参数czs(k)=1.5*z(k-l)-0.7*z(k-2)+u(k-l)+0.5*u(k-2);%系统在M序列的输入下不考虑扰动时的输出响应zm(k)=-z(k-1),-z(k-2),u(k-1),u(k-2)*cl(l);c1(2);cl;cl(4);%模型在M序列的输入下不考虑扰动时的的输出响应zmd(k)=h1'*cl;%模型在M序列的输入下的的输出响应e1=c1-c0;e2=e1./c0;%求参数的相对变化e(:,k)=e2;cO=cl;%给下一次用c(:,k)=cl;%把辨识参数c列向量加入辨识参数矩阵p1=p0-k1*k1'*h1'*p0*h1+1;%findp(k)pO=pl;%给下次用ife2<=Ebreak;%若收敛情况满足要求，终止计算end%判断结束end%循环结束c,e,%显示被辨识参数及参数收敛情况z,zmd%显示输出采样值、系统实际输出值、模型输出值%分离变量al=c(l,;a2=c(2,;bl=c(3,;b2=c(4,;%分离出al、a2、bl、b2d1=c(5,;d2=c(6,;d3=c(7,;%分离出d1、d2、d3eal=e(l,;ea2=e(2,;ebl=e(3,;eb2=e(4,;%分离出al、a2、bl、b2的收敛情况ed1=e(5,;ed2=e(6,;ed3=e(7,;%分离出d1、d2、d3的收敛情况figure(2);%画第二个图形i=l:60;Plot(i,a1,'r',i,a2,'r:',i,b1,'b',i,b2,'b:',i,d1,'g',i,d2,'g:',i,d3,'g+')%画出各个被辨识参数title('ParameterIdentificationwithRecursiveLeastSquaresMethod')%标题figure(3);%画出第三个图形i=l:60;Plot(i,ea1,'r',i,ea2,'r:',i,eb1,'b',i,eb2,'b:',i,ed1,'g',i,ed2,'g:',i,ed2,'r+')%画出各个参数收敛情况title('IdentificationError')%标题%response%响应figure(4);%画出第四个图形subplot(4,l,l);%画出第四个图形中的四个子图的第一个子图i=l:60;plot(i,zs(i),'r')%画出被辨识系统在没有噪声情况下的实际输出响应subplot(4,l,2);%画出第四个图形中的四个子图的第二个子图i=l:60;plot(i,z(i),'g')%画出被辨识系统的采样输出响应subplot(4,l,3);%画出第四个图形中的四个子图的第三个子图i=1:60;plot(i,zmd(i),'b')%画出模型含有噪声的输出响应subplot(4,1,4);%画出第四个图形中的四个子图的第四个子图i=1:60;plot(i,zm(i),'b')%画出模型去除噪声后的输出响应【8】梯度校正的最小二乘辨识程序%FLch4GAeg1clearu=-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),holdonk=1:20plot(k,y)%给出初始值h1=-1,0,0'h2=-1,-1,0'g=0,0,0'I=1,0,0;0,1/2,0;0,0,1/4;h=h1,h2,zeros(3,16);%计算输入样本数据h(k)fork=3:18h(:,k)=u(k),u(k-1),u(k-2)'end%计算出权矩阵R(k)和估计值gfork=1:18a=h(l,k)人2+(h(2,k)人2)/2+(h(3,k)人2)/4;%按照式(4.45)开始计算权矩阵a1=1/a;R=a1*1;%按照式(4.45)计算出权矩阵g(:,k+l)=g(:,k)+R*h(:,k)*(y(k+l)-h(:,k)'*g(:,k);%按照式(4.44)计算脉冲响应估计值end%画出图形g1=g(1,;g2=g(2,;g3=g(3,;figure(2)k=1:19;subplot(121)plot(k,g1,'r',k,g2,'g',k,g3,'b'),gridon%计算模型输出值ym及系统输出与模型输出之间的误差Eyfork=1:18ym(k)=h(:,k)'*g(:,k);Ey(k)=y(k+1)-ym(k);endk=1:18;subplot(122)plot(k,Ey),gridong,ym,Ey%显示脉冲响应估计值、模型输出值及系统输出与模型输出之间的误差figure(3)%画出脉冲响应曲线x=0:1:3;y=0,g(1,18),g(2,18),g(3,18);xi=linspace(0,3);yi=interpl(x,y,xi,'cubic');%三次插值plot(x,y,'o',xi,yi,'m'),gridon%画出脉冲响应估计值及其三次插值曲线【9】递推的极大似然辨识程序clear%清零a(1)=1;b(1)=0;c(1)=1;d(1)=0;u(1)=d(1);z(1)=0;z(2)=0;%初始化fori=2:1200%产生m序列u(i)a(i)=xor(c(i-1),d(i-1);b(i)=a(i-1);c(i)=b(i-1);d(i)=c(i-1);u(i)=d(i);endu;%若取去；可以在程序运行中观测到m序列v=randn(1200,1);%产生正态分布随机数V=0;%计算噪声方差fori=1:1200V=V+v(i)*v(i);endV1=V/1200;fork=3:1200%根据v和u计算zz(k)=1.2*z(k-1)-0.6*z(k-2)+u(k-1)+0.5*u(k-2)+v(k)-v(k-1)+0.2*v(k-2);endo1=0.001*ones(6,1);p0=eye(6,6);%赋初值zf(1)=0.1;zf(2)=0.1;vf(2)=0.1;vf(1)=0.1;uf(2)=0.1;uf(1)=0.1;%迭代计算参数值和误差值fork=3:1200h=-z(k-1);-z(k-2);u(k-1);u(k-2);v(k-1);v(k-2);hf=h;K=p0*hf*inv(hf'*p0*hf+1);p=eye(6,6)-K*hf'*p0;v(k)=z(k)-h'*o1;o=o1+K*v(k);p0=p;o1=o;a1(k)=o(1);a2(k)=o(2);b1(k)=o(3);b2(k)=o(4);d1(k)=o(5);d2(k)=o(6);e1(k)=abs(a1(k)+1.2);e2(k)=abs(a2(k)-0.6);e3(k)=abs(b1(k)-1.0);e4(k)=abs(b2(k)-0.5);e5(k)=abs(d1(k)+1.0);e6(k)=abs(d2(k)-0.2);zf(k)=z(k)-d1(k)*zf(k-1)-d2(k)*zf(k-2);uf(k)=u(k)-d1(k)*uf(k-1)-d2(k)*uf(k-2);vf(k)=v(k)-d1(k)*vf(k-1)-d2(k)*vf(k-2);hf=-zf(k-1);-zf(k-2);uf(k-1);uf(k-2);vf(k-1);vf(k-2);end01%若取去；可以在程序运行中观测到参数V1%绘图subpl0t(4,1,1)k=1:1200;pl0t(k,a1,'k:',k,a2,'b',k,b1,'r',k,b2,'m:',k,d1,'g',k,d2,'k');xlabel('k')ylabel('parameter')legend('a1=-1.2,','a2=0.6','b1=1.0','b2=0.5','d1=-1.0','d2=0.2');%图标炷title('Theparameteridendificati0n0ftheRML');endsubpl0t(4,1,2)k=1:1200;pl0t(k,e1,'k',k,e2,'b',k,e3,'r',k,e4,'m',k,e5,'g',k,e6,'k');xlabel('k');ylabel('error');%title('误差曲线')endsubpl0t(4,1,3)k=1:1200;pl0t(k,u);xlabel('k');ylabel('input');%title('系统输入信号')endsubpl0t(4,1,4)k=1:1200;pl0t(k,v);xlabel('k');ylabel('randomnoise');')end