差分进化算法(differential evolution)的Python实现
演化算法被归类为一组用于生物进化全局优化的算法，并且基于元启发式搜索方法。可能的解决方案通常跨越问题域上的n维向量空间，并且我们模拟几个总体粒子以达到全局最优。
基本形式中的优化问题包括通过根据为算法提供的程序指令从可能的解元素（矢量）池中选择值来解决最大化或最小化实际函数的任务。进化方法通常遵循具有不同变化的特定策略，以从种群集合中选择候选元素，并应用交叉和/或突变来修改元素，同时尝试提高修饰元素的质量。
这些算法也可以应用于一些有趣的应用，并且在优化NP难题方面也表现出色，包括旅行推销员问题，作业车间调度，图着色，同时也在诸如信号和系统，机械工程和解决数学优化问题。
一种属于演化算法家族的算法是差分进化（DE）算法。在这篇文章中，我们将讨论差分进化算法的一些特性，同时在Python中实现它以优化几个测试函数。
差分进化算法
DE接近一个优化问题，迭代地尝试改进一组给定质量指标（成本函数）的候选解决方案。这些算法集属于元启发式算法，因为它们对被优化的问题做出很少或根本不做任何假设，并且可以搜索可能的解决方案元素的非常大的空间。该算法涉及维持候选解决方案的群体经受重组，评估和选择的迭代。创建新的候选解决方案需要使用称为群体差分权重的参数F对选定元素应用线性运算以生成向量元素，然后基于参数交叉概率CR随机应用交叉。
该算法遵循以下列出的步骤：
实现算法
为了实现上述方法，代码结构如下：
其中，differential_evolution.py我们将执行该算法的主文件。helpers目录由helper类和用于几个操作的函数组成，例如处理与候选元素（point.py）有关的点对象和向量操作，处理所有这些点的集合以及构建population（collection.py），测试要使用的函数的方法/成本函数来测试算法的效率（test_functions.py）。
构建点类的Python实现
# helpers/point.py
import numpy as np
import scipy as sp
class Point:
def __init__(self, dim=2, upper_limit=10, lower_limit=-10, objective=None):
self.dim = dim
self.coords = np.zeros((self.dim,))
self.z = None
self.range_upper_limit = upper_limit
self.range_lower_limit = lower_limit
self.objective = objective
self.evaluate_point()
def generate_random_point(self):
self.coords = np.random.uniform(self.range_lower_limit, self.range_upper_limit, (self.dim,))
self.evaluate_point()
def evaluate_point(self):
# self.z = evaluate(self.coords)
self.z = self.objective.evaluate(self.coords)
在这里，我们初始化Point类，dim它是矢量的尺寸大小，lower_limit并upper_limit指定矢量的每个坐标的域。self.z是该点的目标函数值，与每个实例相关联，以便根据其目标函数值对其进行排序。该evaluate_point函数针对测试函数上的给定点运行目标函数。该Point创建矢量对象的实例，以表示种群中的每个个体。Population类中定义了个体的集合。
Population 类的Python实现
# helpers/population.py
import copy
import numpy as np
from matplotlib import pyplot as plt
from point import Point
from matplotlib import pyplot as plt
class Population:
def __init__(self, dim=2, num_points=50, upper_limit=10, lower_limit=-10, init_generate=True, objective=None):
self.points = []
self.num_points = num_points
self.init_generate = init_generate
self.dim = dim
self.range_upper_limit = upper_limit
self.range_lower_limit = lower_limit
self.objective = objective
# If initial generation parameter is true, then generate collection
if self.init_generate == True:
for ix in xrange(num_points):
new_point = Point(dim=dim, upper_limit=self.range_upper_limit,
lower_limit=self.range_lower_limit, objective=self.objective)
new_point.generate_random_point()
self.points.append(new_point)
def get_average_objective(self):
avg = 0.0
for px in self.points:
avg += px.z
avg = avg/float(self.num_points)
return avg
本Population类包含一组作用在种群中的个体点类的实例。个体存储在self.points列表中。该类的参数num_points包含关于种群大小的信息dim，upper_limit 并且lower_limit如上所述。作为一个可选参数，init_generate控制初始种群的生成，并objective引用Function该类的一个对象，并且是目标函数。如果设置为False，则初始总体将为空，并且需要通过算法的主程序添加元素。该get_average_objective函数返回种群的平均评估目标值。
目标函数的Python实现
# helpers/test_functions.py
import numpy as np
class Function:
def __init__(self, func=None):
self.objectives = {
'sphere': self.sphere,
'ackley': self.ackley,
'rosenbrock': self.rosenbrock,
'rastrigin': self.rastrigin,
}
if func is None:
self.func_name = 'sphere'
self.func = self.objectives[self.func_name]
else:
if type(func) == str:
self.func_name = func
self.func = self.objectives[self.func_name]
else:
self.func = func
self.func_name = func.func_name
def evaluate(self, point):
return self.func(point)
def sphere(self, x):
d = x.shape[0]
f = 0.0
for dx in xrange(d):
f += x[dx] ** 2
return f
def ackley(self, x):
z1, z2 = 0, 0
for i in xrange(len(x)):
z1 += x[i] ** 2
z2 += np.cos(2.0 * np.pi * x[i])
return (-20.0 * np.exp(-0.2 * np.sqrt(z1 / len(x)))) - np.exp(z2 / len(x)) + np.e + 20.0
def rosenbrock(self, x):
v = 0
for i in xrange(len(x) - 1):
v += 100 * (x[i + 1] - x[i] ** 2) ** 2 + (x[i] - 1) ** 2
return v
def rastrigin(self, x):
v = 0
for i in range(len(x)):
v += (x[i] ** 2) - (10 * np.cos(2 * np.pi * x[i]))
return (10 * len(x)) + v
在test_functions.py包含执行Function类，它创建了一个目标函数对象。构造函数的参数func可以是字符串，也可以是函数。如果None，它会存储self.func中的sphere函数，否则应检查字符串值。对于一个字符串，它将为该类指定具有相同名称的函数（存储在字典下self.objectives）。对于一个函数，这个假定函数接受一个numpy的ndarray作为输入并返回一个标量作为目标函数值。
在默认情况下实现的目标函数目前包括sphere，ackley，rosenbrock，和rastrigin函数。可以在这里找到optomization测试功能列表。这些都是在多维矢量空间中定义的，并且具有单模式或多模式属性。例如，球函数是单模凸函数，而rastrigin函数是多模非凸函数。显示rastrigin了三维空间中函数的表示（垂直轴是目标函数的值）：
3-D中的Rastrigin函数
差分进化算法类的Python实现
# differential_evolution.py
import copy
import random
import time
from helpers.population import Population
from helpers import get_best_point
from helpers.test_functions import Function
class DifferentialEvolution(object):
def __init__(self, num_iterations=10, CR=0.4, F=0.48, dim=2, population_size=10, print_status=False, func=None):
random.seed()
self.print_status = print_status
self.num_iterations = num_iterations
self.iteration = 0
self.CR = CR
self.F = F
self.population_size = population_size
self.func = Function(func=func)
self.population = Population(dim=dim, num_points=self.population_size, objective=self.func)
def iterate(self):
for ix in xrange(self.population.num_points):
x = self.population.points[ix]
[a, b, c] = random.sample(self.population.points, 3)
while x == a or x == b or x == c:
[a, b, c] = random.sample(self.population.points, 3)
R = random.random() * x.dim
y = copy.deepcopy(x)
for iy in xrange(x.dim):
ri = random.random()
if ri < self.CR or iy == R:
y.coords[iy] = a.coords[iy] + self.F * (b.coords[iy] - c.coords[iy])
y.evaluate_point()
if y.z < x.z:
self.population.points[ix] = y
self.iteration += 1
def simulate(self):
pnt = get_best_point(self.population.points)
print("Initial best value: " + str(pnt.z))
while self.iteration < self.num_iterations:
if self.print_status == True and self.iteration%50 == 0:
pnt = get_best_point(self.population.points)
print pnt.z, self.population.get_average_objective()
self.iterate()
pnt = get_best_point(self.population.points)
print("Final best value: " + str(pnt.z))
return pnt.z
在这里，在这个DifferentialEvolution类中，初始化参数是：
基本上有两个成员函数，self.iterate并且self.simulate。该self.iterate函数运行Differential Evolution过程的一次迭代，通过对群体中的每个个体应用变换操作和交叉，并且self.simulate函数调用迭代函数，直到满足停止条件，然后打印目标函数的最佳值。
演示
现在我们已经为差分进化算法的所有必需类实现了一个实现，我们可以编写一个小脚本来测试所有内容并查看结果。
# demo.py
from differential_evolution import DifferentialEvolution
import datetime
import numpy as np
from matplotlib import pyplot as plt
if __name__ == '__main__':
number_of_runs = 5
val = 0
print_time = True
for i in xrange(number_of_runs):
start = datetime.datetime.now()
de = DifferentialEvolution(num_iterations=200, dim=10, CR=0.4, F=0.48, population_size=75, print_status=False, func='sphere')
val += de.simulate()
if print_time:
print "Time taken:", datetime.datetime.now() - start
print '-'*80
print "Final average of all runs:", val / number_of_runs
这个脚本初始化变量number_of_runs，val和print_time。number_of_runs用于启动算法的几次运行，最后，优化后的目标函数的平均结果在这些运行之后返回。val 存储每次运行的优化目标函数值，并稍后用于计算平均值。print_time是一个布尔值，用于控制是否应为每次运行打印计算时间。
上述代码的输出（即使用差分进化算法优化球体测试函数）在50维（50维矢量空间）上运行，每次运行200次迭代后会产生以下输出：
# Output
Initial best value: 1285.50913073
Final best value: 0.0258755727525
Time taken: 0:00:05.931056
Initial best value: 1218.54112743
Final best value: 0.0323126608382
Time taken: 0:00:05.560921
Initial best value: 1253.1145944
Final best value: 0.0340955810298
Time taken: 0:00:06.081233
Initial best value: 1298.5615981
Final best value: 0.0439433666035
Time taken: 0:00:04.511034
Initial best value: 1228.13894559
Final best value: 0.0405344973595
Time taken: 0:00:05.081286
--------------------------------------------------------------------
Final average of all runs: 0.0353523357167
目标函数值与50D中球体测试函数的迭代次数和50D中的Rastrigin测试函数的关系图如下:
