0

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?

NM_
  • 1,887
  • 3
  • 12
  • 27

2 Answers2

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"
G. Grothendieck
  • 254,981
  • 17
  • 203
  • 341
  • Can you comment on how this approach will ensure that the parameters are non-negative integers and that they sum to k? In this particular case it seems to find an integer solution, but what guarantee is there? – user2974951 Oct 26 '21 at 11:17
  • It does not guarantee that the result is integer but in this particular case it was. – G. Grothendieck Oct 26 '21 at 11:17
  • So this doesn't exactly solve the OP's problem? In general. – user2974951 Oct 26 '21 at 11:18
  • I am not sure whether it would solve it in general or not. There do exist conditions under which linear program relaxations will provide an integer solution. At any rate, It did solve it in this case as the problem without the integer constraint satisfies the constraint so the integer constraint was not active. I have added a second approach just in case. – G. Grothendieck Oct 26 '21 at 11:45
1

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
Enrico Schumann
  • 1,278
  • 7
  • 8