📄 spso-bp-gann.m
字号:
%**************************************************************
% 胜华炼油厂常压塔产品质量软测量模型————PSOABPNN模型
%**************************************************************
clear
clc
% 装载数据
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';
%***********************************************
% 微粒群优化算法训练神经网络
%***********************************************
% 神经网络结构初始化(6-12-1结构)
[h,l]=size(PP);
II=l; % 输入层节点数
JJ=12; % 中间层节点数
% 初始化
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; % 位置初始化---连接权值
popul1=rand(JJ+1,PopSize)*6.0-3.0; % ----阈值
popul=[popul0',popul1']';
clear popul0; clear popul1;
vel=rand(dim,PopSize)*4.0-2.0; % 速度初始化
% 计算初始适应值
for m=1:PopSize, % 微粒群个数
e(m)=0;
for k=1:h % 数据组数
for j=1:JJ % 中间层节点数
net1(j)=0;
for i=1:II % 输入层节点数
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型函数---tansig 函数
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 函数
e(m)=e(m)+(TT(k)-y(k))^2;
end
ee(m)=e(m)/2;
end
ibestpos=popul; % 个体最好位置初始化
ibestfit=ee; % 各个体的适应值
[fbestpart,g]=min(ibestfit); % 找全局最好的适应值
gbestfit=fbestpart; % 全局最好的适应值
gbestpos=popul(:,g); % 全局最好的适应值对应的个体
g1=g;
change=0; % 标识符
stepcounter1=2; % 连续n步不动就重新启动
for k=1:stepcounter1 % 计数器清零
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; % 按余弦函数关系（减小）
for m=1:PopSize, % 将全局最好的适应值对应的个体展开
A(:,m)=gbestpos;
end
R1=rand(dim,PopSize); % 产生随机数
R2=rand(dim,PopSize);
if change==1,
popul(:,g1)=gbestpos;
end
vel=0.5*(w*vel+c1*R1.*(ibestpos-popul)+c2*R2.*(A-popul)); % 速度计算
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; % 位置计算
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, % 微粒群个数
e(m)=0;
for k=1:h % 数据组数
for j=1:JJ % 中间层节点数
net1(j)=0;
for i=1:II % 输入层节点数
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 函数
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 函数
e(m)=e(m)+(TT(k)-y(k))^2;
end
ee(m)=e(m)/2;
end
for m=1:PopSize, % 更新个体历史最好位置
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; % 历史全局最优适应度轨迹
if abs(gbestfit)<ErrGoal % 判断是否迭代结束
success=1;
end
end
gbestpos; % 输出最好个体及最好适应值
gbestfit;
iter;
% PSONN训练结束
%******************************************************
%******************************************************
%%% BP算法训练神经网络
%******************************************************
S1=JJ; % 中间层节点数
P=PP'; % 输入样本数据
T=TT'; % 输出样本数据
[w1,b1,w2,b2]=initff(P,S1,'tansig',T,'tansig'); % 初始化网络结构
df=100; % 训练参数设置
me=MaxIt; % 最大迭代次数
eg=ErrGoal; % 收敛误差限
lr=0.01; % % 学习速率
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训练过程
[w1,b1,w2,b2,ep,tr]=trainbp(w1,b1,'tansig',w2,b2,'tansig',P,T,tp); % BP训练过程
plottr(tr,eg); % 画训练过程曲线
yBP=simuff(P,w1,b1,'tansig',w2,b2,'tansig'); % 反算（仿真）输出结果
% BPNN训练结束
%****************************************************
%****************************************************
% 遗传算法求解
%****************************************************
% 参数设定
Size=PopSize; % 群体总个体数量
G=10; %MaxIt; % 允许最大迭代次数
CodeL=16; % 每个变量拥有的染色体数量
umax1=5.0; % 连接权值的最大值
umin1=-5.0; % 连接权值的最小值
umax2=2.5; % 阈值的最大值
umin2=-2.5; % 阈值的最小值
% 对所求参数进行初始二进制编码
E=round(rand(Size,((II+2)*JJ+1)*CodeL)); % 初始化个体代码（染色体）
% 遗传神经网络训练开始
bf1=0;
zj=0;
GAee00=100;
for k=1:G
k
% ********************** 对各个体进行解码、求适应值*********
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; % 连接权值解码
else
wb(ii,m)=(umax2-umin2)*yoy(ii,m)/(2^CodeL-1)+umin2; % 阈值解码
end
end
% 将解码后的值代入目标函数求适应值（一定要为正数） m
ee0(m)=0;
for kk=1:h % 数据组数
for j=1:JJ % 中间层节点数
net1(j)=0;
for i=1:II % 输入层节点数
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 函数
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 函数
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; % 保存历史最优
Bestwb=wb(:,Indexfi(Size));
end
end
if min(ee)<GAee00
GAee00=min(ee);
end
% ********* 选择函数 ****************
Ji=1./F;
BestJ(k)=min(Ji);
fi=F; % 适应值函数
[Oderfi,Indexfi]=sort(fi); % 从大到小排序
Bestfi=Oderfi(Size);
BestS=E(Indexfi(Size),:);
%****** 选择、产生下一代进行交叉操作 ******
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; % 交叉率设定(人为调试)
n=ceil(((II+2)*JJ+1)*CodeL*rand); % 随机产生交叉操作点
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; % 变异率设定(人为调试)
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
%*********** 保存最优个体 **************
TempE(Size,:)=BestS;
E=TempE;
end
% ****** 输出结果 **************
GAbestx=Bestwb; % 输出最优变量
%Bestbp=tr(MaxIt) % 输出最优适应值-BPNN
GABest_Value=bf1 % 输出最优适应值-GANN
gbestfit % 输出最优适应值-PSONN
GAee00
% 遗传算法求解结束
% ***********************************************
%************************************************
% 输出结果
%************************************************
% PSOANN计算最优结果（粗汽油干点值）
for k=1:h % 数据组数
for j=1:JJ % 中间层节点数
net1(j)=0;
for i=1:II % 输入层节点数
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函数
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函数
end
% GANN计算最优结果（粗汽油干点值）
for k=1:h % 数据组数
for j=1:JJ % 中间层节点数
net1(j)=0;
for i=1:II % 输入层节点数
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函数
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函数
end
%******************************************************
% 输出、显示、比较结果
%******************************************************
% 训练结果对比
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轴'); % 坐标标注
ylabel('y轴');
title('函数对比图：');
legend('实际值','PSONN预测值','GANN预测值','BPNN预测值');
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误差的绝对值平均值
SSE2=sqrt(SSE2/h) % GANN均方差
ABSE2=ABSE2/h % GANN误差的绝对值平均值
SSE3=sqrt(SSE3/h) % BPNN均方差
ABSE3=ABSE3/h % BPNN误差的绝对值平均值
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -