38 KiB
TP1 - Branch & Bound, Cut-Generation for the TSP¶
1. The minilp python package¶
The minilp python package is a small python package that allows us to easily model (integer) linear program. The package comes with an interface to common linear programming solvers (cplex, scipy) but no integer linear programming solver.
The minilp module has been implemented for these practical sessions so you will not find any relevant documentation on Google.
To get help on the module and its components, you can use the built-in help function:
help(minilp)
import minilp
lp = minilp.problem("My first LP problem")
# Create two continuous variables within [0, 4]:
x1, x2 = lp.continuous_var_list(2, 0, 4)
# Add constraints:
lp.add_constraint(-3 * x1 + 4 * x2 <= 7)
lp.add_constraint(2 * x2 <= 5)
lp.add_constraint(6 * x1 + 4 * x2 <= 25)
lp.add_constraint(2 * x1 - x2 <= 6)
# Set the objective function:
lp.set_objective("max", x1 + 2 * x2)
# Solve the problem:
res = lp.lp_solve()
print(res)
print("x1 = {:.4f}, x2 = {:.4f}".format(res.get_value(x1), res.get_value(x2)))
The minilp package also allows you to modelise simple integer linear programs.
There are also binary_var and binary_var_list method to create binary variable (integer variable constrained to 0 and 1).
The _list methods returns standard python list object, so you can combine them to create multi-dimensional lists of
minilp variables.
import minilp
lp = minilp.problem("My first ILP problem")
# Create two integer variables within [0, 4]:
x1, x2 = lp.integer_var_list(2, 0, 4)
# Add constraints:
lp.add_constraint(-3 * x1 + 4 * x2 <= 7)
lp.add_constraint(2 * x2 <= 5)
lp.add_constraint(6 * x1 + 4 * x2 <= 25)
lp.add_constraint(2 * x1 - x2 <= 6)
# Set the objective function:
lp.maximize(x1 + 2 * x2)
While minilp allows you to model (mixed-)integer linear program, it does not provide a integer linear program solver — the lp_solve method will always solve the linear relaxation of the problem.
You can use a different solver by passing a different object to the lp_solve method. Available solvers are:
minilp.solvers.scipy— A solver based on thescipymodule, wildly available.minilp.solvers.docplex— A solver based on thedocplexmodule, which requires a valid CPLEX installation.
The default solver used is docplex if available, otherwize scipy.
res = lp.lp_solve()
print(res)
print("x1 = {:.4f}, x2 = {:.4f}".format(res.get_value(x1), res.get_value(x2)))
The minilp package allows you to modelise <=, >= or == (in)equalities. You can create linear expression by simply adding, substracting or multiplying values (int or float) and variables or existing expressions. You can use the standard python sum to sum a bunch of expressions or variables, and the minilp.dot function to compute the dot product of two vectors.
Exercice: Complete the following code to create a simple model for the knapsack problem.
Make your program as generic as possible, i.e., if N or K changes, you should not have to
modify the code of your problem.
import minilp
N = 5
p = [1, 4, 5, 3, 5] # profits
w = [3, 4, 3, 5, 9] # weights
K = 10 # capacity
assert N == len(w) and N == len(p)
# A simple knapsack
kp = minilp.problem("Simple knapsack")
# TODO: Create variables, add constraints and set the objective.
x = ...
# We can solve the linear relaxation:
res = kp.lp_solve()
print(res)
print(res.get_values(x))
2. Generic Branch & Bound¶
The purpose of the following section is to implement a generic branch-and-bound procedure based on a minilp.problem instance.
2.1. Node structure¶
Do not overthink this section!
The three methods you have to implement in this section (create_root_node, create_root_node and is_empty) are
all one-liner and can be implemented with a single return statement.
We are going to use a simple list to represent the current set of leaf nodes in the branch-and-bound tree. For each node, we are only going to store the list of separation constraints of the node (and not an actual problem).
For instance, assume a problem with variables $x_1$ and $x_2$, we two separation constraints, our node could be created with:
node = [x1 <= 2, x2 >= 4]
import minilp
from typing import List, Optional
# minilp.cons is the Python type representing minilp constraints:
Node = List[minilp.cons]
Question: What does the root node of the branching tree contain?
Exercice: Implement the function create_root_node() that creates the root node of the branching tree.
def create_root_node() -> Node:
"""
Creates the root node for a branch-and-bound procedure.
Returns:
The root node for a branch-and-bound procedure.
"""
return ... # TODO
For the sake of simplicity, we are going to process node in a first-in/first-out order.
Exercice: Implement the function get_next_node(nodes) that extract the next node to process from the tree.
You can check the list.pop method from python.
def extract_next_node(nodes: List[Node]) -> Node:
"""
Extracts the next node to process from the given list of nodes. The node
is removed from the given list before being returned.
Args:
nodes: Current list of nodes to extract the next node from.
Return:
The next node process.
"""
return ... # TODO
Exercice: Implement the function is_empty(nodes) that returns True if there are no more nodes to process in the list.
def is_empty(nodes: List[Node]):
"""
Checks if there are any nodes remaining to process in the given list.
Args:
- nodes: The list of nodes to check.
Returns:
True if there are no more nodes, False otherwise.
"""
return ... # TODO
2.2. Comparing minilp.result¶
The minilp.problem.lp_solve method returns an instance of minilp.result. The following method compare two minilp.result objects for a given
problem and returns True if the left one is better.
Don't forget to execute the cell even if you do not have to modify it in order to have access to the compare_solution function.
def compare_solutions(
problem: minilp.problem, l: minilp.result, r: minilp.result
) -> bool:
"""
Compares the two given solutions, returning True if the left
one is better than the right one for the current problem.
Args:
problem: The problem for which the solution are compared.
l, r: The two solutions (minilp.result) to compare.
Returns:
True if the left solution is better than the right one, or
if the right solution has no solution.
"""
if problem.isnan(r.objective):
return True
if problem.sense == "min":
return l.objective < r.objective
return l.objective > r.objective
2.3. Finding bound variable¶
In order to execute the branch-and-bound procedure, we must be able to find a variable to separate the problem on. For this notebook, we are always going to separate on the first non-integer variable.
You can access the list of variables in the problem with minilp.problem.variables and their values in the
solution with minilp.result.get_value or minilp.result.get_values.
There are variables in the problem that do not have to be integers in the solution (e.g., transport variables in the warehouse allocation problem). These variables should not be checked against when looking for a non-integer value.
You can retrieve the type of a variable (int or float) using minilp.expr.var.category.
Exercire: Implement the get_first_non_integral(problem, result) method that, given a problem and solution, returns
the first variable of the problem that is not integral in the result (and should be), or None if there is no such variable
In order to solve the linear relaxation of the problem, multiple matrix operations have to be performed. These
operations are numerically imprecise, thus it is common to find near-integral values (1.0000001 or 2.9999999).
Such values should be considered integral for the purpose of our algorithms, as long as their distance to the
nearest integer is less than eps (default to 1e-6).
def find_first_non_integral(
problem: minilp.problem, result: minilp.result, eps: float = 1e-6
) -> Optional[minilp.var]:
"""
Retrieves the first integer variable in the given problem
whose value is not integral in the given solution.
Args:
problem: The problem to find a variable from.
result: A solution of the problem.
eps: The maximum allowed distance to consider a variable integral.
Returns:
The first variable (minilp.expr.var) whose value is not integral,
or None if no such variable exists.
"""
... # TODO
2.4. Relaxation, separation, iteration¶
Exercice: Implement the iterate method below that performs a single iteration of the branch-and-bound algorithm, i.e., extract a node,
solve its relaxation, and then update the current best solution (return) or separate the problem.
You can use the minilp.problem.add_constraints and minilp.problem.del_constraints methods to add or remove constraints from
a minilp.problem.
Do not forget to return the minilp.result you found if it is your new best integer solution!
Again, do not overthink this section! While iterate is the core function of the branch-and-bound algorith, it is a
pretty simply and short function that can be implemented in about 10 lines.
import math
def iterate(
problem: minilp.problem,
current_best: minilp.result,
nodes: List[Node],
solver: minilp.solver = minilp.solvers.get_default_solver(),
eps: float = 1e-6,
) -> Optional[minilp.result]:
"""
Performs an iteration of the branch-and-bound algorithm.
Args:
problem: Problem to perform an iteration for.
current_best: Current best known feasible (integral) solution. Instance of minilp.result.
nodes: Current list of nodes to update.
solver: Solver to use to solve the linear relaxation.
eps: The maximum allowed distance to consider a variable integral.
Returns:
The new best solution if the current one was improved, or None if no new solution
was found.
"""
... # TODO
2.5. The branch-and-bound algorithm¶
The cell below defines the global branch_and_bound procedure.
import datetime as dt
def log_solution(
niterations: int, nodes: List[Node], res: minilp.result, new_best: bool = False
):
"""
Logs the given solution with time information.
Args:
niterations: Number of iterations.
nodes: List of nodes.
res: Solution (minilp.result) to log.
new_best: Indicates if this solution is new best solution (integer).
"""
print(
"{} {:5d} {:5d} {:9g}{}".format(
dt.datetime.now().strftime("%T"),
niterations,
len(nodes),
res.objective,
"*" if new_best else "",
)
)
def branch_and_bound(
problem: minilp.problem,
lp_solver: minilp.solver = minilp.solvers.get_default_solver(),
eps: float = 1e-6,
log_frequency: int = 10,
) -> minilp.result:
"""
Applies a branch-and-bound algorithm to solve the given problem.
Args:
problem: A minilp.problem instance corresponding to the problem to solve.
lp_solver: Solver to use to solve the linear relaxation.
eps: The maximum allowed distance to consider a variable integral.
log_frequency: Number of iterations between two log (not including exceptional log).
Returns:
A minilp.result containing the result of running the branch-and-bound
algorithm.
"""
print(
"B&B using {} to solve linear relaxation".format(lp_solver.__class__.__name__)
)
# Insert the first node in the list:
nodes = []
nodes.append(create_root_node())
# Current best result (unknown):
current_best = minilp.result()
# Counter for iterations:
nb_iterations = 0
while not is_empty(nodes):
# Increment counter and log.
if nb_iterations % log_frequency == 0:
log_solution(nb_iterations, nodes, current_best)
nb_iterations += 1
# Iterate:
new_best = iterate(problem, current_best, nodes, lp_solver, eps)
# Check if we have a new current best:
if new_best is not None:
current_best = new_best
log_solution(nb_iterations, nodes, current_best, True)
# Return the best solution found (if any).
return current_best
Exercice: Use the branch_and_bound method to solve the knapsack instance defined at the beginning of the notebook.
import minilp
N = 5
p = [1, 4, 5, 3, 5] # profits
w = [3, 4, 3, 5, 9] # weights
K = 10 # capacity
assert N == len(w) and N == len(p)
# A simple knapsack
kp = minilp.problem("Simple knapsack")
# TODO: Create variables, add constraints and set the objective.
x = kp.binary_var_list(N)
kp.add_constraint(kp.dot(x, w) <= K)
kp.maximize(kp.dot(x, p))
# We can solve the linear relaxation:
res = branch_and_bound(kp)
print(res)
print(res.get_values(x))
Exercice: Create other instances of the knapsack problem to reach the "limits" of your implementation — What is the largest instance you can solve in e.g. less than 5 seconds?
You can use the numpy.random module to generate arrays of random integer or floating point values.
import minilp
import numpy as np
N = 20
p = np.random.randint(50, 100, size=N)
w = np.random.randint(20, 40, size=N)
K = np.sum(w) // (1.5 + np.random.rand() * 1)
assert N == len(w) and N == len(p)
print("Knapsack problem with {} items and a capacity of {}.".format(N, K))
print(" Profits: {}".format(p))
print(" Weights: {}".format(w))
# A simple knapsack
kp = minilp.problem("Simple knapsack")
# TODO: Create variables, add constraints and set the objective.
x = kp.binary_var_list(N)
kp.add_constraint(kp.dot(x, w) <= K)
kp.maximize(kp.dot(x, p))
# We can solve the linear relaxation:
res = branch_and_bound(kp, log_frequency=5)
print(res)
print(res.get_values(x))
3. The Travelling Salesman Problemn (TSP)¶
Given a list of $n$ cities and the distances $c_{ij}$ between each pair of cities, you want to find the shortest circuit that visits each city exactly once and comes back to the first visited city.
The tsp.data packages contains grid of distances of various sizes (5, 6, 7, 8, 9, 10, 15, 17, 26, 42).
The goal of this section is to implement a cut-generation algorithm for the travelling salesman problem, using the minilp python package.
3.1. Creating a model for the TSP¶
3.1.1. Relaxation of the TSP¶
We call "TSP relax" the relaxation of the TSP problem that do not include constraints to eliminate subtours. The model is given below:
- $x_{ij}\in\{0, 1\}$ — Binary variable indicating if we go directly from city $i$ to city $j$.
$ \begin{align} \text{min.} \quad & \sum_{i = 1}^{n} \sum_{j=1}^{n} c_{ij}x_{ij} & \\ \text{s.t.} \quad & \sum_{\substack{j = 1\\ j \neq i}}^{n} x_{ij} = 1, & \forall i\in \left\{1,~\ldots,~n\right\}\\ & \sum_{\substack{i = 1\\i \neq j}}^{n} x_{ij} = 1, & \forall j\in \left\{1,~\ldots,~n\right\}\\ & x_{ij} \in\{0,1\}, & \forall i,j \in \{1,\ldots,n\} \end{align} $
Exercice: Create the tsp_relax function that creates a minilp.problem instance corresponding to
the relaxed TSP problem.
from typing import Sequence, Tuple
import minilp
def tsp_relax(
distances: Sequence[Sequence[float]], name: str = "TSP relax"
) -> Tuple[minilp.problem, Sequence[Sequence[minilp.var]]]:
"""
Create a 'relaxed' model for the TSP. A relaxed includes all the standard
constraints of the TPS, but not the subtours constraints.
Args:
distances: The matrix of distances (diagonal is zeros).
name: The name of the model.
Returns:
A tuple (model, vars) where model is the TSP model and vars is the matrix
of x variables.
"""
N = len(distances)
tsp = minilp.problem(name)
... # TODO
# Returns both the problem and the variables:
return tsp, x
You can visualize a minilp model by using repr(model) or simply writting the variable name at the end of a cell.
import tsp.data as data
tsp, x = tsp_relax(data.grid5)
tsp
3.1.2. MTZ formulation for the TSP¶
We call "TSP relax" the relaxation of the TSP problem that do not include constraints to eliminate subtours. The model is given below:
- $x_{ij}\in\{0, 1\}$ — Binary variable indicating if we go directly from city $i$ to city $j$.
- $u_{i} \in\{1,\ldots,n\}$ — Subtour elimination variables: $u_i$ is the position of city $i$ in the tour.
$ \begin{align} \text{min.} \quad & \sum_{i = 1}^{n} \sum_{j=1}^{n} c_{ij}x_{ij} & \\ \text{s.t.} \quad & \sum_{\substack{j = 1\\ j \neq i}}^{n} x_{ij} = 1, & \forall i\in \left\{1,~\ldots,~n\right\}\\ & \sum_{\substack{i = 1\\i \neq j}}^{n} x_{ij} = 1, & \forall j\in \left\{1,~\ldots,~n\right\}\\ & u_1 = 1 & \\ & 2 \leq u_i \leq n, & \forall i\in\{2,\ldots,n\} \\ & u_i - u_j +1 \leq (n-1)(1 - x_{ij}), &\forall i, j \in \{2,\ldots,n\}\\ & x_{ij} \in\{0,1\},\ u_{i}\in\mathbb{N} & \forall i,j \in \{1,\ldots,n\} \end{align} $
Exercice: Create the tsp_mtz function that creates a minilp.problem instance corresponding to
the MTZ formulation of the TSP.
Use the tsp_relax method to initialize the problem with the basic TSP constraints and retrieve the $x$ variables.
import minilp
def tsp_mtz(
distances: Sequence[Sequence[float]], name: str = "TSP MTZ"
) -> Tuple[minilp.problem, Tuple[Sequence[Sequence[minilp.var]], Sequence[minilp.var]]]:
"""
Create a MTZ model for the TSP.
Args:
distances: The matrix of distances (diagonal is zeros).
name: The name of the model.
Returns:
A tuple (model, (x, u)) where model is the TSP model and (x, u) are the
variables of the model (x is a 2D array, u is a 1D array).
"""
N = len(distances)
# Use tsp_relax to create the basic model:
tsp, x = tsp_relax(distances, name)
... # TODO
# Returns both the problem and the variables:
return tsp, (x, u)
3.1.3. Flow formulation for the TSP¶
We call "TSP relax" the relaxation of the TSP problem that do not include constraints to eliminate subtours. The model is given below:
- $x_{ij}\in\{0, 1\}$ — Binary variable indicating if we go directly from city $i$ to city $j$.
- $y_{ij} \in\mathbb{R}_*^+$ — Subtour elimination variables: $y_{ij}$ is the flow on arc $(i,j)$, each city produces one unit of flow.
$ \begin{align} \text{min.} \quad & \sum_{i = 1}^{n} \sum_{j=1}^{n} c_{ij}x_{ij} & \\ \text{s.t.} \quad & \sum_{\substack{j = 1\\ j \neq i}}^{n} x_{ij} = 1, & \forall i\in \left\{1,~\ldots,~n\right\}\\ & \sum_{\substack{i = 1\\i \neq j}}^{n} x_{ij} = 1, & \forall j\in \left\{1,~\ldots,~n\right\}\\ & \sum_{j=2}^{n} y_{1j} = 1 & \\ & \sum_{j=1}^{n} y_{ij} = \sum_{j=1}^{n} y_{ji} + 1, & \forall i\in\{2,\ldots,n\} \\ & y_{ij} \leq n x_{ij}, &\forall i, j \in \{1,\ldots,n\}\\ & x_{ij} \in\{0,1\}, y_{ij}\in\mathbb{R}_*^{+} & \forall i,j \in \{1,\ldots,n\} \end{align} $
Exercice: Create the tsp_flow function that creates a minilp.problem instance corresponding to
the flow formulation of the TSP.
Use the tsp_relax method to initialize the problem with the basic TSP constraints and retrieve the $x$ variables.
import minilp
def tsp_flow(
distances: Sequence[Sequence[float]], name: str = "TSP Flow"
) -> Tuple[
minilp.problem,
Tuple[Sequence[Sequence[minilp.var]], Sequence[Sequence[minilp.var]]],
]:
"""
Create a Flow model for the TSP.
Args:
distances: The matrix of distances (diagonal is zeros).
name: The name of the model.
Returns:
A tuple (model, (x, y)) where model is the TSP model and (x, y) are the
variables of the model (x and y are 2D arrays).
"""
N = len(distances)
# Use tsp_relax to create the basic model:
tsp, x = tsp_relax(distances, name)
... # TODO
# Returns both the problem and the variables:
return tsp, (x, y)
3.2. Solving small instances of the TSP problem¶
Exercice: Using the tsp_mtz, tsp_flow and branch_and_bound functions you implemented, solve the small TSP instances found in tsp.data.
Question: How large are the instances you are able to solve in a reasonable amount of time?
import tsp.data as data
distances = data.grid7
print("=== MTZ ===")
tsp, (x, u) = tsp_mtz(distances)
r = branch_and_bound(tsp)
print(r)
print("=== Flow ===")
tsp, (x, y) = tsp_flow(distances)
r = branch_and_bound(tsp, eps=1e-8, log_freqeuncy=2)
print(r)
3.3. Cut generation algorithm for the TSP¶
We are going to implement a cut-generation procedure using the following TSP formulation:
- $x_{ij}\in\{0, 1\}$ — Binary variable indicating if we go directly from city $i$ to city $j$.
$ \begin{align} \text{min.} \quad & \sum_{i = 1}^{n} \sum_{j=1}^{n} c_{ij}x_{ij} & \\ \text{s.t.} \quad & \sum_{\substack{j = 1\\ j \neq i}}^{n} x_{ij} = 1, & \forall i\in \left\{1,~\ldots,~n\right\}\label{tsp2:leave}\\ & \sum_{\substack{i = 1\\i \neq j}}^{n} x_{ij} = 1, & \forall j\in \left\{1,~\ldots,~n\right\}\label{tsp2:enter}\\ & \sum_{i\in S} \sum_{j \notin S} x_{ij} \geq 1 & \mathcal{S} \subset \{1, \ldots, n\},\ S \ne \emptyset \label{tsp4:subtour:0} \\ & x_{ij} \in\{0,1\}, & \forall i,j \in \{1,\ldots,n\} \end{align} $
Where $\mathcal{S}$ is the set of subtours in the graph. The generation procedure is as follow:
- We start with an empty set $\mathcal{S} = \emptyset$.
- We solve the problem (completely) using our
branch_and_boundprocedure. - We find all the subtours in the solution:
- If there is a single subtours, we have found the optimal solution.
- Otherwize, we add the subtours to the set $\mathcal{S}$ and we got back to 2.
3.3.1. Generating subtour constraints¶
Exercice: Implement the add_subtour_constraints method that, given a list of subtours (a list of list of integers, e.g., [1, 2, 4] is a subtour going through nodes 1, 2 and 4),
create subtour constraints and add them to the given problem.
def add_subtour_constraints(
tsp: minilp.problem,
x: Sequence[Sequence[minilp.var]],
subtours: Sequence[Sequence[int]],
) -> Sequence[minilp.cons]:
"""
Generates and adds subtours constraints for the given problem.
Args:
tsp: The current TSP problem.
x: The x variables (2D-array) of the TSP formulation.
subtours: The list of subtours. A subtour is a list of int containing the ID of the node.
Returns:
The generated constraints.
"""
... # TODO
# We can test the method by generated a relaxed-TSP and asking constraints for random subtours:
tsp, x = tsp_relax(data.grid5)
add_subtour_constraints(tsp, x, [[1, 2, 3], [0, 4]])
3.3.2. Finding subtours in a solution¶
Exercice: Implement the method find_subtours that given a 2D-array x of integer (0 or 1), returns a list of subtours in x.
def find_subtours(x: Sequence[Sequence[bool]]) -> Sequence[Sequence[int]]:
"""
Extracts subtours from the given 2D-array.
Args:
x: A two-dimensional array corresponding to the x variable in the TSP formulation, where
x[i][j] is 1 if arc (i, j) is used.
Returns:
A list of subtours, where each subtour is a list.
"""
... # TODO
# We can check the method by using a custom x array corresponding
# to two subtours (0 -> 2 -> 0) and (1 -> 4 -> 3 -> 1):
find_subtours(
[
[0, 0, 1, 0, 0], # 0 -> 2
[0, 0, 0, 0, 1], # 1 -> 4
[1, 0, 0, 0, 0], # 2 -> 0
[0, 1, 0, 0, 0], # 3 -> 1
[0, 0, 0, 1, 0], # 4 -> 3
]
)
3.3.3. Branch-and-cut for the TSP¶
Exercice: Complete the tsp_branch_and_cut method below.
from typing import Callable
def tsp_branch_and_cut(
distances: Sequence[Sequence[float]],
ilp_solver: Callable[[minilp.problem], minilp.result] = branch_and_bound,
):
"""
Solves the given TSP instance using a branch-and-cut with the given solver.
Args:
distances: Transport costs for the TSP.
ilp_solver: Function that can be called as ilp_solver(ilp) to solve integer linear program.
Returns:
A minilp.result solution for the given TSP instance.
"""
# Create the relaxation:
N = len(distances)
tsp, x = tsp_relax(distances, "TSP Branch & Cut")
while True:
# Solve the problem:
res = ilp_solver(tsp)
... # TODO
return res
Exercice: Test your tsp_branch_and_cut implementation.
import tsp.data as data
res = tsp_branch_and_cut(data.grid5)