前言

自從上三篇博客詳細講解了Python遺傳和進化算法工具箱及其在帶約束的單目標函數(shù)值優(yōu)化中的應(yīng)用、利用遺傳算法求解有向圖的最短路徑、利用進化算法優(yōu)化SVM參數(shù)之后，這篇不再局限于單一的進化算法工具箱的講解，相反，這次來個橫向?qū)Ρ?#xff0c;比較目前最流行的幾個python進化算法工具箱/框架在求解多目標問題上的表現(xiàn)。

正文

1 多目標優(yōu)化概念

首先大致講一下多目標優(yōu)化：

在生活中的優(yōu)化問題，往往不只有一個優(yōu)化目標，并且往往無法同時滿足所有的目標都最優(yōu)。例如工人的工資與企業(yè)的利潤。那么多目標優(yōu)化里面什么解才算是優(yōu)秀的？我們一般采用“帕累托最優(yōu)”來衡量解是否優(yōu)秀，其定義我這里摘錄百度百科的一段話：

帕累托最優(yōu)（Pareto Optimality），是指資源分配的一種理想狀態(tài)。假定固有的一群人和可分配的資源，從一種分配狀態(tài)到另一種狀態(tài)的變化中，在沒有使任何人境況變壞的前提下，使得至少一個人變得更好。帕累托最優(yōu)狀態(tài)就是不可能再有更多的帕累托改進的余地；換句話說，帕累托改進是達到帕累托最優(yōu)的路徑和方法。 帕累托最優(yōu)是公平與效率的“理想王國”。是由帕累托提出的。

這段話好像讓人看著依舊有點懵逼，下面直接摘錄一段學(xué)術(shù)性的定義：

【摘自：http://www.geatpy.com】

因此，只要找到一組解，其對應(yīng)的待優(yōu)化目標函數(shù)值的點均落在上面的黑色加粗線上，那么就是我們想要的“帕累托最優(yōu)解”了。此外，假如帕累托最優(yōu)解個數(shù)不可數(shù)，那么我們只需找到上面黑色加粗線上的若干個點即可，并且這些點越分散、分布得越均勻，說明算法的效果越好。

2 多目標進化優(yōu)化算法

多目標進化優(yōu)化算法即利用進化算法結(jié)合多目標優(yōu)化策略來求解多目標優(yōu)化問題。經(jīng)典而久經(jīng)不衰的多目標優(yōu)化算法有：NSGA2、NSGA3、MOEA/D等。其中NSGA2和NSGA3是基于支配的MOEA(Multi-objective evolutionary algorithm)，而MOEA/D是基于分解的MOEA。

前兩者（NSGA2、NSGA3）通過非支配排序（后面馬上講到）來篩選出一堆解中的“非支配解”，并且通過種群的不斷進化，使得種群個體對應(yīng)的解對應(yīng)的目標函數(shù)值的點不斷逼近上圖的黑色加粗線。具體算法就不作詳細闡述了，可詳見以下參考文獻，或看下方的代碼實戰(zhàn)部分：

Deb K , Pratap A , Agarwal S , et al. A fast and elitist multiobjective genetic algorithm: NSGA-II[J]. IEEE Transactions on Evolutionary Computation, 2002, 6(2):0-197.

Deb K , Jain H . An Evolutionary Many-Objective Optimization Algorithm Using Reference-Point-Based Nondominated Sorting Approach, Part I: Solving Problems With Box Constraints[J]. IEEE Transactions on Evolutionary Computation, 2014, 18(4):577-601.

后者（MOEA/D）通過線性或非線性的方式將原多目標問題中各個目標進行聚合，即將多個目標聚合成一個目標，然后利用種群進化不斷逼近全局帕累托最優(yōu)解。這里可能有人會有疑問：“為什么MOEA/D是基于分解的MOEA，但過程中需要對各個目標進行聚合？那不就不叫分解了嗎？”答案很簡單：分解是指將多目標優(yōu)化問題分解為一組單目標子問題或多個多目標子問題，利用子問題之間的鄰域關(guān)系，通過協(xié)作來同時優(yōu)化所有的子問題，從而不斷逼近全局帕累托最優(yōu)解；而聚合是指將多個目標聚合成一個目標，因此MOEA/D里面有“分解”和“聚合”兩個步驟，分解是確定鄰域關(guān)系，聚合是用來方便比較解的優(yōu)劣，兩者并不是矛盾的。具體算法就不作詳細闡述了，可詳見以下參考文獻，或看下方的代碼實戰(zhàn)部分：

Qingfu Zhang, Hui Li. MOEA/D: A Multiobjective Evolutionary Algorithm Based on Decomposition[M]. IEEE Press, 2007.

那么存在只利用聚合而沒有分解這一步來進化優(yōu)化的算法嗎？答案是存在的，比如多目標優(yōu)化自適應(yīng)權(quán)重法（awGA），代碼詳見鏈接：

https://github.com/geatpy-dev/geatpy/blob/master/geatpy/templates/moeas/awGA/moea_awGA_templet.py

這個算法就是通過自適應(yīng)地生成一個權(quán)重向量，來將所有的優(yōu)化目標聚合成單一的優(yōu)化目標，然后進行進化優(yōu)化，當(dāng)然這樣效果自然比不上MOEA/D。

有些單目標優(yōu)化學(xué)得比較溜的讀者可能會疑問：”我找一組固定的權(quán)重，把各個優(yōu)化目標加權(quán)聚合成一個目標，再用單目標優(yōu)化的方法進行優(yōu)化不就完事了嗎？“答案非常簡單：如果有理論能證明所找的這組權(quán)重是最合理、最適合當(dāng)前的優(yōu)化模型的，那么用單目標優(yōu)化的方法來解決多目標優(yōu)化問題當(dāng)然是好事；相反，假如沒有依據(jù)地隨便設(shè)置一組權(quán)重，那么肯定不能用這種方法。

3 代碼實戰(zhàn)

下面以一個雙目標優(yōu)化測試函數(shù)ZDT1和一個三目標優(yōu)化測試函數(shù)DTLZ1為例，橫向?qū)Ρ萪eap、pymoo和geatpy三款進化算法代碼包的NSGA2、NSGA3和MOEA/D算法的表現(xiàn)，版本分別為1.3、0.4.0、2.5.0，測試代碼均為三款代碼包官網(wǎng)給出的案例（在代碼組織結(jié)構(gòu)上稍作修改以方便本文顯示）。

【注】：由于在計算ZDT1和DTLZ1時deap默認采用的是循環(huán)來計算種群中每個個體的目標函數(shù)值，而pymoo和geatpy均為利用numpy來矩陣化計算的，為了統(tǒng)一個體評價時間，避免前者帶來的性能下降，這里將deap也改用與后兩者相同的方法。

實驗設(shè)備：cpu i5 9600k，16g ddr4內(nèi)存，windows10，Python 3.7 x64。

相關(guān)庫的安裝：

pip install deap pip install pymoo pip install geatpy

ZDT1的模型為：

其中M為待優(yōu)化的目標個數(shù)，y為決策變量。

3.1 用NSGA2 優(yōu)化 ZDT1

下面用NSGA2算法來優(yōu)化上面的ZDT1，實驗參數(shù)為：【種群個體數(shù)40，進化代數(shù)200，其他相關(guān)參數(shù)均設(shè)為一樣】，代碼1、代碼2、代碼3分別為用deap、pymoo和geatpy的nsga2來優(yōu)化ZDT1：

代碼 1. deap_nsga2.py

import array import random from deap import base from deap import creator from deap import tools import numpy as np import matplotlib.pyplot as plt import timedef ZDT1(Vars, NDIM): # 目標函數(shù)ObjV1 = Vars[:, 0]gx = 1 + 9 * np.sum(Vars[:, 1:], 1) / (NDIM - 1)hx = 1 - np.sqrt(np.abs(ObjV1) / gx) # 取絕對值是為了避免浮點數(shù)精度異常帶來的影響ObjV2 = gx * hxreturn np.array([ObjV1, ObjV2]).Tcreator.create("FitnessMin", base.Fitness, weights=(-1.0, -1.0)) creator.create("Individual", array.array, typecode='d', fitness=creator.FitnessMin) toolbox = base.Toolbox() # Problem definition BOUND_LOW, BOUND_UP = 0.0, 1.0 NDIM = 30 def uniform(low, up, size=None):try:return [random.uniform(a, b) for a, b in zip(low, up)]except TypeError:return [random.uniform(a, b) for a, b in zip([low] * size, [up] * size)]toolbox.register("attr_float", uniform, BOUND_LOW, BOUND_UP, NDIM) toolbox.register("individual", tools.initIterate, creator.Individual, toolbox.attr_float) toolbox.register("population", tools.initRepeat, list, toolbox.individual) toolbox.register("mate", tools.cxSimulatedBinaryBounded, low=BOUND_LOW, up=BOUND_UP, eta=20.0) toolbox.register("mutate", tools.mutPolynomialBounded, low=BOUND_LOW, up=BOUND_UP, eta=20.0, indpb=1.0/NDIM) toolbox.register("select", tools.selNSGA2)def main(seed=None):random.seed(seed)NGEN = 300MU = 40CXPB = 1.0pop = toolbox.population(n=MU)# Evaluate the individuals with an invalid fitnessinvalid_ind = [ind for ind in pop if not ind.fitness.valid]ObjV = ZDT1(np.array(invalid_ind), NDIM)for ind, i in zip(invalid_ind, range(MU)):ind.fitness.values = ObjV[i, :]# This is just to assign the crowding distance to the individuals# no actual selection is donepop = toolbox.select(pop, len(pop))# Begin the generational processfor gen in range(1, NGEN):# Vary the populationoffspring = tools.selTournamentDCD(pop, len(pop))offspring = [toolbox.clone(ind) for ind in offspring] for ind1, ind2 in zip(offspring[::2], offspring[1::2]):if random.random() <= CXPB:toolbox.mate(ind1, ind2)toolbox.mutate(ind1)toolbox.mutate(ind2)del ind1.fitness.values, ind2.fitness.values# Evaluate the individuals with an invalid fitnessinvalid_ind = [ind for ind in offspring if not ind.fitness.valid]ObjV = ZDT1(np.array([ind for ind in offspring if not ind.fitness.valid]), NDIM)for ind, i in zip(invalid_ind, range(MU)):ind.fitness.values = ObjV[i, :]# Select the next generation populationpop = toolbox.select(pop + offspring, MU)return pop if __name__ == "__main__":start = time.time()pop = main()end = time.time()p = np.array([ind.fitness.values for ind in pop])plt.scatter(p[:, 0], p[:, 1], marker="o", s=10)plt.grid(True)plt.show()print('耗時：', end - start, '秒')

上述代碼1利用deap提供的進化算法框架來實現(xiàn)NSGA2算法。deap的主要特點是輕量化和支持擴展。整個deap內(nèi)部的代碼量很少，可以通過”函數(shù)注冊“來擴展模塊，但由此帶來的結(jié)果便是需要自己寫大量的代碼。比如要在deap上寫一個基因表達式編程(GEP)，以geppy這個GEP框架為例，它除了使用到了deap中的”函數(shù)注冊“功能（方便模塊調(diào)用）外，幾乎等于重寫了一遍deap。

代碼2：pymoo_nsga2.py

from pymoo.algorithms.nsga2 import NSGA2 from pymoo.factory import get_problem from pymoo.optimize import minimize import matplotlib.pyplot as plt import time problem = get_problem("ZDT1") algorithm = NSGA2(pop_size=40, elimate_duplicates=False) start = time.time() res = minimize(problem,algorithm,('n_gen', 300),verbose=False) end = time.time() plt.scatter(res.F[:, 0], res.F[:, 1], marker="o", s=10) plt.grid(True) plt.show() print('耗時：', end - start, '秒')

上述代碼2調(diào)用了pymoo內(nèi)置的NSGA2算法模塊，不過該模塊與pymoo框架深度耦合，無法直觀地看到NSGA2的執(zhí)行過程，非pymoo開發(fā)者想要在上面擴展自己設(shè)計的算法會遇到較大的困難。

代碼3：geatpy_nsga2.py

# -*- coding: utf-8 -*- # -*- coding: utf-8 -*- import numpy as np import geatpy as eaclass ZDT1(ea.Problem): # 繼承Problem父類def __init__(self):name = 'ZDT1' # 初始化name（函數(shù)名稱，可以隨意設(shè)置）M = 2 # 初始化M（目標維數(shù)）maxormins = [1] * M # 初始化maxormins（目標最小最大化標記列表，1：最小化該目標；-1：最大化該目標）Dim = 30 # 初始化Dim（決策變量維數(shù)）varTypes = [0] * Dim # 初始化varTypes（決策變量的類型，0：實數(shù)；1：整數(shù)）lb = [0] * Dim # 決策變量下界ub = [1] * Dim # 決策變量上界lbin = [1] * Dim # 決策變量下邊界（0表示不包含該變量的下邊界，1表示包含）ubin = [1] * Dim # 決策變量上邊界（0表示不包含該變量的上邊界，1表示包含）# 調(diào)用父類構(gòu)造方法完成實例化ea.Problem.__init__(self, name, M, maxormins, Dim, varTypes, lb, ub, lbin, ubin)def aimFunc(self, pop): # 目標函數(shù)Vars = pop.Phen # 得到?jīng)Q策變量矩陣ObjV1 = Vars[:, 0]gx = 1 + 9 * np.sum(Vars[:, 1:], 1) / (self.Dim - 1)hx = 1 - np.sqrt(np.abs(ObjV1) / gx) # 取絕對值是為了避免浮點數(shù)精度異常帶來的影響ObjV2 = gx * hxpop.ObjV = np.array([ObjV1, ObjV2]).T # 把結(jié)果賦值給ObjVif __name__ == "__main__":"""================================實例化問題對象============================="""problem = ZDT1() # 生成問題對象"""==================================種群設(shè)置================================"""Encoding = 'RI' # 編碼方式NIND = 40 # 種群規(guī)模Field = ea.crtfld(Encoding, problem.varTypes, problem.ranges, problem.borders) # 創(chuàng)建區(qū)域描述器population = ea.Population(Encoding, Field, NIND) # 實例化種群對象（此時種群還沒被初始化，僅僅是完成種群對象的實例化）"""================================算法參數(shù)設(shè)置==============================="""myAlgorithm = ea.moea_NSGA2_templet(problem, population) # 實例化一個算法模板對象myAlgorithm.MAXGEN = 300 # 最大進化代數(shù)myAlgorithm.logTras = 0 # 設(shè)置每多少代記錄日志，若設(shè)置成0則表示不記錄日志myAlgorithm.verbose = True # 設(shè)置是否打印輸出日志信息myAlgorithm.drawing = 1 # 設(shè)置繪圖方式（0：不繪圖；1：繪制結(jié)果圖；2：繪制目標空間過程動畫；3：繪制決策空間過程動畫）"""==========================調(diào)用算法模板進行種群進化=========================調(diào)用run執(zhí)行算法模板，得到帕累托最優(yōu)解集NDSet以及最后一代種群。NDSet是一個種群類Population的對象。NDSet.ObjV為最優(yōu)解個體的目標函數(shù)值；NDSet.Phen為對應(yīng)的決策變量值。詳見Population.py中關(guān)于種群類的定義。"""[NDSet, population] = myAlgorithm.run() # 執(zhí)行算法模板，得到非支配種群以及最后一代種群NDSet.save() # 把非支配種群的信息保存到文件中"""==================================輸出結(jié)果=============================="""print('用時：%f 秒' % myAlgorithm.passTime)

上面代碼3調(diào)用了Geatpy內(nèi)置的NSGA2算法模板，見鏈接：

https://github.com/geatpy-dev/geatpy/blob/master/geatpy/templates/moeas/nsga2/moea_NSGA2_templet.py

該算法模板代碼與Geatpy提供的進化算法框架的耦合度很小，可以清晰地看到NSGA2算法的執(zhí)行過程。只需遵循Geatpy中設(shè)計的種群數(shù)據(jù)結(jié)構(gòu)（如染色體用什么表示、目標函數(shù)值用什么表示、約束用什么表示），就能很容易在上面擴展自己設(shè)計的新進化算法。

此外由代碼3可以看到使用Geatpy時可以更加專注于待優(yōu)化模型的建立，即在做進化算法的應(yīng)用時，不用管進化算法的細節(jié)，而專注于把待優(yōu)化模型寫成一個自定義問題類。

代碼1、2、3的運行結(jié)果如下圖：

3.2 用NSGA3 優(yōu)化 DTLZ1

下面用NSGA3算法來優(yōu)化上面的DTLZ1，【種群個體數(shù)91，進化代數(shù)500，其他相關(guān)參數(shù)均設(shè)為一樣】，代碼4、代碼5、代碼6分別為用deap、pymoo和geatpy的nsga3來優(yōu)化DTLZ1：

代碼4：deap_nsga3.py

from math import factorial import random import matplotlib.pyplot as plt import numpy as np from deap import algorithms from deap import base from deap import creator from deap import tools import matplotlib.pyplot as plt import mpl_toolkits.mplot3d as Axes3d import time# Problem definition NOBJ = 3 K = 5 NDIM = NOBJ + K - 1 P = 12 H = factorial(NOBJ + P - 1) / (factorial(P) * factorial(NOBJ - 1)) BOUND_LOW, BOUND_UP = 0.0, 1.0def DTLZ1(Vars, NDIM): # 目標函數(shù)XM = Vars[:,(NOBJ-1):]g = 100 * (NDIM - NOBJ + 1 + np.sum(((XM - 0.5)**2 - np.cos(20 * np.pi * (XM - 0.5))), 1, keepdims = True))ones_metrix = np.ones((Vars.shape[0], 1))f = 0.5 * np.hstack([np.fliplr(np.cumprod(Vars[:,:NOBJ-1], 1)), ones_metrix]) * np.hstack([ones_metrix, 1 - Vars[:, range(NOBJ - 2, -1, -1)]]) * (1 + g)return f## # Algorithm parameters MU = int(H) NGEN = 500 CXPB = 1.0 MUTPB = 1.0/NDIM ## # Create uniform reference point ref_points = tools.uniform_reference_points(NOBJ, P) # Create classes creator.create("FitnessMin", base.Fitness, weights=(-1.0,) * NOBJ) creator.create("Individual", list, fitness=creator.FitnessMin) # Toolbox initialization def uniform(low, up, size=None):try:return [random.uniform(a, b) for a, b in zip(low, up)]except TypeError:return [random.uniform(a, b) for a, b in zip([low] * size, [up] * size)] toolbox = base.Toolbox() toolbox.register("attr_float", uniform, BOUND_LOW, BOUND_UP, NDIM) toolbox.register("individual", tools.initIterate, creator.Individual, toolbox.attr_float) toolbox.register("population", tools.initRepeat, list, toolbox.individual) toolbox.register("mate", tools.cxSimulatedBinaryBounded, low=BOUND_LOW, up=BOUND_UP, eta=30.0) toolbox.register("mutate", tools.mutPolynomialBounded, low=BOUND_LOW, up=BOUND_UP, eta=20.0, indpb=1.0/NDIM) toolbox.register("select", tools.selNSGA3, ref_points=ref_points) ## def main(seed=None):random.seed(seed)# Initialize statistics objectpop = toolbox.population(n=MU)# Evaluate the individuals with an invalid fitnessinvalid_ind = [ind for ind in pop if not ind.fitness.valid]ObjV = DTLZ1(np.array(invalid_ind), NDIM)for ind, i in zip(invalid_ind, range(MU)):ind.fitness.values = ObjV[i, :]# Begin the generational processfor gen in range(1, NGEN):offspring = algorithms.varAnd(pop, toolbox, CXPB, MUTPB)# Evaluate the individuals with an invalid fitnessinvalid_ind = [ind for ind in offspring if not ind.fitness.valid]ObjV = DTLZ1(np.array([ind for ind in offspring if not ind.fitness.valid]), NDIM)for ind, i in zip(invalid_ind, range(MU)):ind.fitness.values = ObjV[i, :]# Select the next generation population from parents and offspringpop = toolbox.select(pop + offspring, MU)return popif __name__ == "__main__":start = time.time()pop = main()end = time.time()fig = plt.figure()ax = fig.add_subplot(111, projection="3d")p = np.array([ind.fitness.values for ind in pop])ax.scatter(p[:, 0], p[:, 1], p[:, 2], marker="o", s=10)ax.view_init(elev=30, azim=45)plt.grid(True)plt.show()print('耗時：', end - start, '秒')

代碼5：pymoo_nsga3.py

from pymoo.algorithms.nsga3 import NSGA3 from pymoo.factory import get_problem, get_reference_directions from pymoo.optimize import minimize from pymoo.visualization.scatter import Scatter import time# create the reference directions to be used for the optimization M = 3 ref_dirs = get_reference_directions("das-dennis", M, n_partitions=12) N = ref_dirs.shape[0] # create the algorithm object algorithm = NSGA3(pop_size=N,ref_dirs=ref_dirs) start = time.time() # execute the optimization res = minimize(get_problem("dtlz1", n_obj = M),algorithm,termination=('n_gen', 500)) end = time.time() Scatter().add(res.F).show() print('耗時：', end - start, '秒')

代碼6：geatpy_nsga3.py

# -*- coding: utf-8 -*- import numpy as np import geatpy as eaclass DTLZ1(ea.Problem): # 繼承Problem父類def __init__(self, M):name = 'DTLZ1' # 初始化name（函數(shù)名稱，可以隨意設(shè)置）maxormins = [1] * M # 初始化maxormins（目標最小最大化標記列表，1：最小化該目標；-1：最大化該目標）Dim = M + 4 # 初始化Dim（決策變量維數(shù)）varTypes = np.array([0] * Dim) # 初始化varTypes（決策變量的類型，0：實數(shù)；1：整數(shù)）lb = [0] * Dim # 決策變量下界ub = [1] * Dim # 決策變量上界lbin = [1] * Dim # 決策變量下邊界（0表示不包含該變量的下邊界，1表示包含）ubin = [1] * Dim # 決策變量上邊界（0表示不包含該變量的上邊界，1表示包含）# 調(diào)用父類構(gòu)造方法完成實例化ea.Problem.__init__(self, name, M, maxormins, Dim, varTypes, lb, ub, lbin, ubin)def aimFunc(self, pop): # 目標函數(shù)Vars = pop.Phen # 得到?jīng)Q策變量矩陣XM = Vars[:,(self.M-1):]g = 100 * (self.Dim - self.M + 1 + np.sum(((XM - 0.5)**2 - np.cos(20 * np.pi * (XM - 0.5))), 1, keepdims = True))ones_metrix = np.ones((Vars.shape[0], 1))f = 0.5 * np.hstack([np.fliplr(np.cumprod(Vars[:,:self.M-1], 1)), ones_metrix]) * np.hstack([ones_metrix, 1 - Vars[:, range(self.M - 2, -1, -1)]]) * (1 + g)pop.ObjV = f # 把求得的目標函數(shù)值賦值給種群pop的ObjVif __name__ == "__main__":"""================================實例化問題對象============================="""problem = DTLZ1(3) # 生成問題對象"""==================================種群設(shè)置================================"""Encoding = 'RI' # 編碼方式NIND = 91 # 種群規(guī)模Field = ea.crtfld(Encoding, problem.varTypes, problem.ranges, problem.borders) # 創(chuàng)建區(qū)域描述器population = ea.Population(Encoding, Field, NIND) # 實例化種群對象（此時種群還沒被初始化，僅僅是完成種群對象的實例化）"""================================算法參數(shù)設(shè)置==============================="""myAlgorithm = ea.moea_NSGA3_templet(problem, population) # 實例化一個算法模板對象`myAlgorithm.MAXGEN = 500 # 最大進化代數(shù)myAlgorithm.logTras = 0 # 設(shè)置每多少代記錄日志，若設(shè)置成0則表示不記錄日志myAlgorithm.verbose = True # 設(shè)置是否打印輸出日志信息myAlgorithm.drawing = 1 # 設(shè)置繪圖方式（0：不繪圖；1：繪制結(jié)果圖；2：繪制目標空間過程動畫；3：繪制決策空間過程動畫）"""==========================調(diào)用算法模板進行種群進化=========================調(diào)用run執(zhí)行算法模板，得到帕累托最優(yōu)解集NDSet以及最后一代種群。NDSet是一個種群類Population的對象。NDSet.ObjV為最優(yōu)解個體的目標函數(shù)值；NDSet.Phen為對應(yīng)的決策變量值。詳見Population.py中關(guān)于種群類的定義。"""[NDSet, population] = myAlgorithm.run() # 執(zhí)行算法模板，得到非支配種群以及最后一代種群NDSet.save() # 把非支配種群的信息保存到文件中"""==================================輸出結(jié)果=============================="""print('用時：%f 秒' % myAlgorithm.passTime)

上面代碼6調(diào)用了Geatpy內(nèi)置的NSGA3算法模板，見鏈接：

https://github.com/geatpy-dev/geatpy/blob/master/geatpy/templates/moeas/nsga3/moea_NSGA3_templet.py

代碼4、5、6的運行結(jié)果如下圖：

3.3 用MOEA/D 優(yōu)化 ZDT1

下面用MOEA/D算法來優(yōu)化上面的ZDT1，【種群個體數(shù)40，進化代數(shù)300，其他相關(guān)參數(shù)均設(shè)為一樣】，由于deap并無提供MOEA/D的代碼，因此這里只比較pymoo和geatpy。

值得注意的是：理論上MOEA/D會比NSGA2快，但由于python語言的限制，MOEA/D的算法設(shè)計在python上的執(zhí)行效率會大打折扣。這種情況同樣存在于Matlab中。如果是純C/C++或者Java等，MOEA/D的高效率才能夠真正展現(xiàn)出來。當(dāng)然，假如在C/C++或Java中采用細粒度的并行來執(zhí)行進化，那么MOEA/D的執(zhí)行效率是遠遠比不上NSGA2的，這是因為MOEA/D的算法設(shè)計天然地不支持細粒度的并行。

下面的代碼7、代碼8分別為用pymoo和geatpy的MOEA/D來優(yōu)化ZDT1的代碼：

代碼7：pymoo_moead.py

from pymoo.algorithms.moead import MOEAD from pymoo.factory import get_problem, get_reference_directions from pymoo.optimize import minimize import matplotlib.pyplot as plt import time problem = get_problem("ZDT1") algorithm = MOEAD(get_reference_directions("das-dennis", 2, n_partitions=40),n_neighbors=40,decomposition="tchebi",prob_neighbor_mating=1 ) start = time.time() res = minimize(problem, algorithm, termination=('n_gen', 300)) end = time.time() plt.scatter(res.F[:, 0], res.F[:, 1], marker="o", s=10) plt.grid(True) plt.show() print('耗時：', end - start, '秒')

代碼8：geatpy_moead.py

# -*- coding: utf-8 -*- import numpy as np import geatpy as eaclass ZDT1(ea.Problem): # 繼承Problem父類def __init__(self):name = 'ZDT1' # 初始化name（函數(shù)名稱，可以隨意設(shè)置）M = 2 # 初始化M（目標維數(shù)）maxormins = [1] * M # 初始化maxormins（目標最小最大化標記列表，1：最小化該目標；-1：最大化該目標）Dim = 30 # 初始化Dim（決策變量維數(shù)）varTypes = [0] * Dim # 初始化varTypes（決策變量的類型，0：實數(shù)；1：整數(shù)）lb = [0] * Dim # 決策變量下界ub = [1] * Dim # 決策變量上界lbin = [1] * Dim # 決策變量下邊界（0表示不包含該變量的下邊界，1表示包含）ubin = [1] * Dim # 決策變量上邊界（0表示不包含該變量的上邊界，1表示包含）# 調(diào)用父類構(gòu)造方法完成實例化ea.Problem.__init__(self, name, M, maxormins, Dim, varTypes, lb, ub, lbin, ubin)def aimFunc(self, pop): # 目標函數(shù)Vars = pop.Phen # 得到?jīng)Q策變量矩陣ObjV1 = Vars[:, 0]gx = 1 + 9 * np.sum(Vars[:, 1:], 1) / (self.Dim - 1)hx = 1 - np.sqrt(np.abs(ObjV1) / gx) # 取絕對值是為了避免浮點數(shù)精度異常帶來的影響ObjV2 = gx * hxpop.ObjV = np.array([ObjV1, ObjV2]).T # 把結(jié)果賦值給ObjVif __name__ == "__main__":"""================================實例化問題對象============================="""problem = ZDT1() # 生成問題對象"""==================================種群設(shè)置================================"""Encoding = 'RI' # 編碼方式NIND = 40 # 種群規(guī)模Field = ea.crtfld(Encoding, problem.varTypes, problem.ranges, problem.borders) # 創(chuàng)建區(qū)域描述器population = ea.Population(Encoding, Field, NIND) # 實例化種群對象（此時種群還沒被初始化，僅僅是完成種群對象的實例化）"""================================算法參數(shù)設(shè)置==============================="""myAlgorithm = ea.moea_MOEAD_templet(problem, population) # 實例化一個算法模板對象`myAlgorithm.Ps = 1.0 # (Probability of Selection)表示進化時有多大的概率只從鄰域中選擇個體參與進化myAlgorithm.MAXGEN = 300 # 最大進化代數(shù)myAlgorithm.logTras = 0 # 設(shè)置每多少代記錄日志，若設(shè)置成0則表示不記錄日志myAlgorithm.verbose = True # 設(shè)置是否打印輸出日志信息myAlgorithm.drawing = 1 # 設(shè)置繪圖方式（0：不繪圖；1：繪制結(jié)果圖；2：繪制目標空間過程動畫；3：繪制決策空間過程動畫）"""==========================調(diào)用算法模板進行種群進化=========================調(diào)用run執(zhí)行算法模板，得到帕累托最優(yōu)解集NDSet以及最后一代種群。NDSet是一個種群類Population的對象。NDSet.ObjV為最優(yōu)解個體的目標函數(shù)值；NDSet.Phen為對應(yīng)的決策變量值。詳見Population.py中關(guān)于種群類的定義。"""[NDSet, population] = myAlgorithm.run() # 執(zhí)行算法模板，得到非支配種群以及最后一代種群NDSet.save() # 把非支配種群的信息保存到文件中"""==================================輸出結(jié)果=============================="""print('用時：%f 秒' % myAlgorithm.passTime)

上面代碼8調(diào)用了Geatpy內(nèi)置的MOEA/D算法模板，見鏈接：

https://github.com/geatpy-dev/geatpy/blob/master/geatpy/templates/moeas/moead/moea_MOEAD_templet.py

代碼7、8的運行結(jié)果如下：

這里比較奇怪的是pymoo的基于切比雪夫法的MOEA/D經(jīng)常會出現(xiàn)一個很差的解，如上面左圖的最左邊的點，應(yīng)該是該版本（0.4.0）有BUG。

3.3 更多的比較

下面加大計算規(guī)模，比較deap、pymoo和geatpy在大規(guī)模種群下的表現(xiàn)。

先將代碼1、2、3中種群規(guī)模調(diào)整為500，進化代數(shù)依舊是300，運行結(jié)果如下：

再將代碼2、3、4中種群規(guī)模調(diào)整為496，進化代數(shù)依舊是500，運行結(jié)果如下：

更多比較結(jié)果見下表，這里增添4個經(jīng)典的EA框架/平臺：Pygmo、Pagmo、JMetal和Platemo來共同參與對比，各實驗參數(shù)均保持一致。其中Pygmo的底層是Pagmo，頂層用Python代碼驅(qū)動；Pagmo是純C++的EA框架；JMetal是Java的EA框架；Platemo是Matlab的MOEA實驗與測試平臺。

此外需要注明的是本次實驗中除了Pymoo和Platemo在大規(guī)模計算下會自動使用CPU并行計算外，其余框架的測試均不會自動使用CPU進行并行計算。其中Pagmo可以利用C++的Eigen矩陣庫來加速進化計算，但眾所周知Eigen速度不及Matlab的MKL矩陣庫，且鑒于實驗環(huán)境配置復(fù)雜，這里就沒有使用Pagmo的矩陣計算了。

在上面的實驗中沒有計算多目標優(yōu)化的各項評價指標，但從圖便可大致推斷效果的好壞。其中性能大幅領(lǐng)先于其他的是國產(chǎn)的Geatpy和Platemo，這是因為Geatpy的種群進化過程采用的是Geatpy團隊自主研發(fā)的超高性能矩陣庫進行計算（注：由于Geatpy尚未提供一個統(tǒng)一的標簽來設(shè)置是否采用內(nèi)核并行，上面的實驗為了圖方便均沒使用Geatpy的內(nèi)核并行）；而Platemo利用Matlab自帶的高性能矩陣庫MKL進行計算（在大規(guī)模下會智能地采用CPU并行計算）。

總結(jié)

本文比較了幾種經(jīng)典的進化算法框架/工具箱/平臺在多目標優(yōu)化上的表現(xiàn)。由于都具有權(quán)威性和各自的特色，這里我不作好壞對比，即不評價孰優(yōu)孰劣，也不強行說推薦一定要用哪個，一切以自己習(xí)慣的編程語言為主。

不過如果讀者有將進化計算與機器學(xué)習(xí)、神經(jīng)網(wǎng)絡(luò)、深度學(xué)習(xí)結(jié)合的需求，我在此建議使用Python的進化算法框架。用Python的話編程效率高，并且Python作為一種膠水語言，可以很容易嵌套進各種實際應(yīng)用項目中，同時由上面的實驗也可看出某些Python進化算法框架/工具箱也能達到乃至超越Java、C/C++的進化算法框架的執(zhí)行速度（如Geatpy）。

最后有句話不得不說，尤其是中興事件和華為事件之后，美國技術(shù)變得越來越不可靠了。我個人建議大家在以后的學(xué)習(xí)、工作、研究中逐步減小對matlab的依賴。因為matlab是閉源的，我們能否使用全靠美國總統(tǒng)一張嘴。日后哪天上臺一個更加激進的美國總統(tǒng)，下令全面封禁我們使用matlab，那就麻煩大了。而Python、C/C++等是開源的，不懼怕被封禁。因此更加推薦使用它們來開發(fā)更多國產(chǎn)平臺、框架。

總結(jié)

以上是生活随笔為你收集整理的算法代码_Python进化算法之多目标优化与代码实战的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。