R Binary Integer Optimization with Groups

Question

I am trying to get Rsolnp to constrain my parameters to binary integers or to decimals that are nearly the same (.999 is close enough to 1 for example).

I have three vectors of equal length (52), each of which will get multiplied by my binary parameter vector in my objective function.

library(Rsolnp)
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)

x is my parameter vector and below if my objective function.

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)
}

I essentially want 1 paramater in each group of 4 equal to 1 and the rest 0 and I'm not sure how to create the correct constraints for this but I believe I need to use these sum constraints in combination with another type of constraint as well. Here are my constraints:

eqn1=function(x){
  z1=sum(x[1:4])
  z2=sum(x[5:8])
  z3=sum(x[9:12])
  z4=sum(x[13:16])
  z5=sum(x[17:20])
  z6=sum(x[21:24])
  z7=sum(x[25:28])
  z8=sum(x[29:32])
  z9=sum(x[33:36])
  z10=sum(x[37:40])
  z11=sum(x[41:44])
  z12=sum(x[45:48])
  z13=sum(x[49:52])
  return(c(z1,z2,z3,z4,z5,z6,z7,z8,z9,z10,z11,z12,z13))
}

And finally, here is my function call:

opti<-solnp(pars=rep(1,52), fun = objective_function, eqfun = eqn1, eqB = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1), LB=rep(0,52))

Calling opti$pars returns my solution vector:

 [1] 7.199319e-01 2.800680e-01 6.015388e-08 4.886578e-10 5.540961e-01 4.459036e-01 2.906853e-07 4.635970e-08 5.389325e-01
[10] 4.610672e-01 2.979195e-07 3.651954e-08 6.228346e-01 3.771652e-01 1.980380e-07 3.348488e-09 5.389318e-01 4.610679e-01
[19] 2.979195e-07 3.651954e-08 5.820231e-01 4.179766e-01 2.099869e-07 2.624076e-08 5.389317e-01 4.610680e-01 2.979195e-07
[28] 3.651954e-08 6.499878e-01 3.500120e-01 1.959133e-07 1.059012e-08 6.249098e-01 3.750900e-01 2.588037e-07 1.752927e-08
[37] 6.249106e-01 3.750892e-01 2.588037e-07 1.752927e-08 6.095743e-01 3.904254e-01 2.741968e-07 2.233806e-08 6.095743e-01
[46] 3.904254e-01 2.741968e-07 2.233806e-08 5.679608e-01 4.320385e-01 6.821224e-07 3.997882e-08

As one can see the weight is getting split between multiple variables in each group of 4 instead of being forced onto just 1 with the rest being 0.

If this is not possible with this package could someone show me how to convert my objective function to work with other optimization packages? From what I have seen, they require the objective function to be converted to a vector of coefficients. Any help is appreciated. Thanks!

Looks like a problem with a nonlinear objective function, linear constraints and binary variables. This is called an MINLP problem (Mixed-Integer Non-Linear Programming). — Erwin Kalvelagen
– Erwin Kalvelagen, Commented Aug 25, 2020 at 19:31
@Christoph You should be able to run the code provided. Ideally my solution vector would look something like 0,1,0,0,1,0,0,0... with one in every four paramaters being equal to 1. — throwawaydisplayname
– throwawaydisplayname, Commented Aug 25, 2020 at 20:54

Erwin Kalvelagen · Accepted Answer · 2020-08-25 22:30:50Z

I tried with a few solvers. With MINLP solvers Couenne and Baron we can solve this directly. With Gurobi we need to decompose the objective into two quadratic parts. All these solvers give:

----    119 VARIABLE x.L  

i1  1.000,    i5  1.000,    i9  1.000,    i14 1.000,    i17 1.000,    i21 1.000,    i25 1.000,    i29 1.000
i34 1.000,    i38 1.000,    i41 1.000,    i46 1.000,    i49 1.000


----    119 VARIABLE z.L                   =     -889.346  obj

Zeroes are not printed here.

I used GAMS (commercial) but if you want to use free tools you can use Pyomo(Python) + Couenne. I am not sure about MINLP solvers for R, but Gurobi can be used from R.

Note that the group constraint is simply:

groups(g)..  sum(group(g,i),x(i)) =e= 1;

where g are the groups and group(g,i) is a 2d set with the mapping between groups and items.

For Gurobi you need to do something like (in pseudo code):

z1 = 1166 *  sum(i,x(i)*a(i)) / 2000        (linear)
z2 = ((sum(i, x(i)*b(i))) / 2100) + .05     (linear) 
z3 = ((sum(i, x(i)*c(i)))/1500) + 1.5       (linear) 
z23 = z2*z3                                 (non-convex quadratic)
obj = -z1*z23                               (non-convex quadratic)

and tell Gurobi to use the nonconvex MIQCP solver.

Sorry, no R code for this. But it may give you something to think about.

Alex Fleischer · Accepted Answer · 2020-08-30 07:22:49Z

within CPLEX you may try mathematical programming as Paul wrote, but you may also use Constraint Programming.

In OPL (CPLEX modeling language)

using CP;

execute
{
  cp.param.timelimit=5; // time limit 5 seconds
  
}

int n=52;
range r=1..n;

int a[r]=[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];
int b[r]=[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];
int c[r]=[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];

// decision variable 
 dvar boolean x[r];
 
 // objective
 dexpr float obj=
 
  -(1166 *  sum(i in r) (x[i]*a[i]) / 2000) *
    (((sum(i in r) (x[i]* b[i])) / 2100) + .05) *
    (((sum(i in r) (x[i]*c[i]))/1500) + 1.5);
    
minimize obj;

subject to
{
  // one and only one out of 4 is true
  forall(i in 1..n div 4) count(all(j in 1+(i-1)*4..4+(i-1)*4)x[j],1)==1;
  
}

gives

// solution with objective -889.3463
x = [1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1
         0 0 0 0 1 0 0 1 0 0 0 0 1 0 0 0 1 0 0 0 1 0 0 1 0 0 0 1 0 0 0 1 0
         0 0];

within 5 seconds

NB: You could call OPL CPLEX from R or rely on any other CPLEX API

And in python you can write the same

from docplex.cp.model import CpoModel

n=52

r=range(0,n)

a =[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 =[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 =[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]

mdl = CpoModel(name='x')

#decision variables
mdl.x = {i: mdl.integer_var(0,n,name="x"+str(i+1)) for i in r}

mdl.minimize(-1166 *  sum(mdl.x[i]*a[i] / 2000 for i in r) \
             *((sum(mdl.x[i]* b[i] / 2100  for i in r) +0.05)) \
             *((sum(mdl.x[i]*c[i]/1500  for i in r) +1.5)) )

for i in range(0,n // 4):
    mdl.add(1==sum( mdl.x[j] for j in range(i*4+0,i*4+4)))

msol=mdl.solve(TimeLimit=5)

# Dislay solution
for i in r:
    if (msol[mdl.x[i]]==1):
        print(i+1," ")

and that gives

! Best objective         : -889.3464 
 
1  
5  
9  
13  
17  
22  
25  
30  
34  
38  
41  
45  
49

@ Alex Fleischer : Great answer! Can you please take a look at my question? stackoverflow.com/questions/69319495/… thanks!

prubin · Accepted Answer · 2020-08-28 23:21:38Z

1

I set up an R notebook to solve (or try to solve) the problem as a mixed integer linear program, using CPLEX as the MIP solver and the Rcplex package as the interface to it. The results were unspectacular. After five minutes of grinding, CPLEX had a solution somewhat inferior to what Erwin got (-886.8748 v. his -889.346) with a gap over 146% (which, given Erwin's result, is mostly just the upper bound converging very slowly). I'm happy to share the notebook, which shows the linearization, but to use it you would need to have CPLEX installed.

Update: I have a second notebook, using the GA genetic algorithm package, that consistently gets close to Erwin's solution (and occasionally hits it) in under five seconds. The results are random, so rerunning may do better (or worse), and there is no proof of optimality.

edited Aug 28, 2020 at 23:21

answered Aug 28, 2020 at 22:18

prubin

3762 silver badges16 bronze badges

4 Comments

throwawaydisplayname Over a year ago

I'm currently working on doing this in python but if you could share both of those notebooks with me that would be awesome @prubin

Alex Fleischer Over a year ago

I updated my answer with the cplex python code I used

prubin Over a year ago

First notebook (using CPLEX): rubin.msu.domains/blog/groupselect.nb.html.

prubin Over a year ago

Second notebook (using GA): rubin.msu.domains/blog/groupselect2.nb.html.

Collectives™ on Stack Overflow

R Binary Integer Optimization with Groups

3 Answers 3

Comments

2 Comments

4 Comments

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

3 Answers 3

Comments

2 Comments

4 Comments

Your Answer

Sign up or log in

Post as a guest

Linked

Related