GPC-matlab程序

%%

%GPC控制器函数

function [d_unow] = Control(A,B,SP,Zk,d_u,lanbt,Error)

%A表示差分方程右边输出变量的系数（从y(t-1),y(t-2),......,y(t-model_order1)）

%B表示差分方程右边输入变量的系数（从u(t-1),u(t-2),......,u(t-model_order2)）

%%

a = A(2:end); %差分方程右侧输出变量系数

b = B; %差分方程右侧输入变量系数

N0 = 1; %最小预测时域长度，一般取N0=1

N1 = 25; %最大预测时域长度

Nu = 25; %控制时域长度

% b = [zeros(1,9),B(10:end)];

E = zeros(N1,N1); %存储误差

% lanbt = 1; %加权系数

alpha = 0.19;

na = length(a);

nb = length(b)-1;

aa = zeros(1,na+1); if na>=2

for i = 2:na

aa(i) = (a(i)-a(i-1)); end

end

aa(1) = (a(1)-1);

%定义一个与输出变量阶数相等的数组

%计算差分方程右侧输出变量相邻阶次的系数差值

aa(na+1) = -a(na);

%%

%定义Diophantine方程求解过程多项式，参考王伟编写的《广义预测控制》2.1节广义预测控制的基本方法和2.2节Diophantine方程的递推求解

F = zeros(N1,na+1); %维数=预测时域长度*输出变量阶数

F(1,:) = -aa; %F初值

e = zeros(1,N1); %定义误差缓存数组

% f = zeros(1,na+1);

H = zeros(N1,nb); %维数=预测时域长度*（输入变量阶数-1）

H(1,:) = b(2:end); %H初值

g = zeros(1,N1);

e(1) = 1; %初始E=1

g(1) = e(1)*b(1); %初始G=e(1)*b(1)

%%

%求解diophantine方程

%递推求解F，E，H，G，参见王伟编写书的2.2节内容

for j = 1:N1-1

for i = 1:na

F(j+1,i) = F(j,i+1)-aa(i)*F(j,1);

end

F(j+1,na+1) = -aa(na+1)*F(j,1);

e(j+1) = F(j,1);

if nb>1

for i = 1:nb-1

H(j+1,i) = e(j+1)*b(i+1)+H(j,i+1);

end

end

H(j+1,nb) = e(j+1)*b(nb+1);

g(j+1) = e(j+1)*b(1)+H(j,1);

end

%计算矩阵G

G = zeros(N1,Nu); %定义G

for i = 1:Nu

G(i:end,i) = g(1:N1-i+1);

G(1:N0,:) = 0; %预测时域从N0~N1，固小于N0的G矩阵部分都取0，即无法预测

F(1:N0,:) = 0; %道理同G

H(1:N0,:) = 0; %同上

end

%计算（G'*G+lanbt*eye(size(G'*G))）的逆矩阵

aaa = size(G'*G);

R = diag([1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1]); %N1*N1(24*24)单位对角阵

% R = diag([1,1,1,1,1,1,1,1,1,1,1,1,1,1,1]); %N1*N1单位对角阵

Q = inv(G'*R*G+lanbt*eye(aaa(1)))*G'*R;

p = Q(1,:);

yr = zeros(N1+1,1);

yr(1) = Zk(1);

for i = 1:N1

yr(i+1) = alpha*yr(i)+(1-alpha)*SP;

end

% f = p(1)*F(1,:);

% for i = 2:N1

% f = f+ p(i)*F(i,:);

% end

yr1 = yr(2:end);

yr1=[SP;yr1(1:24)];

% yr1 = [yr11;yr12;yr13;yr14;yr15;yr16;yr17;yr18;yr19;yr20;yr1(1:15)];

% Error = 1.5782;

a1 = p*yr1;

a2 = p*F*Zk;

a3 = p*H*d_u;

d_unow = a1-a2-a3-Error;