gpc.m

隐匿广义预测自校正控制算法,基于MATLAB开发环境

clear
disp('广义预测控制算法')
nn = input('时域长度nn=');
n = input('预测长度n=');
m = input('控制长度m=');
t0 = input('控制加权系数=');
a = input('柔化系数a=');
disp('最小二乘公式初始值')
t1 = 1;
d1 = input('(n+1)阶方阵P的形式：0-对角阵，1-方阵：');
d2 = input('(n+1)阶方阵P的初始值：0-自动赋值，1-键盘输入：');
if(d1==1)
    if(d2==1)
        P = input('在方括号[]中，输入(n+1)阶方阵P的值:');
    else
        P = (1e+5)*ones(n+1);
    end
else
    if(d2==1)
        PP = input('在方括号[]中，输入(n+1)阶对角方阵P的值:');
        P = diag(PP);
    else
        P = (1e+5)*eye(n+1);
    end
end
%初始化参数
uuu = 0; yyy = 0;
uu = zeros(n,1); u = zeros(m,1);
yy = zeros(n,1); y1 = zeros(n,1);
Q = zeros(n+1,1); Q(1,1) = 1; Q(n+1,1) = 1;
%产生周期为100,时间为T,幅值为1的方波信号的给定值;
T = 300; [yr0,t] = gensig('square',100,T,1);
d3 = input('输出曲线是否去掉前100步:0-不,1-去掉');
nm = length(t);%确定循环次数
for ij = 2:nm
    yr = yr0(ij) + 1;%产生周期为100,时间为T,幅值在1和2之间变化的方波信号的给定值
    %根据系统模型,计算k时刻的输出值y(k)
    y = 2.5*yy(n,1)-2.2*yy(n-1,1)+0.7*yy(n-2,1)+(uu(n,1)+1.5*uu(n-1,1))/12.5;
    %(nn=6;n=6;m=2;t0=0.8;a=0.5;t1=1)
    %产生均匀分布的白噪声
    a9 = 0;
    for i=1:1
        a9 = a9 + rand;
    end
    a8 = 0.01*(a9-6);
    %保存k时刻及以前的n个输出值y(k),y(k-1),...,y(k-n),以供模型运算
    for i=1:1
        yy(i,1) = yy(i+1,1);
    end
    yy(n,1) = y;
    yyy = [yyy;y];%保存各k时刻的nm个输出量以便绘图
    %根据最小二乘公式,由y(k)计算G阵的各元素值g0,g1,...,gn
    for i=1:n
        X(1,i) = uu(i,1);
    end
    X(1,n+1) = 1;
    K = P*X'*inv(t1+X*P*X');
    P = (eye(n+1)-K*X)*P/t1;
    Q = Q+K*(y-X*Q);
    %根据元素值g0,g1,...,gn,求G阵
    for j=1:m
        for i=n:-1:j
            i1 = n-i+j;
            G(i1,j) = Q(i,1);
        end
    end
    %根据y1(y1为上一时刻的向量),求n维f向量f(k+1),...,f(k+n)
    e = y-y1(1,1);
    f(1:n-1,1) = y1(2:n,1)+e;
    f(n,1) = y1(n,1)+e;
    y1 = f+G*u;
    %由当前k时刻的输出值y(k)和给定值yr,求k时刻以后的n个参考轨迹w(k+1),...,w(k+n)
    w = a*y+(1-a)*yr;
    for i=2:n
        w = [w;a^i*y+(1-a^i)*yr];
    end
    %计算k时刻及以后的m个控制增量Du(k),...Du(k+m)
    u = inv(G'*G+t0*eye(m))*G'*(w-f);
    %保存k时刻及以后的n个控制增量,以模型运算
    for i=1:n-1
        uu(i,1) = uu(i+1,1);
    end
    uu(n,1) = u(1,1);
    uuu = [uuu,u(1,1)];%保存各k时刻的nm个控制增量以便绘图
    %控制量限幅
    if(u(1,1)>0.5*yr)
        u(1,1) = 0.5*yr;
    end
    if(u(1,1)<-0.5*yr)
        u(1,1) = -0.5*yr;
    end
end
%绘制给定值、输出值和控制增量曲线
if(d3==1)
    %绘制去掉前100的给定值、输出值和控制增量曲线
    yyy1(1:(T-100),1) = yyy(101:T,1);
    uuu1(1:(T-100),1) = uuu(101:T,1);
    t1(1:(T-100),1) = t(101:T,1);
    yr01(1:(T-100),1) = yr0(101:T,1);
    subplot(2,1,1): plot(t1,(yr01+1),t1,yyy1);
    axis([0,nm-100,0,2.5]);
    xlabel('t');
    ylabel('yr,y');
    subplot(2,1,2): plot(t1,uuu1);
    axis([0,nm-100,-1.5,1.5]);
    xlabel('t');
    ylabel('u');
else
    %绘制完整的给定值、输出值和控制增量曲线
    subplot(2,1,1): plot(t,(yr0+1),t,yyy);
    axis([0,nm,0,2.5]);
    xlabel('t');
    ylabel('yr,y');
    subplot(2,1,2): plot(t,uuu);
    axis([0,nm,-1.5,1.5]);
    xlabel('t');
    ylabel('u');
end