Optimise a function with two parameters and constraints in R using solnp function?

Question

I have the following function and would like to optimise parameters x and y in order to maximise the output. I am trying to do that using optimise.

EDIT

Thanks to the comments here I am able to optimise using Rsolnp::solnp. But now I get error when I add upper and lower bounds.

See example below

Ojective function:

my_fn <- function(x) { 
z <- ifelse(x[1] < 1000000, 0, 
     ifelse(x[1] < 19500000, x[1] * 0.0175,
     ifelse(x[1] < 20500000, x[1] * 0.0190,
     ifelse(x[1] < 21500000, x[1] * 0.0220,
     ifelse(x[1] < 22500000, x[1] * 0.0270,
            x[1] * 0.0340))))) +

    ifelse(x[2] < 750000, 0, 
    ifelse(x[2] < 15750000, x[2] * 0.0100,
    ifelse(x[2] < 16500000, x[2] * 0.0150,
    ifelse(x[2] < 17250000, x[2] * 0.0275,
    ifelse(x[2] < 18000000, x[2] * 0.0375,
    ifelse(x[2] < 18750000, x[2] * 0.0450,
           x[2] * 0.0550)))))) + 

return(-z) # because I want to maximise
}

Equality function:

my_fn_equal <- function(x) {
        z1 = x[1] + x[2]
        return(z1)
    }

I am using the following:

library(Rsolnp)

solnp(pars = c(1000000, 750000),   # initial values
      fun = my_fn,                 # objective function
      eqfun = my_fn_equal,         # equality function
      eqB = c(17000000)            # euqlity constraint
)

And it works fine. But once I add lower and upper bounds I start to get errors:

solnp(pars = c(1000000, 750000),        # initial values
      fun = my_fn,                       # objective function
      eqfun = my_fn_equal,               # equality function
      eqB = c(17000000),                 # euqlity constraint
      LB = c(1000000, 750000),           # lower bound for parameters 
      UB = c(16000000, 16000000)         # upper bound for parameters (I chose 16M randomly)
            )

Result with Problem Inverting Hessian:

solnp-->The linearized problem has no feasible
solnp-->solution.  The problem may not be feasible.

Iter: 1 fn: -25000.0000  Pars:  1000000.00000  750000.00000
solnp--> Solution not reliable....Problem Inverting Hessian.
$pars
[1] 1000000  750000

$convergence
[1] 2

$values
[1] -25000 -25000

$lagrange
[1] 0

$hessian
     [,1] [,2]
[1,]    1    0
[2,]    0    1

$ineqx0
NULL

$nfuneval
[1] 7

$outer.iter
[1] 1

$elapsed
Time difference of 0.1664431 secs

$vscale
[1]    25000 15250000        1        1

The constraints are:

# x[1] >= 1000000
# x[2] >= 750000
# x[1] + x[2] = 17000000

Why do I get these errors when I select bounds? What do they mean?

Also, it seems that everytime I change my initial values, the optimised values change. Why don't I get the same parameters that maximise my output everytime?

Thanks!

Answer 1

Here is how to solve the problem with as a MIP (Mixed Integer Program).

Assuming we have two vendors A and B. Each vendor provides a cash amount as bonus depending on the amount we give them (spend). The maximum spend allowed is $30 Million and we would like to maximise the total bonus. There are minimum and maximum spend constraints for each vendor.

In the real scenario there are 7 vendors wiht different constraints and bonus rules. But for simplicity I am using only 2 vendors.

Data Prep

library(tidyverse)

# run spend & bonus for A
A_value <- data.frame()
for (i in seq(0, 35000000, 100000)) {
  A_value_temp <- data.frame(vendor = 'A',
                                   spend = i,
                                   bonus = ifelse(i < 1000000, 0,
                                                  ifelse(i < 19500000, i * 0.0175,
                                                         ifelse(i < 20500000, i * 0.0190,
                                                                ifelse(i < 21500000, i * 0.0220,
                                                                       ifelse(i < 22500000, i * 0.0270,
                                                                              i * 0.0340)))))
  )

  A_value <- rbind(A_value, A_value_temp)
}

# run spend & bonus for B
B_value <- data.frame()
for (i in seq(0, 35000000, 100000)) {
  B_value_temp <- data.frame(vendor = 'B',
                               spend = i,
                               bonus = ifelse(i < 750000, 0,
                                              ifelse(i < 15750000, i * 0.0100,
                                                     ifelse(i < 16500000, i * 0.0150,
                                                            ifelse(i < 17250000, i * 0.0275,
                                                                   ifelse(i < 18000000, i * 0.0375,
                                                                          ifelse(i < 18750000, i * 0.0450,
                                                                                 i * 0.0550))))))
  )

  B_value <- rbind(B_value, B_value_temp)
}



# create control data frame
control <- rbind.data.frame(A_value, B_value)


# Build spend matrices
spend_matrix <- as.matrix(cbind(control %>%
                                      filter(vendor == "A") %>%
                                      mutate(A = spend) %>%
                                      select(A) ,

                                    control %>%
                                      filter(vendor == "B") %>%
                                      mutate(B = spend) %>%
                                      select(B)
                                    )
                          )

# Build bonus matrices
bonus_matrix <- as.matrix(cbind(control %>%
                                  filter(vendor == "A") %>%
                                  mutate(A = bonus) %>%
                                  select(A) ,

                                control %>%
                                  filter(vendor == "B") %>%
                                  mutate(B = bonus) %>%
                                  select(B)
                                )
                          )

# Check spend and bonus matrices
> tail(spend_matrix)
              A        B
[346,] 34500000 34500000
[347,] 34600000 34600000
[348,] 34700000 34700000
[349,] 34800000 34800000
[350,] 34900000 34900000
[351,] 35000000 35000000
> tail(bonus_matrix)
             A       B
[346,] 1173000 1897500
[347,] 1176400 1903000
[348,] 1179800 1908500
[349,] 1183200 1914000
[350,] 1186600 1919500
[351,] 1190000 1925000

Run MIP Optimisation

# Optimisation
library(ompr)
library(ROI)
library(ROI.plugin.glpk)
library(ompr.roi)

# Spend constraint: max spend allowed
max_spend <- 30000000

# Number of scenarios (different scenarios of spend)
scenarios <- nrow(spend_matrix)

# Number of vendors
vendors <- unique(control$vendor)


# MIP model
MIP_sol <- MIPModel() %>%
  add_variable(x[i,j], i=1:scenarios , j=1:length(vendors),  type="binary") %>%

  # Every vendor must have spend constraint
  add_constraint(sum_expr(x[i,j], i=1:scenarios) == 1, j=1:length(vendors))  %>%

  # Minimun spend constraints
  add_constraint(sum_expr(spend_matrix[i,j] * x[i,j], i=1:scenarios, j=1) >= 2500000) %>%  # Vendor A
  add_constraint(sum_expr(spend_matrix[i,j] * x[i,j], i=1:scenarios, j=2) >= 3500000) %>%  # Vendor B

  # Maximum spend constraints
  add_constraint(sum_expr(spend_matrix[i,j] * x[i,j], i=1:scenarios, j=1) <= 26000000) %>%  # Vendor A
  add_constraint(sum_expr(spend_matrix[i,j] * x[i,j], i=1:scenarios, j=2) <= 19000000) %>%  # Vendor B

  # Equality sum constraint
  add_constraint(sum_expr(spend_matrix[i,j] * x[i,j], i=1:scenarios, j=1:length(vendors)) <= max_spend) %>%   

  # Set objective: max bonus
  set_objective(sum_expr(bonus_matrix[i,j] * x[i,j], i=1:scenarios, j=1:length(vendors)), "max") %>%

  # solve model
  solve_model(with_ROI(solver = "glpk", verbose = TRUE))


> print(MIP_sol)
Status: optimal
Objective value: 1237500

Solution

# Solution
cat("Status:",solver_status(MIP_sol),"\n")
cat("Objective:",objective_value(MIP_sol),"\n")


optim_solution <- get_solution(MIP_sol, x[i, j]) %>%
  filter(value > 0) %>%
  mutate(vendor = vendors[j],
         scenario = i,
         spend = diag(spend_matrix[i,j]),
         bonus = diag(bonus_matrix[i,j])) %>%
  select(scenario, vendor, spend, bonus)

> print(optim_solution)
  scenario vendor    spend   bonus
1      111      A 11000000  192500
2      191      B 19000000 1045000