diff --git a/README.md b/README.md index 8d5303ca..07c42ea3 100644 --- a/README.md +++ b/README.md @@ -103,6 +103,10 @@ This is the same [Boltzmann Wealth](https://github.com/projectmesa/mesa-examples This model is based on the NetLogo [Virus on a Network](https://ccl.northwestern.edu/netlogo/models/VirusonaNetwork) model. +### [Ant System for Traveling Salesman Problem](https://github.com/projectmesa/mesa-examples/tree/main/examples/aco_tsp) + +This is based on Dorigo's Ant System "Swarm Intelligence" algorithm for generating solutions for the Traveling Salesman Problem. + ## Visualization Examples ### [Boltzmann Wealth Model (Experimental)](https://github.com/projectmesa/mesa-examples/tree/main/examples/boltzmann_wealth_model_experimental) diff --git a/examples/aco_tsp/README.md b/examples/aco_tsp/README.md new file mode 100644 index 00000000..a810738d --- /dev/null +++ b/examples/aco_tsp/README.md @@ -0,0 +1,55 @@ +Ant System for the Traveling Salesman Problem +======================== + +This is an implementation of the Ant System (AS) algorithm for solving the Traveling Salesman Problem (TSP). This example uses Mesa's Network Grid to model a TSP by representing cities as nodes and the possible paths between them as edges. Ants are then modeled as Mesa Agents that generate solutions by traversing the network using a "swarm intelligence" algorithm. + +When an ant is choosing its next city, it consider both a pheromone trail laid down by previous ants and a greedy heuristic based on city proximity. Pheromone evaporates over time and the strength of new pheromone trail laid by an ant is proportional to the quality of its TSP solution. This produces an emergent solution as the pheromone trail is continually updated and guides ants to high quality solutions as they are discovered. + +As this model runs, more pheromone will be laid on better solutions, and less traveled paths will have their pheromone evaporate. Ants will therefore reinforce good paths and abandon bad ones. Since decisions are ultimately samples from a weighted probability distribution, ants will sometimes explore unlikely paths, which might lead to new strong solutions that will be reflected in the updated pheromone levels. + +Here, we plot the best solution per iteration, the best solution so far in all iterations, and a graph representation where the edge width is proportional to the pheromone quantity. You will quickly see most of the edges in the fully connected graph disappear and a subset of the paths emerge as reasonable candidates in the final TSP solution. + +## How to run +To launch the interactive visualization, run `solara run app.py` in this directory. Tune the $\alpha$ and $\beta$ parameters to modify how much the pheromone and city proximity influence the ants' decisions, respectively. See the Algorithm details section for more. + +Alternatively, to run for a fixed number of iterations, run `python run_tsp.py` from this directory (and update that file with the parameters you want). + +## Algorithm details +Each agent/ant is initialized to a random city and constructs a solution by choosing a sequence of cities until all are visited, but none are visited more than once. Ants then deposit a "pheromone" signal on each path in their solution that is proportional to 1/d, where d is the final distance of the solution. This means shorter paths are given more pheromone. + +When an ant is on city $i$ and deciding which city to choose next, it samples randomly using the following probabilities of transition from city $i$ to $j$: + +$$ +p_{ij}^k = \frac{\tau_{ij}^\alpha \eta_{ij}^\beta}{\sum_{l \in J_i^k} \tau_{il}^\alpha \eta_{il}^\beta} +$$ + +where: +- $\tau_{ij}$ is the amount of path pheromone +- $\eta_{ij}$ the a greedy heuristic of desireability + - In this case, $\eta_{ij} = 1/d_{ij}$, where $d_{ij}$ is the distance between + cities +- $\alpha$ is a hyperparameter setting the importance of the pheromone +- $\beta$ a hyperparameter for setting the importance of the greedy heuristic +- And the denominator sum is over $J_i^k$, which is the set of cities not yet + visited by ant $k$. + +In other words, $\alpha$ and $\beta$ are tuned to set the relative importance of the phermone trail left by prior ants, and the greedy heuristic of 1-over-distance. + +## Data collection +The following data is collected and can be used for further analysis: +- Agent-level (individual ants, reset after each iteration) + - `tsp_distance`: TSP solution distance + - `tsp_solution`: TSP solution path +- Model-level (collection of ants over many iterations) + - `num_steps`: number of algorithm iterations, where one step means each ant generates a full TSP solution and the pheromone trail is updated + - `best_distance`: the distance of the best path found in all iterations + - This is the best solution yet and can only stay flat or improve over time + - `best_distance_iter`: the distance of the best path of all ants in a single iteration + - This changes over time as the ant colony explores different solutions and can be used to understand the explore/exploit trade-off. E.g., if the colony quickly finds a good solution, but then this value trends upward and stays high, then this suggests the ants are stuck re-inforcing a suboptimal solution. + - `best_path`: the best path found in all iterations + +## References +- Original paper: Dorigo, M., Maniezzo, V., & Colorni, A. (1996). Ant system: optimization by a +colony of cooperating agents. IEEE transactions on systems, man, and cybernetics, +part b (cybernetics), 26(1), 29-41. +- Video series of this code being implemented: https://www.youtube.com/playlist?list=PLSgGvve8UweGk2TLSO-q5OSH59Q00ZxCQ \ No newline at end of file diff --git a/examples/aco_tsp/aco_tsp/__init__.py b/examples/aco_tsp/aco_tsp/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/examples/aco_tsp/aco_tsp/data/kroA100.tsp b/examples/aco_tsp/aco_tsp/data/kroA100.tsp new file mode 100644 index 00000000..b9b48209 --- /dev/null +++ b/examples/aco_tsp/aco_tsp/data/kroA100.tsp @@ -0,0 +1,107 @@ +NAME: kroA100 +TYPE: TSP +COMMENT: 100-city problem A (Krolak/Felts/Nelson) +DIMENSION: 100 +EDGE_WEIGHT_TYPE : EUC_2D +NODE_COORD_SECTION +1 1380 939 +2 2848 96 +3 3510 1671 +4 457 334 +5 3888 666 +6 984 965 +7 2721 1482 +8 1286 525 +9 2716 1432 +10 738 1325 +11 1251 1832 +12 2728 1698 +13 3815 169 +14 3683 1533 +15 1247 1945 +16 123 862 +17 1234 1946 +18 252 1240 +19 611 673 +20 2576 1676 +21 928 1700 +22 53 857 +23 1807 1711 +24 274 1420 +25 2574 946 +26 178 24 +27 2678 1825 +28 1795 962 +29 3384 1498 +30 3520 1079 +31 1256 61 +32 1424 1728 +33 3913 192 +34 3085 1528 +35 2573 1969 +36 463 1670 +37 3875 598 +38 298 1513 +39 3479 821 +40 2542 236 +41 3955 1743 +42 1323 280 +43 3447 1830 +44 2936 337 +45 1621 1830 +46 3373 1646 +47 1393 1368 +48 3874 1318 +49 938 955 +50 3022 474 +51 2482 1183 +52 3854 923 +53 376 825 +54 2519 135 +55 2945 1622 +56 953 268 +57 2628 1479 +58 2097 981 +59 890 1846 +60 2139 1806 +61 2421 1007 +62 2290 1810 +63 1115 1052 +64 2588 302 +65 327 265 +66 241 341 +67 1917 687 +68 2991 792 +69 2573 599 +70 19 674 +71 3911 1673 +72 872 1559 +73 2863 558 +74 929 1766 +75 839 620 +76 3893 102 +77 2178 1619 +78 3822 899 +79 378 1048 +80 1178 100 +81 2599 901 +82 3416 143 +83 2961 1605 +84 611 1384 +85 3113 885 +86 2597 1830 +87 2586 1286 +88 161 906 +89 1429 134 +90 742 1025 +91 1625 1651 +92 1187 706 +93 1787 1009 +94 22 987 +95 3640 43 +96 3756 882 +97 776 392 +98 1724 1642 +99 198 1810 +100 3950 1558 +EOF \ No newline at end of file diff --git a/examples/aco_tsp/aco_tsp/model.py b/examples/aco_tsp/aco_tsp/model.py new file mode 100644 index 00000000..1683aee9 --- /dev/null +++ b/examples/aco_tsp/aco_tsp/model.py @@ -0,0 +1,251 @@ +from dataclasses import dataclass + +import mesa +import networkx as nx +import numpy as np + + +@dataclass +class NodeCoordinates: + city: int + x: float + y: float + + @classmethod + def from_line(cls, line: str): + city, x, y = line.split() + return cls(int(city), float(x), float(y)) + + +class TSPGraph: + def __init__(self, g: nx.Graph, pheromone_init: float = 1e-6): + self.g = g + self.pheromone_init = pheromone_init + self._add_edge_properties() + + @property + def pos(self): + return {k: v["pos"] for k, v in dict(self.g.nodes.data()).items()} + + @property + def cities(self): + return list(self.g.nodes) + + @property + def num_cities(self): + return len(self.g.nodes) + + def _add_edge_properties(self): + for u, v in self.g.edges(): + u_x, u_y = self.g.nodes[u]["pos"] + v_x, v_y = self.g.nodes[v]["pos"] + self.g[u][v]["distance"] = ((u_x - v_x) ** 2 + (u_y - v_y) ** 2) ** 0.5 + self.g[u][v]["visibility"] = 1 / self.g[u][v]["distance"] + self.g[u][v]["pheromone"] = self.pheromone_init + + @classmethod + def from_random(cls, num_cities: int, seed: int = 0) -> "TSPGraph": + g = nx.random_geometric_graph(num_cities, 2.0, seed=seed).to_directed() + + return cls(g) + + @classmethod + def from_tsp_file(cls, file_path: str) -> "TSPGraph": + with open(file_path) as f: + lines = f.readlines() + # Skip lines until reach the text "NODE_COORD_SECTION" + while lines.pop(0).strip() != "NODE_COORD_SECTION": + pass + + g = nx.Graph() + for line in lines: + if line.strip() == "EOF": + break + node_coordinate = NodeCoordinates.from_line(line) + + g.add_node( + node_coordinate.city, pos=(node_coordinate.x, node_coordinate.y) + ) + + # Add edges between all nodes to make a complete graph + for u in g.nodes(): + for v in g.nodes(): + if u == v: + continue + g.add_edge(u, v) + + return cls(g) + + +class AntTSP(mesa.Agent): + """ + An agent + """ + + def __init__(self, unique_id, model, alpha: float = 1.0, beta: float = 5.0): + """ + Customize the agent + """ + self.unique_id = unique_id + self.alpha = alpha + self.beta = beta + super().__init__(unique_id, model) + self._cities_visited = [] + self._traveled_distance = 0 + self.tsp_solution = [] + self.tsp_distance = 0 + + def calculate_pheromone_delta(self, q: float = 100): + results = {} + for idx, start_city in enumerate(self.tsp_solution[:-1]): + end_city = self.tsp_solution[idx + 1] + results[(start_city, end_city)] = q / self.tsp_distance + + return results + + def decide_next_city(self): + # Random + # new_city = self.random.choice(list(self.model.all_cities - set(self.cities_visited))) + # Choose closest city not yet visited + g = self.model.grid.G + current_city = self.pos + neighbors = list(g.neighbors(current_city)) + candidates = [n for n in neighbors if n not in self._cities_visited] + if len(candidates) == 0: + return current_city + + # p_ij(t) = 1/Z*[(tau_ij)**alpha * (1/distance)**beta] + results = [] + for city in candidates: + val = ( + (g[current_city][city]["pheromone"]) ** self.alpha + * (g[current_city][city]["visibility"]) ** self.beta + ) + results.append(val) + + results = np.array(results) + norm = results.sum() + results /= norm + + new_city = self.model.random.choices(candidates, weights=results)[0] + + return new_city + + def step(self): + """ + Modify this method to change what an individual agent will do during each step. + Can include logic based on neighbors states. + """ + g = self.model.grid.G + for idx in range(self.model.num_cities - 1): + # Pick a random city that isn't in the list of cities visited + current_city = self.pos + new_city = self.decide_next_city() + self._cities_visited.append(new_city) + self.model.grid.move_agent(self, new_city) + self._traveled_distance += g[current_city][new_city]["distance"] + + self.tsp_solution = self._cities_visited.copy() + self.tsp_distance = self._traveled_distance + self._cities_visited = [] + self._traveled_distance = 0 + + +class AcoTspModel(mesa.Model): + """ + The model class holds the model-level attributes, manages the agents, and generally handles + the global level of our model. + + There is only one model-level parameter: how many agents the model contains. When a new model + is started, we want it to populate itself with the given number of agents. + + The scheduler is a special model component which controls the order in which agents are activated. + """ + + def __init__( + self, + num_agents: int = 20, + tsp_graph: TSPGraph = TSPGraph.from_random(20), + max_steps: int = int(1e6), + ant_alpha: float = 1.0, + ant_beta: float = 5.0, + ): + super().__init__() + self.num_agents = num_agents + self.tsp_graph = tsp_graph + self.num_cities = tsp_graph.num_cities + self.all_cities = set(range(self.num_cities)) + self.max_steps = max_steps + self.schedule = mesa.time.RandomActivation(self) + self.grid = mesa.space.NetworkGrid(tsp_graph.g) + + for i in range(self.num_agents): + agent = AntTSP(unique_id=i, model=self, alpha=ant_alpha, beta=ant_beta) + self.schedule.add(agent) + + city = tsp_graph.cities[self.random.randrange(self.num_cities)] + self.grid.place_agent(agent, city) + agent._cities_visited.append(city) + + self.num_steps = 0 + self.best_path = None + self.best_distance = float("inf") + self.best_distance_iter = float("inf") + # Re-initialize pheromone levels + tsp_graph._add_edge_properties() + + self.datacollector = mesa.datacollection.DataCollector( + model_reporters={ + "num_steps": "num_steps", + "best_distance": "best_distance", + "best_distance_iter": "best_distance_iter", + "best_path": "best_path", + }, + agent_reporters={ + "tsp_distance": "tsp_distance", + "tsp_solution": "tsp_solution", + }, + ) + + self.running = True + + def update_pheromone(self, q: float = 100, ro: float = 0.5): + # tau_ij(t+1) = (1-ro)*tau_ij(t) + delta_tau_ij(t) + # delta_tau_ij(t) = sum_k^M {Q/L^k} * I[i,j \in T^k] + delta_tau_ij = {} + for k, agent in enumerate(self.schedule.agents): + delta_tau_ij[k] = agent.calculate_pheromone_delta(q) + + for i, j in self.grid.G.edges(): + # Evaporate + tau_ij = (1 - ro) * self.grid.G[i][j]["pheromone"] + # Add ant's contribution + for k, delta_tau_ij_k in delta_tau_ij.items(): + tau_ij += delta_tau_ij_k.get((i, j), 0.0) + + self.grid.G[i][j]["pheromone"] = tau_ij + + def step(self): + """ + A model step. Used for collecting data and advancing the schedule + """ + self.datacollector.collect(self) + self.schedule.step() + self.num_steps += 1 + self.update_pheromone() + + # Check len of cities visited by an agent + best_instance_iter = float("inf") + for agent in self.schedule.agents: + # Check for best path + if agent.tsp_distance < self.best_distance: + self.best_distance = agent.tsp_distance + self.best_path = agent.tsp_solution + + if agent.tsp_distance < best_instance_iter: + best_instance_iter = agent.tsp_distance + + self.best_distance_iter = best_instance_iter + + if self.num_steps >= self.max_steps: + self.running = False diff --git a/examples/aco_tsp/app.py b/examples/aco_tsp/app.py new file mode 100644 index 00000000..068b2826 --- /dev/null +++ b/examples/aco_tsp/app.py @@ -0,0 +1,76 @@ +""" +Configure visualization elements and instantiate a server +""" + +import networkx as nx +import solara +from aco_tsp.model import AcoTspModel, TSPGraph +from matplotlib.figure import Figure +from mesa.visualization import SolaraViz + + +def circle_portrayal_example(agent): + return {"node_size": 20, "width": 0.1} + + +tsp_graph = TSPGraph.from_tsp_file("aco_tsp/data/kroA100.tsp") +model_params = { + "num_agents": tsp_graph.num_cities, + "tsp_graph": tsp_graph, + "ant_alpha": { + "type": "SliderFloat", + "value": 1.0, + "label": "Alpha: pheromone exponent", + "min": 0.0, + "max": 10.0, + "step": 0.1, + }, + "ant_beta": { + "type": "SliderFloat", + "value": 5.0, + "label": "Beta: heuristic exponent", + "min": 0.0, + "max": 10.0, + "step": 0.1, + }, +} + + +def make_graph(model): + fig = Figure() + ax = fig.subplots() + ax.set_title("Cities and pheromone trails") + graph = model.grid.G + pos = model.tsp_graph.pos + weights = [graph[u][v]["pheromone"] for u, v in graph.edges()] + # normalize the weights + weights = [w / max(weights) for w in weights] + + nx.draw( + graph, + ax=ax, + pos=pos, + node_size=10, + width=weights, + edge_color="gray", + ) + + solara.FigureMatplotlib(fig) + + +def ant_level_distances(model): + # ant_distances = model.datacollector.get_agent_vars_dataframe() + # Plot so that the step index is the x-axis, there's a line for each agent, + # and the y-axis is the distance traveled + # ant_distances['tsp_distance'].unstack(level=1).plot(ax=ax) + pass + + +page = SolaraViz( + AcoTspModel, + model_params, + space_drawer=None, + measures=["best_distance_iter", "best_distance", make_graph], + agent_portrayal=circle_portrayal_example, + play_interval=1, +) diff --git a/examples/aco_tsp/run_tsp.py b/examples/aco_tsp/run_tsp.py new file mode 100644 index 00000000..8b33e5ce --- /dev/null +++ b/examples/aco_tsp/run_tsp.py @@ -0,0 +1,47 @@ +from collections import defaultdict + +import matplotlib.pyplot as plt +from aco_tsp.model import AcoTspModel, TSPGraph + + +def main(): + # tsp_graph = TSPGraph.from_random(num_cities=20, seed=1) + tsp_graph = TSPGraph.from_tsp_file("aco_tsp/data/kroA100.tsp") + model_params = { + "num_agents": tsp_graph.num_cities, + "tsp_graph": tsp_graph, + } + number_of_episodes = 50 + + results = defaultdict(list) + + best_path = None + best_distance = float("inf") + + model = AcoTspModel(**model_params) + + for e in range(number_of_episodes): + # model = AcoTspModel(**model_params) + model.step() + results["best_distance"].append(model.best_distance) + results["best_path"].append(model.best_path) + print( + f"Episode={e + 1}; Min. distance={model.best_distance:.2f}; pheromone_1_8={model.grid.G[17][15]['pheromone']:.4f}" + ) + if model.best_distance < best_distance: + best_distance = model.best_distance + best_path = model.best_path + print(f"New best distance: distance={best_distance:.2f}") + + print(f"Best distance: {best_distance:.2f}") + print(f"Best path: {best_path}") + # print(model.datacollector.get_model_vars_dataframe()) + + _, ax = plt.subplots() + ax.plot(results["best_distance"]) + ax.set(xlabel="Episode", ylabel="Best distance", title="Best distance per episode") + plt.show() + + +if __name__ == "__main__": + main()