系统辨识实验报告Word下载.docx
《系统辨识实验报告Word下载.docx》由会员分享，可在线阅读，更多相关《系统辨识实验报告Word下载.docx（26页珍藏版）》请在冰豆网上搜索。
RLS
-0.7824
三阶
-0.4381
-0.1228
0.0407
-0.0078
0.5652
0.5106
0.1160
以下给出RLS中的参数估计过程曲线和误差曲线
程序清单：
LS（二阶）：
M=UY（:
1）;
z=UY（:
2）;
H=zeros（,5）;
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）;
Estimate=inv（H'
*H）*H'
*（z（3:
201））
RLS（二阶）：
clc
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];
K=[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）;
i-1）=Theta（:
i-2）+K*（z（i）-h'
*Theta（:
i-2））;
P=（eye（5）-K*h'
）*P;
i-1）=[P（1,1）,P（2,2）,P（3,3）,P（4,4）,P（5,5）]'
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,:
估计方差变化过程'
同理可以写出三阶的LS以及RLS算法，此处略去。
2.在t=250出加入一个阶跃扰动，令扰动值为0.05
利用RLS和带遗忘因子的RLS计算结果如下表。
遗忘因子
1
-0.4791
-0.0917
0.0349
-0.0068
0.5659
0.4882
0.1033
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求解，得到如下的结果
参数
d1
d2
ELS
-0.7424
0.0819
0.2346
0.1039
绘制出ELS算法下参数的变化曲线以及误差曲线如下：
实验程序：
GLS
%%%%%%%%%%%%%%%%广义最小二乘法辨识%%%%%%%%%%%
n=2;
N=799;
%y=y1;
%u=u1;
y=UY（:
u=UY（:
%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:
phi4=u（n-1:
phi=[phi1phi2phi3phi4];
theta=inv（phi'
*phi）*phi'
*Y
%%%%%%以上为最小二乘计算Theta估计值%%%%%%%%%
f0=[10;
while1
e（n+1:
n+N,1）=y（n+1:
n+N）-phi*theta;
omega（:
1）=-e（n:
2）=-e（n-1:
f=inv（omega'
*omega）*omega'
*e（n+1:
n+N,1）;
if（f-f0）'
*（f-f0）<
0.0001
break
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）]'
yy（1,1）=y（1,1）;
uu（2,1）=u（2,1）;
%%%%%%%%%%%以上为数据滤波%%%%%%%%
phi（:
1）=-yy（n:
2）=-yy（n-1:
3）=uu（n:
4）=uu（n-1:
*yy（n+1:
f
theta
%增广最小二乘的递推算法
%===========产生均值为0，方差为1的高斯白噪声=========
v
（1）=0;
v
（2）=0;
%递推求解
P=100*eye（6）;
Pstore=zeros（6,800）;
1）=[P（1,1）,P（2,2）,P（3,3）,P（4,4）,P（5,5）,P（6,6）];
Theta=zeros（6,401）;
1）=[3;
3;
3];
801
M（i-2）;
v（i-1）;
v（i-2）];
v（i）=z（i）-h'
i-2）;
P=（eye（6）-K*h'
i-1）=[P（1,1）,P（2,2）,P（3,3）,P（4,4）,P（5,5）,P（6,6）];
%=======================================================================
disp（'
参数a1、a2、b1、b2、d1、d2估计结果：
'
800）
800;
）,i,Theta（6,:
）,i,Pstore（6,:
大作业程序
第一题
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;
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:
13,GG,'
filled'
）;
gridon;
脉冲响应序列'
求系统的脉冲传递函数
g=[-0.00360.29220.36830.26730.17050.08200.0489]'
a=[];
1:
6
a=[a;
g（k）,g（k+1）];
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）;
B=D*C;
b1=B
（1）
b2=B
（2）
zuixiao
相关二步法
P=100*eye（4）;
Theta=zeros（4,）;
2）=[3;
3）=[3;
4）=[3;
fori=5:
hstar=[M（i-1）;
M（i-3）;
M（i-4）];
K=P*hstar*inv（h'
*P*hstar+1）;
i）=Theta（:
i-1）+K*（z（i）-h'
i-1））;
P=（eye（4）-K*h'
相关二步法参数过度过程'
gridon
最小二乘递推算法相关参数
clc;
lamt=1;
c0=[0.0010.0010.0010.001]'
p0=10^4*eye（4,4）;
c=[c0,zeros（4,）];
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,:
plot（i,a1,'
r'
i,a2,'
b'
i,b1,'
y'
i,b2,'
black'
）;
递推最小二乘法参数辨识'
带遗忘因子的最小二乘法
将上述程序中的lamt=0.9；
辨识系统阶次
方法一、利用残差最小辨识
Data=UY;
L=length（Data）;
Jn=zeros（1,10）;
t=zeros（1,10）;
rm=zeros（10,10）;
forn=1:
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）;
Y=Data（n+1:
n+N,2）;
thita=inv（FIA'
*FIA）*FIA'
*Y;
Jn0=0;
fork=n+1:
forj=1:
n
du（j）=Data（k-j,1）;
dy（j）=Data（k+1-j,2）;
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;
Jn（n）=Jn0;
form=0:
9
form2=n+1:
L-m
r1=r1+E1（m2）*E1（m2+m）/（L-m-n）;
r0=r0+E1（m2）^2/（L-m-n）;
rm（n,m+1）=r1/r0;
subplot（2,1,1）;
plot（1:
10,Jn,'
残差平方和-jn曲线图'
subplot（2,1,2）;
10,t,'
g'
F检验结果图'
figure
（2）;
9,rm（1,:
）,'
）,holdon;
9,rm（2,:
9,rm（3,:
残差白色丁阶结果图'
方法二、AIC方法订阶次
Data=UY;
U=Data（:
Y=Data（:
forn=1:
5
N=201-n
yd（1:
N,1）=Y（n+1:
fori=1:
forl=1:
2*n
FIA（i,l）=-Y（n+i-l）;
else
FIA（i,l）=U（2*n+i-l）;
thita=inv（FIA'
*yd;
omiga=（yd-FIA*thita）'
*（yd-FIA*thita）/N;
AIC（n）=N*log（omiga）+4*n;
plot（AIC,'
-'
AI订阶法'
xlabel（'
n'
ylabel（'
AIC'
第三题
广义最小二乘法：
phi3=u（n+1:
phi4=u（n:
phi5=u（n-1:
n+N-2）
phi=[phi1phi2phi3phi4phi5];
*Y;
3）=uu（n+1:
4）=uu（n:
5）=uu（n-1:
增广最小二乘法：
程序：
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）;
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;
thita=thita（:
第四题
%使用夏氏修正法，对2阶系统进行参数辨识
n=2;
N=L-n;
%首先使用最小二乘法，求系统参数向量初值
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;
夏氏修正法'
Xsa1=Xsthita
（1）,Xsa2=Xsthita
（2）,Xsb1=Xsthita（3）,Xsb2
升级会员 