系统辨识大作业.docx
《系统辨识大作业.docx》由会员分享，可在线阅读，更多相关《系统辨识大作业.docx（20页珍藏版）》请在冰豆网上搜索。
系统辨识大作业
系统辨识大作业
第一题
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；
利用最小二乘的求得
作图验证如下：
第二题，搭建的模型如图（m序列发生器采用三个移位寄存器，其他参数同实验二）
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=zeros（4,400）;%ÔöÒæ¾ØÕó
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
得到相关系数为ans=
-0.7982
0.1377
0.5631
0.3149
最小二乘递推算法相关参数
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;%¸üÐÂC1
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
得到相关系数为：
ans=
-0.7971
0.1423
0.5616
0.3127
带遗忘因子的最小二乘法
将上述程序中的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;
%°´AkkaikeÐÅÏ¢×¼Ôò¶¨½×
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=800;%nΪÂö³åÊýÄ¿
M=UY（1:
800,1）;
%===========²úÉú¾ùֵΪ0£¬·½²îΪ1µÄ¸ß˹°×ÔëÉù=========
v=randn（1,800）;
e=[];
e
（1）=v
（1）;
e
（2）=v
（2）;
fori=3:
800
e（i）=0*e（i-1）+0*e（i-2）+v（i）;
end
%==============²úÉú¹Û²âÐòÁÐz=================
z=UY（1:
800,2）;
%±ä»»ºóµÄ¹Û²âÐòÁÐ
zf=[];
zf
（1）=-1;
zf
（2）=0;
fori=3:
800
zf（i）=z（i）-0*z（i-1）-0*z（i-2）;
end
%±ä»»ºóµÄÊäÈëÐòÁÐ
uf=[];
uf
（1）=M
（1）;
uf
（2）=M
（2）;
fori=3:
800
uf（i）=M（i）-0*M（i-1）-0*M（i-2）;
end
%¸³³õÖµ
P=100*eye（4）;%¹À¼Æ·½²î
Theta=zeros（4,800）;%²ÎÊýµÄ¹À¼ÆÖµ£¬´æ·ÅÖмä¹ý³Ì¹ÀÖµ
Theta（:
2）=[3;3;3;3];
K=[10;10;10;10];%ÔöÒæ
PE=10*eye
（2）;
ThetaE=zeros（2,800）;
ThetaE（:
2）=[0.5;0.3];
KE=[10;10];
fori=3:
800
h=[-zf（i-1）;-zf（i-2）;uf（i-1）;uf（i-2）];
K=P*h*inv（h'*P*h+1）;
Theta（:
i）=Theta（:
i-1）+K*（z（i）-h'*Theta（:
i-1））;
P=（eye（4）-K*h'）*P;
he=[-e（i-1）;-e（i-2）];
KE=PE*he*inv（1+he'*PE*he）;
ThetaE（:
i）=ThetaE（:
i-1）+KE*（e（i）-he'*ThetaE（:
i-1））;
PE=（eye
（2）-KE*he'）*PE;
end
%=====================Êä³ö½á¹û¼°×÷ͼ=========================
i=1:
800;
figure
（1）;
plot（i,Theta（1,:
）,i,Theta（2,:
）,i,Theta（3,:
）,i,Theta（4,:
））
title（'´ý¹À²ÎÊý¹ý¶É¹ý³Ì'）
figure（3）;
plot（i,ThetaE（1,:
）,
升级会员 