Group Selection

This notebook contains a mixed integer linear programming (MILP) model for a problem posed on Stack Overflow. The problem involves 52 binary variables in 13 groups of four each. Exactly one variable from each group must be set to 1, with the objective of minimizing a cubic function of the variables. In a blog post, Erwin Kalvelagen compared MINLP and MIQCP models, using GAMS for the modeling and four different solvers. Here I will linearize the problem and solve it using CPLEX.

Set up

We will need the Rcplex package to interface with CPLEX and the Matrix package for sparse matrices (given the size of the resulting MILP model).

library(Matrix)
library(Rcplex)

Next, we load coefficient vectors a, b and c from the original SO post.

a <- c(251, 179, 215, 251, 63, 45, 54, 63, 47, 34, 40, 47, 141, 101, 121, 141, 47, 34, 40, 47, 94, 67, 81, 94, 47, 34, 40, 47, 157, 108, 133, 157, 126, 85, 106, 126, 126, 85, 106, 126, 110, 74, 92, 110, 110, 74, 92, 110, 63, 40, 52, 63)
b <- c(179, 251, 215, 0, 45, 63, 54, 0, 34, 47, 40, 0, 101, 141, 121, 0, 34, 47, 40, 0, 67, 94, 81, 0, 34, 47, 40, 0, 108, 157, 133, 0, 85, 126, 106, 0, 85, 126, 106, 0, 74, 110, 92, 0, 74, 110, 92, 0, 40, 63, 52, 0)
c <- c(179, 179, 118, 179, 45, 45, 30, 45, 34, 34, 22, 34, 101, 101, 67, 101, 34, 34, 22, 34, 67, 67, 44, 67, 34, 34, 22, 34, 108, 108, 71, 108, 85, 85, 56, 85, 85, 85, 56, 85, 74, 74, 49, 74, 74, 74, 49, 74, 40, 40, 27, 40)

Variables

Let \(y_{ij} = x_i * x_j\) and \(z_{ijk} = x_i * x_j * x_k = x_i * y_{jk}\) for all combinations of \(i,j,k\). Although \(y_{ij}\) and \(z_{ijk}\) can only be 0 or 1, they need not be declared binary. The constraints that define them will restrict them to being 0 or 1.

Note that, due to \(x\) being binary, \(y_{ii} = x_i^2 = x_i\) and \(z_{iii} = x_i^3 = x_i\). Also, \(y_{ij} = y_{ji}\) and \(z_{ijk} = z_{ikj} = \dots = z_{kji}\) since multiplication is commutative. Moreover, \(z_{iij} = z_{iji} = z_{jii} = y_{ij}\) for all \(i\neq j\). To exploit all that requires getting a bit “creative” about indexing variables.

Since Rcplex works with matrix representation, we need to “flatten” the variable matrices in a single vector. The variables will be listed in the order \(x\), \(y\), \(z\) with \(y\) and \(z\) sequenced so that the last index increases fastest and the first index slowest. To facilitate moving between algebraic and matrix representations, we will use a set of arrays that maps tuples of \(i\), \(j\) and \(k\) to 1-based indices in the master variable vector.

For convenience, we assign a label to the number of binary variables.

nx <- length(a)      # number of x variables (52)

# Index array for subscripts of x (which do not need modification).
ix <- 1:nx
# Index array for subscripts of y (which need to be mapped after the x entries).
t <- nx + 1  # next variable index
iy <- matrix(0, nrow = nx, ncol = nx)
for (i in 1:nx) {
  for (j in i:nx) {
    if (i == j) {
      iy[i, j] <- ix[i]
    } else {
      iy[i, j] <- t
      iy[j, i] <- t
      t <- t + 1
    }
  }
}
# Index array for subscripts of z (which need to be mapped after the y entries).
iz <- array(0, dim = c(nx, nx, nx))
for (i in 1:nx) {
  for (j in i:nx) {
    for (k in j: nx) {
      if (i == j) {
        if (j == k) {
          # All indices the same: reduces to an x variable.
          q <- ix[i]
        } else {
          # First two indices the same: reduces to a y variable.
          q <- iy[i, k]
        }
      } else if (j == k || i == k) {
        # Third index matches one of the others: reduces to a y variable.
        q <- iy[i, j]
      } else {
        # All indices are distinct.
        q <- t
        t <- t + 1
      }
      # All permutations of the indices get the same index.
      iz[i, j, k] <- q
      iz[i, k, j] <- q
      iz[j, i, k] <- q
      iz[j, k, i] <- q
      iz[k, i, j] <- q
      iz[k, j, i] <- q
    }
  }
}

For convenience, we’ll assign a name to the total number of distinct variables.

nvar <- nx + choose(nx, 2) + choose(nx, 3)

One of the arguments to Rcplex will be a vector indicating the type (continuous, integer, binary) of each variable. We can use scalars for the lower and upper bounds of the variables, since all variables are between 0 and 1.

# Start with all continuous variables.
vtype <- rep.int("C", nvar)
# The x variables (only) are binary.
vtype[1:nx] <- "B"
lb <- 0  # lower bound of all variables
ub <- 1  # upper bound of all variables

Objective function

The original post minimized the negative of a product. For a bit more simplicity, we will maximize the product without negating it. The original objective function has the following form (where the \(x_i\) are the binary variables): \[ \left(\alpha\sum_{i}a_{i}x_{i}\right)\times\left(\beta\sum_{j}b_{j}x_{j}+\beta_{0}\right)\times\left(\gamma\sum_{k}c_{k}x_{k}+\gamma_{0}\right) \] where \(\alpha = 1166/2000\), \(\beta = 1/2100\), \(\beta_0 = 0.05\), \(\gamma = 1/1500\) and \(\gamma_0 = 1.5\). We store those values here.

alpha <- 1166/2000
beta <- 1/2100
beta0 <- 0.05
gamma <- 1/1500
gamma0 <- 1.5

We can rewrite the objective function in linear form as follows: \[ \sum_{i,j,k}\alpha\beta\gamma a_{i}b_{j}c_{k}z_{ijk}+\sum_{i,j}\left(\alpha\beta\gamma_{0}a_{i}b_{j}+\alpha\beta_{0}\gamma a_{i}c_{j}\right)y_{ij}+\sum_{i}\alpha\beta_{0}\gamma_{0}a_{i}x_{i}. \]

To cut down on repetitive multiplications, we can assign some of the products to variables.

abg <- alpha * beta * gamma
abg0 <- alpha * beta * gamma0
ab0g <- alpha * beta0 * gamma
ab0g0 <- alpha * beta0 * gamma0

Now we build the objective coefficient vector.

# Start with a vector of zeros.
obj <- rep.int(0, nvar)
# Fill in the x coefficients.
obj[ix[1]:ix[nx]] <- ab0g0 * a
# Fill in the y coefficients using a loop. Note that we add to the existing values because the coefficient of y_{ii} adds to that of x_i and the coefficient of y_{ij} with j < i adds to that of y_{ji}.
for (i in 1:nx) {
  for (j in 1:nx) {
    obj[iy[i, j]] <- obj[iy[i, j]] + abg0 * a[i] * b[j] + ab0g * a[i] * c[j]
  }
}
# Fill in the z coefficients similarly.
for (i in 1:nx) {
  for (j in 1:nx) {
    for (k in 1:nx) {
      obj[iz[i, j, k]] <- obj[iz[i, j, k]] + abg * a[i] * b[j] * c[k]
    }
  }
}

Constraints

We will use a sparse matrix for the constraint coefficient matrix. To do that, we create three vectors of equal length, one (rindex) to hold row indices of nonzero entries, one (cindex) to hold column indices, and one (entry) to hold the entries themselves.

rindex <- c()
cindex <- c()
entry <- c()

We will use vector rhs for the right hand sides and vector sense for the constraint directions.

rhs <- c()
sense <- c()

The constraint on \(x\) is that in each consecutive group of four entries, exactly one must be 1 and the other three 0. So \(x_i + x_{i+1} + x_{i+2} + x_{i+3} = 1\) for \(i=1,5,9,\dots,49\).

r <- 0 # last row index used
for (i in seq(from = 1, to = nx, by = 4)) {
  r <- r + 1 # row index for the constraint
  # Create four entries (same row, consecutive columns) with value 1.
  rindex <- c(rindex, rep.int(r, 4))
  cindex <- c(cindex, ix[i]:ix[i + 3])
  entry <- c(entry, rep.int(1, 4))
  # The right hand side is 1.
  rhs <- c(rhs, 1)
  # The constraint is an equation.
  sense <- c(sense, "E")
}

To define the \(y\) variables, we assert the following constraints for each \(i\) and \(j\) with \(i<j\): \[ y_{ij} \le x_{i}\\ y_{ij} \le x_{j}\\ y_{ij} \ge x_{i}+x_{j}-1. \]

# Create vectors to hold the indices and coefficients. Each (i,j) pair with i < j generates three constraints with a total of seven coefficients.
ncoef <- 7 * choose(nx, 2)
rtemp <- rep.int(0, ncoef)
ctemp <- rep.int(0, ncoef)
etemp <- rep.int(0, ncoef)
t <- 1  # index for next entry
for (i in 1:(nx - 1)) {
  for (j in (i + 1):nx) {
    r <- r + 1
    rtemp[t:(t + 1)] <- r
    ctemp[t] <- iy[i, j]
    ctemp[t + 1] <- ix[i]
    etemp[t] <- 1
    etemp[t + 1] <- -1
    t <- t + 2
    r <- r + 1
    rtemp[t:(t + 1)] <- r
    ctemp[t] <- iy[i, j]
    ctemp[t + 1] <- ix[j]
    etemp[t] <- 1
    etemp[t + 1] <- -1
    t <- t + 2
    r <- r + 1
    rtemp[t:(t + 2)] <- r
    ctemp[t] <- iy[i, j]
    ctemp[t + 1] <- ix[i]
    ctemp[t + 2] <- ix[j]
    etemp[t] <- 1
    etemp[t + 1] <- -1
    etemp[t + 2] <- -1
    t <- t + 3
  }
}
# Merge the coefficients into the existing arrays.
rindex <- c(rindex, rtemp)
cindex <- c(cindex, ctemp)
entry <- c(entry, etemp)
# Each (i, j) entry has three constraints with the pattern "L", "L", "G".
npair <- choose(nx, 2)  # number of i, j pairs with i < j
sense <- c(sense, rep.int(c("L", "L", "G"), npair))
# Right hand sides for the three constraints are 0, 0, -1.
rhs <- c(rhs, rep.int(c(0, 0, -1), npair))

To define the \(z\) variables, we assert the following constraints for all combinations of \(i < j < k\): \[ z_{ijk} \le x_i \\ z_{ijk} \le y_{jk} \\ z_{ijk} \ge x_i + y_{jk} - 1. \]

# Create vectors to hold the indices and coefficients. Each (i,j,k) triplet generates three constraints with a total of seven coefficients.
ncoef <- 7 * nx^2
rtemp <- rep.int(0, ncoef)
ctemp <- rep.int(0, ncoef)
etemp <- rep.int(0, ncoef)
t <- 1  # index for next entry
for (i in 1:(nx- 2)) {
  for (j in (i + 1):(nx - 1)) {
    for (k in (j + 1):nx) {
      r <- r + 1
      rtemp[t:(t + 1)] <- r
      ctemp[t] <- iz[i, j, k]
      ctemp[t + 1] <- ix[i]
      etemp[t] <- 1
      etemp[t + 1] <- -1
      t <- t + 2
      r <- r + 1
      rtemp[t:(t + 1)] <- r
      ctemp[t] <- iz[i, j, k]
      ctemp[t + 1] <- iy[j, k]
      etemp[t] <- 1
      etemp[t + 1] <- -1
      t <- t + 2
      r <- r + 1
      rtemp[t:(t + 2)] <- r
      ctemp[t] <- iz[i, j, k]
      ctemp[t + 1] <- ix[i]
      ctemp[t + 2] <- iy[j, k]
      etemp[t] <- 1
      etemp[t + 1] <- -1
      etemp[t + 2] <- -1
      t <- t + 3
    }
  }
}
# Merge the coefficients into the existing arrays.
rindex <- c(rindex, rtemp)
cindex <- c(cindex, ctemp)
entry <- c(entry, etemp)
# Each (i, j, k) entry has three constraints with the pattern "L", "L", "G".
ntriples <- choose(nx, 3)
sense <- c(sense, rep.int(c("L", "L", "G"), ntriples))
# Right hand sides for the three constraints are 0, 0, -1.
rhs <- c(rhs, rep.int(c(0, 0, -1), ntriples))

Solution

We need to convert the constraint coefficients into a sparse matrix.

A <- sparseMatrix(i = rindex, j = cindex, x = entry)

We also need to set a time limit, since it is difficult to stop the CPLEX solver once it starts.

tl <- 300 # time limit in seconds

Now we call Rcplex to solve (hopefully) the model.

solution <- Rcplex(obj, A, rhs, lb = 0, ub = 1, objsense = "max", sense = sense, vtype = vtype, control = list(tilim = tl))

Once the solver is done, we can release it.

Rcplex.close()

We can inspect the results from CPLEX.

cat("The objective value of the incumbent solution is ", solution$obj, ".\n", "Incumbent solution:\n")
which(solution$xopt[1:nx] == 1)

Finally, we can use the original poster’s objective function (with the leading minus sign removed) to compare results and confirm that we have the correct objective value for the MIP solution.

objective_function <- function(x) {
  (1166 *  sum(x[1:52] * a) / 2000) *
  (((sum(x[1:52] * b)) / 2100) + .05) *
  (((sum(x[1:52] * c))/1500) + 1.5)
}
objective_function(solution$xopt)