系统辨识实验报告 2.docx
《系统辨识实验报告 2.docx》由会员分享，可在线阅读，更多相关《系统辨识实验报告 2.docx（29页珍藏版）》请在冰点文库上搜索。
系统辨识实验报告2
自动化09-3宋佳瑛09051304
系统辨识实验报告
实验一：
系统辨识的经典方法
实验目的：
掌握系统的数学模型与系统的输入，输出信号之间的关系，掌握经典辨识的实验测试方法和数据处理方法。
熟悉matlab/Simulink环境。
实验内容：
1.用系统阶跃响应法测试给定系统的数学模型。
在系统没有噪声干扰的条件下通过测试系统的阶跃响应获得系统的一阶加纯滞后或二阶加纯滞后模型，对模型进行验证。
2.在被辨识的系统加入噪声干扰，重复上述1的实验过程。
1.没有噪声
搭建对象
测试对象流程图
实验结果为：
2、加入噪声干扰
搭建对象
实验结果：
加入噪声干扰之后水箱输出不平稳，有波动。
实验二：
相关分析法
搭建对象：
处理程序：
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之间，设计一个阶跃扰动，用最小二乘法和带遗忘因子的最小二乘法实现，对这两种算法的特点进行说明；
实验内容结果与程序代码：
利用LS和RLS得到的二阶，三阶参数
算法
阶次
A1
A2
A3
B0
B1
B2
B3
LS
二阶
-0.7842
0.1373
-0.0036
0.5668
0.3157
RLS
二阶
-0.7824
0.1373
-0.0036
0.5668
0.3157
LS
三阶
-0.4381
-0.1228
0.0407
-0.0078
0.5652
0.5106
0.1160
RLS
三阶
-0.4381
-0.1228
0.0407
-0.0078
0.5652
0.5106
0.1160
以下给出RLS中的参数估计过程曲线和误差曲线
程序清单：
LS（二阶）：
M=UY（:
1）;
z=UY（:
2）;
H=zeros（199,5）;
fori=1:
199
H（i,1）=-z（i+1）;
H（i,2）=-z（i）;
H（i,3）=M（i+2）;
H（i,4）=M（i+1）;
H（i,5）=M（i）;
end
Estimate=inv（H'*H）*H'*（z（3:
201））
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=[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（'估计方差变化过程'）
同理可以写出三阶的LS以及RLS算法，此处略去。
2.在t=250出加入一个阶跃扰动，令扰动值为0.05
利用RLS和带遗忘因子的RLS计算结果如下表。
算法
遗忘因子
A1
A2
A3
B0
B1
B2
B3
RLS
1
-0.4791
-0.0917
0.0349
-0.0068
0.5659
0.4882
0.1033
RLS
0.98
-0.5498
-0.0390
0.0270
-0.0069
0.5628
0.4515
0.0777
RLS的参数变化过程曲线和误差曲线如下：
带0.98遗忘因子的参数过度曲线和误差曲线
根据上面两种方法所得到的误差曲线和参数过渡过程曲线，我们可以看出来利用最小二乘法得到的参数最终趋于稳定，为利用带遗忘因子的最小二乘算法，曲线参数最终还是有小幅度震荡。
由此可以看出两种算法的一些特点与区别。
最小二乘法：
递推算法没获得一次新的观测数据就修正一次参数估计值，随着时间的推移，便能获得满意的辨识结果。
带遗忘因子的最小二乘法。
其本质还是最小二乘法，只不过加强新的数据提供的信息量，逐渐削弱老的数据，防止数据饱和。
2.参照index3，设计符合GLS和ELS的对象模型，改写参照程序，实现相应的算法。
GLS与ELS的实现
参照index3，设计符合GLS和ELS的对象模型，改写参照程序，实现相应的算法。
利用GLS算法求出参数如下：
a1
a2
b1
b2
-0.7411
0.0811
0.2364
0.1040
即y（k）=0.7411y（k—1）-0.0811y（k-2）+0.2364u（k-1）+0.1040u（k-2）;
利用ELS求解，得到如下的结果
参数
a1
a2
b1
b2
d1
d2
ELS
-0.7424
0.0819
0.2346
0.1039
0
0
绘制出ELS算法下参数的变化曲线以及误差曲线如下：
实验程序：
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,Pstore（3,:
）,i,Pstore（4,:
）,i,Pstore（5,:
）,i,Pstore（6,:
））
title（'估计方差变化过程'）
大作业程序
第一题
x=[1.012.033.024.0156.027.038.049.0310];
y=[9.64.11.30.40.050.10.71.73.89.1];
[ma,na]=size（x）;
sumx=0;sumy=0;sumxy=0;sumx2=0;sumx2y=0;sumx3=0;sumx4=0;
fork=1:
na
sumx=sumx+x（k）;
sumy=sumy+y（k）
sumxy=sumxy+x（k）*y（k）;
sumx2=sumx2+x（k）*x（k）;
sumx2y=sumx2y+x（k）*x（k）*y（k）;
sumx3=sumx3+x（k）^3;
sumx4=sumx4+x（k）^4;
end
A=[na,sumx,sumx2;sumx,sumx2,sumx3;sumx2,sumx3,sumx4];
B=[sumy;sumxy;sumx2y];
C=A\B；
第二题
aa=5;NNPP=7;ts=2;
RR=ones（7）+eye（7）;
UU=[UY（15:
21,1）';UY（14:
20,1）';UY（13:
19,1）';UY（12:
18,1）';UY（11:
17,1）';
UY（10:
16,1）';UY（9:
15,1）'];
YY=[UY（8:
14,2）];
GG=（RR*UU*YY+4.4474）/[aa*aa*（NNPP+1）*ts];
plot（0:
2:
13,GG）;
holdon
stem（0:
2:
13,GG,'filled'）;
gridon;
title（'脉冲响应序列'）
求系统的脉冲传递函数
g=[-0.00360.29220.36830.26730.17050.08200.0489]';
a=[];
fork=1:
1:
6
a=[a;g（k）,g（k+1）];
end
G=[a;g（7）,g
（1）];
G1=[g（3:
7,1）;g
（1）;g
（2）];
A=G\（-G1）;
a1=A
（2）
a2=A
（1）
D=[10;a11];
C=[g
（1）;g
（2）];
B=D*C;
b1=B
（1）
b2=B
（2）
zuixiao
相关二步法
z=UY（:
2）;
M=UY（:
1）;
P=100*eye（4）;
Theta=zeros（4,196）;
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:
196
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;
end
i=1:
196;
plot（i,Theta（1,:
）,i,Theta（2,:
）,i,Theta（3,:
）,i,Theta（4,:
））
title（'相关二步法参数过度过程'）;
gridon
最小二乘递推算法相关参数
clc;
lamt=1;
z=UY（:
2）;
u=UY（:
1）;
c0=[0.0010.0010.0010.001]';
p0=10^4*eye（4,4）;
c=[c0,zeros（4,198）];
fork=3:
200
h1=[-z（k-1）,-z（k-2）,u（k-1）,u（k-2）]';
x=h1'*p0*h1+1*lamt;
x1=inv（x）;
k1=p0*h1*x1;
d1=z（k）-h1'*c0;
c1=c0+k1*d1;
p1=1/lamt*（eye（4）-k1*h1'）*p0;
c（:
k-1）=c1;
c0=c1;
p0=p1;
end
a1=c（1,:
）;a2=c（2,:
）;b1=c（3,:
）;b2=c（4,:
）;
i=1:
199;
plot（i,a1,'r',i,a2,'b',i,b1,'y',i,b2,'black'）;
title（'递推最小二乘法参数辨识'）;
gridon
带遗忘因子的最小二乘法
将上述程序中的lamt=0.9；
辨识系统阶次
方法一、利用残差最小辨识
Data=UY;
L=length（Data）;
Jn=zeros（1,10）;
t=zeros（1,10）;
rm=zeros（10,10）;
forn=1:
1:
10;
N=L-n;
FIA=zeros（N,2*n）;
du=zeros（n,1）;
dy=zeros（n+1,1）;
r1=0;r0=0;
fori=1:
N
forl=1:
n*2
if（l<=n）
FIA（i,l）=-Data（n+i-l,2）;
elseif（n+i-l+2>0）
FIA（i,l）=Data（n+i-l+2,1）;
end
end
end
Y=Data（n+1:
n+N,2）;
thita=inv（FIA'*FIA）*FIA'*Y;
Jn0=0;
fork=n+1:
n+N
forj=1:
n
du（j）=Data（k-j,1）;
dy（j）=Data（k+1-j,2）;
end
dy（n+1）=Data（1,2）;
E1（k）=[1,thita（1:
n）']*dy-thita（n+1:
2*n）'*du;
Jn1=Jn0+E1（k）^2;
t（n）=abs（（Jn0-Jn1）/Jn1*（N-2*n-2）/2）;
Jn0=Jn1;
end
Jn（n）=Jn0;
form=0:
1:
9
form2=n+1:
1:
L-m
r1=r1+E1（m2）*E1（m2+m）/（L-m-n）;
r0=r0+E1（m2）^2/（L-m-n）;
end
rm（n,m+1）=r1/r0;
end
end
subplot（2,1,1）;
plot（1:
10,Jn,'r'）;
title（'残差平方和-jn曲线图'）;
subplot（2,1,2）;
plot（1:
1:
10,t,'g'）;
title（'F检验结果图'）;
figure
（2）;
plot（0:
9,rm（1,:
）,'g'）,holdon;
plot（0:
9,rm（2,:
）,'b'）,holdon;
plot（0:
9,rm（3,:
）,'r'）;
title（'残差白色丁阶结果图'）;
方法二、AIC方法订阶次
Data=UY;
L=length（Data）;
U=Data（:
1）;
Y=Data（:
2）;
forn=1:
1:
5
N=201-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（'AI订阶法'）;
xlabel（'n'）;
ylabel（'AIC'）;
gridon
第三题
广义最小二乘法：
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+1:
n+N）;
phi4=u（n:
n+N-1）;
phi5=u（n-1:
n+N-2）
phi=[phi1phi2phi3phi4phi5];
theta=inv（phi'*phi）*phi'*Y;
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;
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+1:
n+N）;
phi（:
4）=uu（n:
n+N-1）;
phi（:
5）=uu（n-1:
n+N-2）
theta=inv（phi'*phi）*phi'*yy（n+1:
n+N）;
end
f
theta
增广最小二乘法：
程序：
y=UY（:
2）;
u=UY（:
1）;
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
thita=thita（:
800）
第四题
Data=UY;
%使用夏氏修正法，对2阶系统进行参数辨识
n=2;L=length（Data）;N=L-n;
U=Data（:
1）;
Y=Data（:
2）;
%首先使用最小二乘法，求系统参数向量初值
glOL=[-Y（2:
L-1）,-Y（1:
L-2）,U（2:
L-1）,U（1:
L-2）];
Zgl1=Data（3:
L,2）;
Sgl1=glOL'*glOL;Sgl2=inv（Sgl1）;Sgl3=glOL'*Zgl1;
Xsthita=Sgl2*Sgl3;%计算参数矩阵
%得到初值后，以下使用夏氏修正法进行参数辨识：
thitab0=0;%设偏差项的偏差初值为0
Fa=Sgl2*glOL';
M=eye（N）-glOL*Sgl2*glOL';
F=eye（N）;
i=0;
if（F>=10^-6*eye（N））
E1=Zgl1-glOL*Xsthita;%计算残差E
omiga（2:
N,1）=-E1（1:
N-1）;
omiga（3:
N,2）=-E1（1:
N-2）;
D=omiga'*M*omiga;
Fx=inv（D）*omiga'*M*Zgl1;
thitab=Fa*omiga*Fx;
Xsthita=Xsthita-thitab;
F=thitab-thitab0;
thitab0=thitab;
end
disp（'夏氏修正法'）
Xsa1=Xsthita
（1）,Xsa2=Xsthita（2
升级会员 