www.pudn.com > zishiying.rar > zishiying.m, change:2010-07-06,size:3690b


clear  
%disp('单变量系统的测控制算法'); 
%disp('广义预测控制算法的初始植'); 
n=4;%input('预测长度n='); 
m=3;%input('控制长度m='); 
t0=0.5;%input('控制加权='); 
a=0.4;%input('柔化系数a='); 
p=0.98; 
A=[1;-1.5;0.7];B=[1;1.5];A1=[2.5,-2.2,0.7];B1=[1.5]; 
%参数初始化 
T=50; 
na=length(A)-1; 
nb=length(B)-1; 
A2=zeros(na+1,1); 
A2(1,1)=1; 
B2=zeros(nb+1,1); 
Q0=zeros(na+nb+1,1); 
for l=1:1:na+nb+1 
    if(l<na+1) 
     Q0(l,1)=A(l+1,1); 
    end 
   if (l>na)    
   Q0(l,1)=B(l-na,1); 
   end 
end 
Q1=zeros(na+nb+1,1); 
Q2=zeros(na+nb+1,1); 
P1=(1e+5)*eye(na+nb+1); 
P2=(1e+5)*eye(na+nb+1); 
uu=zeros(n,1); 
yy=zeros(n,1); 
y1=zeros(na+1,1); 
u1=zeros(nb,1); 
G=zeros(n,m); 
u=zeros(m,1); 
t=zeros(T,1); 
t(1,1)=1; 
t(2,1)=2; 
uuu=zeros(T,1); 
yyy=zeros(T,1); 
uuu(1,1)=0; 
uuu(2,1)=0; 
yyy(1,1)=0; 
yyy(2,1)=0; 
yyr=zeros(T,1); 
yyr(1,1)=1.5; 
yyr(2,1)=1.5; 
yd=zeros(n,1); 
 dy=zeros(na,1); 
 ddu=zeros(nb+1,1); 
%求F,H,G,E, 
% 开始循环 
for ij=3:T 
     t(ij,1)=ij; 
for s=1:1:na 
    if(ij-s-1>0) 
        dy(s,1)=-(yyy(ij-s,1)-yyy(ij-s-1,1)); 
    end 
end 
for s1=1:1:nb+1 
    if(ij-t-1>0) 
        ddu(s1,1)=uuu(ij-s1,1)-uuu(ij-s1-1,1); 
    end 
end 
%DDU=ddu 
X=[dy;ddu]; 
ddy=X'*Q0; 
ss=ddy-X'*Q1; 
Q2=Q1+(P1*X*ss)/(p+(X'*P1*X)); 
P2=(P1-(P1*X*X'*P1)/(p+(X'*P1*X)))/p; 
Q1=Q2; 
%q2=Q2 
P1=P2; 
for f=2:1:na+1 
    A2(f,1)=Q1(f-1,1); 
end 
for f1=1:1:nb+1 
    B2(f1,1)=Q1(f1+na,1); 
end 
if(nb==1) 
    B1=B2(2,1); 
end 
if(nb>1) 
    for k=1:1:nb 
    B1(1,k)=B2(k+1,1); 
    end 
end 
%b1=B1 
for m1=1:1:na 
A1(1,m1)=A2(m1,1)-A2(m1+1,1); 
end 
A1(1,na+1)=A2(na+1,1); 
%a1=A1 
F1=zeros(1,na+1); 
F2=zeros(1,na+1); 
H1=zeros(1,nb); 
H2=zeros(1,nb); 
e=1; 
E=e; 
g=E(1,1)*B(1,1); 
G1=g; 
for i=1:na+1 
    F1(1,i)=A1(1,i); 
end 
if (nb==1) 
    H1(1,1)=B1(1,1); 
end 
if (nb>1) 
    for j=1:nb 
    H1(1,j)=B1(1,j); 
    end 
end 
F=F1; 
H=H1; 
for kh=2:n 
    e=F1(1,1); 
    E=[E;e]; 
    g=e*B(1,1)+H1(1,1); 
    G1=[G1;g]; 
    for k=1:na 
        F2(1,k)=F1(1,k+1)+A1(1,k)*e; 
    end 
   F2(1,na+1) =A1(1,na+1)*e; 
   F1=F2; 
   F=[F;F2]; 
   if  (nb==1) 
       H2(1,1)=B1(1,1)*e; 
   end 
   if(nb==2) 
       H2(1,1)=H1(1,2)+B1(1,1)*e; 
   end 
   if(nb>2) 
    for h=1:nb-1 
        H2(1,h)=H1(1,h+1)+B1(1,h)*e; 
    end 
   H2(1,nb) =B1(1,nb)*e; 
   end 
   H=[H;H2]; 
   H1=H2; 
end 
for s2=1:m     %G=zeros(n,m); 
     for l1=n-s2+1:-1:1 
        il=l1+s2-1; 
        G(il,s2)=G1(l1,1); 
    end 
end 
  yr=1.5;a9=0; 
    for i=1:1 
        a9=a9+rand;%noises 
    end 
    a8=0.04*a9; 
    yyr(ij,1)=yr+a8;  
  y=1.5*yy(n,1)-0.7*yy(n-1,1)+uu(n,1)+1.5*uu(n-1,1); 
  yyy(ij,1)=y; 
    
    for p1=na+1:-1:2 
         y1(p1)=y1(p1-1); 
     end 
      y1(1,1)=y; 
         
     for v=1:1:n-1 
          yy(v,1)= yy(v+1,1); 
      end 
      yy(n,1)=y; 
      if(nb==1) 
        u1(1,1)= uuu(ij-1,1)- uuu(ij-2,1);   
      end 
      if(nb==2) 
          u1(1,1)=u1(2,1);  
           u1(2,1)= uuu(ij-1,1)- uuu(ij-2,1);   
      end 
      if(nb>2) 
     for q=nb:-1:2 
         u1(q,1)=u1(q-1,1); 
     end 
     u1(1,1)= uuu(ij-1,1)- uuu(ij-2,1);    
      end 
   y0=F*y1+H*u1; 
   yd(1,1)=a*y+(1-a)*yr; 
   for i=2:n 
       yd(i)=a^i*yd(i-1,1)+(1-a^i)*yr; 
   end 
   u=inv(G'*G+t0*eye(m))*G'*(yd-y0);%只是取第一个控制增量 
 if(u(1,1)>0.2) 
      u(1,1)=0.2; 
 end 
  if(u(1,1)<-0.2) 
       u(1,1)=-0.2; 
  end  
   for j=1:n-1 
         uu(j,1)= uu(j+1,1); 
  end 
    uu(n,1)= uuu(ij-1,1)+u(1,1); 
     uuu(ij)= uu(n,1);   
     ddu(2,1)=ddu(1,1); 
     ddu(1,1)=u(1,1); 
end 
    %plot(t,yyr,t,yyy,'r'); 
      subplot(2,1,1):plot(t,yyr,t,yyy,'r'); 
    xlabel('t');ylabel('yr,y') 
      subplot(2,1,2):plot(t,uuu); 
    xlabel('t');ylabel('u')