640 KiB
Algorithm project¶
Project context ¶
ADEME has launched a call to promote new mobility solutions. CesiCDP, with its team of 4 people, is responding to this call by focusing on the management of delivery rounds. The aim is to minimize the total duration of the round, taking into account traffic and other constraints. Attractive financing is at stake to develop CesiCDP's business. Adding further constraints to the problem will make the study more realistic and more difficult to solve, but it will demonstrate the team's expertise.
Our aim is to develop an algorithm that will enable us to pass through all the delivery points in an optimized time.
Selected constraint¶
The constraint we chose was the following:
- To have several trucks available simultaneously to make deliveries.
Formulation of the problem ¶
Consider a graph $G=(V,E)$, where $V$ is the set of cities (or delivery points) and $E$ is the set of roads between cities. There are $k$ trucks available to make deliveries.
The problem is to find a route for each truck, so that all deliveries are made in the shortest possible time, both to and from the depot.
The problem we have with the above constraints is the VRP (Vehicle Routing Problem).
Problem constraints ¶
List of problem constraints:
- All customers must be served
- A customer can only be served by one vehicle.
- When leaving a customer, a vehicle can only go to one other customer.
We will therefore assign each customer to a route served by a single vehicle.
Demonstrating the complexity of the vehicle routing problem (VRP) ¶
Introduction¶
The Vehicle Routing Problem (VRP) is an extension of the Traveling Salesman Problem (TSP), and is known to be an NP-hard problem. We will demonstrate the complexity of VRP based on the Hamiltonian chain problem, which is known to be NP-complete.
Preliminary definitions¶
Hamiltonian chain problem**: the problem consists in determining whether a given undirected graph has a Hamiltonian chain, i.e. a path that visits each vertex exactly once.
Traveling Salesman Problem (TSP)**: the problem consists in finding the shortest possible circuit that visits each city in a given set of cities and returns to the original city.
Vehicle Routing Problem (VRP)**: the problem consists in delivering a series of customers with several vehicles while minimizing the total cost, such as delivery time or distance traveled.
Proof of the complexity of the TSP.¶
The TSP is an extension of the Hamiltonian chain problem. In fact, a special case of the TSP is the Hamiltonian chain problem, in which all edges have the same weight (or cost). In this case, finding the shortest circuit in the TSP is equivalent to finding a Hamiltonian chain in the graph. Since the Hamiltonian chain problem is NP-complete, the TSP must be at least as difficult, so the TSP is NP-complete.
Proof of the complexity of the VRP.¶
The VRP is an extension of the TSP. In fact, a special case of the VRP is the TSP where there is only one vehicle available to make deliveries. In this case, finding the solution to the VRP is the same as finding the solution to the TSP.
Since the TSP is NP-complete, the VRP must be at least as difficult, so the VRP is NP-difficult. Furthermore, VRP introduces additional constraints, such as multiple vehicles and potentially vehicle capacities and time windows, making it even more complex.
Conclusion¶
In conclusion, we have shown that VRP is an NP-hard problem by reducing it to the NP-complete TSP problem, which in turn can be reduced to the NP-complete Hamiltonian chain problem. This demonstration highlights the difficulty of solving the VRP, particularly for large instances. In practice, (meta)heuristic and approximate methods are often used to solve the VRP, as we shall see later.
Mathematical modeling ¶
Set and parameters¶
$V=\{0,1,2,...,n_v\}$ : the set of cities, where 0 is the base (or depot), $1,2,...,n_v$ are the delivery cities. $n_v+1$ will be the depot for the return.
$K=\{1,2,...,k\}$ : all trucks.
$E$ represents the arcs between two customers $i,j \in V$
$G=(V,E)$ : the graph with V as vertices (cities) and E as edges
$d_{ij}$ : distance (or travel time) from city i to city $j$ (travel cost)
$M$ : a great constant.
Decision variables¶
$x_{ijk}$ : binary variable worth 1 if truck $k$ moves from city $i$ to city $j$, and 0 otherwise.
Objective function¶
$$\min \sum_{k∈K} \sum_{i∈V} \sum_{j∈V} d_{ij} x_{ijk} $$
VRP constraints¶
Each city is visited once and only once: $$\sum_{k \in K} \sum_{j \in V} x_{ijk} = 1, \forall i \in V, i \ne 0$$
If a truck visits a city, it must leave it: $$\sum_{i \in V} x_{ijk} = \sum_{j \in V} x_{ijk}, \forall k \in K, \forall i \in V, \forall j \in V $$
Sub-tour elimination constraint (to ensure that each truck completes a full tour) : $$\sum_{i \in S, j \notin S} x_{ijk} \geq 1, \forall k \in K, \forall \; subset \; S \; de \; V, 0 \in S, S \ne V $$
VRPTW constraints¶
Each city is visited once and only once: $$\sum_{k \in K} \sum_{j \in V}x_{ijk} = 1, \forall i \in V, i \ne 0 $$
If a truck visits a city, it must leave it: ∑𝑖∈𝑉𝑥𝑖𝑗𝑘=∑𝑗∈𝑉𝑥𝑖𝑗𝑘, ∀𝑘∈𝐾, ∀𝑖∈𝑉, ∀𝑗∈𝑉
Sub-tour elimination constraint (to ensure that each truck completes a full tour): ∑𝑖∈𝑆,𝑗∉𝑆𝑥𝑖𝑗𝑘≥1, ∀𝑘∈𝐾, ∀𝑠𝑢𝑏𝑠𝑒𝑡𝑆𝑑𝑒𝑉, 0∈𝑆, 𝑆≠𝑉
Each truck arrives within the time window of each city: a_ik + t_ij ≤ a_jk + M(1 - x_ijk), ∀k∈K, ∀(i, j)∈A, ∀a_ik ∈ [e_i, l_i]
Each truck leaves the city after the service time: a_ik + s_i = d_ik, ∀k∈K, ∀i∈V, ∀a_ik ∈ [e_i, l_i]
Note: Here, a_ik is the arrival time of vehicle k at node i, t_ij is the travel time from node i to node j, M is a large positive number, x_ijk is a binary variable that is 1 if vehicle k travels from node i to node j and 0 otherwise, e_i is the earliest time that a service can start at node i, l_i is the latest time that a service can start at node i, and s_i is the service time at node i.
These constraints take into account not just the need to visit each city once, but also the requirement to arrive and depart within specified time windows, which makes the problem more complex than the standard VRP.
Resolution algorithm¶
As part of our project, we plan to adopt three different heuristic algorithms to solve our inter-city route optimization problem. Our aim is to compare them in order to determine which one will provide us with the route closest to the optimum. The three algorithms we will be implementing are as follows:
- Ant colony algorithm
- Simulated annealing algorithm
- The particle swarm algorithm
We will then examine their implementation and operation to evaluate their performance.
Ant colony algorithm¶
The ant colony algorithm is an approach inspired by the behavior of ants in the search for optimal paths. Using heuristics and pheromones, this algorithm enables optimization problems to be solved efficiently.
# Variables to edit: nb_ville = 100 max_coords = 1000 nb_truck = 4 max_time = 4 nb_ants = 10 alpha = 1 beta = 6 evaporation = 0.5 intensification = 2
# Import necessary libraries import matplotlib.pyplot as plt import random, time from tests.libs.clustering import split_tour_across_clusters from tests.libs.aco import AntColony, total_distance # Function to generate city coordinates def generate_cities(nb, max_coords=1000): # Generate a list of nb cities, each with random coordinates within the range of max_coords return [random.sample(range(max_coords), 2) for _ in range(nb)] # Calculate maximum time per cluster max_time_per_cluster = max_time / nb_truck # Track the time taken to generate cities start_time_generate = time.time() # Generate cities cities = generate_cities(nb_ville, max_coords) # Set the first city's coordinates at the center cities[0] = [max_coords/2, max_coords/2] stop_time_generate = time.time() # Track the time taken to split cities across clusters start_time_split = time.time() # Split cities into clusters clusters = split_tour_across_clusters(cities, nb_truck) stop_time_split = time.time() # Print out the time taken for each process print("\n---- TIME ----") print("Generate cities time: ", stop_time_generate - start_time_generate) print("Split cities time: ", stop_time_split - start_time_split) # Initialize a new figure for displaying paths plt.figure() plt.title("Ant Colony Optimization") # List of colors for plotting colors = [ '#1f77b4', # Bleu moyen '#ff7f0e', # Orange '#2ca02c', # Vert '#d62728', # Rouge '#9467bd', # Violet '#8c564b', # Marron '#e377c2', # Rose '#7f7f7f', # Gris '#bcbd22', # Vert olive '#17becf', # Turquoise '#1b9e77', # Vert Teal '#d95f02', # Orange foncé '#7570b3', # Violet moyen '#e7298a', # Fuchsia '#66a61e', # Vert pomme '#e6ab02', # Jaune or '#a6761d', # Bronze '#666666', # Gris foncé '#f781bf', # Rose clair '#999999', # Gris moyen ] # Initialize a list for storing the best routes best_routes = [] print("\n---- CLUSTER ----") # Iterate through each cluster for i, cluster_indices in enumerate(clusters.values()): # Select a color for the cluster color = colors[i % len(colors)] # Retrieve city coordinates for the cluster cluster_cities = [cities[index] for index in cluster_indices] # Run the AntColony algorithm and store the best route ant_colony = AntColony(cluster_cities, n_ants=nb_ants, max_time=max_time_per_cluster, alpha=alpha, beta=beta, evaporation=evaporation, intensification=intensification) best_route = ant_colony.run() best_routes.append((best_route, color)) # Print the total distance for the cluster print("Total distance for cluster {} ({} cities) : {}".format(i+1, len(cluster_cities), total_distance(best_route))) # Calculate and print total distance for all clusters full_total_distance = 0 for route, color in best_routes: full_total_distance += total_distance(route) print("Total distance for all clusters: ", full_total_distance) # Plot each route for i, (route, color) in enumerate(best_routes): x = [city[0] for city in route] y = [city[1] for city in route] # Close the loop of the route x.append(x[0]) y.append(y[0]) # Plot the route plt.plot(x, y, color=color, marker='o', linestyle='-', label=f"Cluster {i}") # Display the plot plt.show()
---- TIME ---- Generate cities time: 0.0 Split cities time: 0.04548168182373047 ---- CLUSTER ---- Total distance for cluster 1 (24 cities) : 2315.875564700583 Total distance for cluster 2 (27 cities) : 2161.8716667459307 Total distance for cluster 3 (28 cities) : 2309.9191397305735 Total distance for cluster 4 (28 cities) : 2261.7683336604277 Total distance for all clusters: 9049.434704837515
Simulated annealing algorithm¶
The simulated annealing algorithm is an optimization technique inspired by the annealing process used in metallurgy. It allows probabilistic exploration of the solution space, occasionally accepting moves that may initially seem less favorable. This approach makes it possible to find global solutions and avoid getting stuck in local optima.
# Variables to edit: nb_ville = 100 max_coords = 1000 nb_truck = 4 temperature = 10000 cooling_rate = 0.9999 temperature_ok = 0.001
# Import necessary libraries import matplotlib.pyplot as plt import random, time from tests.libs.clustering import split_tour_across_clusters from tests.libs.simulated_annealing import SimulatedAnnealing, total_distance # Function to generate city coordinates def generate_cities(nb, max_coords=1000): # Generate a list of nb cities, each with random coordinates within the range of max_coords return [random.sample(range(max_coords), 2) for _ in range(nb)] # Track the time taken to generate cities start_time_generate = time.time() # Generate cities cities = generate_cities(nb_ville, max_coords) # Set the first city's coordinates at the center cities[0] = [max_coords/2, max_coords/2] stop_time_generate = time.time() # Track the time taken to split cities across clusters start_time_split = time.time() # Split cities into clusters clusters = split_tour_across_clusters(cities, nb_truck) stop_time_split = time.time() # Print out the time taken for each process print("\n---- TIME ----") print("Generate cities time: ", stop_time_generate - start_time_generate) print("Split cities time: ", stop_time_split - start_time_split) # Initialize a new figure for displaying paths plt.figure() plt.title("Simulated Annealing") # List of colors for plotting colors = [ '#1f77b4', # Bleu moyen '#ff7f0e', # Orange '#2ca02c', # Vert '#d62728', # Rouge '#9467bd', # Violet '#8c564b', # Marron '#e377c2', # Rose '#7f7f7f', # Gris '#bcbd22', # Vert olive '#17becf', # Turquoise '#1b9e77', # Vert Teal '#d95f02', # Orange foncé '#7570b3', # Violet moyen '#e7298a', # Fuchsia '#66a61e', # Vert pomme '#e6ab02', # Jaune or '#a6761d', # Bronze '#666666', # Gris foncé '#f781bf', # Rose clair '#999999', # Gris moyen ] # Initialize a list for storing the best routes best_routes = [] # Iterate through each cluster for i, cluster_indices in enumerate(clusters.values()): # Select a color for the cluster color = colors[i % len(colors)] # Retrieve city coordinates for the cluster cluster_cities = [cities[index] for index in cluster_indices] # Run the SimulatedAnnealing algorithm and store the best route simulated_annealing = SimulatedAnnealing(cluster_cities, temperature=10000, cooling_rate=0.999, temperature_ok=0.01) best_route = simulated_annealing.run() best_routes.append((best_route, color)) # Print the total distance for the cluster print("Total distance for cluster {} ({} cities) : {}".format(i+1, len(cluster_cities), total_distance(best_route))) # Calculate and print total distance for all clusters full_total_distance = 0 for route, color in best_routes: full_total_distance += total_distance(route) print("Total distance for all clusters: ", full_total_distance) # Plot each route for i, (route, color) in enumerate(best_routes): x = [city[0] for city in route] y = [city[1] for city in route] # Close the loop of the route x.append(x[0]) y.append(y[0]) # Plot the route plt.plot(x, y, color=color, marker='o', linestyle='-', label=f"Cluster {i}") # Add legend to the plot plt.legend(loc="best") # Display the plot plt.show()
---- TIME ---- Generate cities time: 0.0009975433349609375 Split cities time: 0.04487872123718262 Total distance for cluster 1 (27 cities) : 2664.6125913825863 Total distance for cluster 2 (27 cities) : 2419.498935358293 Total distance for cluster 3 (27 cities) : 2792.878477974977 Total distance for cluster 4 (27 cities) : 2914.2725993190484 Total distance for all clusters: 10791.262604034904
Particle swarm algorithm¶
The particle swarm algorithm is an optimization method based on the collective behavior of animal swarms, such as birds or fish. It uses particles that move through the search space and communicate with each other to find optimal solutions. This algorithm is effective in exploring the solution space and finding local and global optima.
# Variables to edit: num_cities = 50 num_trucks = 5 num_particles = 50 max_iterations = 300 inertia_weight = 2 cognitive_weight = 1 social_weight = 1.5 incompleteness_probability = 0.8 start_time = 8 end_time = 18 infinite_distance_value = 1e6
# Importing necessary libraries import random import numpy as np import networkx as nx import matplotlib.pyplot as plt import matplotlib from IPython.display import display, Markdown import matplotlib.cm as cm from tests.libs.pso import Particle # Function to generate a matrix of distances between cities def generate_distance_matrix(num_cities, incompleteness_probability): # Initialize an empty matrix with zeros distances = np.zeros((num_cities, num_cities)) # Iterate over each city for i in range(num_cities): # For each city, iterate over every other city to calculate distance for j in range(i+1, num_cities): # With certain probability, some distances will be missing (infinite) if random.random() < incompleteness_probability: distances[i, j] = float('inf') distances[j, i] = float('inf') else: # Assign random integer value to the distances between i and j distances[i, j] = random.randint(1, 10) distances[j, i] = distances[i, j] return distances # Function to generate random time windows for each city def generate_time_windows(num_cities, start_time, end_time): time_windows = [] # For each city, generate a random start and end time for _ in range(num_cities): start = random.randint(start_time, end_time - 1) end = random.randint(start + 1, end_time) time_windows.append((start, end)) return time_windows # Function to generate a graph representation of the cities and the distances between them def generate_city_graph(distances, time_windows): num_cities = distances.shape[0] # Create an empty graph G = nx.Graph() # Add cities as nodes with their corresponding time windows for i in range(num_cities): G.add_node(i+1, time_window=time_windows[i]) # Add edges between cities with their corresponding distances for u in range(num_cities): for v in range(u+1, num_cities): if distances[u, v] != float('inf'): G.add_edge(u+1, v+1, distance=distances[u, v]) return G # Generate the distances and time windows distances = generate_distance_matrix(num_cities, incompleteness_probability) time_windows = generate_time_windows(num_cities, start_time, end_time) # Generate the city graph G = generate_city_graph(distances, time_windows) # Initialize the particles and the global best position and cost particles = [] global_best_position = np.random.permutation(range(1, num_cities+1)) global_best_cost = float('inf') for _ in range(num_particles): particle = Particle(num_cities, num_trucks, distances, time_windows, infinite_distance_value) particles.append(particle) particle_cost = particle.evaluate_cost() if particle_cost < global_best_cost: global_best_cost = particle_cost global_best_position = particle.position.copy() # Main loop for the particle swarm optimization algorithm iteration = 0 best_costs = [] while iteration < max_iterations: for particle in particles: particle.update_velocity(global_best_position, inertia_weight, cognitive_weight, social_weight) particle.update_position() particle_cost = particle.evaluate_cost() if particle_cost < global_best_cost: global_best_cost = particle_cost global_best_position = particle.position.copy() best_costs.append(global_best_cost) iteration += 1 # Display the best solution found best_position_str = ", ".join(str(node) for node in global_best_position) best_cost_str = str(global_best_cost) markdown_text = f"### Best solution found:\n\n- **City Positions**: {best_position_str}\n- **Total Cost**: {best_cost_str}" display(Markdown(markdown_text)) # ---------------- Plot the best cost for each iteration ---------------------------------------------------------------------- # Compute the layout positions for the nodes using the spring layout algorithm pos = nx.spring_layout(G) # Create labels for each node in the graph labels = {node: str(node) for node in G.nodes()} # Draw the nodes of the graph using the computed positions nx.draw_networkx_nodes(G, pos) # Draw the edges of the graph using the computed positions nx.draw_networkx_edges(G, pos) # Draw the labels for the nodes using the computed positions and labels nx.draw_networkx_labels(G, pos, labels) # Set the title of the plot plt.title("Graph of cities with time windows") # Turn off the axis display plt.axis("off") # Display the plot plt.show() # ---------------- Plot Truck Paths ---------------------------------------------------------------------- # Initialize an empty list for truck paths truck_paths = [] # Create a colormap for the trucks cmap = matplotlib.colormaps['gist_ncar'] colors = [cmap(i/num_trucks) for i in range(num_trucks)] # Split the global best position into separate paths for each truck for truck in range(num_trucks): path = [] for i in range(len(global_best_position)): if i % num_trucks == truck: city = global_best_position[i] path.append(int(city)) path.append(path[0]) # Making the path circular by adding the first city to the end of the path truck_paths.append(path) # Plot the truck paths on a graph plt.figure(figsize=(10, 6)) plt.title("Truck Paths") plt.axis("off") for i in range(num_trucks): path = truck_paths[i] path_str = " -> ".join(str(city) for city in path) text = f"Truck {i+1}: {path_str}" display(Markdown(text)) edges = [(path[j], path[j+1]) for j in range(len(path)-1)] nx.draw_networkx_edges(G, pos, edgelist=edges, edge_color=colors[i], label=f"Truck {i+1}") nx.draw_networkx_labels(G, pos, labels) plt.legend() plt.show() # ---------------- Plot path cost for Each Truck ---------------------------------------------------------------------- # Initialize an empty list for the total cost of each path path_costs = [] # Calculate the total cost of each path for path in truck_paths: path_cost = 0 for i in range(len(path) - 1): u, v = path[i] - 1, path[i + 1] - 1 # Adjust the indices # If the distance is infinite, use a replacement value. distance = distances[u, v] if distances[u, v] != float('inf') else infinite_distance_value path_cost += distance path_costs.append(path_cost) # Plot a histogram of the path costs plt.figure(figsize=(10, 6)) plt.bar(range(1, num_trucks + 1), path_costs, color=cm.rainbow(np.linspace(0, 1, num_trucks))) plt.xlabel("Truck") plt.ylabel("Path Cost") plt.title("Path Cost for Each Truck") plt.show() # ---------------- Plot Evolution of the Best Cost Over Iterations ---------------------------------------------------------------------- # Plot a graph of the evolution of the best cost plt.figure(figsize=(10, 6)) plt.plot(range(1, max_iterations + 1), best_costs) plt.title("Evolution of the Best Cost Over Iterations") plt.xlabel("Iteration") plt.ylabel("Best Cost") plt.grid(True) plt.show()
Best solution found:¶
- City Positions: 22, 18, 5, 15, 27, 4, 36, 33, 16, 45, 17, 24, 50, 20, 48, 7, 1, 8, 3, 35, 40, 6, 25, 11, 43, 28, 34, 30, 38, 29, 26, 23, 13, 39, 47, 12, 2, 31, 37, 9, 49, 44, 19, 10, 41, 42, 32, 21, 46, 14
- Total Cost: 86000793.0
Truck 1: 22 -> 4 -> 17 -> 7 -> 40 -> 28 -> 26 -> 12 -> 49 -> 42 -> 22
Truck 2: 18 -> 36 -> 24 -> 1 -> 6 -> 34 -> 23 -> 2 -> 44 -> 32 -> 18
Truck 3: 5 -> 33 -> 50 -> 8 -> 25 -> 30 -> 13 -> 31 -> 19 -> 21 -> 5
Truck 4: 15 -> 16 -> 20 -> 3 -> 11 -> 38 -> 39 -> 37 -> 10 -> 46 -> 15
Truck 5: 27 -> 45 -> 48 -> 35 -> 43 -> 29 -> 47 -> 9 -> 41 -> 14 -> 27
Test de la variation du coût optimal en faisant varier le paramètre "inertia weight"¶
# List of inertia values to test inertia_values = np.linspace(0.1, 5, 10) # Initialize list to store optimal cost for each inertia value optimal_costs = [] # Loop over each inertia value for inertia_weight in inertia_values: # Initialize particles and the global best cost particles = [] global_best_position = np.random.permutation(range(1, num_cities+1)) global_best_cost = float('inf') # Initialize particles and find the best cost for _ in range(num_particles): particle = Particle(num_cities, num_trucks, distances, time_windows, infinite_distance_value) particles.append(particle) particle_cost = particle.evaluate_cost() if particle_cost < global_best_cost: global_best_cost = particle_cost global_best_position = particle.position.copy() # Iterate to find the optimal solution iteration = 0 best_costs = [] while iteration < max_iterations: for particle in particles: particle.update_velocity(global_best_position, inertia_weight, cognitive_weight, social_weight) particle.update_position() particle_cost = particle.evaluate_cost() if particle_cost < global_best_cost: global_best_cost = particle_cost global_best_position = particle.position.copy() best_costs.append(global_best_cost) iteration += 1 # Store the optimal cost for this inertia value optimal_costs.append(global_best_cost) # Plot the evolution of the optimal cost as a function of inertia plt.figure(figsize=(10, 6)) plt.plot(inertia_values, optimal_costs) plt.title("Evolution of Optimal Cost as a Function of Inertia") plt.xlabel("Inertia") plt.ylabel("Optimal Cost") plt.show()
Test de la variation du coût optimal en faisant varier le paramètre "cognitive weight"¶
# List of cognitive weight values to test cognitive_values = np.linspace(0.1, 5, 10) # Initialize list to store optimal cost for each cognitive weight value optimal_costs = [] # Loop over each cognitive weight value for cognitive_weight in cognitive_values: # Initialize particles and the global best cost particles = [] global_best_position = np.random.permutation(range(1, num_cities+1)) global_best_cost = float('inf') # Initialize particles and find the best cost for _ in range(num_particles): particle = Particle(num_cities, num_trucks, distances, time_windows, infinite_distance_value) particles.append(particle) particle_cost = particle.evaluate_cost() if particle_cost < global_best_cost: global_best_cost = particle_cost global_best_position = particle.position.copy() # Iterate to find the optimal solution iteration = 0 best_costs = [] while iteration < max_iterations: for particle in particles: particle.update_velocity(global_best_position, inertia_weight, cognitive_weight, social_weight) particle.update_position() particle_cost = particle.evaluate_cost() if particle_cost < global_best_cost: global_best_cost = particle_cost global_best_position = particle.position.copy() best_costs.append(global_best_cost) iteration += 1 # Store the optimal cost for this cognitive weight value optimal_costs.append(global_best_cost) # Plot the evolution of the optimal cost as a function of cognitive weight plt.figure(figsize=(10, 6)) plt.plot(cognitive_values, optimal_costs) plt.title("Evolution of Optimal Cost as a Function of Cognitive Weight") plt.xlabel("Cognitive Weight") plt.ylabel("Optimal Cost") plt.show()