製造業出身のデータサイエンティストがお送りする記事
今回は多目的最適化手法の中で、NSGA-Ⅱを実装（サンプルコード使用）しました。

##はじめに
先日、ご紹介した多目的最適化（NSGA-Ⅱ）を実装したのでご紹介します。
多目的最適化、NSGA-Ⅱとは何か？と言うことについては上記の記事に書いておりますので、そちらをご参照ください。

##使用するライブラリー（deap）
今回は遺伝的アルゴリズムライブラリdeapを使って実装したいと思います。
他にもライブラリーはあるのですが、今回は開発が盛んに行っているらしいdeapを使用しようと思います。

##NSGA-Ⅱの実装
今回はdeapのチュートリアルにNSGA-Ⅱのサンプルがありましたので、そちらを活用しました。

pythonのコードは下記の通りです。

# 必要なライブラリーのインポート
import array
import random
import json

import numpy as np
import matplotlib.pyplot as plt

from math import sqrt

from deap import algorithms
from deap import base
from deap import benchmarks
from deap.benchmarks.tools import diversity, convergence, hypervolume
from deap import creator
from deap import tools

最初に遺伝的アルゴリズムの設計を行っていきます。

# 適合度を最小化することで最適化されるような適合度クラスの作成
creator.create("FitnessMin", base.Fitness, weights=(-1.0, -1.0))
# 個体クラスIndividualを作成
creator.create("Individual", array.array, typecode='d', fitness=creator.FitnessMin)
# Toolboxの作成
toolbox = base.Toolbox()
# 遺伝子が取り得る値の範囲を指定
BOUND_LOW, BOUND_UP = 0.0, 1.0
# 1つの個体内の遺伝子の数を指定
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)]

# 遺伝子を生成する関数"attr_gene"を登録
toolbox.register("attr_float", uniform, BOUND_LOW, BOUND_UP, NDIM)
# 個体を生成する関数”individual"を登録
toolbox.register("individual", tools.initIterate, creator.Individual, toolbox.attr_float)
# 個体集団を生成する関数"population"を登録
toolbox.register("population", tools.initRepeat, list, toolbox.individual)
# 評価関数"evaluate"を登録
toolbox.register("evaluate", benchmarks.zdt1)
# 交叉を行う関数"mate"を登録
toolbox.register("mate", tools.cxSimulatedBinaryBounded, low=BOUND_LOW, up=BOUND_UP, eta=20.0)
# 変異を行う関数"mutate"を登録
toolbox.register("mutate", tools.mutPolynomialBounded, low=BOUND_LOW, up=BOUND_UP, eta=20.0, indpb=1.0/NDIM)
# 個体選択法"select"を登録
toolbox.register("select", tools.selNSGA2)

deapにおいて個体は、list型を拡張したデータ型を作成することで個体として扱うことができます。
各コードで何を実施しているのかはコメント文で分かるように記載しているため、コード全体が見にくくなっている部分はご容赦願います。

今回は、目的関数が2つであり、どちらも最小化する設定で問題設定を行っております。deapでは最小化「-1」、最大化「1」を使用します。

次に実際の進化計算部分になります。

def main():
    random.seed(1)

    NGEN = 250 # 繰り返し世代数
    MU = 100 # 集団内の個体数
    CXPB = 0.9 # 交叉率
    
    # 世代ループ中のログに何を出力するかの設定
    stats = tools.Statistics(lambda ind: ind.fitness.values)
    stats.register("min", np.min, axis=0)
    stats.register("max", np.max, axis=0)
    
    logbook = tools.Logbook()
    logbook.header = "gen", "evals", "std", "min", "avg", "max"
    
    # 第一世代の生成
    pop = toolbox.population(n=MU)
    pop_init = pop[:]
    invalid_ind = [ind for ind in pop if not ind.fitness.valid]
    fitnesses = toolbox.map(toolbox.evaluate, invalid_ind)
    for ind, fit in zip(invalid_ind, fitnesses):
        ind.fitness.values = fit
    pop = toolbox.select(pop, len(pop))
    
    record = stats.compile(pop)
    logbook.record(gen=0, evals=len(invalid_ind), **record)
    print(logbook.stream)

    # 最適計算の実行
    for gen in range(1, NGEN):
        # 子母集団生成
        offspring = 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
        
        # 適応度を削除した個体について適応度の再評価を行う
        invalid_ind = [ind for ind in offspring if not ind.fitness.valid]
        fitnesses = toolbox.map(toolbox.evaluate, invalid_ind)
        for ind, fit in zip(invalid_ind, fitnesses):
            ind.fitness.values = fit

        # 次世代を選択
        pop = toolbox.select(pop + offspring, MU)
        record = stats.compile(pop)
        logbook.record(gen=gen, evals=len(invalid_ind), **record)
        print(logbook.stream)

    # 最終世代のハイパーボリュームを出力
    print("Final population hypervolume is %f" % hypervolume(pop, [11.0, 11.0]))

    return pop, pop_init, logbook

一般的な遺伝的アルゴリズムでは、交叉率0.9（CXPB=0.9）は高いと感じるかと思います。しかし、NSGA-Ⅱでは親母集団を保存しているため、交叉率は高くて問題ないという認識です。個人的には交叉率は1.0でも良いのではないかと思っております。

あとは、このmain関数を呼び出します。

if __name__ == '__main__':

    pop, pop_init, stats = main()

今回は、初期サンプルと最終世代の結果を可視化してみたいと思います。

fitnesses_init = np.array([list(pop_init[i].fitness.values) for i in range(len(pop_init))])
fitnesses = np.array([list(pop[i].fitness.values) for i in range(len(pop))])
plt.plot(fitnesses_init[:,0], fitnesses_init[:,1], "b.", label="Initial")
plt.plot(fitnesses[:,0], fitnesses[:,1], "r.", label="Optimized" )
plt.legend(loc="upper right")
plt.title("fitnesses")
plt.xlabel("f1")
plt.ylabel("f2")
plt.grid(True)

##さいごに
最後まで読んで頂き、ありがとうございました。
今回は、多目的最適化のNSGA-Ⅱについてサンプルコードを確認しました。

訂正要望がありましたら、ご連絡頂けますと幸いです。