Introduction to Design and Implementation of Metaheuristics

メタヒューリスティクスの設計と実装⼤阪⼤学⼤学院情報科学研究科梅⾕俊治 2023年1⽉19⽇

はじめに • 実務に現れる組合せ最適化問題には，汎⽤の数理最適化ソルバーでは対応できない問題が少なくありません． • これらの問題では，貪欲法や局所探索法にさまざまな⼿法を組み合わせたメタヒューリスティクスを開発します． • しかし，メタヒューリスティクスを詳細に解説した書籍は⾮常に少なく，実装のノウハウを習得することは容易ではありません．
• 本講演では，巡回セールスマン問題と⼀般化割当問題を通じて，局所探索法とメタヒューリスティクスの実装を解説します． 2 注意1︓本講演では多くの聴衆が親しんでいるPythonでサンプルコードを記述していますが，Pythonによる実装を推奨しているわけではありません．Pythonで実装したコードは実⾏速度が遅いので，C/C++など実⾏速度の速いプログラミング⾔語での実装をおすすめします．注意2︓本講演で紹介する実装は⼀例であり，特に性能が優れているわけでも実装が洗練されているわけでもないことに注意して下さい．注意3︓本講演で紹介するプログラムは説明のため簡略化しています．実際のプログラムとは異なり，そのままでは動作しないことに注意して下さい．

講演の概要 • 局所探索法の設計と実装 ü 巡回セールスマン問題(TSP) ü TSP: 最近傍法(NN), 2-opt法，Or-opt法，3-opt法 ü
TSP: 局所探索法の⾼速化 ü ⼀般化割当問題(GAP) ü GAP: 制約付き問題に対するペナルティ関数法 ü GAP: シフト近傍⾛査，スワップ近傍⾛査 • メタヒューリスティクスの設計と実装 ü TSP: ランダム多スタート局所探索法(MLS,GRASP) ü TSP: 反復局所探索法(ILS) ü TSP: 遺伝的アルゴリズム(GA,MA) ü TSP: アニーリング法(SA) ü TSP: タブー探索法(TS) ü TSP: 誘導局所探索法(GLS)，Fast Local Search(FLS) ü GAP: 制約付き問題に対する重み付き局所探索法(WLS) • 数値実験(TSP, GAP) 3

巡回セールスマン問題(TSP) • 都市の集合と各2都市間の移動コスト(距離など)が与えられたとき，全ての都市をちょうど⼀度ずつ訪問し，出発した都市に戻る巡回路の中で総移動コストが最⼩となるものを求める問題． • ユークリッドTSP︓都市はd次元空間内の点，都市間の移動コストは2点間のユークリッド距離で与えられる． • △TSP︓都市間の移動コストが三⾓不等式を満たす．
• 対称TSP︓都市間の移動コストは向きによらず同じ． • ⾮対称TSP︓都市間の移動コストが向きにより異なっても良い． 5 ⼀般特殊 u v w dist(u,w) <= dist(u,v) + dist(v,w) 三⾓不等式基板のドリル⽳あけ計画，配送計画などの応⽤を持つ

巡回セールスマン問題の局所探索法の実装 • ⼊⼒データの読込み • 初期解の⽣成，解の評価，解の出⼒ • 近傍⾛査 6 ⼊⼒データの読込み最近近傍法
2-opt法 Or-opt法 3-opt法実⾏結果の出⼒最近近傍法⼊⼒データの読込み最近近傍法実⾏結果の出⼒ 2-opt法 2-opt法⼊⼒データの読込み最近近傍法実⾏結果の出⼒局所探索法改善解が得られたら 2-opt法にもどる

⼊⼒データの読込み • TSPLIB*から最適値が既知の問題例40題(都市数100〜783) • 各都市(u,v)間のユークリッド距離を整数に丸めた値を距離とする． 7 *http://comopt.ifi.uni-heidelberg.de/software/TSPLIB95/ NAME : eil51
COMMENT : 51-city problem (Christofides/Eilon) TYPE : TSP DIMENSION : 51 EDGE_WEIGHT_TYPE : EUC_2D NODE_COORD_SECTION 1 37 52 2 49 49 3 52 64 (中略) 48 25 55 49 48 28 50 56 37 51 30 40 EOF def calc_dist(u,v): xd = float(coord[u][0] - coord[v][0]) yd = float(coord[u][1] - coord[v][1]) return int(math.sqrt(xd * xd + yd * yd) + 0.5) TSPLIBの問題例局所探索法の実⾏結果都市(u,v)間の距離計算四捨五⼊している各都市の番号，x座標，y座標都市数

最近近傍法* • 出発する都市を適当に選ぶ． • いま訪問している都市から最も近い未訪問の都市に移動する． • 全ての都市を訪問したら出発した都市に戻る． 8 *データマイニングにおけるk-近傍法と名前が似ているので注意すること．

最近近傍法の実装 • 巡回路を配列tourに格納する*．tour[i]=vはi番⽬に都市vを訪問することを表す． • i-1番⽬までに訪問する都市が決定済みとして，都市tour[i-1]から最も近い未訪問の都市を⾛査する． 9 *都市数が⾮常に多い⼤規模な問題例では，two-level tree
tourと呼ばれるデータ構造に格納する． *2都市間の距離をメモリに記憶すると計算効率が改善しますが，今回は実施していません． *説明のためコードを簡略化していますので，そのままでは動作しないことに注意して下さい． # initialize tour tour = [i for i in range(num_node)] # nearest neighbor for i in range(1,num_node): # find nearest unvisited node min_dist = float('inf') arg_min_dist = None for j in range(i,num_node): dist = calc_dist(tour[i-1],tour[j]) if dist < min_dist: min_dist = dist arg_min_dist = j # set nearest unvisited node tour[i], tour[arg_min_dist] = tour[arg_min_dist], tour[i] 0 1 2 3 4 5 6 7 0 1 2 3 4 5 6 7 0 1 3 2 4 5 6 7 tour[i-1]に最も近い都市を⾛査 tour[i]とtour[arg_min_dist]を⼊れ替え最近近傍法の⼿続き*

k-opt法，Or-opt法 • k-opt法︓現在の巡回路から辺をk本交換する． • Or-opt法︓部分路を別の場所に挿⼊する． 10 *時間枠制約付き配送計画問題ではOr-opt法を使うことが多い． 3-opt近傍の例 Or-opt近傍の例 2-opt近傍の例

2-opt近傍⾛査の実装 • (u,next_u),(v,next_v)を(u,v),(next_u,next_v)につなぎ替える． • 改善解が得られたら始めから近傍を⾛査し直す(即時移動戦略)*． 11 *いつも始めの都市から⾛査すると探索に偏りが⽣じるので，tour[i]の次の都市から⾛査する，もしくは，⾛査順をランダムシャッフルすると性能が向上することがあります． *nbhdはリストではなくジェネレータ式を⽤いることに注意して下さい．C⾔語であれば2重のfor⽂を書いて，goto⽂を⽤いてネストしたループから脱出させたりします．
restart = True while restart: restart = False nbhd = ((i,j) for i in range(num_node) for j in range(i+2,num_node)) for i,j in nbhd: # evaluate difference delta = eval_diff(tour,i,j) if delta < 0: # change tour change_tour(tour,i,j) restart = True break 2-opt近傍⾛査の⼿続き* tour[i] u tour[i+1] next_u v tour[j] next_v tour[j+1] tour[i] u tour[i+1] next_u v tour[j] next_v tour[j+1] リストではなくジェネレータ式改善解が得られたら始めから近傍を⾛査し直す

2-opt近傍⾛査の実装 • (u,next_u),(v,next_v)を(u,v),(next_u,next_v)につなぎ替える． • 巡回路の変化量はつなぎ替えた辺の⻑さの差分を計算すればOK． • 巡回路の更新はtour[i+1]からtour[j]の部分を反転する． 12 *巡回路の更新の計算⼿間はO(n)なので，巡回路を双⽅向リストに格納すれば計算効率が⼤幅に改善しそうですが，巡回路が更新される頻度は⾮常に少ないので思ったほどには改善しません．
def eval_diff(tour,i,j): u, next_u = tour[i],tour[(i+1) % num_node] v, next_v = tour[j],tour[(j+1) % num_node] cur = calc_dist(u,next_u) + calc_dist(v,next_v) new = calc_dist(u,v) + calc_dist(next_u, next_v) return new - cur 巡回路⻑の差分計算 def change_tour(tour,i,j): # reverse sub-path [i+1,...,j] tour[i+1:j+1] = list(reversed(tour[i+1:j+1])) # update objective value obj += eval_diff(tour,i,j) 巡回路の更新* tour[i] u v tour[j] tour[i+1] next_u next_v tour[j+1] next_u v v next_u tour[i+1]からtour[j]の部分を反転する tour[n-1]の次はtour[0] であることに注意する．

Or-opt近傍⾛査の実装 • (prev_p,head_p),(tail_p,next_p),(v,next_v)を (prev_p,next_p),(v,head_p),(tail_p,next_v)につなぎ替える． • パラメータsで部分路(head_p,...,tail_p)の⻑さを指定する (sizeが最⼤値) 13 *実際には部分路を逆向きにつなぎ替える近傍操作も同時にチェックしますがここでは省略します． restart
= True while restart: restart = False nbhd = ((s,i,j) for s in range(1,size+1) for i in range(num_node) for j in range(i+s,i+num_node-1)) for s,i,j in nbhd: # evaluate difference delta = eval_diff(tour,s,i,j) if delta < 0: # change tour change_tour(tour,s,i,j) restart = True break Or-opt近傍⾛査の⼿続き tour[i] head_p tour[i+s-1] tail_p prev_p next_p tour[j+1] next_v head_p tail_p prev_p next_p next_v v tour[j] v

Or-opt近傍⾛査の実装 • (tour[i],...,tour[i+s-1])と(tour[i+s],...,tour[j])を⼊れ替える． 14 def eval_diff(tour,s,i,j): head_p, tail_p = tour[i],
tour[(i+s-1) % num_node] prev_p, next_p = tour[(i-1) % num_node], tour[(i+s) % num_node] v, next_v = tour[j % num_node], tour[(j+1) % num_node] cur = calc_dist(prev_p,head_p) + calc_dist(tail_p,next_p) + calc_dist(v,next_v) new = calc_dist(prev_p,next_p) + calc_dist(v,head_p) + calc_dist(tail_p,next_v) return new - cur 巡回路⻑の差分計算 def change_tour(tour,s,i,j): # get sub-path [i,...,i+s-1] subpath = [] for h in range(s): subpath.append(tour[(i+h) % num_node]) # move sub-path [i,...,i+s-1] to j+1 for h in range(i+s,j+1): tour[(h-s) % num_node] = tour[h % num_node] for h in range(s): tour[(j+1-s+h) % num_node] = subpath[h] # update objective value obj += eval_diff(tour,s,i,j) 巡回路の更新 i i+s-1 tour[i]からtour[i+s-1]の部分を切り出す j j+1 i+s i j-s i j-s j-s+1 j j+1

3-opt近傍⾛査の実装 • (u,next_u),(v,next_v),(w,next_w)のつなぎ替えが4通りある． type1: (u,next_v),(w,next_u),(v,next_w) type2: (u,w),(next_v,next_u),(v,next_w) type3: (u,next_v),(w,v),(next_u,next_w) type4:
(u,v),(next_u,w),(next_v,next_w) 15 u tour[i] next_u tour[i+1] v tour[j] tour[j+1] next_v w tour[k] next_w tour[k+1] u next_u v next_v w next_w u next_u v next_v w next_w u next_u v next_v w next_w u v next_v w next_w next_u type1 type2 type3 type4

3-opt近傍⾛査の実装 • (u,next_u),(v,next_v),(w,next_w)のつなぎ替えが4通りある． type1: (u,next_v),(w,next_u),(v,next_w) type2: (u,w),(next_v,next_u),(v,next_w) type3: (u,next_v),(w,v),(next_u,next_w) type4:
(u,v),(next_u,w),(next_v,next_w) 16 restart = True while restart: restart = False nbhd = ((i,j,k) for i in range(num_node) for j in range(i+2,num_node) for k in range(j+2,num_node)) for i,j,k in nbhd: # evaluate difference delta, oper = eval_diff(tour,i,j,k) if delta < 0: # change tour change_tour(tour,i,j,k,oper) restart = True break 3-opt近傍⾛査の⼿続き operでtype1-4を指定 def eval_diff(tour, i, j, k): best, arg_best = float('inf'), None u, next_u = tour[i], tour[(i+1) % num_node] v, next_v = tour[j], tour[(j+1) % num_node] w, next_w = tour[k], tour[(k+1) % num_node] cur = calc_dist(u,next_u) + calc_dist(v,next_v) + calc_dist(w,next_w) new = calc_dist(u,next_v) + calc_dist(v,next_w) + calc_dist(w,next_u) if new - cur < best: best, arg_best = new - cur, 'type1' new = calc_dist(u,w) + calc_dist(next_v,next_u) + calc_dist(v,next_w) if new - cur < best: best, arg_best = new - cur, 'type2' new = calc_dist(u,next_v) + calc_dist(w,v) + calc_dist(next_u,next_w) if new - cur < best: best, arg_best = new - cur, 'type3' new = calc_dist(v,u) + calc_dist(next_w,next_v) + calc_dist(w,next_u) if new - cur < best: best, arg_best = new - cur, 'type4' return best, arg_best 巡回路⻑の差分計算 type1-4で最も良かった近傍操作を返す

3-opt近傍⾛査の実装 • type1: (tour[i+1],...,tour[j])と(tour[j+1],...,tour[k])を⼊れ替え． 17 *巡回路の⻑さの計算⼿間はO(n)なので，差分計算で置き換えれば計算効率が⼤幅に改善しそうですが，巡回路が更新される頻度は⾮常に少ないので思ったほどには改善しません． i i+1
k+1 k j j+1 u next_u v next_v w next_w i k+1 巡回路の更新* tour[i] tour[i+1] tour[j] tour[j+1] tour[k] tour[k+1] def change_tour(tour,i,j,k,oper): if oper == 'type1’: tour[i+1:k+1] = tour[j+1:k+1] + tour[i+1:j+1] elif oper == 'type2': tour[i+1:k+1] = list(reversed(tour[j+1:k+1])) + tour[i+1:j+1] elif oper == 'type3': tour[i+1:k+1] = tour[j+1:k+1] + list(reversed(tour[i+1:j+1])) elif oper == 'type4': tour[i+1:k+1] = list(reversed(tour[i+1:j+1])) + list(reversed(tour[j+1:k+1])) # update objective value obj = tour_length(tour)

3-opt近傍⾛査の実装 • type2: (tour[j+1],...,tour[k])を反転し(tour[i+1],...,tour[j])と⼊れ替え． 18 i i+1 k+1 k
j j+1 i k+1 u next_u v next_v w next_w tour[i] tour[i+1] tour[j] tour[j+1] tour[k] tour[k+1] 反転する def change_tour(tour,i,j,k,oper): if oper == 'type1’: tour[i+1:k+1] = tour[j+1:k+1] + tour[i+1:j+1] elif oper == 'type2': tour[i+1:k+1] = list(reversed(tour[j+1:k+1])) + tour[i+1:j+1] elif oper == 'type3': tour[i+1:k+1] = tour[j+1:k+1] + list(reversed(tour[i+1:j+1])) elif oper == 'type4': tour[i+1:k+1] = list(reversed(tour[i+1:j+1])) + list(reversed(tour[j+1:k+1])) # update objective value obj = tour_length(tour) 巡回路の更新

3-opt近傍⾛査の実装 • type3: (tour[i+1],...,tour[j])を反転し(tour[j+1],...,tour[k]) と⼊れ替え． 19 i i+1 k+1 k
j j+1 i k+1 def change_tour(tour,i,j,k,oper): if oper == 'type1’: tour[i+1:k+1] = tour[j+1:k+1] + tour[i+1:j+1] elif oper == 'type2': tour[i+1:k+1] = list(reversed(tour[j+1:k+1])) + tour[i+1:j+1] elif oper == 'type3': tour[i+1:k+1] = tour[j+1:k+1] + list(reversed(tour[i+1:j+1])) elif oper == 'type4': tour[i+1:k+1] = list(reversed(tour[i+1:j+1])) + list(reversed(tour[j+1:k+1])) # update objective value obj = tour_length(tour) 巡回路の更新 u next_u v next_v w next_w tour[i] tour[i+1] tour[j] tour[j+1] tour[k] tour[k+1] 反転する

3-opt近傍⾛査の実装 • type4: (tour[i+1],...,tour[j])と(tour[j+1],...,tour[k])をそれぞれ反転． 20 i i+1 k+1 k
j j+1 i i+1 k+1 k j j+1 def change_tour(tour,i,j,k,oper): if oper == 'type1’: tour[i+1:k+1] = tour[j+1:k+1] + tour[i+1:j+1] elif oper == 'type2': tour[i+1:k+1] = list(reversed(tour[j+1:k+1])) + tour[i+1:j+1] elif oper == 'type3': tour[i+1:k+1] = tour[j+1:k+1] + list(reversed(tour[i+1:j+1])) elif oper == 'type4': tour[i+1:k+1] = list(reversed(tour[i+1:j+1])) + list(reversed(tour[j+1:k+1])) # update objective value obj = tour_length(tour) 巡回路の更新反転する反転する u v next_v w next_w next_u tour[i] tour[i+1] tour[j] tour[j+1] tour[k] tour[k+1]

局所探索法の⾼速化 • Or-opt法と3-opt法は，2-opt法より計算⼿間が⼤きいので，改善解が得られた時点で即座に2-opt法に戻る実装をすることも多い． 21 def or_opt_search(tour, size=OR_OPT_SIZE): nbhd =
((s,i,j) for s in range(1,size+1) for i in range(num_node) for j in range(i+s,i+num_node-1)) for s,i,j in nbhd: # evaluate difference delta = eval_diff(tour,s,i,j) if delta < 0: # change tour change_tour(tour,s,i,j) return True Or-opt近傍⾛査の⼿続き # nearest neighbor nearest_neighbor(tour) # local search while True: # 2-opt neighborhood search two_opt_search(tour) # Or-opt neighborhood search if or_opt_search(tour): continue # 3-opt neighborhood search if three_opt_search(tour): continue break 局所探索法の⼿続き改善解が得られたら即座に2-opt近傍の⾛査に戻る改善解が得られたら即座に Or-opt近傍の⾛査を終了する

隣接リストを⽤いた局所探索法の⾼速化 • 最適巡回路は近隣の都市間を繋ぐ枝を含む場合が多い． • 各都市に対して近隣にある都市を列挙して隣接リストに保持する． • 平⾯上の巡回セールスマン問題であれば，k-d⽊などのデータ構造を⽤いて近隣にある都市を⾼速に検索できる． 22 1
2 3 4 5 6 7 8 9 10 11 12 13 15 16 14 領域分割を繰返して都市をクラスタに分割する 1,2,4,5 3,16,15,14 6,7,10,11 8,9,13,12 cutdim=x cutval=x[14] cutdim=y cutval=y[2] cutdim=y cutval=y[6] ⽣成されたk-d⽊ 2,3,5,16 *久保幹雄「巡回セールスマン問題に対するLin-Kernighan近傍を⽤いたGuided Local Search」より引⽤ https://www.logopt.com/mikiokubo/presen/gls-tsp.ppt

隣接リストを⽤いた局所探索法の実装 • 隣接リストneighborは前処理で作成する． • 隣接リストを⽤いた近傍探索では，⾛査順が巡回路の訪問順に沿わないため実装に若⼲の修正が必要になる． • 巡回路の更新のため，配列tourを逆参照する配列posを⽤意する． pos[v]=iは都市vをi番⽬に訪問することを表す． 23
def gen_neighbor(): neighbor = [[] for _ in range(num_node)] for i in range(num_node): temp = [(calc_dist(i,j),j) for j in range(num_node) if j != i] temp.sort(key=lambda x: x[0]) neighbor[i] = [temp[h][1] for h in range(min(NB_LIST_SIZE,num_node))] 隣接リストの作成 pos = [i for i in range(len(tour))] for i in range(len(tour)): pos[tour[i]] = i 配列posの作成 0 1 2 3 4 5 6 7 0 1 3 2 6 7 5 4 pos 0 v 1 2 3 4 5 6 7

隣接リストを⽤いた2-opt近傍⾛査の実装 • 巡回路のi,j番⽬の都市ではなく，都市u,vで近傍を⾛査する． • 巡回路の更新時に配列posもあわせて更新する． • Or-opt法も同様の修正で実装できる． 24 restart =
True while restart: restart = False nbhd = ((u,v) for u in tour for v in neighbor[u]) for u,v in nbhd: # evaluate difference delta = eval_diff(tour,u,v) if delta < 0: # change tour change_tour(tour,pos,u,v) restart = True break 2-opt法の⼿続き def eval_diff(tour,u,v): next_u, next_v = tour_next(u), tour_next(v) cur = calc_dist(u,next_u) + calc_dist(v,next_v) new = calc_dist(u,v) + calc_dist(next_u,next_v) return new - cur 巡回路⻑の差分計算 def tour_next(v): return tour[(pos[v]+1) % num_node] 都市vの次に訪問する都市 def change_tour(tour,pos,u,v): if pos[u] < pos[v]: i,j = pos[u],pos[v] else: i,j = pos[v],pos[u] # reverse sub-path [i+1,...,j] tour[i+1:j+1] = list(reversed(tour[i+1:j+1])) # update positions for h in range(num_node): pos[tour[h]] = h # update objective value obj = tour_length(tour) 巡回路の更新

隣接リストを⽤いた3-opt近傍⾛査の実装 • type1,3,4とtype2で近傍⾛査を分ける必要が⽣じる︖ 25 restart = True while restart: restart
= False # scan 3-opt neighborhood (type1,3,4) nbhd = ((u,v,w) for u in tour for v in neighbor[u] for w in neighbor[tour_next(u)]) for u,v,w in nbhd: # evaluate difference delta, oper = eval_diff_type134(tour,u,v,w) if delta < 0: # change tour change_tour(tour,u,v,w,oper) restart = True break # scan 3-opt neighborhood (type2) nbhd = ((u,v,w) for u in tour for v in neighbor[tour_next(u)] for w in neighbor[u]) for u,v,w in nbhd: # evaluate difference delta, oper = eval_diff_type2(tour,u,v,w) if delta < 0: # change tour change_tour(tour,u,v,w,oper) restart = True break 3-opt法の⼿続き type1 u next_u type2 u next_u u next_u type3 u next_u type4 v w v w w v v w type1,3,4では都市vが都市uに隣接 type2では都市wが都市uに隣接 prev_v next_w prev_v next_w prev_v prev_w next_v next_w

隣接リストを⽤いた3-opt近傍⾛査の実装 26 def eval_diff_type134(tour, pos, u, v, w): best, arg_best
= float('inf'), None cur = calc_dist(u,tour_next(u)) + calc_dist(tour_prev(v),v) + calc_dist(w,tour_next(w)) new = calc_dist(u,v) + calc_dist(tour_prev(v),tour_next(w)) + calc_dist(w,tour_next(u)) if new - cur < best and pos[v] >= pos[u]+3 and pos[w] >= pos[v]+1: best, arg_best = new - cur, 'type1' cur = calc_dist(u,tour_next(u)) + calc_dist(tour_prev(v),v) + calc_dist(tour_prev(w),w) new = calc_dist(u,v) + calc_dist(tour_prev(w),tour_prev(v)) + calc_dist(tour_next(u),w) if new - cur < best and pos[v] >= pos[u]+3 and pos[w] >= pos[v]+2: best, arg_best = new - cur, 'type3' cur = calc_dist(u,tour_next(u)) + calc_dist(v,tour_next(v)) + calc_dist(w,tour_next(w)) new = calc_dist(u,v) + calc_dist(tour_next(u),w) + calc_dist(tour_next(v),tour_next(w)) if new - cur < best and pos[v] >= pos[u]+2 and pos[w] >= pos[v]+2: best, arg_best = new - cur, 'type4' return best, arg_best 巡回路⻑の差分計算(type1,3,4) def eval_diff_type2(tsp, pos, u, v, w): cur = calc_dist(u,tour_next(u)) + calc_dist(tour_prev(v),v) + calc_dist(w,tour_next(w)) new = calc_dist(u,w) + calc_dist(v,tour_next(u)) + calc_dist(tour_prev(v),tour_next(w)) if pos[v] >= pos[u]+3 and pos[w] >= pos[v]+1: return new - cur, 'type2' else: return float('inf'), None 巡回路⻑の差分計算(type2) 巡回路上でuはv より3つ以上先⾏， vはwに先⾏する type1,3,4で最も良かった近傍操作を返す

隣接リストを⽤いた3-opt近傍⾛査の実装 • type1: (tour[i+1],...,tour[j-1])と(tour[j],...,tour[k])を⼊れ替え． 27 def change_tour(tour,pos,u,v,w,oper): i,j,k =
pos[u],pos[v],pos[w] if oper == 'type1': tour[i+1:k+1] = tour[j:k+1] + tour[i+1:j] elif oper == 'type2': tour[i+1:k+1] = list(reversed(tour[j:k+1])) + tour[i+1:j] elif oper == 'type3': tour[i+1:k] = tour[j:k] + list(reversed(tour[i+1:j])) elif oper == 'type4': tour[i+1:k+1] = list(reversed(tour[i+1:j+1])) + list(reversed(tour[j+1:k+1])) # update positions for h in range(num_node): pos[tour[h]] = h # update objective value obj = tour_length(tour) 巡回路の更新 i i+1 k+1 k j-1 j u next_u prev_v v w next_w i k+1 tour[i] tour[j-1] tour[j] tour[k] tour[k+1] tour[i+1]

隣接リストを⽤いた3-opt近傍⾛査の実装 • type2: (tour[j],...,tour[k])を反転して(tour[i+1],...,tour[j-1]) と⼊れ替え． • type3: (tour[i+1],...,tour[j-1])を反転して(tour[j],..., tour[k-1])と⼊れ替え． 28
i i+1 k+1 k j-1 j i k+1 u next_u prev_v v w next_w tour[i] 反転する tour[i+1] tour[j-1] tour[j] tour[k] tour[k+1] i i+1 k k-1 j-1 j i k u next_u prev_v v prev_w w tour[i] tour[i+1] tour[j-1] tour[j] tour[k-1] tour[k] 反転する

⼀般化割当問題(GAP) • 与えられたn個の仕事をm個のエージェントに割り当てたとき，割当にともなう総費⽤を最⼩化する問題． • 仕事jをエージェントiに割り当てたときのコストc[i,j]と資源の要求量a[i,j]，各エージェントiの利⽤可能資源量b[i]が与えられる． • 実⾏可能解の有無を判定する決定問題がNP完全で，実⾏可能解を 1つ求めることすら容易ではない．
30 1 1 2 2 3 c[1,1]=2 a[1,1]=2 4,2 b[1]=4 2,2 3,1 5,2 問題例実⾏可能解 1 2 3 4,2 2,3 3,1 3 2 1 2,3 3 3 1 2

⼀般化割当問題(GAP) • 与えられたn個の仕事をm個のエージェントに割り当てたとき，割当にともなう総費⽤を最⼩化する問題． • 仕事jをエージェントiに割り当てたときのコストc[i,j]と資源の要求量a[i,j]，各エージェントiの利⽤可能資源量b[i]が与えられる． • 実⾏可能解の有無を判定する決定問題がNP完全で，実⾏可能解を 1つ求めることすら容易ではない．
31 エージェントiの容量制約仕事jの割当制約⼀般化割当問題の定式化

制約付き問題に対するペナルティ関数法 • 実⾏可能解を1つ求めることも容易ではないので，制約条件を緩和して実⾏不可能解も探索空間に含める． • エージェントiの容量制約に対する違反度p[i]に重みw[i]をかけた値をペナルティとして⽬的関数に加える． • ペナルティ重みw[i]の調整は容易ではなく，適当な値を与えて局所探索法を1回適⽤するだけは質の⾼い実⾏可能解は得られない*．
32 エージェントiの容量制約の違反量p[i]を追加 *ペナルティ重みの適応的な調整は，後述の重み付き局所探索法で紹介します．

制約付き問題に対するペナルティ関数法 • 実⾏可能解を1つ求めることも容易ではないので，制約条件を緩和して実⾏不可能解も探索空間に含める． • エージェントiの容量制約に対する違反度p[i]に重みw[i]をかけた値をペナルティとして⽬的関数に加える． • ペナルティ重みw[i]の調整は容易ではなく，適当な値を与えて局所探索法を1回適⽤するだけは質の⾼い実⾏可能解は得られない*．
33 *ペナルティ重みの適応的な調整は，後述の重み付き局所探索法で紹介します．実⾏不可能領域実⾏不可能領域実⾏不可能領域最⼩化最⼩化最⼩化適切なペナルティ重みw[i]を決定することは困難ペナルティ重みが⼤きいペナルティ重みが⼩さいペナルティ重みがちょうど良い

⼀般化割当問題の局所探索法の実装 • ⼊⼒データの読込み • 初期解の⽣成，解の評価，解の出⼒ • 近傍⾛査 34 シフト近傍⾛査スワップ近傍⾛査
⼊⼒データの読込みランダム割当実⾏結果の出⼒局所探索法 # generate initial solution for j in range(num_job): best_sol[j] = random.randrange(0,num_agent) # initialize current solution sol = best_sol[:] # local search while True: # shift neighborhood search shift_nb_search(best_sol,sol) # swap neighborhood search if swap_nb_search(best_sol,sol): continue break 局所探索法の⼿続き改善解が得られたら即座にシフト近傍の⾛査に戻る仕事をランダムに割り当て* *貪欲法などを⽤いた初期解の⽣成も提案されている．いずれの⽅法にせよ，必ずしも実⾏可能解が得られるわけではないので，ここでは簡単なランダム割当を⽤いた． sol i j sol[j]=iは仕事jがエージェントi に割り当てられることを表す

シフト近傍，スワップ近傍 • シフト近傍︓仕事jの割当先のエージェントをi1からi2に移動する． • スワップ近傍︓仕事j1とj2の割当先のエージェントを交換する． • シフト近傍よりスワップ近傍の⽅が容量制約を違反しづらい． • シフト近傍の⼤きさO(mn) <
スワップ近傍の⼤きさO(n2)． 35 シフト近傍の例 j j スワップ近傍の例 j i1 i2 x[i1,j]=1→0 x[i2,j]=0→1 i1 i2 j1 j1 j2 j2 i1 i2 j1 j1 j2 j2 x[i1,j1]=1 x[i2,j2]=1 x[i1,j1]=1→0 x[i2,j2]=1→0 x[i2,j1]=0→1 x[i1,j2]=0→1 i1 j sol i2 j sol i1 i2 sol j1 j2 i2 i1 sol j1 j2

シフト近傍⾛査の実装 • 最良の実⾏可能解(暫定解)best_solと現在の解solを更新する． • 各エージェントiの使⽤資源量used[i]を保持すると，差分をO(1) の⼿間で計算できる． 36 シフト近傍⾛査の⼿続き j sol[j]
i x[sol[j],j]=1 x[i,j]=0 j sol[j] i x[sol[j],j]=0 x[i,j]=1 restart = True while restart: restart = False nbhd = ((j,i) for j in range(num_job) for i in range(num_agent) if i != sol[j]) for j,i in nbhd: # calculate difference delta_obj, delta_plt = calc_diff(sol,j,i) if (plt + delta_plt < NUM_EPSILON and (best_plt > NUM_EPSILON or obj + delta_obj < best_obj - NUM_EPSILON)): # update incumbent solution best_sol = sol[:] update_sol(best_sol,j,i) if delta_obj + delta_plt < -NUM_EPSILON: # update current solution update_sol(sol,j,i) restart = True break 最良の実⾏可能解も更新する* *best_solが実⾏可能解でない場合の対応が必要なことに注意する． j j

シフト近傍⾛査の実装 • 最良の実⾏可能解(暫定解)best_solと現在の解solを更新する． • 各エージェントiの使⽤資源量used[i]を保持すると，差分をO(1) の⼿間で計算できる． 37 def calc_diff(sol,j,i,flag): i1,i2
= sol[j],i delta_obj = cost[i2][j] - cost[i1][j] cur_plt = max(0, used[i1] - cap[i1]) new_plt = max(0, used[i1] - res[i1][j] - cap[i1]) delta_plt = wt[i1] * (new_plt - cur_plt) cur_plt = max(0, used[i2] - cap[i2]) new_plt = max(0, used[i2] + res[i2][j] - cap[i2]) delta_plt += wt[i2] * (new_plt - cur_plt) return delta_obj, delta_plt 評価関数の差分計算 def update_sol(sol,j,i): i1,i2 = sol[j],i sol[j] = i2 used[i1] -= res[i1][j] used[i2] += res[i2][j] obj += cost[i2][j] - cost[i1][j] plt = 0.0 for i in range(num_agent): plt += wt[i] * max(0, used[i] – cap[i]) 割当の更新 used[i1] cap[i1] cap[i2] used[i2] j j i1 i2 x[i1,j]=1 x[i2,j]=0 j x[i1,j]=0 x[i2,j]=1 j i1 i2 used[i1] -res[i1][j] used[i2] +res[i2][j]

スワップ近傍⾛査の実装 38 nbhd = ((j1,j2) for j1 in range(num_job) for
j2 in range(j1+1,num_job) if sol[j2] != sol[j1]) for j1,j2 in nbhd: # calculate difference delta_obj, delta_plt = calc_diff(sol,j1,j2) if (plt + delta_plt < NUM_EPSILON and (best_plt > NUM_EPSILON or obj + delta_obj < best_obj - NUM_EPSILON)): # update incumbent solution best_sol = sol[:] update_sol(best_sol,j1,j2) if delta_obj + delta_plt < -NUM_EPSILON: # update current solution update_sol(sol,j1,j2) return True return False スワップ近傍⾛査の⼿続き i1=sol[j1] i2=sol[j2] j1 j1 j2 j2 x[i1,j1]=1 x[i2,j2]=1 i1=sol[j1] i2=sol[j2] j1 j1 j2 j2 x[i1,j1]=0 x[i2,j2]=0 x[i2,j1]=1 x[i1,j2]=1 改善解が得られたら即座にスワップ近傍の⾛査を終了する最良の実⾏可能解も更新する* *best_solが実⾏可能解でない場合の対応が必要なことに注意する．

スワップ近傍⾛査の実装 39 def calc_diff(sol,j1,j2,flag): i1,i2 = sol[j1],sol[j2] delta_obj = cost[i2][j1]
+ cost[i1][j2] - cost[i1][j1] - cost[i2][j2] cur_plt = max(0, used[i1] - cap[i1]) new_plt = max(0, used[i1] - res[i1][j1] + res[i1][j2] - cap[i1]) delta_plt = wt[i1] * (new_plt - cur_plt) cur_plt = max(0, used[i2] - cap[i2]) new_plt = max(0, used[i2] - res[i2][j2] + res[i2][j1] - cap[i2]) delta_plt += wt[i2] * (new_plt - cur_plt) return delta_obj, delta_plt 評価関数の差分計算 def update_sol(sol,j1,j2): i1,i2 = sol[j1],sol[j2] sol[j1], sol[j2] = i2, i1 used[i1] += res[i1][j2] - res[i1][j1] used[i2] += res[i2][j1] - res[i2][j2] obj += cost[i2][j1] + cost[i1][j2] - cost[i1][j1] - cost[i2][j2] plt = 0.0 for i in range(num_agent): plt += wt[i] * max(0, used[i] – cap[i]) 割当の更新 i1 i2 j1 j1 j2 j2 x[i1,j1]=1 x[i2,j2]=1 i1 i2 j1 j1 j2 j2 x[i1,j1]=0 x[i2,j2]=0 x[i2,j1]=1 x[i1,j2]=1 cap[i1] cap[i2] used[i1] used[i2] used[i1] -res[i1][j1] +res[i1][j2] used[i2] -res[i2][j2] +res[i2][j1]

メタヒューリスティクスの実装 • 貪欲法や局所探索法に様々な⼿法を組合せて探索の集中化と多様化をバランス良く実現する． • 複数の初期解に対して局所探索法を適⽤する． → ランダム多スタート局所探索法，GRASP，反復局所探索法，遺伝的アルゴリズム • 改悪解への移動を許容して局所最適解に留まらず探索を継続する．
→ アニーリング法，タブー探索法 • 評価関数を適応的に変形して局所最適解に留まらず探索を継続する． → 誘導局所探索法，重み付き局所探索法 41 ⼊⼒データの読込み初期解の⽣成局所探索法実⾏結果の出⼒反復局所探索法局所探索法⼊⼒データの読込み初期解の⽣成局所探索法実⾏結果の出⼒初期解の⽣成と局所探索法を繰り返し適⽤する

ランダム多スタート局所探索法 • ランダムに⽣成した複数の初期解から局所探索法を実⾏し，得られた解の中で最良のものを出⼒する． • 各試⾏が独⽴なため，局所探索法の適⽤回数を増やしてもより良い解が得られる効果は急速に⼩さくなる． • ランダム性を加えた貪欲法を⽤いて初期解を⽣成する⼿法は GRASP(Greedy
Randomized Adaptive Search Procedure)と呼ぶ． 42 複数の初期解に局所探索法を適⽤解空間⽬的関数値

ランダム多スタート局所探索法の実装 • 2回⽬以降の局所探索法はランダム順列で初期解を作成． • 暫定解を保持して局所探索法が終わる度に更新する． • 局所最適解に到達するまでの時間が⻑く，局所探索法の反復回数が少ない． 43 #
nearest neighbor nearest_neighbor(tour) obj = tour_length(tour) start_time = cur_time = time.time() while cur_time - start_time < time_limit: # generate random tour random.shuffle(cur_tour) # local search local_search(cur_tour) cur_obj = tour_length(cur_tour) # update best tour if cur_obj < obj: tour = cur_tour[:] obj = cur_obj cur_time = time.time() ランダム多スタート法の⼿続きランダム順列を⽣成

GRASPの実装 • i-1番⽬までに訪問する都市が決定済みとして，都市tour[i-1]から TOP_K_SIZE番⽬までに近い未訪問の都市をランダムに選んで⾛査する． • ヒープを⽤いて都市tour[i-1]からTOP_K_SIZE番⽬までに近い未訪問の都市を格納する．ヒープの根に最も距離の遠い都市を格納する． 44 def random_nearest_neighbor(tour):
# initialize tour tour = [i for i in range(num_node)] # randomized k-nearest neighbor for i in range(1,num_node): # initialize heap top_k = [] # find top-k nearest unvisited node for j in range(i,num_node): dist = calc_dist(tour[i-1],tour[j]) if len(top_k) < TOP_K_SIZE: heapq.heappush(top_k, (-dist, j)) elif -dist > top_k[0][0]: heapq.heapreplace(top_k, (-dist, j)) # set randomly selected unvisited node from top-k idx = top_k[random.randrange(min(TOP_K_SIZE,len(top_k)))][1] tour[i], tour[idx] = tour[idx], tour[i] # clear heap top_k.clear() GRASPの⼿続きヒープ内の最も遠い都市と⼊替えヒープ内の都市からランダムに選択 TOP_SIZE_K個以下なら都市を追加* TOP_SIZE_K個ならばなら最も遠い都市と⼊替え* ヒープ *追加・⼊替えられた要素は直後に⾏われるヒープの更新により位置が変わることに注意．

反復局所探索法 • 過去の探索で得られた良い解にランダムな変形を加えたものを初期解として局所探索法を適⽤する⼿続きを繰り返す． • 逆戻りを防ぐために，局所探索法とは異なる近傍操作を⽤いて解に変形を加えることが望ましい． • 変形時に⽤いる近傍の⼤きさを適応的に変化させる可変近傍探索法も知られる．
45 複数の初期解に局所探索法を適⽤過去の探索で得られた良い解に変形を加える解空間⽬的関数値局所探索では2-opt近傍変形にdouble bridge近傍

反復局所探索法の実装 • 2回⽬以降の局所探索法は暫定解に摂動を加えて初期解を作成． • double-bridge近傍操作を適⽤して暫定解に摂動を加える*． 46 # nearest neighbor nearest_neighbor(tour)
obj = tour_length(tour) start_time = cur_time = time.time() while cur_time - start_time < time_limit: # kick incumbent tour cur_tour = tour[:] kick_tour(cur_tour) # local search local_search(cur_tour) cur_obj = tour_length(cur_tour) # update best tour if cur_obj < obj: tour = cur_tour[:] obj = cur_obj cur_time = time.time() 反復局所探索法の⼿続き def kick_tour(tour): i = random.randrange(0,len(tour)-7) j = random.randrange(i+2,len(tour)-5) k = random.randrange(j+2,len(tour)-3) l = random.randrange(k+2,len(tour)-1) tour[i+1:l+1] = tour[k+1:l+1] + tour[j+1:k+1] + tour[i+1:j+1] 巡回路の摂動 tour[i] tour[i+1] tour[j] tour[j+1] tour[k] tour[k+1] tour[l] tour[l+1] i i+1 j j+1 k k+1 l l+1 i l+1 *2-opt近傍操作を適⽤すると直後の局所探索法ですぐに巻き戻ってしまう．

反復局所探索法の実装 • 2回⽬以降の局所探索法では，局所最適解に到達するまでの時間が短く，局所探索法の反復回数が多い*． • 常に暫定解を初期解の種とするのではなく，アニーリング法のアイデアを⽤いて初期解の種となる解を更新すると性能を改善できる*． 47 *ランダム多スタート局所探索法の数倍程度に反復回数が増える． *単に現在の解を摂動すると暫定解から⼤きく離れてしまい質の良い解が得られない．
*ここでは省略しているが，実際には暫定解tourも更新している． seed_tour = tour[:] seed_obj = tour_length(seed_tour) while cur_time - start_time < time_limit: # kick incumbent tour cur_tour = seed_tour[:] kick_tour(cur_tour) # local search local_search(cur_tour) cur_obj = tour_length(cur_tour) # set temperature temp = set_temp(cur_time) # update seed tour delta = cur_obj - seed_obj if delta < 0 or math.exp(-delta/temp) > random.random(): seed_tour = cur_tour[:] cur_time = time.time() 反復局所探索法の⼿続き(修正版)* 改悪でも確率 exp(-delta/temp)で seed_tourを更新

遺伝的アルゴリズム • 複数の解を集団として保持することで探索の多様化を実現する． • 交叉や突然変異などの操作を適⽤して新たな解を⽣成し，選択により次世代の集団を決定する⼀連の⼿続きを繰り返し適⽤する． • より⾃由な枠組みに基づいて複数の解から新たな解を⽣成する散布局所探索法やパス再結合法などが知られる． 48
複数の初期解に局所探索法を適⽤過去の探索で得られた良い解に変形を加える集団新たな集団次世代の集団交叉突然変異選択新しい個体を局所探索法で改善する⽅法は，遺伝的局所探索法もしくはミームアルゴリズムと呼ばれる．

遺伝的アルゴリズムの実装 • 交叉︓巡回路tour1から適当な⻑さの部分路subpathをランダムに切り取り巡回路tour2のランダムな位置に挿⼊する． • 巡回路tour2から部分路subpathに含まれる都市は取り除き，含まれない都市は元の巡回路の順序を保つように並べる． 49 1 2
3 4 5 6 7 8 6 3 8 2 7 1 4 5 6 8 7 2 3 4 1 5 def crossover(tour1, tour2, new_tour): # cut subpath from tour1 length = int(random.uniform(MIN_PATH_RATIO,MAX_PATH_RATIO) * num_node) head = random.randint(0,num_node-length) subpath = tour1[head:head+length] # generate a new_tour new_tour = [None for _ in range(num_node)] head = random.randint(0,num_node-length) new_tour[head:head+length] = subpath k = 0 for i in range(num_node): if tour2[i] not in subpath: while new_tour[k] is not None: k += 1 new_tour[k] = tour2[i] 交叉の⼿続き tour1 tour2 new_tour 順序交叉の変形 Or-opt操作の⼀般化

遺伝的アルゴリズムの実装 • 複製選択︓巡回路⻑を無視して集団から2個体をランダムに⾮復元抽出する． • ⽣存選択︓両親と⼦から最良の1個体とランダムに選択した1個体を次世代に残す(家族内選択)． 50 *遺伝的アルゴリズムの提案からまもない時期に局所探索法を組み込んだ実装も提案されている． *家族内選択では⼀世代で2個体しか⼊替えが⽣じない．⼀般的な世代交代とはイメージが異なる．
# initialize population population = [[] for _ in range(POP_SIZE)] init_population(population) while cur_time - start_time < time_limit: # cross-over tours to generate new tour tour1 = population.pop(random.randrange(len(population))) tour2 = population.pop(random.randrange(len(population))) cur_tour = [i for i in range(num_node)] crossover(tour1, tour2, cur_tour) # local search local_search(cur_tour) cur_obj = tour_length(cur_tour) # update best tour if cur_obj < obj: tour = cur_tour[:] obj = cur_obj # update population update_population(population, tour1, tour2, cur_tour) cur_time = time.time() 遺伝的アルゴリズムの⼿続き* def update_population(population, tour1, tour2, cur_tour): cand = [tour1, tour2, cur_tour] best_obj, arg_best = float('inf'), None for k in range(len(cand)): if cand[k].obj < best_obj: best_obj, arg_best = cand[k].obj, k population.append(cand[arg_best]) cand.pop(arg_best) arg_rand = random.randrange(len(cand)) population.append(cand[arg_rand]) 集団の更新 def init_population(population): for k in range(POP_SIZE): # generate random tour population[k] = [i for i in range(num_node)] random.shuffle(population[k]) # local search local_search(population[k]) 初期集団の⽣成各個体に局所探索法を適⽤

アニーリング法 • 現在の解の近傍からランダムに選んだ解が改善解ならば移動し，改悪解ならば⽬的関数の値の変化量に応じた確率で移動することで，質の悪い局所最適解から脱出し，より良い解を探索する． • 物理現象の焼きなましにアイデアを借りて，温度と呼ばれるパラメータにより改悪解に移動する確率を適応的に調整する． •
探索の開始時にはランダムな移動が⽣じ易いように温度を⾼めに設定し，探索が進むにつれて徐々に温度を低くする． 51 改悪解への移動を許容して局所探索法を再開解空間⽬的関数値アニーリング法の解の移動確率

アニーリング法の実装 • 複数の近傍を順に⾛査し解候補をランダムサンプリングする． • 改悪でも確率で現在の解cur_tourを更新する. • とすれば，t=0で , t=1で
となる幾何冷却スケジューリングを実現できる*． 52 # nearest neighbor nearest_neighbor(tour) obj = tour_length(tour) start_time = cur_time = time.time() while cur_time - start_time < time_limit: # set temperature temp = set_temp(cur_time) # random neighborhood operation nb_set = ['two_opt','or_opt','three_opt_type134', 'three_opt_type2'] for nb in nb_set: if eval(nb)(cur_tour,temp): break # update best tour cur_obj = tour_length(cur_tour) if cur_obj < obj: tour = cur_tour[:] cur_time = time.time() アニーリング法の⼿続き u = random.choice(tour) v = random.choice(neighbor[u]) # evaluate difference delta = eval_diff(tour, u, v) if delta < 0 or math.exp(-delta/temp) > random.random(): # change tour change_tour(tour, u, v) return True else: return False 2-opt法の⼿続き改悪でも確率 exp(-delta/temp)で現在の解を更新 def set_temp(cut_time): t = (cur_time – start_time) / time_limit return math.pow(max_temp, 1.0 – t) * math.pow(min_temp, t) 温度パラメータの設定開始時start_timeにmax_temp, 終了時time_limitにmin_tempとなる *線形冷却スケジューリングはやや性能が落ちる．

アニーリング法の実装 • なら現在の解よりだけ改悪となる解をの確率で受理． • 他⼿法とは挙動が異なり，開始直後は暫定解から離れた解を探索し，中盤から終盤にかけて徐々に暫定解を更新する*． • 開始前に2-opt近傍内の解候補を適当な数ランダムサンプリングし，
, (>0)とする*． 53 def init_temp(tour): num_sample = NB_LIST_SIZE * num_node min_delta, avg_delta = float('inf'), 0.0 cnt = 0 for k in range(num_sample): u = random.choice(tour) v = random.choice(neighbor[u]) next_u, next_v = tour_next(u), tour_next(v) delta = calc_dist(u,v) + calc_dist(next_u,next_v) – calc_dist(u,next_u) – calc_dist(v,next_v) if delta > 0: avg_delta += delta cnt += 1 if delta < min_delta: min_delta = delta avg_delta /= float(cnt) return min_delta, avg_delta 温度パラメータの最⼤・最⼩値の設定 *終了直前に暫定解を更新することが多いので，time_limitより少し早い時間にとなるように冷却スケジュールを修正した⽅が良いかも知れない． *例えばとすると性能が悪くなる．求解性能が温度設定に⼤きく左右されることに注意． delta > 0となる解候補のみサンプリング

タブー探索法 • 現在の解を除く近傍内の最良解に移動する⼿続きを繰り返すことで局所最適解から脱出する． • 最近に探索した解に戻る巡回を防ぐために，タブーリストと呼ばれる移動を禁⽌する解の集合を⽤意し，現在の解とタブーリストTを除く近傍内の最良解に移動する． • 最近に探索した解そのものではなく，近傍操作により変化した変数
(とその値の組)をタブーリストに記憶することが多い． 54 改悪解への移動を許容して局所探索法を再開探索の履歴に応じて探索⽅向を制御解空間⽬的関数値 T[1]={} T[2]={1} T[3]={1,2} T[4]={1,2,3} T[5]={2,3,4} T[6]={3,4,5} tenure=3 tabu list t=1 t=6

タブー探索法の実装 • 規則1︓現在の解cur_tourの更新時に追加された辺(u,v)を記憶し， (u,v)を削除する操作を⼀定期間禁⽌する． • 規則2︓現在の解cur_tourの更新時に削除された辺(u,v)を記憶し， (u,v)を追加する操作を⼀定期間禁⽌する． • tenure回の反復後に操作を許可する旨を辞書moveに設定する*． 55
*無限ループが⽣じづらいようにtenureを[1,sqrt(num_node)]の範囲の⼀様乱数に設定する． def set_tabu(move, u, v): tenure = random.randint(1,int(math.sqrt(num_node))+1) if u <= v: move[(u,v)] = cnt + tenure else: move[(v,u)] = cnt + tenure タブー期間の設定(規則1) def get_tabu(move, u, v): if u > v: u,v = v,u if (u,v) in move and cnt <= move[(u,v)]: return True else: return False タブー期間の確認(規則1) def update_tabu(u, v): set_tabu(u,v) set_tabu(tour_next(u), tour_next(v)) def check_tabu(u, v): if get_tabu(u,tour_next(u)) or get_tabu(v,tour_next(v)): return True else: return False 2-opt法のタブー設定(規則1) 2-opt法のタブー確認(規則1) (u,v), (next_u,next_v) の追加を記憶 (u,next_u),(v,next_v) の削除を禁⽌

タブー探索法の実装 • 規則1︓現在の解cur_tourの更新時に追加された辺(u,v)を記憶し， (u,v)を削除する操作を⼀定期間禁⽌する． • 規則2︓現在の解cur_tourの更新時に削除された辺(u,v)を記憶し， (u,v)を追加する操作を⼀定期間禁⽌する． • tenure回の反復後に操作を許可する旨を辞書moveに設定する*． 56
*無限ループが⽣じづらいようにtenureを[1,NB_LIST_SIZE*sqrt(num_node)]の範囲の⼀様乱数に設定する． def set_tabu(move, u, v): tenure = random.randint(1,int(NB_LIST_SIZE * math.sqrt(num_node))+1) if u <= v: move[(u,v)] = cnt + tenure else: move[(v,u)] = cnt + tenure タブー期間の設定(規則2) def get_tabu(move, u, v): if u > v: u,v = v,u if (u,v) in move and cnt <= move[(u,v)]: return True else: return False タブー期間の確認(規則2) def update_tabu(u, v): set_tabu(u,tour_next(u)) set_tabu(v, tour_next(v)) def check_tabu(u, v): if get_tabu(u,v) or get_tabu(tour_next(u),tour_next(v)): return True else: return False 2-opt法のタブー設定(規則2) 2-opt法のタブー確認(規則2) (u,next_u), (v,next_v) の削除を記憶 (u,v),(next_u,next_v) の追加を禁⽌

タブー探索法の実装 • 特別選択基準︓暫定解を更新するならタブー期間でも操作を実⾏する． • 複数の近傍を順に⾛査するが改善解が得られたら即座に2-optに戻る． 57 min_delta, arg_min_delta = float('inf'),
None nbhd = ((u,v) for u in cur_tour for v in neighbor[u]) for u,v in nbhd: asp_flag = False # evaluate difference delta = eval_diff(cur_tour, u, v) # update incumbent tour and check aspiration if cur_obj + delta < tour_obj: asp_flag = True tour = cur_tour[:] obj = cur_obj # update best solution in the neighborhood if delta < min_delta and (not check_tabu(u,v) or asp_flag): min_delta, arg_min_delta = delta, (u,v) if arg_min_delta is not None: # update tabu-list update_tabu(*arg_min_delta) # update current tour change_tour(cur_tour, *arg_min_delta) cnt += 1 if min_delta < 0: return True else: return False 2-opt法の⼿続き start_time = cur_time = time.time() while cur_time - start_time < time_limit: while True: # 2-opt neighborhood search if two_opt_search(tour, cur_tour): continue # Or-opt neighborhood search if or_opt_search(tour, cur_tour): continue # 3-opt neighborhood search if three_opt_search(tour, cur_tour): continue break cur_time = time.time() タブー探索法の⼿続き改善解が得られたら即座に2-opt近傍の⾛査に戻る* (u,v)がタブーリストに含まれない or 特別選択基準を満たす *同じ2-opt近傍の⾛査でもタブーリストが変化すると同じ結果にはならないことに注意する．

誘導局所探索法 • ⽬的関数とは異なる評価関数を適応的に変化させつつ局所探索法を繰り返し適⽤する． • 直前の局所探索法で得られた局所最適解の構成要素の中でコストが⼤きな要素にペナルティを加えて評価関数を更新することで，局所最適解から脱出する． • ⽬的関数に制約条件に対するペナルティ関数を加えた評価関数を⽤
いて，制約条件のペナルティ重みの更新と局所探索法の適⽤を交互に適⽤する重み付け法も知られる． 58 ⽬的関数と異なる評価関数の下で局所探索法を実⾏探索の履歴に応じて評価関数を適応的に変形解空間評価関数値評価関数値解空間

誘導局所探索法の実装 • 局所探索法では，各都市間の距離dist[u,v]ではなくペナルティ付き距離pdist[u,v]=dist[u,v] + alpha * penalty[u,v]で評価する． • ペナルティpenalty[u,v]の初期値は0とし，局所探索法の後に巡回路
の中でdist[v,next_v]/(1+penalty[v,next_v])が最⼤の(v,next_v) についてpenalty[v,next_v]を1増加する． 59 while cur_time - start_time < time_limit: # set coefficient of penalized distance alpha = PENALTY_RATIO * float(tour_length(tour)) / float(num_node) # local search algorithm local_search(tour, cur_tour) # update penalty update_penalty(cur_tour) cur_time = time.time() 誘導局所探索法の⼿続き def update_penalty(tour): max_val, arg_max_val = 0.0, None for i in range(num_node): v, next_v = tour[i], tour[(i+1) % num_node] val = float(calc_dist(v,next_v)) / float(1.0 + penalty[v,next_v]) if val > max_val: max_val, arg_max_val = val, (v,next_v) v, next_v = arg_max_val penalty[v,next_v] += 1.0 penalty[next_v,v] += 1.0 ペナルティ重み更新 alphaは巡回路に含まれる辺の平均距離*PENALTY_RATIO 距離が⻑くペナルティ値が⼩さい辺を選ぶ

誘導局所探索法の実装 • 近傍内の解候補を距離dist[u,v]でも評価し暫定解tourを更新する． • ペナルティ値の更新後の局所探索法による巡回路の変化が⼩さい • 直前の局所探索法から変化した部分にのみ近傍⾛査を集中するFast Local Searchを併⽤する*． 60
restart = True while restart: restart = False nbhd = ((u,v) for u in tour for v in neighbor[u]) for u,v in nbhd: # evaluate difference in original distance delta = eval_diff(cur_tour, u, v, 'dist') if cur_obj + delta < obj - NUM_EPSILON: # update incumbent tour tour = cur_tour[:] change_tour(tour, u, v) obj = calc_length(tour) # evaluate difference in penalized cost delta = eval_diff(cur_tour, u, v, 'pdist') if delta < -NUM_EPSILON: # change current tour change_tour(cur_tour, u, v) cur_obj = calc_plength(tour) restart = True break 2-opt近傍の⾛査 def calc_pdist(u,v): return calc_dist(u,v) + alpha * penalty[u,v] def eval_diff(tour, u, v, flag): next_u, next_v = tour_next(u), tour_next(v) if flag == 'pdist': cur = calc_pdist(u,next_u) + calc_pdist(v,next_v) new = calc_pdist(u,v) + calc_pdist(next_u,next_v) return new - cur else: cur = calc_dist(u,next_u) + calc_dist(v,next_v) new = calc_dist(u,v) + calc_dist(next_u,next_v) return new - cur ペナルティ付き距離の計算巡回路⻑の差分計算現在の解だけではなく暫定解も更新する *Fast Local Searchは巡回セールスマン問題ではdon't look bitとも呼ばれる．

Fast Local Searchの実装* • 直前の局所探索法から変化した部分にのみ近傍⾛査を集中するため，フラグactive[v]=Trueとなる都市vを起点とする近傍のみ⾛査する． • 局所探索法の後に，摂動やペナルティ値の更新が⽣じた辺(u,v)とその隣接辺のフラグactiveはTrue，それ以外はFalseに設定する． •
近傍⾛査で改善解が得られなければactive[v]=Falseとする． 61 *誘導局所探索法に限らず，反復局所探索法など⼩さな摂動を適⽤した後の局所探索法で効果的に働く． restart = True while restart: restart = False for u in tour: if active[u]: # 2-opt neighborhood search if two_opt_search(tour, cur_tour, u): restart = True break # Or-opt neighborhood search if or_opt_search(tour, cur_tour, u): restart = True break # 3-opt neighborhood search if three_opt_search(tour, cur_tour, u): restart = True break # inactivate node u active[u] = False 局所探索法の⼿続き def update_penalty(tour): max_val, arg_max_val = 0.0, None for i in range(num_node): v, next_v = tour[i], tour[(i+1) % num_node] val = float(calc_dist(v,next_v)) / float(1.0 + penalty[v,next_v]) if val > max_val: max_val, arg_max_val = val, (v,next_v) v, next_v = arg_max_val penalty[v,next_v] += 1.0 penalty[next_v,v] += 1.0 # set active nodes for v in tour: active[v] = False active[v] = active[next_v] = True for u in neighbor[v]: active[u] = True for u in neighbor[next_v]: active[u] = True ペナルティ重み更新(修正版) 隣接する都市のフラグ activeもTrueにする

Fast Local Searchの実装 • 局所探索法の近傍⾛査で改善解が得られた際に，更新が⽣じた辺 (u,v)のフラグactiveはTrueに設定する． 62 for v in
neighbor[u]: # evaluate difference (original) delta = eval_diff(tour, u, v, 'dist') if cur_obj + delta < obj - NUM_EPSILON: # update incumbent tour tour = cur_tour[:] change_tour(tour, u, v) obj = calc_length(tour) # evaluate difference (penalized) delta = eval_diff(tour, u, v, 'pdist') if delta < -NUM_EPSILON: # activate nodes active[u] = active[tour_next(u)] = True active[v] = active[tour_next(v)] = True # change current tour change_tour(tour, u, v) cur_obj = calc_plength(tour) return True return False 2-opt近傍の⾛査 Or-opt近傍⾛査のフラグ更新 def activate_node(tour, s, u, v): i = pos[u] head_p, tail_p = u, tour[(i+s-1) % num_node] prev_p, next_p = tour[(i-1) % num_node], tour[(i+s) % num_node] active[head_p] = active[tail_p] = True active[prev_p] = active[next_p] = True active[v] = active[tour_next(v)] = True 3-opt近傍⾛査のフラグ更新 def activate_node(tour, u, v, w, oper): active[u] = active[v] = active[w] = True if oper == 'type1' or oper == 'type2': active[tour_next(u)] = True active[tour_prev(v)] = True active[tour_next(w)] = True elif oper == 'type3': active[tour_next(u)] = True active[tour_prev(v)] = True active[tour_prev(w)] = True elif oper == 'type4': active[tour_next(u)] = True active[tour_next(v)] = True active[tour_next(w)] = True

重み付き局所探索法の実装 • 制約付き問題ではペナルティ関数法を導⼊しても，適切なペナルティ重みw[i]を決定することは容易ではない． • ペナルティ重みw[i]の更新と局所探索法を交互に繰り返し適⽤する． • ならば下式に従ってペナルティ重みを増加する．そうでなければ，⼀律に倍して減少する．
• 評価関数を変形するため局所最適解から脱出する効果も期待できる． 63 実⾏不可能領域実⾏可能領域有望な解が含まれる領域更新式1: 更新式2*: *更新式2では，ペナルティ重み更新後にとなる．

ペナルティ重みの⾃動調整 • ペナルティ重みw[i]の初期値はに設定する． • ペナルティ重み減少の直後にペナルティ重みを増加する⽅法もある． 64 # initialize current
solution sol = best_sol[:] # initialize penalty weight init_weight(wt) # weighting local search start_time = cur_time = time.time() while cur_time - start_time < time_limit: # local search algorithm local_search(best_sol,sol) # update penalty weight update_weight(sol,best_obj) cur_time = time.time() 重み付き局所探索法の⼿続き def init_weight(wt): for i in range(num_agent): wt[i] = -float('inf') for j in range(num_job): if cost[i][j] > wt[i]: wt[i] = cost[i][j] # update weighted penalty plt = 0.0 for i in range(num_agent): plt += wt[i] * max(0, used[i] – cap[i]) ペナルティ重みの初期化 def update_weight(sol, best_obj): if obj + plt > best_obj - NUM_EPSILON: # decrease penalty weight for i in range(num_agent): wt[i] = max(NUM_EPSILON, (1.0 - ALPHA) * wt[i]) else: # increase penalty weight max_plt = 0.0 for i in range(num_agent): if used[i] - cap[i] > max_plt: max_plt = used[i] - (gap.cap)[i] for i in range(num_agent): wt[i] *= 1.0 + DELTA * max(0, used[i] - cap[i]) / max_plt # update weighted penalty plt = 0.0 for i in range(num_agent): plt += wt[i] * max(0, used[i] – cap[i]) ペナルティ重みの更新更新式1

巡回セールスマン問題の数値実験* • TSPLIBから最適値が既知の問題例40題(都市数100〜783)． • Gurobi9.1.0(MTZ制約を⽤いた定式化*)，LocalSolver10.1と⽐較． • Gurobiは都市数150を超えるとNNと同程度の解すら難しい． • 隣接リストを導⼊すると計算速度は⼤幅に向上する． •
FLSの初回実⾏時の計算時間はLSとあまり変わらない． 66 *実⾏環境はMac mini(2018), CPU: Intel Core i7 (6core), Mem: 64GB *サンプルコードはPyPy3で実⾏．通常のPython3よりも約10倍ほど⾼速です． *部分巡回路除去制約を逐次的に追加する分枝カット法では1つも実⾏可能解が⾒つかりませんでした． *1000都市程度までなら最適解は巡回セールスマン問題専⽤のソルバーConcordeで求められます． Gurobi LocalSolver NN LS LS-NL (NL=20) LS-NL (NL=5) FLS (NL=20) 最適値からのギャップ* 75.72% 3.88% 24.58% 2.93% 3.58% 5.53% 4.00% 実⾏時間 60.00s 60.00s 0.01s 21.87s 0.63s 0.21s 0.40s 実⾏可能解 30/40 40/40 40/40 40/40 40/40 40/40 40/40 NN:最近近傍法，LS:局所探索法，LS-NL:隣接リストを⽤いた局所探索法， FLS:Fast Local Search

巡回セールスマン問題の数値実験 • メタヒューリスティクスの実装には全て隣接リストを導⼊*． • MLSは初期解の質が悪いのでLS1回当たりの実⾏時間が⻑い． • ILS-FLSはLSの実⾏回数は増加するが性能はあまり改善しない． • ILSにSAのアイデアを導⼊して初期解の種となる解を更新する修正版 ILS+は性能が良い．
• ILSの摂動に2-opt近傍操作を⽤いたILS'は性能があまり良くない． 67 LS MLS GRASP ILS ILS-FLS ILS+ ILS' 最適値からのギャップ 2.93% 2.27% 1.04% 0.34% 0.33% 0.16% 1.60% LSの実⾏回数 1 1630.9 2033.7 4491.1 18256.6 24971.9 76239.6 MLS: ランダム多スタート局所探索法，GRASP: GRASP ILS: 反復局所探索法，ILS+: 反復局所探索法(修正版)，ILS': 反復局所探索法(2-opt近傍操作)

巡回セールスマン問題の数値実験 • MAはFLSを導⼊すると若⼲性能が改善する． • SAは温度設定により性能が変わる．初期温度は低めの⽅が良い． • JohnsonによるSAの実装(SA')やParallel tempering(PT,レプリカ法) などの実装はあまり性能が良くない． •
TSは単純な実装では性能が良くない． • GLSにFLSを導⼊するとLSの実⾏回数が増加し性能も改善する． 68 *隣接リストの⻑さはSAのみNL=20，それ以外はNL=5． *SA'の終了条件は実⾏時間ではない．平均実⾏時間は2.05秒． MA: Memetic algorithm， SA: アニーリング法，SA': アニーリング法(Johnson)，PT: Parallel tempering TS1: タブー探索法(規則1)，TS2: タブー探索法(規則2) GLS: 誘導局所探索法 MA-FLS SA SA' PT TS1 TS2 GLS GLS-FLS 最適値からのギャップ 0.45% 0.30% 2.64% 1.49% 3.00% 2.53% 1.30% 0.61% LSの実⾏回数 5929.3 1294.9 17072.3

数値実験 • rat783に対してアルゴリズムを適⽤した際の，時間経過にともなう現在の解(⻘)と最良解(橙)の⽬的関数値の推移を描画*． 69 ILS+ ILS-FLS GRASP ILS' MLS
*縦軸はlog scaleにしています．

数値実験 • rat783に対してアルゴリズムを適⽤した際の，時間経過にともなう現在の解(⻘)と最良解(橙)の⽬的関数値の推移を描画． • 全アルゴリズムの最良解の⽬的関数値の推移を⽐較． SA GLS-FLS TS1 TS2
MA-FLS 70

⼀般化割当問題の数値実験 • ベンチマーク問題集TypeA〜Eの57題(エージェント数m=5〜80，仕事数n=100〜1600) • Gurobi9.1.0，Python-MIP(CBC)と⽐較*． • 緩和問題から得られる下界が⾮常に良く，GurobiでもPython-MIP (CBC)でも短時間で最適解が求まることが多い． •
WLSは両⽅の更新式で良い実⾏可能解が得られる(最適解は難しい)． 71 Gurobi Python- MIP(CBC) WLS (更新式1) WLS+ (更新式1) WLS+ (更新式2) 最良値からのギャップ 0.04% 0.05% 0.43% 0.38% 0.40% LSの実⾏回数 14222.8 44549.5 39448.5 WLS:重み付き局所探索法(更新式1,2) WLS+:重み付き局所探索法(修正版*) *Gurobiの実⾏時間は60秒，それ以外はn<=400は60秒，n=900は300秒，n=1600は600秒． *少なくとも⼀⽅の仕事が容量制約を違反するエージェントに含まれるスワップ近傍のみに⾛査を制限．

まとめ • 局所探索法の設計と実装 ü 巡回セールスマン問題(TSP) ü TSP: 最近傍法(NN), 2-opt法，Or-opt法，3-opt法 ü
TSP: 局所探索法の⾼速化 ü ⼀般化割当問題(GAP) ü GAP: 制約付き問題に対するペナルティ関数法 ü GAP: シフト近傍⾛査，スワップ近傍⾛査 • メタヒューリスティクスの設計と実装 ü TSP: ランダム多スタート局所探索法(MLS,GRASP) ü TSP: 反復局所探索法(ILS) ü TSP: 遺伝的アルゴリズム(GA,MA) ü TSP: アニーリング法(SA) ü TSP: タブー探索法(TS) ü GAP: 誘導局所探索法(GLS)，Fast Local Search(FLS) ü GAP: 制約付き問題に対する重み付き局所探索法(WLS) • 数値実験 72 現実問題の解決にメタヒューリスティクスを活⽤して下さい︕

参考⽂献 • 柳浦睦憲，茨⽊俊秀，組合せ最適化―メタ戦略を中⼼として―，朝倉書店，2001． • 久保幹雄，J.P.ペドロソ，メタヒューリスティクスの数理，共⽴出版，2012． • 梅⾕俊治，メタヒューリスティクス事始め︓まずは局所探索法から，オペレーションズ・リサーチ，58 (2013)，689-694． •
今堀慎治，柳浦睦憲，概説メタヒューリスティクス，オペレーションズ・リサーチ， 58(2013)，695-702． • 橋本英樹，野々部宏司，⼊⾨タブー探索法，オペレーションズ・リサーチ，58(2013), 703-707. • 永⽥裕⼀，多点探索型アルゴリズムの基礎と最前線，オペレーションズ・リサーチ， 58(2013), 708-715. • 野々部宏司，メタヒューリスティックアルゴリズム実装の⼿順と留意点，オペレーションズ・リサーチ，58(2013)，716-722． • 久保幹雄，巡回セールスマン問題への招待II，オペレーションズ・リサーチ，39(1994)， 91-96. • 久保幹雄，巡回セールスマン問題への招待III，オペレーションズ・リサーチ，39(1994)， 156-162. • 岩⽥陽⼀，Introduction to Heuristics Contest解説，AtCoder, 2020. 73

参考⽂献 • E.Aarts, K.Lenstra (eds.), Local Search in combinatorial optimization,
Princeton University Press, 1997. • M.Gendreau, J.-Y.Potvin(eds.), Handbook of Metaheuristics (3rd ed.), Springer, 2019. • D.L.Applegate, R.E.Bixby, V.Chvatal and W.J.Cook, The Traveling Salesman Problem: A Computational Study, Princeton University Press, 2006. • E.L.Lawler, J.K.Lenstra, A.H.G.Rinnooy Kan and D.B.Shmoys (eds.), The Traveling Salesman Problem, Wiley, 1985. • TSPサンプルプログラム https://github.com/shunji-umetani/tsp-solver • GAPサンプルプログラム https://github.com/shunji-umetani/gap-solver 74

Introduction to Design and Implementation of Metaheuristics

Introduction to Design and Implementation of Metaheuristics

More Decks by Shunji Umetani

Other Decks in Technology

Featured

Transcript