test.m

Clark等人提出的广义预测控制自校正控制器是一种基于参数模型的预测控制算法

clear all;close all;%程序开始，清空工作间关闭窗口

s=tf('s');
Ts=12;%Ts为采样周期(步长）
T_final = 125*Ts; %T_final 仿真时间 
g11=0.162*(275*s+1)/((168*s+1)^2*(11.5*s+1));
figure(1);
step(g11);
%g12=-0.081*(0.01+0.99/(97*s+1));
%g21=2.116*(457*s+1)/((221*s+1)^2*(21.8*s+1));
%g22=1.483*s*(150*s+1)/((632*s^2+40*s+1)*(2.7*s+1));
%G=[g11,g12;g21,g22];
%step(G);
G11=c2d(g11,Ts);
[num,den]=tfdata(G11,'v');

%initialization
na = length(den)-1;
nb = length(num)-1;
k = 1;%从第k 步开始预测
st = 1;%系统设定值setpoint在初始时读入
p = 20; %根据需要选取预测步长P
M = 8; %控制步长为M
aifa = 0.05 ; %输入柔化因子aifa
lambda = 10; %输入增量参数lambda
b = 1;%wug's的阶梯因子
h=[];

%the first line of F and G
F(1,1:na) = den(1,1:na)-den(1,2:na+1);
F(1,na+1) = den(1,na+1);%F（1，：）(q^-1)=(1-a1)+(a1-a2)*q^-1+...+ana*q^-na
E(1,1)=1;
G_ef(1,:) = num;

%solve the diophantus equation
for j= 2:k+p-1
    %get G
    G_ef(j,1:j-1) = G_ef(j-1,1:j-1);
    for i= j:nb+j-1
        G_ef(j,i) = G_ef(j-1,i) + F(j-1,1)*num(i-j+1);
    end
    G_ef(j,nb+j) = F(j-1,1) * num(nb+1);
    %get F
    for i= 1:na
        F(j,i) = F(j-1,i+1) + F(j-1,1)*F(1,i);
    end
    F(j,na+1) = F(j-1,1) * den(na+1);
end
%*-------------finished solving the diophantus equation-----------*

%calculate wug's G2
for i = 1:p%构造讲义中的 G
for j = 1:i
G(i,j) = G_ef(i,i-j+1);
end
end

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

for i=1:M
    %h=[h;b^(i-1)];
    h=[h;b^(i-1)];
end
G2=G1*h;%引入阶梯因子后的G
% end calculate G2

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

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)
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) - den(1,2:na+1) * delt_ypast(1:na,1) + num(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,ypre1为
    
    %Tender the setpoint track vector
    W(1,1) = yt;
    for i = 2:p+1
        W(i,1)=aifa*W(i-1,1)+(1-aifa)*st;
    end
    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];
figure(2);
plot(t,outY),hold on,
figure(3);
plot(t,outU,'r')