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')

```