Dynamic routing example

Dynamic_routing_example.ipynb Open In Colab Kaggle Gradient Open In SageMaker Studio Lab Hits

Description: Example of interactive optimization with GUI using AMPL and Google Maps

Tags: amplpy,gui

Notebook author: Christian Valente <christian.valente@gmail.com>

# Install dependencies
%pip install -q amplpy pandas

# Google Colab & Kaggle integration
from amplpy import AMPL, ampl_notebook

ampl = ampl_notebook(
    modules=["highs"],  # modules to install
    license_uuid="default",  # license to use
)  # instantiate AMPL object and register magics

Introduction

This example shows how to create a simple GUI for a routing problem that allows to:

  • visualize inputs and solutions on a map

  • dynamically add points by clicking on the map itself

  • show each step of the row generation methodology applied to solve it

Insert your Google Maps API key

Follow the instructions to generate a Google Maps key for development, then insert it in the placeholder below.

API_KEY = ""  # Fill in with your API key

Install and import dependencies

if not API_KEY:
    raise ValueError("Please insert your Google maps API key")
%pip install -q gmaps
import gmaps

gmaps.configure(api_key=API_KEY)
# Activate custom widgets manager
from google.colab import output

output.enable_custom_widget_manager()
# Import display and ipywidgets
from IPython.display import display
import ipywidgets as widgets

# Import needed data/maths libs
import pandas as pd
import math
import collections
import sys
# collections.Iterable after python 3.9
if sys.version_info >= (3, 10):
    collections.Iterable = collections.abc.Iterable

AMPL model: TSP model with subtour elimination

The standard integer programming formulation of the Travelling Salesman Problem is as follows:

Given a set of cities $S={1,..,n}$ and distances $c_{ij}$ ($i,j \in S$), we define:

  • $x_{ij} \in {0,1} \forall i \in S, j \in S, i \ne j $ a variable denoting if path from $i$ to $j$ or from $j$ to $i$ is part of ours solution.

Our aim is to minimize the total tour length:

  • $\sum_{i=1}^{n} \sum_{j=1,j \neq i}^{n} c_{ij}x{ij}$

while visiting all cities once, which is enforced with these two constraints:

  • $\sum_{i=1,i \neq j} x_{ij} \forall j \in S$

and

  • $\sum_{j=1,j \neq i} x_{ij} \forall i \in S$

This formulation notably allows for subtours that can be eliminated with a suitable subtour-elimination constraint; the number of possible subtours grows exponentially with the number of cities, so one approach is to eliminate only the subtours that result from actually obtained solutions.

The following AMPL model expresses the problem above with a slight semplification: we suppose the distance/cost matrix is symmetrical (aka the distances between each pair of cities is the same both ways), so we can consider only half the possible pairs for $X$. So we define:

  • $x_{ij} \in {0,1} \forall i \in S, j \in S, i > j $

and we change the formulation accordingly.

%%ampl_eval
# List of nodes
set NODES ordered;
# Longitude and latitude of each node; not used in the model,
# just for data retrieval
param longitude {NODES};
param latitude {NODES};
# Set of valid city pairs (segments)
set PAIRS := {i in NODES, j in NODES: ord(i) < ord(j)};
# Distance of each valid route
param distance {(i,j) in PAIRS};

# If set to 1, the relative segment is part of the solution
var X{PAIRS} binary;
# Minimize the total tour length
minimize Tour_Length: sum {(i,j) in PAIRS} distance[i,j] * X[i,j];

# Each node must have two segments active, showing that we get in
# and out of it
subject to Visit_All {i in NODES}:
   sum {(i, j) in PAIRS} X[i,j] + sum {(j, i) in PAIRS} X[j,i] = 2;

# The following are used in the subtour elimination procedure
param nSubtours >= 0 integer, default 0;
set SUB {1..nSubtours} within NODES;
# Subtour elimination constraint
# For each stored subtour, make sure that there is at least one segment
# leading into the subtour itself (so that it is not a subtour any longer)
subject to Subtour_Elimination {k in 1..nSubtours}:
  sum {i in SUB[k], j in NODES diff SUB[k]} 
      if (i, j) in PAIRS then X[i, j] else X[j, i] >= 2;

AMPL solution method

To solve the TSP model, we need to make sure it has no subtours; since enumerating them all beforehand is basically impossible, we choose a reactive approach. This function solves the model (with the current set of subtours cut out) and checks if the obtained solution has any subtour. If it does, it displays it and return False, otherwise the solution is a valid path, simply return True.

def solve_ampl(map) -> bool:
    nSubtoursParam = ampl.get_parameter("nSubtours")
    SubtoursSet = ampl.get_set("SUB")
    allsubtours = [
        ampl.get_set("SUB")[i].to_list() for i in range(1, int(nSubtoursParam.value()))
    ]

    # Solve the model
    ampl.solve()

    # Get solution
    print("Solved, getting solution")
    length = ampl.get_objective("Tour_Length").value()
    print(ampl.get_nodes())
    print(ampl.get_arcs())

    subtours = find_subtours(ampl.get_arcs(), ampl.get_nodes())
    # If we have only one tour, the solution is valid
    if len(subtours) <= 1:
        print(f"\nFound only {len(subtours)} subtours, returning")
        map.log(f"Length={length}. No subtour found, shortest valid path.")
        return True

    # Else:
    print(f"\nFound {len(subtours)} subtours, plotting them and adding cuts")
    map.log(f"Length={length}. {len(subtours)} subtours, adding cuts")

    # Plot the subtours
    coords = ampl.get_coords()
    for s in subtours:
        vt = s.copy()
        vt.append(vt[0])
        vt = as_steps(vt)
        vt = [(coords[i], coords[j]) for i, j in vt]
        map.add_arcs(vt)

    # Add the current tours to the list of subtours
    allsubtours.extend(subtours)
    # And add those to the constraints by assigning the values to
    # the parameter and the set
    nSubtoursParam.set(len(allsubtours))
    for i, tour in enumerate(allsubtours):
        print(f"Adding cut: {tour}")
        SubtoursSet[i + 1].set_values(tour)
    return False

“Glue” functions

These functions help glueing the GUI and python data processing part with the underlying AMPL object.

SOLVER = "highs"
from amplpy import DataFrame as aDF


def data_2_AMPL(coordinates: pd.DataFrame, distances: pd.DataFrame):
    # Convert the dataframes to the specific dataframes expected by AMPL and assings
    # those values to the set of NODES, the params longitude and latitude, and the parameter distances
    ampl.eval("reset data;")
    # Convert to amplpy dataframes
    coordinates = aDF.from_pandas(coordinates, index_names=["POINT"])
    distances = aDF.from_dict(
        distances, index_names=["orig", "dest"], column_names=["distance"]
    )
    ampl.set_data(coordinates, set_name="NODES")
    ampl.set_data(distances)


def solve(map: MyMap, alreadySolving: bool):
    # Main solve procedure; if we're not in the middle of a solution process,
    # reads the data from the map and pass it to AMPL then solve it.
    # Otherwise just continue the ongoing solution process
    if not alreadySolving:
        map.set_solving_status("solving")
        nodes = map.get_points()
        if len(nodes) <= 2:
            map.log("Trivial problem, add more points to solve.")
            map.set_solving_status("finished")
            return
        print(f"Solving for {len(nodes)} points")
        (coords, distances) = define_data(nodes)
        data_2_AMPL(coords, distances)
        ampl.option["solver"] = SOLVER
    continue_solve(map)


def continue_solve(map: MyMap):
    # Continue the ongoing solution process
    # Clear the previous solution from the map
    map.clean(markers=False)
    # Call AMPL to solve the current problem instance
    status = solve_ampl(map)
    # Status==True if the solution has no subtours
    if status:
        # If the solution is valid (no subtours), "inform" the map object and
        # display the solution
        map.set_solving_status("finished")
        arcs = ampl.get_arcs()
        coords = ampl.get_coords()
        tours = [(coords[i], coords[j]) for i, j in arcs]
        print(f"Solution has {len(tours)} arcs;")
        print(arcs)
        map.add_arcs(tours)

Show the UI

These commands load a sample dataset and add it to the map

# Create the GUI and use the solve function defined here to actually solve
# the model
m = MyMap(solve)
# Get sample data
(coord, dist) = sample_data()
# Add it to the map
m.add_points(coord[["latitude", "longitude"]].to_records(index=False))
# Show the GUI
m.render()

previous

Dual-Donor Organ Exchange problem

next

Efficient Frontier with Google Sheets