遗传算法（2）

在遗传算法中我们再举一个求极大值的例子。这种例子也是比较多见的，只要我们把一些数据关系描述成函数之后就会有一些求极大值或者极小值的问题。

其实极大值和极小值是一类问题，就是极值问题，解题思路也是一般无二。那我们信手拈来看一个函数。我们假设在空间里有一个函数z=ysin(x)+xcos(y)，图像如下所示。

如果觉得这个图的意义找不到对应的物理解释也没关系，我们就假设这是一个地区当中的地形图，而且地形每个点的(x,y,z)三维坐标是可以用z=ysin(x)+xcos(y)这个函数来进行描述的，山峦起伏是不是？如果我们想求这个函数在x位于[-10,10]和y位于[10,10]之间的最大值怎么求呢？

对于微积分学得好的朋友可能就根本不是问题了。把z对y求偏导数，把z对x求偏导数，这样能够求出满足两个偏导数同时为零的x,y解，很可能有多个，把所有这一组一组的x,y代入z=ysin(x)+xcos(y)，求出的(x,y,z)就是驻点。每个驻点的z值比大小，最大的就是要求的解了。没学过微积分的朋友也先别着急，我们今天介绍的不是这种微积分领域常用的办法，还是考虑用遗传算法的思路来做。

在刚刚这个山峦叠嶂的小小区域，我们再冒充一把上帝。

可以这么想象，我在这个地区无规律地区撒下我创造的很多子民，有的在谷底，有的在半山腰，有的可能已经在山顶或者山顶附近。那么下面让他们一代一代生生不息地繁殖，凡是能爬得更高的就留下，爬不上去的……那就让他们在另一个平行世界生活吧，我们就不观察了。有了这个思路我们就还是按部就班地来试着做做看。

1、基因编码

首先我们遇到的基因编码问题就和背包问题不一样。刚刚我们说的背包问题属于离散型的自变量，一个物品在背包里要么有要么没有，即便是128个物品，我们也知道一共是有2128种情况，基因编码我们最多用128bit也就够了。

但是这个例子中，我们有两个自变量x,y，这两个自变量还都是实数。根据实数稠密性原理，在[-10,10]之间是有无穷多个实数的，我们没有办法用类似背包问题中的基因编码的办法穷举所有的实数可能性。

说到这里我们还是不要惊慌，不知道大家注意过没有，电子计算机在做所有实数运算的时候其实是没有那么精确的——因为本身用电子计算机来计算实数就是一个用离散解决连续，用有穷解决无穷的方案，本身就有天生的不可逾越的局限性。我们不管用64bit来描述一个实数还是用128bit描述一个实数都太有限了，所以我们需要一个精度限制，就是我们平时说的有效数字。问题是我们应该取多少位有效数字？

理想来说我们对于有效数字应该是尽可能多取，有多少取多少，但是在实际生活中却不这么做。原因也很简单，取多有效数字本来是为了提高精确度，降低误差与成本的过程。然而取多有效数字同样需要更多的成本，而多出的有效数字的增长对提高收益如果没有明显的好处，那显然取太多有效数字反而是不划算的。

我们在平时购物的时候取到元小数点后2位就够了，因为再小的货币单位也不发行，再往后标记更小的标价单位毫无意义。一些大宗交易的物品单价通常会取的比较多，比如像小米手机、iphone手机中的电子元件——电容、放大器之类的，很多出厂单价才0.005元一个甚至更低。但是由于购买量巨大，一批就是以千万个原件作为单位的购买量，1000万个电子元件那就是50000人民币了，这种情况下取4位小数就有意义。再比如生产CPU用的高纯硅通常是99.999999999%（11个9）的，纯度低了无法支持22nm级别的生产工艺。太阳能板硅板就不用这么精细了，一般4个9就能满足，如果超过6个9只能是增加成本。

所以对于那些我们无法感知精确的，感知了也无法把控的，把控了也对提高效益无意义的，这些场合下的精度提升就没必要了，在生产生活中大家酌情进行取舍就好了。

再回来说我们刚刚这个例子，我们如果允许一定的误差，只要误差足够小其实就已经能够满足模型的需求了。

假设从-10到10之间指的是一个20公里的地带，那么取精度为1米可以理解为[-10.000,10.000]这个区间范围，取精度0.1米就可以理解为[-10.0000,10.0000]的范围，我们就假设1米的精度已经能够满足误差要求了，那么我们取[-10.000,10.000]作为取值范围，这样无限就变成有限了，这里实际是20001种取值。

如果我们还是要用二进制来进行基因编码应该选用多少位呢？214是16384，215是32768，这个20001种情形。那么我们就取15bit作为基因描述吧，x和y都用15bit的基因——因为这已经能够覆盖所有的设想样本。

那么下一个问题就是把一个15bit的二进制数字映射到[-10.000,10.000]这个区间里去，难吗？一点都不难，我们用第一bit作为正负数的标识，剩下14bit作为具体数字的标识，这样这个15bit的数字就成为了[-16383,0]和[0,16383]两个区间取并集了，注意集合里的元素都是整数。再把每个数字÷16383÷1000，取值范围就被压缩在[-10.000,10.000]之间了。

注意在这里为了让区间包括两边的边界点，我们做一个小小的变换：

我们定义F(x)的内容如下：x是自变量，二进制数字，首位为0代表正数，首位为1代表负数，函数值为x对应的十进制数字。那么我们就完成了编码到整个定义区间上的映射F(x)，最后表示二进制到实数值的函数H(x)如下定义：

细心的朋友可能看到这个地方有问题了，这个地方有两个问题。

问题1：F(x)>=0和F(x)<=0，是有重叠部分的，F(x)究竟算谁的。那就补充说明一下，这个写法确实是不够严谨，具体意义是指x的首位是0还是1的情况，是0那就算F(x)>=0，是1那就算F(x)<=0，因为确实后面14bit都是0的情况下，不管首位是0还是1代表值都是0。

问题2:10*(F(x)+1)/16384和10*(F(x)-1)/16384分别构成了x正负两个部分，但是区间却不是我们原来说的[-10.000,10.000]了，变成了[-10.000,0)和(0,10.000]，0从自变量范围里被拿掉了。我承认这是为了计算方便而已，不过我确信0肯定不是我们要的解，不管是x=0还是y=0，所以从解里面把所有x=0的情况以及所有y=0的情况都去除了，这个变换不影响最终求解。如果你觉得不放心，那么可以尝试重新构造这个映射关系，把[-10.000,10.000]所有的点都覆盖到。

最后距离说明映射关系的具体值：

(000000100100100)2=292，对应0.179，(111000100110110)2=-12598，对应-7.690。

怎么样，容易吧，用现在流行的词儿来说那就是——“完美”。

2、设计初始群体

冒充上帝真不是一件容易的事情，因为设计初始群体是我们的第二个要解决的大问题。

我们倾斜一下看这个小世界，我们要设置一些初始化的子民让他们一代一代繁衍后代，能爬得更高的就继续观察，爬不高的就不管了。在我们刚刚设计的基因里其实有两条基因，一条是x一条是y，这两条基因各有15个基因信息点，也就是215个可能性，我们随机产生8组。

3、适应度函数

适应度在这个场景里倒是不难设计，就用z=ysin(x)+xcos(y)就可以，z就是适应度。

4、产生下一代

在这个场景里，我们在每一代都是可以让1个染色体中的基因X之间和基因Y之间进行组合。

如表格所示，如果染色体1和2进行结合，那么

1染色体的X基因的前7位和2染色体的X基因的后8位将结合；

1染色体的Y基因的前7位和2染色体的Y基因的后8位将结合；

2染色体的X基因的前7位和1染色体的X基因的后8位将结合；

2染色体的Y基因的前7位和1染色体的Y基因的后8位将结合。

由此形成XA、YA、XB、YB四个后代基因，XA和YA将成为新的一组染色体，XB和YB将成为新的一组染色体，形成两个完整的后代基因染色体组。

如果是8组作为初始种群的大小的话，这样就有C8,2=28种组合方式，而每一种组合产生2个后代，那么实际上第一代以后就产生56个个体。

这56个个体的适应度我们可以拿来排序，只取出排名前8的个体。

这里同样可以允许一定的基因突变性，我们在8个已遴选的个体中，随机找到两个个体的，让这两个个体其中一个x染色体发生变异而让另一个y染色体发生变异（这个例子碰巧了由x和y两个自变量构成的染色体跟人类的性染色体同名了）。变异也是随机改变x染色体中的某一位，或随机改变y染色体中的某一位。

这之后再进行两两重组的计算，产生下一代。

这里注意两点：

1、断开点的位置

理论上讲断开点的位置是任意的，只不过，断开点靠左对数值影响变化大，自变量“跳跃”范围也就大；断开点靠右对数值影响变化小，自变量“跳跃”范围也就小。

2、基因变异的位置

和断开点位置的影响是完全一样的，同样是变异点靠左自变量“跳跃”范围大，变异点靠右自变量“跳跃”范围小。

5、迭代计算

迭代计算的过程我们直接还是看代码和执行过程就好了。请注意，由于其中有很多随机的因素，所以我们计算的过程结果可能会不一致，迭代的代数也可能不一致，但是最终结果应该都是一样的。

#coding=utf-8

import random

import math

import numpy as np

#染色体长度

CHROMOSOME_SIZE = 15

#判断退出

def is_finished(last_three):
    s =sorted(last_three)
    if s[0]and s[2] - s[0] < 0.01 * s[0]:
       return True
    else:
       return False

#初始染色体样态

def init():
   chromosome_state1 = ['000000100101001', '101010101010101']
   chromosome_state2 = ['011000100101100', '001100110011001']
   chromosome_state3 = ['001000100100101', '101010101010101']
   chromosome_state4 = ['000110100100100', '110011001100110']
   chromosome_state5 = ['100000100100101', '101010101010101']
   chromosome_state6 = ['101000100100100', '111100001111000']
    chromosome_state7= ['101010100110100', '101010101010101']
   chromosome_state8 = ['100110101101000', '000011110000111']
   chromosome_states = [chromosome_state1,
                    chromosome_state2,
                    chromosome_state3,
                     chromosome_state4,
                    chromosome_state5,
                    chromosome_state6,
                    chromosome_state7,
                    chromosome_state8]
    returnchromosome_states

#计算适应度

def fitness(chromosome_states):
    fitnesses= []
    forchromosome_state in chromosome_states:
        ifchromosome_state[0][0] == '1':
            x= 10 * (-float(int(chromosome_state[0][1:], 2) - 1)/16384)
        else:
            x= 10 * (float(int(chromosome_state[0], 2) + 1)/16384)
        ifchromosome_state[1][0] == '1':
            y= 10 * (-float(int(chromosome_state[1][1:], 2) - 1)/16384)
        else:
            y= 10 * (float(int(chromosome_state[1], 2) + 1)/16384)
        z = y* math.sin(x) + x * math.cos(y)
       fitnesses.append(z)
    returnfitnesses

#筛选

def filter(chromosome_states, fitnesses):
    #top 8 对应的索引值
   chromosome_states_new = []
   top1_fitness_index = 0
    for i innp.argsort(fitnesses)[::-1][:8].tolist():
       chromosome_states_new.append(chromosome_states[i])
       top1_fitness_index = i
    returnchromosome_states_new, top1_fitness_index

#产生下一代

def crossover(chromosome_states):
   chromosome_states_new = []
    whilechromosome_states:
       chromosome_state = chromosome_states.pop(0)
        for vin chromosome_states:
           pos = random.choice(range(8, CHROMOSOME_SIZE - 1))
           chromosome_states_new.append([chromosome_state[0][:pos] + v[0][pos:],chromosome_state[1][:pos] + v[1][pos:]])
           chromosome_states_new.append([v[0][:pos] + chromosome_state[1][pos:],v[0][:pos] + chromosome_state[1][pos:]])
    returnchromosome_states_new

#基因突变

def mutation(chromosome_states):
    n =int(5.0 / 100 * len(chromosome_states))
    while n> 0:
        n -=1
       chromosome_state = random.choice(chromosome_states)
        index= chromosome_states.index(chromosome_state)
        pos =random.choice(range(len(chromosome_state)))
        x =chromosome_state[0][:pos] + str(int(not int(chromosome_state[0][pos]))) +chromosome_state[0][pos+1:]
        y =chromosome_state[1][:pos] + str(int(not int(chromosome_state[1][pos]))) +chromosome_state[1][pos+1:]
       chromosome_states[index] = [x, y]
if __name__ == '__main__':
   chromosome_states = init()
    last_three= [0] * 3
    last_num= 0
    n = 100
    while n> 0:
        n -=1
       chromosome_states = crossover(chromosome_states)
       mutation(chromosome_states)
       fitnesses = fitness(chromosome_states)
       chromosome_states, top1_fitness_index = filter(chromosome_states,fitnesses)
       last_three[last_num] = fitnesses[top1_fitness_index]
        ifis_finished(last_three):
           break
        iflast_num >= 2:
           last_num = 0
        else:
           last_num += 1

#x: 7.69897460938 y:7.69897460938 z:8.79528923825

#x: -8.35693359375 y:-8.35693359375 z:11.3501994249

这里我们最后收敛条件是判断每一代的适应函数最大值做以记录，判断最近3代的最大值如何变化。如果最近3代的适应函数最大值相比较，第一大的（最大的）比第三大的（最小的）增益小于1%，那么就判断为收敛。这种收敛的速度非常快，大概3到4代就可以收敛完毕求出最优解。当然，这个部分我们同样可以再讨论采用更为考究一些的计算方式，这个作为有兴趣的朋友继续研究的提示好了。

请注意这里有一点，大家在执行这段代码中可能会发现这么一个问题。有的时候z值会收敛到8.8附近，此时x和y都在7.7附近，就像我们例子中给的这个数字；有的时候z值会收敛到11.35附近，此时x和y都在-8.35附近。从这两种不同的解上来看，我们就可以知道后者应该是正确的而绝非前者。

在图上我们看一下，圆圈圈住的位置就是(7.7, 7.7)附近的点，而实际上正确的解在方框圈住的位置(-8.35, -8.35)。为什么会出现这种现象呢？其实也不难理解，就是在繁殖下一代的时候又出现一次播撒的情况，但是播撒不够均匀，导致部分爬到“圆圈山”位置的种群个体表现格外良好，其他种群比如在“方框山”半山腰的这些就没能被遴选——凡是第9名开外的都算淘汰。这样是有相当的概率会收敛到8.8附近的“圆圈山”的。这一类的问题可能我们以后在写遗传算法中也同样会遇到，请大家注意。

怎么破呢，我觉得可以考虑以下两个方法。

方法一、初始种群扩大化

初始种群可以不是8个，可以是16个、32个，或者更多，总之就是让一开始就使得第一代有足够的机会爬到“方框山”的比较好的位置去。

方法二、每一代遴选增加名额

每一代现在我们只是留下8个种群个体，同样的，我们可以留下16个、32个，也一样会让更多爬到“方框山”较高位置的对象存活下来。

这两种方法都采用之后是可以让找到最优解的概率大大增加的。