?? spso-bp-gann.m
字號(hào):
%**************************************************************
% 勝華煉油廠常壓塔產(chǎn)品質(zhì)量軟測(cè)量模型————PSOABPNN模型
%**************************************************************
clear
clc
% 裝載數(shù)據(jù)
PP=-1:0.1:1;
PP=PP';
TT=[-0.9602 -0.577 -0.0729 0.3771 0.6405 0.66 0.4609...
0.1336 -0.2013 -0.4344 -0.5 -0.393 -0.1647 0.0988...
0.3072 0.396 0.3449 0.1816 -0.0312 -0.2189 -0.3201];
TT=TT';
%***********************************************
% 微粒群優(yōu)化算法訓(xùn)練神經(jīng)網(wǎng)絡(luò)
%***********************************************
% 神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)初始化(6-12-1結(jié)構(gòu))
[h,l]=size(PP);
II=l; % 輸入層節(jié)點(diǎn)數(shù)
JJ=12; % 中間層節(jié)點(diǎn)數(shù)
% 初始化
success=0; PopSize=20; MaxIt=100; iter=0; ErrGoal=0.001;
dim=(II+1)*JJ+1; maxw=1.8; minw=0.01;
c1=2; c2=2; w=maxw;
popul0=rand((II+1)*JJ,PopSize)*16.0-8.0; % 位置初始化---連接權(quán)值
popul1=rand(JJ+1,PopSize)*6.0-3.0; % ----閾值
popul=[popul0',popul1']';
clear popul0; clear popul1;
vel=rand(dim,PopSize)*4.0-2.0; % 速度初始化
% 計(jì)算初始適應(yīng)值
for m=1:PopSize, % 微粒群個(gè)數(shù)
e(m)=0;
for k=1:h % 數(shù)據(jù)組數(shù)
for j=1:JJ % 中間層節(jié)點(diǎn)數(shù)
net1(j)=0;
for i=1:II % 輸入層節(jié)點(diǎn)數(shù)
net1(j)=net1(j)+PP(k,i)*popul(j+(i-1)*JJ,m);
end
net1(j)=net1(j)+popul((II+1)*JJ+j,m);
o1(j)=(1-exp(-net1(j)))/(1+exp(-net1(j))); % S型函數(shù)---tansig 函數(shù)
end
net2=0;
for j=1:JJ
net2=net2+o1(j)*popul(j+II*JJ,m);
end
net2=net2+popul((II+2)*JJ+1,m);
y(k)=(1-exp(-net2))/(1+exp(-net2));
% y(k)=1/(1+exp(-net2)); % logsig 函數(shù)
e(m)=e(m)+(TT(k)-y(k))^2;
end
ee(m)=e(m)/2;
end
ibestpos=popul; % 個(gè)體最好位置初始化
ibestfit=ee; % 各個(gè)體的適應(yīng)值
[fbestpart,g]=min(ibestfit); % 找全局最好的適應(yīng)值
gbestfit=fbestpart; % 全局最好的適應(yīng)值
gbestpos=popul(:,g); % 全局最好的適應(yīng)值對(duì)應(yīng)的個(gè)體
g1=g;
change=0; % 標(biāo)識(shí)符
stepcounter1=2; % 連續(xù)n步不動(dòng)就重新啟動(dòng)
for k=1:stepcounter1 % 計(jì)數(shù)器清零
for i=1:PopSize
velsign(k,i)=1;
velsign0(k,i)=1;
end
end
while(success==0)&(iter<MaxIt), % 迭代開始
iter=iter+1
w=minw+(maxw-minw)*(1+cos((iter-1)*pi/(MaxIt-1)))/2.0; % 按余弦函數(shù)關(guān)系（減小）
for m=1:PopSize, % 將全局最好的適應(yīng)值對(duì)應(yīng)的個(gè)體展開
A(:,m)=gbestpos;
end
R1=rand(dim,PopSize); % 產(chǎn)生隨機(jī)數(shù)
R2=rand(dim,PopSize);
if change==1,
popul(:,g1)=gbestpos;
end
vel=0.5*(w*vel+c1*R1.*(ibestpos-popul)+c2*R2.*(A-popul)); % 速度計(jì)算
clear A;
for d=1:dim % 速度限幅處理
for m=1:PopSize
if vel(d,m)>4;,
vel(d,m)=4;
end
if vel(d,m)<-4,
vel(d,m)=-4;
end
end
end
for i=1:PopSize
distance1=0;
for d=1:dim
distance1=distance1+vel(d,i).^2;
end
distance1=sqrt(distance1);
if distance1<0.000005
velsign(stepcounter1,i)=0;
elseif distance1>1000
velsign0(stepcounter1,i)=0;
else
velsign(stepcounter1,i)=1;
velsign0(stepcounter1,i)=1;
end
for k=1:(stepcounter1-1)
velsign(k,i)=velsign(k+1,i);
velsign0(k,i)=velsign0(k+1,i);
end
distance1=0;
distance0=0;
for k=1:stepcounter1
distance1=distance1+velsign(k,i);
distance0=distance0+velsign0(k,i);
end
if distance1<0.5
vel(:,i)=rands(dim,1);
end
if distance0<0.5
vel(:,i)=rands(dim,1);
end
end
popul=popul+vel; % 位置計(jì)算
for m=1:PopSize % 位置限
for j=1:JJ % b1
if popul((II+1)*JJ+j,m)>2.8
popul((II+1)*JJ+j,m)=2.8*rand; %rand*1.5;
elseif popul((II+1)*JJ+j,m)<-2.8
popul((II+1)*JJ+j,m)=-2.8*rand; %rand*(-1.5);
end
end
if popul((II+2)*JJ+1,m)>2.5, % b2
popul((II+2)*JJ+1,m)=2.5*rand;
elseif popul((II+2)*JJ+1,m)<-2.5
popul((II+2)*JJ+1,m)=-2.5*rand;
end
end
for m=1:PopSize % w
for j=1:JJ
for i=1:II+1
if popul((i-1)*JJ+j,m)>8
popul((i-1)*JJ+j,m)=8*rand;
elseif popul((i-1)*JJ+j,m)<-8
popul((i-1)*JJ+j,m)=-8*rand;
end
end
end
end
for m=1:PopSize, % 微粒群個(gè)數(shù)
e(m)=0;
for k=1:h % 數(shù)據(jù)組數(shù)
for j=1:JJ % 中間層節(jié)點(diǎn)數(shù)
net1(j)=0;
for i=1:II % 輸入層節(jié)點(diǎn)數(shù)
net1(j)=net1(j)+PP(k,i)*popul(j+(i-1)*JJ,m);
end
net1(j)=net1(j)+popul((II+1)*JJ+j,m);
o1(j)=(1-exp(-net1(j)))/(1+exp(-net1(j))); % tansig 函數(shù)
end
net2=0;
for j=1:JJ
net2=net2+o1(j)*popul(j+II*JJ,m);
end
net2=net2+popul((II+2)*JJ+1,m);
y(k)=(1-exp(-net2))/(1+exp(-net2));
% y(k)=1/(1+exp(-net2)); % logsig 函數(shù)
e(m)=e(m)+(TT(k)-y(k))^2;
end
ee(m)=e(m)/2;
end
for m=1:PopSize, % 更新個(gè)體歷史最好位置
if ee(m)<ibestfit(m),
ibestfit(m)=ee(m);
ibestpos(:,m)=popul(:,m);
end
end
[fbestpart,g]=min(ee); % 更新全局歷史最好位置
if fbestpart<gbestfit;
gbestfit=fbestpart;
gbestpos=popul(:,g);
g1=g;
change=0;
else
change=1;
end
[fbestpart,g]=max(gbestfit); % 更新全局歷史最壞位置----sPSOA
popul(:,g)=gbestpos;
seg0(iter)=gbestfit; % 歷史全局最優(yōu)適應(yīng)度軌跡
if abs(gbestfit)<ErrGoal % 判斷是否迭代結(jié)束
success=1;
end
end
gbestpos; % 輸出最好個(gè)體及最好適應(yīng)值
gbestfit;
iter;
% PSONN訓(xùn)練結(jié)束
%******************************************************
%******************************************************
%%% BP算法訓(xùn)練神經(jīng)網(wǎng)絡(luò)
%******************************************************
S1=JJ; % 中間層節(jié)點(diǎn)數(shù)
P=PP'; % 輸入樣本數(shù)據(jù)
T=TT'; % 輸出樣本數(shù)據(jù)
[w1,b1,w2,b2]=initff(P,S1,'tansig',T,'tansig'); % 初始化網(wǎng)絡(luò)結(jié)構(gòu)
df=100; % 訓(xùn)練參數(shù)設(shè)置
me=MaxIt; % 最大迭代次數(shù)
eg=ErrGoal; % 收斂誤差限
lr=0.01; % % 學(xué)習(xí)速率
tp=[df me eg lr];
%tp=[df me eg];
%[w1,b1,w2,b2,ep,tr]=trainlm(w1,b1,'tansig',w2,b2,'logsig',P,T,tp); % BP訓(xùn)練過程
[w1,b1,w2,b2,ep,tr]=trainbp(w1,b1,'tansig',w2,b2,'tansig',P,T,tp); % BP訓(xùn)練過程
plottr(tr,eg); % 畫訓(xùn)練過程曲線
yBP=simuff(P,w1,b1,'tansig',w2,b2,'tansig'); % 反算（仿真）輸出結(jié)果
% BPNN訓(xùn)練結(jié)束
%****************************************************
%****************************************************
% 遺傳算法求解
%****************************************************
% 參數(shù)設(shè)定
Size=PopSize; % 群體總個(gè)體數(shù)量
G=10; %MaxIt; % 允許最大迭代次數(shù)
CodeL=16; % 每個(gè)變量擁有的染色體數(shù)量
umax1=5.0; % 連接權(quán)值的最大值
umin1=-5.0; % 連接權(quán)值的最小值
umax2=2.5; % 閾值的最大值
umin2=-2.5; % 閾值的最小值
% 對(duì)所求參數(shù)進(jìn)行初始二進(jìn)制編碼
E=round(rand(Size,((II+2)*JJ+1)*CodeL)); % 初始化個(gè)體代碼（染色體）
% 遺傳神經(jīng)網(wǎng)絡(luò)訓(xùn)練開始
bf1=0;
zj=0;
GAee00=100;
for k=1:G
k
% ********************** 對(duì)各個(gè)體進(jìn)行解碼、求適應(yīng)值*********
for m=1:1:Size
mm=E(m,:); % 取要解碼的染色體
for ii=1:(II+2)*JJ+1 % 清零
yoy(ii,m)=0;
end
% ********* 解碼操作 ******************
for ii=1:(II+2)*JJ+1
m1=mm((ii-1)*CodeL+1:1:ii*CodeL);
for i=1:1:CodeL
yoy(ii,m)=yoy(ii,m)+m1(i)*2^(i-1);
end
if ii<((II+1)*JJ+1)
wb(ii,m)=(umax1-umin1)*yoy(ii,m)/(2^CodeL-1)+umin1; % 連接權(quán)值解碼
else
wb(ii,m)=(umax2-umin2)*yoy(ii,m)/(2^CodeL-1)+umin2; % 閾值解碼
end
end
% 將解碼后的值代入目標(biāo)函數(shù)求適應(yīng)值（一定要為正數(shù)） m
ee0(m)=0;
for kk=1:h % 數(shù)據(jù)組數(shù)
for j=1:JJ % 中間層節(jié)點(diǎn)數(shù)
net1(j)=0;
for i=1:II % 輸入層節(jié)點(diǎn)數(shù)
net1(j)=net1(j)+PP(kk,i)*wb(j+(i-1)*JJ,m);
end
net1(j)=net1(j)+wb((II+1)*JJ+j,m);
o1(j)=(1-exp(-net1(j)))/(1+exp(-net1(j))); % tansig 函數(shù)
end
net2=0;
for j=1:JJ
net2=net2+o1(j)*wb(j+II*JJ,m);
end
net2=net2+wb((II+2)*JJ+1,m);
y(kk)=(1-exp(-net2))/(1+exp(-net2));
% y(kk)=1/(1+exp(-net2)); % logsig 函數(shù)
ee0(m)=ee0(m)+(TT(kk)-y(kk))^2;
end
ee(m)=ee0(m)/2;
F(m)=1./ee(m);
end
[Oderfi,Indexfi]=sort(F); % 從大到小排序
Bestfi=Oderfi(Size);
if k>1
if Bestfi>bf1
bf1=Bestfi; % 保存歷史最優(yōu)
Bestwb=wb(:,Indexfi(Size));
end
end
if min(ee)<GAee00
GAee00=min(ee);
end
% ********* 選擇函數(shù) ****************
Ji=1./F;
BestJ(k)=min(Ji);
fi=F; % 適應(yīng)值函數(shù)
[Oderfi,Indexfi]=sort(fi); % 從大到小排序
Bestfi=Oderfi(Size);
BestS=E(Indexfi(Size),:);
%****** 選擇、產(chǎn)生下一代進(jìn)行交叉操作 ******
fi_sum=sum(fi);
fi_Size=(Oderfi/fi_sum)*Size;
fi_S=floor(fi_Size);
kk=1;
for i=1:1:Size
for j=1:1:fi_S(i)
TempE(kk,:)=E(Indexfi(i),:);
kk=kk+1;
end
end
%************ 交叉操作 ************
pc=0.75; % 交叉率設(shè)定(人為調(diào)試)
n=ceil(((II+2)*JJ+1)*CodeL*rand); % 隨機(jī)產(chǎn)生交叉操作點(diǎn)
for i=1:2:(Size-1)
temp=rand;
if pc>temp
for j=n:1:(((II+2)*JJ+1)*CodeL)
TempE(i,j)=E(i+1,j);
TempE(i+1,j)=E(i,j);
end
end
end
TempE(Size,:)=BestS;
E=TempE;
%************ 變異操作 **************
pm=0.1; % 變異率設(shè)定(人為調(diào)試)
for i=1:Size
for j=1:(((II+2)*JJ+1)*CodeL)
temp=rand;
if pm>temp
if TempE(i,j)==0
TempE(i,j)=1;
else
TempE(i,j)=0;
end
end
end
end
%*********** 保存最優(yōu)個(gè)體 **************
TempE(Size,:)=BestS;
E=TempE;
end
% ****** 輸出結(jié)果 **************
GAbestx=Bestwb; % 輸出最優(yōu)變量
%Bestbp=tr(MaxIt) % 輸出最優(yōu)適應(yīng)值-BPNN
GABest_Value=bf1 % 輸出最優(yōu)適應(yīng)值-GANN
gbestfit % 輸出最優(yōu)適應(yīng)值-PSONN
GAee00
% 遺傳算法求解結(jié)束
% ***********************************************
%************************************************
% 輸出結(jié)果
%************************************************
% PSOANN計(jì)算最優(yōu)結(jié)果（粗汽油干點(diǎn)值）
for k=1:h % 數(shù)據(jù)組數(shù)
for j=1:JJ % 中間層節(jié)點(diǎn)數(shù)
net1(j)=0;
for i=1:II % 輸入層節(jié)點(diǎn)數(shù)
net1(j)=net1(j)+PP(k,i)*gbestpos(j+(i-1)*JJ);
end
net1(j)=net1(j)+gbestpos((II+1)*JJ+j);
o1(j)=(1-exp(-net1(j)))/(1+exp(-net1(j))); % tansig函數(shù)
end
net2=0;
for j=1:JJ
net2=net2+o1(j)*gbestpos(j+II*JJ);
end
net2=net2+gbestpos((II+2)*JJ+1);
yPSO(k)=(1-exp(-net2))/(1+exp(-net2));
% yPSO(k)=1/(1+exp(-net2)); % logsig函數(shù)
end
% GANN計(jì)算最優(yōu)結(jié)果（粗汽油干點(diǎn)值）
for k=1:h % 數(shù)據(jù)組數(shù)
for j=1:JJ % 中間層節(jié)點(diǎn)數(shù)
net1(j)=0;
for i=1:II % 輸入層節(jié)點(diǎn)數(shù)
net1(j)=net1(j)+PP(k,i)*Bestwb(j+(i-1)*JJ);
end
net1(j)=net1(j)+Bestwb((II+1)*JJ+j);
o1(j)=(1-exp(-net1(j)))/(1+exp(-net1(j))); % tansig函數(shù)
end
net2=0;
for j=1:JJ
net2=net2+o1(j)*Bestwb(j+II*JJ);
end
net2=net2+Bestwb((II+2)*JJ+1);
yGA(k)=(1-exp(-net2))/(1+exp(-net2));
% yGA(k)=1/(1+exp(-net2)); % logsig函數(shù)
end
%******************************************************
% 輸出、顯示、比較結(jié)果
%******************************************************
% 訓(xùn)練結(jié)果對(duì)比
figure
plot(1:h,TT(1:h),'-k',1:h,yPSO(1:h),'--b',1:h,yGA(1:h),'-.r',1:h,yBP(1:h),':g')
xlabel('x軸'); % 坐標(biāo)標(biāo)注
ylabel('y軸');
title('函數(shù)對(duì)比圖：');
legend('實(shí)際值','PSONN預(yù)測(cè)值','GANN預(yù)測(cè)值','BPNN預(yù)測(cè)值');
SSE1=0;
ABSE1=0;
SSE2=0;
ABSE2=0;
SSE3=0;
ABSE3=0;
for k=1:h
SSE1=SSE1+(yPSO(k)-TT(k))^2;
ABSE1=ABSE1+abs(yPSO(k)-TT(k));
SSE2=SSE2+(yGA(k)-TT(k))^2;
ABSE2=ABSE2+abs(yGA(k)-TT(k));
SSE3=SSE3+(yBP(k)-TT(k))^2;
ABSE3=ABSE3+abs(yBP(k)-TT(k));
end
SSE1=sqrt(SSE1/h) % PSONN均方差
ABSE1=ABSE1/h % PSONN誤差的絕對(duì)值平均值
SSE2=sqrt(SSE2/h) % GANN均方差
ABSE2=ABSE2/h % GANN誤差的絕對(duì)值平均值
SSE3=sqrt(SSE3/h) % BPNN均方差
ABSE3=ABSE3/h % BPNN誤差的絕對(duì)值平均值
?? 快捷鍵說明
復(fù)制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號(hào)
Ctrl + =
減小字號(hào)
Ctrl + -