How to use R to solve a program with piecewise constant objective function?

Question

I would like to solve a minimization problem using R that has a piecewise constant objective function. The idea is that for lower values of my (integer) decision variable x higher penalty costs are incurred than for higher values. I want to minimize total penalty costs, given some constraints.

So, my program looks as follows:

min   P(x)
s.t.  A x <= b
        x >= 0

I have a problem programming the piecewise constant objective function P(x), where P(x) is the sum of all elements of the vector x. I know it cannot be used in combination with the function lp() from the linprog library. However, I cannot find a way to do this without specifying an awful lot of extra variables. Furthermore, an extensive internet search did not provide any helpful ideas.

Let me give an example of how this function P looks like

[,1] [,2] [,3] [,4] [,5] [,6]
  21   11    9    9    0    0
  45   28   17   17    6    0
  14    0    0    0    0    0
  17   11   11    5    5    0
  26   11    0    0    0    0
  38   18   18   13    0    0

This should be read in the following way: if x1=2 a penalty cost of 11 is incurred. If x6=4, a penalty cost of 13 is incurred. That is, for x=c(2, ..., 4), we have that P=c(11, ..., 13) and total penalty costs (objective value) is sum(11, ..., 13).

My matrix A (it is totally unimodular) and vector b look as follows .

A <- matrix(c(1,0,1,0,0,1,1,0,0,1,0,0,0,1,0,1,1,0,0,1,0,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1),nrow=6) b <- c(4,5,1,5,2,4).

So, my question is:

How can I find the minimum of a piecewise constant objective function using R?

it seems as if you're saying that for x = c(2, , , , , ,4), the value of P is c(11, , , , , 13) so P(x) is a vector-valued function. How would you calculate "total penalty costs" as a scalar? — WaltS, Apr 24 '15 at 13:56
So you need a scalar to minimize. So how do you weight the elements of P to make a scalar? — WaltS, Apr 24 '15 at 14:01
Oh, I am sorry to not have included this in my question originally. I simply take the sum of the elements of P. — Braamstruik, Apr 24 '15 at 14:13

Kota Mori · Answer 1 · 2015-12-15T15:44:33.447

I got two approaches for your question. Your problem seems to be an integer programming where x1 ... x6 may take a value in {1, 2, 3, 4, 5, 6}.

Pmat <- matrix(c(25, 11, 9, 9, 0, 0,
                 45, 28, 17, 17, 6, 0,
                 14, 0, 0, 0, 0, 0,
                 17, 11, 11, 5, 5, 0,
                 26, 11, 0, 0, 0, 0,
                 38, 18, 18, 13, 0, 0),
                byrow = TRUE, nrow = 6, ncol = 6)
A <- matrix(c(1,0,1,0,0,1,1,0,0,1,0,0,0,1,0,1,1,0,0,1,0,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1),nrow=6)
b <- c(4,5,1,5,2,4)

Brute force

Since you have 6 such variables, there are only 6^6 = 46656 candidates. Hence it won't take that long time to solve this in brute force; for each combination check if it satisfies the constraint and the value.

d <- expand.grid(x1 = 1:6, x2 = 1:6, x3 = 1:6, x4 = 1:6, x5 =1:6, x6 = 1:6)
condition <- logical(nrow(d))
value <- numeric(nrow(d))
for (i in seq(nrow(d))) {
  x <- as.matrix(as.numeric(d[i, 1:6]), ncol = 1)
  d$condition[i] <- all(A %*% x <= b)
  d$value[i] <- sum(Pmat[cbind(1:6, x)])
}
d <- d[order(d$value, decreasing = FALSE), ]
d <- d[order(d$condition, decreasing = TRUE), ]
print(head(d))

This gets

#     x1 x2 x3 x4 x5 x6 condition value
#7789  1  3  1  1  1  2      TRUE   117
#7783  1  2  1  1  1  2      TRUE   128
#229   1  3  1  2  1  1      TRUE   131
#13    1  3  1  1  1  1      TRUE   137
#223   1  2  1  2  1  1      TRUE   142
#7777  1  1  1  1  1  2      TRUE   145

So the solution is [1, 3, 1, 1, 1, 2] with the value 117.

LP formulation

The brute force is not efficient and you obviously are seeking a way to solve using a LP solver. In order to formulate the LP for your problem, I would redefine a set of binary variables as below (this may be what you call "an awful lot of extra variables", though):

y_{ij} = 1(x_i = j)  for each i = 1,...,6 and j = 1,...,6

Here, 1() is a function that takes 1 if and only if the statement inside is true.

Then vectorize these variables as

y = [y_{11}, y_{12}, ..., y_{16}, y_{21}, ..., y_{66}]

This is the new control variable for the LP. Note that the size of y is 36.

You then need to vectorize the P matrix into the weight vector for y. You also need to convert the matrix A. For example, if the original constraint is

x_1 + x_2 <= b

Then this is equivalent to

  y_{11} + 2*y_{12} + 3*y_{13} + 4*y_{14} + 5*y_{15} + 6*y_{16} 
+ y_{21} + 2*y_{22} + 3*y_{23} + 4*y_{24} + 5*y_{25} + 6*y_{26}  <= b

The Kronecker product is convenient for this conversion.

Finally, for each i, there must be only one j such that y_{ij} = 1. Hence you need

y_{i1} + y_{i2} + ... + y_{i6} = 1  for each i.

The Kronecker product is helpful for this too.

library(lpSolve)
Pvec <- as.numeric(t(Pmat))
A1 <- kronecker(A, matrix(1:6, nrow = 1))
b1 <- b
A2 <- kronecker(diag(6), matrix(rep(1, 6), nrow = 1))
b2 <- rep(1, 6)
o <- lp(direction = "min",
        objective.in = Pvec, const.mat = rbind(A1, A2),
        const.dir = c(rep("<=", 6), rep("=", 6)), const.rhs = c(b1, b2),
        all.bin = TRUE)
# show the variable names
tmp <- expand.grid(1:6, 1:6)
vname <- paste("x", tmp[,2], tmp[,1], sep = "")
names(o$solution) <- vname
print(o)
print(o$solution)

This yields the same solution as the brute force.

Success: the objective function is 117 

x11 x12 x13 x14 x15 x16 x21 x22 x23 x24 x25 x26 x31 x32 x33 x34 x35 x36 
  1   0   0   0   0   0   0   0   1   0   0   0   1   0   0   0   0   0 
x41 x42 x43 x44 x45 x46 x51 x52 x53 x54 x55 x56 x61 x62 x63 x64 x65 x66 
  1   0   0   0   0   0   1   0   0   0   0   0   0   1   0   0   0   0

How to use R to solve a program with piecewise constant objective function?

1 Answers1

Brute force

LP formulation