自校正GPC matlab仿真

自校正广义预测控制 matlab仿真 自校正GPC

%广义预测控制_自校正GPC算法仿真
%被控模型为y（k）=1.5*y（k-1）-0.7*y（k-2）+u（k-1）+0.5*u（k-2）+e（k），其中e（k）为噪声
clc;
clear;
%仿真步数
tim=50;
%预辨识次数
tim0=10;
%预测时域和控制时域
N=5;
Nu=4;
%扰动幅度
ampe=1;
%系数λ
lmd=0.5;
%y（k）的最终给定值yr
yr=50;
%柔化因子a
a=0.5;
%β
b=1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%参数预辨识
%将系统在正弦信号控制作用下运行tim0步，测的一组系统的运行数据
e=ampe*（1-2*rand（1tim0））;
yy（1）=0;yy（2）=0;
uu（1）=0;uu（2）=0;
for k=3:tim0
    uu（k）=6+sin（k）;
end
for k=3:tim0
    yy（k）=1.5*yy（k-1）-0.7*yy（k-2）+uu（k-1）+0.5*uu（k-2）+e（k）;
end
%根据所得运行数据，使用递推最小二乘法辨识系统模型参数
A=[1-11];
B=[11];
na=length（A）-1;
nb=length（B）-1;
a1（1）=-1;a2（1）=1;b0（1）=1;b1（1）=1;
THETA=[-A（na:-1:1）B（nb+1:-1:1）]‘;
P=b^2*eye（na+nb+1）;
for k=3:tim0
    PHI=[yy（k-1:-1:k-na）uu（k-1:-1:k-nb-1）]‘;
    THETA=THETA+P*PHI*（yy（k）-THETA‘*PHI）/（1+PHI‘*P*PHI）;
    P=P-P*PHI*PHI‘*P/（1+PHI‘*P*PHI）;
end
A（na:-1:1）=-THETA（1:na）;
B（nb+1:-1:1）=THETA（na+1:na+nb+1）;
a1（2）=A（2）;a2（2）=A（1）;b0（2）=B（2）;b1（2）=B（1）;

%至此预辨识结束
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%一些计算的中间量
Z=[10];%表示z^（-1）
Zj=cell（1N-1）;%Zj中存储的为z^（-j）次
for j=1:N-1
    Zj{j}（1）=1;
    Zj{j}（2:j+1）=0;
end
delta=[-11];%△

%系统状态的初值
y（1）=0;y（2）=0;
u（1）=0;u（2）=0;
e=ampe*（1-2*rand（1tim））;
%下面开始迭代进行自校正GPC控制
for k=3:tim
    y（k）=1.5*y（k-1）-0.7*y（k-2）+u（k-1）+0.5*u（k-2）+e（k）;
    %根据新采集到的系统运行数据进行一次递推最小二乘法辨识
    PHI=[y（k-1:-1:k-na）u（k-1:-1:k-nb-1）]‘;
    THETA=THETA+P*PHI*（y（k）-THETA‘*PHI）/（1+PHI‘*P*PHI）;
    P=P-P*PHI*PHI‘*P/（1+PHI‘*P*PHI）;
    %一步辨识结束，得到参数
    A（na:-1:1）=-THETA（1:na）;
    B（nb+1:-1:1）=THETA（na+1:na+nb+1）;
    a1（k）=A（2）;a2（k）=A（1）;b0（k）=B（2）;b1（k）=B（1）;
    Adelta=conv（Adelta）;%A△
    %迭代求第一个丢番图方程中的EF
    E=cell（1N）;
    F=cell（1N）;
    E{1}=[1];
    F{1}=deconv（polyadd（1-Adelta）Z）;
    for j=1:N-1
        E{j+1}=polyadd（E{j}F{j}（na+1）*Zj{j}）;
        F{j+1}=deconv（polyadd（F{j}-（F{j}（na+1）*Adelta））Z）;
    end
    %求第二个丢番图方程中的GH
    G=cell（1N）;
    H=cell（1N）;
    for j=1:N
        EB=conv（E{j}B）;
        H{j}=EB（1:length（EB）-j）;
        G{j}=EB（length（EB）-j+1:length（EB））;
    end
    %将EFGH左右翻转，因为MATLAB中的多项式系数是从高次到低次排列，而GPC算法中是相反的，这样翻转过来使用比较方便
    for j=1:N
        E{j}=fliplr（E{j}）;
        F{j}=fliplr（F{j}）;
        G{j}=fliplr（G{j}）;
        H{j}=fliplr（H{j}）;
    end
    %得出后边计算要用到的系数矩阵F_H_G_
    for j=1:N
        F_（j1:na+1）=F{j};
        H_（j1:nb）=H{j};
    end
    for i=1:Nu
        G_（1:Ni）=[zeros（1i-1）G{N}（1:N+1-i）]‘;
    end
    %△U（k）=K（k）（Yd（k）-Y0（k））
    K=（G_‘*G_+lmd*eye（Nu））\G_‘;
    %下面利用上面计算所得结果进行一步自校正GPC控制    
    %求Y0（k）
    deltau（k-1）=u（k-1）-u（k-2）;
    olddeltaU=[deltau（k-1:-1:k-nb）]‘;
    oldY=[y（k:-1:k-na）]‘;
    Y0=F_*oldY+H_*olddeltaU;
    %求Yd（k）
    %计算柔化轨迹
    yd（k）=y（k）;
    for j=1:N
        yd（k+j）=a*yd（k+j-1）+（1-a）*yr;
    end
    Yd=[yd（k+1:k+N）]‘;
    dd（k+1）=yd（k+1）;
    %求最终用到的△U，将△U的第一个u（k）作为当前时刻的输入改变量
    deltaU=K*（Yd-Y0）;
    u（k）=u（k-1）+deltaU（1）;
end

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     文件        3619  2012-10-28 15:34  GPC2\GPC_zijiaozheng.m
     文件         275  2012-09-24 09:21  GPC2\polyadd.m
     目录           0  2012-12-06 18:09  GPC2\