数学建模系列——GA原理与实战

——————————————————————————

    这次来写一下遗传算法。遗传算法是一种仿生算法，最主要的思想是物竞天择，适者生存。这个算法很好的模拟了生物的进化过程，保留好的物种，同样一个物种中的佼佼者才会幸运的存活下来。转换到数学问题中，这个思想就可以很好的转化为优化问题，在求解方程组的时候，好的解视为好的物种被保留，坏的解视为坏的物种而淘汰，设置好进化次数以后开始迭代，记录下这些解里面最好的那个，就是要求解的方程组的解。来介绍一下遗传算法的流程。

1. 1. 编码

    对于任意一个方程， y=kx+b y = k x + b ，随机产生它的解，对应为设置好的种群里面的 x x 个生物，对其进行编码，编为二进制代码。如第一个个体是方程组的一个解 x=1 x = 1 ，将其编码为 010101001010 010101001010 这样的形式，同理，第二个物种 x=4 x = 4 编码为 111000101010 111000101010 这样的形式，编码长度相同，但每个个体之间的编码不同。

2. 2. 个体适应度评价

    怎么评价那个物种是好物种，那个物种是坏物种呢？得找个评价方式。以方程 y=x3+sin(46πx+x2)−arctan(48.6x2−sinx) y = x 3 + sin ⁡ ( 46 π x + x 2 ) − arctan ⁡ ( 48.6 x 2 − sin ⁡ x ) 这样的奇葩方程为例，将每一个解带入原方程求得对应的 yi y i ， yi y i 记为适应度。计算适应度的总和 F F ， F=y1+y2+y3+…+yn F = y 1 + y 2 + y 3 + … + y n 。

    下一步计算每个个体有几率生存下去的概率，见下面的公式，适应度越强，生存下去的几率就越大。

Pi=yiF P i = y i F

3. 3. 新种群复制

    产生和个体一样的十个在 0 0 , 1 1 之间的随机数，比如像下面那样的，假设有十个个体。然后用这十个随机数去选择遗传下去的个体，并不是选择适应性最强的个提取遗传，为什么呢，我感觉他是考虑了天灾人祸等因素，并不是所有的强者都能活下去，也有可能有幸运的弱者会继续活下去。但是适应度还是占了一定比例的。对了，还得求一个遗传累积概率，第一个个体遗传的概率为 y1 y 1 ，第二个个体遗传的概率为 y1+y2 y 1 + y 2 ，第 3 3 个个体的遗传概率为 y1+y2+y3 y 1 + y 2 + y 3 ，以此类推。用随机数的第一个数去比较，比如，假设第2个随机数 0.79 0.79 大于前7个个体的累积和，那么第8个个体被保留(个体适应度同样的占一定比例)，最后一个随机数 0.41 0.41 大于前四个个体的累积和，第五个个体被保留，以此类推，选择出被保留的个体去交配。

0.09978328867338249,0.7927769683560828,0.03785578076377627,0.7989147320211916,
0.7511181280134576,0.050985575043039244,0.24741490093940222,0.7897723998715,
0.6745958293779668,0.4123272416111048,

3. 3. 新种群交配

    产生和个体数目一样的在0,1之间随机数，设大于那个随机数大于0.7的个体才有机会去交配，剩下的个体守寡或者打光棍。

    那么怎么交配呢，前面提到每一个解都被编码为 1100 1100 , 0011 0011 这样的二进制代码形式，好，现在让他交配，把这两个个体视为染色体，高中的那个什么染色体互换就派上用场了，两人各取一半相结合，所以上面两个个体的交配结果就是 1111 1111 , 0000 0000 ，当然也可以不取一半，交配成 1101 1101 , 0010 0010 也是可以的。

    交配完了在模仿一下大自然，基因突变。有时候基因突变是好的，有时候是坏的，但是发生基因突变的概率很低。怎么突变呢， 11110000 11110000 以极低的概率突变成 11110001 11110001 这就是突变。

4. 4. 反复迭代

    不关你是生下来的个体，还是基因突变的个体，都要去面临下一轮的选择和淘汰，重复第一步，再次选择，评价，复制，交配，产生下一代，循环往复，记录每次交配后产生的最优解，多次迭代，保留最优解，也就是你函数的最优解了。

5. 5. 程序实例

5.1 5.1 主函数

from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
import math
import utils 
from selection import selection
from crossover import crossover
from mutation import mutation

#max the function
def aimFunction(x):
    y = x ** 3 - 60*x ** 2 + 900*x + 100 
    return y

x = [i/100 for i in range(900)]  # 0开始，0.01的间距，900个数
y = [0 for i in range(900)]      #0开始，900个0
for i in range(900):
    y[i]=aimFunction(x[i])
#plt.plot(x,y)

#initiate population,产生01编码的种群
population=[]
for i in range(10):
    entity=''
    for j in range(17):  
        entity = entity + str(np.random.randint(0, 2))
    population.append(entity)

maxvalue = 0
for i in range(1, 50):
    t=[]
    #迭代次数
    for i in range(1000):
        #selection
        value = utils.fitness(population, aimFunction)  #some solution
        #轮盘赌选择
        population_new = selection(population, value)
        #crossover，0.7是交配率
        offspring = crossover(population_new, 0.7)
        #mutation 0.02是变异概率
        population = mutation(offspring, 0.01)
        result = []
        for j in range(len(population)):
            result.append(aimFunction(utils.decode(population[j])))
        t.append(max(result))
    if maxvalue <= max(t):
        maxvalue = max(t) 

print("result = " + str(maxvalue))

5.2 5.2 交叉函数(交配)

from __future__ import division
import numpy as np

# 交叉
def crossover(population_new, pc):

    half = int(len(population_new)/2)
    father = population_new[:half]
    mother = population_new[half:]
    # 打乱顺序
    #np.random.shuffle(father)
    #np.random.shuffle(mother)
    offspring = []
    for i in range(half):      
        if np.random.uniform(0,1) <= pc:
            copint = np.random.randint(0,int(len(father[i])/2))
            son = father[i][:copint] + (mother[i][copint:])
            daughter = mother[i][:copint] + (father[i][copint:])
        else:
            son = father[i]
            daughter = mother[i]
        offspring.append(son)
        offspring.append(daughter)
    return offspring

5.3 5.3 变异函数

from __future__ import division
import numpy as np

def mutation(offspring, pm):
    for i in range(len(offspring)):
        if np.random.uniform(0, 1) <= pm:
            position = np.random.randint(0, len(offspring[i]))
            #0,1替换产生变异
            if position != 0:
                if offspring[i][position]=='1':
                    offspring[i] = offspring[i][:position-1]+ '0' + offspring[i][position:]
                else:
                    offspring[i] = offspring[i][:position-1]+'1'+offspring[i][position:]
            else:
                if offspring[i][position] == '1':
                    offspring[i] = '0' + offspring[i][1:]
                else:
                    offspring[i] = '1' + offspring[i][1:]
    return offspring

5.4 5.4 选择函数

from __future__ import division
import numpy as np
def selection(population,value):

    #轮盘赌选择
    #累积求和
    fitness_sum=[]
    for i in range(len(value)):
        if i == 0:
            fitness_sum.append(value[i])
        else:
            fitness_sum.append(fitness_sum[i-1]+value[i])

    #个体适应度
    for i in range(len(fitness_sum)):
        fitness_sum[i] /= sum(value)

    #select new population
    population_new=[]
    for i in range(len(value)):
        rand=np.random.uniform(0,1)
        for j in range(len(value)):
            if j == 0:
                if 0 < rand and rand <= fitness_sum[j]:
                    population_new.append(population[j]) 
            else:
                if fitness_sum[j-1] < rand and rand <= fitness_sum[j]:
                    population_new.append(population[j])             
    return population_new

5.5 5.5 适应度函数

from __future__ import division

def decode(x):
    y=0+int(x,2)/(2**17-1)*9
    return y

def fitness(population,aimFunction):
    value=[]
    for i in range(len(population)):
        value.append(aimFunction(decode(population[i])))
        # weed out negative value
        if value[i] < 0:
            value[i] = 0
    return value

6. 6. 求解函数

y=x+5sin(5x)+arctan(7x) y = x + 5 sin ⁡ ( 5 x ) + arctan ⁡ ( 7 x )

result = 14.393654613558743