优化算法 | 遗传算法（附Python代码）

今天为各位讲解遗传算法（GA），相信各位对GA都已非常熟悉，之前我们也出过很多关于GA的教程：

遗传算法（GA）求解带时间窗的车辆路径（VRPTW）问题MATLAB代码
遗传算法（GA）求解容量受限的车辆路径（CVRP）问题MATLAB代码
遗传算法求解0-1背包问题（附matlab源代码）
遗传算法（GA）求解旅行商问题（TSP）附MATLAB代码
优化算法 | 遗传算法的变形
机器学习 | 基于遗传算法的BP神经网络优化算法（附MATLAB代码）
多种群遗传算法的函数优化算法（附MATLAB代码）
基于量子遗传算法的函数寻优算法
遗传算法求解车间调度问题（附MATLAB代码）

之前的教程基本上都是用GA求解离散优化问题，鲜有涉及连续优化问题，因此今天我们主要讲解使用GA求解连续优化问题。

我们首先给出求解问题，目标是求解下述函数最小值：

$f(x)=\sum_{i=1}^{5} x_{i}^{2} \quad-100 \leq x_{i} \leq 100$

▎GA算法求解步骤

01 | 种群初始化

因为求解问题的变量数目为5，所以染色体长度为5。假设种群数目为npop，则随机生成npop个长度为5的个体，并且每个个体在每个维度上的取值应该在-100~100之间。种群初始化的代码如下：

# Initialize Population
pop = empty_individual.repeat(npop)
for i in range(npop):pop[i].position = np.random.uniform(varmin, varmax, nvar)pop[i].cost = costfunc(pop[i].position)if pop[i].cost < bestsol.cost:bestsol = pop[i].deepcopy()

假设npop=5，则初始化种群结果如下，其中position是个体位置，cost是个体目标函数值：

[struct({'position': array([-31.01033353,94.77288865,-91.01647726,1.32204408,-49.62896178]), 'cost': 20692.321989577737}),
struct({'position': array([-69.29141662,63.88025843,69.29075941,34.80789342,81.93170175]), 'cost': 21607.590370704376}),
struct({'position': array([-64.22127924,0.08618959,-56.55556517,-10.68613122,-61.14429882]), 'cost': 11175.730766415392}),
struct({'position': array([ 80.3919982,-56.0863127,-27.38791656,-39.56431118,-31.46714348]), 'cost': 12914.161657029803}),
struct({'position': array([ 27.41598459,24.68736174,-31.05651977,-48.98385997,-43.34750546]), 'cost': 6604.034227715052})]

02 | 适应度值计算

这份代码中的适应度值就是目标函数值，具体代码如下：

# Sphere Test Function
def sphere(x):return sum(x**2)

03 | 选择操作

选择操作就是从当前种群中选择出若干个个体，这份代码使用的是轮盘赌选择策略，也就是适应度值大的个体有更大的概率被选中，具体代码如下：

def roulette_wheel_selection(p):c = np.cumsum(p)r = sum(p)*np.random.rand()ind = np.argwhere(r <= c)return ind[0][0]

假设输入参数p=[1,2,3,4,5]，接下来依次看每一行代码代表什么含义。

cumsum函数是累加函数，则c=[1,3,6,10,15]。

sum函数是求和函数，sum(p)=15，np.random.rand()函数返回一个[0,1)之间的随机数，不包括1，则r的取值范围为[0,15)之间的一个数，不包括15，假设此时r=5。

np.argwhere(a)函数返回非0元素的索引，其中a是要索引数组的条件，并且该函数输出的是一列元素。因此np.argwhere(r <= c)返回的是r <= c的索引，输出结果为

def roulette_wheel_selection(p):c = np.cumsum(p)r = sum(p)*np.random.rand()ind = np.argwhere(r <= c)return ind[0][0]

最终返回ind[0][0]=2，说明r=5落在3~6之间，即落在p=3所属的区间，轮盘赌示意图如下图所示。

04 | 交叉操作

交叉操作就是对选择出的两个个体进行处理，使这两个个体能够互相继承彼此良好的基因片段，具体代码如下：

def crossover(p1, p2, gamma=0.1):c1 = p1.deepcopy()c2 = p1.deepcopy()alpha = np.random.uniform(-gamma, 1+gamma, *c1.position.shape)c1.position = alpha*p1.position + (1-alpha)*p2.positionc2.position = alpha*p2.position + (1-alpha)*p1.positionreturn c1, c2

np.random.uniform(low,high,size)函数返回size相同维度且取值范围在[low,high)之间的随机数，如

np.random.uniform(low = 0, high = 1, size = 4)

返回4个[0,1)之间的随机数，具体返回结果如下

[0.25773193 0.90374664 0.51513421 0.60574244]

下面这两行代码实际上表示产生两个新个体的计算公式。

c1.position = alpha*p1.position + (1-alpha)*p2.position
c2.position = alpha*p2.position + (1-alpha)*p1.position

05 | 变异操作

变异操作是对某一个个体进行处理，更新该个体的位置，具体代码如下，其中mu是变异概率：

def mutate(x, mu, sigma):y = x.deepcopy()flag = np.random.rand(*x.position.shape) <= muind = np.argwhere(flag)y.position[ind] += sigma*np.random.randn(*ind.shape)return y

flag = np.random.rand(*x.position.shape) <= mu这行代码的含义是判断随机生成的与个体长度相等的数字中是否小于等于变异概率mu，返回结果示例为[False True False False False]

ind = np.argwhere(flag)表示找出flag为True的索引，即找出待变异的基因位，返回结果为ind=[[1]]。

y.position[ind] += sigma*np.random.randn(*ind.shape)表示更新待变异的基因位。

▎GA算法流程图

▎GA算法实例验证

实例验证函数如下，目标是求解下述函数最小值：

$f(x)=\sum_{i=1}^{5} x_{i}^{2} \quad-100 \leq x_{i} \leq 100$

ga.py代码如下：

import numpy as np
from ypstruct import structuredef run(problem, params):# Problem Informationcostfunc = problem.costfuncnvar = problem.nvarvarmin = problem.varminvarmax = problem.varmax# Parametersmaxit = params.maxitnpop = params.npopbeta = params.betapc = params.pcnc = int(np.round(pc*npop/2)*2)gamma = params.gammamu = params.musigma = params.sigma# Empty Individual Templateempty_individual = structure()empty_individual.position = Noneempty_individual.cost = None# Best Solution Ever Foundbestsol = empty_individual.deepcopy()bestsol.cost = np.inf# Initialize Populationpop = empty_individual.repeat(npop)for i in range(npop):pop[i].position = np.random.uniform(varmin, varmax, nvar)pop[i].cost = costfunc(pop[i].position)if pop[i].cost < bestsol.cost:bestsol = pop[i].deepcopy()# Best Cost of Iterationsbestcost = np.empty(maxit)# Main Loopfor it in range(maxit):costs = np.array([x.cost for x in pop])avg_cost = np.mean(costs)if avg_cost != 0:costs = costs/avg_costprobs = np.exp(-beta*costs)popc = []for _ in range(nc//2):# Select Parents#q = np.random.permutation(npop)#p1 = pop[q[0]]#p2 = pop[q[1]]# Perform Roulette Wheel Selectionp1 = pop[roulette_wheel_selection(probs)]p2 = pop[roulette_wheel_selection(probs)]# Perform Crossoverc1, c2 = crossover(p1, p2, gamma)# Perform Mutationc1 = mutate(c1, mu, sigma)c2 = mutate(c2, mu, sigma)# Apply Boundsapply_bound(c1, varmin, varmax)apply_bound(c2, varmin, varmax)# Evaluate First Offspringc1.cost = costfunc(c1.position)if c1.cost < bestsol.cost:bestsol = c1.deepcopy()# Evaluate Second Offspringc2.cost = costfunc(c2.position)if c2.cost < bestsol.cost:bestsol = c2.deepcopy()# Add Offsprings to popcpopc.append(c1)popc.append(c2)# Merge, Sort and Selectpop += popcpop = sorted(pop, key=lambda x: x.cost)     #按照目标函数值从小到大排序pop = pop[0:npop]# Store Best Costbestcost[it] = bestsol.cost# Show Iteration Informationprint("Iteration {}: Best Cost = {}".format(it, bestcost[it]))# Outputout = structure()out.pop = popout.bestsol = bestsolout.bestcost = bestcostreturn outdef crossover(p1, p2, gamma=0.1):c1 = p1.deepcopy()c2 = p1.deepcopy()alpha = np.random.uniform(-gamma, 1+gamma, *c1.position.shape)c1.position = alpha*p1.position + (1-alpha)*p2.positionc2.position = alpha*p2.position + (1-alpha)*p1.positionreturn c1, c2def mutate(x, mu, sigma):y = x.deepcopy()flag = np.random.rand(*x.position.shape) <= muind = np.argwhere(flag)y.position[ind] += sigma*np.random.randn(*ind.shape)return ydef apply_bound(x, varmin, varmax):x.position = np.maximum(x.position, varmin)x.position = np.minimum(x.position, varmax)def roulette_wheel_selection(p):c = np.cumsum(p)r = sum(p)*np.random.rand()ind = np.argwhere(r <= c)return ind[0][0]
function [f] = Sphere(x)
f= sum(x.^2);
end

app.py代码如下：

import numpy as np
import matplotlib.pyplot as plt
from ypstruct import structure
import ga# Sphere Test Function
def sphere(x):return sum(x**2)# Problem Definition
problem = structure()
problem.costfunc = sphere
problem.nvar = 5
problem.varmin = [-100, -100, -100, -100,  -100]
problem.varmax = [ 100,  100,  100,  100, 100]# GA Parameters
params = structure()
params.maxit = 500
params.npop = 5
params.beta = 1
params.pc = 1
params.gamma = 0.1
params.mu = 0.01
params.sigma = 0.1# Run GA
out = ga.run(problem, params)# Results
plt.plot(out.bestcost)
print(out.bestsol)
# plt.semilogy(out.bestcost)
plt.xlim(0, params.maxit)
plt.xlabel('Iterations')
plt.ylabel('Best Cost')
plt.title('Genetic Algorithm (GA)')
plt.grid(True)
plt.show()

求解结果如下：

最优解：struct({'position': array([-1.32382941e-04, -2.07103364e-05, -1.31247597e-03,  1.21645087e-04,2.23626596e-08]), 'cost': 1.7553448711624725e-06})

▎参考

[1]Mostapha Kalami Heris, Practical Genetic Algorithms in Python and MATLAB – Video Tutorial (URL:

https://yarpiz.com/632/ypga191215-practical-genetic-algorithms-in-python-and-matlab), Yarpiz, 2020.

咱们下期再见

▎近期你可能错过了的好文章：

新书上架 | 《MATLAB智能优化算法：从写代码到算法思想》

优化算法 | 灰狼优化算法（文末有福利）

优化算法 | 鲸鱼优化算法

遗传算法（GA）求解带时间窗的车辆路径（VRPTW）问题MATLAB代码

粒子群优化算法（PSO）求解带时间窗的车辆路径问题（VRPTW）MATLAB代码