Optimizing Non-Linear Function using nloptr - Unexpected Results

0
votes

When I try to maximize my objective function using nloptr, it is terminating on the initial values that I have set, which I think is very unlikely to return the maximum.

For some background, I am trying to allocate 80M of investment (x) between 114 products in a way that maximizes revenue (y). Predicted revenue for each product and investment level is based off a two-part model i.e. it is the result of multiplying the predicted probability of any revenue being generated, by the predicted level of revenue, conditional on revenue being generated. Each line of the objective function represents the fitted model for each product.

library(nloptr)

# objective function
eval_f0 <- function(x){
  return(-(plogis(0.872 + 0.0000471*x[1]-1.04483)*exp(8.459+0.468*log(x[1])+1.07743) +
             plogis(0.872 + 0.0000471*x[2]-1.04483)*exp(8.459+0.468*log(x[2])+1.07743) +
             plogis(0.872 + 0.0000471*x[3]-1.04483)*exp(8.459+0.468*log(x[3])+1.07743) +
             plogis(0.872 + 0.0000471*x[4]-1.04483)*exp(8.459+0.468*log(x[4])+1.07743) +
             plogis(0.872 + 0.0000471*x[5]-1.04483)*exp(8.459+0.468*log(x[5])+1.07743) +
             plogis(0.872 + 0.0000471*x[6]-1.04483)*exp(8.459+0.468*log(x[6])+1.07743) +
             plogis(0.872 + 0.0000471*x[7]-1.04483)*exp(8.459+0.468*log(x[7])+1.07743) +
             plogis(0.872 + 0.0000471*x[8]-1.04483)*exp(8.459+0.468*log(x[8])+1.07743) +
             plogis(0.872 + 0.0000471*x[9]-1.04483)*exp(8.459+0.468*log(x[9])+1.07743) +
             plogis(0.872 + 0.0000471*x[10]-1.04483)*exp(8.459+0.468*log(x[10])+1.07743) +
             plogis(0.872 + 0.0000471*x[11]-1.04483)*exp(8.459+0.468*log(x[11])+1.07743) +
             plogis(0.872 + 0.0000471*x[12]-1.04483)*exp(8.459+0.468*log(x[12])+1.07743) +
             plogis(0.872 + 0.0000471*x[13]-1.04483)*exp(8.459+0.468*log(x[13])+1.07743) +
             plogis(0.872 + 0.0000471*x[14]-1.04483)*exp(8.459+0.468*log(x[14])+1.07743) +
             plogis(0.872 + 0.0000471*x[15]-1.04483)*exp(8.459+0.468*log(x[15])+1.07743) +
             plogis(0.872 + 0.0000471*x[16]-1.04483)*exp(8.459+0.468*log(x[16])+1.07743) +
             plogis(0.872 + 0.0000471*x[17]-1.04483)*exp(8.459+0.468*log(x[17])+1.07743) +
             plogis(0.872 + 0.0000471*x[18]-1.04483)*exp(8.459+0.468*log(x[18])+1.07743) +
             plogis(0.872 + 0.0000471*x[19]-1.04483)*exp(8.459+0.468*log(x[19])+1.07743) +
             plogis(0.872 + 0.0000471*x[20]-1.04483)*exp(8.459+0.468*log(x[20])+1.07743) +
             plogis(0.872 + 0.0000471*x[21]-1.04483)*exp(8.459+0.468*log(x[21])+1.07743) +
             plogis(0.872 + 0.0000471*x[22]-1.04483)*exp(8.459+0.468*log(x[22])+1.07743) +
             plogis(0.872 + 0.0000471*x[23]-1.04483)*exp(8.459+0.468*log(x[23])+1.07743) +
             plogis(0.872 + 0.0000471*x[24]-1.04483)*exp(8.459+0.468*log(x[24])+1.07743) +
             plogis(0.872 + 0.0000471*x[25]-1.04483)*exp(8.459+0.468*log(x[25])+1.07743) +
             plogis(0.872 + 0.0000471*x[26]-1.04483)*exp(8.459+0.468*log(x[26])+1.07743) +
             plogis(0.872 + 0.0000471*x[27]-1.04483)*exp(8.459+0.468*log(x[27])+1.07743) +
             plogis(0.872 + 0.0000471*x[28]-1.04483)*exp(8.459+0.468*log(x[28])+1.07743) +
             plogis(0.872 + 0.0000471*x[29]-1.04483)*exp(8.459+0.468*log(x[29])+1.07743) +
             plogis(0.872 + 0.0000471*x[30]-1.04483)*exp(8.459+0.468*log(x[30])+1.07743) +
             plogis(0.872 + 0.0000471*x[31]-1.04483)*exp(8.459+0.468*log(x[31])+1.07743) +
             plogis(0.872 + 0.0000471*x[32]-1.04483)*exp(8.459+0.468*log(x[32])+1.07743) +
             plogis(0.872 + 0.0000471*x[33]-1.04483)*exp(8.459+0.468*log(x[33])+1.07743) +
             plogis(0.872 + 0.0000471*x[34]-1.04483)*exp(8.459+0.468*log(x[34])+1.07743) +
             plogis(0.872 + 0.0000471*x[35]-1.04483+0.289)*exp(8.459+0.468*log(x[35])+1.07743-0.83) +
             plogis(0.872 + 0.0000471*x[36]-1.04483+0.289)*exp(8.459+0.468*log(x[36])+1.07743-0.83) +
             plogis(0.872 + 0.0000471*x[37]-1.04483+0.289)*exp(8.459+0.468*log(x[37])+1.07743-0.83) +
             plogis(0.872 + 0.0000471*x[38]-1.04483+0.289)*exp(8.459+0.468*log(x[38])+1.07743-0.83) +
             plogis(0.872 + 0.0000471*x[39]-1.04483+0.289)*exp(8.459+0.468*log(x[39])+1.07743-0.83) +
             plogis(0.872 + 0.0000471*x[40]-1.04483+0.289)*exp(8.459+0.468*log(x[40])+1.07743-0.83) +
             plogis(0.872 + 0.0000471*x[41]-1.04483+0.289)*exp(8.459+0.468*log(x[41])+1.07743-0.83) +
             plogis(0.872 + 0.0000471*x[42]-1.04483+0.289)*exp(8.459+0.468*log(x[42])+1.07743-0.83) +
             plogis(0.872 + 0.0000471*x[43]-1.04483+0.289)*exp(8.459+0.468*log(x[43])+1.07743-0.83) +
             plogis(0.872 + 0.0000471*x[44]-1.04483+0.289)*exp(8.459+0.468*log(x[44])+1.07743-0.83) +
             plogis(0.872 + 0.0000471*x[45]-1.04483+0.289)*exp(8.459+0.468*log(x[45])+1.07743-0.83) +
             plogis(0.872 + 0.0000471*x[46]-1.04483+0.289)*exp(8.459+0.468*log(x[46])+1.07743-0.83) +
             plogis(0.872 + 0.0000471*x[47]-1.04483+0.289)*exp(8.459+0.468*log(x[47])+1.07743-0.83) +
             plogis(0.872 + 0.0000471*x[48]-1.04483+0.289)*exp(8.459+0.468*log(x[48])+1.07743-0.83) +
             plogis(0.872 + 0.0000471*x[49]-1.04483+0.289)*exp(8.459+0.468*log(x[49])+1.07743-0.83) +
             plogis(0.872 + 0.0000471*x[50]-1.04483+0.289)*exp(8.459+0.468*log(x[50])+1.07743-0.83) +
             plogis(0.872 + 0.0000471*x[51]-1.04483+0.289)*exp(8.459+0.468*log(x[51])+1.07743-0.83) +
             plogis(0.872 + 0.0000471*x[52]-1.04483+0.289)*exp(8.459+0.468*log(x[52])+1.07743-0.83) +
             plogis(0.872 + 0.0000471*x[53]-1.04483+0.289)*exp(8.459+0.468*log(x[53])+1.07743-0.83) +
             plogis(0.872 + 0.0000471*x[54]-1.04483+0.289)*exp(8.459+0.468*log(x[54])+1.07743-0.83) +
             plogis(0.872 + 0.0000471*x[55]-1.04483+0.289)*exp(8.459+0.468*log(x[55])+1.07743-0.83) +
             plogis(0.872 + 0.0000471*x[56]-1.04483+0.289)*exp(8.459+0.468*log(x[56])+1.07743-0.83) +
             plogis(0.872 + 0.0000471*x[57]-1.04483+0.289)*exp(8.459+0.468*log(x[57])+1.07743-0.83) +
             plogis(0.872 + 0.0000471*x[58]-1.04483+0.289)*exp(8.459+0.468*log(x[58])+1.07743-0.83) +
             plogis(0.872 + 0.0000471*x[59]-1.04483+0.289)*exp(8.459+0.468*log(x[59])+1.07743-0.83) +
             plogis(0.872 + 0.0000471*x[60]-1.04483+0.289)*exp(8.459+0.468*log(x[60])+1.07743-0.83) +
             plogis(0.872 + 0.0000471*x[61]-1.04483+0.289)*exp(8.459+0.468*log(x[61])+1.07743-0.83) +
             plogis(0.872 + 0.0000471*x[62]-1.04483+0.289)*exp(8.459+0.468*log(x[62])+1.07743-0.83) +
             plogis(0.872 + 0.0000471*x[63]-1.04483+0.289)*exp(8.459+0.468*log(x[63])+1.07743-0.83) +
             plogis(0.872 + 0.0000471*x[64]-1.04483+0.289)*exp(8.459+0.468*log(x[64])+1.07743-0.83) +
             plogis(0.872 + 0.0000471*x[65]-1.04483+0.289)*exp(8.459+0.468*log(x[65])+1.07743-0.83) +
             plogis(0.872 + 0.0000471*x[66]-1.04483+0.289)*exp(8.459+0.468*log(x[66])+1.07743-0.83) +
             plogis(0.872 + 0.0000471*x[67]-1.04483+0.289)*exp(8.459+0.468*log(x[67])+1.07743-0.83) +
             plogis(0.872 + 0.0000471*x[68]-1.04483+0.289)*exp(8.459+0.468*log(x[68])+1.07743-0.83) +
             plogis(0.872 + 0.0000471*x[69]-1.04483+1.081)*exp(8.459+0.468*log(x[69])+1.07743-0.637) +
             plogis(0.872 + 0.0000471*x[70]-1.04483+1.081)*exp(8.459+0.468*log(x[70])+1.07743-0.637) +
             plogis(0.872 + 0.0000471*x[71]-1.04483+1.081)*exp(8.459+0.468*log(x[71])+1.07743-0.637) +
             plogis(0.872 + 0.0000471*x[72]-1.04483+1.081)*exp(8.459+0.468*log(x[72])+1.07743-0.637) +
             plogis(0.872 + 0.0000471*x[73]-1.04483+1.081)*exp(8.459+0.468*log(x[73])+1.07743-0.637) +
             plogis(0.872 + 0.0000471*x[74]-1.04483+1.081)*exp(8.459+0.468*log(x[74])+1.07743-0.637) +
             plogis(0.872 + 0.0000471*x[75]-1.04483+1.081)*exp(8.459+0.468*log(x[75])+1.07743-0.637) +
             plogis(0.872 + 0.0000471*x[76]-1.04483+1.081)*exp(8.459+0.468*log(x[76])+1.07743-0.637) +
             plogis(0.872 + 0.0000471*x[77]-1.04483+1.081)*exp(8.459+0.468*log(x[77])+1.07743-0.637) +
             plogis(0.872 + 0.0000471*x[78]-1.04483+1.081)*exp(8.459+0.468*log(x[78])+1.07743-0.637) +
             plogis(0.872 + 0.0000471*x[79]-1.04483+1.081)*exp(8.459+0.468*log(x[79])+1.07743-0.637) +
             plogis(0.872 + 0.0000471*x[80]-1.04483+1.081)*exp(8.459+0.468*log(x[80])+1.07743-0.637) +
             plogis(0.872 + 0.0000471*x[81]-1.04483+1.081)*exp(8.459+0.468*log(x[81])+1.07743-0.637) +
             plogis(0.872 + 0.0000471*x[82]-1.04483+1.081)*exp(8.459+0.468*log(x[82])+1.07743-0.637) +
             plogis(0.872 + 0.0000471*x[83]-1.04483+1.081)*exp(8.459+0.468*log(x[83])+1.07743-0.637) +
             plogis(0.872 + 0.0000471*x[84]-1.04483+1.081)*exp(8.459+0.468*log(x[84])+1.07743-0.637) +
             plogis(0.872 + 0.0000471*x[85]-1.04483+1.081)*exp(8.459+0.468*log(x[85])+1.07743-0.637) +
             plogis(0.872 + 0.0000471*x[86]-1.04483+1.081)*exp(8.459+0.468*log(x[86])+1.07743-0.637) +
             plogis(0.872 + 0.0000471*x[87]-1.04483+1.081)*exp(8.459+0.468*log(x[87])+1.07743-0.637) +
             plogis(0.872 + 0.0000471*x[88]-1.04483+1.081)*exp(8.459+0.468*log(x[88])+1.07743-0.637) +
             plogis(0.872 + 0.0000471*x[89]-1.04483+1.081)*exp(8.459+0.468*log(x[89])+1.07743-0.637) +
             plogis(0.872 + 0.0000471*x[90]-1.04483+1.081)*exp(8.459+0.468*log(x[90])+1.07743-0.637) +
             plogis(0.872 + 0.0000471*x[91]-1.04483+1.081)*exp(8.459+0.468*log(x[91])+1.07743-0.637) +
             plogis(0.872 + 0.0000471*x[92]-1.04483+1.081)*exp(8.459+0.468*log(x[92])+1.07743-0.637) +
             plogis(0.872 + 0.0000471*x[93]-1.04483+1.081)*exp(8.459+0.468*log(x[93])+1.07743-0.637) +
             plogis(0.872 + 0.0000471*x[94]-1.04483+1.081)*exp(8.459+0.468*log(x[94])+1.07743-0.637) +
             plogis(0.872 + 0.0000471*x[95]-1.04483+1.081)*exp(8.459+0.468*log(x[95])+1.07743-0.637) +
             plogis(0.872 + 0.0000471*x[96]-1.04483+1.081)*exp(8.459+0.468*log(x[96])+1.07743-0.637) +
             plogis(0.872 + 0.0000471*x[97]-1.04483+1.081)*exp(8.459+0.468*log(x[97])+1.07743-0.637) +
             plogis(0.872 + 0.0000471*x[98]-1.04483+1.081)*exp(8.459+0.468*log(x[98])+1.07743-0.637) +
             plogis(0.872 + 0.0000471*x[99]-1.04483+1.081)*exp(8.459+0.468*log(x[99])+1.07743-0.637) +
             plogis(0.872 + 0.0000471*x[100]-1.04483+1.081)*exp(8.459+0.468*log(x[100])+1.07743-0.637) +
             plogis(0.872 + 0.0000471*x[101]-1.04483+1.081)*exp(8.459+0.468*log(x[101])+1.07743-0.637) +
             plogis(0.872 + 0.0000471*x[102]-1.04483+1.081)*exp(8.459+0.468*log(x[102])+1.07743-0.637) +
             plogis(0.872 + 0.0000471*x[103]-1.04483-1.978)*exp(8.459+0.468*log(x[103])+1.07743+0.336) +
             plogis(0.872 + 0.0000471*x[104]-1.04483-1.978)*exp(8.459+0.468*log(x[104])+1.07743+0.336) +
             plogis(0.872 + 0.0000471*x[105]-1.04483-1.978)*exp(8.459+0.468*log(x[105])+1.07743+0.336) +
             plogis(0.872 + 0.0000471*x[106]-1.04483-1.978)*exp(8.459+0.468*log(x[106])+1.07743+0.336) +
             plogis(0.872 + 0.0000471*x[107]-1.04483+0.289-1.978)*exp(8.459+0.468*log(x[107])+1.07743-0.83+0.336) +
             plogis(0.872 + 0.0000471*x[108]-1.04483+0.289-1.978)*exp(8.459+0.468*log(x[108])+1.07743-0.83+0.336) +
             plogis(0.872 + 0.0000471*x[109]-1.04483+0.289-1.978)*exp(8.459+0.468*log(x[109])+1.07743-0.83+0.336) +
             plogis(0.872 + 0.0000471*x[110]-1.04483+0.289-1.978)*exp(8.459+0.468*log(x[110])+1.07743-0.83+0.336) +
             plogis(0.872 + 0.0000471*x[111]-1.04483+1.081-1.978)*exp(8.459+0.468*log(x[111])+1.07743-0.637+0.336) +
             plogis(0.872 + 0.0000471*x[112]-1.04483+1.081-1.978)*exp(8.459+0.468*log(x[112])+1.07743-0.637+0.336) +
             plogis(0.872 + 0.0000471*x[113]-1.04483+1.081-1.978)*exp(8.459+0.468*log(x[113])+1.07743-0.637+0.336) +
             plogis(0.872 + 0.0000471*x[114]-1.04483+1.081-1.978)*exp(8.459+0.468*log(x[114])+1.07743-0.637+0.336)) )
}

# constraint function
eval_g0 <- function(x)
{
  return(x[1]+x[2]+x[3]+x[4]+x[5]+x[6]+x[7]+x[8]+x[9]+x[10]+x[11]+x[12]+x[13]+x[14]+x[15]+x[16]+x[17]+x[18]+x[19]+x[20]+x[21]+x[22]+x[23]+x[24]+x[25]+x[26]+x[27]+x[28]+x[29]+x[30]+x[31]+x[32]+x[33]+x[34]+x[35]+x[36]+x[37]+x[38]+x[39]+x[40]+x[41]+x[42]+x[43]+x[44]+x[45]+x[46]+x[47]+x[48]+x[49]+x[50]+x[51]+x[52]+x[53]+x[54]+x[55]+x[56]+x[57]+x[58]+x[59]+x[60]+x[61]+x[62]+x[63]+x[64]+x[65]+x[66]+x[67]+x[68]+x[69]+x[70]+x[71]+x[72]+x[73]+x[74]+x[75]+x[76]+x[77]+x[78]+x[79]+x[80]+x[81]+x[82]+x[83]+x[84]+x[85]+x[86]+x[87]+x[88]+x[89]+x[90]+x[91]+x[92]+x[93]+x[94]+x[95]+x[96]+x[97]+x[98]+x[99]+x[100]+x[101]+x[102]+x[103]+x[104]+x[105]+x[106]+x[107]+x[108]+x[109]+x[110]+x[111]+x[112]+x[113]+x[114]-80000)
}

# Set optimization options.
opts <- list( "algorithm"= "NLOPT_GN_ISRES",
              "xtol_rel"= 1.0e-15,
              "maxeval"= 160000,
              "print_level" = 0 )

# Solve using NLOPT_LN_COBYLA without gradient information
res1 <- nloptr( x0=c(100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100),
                eval_f=eval_f0,
                lb = c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
                ub = c(90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,90000,60000,60000,60000,60000,60000,60000,60000,60000,60000,60000,60000,60000),
                eval_g_ineq = eval_g0,
                opts = opts)
print(res1)

The results are below:

Minimization using NLopt version 2.4.2 

NLopt solver status: 5 ( NLOPT_MAXEVAL_REACHED: Optimization stopped because maxeval (above) was 
reached. )

Number of Iterations....: 160000 
Termination conditions:  xtol_rel: 0.000000000000001    maxeval: 160000 
Number of inequality constraints:  1 
Number of equality constraints:    0 
Current value of objective function:  -4536536.34788097 
Current value of controls: 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 
100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 
100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 
100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 
100 100 100 100 100 100 100 100 100 100 100 100 100 100

Could anyone take a look at my code and let me know whether I am doing anything wrong? Do I need to specify more iterations?

Edit: As @AirSquid suggested, I tried cutting the number of options down to 6 instead and that seems to have worked. So I need to understand why it doesn't work for a larger number of options

optimizationnonlinear-optimizationnlopt
Not an expert with nlopt, but a couple things to look at. Not sure what plogis is. Have you tested it? If it is returning junk or a constant, that might be the problem. Have you tried chopping your data down to something very small to see if you get believable results (maybe 4 options to start instead of 100+). Also, I think you can simplify your exponential expression and get rid of the exponent and the log w/ some algebra... - AirSquid
Retract the portion of the comment about plogis ... just looked it up - AirSquid
@AirSquid I tried cutting it down to 6 options instead and that does actually seem to have worked. Does that imply I just need to increase the iterations? The first part of my two-part model is a logistic GLM regression, and so plogis is used to obtain the probability from the log-odds that is estimated. I can look at simplifying the exponential expression - do you think that would help? - user16993967
too difficult to T/S without the whole thing, but you are on the right track. Many of your items have the same coefficient set, so it isn't clear why you set your problem up this way. x[1] - x[34] are indistinguishable in your model, so you can simplify that way too and group them. Start small, get it working, grow it. And simplify the algebra. - AirSquid
Thanks @AirSquid. The reason why they have the same coefficients is because there are multiple products which come from the same 'group' and therefore have the same coefficients. I want to allow for the option of allocating to five products from Group 1, 10 products from Group 2 etc. - user16993967