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 thescipy
module, wildly available.minilp.solvers.docplex
— A solver based on thedocplex
module, 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_bound
procedure. - 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)