Optimization in R with constraint on the sum and type of optimization parameters

Question

I have written a function below that I would like to optimize

my_function = function(param, q, m){
  out = sum(-param*q + param*m)
  return(-out)
}

I am able to run the function and obtain an optimized result

> init = c(0,0,0)
> q = c(0.6, 0.14, 0.18)
> m = c(0, 2.5 , 4.2)
>
> nlminb(init, my_function, q=q, m=m, lower=c(0,0,0), upper=c(3,3,3))
$par
[1] 0 3 3

$objective
[1] -19.14

$convergence
[1] 0

$iterations
[1] 3

$evaluations
function gradient 
       4        9 

$message
[1] "both X-convergence and relative convergence (5)"

I would like to introduce the following constraints but I'm not sure how to do this

• The output parameters should be non-negative integers
• The parameters should sum up to some value k

Can someone inform me on how I can achieve this please?

Answer 1

1) Define a function proj such that for any input vector x the output vector y satisfies sum(y) = k. Then we have the following.

Note that this is a relaxation of the original problem where we have not applied the integer constraint; however, if the relaxed problem satisfies the constraint then it must be the solution to the original problem as well.

proj <- function(x, k = 3) k * x / sum(x)
obj <- function(x, ...) my_function(proj(x), ...)
out <- nlminb(c(1, 1, 1), obj, q = q, m = m, lower = 0)
str(out)
## List of 6
##  $ par        : num [1:3] 0 0 5.05
##  $ objective  : num -12.1
##  $ convergence: int 0
##  $ iterations : int 4
##  $ evaluations: Named int [1:2] 5 12
##   ..- attr(*, "names")= chr [1:2] "function" "gradient"
##  $ message    : chr "both X-convergence and relative convergence (5)"

proj(out$par) # solution
## [1] 0 0 3

2) Another approach is to use integer programming. This one does explicitly impose the integer constraint.

library(lpSolve)
res <- lp("min", q-m, t(rep(1, 3)), "=", 3, all.int = TRUE)
str(res)

giving the following (res$solution is the solution).

List of 28
 $ direction       : int 0
 $ x.count         : int 3
 $ objective       : num [1:3] 0.6 -2.36 -4.02
 $ const.count     : int 1
 $ constraints     : num [1:5, 1] 1 1 1 3 3
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:5] "" "" "" "const.dir.num" ...
  .. ..$ : NULL
 $ int.count       : int 3
 $ int.vec         : int [1:3] 1 2 3
 $ bin.count       : int 0
 $ binary.vec      : int 0
 $ num.bin.solns   : int 1
 $ objval          : num -12.1
 $ solution        : num [1:3] 0 0 3
 $ presolve        : int 0
 $ compute.sens    : int 0
 $ sens.coef.from  : num 0
 $ sens.coef.to    : num 0
 $ duals           : num 0
 $ duals.from      : num 0
 $ duals.to        : num 0
 $ scale           : int 196
 $ use.dense       : int 0
 $ dense.col       : int 0
 $ dense.val       : num 0
 $ dense.const.nrow: int 0
 $ dense.ctr       : num 0
 $ use.rw          : int 0
 $ tmp             : chr "Nobody will ever look at this"
 $ status          : int 0
 - attr(*, "class")= chr "lp"

Answer 2

You could try a brute-force grid search:

my_function <- function(param, q, m){
  out <- sum(-param*q + param*m)
  -out
}
q <- c(0.6, 0.14, 0.18)
m <- c(0, 2.5 , 4.2)

library("NMOF")
ans <- gridSearch(fun = my_function,
                  lower = c(0, 0, 0),
                  upper = c(3, 3, 3),
                  n = 4,  ## 4 levels from lower to upper: 0,1,2,3
                  q = q, m = m)

The answer is a list of all the possible combinations and their objective-function values:

ans
## $minfun
## [1] -19.14
## 
## $minlevels
## [1] 0 3 3
## 
## $values
##  [1]   0.00   0.60   1.20   1.80  -2.36  -1.76  -1.16  -0.56  -4.72  -4.12
## [11]  -3.52  -2.92  -7.08  -6.48  -5.88  -5.28  -4.02  -3.42  -2.82  -2.22
## [21]  -6.38  -5.78  -5.18  -4.58  -8.74  -8.14  -7.54  -6.94 -11.10 -10.50
## [31]  -9.90  -9.30  -8.04  -7.44  -6.84  -6.24 -10.40  -9.80  -9.20  -8.60
## [41] -12.76 -12.16 -11.56 -10.96 -15.12 -14.52 -13.92 -13.32 -12.06 -11.46
## [51] -10.86 -10.26 -14.42 -13.82 -13.22 -12.62 -16.78 -16.18 -15.58 -14.98
## [61] -19.14 -18.54 -17.94 -17.34
## 
## $levels
## $levels[[1]]
## [1] 0 0 0
## 
## $levels[[2]]
## [1] 1 0 0
## 
## $levels[[3]]
## [1] 2 0 0
## 
## .....
## 
## $levels[[64]]
## [1] 3 3 3

The levels are non-negative integers, but their sum is unconstrained. To add a sum constraint, either check in the objective function and return a large value if the particular solution violates the constraint (i.e. the solution gets marked as bad). Or filter the results; for instance, suppose the sum should be 2:

valid <- sapply(ans$levels, sum) == 2
ans$values[valid]
## [1]  1.20 -1.76 -4.72 -3.42 -6.38 -8.04
ans$levels[valid]
## [[1]]
## [1] 2 0 0
## 
## [[2]]
## [1] 1 1 0
## 
## [[3]]
## [1] 0 2 0
## 
## [[4]]
## [1] 1 0 1
## 
## [[5]]
## [1] 0 1 1
## 
## [[6]]
## [1] 0 0 2

Then keep only the best of the valid solutions.

best <- which.min(ans$values[valid])
ans$values[valid][best]
## [1] -8.04
ans$levels[valid][best]
## [[1]]
## [1] 0 0 2