系统辨识及其matlab仿真

1、【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;for k=1:100x1=mod (x2,M);v1=x1

2、/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('(0,1) 均匀分布的随机序列 ')【2】白噪声产生程序A=6; x0=1; M=255; f=2; N=100;for k=1:Nx2=A*x0;x1=mod (x2,M);v1=x1/256;v(:,k)=(v1-0.5)*f;v0=v1;endv2=vk1=k;%grapherk=1:k1;plot(k,v,k,v,&#

3、39;r');xlabel('k'), ylabel('v');title('(-1,+1) 均匀分布的白噪声 ')【 3 】 M 序列产生程序Yi 为移位寄存器各级输出X1=1;X2=0;X3=1;X4=0; %移位寄存器输入 Xi 初 T 态(0101 )m=60; % 置 M 序列总长度for i=1:m %1#Y4=X4; Y3=X3; Y2=X2; Y1=X1;X4=Y3; X3=Y2; X2=Y1;X1=xor(Y3,Y4); % 异或运算if Y4=0U(i)=-1;else endendM=U%绘图i1=ik=1:1:i1

4、;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); % 定义输出观测值的长度for k=3:16z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2); % 用理想输出值作为观测值endsubplot(

5、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,z% 显示输入信号和输出观测信号%L=14% 数据长度HL=-z(2) -z(1) u(2) u(1);-z(3) -z(

6、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) % 给样

7、本矩阵 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 赋 值%calculating parameters% 计算参数c1=HL'*HL; c2=inv(c1); c3=HL'*ZL; c=c2*c3 % 计算并显示%DISPLAY PARAMETERSa1=c(1), a2=c(2), b1=c(3), b2=c(4) % 从 中分离出并显示 a1 、a2、 b1 、 b2%End【5】实际压力系统的最小二乘辨识程序%FLch3LSe

8、g2clear% 工作间清零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 之间的关系for i=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

9、(V(4),1;-log(V(5),1;-log(V(6),1 %给 HL 赋值%calculat ing parameters% 计算参数c1=HL'*HL; c2=i nv(c1); c3=HL'*ZL; c4=c2*c3% 计算%Separati on of Parameters% 分离变量alpha=c4(1) % 为c4的第1个元素beita=exp(c4(2) %为以自然数为底的c4的第2个元素的指数【6】递推的最小二乘辨识程序%FLch3RLSeg3clear%清理工作间变量L=15;% M序列的周期y1=1；y2=l；y3=1;y4=0;%四个移位积存器的输岀初

10、始值for i=1乞;开始循环，长度为Lx1=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)=-0.03;% 如果 M 序列的值为 "1" 时,辨识的输入信号取 “-0.03 ”else u(i

11、)=0.03;% 当 M 序列的值为 "0" 时 ,辨识的输入信号取 “0.03 ”end% 小循环结束y1=x1;y2=x2;y3=x3;y4=x4;% 为下一次的输入信号做准备end% 大循环结束，产生输入信号 ufigure(1);% 第 1 个图形stem(u),grid on% 以径的形式显示出输入信号并给图形加上网格z(2)=0;z(1)=0;% 取 z 的前两个初始值为零for k=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 递推

12、最小二乘辨识c0=0.001 0.001 0.001 0.001'% 直接给出被辨识参数的初始值 ,即一个充分小的实向量p0=10A6*eye(4,4);%直接给岀初始状态P0，即一个充分大的实数单位矩阵E=0.000000005;% 相对误差 E=0.000000005c=c0,zeros(4,14);% 被辨识参数矩阵的初始值及大小e=zeros(4,15);% 相对误差的初始值及大小for k=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=p

13、0*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;% 给下次用if e2<=E break;% 若参数收敛满足要求，终止计算end% 小循

14、环结束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 到 15图形标题plot(i,a1,'r',i,a2,':',i,b1,'g',i,b2,':') % 画出 a1，a2，b1 ，b2 的各次辨识结果 title('Parameter IdentifiCation with

15、 ReCursive Least Squares Method')% figure(3); % 第 3 个图形i=1:15; % 横坐标从 1 到 15plot(i,ea1,'r',i,ea2,'g',i,eb1,'b',i,eb2,'r:') % 画岀 al , a2 , bl , b2 的各次辨识结果的收敛情况 ti tle('lde ntificatio n Precisio n') % 图形标题【7】增广的最小二乘辨识程序%FLch3ELSeg4clearL=60;%四位移位积存器产生的 M序列的周

16、期y1=1y2=1;y3=1;y4=0;%四个移位积存器的输岀初始值x1=xor(y3,y4);%第一个移位积存器的输入信号x2=y1;%第二个移位积存器的输入信号x3=y2;%第三个移位积存器的输入信号x4=y3;%第四个移位积存器的输入信号y(i)=y4;%第四个移位积存器的输岀信号，幅值"0"和"1"if y(i)>0.5,u(i)=-1;%M 序列的值为"1"时，辨识的输入信号取-1else u(i)=1;%M 序列的值为"0"时，辨识的输入信号取1 ”endy1=x1;y2=x2;y3=x3;y4=

17、x4;%为下一次的输入信号作准备endfigure(1);%画第一个图形subplot(2,1,1); %画第一个图形的第一个子图stem(u),grid on%画岀 M 序列输入信号v=randn(1,60); %产生一组 60 个正态分布的随机噪声subplot(2,1,2); %画第一个图形的第二个子图plot(v),grid on;%画岀随机噪声信号R=corrcoef(u,v);%计算输入信号与随机噪声信号的相关系数r=R(1,2);% 取出互相关系数u% 显示输入型号v% 显示噪声型号z=zeros(7,60);zs=zeros(7,60);zm=zeros(7,60);zmd=z

18、eros(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.001 0.001 0.001 0.001 0.001 0.001 0.001'%直接给出被辨识参数的初始值 ,即一个充分小的实向量p0=10A6*eye(7,7);%直接给岀初始状态P0，即一个充分大的实数单位矩阵E=5.0e-15;% 取相对误差 Ec

19、=c0,zeros(7,59);% 被辨识参数矩阵的初始值及大小e=zeros(7,60);% 相对误差的初始值及大小for k=3:60; % 开始求 K z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2)+v(k)-v(k-1)+0.2*v(k-2);% 系统在 M 序列输入下的 输出采样信号h1=-z(k-1),-z(k-2),u(k-1),u(k-2),v(k),v(k-1),v(k-2)'%为求 K(k) 作准备x=h1'*p0*h1+1; x1=inv(x); k1=p0*h1*x1; %Kd1=z(k)-h1'*c0;

20、 c1=c0+k1*d1;% 辨识参数 czs(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2);% 系统在 M 序列的输入下不考虑扰动时的输出响 应zm(k)=-z(k-1),-z(k-2),u(k-1),u(k-2)*c1(1);c1(2);c1(3);c1(4);% 模型在 M 序列的输入下不考虑扰 动时的的输出响应zmd(k)=h1'*c1;% 模型在 M 序列的输入下的的输出响应e1=c1-c0;e2=e1./c0; % 求参数的相对变化e(:,k)=e2;c0=c1;% 给下一次用c(:,k)=c1;% 把辨识参数 c 列向量加入辨识参数

21、矩阵p1=p0-k1*k1'*h1'*p0*h1+1;%find p(k)p0=p1;% 给下次用if e2<=E break;% 若收敛情况满足要求，终止计算end% 判断结束end% 循环结束c, e, % 显示被辨识参数及参数收敛情况z,zmd % 显示输出采样值、系统实际输出值、模型输出值%分离变量a1=c(1,; a2=c(2,; b1=c(3,; b2=c(4,;% 分离出 a1 、 a2 、 b1 、 b2d1=c(5,; d2=c(6,; d3=c(7,; % 分离出 d1、 d2 、 d3ea1=e(1,; ea2=e(2,; eb1=e(3,; eb2

22、=e(4,; %分离出 a1、 a2、 b1 、 b2 的收敛情况ed1=e(5,; ed2=e(6,; ed3=e(7,; %分离出 d1 、 d2 、d3 的收敛情况figure(2);% 画第二个图形i=1:60;plot(i,a1,'r',i,a2,'r:',i,b1,'b',i,b2,'b:',i,d1,'g',i,d2,'g:',i,d3,'g+')% 画出各个被辨识参数title('Parameter Identification with Recursive

23、Least Squares Method')%标题figure(3);% 画出第三个图形i=1:60;plot(i,ea1,'r',i,ea2,'r:',i,eb1,'b',i,eb2,'b:',i,ed1,'g',i,ed2,'g:',i,ed2,'r+')% 画出各个参数收敛情况title('Identification Error')%标题%response% 响应figure(4);% 画出第四个图形subplot(4,1,1); % 画出第四个图形中

24、的四个子图的第一个子图i=1:60;plot(i,zs(i),'r')% 画出被辨识系统在没有噪声情况下的实际输出响应subplot(4,1,2); % 画出第四个图形中的四个子图的第二个子图i=1:60;plot(i,z(i),'g')% 画出被辨识系统的采样输出响应subplot(4,1,3); % 画出第四个图形中的四个子图的第三个子图i=1:60;plot(i,zmd(i),'b')%画出模型含有噪声的输出响应subplot(4,1,4); %画出第四个图形中的四个子图的第四个子图i=1:60;plot(i,zm(i),'b

25、9;)%画出模型去除噪声后的输出响应【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), hold onk=1:20plot(k,y)%给出初始值h1=-1,0,0'h2=-1,-1,0'g=0,0,0'

26、I=1,0,0;0,1/2,0;0,0,1/4;h=h1,h2,zeros(3,16);%计算输入样本数据 h(k)h(:,k)=u(k),u(k-1),u(k-2)'end%计算出权矩阵 R(k) 和估计值 gfor k=1:18a=h(1,k)A2+(h(2,k)A2)/2+(h(3,kF2)/4%按照式(4.45)开始计算权矩阵a1=1/a;R=a1*I;% 按照式( 4.45 )计算出权矩阵g(:,k+1)=g(:,k)+R*h(:,k)*(y(k+1)-h(:,k)'*g(:,k);% 按照式( 4.44)计算脉冲响应估计值end%画出图形g1=g(1,; g2=g(

27、2,; g3=g(3,;figure(2)k=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

28、=0,g(1,18),g(2,18),g(3,18);xi=linspace(0,3); yi=interp1(x,y,xi,'cubic');% 三次插值plot(x,y,'o',xi,yi,'m'),grid on% 画出脉冲响应估计值及其三次插值曲线【9】递推的极大似然辨识程序clear % 清零a(1)=1;b(1)=0;c(1)=1;d(1)=0;u(1)=d(1);z(1)=0;z(2)=0; % 初始化for i=2:1200 % 产生 m 序列 u(i)a(i)=xor(c(i-1),d(i-1);b(i)=a(i-1);c(i)

29、=b(i-1);d(i)=c(i-1);u(i)=d(i);endu; %若取去；'可以在程序运行中观测到 m 序列v=randn(1200,1); % 产生正态分布随机数V=0; % 计算噪声方差for i=1:1200V=V+v(i)*v(i);endV1=V/1200;for k=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

30、(2)=0.1;vf(1)=0.1;uf(2)=0.1;uf(1)=0.1;%迭代计算参数值和误差值for k=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