bp_pso.m

用pso优化bp初始权值

%混合bp算法，即pso_bp算法：
%算法思路：利用pso的全局搜索能力对bp网络的权值进行优化，再利用bp算法对其进行局部寻优。
%pso算法部分如下：
%1.1 初始化格式
clear all
clc;
format long;
%1.2 初始化微粒群的参数
wmax=1.2;
wmin=0.6;
N = 40 ; %微粒数目
%N=60;
D = 16 ; %每个微粒的空间维数
a=-1;
b=1;
c1=2.5; %学习因子1
%c1=2;
c2=2.5;  %学习因子2
%c2=2;
%w=1.2; %惯性权重
MaxDT=20 ; %最大迭代次数
%rand('state',sum(100*clock));
%x=a+(b-a)*rand(N,D);%[-1,1]
%v=a+(b-a)*rand(N,D);%[-1,1]
x=randn(N,D);  %[0,1]
v=randn(N,D);  %[0,1]

%初始化神经网络
P=-0.5:0.05:0.5;
T=(sin(2*3.14*P)).^2.*exp(-4*(P+0.5));
I=1;H=5;O=1;%1-5-1结构
net = newff([-0.5 0.5],[H O],{'tansig' 'purelin'});

%计算各个粒子的适应度值
for i=1:N
    p(i) = fitcall(x(i,:),net,I,H,O,D,P,T); %局部最优值
    y(i,:) = x(i,:);
end
pg = x(1,:); %pg为全局最优值
for i=2:N
    if fitcall(x(i,:),net,I,H,O,D,P,T)<fitcall(pg,net,I,H,O,D,P,T)
        pg = x(i,:);
    end
end

%迭代公式
for t = 1:MaxDT
    w = wmin + (wmax-wmin)*(1+cos((t-1)*pi/(MaxDT-1)))/2;
    for i = 1:N
        v(i,:)=w*v(i,:)+c1*rand*(y(i,:)-x(i,:))+c2*rand*(pg-x(i,:));
        x(i,:)=x(i,:)+v(i,:);
        if fitcall(x(i,:),net,I,H,O,D,P,T)<p(i)
            p(i) = fitcall(x(i,:),net,I,H,O,D,P,T);
            y(i,:) = x(i,:);
        end
        if fitcall(x(i,:),net,I,H,O,D,P,T)<fitcall(pg,net,I,H,O,D,P,T)
            pg = y(i,:);
        end
    end

    Pbest(t) = fitcall(pg,net,I,H,O,D,P,T);
end

%bp网络
%clear all;
%clc;
%p = -0.5:0.05:0.5;
%t = (sin(2*3.14*p)).^2.*exp(-4*(p+0.5));
%net=newff(minmax(p),[5 1],{'transig' 'purelin'},'trainlm');
%net = newff([-0.5 0.5],[5 1],{'tansig' 'purelin'});
%net = initnw(net,1);
for t = 1:H
    x2iw(t,:) = pg(1,((t-1)*I+1):t*I);%wi为5*1的矩阵
end
for r = 1:O
    x2lw(r,:) = pg(1,(I*H+1):(I*H+H));%wo为1*5的矩阵
end

x2b=pg(1,((I+1)*H+1):D);
x2b1=x2b(1:H).';
x2b2=x2b(H+1:H+O).';
net.iw{1,1}=x2iw(:,1);
net.lw{2,1}=x2lw;
net.b{1,1}=x2b1;
net.b{2,1}=x2b2;

net.trainParam.epoch = 500; %训练次数为500次；
net.trainParam.goal = 0.0001; %训练目标为误差小于0.001;
net.trainParam.lr = 0.01; %学习速率为0.01;
net = train(net,P,T);
y1 = sim(net,P);
figure(1);
%plot(p,y1,'-b')
%plot(p,t,'-r',p,y1,'-b','LineWidth',2);
plot(P,T,'r*');
hold on;
plot(P,y1,'-b');
axis([-0.5 0.5 -0.05 0.45]);
title('非线性函数逼近');
xlabel('输出');
ylabel('时间');