系统辨识报告.docx
《系统辨识报告.docx》由会员分享，可在线阅读，更多相关《系统辨识报告.docx（34页珍藏版）》请在冰点文库上搜索。
系统辨识报告
实验一：
系统辨识的经典方法
实验目的：
掌握系统的数学模型与系统的输入，输出信号之间的关系，掌握经典辨识的实验测试方法和数据处理方法。
熟悉matlab/Simulink环境。
实验内容：
1.用系统阶跃响应法测试给定系统的数学模型。
在系统没有噪声干扰的条件下通过测试系统的阶跃响应获得系统的一阶加纯滞后或二阶加纯滞后模型，对模型进行验证。
2.在被辨识的系统加入噪声干扰，重复上述1的实验过程。
实验内容：
利用非线性水槽模型搭建单水槽系统模型，如下图
由上图及其计算的H：
H为一二元数组，分别表示第一个、第二个水箱的液位。
二．实验方法：
运用阶跃响应法：
第一个水箱的参数辨识（一阶）：
一阶惯性环节的传递函数为：
其中：
将H归一化：
在H中查得
时对应T=2.3
故模型为
第二个水箱的参数辨识（二阶）：
二阶系统的传递函数为：
其中：
在H中可得：
故有：
由公式
可求出：
故：
实验二：
相关分析法
搭建对象：
处理程序：
fori=1:
15
m（i,:
）=UY（32-i:
46-i,1）;
end
y=UY（31:
45,2）;
gg=ones（15）+eye（15）;
g=1/（25*16*2）*gg*m*y;
plot（g）;
holdon;
stem（g）;
实验结果：
相关分析法
最小二乘法建模：
二、三次实验
本次实验要完成的内容：
1．参照index2，设计对象，从workspace空间获取数据，取二阶，三阶对象实现最小二乘法的一次完成算法和最小二乘法的递推算法（LSandRLS）;
2．对设计好的对象，在时间为200-300之间，设计一个阶跃扰动，用最小二乘法和带遗忘因子的最小二乘法实现，对这两种算法的特点进行说明；
实验三最小二乘法参数估计
一．实验内容
（1）参照index2，设计对象，从workspace空间获取数据，取二阶，三阶对象实现最小二乘法的一次完成算法和最小二乘法的递推算法（LSandRLS）;
（2）对设计好的对象，在时间为200-300之间，设计一个阶跃扰动，用最小二乘法和带遗忘因子的最小二乘法实现，对这两种算法的特点进行说明
（3）参照index3，设计符合GLS和ELS的对象模型，改写参照程序，实现相应的.算法。
二．实验源程序
内容
（1）：
（a）%最小二乘法的一次完成算法
M=UY（:
1）;
z=UY（:
2）;
H=zeros（100,4）;
fori=1:
100
H（i,1）=-z（i+1）;
H（i,2）=-z（i）;
H（i,3）=M（i+1）;
H（i,4）=M（i）;
end
Estimate=inv（H'*H）*H'*（z（3:
102））
Estimate=
-0.7867
0.1388
0.0007
0.5708
0.3116
（b）%最小二乘的递推算法
M=UY（:
1）;
z=UY（:
2）;
P=1000*eye（5）;%¹À¼Æ·½²î
%Pstore=zeros（5,200）;
%Pstore（:
1）=[P（1,1）,P（2,2）,P（3,3）,P（4,4）,P（5,5）]';
Theta=zeros（5,200）;%²ÎÊýµÄ¹À¼ÆÖµ£¬´æ·ÅÖмä¹ý³Ì¹ÀÖµ
Theta（:
1）=[0;0;0;0;0];
K=zeros（4,400）;%ÔöÒæ¾ØÕó
K=[10;10;10;10;10];
fori=3:
201
h=[-z（i-1）;-z（i-2）;M（i）;M（i-1）;M（i-2）];
K=P*h*inv（h'*P*h+1）;
Theta（:
i-1）=Theta（:
i-2）+K*（z（i）-h'*Theta（:
i-2））;
P=（eye（5）-K*h'）*P;
%Pstore（:
i-1）=[P（1,1）,P（2,2）,P（3,3）,P（4,4）,P（5,5）]';
end
i=1:
200;
figure
（1）
plot（i,Theta（1,:
）,i,Theta（2,:
）,i,Theta（3,:
）,i,Theta（4,:
）,i,Theta（5,:
））
title（'´ý¹À²ÎÊý¹ý¶É¹ý³Ì'）
grid
内容
（2）：
（b）%带遗忘因子的最小二乘法
M=UY（:
1）;
z=UY（:
2）;
P=1000*eye（5）;%¹À¼Æ·½²î
Theta=zeros（5,200）;%²ÎÊýµÄ¹À¼ÆÖµ£¬´æ·ÅÖмä¹ý³Ì¹ÀÖµ
Theta（:
1）=[0;0;0;0;0];
K=zeros（4,400）;%ÔöÒæ¾ØÕó
K=[10;10;10;10;10];
lamda=0.99;%ÒÅÍüÒò×Ó
fori=3:
201
h=[-z（i-1）;-z（i-2）;M（i）;M（i-1）;M（i-2）];
K=P*h*inv（h'*P*h+lamda）;
Theta（:
i-1）=Theta（:
i-2）+K*（z（i）-h'*Theta（:
i-2））;
P=（eye（5）-K*h'）*P/lamda;
end
i=1:
200;
figure
（1）
plot（i,Theta（1,:
）,i,Theta（2,:
）,i,Theta（3,:
）,i,Theta（4,:
）,i,Theta（5,:
）
title（'´ý¹À²ÎÊý¹ý¶É¹ý³Ì'）
grid
内容（3）：
（a）%广义最小二乘法辨识GLS:
clc
n=2;
N=799;
u=UY（:
1）;
y=UY（:
2）;
Y=y（n+1:
n+N）;
phi1=-y（n:
n+N-1）;
phi2=-y（n-1:
n+N-2）;
phi3=u（n:
n+N-1）;
phi4=u（n-1:
n+N-2）;
phi=[phi1phi2phi3phi4];
theta=inv（phi'*phi）*phi'*Y;
%%%%%%ÒÔÉÏΪ×îС¶þ³Ë¼ÆËãTheta¹À¼ÆÖµ%%%%%%%%%
f0=[10;10];
while1
e（n+1:
n+N,1）=y（n+1:
n+N）-phi*theta;
omega（:
1）=-e（n:
n+N-1）;
omega（:
2）=-e（n-1:
n+N-2）;
f=inv（omega'*omega）*omega'*e（n+1:
n+N,1）;
if（f-f0）'*（f-f0）<0.0001
break
end
f0=f;
%%%%%%ÒÔÉÏΪ×îС¶þ³Ë¼ÆËãf¹À¼ÆÖµ%%%%%%%%%
fori=n+1:
n+N
yy（i,1）=[y（i）,y（i-1）,y（i-2）]*[1,f
（1）,f
（2）]';
uu（i,1）=[u（i）,u（i-1）,u（i-2）]*[1,f
（1）,f
（2）]';
end
yy（1,1）=y（1,1）;
uu（2,1）=u（2,1）;
%%%%%%%%%%%ÒÔÉÏΪÊý¾ÝÂ˲¨%%%%%%%%
phi（:
1）=-yy（n:
n+N-1）;
phi（:
2）=-yy（n-1:
n+N-2）;
phi（:
3）=uu（n:
n+N-1）;
phi（:
4）=uu（n-1:
n+N-2）;
theta=inv（phi'*phi）*phi'*yy（n+1:
n+N）;
end
f
theta
（b）%增广最小二乘法的递推算法（EM）：
u=UY（:
1）;
y=UY（:
2）;
n=2;
N=799;
thitak=ones（7,1）*0.001;
p1=eye（7,7）*（1.0e6）;p2=zeros（7,7）;
K=zeros（7,1）;e=zeros（801,1）;I=eye（7,7）;
fori=3:
801
phi1=-y（i-1）;
phi2=-y（i-2）;
phi3=u（i）;
phi4=u（i-1）;
phi5=u（i-2）;
phi6=e（i-1）;
phi7=e（i-2）;
phi=[phi1phi2phi3phi4phi5phi6phi7];
K=p1*phi'/（phi*p1*phi'+1）;
p2=（I-K*phi）*p1;
thita（:
i）=thitak+K*（y（i,1）-phi*thitak）;
p1=p2;
thitak=thita（:
i）;
e（i）=y（i）-phi*thitak;
end
thitak
thitak=
-0.7426
0.0821
0.0000
0.2364
0.1038
0.0010
0.0010
三.实验结果
内容
（1）
（b）最小二乘的递推算法的实验结果
内容
（2）
（a）用最小二乘法得出的辨识结果
（b）用带遗忘因子的最小二乘法得到的结果
当λ=0.99时：
当λ=0.98时：
当λ=0.97时：
当λ=0.95时：
当λ=0.93时：
当λ=0.92时：
内容（3）：
（a）广义最小二乘法辨识GLS得到的结果:
f=
-0.3506
0.2873
theta=
-0.7411
0.0809
0.2366
0.1041
则系统的差分方程为
（b）增广最小二乘法的递推算法（EM）得到的结果：
四．实验结果分析
比较以上六个图可以知道：
加入扰动时，最小二乘法的递推算法不能将扰动引起的干扰消除，其辨识出来的参数不能迅速复原；而渐消记忆的递推算法在扰动加入后，能较快复原，且随着遗忘因子的减小，其效果越明显；
五．附录：
系统辨识大作业源程序清单
作业一：
X=[1.01;2.03;3.02;4.01;5;6.02;7.03;8.04;9.03;10];
Y=[9.6;4.1;1.3;0.4;0.05;0.1;0.7;1.7;3.8;9.1];
X0=ones（10,1）;
A=[X.^2,X,X0];
E=inv（A'*A）*A'*Y;
YE=A*E;
Dif=abs（Y）-abs（YE）;
plot（X,Y,'-rs',X,YE,'--gs',X,Dif,'.'）;
legend（'观测值Y','拟合值YE','差值Y-YE'）;
xlabel（'x'）;
ylabel（'y'）;
title（'作业一y=ax^2+bx+c参数求解'）;
作业二：
%利用相关法求解系统的脉冲响应。
Np=15;
a=5;
delta=2;
C=4*Np+1;
fori=1:
Np
C=C-1;
forj=1:
Np
U（i,j）=UY（C+j-1,1）;
end
Y（i）=UY（i+4*Np-1,2）;
end
R=ones（15）+eye（15）;
G=R*U*Y'/（a^2*（Np+1）*delta）;
t=0:
2:
2*Np-2;
stem（t,G,'gs'）
holdon
plot（t,G,'--'）
%利用脉冲相应求解系统G（z）
n=2;%系统阶次是2
ga=[G（2+i）G（3+i）;G（3+i）G（4+i）];
y=-G（4+i:
5+i）;
az=inv（ga）*y;
r=[100;az
（2）10;az
（1）az
（2）1];
gb=G（1+i:
3+i）;
bz=r*gb;
Gz=tf（bz',[1az'],2）;
%利用脉冲相应求解系统G（s）
as=ga\（-G（1:
2））;
x=roots（[as
（2）as
（1）1]）;
s=log（x）;
cs=[11;x
（1）x
（2）]\G（1:
2）;
den=conv（[1-s
（1）],[1-s
（2）]）;num=cs
（1）*cs
（2）;
Gs=tf（num,den）
%利用相关最小二乘法求系统的参数
clc
z=UY（:
2）;
M=UY（:
1）;
%递推求解
P=100*eye（4）;%估计方差
Pstore=zeros（4,301）;
Pstore（:
1）=[P（1,1）,P（2,2）,P（3,3）,P（4,4）];
Theta=zeros（4,301）;%参数的估计值，存放中间过程估计值
Theta（:
1）=[3;3;3;3];
Theta（:
2）=[3;3;3;3];
Theta（:
3）=[3;3;3;3];
Theta（:
4）=[3;3;3;3];
K=[10;10;10;10];
fori=5:
301
h=[-z（i-1）;-z（i-2）;M（i-1）;M（i-2）];
hstar=[M（i-1）;M（i-2）;M（i-3）;M（i-4）];%辅助变量
K=P*hstar*inv（h'*P*hstar+1）;
Theta（:
i）=Theta（:
i-1）+K*（z（i）-h'*Theta（:
i-1））;
P=（eye（4）-K*h'）*P;
Pstore（:
i-1）=[P（1,1）,P（2,2）,P（3,3）,P（4,4）];
end
%==================结果输出及作图===================
disp（'参数a1a2b1b2的估计结果：
'）
Theta（:
301）
i=1:
301;
figure
（1）
plot（i,Theta（1,:
）,i,Theta（2,:
）,i,Theta（3,:
）,i,Theta（4,:
））
title（'参数过渡过程'）
figure
（2）
plot（i,Pstore（1,:
）,i,Pstore（2,:
）,i,Pstore（3,:
）,i,Pstore（4,:
））
title（'估计方差变化过程'）
%最小二乘递推算法（RLS）求解模型参数
clc
M=UY（:
1）;
z=UY（:
2）;
P=100*eye（5）;%估计方差
Pstore=zeros（5,200）;
Pstore（:
1）=[P（1,1）,P（2,2）,P（3,3）,P（4,4）,P（5,5）]';
Theta=zeros（5,200）;%参数的估计值，存放中间过程估值
Theta（:
1）=[0;0;0;0;0];
%K=zeros（4,400）;%增益矩阵
K=[10;10;10;10;10;10;10];
fori=3:
201
h=[-z（i-1）;-z（i-2）;M（i）;M（i-1）;M（i-2）];
K=P*h*inv（h'*P*h+1）;
Theta（:
i-1）=Theta（:
i-2）+K*（z（i）-h'*Theta（:
i-2））;
P=（eye（5）-K*h'）*P;
Pstore（:
i-1）=[P（1,1）,P（2,2）,P（3,3）,P（4,4）,P（5,5）]';
end
i=1:
200;
figure
（1）
plot（i,Theta（1,:
）,i,Theta（2,:
）,i,Theta（3,:
）,i,Theta（4,:
）,i,Theta（5,:
））
title（'待估参数过渡过程'）
figure
（2）
plot（i,Pstore（1,:
）,i,Pstore（2,:
）,i,Pstore（3,:
）,i,Pstore（4,:
）,i,Pstore（5,:
））
title（'估计方差变化过程'）
Theta
%带遗忘因子的RLS算法
clc
M=UY（:
1）;
z=UY（:
2）;
larmda=0.99;%遗忘因子
P=100*eye（5）;%估计方差
Pstore=zeros（5,200）;
Pstore（:
1）=[P（1,1）,P（2,2）,P（3,3）,P（4,4）,P（5,5）]';
Theta=zeros（5,200）;%参数的估计值，存放中间过程估值
Theta（:
1）=[0;0;0;0;0];
K=[10;10;10;10;10;10;10];
fori=3:
201
h=[-z（i-1）;-z（i-2）;M（i）;M（i-1）;M（i-2）];
K=P*h*inv（h'*P*h+larmda）;
Theta（:
i-1）=Theta（:
i-2）+K*（z（i）-h'*Theta（:
i-2））;
P=（eye（5）-K*h'）*P/larmda;
Pstore（:
i-1）=[P（1,1）,P（2,2）,P（3,3）,P（4,4）,P（5,5）]';
end
i=1:
200;
figure
（1）
plot（i,Theta（1,:
）,i,Theta（2,:
）,i,Theta（3,:
）,i,Theta（4,:
）,i,Theta（5,:
））
title（'待估参数过渡过程'）
figure
（2）
plot（i,Pstore（1,:
）,i,Pstore（2,:
）,i,Pstore（3,:
）,i,Pstore（4,:
）,i,Pstore（5,:
））
title（'估计方差变化过程'）
Theta
%%%%%%%%%%%%%%阶次辨识（残差最小）%%%%%%%%%%%
clc
y=UY（:
2）;
u=UY（:
1）;
%subplot（2,2,1）;
forn=1:
5%对各阶次计算Jn
N=301-n;
Y=y（n+1:
n+N,1）;
forj=1:
N
fork=1:
n
phi（j,k）=-y（n+j-k）;
phi（j,n+k）=u（n+j-k）;
end
end
theta=inv（phi'*phi）*phi'*Y;
e（n+1:
n+N,1）=y（n+1:
n+N）-phi*theta;
J（n,1）=e'*e;%计算Jn
phi=0;
m（n）=n;
end
plot（m,J,'r-o'）;
title（'数据阶次曲线图'）;
holdon
gridon
%************阶次计算AIC***************%
Data=UY;
%按Akkaike信息准则定阶
L=length（Data）;
U=Data（:
1）;
Y=Data（:
2）;
forn=1:
1:
5
N=301-n
yd（1:
N,1）=Y（n+1:
n+N）;
fori=1:
N
forl=1:
1:
2*n
if（l<=n）
FIA（i,l）=-Y（n+i-l）;
else
FIA（i,l）=U（2*n+i-l）;
end
end
end
thita=inv（FIA'*FIA）*FIA'*yd;
omiga=（yd-FIA*thita）'*（yd-FIA*thita）/N;
AIC（n）=N*log（omiga）+4*n;
end
plot（AIC,'-'）;
title（'AIC随阶次的变化曲线'）;
xlabel（'n'）;
ylabel（'AIC'）;
作业三:
广义最小二乘法（GLS）：
%广义最小二乘法辨识
clc
n=2;
N=799;
%y=y1;
%u=u1;
y=UY（:
2）;
u=UY（:
1）;
%y=y3;
%u=u3;
Y=y（n+1:
n+N）;
phi1=-y（n:
n+N-1）;
phi2=-y（n-1:
n+N-2）;
phi3=u（n:
n+N-1）;
phi4=u（n-1:
n+N-2）;
phi=[phi1phi2phi3phi4];
theta=inv（phi'*phi）*phi'*Y
%%%%%%以上为最小二乘计算Theta估计值%%%%%%%%%
f0=[10;10];
while1
e（n+1:
n+N,1）=y（n+1:
n+N）-phi*theta;
omega（:
1）=-e（n:
n+N-1）;
omega（:
2）=-e（n-1:
n+N-2）;
f=inv（omega'*omega）*omega'*e（n+1:
n+N,1）;
if（f-f0）'*（f-f0）<0.0001
break
end
f0=f;
%%%%%%以上为最小二乘计算f估计值%%%%%%%%%
fori=n+1:
n+N
yy（i,1）=[y（i）,y（i-1）,y（i-2）]*[1,f
（1）,f
（2）]';
uu（i,1）=[u（i）,u（i-1）,u（i-2）]*[1,f
（1）,f
（2）]';
end
yy（1,1）=y（1,1）;
uu（2,1）=u（2,1）;
%%%%%%%%%%%以上为数据滤波%%%%%%%%
phi（:
1）=-yy（n:
n+N-1）;
phi（:
2）=-yy（n-1:
n+N-2）;
phi（:
3）=uu（n:
n+N-1）;
phi（:
4）=uu（n-1:
n+N-2）;
theta=inv（phi'*phi）*phi'*yy（n+1:
n+N）;
end
f
theta
ELS：
%增广最小二乘的递推算法
%===========产生均值为0，方差为1的高斯白噪声=========
v
（1）=0;
v
（2）=0;
z=UY（:
2）;
M=UY（:
1）;
%递推求解
P=100*eye（6）;%估计方差
Pstore=zeros（6,800）;
Pstore（:
1）=[P（1,1）,P（2,2）,P（3,3）,P（4,4）,P（5,5）,P（6,6）];
Theta=zeros（6,401）;%参数的估计值，存放中间过程估值
Theta（:
1）=[3;3;3;3;3;3];
K=[10;10;10;10;10;10];
fori=3:
801
h=[-z（i-1）;-z（i-2）;M（i-1）;M（i-2）;v（i-1）;v（i-2）];
K=P*h*inv（h'*P*h+1）;
Theta（:
i-1）=Theta（:
i-2）+K*（z（i）-h'*Theta（:
i-2））;
v（i）=z（i）-h'*Theta（:
i-2）;
P=（eye（6）-K*h'）*P;
Pstore（:
i-1）=[P（1,1）,P（2,2）,P（3,3）,P（4,4）,P（5,5）,P（6,6）];
end
%=======================================================================
disp（'参数a1、a2、b1、b2、d1、d2估计结果：
'）
Theta（:
800）
i=1:
800;
figure
（1）
plot（i,Theta（1,:
）,i,Theta（2,:
）,i,Theta（3,:
）,i,Theta（4,:
）,i,Theta（5,:
）,i,Theta（6,:
））
title（'待估参数过渡过程'）
figure
（2）
plot（i,Pstore（1,:
）,i,Pstore（2,:
）,i,Ps
升级会员 