遗传算法优化的BP神经网络建模

遗传算法优化的BP神经网络建模。

目标：

   对y=x1^2+x2^2非线性系统进行建模，用1500组数据对网络进行构建网络，500组数据测试网络。由于BP神经网络初始神经元之间的权值和阈值一般随机选择，因此容易陷入局部最小值。本方法使用遗传算法优化初始神经元之间的权值和阈值，并对比使用遗传算法前后的效果。

步骤：

未经遗传算法优化的BP神经网络建模

1、 随机生成2000组两维随机数(x1,x2)，并计算对应的输出y=x1^2+x2^2，前1500组数据作为训练数据input_train，后500组数据作为测试数据input_test。并将数据存储在data中待遗传算法中使用相同的数据。

2、  数据预处理：归一化处理。

3、  构建BP神经网络的隐层数，次数，步长，目标。

4、 使用训练数据input_train训练BP神经网络net。

5、 用测试数据input_test测试神经网络，并将预测的数据反归一化处理。

6、  分析预测数据与期望数据之间的误差。

遗传算法优化的BP神经网络建模

1、  读取前面步骤中保存的数据data；

2、  对数据进行归一化处理；

3、  设置隐层数目；

4、  初始化进化次数，种群规模，交叉概率，变异概率

5、 对种群进行实数编码，并将预测数据与期望数据之间的误差作为适应度函数；

6、 循环进行选择、交叉、变异、计算适应度操作，直到达到进化次数，得到最优的初始权值和阈值；

7、  将得到最佳初始权值和阈值来构建BP神经网络；

8、 使用训练数据input_train训练BP神经网络net;

9、 用测试数据input_test测试神经网络，并将预测的数据反归一化处理；

10、             分析预测数据与期望数据之间的误差。

算法流程图如下：

运行后使用遗传算法改进前后误差的对比图：

程序：

1、未经遗传算法优化的BP神经网络建模

clear;

clc;

%%%%%%%%%%%%%输入参数%%%%%%%%%%%%%%

N=2000;           %数据总个数

M=1500;           %训练数据

%%%%%%%%%%%%%训练数据%%%%%%%%%%%%%%

for i=1:N

   input(i,1)=-5+rand*10;

   input(i,2)=-5+rand*10;

end

output=input(:,1).^2+input(:,2).^2;

save data input output

 

load data.mat

 

%从1到N随机排序

k=rand(1,N);

[m,n]=sort(k);

%找出训练数据和预测数据

input_train=input(n(1:M),:)';

output_train=output(n(1:M),:)';

input_test=input(n((M+1):N),:)';

output_test=output(n((M+1):N),:)';

%数据归一化

[inputn,inputs]=mapminmax(input_train);

[outputn,outputs]=mapminmax(output_train);

%构建BP神经网络

net=newff(inputn,outputn,5);

net.trainParam.epochs=100;

net.trainParam.lr=0.1;

net.trainParam.goal=0.0000004;

%BP神经网络训练

net=train(net,inputn,outputn);

%测试样本归一化

inputn_test=mapminmax('apply',input_test,inputs);

%BP神经网络预测

an=sim(net,inputn_test);

%%网络得到数据反归一化

BPoutput=mapminmax('reverse',an,outputs);

 

figure(1)

%plot(BPoutput,':og');

scatter(1:(N-M),BPoutput,'rx');

hold on;

%plot(output_test,'-*');

scatter(1:(N-M),output_test,'o');

legend('预测输出','期望输出','fontsize',12);

title('BP网络预测输出','fontsize',12);

xlabel('样本','fontsize',12);

xlabel('优化前输出的误差','fontsize',12);

 

figure(2)

error=BPoutput-output_test;

plot(1:(N-M),error);

xlabel('样本','fontsize',12);

ylabel('优化前输出的误差','fontsize',12);

%save net net inputs outputs

2、遗传算法优化的BP神经网络建模

（1）主程序

%清空环境变量

clc

clear

 

%读取数据

load data.mat

 

%节点个数

inputnum=2;

hiddennum=5;

outputnum=1;

 

%训练数据和预测数据

input_train=input(1:1500,:)';

input_test=input(1501:2000,:)';

output_train=output(1:1500)';

output_test=output(1501:2000)';

 

%选连样本输入输出数据归一化

[inputn,inputps]=mapminmax(input_train);

[outputn,outputps]=mapminmax(output_train);

 

%构建网络

net=newff(inputn,outputn,hiddennum);

 

%% 遗传算法参数初始化

maxgen=10;                        %进化代数，即迭代次数

sizepop=30;                       %种群规模

pcross=[0.3];                      %交叉概率选择，0和1之间

pmutation=[0.1];                   %变异概率选择，0和1之间

 

%节点总数

numsum=inputnum*hiddennum+hiddennum+hiddennum*outputnum+outputnum;

 

lenchrom=ones(1,numsum);       

bound=[-3*ones(numsum,1)3*ones(numsum,1)];   %数据范围

 

%------------------------------------------------------种群初始化------------------------------%------------------

 

--------

individuals=struct('fitness',zeros(1,sizepop),'chrom',[]);  %将种群信息定义为一个结构体

%avgfitness=[];                     %每一代种群的平均适应度

bestfitness=[];                    %每一代种群的最佳适应度

bestchrom=[];                      %适应度最好的染色体

%初始化种群

for i=1:sizepop

   %随机产生一个种群

   individuals.chrom(i,:)=Code(lenchrom,bound);   %编码

   x=individuals.chrom(i,:);

    %计算适应度

   individuals.fitness(i)=fun(x,inputnum,hiddennum,outputnum,net,inputn,outputn);  %染色体的适应度

end

 

%找最好的染色体

[bestfitnessbestindex]=min(individuals.fitness);

bestchrom=individuals.chrom(bestindex,:); %最好的染色体

%avgfitness=sum(individuals.fitness)/sizepop;%染色体的平均适应度

% 记录每一代进化中最好的适应度和平均适应度

%trace=[avgfitness bestfitness];

 

%% 迭代求解最佳初始阀值和权值

% 进化开始

for i=1:maxgen

    i

    % 选择

   individuals=Select(individuals,sizepop);

   %avgfitness=sum(individuals.fitness)/sizepop;

    %交叉

   individuals.chrom=Cross(pcross,lenchrom,individuals.chrom,sizepop,bound);

    % 变异

   individuals.chrom=Mutation(pmutation,lenchrom,individuals.chrom,sizepop,i,maxgen,bound);

   

    % 计算适应度

    forj=1:sizepop

       x=individuals.chrom(j,:); %解码

       individuals.fitness(j)=fun(x,inputnum,hiddennum,outputnum,net,inputn,outputn);  

    end

   

  %找到最小和最大适应度的染色体及它们在种群中的位置

   [newbestfitness,newbestindex]=min(individuals.fitness);

   [worestfitness,worestindex]=max(individuals.fitness);

    %代替上一次进化中最好的染色体

    ifbestfitness>newbestfitness

       bestfitness=newbestfitness;

       bestchrom=individuals.chrom(newbestindex,:);

    end

   individuals.chrom(worestindex,:)=bestchrom;

   individuals.fitness(worestindex)=bestfitness;

   

   %avgfitness=sum(individuals.fitness)/sizepop;

   

   %trace=[trace;avgfitness bestfitness]; %记录每一代进化中最好的适应度和平均适应度

 

end

%% 遗传算法结果分析

 %figure(3)

%[r c]=size(trace);

%plot([1:r]',trace(:,2),'b--');

%title(['适应度曲线  ' '终止代数＝'num2str(maxgen)]);

%xlabel('进化代数');ylabel('适应度');

%legend('平均适应度','最佳适应度');

disp('适应度                  变量');

x=bestchrom;

 

%% 把最优初始阀值权值赋予网络预测

% %用遗传算法优化的BP网络进行值预测

w1=x(1:inputnum*hiddennum);

B1=x(inputnum*hiddennum+1:inputnum*hiddennum+hiddennum);

w2=x(inputnum*hiddennum+hiddennum+1:inputnum*hiddennum+hiddennum+hiddennum*outputnum);

B2=x

 

(inputnum*hiddennum+hiddennum+hiddennum*outputnum+1:inputnum*hiddennum+hiddennum+hiddennum*outputnum+outputnum);

 

net.iw{1,1}=reshape(w1,hiddennum,inputnum);

net.lw{2,1}=reshape(w2,outputnum,hiddennum);

net.b{1}=reshape(B1,hiddennum,1);

net.b{2}=B2;

 

%% BP网络训练

%网络进化参数

net.trainParam.epochs=100;

net.trainParam.lr=0.1;

%net.trainParam.goal=0.00001;

 

%网络训练

[net,per2]=train(net,inputn,outputn);

 

%% BP网络预测

%数据归一化

inputn_test=mapminmax('apply',input_test,inputps);

an=sim(net,inputn_test);

test_simu=mapminmax('reverse',an,outputps);

error=test_simu-output_test;

 

%figure(4);

hold on;plot(1:500,error,'r');

legend('优化前的误差','优化后的误差','fontsize',12)

 

(2)编码子程序code.m

function ret=Code(lenchrom,bound)

%本函数将变量编码成染色体，用于随机初始化一个种群

% lenchrom  input : 染色体长度

%bound     input : 变量的取值范围

%ret       output: 染色体的编码值

flag=0;

while flag==0

   pick=rand(1,length(lenchrom));

   ret=bound(:,1)'+(bound(:,2)-bound(:,1))'.*pick;%线性插值，编码结果以实数向量存入ret中

   flag=test(lenchrom,bound,ret);    %检验染色体的可行性

end

 

(3)适应度函数fun.m

function error =fun(x,inputnum,hiddennum,outputnum,net,inputn,outputn)

%该函数用来计算适应度值

%x         input    个体

%inputnum  input    输入层节点数

%outputnum input    隐含层节点数

%net       input    网络

%inputn    input    训练输入数据

%outputn   input    训练输出数据

%error     output   个体适应度值

%提取

w1=x(1:inputnum*hiddennum);

B1=x(inputnum*hiddennum+1:inputnum*hiddennum+hiddennum);

w2=x(inputnum*hiddennum+hiddennum+1:inputnum*hiddennum+hiddennum+hiddennum*outputnum);

B2=x(inputnum*hiddennum+hiddennum+hiddennum*outputnum+1:inputnum*hiddennum+hiddennum+hiddennum*outputnum+outputnum);

net=newff(inputn,outputn,hiddennum);

%网络进化参数

net.trainParam.epochs=20;

net.trainParam.lr=0.1;

net.trainParam.goal=0.00001;

net.trainParam.show=100;

net.trainParam.showWindow=0;

%网络权值赋值

net.iw{1,1}=reshape(w1,hiddennum,inputnum);

net.lw{2,1}=reshape(w2,outputnum,hiddennum);

net.b{1}=reshape(B1,hiddennum,1);

net.b{2}=B2;

%网络训练

net=train(net,inputn,outputn);

an=sim(net,inputn);

error=sum(abs(an-outputn));

 

(4)选择操作Select.m

function ret=select(individuals,sizepop)

% 该函数用于进行选择操作

% individualsinput   种群信息

%sizepop    input   种群规模

%ret        output   选择后的新种群

 

%求适应度值倒数  

[a bestch]=min(individuals.fitness);

%b=individuals.chrom(bestch);

%c=individuals.fitness(bestch);

fitness1=10./individuals.fitness;%individuals.fitness为个体适应度值

 

%个体选择概率

sumfitness=sum(fitness1);

sumf=fitness1./sumfitness;

 

%采用轮盘赌法选择新个体

index=[];

for i=1:sizepop  %sizepop为种群数

   pick=rand;

    whilepick==0   

       pick=rand;       

    end

    fori=1:sizepop   

       pick=pick-sumf(i);       

       ifpick<0       

           index=[indexi];           

           break; 

       end

    end

end

%index=[index bestch];

%新种群

individuals.chrom=individuals.chrom(index,:);  %individuals.chrom为种群中个体

individuals.fitness=individuals.fitness(index);

%individuals.chrom=[individuals.chrom;b];

%individuals.fitness=[individuals.fitness;c];

ret=individuals;

 

(5)交叉操作cross.m

functionret=Cross(pcross,lenchrom,chrom,sizepop,bound)

%本函数完成交叉操作

%pcorss               input  : 交叉概率

%lenchrom             input  : 染色体的长度

%chrom    input  : 染色体群

%sizepop              input  : 种群规模

%ret                  output : 交叉后的染色体

 for i=1:sizepop %每一轮for循环中，可能会进行一次交叉操作，染色体是随机选择的，交叉位置也是随机选择的，%但该轮for循环中是否进行交叉操作则由交叉概率决定（continue控制）

    % 随机选择两个染色体进行交叉

    pick=rand(1,2);

    while prod(pick)==0

        pick=rand(1,2);

    end

    index=ceil(pick.*sizepop);

    % 交叉概率决定是否进行交叉

    pick=rand;

    while pick==0

        pick=rand;

    end

    if pick>pcross

        continue;

    end

    flag=0;

    while flag==0

        % 随机选择交叉位

        pick=rand;

        while pick==0

            pick=rand;

        end

        pos=ceil(pick.*sum(lenchrom));%随机选择进行交叉的位置，即选择第几个变量进行交叉，注意：两个染色体交叉的位置相同

        pick=rand; %交叉开始

        v1=chrom(index(1),pos);

        v2=chrom(index(2),pos);

        chrom(index(1),pos)=pick*v2+(1-pick)*v1;

        chrom(index(2),pos)=pick*v1+(1-pick)*v2; %交叉结束

        flag1=test(lenchrom,bound,chrom(index(1),:)); %检验染色体1的可行性

        flag2=test(lenchrom,bound,chrom(index(2),:)); %检验染色体2的可行性

        if   flag1*flag2==0

            flag=0;

        else flag=1;

        end   %如果两个染色体不是都可行，则重新交叉

    end

 end

ret=chrom;

 

(6)变异操作Mutation.m

functionret=Mutation(pmutation,lenchrom,chrom,sizepop,num,maxgen,bound)

% 本函数完成变异操作

%pcorss               input  : 变异概率

%lenchrom             input  : 染色体长度

%chrom    input  : 染色体群

%sizepop              input  : 种群规模

%opts                 input  : 变异方法的选择

%pop                  input  : 当前种群的进化代数和最大的进化代数信息

%bound                input  : 每个个体的上届和下届

%maxgen               input  ：最大迭代次数

%num                  input  : 当前迭代次数

%ret                  output : 变异后的染色体

for i=1:sizepop  %每一轮for循环中，可能会进行一次变异操作，染色体是随机选择的，变异位置也是随机选择的，

   %但该轮for循环中是否进行变异操作则由变异概率决定（continue控制）

    %随机选择一个染色体进行变异

   pick=rand;

    whilepick==0

       pick=rand;

    end

   index=ceil(pick*sizepop);

    %变异概率决定该轮循环是否进行变异

   pick=rand;

    ifpick>pmutation

       continue;

    end

    flag=0;

    whileflag==0

       % 变异位置

       pick=rand;

       whilepick==0     

           pick=rand;

       end

       pos=ceil(pick*sum(lenchrom)); %随机选择了染色体变异的位置，即选择了第pos个变量进行变异

       pick=rand;%变异开始    

       fg=(rand*(1-num/maxgen))^2;

       if pick>0.5

           chrom(i,pos)=chrom(i,pos)+(bound(pos,2)-chrom(i,pos))*fg;

       else

           chrom(i,pos)=chrom(i,pos)-(chrom(i,pos)-bound(pos,1))*fg;

       end   %变异结束

       flag=test(lenchrom,bound,chrom(i,:));    %检验染色体的可行性

    end

end

ret=chrom;