这里是N4SID的,想改成pem的算法,然后对比两种算法
clear all;
N=1000;
u1=idinput(N,'rgs',[0 1]);
u2=idinput(N,'rgs',[0 1]);
u3=idinput(N,'rgs',[0 1]);
y1=zeros(N,1);y2=zeros(N,1);y3=zeros(N,1);
m3=[u1,u2,u3];
e = randn(1000,3).*0.001;
A=[0 1 0;0 0 1;-0.1 -0.6 1.5];
B=[0.84 0.42 0.21;-0.68 0.22 -0.15;-0.66 0.3 -0.15];
C=[1 0 0;0 1 0; 0 0 1];
D=[0 0 0;0 0 0; 0 0 0];
K =[1 0 0;0 1 0; 0 0 1];
x0 = [0;0;0];
th0=idss(A,B,C,D,K,x0,'ts',0.4);
m0.As = [0 NaN 0;0 0 NaN;NaN NaN NaN];
y=sim(th0,[m3 e]);
y1=y(:,1);
y2=y(:,2);
y3=y(:,3);
z = iddata([y1 y2 y3],[ u1 u2 u3]);
m = n4sid(z,[1:10],'foc','sim','ssp','can','ts',0.4);
ffplot(m,'r-',th0,'b-.');
clear all;
N=1000;
u1=idinput(N,'rgs',[0 1]);
u2=idinput(N,'rgs',[0 1]);
u3=idinput(N,'rgs',[0 1]);
y1=zeros(N,1);y2=zeros(N,1);y3=zeros(N,1);
m3=[u1,u2,u3];
e = randn(1000,3).*0.001;
A=[0 1 0;0 0 1;-0.1 -0.6 1.5];
B=[0.84 0.42 0.21;-0.68 0.22 -0.15;-0.66 0.3 -0.15];
C=[1 0 0;0 1 0; 0 0 1];
D=[0 0 0;0 0 0; 0 0 0];
K =[1 0 0;0 1 0; 0 0 1];
x0 = [0;0;0];
th0=idss(A,B,C,D,K,x0,'ts',0.4);
m0.As = [0 NaN 0;0 0 NaN;NaN NaN NaN];
y=sim(th0,[m3 e]);
y1=y(:,1);
y2=y(:,2);
y3=y(:,3);
z = iddata([y1 y2 y3],[ u1 u2 u3]);
m = n4sid(z,[1:10],'foc','sim','ssp','can','ts',0.4);
ffplot(m,'r-',th0,'b-.');