多变量系统辨识matlab程序

多变量系统辨识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')
hold on
xi=y';
save sub.txt xi –ascii
程序5
clear
%CRA模型基于模型阶次递增的辨识。

clc
close all
z=load('sub.txt');
u1=load('prbs1.txt');
u2=load('prbs2.txt');
for i=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);
for i=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);
%**************************
for n=2:1:10
for i=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;
for i=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')
hold on
e=z-y';
save error.txt e –ascii
程序7
clear
%子模型基于模型阶次递增的辨识。

clc
close all
z=load('sub.txt');
u1=load('prbs1.txt');
u2=load('prbs2.txt');
e=load('error.txt');
for i=1:1:100
end
theta=(1e-3)*ones(4,1);
P=(1e8)*eye(4);
for i=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);
%**************************
for n=2:1:10
for i=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=[H1 H2];
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;
for i=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')
hold on
save ssub1.txt y1 -ascii
save ssub2.txt y2 –ascii
程序9
clear
%子子模型1的基于阶次递增的辨识。

clc
% close all
z=load('ssub1.txt');
u=load('prbs1.txt');
for i=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);
for i=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;
%**************************
for n=2:1:10
for i=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=[H1 H2];
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
close all
z=load('ssub2.txt');
u=load('prbs2.txt');
for i=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);
for i=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。