gpc2.m

预测控制是适用于含有大滞后环节、积分环节的复杂工业系统控制。GPC为广义预测控制

%gpc1(Aplant,Bplant,q,p,aifa).m
clear all;close  all;   %程序开始，清空工作间关闭窗口
%创建传递函数形式的模型，采样周期指定为Ts,可由具体要求而定
 
Amodel = [1 -1.4 0.49 -0.343];             %A(q^-1)=1+a1*q^-1+a2*q^-2+......+ana*q^ -na 
Bmodel = [-0.2 1];             %B(q^-1)=bo+b1*q^-1+b2*q^-2+......+bnb*q^ -nb
%Amodel = [1 -0.9355];
%Bmodel = [0.06449];
na = length(Amodel)-1;        
nb = length(Bmodel)-1;              %控制对象的传递函数分子A,分母B的阶次分别为na ,nb
k = 4;                              %从第k 步开始预测   
p = 14;                             %根据需要选取预测步长P
SPt = 10;                           %系统设定值setpoint在初始时读入
aifa = 0.99 ;                       %输入柔化因子aifa 
Ts = 1;                             %Ts为采样周期
T_final = 500;                      %T_final  仿真时间
lambda = 0;                         %读入参数lambda
M = 6;                              %控制步长为M
Aplant = Amodel;
Bplant = Bmodel;
            
F(1,1:na) = Amodel(1,1:na)-Amodel(1,2:na+1);              
F(1,na+1) = Amodel(1,na+1);         %F1(q^-1)=(1-a1)+(a1-a2)*q^-1+...+ana*q^-na
E(1,1)=1;                           % F1=poly2sym(F1)   构造F1,E1
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));
    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
    for j = 1:i 
        G(i,j) = G_ef(i,i-j+1);
    end
end                              %构造讲义中的 G                             

for i= 1:p
        G_poly(i,1:nb+k-1) = G_ef(k+i-1,i+1:nb+k+i-1);
end                                 
%构造计算y^1时所需的矩阵G_poly
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)
outU = [];
outY = [];
for  I =  0:Ts:T_final 
    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);
    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
    
    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)*SPt;
    end                                    %未来t+j时刻系统的柔化设定值记为w(j,1)
    deltU1 = inv(G1'*G1+lambda*eye(M))*G1'*(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,plot(t,outU,'r')