danzhidegpc.m

阶梯式广义预测控制的程序,单值广义预测控制,有纯滞后的阶梯式广义预测控制,一般有纯滞后的GPC及带迟滞的一阶环节辨识

%两次都为单值预测的充分利用信息
clear;
Amodel = [1 -1.474 0.474]; %A(q^-1)=1+a1*q^-1+a2*q^-2+......+ana*q^ -na
Bmodel = [0.445 0.163 0.029]; %B(q^-1)=bo+b1*q^-1+b2*q^-2+......+bnb*q^ -nb,nb=1
na = length(Amodel)-1;
nb = length(Bmodel)-1; %控制对象的传递函数分子A,分母B的阶次分别为na ,nb

p = 3; %根据需要选取预测步长P
M = 3; %控制步长为M 
%landa =0.8; %控制量加权系数landa
aifa =0; %输入柔化因子aifa

st = 1; %系统设定值setpoint在初始时读入
Ts = 1; %Ts为采样周期
T_final =50; %T_final 仿真时间
Aplant = Amodel;%实际对象的传函系数，分母，分子
Bplant = Bmodel;

S(1,1:na) = Amodel(1,1:na)-Amodel(1,2:na+1);
S(1,na+1) = Amodel(1,na+1);%S1(q^-1)=(1-a1)+(a1-a2)*q^-1+...+ana*q^-na
R(1,1)=1; %构造R1,S1
G_BR(1,:) = Bmodel;

for j= 2:p
R(1,j) = S(j-1,1);
S(j,1:na) = S(j-1,2:na+1) - R(1,j)*(Amodel(2:na+1)-Amodel(1:na));
S(j,na+1) =S(j-1,1)*Amodel(1,na+1);
G_BR(j,1:nb+j) = conv(R,Bmodel);
end
%以上为构造Sj(q^-1)、Rj(q^-1);Rj(q^-1)为j-1次首一多项式,Sj(q^-1)为na次
%构造G_BR,为nb+j-1次(含有nb+j项),得到g0,....,gnb+j-1
ypast(:,1) = zeros(na+1,1); %初始值ypast(na+1项)，从y(t-1)开始至y(t-na-1),列向量
Upast(:,1) = zeros(nb+2,1); %预测前已知的控制量u(nb+2项)，从u(t-1)开始至u(t-nb-1),列向量u(t-1)*z^(-nb),
                            %delt u(t-nb-1)就需要u(t-nb-2)，共nb+2项
outU = [];
outY = [];                            
for I = 0:Ts:T_final %开始循环
    delt_ypast = ypast(1:na,1)-ypast(2:na+1,1);%delty(t-1)---delty(t-na)
    deltU = Upast(1:nb+1,1) - Upast(2:nb+2,1);%产生已知输入控制信号deltu(t-1)---delt(t-nb-1)
%采样，得到仿真对象的当前输出值yt
yt = ypast(1,1) - Aplant(1,2:na+1) * delt_ypast(1:na,1) + Bplant(1,:) * deltU(1:nb+1,1);
outY = [outY;yt];%列向量下增
ypast = [yt;ypast];
ypast = ypast(1:na+1,1);%ypast为y(t)---y(t-na)
W(1,1)=yt;
for i = 2:p+1
W(i,1)=aifa*W(i-1,1)+(1-aifa)*st;
end %未来t+j时刻系统的柔化设定值记为w(j,1)
deltut=((W(p+1,1)-S(p,:)*ypast(:,1))-G_BR(1,p+1:end)*deltU(1:nb,1))/G_BR(p,p);
ut = Upast(1,1)+deltut;%无补偿的计算输入
Upast = [ut;Upast];
Upast = Upast(1:nb+2,1);
deltU = Upast(1:nb+1,1) - Upast(2:nb+2,1);
y1=G_BR(1,1:nb+1)*deltU+S(1,:)*ypast;%做为t+1时刻的已知"yt"的预测值
ypasth=ypast;
ypasth = [y1;ypast];
ypasth = ypasth(1:na+1,1);
deltut1=((st-S(p,:)*ypasth(:,1))-G_BR(1,p+1:end)*deltU(1:nb,1))/G_BR(p,p);
uth=ut+0.4*deltut1;%充分利用信息后的控制量 
Upast(1,1)=uth;
outU = [outU;uth];
%Upast = [ut;Upast];
%Upast = Upast(1:nb+2,1);
%outU = [outU;ut];
end
t=[0:Ts:T_final];
figure(1);hold on
%axis([0,50,0,1.1]);grid on
plot(t,outY)
xlabel('t/s');
ylabel('y');
%plot(t,outU);%gtext('基本');
%figure(2);
%axis([0,50,-0.1,0.6]);
hold on;
%plot(t,outU)
xlabel('t/s');
ylabel('y');