广义预测控制 GPC

广义预测控制（GPC）

GPC算法仿真

被控对象模型

动态矩阵控制算法的编程原理

(1)设置GPC参数,例如采样周期，预测时域，控制时域，截断步长等。

(2)建立系统阶跃响应模型

(3)设置初始时刻参数，例如系统的初始时刻值，柔化系数等。

(4)计算参考轨迹

(5)计算控制作用增量

(6)实施GPC控制

(7)输出结果，绘制曲线

GPC算法：

1.初选控制参数：Q、R、P、M、 ysp 、α、Â(z-1)

2.采集输入、输出样本{∆u(k),∆y(k)}

3.用RLS算法估计参数

4.递推求解Diophantine方程，得到

5.计算F(k)

6.在线计算控制器参数d T

7.得到控制增量∆u(k)和控制输入u(k) ＝u(k－1) ＋∆u(k)

8.k+1 →k，进入下一周期预测计算和滚动优化

GPC程序：

%Clarke广义预测控制（C=1）（对象参数已知）

%N1=d、N、Nu取不同的值

clear all；close all;

a=cell(1,2) ;b=cell(1,2) ;c=cell(1,1);d=cell(1,1);%对象参数syms k;

k=length(k);

if (0<=k<=150)

a=[1 0.9234]; b=[7.2402 0.9485]; c=1; d=1;

elseif (150

a=[1 0.8981]; b=[9.9901 0.14142]; c=1; d=1; elseif (300

a=[1 0.8838]; b=[9.6041 0.34067]; c=1; d=1;

else (450

a=[1 0.9234]; b=[7.2402 0.9485]; c=1; d=1;

end

na=length(a)-1;b=[zeros(1,d-1) b];nb =length(b)-1;%na、nb为多项式A、B阶次（因d！=1,对b添0）

aa=conv(a,[1 -1]);naa=na+1;%aa的阶次

N1=d;N=15;Nu=5;%最小输出长度、预测长度、控制长度

gamma=1*eye(Nu);alpha=0.11;%控制加权矩阵、输出柔化系数

L=600;%控制步数

uk=zeros(d+nb,1);%输入初值：uk(i)表示u(k-i)

duk=zeros(d+nb,1);%控制增量初值

yk=zeros(naa,1);%输出初值

w=10*[ones(L/4,1);-ones(L/4,1);ones(L/4,1);-ones(L/4+d,1)]; %设定值

xi=sqrt(0.01)*randn(L,1);%白噪声序列

%求解多步Diophantine方程并构建F1、F2、G

[E,F,G]=multidiophantine(aa,b,c,N);

G=G(N1: N, : );

F1=zeros(N-N1+1,Nu); F2=zeros(N-N1+1,nb);

for i=1:N-N1+1

for j=1:min(i,Nu); F1(i,j)=F(i+N1-1,i+N1-1-j+1);end for j=1:nb; F2(i,j)=F(i+N1-1,i+N1-1+j);end

end

for k=1:L

if (1<=k<=150)

time(k)=k;

a=[1 0.9234]; b=[7.2402 0.9485]; c=1; d=1;

y(k)=-aa(2:naa+1)*yk+b*duk(1:nb+1)+xi(k);%采集输出数据Yk=[y(k);yk(1:na)];%构建向量Y(k)

dUk=duk(1:nb);%构建向量△U(k-j)

elseif (150

time(k)=k;

a=[1 0.8981]; b=[9.9901 0.14142]; c=1; d=1;

y(k)=-aa(2:naa+1)*yk+b*duk(1:nb+1)+xi(k);%采集输出数据Yk=[y(k);yk(1:na)];%构建向量Y(k)

dUk=duk(1:nb);%构建向量△U(k-j)

elseif (300

time(k)=k;

a=[1 0.8838]; b=[9.6041 0.34067]; c=1; d=1; y(k)=-aa(2:naa+1)*yk+b*duk(1:nb+1)+xi(k);%采集输出数据Yk=[y(k);yk(1:na)];%构建向量Y(k)

dUk=duk(1:nb);%构建向量△U(k-j)

else (450

time(k)=k;

a=[1 0.9234]; b=[7.2402 0.9485]; c=1; d=1; y(k)=-aa(2:naa+1)*yk+b*duk(1:nb+1)+xi(k);%采集输出数据Yk=[y(k);yk(1:na)];%构建向量Y(k)

dUk=duk(1:nb);%构建向量△U(k-j)

end

%参考轨迹

yr(k)=y(k);

for i=1:N

yr(k+i)=alpha*yr(k+i-1)+(1-alpha)*w(k+d); end

Yr=[yr(k+N1:k+N)]';%构建向量Yk(k)

%求控制量

dU=inv(F1'*F1+gamma)*F1'*(Yr-F2*dUk-G*Yk); %ΔU du(k)=dU(1); u(k)=uk(1)+du(k);

%更新数据

for i=1+nb:-1:2

uk(i)=uk(i-1);

duk(i)=duk(i-1);

end

uk(1)=u(k);

duk(1)=du(k);