%廣義預(yù)測控制_自校正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;
%仿真步數(shù)
tim=50;
%預(yù)辨識次數(shù)
tim0=10;
%預(yù)測時域和控制時域
N=5;
Nu=4;
%擾動幅度
ampe=1;
%系數(shù)λ
lmd=0.5;
%y（k）的最終給定值yr
yr=50;
%柔化因子a
a=0.5;
%β
b=1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%參數(shù)預(yù)辨識
%將系統(tǒng)在正弦信號控制作用下運行tim0步，測的一組系統(tǒng)的運行數(shù)據(jù)
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
%根據(jù)所得運行數(shù)據(jù)，使用遞推最小二乘法辨識系統(tǒng)模型參數(shù)
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）;

%至此預(yù)辨識結(jié)束
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%一些計算的中間量
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];%△

%系統(tǒng)狀態(tài)的初值
y（1）=0;y（2）=0;
u（1）=0;u（2）=0;
e=ampe*（1-2*rand（1tim））;
%下面開始迭代進(jìn)行自校正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）;
????%根據(jù)新采集到的系統(tǒng)運行數(shù)據(jù)進(jìn)行一次遞推最小二乘法辨識
????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）;
????%一步辨識結(jié)束，得到參數(shù)
????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左右翻轉(zhuǎn)，因為MATLAB中的多項式系數(shù)是從高次到低次排列，而GPC算法中是相反的，這樣翻轉(zhuǎn)過來使用比較方便
????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
????%得出后邊計算要用到的系數(shù)矩陣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_‘;
????%下面利用上面計算所得結(jié)果進(jìn)行一步自校正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）作為當(dāng)前時刻的輸入改變量
????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\