算法流程：
实现：
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‘breakif stay_counter > self.max_stay_counter:
stop_code = ‘Stay unchanged in the last {stay_counter} iterations‘.format(stay_counter=stay_counter)
breakreturn 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 SAfrom 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 resultimport 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 ~ = 111000mdef 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 routinefrom 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 animationimport 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问题求解的过程
原文：https://www.cnblogs.com/randy-lo/p/13286020.html
【说明】：本文章由站长整理发布，文章内容不代表本站观点，如文中有侵权行为，请与本站客服联系(QQ:254677821)！