遗传算法求数值函数的最值

遗传算法求数值函数的最值

0. 引言

设有函数：

f(x) = x + 10*sin(5*x) + 7*cos(4*x);

其图像容易画出，如下所示：

先要求求该函数的最大值，读者可能已经有了很多种思路，本文介绍遗传算法是如何解决此类问题的。

1. 遗传算法简介

如果不关心算法的实现细节的话，遗传算法可以使用如下的流程描述。

这基本是借鉴生物种群的自然演化规律而抽象得到的流程图。下面分别解释流程图中的各个步骤。

编码

众所周知的，生物中都有保存其遗传信息的物质——染色体，染色体中的记录的遗传信息——基因决定了一个生物个体的表现性状，可以简单地认为一个染色体唯一确定一个生物个体。

染色体的本质是一种蛋白质，其结构非常复杂，用计算机模拟是不现实的。因此遗传算法简单地把染色体视为一串二进制位组成的序列（即该序列中的所有元素取自集合{0,1}）。

假设有命名为A的生物，规定其染色体序列长度为8，这时可以取得一个个体：

[0,1,0,0,1,0,1,1]

初始化种群

显然一个生物种群不可能只包含一个个体，因此可以规定种群的大小，为方便说明原理，这是规定其大小为4。初始种群可以随机生成。如下：

[1,1,1,1,1,1,0,0]
[1,1,1,1,0,1,1,0]
[1,0,1,0,0,0,0,1]
[0,0,0,1,0,1,0,1]

评估种群中的个体适应度

生物种群的存续与兴衰决定于其对于环境的适应能力。对于使用二进制位序列表示的个体，总可以找出一个函数来评估这个个体对环境的适应能力。可以规定函数：

f(x) = int(x)
# x是输入的一个个体，如[0,1,0,0,1,0,1,1]
# 返回值是该二进制序列表示的整数

此外还可以规定上述函数的返回值（一个整数）越大，个体的适应度越高。显然适应度最高的个体是：

[1,1,1,1,1,1,1,1]

选择

选择操作用于选出一个种群中存活概率更高的个体，这些个体组成一个新的种群。显然个体适应度越高的个体存活概率更高。因此可以使用赌轮盘算法选出新种群。赌轮盘算法的原理不是本文的重点，因此略去，其实现代码见后文。

交叉

交叉操作也是对生物遗传过程的抽象。设有两个生物个体，这两个个体的生殖繁衍过程中，其染色体会进行交叉交换，生成新的染色体，繁衍得来的新的个体携带新的染色体。

设有两个生物个体如下：

[1,1,1,1,1,1,0,0]
[1,1,1,1,0,1,1,0]

设对这两个生物个体有一定的概率发生交叉操作，而交叉操作必然发生在二进制位序列中的某一位置。这一位置的选择也是随机的。假设从第5位（从1开始计数）开始交叉，其步骤如下：

[1,1,1,1,1,1,0,0]
[1,1,1,1,0,1,1,0]
->
[1,1,1,1,1],[1,0,0]
[1,1,1,1,0],[1,1,0]
->
[1,1,1,1,1],[1,1,0]
[1,1,1,1,0],[1,0,0]
->
[1,1,1,1,1,1,1,0]
[1,1,1,1,0,1,0,0]

变异

生物在自然环境中受环境的影响，有可能发生遗传信息的永久改变，也就是变异。遗传算法规定，二进制位序列的变异指的是在一定的概率下，随机选择该序列中的一位取反。假设在第2位发生了变异，步骤如下：

[0,1,0,0,1,0,1,1]
->
[0],1,[0,0,1,0,1,1]
->
[0],0,[0,0,1,0,1,1]
->
[0,0,0,0,1,0,1,1]

小结

上述的各个步骤按照流程图所示依次执行，直到一个种群中所有的个体都相同（即收敛），如：

[1,0,1,0,0,0,0,1]
[1,0,1,0,0,0,0,1]
[1,0,1,0,0,0,0,1]
[1,0,1,0,0,0,0,1]

或者主循环的迭代次数超过规定次数，算法停止，输出结果。

Matlab入门

本文中的算法使用Matlab描述，因此需要对Matlab有一定的了解。但是帮读者入门Matlab不是我们的任务，因此推荐以下手册，可以帮助读者快速掌握Matlab：
https://ww2.mathworks.cn/help/pdf_doc/matlab/getstart_zh_CN.pdf

Matlab是收费软件，而且对于运行该软件的性能有一定要求，因此可以使用开源免费且占用资源较少的Octave作为替代品，其官网如下：
https://www.gnu.org/software/octave/

Octave使用的语法规则与Matlab基本相同，不兼容之处见下：
https://en.wikibooks.org/wiki/MATLAB_Programming/Differences_between_Octave_and_MATLAB

本文所附的代码在Octave上运行通过，并未在Matlab上进行测试，如果有需要使用Matlab运行代码，需要自己检查兼容性。

算法实现

function retval = main (input1, input2)
  % 主函数，input1,input2,以及retval都没有实际作用
  clear all;
  close all;
  clc;
  
  NP = 50; % 种群大小
  L = 16; % 染色体序列长度
  Pc = 0.8; % 交叉概率
  Pm = 0.1; % 变异概率
  G = 50; % 主循环迭代次数
  Xs = 10; % 所求函数的定义域的最大值
  Xx = 0; % 所求函数的定义域的最小值
  S = 10; % 什么意思自己悟吧，我只是不想硬编码成10
  
  pop = initPop(NP,L);
  
  for i = 1:G
    dpop = decode(pop);
    xpop = getX(dpop,L,Xx,Xs);
    fit = fitValue(xpop);
    if all(fit == fit(1)) % 如果收敛，停止主循环
      break;
    endif
    [maxFit,maxIndex] = max(fit);
    maxX(i) = getX(decode(pop(maxIndex,:)),L,Xx,Xs);
    maxY(i) = maxFit;
    pop = select(pop,fit,xpop);    
    pop = crossOver(pop,Pc);    
    pop = mutation(pop,Pm);
  endfor
  
  
  x = Xx:(Xs-Xx)/(2^L-1)/S:Xs;
  y = targetFunc(x);  
  subplot(2,1,1);
  plot(x,y); % 绘制函数图像
  hold on;
  plot(maxX,maxY,'r*'); % 绘制历次迭代生成的最大值点
  
  [px,py] = size(maxY);
  x = 1:py;
  y = maxY;  
  subplot(2,1,2);
  plot(x,y); % 绘制历次生成的最大值点关于迭代次数的折线图
  hold off;
  
endfunction


function pop = initPop(NP,L)
  pop = round(rand(NP,L));
endfunction

function res = targetFunc(x)
  res = x + 10*sin(5*x) + 7*cos(4*x);
endfunction

function res = decode(pop)
  % 用于将二进制位序列转化成整数
  [px,py] = size(pop);
  res = ones(px,1);
  for j = 1:py
    res(:,j) = 2.^(py-j)*pop(:,j);
  endfor
  res = sum(res,2);
endfunction

function res = getX(dpop,L,Xx,Xs)
  % 缩放二进制位序列表示成的整数到
  % 定义域中
  res = dpop./(2^L-1).*(Xs-Xx);
endfunction

function res = fitValue(xpop)
  res = targetFunc(xpop);
endfunction

function newpop = select(pop,fit,xpop)
  % 赌轮盘算法的实现
  [maxFit,maxIndex] = max(fit);
  [minFit,minIndex] = min(fit);
  [px,py] = size(pop);
  newpop = ones(px,py);
  fit = (fit-minFit)./(maxFit-minFit);
  fit = fit./sum(fit);
  cumfit = cumsum(fit);
  ms = sort(rand(px,1));
  fiti = 1;
  newi = 1;
  while newi <= px
    if ms(newi) < cumfit(fiti)
      newpop(newi,:) = pop(fiti,:);
      newi = newi + 1;
    else
      fiti = fiti + 1;
    endif
  endwhile
endfunction

function newpop = crossOver(pop,Pc)
  % 交叉操作
  [px,py] = size(pop);
  newpop = ones(size(pop));
  for i = 1:2:px-1
    if rand<Pc
      cpoint = ceil(rand*py);
      newpop(i,:) = [pop(i,1:cpoint),pop(i+1,cpoint+1:py)];
      newpop(i+1,:) = [pop(i+1,1:cpoint),pop(i,cpoint+1:py)];
     else
      newpop(i,:) = pop(i,:);
      newpop(i+1,:) = pop(i+1,:);
    endif
  endfor
endfunction

function newpop = mutation(pop,Pm)
  % 变异操作
  [px,py] = size(pop);
  newpop = ones(px,py);
  for i = 1:px
    if rand<Pm
      mpoint=ceil(rand*py);
      newpop(i,mpoint) = 1 - pop(i,mpoint);
      newpop(i,:) = pop(i,:);
    else
      newpop(i,:) = pop(i,:);
    endif
  endfor
endfunction

绘图结果如下：