-
Notifications
You must be signed in to change notification settings - Fork 1
/
SA.py
164 lines (136 loc) · 5.13 KB
/
SA.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
import xml.etree.ElementTree as ET
import random
import math
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
number_of_dest = 10
initial_temp = 1000
minimum_temp = 0
cooling_rate = 1
num_iterations = 100
def main():
map_name = 'map_uwaterloo.osm.xml'
osm = ET.parse(map_name)
root = osm.getroot()
# get the boundary of the map
bounds = root.find('bounds').attrib
x_max = lon_to_x(bounds['maxlon'], bounds['minlon'])
y_max = lat_to_y(bounds['maxlat'], bounds['minlat'])
z_max = 100
nodes = generate_nodes(x_max, y_max, z_max, number_of_dest + 1)
home = nodes[0]
destinations = nodes[1:]
# print(f'home: {home}')
# print(f'destinations: {destinations}')
neighbourhood = build_neighbourhood(destinations)
cost_dict = build_cost_dict(nodes)
run_sa(nodes, cost_dict, neighbourhood)
def run_sa(nodes, cost_dict, neighbourhood):
print('running SA......')
print('======================================================================')
# Set current temperature to an initial temperature t = t0
current_temp = initial_temp
# Set current solution to an initial solution s=s0
path = nodes
cost = 0
for i in range(len(path)):
cost += cost_dict[(path[i-1], path[i])]
print(f'initial path: {path}')
print(f'initial cost: {cost}')
print('======================================================================')
# plot the initial path
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
xs = [node[0] for node in path]
ys = [node[1] for node in path]
zs = [node[2] for node in path]
ax.plot(xs, ys, zs)
plt.show()
# do SA
while current_temp > minimum_temp:
for i in range(num_iterations):
# Select a solution si from the neighborhood N (s)
new_path, new_cost = select_a_solution(
neighbourhood,
# use the current temp and current iteration as an unique seed
current_temp * num_iterations + i,
path,
cost_dict
)
# Calculate change in cost C according to the new solution
cost_diff = new_cost - cost
# If C < 0 then accept new solution (it is improving) -> s= si
if cost_diff < 0:
path = new_path
cost = new_cost
# Else generate a random number x in the range (0,1)
else:
x = random.random()
# If x<e c/t then accept new solution s=si
if x < calculate_acceptance_probability(current_temp, cost_diff):
path = new_path
cost = new_cost
# Decrease t using alpha
current_temp = decrement_temp(current_temp)
# plot the final path
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
xs = [node[0] for node in path]
ys = [node[1] for node in path]
zs = [node[2] for node in path]
ax.plot(xs, ys, zs)
plt.show()
print(f'final path: {path}')
print(f'final cost: {cost}')
def select_a_solution(neighbourhood, seed, current_path, cost_dict):
# use the current temp and current iteration as an unique seed
random.seed(seed)
# get a random neighbour from the neighbourhood
neighbour = neighbourhood[random.randint(0, len(neighbourhood)-1)]
# get the index of the first swapped node
swap_first = current_path.index(neighbour[0])
# get the index of the second swapped node
swap_second = current_path.index(neighbour[1])
# make the swap
new_path = current_path
temp = current_path[swap_first]
new_path[swap_first] = current_path[swap_second]
new_path[swap_second] = temp
# calculate the new cost
new_cost = 0
for i in range(len(new_path)):
new_cost += cost_dict[(new_path[i - 1], new_path[i])]
return new_path, new_cost
def calculate_acceptance_probability(current_temp, cost_diff):
return math.exp(- cost_diff / current_temp)
# Linear rule
def decrement_temp(current_temp):
return current_temp - cooling_rate
def lon_to_x(lon, min_lon):
return int((float(lon) - float(min_lon)) * 10000)
def lat_to_y(lat, min_lat):
return int((float(lat) - float(min_lat)) * 10000)
def generate_nodes(x_max, y_max, z_max, number):
destinations = []
for i in range(number):
destinations.append(
(int(x_max * random.random()), int(y_max * random.random()), int(z_max * random.random())))
return destinations
def calculate_cost(a, b):
return math.sqrt((a[0] - b[0]) ** 2 + (a[1] - b[1]) ** 2 + (a[2] - b[2]) ** 2)
def build_neighbourhood(destinations):
neighbourhood = []
# while len(neighbourhood) <= neighbourhood_size:
for i in range(len(destinations)):
for j in range(i+1, len(destinations)):
neighbourhood.append((destinations[i], destinations[j]))
return neighbourhood
def build_cost_dict(nodes):
cost_dict = {}
for node1 in nodes:
for node2 in nodes:
if node1 != node2:
cost_dict[(node1, node2)] = calculate_cost(node1, node2)
return cost_dict
if __name__ == '__main__':
main()