yinshigpc.m

yinshiGPC model predictive control

clear
nn=6;   %nn=input('时域长度nn=');
n=6;  %n=input('预测长度n=');
m=2 ;  %m=input('控制长度m=');
t0=0.8;   %t0=input('控制加权系数=');
a=0.3;     %a=input('柔化系数a=');
disp('最小二乘公式初始值')
t1=1;  %t1=1;%input('遗忘因子')
P=(1e+5)*eye(n+1);   
%参数初始值
uuu=0;yyy=0;
uu=zeros(n,1);u=zeros(m,1);
yy=zeros(n,1);y1=zeros(n,1);
Q=zeros(n+1,1);Q(1,1)=1;Q(n+1,1)=1;
%产生周期为100，时间为T，幅值为1的方波信号的给定值
T=300;[yr0,t]=gensig('square',100,T,1);
nm=length(t);%确定循环次数
for ij=2:nm;
yr=yr0(ij)+1;%产生周期为100，时间为T，幅值在1和2之间变化的防拨信号的给定值
a9=0;
for i=1:1
    a9=a9+rand;
end
a8=0.01*(a9-6);
%根据系统模型，计算k时刻的输出值y(k)
y=1.496585*yy(n,1)-0.496585*yy(n-1,1)+0.5*uu(n-1,1);
%(nn=6;n=6;m=2;t0=0.8;a=0.3;t1=1)
%产生均匀分布白噪声

for i=1:n-1
 yy(i,1)=yy(i+1,1);
end
yy(n,1)=y;
yyy=[yyy;y];
for i=1:n
    X(1,i)=uu(i,1);
end
X(1,n+1)=1;
K=P*X'*inv(t1+X*P*X');
P=(eye(n+1)-K*X)*P/t1;
Q=Q+K*(y-X*Q);
%根据元素值g0,g1,...gn,求G矩阵
for j=1:m
    for i=n:-1:j
        i1=n-i+j;
        G(i1,j)=Q(i,1);
    end
end
%求nn维y0向量(y1为上一时刻y0的向量) 
y0(1:nn-1,1)=y1(2:nn,1); y0(nn,1)=y1(nn,1);
y0=y0+(y-y1(1,1));
for i=1:n
    y1(i,1)=y0(i,1)+u(1,1)*Q(n+1-i,1);
end 
for i=n+1:nn
    y1(i,1)=y1(n,1);
end
%根据y0,求n维f向量f(k+1),...,f(k+n)
f(1:n,1)=y0(1:n,1);
%由当前k时刻的输出值y(k)和给定值yr，求k时刻后的n个参考轨迹w(k+1),...,w(k+n)
w=[a*y+(1-a)*yr];
for i=2:n
    w=[w;a^i*y+(1-a^i)*yr];
end
%计算k时刻以后的m个控制增量Du(k),...,Du(k+m)
u=inv(G'*G+t0*eye(m))*G'*(w-f);
%保存k时刻及以后n个控制增量，以供模型运算
for i=1:n-1
    uu(i,1)=uu(i+1,1);
end
uu(n,1)=u(1,1);
uuu=[uuu;u(1,1)];%保存各k时刻的nm个控制增量以便绘图
%控制量限幅
if (u(1,1)>1)
    u(1,1)=1;
end
if(u(1,1)<-1)
    u(1,1)=-1;
end
end
%绘制完整的给定值，输出值和控制增量曲线
subplot(2,1,1):plot(t,(yr0+1),t,yyy);
axis([0,nm,0,4]);  xlabel('t');ylabel('yr,y')
subplot(2,1,2):plot(t,uuu);
axis([0,nm,-1.5,1.5]);xlabel('t');ylabel('u')