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



参考文献：
[1]系统辨识及其matlab仿真。侯媛彬，汪梅，王立琦。科学出版社2004。
[2]过程辨识。方崇智。清华大学出版社，1987。
[3]系统辨识讲义。伍清河，北京理工大学，2005。（见
http://www.

主要内容：
【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、

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

%FLch3LSeg2
clear%工作间清零
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,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;%若参数收敛满足要求，终止计算

【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') %图形标题
。。。。。。。

，飘过