基础问题求解
该算法包寻找最优解的基本分两步：第一步，搭建问题框架，编写目标函数、约束函数等；第二步；问题求解。简单介绍一个单目标问题求解的样例。
编写问题框架
#定义问题
import numpy as np
import geatpy as ea
"""
max f = x * np.sin(10 * np.pi * x) + 2.0
s.t.
-1 <= x <= 2
"""
class MyProblem(ea.Problem): # 继承Problem父类
def __init__(self):
name = 'MyProblem' # 初始化name（函数名称，可以随意设置）
M = 1 # 初始化M（目标维数）
maxormins = [-1] # 初始化maxormins（目标最小最大化标记列表，1：最小化该目标；-1：最大化该目标）
Dim = 1 # 初始化Dim（决策变量维数）
varTypes = [0] * Dim # 初始化varTypes（决策变量的类型，元素为0表示对应的变量是连续的；1表示是离散的）
lb = [-1] # 决策变量下界
ub = [2] # 决策变量上界
lbin = [1] * Dim # 决策变量下边界（0表示不包含该变量的下边界，1表示包含）
ubin = [1] * Dim # 决策变量上边界（0表示不包含该变量的上边界，1表示包含）
# 调用父类构造方法完成实例化
ea.Problem.__init__(self, name, M, maxormins, Dim, varTypes, lb, ub, lbin, ubin)
def evalVars(self, x): # 目标函数
f = x * np.sin(10 * np.pi * x) + 2.0
return f
问题求解
"""
该案例展示了一个简单的连续型决策变量最大化目标的单目标优化问题的求解。问题的定义详见MyProblem.py。
"""
if __name__ == '__main__':
# 实例化问题对象
problem = MyProblem()
# 构建算法
algorithm = ea.soea_SEGA_templet(problem,
ea.Population(Encoding='RI', NIND=10),#编码规则、种群个数
MAXGEN=25, # 最大进化代数。
logTras=5, # 表示每隔多少代记录一次日志信息，0表示不记录。
trappedValue=1e-6, # 单目标优化陷入停滞的判断阈值。
maxTrappedCount=10) # 进化停滞计数器最大上限值。
# 求解
res = ea.optimize(algorithm, verbose=True, drawing=1, outputMsg=True, drawLog=False, saveFlag=True)
print(res)
其结果很清晰地展示了求解过程，每一代的平均值、最大值、最小值等。以上是一个简单的函数最大值求解案例。
import numpy as np
import geatpy as ea
"""
max f = 4*x1 + 2*x2 + x3
s.t.
2*x1 + x2 - 1 <= 0
x1 + 2*x3 - 2 <= 0
x1 + x2 + x3 - 1 == 0
0 <= x1,x2 <= 1
0 < x3 < 2
"""
class MyProblem(ea.Problem): # 继承Problem父类
def __init__(self):
name = 'MyProblem' # 初始化name（函数名称，可以随意设置）
M = 1 # 初始化M（目标维数）
maxormins = [-1] # 初始化maxormins（目标最小最大化标记列表，1：最小化该目标；-1：最大化该目标）
Dim = 3 # 初始化Dim（决策变量维数）
varTypes = [0] * Dim # 初始化varTypes（决策变量的类型，元素为0表示对应的变量是连续的；1表示是离散的）
lb = [0, 0, 0] # 决策变量下界
ub = [1, 1, 2] # 决策变量上界
lbin = [1, 1, 0] # 决策变量下边界（0表示不包含该变量的下边界，1表示包含）
ubin = [1, 1, 0] # 决策变量上边界（0表示不包含该变量的上边界，1表示包含）
# 调用父类构造方法完成实例化
ea.Problem.__init__(self, name, M, maxormins, Dim, varTypes, lb, ub, lbin, ubin)
def evalVars(self, Vars): # 目标函数
x1 = Vars[:, [0]]
x2 = Vars[:, [1]]
x3 = Vars[:, [2]]
f = 4 * x1 + 2 * x2 + x3
# 采用可行性法则处理约束
CV = np.hstack([2 * x1 + x2 - 1,
x1 + 2 * x3 - 2,
np.abs(x1 + x2 + x3 - 1)])
return f, CV
def calReferObjV(self): # 设定目标数参考值（本问题目标函数参考值设定为理论最优值）
referenceObjV = np.array([[2.5]])
return referenceObjV
"""
该案例展示了一个带等式约束的连续型决策变量最大化目标的单目标优化问题的求解。
"""
if __name__ == '__main__':
# 实例化问题对象
problem = MyProblem()
# 构建算法
algorithm = ea.soea_DE_rand_1_bin_templet(problem,
ea.Population(Encoding='RI', NIND=100),
MAXGEN=500, # 最大进化代数。
logTras=100) # 表示每隔多少代记录一次日志信息，0表示不记录。
algorithm.mutOper.F = 0.5 # 差分进化中的参数F
algorithm.recOper.XOVR = 0.7 # 重组概率
# 求解
res = ea.optimize(algorithm, verbose=True, drawing=1, outputMsg=True, drawLog=False, saveFlag=True)
print(res)
TSP问题求解
Geatpy提供了TSP问题求解的建模方法，同样是先建模，再求解的步骤。
import numpy as np
import geatpy as ea
import matplotlib.pyplot as plt
"""
有十座城市：A, B, C, D, E, F, G, H, I, J，坐标如下：
X Y
[[0.4, 0.4439],
[0.2439,0.1463],
[0.1707,0.2293],
[0.2293,0.761],
[0.5171,0.9414],
[0.8732,0.6536],
[0.6878,0.5219],
[0.8488,0.3609],
[0.6683,0.2536],
[0.6195,0.2634]]
某旅行者从A城市出发，想逛遍所有城市，并且每座城市去且只去一次，最后要返回出发地，
而且需要从G地拿重要文件到D地，另外要从F地把公司的车开到E地，那么他应该如何设计行程方案，才能用
最短的路程来满足他的旅行需求？
分析：在这个案例中，旅行者从A地出发，把其他城市走遍一次后回到A地，因此我们只需要考虑中间途
径的9个城市的访问顺序即可。这9个城市需要排列组合选出满足约束条件的最优的排列顺序作为最终的路线方案。
"""
class MyProblem(ea.Problem): # 继承Problem父类
def __init__(self):
name = 'MyProblem' # 初始化name（函数名称，可以随意设置）
M = 1 # 初始化M（目标维数）
maxormins = [1] # 初始化maxormins（目标最小最大化标记列表，1：最小化该目标；-1：最大化该目标）
Dim = 9 # 初始化Dim（决策变量维数）
varTypes = [1] * Dim # 初始化varTypes（决策变量的类型，元素为0表示对应的变量是连续的；1表示是离散的）
lb = [1] * Dim # 决策变量下界
ub = [9] * Dim # 决策变量上界
lbin = [1] * Dim # 决策变量下边界（0表示不包含该变量的下边界，1表示包含）
ubin = [1] * Dim # 决策变量上边界（0表示不包含该变量的上边界，1表示包含）
# 调用父类构造方法完成实例化
ea.Problem.__init__(self, name, M, maxormins, Dim, varTypes, lb, ub, lbin, ubin)
# 新增一个属性存储旅行地坐标
self.places = np.array([[0.4, 0.4439],
[0.2439, 0.1463],
[0.1707, 0.2293],
[0.2293, 0.761],
[0.5171, 0.9414],
[0.8732, 0.6536],
[0.6878, 0.5219],
[0.8488, 0.3609],
[0.6683, 0.2536],
[0.6195, 0.2634]])
def evalVars(self, x): # 目标函数
# 添加从0地出发且最后回到出发地
X = np.hstack([np.zeros((x.shape[0], 1)), x, np.zeros((x.shape[0], 1))]).astype(int)
ObjV = [] # 存储所有种群个体对应的总路程
for i in range(X.shape[0]):
journey = self.places[X[i], :] # 按既定顺序到达的地点坐标
distance = np.sum(np.sqrt(np.sum(np.diff(journey.T) ** 2, 0))) # 计算总路程
ObjV.append(distance)
f = np.array([ObjV]).T
# 找到违反约束条件的个体在种群中的索引，保存在向量exIdx中（如：若0、2、4号个体违反约束条件，则编程找出他们来）
exIdx1 = np.where(np.where(x == 3)[1] - np.where(x == 6)[1] < 0)[0]
exIdx2 = np.where(np.where(x == 4)[1] - np.where(x == 5)[1] < 0)[0]
exIdx = np.unique(np.hstack([exIdx1, exIdx2]))
CV = np.zeros((x.shape[0], 1))
CV[exIdx] = 1 # 把求得的违反约束程度矩阵赋值给种群pop的CV
return f, CV
"""
该案例展示了一个带约束的单目标旅行商问题的求解。
"""
if __name__ == '__main__':
# 实例化问题对象
problem = MyProblem()
# 构建算法
algorithm = ea.soea_SEGA_templet(problem,
ea.Population(Encoding='P', NIND=50),
MAXGEN=200, # 最大进化代数
logTras=50) # 表示每隔多少代记录一次日志信息，0表示不记录。
algorithm.mutOper.Pm = 0.5 # 变异概率
# 求解
res = ea.optimize(algorithm, verbose=True, drawing=1, outputMsg=True, drawLog=False, saveFlag=True)
# 绘制路线图
if res['success']:
print('最短路程为：%s' % res['ObjV'][0][0])
print('最佳路线为：')
best_journey = np.hstack([0, res['Vars'][0, :], 0])
for i in range(len(best_journey)):
print(int(best_journey[i]), end=' ')
print()
# 绘图
plt.figure()
plt.plot(problem.places[best_journey.astype(int), 0], problem.places[best_journey.astype(int), 1], c='black')
plt.plot(problem.places[best_journey.astype(int), 0], problem.places[best_journey.astype(int), 1], 'o',
c='black')
for i in range(len(best_journey)):
plt.text(problem.places[int(best_journey[i]), 0], problem.places[int(best_journey[i]), 1],
chr(int(best_journey[i]) + 65), fontsize=20)
plt.grid(True)
plt.xlabel('x')
plt.ylabel('y')
plt.savefig('roadmap.svg', dpi=600, bbox_inches='tight')
plt.show()
else:
print('没找到可行解。')
用遗传算法聚类
geatpy也可以解决聚类问题。
"""
该案例展示了如何利用进化算法进行仿k-means聚类（可称之为EA-KMeans算法）。问题的定义详见MyProblem.py。
本案例采用与k-means类似的聚类方法，采用展开的聚类中心点坐标作为染色体的编码，基本流程大致如下：
1) 初始化种群染色体。
2) 迭代进化（循环第3步至第6步），直到满足终止条件。
3) 重组变异，然后根据得到的新染色体计算出对应的聚类中心点。
4) 计算各数据点到聚类中心点的欧式距离。
5) 把与各中心点关联的数据点的坐标平均值作为新的中心点，并以此更新种群的染色体。
6) 把各中心点到与其关联的数据点之间的距离之和作为待优化的目标函数值。
注意：导入的数据是以列为特征的，即每一列代表一个特征（如第一列代表x，第二列代表y......）。
"""
import matplotlib.pyplot as plt
import numpy as np
import geatpy as ea
class MyProblem(ea.Problem): # 继承Problem父类
def __init__(self):
# 目标函数计算中用到的一些数据
self.datas = np.loadtxt('data.csv', delimiter=',') # 读取数据
self.k = 4 # 分类数目
# 问题类设置
name = 'MyProblem' # 初始化name（函数名称，可以随意设置）
M = 1 # 初始化M（目标维数）
maxormins = [1] # 初始化maxormins（目标最小最大化标记列表，1：最小化该目标；-1：最大化该目标）
Dim = self.datas.shape[1] * self.k # 初始化Dim
varTypes = [0] * Dim # 初始化varTypes（决策变量的类型，元素为0表示对应的变量是连续的；1表示是离散的）
lb = list(np.min(self.datas, 0)) * self.k # 决策变量下界
ub = list(np.max(self.datas, 0)) * self.k # 决策变量上界
lbin = [1] * Dim # 决策变量下边界（0表示不包含该变量的下边界，1表示包含）
ubin = [1] * Dim # 决策变量上边界（0表示不包含该变量的上边界，1表示包含）
# 调用父类构造方法完成实例化
ea.Problem.__init__(self, name, M, maxormins, Dim, varTypes, lb, ub, lbin, ubin)
def aimFunc(self, pop): # 目标函数
centers = pop.Phen.reshape(int(pop.sizes * self.k), int(pop.Phen.shape[1] / self.k)) # 得到聚类中心
dis = ea.cdist(centers, self.datas, 'euclidean') # 计算距离
dis_split = dis.reshape(pop.sizes, self.k, self.datas.shape[0]) # 分割距离矩阵，把各个聚类中心到各个点之间的距离的数据分开
labels = np.argmin(dis_split, 1)[0] # 得到聚类标签值
uni_labels = np.unique(labels)
for i in range(len(uni_labels)):
centers[uni_labels[i], :] = np.mean(self.datas[np.where(labels == uni_labels[i])[0], :], 0)
# 直接修改染色体为已知的更优值，加快收敛
pop.Chrom = centers.reshape(pop.sizes, self.k * centers.shape[1])
pop.Phen = pop.decoding() # 染色体解码（要同步修改Phen，否则后面会导致数据不一致）
dis = ea.cdist(centers, self.datas, 'euclidean')
dis_split = dis.reshape(pop.sizes, self.k, self.datas.shape[0])
pop.ObjV = np.sum(np.min(dis_split, 1), 1, keepdims=True) # 计算个体的目标函数值
def draw(self, centers): # 绘制聚类效果图
dis = ea.cdist(centers, self.datas, 'euclidean')
dis_split = dis.reshape(1, self.k, self.datas.shape[0])
labels = np.argmin(dis_split, 1)[0]
colors = ['r', 'g', 'b', 'y']
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
for i in range(self.k):
idx = np.where(labels == i)[0] # 找到同一类的点的下标
datas = self.datas[idx, :]
ax.scatter(datas[:, 0], datas[:, 1], datas[:, 2], c=colors[i])
if __name__ == '__main__':
# 实例化问题对象
problem = MyProblem()
# 构建算法
algorithm = ea.soea_DE_rand_1_bin_templet(problem,
ea.Population(Encoding='RI', NIND=10),
MAXGEN=20, # 最大进化代数。
logTras=5, # 表示每隔多少代记录一次日志信息，0表示不记录。
trappedValue=1e-4, # 单目标优化陷入停滞的判断阈值。
maxTrappedCount=20) # 进化停滞计数器最大上限值。
# 求解
res = ea.optimize(algorithm, verbose=True, drawing=1, outputMsg=True, drawLog=False, saveFlag=True)
# 检验结果
if res['success']:
print('最优的聚类中心为：')
Vars = res['Vars'][0, :]
centers = Vars.reshape(problem.k, int(len(Vars) / problem.k)) # 得到最优的聚类中心
print(centers)
"""=================================检验结果==============================="""
problem.draw(centers)
微信公众号后台回复
加群：加入全球华人OR|AI|DS社区硕博微信学术群
资料：免费获得大量运筹学相关学习资料
人才库：加入运筹精英人才库，获得独家职位推荐
电子书：免费获取平台小编独家创作的优化理论、运筹实践和数据科学电子书，持续更新中ing...
加入我们：加入「运筹OR帷幄」，参与内容创作平台运营
知识星球：加入「运筹OR帷幄」数据算法社区，免费参与每周「领读计划」、「行业inTalk」、「OR会客厅」等直播活动，与数百位签约大V进行在线交流
文章须知
文章作者：用户007
微信编辑：疑疑
文章转载自『Python学习杂记』公众号，原文链接：geatpy遗传算法包介绍（一）
关注我们
FOLLOW US
