rnonlinear-optimizationnlopt

How to construct an objective function with n terms for optimisation in R using nloptr?


I have updated this question to specify the problem more clearly.

I've got a dataset of the form:

#sample dataset
dt <- read.table(text = 'Plant Type Ownership Capacity VC TargetPLF
1      T Base       Pub 1000 1.50 0.80  
2      S Base       Pub 500  1.62 0.80
3      Y Base       Pub 500  2.43 0.75
4      K Base       Pub 500  1.93 0.80
5      J Base       Pub 500  3.40 0.80', header =T)

demand <- c(1500,2000,2200,3000,2800)

My objective function looks like this:

eval_f <- function(x)   #objective function
{
    return((1.50*1000*x[1] + 1.62*500*x[2] + 2.43*500*x[3] + 1.93*500*x[4] + 3.40*500*x[5]) / (1000*x[1] + 500*x[2] + 500*x[3] + 500*x[4] + 500*x[5]))
}

Other arguments:

x0 <- rep(0.45,5) #initialisation values

local_opts <- list( "algorithm" = "NLOPT_LD_MMA", "xtol_rel" = 1.0e-15 )  #options for optimisation
opts = list( "algorithm"= "NLOPT_GN_ISRES","xtol_rel"= 1.0e-15,"maxeval"= 160000,"local_opts" = local_opts,"print_level" = 0 ) #optimisation algorithm

# constraint function
eval_g0_eq <- function(x) {
return(1000*x[1] + 500*x[2] + 500*x[3] + 500*x[4] + 500*x[5] - 1500)
}

lb = c(rep(0.40,5))   #lower bounds of controls
ub = dt$Target        #upper bounds of controls

nloptr(x0,eval_f,lb=lb,ub=ub,eval_g_eq = eval_g0_eq,opts=opts)  #optimisation

How can I conduct this optimisation 5 times, with a different constraint function each time?

So in the second iteration, the constraint function should take the second element from the vector 'demand' and become this:

# constraint function for second iteration
eval_g0_eq <- function(x) {
return(1000*x[1] + 500*x[2] + 500*x[3] + 500*x[4] + 500*x[5] - 2000)
}
nloptr(x0,eval_f,lb=lb,ub=ub,eval_g_eq = eval_g0_eq,opts=opts)  #optimisation

Further, is there a way to construct eval_f and eval_g0_eq using the values in dt, and not manually as done now?


Solution

  • I am not aware of this nlopt function or library where it comes from, but I think you do not need to create objective function or constraint function repeatedly. You may use the R's vectorised calculation potential to do this.

    dt <- read.table(text = 'Plant Type Ownership Capacity VC TargetPLF
    1      T Base       Pub 1000 1.50 0.80  
    2      S Base       Pub 500  1.62 0.80
    3      Y Base       Pub 500  2.43 0.75
    4      K Base       Pub 500  1.93 0.80
    5      J Base       Pub 500  3.40 0.80', header =T)
    
    demand <- c(1500,2000,2200,3000,2800)
    
    library(nloptr)
    dt
    #>   Plant Type Ownership Capacity   VC TargetPLF
    #> 1     T Base       Pub     1000 1.50      0.80
    #> 2     S Base       Pub      500 1.62      0.80
    #> 3     Y Base       Pub      500 2.43      0.75
    #> 4     K Base       Pub      500 1.93      0.80
    #> 5     J Base       Pub      500 3.40      0.80
    

    The below code will do iterations equal to number of rows in dt

    lapply(seq(nrow(dt)), function(.x) nloptr(x0 = rep(0.45,5), 
                                          eval_f = function(x) sum(dt$VC * dt$Capacity * x)/sum(dt$Capacity * x), 
                                          lb=c(rep(0.40,5)), 
                                          ub = dt$Target, 
                                          eval_g_eq = function(x) sum(dt$Capacity * x) - demand[.x] , 
                                          opts = list( "algorithm"= "NLOPT_GN_ISRES",
                                                       "xtol_rel"= 1.0e-15,
                                                       "maxeval"= 160000,
                                                       "local_opts" = list("algorithm" = "NLOPT_LD_MMA", 
                                                                           "xtol_rel" = 1.0e-15 ),
                                                       "print_level" = 0 ) 
                                          )
           )
    

    Its output will be

    #> [[1]]
    #> 
    #> Call:
    #> nloptr(x0 = rep(0.45, 5), eval_f = function(x) sum(dt$VC * dt$Capacity * 
    #>     x)/sum(dt$Capacity * x), lb = c(rep(0.4, 5)), ub = dt$Target, 
    #>     eval_g_eq = function(x) sum(dt$Capacity * x) - demand[.x], 
    #>     opts = list(algorithm = "NLOPT_GN_ISRES", xtol_rel = 1e-15, 
    #>         maxeval = 160000, local_opts = list(algorithm = "NLOPT_LD_MMA", 
    #>             xtol_rel = 1e-15), print_level = 0))
    #> 
    #> 
    #> 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: 1e-15 maxeval: 160000 
    #> Number of inequality constraints:  0 
    #> Number of equality constraints:    1 
    #> Current value of objective function:  1.9506666666645 
    #> Current value of controls: 0.7 0.4 0.4 0.4 0.4
    #> 
    #> 
    #> 
    #> [[2]]
    #> 
    #> Call:
    #> nloptr(x0 = rep(0.45, 5), eval_f = function(x) sum(dt$VC * dt$Capacity * 
    #>     x)/sum(dt$Capacity * x), lb = c(rep(0.4, 5)), ub = dt$Target, 
    #>     eval_g_eq = function(x) sum(dt$Capacity * x) - demand[.x], 
    #>     opts = list(algorithm = "NLOPT_GN_ISRES", xtol_rel = 1e-15, 
    #>         maxeval = 160000, local_opts = list(algorithm = "NLOPT_LD_MMA", 
    #>             xtol_rel = 1e-15), print_level = 0))
    #> 
    #> 
    #> 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: 1e-15 maxeval: 160000 
    #> Number of inequality constraints:  0 
    #> Number of equality constraints:    1 
    #> Current value of objective function:  1.893 
    #> Current value of controls: 0.8 0.8 0.4 0.8 0.4
    #> 
    #> 
    #> 
    #> [[3]]
    #> 
    #> Call:
    #> nloptr(x0 = rep(0.45, 5), eval_f = function(x) sum(dt$VC * dt$Capacity * 
    #>     x)/sum(dt$Capacity * x), lb = c(rep(0.4, 5)), ub = dt$Target, 
    #>     eval_g_eq = function(x) sum(dt$Capacity * x) - demand[.x], 
    #>     opts = list(algorithm = "NLOPT_GN_ISRES", xtol_rel = 1e-15, 
    #>         maxeval = 160000, local_opts = list(algorithm = "NLOPT_LD_MMA", 
    #>             xtol_rel = 1e-15), print_level = 0))
    #> 
    #> 
    #> 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: 1e-15 maxeval: 160000 
    #> Number of inequality constraints:  0 
    #> Number of equality constraints:    1 
    #> Current value of objective function:  1.92181376479112 
    #> Current value of controls: 0.7883007 0.7815353 0.583928 0.7845948 0.4043761
    #> 
    #> 
    #> 
    #> [[4]]
    #> 
    #> Call:
    #> nloptr(x0 = rep(0.45, 5), eval_f = function(x) sum(dt$VC * dt$Capacity * 
    #>     x)/sum(dt$Capacity * x), lb = c(rep(0.4, 5)), ub = dt$Target, 
    #>     eval_g_eq = function(x) sum(dt$Capacity * x) - demand[.x], 
    #>     opts = list(algorithm = "NLOPT_GN_ISRES", xtol_rel = 1e-15, 
    #>         maxeval = 160000, local_opts = list(algorithm = "NLOPT_LD_MMA", 
    #>             xtol_rel = 1e-15), print_level = 0))
    #> 
    #> 
    #> 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: 1e-15 maxeval: 160000 
    #> Number of inequality constraints:  0 
    #> Number of equality constraints:    1 
    #> Current value of objective function:  1.93672658072513 
    #> Current value of controls: 0.7870349 0.7534322 0.6660887 0.770566 0.4118776
    #> 
    #> 
    #> 
    #> [[5]]
    #> 
    #> Call:
    #> nloptr(x0 = rep(0.45, 5), eval_f = function(x) sum(dt$VC * dt$Capacity * 
    #>     x)/sum(dt$Capacity * x), lb = c(rep(0.4, 5)), ub = dt$Target, 
    #>     eval_g_eq = function(x) sum(dt$Capacity * x) - demand[.x], 
    #>     opts = list(algorithm = "NLOPT_GN_ISRES", xtol_rel = 1e-15, 
    #>         maxeval = 160000, local_opts = list(algorithm = "NLOPT_LD_MMA", 
    #>             xtol_rel = 1e-15), print_level = 0))
    #> 
    #> 
    #> 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: 1e-15 maxeval: 160000 
    #> Number of inequality constraints:  0 
    #> Number of equality constraints:    1 
    #> Current value of objective function:  1.93152444911471 
    #> Current value of controls: 0.795248 0.7687759 0.6712354 0.7692782 0.4034175