前回の「遺伝的アルゴリズムによるレイリー減衰のパラメータ探査」で試した、一定減衰に相当するレイリー減衰のパラメータ探査を、今回はPythonの進化計算ライブラリのDEAPを使ってやってみます。
プログラムは、データサイエンス教育総合研究所のサイトを参考にさせていただきました。参考にしたというより、目的関数のところを一定減衰Shakeとレイリー減衰応答解析の差分2乗和とする関数に置き換えただけです。
前回は、遺伝的アルゴリズムをプログラミングしてパラメータを探しましたが、DEAPを使うとアルゴリズムの詳細をプログラミングしなくてもパラメータ探査できます。
実行結果は以下のとおり。
最も良い個体は [2.572249424150136, 6.2771604118935755]で、そのときの目的関数の値は (7.337222985431525,)
上記の個体の数値がレイリー減衰パラメータである二つの振動数になります。
このレイリー減衰のパラメータで応答計算し、評価位置で応答スペクトルを比較したものが下図です。
一定減衰とレイリー減衰の応答性状がよく合致しています。
参考までに、全く的外れなレーリー減衰のパラメータとして、0.02Hzと50Hzを指定した場合の比較は下図になります。
以下のプログラムでは、この0.02Hzから50Hzを遺伝子の範囲としています。
#個体の各遺伝子を決めるために使用
import random
import numpy as np
import pandas as pd
#DEAPの中にある必要なモジュールをインポート
from deap import base
from deap import creator
from deap import tools
from deap import algorithms
from FreeField import FreeField
#最小化問題として設定(-1.0で最小化、1.0で最大化問題)
#面食らうかもしれませんが、目的が最小化の場合は以下のように設定します。
creator.create("FitnessMin", base.Fitness, weights=(-1.0,))
#個体の定義(list型と指定しただけで、中身の遺伝子は後で入れる)
#これも面食らうかもしれませんが、個体を準備したと思ってください。
creator.create("Individual", list, fitness=creator.FitnessMin)
#目的関数の定義。必ずreturnの後に,をつける
#レイリー減衰の応答解析とターゲットとなる一定減衰Shakeの差分2乗和を計算する
def resp(individual):
#一定減衰Shakeのターゲットスペクトルの読み込み
df_target_spec = pd.read_csv("shake.csv")
#レイリー減衰応答解析のスペクトル計算
FreeField_spec = np.array([FreeField(individual).resp()])
object = np.sum((df_target_spec["shake"] - FreeField_spec[0]) ** 2)
return object,
#各種関数の設定を行います
#交叉、選択、突然変異などには、DEAPのToolbox内にある関数を利用
toolbox = base.Toolbox()
#random.uniformの別名をattribute関数として設定。各個体の遺伝子の中身を決める関数(各遺伝子は0.02Hz~50Hzのランダムな値)
toolbox.register("attribute", random.uniform, 0.02,50)
#individualという関数を設定。それぞれの個体に含まれる2個の遺伝子をattributeにより決めるよ、ということ。
toolbox.register("individual", tools.initRepeat, creator.Individual, toolbox.attribute, 2)
#集団の個体数を設定するための関数を準備
toolbox.register("population", tools.initRepeat, list, toolbox.individual)
#トーナメント方式で次世代に子を残す親を選択(tornsizeは各トーナメントに参加する個体の数)
toolbox.register("select", tools.selTournament, tournsize=5)
#交叉関数の設定。ブレンド交叉という手法を採用
toolbox.register("mate", tools.cxBlend,alpha=0.2)
#突然変異関数の設定。indpbは各遺伝子が突然変異を起こす確率。muとsigmaは変異の平均と標準偏差
toolbox.register("mutate", tools.mutGaussian, mu=[0.0, 0.0], sigma=[20.0, 20.0], indpb=0.2)
#評価したい関数の設定(目的関数のこと)
toolbox.register("evaluate", resp)
#以下でパラメータの設定
#今回は最も単純な遺伝的アルゴリズムの手法を採用
#乱数を固定
random.seed(64)
#何世代まで行うか
NGEN = 30
#集団の個体数
POP = 10
#交叉確率
CXPB = 0.9
#個体が突然変異を起こす確率
MUTPB = 0.1
#集団は80個体という情報の設定
pop = toolbox.population(n=POP)
#集団内の個体それぞれの適応度(目的関数の値)を計算
for individual in pop:
individual.fitness.values = toolbox.evaluate(individual)
#パレート曲線上の個体(つまり、良い結果の個体)をhofという変数に格納
hof = tools.ParetoFront()
#今回は最も単純なSimple GAという進化戦略を採用
algorithms.eaSimple(pop, toolbox, cxpb=CXPB, mutpb=MUTPB, ngen=NGEN, halloffame=hof)
#最終的な集団(pop)からベストな個体を1体選出する関数
best_ind = tools.selBest(pop, 1)[0]
#結果表示
print("最も良い個体は %sで、そのときの目的関数の値は %s" % (best_ind, best_ind.fitness.values))
以下、関連リンク
・DEAPの本家サイト