遗传算法解决带约束的单目标优化问题

目录

【第一步】：通过继承 Problem 问题类完成对问题模型的描述。

【第二步】：编写执行脚本“main.py”调用算法模板进行求解。

 代码解析：

---

关于下述所涉及的geatpy库的学习，可以通过 Geatpy库学习流程

待优化的问题模型如下：

        这是一个带不等式约束和等式约束的单目标最大化优化问题，存在多个局部最优解， 对进化算法具有一定的挑战性。全局最优解为：f (0.5, 0, 0.5) = 2.5。这里拟采用差分进 化算法“DE/best/1/L”来求解该问题 (算法模板源码详见“soea_DE_best_1_L_templet.py” )，此时只需要进行编写问题子类和编写执行脚本两个步骤即可完成问题的求解。

【第一步】：通过继承 Problem 问题类完成对问题模型的描述。

        编写“MyProblem.py”文件：

# -*- coding: utf-8 -*-
"""MyProblem.py"""
import numpy as np
import geatpy as ea
class MyProblem(ea.Problem): # 继承Problem父类
    def __init__(self):
        name = 'MyProblem' # 初始化name（函数名称，可以随意设置）
        M = 1 # 初始化M（目标维数）
        maxormins = [-1] # 初始化目标最小最大化标记列表，1：min；-1：max
        Dim = 3 # 初始化Dim（决策变量维数）
        varTypes = [0] * Dim # 初始化决策变量类型，0：连续；1：离散
        lb = [0,0,0] # 决策变量下界
        ub = [1,1,2] # 决策变量上界
        lbin = [1,1,0] # 决策变量下边界
        ubin = [1,1,0] # 决策变量上边界
        # 调用父类构造方法完成实例化
        ea.Problem.__init__(self, name, M, maxormins, Dim, varTypes, lb,
        ub, lbin, ubin)
    def aimFunc(self, pop): # 目标函数，pop为传入的种群对象
        Vars = pop.Phen # 得到决策变量矩阵
        x1 = Vars[:, [0]] # 取出第一列得到所有个体的x1组成的列向量
        x2 = Vars[:, [1]] # 取出第二列得到所有个体的x2组成的列向量
        x3 = Vars[:, [2]] # 取出第三列得到所有个体的x3组成的列向量
        # 计算目标函数值，赋值给pop种群对象的ObjV属性
        pop.ObjV = 4*x1 + 2*x2 + x3
        # 采用可行性法则处理约束，生成种群个体违反约束程度矩阵
        pop.CV = np.hstack([2*x1 + x2 - 1, # 第一个约束
        x1 + 2*x3 - 2, # 第二个约束
        np.abs(x1 + x2 + x3 - 1)]) # 第三个约束

        代码解析：

        语句 class MyProblem(ea.Problem) 意思是定义一个名为 MyProblem 的类，括号中的“ea.Problem”指明了该类是继承了 Problem 类的。

        （特别注意：Python 的面向对象和 Java 等语言的面向对象有个很大的不同之处在于： 子类实例化时并不会自动调用父类的构造方法。这意味着子类单纯编写 __init__() 构造 方法后，仍然无法继承父类所拥有的属性。而父类 Problem 中定义了跟问题有关的一些 重要属性（详见“Problem.py”），因此自定义的问题子类必须在构造方法中显式调用父 类的构造方法。）

        上面的 maxormins 被定义为一个 list 列表，在调用 Problem 父类的构造方法后，问 题类的成员属性 maxormins 将被转化为一个 Numpy array 行向量。

        问题类的“aimFunc()”函数即为待优化的目标函数。上面的目标函数传入了名为 “pop”的参数，它是一个 Population 类的对象，代表着一个种群。Population 类拥有 Chrom、 Phen、ObjV、CV、FitnV 等重要属性，分别指代种群的染色体矩阵、染色体解码后得到的 表现型矩阵、目标函数值矩阵、违反约束程度矩阵、适应度列向量（详见 Population.py 关于种群类的定义）。这些属性的每一行都对应着一个个体。这里的“表现型”矩阵实质上是等价于决策变量矩阵的。因为染色体解码后对应的是决策变量。

        值得注意的是，“决 策变量矩阵”是整个种群所有个体的决策变量组成的矩阵，其每一行对应一个个体，每 一列对应一个决策变量。在计算目标函数后，赋值给种群对象 pop 的 ObjV 属性，至此 完成目标函数值的计算。

        本案例的问题包含不等式约束和等式约束，在 Geatpy 中有两种处理约束条件的方 法：罚函数法和可行性法则。上面的代码展示的即为采用可行性法则处理约束的用法， 它需要计算每个个体违反约束的程度，并把结果保存在种群类的 CV 矩阵中。CV 矩阵 的每一行对应一个个体、每一列对应一个约束条件（可以是等式约束也可以是不等式约 束，关于 CV 的数据结构详见前面的“Geatpy 数据结构”章节），CV 矩阵中元素小于或 等于 0 表示对应个体满足对应的约束条件，否则是违反对应的约束条件，大于 0 的值越大，表示违反约束的程度越高。

        在上面的例子中，第一个约束是 2x1 + x2 ≤ 1，因此可以用 2x1 + x2 − 1 来衡量违 反该约束的程度。由于 x1、x2 和 x3 都是列向量，因此求得的 2x1 + x2 − 1 也是列向 量。第二、三个约束与之类似，最后把三个列向量（对应三个约束条件）用 Numpy 的 “hstack()”函数拼合在一起，得到 CV 矩阵。

        采用可行性法则处理约束只需要关注如何生成 CV 矩阵，不需要关注算法具体是如 何根据 CV 矩阵处理约束条件的。因为不同的算法根据 CV 矩阵处理约束的方式往往会 有所不同，可详见 Geatpy 内置算法模板注释中所标注的算法参考文献。

        若要使用罚函数法，则不需要生成 CV 矩阵，最简单的方法是利用 Numpy 的“where” 语句把违反约束条件的个体索引找到，并根据该索引对种群的对应位置上的目标函数值 加以惩罚即可。因此若要采用这种方法，目标函数“aimFunc()”可如下定义：

        

def aimFunc(self, pop): # 目标函数，pop为传入的种群对象
    Vars = pop.Phen # 得到决策变量矩阵
    x1 = Vars[:, [0]] # 取出第一列得到所有个体的x1组成的列向量
    x2 = Vars[:, [1]] # 取出第二列得到所有个体的x2组成的列向量
    x3 = Vars[:, [2]] # 取出第三列得到所有个体的x3组成的列向量
    f = 4*x1 + 2*x2 + x3 # 计算目标函数值
    # 采用罚函数法处理约束
    exIdx1 = np.where(2*x1+x2>1)[0] # 找到违反约束条件1的个体索引
    exIdx2 = np.where(x1+2*x3>2)[0] # 找到违反约束条件2的个体索引
    exIdx3 = np.where(x1+x2+x3!=1)[0] # 找到违反约束条件3的个体索引
    exIdx = np.unique(np.hstack([exIdx1,exIdx2,exIdx3])) # 合并索引
    alpha = 2 # 惩罚缩放因子
    beta = 1 # 惩罚最小偏移量
    f[exIdx] += self.maxormins[0]*alpha *(np.max(f)-np.min(f)+beta)
    pop.ObjV = f # 把目标函数值矩阵赋值给种群的ObjV属性

        但本案例中包含一个等式约束，用这种简单的惩罚方法难以找到可行解（最坏结果 会是一个可行解都没有找到），读者可以自行尝试。

【第二步】：编写执行脚本“main.py”调用算法模板进行求解。

# -*- coding: utf-8 -*-
"""main.py"""
import numpy as np
import geatpy as ea # import geatpy
from MyProblem import MyProblem # 导入自定义问题接口
"""============================实例化问题对象========================"""
problem = MyProblem() # 实例化问题对象
"""==============================种群设置==========================="""
Encoding = 'RI' # 编码方式
NIND = 50 # 种群规模
Field = ea.crtfld(Encoding, problem.varTypes, problem.ranges,
problem.borders) # 创建区域描述器
population = ea.Population(Encoding, Field, NIND)
# 实例化种群对象（此时种群还没被真正初始化，仅仅是生成一个种群对象）
"""===========================算法参数设置=========================="""
myAlgorithm = ea.soea_DE_best_1_L_templet(problem, population) #
实例化一个算法模板对象
myAlgorithm.MAXGEN = 1000 # 最大遗传代数
myAlgorithm.mutOper.F = 0.5 # 设置差分进化的变异缩放因子
myAlgorithm.recOper.XOVR = 0.5 # 设置交叉概率
myAlgorithm.drawing = 1 # 设置绘图方式
"""=====================调用算法模板进行种群进化====================="""
[bst_trace, population] = myAlgorithm.run() # 执行算法模板
# print(population)
# print(obj_trace)
# # 输出结果
# best_gen = np.argmax(obj_trace[:, 1]) # 记录最优种群是在哪一代
# best_ObjV = obj_trace[best_gen, 1]
# print('最优的目标函数值为：%s'%(best_ObjV))
# print('最优的决策变量值为：')
# for i in range(var_trace.shape[1]):
#     print(var_trace[best_gen, i])
# print('有效进化代数：%s'%(obj_trace.shape[0]))
# print('最优的一代是第 %s 代'%(best_gen + 1))
# print('评价次数：%s'%(myAlgorithm.evalsNum))
# print('时间已过 %s 秒'%(myAlgorithm.passTime))

 代码解析：

        在“main.py”执行脚本中，一开始需要实例化一个问题对象。然后是种群对象的实例化。在实例化种群对象前，需要设定种群的编码方式 Encoding、种群规模 NIND，并且生成区域描述器 Field（或称译码矩阵），因为种群类的构造方法中需要至 少用到这三个参数（详见“Population.py”中种群类的构造方法）。

        值得注意的是，此时种群并未真正初始化（它还没有染色体），它仅仅是被实例化而已。种群的初始化工作 在具体的算法模板中完成。 在完成了问题类对象和种群对象的实例化后，将其传入算法模板类的构造方法来实例化一个算法模板对象。

        本案例中实例化的是“soea_DE_best_1_L_templet”算法模板。 它实现的是“DE/best/1/L”差分进化算法。里面的进化算法具体是如何操作的，可详见 “soea_DE_best_1_L_templet.py”。

        完成算法模板的实例化后，可以在运行算法模板之前对算法模板的一些参数进行设置，比如交叉概率、最大进化代数、绘图方式等等。也可以不进行设置，此时这些参数的值将采用默认的设置（详见算法模板源码以及“Algorithm.py”）。

        在算法模板类中有一个名为“drawing”的参数用于设置绘图的方式。它有三个可以选择的值：

        • drawing = 0 不绘图

        • drawing = 1 绘制结果图

        • drawing = 2 绘制进化过程动态图

        默认情况下 drawing 的值为 1，此时在执行进化算法后，将会对结果进行绘图。

        在完成对进化算法模板的参数设置后，便可以调用其“run()”方法来执行算法模板。 后面便是返回 bst_trace、population。这个 population 为最后一代的种群（可详见算法模板中的返回参数），而bst_trace 记录着种群最优个体的目标函数值和对应的决策变量值。