SA
代码直接用就行，成功率极高
import random
import numpy as np
import matplotlib.pyplot as plt
# randomly generate the map with constraint of [-100, 100]
def gen_cities(city_num, random_state=True):
if random_state:
cities = (np.random.uniform(0, 2, (city_num, 2))-1)*200
np.savetxt('cities.txt', cities)
else:
cities = np.loadtxt('cities.txt')
return cities
# get distance matrix
def get_distance_mateix(cities):
distances = np.zeros((len(cities), len(cities)))
for i in range(len(cities)):
for j in range(len(cities)):
distances[i][j] = np.sqrt((cities[i][0]-cities[j][0])**2 + (cities[i][1]-cities[j][1])**2)
return distances
# plot cities on a graph
def plot_map(cities, path):
plt.figure(figsize=(10,10), dpi=100)
plt.scatter(cities[:,0], cities[:,1], s=150)
for i in range(len(cities)-1):
plt.plot([cities[path[i]][0], cities[path[i+1]][0]], [cities[path[i]][1], cities[path[i+1]][1]], c='r')
plt.plot([cities[path[0]][0], cities[path[len(cities)-1]][0]], [cities[path[0]][1], cities[path[len(cities)-1]][1]], c='r')
plt.savefig('cities.png')
# find fittest path called every generation
def find_fittest(routes_length, pop_size):
key = 10000000000
fittest = 0
for i in range(pop_size):
if routes_length[i] < key:
key = routes_length[i]
fittest = i
return fittest
def SA(cities, distances):
num = cities.shape[0]
"""***********退火系数**********"""
delta = 0.99
t = 100
tk = 1
solution_new = np.arange(num)
solution_current = solution_new.copy()
value_current = 99999
solution_best = solution_new.copy()
value_best = 99999
peak_length = []
peak = []
hillside = []
hillside_length = []
while t > tk:
for i in np.arange(1000):
if np.random.rand() > 0.5:
while True:
loc1 = int(np.ceil(np.random.rand() * (num - 1)))
loc2 = int(np.ceil(np.random.rand() * (num - 1)))
if loc1 != loc2:
break
solution_new[loc1], solution_new[loc2] = solution_new[loc2], solution_new[loc1]
else:
while True:
loc1 = int(np.ceil(np.random.rand() * (num - 1)))
loc2 = int(np.ceil(np.random.rand() * (num - 1)))
loc3 = int(np.ceil(np.random.rand() * (num - 1)))
if (loc1 != loc2) & (loc2 != loc3) & (loc1 != loc3):
break
if loc1 > loc2:
loc1, loc2 = loc2, loc1
if loc2 > loc3:
loc2, loc3 = loc3, loc2
if loc1 > loc2:
loc1, loc2 = loc2, loc1
loc1_loc2 = solution_new[loc1:loc2].copy()
solution_new[loc1:loc3 - loc2 + 1 + loc1] = solution_new[loc2:loc3 + 1].copy()
solution_new[loc3 - loc2 + 1 + loc1:loc3 + 1] = loc1_loc2.copy()
value_new = 0
for j in range(num - 1):
value_new += distances[solution_new[j]][solution_new[j + 1]]
value_new += distances[solution_new[0]][solution_new[num - 1]]
if value_new < value_current:
value_current = value_new
solution_current = solution_new.copy()
hillside_length.append(value_new)
hillside.append(solution_new)
if value_new < value_best:
value_best = value_new
solution_best = solution_new.copy()
peak_length.append(value_best)
peak.append(solution_best)
else:
if np.random.rand() < np.exp(-(value_new - value_current) / t):
value_current = value_new
solution_current = solution_new.copy()
hillside_length.append(value_new)
hillside.append(solution_new)
else:
solution_new = solution_current.copy()
t = delta * t
return peak_length, peak, hillside_length, hillside
def main():
cities = gen_cities(city_num=50, random_state=True)
distances = get_distance_mateix(cities)
best_lengths, beat_paths, hillside_length, hillside = SA(cities, distances)
index = find_fittest(best_lengths, len(best_lengths))
best_path = beat_paths[index]
print("Best path is:", best_path, "with length", best_lengths[index], '\n')
for i in range(len(hillside_length)):
print('length=', hillside_length[i], '\n')
for i in range(len(best_lengths)):
print('peak_length=', best_lengths[i], '\n')
plot_map(cities, beat_paths[index])
if __name__=="__main__":
main()
plt.show()
上面分别为50点和100点的运行图，都是第一次生成就直接找到全局最优解，后面经过几十次测试大概感觉全局最优解的概率在百分之80到90之间，点少的情况太简单就不展示了，可以接近百分之百的全局最优解（刚好围成一圈并且没有线段交叉）
我也在代码中添加了记录爬山过程的部分和每一个山峰（极大值点），感觉看着结果时增时降也是很有趣的，就好像自己在爬山一样（不过爬山过程太长了，爬的实在是太累了）
我发现在二维的低维空间下退火算法相对与遗传算法有碾压性的优势：
- 扰动方式简单，退火算法中我仅仅采用二点交换与三角交换进行扰动，不需要有序交叉或者单点变异，极大程度地降低了结果地不稳定性。
- 资源要求少，退火算法的初始种群仅有1个，但是仍然能够以出色的概率找到全局最优解，本质上是节省了非常多的时间在搜索与筛选上面的
- 效果好，遗传算法在8点以前可以做到大于二分之一的概率找到全局最优解，但是超过9点之后出现全局最优解的概率几乎是零（我认为是路径变复杂后种群数量不足导致的）
SA+GA
后来老师说可以用遗传算法优化退火算法的初始种群，于是我就做了一个缝合版
import random
import numpy as np
import matplotlib.pyplot as plt
# randomly generate the map with constraint of [-100, 100]
def gen_cities(city_num, random_state=True):
if random_state:
cities = (np.random.rand(city_num,2)-0.5)*200
np.savetxt('cities.txt', cities)
else:
cities = np.loadtxt('cities.txt')
return cities
# get distance matrix
def get_distance_mateix(cities):
distances = np.zeros((len(cities), len(cities)))
for i in range(len(cities)):
for j in range(len(cities)):
distances[i][j] = np.sqrt((cities[i][0]-cities[j][0])**2 + (cities[i][1]-cities[j][1])**2)
return distances
# plot cities on a graph
def plot_map(cities, path):
plt.figure(figsize=(10,10), dpi=100)
plt.scatter(cities[:,0], cities[:,1], s=150)
for i in range(len(cities)-1):
plt.plot([cities[path[i]][0], cities[path[i+1]][0]], [cities[path[i]][1], cities[path[i+1]][1]], c='r')
plt.plot([cities[path[0]][0], cities[path[len(cities)-1]][0]], [cities[path[0]][1], cities[path[len(cities)-1]][1]], c='r')
plt.savefig('cities.png')
# calculates distance between 2 cities
def calc_distance(distances, city1, city2):
return distances[city1][city2] # ord('A')=65
# creates a random route
def create_route(cities):
shuffled = random.sample(range(len(cities)), len(cities))
return shuffled
# calculates length of an route
def calc_route_length(population, distances, cities, routes_length, fitness):
pop_size = len(population)
for i in range(pop_size):
route_l = 0
for j in range(1, len(cities)):
route_l = route_l + calc_distance(distances, population[i][j - 1], population[i][j])
# route_l = route_l + calc_distance(population[i][len(cities)-1], population[i][1]) calculate distance from last to first
routes_length[i] = route_l
fitness[i] = 1 / routes_length[i]
return routes_length, fitness
# pop_size = len(fitness)
# for i in range(pop_size):
# route_l = 0
# for j in range(0, len(cities)-1):
# route_l += distances[population[i][j]][population[i][j + 1]]
# route_l += distances[population[i][0]][population[i][len(cities)-1]]
# routes_length[i] = route_l
# fitness[i] = 1 / route_l
#
# print(routes_length)
#
# return routes_length, fitness
# creates starting population
def create_population(pop_size, cities):
population = []
for i in range(pop_size):
population.append(create_route(cities))
return population
def swap_mutation(ind, cities, population):
a, b = random.sample(range(len(cities)), 2)
population[ind][a], population[ind][b] = population[ind][b], population[ind][a]
unique_cities = set()
for city in population[ind]:
if city in unique_cities:
population[ind] = create_route(cities)
break
unique_cities.add(city)
return population
# print("Mutated path: ", population[ind])
"""*********这玩意是不是叫partially_mapped_crossover**********"""
def partially_matched_crossover(ind1, ind2, cities):
size = len(ind1)
site1 = random.randint(0, size - 1)
site2 = random.randint(site1 + 1, size)
child1 = [-1] * size
child2 = [-1] * size
child1[site1:site2] = ind1[site1:site2]
child2[site1:site2] = ind2[site1:site2]
for i in range(size):
if i < site1 or i >= site2:
for j in range(size):
if ind2[j] not in child1:
child1[i] = ind2[j]
break
for j in range(size):
if ind1[j] not in child2:
child2[i] = ind1[j]
break
return child1, child2
"""********适应度缩放选择 减缓了轮盘选择的收敛 增加多样性*********"""
# function that picks a parent Fitness Proportionate Selection
def roulette_wheel_selection(pop_size, fitness):
maximum = fitness[pop_size-1]
minimum = fitness[0]
a = (maximum - minimum) / (maximum - minimum + pop_size)
b = maximum - a * maximum
sum_fitness = sum(fitness)
cumulative_probability = 0
for i in range(pop_size):
fitness[i] = (a * fitness[i] + b) / sum_fitness
cumulative_probability += fitness[i]
selected = random.uniform(0, 1)
index = pop_size - 1
for i in range(pop_size - 1):
if fitness[i] < selected :
fitness[i + 1] += fitness[i]
else:
index = i
break
return index
"""****************退火********************"""
def tuihuo(route, distances):
delta = 0.99 #原来为0.99
t = 100
tk = 1
solution_new = route
solution_current = solution_new.copy()
value_current = 99999
solution_best = solution_new.copy()
value_best = 99999
"""记录爬山过程
peak_length = []
peak = []
hillside = []
hillside_length = []
"""
while t > tk:
for i in np.arange(1000):
if np.random.rand() > 0.5:
while True:
loc1 = int(np.ceil(np.random.rand() * (num - 1)))
loc2 = int(np.ceil(np.random.rand() * (num - 1)))
if loc1 != loc2:
break
solution_new[loc1], solution_new[loc2] = solution_new[loc2], solution_new[loc1]
else:
while True:
loc1 = int(np.ceil(np.random.rand() * (num - 1)))
loc2 = int(np.ceil(np.random.rand() * (num - 1)))
loc3 = int(np.ceil(np.random.rand() * (num - 1)))
if (loc1 != loc2) & (loc2 != loc3) & (loc1 != loc3):
break
if loc1 > loc2:
loc1, loc2 = loc2, loc1
if loc2 > loc3:
loc2, loc3 = loc3, loc2
if loc1 > loc2:
loc1, loc2 = loc2, loc1
loc1_loc2 = solution_new[loc1:loc2].copy()
solution_new[loc1:loc3 - loc2 + 1 + loc1] = solution_new[loc2:loc3 + 1].copy()
solution_new[loc3 - loc2 + 1 + loc1:loc3 + 1] = loc1_loc2.copy()
value_new = 0
for j in range(num - 1):
value_new += distances[solution_new[j]][solution_new[j + 1]]
value_new += distances[solution_new[0]][solution_new[num - 1]]
if value_new < value_current:
value_current = value_new
solution_current = solution_new.copy()
# hillside_length.append(value_new)
# hillside.append(solution_new)
if value_new < value_best:
value_best = value_new
solution_best = solution_new.copy()
# peak_length.append(value_best)
# peak.append(solution_best)
else:
if np.random.rand() < np.exp(-(value_new - value_current) / t):
value_current = value_new
solution_current = solution_new.copy()
# hillside_length.append(value_new)
# hillside.append(solution_new)
else:
solution_new = solution_current.copy()
t = delta * t
return solution_best
# find fittest path called every generation
def find_fittest(routes_length, pop_size):
key = 10000000000
fittest = 0
for i in range(pop_size):
if routes_length[i] < key:
key = routes_length[i]
fittest = i
return fittest
# sorts parallely the lists
#def sort_alongside(routes_length, population):
# routes_length, population = (list(i) for i in zip(*sorted(zip(routes_length, population))))
# initialize algorithm
def main():
pop_size = 30
mutate_prob = 0.1
n_generations = 30
routes_length = [0]*pop_size
fitness = [0]*pop_size
best_path = 10000000000
cities = gen_cities(city_num=8,random_state=True)
best_route = [0] * cities.shape[0]
distances = get_distance_mateix(cities)
population = create_population(pop_size=pop_size, cities=cities)
# print("Population initialization:", "\n", population)
routes_length, fitness = calc_route_length(population, distances, cities, routes_length, fitness)
# print("Population's paths length:", "\n", routes_length)
for j in range(n_generations):
for i in range(0, pop_size, 2):
# pick parents for crossover
parent1 = roulette_wheel_selection(pop_size, fitness)
parent2 = roulette_wheel_selection(pop_size, fitness)
# always pick different parents (not necessary)
while True:
if parent1 == parent2:
parent2 = roulette_wheel_selection(pop_size, fitness)
else:
break
# update population
population[i], population[i + 1] = partially_matched_crossover(population[parent1], population[parent2], cities)
# calculate lengths for updated generation
routes_length, fitness = calc_route_length(population, distances, cities, routes_length, fitness)
# pick the paths for mutation based on a probability
for i in range(pop_size):
rand = random.uniform(0, 1)
if rand < mutate_prob:
population = swap_mutation(i, cities, population)
if i % 5 == 0:
for i in range(pop_size):
population[i] = tuihuo(population[i], distances)
# calculate lengths after mutation
routes_length, fitness = calc_route_length(population, distances, cities, routes_length, fitness)
# find best path overall
if routes_length[find_fittest(routes_length, pop_size)] < best_path:
index = find_fittest(routes_length, pop_size)
best_path = routes_length[index]
best_route = population[index]
# print("Best route of generation", j+1, ": ", population[find_fittest()], "\n" "Route length: ",
# routes_length[find_fittest()])
# print("Population of generation", j+1, ": \n", population)
# print("Routes lengths:", routes_length, "\n")
index = find_fittest(routes_length, pop_size)
print("Best path is:", best_route, "with length", best_path)
plot_map(cities, population[index])
if __name__=="__main__":
main()
plt.show()不幸的是，导致产出结果并不是很好，代码部分应该还有很大不足的地方，所以就不展示图片了。
标签：length,solution,loc2,算法,模拟退火,new,SA,cities,population From： https://blog.csdn.net/m0_61546651/article/details/136969212