系统辨识及其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; for k=1:100 x2=A*x0; x1=mod (x2,M); v1=x1/256; v(:,k)=v1; x0=x1; v0=v1; end v2=v k1=k; %grapher k=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:N x2=A*x0; x1=mod (x2,M); v1=x1/256; v(:,k)=(v1-0.5)*f; x0=x1;

文档大全

实用标准文案

v0=v1; end v2=v k1=k; %grapher k=1:k1; plot(k,v,k,v,'r');

xlabel('k'), ylabel('v');title('(-1,+1)均匀分布的白噪声')

【3】M序列产生程序

X1=1;X2=0;X3=1;X4=0; %移位寄存器输入Xi初T态（0101）， Yi为移位寄存器各级输出 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==0 U(i)=-1; else U(i)=Y4; end end M=U %绘图 i1=i k=1:1:i1; plot(k,U,k,U,'rx') xlabel('k') ylabel('M序列')

title('移位寄存器产生的M序列')

【4】二阶系统一次性完成最小二乘辨识程序 %FLch3LSeg1

u=[-1,1,-1,1,1,1,1,-1,-1,-1,1,-1,-1,1,1]; %系统辨识的输入信号为一个周期的M序列 z=zeros(1,16); %定义输出观测值的长度 for k=3:16

z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2); %用理想输出值作为观测值 end

subplot(3,1,1) %画三行一列图形窗口中的第一个图形 stem(u) %画出输入信号u的经线图形

subplot(3,1,2) %画三行一列图形窗口中的第二个图形 i=1:1:16; %横坐标范围是1到16，步长为1

plot(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(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赋值

%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 、a2、 b1、 b2 %End

【5】实际压力系统的最小二乘辨识程序

%FLch3LSeg2 clear%工作间清零

V=[.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,1][alpha,log(beita)]'=HL*sita%注释P、V之间的关系 for i=1:6;%循环变量的取值为从1到6 Z(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赋值 %calculating parameters%计算参数

c1=HL'*HL; c2=inv(c1); c3=HL'*ZL; c4=c2*c3%计算 %Separation of Parameters%分离变量 alpha=c4(1) % 为c4的第1个元素

beita=exp(c4(2)) % 为以自然数为底的c4的第2个元素的指数

【6】递推的最小二乘辨识程序

%FLch3RLSeg3 clear%清理工作间变量 L=15;% M序列的周期

y1=1;y2=1;y3=1;y4=0;%四个移位积存器的输出初始值 for i=1

;%开始循环，长度为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)=-0.03;%如果M序列的值为\"1\"时,辨识的输入信号取“-0.03” else u(i)=0.03;%当M序列的值为\"0\"时,辨识的输入信号取“0.03” end%小循环结束

y1=x1;y2=x2;y3=x3;y4=x4;%为下一次的输入信号做准备 end%大循环结束，产生输入信号u figure(1);%第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=10^6*eye(4,4);%直接给出初始状态P0，即一个充分大的实数单位矩阵 E=0.000000005;%相对误差E=0.000000005 c=[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; %把当前相对变化的列向量加入误差矩阵的最后一列 c0=c1;%新获得的参数作为下一次递推的旧参数

c(:,k)=c1;%把辨识参数c 列向量加入辨识参数矩阵的最后一列 p1=p0-k1*k1'*[h1'*p0*h1+1];%求出 p(k)的值 p0=p1;%给下次用

if e2<=E break;%若参数收敛满足要求，终止计算 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到15

plot(i,a1,'r',i,a2,':',i,b1,'g',i,b2,':') %画出a1，a2，b1，b2的各次辨识结果 title('Parameter Identification with Recursive Least Squares Method')%图形标题 figure(3); %第3个图形 i=1:15; %横坐标从1到15

文档大全

实用标准文案

plot(i,ea1,'r',i,ea2,'g',i,eb1,'b',i,eb2,'r:') %画出a1，a2，b1，b2的各次辨识结果的收敛情况 title('Identification Precision') %图形标题

【7】增广的最小二乘辨识程序

%FLch3ELSeg4 clear

L=60;%四位移位积存器产生的M序列的周期

y1=1;y2=1;y3=1;y4=0;%四个移位积存器的输出初始值 for i=1

;

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\"时,辨识的输入信号取“-1” else u(i)=1;%M序列的值为\"0\"时,辨识的输入信号取“1” end

y1=x1;y2=x2;y3=x3;y4=x4;%为下一次的输入信号作准备 end

figure(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=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.001 0.001 0.001 0.001 0.001 0.001 0.001]';%直接给出被辨识参数的初始值,即一个充分小的实向量

p0=10^6*eye(7,7);%直接给出初始状态P0，即一个充分大的实数单位矩阵 E=5.0e-15;%取相对误差E

c=[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; %K d1=z(k)-h1'*c0; c1=c0+k1*d1;%辨识参数c

zs(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 列向量加入辨识参数矩阵 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、 b2 d1=c(5,; d2=c(6,; d3=c(7,; %分离出d1、 d2、 d3

ea1=e(1,; ea2=e(2,; eb1=e(3,; eb2=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 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); %画出第四个图形中的四个子图的第一个子图 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')%画出模型去除噪声后的输出响应

【8】梯度校正的最小二乘辨识程序

%FLch4GAeg1 clear

u=[-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 on k=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;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:18

a=h(1,k)^2+(h(2,k)^2)/2+(h(3,k)^2)/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(2,; g3=g(3,; figure(2) k=1:19; subplot(121)

plot(k,g1,'r',k,g2,'g',k,g3,'b'),grid on

%计算模型输出值ym及系统输出与模型输出之间的误差Ey for k=1:18

ym(k)=h(:,k)'*g(:,k); Ey(k)=y(k+1)-ym(k); end k=1:18; subplot(122) plot(k,Ey),grid on

文档大全

实用标准文案

g, 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,'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)=b(i-1); d(i)=c(i-1); u(i)=d(i); end

u; %若取去‘；’可以在程序运行中观测到m序列

v=randn(1200,1); %产生正态分布随机数 V=0; %计算噪声方差 for i=1:1200 V=V+v(i)*v(i); end

V1=V/1200;

for k=3:1200 %根据v和u计算z

z(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); end

o1=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; %迭代计算参数值和误差值 for k=3:1200

h=[-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)]; end

o1 %若取去‘；’可以在程序运行中观测到参数 V1 %绘图 subplot(4,1,1) k=1:1200;

plot(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('The parameter idendification of the RML'); end

subplot(4,1,2) k=1:1200;

plot(k,e1,'k',k,e2,'b',k,e3,'r',k,e4,'m',k,e5,'g',k,e6,'k'); xlabel('k'); ylabel('error'); %title('误差曲线') end

subplot(4,1,3) k=1:1200; plot(k,u);

xlabel('k'); ylabel('input'); %title('系统输入信号') end

subplot(4,1,4) k=1:1200; plot(k,v);

文档大全

实用标准文案

xlabel('k'); ylabel('random noise'); %title('系统所加的随机噪声') end

文档大全