Greedy algorithms - Part 1

Greedy Algorithms, part 2 →
Back to the IO book page

This page is the first of a series of short tutorials on topics covered by the upcoming book Intelligent Optimization - Optimization plus Machine Learning for Reliable Artificial Intelligence by Roberto Battiti, Kevin Tierney and Mauro Brunato.

This tutorial is released under the Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International license (CC BY-NC-SA 4.0)

Prerequisites

In order to be able to experiment this notebook locally, you must have at least one of the following:

To set up all necessary packages in your environment, the following installation command only needs to be run once:

!pip install pyvrp matplotlib

What are greedy algorithms?

The Greedy method is one of the basic “building blocks” for the development of optimization schemes. In greedy algorithms “better an egg today than a hen tomorrow”.

A solution is built through a sequence of decisions, each step involving a small part of it. The algorithm always makes the choice that looks best at the moment. A greedy algorithm is therefore shortsighted: it can make commitments to certain choices too early, which may prevent it from finding the best overall solution later. Only for specific problems like the Minimum Spanning Tree is the greedy solution guaranteed to be optimal. In general, greedy tends to be a simple and useful method to get an initial solution.

The Vehicle Routing Problem (VRP)

In this tutorial we explain the greedy algorithm in the context of the Vehicle Routing Problem (VRP).

A VRP instance defines a distribution network and we are asked to plan it operations in order to minimize a cost. As a minimum, we are given:

  1. A set of depots, with optional capacities.

  2. A set of clients with specific demands (e.g., quantity of product they need to receive).

  3. A set of vehicles that can deliver the product, with optional capacities.

In order to define costs, the cost of traveling between depots and customers are given, either directly or by defining locations.

In the following we assume that depots and customers live in a Cartesian plane and that the cost for moving between them is given by their Euclidean distance. The vehicles obey a number of constraints:

  • they are assigned to specific depots,
  • they have a maximum capacity,
  • they must go back to the same depot they left from.

Another common constraint is:

  • every vehicle makes one trip, possibly serving multiple customers and moving back to the depot; once it reaches the depot, it cannot reload.

In order to simulate multiple trips from the same vehicle, we simply assign more vehicles to the depot.

Building a problem instance

First of all, let us define a few parameters:

# Number of depots serving the area
n_depots = 4
# Number of customers to be served
n_clients = 100
# Client demands are random integers, from 1 to a maximum
client_maximum_demand = 10
# Number of vehicles per depot: we set a large number to avoid running out
# (we want our instance to be feasible, no matter what)
n_vehicles_per_depot = n_clients
# Maximum amount of goods that can be loaded into a vehicle
vehicle_capacity = 100
# The pyvrp package provides a CP solver, we want to limit its runtime
maximum_runtime = 1 # seconds
# The RNG seed, for uniformity and reproducibility;
# feel free to change it to explore different instances
seed = 1

To manage our VRP instance we use the PyVRP package. For this tutorial’s purposes, we are interested in creating instances and solutions from scratch; however, the PyVRP package will be useful to check the correctness and completeness of our solutions and to “manage our expectations”.

import pyvrp

model = pyvrp.Model()

Now we add depots to the model. Every depot has \((x,y)\) coordinates inside the \([0,100]\times[0,100]\) rectangle; furthermore, every depot has a number of vehicles, each having a specific capacity:

import random
random.seed(seed)

for _ in range(n_depots):
    depot = model.add_depot(
        x=random.uniform(0, 100),
        y=random.uniform(0, 100)
    )
    model.add_vehicle_type(
        n_vehicles_per_depot,
        capacity=vehicle_capacity,
        start_depot=depot,
        end_depot=depot,
        reload_depots=[depot]
    )

We do the same for clients; each client has a position and a specific demand (the delivery parameter):

for coord in range(n_clients):
    model.add_client(
        x=random.uniform(0, 100),
        y=random.uniform(0, 100),
        delivery=random.randint(1, client_maximum_demand)
    )

To complete the problem definition, we set up travel distances between all model locations (depots and clients). Note that edges between depots are not required, but they do no harm.

Observe that the library requires distances to be integers; to try preserving as much detail as possible, we multiply all distances by 1000 before rounding.

import math
from itertools import permutations

# Iterate through all pairs of different locations:
for frm, to in permutations(model.locations, 2):
    if type(frm) == pyvrp.Depot and type(to) == pyvrp.Depot: continue
    dx = to.x - frm.x
    dy = to.y - frm.y
    distance = int(.5 + 1000 * math.sqrt(dx**2 + dy**2))
    model.add_edge(frm, to, distance)

Let us have a look at our instance. The plot_coordinates() function is fine, but it mixes the legend with the plot, therefore we declare a custom ax:

import matplotlib.pyplot as plt
import pyvrp.plotting

ax = plt.axes()
pyvrp.plotting.plot_coordinates(model.data(), ax=ax)
_=ax.legend(bbox_to_anchor=(1,1))

Invoking the package’s solver

Before we work on our solution, let us fire up the package’s internal solver to see how it works. The solve() method prints out a summary of its execution and returns an object containing the optimization trace by the Constraint Programming solver and the final results.

cp_solution = model.solve(stop=pyvrp.stop.MaxRuntime(maximum_runtime))
PyVRP v0.13.3

Solving an instance with:
    4 depots
    100 clients
    400 vehicles (4 vehicle types)

    Iters    Time |      Current OK    Candidate OK         Best OK

Search terminated in 1.00s after 1494 iterations.
Best-found solution has cost 897683.

Solution results
================
    # routes: 4
     # trips: 6
   # clients: 100
   objective: 897683
    distance: 897683
    duration: 0
# iterations: 1494
    run-time: 1.00 seconds

The evolution of the optimization process can be plotted out by a helper function:

pyvrp.plotting.plot_objectives(cp_solution)

We can also plot out the best solution found (as before, we create a custom axis to take the legend off the plot):

ax = plt.axes()
pyvrp.plotting.plot_solution(cp_solution.best, model.data(), ax=ax)
_=ax.legend(bbox_to_anchor=(1,1))

A greedy solution

Let us work on a custom solution, beginning with a greedy one. Here we assume that every depot has enough vehicles to serve all routes that we come up with.

Our greedy setup will work as follows, where we maintain a list of unserved clients (initially, all of them):

  • As long as there are unserved clients:
    • locate the closest (depot, unserved client) pair;
    • initiate a new route from the chosen depot, load the vehicle with the maximum amount of goods
    • while the vehicle is not empty:
      • locate the closest unserved client to the vehicle whose demand doesn’t exceed the vehicle’s current load
      • if no such clients are available: exit the inner loop
      • move the vehicle to the client (add he new client to the current route)
      • unload the required amount of goods
      • mark the client as served
    • add the depot as the final step of the route.

The pyvrp package define classes for clients, depots, vehicle types and routes, but many functions refer to them using their numeric position in the respective arrays. Moreover in route lists and in the distance matrix depots and clients share a common indexing scheme (depots fitst, clients later) and are referred to referred to as “locations”.

In the following code, we will use this naming convention:

  • an index is an integer referring to a specific object (i.e., the index of clients, depots, routes in their respective arrays in the problem object);
  • a location (short for “location index”) is the unified enumeration of all locations (i.e., the index of depots and clients inside the distance matrix).
    • For depots, index and location are the same; for clients, they differ by the number of depots.
  • a position refers to the order of a client within a route.

The following function locates the closest (depot, client) pair where the client hasn’t been visited yet.

  • Arguments:
    • data contains all instance data, in particular the matrix with all mutual distances. The matrix has as many rows and columns as the number of locations (depots plus clients), with depots appearing before clients.
    • unserved clients is the set of all (zero-based) client indexes that haven’t been included in a route yet.
  • Return value:
    • the zero-based index of the chosen depot, None if there arent’any.
from typing import Optional

def closest_depot_to_unserved_clients(
    data: pyvrp.ProblemData,
    unserved_clients:set[int]
) -> Optional[int]:
    # Extract the distance matrix from the problem instance
    distance_matrix = data.distance_matrix(0)
    # Remember the smallest distance
    best_distance: Optional[int] = None
    depot_index: Optional[int] = None
    # Iterate on all depots
    for i in range(data.num_depots):
        # Iterate on still unassigned clients
        for j in unserved_clients:
            # get the distance; observe that the location index of
            # client j in the distance matrix is data.num_depots+j
            d = distance_matrix[i, data.num_depots + j]
            # if this is the first distance or the new distance
            # is better, remember it
            if best_distance is None or d < best_distance:
                best_distance = d
                depot_index = i
    # return the best distance (or none)
    return depot_index

A similar function to return the next step in the route (if any).

  • Arguments:
    • data, unserved_clients: same as above;
    • residual_capacity: the quantity of goods currently in the vehicle (there’s no point in visiting a customer if the vehicle doesn’t have enough residual capacity to fulfill the delivery);
    • current_location_index: index of the current location (zero based, depots first).
  • Return value:
    • zero-based index of the next client
def closest_client(
            data: pyvrp.ProblemData,
            unserved_clients: set[int],
            residual_capacity: int,
            current_location_index: int
        ) -> Optional[int]:
    distance_matrix = data.distance_matrix(0)
    best_distance = None
    client_index = None
    for j in unserved_clients:
        if data.clients()[j].delivery[0] <= residual_capacity:
            d = distance_matrix[current_location_index, data.num_depots+j]
            if best_distance is None or d < best_distance:
                best_distance = d
                client_index = j
    return client_index

A route can be greedily built by selecting a depot, then extending a list of visited nodes.

  • Arguments:
    • data, unserved_clients: as before;
    • depot_index: the index of the depot that manages the route;
    • added_client_callback: an optional function to be called every time a new client is added to the current route; the callback will be used to animate the greedy construction - but you can safely ignore its existence.
  • Returns:
    • the Route object.

Observe that, contrary to the previous functions, this function actually modifies the unserved_clients set by removing clients as soon as they are visited and added to the route.

from typing import Callable

def build_greedy_route(
    data: pyvrp.ProblemData,
    unserved_clients: set[int],
    depot_index: int,
    added_client_callback: Optional[Callable[[pyvrp.Route],None]] = None
) -> pyvrp.Route:
    # identify the (only) vehicle type associated to the chosen depot
    vehicle_type_index = [
        i for i in range(data.num_vehicle_types)
        if data.vehicle_types()[i].start_depot == depot_index
    ][0]
    vehicle = data.vehicle_types()[vehicle_type_index]
    # initialize the list of visited client indexes
    visited_clients = []
    # "load" the vehicle to full capacity
    residual_load = vehicle.capacity[0]
    # Index of the current location,
    # to be updated as the greedy construction proceeds
    current_location_index = depot_index
    # Exit condisions will be checked within the loop
    while True:
        # Find the closest index to the current location
        # compatible with the current vehicle load
        client_index = closest_client(
            data,
            unserved_clients,
            residual_load,
            current_location_index
        )
        # If no client is found, it's time to go back home
        if client_index is None:
            # Embed the information in a Route object and return it
            return pyvrp.Route(data, visited_clients, vehicle_type_index)
        # otherwise, add the new location index to the visits list
        # (the index must include depots)
        visited_clients.append(data.num_depots + client_index)
        # move the current location to the new client
        # (location index includes depots)
        current_location_index = data.num_depots + client_index
        # remove the client's delivery from the current vehicle load
        residual_load -= data.clients()[client_index].delivery[0]
        # remove the client from the unserved client set
        unserved_clients.remove(client_index)
        # Now notify the callback that a new client has been added
        if added_client_callback is not None:
            added_client_callback(
                pyvrp.Route(data, visited_clients, vehicle_type_index)
            )
        # repeat until no more clients can be added

The greedy solution is built by repeatedly adding greedy routes until all clients are exhausted-

  • Arguments:
    • data the problem instance data;
    • update_callback: as in build_greedy_route, an optional function to be called every time a new client is added to the current route, used to animate the greedy construction - again, you can safely ignore its existence.
  • Returns:
    • the full solution
def build_greedy_solution(
    data: pyvrp.ProblemData,
    update_callback: Optional[Callable[[pyvrp.Solution], None]] = None
) -> pyvrp.Solution:
    # list of greedy routes,
    # used to build the solution object after it's filled up
    greedy_routes: list[pyvrp.Route] = []
    # initialize the set of unserved client indexes to the whole client set
    unserved_clients = set(range(data.num_clients))
    # repeat untile there are no more clients
    while len(unserved_clients) > 0:
        # find the depot with the closest unserved client
        depot_index = closest_depot_to_unserved_clients(data, unserved_clients)
        # since there are unserved clients, we are sure to find a depot
        assert depot_index is not None
        # get the next greedy route in terms of vehicle type
        # (encoding the chosen depot) and a list of visited client indexes
        route = build_greedy_route(
            data,
            unserved_clients,
            depot_index,
            # if the callback is defined, create a partial solution with the
            # current routes and invoke it
            None if update_callback is None
            else lambda partial_route:
                update_callback(
                    pyvrp.Solution(data, greedy_routes + [partial_route])
                )
        )
        # create the corresponding Route object
        # and append it to the list of routes
        greedy_routes.append(route)
    # Create the solution object and return it
    return pyvrp.Solution(data, greedy_routes)

Testing the algorithm

Let us actually build the solution. Before calling the greedy constructor, we define a function that plots the partial solution while it’s being built, using the callback parameters defined in the previous functions.

# import packages for real-time visualization
from IPython.display import display, clear_output
from time import sleep

def solution_callback(partial_solution: pyvrp.Solution) -> None:
    # clear the current output
    clear_output(wait=True)
    # forget the current plot
    plt.close()
    # plot the solution on the current axis
    ax = plt.axes()
    pyvrp.plotting.plot_solution(
        partial_solution, model.data(), ax=ax, plot_clients=True
    )
    # No legend
    ax.legend([])
    # Actually show the plot and wait a short time
    display(plt.gcf())
    sleep(.01)

# Now build the greedy solution
greedy_solution = build_greedy_solution(model.data(), solution_callback)
# Result already displayed
plt.close()

You may notice that from time to time the greedy construction seems to disregard a closer client in order to reach a farther one: this is not a bug - remember that the vehicle has a limited capacity and will only visit clients whose demand can be satisfied.

The algorithm performance

Let us see how much our greedy solution is worse than the CP optimum:

best_distance = cp_solution.best.distance()
greedy_distance = greedy_solution.distance()
improvement = (greedy_distance - best_distance) / best_distance
print(f'''Total distance = {greedy_solution.distance()}
CP optimum: {cp_solution.best.distance()}
Greedy is worse by {int(100*improvement+.5)}%''')
Total distance = 1434922
CP optimum: 897683
Greedy is worse by 60%

Observe that the solution is far from being optimal: while most routes start from clients in the vicinity of a depot, they tend to wander far from it, even going close to other depots, until the load of the vehicle is spent, followed by a long journey home:

ax = plt.axes()
pyvrp.plotting.plot_solution(greedy_solution, model.data(), ax=ax)
_=ax.legend(bbox_to_anchor=(1,1))

What next?

As we have just seen, there is quite a large gap between our simple greedy solution and the result obtained with the specialized Constraint Programming optimizer implemented in the pyvrp package.

In the upcoming tutorials we shall explore refinements of our algorithm and other improvements achieving better and better results.