17 KiB
TP2 - Modeling using docplex¶
1. The docplex python package¶
DOcplex is a python package developped by IBM — It provides easy-to-use API for IBM solvers Cplex and Cpoptimizer.
DOcplex documentation for mathematical programming can be found here: http://ibmdecisionoptimization.github.io/docplex-doc/#mathematical-programming-modeling-for-python-using-docplex-mp-docplex-mp
from docplex.mp.model import Model
import tsp.data as data
distances = data.grid42
N = len(distances)
tsp = Model("TSP")
tsp.log_output = True
... # TODO
solution = tsp.solve()
print("z* =", solution.objective_value)
The largest set of distances contains 42 nodes, and should be easily solved by docplex.
2.2. Generating random TSP instances¶
Question: What method could you implement to generate a realistic set of distances for $n$ customers?
Exercice: Implement the method proposed above and test it.
import numpy as np
import scipy.spatial.distance
def generate_distances(n: int):
... # TODO
from docplex.mp.model import Model
distances = generate_distances(50)
print(distances)
N = len(distances)
tsp = Model("TSP")
tsp.log_output = True
... # TODO: Copy your model from the first question here.
solution = tsp.solve()
print("z* =", solution.objective_value)
3. Solving Warehouse Allocation using Benders decomposition with docplex¶
3.1. The warehouse problem¶
A company needs to supply a set of $n$ clients and needs to open new warehouses (from a set of $m$ possible warehouses). Opening a warehouse $j$ costs $f_j$ and supplying a client $i$ from a warehouse $j$ costs $c_{ij}$ per supply unit. Which warehouses should be opened in order to satisfy all clients while minimizing the total cost?
3.2. Solving the warehouse problem with a single ILP¶
- $y_{j} = 1$ if and only if warehouse $j$ is opened.
- $x_{ij}$ is the fraction supplied from warehouse $j$ to customer $i$.
$ \begin{align} \text{min.} \quad & \sum_{i=1}^{n} \sum_{j=1}^{m} c_{ij} x_{ij} + \sum_{j=1}^{m} f_{j} y_{j} & \\ \text{s.t.} \quad & \sum_{j=1}^{m} x_{ij} = 1, & \forall i \in\{1,\ldots,n\}\\ & x_{ij} \leq y_{j}, & \forall i\in\{1,\ldots,n\},\forall j\in\{1,\ldots,m\}\\ & y_{j} \in \left\{0,~1\right\}, & \forall j \in \left\{1,~\ldots,~m\right\}\\ & x_{ij} \geq 0, & \forall i \in \left\{1,~\ldots,~n\right\}, \forall j \in \left\{1,~\ldots,~m\right\} \end{align} $
Exercice: Implement the ILP model for the warehouse allocation problem and test it on the given instance.
from docplex.mp.model import Model
# We will start with a small instances with 3 warehouses and 3 clients:
N = 3
M = 3
# Opening and distribution costs:
f = [20, 20, 20]
c = [[15, 1, 2], [1, 16, 3], [4, 1, 17]]
wa = Model("Warehouse Allocation")
wa.log_output = True
... # TODO: Model for the warehouse allocation.
solution = wa.solve()
print("z* =", solution.objective_value)
3.3. Benders' decomposition for the Warehouse Allocation problem¶
We are going to use Benders' decomposition to solve the Warehouse Allocation problem.
Dual subproblem¶
$ \begin{align*} \text{max.} \quad & \sum_{i=1}^{n} v_{i} - \sum_{i=1}^{n}\sum_{j=1}^{m} \bar{y}_{j} u_{ij} & \\ \text{s.t.} \quad & v_{i} - u_{ij} \leq c_{ij}, & \forall i\in\{1,\ldots,n\},\forall j\in\{1,\ldots,m\}\\ & v_{i} \in\mathbb{R},\ u_{ij} \geq 0 & \forall i \in\{1,\ldots,n\}, \forall j\in\{1,\ldots,m\} \end{align*} $
Master problem¶
$ \begin{align*} \text{min.} \quad & \sum_{j=1}^{m} f_j y_j + z & \\ \text{s.t.} \quad & z \geq \sum_{i=1}^{n}v_i^p - \sum_{i=1}^{n} \sum_{j=1}^{m} u_{ij}^p y_j, & \forall p\in l_1\\ & 0 \geq \sum_{i=1}^{n}v_i^r - \sum_{i=1}^{n} \sum_{j=1}^{n} u_{ij}^r y_j, & \forall r\in l_2\\ & y_{j} \in\{0,1\}, & \forall j\in\{1,\ldots,m\} \end{align*} $
Exercice: Implement the method create_master_problem that creates the initial master problem (without feasibility or optimality constraints) for the warehouse allocation problem.
You can use print(m.export_as_lp_string()) to display a textual representation of a docplex model m.
from docplex.mp.model import Model
from docplex.mp.linear import Var
from docplex.mp.constr import AbstractConstraint
from typing import List, Sequence, Tuple
def create_master_problem(
N: int, M: int, f: Sequence[float], c: Sequence[Sequence[float]]
) -> Tuple[Model, Var, Sequence[Var]]:
"""
Creates the initial Benders master problem for the Warehouse Allocation problem.
Args:
N: Number of clients.
M: Number of warehouses.
f: Array-like containing costs of opening warehouses.
c: 2D-array like containing transport costs from client to warehouses.
Returns:
A 3-tuple containing the docplex problem, the z variable and the y variables.
"""
wa = Model("Warehouse Allocation - Benders master problem")
... # TODO
return wa, z, y
# Check your method:
wa, z, y = create_master_problem(N, M, f, c)
print(wa.export_as_lp_string())
Exercice: Implement the method add_optimality_constraints that add optimality constraints to the Benders master problem.
def add_optimality_constraint(
N: int,
M: int,
wa: Model,
z: Var,
y: Sequence[Var],
v: Sequence[float],
u: Sequence[Sequence[float]],
) -> List[AbstractConstraint]:
"""
Adds an optimality constraints to the given Warehouse Allocation model
using the given optimal values from the Benders dual subproblem.
Args:
N: Number of clients.
M: Number of warehouses.
wa: The Benders master problem (docplex.mp.model.Model).
z: The z variable of the master problem.
y: The y variables of the master problem.
v: The optimal values for the v variables of the Benders dual subproblem.
u: The optimal values for the u variables of the Benders dual subproblem.
Return: The optimality constraint added.
"""
return ... # TODO
Exercice: Implement the method add_feasibility_constraints that add feasibility constraints to the Benders master problem.
def add_feasibility_constraints(
N: int,
M: int,
wa: Model,
z: Var,
y: Sequence[Var],
v: Sequence[float],
u: Sequence[Sequence[float]],
) -> List[AbstractConstraint]:
"""
Adds an optimality constraints to the given Warehouse Allocation model
using the given optimal values from the Benders dual subproblem.
Args:
- N: Number of clients.
- M: Number of warehouses.
- wa: The Benders master problem (docplex.mp.model.Model).
- z: The z variable of the master problem.
- y: The y variables of the master problem.
- v: The extreme rays for the v variables of the Benders dual subproblem.
- u: The extreme rays for the u variables of the Benders dual subproblem.
Returns:
The feasibility constraint added.
"""
# TODO:
return
Exercice: Implement the method create_dual_subproblem that, given a solution y of the master problem, create the corresponding Benders dual subproblem.
$ \begin{align*} \text{max.} \quad & \sum_{i=1}^{n} v_{i} - \sum_{i=1}^{n}\sum_{j=1}^{m} \bar{y}_{j} u_{ij} & \\ \text{s.t.} \quad & v_{i} - u_{ij} \leq c_{ij}, & \forall i\in\{1,\ldots,n\},\forall j\in\{1,\ldots,m\}\\ & v_{i} \in\mathbb{R},\ u_{ij} \geq 0 & \forall i \in\{1,\ldots,n\}, \forall j\in\{1,\ldots,m\} \end{align*} $
from docplex.mp.model import Model
def create_dual_subproblem(
N: int, M: int, f: Sequence[float], c: Sequence[Sequence[float]], y: Sequence[int]
) -> Tuple[Model, Sequence[Var], Sequence[Sequence[Var]]]:
"""
Creates a Benders dual subproblem for the Warehouse Allocation problem corresponding
to the given master solution.
Args:
N: Number of clients.
M: Number of warehouses.
f: Array-like containing costs of opening warehouses.
c: 2D-array like containing transport costs from client to warehouses.
y: Values of the y variables from the Benders master problem.
Returns:
A 3-tuple containing the docplex problem, the v variable and the u variables.
"""
dsp = Model("Warehouse Allocation - Benders dual subproblem")
# We disable pre-solve to be able to retrieve a meaningful status in the main
# algorithm:
dsp.parameters.preprocessing.presolve.set(0)
... # TODO
return dsp, v, u
# Check your method (assuming y = [1 1 1 ... 1]):
dsp, v, u = create_dual_subproblem(N, M, f, c, [1] * M)
print(dsp.export_as_lp_string())
Exercice: Using the methods you implemented, write the Benders decomposition algorithm for the warehouse allocation problem.
The get_extreme_rays function can be used to retrieve the extreme rays associated with an unbounded solution of the dual subproblem.
You can use model.get_solve_status() to obtain the status of the resolution and compare it to members of JobSolveStatus:
if model.get_solve_status() == JobSolveStatus.OPTIMAL_SOLUTION:
pass
from docplex.mp.model import Model
from docplex.util.status import JobSolveStatus
def get_extreme_rays(
N: int, M: int, model: Model, v: Sequence[Var], u: Sequence[Sequence[Var]]
) -> Tuple[Sequence[float], Sequence[Sequence[float]]]:
"""
Retrieves the extreme rays associated to the dual subproblem.
Args:
N: Number of clients.
M: Number of warehouses.
model: The Benders dual subproblem model (docplex.mp.model.Model).
v: 1D array containing the v variables of the subproblem.
u: Either a 2D array of a tuple-index dictionary containing the u variables
of the subproblem.
Returns:
A 2-tuple containing the list of extreme rays correspondig to v,
and the 2D-list of extreme rays corresponding to u.
"""
ray = model.get_engine().get_cplex().solution.advanced.get_ray()
if isinstance(u, dict):
def get_uij(i, j):
return u[i, j]
else:
def get_uij(i, j):
return u[i][j]
return (
[ray[v[i].index] for i in range(N)],
[[ray[get_uij(i, j).index] for j in range(M)] for i in range(N)],
)
# We will start with a small instances with 3 warehouses and 3 clients:
N = 3
M = 3
# Opening and distribution costs:
f = [20, 20, 20]
c = [[15, 1, 2], [1, 16, 3], [4, 1, 17]]
# We stop iterating if the new solution is less than epsilon
# better than the previous one:
epsilon = 1e-6
wa, z, y = create_master_problem(N, M, f, c)
n = 0
while True:
# Print iteration:
n = n + 1
print("Iteration {}".format(n))
... # TODO
print("Done.")
3.4. Generating instances for the Warehouse Allocation problem¶
Exercice: Using the TSP instances contained in tsp.data or the generate_distances method, create instances for the warehouse allocation problem with randomized opening costs.