Scippy

GCG

Branch-and-Price & Column Generation for Everyone

The Magic of Genericity

Table of Contents

The Magic of Genericity

In this use case, we make GCG do the task it has been made for: solving structured mixed-integer programs. We will not go into any details concerning decomposition, reformulation or branch-and-price-and-cut, but only explain where one can find these steps in the output log and - just very superficially - what they mean.

Introduction

This is Andrew. He lives near Edinburgh and has just founded his Start-up "Edin-OR" which specializes on optimizing planning and decision making processes within companies. Luckily, he has a friend over at Royal Mail, the local postal services company, who might help him with getting his first real job, scheduling the vehicles that deliver the mail. Since the management is doubtful of Andrew's capacilities, they ask him to begin with artificial data that is, however, very close to theirs. As a start, they suggest instances first used by Christofides et al. [2].
As always, Andrew starts with the exact problem definition.

Problem Definition


Sources: motoringresearch.com and electrek.co
The Royal Mail has a fleet consisting of about 50.000 vehicles and a subset of them operates in the Edinburgh region. Soon getting new, more efficient and environmentally-friendly electric vehicles, the local Royal Mail also wants to optimize the routes that they (and older vehicles) drive to deliver mail.

The Capacitated Vehicle Routing Problem

The goal is to deliver mail to all their customers using a given number of vehicles with limited capacity each. These will each drive their tour, starting from the depot and returning after their shift (after having delivered the mail to all recipients). Doing this, they don't have any time constraints, but instead just the capacity of their own vehicle. This is called the Capacitated Vehicle Routing Problem and has already been known for quite some time. A standard model is defined as following:

\begin{align} \min\quad & \sum_{(i,j) \in E} c_{ij}x_{ij}^{k}\\ s.t. \quad & \sum_{(i,j)\in \delta^{+} (i)} x_{ij}^{k}=\sum_{(j,i)\in \delta^{-} (i)} x_{ji}^{k}= \begin{cases} \begin{aligned} 1 &\qquad i=d,\\ z_{i}^{k} &\qquad i\neq d \end{aligned} \end{cases} \qquad && k \in \{1,\dots,K\}\\ & \sum_{k\in\{1,\dots,K\}} z_{i}^{k} = 1 && i\in V\\ & \sum_{(i,j)\in\delta^+(S)}x_{ij}\geq 1 \qquad && \emptyset \neq S \subsetneq V \\ & x_{ij}^{k} \in \{0,1\} && k\in\{1,\dots,K\},\ (i,j)\in E\\ & z_{i}^{k} \in \{0,1\} && k\in\{1,\dots,K\},\ i\in V\\ \end{align}

Andrew tries to choose an instance of medium difficulty to start with, even though they might become much bigger once he gets real data. The instance is called E021-06m. It features 21 customers to visit and 6 vehicles, each with a capacity of 4000. Note that the model by Christofides et al. is not exactly as the one above, since they refined it a bit further (see [2] for more information).
Andrew could also specify the problem himself using e.g. ZIMPL (see second use case).

Searching for a Solver

With the problem at hand, Andrew starts searching for solvers of mixed integer programs. He quickly comes across the most common ones like Gurobi and SCIP. The latter has the advantage, that it is source-open and can be used without having to purchase a license, which is a big advantage because Andrew just wants do to a bit of testing and furthermore does not have that much capital yet (note that for commercial use of GCG and SCIP, the SCIP team has to be contacted first, since they require a different license). Thus, he tries to solve his instance with SCIP. After more than two hours he aborts the solving:

...
time | node | left |LP iter|LP it/n|mem/heur|mdpt |vars |cons |rows |cuts |sepa|confs|strbr|dualbound |primalbound |gap |compl.
7570s| 1014k|140418 | 39911k| 39.2 | 3302M | 73 |2682 |9236 |2761 | 757k| 3 | 863k| 114k|3.634914e+02 |4.308847e+02 |18.54%|22.30%
7571s| 1014k|140427 | 39914k| 39.2 | 3302M | 73 |2682 |9238 |2761 | 757k| 0 | 863k| 114k|3.634914e+02 |4.308847e+02 |18.54%|22.30%
7572s| 1014k|140440 | 39918k| 39.2 | 3302M | 73 |2682 |9239 |2740 | 757k| 1 | 863k| 114k|3.634927e+02 |4.308847e+02 |18.54%|22.30%
^Cpressed CTRL-C 1 times (5 times for forcing termination)
SCIP Status : solving was interrupted [user interrupt]
Solving Time (sec) : 7572.69
Solving Nodes : 1014739 (total of 1017553 nodes in 2 runs)
Primal Bound : +4.30884678000000e+02 (31 solutions)
Dual Bound : +3.63492688685070e+02
Gap : 18.54 %
SCIP>

He thinks that he could do better, since Andrew knows that the problem at hand is a structured integer program. Thus, he starts searching for state-of-the-art methods for solving those. He encounters many proprietary approaches, where scientists had implemented their own customized solvers for their own vehicle routing and many other problems, but since Andrew has not yet gotten the hang of all this "decomposing", "pricing" or "cutting" and he also wants to reuse code as efficiently and often as possible, he starts looking for generic solvers that would do all of that automatically for him and finally finds "Generic Column Generation", which seems to be fitting. He is happy about the interface of GCG being just like those of all other commonly used solvers, but since the search for existing approaches had taken him so long, there is not much time left to solve the vehicle routing problem the management gave him, so he just goes ahead and installs GCG.

First Time: Using GCG

After the installation, Andrew fires up GCG with a simple

./bin/gcg

The output he is getting shows him the software versions and some copyright information. If he were to add more external code, which could improve solving speed, they would appear under "External Codes". As LP solver, he again can rely on the source-open solver "SoPlex", which, in contrast to "CPLEX" did not cost him anything. It is part of the whole SCIP Optimization Suite. Finally, he can see that the default settings are being used, which is fine, since he will not do any fine-tuning.

GCG version 3.1.0 [GitHash: 90d8bdf3]
Copyright (C) 2010-2020 Operations Research, RWTH Aachen University
Konrad-Zuse-Zentrum fuer Informationstechnik Berlin (ZIB)
SCIP version 7.0.0 [precision: 8 byte] [memory: block] [mode: optimized] [LP solver: SoPlex 5.0.0] [GitHash: 0bc4dc9c6]
Copyright (C) 2002-2020 Konrad-Zuse-Zentrum fuer Informationstechnik Berlin (ZIB)
External codes:
Readline 8.0 GNU library for command line editing (gnu.org/s/readline)
SoPlex 5.0.0 Linear Programming Solver developed at Zuse Institute Berlin (soplex.zib.de) [GitHash: 6535a3c8]
CppAD 20180000.0 Algorithmic Differentiation of C++ algorithms developed by B. Bell (www.coin-or.org/CppAD)
ZLIB 1.2.11 General purpose compression library by J. Gailly and M. Adler (zlib.net)
user parameter file <gcg.set> not found - using default parameters

To continue, Andrew now wants to read in his problem. This can be done with a simple read and then the problem file. GCG will signalize that the problem was read successfully and already give some statistics: E021-06m only contains a few continuous variables and quite some integer variables and 2600 inequalities that are constraining the problem.

GCG> read E021-06m.cvrp.lp.gz
read problem <E021-06m.cvrp.lp.gz>
============
original problem has 2766 variables (0 bin, 2646 int, 0 impl, 120 cont) and 2600 constraints

Finally, after looking at the Getting started Guide in GCG's documentation, he proceeds with optimizing his problem.

GCG> optimize
...
Presolving Time: 0.75
...
Detection Time: 0.24
...
A Dantzig-Wolfe reformulation is applied to solve the original problem.
...
time | node | left |SLP iter|MLP iter|LP it/n| mem |mdpt |ovars|mvars|ocons|mcons|mcuts| dualbound | primalbound | deg | gap
p 0.9s| 1 | 0 | 0 | 0 | - | 23M| 0 |2682 | 0 |2588 | 0 | 0 | 0.000000e+00 | 6.217321e+02 | -- | Inf
1.0s| 1 | 0 | 0 | 0 | - | 23M| 0 |2682 | 0 |2588 | 0 | 0 | 3.242201e+02 | 6.217321e+02 | -- | 91.76%
...
SCIP Status : problem is solved [optimal solution found]
Solving Time (sec) : 39.13
Solving Nodes : 1
Primal Bound : +4.30884678000000e+02 (2 solutions)
Dual Bound : +4.30884678000000e+02
Gap : 0.00 %

Having already worked with mixed integer programming solvers before, Andrew sees that there are many similarities to them. The optimal solution can be read at the bottom ("Primal Bound"), along with the time needed to solve it. If Andrew aborted before GCG was able to find an optimal solution, he would still be able to see the two bounds:

SCIP Status : solving was interrupted
Solving Time (sec) : 8.23
Solving Nodes : 1
Primal Bound : +4.30884678000000e+02 (2 solutions)
Dual Bound : +4.18548834278730e+02
Gap : 2.95 %

And would know that the objective function value of the optimal solution lies somewhere between 419 and 430. This result, which was obtained in just a few seconds, is already stronger than the two bounds that a non BP&C solver like SCIP calculated in over 2 hours (see above).

Understanding the Output

For many, this is not very relevant (they can just skip to the section Getting the Solution), but since Andrew wants to understand what is happening (at least someday), he takes a brief look at what was printed there. Before GCG starts with the "real" solving process (the Branch and Price and Cutting), there are a few steps that are executed and that distinguish GCG from other solvers.

Presolving

More details and the corresponding theory and implementation can be found on the page The GCG Presolving.

During the presolving, trivially fulfilled or unsatisfied constraints are indentified as well as redundant variables and more. It can happen that a problem is already solved after presolving. The full log of a presolving process looks as following:

presolving:
(round 1, fast) 0 del vars, 120 del conss, 0 add conss, 120 chg bounds, 0 chg sides, 6 chg coeffs, 0 upgd conss, 0 impls, 314 clqs
(round 2, exhaustive) 0 del vars, 120 del conss, 0 add conss, 120 chg bounds, 0 chg sides, 6 chg coeffs, 284 upgd conss, 0 impls, 314 clqs
(round 3, medium) 0 del vars, 126 del conss, 114 add conss, 120 chg bounds, 0 chg sides, 6 chg coeffs, 284 upgd conss, 0 impls, 440 clqs
(round 4, exhaustive) 26 del vars, 126 del conss, 114 add conss, 120 chg bounds, 0 chg sides, 12 chg coeffs, 284 upgd conss, 0 impls, 6230 clqs
(round 5, exhaustive) 62 del vars, 126 del conss, 114 add conss, 120 chg bounds, 0 chg sides, 12 chg coeffs, 284 upgd conss, 0 impls, 9622 clqs
(0.5s) probing: 1000/2646 (37.8%) - 84 fixings, 0 aggregations, 356518 implications, 0 bound changes
(0.6s) probing: 1190/2646 (45.0%) - 84 fixings, 0 aggregations, 451993 implications, 0 bound changes
(0.6s) probing aborted: 1000/1000 successive useless probings
(round 6, exhaustive) 84 del vars, 126 del conss, 114 add conss, 120 chg bounds, 0 chg sides, 12 chg coeffs, 284 upgd conss, 1490 impls, 9793 clqs
(0.7s) probing: 1290/2646 (48.8%) - 84 fixings, 0 aggregations, 503003 implications, 0 bound changes
(0.7s) probing aborted: 1000/1000 successive useless probings
presolving (7 rounds: 7 fast, 6 medium, 5 exhaustive):
84 deleted vars, 126 deleted constraints, 114 added constraints, 120 tightened bounds, 0 added holes, 0 changed sides, 12 changed coefficients
1690 implications, 9995 cliques
presolved problem has 2682 variables (2562 bin, 0 int, 0 impl, 120 cont) and 2588 constraints
6 constraints of type <knapsack>
386 constraints of type <setppc>
2196 constraints of type <linear>
Presolving Time: 0.70

Andrew easily identifies that the presolving happens in rounds (1-7), where each round has different settings and different outcomes. Most importantly, in the line

presolving (7 rounds: 7 fast, 6 medium, 5 exhaustive):
84 deleted vars, 126 deleted constraints, 114 added constraints, 120 tightened bounds, 0 added holes, 0 changed sides, 12 changed coefficients

he can see how variables and constraints were deleted and thus do not need to be respected anymore and the "bounds were tightened", meaning that these minor changes to his instance lead to a potentially faster solving of his problem.

Detection

More details and the corresponding theory and implementation can be found on the page The GCG Detection.

During the detection, one of GCG's main features, the solver will find different structures in your instance and will then "decompose" the problem such that it can be reformulated (see next section).

Consclassifier "nonzeros" yields a classification with 6 different constraint classes
Consclassifier "constypes" yields a classification with 4 different constraint classes
Consclassifier "constypes according to miplib" yields a classification with 5 different constraint classes
Varclassifier "vartypes" yields a classification with 2 different variable classes
Varclassifier "varobjvals" yields a classification with 179 different variable classes
Varclassifier "varobjvalsigns" yields a classification with 2 different variable classes
the current varclass distribution includes 179 classes but only 18 are allowed for GCGconshdlrDecompCalcCandidatesNBlocks()
in dec_consclass: there are 3 different constraint classes
the current constraint classifier "nonzeros" consists of 6 different classes
the current constraint classifier "constypes" consists of 4 different classes
the current constraint classifier "constypes according to miplib" consists of 5 different classes
dec_consclass found 109 new partialdecs
dec_densemasterconss found 1 new partialdec
dec_neighborhoodmaster found 1 new partialdec
the current varclass distribution includes 179 classes but only 9 are allowed for propagatePartialdec() of var class detector
POSTPROCESSING of decompositions. Added 14 new decomps.
Found 106 finished decompositions.
Measured running time per detector:
Detector postprocess worked on 14 finished decompositions and took a total time of 0.007
Detector consclass worked on 101 finished decompositions and took a total time of 0.016
Detector connectedbase worked on 105 finished decompositions and took a total time of 0.064
Detector varclass worked on 4 finished decompositions and took a total time of 0.001
Detection Time: 0.26

Reformulation

More details and the corresponding theory and implementation can be found on the page The GCG Pricing.

After decomposing the problem, GCG will apply a Dantzig-Wolfe (default) or Benders reformulation to the problem. With that formulation, the problem can in some cases be solved far more efficiently using Column Generation.

A Dantzig-Wolfe reformulation is applied to solve the original problem.
Chosen structure has 6 blocks and 20 linking constraints.
This decomposition has a maxwhite score of 0.826893.
Warning: Discretization with continuous variables is only an experimental feature.
Master problem is a set partitioning problem.
Matrix has 6 blocks, using 1 aggregated pricing problem.

Branch-and-Price-and-Cut

More details and the corresponding theory and implementation can be found on the page The GCG Pricing.

As a BP&C solver, GCG will apply Column Generation and apply cuts in every node of the branch-and-bound tree. Concerning the log, this is the part that Andrew already knows quite well. While most columns are not that important to him yet, there is one that is particularly, and that is the gap. It tells him how close GCG currently is to the optimal solution. It can also be read at the "dualbound" and "primalbound" columns.

time | node | left |SLP iter|MLP iter|LP it/n| mem |mdpt |ovars|mvars|ocons|mcons|mcuts| dualbound | primalbound | deg | gap
p 0.9s| 1 | 0 | 0 | 0 | - | 23M| 0 |2682 | 0 |2588 | 0 | 0 | 0.000000e+00 | 6.217321e+02 | -- | Inf
1.0s| 1 | 0 | 0 | 0 | - | 23M| 0 |2682 | 0 |2588 | 0 | 0 | 3.242201e+02 | 6.217321e+02 | -- | 91.76%
1.0s| 1 | 0 | 0 | 0 | - | 23M| 0 |2682 | 6 |2589 | 22 | 0 | 3.242201e+02 | 6.217321e+02 | 0.00%| 91.76%
1.0s| 1 | 0 | 0 | 0 | - | 23M| 0 |2682 | 6 |2589 | 22 | 0 | 3.242201e+02 | 6.217321e+02 | 0.00%| 91.76%
Starting reduced cost pricing...
*r14.2s| 1 | 0 | 19718 | 145 | - | 25M| 0 |2682 | 118 |2589 | 22 | 0 | 3.761693e+02 | 4.308847e+02 | 37.05%| 14.55%
14.2s| 1 | 0 | 19718 | 145 | - | 25M| 0 |2682 | 118 |2589 | 22 | 0 | 3.761693e+02 | 4.308847e+02 | 37.05%| 14.55%
39.1s| 1 | 0 | 90613 | 174 | - | 25M| 0 |2682 | 146 |2589 | 22 | 0 | 4.308847e+02 | 4.308847e+02 | 45.44%| 0.00%
39.1s| 1 | 0 | 90613 | 174 | - | 25M| 0 |2682 | 118 |2589 | 22 | 0 | 4.308847e+02 | 4.308847e+02 | -- | 0.00%

Getting the Solution

Apart from the optimal objective function value, Andrew also wants to see which vehicle is routed where. To do this, he can simply (after having solved the problem) type display solution and will be able to see what values each variable has in the optimal solution.

objective value: 430.884678
x#1#3#1 1 (obj:8.062258)
x#1#6#2 1 (obj:24.738634)
x#1#8#3 1 (obj:17.262677)
...

Depending on the definition of the model, this can have different meanings, but here, the most intuitive applies: if x#1#3#1 is 1, then \(x^1_{31}=1\), i.e. vehicle \(1\) drives from customer \(3\) to customer \(1\). Since driving there is not free, this comes at a cost of about 8.06 units for the objective function value.

Disclaimer

With this (not even that simple!) VRP instance, we were quite fortunate. SCIP does not seem to solve in reasonable time at all, and even commercial solvers such as gurobi are not faster than GCG. However, this is not always the case. It can happen that GCG performs significantly slower than other solvers not using Branch-and-Price. Thus, the main takeaway should be that you should always try column generation, since it can in many cases still be very beneficial to your solving performance. If you want to understand what the main factor of success of Branch-and-Price is, you should check out the page on Pricing Solvers.