背景

对于未知的非线性函数，仅通过函数的输入输出数据难以准确寻找函数极值。这类问题可以通过神经网络结合遗传算法求解，利用神经网络的非线性拟合能力和遗传算法的非线性寻优能力寻找函数极值。本文用神经网络遗传算法寻优如下非线性函数极值，函数表达式为

函数图形如下图1所示。

图1 非线性函数图形

从函数方程和图形可以看出，该函数的全局最小值为0，对应的坐标为(0,0)。虽然从函数方程和图形中很容易找出函数极值及极值对应坐标，但是在函数方程未知的情况下函数极值及极值对应坐标就很难找到。

模型建立

神经网络遗传算法函数极值寻优主要分为BP神经网络训练拟合和遗传算法极值寻优两步，算法流程如下图2所示。

图2 算法流程图

神经网络训练拟合根据寻优函数的特点构建合适的BP神经网络，用非线性函数的输出数据训练BP网络，训练后的BP神经网络就可以预测函数输出。遗传算法极值寻优把训练后的BP神经网络预测结果作为个体适应度值，通过选择、交叉和变异操作寻找函数的全局最优值及对应输入值。
本文根据非线性函数有2个输入参数、1个输出参数，确定BP神经网络结构为2-5-1.取函数的4 000组输入输出数据，从中随机选取3 900组数据训练网络，100组数据测试网络性能，网络训练好后用于预测非线性函数输出。
遗传算法中个体采用实数编码，由于寻优函数只有2个输入参数，所以个体长度为2。个体适应度值为BP神经网络预测值，适应度值越小。交叉概率为0.4，变异概率为0.2。

MATLAB实现

1、BP神经网络训练拟合

用函数输入输出数据训练BP神经网络，使训练后的网络能够拟合非线性函数输出，保存训练好的网络用语计算个体适应度值。根据非线性函数方程随机得到该函数的4 000组输入输出数据，存储于data.mat中，其中input为函数输入数据，output为函数对应输出数据，从中随机抽取3 900组训练数据训练网络，100组测试数据测试网络拟合性能。最后保存训练好的网络。

%% 清空环境变量
clc
clear

tic
%% 训练数据预测数据提取及归一化
%下载输入输出数据
load data1 input output

%从1到2000间随机排序
k=rand(1,4000);
[m,n]=sort(k);

%找出训练数据和预测数据
input_train=input(n(1:3900),:)';
output_train=output(n(1:3900),:)';
input_test=input(n(3901:4000),:)';
output_test=output(n(3901:4000),:)';

%选连样本输入输出数据归一化
[inputn,inputps]=mapminmax(input_train);
[outputn,outputps]=mapminmax(output_train);

%% BP网络训练
% %初始化网络结构
net=newff(inputn,outputn,5);

net.trainParam.epochs=100;
net.trainParam.lr=0.1;
net.trainParam.goal=0.0000004;

%网络训练
net=train(net,inputn,outputn);

%% BP网络预测
%预测数据归一化
inputn_test=mapminmax('apply',input_test,inputps);
 
%网络预测输出
an=sim(net,inputn_test);
 
%网络输出反归一化
BPoutput=mapminmax('reverse',an,outputps);

%% 结果分析

figure(1)
plot(BPoutput,':og')
hold on
plot(output_test,'-*');
legend('预测输出','期望输出','fontsize',12)
title('BP网络预测输出','fontsize',12)
xlabel('样本','fontsize',12)
ylabel('输出','fontsize',12)

%预测误差
error=BPoutput-output_test;

figure(2)
plot(error,'-*')
title('神经网络预测误差')

figure(3)
plot((output_test-BPoutput)./BPoutput,'-*');
title('神经网络预测误差百分比')

errorsum=sum(abs(error))

toc

2、适应度函数

把训练好的BP神经网络预测输出作为个体适应度值。

function fitness = fun(x)
% 函数功能：计算该个体对应适应度值
% x           input     个体
% fitness     output    个体适应度值

%
load data net inputps outputps

%数据归一化
x=x';
inputn_test=mapminmax('apply',x,inputps);
 
%网络预测输出
an=sim(net,inputn_test);
 
%网络输出反归一化
fitness=mapminmax('reverse',an,outputps);

3、遗传算法主函数

%% 清空环境变量
clc
clear

%% 初始化遗传算法参数
%初始化参数
maxgen=100;                         %进化代数，即迭代次数
sizepop=20;                        %种群规模
pcross=[0.4];                       %交叉概率选择，0和1之间
pmutation=[0.2];                    %变异概率选择，0和1之间

lenchrom=[1 1];          %每个变量的字串长度，如果是浮点变量，则长度都为1
bound=[-5 5;-5 5];  %数据范围


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);   %染色体的适应度
end
%找最好的染色体
[bestfitness bestindex]=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);
    
    % 计算适应度 
    for j=1:sizepop
        x=individuals.chrom(j,:); %解码
        individuals.fitness(j)=fun(x);   
    end
    
  %找到最小和最大适应度的染色体及它们在种群中的位置
    [newbestfitness,newbestindex]=min(individuals.fitness);
    [worestfitness,worestindex]=max(individuals.fitness);
    % 代替上一次进化中最好的染色体
    if bestfitness>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();
[r c]=size(trace);
plot([1:r]',trace(:,2),'r-');
title('适应度曲线','fontsize',12);
xlabel('进化代数','fontsize',12);
ylabel('适应度','fontsize',12);
axis([0,100,0,1])

% 窗口显示
disp('适应度');
x=bestchrom;
disp([bestfitness]);
disp('变量');
disp([x]);

4、结果分析

BP神经网络拟合结果分析
本文中个体的适应度值为BP神经网络预测值，因此BP神经网络预测精度对于最优位置的寻找具有非常重要的意义。由于寻优非线性函数有2个输入参数、1个输出参数，所以构建的BP神经网络的结构为2-5-1。共取非线性函数4 000组输入输出数据，从中随机选择3 900组数据训练BP神经网络，100组数据作为测试数据测试BP神经网络拟合性能，BP神经网络预测输出和期望输出对比如下图3所示。

图3 BP神经网络预测输出

从BP神经网络预测结果可以看出，BP神经网络可以准确预测非线性函数输出，可以把网络预测近似看成函数实际输出。

遗传算法寻优结果分析BP神经网络训练结束后，可以利用遗传算法寻找该非线性函数的最小值。遗传算法的迭代次数是100次，种群规模是20，交叉概率为0.4，变异概率为0.2，采用浮点数编码，个体长度为21，优化过程中最优个体适应度值变化曲线如下图4所示。

图4 适应度变化曲线

讨论

本文所使用的方法有比较重要的工程应用价值，比如对于某项试验来说，试验目的是获取到最大试验结果对应的实验条件，但是由于时间和经费限制，该试验只能进行有限次，可能单靠试验结果找不到最优的试验条件。这时可以在已知试验数据的基础上，通过本文介绍的神经网络遗传算法寻找最优试验条件。
思路就是先根据试验条件数和试验结果数确定BP神经网络结构；然后把试验条件作为输入数据，试验结果作为输出数据训练BP网络，使得训练后的网络可以预测一定试验条件下的试验结果；最后把试验条件作为遗传算法中的种群个体，把网络预测的试验结果作为个体适应度值，通过遗传算法推导最优试验结果及其对应试验条件。

本文训练数据及程序源代码见评论区

获取更多干货文章，请关注微信公众号【科学计算与编程】