多变量系统辨识matlab程序.docx
《多变量系统辨识matlab程序.docx》由会员分享，可在线阅读，更多相关《多变量系统辨识matlab程序.docx（9页珍藏版）》请在冰豆网上搜索。
多变量系统辨识matlab程序
多变量系统辨识matlab程序
y（i）=0.05129*u1_1+0.0418；u1_3=u1_2;u1_2=u1_1;u1_1；u2_3=u2_2;u2_2=u2_1;u2_1；y_3=y_2;y_2=y_1;y_1=y（i）；r_3=r_2;r_2=r_1;r_1=r（i）；end；plot（time,y,'b'）；holdon；xi=y';；savesub.
y（i）=0.05129*u1_1+0.0418*u2_1+0.6386*y_1+0.06268*u1_2+0.0346*u2_2-0.1179*y_2-0.004184*u1_3-0.00218*u2_3+0.006738*y_3+0.091*r_1-0.114*r_2+0.0509*r_3;
u1_3=u1_2;u1_2=u1_1;u1_1=u1（i）;
u2_3=u2_2;u2_2=u2_1;u2_1=u2（i）;
y_3=y_2;y_2=y_1;y_1=y（i）;
r_3=r_2;r_2=r_1;r_1=r（i）;
end
plot（time,y,'b'）
holdon
xi=y';
savesub.txtxi–ascii
程序5
clear
%CRA模型基于模型阶次递增的辨识。
clc
closeall
z=load（'sub.txt'）;
u1=load（'prbs1.txt'）;
u2=load（'prbs2.txt'）;
fori=1:
1:
100
H（i,:
）=[u1（20+i-1）u2（20+i-1）-1*z（20+i-1）];
end
theta=（1e-3）*ones（3,1）;
P=（1e8）*eye（3）;
fori=1:
1:
100
K=P*H（i,:
）'./（H（i,:
）*P*H（i,:
）'+1）;
theta=theta+K*（z（i+20）-H（i,:
）*theta）;
P=（eye（3）-K*H（i,:
））*P;
end
theta1=theta
H1=H;
J
（1）=（z（21:
120）-H1*theta1）'*（z（21:
120）-H1*theta1）;
ZZ=inv（H1'*H1）;
%**************************
forn=2:
1:
10
fori=1:
1:
100
H2（i,:
）=[u1（20+i-n）u2（20+i-n）-1*z（20+i-n）];
end
B=inv（H2'*H2-H2'*H1*ZZ*H1'*H2）;
A=ZZ*H1'*H2*B;
theta2=B*H2'*（z（21:
120）-H1*theta1）;
theta1=theta1-A*H2'*（z（21:
120）-H1*theta1）;
theta1=[theta1;theta2]
ZZ1=[ZZ+A*H2'*H1*ZZ-A];
ZZ2=[-A'B];
ZZ=[ZZ1;ZZ2];
J（n）=（z（21:
120）-H1*theta1）'*（z（21:
120）-H1*theta1）;
F（n-1）=（（J（n-1）-J（n））/2）/（（J（n））/（100-2*n））;
time（n-1）=n;
TEST（n-1）=3;
end
plot（time,F,'r-*',time,TEST）
title（'F统计值随系统阶次的变化'）
xlabel（'系统阶次'）
ylabel（'F统计值'）
legend（'F（2（n_2-n_1）,100-2n_2）','F（2,100）'）
程序6
clear
%****************CAR模型最佳辨识的验证，同时获取CARMA模型的残差序列，存于error.txt中。
clc
u1=load（'prbs1.txt'）;
u2=load（'prbs2.txt'）;
z=load（'sub.txt'）
u1_6=0;u1_5=0;u1_4=0;u1_1=0;u1_2=0;u1_3=0;
u2_6=0;u2_5=0;u2_4=0;u2_1=0;u2_2=0;u2_3=0;
y_6=0;y_5=0;y_4=0;y_1=0;y_2=0;y_3=0;
r_1=0;r_2=0;r_3=0;
fori=1:
1:
120
time（i）=i;
y（i）=0.0496*u1_1+0.0417*u2_1-0.6724*y_1+0.1300*u1_2+0.0902*u2_2-0.4219*y_2+0.1352*u1_3+0.0911*u2_3-0.1887*y_3+0.1032*u1_4+0.0707*u2_4+-0.0188*y_4+0.0639*u1_5+0.0401*u2_5+00.1125*y_5+0.0210*u1_6+0.0132*u2_6-0.0101*y_6;
u1_6=u1_5;u1_5=u1_4;u1_4=u1_3;u1_3=u1_2;u1_2=u1_1;u1_1=u1（i）;
u2_6=u2_5;u2_5=u2_4;u2_4=u2_3;u2_3=u2_2;u2_2=u2_1;u2_1=u2（i）;
y_6=y_5;y_5=y_4;y_4=y_3;y_3=y_2;y_2=y_1;y_1=y（i）;
end
plot（time,y,'g'）
holdon
e=z-y';
saveerror.txte–ascii
程序7
clear
%子模型基于模型阶次递增的辨识。
clc
closeall
z=load（'sub.txt'）;
u1=load（'prbs1.txt'）;
u2=load（'prbs2.txt'）;
e=load（'error.txt'）;
fori=1:
1:
100
end
theta=（1e-3）*ones（4,1）;
P=（1e8）*eye（4）;
fori=1:
1:
100
K=P*H（i,:
）'./（H（i,:
）*P*H（i,:
）'+1）;
theta=theta+K*（z（i+20）-H（i,:
）*theta）;
P=（eye（4）-K*H（i,:
））*P;
end
theta1=theta
H1=H;
J
（1）=（z（21:
120）-H1*theta1）'*（z（21:
120）-H1*theta1）;
ZZ=inv（H1'*H1）;
%**************************
forn=2:
1:
10
fori=1:
1:
100
H2（i,:
）=[u1（20+i-n）u2（20+i-n）-1*z（20+i-n）e（20+i-n）];
end
B=inv（H2'*H2-H2'*H1*ZZ*H1'*H2）;
A=ZZ*H1'*H2*B;
theta2=B*H2'*（z（21:
120）-H1*theta1）;
theta1=theta1-A*H2'*（z（21:
120）-H1*theta1）;
theta1=[theta1;theta2]
ZZ1=[ZZ+A*H2'*H1*ZZ-A];
ZZ2=[-A'B];
ZZ=[ZZ1;ZZ2];
H1=[H1H2];
J（n）=（z（21:
120）-H1*theta1）'*（z（21:
120）-H1*theta1）;
F（n-1）=（（J（n-1）-J（n））/3）/（（J（n））/（100-3*n））;
time（n-1）=n;
TEST（n-1）=3;
end
plot（time,F,'r-*',time,TEST）
title（'F统计值随系统阶次的变化'）
xlabel（'系统阶次'）
ylabel（'F统计值'）
legend（'F（2（n_2-n_1）,100-2n_2）','F（2,100）'）
程序8
clear
%****************CARMA模型最佳辨识的验证，同时获取子子模型的输出，分别存于ssub1.txt与ssub2.txt中。
clc
u1=load（'prbs1.txt'）;
u2=load（'prbs2.txt'）;
z=load（'sub.txt'）;
e=load（'error.txt'）;
u1_1=0;u1_2=0;u1_3=0;u1_4=0;
y_1=0;y_2=0;y_3=0;y_4=0;
e_4=0;e_1=0;e_2=0;e_3=0;
y1_1=0;y1_2=0;y1_3=0;y1_4=0;
y2_1=0;y2_2=0;y2_3=0;y2_4=0;
r_1=0;r_2=0;r_3=0;
fori=1:
1:
120
time（i）=i;
y（i）=0.0507*u1_1+0.0426*u2_1-0.8329*y_1+0.2335*e_1+0.1383*u1_2+0.0974*u2_2-0.2791*y_2-0.1006*e_2+0.1437*u1_3+0.0956*u2_3+0.4623*y_3-0.7050*e_3+0.0654*u1_4+0.0400*u2_4-0.0736*y_4-0.0502*e_4;
y1（i）=0.0507*u1_1-0.8329*y1_1+0.1383*u1_2-0.2791*y1_2+0.1437*u1_3+0.4623*y1_3+0.0654*u1_4-0.0736*y1_4;
y2（i）=0.0426*u2_1-0.8329*y2_1+0.0974*u2_2-0.2791*y2_2+0.0956*u2_3+0.4623*y2_3+0.0400*u2_4-0.0736*y2_4;
u1_4=u1_3;u1_3=u1_2;u1_2=u1_1;u1_1=u1（i）;
u2_4=u2_3;u2_3=u2_2;u2_2=u2_1;u2_1=u2（i）;
y_4=y_3;y_3=y_2;y_2=y_1;y_1=y（i）;
e_4=e_3;e_3=e_2;e_2=e_1;e_1=e（i）;
y1_4=y1_3;y1_3=y1_2;y1_2=y1_1;y1_1=y1（i）;
y2_4=y2_3;y2_3=y2_2;y2_2=y2_1;y2_1=y2（i）;
end
plot（time,y,'m',time,y1+y2,'k'）
holdon
savessub1.txty1-ascii
savessub2.txty2–ascii
程序9
clear
%子子模型1的基于阶次递增的辨识。
clc
%closeall
z=load（'ssub1.txt'）;
u=load（'prbs1.txt'）;
fori=1:
1:
100
H（i,:
）=[u（20+i-1）-1*z（20+i-1）];
end
theta=（1e-3）*ones（2,1）;
P=（1e8）*eye
（2）;
fori=1:
1:
100
K=P*H（i,:
）'./（H（i,:
）*P*H（i,:
）'+1）;
theta=theta+K*（z（i+20）-H（i,:
）*theta）;
P=（eye
（2）-K*H（i,:
））*P;
end
theta1=theta
sigmaF
（1）=（z（21:
120）'-H1*theta1）'*（z（21:
120）'-H1*theta1）/100;AIC
（1）=100*log10（sigmaF
（1））+2*（1+1）;
ZZ=inv（H1'*H1）;
ZZ=inv（H1'*H1）;
THETA=zeros（10,20）;
%THETA（2,:
）=theta1;
%**************************
forn=2:
1:
10
fori=1:
1:
100
H2（i,:
）=[u（20+i-n）-1*z（20+i-n）];
end
B=inv（H2'*H2-H2'*H1*ZZ*H1'*H2）;
A=ZZ*H1'*H2*B;
theta2=B*H2'*（z（21:
120）'-H1*theta1）;
theta1=theta1-A*H2'*（z（21:
120）'-H1*theta1）;
theta1=[theta1;theta2]
ZZ1=[ZZ+A*H2'*H1*ZZ-A];
ZZ2=[-A'B];
ZZ=[ZZ1;ZZ2];
H1=[H1H2];
time（n-1）=n;
sigmaF（n）=（z（21:
120）'-H1*theta1）'*（z（21:
120）'-H1*theta1）/100;AIC（n-1）=100*log10（sigmaF（n））+2*（n+n）;
end
plot（time,AIC,'r-*'）
title（'信息准则AIC值随系统阶次的变化'）
xlabel（'系统阶次'）
ylabel（'信息准则AIC值'）
程序10
clear
%子子模型2的基于阶次递增的辨识。
clc
closeall
z=load（'ssub2.txt'）;
u=load（'prbs2.txt'）;
fori=1:
1:
100
H（i,:
）=[u（20+i-1）-1*z（20+i-1）];
end
theta=（1e-3）*ones（2,1）;
P=（1e8）*eye
（2）;
fori=1:
1:
100
K=P*H（i,:
）'./（H（i,:
）*P*H（i,:
）'+1）;
theta=theta+K*（z（i+20）-H（i,:
）*theta）;
P=（eye
（2）-K*H（i,:
））*P;
end
theta1=theta
升级会员 