original.m

用m文件实现的广义预测基本算法

M=4;P=6;%N1=1;N2=P;
A=[1 0.36 0.24];
A1=[1 -0.64 -0.12 -0.24];
B=[0.5 0.25 ];
% A=[1 -0.496585];
% A1=[1 -1.496585 0.496585];
% B=[0 0.5 -0.5];
la=length(A);
lb=length(B);
na=la-1;nb=lb-1;
%S(1,:)=[2.2 -1.7 0.5]; 
S=zeros(P,na+1);
S(1,:)=[0.64 0.12 0.24];
% S(1,:)=[1.496585 -0.496585];
for j=1:P-1 
    for i=1:na  
       S(j+1,i)=S(j,i+1)-A1(i+1)*S(j,1); 
    end
    S(j+1,na+1)=-S(j,1)*A1(na+2);
end
% R=zeros(P,P);
% R(1,1)=1;
% for j=1:P-1
%     for i=1:j
%        R(j+1,i)=R(j,i);
%     end
%     R(j+1,j+1)=S(j,1);
% end
G=zeros(P,P+nb);
for j=1:nb+1 
G(1,j)=B(j);
end
for j=1:P-1
    for i=1:j+nb
         if i<j
           bij=0;
         else 
           bij=B(i-j+1);
         end
         G(j+1,i+1)=G(j,i+1)+S(j,1)*bij;
     end
end
for j=2:P
    G(j,1)=B(1);
end
% G1=zeros(P,M);
% for j=1:P
%     for i=j:M
%         G1(i,j)=G(i-j+1,i-j+1);
%     end
% end
G1=zeros(P,P);
for i=1:P
    for j=1:i
        G1(i,j)=G(i-j+1,i-j+1);
    end
end
G11=zeros(P,M);
for i=1:P
    for j=1:M
        G11(i,j)=G1(i,j);
    end
end
G2=zeros(P,nb);
for i=1:P
    for j=1:nb
        G2(i,j)=G(i,nb+1+i-j);
    end
end
u1=zeros(P,1); 
u0=zeros(nb,1);
y0=zeros(na+1,1);
y=[];
u01=zeros(na+2,1);
u=[];
for j=1:na+1
u(j)=0;
end
u=u';
yr=[];yr(1)=0;
yr1=zeros(P,1);
sy=zeros(P,1);
d1=zeros(1,M);
d1(1,1)=1;%
f=zeros(P,1);
% w=2;r=2.4;a=0.8;
w=1;r=0.5;a=0.1;
for t=1:200 
for j=1:P
yr(t+j)=a*yr(t+j-1)+(1-a)*w;
yr1(j)=yr(t+j);
end    
for j=1:P
    syy=0;
    for i=1:na+1
    syy=syy+S(j,i)*y0(i);
    end
    sy(j)=syy;
end
f=G2*u0+sy;
lamda=r*eye(M);
lamda(1,1)=1;
u1=inv(G11'*eye(P)*G11+lamda)*G11'*eye(P)*(yr1-f);
ut1=d1*u1;
ut=ut1+u01(na+2);
for i=1:na+1
    u01(i)=u01(i+1);
end
u01(na+2)=ut;
u=[u;ut];
for i=1:nb-1
u0(i)=u0(i+1);
end
u0(nb)=ut1; 
if t<30 
    yyy1=0;
    yyy2=0;  
    for i=1:lb
        yyy1=yyy1+B(i)*u(nb+1+t-i);
    end
    for i=2:la
        yyy2=yyy2+A(i)*y0(i);
    end
    y1=yyy1-yyy2;
else
    yyy1=0;
    yyy2=0; 
    for i=1:lb
        yyy1=yyy1+B(i)*u(nb+1+t-i);
    end
    for i=2:la
        yyy2=yyy2+A(i)*y0(i)-0.1;
    end
    y1=yyy1-yyy2;
end
y=[y;y1];
for i=la:-1:2
    y0(i)=y0(i-1);
end
y0(1)=y1; 
end 
t=1:200;
figure;
plot(t,y);
grid on