写在前面：随着研究问题的复杂化，有的时候会有一些寻找最优参数的科研任务，这就需要我们会使用一些优化方法，而遗传算法则是众多优化算法中最基础和最常见的方法。本次分享依旧属于科普贴，拟分两期给初次接触优化算法的同学介绍遗传算法，第一期借助一个遗传算法的案例代码，对遗传算法的原理进行简单介绍；第二期借助GA工具箱，对遗传算法程序进行优化与简化，以提高遗传算法运行效率和查错效率。

遗传算法

我们都知道进化论的基本规则就是“自然选择，适者生存”，每个个体通过交配繁衍出后代，后代在成长过程中经过自然法则的筛选而淘汰掉部分个体，而那些拥有优良性状的个体得以存活并且通过繁衍使得优秀遗传信息得以保留和扩散。

遗传算法（Genetic Algorithm, GA）正是借鉴了这种思想。我们把生物的进化过程抽象成一个数学模型，生物进化指的是初始种群里拥有各种各样的同种生物，经过繁衍和自然选择得到的最终种群里，每个个体或者说绝大部分是最适应环境的优秀个体。那么对应的数学模型则可简化成，初始集合中包含一定数量的随机数，而这些随机数通过一定规则反复迭代计算，然后设计一种个体的评价体系，把不符合规则或者评分比较低的个体淘汰，得到的最终集合里绝大部分是最符合规则的数。

有了生物进化和对应的抽象数学模型之后，我们用一个图示来简要说明遗传算法中的一些术语和流程。

下面我们用一个小案例以及代码，来详细展示遗传算法的流程以及代码的书写。

问题：（二元函数在给定区间的最值）

原函数在区间上的图像如图所示：

syms x y

x=-1:0.01:12;

y=4:0.01:6;

[X,Y]=meshgrid(x,y);

f=X.*sin(pi.*X)+Y.*sin(exp(Y));

mesh(X,Y,f);

可以看出，函数是一个在x方向上震荡，同时在y方向上震荡的图形，并且该函数是一个超越函数，用常规求最大值方法比如求导或者离散方法求解起来十分复杂。

下面我们用遗传算法来求解该问题。

首先，我们需要设置初始种群的大小还有每个个体遗传信息的编码（coding）与解码（decoding）。初始种群大小可以根据实际问题自由设定，较大种群包含的遗传信息更多，但是迭代起来更慢。仿照DNA的结构，DNA含有四种碱基，这四种碱基可以对应生成人体所有的蛋白质，所以一个DNA序列可以认为是一种四进制编码。数学中常见的进制编码有二进制编码，八进制编码，十进制编码，十六进制编码，计算机的底层存储通常是二进制编码，所以这里我们选用二进制编码比较合适，这样计算起来会快很多。那么我们需要进行编码与解码的过程，编码则是把性状变为遗传信息，而解码则是把遗传信息变为性状。以x变量为例，x变量的范围为[-1,12]，区间长度为13，我们需要的精度为小数点后四位，那么，通过进制数之间的转化

那么x变量的二进制编码长度至少为17位，同理若是y变量的精度为小数点后四位，那么y变量所需要的二进制编码长度至少为15位

那么同时包含变量x，y的编码总共需要的二进制编码长度为32位。这32位二进制编码则是生物学中的遗传信息，而对应的性状则是变量x，y的十进制值，所以解码过程则是把二进制信息转换成十进制性状。由此可以看出，遗传算法中的编码与解码，实际上就是数学中不同进制数之间的相互转换。解码过程不再详细赘述，这里只给出一个公式：，x的二进制编码为，则长度为n+1位，那么十进制数为

设置好了种群大小和编码与解码过程，我们可以写出种群设置代码。

%% 初始种群

M = 40;N =32;

pop = round(rand(M,N)); %随机0-1矩阵生成，产生初始种群遗传信息

x1 = decode_x1(pop(:,1:17)); %对x1和x2进行解码

x2 = decode_x2(pop(:,18:end));

子函数decode_x1和decode_x2类似，主要实现二进制数和十进制数之间的转换

function [x1]=decode_x1(code)

[M,N] = size(code);

sum = zeros(N);

for i=1:M

   for j=1:N

   sum(i)=sum(i)+code(i,j)*2^(N-j);

   end

 x2(i) = 4.0 + sum(i)*(12.0 - (-4.0))/(2^N - 1);
end

得到初始种群的性状之后，我们要对初始种群进行定量评价，来确定种群中个体之间的适应性（fitness），通常就是我们研究问题的目标，

则有代码

fitness = x1.*sin(pi*x1) +x2.*sin(exp(x2)); %适应度

那么fitness大的个体就是我们所需要寻找的，小的个体则被淘汰，所以后续代码则实现如何产生并选择fitness大的个体。

我们设计一个迭代规则，主要是为了实现优秀亲代如何产生子代的过程，其中包括选择，交互和变异。

首先看看选择（selection）过程，这里是为了选择出优秀亲代，让优秀亲代之间产生子代。通常所使用的方法是赌轮盘选择，该方法简单，容易编写代码，但是选择误差较大。

比如我们有三条染色体，他们所对应的适应度评分分别为5，15，20，那么累积总适应度为，所以，每个个体被选中的概率为

那么想象有这样一个轮盘，中间有三个区域对应这三种个体，转一次轮盘，指针所指个体即为所选中个体。

为了实现赌轮盘的选择规则，我们写出如下赌轮盘子函数：

%赌轮盘子函数

function [B]=seclection(fitness,A)

   [M,N]=size(A);

   fit_sum = sum(fitness);

   for i = 1:M

     probability(i) = fitness(i)/fit_sum;

   end

   for i = 2:M

     probability(i) =probability(i) + probability(i-1);
  
   end

   for j = 1:M

     p = rand();

     i = 1;

     while p > probability(i)

       i = i+1;

     end
  
     B(j,:) = A(i,:);
    
   end

end

然后，来看看优秀亲代之间的交互（crossover）过程。通过生物学常识，类比于染色体的重组与交互，我们知道亲代中父本和母本遗传信息的交互生成了子代的遗传信息。交互过程中亲代的选择与遗传信息的具体结合方式有很多种，可以自由设计，这里我们采取一种最简单的方式，即种群编号相邻的两个个体为亲代的父本和母本，并且随机生成一个二进制编码点位，该点位到编码最后一位即为交互片段。为了实现该过程，我们可以写出交互子函数：

%交互子函数

function [newpop]=crossover(pc,A)

  newpop=A;

  [M,N]=size(A);

   for i= 1:2:M-1

     if rand < pc

       p = ceil(rand()*N);

         for j= p:N

           ch = A(i,j);

           newpop(i,j) = A(i+1,j);

           newpop(i+1,j) = ch;

         end

       end

    end

end

这里的pc则是交互概率，根据实际问题提前设定好pc的值，本例中pc选为0.65。

最后则是变异（mutation）过程，变异过程最主要的意义在于避免在迭代中陷入局部最大值附近，也就是函数极大值附近，就像本例中，x方向和y方向上反复震荡，有许多的峰值，如果没有变异过程，那么该算法所得到的结果容易进入极值点附近，而非最大值附近。正常情况下遗传信息的变异很复杂，包含有替换，缺失，增添三类，为了简化变异过程，针对二进制编码的特殊性，变异过程简易设计为替换——在某概率下，原编码中的‘0’变为‘1’，‘1’变为‘0’。那么为了实现变异过程，我们可以写出变异子函数：

%变异子函数

function [newpop]=mutation(pm,A)

  [M,N]=size(A);

  newpop=A;

  W = rand(M,N)<pm;

   for i=1:M

     for j=1:N

       if W(i,j) ==1

         if A(i,j)~=1

           newpop(i,j)= 1;

         else

           newpop(i,j) = 0;

         end

        end

       end

    end

end

其中pm则是变异概率，同样的，根据实际问题提前设置好pm的值，在本例中，pm值选为0.08。

**写好了选择，交互和变异过程的子函数，我们可以来写遗传算法迭代主体代码了。**

%% 迭代主题代码

gen = 1;

while gen < generation

   [B] = seclection(fitness,pop); %赌轮盘子函数

    [newpop] = crossover(pc,B); %交互子函数

    [B] = mutation(pm,newpop); %变异子函数

    pop = B;

   x1 = decode_x1(pop(:,1:17));

   x2 = decode_x2(pop(:,18:end));

   fitness = x1.*sin(pi*x1) +x2.*sin(exp(x2)); %该代子代适应度计算

   [fit_best,index] = max(fitness);%该代子代中最优个体适应度以及对应的遗传信息

  %所有子代中最优个体以及对应性状（变量值）寻找

   if fit_best >= maxdata

     maxdata = fit_best;

     verter = pop(index,:);

     x_1 = decode_x1(verter(1:17));

     x_2 = decode_x2(verter(18:end));

   end

  num(gen) = maxdata; %每一代子代中最优个体集合

   gen = gen + 1;

end

遗传算法预先可以设置迭代停止条件，比如本代最优个体与所有子代中最优个体之间差值小于某个范围则停止迭代，这里我们采用的是最简单的停止方法，当代数进行到一定值时，停止迭代，即generation，预先设置为200代。初始最优适应度maxdata可以任意选择，在本例中，maxdata选为-2。

剩下的，则是对运行程序得到的结果进行处理，遗传算法通常需要的结果大致有两个：

1）最优个体（变量x和y的值）以及最优适应度（二元函数的最大值）；

2）遗传算法每一代适应度变化过程。

本例中，结果处理的代码为

%% 结果与画图

disp(sprintf('max f=£º%.6f',num(gen-1)));%输出最优解

disp(sprintf('x_1=£º%.5f x_2=£º%.5f',x_1,x_2));%最优解对应的变量值

figure(1)  

plot(num,'k');%画图

所得到结果为

这里变量的结果精度为小数点后五位，并不是原先设定的四位，原因是在解码时候产生的。

小结一下：

遗传算法模拟了生物的进化过程，看起来很麻烦，但实际上理清生物进化的几个过程之后，分步写出对应各个过程的子函数，遗传算法代码也是比较容易书写的。虽是如此，有的同学仍然觉着遗传算法代码好长啊，找bug时候比较困难，没关系，下一期会介绍遗传算法工具箱的使用方法，十几行的代码就可以写完一个遗传算法程序。

最后，为了加强遗传算法的应用，我们来看一个实际的优化问题。姜锐老师研究传统车辆的行驶之后提出了全速度差模型

从模型上来看，这是一个确定性模型，并且有两个不同的敏感度系数和 ，我们把全速度差模型设计成自动车的上层模型，那么一个自动车车队每辆车都可以认为是同质的，不存在传统车辆中不同司机对于敏感度系数和 的影响，模型中参数变量只有和 。

那么，我们如何用遗传算法去寻找这样的参数和 ：

在保证自动车车队满足系统稳定性的前提之下，使得车队的总能耗最小？

我们是一个有灵魂的团队，坚持探索，致力于分享交流学习经验。

想获取更多交通建模，论文写作，开源资料等科研信息的小伙伴就请关注

微信公众号【交通科研Lab】 （所有信息均在公众号第一时间发布）

交通科研Lab.png