广义预测控制 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<k<=300)
a=[1 0.8981]; b=[9.9901 0.14142]; c=1; d=1; elseif (300<k<=450)
a=[1 0.8838]; b=[9.6041 0.34067]; c=1; d=1;
else (450<k<=600)
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<k<=300)
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<k<=450)
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<k<=L)
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);
for i=naa:-1:2
yk(i)=yk(i-1);
end
yk(1)=y(k);
end
subplot(2,1,1);
plot(time,w(1:L),'m:',time,y);
xlabel('k');ylabel('w(k)、y(k)');
legend('w(k)','y(k)');
subplot(2,1,2);
plot(time,u);
xlabel('k');ylabel('u(k)');
function[E,F,G]=multidiophantine(a,b,c,N)
%********************************************************** *
%功能：多步Diophanine方程的求解
%调用格式：[E,F,G]=sindiophantine(a,b,c,N)(注：d=1)
%输入参数：多项式A,B,C系数向量及预测步数（共4个）
%输出参数：Diophanine方程的解E,F,G（共3个）
%********************************************************** ***
na=length(a)-1;nb =length(b)-1;nc=length(c)-1;%A、B、C的阶次
%E、F、G的初值
E=zeros(N);E(1,1)=1;F(1,:)=conv(b,E(1,:));
if na>=nc
G(1,:)=[c(2:nc+1) zeros(1,na-nc)]-a(2:na+1);%令c(nc+2)=c(nc+3)=...=0
else
G(1,:)=c(2:nc+1) -[a(2:na+1)-zeros(1,nc-na)];%令a(nc+2)=a(nc+3)=...=0
end
%求E、F、G
for j=1:N-1
for i=1:j
E(j+1,i)=E(j,i);
end
E(j+1,j+1)=G(j,1);
for i=2:na
G(j+1,i-1)=G(j,i)-G(j,1)*a(i);
end
G(j+1,na)=-G(j,1)*a(na+1);
F(j+1,: )=conv(b,E(j+1,:));
end
仿真结果
N=15 Nu=5 alpha=0.11
N=10 Nu=5 alpha=0.11
N=15 Nu=3 alpha=0.11
N=15 Nu=3 alpha=0.31
结论
可以得出,当保持其他参数不变而改变一或几个变量时会有不同的情形。

当预测步长越大，系统越稳定，预测步长越小，系统的快速性会变好，但对系统响应没有多大影响。

当预测步长不变时，随着控制步长的减小，系统的稳定性和鲁棒性会更好，但可能跟踪性能会变差。

随着控制步长的增大，系统的灵敏性更好，动态性能会更好。

柔化系数越大参考轨迹柔性变好，鲁棒性变好，而系统的快速性下降。

心得体会
通过这次对GPC的仿真，让我更好的掌握了Matlab的运用方法，更熟悉了Matlab的一些技巧和仿真步骤。

虽然在先前的几次仿真中遇到或失败或者错误之类的，但经过同学之间的交流和师兄的指导和查阅相关资料，这些问题都得到了一一克服，我觉得这是对我对书本知识运用的一次很好的历练。