2

I have created an optimization model in Excel using the standard Solver and now would like to make a similar model in R as this will allow me to make larger models like this. Unfortunately, I have a bit of a hard time finding good examples that I can model my concept on. Therefore, I would like to ask you whether anybody is able to give me some hints on how to make a similar model in R.

I have uploaded my Excel sheet to http://dl.dropbox.com/u/9641130/R/Positioning%20Optimization%20R.xlsx

The basic idea is that cell B3 is maximized by changing max 8 cells in range E10:L19 to one. The B3 cell includes a sumproduct() of the range E10:L19 and a number of similar ranges.

I am looking forward to see some hints on how to build a similar model in R.

Thanks! Jochem

========

Update following Chase's suggestion

I would like to clarify my question a bit with some repoducable R code. This is roughly the same model as you will find in the Excel code above.

The initial set of matrices:

A <- as.matrix(structure(list(X0 = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L), X0.1 = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L), X0.2 = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L), X0.3 = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L), X0.4 = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L), X0.5 = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L), X0.6 = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L), X0.7 = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L)), .Names = c("X0", "X0.1", "X0.2", "X0.3", "X0.4", "X0.5", "X0.6", "X0.7"), class = "data.frame", row.names = c(NA, -9L)))
B <- as.matrix(structure(list(X1 = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), X1.1 = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), X1.2 = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), X1.3 = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), X.100000 = c(-100000L, -100000L, -100000L, -100000L, 1L, 1L, 1L, 1L, 1L), X.100000.1 = c(-100000L, -100000L, -100000L, -100000L, 1L, 1L, 1L, 1L, 1L), X.100000.2 = c(-100000L, -100000L, -100000L, -100000L, 1L, 1L, 1L, 1L, 1L), X.100000.3 = c(-100000L, -100000L, -100000L, -100000L, 1L, 1L, 1L, 1L, 1L)), .Names = c("X1", "X1.1", "X1.2", "X1.3", "X.100000", "X.100000.1", "X.100000.2", "X.100000.3"), class = "data.frame", row.names = c(NA, -9L)))
C <- as.matrix(structure(list(X1 = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), X1.1 = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), X1.2 = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), X1.3 = c(1L, 1L, 1L, 1L, 1L, -100000L, 1L, 1L, 1L), X1.4 = c(1L, 1L, 1L, 1L, 1L, -100000L, 1L, 1L, 1L), X1.5 = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), X1.6 = c(1L, 1L, 1L, 1L, 1L, -100000L, 1L, 1L, -100000L), X1.7 = c(1L, 1L, 1L, 1L, -100000L, -100000L, 1L, 1L, -100000L)), .Names = c("X1", "X1.1", "X1.2", "X1.3", "X1.4", "X1.5", "X1.6", "X1.7"), class = "data.frame", row.names = c(NA, -9L)))
D <- as.matrix(structure(list(X775 = c(385L, 1233L, 1067L, 5L, 730L, 1123L, 837L, 5L, 3087L), X704 = c(625L, 1338L, 804L, 110L, 659L, 1363L, 942L, -165L, 3350L), X704.1 = c(625L, 1338L, 804L, 110L, 659L, 1363L, 942L, -165L, 3350L), X944 = c(625L, 1263L, 898L, 35L, 899L, 1363L, 867L, -65L, 3110L), X775.1 = c(385L, 1233L, 1067L, 5L, 730L, 1123L, 837L, 5L, 3087L), X775.2 = c(385L, 1233L, 1067L, 5L, 730L, 1123L, 837L, 5L, 3087L), X944.1 = c(625L, 1263L, 898L, 35L, 899L, 1363L, 867L, -65L, 3110L), X944.2 = c(625L, 1263L, 898L, 35L, 899L, 1363L, 867L, -65L, 3110L)), .Names = c("X775", "X704", "X704.1", "X944", "X775.1", "X775.2", "X944.1", "X944.2"), class = "data.frame", row.names = c(NA, -9L)))

The result of the function sum(A*B*C*D) is currently 0. This is logical since in matrix A all cells have a value of 0. However, I would like to know with what formula I can maximize the value of the function sum(A*B*C*D). 

sum(A*B*C*D)
[1] 0

I want to do this by changing the values in Matrix A from 0 to 1. Moreover, the following constrains should be taken in consideration. 1. Each row can only include one cell with the value 1. 2. Each column can only include one cell with the value 1; this means that we can place a maximum of eight times the value 1 in Matrix A.

Does anybody have a suggestion on how to accomplis this?

Jochem
  • 3,295
  • 4
  • 30
  • 55
  • Don't have time to look into this now, but I asked a question about optimization with constraints [here](http://stackoverflow.com/questions/9817001/optimization-with-constraints) that may be of interest to you. If there are no constraints, then the code should be a bit more straight forward. If it's possible to include a snippet of data and show the expected output in the question, people will likely be more inclined to help, as downloading a xlsm may scare folks. Good tips on reproducible questions [here](http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example) – Chase May 29 '12 at 20:42
  • @Chase Thank you, I have updated the question a bit so there is no need to open the Excel file to answer the question. Thank you for your comments. So far I have been unable to work it out on the basis of the article that you linked in your comment. – Jochem May 30 '12 at 19:49

1 Answers1

1

Your decision variables (the cell values in A) are booleans (0 or 1) and your objective and constraints are linear functions of those variables, which puts you in a category of optimization problems called Mixed-Integer Linear Programmings. These problems can be solved using the Rglpk package for example. Here is my solution:

n1 <- nrow(A)
n2 <- ncol(A)

# objective coefficients
obj <- as.vector(B*C*D) # objective

# matrix of constraints weights
mat <- matrix(0, n1+n2, n1*n2)
for (i in 1:n1) {
   mat[i, ] <- as.numeric(row(A) == i)
}
for (j in 1:n2) {
   mat[n1+j, ] <- as.numeric(col(A) == j)
}

dir   <- rep("<=", n1+n2) # constraint directions
rhs   <- rep(1, n1+n2)    # constraints upper-bounds
types <- rep("B", n1*n2)  # variable types (boolean)

library(Rglpk)
opt <- Rglpk_solve_LP(obj, mat, dir, rhs, types,
                      max = TRUE, verbose = TRUE)
opt

# $optimum
# [1] 9950

# $solution
#  [1] 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0
# [39] 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0
# [77] 0 1 0 0

# $status
# [1] 0

opt.A <- matrix(opt$solution, n1, n2)
opt.A

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

Is 9950 the optimal value you were also getting for sum(A*B*C*D)? (I do not have Excel on this computer...)

flodel
  • 87,577
  • 21
  • 185
  • 223
  • Well the csv files and the excel sheet were not identical. However, I have recalculated the optimization in Excel (and uploaded the new sheet) and the result is 9950 as well; it only takes longer to compute than R ;-). Anyway, thank you for your solution! The code seems to work, I need to study it a little more in order to fully understand it! Cheers! – Jochem May 31 '12 at 06:26
  • Hi flodel, I have studied the code you provided a bit and I have some difficulty to understand what and why is happening where "matrix of constraints weights" are defined. Why lining up the rows and column in one matrix? – Jochem May 31 '12 at 19:35
  • @Jochem: constraints are defined in matrix form, via `mat`, `dir` and `rhs`. Take row `i` of `mat` for example: the constraint reads `sum(mat[i, ] * x) dir[i] rhs[i]`, or `sum(mat[i, ] * x) <= 1`. The first `n1` rows of `mat` are representing your "the sum per row is <= 1" constraints and the next `n2` rows are representing your "the sum per column is <= 1". – flodel Jun 01 '12 at 02:50
  • In my comment above, `x` is what you are solving for, it is a _vector_ of boolean variables (0 or 1) corresponding to the values in your matrix `A`. `x` has a length `n1*n2`; the first `n1` values are for the first column of `A`, the next `n1` values are for the second column, etc. Hope this puts you on the right track. – flodel Jun 01 '12 at 02:53