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!
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