Breaking symmetry in linear programming formulation

Jose M Sallan 2022-01-12 8 min read

An integer linear program (ILP) is symmetric if its variables can be permuted without changing the structure of the problem. Even for relatively modestly sized problems, ILPs with large symmetry groups are difficult to solve using traditional branch-and- bound or branch-and-cut algorithms (Margot, 2010). This is because once the algorithm finds an optimal solution, it needs to explore all its permutations until it knows for sure that has found the optimum. Finding an optimal solution for symmetrical problems can be easy, bot proving that it is optimal is not. We can tackle symmetry with slight perturbations of cost coefficients to break symmetry or fixing variables.

Here I will present a formulation of the bin packing problem with symmetry, to examine how can we eliminate symmetry with an alternative formulation. I will be using Rglpk to solve linear programming and use rbenchmark to compare times of execution.

library(Rglpk)
library(rbenchmark)

A formulation of the bin packing problem

In the bin packing problem (BPP), items of different volumes must be packed into equal bins of a fixed given volume in a way that minimizes the number of bins used. An instance of the bin packing problem is defined by the volume \(v_i\) of each of the \(n\) elements to pack, and by the volume \(V\) of each of the bins used.

To examine the performance of different BPP formulations, I will use the following instance:

vol <- c(0.61, 0.96, 0.95, 0.91, 0.13, 0.53, 0.53, 0.05, 0.65, 0.66)
V <- 1

A straigthforward formulation of the BPP can be obtained defining variables:

  • binaries \(x_{ij}\) equal to one if element \(i\) goes into bin \(j\) and zero otherwise.
  • binaries \(y_j\) equal to one if we use bin \(j\) and zero otherwise.

As we don’t know how many bins we will use, we set \(i = 1, \dots, n\) and \(j = 1, \dots, n\).

The objective function is the number of bins to use:

\[ \text{MIN } \sum_{j=1}^n y_j \]

As constraints, we first establish that each element must be in one bin:

\[ \sum_{i=1}^n x_{ij} = 1, \ \ \ i = 1, \dots, n \]

The second set of constraints makes that we can use each bin only when its binary variable is equal to one, taking into account that bins have capacity \(V\)

\[ \sum_{i=1}^n v_i x_{ij} \leq Vy_j, \ \ \ j = 1, \dots, n \]

To break this symmetry somewhat (although not completely) we can force that bin \(j\) cannot be activated if bin \(j-1\) has not been activated.

\[ y_{j} \leq y_{j-1} \ \ \ j = 2, \dots, n \]

The implementation of this formulation in MathProg for the instance is:

bpp1 <- 
"/* number of elements in bin */
param n; 

/* volume of bin */
param V; 

/* volume of each element */
param volume{1..n};

/* equals one if element i is in bin j */
var x{i in 1..n, j in 1..n}, binary;

/* equals one if we are using bin j */
var y{j in 1..n}, binary;

minimize bins: sum{j in 1..n} y[j];

/* each element is only in one bin */
s.t. element{i in 1..n}: sum{j in 1..n} x[i, j] = 1;

/* bin j is active if y[j] = 1 */
s.t. binscap{j in 1..n}: sum{i in 1..n} volume[i] * x[i, j] - V * y[j] <= 0;

/* to break symmetry, y[i] can be active only if y[i-1] is */
s.t. order{j in 2..n}: y[j] - y[j-1] <= 0;

data;

param n := 10;

param V := 1;

param volume := 1 0.61
                2 0.96
                3 0.95
                4 0.91
                5 0.13
                6 0.53
                7 0.53
                8 0.05
                9 0.65
                10 0.66;

end;"

The solution for this formulation is:

##  x[1,4]  x[2,7]  x[3,8]  x[4,2]  x[5,1]  x[6,5]  x[7,1]  x[8,2]  x[9,6] x[10,3] 
##       1       1       1       1       1       1       1       1       1       1 
##    y[1]    y[2]    y[3]    y[4]    y[5]    y[6]    y[7]    y[8] 
##       1       1       1       1       1       1       1       1

The same solution presented as data frame:

element bin
1 4
2 7
3 8
4 2
5 1
6 5
7 1
8 2
9 6
10 3

Thanks to the order constraints we are using only bins 1 to 8, but some symmetry persists, because bins can be relabelled. This second table corresponds with the same solution:

element bin
1 5
2 8
3 1
4 3
5 2
6 6
7 2
8 3
9 7
10 4

An alternative formulation avoiding symmetry

A way of supressing the symmetry of the formulation above is to identify each with the lowest-index item it contains. For instance, bin 1 can contain elements 1 or higher, bin 3 elements 3 or higher, and so on. To achieve this, we define variables \(u_{ij}\) that equal one if element \(j \geq i\) is in the bin with lowest-index element \(i\) and zero otherwise. For this formulation to work, we need that variables \(u_{ij}\) for \(j > i\) can be 1 only if \(u_{ii} = 1\).

For a bin packing problem of \(n=5\), we would have the following scheme of variables:

  • elements in bin 1: \(u_{11}\), \(u_{12}\), \(u_{13}\), \(u_{14}\), \(u_{15}\)
  • elements in bin 2: \(u_{22}\), \(u_{23}\), \(u_{24}\), \(u_{25}\)
  • elements in bin 3: \(u_{33}\), \(u_{34}\), \(u_{35}\)
  • elements in bin 4: \(u_{44}\), \(u_{45}\)
  • elements in bin 5: \(u_{55}\)

Note that these variables are defined for \(j \geq i\), so we have \(n \left( n+1\right)/2\) variables.

Let’s formulate the bin packing problem with these variables. As bin \(i\) is open when \(u_{ii} = 1\), the objective function is:

\[ \text{MIN } \sum_{i=1}^n u_{ii} \]

Each element \(j\) must be in one bin:

\[ \sum_{i=1}^j u_{ij} = 1, \ \ \ j = 1, \dots, n \]

And no elements can be placed in a bin that is not open:

\[ \sum_{j=i}^n v_j u_{ij} \leq Vu_{ii}, \ \ \ i = 1, \dots, n \]

The implementation of this formulation in MathProg for this instance is:

bpp2 <- 
"/* number of elements in bin */
param n;

/* volume of bin */
param V;

/* volume of each element */
param volume{1..n}; 

/* equals one if element j >=i is in the same bin as i */
var u{i in 1..n, j in i..n}, binary; 

minimize bins: sum{i in 1..n} u[i, i];

/* each element in one bin */
s.t. element{j in 1..n}: sum{i in 1..j} u[i, j] = 1; 

/* bin of element i is active if u[i, i] = 1 */
/* variables u[i, j] can be 1 only if u[i, i] = 1 */
s.t. bincap{i in 1..n}: sum{j in i..n} volume[j] * u[i, j] - V * u[i, i] <= 0; 

data;

param n := 10;

param V := 1;

param volume := 1 0.61
                2 0.96
                3 0.95
                4 0.91
                5 0.13
                6 0.53
                7 0.53
                8 0.05
                9 0.65
                10 0.66;

end;"

And the non-zero elements of the solution are:

##  u[1,1]  u[2,2]  u[3,3]  u[4,4]  u[5,5]  u[6,6]  u[7,7]  u[8,8]  u[5,9] u[8,10] 
##       1       1       1       1       1       1       1       1       1       1

Note that now are fixing the index of the label with the element of lower index. Each element is in its own bin, except elements 9 and 10, that are contained in the bins of elements 5 and 8, respectively.

Time of execution

This second formulation can be advantageous in terms of time of execution. To examine that, I will read each of the two formulations into model_bpp1 and model_bpp2, respectively:

writeLines(bpp1, "bpp.mod")
model_bpp1 <- Rglpk_read_file("bpp.mod", type = "MathProg", verbose = FALSE)
unlink("bpp.mod")

writeLines(bpp2, "bpp.mod")
model_bpp2 <- Rglpk_read_file("bpp.mod", type = "MathProg", verbose = FALSE)
unlink("bpp.mod")

Then, I am creating a function to solve each model, without elaborating on the function output:

solve_lpmod <- function(model){
  solve <- Rglpk_solve_LP(obj = model$objective,
                        mat = model$constraints[[1]],
                        dir = model$constraints[[2]],
                        rhs = model$constraints[[3]],
                        types = model$types,
                        bounds = model$bounds,
                        max = model$maximum)
  return(solve)
}

Finally, I am using rbenchmark to compare the performance of both models:

speed_test <- 
benchmark(solve_lpmod(model_bpp1), solve_lpmod(model_bpp2), columns=c('test', 'replications', 'elapsed', 'relative', 'user.self', 'sys.self'), replications = 10, order='elapsed')
speed_test
##                      test replications elapsed relative user.self sys.self
## 2 solve_lpmod(model_bpp2)           10   0.015      1.0     0.015        0
## 1 solve_lpmod(model_bpp1)           10   8.016    534.4     8.014        0

We observe that the second formulation is considerably faster: the second formulation is 534.4 times faster than the second.

Acknowledgements

Thanks to Mari Albareda for proofreading this post, suggesting a satisfactory definition of \(u_{ij}\) variables and further work on this formulation. All remaining errors are my own.

References

Session info

## R version 4.1.2 (2021-11-01)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Debian GNU/Linux 10 (buster)
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.8.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.8.0
## 
## locale:
##  [1] LC_CTYPE=es_ES.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=es_ES.UTF-8        LC_COLLATE=es_ES.UTF-8    
##  [5] LC_MONETARY=es_ES.UTF-8    LC_MESSAGES=es_ES.UTF-8   
##  [7] LC_PAPER=es_ES.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=es_ES.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] rbenchmark_1.0.0 Rglpk_0.6-4      slam_0.1-49      kableExtra_1.3.4
## 
## loaded via a namespace (and not attached):
##  [1] rstudioapi_0.13   knitr_1.33        xml2_1.3.2        magrittr_2.0.1   
##  [5] munsell_0.5.0     rvest_1.0.2       viridisLite_0.4.0 colorspace_2.0-1 
##  [9] R6_2.5.0          rlang_0.4.12      highr_0.9         stringr_1.4.0    
## [13] httr_1.4.2        tools_4.1.2       webshot_0.5.2     xfun_0.23        
## [17] jquerylib_0.1.4   systemfonts_1.0.2 htmltools_0.5.1.1 yaml_2.2.1       
## [21] digest_0.6.27     lifecycle_1.0.0   bookdown_0.24     sass_0.4.0       
## [25] glue_1.4.2        evaluate_0.14     rmarkdown_2.9     blogdown_1.5     
## [29] stringi_1.7.3     compiler_4.1.2    bslib_0.2.5.1     scales_1.1.1     
## [33] svglite_2.0.0     jsonlite_1.7.2