Local Search

← Greedy Algorithms, part 2

Parking a car →

Back to the IO book page

---

This page is the third 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

We need the packages already installed at the beginning of the previous tutorials.

Relocation

Let us now take a greedy solution as seen in the previous notebook and try to improve it by making small changes to it.

Many possible “small changes” have been proposed in the literature, and the pyvrp package provides lots of pre-packaged moves.

In this tutorial, however, we will build and test one such move from scratch.

To start with a simple one, consider the so-called “relocation” move: remove a client from a route and reassign it to a new position in the same route or in another one.

A few considerations:

• We are interested in feasible solutions: while removing a client does not violate constraints (other than making the solution incomplete), the addition to other routes might be impossible due to the vehicle’s capacity constraint.
• If we also allow the creation of new routes with the detached client alone, served from an arbitrary depot, then every feasible solution can be reached from any other by repeated applications of this move.
• This is very important, because we don’t want to cut out any solution because of hte choice of the starting point. However, guiding our moves toward the best solution can be hard.

We will build our local searcher step by step, by implementing a series of helper functions.

Technical note: building a class step by step

For reasons that will become apparent as our work proceeds (mainly related to code reusability), we want to collect all the solver code into a class. However, we want to build it gradually, and keep interleaving our code blocks with text as long as we write it.

Therefore, let us define a very small class with just the constructor code and a few attributes. As soon as we need to insert a new atttribute or method, we shall do it dynamically with the setattr() builtin function.

We create the Searcher class and initialize it with the model (the problem that we want to solve):

from typing import Optional
import numpy as np
import pyvrp

class Searcher:
    # Some elements in the problem will need to be retrieved again and again,
    # so let us make them available directly as class attributes
    data: pyvrp.ProblemData
    num_clients: int
    num_depots: int
    distance_matrix: np.ndarray[int]
    clients: list[pyvrp.Client]
    depots: list[pyvrp.Depot]
    vehicle_types: list[pyvrp.VehicleType]
    # All local search techniques maintain a state, consisting (at least)
    # of the current solution and the best solution found so far.
    current_solution: Optional[pyvrp.Solution]
    current_cost: Optional[int]
    best_solution: Optional[pyvrp.Solution]
    best_cost: Optional[int]
    # Constructor: take the model and extract the problem data
    def __init__(self, model: pyvrp.Model) -> None:
        self.data = model.data()
        self.num_clients = self.data.num_clients
        self.num_depots = self.data.num_depots
        self.distance_matrix = self.data.distance_matrix(0)
        self.clients = self.data.clients()
        self.depots = self.data.depots()
        self.vehicle_types = self.data.vehicle_types()
        # The initial step will be provided later
        self.current_solution = self.best_solution = None
        self.current_cost = self.best_cost = None

As we want to test our local search method from many possible initial solutions (starting from the greedy solution introduced in the previous tutorial), let us prepare a convenience method to “restart” the same searcher from an initial solution.

Arguments:

• self: reference to the Searcher object: remember that this function will be dynamically added as a method to the Searcher class.
• solution: the initial solution.

def restart(self: Searcher, solution: pyvrp.Solution) -> None:
    self.current_solution = self.best_solution = solution
    self.current_cost = self.best_cost = solution.distance_cost()

Let’s dynamically add this function as method to the Searcher class. The builtin function setattr() requires three arguments:

• the class reference
• name of the method (which we can access by the __name__ attribute of the function itself)
• the function reference.

setattr(
    Searcher,
    restart.__name__,
    restart
)

Evaluating a move

First of all, we need to be able to evaluate the cost reduction associated to a given move.

Let us introduce a function that evaluates the cost change when removing a specific client from a given position in a route.

Arguments:

• self: reference to the LocalSearch object: remember that this function will be dynamically added as a method to the Searcher class.
• route: the route from which we are going to remove the client.
• position_in_route: which position in the order of visit the client is currently occupying (0 being the first visited client).

Return value:

• the difference in cost (travelled distance for the route). Because of triangular inequality, the difference will always be negative (the cost after the removal is less than the cost before).

Observe that we are not evaluating a full move here, just the initial part (removing the node): we will also need to add back the client elsewhere.

def removal_cost(self: Searcher, route: pyvrp.Route, position_in_route: int) -> int:
    # Retrieve the array of visits from the route
    visits = route.visits()
    # Get the location (index in the distance matrix)
    # of the client that we are going to remove
    location = visits[position_in_route]
    # To evaluate the change in cost we need to know the location visited before
    # the removed client. Observe that, if the client was in the first (0)
    # position, then the previous location was the start depot.
    previous_location = \
        visits[position_in_route-1] if position_in_route > 0 \
        else route.start_depot()
    # Same for the location visited immediately after.
    next_location = \
        visits[position_in_route+1] if position_in_route < len(visits)-1 \
        else route.end_depot()
    # The change in cost is the distance between the previous and next location,
    # minus the distances of the two locations from the removed client:
    return self.distance_matrix[previous_location, next_location] - (
        self.distance_matrix[previous_location,location]
        + self.distance_matrix[location, next_location]
    )
setattr(Searcher, removal_cost.__name__, removal_cost)

Let us now evaluate the cost increase when we add a client to an existing route.

First of all, we need to check if the insertion of the new client still keeps the route feasible (the sum of the clients’ deliveries must not exceed the vehicle’s capacity).

Argument:

• self
• route
• client_index: the index of the client we want to add (to retrieve its delivery requirement).

Return value:

• True if the client can be added, False otherwise.

We don’t need to know the actual position in the route, because we only care about there being enough delivery for everyone.

def can_insert(self: Searcher, route: pyvrp.Route, client_index: int) -> bool:
    # Get the vehicle type, to know its capacity
    vehicle = self.vehicle_types[route.vehicle_type()]
    # Total delivery required by the route
    route_delivery = route.delivery()[0]
    # Additional delivery requested by the new client
    client_delivery = self.clients[client_index].delivery[0]
    # Check and return
    return route_delivery + client_delivery <= vehicle.capacity[0]
setattr(Searcher, can_insert.__name__, can_insert)

The following function evaluates the cost of adding the client to an existing route.

Arguments:

• self
• route
• client_index: the client to be added
• position_in_route

Return value:

• the difference in cost. None if the addition makes the rute infeasible (excessive delivery requirements wrt the vehicle’s capacity)

def insertion_cost(
    self: Searcher,
    route: pyvrp.Route,
    client_index: int,
    position_in_route: int
) -> Optional[int]:
    # if infeasible, return None
    if not self.can_insert(route, client_index): return None
    # get the sequence of client locations from the route
    visits = route.visits()
    # compute the new client location from its index
    client_location = self.num_depots + client_index
    # identify the location before the insertion point.
    # If the insertion point is 0, then use the start depot
    previous_location = \
        visits[position_in_route-1] if position_in_route > 0 \
        else route.start_depot()
    # same for the location at the insertion point. If the position is beyond
    # the list size, then the next location is the end depot.
    next_location = \
        visits[position_in_route] if position_in_route < len(visits) \
        else route.end_depot()
    # Return the cost difference of adding the location.
    return self.distance_matrix[previous_location,client_location] + (
        self.distance_matrix[client_location,next_location]
        - self.distance_matrix[previous_location,next_location]
    )
setattr(Searcher, insertion_cost.__name__, insertion_cost)

Let us also consider the possibility of creating a new route for the client.

Arguments:

• self
• depot_index: the index of the depot (start and end) for the new route.
• client_index

Return value:

• the cost of the new route (depot to client to depot)

def new_route_cost(self: Searcher, depot_index: int, client_index: int) -> int:
    #compute the client location in the matrix
    # (the depot index and location are the same)
    client_location = self.num_depots + client_index
    # return the round trip cost. Although the matrix should be symmetric,
    # let us not rely on that.
    return \
        self.distance_matrix[depot_index,client_location] \
        + self.distance_matrix[client_location,depot_index]
setattr(Searcher, new_route_cost.__name__, new_route_cost)

Now we have all elements needed to compute the change in cost of a given move.

To simplify things in the following part, let us collect a whole move description in a named tuple:

from typing import NamedTuple

class RelocationMove(NamedTuple):
    # Route index and position from which the client must be removed
    from_route_index: int
    from_position: int
    # Route index and position to which the client must be moved
    to_route_index: Optional[int] = None
    to_position: Optional[int] = None
    # Alternatively, depot index for a new route with the moved client
    to_depot_index: Optional[int] = None

In the above named tuple, we must set a few conventions:

• if to_depot_index is None, then both to_route_index and to_position must be integers;
• if to_depot_index is an integer, then both to_route_index and to_position must be None;
• if from_route_index is equal to to_route_index (the client must be moved to a different place in the same route), then to_position refers to the position in the original route, before the client is removed.

We can finally define a unique function to evaluate the move (as applied to the current solution).

Arguments:

• self
• move: the move that we want to evaluate.

Return value:

• the change in cost, or None if the move leads to an infeasible route as seen in the above functions.

def move_cost(self: Searcher, move: RelocationMove) -> Optional[int]:
    # get the current routes
    routes = self.current_solution.routes()
    # get the route from which we want to remove the client
    from_route = routes[move.from_route_index]
    # from the route, get the visits array with the sequence of visited clients
    from_visits = from_route.visits()
    # get the client index (consider that visits containt the locations)
    client_index = from_visits[move.from_position] - self.num_depots
    # compute the cost change in removing the client
    cost_1 = self.removal_cost(from_route, move.from_position)
    # check if the move requires the creation of a new route
    if move.to_depot_index is not None:
        # if so, evaluate the cost increase by invoking the appropriate function
        cost_2 = self.new_route_cost(move.to_depot_index, client_index)
    else:
        # else, expect both destination route fields to be non empty
        assert move.to_route_index is not None and move.to_position is not None
        # compute the cost of inserting the client in the required position
        # of the new route
        cost_2 = self.insertion_cost(
            routes[move.to_route_index],
            client_index, move.to_position
        )
        # if the insertion is infeasible, return None
        if cost_2 is None: return None
    # Return the overall cost (removal + insertion)
    return cost_1 + cost_2
setattr(Searcher, move_cost.__name__, move_cost)

Performing a move

Up to now, we only evaluated costs without actually moving clients. A soon as a good move is found, however, we can proceed to actually perform it.

As before, let us divide the task into the removal and insertion subtask, starting with the removal.

Observe that we are treating routes as immutable objects: insertion and removal will consist in the creation of new routes.

Arguments (same as in removal_cost()):

• self
• route
• position_in_route

Return values:

• the new route without the removed node (None if the route doesn’t have any client after removal)
• the index of the removed client (to be used for adding the same client elsewhere)

def remove_client(
    self: Searcher,
    route: pyvrp.Route,
    position_in_route: int
) -> tuple[Optional[pyvrp.Route], int]:
    # Get the sequence of visited clients
    visits = route.visits()
    # Extract the client location from the visits array
    client_location = visits.pop(position_in_route)
    # Compute the client index from the location
    client_index = client_location - self.num_depots
    # Create the new route object, unless the visits array is empty
    new_route = \
        pyvrp.Route(self.data, visits, route.vehicle_type()) if len(visits) > 0 \
        else None
    # Return the new route and the removed client index
    return new_route, client_index
setattr(Searcher, remove_client.__name__, remove_client)

Next, let us create a function to insert the previously removed client onto an existing route.

Arguments (same as insertion_cost()):

• self
• route
• client_index
• position_in_route

Returned value:

• the modified route

def insert_client(
    self: Searcher,
    route: pyvrp.Route,
    client_index: int,
    position_in_route: int
) -> pyvrp.Route:
    # get the sequence of visited clients
    visits = route.visits()
    # insert the client location (obtained from its index) in the requested position
    visits.insert(position_in_route, self.num_depots + client_index)
    # create and return the new route
    return pyvrp.Route(self.data, visits, route.vehicle_type())
setattr(Searcher, insert_client.__name__, insert_client)

A new route can be created with a single client.

Arguments (same as in new_route_cost()):

• self
• depot_index: the index of the depot (start and end) for the new route.
• client_index

Return value:

• the new route (depot to client to depot)

def new_route(self: Searcher, depot_index: int, client_index: int) -> pyvrp.Route:
    # get the vehicle type associated to the selected depot
    vehicle_type_index = [
        i for i in range(self.num_depots)
        if self.vehicle_types[i].start_depot == depot_index
    ][0]
    # create and return the new route with the single client (location, not index)
    return pyvrp.Route(
        self.data,
        [self.num_depots + client_index],
        vehicle_type_index
    )
setattr(Searcher, new_route.__name__, new_route)

Finally, we can create a function that applies an arbitrary move to the current solution.

Arguments (same as move_cost()):

• self
• move: the move that we want to evaluate.

Note that we assume that the move is feasible, and we don’t check if it is not: we assume that the move has been selected based on the cost improvement, and the cost functions have already determined its feasibility.

def update_cost_and_best(self: Searcher) -> None:
    self.current_cost = self.current_solution.distance_cost()
    if self.best_cost is None or self.current_cost < self.best_cost:
        self.best_solution = self.current_solution
        self.best_cost = self.current_cost
setattr(Searcher, update_cost_and_best.__name__, update_cost_and_best)

def apply_move(self: Searcher, move: RelocationMove) -> None:
    # get the routes list from the current solution
    routes = self.current_solution.routes()
    # get the route from which the client must be removed
    from_route = routes[move.from_route_index]
    # Remove the client from the requested position
    # (the function returns a new route and the client's index);
    # store the new route in place of the former
    routes[move.from_route_index], client_index = \
        self.remove_client(from_route, move.from_position)
    # Check if the move requires the inertion into an existing route
    if move.to_route_index is not None:
        # In this case, the position should be valid too
        assert move.to_position is not None
        # get the destination route from its index
        to_route = routes[move.to_route_index]
        # the destination position should be reduced by one if the client
        # is being moved forward in the same route
        to_position = \
            move.to_position if move.to_route_index != move.from_route_index \
                or move.to_position < move.from_position \
            else move.to_position-1
        #  perform the insertion replacing the old route with the new one
        routes[move.to_route_index] = self.insert_client(
            to_route,
            client_index,
            to_position
        )
    else:
        # if the creation of a new route was required instead,
        # create the route and append it to the routes array
        routes.append(self.new_route(move.to_depot_index, client_index))
    # the source route might now be empty;
    # in this case, remove the None from the routes list
    if routes[move.from_route_index] is None: routes.pop(move.from_route_index)
    # create the new solution; update cost and best
    self.current_solution = pyvrp.Solution(self.data, routes)
    self.update_cost_and_best()
setattr(Searcher, apply_move.__name__, apply_move)

Performing the local search

Now that we have all functions needed to evaluate and perform the moves, let us proceed by writing a utility iterator to generate all valid moves from the current solution.

Arguments:

• self

Return values:

• an iterator that yields all available (feasible or not) moves

from typing import Iterator

def iterate_moves(self: Searcher) -> Iterator[RelocationMove]:
    # The array of all route lengths
    route_lengths = [len(r.visits()) for r in self.current_solution.routes()]
    # iterate all possible source routes
    for from_route_index, from_route_length in enumerate(route_lengths):
        # for every source route, iterate
        # all possible source positions within the route
        for from_position in range(from_route_length):
            # for every source route and position, iterate all destination routes
            for to_route_index, to_route_length in enumerate(route_lengths):
                # destination positions must include
                # one more position for inserting after the last position
                for to_position in range(to_route_length+1):
                    # exclude a destination that would force us to add the client
                    # to the same position from which it was removed
                    if to_route_index != from_route_index \
                            or to_position < from_position \
                            or to_position > from_position + 1:
                        # yield the move from source to destination route and position
                        yield RelocationMove(
                            from_route_index=from_route_index,
                            from_position=from_position,
                            to_route_index=to_route_index,
                            to_position=to_position
                        )
            # likewise, for every source route and position consider the creation of
            # a singleton route starting and ending from each depot
            for to_depot_index in range(self.num_depots):
                yield RelocationMove(
                    from_route_index=from_route_index,
                    from_position=from_position,
                    to_depot_index=to_depot_index
                )
setattr(Searcher, iterate_moves.__name__, iterate_moves)

A single step in our search consists of finding an improving, feasible move and applying it.

Arguments:

• self
• current_solution
• first_improvement: True if the search for a move should end at the first improving one, or proceed with aanalyzing all moves until the best one is found.

Return value:

• True if there was an improvement, otherwise False.

The first_improvement argument lets us explore the two basic strategies in local search:

• First improvement: choose the first improving feasible move that you find; of course, to avoid bias it is necessary to explore moves in random order.
• Best move: explore all feasible moves and choose the one that offers the best cost reduction; while randomness is less important, it can help in cases where ties are expected.

import random

def search_step(self: Searcher, *, first_improvement: bool) -> bool:
    # we must remember the best improvement found so far,
    # and the associated move (initially none)
    best_cost_delta: Optional[int] = None
    best_move: Optional[RelocationMove] = None
    # let us iterate through all possible moves
    # from the current solution in random order
    move_list = [m for m in self.iterate_moves()]
    random.shuffle(move_list)
    for move in move_list:
        # compute the variation in cost due to the move
        cost_delta = self.move_cost(move)
        # if the move is infeasible (i.e., the relocated client
        # cannot be accepted in the destination route because of
        # delivery constraints) or the cost increases (we are
        # only interested in negative deltas), then skip to the next move
        if cost_delta is None or cost_delta >= 0: continue
        # if this is the dirst acceptable improvement or it is better
        # (more negative) that the previous best
        if best_cost_delta is None or cost_delta < best_cost_delta:
            # remember the new improvement and move.
            best_cost_delta = cost_delta
            best_move = move
            # if we aim for the first improving move, let us stop the search here.
            if first_improvement: break
    # if no improving move has been found then return None
    if best_move is None: return False
    # otherwise, return the new solution after the application of the move
    self.apply_move(best_move)
    return True
setattr(Searcher, search_step.__name__, search_step)

Our local search heuristic will consist of the repeated application of the search step until no improvement is possible.

Let us define a named tuple to hold our search history.

class SearchRecord(NamedTuple):
    step: int
    time: float
    cost: int
SearchHistory = list[SearchRecord]

The local search procedure starts from the current solution (which must have been preser by calling restart()) and performs a series of improving steps until a local minimum is found or the maximum number of steps is performed.

Arguments:

• self
• first_improvement: the move selection strategy
• time_limit: if not None, if a minimum is not found before, maximum epoch before ending the search

Return values:

• best (i.e., last) solution found
• search history

# import the unlimited iterator (to iterate over unlimited restarts)
from itertools import count
# we want to keep track of time
from time import time

def local_search(
        self: Searcher,
        *,
        first_improvement: bool,
        time_limit: Optional[float] = None,
        update_callback: Optional[Callable[[pyvrp.Solution], None]] = None
    ) -> tuple[pyvrp.Solution,SearchHistory]:
    # initialize the search history with info from the initial solution
    # and the initial time of search
    search_history: SearchHistory = [SearchRecord(0, time(), self.current_cost)]
    # count from 1 upward
    for step in count(1):
        # Perform a search step
        improvement = self.search_step(first_improvement=first_improvement)
        # if no improving feasible move, then we are done
        if not improvement: break
        # record the new search state
        t = time()
        search_history.append(SearchRecord(step, t, self.current_cost))
        # invoke the callback, if it has been provided
        if update_callback is not None: update_callback(self.best_solution)
        # End the search if the time limit has been reached
        if time_limit is not None and t > time_limit: break
    # Return the last solution found and the search history
    return self.best_solution, search_history
setattr(Searcher, local_search.__name__, local_search)

Testing and comparing the local search heuristics

Let us now apply our local search function to improve the greedy solution from the previous tutorial.

First of all, let us try a test run with a visualization of the local steps.

We import the model and the solution_callback function from the first tutorial; remember that solution_callback displays the current solution in real time. We also import the greedy solution from the second tutorial as a starting point for the local search.

from Code.greedy_algorithms_1 import model, solution_callback
from Code.greedy_algorithms_2 import build_greedy_closest_depot_solution

Now we can finally instantiate our local search object, which we’ll use to run our tests:

searcher = Searcher(model)

Now we can create the greedy solution, run the local searcher on it and display the solution step by step, displaying the improvements. Let’s choose the best-move strategy for this visual check:

from matplotlib import pyplot as plt
initial_solution = build_greedy_closest_depot_solution(model.data())
searcher.restart(initial_solution)
searcher.local_search(
    first_improvement=False,
    update_callback=solution_callback
)
plt.close()

Observe how the local search iterations try to disentangle and balance the initial greedy solution; while this procedure does not yield the best possible result (i.e., it gets stuck in a local minimum - we will get to this in the next tutorial), the improvement is significant: let us measure it.

Let’s test the best move and first improvement strategies starting from the greedy solution developed in the first tutorial (the one that does not use closest-depot constraints):

from Code.greedy_algorithms_1 import build_greedy_solution
initial_solution = build_greedy_solution(model.data())
searcher.restart(initial_solution)
first_improvement_solution, first_improvement_history = \
    searcher.local_search(first_improvement=True)
searcher.restart(initial_solution)
best_move_solution, best_move_history = \
    searcher.local_search(first_improvement=False)
print('Initial cost:', initial_solution.distance_cost())
print('First improvement solution cost:', first_improvement_solution.distance_cost())
print('best move solution cost:', best_move_solution.distance_cost())

Initial cost: 1434922
First improvement solution cost: 1132036
best move solution cost: 1060825

Let us compare the evolution of the two heuristics in terms of time and number of moves:

import matplotlib.pyplot as plt

_, axes = plt.subplots(figsize=(13,5), ncols=2, sharey=True)
# First axis: times
t0 = first_improvement_history[0].time
axes[0].plot(
    [r.time-t0 for r in first_improvement_history],
    [r.cost for r in first_improvement_history],
    label='First improvement'
)
t0 = best_move_history[0].time
axes[0].plot(
    [r.time-t0 for r in best_move_history],
    [r.cost for r in best_move_history],
    label='Best move'
)
axes[0].set_xlabel('Time (s)')
axes[0].set_ylabel('Cost')
axes[0].grid()
_=axes[0].legend()
# Second axis: steps
axes[1].plot(
    [r.step for r in first_improvement_history],
    [r.cost for r in first_improvement_history]
)
axes[1].plot(
    [r.step for r in best_move_history],
    [r.cost for r in best_move_history]
)
axes[1].set_xlabel('Step')
axes[1].grid()

You can notice from the left plot how the best-move strategy takes a much longer time (it checks all possible improvements before deciding every move), but in the end yields a much better result (although still not optimal). From the right-hand plot we see that the overall number of steps is lower for the best-move strategy.

Let us directly compare the solution values:

import pyvrp.plotting
_, axes = plt.subplots(figsize=(13,4), ncols=3, sharey=True)
pyvrp.plotting.plot_solution(initial_solution, model.data(), ax=axes[0])
axes[0].set_title('Initial (greedy)')
_=axes[0].legend([])
pyvrp.plotting.plot_solution(first_improvement_solution, model.data(), ax=axes[1])
axes[1].set_title('Final (first improvement)')
_=axes[1].legend([])
pyvrp.plotting.plot_solution(best_move_solution, model.data(), ax=axes[2])
axes[2].set_title('Final (best move)')
_=axes[2].legend([])

Let us repeat the same procedure starting from the closest-depot partition of clients:

initial_solution = build_greedy_closest_depot_solution(model.data())
searcher.restart(initial_solution)
first_improvement_solution, first_improvement_history = \
    searcher.local_search(first_improvement=True)
searcher.restart(initial_solution)
best_move_solution, best_move_history = \
    searcher.local_search(first_improvement=False)
print('Initial cost:', initial_solution.distance())
print('First improvement solution cost:', first_improvement_solution.distance())
print('best move solution cost:', best_move_solution.distance())

Initial cost: 1172456
First improvement solution cost: 1018856
best move solution cost: 969294

Again, let us compare the evolution of the local search procedures in terms of time and number of moves:

_, axes = plt.subplots(figsize=(13,5), ncols=2, sharey=True)
# First axis: times
t0 = first_improvement_history[0].time
axes[0].plot(
    [r.time-t0 for r in first_improvement_history],
    [r.cost for r in first_improvement_history],
    label='First improvement'
)
t0 = best_move_history[0].time
axes[0].plot(
    [r.time-t0 for r in best_move_history],
    [r.cost for r in best_move_history],
    label='Best move'
)
axes[0].set_xlabel('Time (s)')
axes[0].set_ylabel('Cost')
axes[0].grid()
_=axes[0].legend()
# Second axis: steps
axes[1].plot(
    [r.step for r in first_improvement_history],
    [r.cost for r in first_improvement_history]
)
axes[1].plot(
    [r.step for r in best_move_history],
    [r.cost for r in best_move_history]
)
axes[1].set_xlabel('Step')
axes[1].grid()

Finally, here is the visual comparison between the final results of the two local search strategies:

_, axes = plt.subplots(figsize=(13,4), ncols=3)
pyvrp.plotting.plot_solution(initial_solution, model.data(), ax=axes[0])
axes[0].set_title('Initial (greedy)')
_=axes[0].legend([])
pyvrp.plotting.plot_solution(first_improvement_solution, model.data(), ax=axes[1])
axes[1].set_title('Final (first improvement)')
_=axes[1].legend([])
pyvrp.plotting.plot_solution(best_move_solution, model.data(), ax=axes[2])
axes[2].set_title('Final (best move)')
_=axes[2].legend([])

The comparison between the two strategies depends on the initial configuration and on the sequence of random moves. Thus, in order to perform better comparisons, we need to be able to repeat the local search move starting from different initial solutions. This will be the topic of our next tutorial on Repeated Local Search.