AMPL Capacitated p-Median Problem with GCG#
Description: Dantzig-Wolfe decomposition for Capacitated p-Median Problem with GCG
Tags: GCG, cpmp, amplpy, dantzig-wolfe decomposition, branch-price-and-cut, highlights
Notebook author: Jurgen Lentz <jurgenlentz26@gmail.com>
# Install dependencies
%pip install -q amplpy
# Google Colab & Kaggle integration
from amplpy import AMPL, ampl_notebook
ampl = ampl_notebook(
modules=["gcg"], # modules to install
license_uuid="default", # license to use
) # instantiate AMPL object and register magics
Capacitated \(p\)-median problem#
The goal of the capacitated \(p\)-median problem (CPMP) is to fix \(p \in \mathbb{N}\) locations as medians from a given number of \(n \in \mathbb{N}\) locations. Furthermore, we are given demands \(q_{i} \geq 0\) and capacities \(Q_{i} \geq 0\) for every location \(i \in \{1,\ldots,n\}\) as well as distances \(d_{ij} \geq 0\) for each pair of locations \(i,j \in \{1,\ldots,n\}\). While fixing the medians, we also assign every location \(i \in \{1,\ldots,n\}\) to exactly one median such that the demands of all locations assigned to the median can be satisfied by the capacity of the median. An assignment of locations to a median is called a cluster. The overall objective is to minimize the total distance between the medians and their assigned locations.
CPMP can be modeled as follows (see Ceselli et al.):
\[\begin{split}
\begin{aligned}
\text{minimize} \quad &\sum_{i = 1}^{n} \sum_{j = 1}^{n} d_{i j} x_{i j} \\
\text{subject to} \quad &\sum_{j = 1}^{n} x_{i j} = 1 \quad \forall i \in \{1,...,n\} \\
&\sum_{i = 1}^{n} q_{i} x_{i j} \leq Q_{j} y_{j} \quad \forall j \in \{1,...,n\} \\
&\sum_{j = 1}^{n} y_{j} = p \\
&x_{i j} \in \{0,1\} \quad \forall i,j \in \{1,...,n\} \\
&y_{j} \in \{0,1\} \quad \forall j \in \{1,...,n\}
\end{aligned}
\end{split}\]
We will solve an CPMP instance using Dantzig-Wolfe decomposition with the GCG solver in AMPL. Therefore, we first model CPMP in AMPL (shown below).
%%ampl_eval
param n;
param p;
set I = 1..n ordered;
set J = 1..n ordered;
param d {I,J};
param w {I};
param C {J};
var x {I,J} binary;
var y {J} binary;
minimize Cost: sum {i in I} sum {j in J} d[i,j] * x[i,j];
subject to Allocate {i in I}:
sum {j in J} x[i,j] = 1;
subject to Capacity {j in J}:
sum {i in I} w[i] * x[i,j] <= C[j] * y[j];
subject to NFacilities:
sum{j in J} y[j] <= p;
We generate a small instance with 5 locations and 2 locations as medians.
ampl.param["n"] = 5
ampl.param["p"] = 2
d = [
[0, 6, 54, 52, 19],
[6, 0, 28, 75, 61],
[54, 28, 0, 91, 40],
[52, 75, 91, 0, 28],
[19, 61, 40, 28, 0],
]
ampl.param["d"] = {(i, j): d[i - 1][j - 1] for i in range(1, 6) for j in range(1, 6)}
ampl.param["w"] = [14, 13, 9, 15, 6]
ampl.param["C"] = [39, 39, 39, 39, 39]
Automatic Mode in GCG with AMPL#
We use AMPL to call the solver GCG to solve our CPMP instance automatically without providing any information about the Dantzig-Wolfe decomposition. Here, GCG detects different decompositions and chooses heuristically the best decomposition. Afterwards, the solver uses a branch-price-and-cut algorithm to solve it to optimality.
ampl.option["solver"] = "gcg"
ampl.option["gcg_options"] = "outlev=1"
ampl.solve()
GCG 4.0.0: tech:outlev = 1
presolving:
(round 1, exhaustive) 0 del vars, 0 del conss, 0 add conss, 0 chg bounds, 0 chg sides, 0 chg coeffs, 10 upgd conss, 0 impls, 5 clqs
(round 2, exhaustive) 0 del vars, 0 del conss, 0 add conss, 0 chg bounds, 0 chg sides, 0 chg coeffs, 11 upgd conss, 0 impls, 30 clqs
(0.0s) probing cycle finished: starting next cycle
presolving (3 rounds: 3 fast, 3 medium, 3 exhaustive):
0 deleted vars, 0 deleted constraints, 0 added constraints, 0 tightened bounds, 0 added holes, 0 changed sides, 0 changed coefficients
0 implications, 30 cliques
presolved problem has 30 variables (30 bin, 0 int, 0 impl, 0 cont) and 11 constraints
6 constraints of type <knapsack>
5 constraints of type <setppc>
transformed objective value is always integral (scale: 1)
Presolving Time: 0.00
Consclassifier "nonzeros" yields a classification with 2 different constraint classes
Consclassifier "constypes" yields a classification with 2 different constraint classes
Consclassifier "constypes according to miplib" yields a classification with 3 different constraint classes
Consclassifier "gamsdomain" yields a classification with 1 different constraint classes
Consclassifier "gamssymbols" yields a classification with 1 different constraint classes
Conspartition "gamssymbols" is not considered since it offers the same structure as "gamsdomain" conspartition
Varclassifier "gamsdomain" yields a classification with 1 different variable classes
Varclassifier "gamssymbols" yields a classification with 1 different variable classes
Varpartition "gamssymbols" is not considered since it offers the same structure as "gamsdomain"
Varclassifier "vartypes" yields a classification with 1 different variable classes
Varpartition "vartypes" is not considered since it offers the same structure as "gamsdomain"
Varclassifier "varobjvals" yields a classification with 10 different variable classes
Varclassifier "varobjvalsigns" yields a classification with 2 different variable classes
Added reduced version of varpartition varobjvals with 9 different variable classes
in dec_consclass: there are 4 different constraint classes
the current constraint classifier "nonzeros" consists of 2 different classes
the current constraint classifier "constypes" consists of 2 different classes
the current constraint classifier "constypes according to miplib" consists of 3 different classes
the current constraint classifier "gamsdomain" consists of 1 different classes
dec_consclass found 14 new partialdecs
dec_densemasterconss found 1 new partialdec
dec_neighborhoodmaster found 1 new partialdec
the current varclass distribution includes 10 classes but only 9 are allowed for propagatePartialdec() of var class detector
POSTPROCESSING of decompositions. Added 0 new decomps.
Found 521 finished decompositions.
Measured running time per detector:
Detector consclass worked on 7 finished decompositions and took a total time of 0.000
Detector neighborhoodmaster worked on 1 finished decompositions and took a total time of 0.000
Detector connectedbase worked on 520 finished decompositions and took a total time of 0.003
Detector varclass worked on 512 finished decompositions and took a total time of 0.002
Detection Time: 0.02
A Dantzig-Wolfe reformulation is applied to solve the original problem.
Chosen structure has 6 blocks and 5 linking constraints.
This decomposition has a maxwhite score of 0.454545.
Matrix has 6 blocks, using 6 pricing problems.
time | node | left |SLP iter|MLP iter|LP it/n| mem |mdpt |ovars|mvars|ocons|mcons|mcuts| dualbound | primalbound | deg | gap
p 0.0s| 1 | 0 | 0 | 0 | - |3798k| 0 | 30 | 0 | 12 | 0 | 0 | 0.000000e+00 | 1.480000e+02 | -- | Inf
p 0.0s| 1 | 0 | 0 | 0 | - |3799k| 0 | 30 | 0 | 12 | 0 | 0 | 0.000000e+00 | 8.800000e+01 | -- | Inf
p 0.0s| 1 | 0 | 0 | 0 | - |3797k| 0 | 30 | 0 | 12 | 0 | 0 | 0.000000e+00 | 6.200000e+01 | -- | Inf
0.0s| 1 | 0 | 0 | 0 | - |3793k| 0 | 30 | 0 | 12 | 0 | 0 | 0.000000e+00 | 6.200000e+01 | -- | Inf
0.0s| 1 | 0 | 0 | 0 | - |3857k| 0 | 30 | 6 | 12 | 12 | 0 | 0.000000e+00 | 6.200000e+01 | 0.00%| Inf
0.0s| 1 | 0 | 0 | 0 | - |3886k| 0 | 30 | 18 | 12 | 12 | 0 | 0.000000e+00 | 6.200000e+01 | 0.00%| Inf
Starting reduced cost pricing...
0.0s| 1 | 0 | 14 | 14 | - |3899k| 0 | 30 | 21 | 12 | 12 | 0 | 0.000000e+00 | 6.200000e+01 | 9.09%| Inf
0.0s| 1 | 0 | 26 | 26 | - |3923k| 0 | 30 | 22 | 12 | 12 | 3 | 0.000000e+00 | 6.200000e+01 | 6.06%| Inf
0.0s| 1 | 0 | 28 | 28 | - |3929k| 0 | 30 | 24 | 12 | 12 | 3 | 3.542308e+01 | 6.200000e+01 | 5.41%| 75.03%
0.0s| 1 | 0 | 31 | 31 | - |3943k| 0 | 30 | 24 | 12 | 12 | 4 | 3.979412e+01 | 6.200000e+01 | 5.59%| 55.80%
0.0s| 1 | 0 | 32 | 32 | - |3954k| 0 | 30 | 24 | 12 | 12 | 5 | 4.047059e+01 | 6.200000e+01 | 4.89%| 53.20%
0.0s| 1 | 0 | 33 | 33 | - |3961k| 0 | 30 | 24 | 12 | 12 | 6 | 4.047059e+01 | 6.200000e+01 | 4.89%| 53.20%
0.0s| 1 | 0 | 33 | 33 | - |3994k| 0 | 30 | 24 | 12 | 12 | 6 | 6.200000e+01 | 6.200000e+01 | -- | 0.00%
SCIP Status : problem is solved [optimal solution found]
Solving Time (sec) : 0.05
Solving Nodes : 1
Primal Bound : +6.20000000000000e+01 (5 solutions)
Dual Bound : +6.20000000000000e+01
Gap : 0.00 %
WARNING: Dual information only available for pure LPs (only continuous variables).
GCG 4.0.0: optimal solution; objective 62
10 simplex iterations
1 branching nodes
Using a custom decomposition in GCG with AMPL#
AMPL allows the users to create their own decomposition and forwards it to GCG using suffixes. Here, we assign the allocation and pmedian constraints to the master/linking constraints and each capacity constraint to a different pricing problem.
%%ampl_eval
suffix master IN, binary;
suffix block IN, integer;
let {i in I}
Allocate[i].master := 1;
let NFacilities.master := 1;
let {j in J}
Capacity[j].block := j;
ampl.solve()
GCG 4.0.0: tech:outlev = 1
added complete decomp for original problem with 5 blocks and 6 masterconss, 0 linkingvars, 0 mastervars, and max white score of 0.363636
presolving:
(round 1, exhaustive) 0 del vars, 0 del conss, 0 add conss, 0 chg bounds, 0 chg sides, 0 chg coeffs, 10 upgd conss, 0 impls, 5 clqs
(round 2, exhaustive) 0 del vars, 0 del conss, 0 add conss, 0 chg bounds, 0 chg sides, 0 chg coeffs, 11 upgd conss, 0 impls, 30 clqs
(0.0s) probing cycle finished: starting next cycle
presolving (3 rounds: 3 fast, 3 medium, 3 exhaustive):
0 deleted vars, 0 deleted constraints, 0 added constraints, 0 tightened bounds, 0 added holes, 0 changed sides, 0 changed coefficients
0 implications, 30 cliques
presolved problem has 30 variables (30 bin, 0 int, 0 impl, 0 cont) and 11 constraints
6 constraints of type <knapsack>
5 constraints of type <setppc>
transformed objective value is always integral (scale: 1)
Presolving Time: 0.00
calculated translation; number of missing constraints: 0; number of other partialdecs: 1
Preexisting decomposition found. Solution process started...
A Dantzig-Wolfe reformulation is applied to solve the original problem.
Chosen structure has 5 blocks and 6 linking constraints.
This decomposition has a maxwhite score of 0.363636.
Matrix has 5 blocks, using 5 pricing problems.
time | node | left |SLP iter|MLP iter|LP it/n| mem |mdpt |ovars|mvars|ocons|mcons|mcuts| dualbound | primalbound | deg | gap
p 0.0s| 1 | 0 | 0 | 0 | - |2766k| 0 | 30 | 0 | 12 | 0 | 0 | 0.000000e+00 | 1.480000e+02 | -- | Inf
p 0.0s| 1 | 0 | 0 | 0 | - |2767k| 0 | 30 | 0 | 12 | 0 | 0 | 0.000000e+00 | 8.800000e+01 | -- | Inf
p 0.0s| 1 | 0 | 0 | 0 | - |2765k| 0 | 30 | 0 | 12 | 0 | 0 | 0.000000e+00 | 6.200000e+01 | -- | Inf
0.0s| 1 | 0 | 0 | 0 | - |2761k| 0 | 30 | 0 | 12 | 0 | 0 | 0.000000e+00 | 6.200000e+01 | -- | Inf
0.0s| 1 | 0 | 0 | 0 | - |2825k| 0 | 30 | 5 | 12 | 12 | 0 | 0.000000e+00 | 6.200000e+01 | 0.00%| Inf
0.0s| 1 | 0 | 0 | 0 | - |2850k| 0 | 30 | 15 | 12 | 12 | 0 | 0.000000e+00 | 6.200000e+01 | 0.00%| Inf
Starting reduced cost pricing...
0.0s| 1 | 0 | 16 | 16 | - |2888k| 0 | 30 | 26 | 12 | 12 | 0 | 6.200000e+01 | 6.200000e+01 | 54.55%| 0.00%
0.0s| 1 | 0 | 16 | 16 | - |2862k| 0 | 30 | 15 | 12 | 12 | 0 | 6.200000e+01 | 6.200000e+01 | -- | 0.00%
SCIP Status : problem is solved [optimal solution found]
Solving Time (sec) : 0.01
Solving Nodes : 1
Primal Bound : +6.20000000000000e+01 (3 solutions)
Dual Bound : +6.20000000000000e+01
Gap : 0.00 %
WARNING: Dual information only available for pure LPs (only continuous variables).
GCG 4.0.0: optimal solution; objective 62
10 simplex iterations
1 branching nodes
Analogously, the user can fix constraints as pricing block or master constraint and variables as pricing block, master or linking variables. It is not needed to provide a complete decomposition. If the user provides a partial decomposition, GCG completes the decomposition by fixing only the left constraints and variables using its detection loop.
NOTE: The index of pricing block constraints/variables starts at 1.