Solving optimal shipment route problem using Christofides–Serdyukov algorithm, Simplex and Reinforcement Learning

Solving optimal routes problem is a classic industry problem that has many practical applications in the supply chain optimization sector. One could look at it as a variant of the Traveling salesman problem , where we have a starting point, for example our product depot, and a set of cities to deliver to, each route between two cities has a cost , for example fuel cost and tolls and taxes. We could also assume that each road has a form of total available capacity, and each city has a certain demand and surplus it can give.

So now our problem is the following : We want to maximize the quanitity of shipped items whilst minimizing distance and cost and delivering to each single city.

it is a challenging combinatorial optimization problem , and in this article we’ll explore how we can solve it using the above mentioned methods.

We’ll start simple and then add more complexity as we go along

1- Data generation

Let us assume our shipping company is in France , our depot is in Paris and we have to deliver to multiple cities across the country. We’ll use networkx library in python to generate some graph data representing this problem.

import networkx as nx
import geopandas as gpd
from shapely.geometry import Point, LineString
import matplotlib.pyplot as plt

# Set the number of nodes (cities)
num_nodes = 10

# Define the city data
city_data = {
    'city': ['Paris', 'Marseille', 'Lyon', 'Toulouse', 'Nice', 'Nantes', 'Montpellier', 'Strasbourg', 'Bordeaux', 'Lille'],
    'lat': [48.8567, 43.2964, 45.7600, 43.6045, 43.7034, 47.2181, 43.6119, 48.5833, 44.8400, 50.6278],
    'lng': [2.3522, 5.3700, 4.8400, 1.4440, 7.2663, -1.5528, 3.8772, 7.7458, -0.5800, 3.0583]
}

# Create a complete graph with num_nodes nodes
G = nx.complete_graph(num_nodes)

# Add nodes to the graph with city names as labels
for i in range(num_nodes):
    G.nodes[i]['city'] = city_data['city'][i]

# Create a GeoDataFrame for the cities with Point geometries
gdf_cities = gpd.GeoDataFrame(city_data, geometry=gpd.points_from_xy(city_data['lng'], city_data['lat']))

# Create a GeoDataFrame for the routes (edges)
routes = []
for edge in G.edges():
    city1 = G.nodes[edge[0]]['city']
    city2 = G.nodes[edge[1]]['city']
    route = LineString([(gdf_cities[gdf_cities['city'] == city1].geometry.x.iloc[0], gdf_cities[gdf_cities['city'] == city1].geometry.y.iloc[0]),
                        (gdf_cities[gdf_cities['city'] == city2].geometry.x.iloc[0], gdf_cities[gdf_cities['city'] == city2].geometry.y.iloc[0])])
    routes.append(route)

gdf_routes = gpd.GeoDataFrame(geometry=routes)

# Create a plot for the map of France
fig, ax = plt.subplots(figsize=(10, 8))
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
france_shapefile = world[world.name == 'France']
france_shapefile.boundary.plot(ax=ax, color='lightgray', linewidth=0.8)
gdf_cities.plot(ax=ax, marker='o', color='red', markersize=100, alpha=0.7)
gdf_routes.plot(ax=ax, color='blue', linewidth=1, alpha=0.6)

# Add labels for the cities
for x, y, label in zip(gdf_cities.geometry.x, gdf_cities.geometry.y, gdf_cities['city']):
    ax.annotate(label, xy=(x, y), xytext=(3, 3), textcoords="offset points", color='black', fontsize=10)

# Set the axis limits to focus on France
ax.set_xlim(-5, 10)
ax.set_ylim(40, 54)

# Show the plot
plt.title('Cities and Routes in France')
plt.show()

Now let us start with a simple problem where we only want to minimize distance cost, we’ll assume the cost is 0.25euros per km , and compute the weightts for each edge as such.

# Compute the distance and add weight to each edge
from geopy.distance import geodesic
import numpy as np
for edge in G.edges():
    city1 = G.nodes[edge[0]]['city']
    city2 = G.nodes[edge[1]]['city']
    coord1 = (city_data['lat'][city_data['city'].index(city1)], city_data['lng'][city_data['city'].index(city1)])
    coord2 = (city_data['lat'][city_data['city'].index(city2)], city_data['lng'][city_data['city'].index(city2)])
    distance_km = geodesic(coord1, coord2).kilometers
    G.edges[edge]['weight'] = np.round(0.25 * distance_km)
# Create a plot for the graph data with assigned weights
fig, ax = plt.subplots(figsize=(10, 8))

# Create a GeoDataFrame for the cities with Point geometries
gdf_cities = gpd.GeoDataFrame(city_data, geometry=gpd.points_from_xy(city_data['lng'], city_data['lat']))

# Draw the cities as nodes
pos = {i: (city_data['lng'][i], city_data['lat'][i]) for i in range(num_nodes)}
nx.draw_networkx_nodes(G, pos, node_size=200, node_color='red', alpha=0.7, ax=ax)

# Draw the routes as edges with weights as labels
edges = nx.draw_networkx_edges(G, pos, edge_color='blue', alpha=0.6, width=1)
edge_labels = nx.get_edge_attributes(G, 'weight')
nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels, font_size=8, ax=ax)

# Add labels for the cities
for x, y, label in zip(gdf_cities.geometry.x, gdf_cities.geometry.y, gdf_cities['city']):
    ax.annotate(label, xy=(x, y), xytext=(3, 3), textcoords="offset points", color='black', fontsize=10)

# Set the axis limits to focus on the cities
ax.set_xlim(min(city_data['lng']) - 1, max(city_data['lng']) + 1)
ax.set_ylim(min(city_data['lat']) - 1, max(city_data['lat']) + 1)

# Show the plot
plt.title('Graph Data with Assigned Weights (Distances in Euros)')
plt.axis('equal')
plt.show()

2- Simple problem: Shotest path minimizing just the distance cost

Alright now we can try and find a solution for this problem, we’ll set our start point as Paris and use the Christofides–Serdyukov algorithm to approximate the solution, this can be done using the library networkx , specifically the function : nx.approximation.traveling_salesman_problem(G, cycle=True)), Then we’ll plot the solution found

# Convert the problem to finding a Hamiltonian cycle starting and ending at the specified start point
start_point = 'Paris'
# Add the start point to the city_data dictionary
start_lat = city_data['lat'][city_data['city'].index(start_point)]
start_lng = city_data['lng'][city_data['city'].index(start_point)]
city_data['city'].append(start_point)
city_data['lat'].append(start_lat)
city_data['lng'].append(start_lng)
# Create a new node representing the start point
new_node_id = num_nodes
G.add_node(new_node_id, city=start_point)

# Connect the new node to all other nodes
for i in range(num_nodes):
    city = G.nodes[i]['city']
    coord = (city_data['lat'][city_data['city'].index(city)], city_data['lng'][city_data['city'].index(city)])
    distance_km = geodesic(coord, (start_lat, start_lng)).kilometers
    G.add_edge(new_node_id, i, weight=0.25 * distance_km)

# Find the Hamiltonian cycle starting and ending at the specified start point
optimal_route = list(nx.approximation.traveling_salesman_problem(G, cycle=True))

# Remove the added node representing the start point
optimal_route.remove(new_node_id)

# Rearrange the route to start from the specified start point
start_index = optimal_route.index(city_data['city'].index(start_point))
optimal_route = optimal_route[start_index+1:] + optimal_route[:start_index+1]


# Remove the added node representing the start point
G.remove_node(new_node_id)

# Rearrange the route to start from the specified start point
start_index = optimal_route.index(new_node_id)
optimal_route = optimal_route[start_index+1:] + optimal_route[:start_index+1]
# Create a GeoDataFrame for the optimal route
route_points = [(city_data['lng'][i], city_data['lat'][i]) for i in optimal_route]
gdf_route = gpd.GeoDataFrame({'geometry': [LineString(route_points + [route_points[0]])]})

# Set the positions of nodes using spring_layout
node_positions = nx.spring_layout(G, seed=42, iterations=100, scale=30)
# Plot the optimal route in red color with arrows
optimal_route_edges = [(optimal_route[i], optimal_route[i + 1]) for i in range(len(optimal_route) - 1)] + [(optimal_route[-1], optimal_route[0])]
optimal_route_positions = {node: node_positions[node-1] if node != 0 else (0, 0) for node in optimal_route}

# Create a plot for the graph data with assigned weights
fig, ax = plt.subplots(figsize=(10, 8))

# Create a GeoDataFrame for the cities with Point geometries
gdf_cities = gpd.GeoDataFrame(city_data, geometry=gpd.points_from_xy(city_data['lng'], city_data['lat']))

# Draw the routes as edges with weights as labels
edges = nx.draw_networkx_edges(G, pos, edge_color='blue', alpha=0.6, width=1)
edge_labels = nx.get_edge_attributes(G, 'weight')
nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels, font_size=8, ax=ax)

# Add labels for the cities
for x, y, label in zip(gdf_cities.geometry.x, gdf_cities.geometry.y, gdf_cities['city']):
    ax.annotate(label, xy=(x, y), xytext=(3, 3), textcoords="offset points", color='black', fontsize=10)

# Plot the optimal route in red color with arrows
for i in range(len(optimal_route) - 1):
    city1 = city_data['city'][optimal_route[i]]
    city2 = city_data['city'][optimal_route[i + 1]]
    x_coords = [city_data['lng'][optimal_route[i]], city_data['lng'][optimal_route[i + 1]]]
    y_coords = [city_data['lat'][optimal_route[i]], city_data['lat'][optimal_route[i + 1]]]
    ax.plot(x_coords, y_coords, color='red', linewidth=2, alpha=0.6, zorder=1)

    # Add arrows to the line segments
    dx = x_coords[1] - x_coords[0]
    dy = y_coords[1] - y_coords[0]
    ax.arrow(x_coords[0], y_coords[0], dx, dy, head_width=0.1, head_length=0.1, fc='red', ec='red')

# Connect the last city to the first city in the optimal route
city1 = city_data['city'][optimal_route[-1]]
city2 = city_data['city'][optimal_route[0]]
x_coords = [city_data['lng'][optimal_route[-1]], city_data['lng'][optimal_route[0]]]
y_coords = [city_data['lat'][optimal_route[-1]], city_data['lat'][optimal_route[0]]]
ax.plot(x_coords, y_coords, color='red', linewidth=2, alpha=0.6, zorder=1)

# Draw the cities as nodes on top of the edges
#nx.draw_networkx_nodes(G, pos, node_size=200, node_color='red', alpha=0.7, ax=ax)

# Set the axis limits to focus on the cities
ax.set_xlim(min(city_data['lng']) - 1, max(city_data['lng']) + 1)
ax.set_ylim(min(city_data['lat']) - 1, max(city_data['lat']) + 1)

# Show the plot
plt.title('Graph Data with Assigned Weights (Distances in Euros) and Optimal Route')
plt.axis('equal')
plt.show()

Okay , so here the solution : Start at Paris -> Nantes -> Bordeaux -> Toulouse -> Montpellier -> Nice -> Marseille -> Lyon -> Strasbourg -> Lille and finally go back to Paris.

3- Adding more complexity : Minimzing the distance cost and transportation cost

Okay , Now let us assume each road has another specific cost like tolls and taxes , the cost per road is no longer linear to the distance in km. Let us try and simulate that by adding some random positive cost to each edge.

import random
# Add random positive cost to each edge
for edge in G.edges():
     G[edge[0]][edge[1]]['cost'] = random.randint(25, 200)
# Define a custom weight function that considers both distance and cost
def custom_weight(u, v, data):
    distance_weight = data['weight']  # Distance weight (previously calculated)
    cost_weight = data['cost']  # Custom cost weight for the edge
    return distance_weight + cost_weight

# Compute the shortest path with the custom weight function
optimal_route = nx.approximation.traveling_salesman_problem(G, cycle=True, weight=custom_weight)

Now let’s see the solution

So now the solution is Start at Paris -> Lille -> Paris -> Nantes -> Nice -> Montpellier -> Toulouse -> Paris -> Strasbourg -> Lyon -> Bordeaux -> Marseille -> Paris, note that we did not add a condition to not go back to a city, but to at least visit each city once.

4- Solving using an RL agent.

To solve using a RL agent, we need to define our environment, meaning the states: Which cities we are currently in (visited cities and cities not visited yet ), the set of possible actions (which city to visit next ) and the reward (negative reward representing the cost ) for each action.

For that we will use the following libraries in python: stable_baselines3 and gym

from stable_baselines3 import DQN
from stable_baselines3.common.env_util import DummyVecEnv
import gym
from gym import spaces
# Define the RL environment as a wrapper for the TSPEnvironment
class TSPGymEnv(gym.Env):
    def __init__(self, G):
        self.G = G
        self.num_nodes = G.number_of_nodes()
        self.reset()

        # Action space
        self.action_space = spaces.Discrete(self.num_nodes)

        # Observation space (current city)
        self.observation_space = spaces.Discrete(self.num_nodes)
        # Initialize the set of visited cities
        self.visited_cities = set()

        # Initialize the starting city as the current city
        self.current_city = 0
    def reset(self):
        self.visited_cities = set()
        self.current_city = random.choice(list(self.G.nodes()))
        self.visited_cities.add(self.current_city)
        return self.current_city



    def step(self, action):
        if self.current_city not in self.visited_cities:
            # Add the current city to the set of visited cities
            self.visited_cities.add(self.current_city)

        next_city = list(self.G[self.current_city])[action - 1 if action != 0 else 0]
        reward = -self.G[self.current_city][next_city]['cost']  # Negative cost as the reward
        done = len(self.visited_cities) == self.num_nodes
        if not done:
            self.current_city = next_city
        return next_city, reward, done, {}


# Create the RL environment
env = DummyVecEnv([lambda: TSPGymEnv(G)])
# Create the DQN model
model = DQN("MlpPolicy", env, verbose=1)

# Train the model
model.learn(total_timesteps=30000)
# Create a list to store the optimal route
optimal_route = [0]

# Continue predicting the next cities until all cities have been visited
while len(optimal_route) < num_nodes:
    # Convert the current_city to a numpy array
    current_city_array = np.array([current_city])

    # Predict the next action using the RL model
    action, _ = model.predict(current_city_array)

    # Get the next city based on the predicted action
    next_city = list(G[current_city])[action[0] - 1 if action[0] != 0 else 0]

    # If the next city has not been visited yet, add it to the optimal route
    if next_city not in optimal_route:
        optimal_route.append(next_city)

    # Update the current_city for the next iteration
    current_city = next_city

This is the solution found by the Agent : Start at Paris -> Lyon -> Toulouse -> Bordeaux -> Lille -> Nice -> Strasbourg -> Nantes -> Marseille -> Montpellier

5-Minimum-cost flow problem

Finally, Let us assume you a total inventory of d items to ship , each city has a required demand f and a possbile surplas it can give , and foreach route there is a maximum capacity of items you can ship on that specific route, and we want to find out the optimal route to follow so we can minimize total cost and also satisfy as much demand as possible. now the problem has really become quite the challenge.

Mathematically this a variation of a Minimum-cost flow problem or a flow network, we can formulate this as optimization problem as follows:

Let G = (V , E ) be a directed graph with a source vertex s ∈ V and a sink vertex t ∈ V. Each edge (u, v) ∈ E has a positive capacity c(u, v), a flow value f(u, v), and a cost a(u, v). The cost of sending flow along an edge (u, v) is given by f(u, v) * a(u, v). The objective is to find the minimum-cost flow of a specified amount d from source s to sink t, subject to the capacity constraints on the edges. The goal is to minimize the total cost of the flow over all edges :

To solve this we’ll use the simplex algorithm in networkx.

The simplex algorithm is a popular algorithm for solving linear programming problems, which are optimization problems where the objective function and the constraints are linear. The algorithm works by moving from one vertex of the feasible region (the set of points that satisfy the constraints) to another, along the edges of the polytope (the geometric shape formed by the constraints), until it reaches the optimal solution or finds that the problem is unbounded12

The algorithm involves adding slack variables to convert the inequalities into equalities, and then performing a series of operations called pivoting, which exchange a basic variable (one that appears in the left-hand side of an equality) with a non-basic variable (one that appears in the right-hand side or the objective function), while maintaining feasibility and improving the objective value34

1: Simplex algorithm — Wikipedia 2: Simplex algorithm — Cornell University Computational Optimization Open Textbook — Optimization Wiki 3: 1 The Simplex Method — Department of Computer Science 4: Simplex Method for Solution of L.P.P (With Examples) | Operation Research

Okay, so first we’ll have to redefine our graph, and create new data.

The simplex solver in networkx only accepts directed graphs, so we’ll create a directed graph representing our cities and constraints, we’ll assume we have a total start inventory of 1000, and to simplfy data generation we’ll assign random demand per city, then we’ll use the networkx nx.network_simplex(G) method to find a solution.

Also note that in this approach, we no longer assume we only have one depot in Paris and one delivery truck going across all cities, but rather every city with a surplus can be seen as a depot and can send a small truck (or multiple trucks) to deliver goods and fulfill some of the demand of the others cities. In a futur article we shall explore how to further tighten these constraints , for example force one delivery truck to go around , and have a sort of sequential route, that is for example we can’t have multiple outbound routes from one node at the same time , the truck must leave node i to node j , before returning to node i for resupply for example, or we can add another constraint on the capacity of the truck itself, but we can also assume that in this scenario the capacities of the roads actually represent the capacities of delivery trucks running on them.

import random
import networkx as nx
import matplotlib.pyplot as plt

# Set the number of nodes (cities) and maximum possible routes between two nodes
num_nodes = 10
max_routes = 3
required_flow = 1000  # Total flow across the graph

# Define the city data
city_data = {
    'city': ['Paris', 'Marseille', 'Lyon', 'Toulouse', 'Nice', 'Nantes', 'Montpellier', 'Strasbourg', 'Bordeaux', 'Lille'],
    'lat': [48.8567, 43.2964, 45.7600, 43.6045, 43.7034, 47.2181, 43.6119, 48.5833, 44.8400, 50.6278],
    'lng': [2.3522, 5.3700, 4.8400, 1.4440, 7.2663, -1.5528, 3.8772, 7.7458, -0.5800, 3.0583]
}

# Create a random graph with num_nodes nodes and random edges
G = nx.DiGraph()  # Convert to directed graph

# Add nodes to the graph with demand attributes
total_demand = 0
# Set the source node (Paris) demand to a negative value

for city in city_data['city']:
    if city_data['city'] != 'Paris':
      demand = random.randint(-required_flow, required_flow)
      G.add_node(city, demand=demand)
      total_demand += demand
    else:
     demand = -required_flow
     G.add_node(city, demand=demand)
     total_demand += demand  
# Balance the total demand to sum to 0
G.nodes[city_data['city'][-1]]['demand'] -= total_demand


# Generate random edges between nodes
for i in range(num_nodes):
    for j in range(i + 1, num_nodes):
        if random.random() < 0.5:
            num_routes = random.randint(1, max_routes)
            for _ in range(num_routes):
                distance = random.randint(50, 500)  # Random distance between 50 to 500 km
                cost = random.randint(10, 100)  # Random cost between 10 to 100 Euros
                capacity = random.randint(100, 500)  # Random capacity between 100 to 500 units
                G.add_edge(city_data['city'][i], city_data['city'][j], distance=distance, cost=cost, capacity=capacity)
                G.add_edge(city_data['city'][j], city_data['city'][i], distance=distance, cost=cost, capacity=capacity)  # Reverse direction

# Solve the problem using Minimum Cost Flow algorithm with required flow constraint
flow_cost, flow_dict = nx.network_simplex(G)

Lets us the found solution now :

Optimal Route: [(‘Paris’, ‘Nantes’), (‘Paris’, ‘Strasbourg’), (‘Paris’, ‘Lille’), (‘Marseille’, ‘Strasbourg’), (‘Marseille’, ‘Lille’), (‘Lyon’, ‘Paris’), (‘Lyon’, ‘Marseille’), (‘Lyon’, ‘Lille’), (‘Toulouse’, ‘Marseille’), (‘Toulouse’, ‘Bordeaux’), (‘Nice’, ‘Marseille’), (‘Nice’, ‘Bordeaux’), (‘Nantes’, ‘Toulouse’), (‘Nantes’, ‘Bordeaux’), (‘Montpellier’, ‘Lyon’), (‘Montpellier’, ‘Lille’), (‘Bordeaux’, ‘Marseille’), (‘Bordeaux’, ‘Lille’)]

Total Distance: 4895 km

Total Cost: 5763 Euros

Total Capacity: 5697 units

Total flow achieved: 711 units

That’s it for now , Next time we’ll see how we can build an RL agent to try and solve this complex problem , as the environment is now way more complexe, for instance the set of actions is no longer just where to go next but also how much can you carry with you whilst respecting all constraints and also trying to minimize the total cost.

I hope you find this article useful for you :)