遺伝的アルゴリズムで関数の最適値を求める（その１）

進化アルゴリズムとは

進化的アルゴリズム（evolutionary algorithm, EA）とは、最適化アルゴリズムの一つで、自然界の淘汰、遺伝子組み換え、突然変異などをモデルとして組み込んだ最適化手法である。

進化的アルゴリズムのは下記の4種類に大分される。
- 遺伝的アルゴリズム
- 遺伝的プログラミング
- 進化戦略
- 進化的プログラミング

今回は遺伝的アルゴリズム（genetic algorithm, GA）を扱う。

使用したライブラリ

下記のDEAPというpythonライブラリを用いる。
https://deap.readthedocs.io/en/master/

コード

目的関数

下記に示すBird Bunctionの最小値を探索する。
$$f(x, y) = sin(x)e^{(1-cos(y))^2}+cos(y)e^{(1-sin(x))^2}+(x-y)^2$$
Bird Functionは多峰性の関数で、(x,y)=(4.70104, 3.15294),(-1.58214, -3.13024)に大域的最適解-106.764537、その周辺に多数の局所解を持つ。

sample_GA.py

import numpy as np
def bird(x):
    x1 = x[0]
    x2 = x[1]
    t = np.sin(x1)*np.exp((1-np.cos(x2))**2) + np.cos(x2)*np.exp((1-np.sin(x1))**2) + (x1-x2)**2 
    return t,

# meshgrid上で計算
c1 = c2 = np.arange(-10, 10, 0.1)
xx, yy = np.meshgrid(c1, c2)
length = np.size(c1)  
z = np.empty((length, length))
for i in range(length):
    for j in range(length):
        xin=np.array([xx[i,j], yy[i,j]])
        z[i,j], = bird(xin)

# 描画
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

fig0 = plt.figure(figsize=(12,12), dpi=80).add_subplot(111, projection='3d')
fig0.plot_surface(xx, yy, z, alpha=1)
plt.xlim(10, -10)            
plt.ylim(-10, 10)

ライブラリのインポート

必要なライブラリをインポートしておく。
DEAPはrandomライブラリを使用するため、1試行ごとに結果が変化する。
これを避けるためには、random.seed(1)とすること。

sample_GA.py

import pandas as pd
import random
from deap import algorithms
from deap import base,creator,tools
# random.seed(1)

初期個体群の生成

初期個体を生成する関数を作っておく。
ここでは、x, yともに-10〜10の範囲で、10×10の格子状に初期個体を生成した。

sample_GA.py

def initPopulation(pcls, ind_init):
    x = np.linspace(-10, 10, 10)
    y = np.linspace(-10, 10, 10)
    inix, iniy = np.meshgrid(x,y)
    contents = np.concatenate([inix.reshape([-1,1]), iniy.reshape([-1,1])], axis=1)
    pcls = [ind_init(c) for c in contents]
    return pcls

toolboxに登録

評価（"evaluate"）、個体群（"population"）、交叉（"mate"）、突然変異（"mutate"）、選択（"select"）等の手法をtoolboxに登録しておく。

例えば、下記のtoolbox.register("mate", tools.cxBlend, alpha=0.2)は、
tool.cxblendという関数を交叉手法として登録します、という意味（alpha=0.2はtools.cxBlendの引数）。

交叉方法にブレンド交叉、突然変異にガウス変異、選択形式は10個体のトーナメントで実施する。
初期個体群を指定したい場合、"population_guess"に登録しておく。

sample_GA.py

creator.create("FitnessMax", base.Fitness, weights=(-1.0,))
creator.create("Individual", np.ndarray, fitness=creator.FitnessMax)

toolbox = base.Toolbox()
toolbox.register("evaluate", bird)
toolbox.register("population_guess", initPopulation, list, creator.Individual)
toolbox.register("mate", tools.cxBlend, alpha=0.2)
toolbox.register("mutate", tools.mutGaussian, mu=[0.0, 0.0], sigma=[2.0, 2.0], indpb=1)
toolbox.register("select", tools.selTournament, tournsize=10)

実行

ここまで設定すれば、あとは簡単なコマンドで動く。

下記の例では各世代の評価個体数、平均値、標準偏差、最大値、最小値をlogbookに記録していく。
hofとは、Holl of Fameの略で、評価が高い順に指定した数だけ保存する事ができる。

algorithms.eaSimpleの引数cxpbは交叉率、mutpbは突然変異率、ngenは世代数を示す。
1世代あたりの個体数はpopで定義したものを引き継ぎ10×10=100個体となる。
初期サンプリングと世代交代時の個体数を別々に指定することも可能。

sample_GA.py

# 初期個体群を生成
pop = toolbox.population_guess()

# 設定
hof = tools.HallOfFame(5, similar=np.array_equal)
stats = tools.Statistics(lambda ind: ind.fitness.values)
stats.register("avg", np.mean)
stats.register("std", np.std)
stats.register("min", np.min)
stats.register("max", np.max)

# 実行
pop, logbook= algorithms.eaSimple(pop, toolbox,
                        cxpb=0.9, mutpb=0.1, ngen=20,
                        stats=stats, halloffame=hof, verbose=True)

# 上位3個体を表示
print("holl of fame.1 : %s" % (hof[0]))
print("holl of fame.2 : %s" % (hof[1]))
print("holl of fame.3 : %s" % (hof[2]))

terminal

gen nevals  avg     std     min         max   
0   100     80.0665 99.3005 -83.8352    414.98
1   74      7.49092 62.2371 -103.255    301.724
2   89      11.499  108.449 -105.693    665.594
3   84      -60.283 53.2721 -106.71     161.645
4   87      -95.3896    32.299  -106.723    40.9594
5   75      -99.3516    28.4967 -106.758    23.0967
6   73      -104.068    15.7185 -106.764    4.33984
7   76      -103.476    18.6955 -106.764    10.0824
8   80      -101.665    22.5172 -106.764    16.8155
9   92      -102.631    20.6472 -106.764    26.921 
10  77      -102.882    19.2791 -106.764    6.34586
11  90      -99.4555    30.4939 -106.764    56.3788
12  89      -100.566    27.1489 -106.764    30.2934
13  79      -100.978    25.2596 -106.764    51.745 
14  78      -98.4071    32.1796 -106.764    85.5625
15  76      -105.728    10.3096 -106.764    -3.14868
16  89      -95.2292    38.3427 -106.764    80.6272 
17  91      -102.44     25.6436 -106.764    96.6545 
18  80      -105.432    11.2501 -106.764    4.33866 
19  83      -102.271    23.504  -106.764    67.6912 
20  79      -103.856    16.5553 -106.764    -3.86946
holl of fame.1 : [4.70021671 3.1529326 ]
holl of fame.2 : [4.70021671 3.15293287]
holl of fame.3 : [4.70021671 3.15293358]

大域的最適解-106.765付近に収束した。
各個体の微分情報を使わないアルゴリズムなので、微小な誤差は出てしまう。

まとめと所感

• 遺伝的アルゴリズムDEAPで関数の最適化を行った
• 乱数を使うため、試行ごとに収束するかどうかが怪しい
• 離散値で多変数の最適化向けか？