はじめに
平野廣美氏による1995年の論文クラスタ平均法を組み込んだ遺伝的アルゴリズムによるジョブショップスケジューリング問題の解法に記載の方法をPythonのDEAPモジュールを使用して実装した。
今回実装を5章に分けて紹介する。この記事は第3章である。
- ジョブショップスケジューリング問題について
- ジョブショップスケジューリング問題の解表現
- ジョブショップスケジューリング問題向きの遺伝的操作（←この記事）
- DEAPライブラリ
- メイン処理、実行結果とまとめ
なおソースコード全体は以下に置いてある。
なお、論文補足のため同氏による2000年初版の書籍 遺伝的アルゴリズムと遺伝的プログラミング オブジェクト指向フレームワークによる構成と応用 での記述および付属ソースコードを参考にしている。ただしクラス構成などはオリジナルとは異なる。オリジナルは洗練されており参照してほしい。
ジョブショップスケジューリング問題向きの遺伝的操作
前章では解表現からガントチャート（機械単位）およびメイクスパンを得る方法を説明した。
この章では乱数による操作を解候補に施して最適な解を探る手順（遺伝的操作）について説明する。
遺伝的アルゴリズムについてはWikipediaなどを手掛かりに調べてみてほしい。
解候補となる個体を数100用意し（この集団を世代とかpopulationなどと呼ぶ）それらに対して操作を行っていく。
なお今回は100個体用意した。
平野氏記載の遺伝的操作は、制約条件を満たさない致死遺伝子が出現しないよう工夫したものである。これら操作を「選択」、「交叉」、「突然変異」、「置換」に分けて記載する。詳細は各文献をみてほしい。処理の手順は次のようになる。
-
populationから2個体を「選択」する。 - 一定の確率で前項の2個体に「交叉処理」を施す。
- 一定の確率で前項の2個体それぞれに「突然変異処理」を施す。
- 前項の2個体を
population内の個体と「置換」する。
選択
ルーレット選択法である。
現世代から遺伝的操作候補となる個体を選択する方法の一つである。
前章で記載の方法で各個体のガントチャートを作成するとガントチャート別に最大完了時刻（メイクスパン）を得ることができる。今回の問題においてこの値が小さな個体がすぐれた個体であり大きな個体が劣った個体である。
ルーレット選択では、すぐれた個体をより高い確率で選択する方法である。
DEAPライブラリでは deap.tools.selRoulette()を使ってルーレット選択法による個体の選択が可能である。
ただし、ルーレット選択法は適応度が大きな個体ほど優れた個体とするものである。今回はメイクスパンの逆数を適応度とすることでルーレット選択を可能にする。
交叉
ルーレット選択法で2個体を選出する。
ある確率でこの2個体を混じらわせ新たな2個体を生み出すことができる。
今回は確率80%で交叉することにした。この値は参考文献のソースコード記載のデフォルト値であるが、論文記載の試行がどのような値であったのかは不明である。
交叉法は各個体の部分遺伝子を入れ替える2点交叉である。致死遺伝子を持たないようJSP向けに工夫したものである。ここではコードの解説のみとする。
交叉処理を行う関数crossoverは、ind1, ind2の2個体を引数にもつ。
def crossover ( ind1, ind2 ) :
"""JSP用の2点交叉処理"""
まずは2点交叉する位置つまり、部分遺伝子の開始位置cxpoint1と終了位置cxpoint2を取得する。
# ind1, ind2の長さは同じ
size = len ( ind1 )
cxpoint1, cxpoint2 = 0, 0
while cxpoint1 == cxpoint2 :
cxpoint1 = random.randrange ( 0, size )
cxpoint2 = random.randrange ( 0, size )
# 2点目が1点目より前なら入れ替える
if cxpoint2 < cxpoint1 :
cxpoint1, cxpoint2 = cxpoint2, cxpoint1
個体1, 2の部分遺伝子をsub1, sub2にセットしまた、交叉後の部分遺伝子格納用にnew_sub1, new_sub2を用意する。
なお各個体はarray.arrayを継承したクラスで実装していることから、部分遺伝子もarray.arrayを使用しているが、
listを使っても大差ないと思う。
# 部分遺伝子を取得
sub1 = array.array ( ind1.typecode, ind1 [ cxpoint1 : cxpoint2 ] )
sub2 = array.array ( ind2.typecode, ind2 [ cxpoint1 : cxpoint2 ] )
new_sub1 = array.array ( ind1.typecode, [ -1 ] * ( cxpoint2 - cxpoint1 ) )
new_sub2 = array.array ( ind2.typecode, [ -1 ] * ( cxpoint2 - cxpoint1 ) )
個体1の部分遺伝子と個体2の部分遺伝子と比較し共通のジョブ番号があれば相手側に複製する。
for s1_idx, s1 in enumerate ( ind1 [ cxpoint1 : cxpoint2 ] ) :
# ind1の部分遺伝子の要素がind2の部分遺伝子にみつからない
if s1 not in sub2 : continue
# みつかったら相手のnew_subに位置を保存してコピー
s2_idx = sub2.index ( s1 )
new_sub1 [ s2_idx ] = s1
new_sub2 [ s1_idx ] = s1
# コピーした要素は-1にしておく
sub1 [ s1_idx ] = -1
sub2 [ s2_idx ] = -1
コピーしなかった遺伝子があれば順序を保って元の個体に保存する。
以下は個体1の部分遺伝子の場合である。
個体2にコピーした遺伝子は-1をセットしてあるので-1以外ならば未コピーの遺伝子である。
コピー先で遺伝子が-1である位置は空いているためそこに未コピー遺伝子をセットする。
# コピーしなかった要素を順序を保存して戻す
for s1 in sub1 :
# -1の場合コピー済み
if s1 == -1 : continue
# new_sub1の未コピー位置を取得
idx = new_sub1.index ( -1 )
new_sub1 [ idx ] = s1
個体2についても同様である。
for s2 in sub2 :
if s2 == -1 : continue
idx = new_sub2.index ( -1 )
new_sub2 [ idx ] = s2
最後に交叉後の部分遺伝子をもとの個体にセットする
# 部分遺伝子を個体にセット
ind1 [ cxpoint1 : cxpoint2 ] = new_sub1
ind2 [ cxpoint1 : cxpoint2 ] = new_sub2
この関数全体をまとめると以下のようになる。
def crossover ( ind1, ind2 ) :
"""JSP用の2点交叉処理"""
# ind1, ind2の長さは同じ
size = len ( ind1 )
cxpoint1, cxpoint2 = 0, 0
while cxpoint1 == cxpoint2 :
cxpoint1 = random.randrange ( 0, size )
cxpoint2 = random.randrange ( 0, size )
# 2点目が1点目より前なら入れ替える
if cxpoint2 < cxpoint1 :
cxpoint1, cxpoint2 = cxpoint2, cxpoint1
# 部分遺伝子を取得
sub1 = array.array ( ind1.typecode, ind1 [ cxpoint1 : cxpoint2 ] )
sub2 = array.array ( ind2.typecode, ind2 [ cxpoint1 : cxpoint2 ] )
new_sub1 = array.array ( ind1.typecode, [ -1 ] * ( cxpoint2 - cxpoint1 ) )
new_sub2 = array.array ( ind2.typecode, [ -1 ] * ( cxpoint2 - cxpoint1 ) )
for s1_idx, s1 in enumerate ( ind1 [ cxpoint1 : cxpoint2 ] ) :
# ind1の部分遺伝子の要素がind2の部分遺伝子にみつからない
if s1 not in sub2 : continue
# みつかったら相手のnew_subに位置を保存してコピー
s2_idx = sub2.index ( s1 )
new_sub1 [ s2_idx ] = s1
new_sub2 [ s1_idx ] = s1
# コピーした要素は-1にしておく
sub1 [ s1_idx ] = -1
sub2 [ s2_idx ] = -1
# コピーしなかった要素を順序を保存して戻す
for s1 in sub1 :
# -1の場合コピー済み
if s1 == -1 : continue
# new_sub1の未コピー位置を取得
idx = new_sub1.index ( -1 )
new_sub1 [ idx ] = s1
for s2 in sub2 :
if s2 == -1 : continue
idx = new_sub2.index ( -1 )
new_sub2 [ idx ] = s2
# 部分遺伝子を個体にセット
ind1 [ cxpoint1 : cxpoint2 ] = new_sub1
ind2 [ cxpoint1 : cxpoint2 ] = new_sub2
return ind1, ind2
突然変異
前節で選択した2個体それぞれを対象とする。
ある確率で個体の遺伝子の値を変更する操作。
今回は確率50%で突然変異することにした。この値は参考文献のソースコード記載のデフォルト値であるが、論文記載の試行がどのような値であったのかは不明である。
致死遺伝子を持たないようJSP向けに工夫し、ランダムに選んだ2か所の値を入れ替える操作としている。ここではコードの解説のみとする。
突然変異を行う関数mutationは、ind個体を引数にもつ。
def mutation ( ind ) :
"""JSP用の突然変異処理"""
# 個体の長さを取得
size = len ( ind )
ランダムなpos1, pos2を選択しそれらの値を交換する。
理由は不明だが、参考文献ではこの操作を2度実施している。
# 理由は不明だがオリジナルソースでは2回実施している
for _ in range ( 2 ) :
pos1 = random.randrange ( 0, size )
pos2 = random.randrange ( 0, size )
ind [ pos1 ], ind [ pos2 ] = ind [ pos2 ], ind [ pos1 ]
なお、今回の個体はarray.arrayを継承したクラスのためスワップ操作を簡易に記述できるが、numpy.arrayを継承したクラスの場合のスワップ操作は一時変数を介する必要がある。
この関数全体をまとめると以下のようになる。
def mutation ( ind ) :
"""JSP用の突然変異処理"""
size = len ( ind )
# 理由は不明だがオリジナルソースでは2回実施している
for _ in range ( 2 ) :
pos1 = random.randrange ( 0, size )
pos2 = random.randrange ( 0, size )
ind [ pos1 ], ind [ pos2 ] = ind [ pos2 ], ind [ pos1 ]
return ind
通常の置換操作
前節、前々節の個体を既存のpopulationの個体と置き換える操作である。
置換する個体を選択する方法として2種類説明する。
通常の置換操作と、平野氏提案のクラスタ平均法（CAM）とを説明する。
通常の置換操作では、単純に既存populationの中で最も適応度が悪い個体（すなわちメイクスパンが最も大きい個体）と遺伝的操作の対象となった個体とを置き換えるものである。
通常の置換対象の個体を選択する関数getArgWorstは個体のリストpopulationと取得する個体数nとを引数に持つ。
def getArgWorst ( population, n ) :
"""
fitnessesの大きい方からn個のインデックスを取得
@param population individualのリスト
@param n 取得する個体数
@return nth worst indexes
"""
まず適応度が悪い個体n個を選択する。
s_inds = getWorst ( population, n )
getWorst()の実装は以下の通りである。
DEAPフレームワークでは適応度(fitness)は、各個体のfitness.values属性に格納している。
次のようにして適応度の降順にした個体群のコピーの状態n個体を取得する。
def getWorst ( population, n ) :
"""
fitnessesの大きい方からn個の個体を取得
@param population individualのリスト
@param n 取得する個体数
@return nth worst individuals
"""
return sorted ( population, key=lambda x: x.fitness.values, reverse=True )[ : n ]
つぎにこれらn個体の元個体群でのindexを取得する。
selected = []
for ind in s_inds :
selected.append ( population.index ( ind ) )
この関数全体をまとめると以下のようになる。
def getWorst ( population, n ) :
"""
fitnessesの大きい方からn個の個体を取得
@param population individualのリスト
@param n 取得する個体数
@return nth worst individuals
"""
return sorted ( population, key=lambda x: x.fitness.values, reverse=True )[ : n ]
def getArgWorst ( population, n ) :
"""
fitnessesの大きい方からn個のインデックスを取得
@param population individualのリスト
@param n 取得する個体数
@return nth worst indexes
"""
s_inds = getWorst ( population, n )
selected = []
for ind in s_inds :
selected.append ( population.index ( ind ) )
return selected
クラスタ平均法（CAM）による置換操作
平野氏が提案したクラスタ平均法（CAM）による置換操作では、各個体をその先頭1遺伝子の値でクラスターに分類し、最大のクラスター内で最も適応度が悪い個体と遺伝的操作対象となった個体とを置き換えるものである。
CAMによる置換対象の個体を選択する関数getArgWorstCAMは個体のリストpopulationと取得する個体数nとを引数に持つ。
def getArgWorstCAM ( population, n ) :
"""クラスタ平均化法(Cluster Averaging Method)により適応度が悪い個体をn個選択
@param population individualのリスト
@param n 個体選択数
@return 選択した個体インデックスリスト
"""
まず各個体の先頭遺伝子の値とキー、このクラスタに属する個体のリストを値とする辞書(clusters)を後述するgetClusters()を用いて取得する。
そして最大のクラスター（max_cluster）およびそのサイズ（max_cluster_sz）そして最小のクラスターサイズ（min_cluster_sz）を得る。
# 各個体の先頭遺伝子別のクラスターを取得
clusters = getClusters ( population )
max_cluster = max ( clusters.values(), key=len )
max_cluster_sz = len ( max_cluster )
min_cluster_sz = len ( min ( clusters.values(), key=len ) )
なおgetClusters()の実装は以下のようになっている。
def getClusters ( population ) :
clusters = {}
for ind in population :
# 各個体の先頭遺伝子だけを取得
if ind [ 0 ] not in clusters :
clusters [ ind [ 0 ] ] = []
clusters [ ind [ 0 ] ].append ( ind )
return clusters
そして、最大クラスターと最小クラスターの大きさの差が40未満（40という数値に説明はないが経験的なものだと考える）の場合クラスターを無視して全体から個体を選択し、40以上の場合最大のクラスターから個体を選択する。
getArgWorst()およびgetWorst()は前節で紹介している。
全体から個体を選択して置換した場合、最大クラスターと最小クラスターとの大きさの差は、大きくも小さくもなり得る。
一方、最大クラスターの個体を置換した場合、最大クラスターの大きさは変わらないか小さくなるかである。つまり最大クラスターと最小クラスターの大きさの差がこれ以上広がることはない。
つまりCAMを採用すれば、クラスタサイズの上限があることで多様なクラスタを保つことができるようになる。
# 最大のクラスターと最小のクラスターと個体数の差が小さければ全体から探す
if ( max_cluster_sz - min_cluster_sz ) < 40 :
selected = getArgWorst ( population, n )
# 最大のクラスターと最小のクラスターと個体数の差が大きければ最大クラスターから探す
else :
# 最大のクラスターの個体の適応度下位n件を取得
s_inds = getWorst ( max_cluster, n )
for ind in s_inds :
selected.append ( population.index ( ind ) )
ここで紹介した関数全体をまとめると以下のようになる。
def getClusters ( population ) :
clusters = {}
for ind in population :
# 各個体の先頭遺伝子だけを取得
if ind [ 0 ] not in clusters :
clusters [ ind [ 0 ] ] = []
clusters [ ind [ 0 ] ].append ( ind )
return clusters
def getArgWorstCAM ( population, n ) :
"""クラスタ平均化法(Cluster Averaging Method)により適応度が悪い個体をn個選択
@param population numpy.array 個体リスト。順序は適応度の小さい順に変わる
@param n 個体選択数
@return 選択した個体インデックスリスト
"""
selected = []
# 各個体の先頭遺伝子別のクラスターを取得
cluster = getClusters ( population )
max_cluster = max ( cluster.values(), key=len )
max_cluster_sz = len ( max_cluster )
min_cluster_sz = len ( min ( cluster.values(), key=len ) )
# 最大のクラスターと最小のクラスターと個体数の差が小さければ全体から探す
if ( max_cluster_sz - min_cluster_sz ) < 40 :
selected = getArgWorst ( population, n )
# 最大のクラスターと最小のクラスターと個体数の差が大きければ最大クラスターから探す
else :
# 最大のクラスターの個体の適応度下位n件を取得
s_inds = getWorst ( max_cluster, n )
for ind in s_inds :
selected.append ( population.index ( ind ) )
return selected
遺伝的操作の呼び出し
上述した「選択」「交叉」「突然変異」「置換」の遺伝的操作を行うtest2()関数は、populationを引数に持つ。
def test2 ( population ) :
""" populationに遺伝的操作を施す """
# 交叉確率、突然変異確率
CXPB, MUTPB = 0.8, 0.5
まずpopulationから２個体を選択する。indsリストに格納する。
ここでgToolboxはglobalに保存したツールボックスオブジェクトである。各操作を抽象化した名前で実行することが可能になる。DEAPライブラリが提供している。
gToolbox.select()には別途選択関数を登録しておりpopulationから2個体を選出する。
gToolbox.cloneは複製操作であり、選出した2個体をそれぞれ複製する。
# idx1, idx2 をルーレット選択し複製
inds = list ( map ( gToolbox.clone, gToolbox.select ( population, 2 ) ) )
一定の確率で交叉処理を実施する。
gToolbox.mate()には別途交叉関数を登録しており引数の2個体に交叉処理を施す。
遺伝的操作を施したことを記録するため、個体の適応度を無効にしておく。
# 交叉確率の割合で交叉処理を実施
if random.random() < CXPB :
gToolbox.mate ( inds [ 0 ], inds [ 1 ] )
# 操作した個体の適応度を無効にする
del inds [ 0 ].fitness.values
del inds [ 1 ].fitness.values
次に、選択した2個体それぞれに対して遺伝的操作を施す。
# 選択した個体それぞれに操作
for ind in inds :
# indに遺伝的操作
一定の確率で突然変異処理を実施する。
gToolbox.mutate()には別途突然変異関数を登録しており引数の個体に処理を施す。
遺伝的操作を施したことを記録するため、個体の適応度を無効にしておく。
# 突然変異確率の割合で突然変異処理を実施
if random.random() < MUTPB :
gToolbox.mutate ( ind )
# 操作した個体の適応度を無効にする
del ind.fitness.values
遺伝的操作を施した個体の適応度を再計算する。
gToolbox.evaluate()には別途評価関数を登録しており引数の個体のメイクスパンを取得する。
# 適応度が無効である個体を再評価
if not ind.fitness.valid :
ind.fitness.values = gToolbox.evaluate ( ind )
最後に適応度が最も悪い個体（メイクスパンが最も大きい）を1個体選択し置換する。
gToolbox.getArgWors()には別途関数を登録しており、引数のpopulationから選択した個体のpopulationでのindexを取得する。
最後に個体を置換する。
# 既存の個体と置換
worst_idx = gToolbox.getArgWorst ( population, 1 )[ 0 ]
population [ worst_idx ] = ind
この関数全体をまとめると以下のようになる。
def test2 ( population ) :
""" populationに遺伝的操作を施す """
# 交叉確率、突然変異確率
CXPB, MUTPB = 0.8, 0.5
# idx1, idx2 をルーレット選択し複製
inds = list ( map ( gToolbox.clone, gToolbox.select ( population, 2 ) ) )
# 交叉確率の割合で交叉処理を実施
if random.random() < CXPB :
gToolbox.mate ( inds [ 0 ], inds [ 1 ] )
# 操作した個体の適応度を無効にする
del inds [ 0 ].fitness.values
del inds [ 1 ].fitness.values
# 選択した個体それぞれに操作
for ind in inds :
# 突然変異確率の割合で突然変異処理を実施
if random.random() < MUTPB :
gToolbox.mutate ( ind )
# 操作した個体の適応度を無効にする
del ind.fitness.values
# 適応度が無効である個体を再評価
if not ind.fitness.valid :
ind.fitness.values = gToolbox.evaluate ( ind )
# 既存の個体と置換
worst_idx = gToolbox.getArgWorst ( population, 1 )[ 0 ]
population [ worst_idx ] = ind
return population
以上の操作をpopulation/2回繰り返す。
def do_generation ( population ) :
for _ in range ( len ( population ) // 2 ) :
test2 ( population )
return population
do_generation()での処理を１世代とする。世代ごとに最良の個体を保存する。3000世代ほど繰り返し世代を通して最良の個体を準最適解とする。
この章で紹介した遺伝的操作はscheduleモジュールに格納している。モジュール全体は以下リンクにある。
また遺伝的操作を実行するtest2()などは以下mainモジュールに格納している。
(つづく）
次の章「DEAPライブラリ」はこちら
前の章「ジョブショップスケジューリング問題の解表現」はこちら