daichunzhihoudejietigpc2.m

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

%有纯滞后的阶梯控制，纯滞后时间为k-1
tic;
%gpc1(Aplant,Bplant,q,p,aifa).m
clear all;close all; %程序开始，清空工作间关闭窗口
%创建传递函数形式的模型，采样周期指定为Ts,可由具体要求而定

%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
Amodel = [1 -1.2 0.36]; 
Bmodel = [0.1 0.2];
na = length(Amodel)-1;
nb = length(Bmodel)-1; %控制对象的传递函数分子A,分母B的阶次分别为na ,nb

k = 4; %从第k 步开始预测(输出滞后k-1,k=N1）
st = 1; %系统设定值setpoint在初始时读入
Ts = 1; %Ts为采样周期(步长）
T_final = 50; %T_final 仿真时间 


p = 3; %根据需要选取预测步长P
M = 3; %控制步长为M
aifa = 0 ; %输入柔化因子aifa
lambda = 0; %输入增量参数lambda
b=0.02;

Aplant = Amodel;
Bplant = Bmodel;
h=[];
% 下面构造F1,E1
F(1,1:na) = Amodel(1,1:na)-Amodel(1,2:na+1);
F(1,na+1) = Amodel(1,na+1);%F（1，：）(q^-1)=(1-a1)+(a1-a2)*q^-1+...+ana*q^-na
E(1,1)=1; 
G_ef(1,:) = Bmodel;

for j= 2:k+p-1
E(1,j) = F(j-1,1);
F(j,1:na) = F(j-1,2:na+1) - E(1,j)*(Amodel(1,2:na+1)-Amodel(1,1:na));%其中E(1,j)=F(j-1,1)
F(j,na+1) = F(j-1,1)*Amodel(1,na+1);
G_ef(j,1:nb+j) = conv(E,Bmodel);
end
%以上为构造Fj(q^-1)、Ej(q^-1);Ej(q^-1)为j-1次首一多项式,Fj(q^-1)为na次
%构造G_ef,为nb+j-1次(含有nb+j项),得到g0,....,gnb+j-1
for i = 1:p%构造讲义中的 G
for j = 1:i
G(i,j) = G_ef(i,i-j+1);
end
end 

for i= 1:p  %构造计算y^1（即G_poly*delte44U+F*yk )时所需的矩阵G_poly
G_poly(i,1:nb+k-1) = G_ef(k-1+i,i+1:nb+k+i-1);
end

G1 = G(1:p,1:M);%根据预测控制时域得到的G1

ypast(:,1) = zeros(na+1,1); %初始值ypast(na+1项)，从y(t-1)开始至y(t-na-1).
Upast(:,1) = zeros(nb+k+1,1); %预测前已知的控制量u(nb+k)，从u(t-1)开始至u(t-nb-k),deltu(t-nb-k+1)时用到u(t-nb-k)
for i=1:M
    %h=[h;b^(i-1)];
    f=0;
 for j=1:i
     f=f+b^(j-1);
 end
    h=[h;f];
end
G2=G1*h;%引入阶梯因子后的G

outU = [];
outY = [];

for I = 0:Ts:T_final %开始循环计算
    if I>25
        st=0.9;
    end
delt_ypast = ypast(1:na,1)-ypast(2:na+1,1);
deltU = Upast(1:nb+k,1) - Upast(2:nb+k+1,1);%产生已知输入控制信号deltu(t)
yt = ypast(1,1) - Aplant(1,2:na+1) * delt_ypast(1:na,1) + Bplant(1,:) * deltU(k:nb+k,1);
%yt = ypast(1,1) - Aplant(1,2:na+1) * delt_ypast(1:na,1) + Bplant(1,:) * deltU(1:nb+1,1);
outY = [outY;yt];
%采样，得到仿真对象的当前输出值yt
ypast = [yt;ypast];
ypast = ypast(1:na+1,1);
for i = 1:p
ypre1(i,1) = F(k+i-1,:) * ypast(:,1) + G_poly(i,:) * deltU(1:nb+k-1,1);
end
%以上为求y^1,ypre1为
%aifa=0.95^(I+1);%输入的柔化因子aifa
%if k == 1
W(1,1) = yt;
%else
%W(1,1) = F(k-1,:) * ypast(:,1) + G_ef(k-1,1:nb+k-1) * deltU(1:nb+k-1,1);
%end
for i = 2:p+1
W(i,1)=aifa*W(i-1,1)+(1-aifa)*st;
end %未来t+j时刻系统的柔化设定值记为w(j,1)
deltU1 = inv(G2'*G2)*G2'*(W(2:p+1,1)-ypre1(:,1)) ;%当前控制律deltU1
ut = Upast(1,1)+deltU1(1,1);
Upast = [ut;Upast];
Upast = Upast(1:nb+k+1,1);
outU = [outU;ut];
end
t=[0:Ts:T_final];
plot(t,outY),hold on,
figure(2);
plot(t,outU,'r')
toc;