算法流程：
實現：
base.py
from abc import ABCMeta, abstractmethod import types class SkoBase(metaclass=ABCMeta): def register(self, operator_name, operator, *args, **kwargs): ''' regeister udf to the class :param operator_name: string :param operator: a function, operator itself :param args: arg of operator :param kwargs: kwargs of operator :return: ''' def operator_wapper(*wrapper_args): return operator(*(wrapper_args + args), **kwargs) setattr(self, operator_name, types.MethodType(operator_wapper, self)) return self # 空函數 class Problem(object): pass
sko目錄下的operaters目錄 mutation.py 用於TSP問題的狀態空間轉移函數的設計
import numpy as np def mutation(self): ''' mutation of 0/1 type chromosome faster than `self.Chrom = (mask + self.Chrom) % 2` :param self: :return: ''' # mask = (np.random.rand(self.size_pop, self.len_chrom) < self.prob_mut) self.Chrom ^= mask return self.Chrom def mutation_TSP_1(self): ''' every gene in every chromosome mutate :param self: :return: ''' for i in range(self.size_pop): for j in range(self.n_dim): if np.random.rand() < self.prob_mut: n = np.random.randint(0, self.len_chrom, 1) self.Chrom[i, j], self.Chrom[i, n] = self.Chrom[i, n], self.Chrom[i, j] return self.Chrom # 交換兩個點 def swap(individual): n1, n2 = np.random.randint(0, individual.shape[0] - 1, 2) if n1 >= n2: n1, n2 = n2, n1 + 1 individual[n1], individual[n2] = individual[n2], individual[n1] return individual # 逆轉序列中指定的一段 def reverse(individual): ''' Reverse n1 to n2 Also called `2-Opt`: removes two random edges, reconnecting them so they cross Karan Bhatia, "Genetic Algorithms and the Traveling Salesman Problem", 1994 https://pdfs.semanticscholar.org/c5dd/3d8e97202f07f2e337a791c3bf81cd0bbb13.pdf ''' n1, n2 = np.random.randint(0, individual.shape[0] - 1, 2) if n1 >= n2: n1, n2 = n2, n1 + 1 individual[n1:n2] = individual[n1:n2][::-1] return individual # 分段交叉 def transpose(individual): # randomly generate n1 < n2 < n3. Notice: not equal n1, n2, n3 = sorted(np.random.randint(0, individual.shape[0] - 2, 3)) n2 += 1 n3 += 2 slice1, slice2, slice3, slice4 = individual[0:n1], individual[n1:n2], individual[n2:n3 + 1], individual[n3 + 1:] individual = np.concatenate([slice1, slice3, slice2, slice4]) return individual def mutation_reverse(self): ''' Reverse :param self: :return: ''' for i in range(self.size_pop): if np.random.rand() < self.prob_mut: self.Chrom[i] = reverse(self.Chrom[i]) return self.Chrom def mutation_swap(self): for i in range(self.size_pop): if np.random.rand() < self.prob_mut: self.Chrom[i] = swap(self.Chrom[i]) return self.Chrom
SA.py 定義了眾多的模擬退火方案
#!/usr/bin/env python3 # -*- coding: utf-8 -*- # @Time : 2019/8/17 # @Author : github.com/guofei9987 import numpy as np from base import SkoBase from sko.operators import mutation class SimulatedAnnealingBase(SkoBase): """ DO SA(Simulated Annealing) Parameters ---------------- func : function The func you want to do optimal n_dim : int number of variables of func x0 : array, shape is n_dim initial solution T_max :float initial temperature T_min : float end temperature L : int num of iteration under every temperature（Long of Chain） Attributes ---------------------- Examples ------------- See https://github.com/guofei9987/scikit-opt/blob/master/examples/demo_sa.py """ #L是單次迭代中馬爾科夫鏈的長度 def __init__(self, func, x0, T_max=100, T_min=1e-7, L=300, max_stay_counter=150, **kwargs): assert T_max > T_min > 0, 'T_max > T_min > 0' self.func = func self.T_max = T_max # initial temperature self.T_min = T_min # end temperature self.L = int(L) # num of iteration under every temperature（also called Long of Chain） self.max_stay_counter = max_stay_counter # stop if best_y stay unchanged over max_stay_counter times self.n_dims = len(x0) self.best_x = np.array(x0) # initial solution self.best_y = self.func(self.best_x) self.T = self.T_max self.iter_cycle = 0 self.best_y_history = [self.best_y] self.best_x_history = [self.best_x] def get_new_x(self, x): u = np.random.uniform(-1, 1, size=self.n_dims)#產生[-1,1]之間的隨機數，維數指定 x_new = x + 20 * np.sign(u) * self.T * ((1 + 1.0 / self.T) ** np.abs(u) - 1.0) return x_new def cool_down(self): self.T = self.T * 0.7 def isclose(self, a, b, rel_tol=1e-09, abs_tol=1e-30): return abs(a - b) <= max(rel_tol * max(abs(a), abs(b)), abs_tol) def run(self): x_current, y_current = self.best_x, self.best_y stay_counter = 0 while True: for i in range(self.L): x_new = self.get_new_x(x_current) y_new = self.func(x_new) # Metropolis df = y_new - y_current if df < 0 or np.exp(-df / self.T) > np.random.rand(): x_current, y_current = x_new, y_new if y_new < self.best_y: self.best_x, self.best_y = x_new, y_new self.iter_cycle += 1 self.cool_down() self.best_y_history.append(self.best_y) self.best_x_history.append(self.best_x) # if best_y stay for max_stay_counter times, stop iteration # 最后一次的best_y值和倒數第二次的best_y值一樣，表示停留在平坦的狀態 if self.isclose(self.best_y_history[-1], self.best_y_history[-2]): stay_counter += 1 else: stay_counter = 0 if self.T < self.T_min: stop_code = 'Cooled to final temperature' break if stay_counter > self.max_stay_counter: stop_code = 'Stay unchanged in the last {stay_counter} iterations'.format(stay_counter=stay_counter) break return self.best_x, self.best_y fit = run class SAFast(SimulatedAnnealingBase): ''' u ~ Uniform(0, 1, size = d) y = sgn(u - 0.5) * T * ((1 + 1/T)**abs(2*u - 1) - 1.0) xc = y * (upper - lower) x_new = x_old + xc c = n * exp(-n * quench) T_new = T0 * exp(-c * k**quench) ''' def __init__(self, func, x0, T_max=100, T_min=1e-7, L=300, max_stay_counter=150, **kwargs): super().__init__(func, x0, T_max, T_min, L, max_stay_counter, **kwargs) self.m, self.n, self.quench = kwargs.get('m', 1), kwargs.get('n', 1), kwargs.get('quench', 1) self.lower, self.upper = kwargs.get('lower', -10), kwargs.get('upper', 10) self.c = self.m * np.exp(-self.n * self.quench) def get_new_x(self, x): r = np.random.uniform(-1, 1, size=self.n_dims) xc = np.sign(r) * self.T * ((1 + 1.0 / self.T) ** np.abs(r) - 1.0) x_new = x + xc * (self.upper - self.lower) return x_new def cool_down(self): self.T = self.T_max * np.exp(-self.c * self.iter_cycle ** self.quench) class SABoltzmann(SimulatedAnnealingBase): ''' std = minimum(sqrt(T) * ones(d), (upper - lower) / (3*learn_rate)) y ~ Normal(0, std, size = d) x_new = x_old + learn_rate * y T_new = T0 / log(1 + k) ''' def __init__(self, func, x0, T_max=100, T_min=1e-7, L=300, max_stay_counter=150, **kwargs): super().__init__(func, x0, T_max, T_min, L, max_stay_counter, **kwargs) self.lower, self.upper = kwargs.get('lower', -10), kwargs.get('upper', 10) self.learn_rate = kwargs.get('learn_rate', 0.5) def get_new_x(self, x): std = min(np.sqrt(self.T), (self.upper - self.lower) / 3.0 / self.learn_rate) * np.ones(self.n_dims) xc = np.random.normal(0, 1.0, size=self.n_dims) x_new = x + xc * std * self.learn_rate return x_new def cool_down(self): self.T = self.T_max / np.log(self.iter_cycle + 1.0) class SACauchy(SimulatedAnnealingBase): ''' u ~ Uniform(-pi/2, pi/2, size=d) xc = learn_rate * T * tan(u) x_new = x_old + xc T_new = T0 / (1 + k) ''' def __init__(self, func, x0, T_max=100, T_min=1e-7, L=300, max_stay_counter=150, **kwargs): super().__init__(func, x0, T_max, T_min, L, max_stay_counter, **kwargs) self.learn_rate = kwargs.get('learn_rate', 0.5) def get_new_x(self, x): u = np.random.uniform(-np.pi / 2, np.pi / 2, size=self.n_dims) xc = self.learn_rate * self.T * np.tan(u) x_new = x + xc return x_new def cool_down(self): self.T = self.T_max / (1 + self.iter_cycle) # SA_fast is the default SA = SAFast class SA_TSP(SimulatedAnnealingBase): def cool_down(self): self.T = self.T_max / (1 + np.log(1 + self.iter_cycle)) def get_new_x(self, x): x_new = x.copy() new_x_strategy = np.random.randint(3) if new_x_strategy == 0: x_new = mutation.swap(x_new) elif new_x_strategy == 1: x_new = mutation.reverse(x_new) elif new_x_strategy == 2: x_new = mutation.transpose(x_new) return x_new
demo_sa.py模擬退火的函數極值搜索案例，使用了一個著名的震盪函數，只能搜索到局部最優解，因為這個函數的極小值群體過於龐大
import numpy as np # demo_func = lambda x: x[0] ** 2 + (x[1] - 0.05) ** 2 + x[2] ** 2 # 震盪函數，有無窮多個極小值點 demo_func = lambda x: 0.5+(np.square(np.sin(np.sqrt(x[0]**2+x[1]**2)))-0.5)/np.square(1+0.001*(x[0]**2+x[1]**2)) # %% Do SA from sko.SA import SA sa = SA(func=demo_func, x0=[1,1], T_max=1, T_min=1e-9, L=300, max_stay_counter=150) best_x, best_y = sa.run() print('best_x:', best_x, 'best_y', best_y) # %% Plot the result import matplotlib.pyplot as plt import pandas as pd plt.plot(pd.DataFrame(sa.best_y_history)) plt.show() # %% from sko.SA import SAFast sa_fast = SAFast(func=demo_func, x0=[1, 1], T_max=1, T_min=1e-9, q=0.99, L=300, max_stay_counter=150) sa_fast.run() print('Fast Simulated Annealing: best_x is ', sa_fast.best_x, 'best_y is ', sa_fast.best_y) # %% from sko.SA import SABoltzmann sa_boltzmann = SABoltzmann(func=demo_func, x0=[1, 1], T_max=1, T_min=1e-9, q=0.99, L=300, max_stay_counter=150) sa_boltzmann.run() print('Boltzmann Simulated Annealing: best_x is ', sa_boltzmann.best_x, 'best_y is ', sa_fast.best_y) # %% from sko.SA import SACauchy sa_cauchy = SACauchy(func=demo_func, x0=[1, 1], T_max=1, T_min=1e-9, q=0.99, L=300, max_stay_counter=150) sa_cauchy.run() print('Cauchy Simulated Annealing: best_x is ', sa_cauchy.best_x, 'best_y is ', sa_cauchy.best_y)
得到的函數值下降曲線為
TSP問題的優化以及動態迭代過程展示
import numpy as np from scipy import spatial import matplotlib.pyplot as plt import sys file_name = sys.argv[1] if len(sys.argv) > 1 else 'data/nctu.csv' points_coordinate = np.loadtxt(file_name, delimiter=',') num_points = points_coordinate.shape[0] distance_matrix = spatial.distance.cdist(points_coordinate, points_coordinate, metric='euclidean') distance_matrix = distance_matrix * 111000 # 1 degree of lat/lon ~ = 111000m def cal_total_distance(routine): '''The objective function. input routine, return total distance. cal_total_distance(np.arange(num_points)) ''' num_points, = routine.shape return sum([distance_matrix[routine[i % num_points], routine[(i + 1) % num_points]] for i in range(num_points)]) # %% from sko.SA import SA_TSP sa_tsp = SA_TSP(func=cal_total_distance, x0=range(num_points), T_max=100, T_min=1, L=10 * num_points) best_points, best_distance = sa_tsp.run() print(best_points, best_distance, cal_total_distance(best_points)) # %% Plot the best routine from matplotlib.ticker import FormatStrFormatter fig, ax = plt.subplots(1, 2) best_points_ = np.concatenate([best_points, [best_points[0]]]) best_points_coordinate = points_coordinate[best_points_, :] ax[0].plot(sa_tsp.best_y_history) ax[0].set_xlabel("Iteration") ax[0].set_ylabel("Distance") ax[1].plot(best_points_coordinate[:, 0], best_points_coordinate[:, 1], marker='o', markerfacecolor='b', color='c', linestyle='-') ax[1].xaxis.set_major_formatter(FormatStrFormatter('%.3f')) ax[1].yaxis.set_major_formatter(FormatStrFormatter('%.3f')) ax[1].set_xlabel("Longitude") ax[1].set_ylabel("Latitude") plt.show() # %% Plot the animation import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation best_x_history = sa_tsp.best_x_history fig2, ax2 = plt.subplots(1, 1) ax2.set_title('title', loc='center') line = ax2.plot(points_coordinate[:, 0], points_coordinate[:, 1], marker='o', markerfacecolor='b', color='c', linestyle='-') ax2.xaxis.set_major_formatter(FormatStrFormatter('%.3f')) ax2.yaxis.set_major_formatter(FormatStrFormatter('%.3f')) ax2.set_xlabel("Longitude") ax2.set_ylabel("Latitude") plt.ion() p = plt.show() def update_scatter(frame): ax2.set_title('iter = ' + str(frame)) points = best_x_history[frame] points = np.concatenate([points, [points[0]]]) point_coordinate = points_coordinate[points, :] plt.setp(line, 'xdata', point_coordinate[:, 0], 'ydata', point_coordinate[:, 1]) return line ani = FuncAnimation(fig2, update_scatter, blit=True, interval=25, frames=len(best_x_history)) plt.show() # ani.save('sa_tsp.gif', writer='pillow')
data/nctu.csv 文件
24.783003,120.997474 24.785737,120.996905 24.784604,120.998911 24.785261,120.998698 24.785266,120.996754 24.784001,120.998107 24.783266,120.99845 24.785178,120.994979 24.784736,120.995445 24.785626,120.995474 24.784542,120.995755 24.788155,120.99547 24.788339,120.995857 24.78886,120.995259 24.788192,120.995671 24.786551,120.995681 24.787233,120.995651 24.786942,120.995354 24.788947,120.994135 24.789596,120.995856 24.78969,120.995223 24.789652,120.994996 24.788933,120.995836 24.789182,120.997927 24.786731,121.001706 24.787541,120.998777 24.786444,120.997757 24.788442,120.999333 24.787059,120.996254 24.788719,120.997473 24.787729,121.000022 24.786509,120.997 24.787374,121.000994 24.787748,120.997273 24.789707,120.99633 24.785103,121.000918 24.787056,120.997282 24.785563,120.998591 24.79033,120.997038 24.785684,120.999245 24.788844,120.999945 24.789328,120.997568 24.786955,120.999126 24.78872,120.999253 24.78777,120.998151 24.788856,120.999068 24.787943,120.996833 24.788764,120.996599 24.784636,120.999428 24.784805,120.99954 24.784957,120.997073 24.783229,120.997842 24.785159,120.998036 24.787014,120.998277 24.787786,121.000889 24.788249,121.000433 24.78445,120.997293 24.786209,120.999259 24.788648,121.000958 24.787658,120.999493 24.784886,121.000451 24.78455,120.997889 24.786806,121.000725 24.786624,120.999088 24.784222,120.996917 24.787919,120.996822 24.788746,121.001269 24.788667,120.998652 24.783467,120.997571 24.788086,121.001596 24.788868,120.998306 24.786342,121.000484 24.785605,120.997364 24.787122,120.998925 24.785791,120.998019 24.786612,120.996584 24.786708,120.999975 24.78581,121.000076
gif迭代過程展示，gif太大無法展示，可以通過上面的Python源碼運行后進行保存即可，根據每次迭代，動態更新的，可以看出TSP問題求解的過程
