jietishigpc.m

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

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

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

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

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
for i = 1:p
for j = 1:i  %构造多步预测中的 G
G(i,j) = G_BR(i,i-j+1);
end
end 

for i= 1:p
F0(i,1:nb) = G_BR(i,i+1:nb+i);%构造计算所需的矩阵F0,因为deltU元素排列与书上不一样
end
G1 = G(1:p,1:M);%根据预测控制时域得到的G1，大于M的输入为0

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项
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);%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)
for i = 1:p
ypre1(i,1) = S(i,:) * ypast(:,1) + F0(i,:) * deltU(1:nb,1);%S中只有sj,i,没有z^(-i)，得列向量
end  %以上为求y^1,ypre1为由已知量求得
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)
%deltU1 = inv(G1'*G1+landa*eye(M))*G1'*(W(2:p+1,1)-ypre1(:,1));%当前控制律向量deltU1
%ut = Upast(1,1)+deltU1(1,1);
%ut=Upast(1,1)+inv(G2'*G2+landa*h'*h)*G2'*(W(2:p+1,1)-ypre1(:,1));
ut=Upast(1,1)+inv(G2'*G2)*G2'*(W(2:p+1,1)-ypre1(:,1));
Upast = [ut;Upast];
Upast = Upast(1:nb+2,1);
outU = [outU;ut];
end
t=[0:Ts:T_final];

figure(2);
hold on;
axis([0,50,0.4,1.2]);
plot(t,outU,'--');gtext('2');
%legend('2-未补偿阶梯式控制',1);
xlabel('t/s');
ylabel('u');

figure(1);hold on
axis([0,50,0.85,1.15]);
plot(t,outY,'--');gtext('2');
%legend('2-未补偿阶梯式控制',1);
xlabel('t/s');
ylabel('y');

toc;