NLP

In this notebook we try to solve a nonlinear optimization model posted on OR Stack Exchange. We deal only with the case where the \(x\) variables are continuous rather than binary. The problem and solution process are discussed in a post on my blog.

The solution approach combines one dimensional line search with solving a linear program. We use OMPR to model the LP and CPLEX (via the ROI plugin) to solve the LPs.

First, we load the required packages.

library(ompr)
library(ompr.roi)
library(ROI)
library(ROI.plugin.cplex)
library(magrittr)

Preparation

Data

The problem author provided an example. We set up the example data here.

p <- c(99.9, 200.2, 300.0)
f <- c(100, 120, 320, 80, 45)
b <- t(matrix(c(10, 100, 20, 39, 4,
                70, 10, 200, 49, 40,
                60, 50, 300, 50, 60), nrow = 3, ncol = 5, byrow = TRUE))
beta <- 1  # placeholder value

Next, we define problem dimensions.

tmax <- 4 # $T$ in the original post
n <- 3

Lastly, we need upper bounds for the nonnegative variables \(q_1\) and \(q_2\).

q1max <- 100
q2max <- 100

Functions

Unless noted otherwise, all functions will take as argument a list containing scalar elements q1 and q2 and vector element x.

We need to define functions for the right-hand sides of the three constraints and their common left-hand side.

lhs <- function(v) drop(p %*% v$x) # LHS of (1)-(3)

rhs1 <- function(v) sum(f / (1 + v$q1)^(0:tmax))  # RHS of (1)

rhs2 <- function(v) sum(drop(b %*% v$x) / (1 + v$q2)^(0:tmax)) # RHS of (2)

Before proceeding, we need to find a value for \(\beta\) that makes the model feasible. To that end, we will pick arbitrary values for \(x\), use equation (1) to find \(q_1\), use equation (2) to find \(q_2\), and then choose \(\beta\) to balance equation (3).

temp <- list(x = c(0.2, 0.1, 0.3)) # arbitrary
junk <- function(q) rhs1(list(q1 = q)) - lhs(temp)
temp$q1 <- uniroot(junk, lower = 0, upper = q1max, tol = 1e-9)$root
junk <- function(q) {
  v <- temp
  v$q2 <- q
  rhs2(v) - lhs(temp)
}
temp$q2 <- uniroot(junk, lower = 0, upper = q2max, tol = 1e-9)$root
beta <- lhs(temp) / rhs1(list(q1 = temp$q1 + temp$q2))
cat("Known feasible solution: q1 = ", temp$q1, ", q2 = ", temp$q2, ", x = ", temp$x, ".\nUsing beta = ", beta, ".\n")

Known feasible solution: q1 =  4.893879 , q2 =  0.3784701 , x =  0.2 0.1 0.3 .
Using beta =  1.01866 .

rm(temp, junk)

Armed with a reasonable value for \(\beta\), we can define the function that computes the right-hand side of equation (3).

rhs3 <- function(v) beta * rhs1(list(q1 = v$q1 + v$q2, q2 = NA, x = NA)) # RHS of (3)

We also need a function to compute the coefficients of \(x\) in the right-hand side of constraint (2), given \(q_2\).

coefs2 <- function(v) drop((1/(1 + v$q2)^(0:tmax)) %*% b)

Given values of \(q_1\) and \(q_2\), we find a suitable vector \(x\in [0,1]^n\) by solving a linear program, using the following function. The return value is a list of all variable values \(q_1\), \(q_2\), \(x\) or NA if no feasible solution for \(x\) can be found.

findX <- function(v) {
  a <- coefs2(v)
  zz <- file("/tmp/junk", open="wt")
  # Divert annoying CPLEX messages to a junk file.
  sink(zz, type = "message")
  sink(zz)
  # Set up and solve the LP.
  sol <- MIPModel() %>% 
           add_variable(x[i], i=1:n, type = "continuous", lb = 0, ub = 1) %>% 
           set_objective(0) %>% 
           add_constraint(sum_expr(p[i] * x[i], i=1:n) == rhs1(v)) %>% 
           add_constraint(sum_expr((p[i] - a[i]) * x[i], i=1:n) == 0) %>% 
           solve_model(with_ROI(solver = "cplex"))
  # Stop diversion and close the temp file.
  sink()
  sink(type = "message")
  close(zz)
  if (sol$status == "infeasible") {
    return(NA)
  } else {
    v$x <- sol %>% get_solution(x[i])
    return(v)
  }
}

Given a value for \(q_1\), we need a function that finds the corresponding value of \(q_2\) using the equality of the right-hand sides of constraints 1 and 3. The return value is a list containing all three of \(q_1\), \(q_2\) and \(x\). If no value of \(q_2\) can be found, the function returns NA.

findQ2 <- function(v) {
  # Define a function of one variable (q2) that computes the difference between
  # the right-hand sides of constraints (1) and (3), given q1.
  r1 = rhs1(v)
  f <- function(q) { 
    v0 <- v
    v0$q2 <- q
    rhs3(v0) - r1
  }
  # Solve for a root of f. If a root is found, add it to the input list v.
  # If not, return NA.
  result <- tryCatch(uniroot(f, c(0, q2max), tol = 1e-5),
                     error = function(e) NA)
  if (is.list(result)) {
    v$q2 <- result$root
  } else {
    return(NA)
  }
  # We now have values for q1 and q2. Solve for x. If a solution is found, add it
  # to v and return the completed list. If not, return NA.
  result <- findX(v)
  if (is.list(result)) {
    # Add the objective value.
    result$obj <- result$q1 + result$q2
    return(result)
  } else {
    return(NA)
  }
}

Solution process

Throughout the process, we will maintain an incumbent solution.

incumbent <- list(obj = -1)  # dummy starting value

For convenience, we will define a function to check for a new incumbent.

checkIncumbent <- function(v) {
  if (v$obj > incumbent$obj) {
    incumbent <<- v
  }
}

Lower limit for \(q_1\)

The first step in the solution process is to establish a lower bound for \(q_1\), using constraint (1). The maximum value of the left-hand side occurs when \(x_i = 1\) for all \(i\). We solve for the corresponding value of \(q_1\) and set that as the domain minimum (q1min).

temp <- function(q) rhs1(list(q1 = q)) - sum(p)
lo <- tryCatch(uniroot(temp, c(0, q1max), tol = 1e-7), error = function(e) NA)
if (is.list(lo)) {
  q1min <- lo$root
} else {
  cat("Unable to find a valid lower bound for q_1.\n")
}
# Clean up.
rm(lo)

The lower limit may not be feasible. We can try a line search with fixed step size until we find a feasible value for \(q_1\).

stepSize <- 0.1   # step size for search
v <- list(q1 = q1min)
temp <- findQ2(v)
while (!is.list(temp)) {
  v$q1 <- v$q1 + stepSize
  if (v$q1 > q1max) {
    break
  } else {
    temp <- findQ2(v)
  }
}
if (is.list(temp)) {
  # We have a feasible solution. Check if it is a new incumbent.
  v <- temp
  checkIncumbent(v)
  # Make it the new lower limit on q1.
  q1min <- v$q1
} else {
  # We failed to find a feasible solution, so give up.
  cat("Cannot find a feasible solution!\n")
}
# Clean up.
rm(v, temp, stepSize)

Before continuing, we need to set a stopping criterion. We will stop when the gap (the width of the interval of uncertainty for the value of \(q_1\)) drops below a specified limit.

gaptol <- 1e-5

Upper limit for \(q_1\)

If the upper limit we arbitrarily set for \(q_1\) is not feasible, we can perform a bisection search to find the largest feasible value for \(q_1\).

v <- findQ2(list(q1 = q1max))
if (!is.list(v)) {
  # The upper limit is not feasible. We need to do a bisection search.
  uncertainty <- c(q1min, q1max)
  while (uncertainty[2] - uncertainty[1] > gaptol) {
    # Try the midpoint.
    mid <- (uncertainty[1] + uncertainty[2]) / 2
    v <- findQ2(list(q1 = mid))
    if (is.list(v)) {
      # The midpoint is feasible and becomes the lower end of the uncertainty interval.
      uncertainty[1] <- mid
      # Check if we stumbled on a new incumbent.
      checkIncumbent(v)
    } else {
      # The midpoint is infeasible and becomes the upper end of the uncertainty interval.
      uncertainty[2] <- mid
    }
  }
  # Set the lower limit of the interval as the upper limit for q1.
  q1max <- uncertainty[1]
  # Clean up
  rm(uncertainty, mid, v)
}
cat("The interval of uncertainty for q1 is [", q1min, ", ", q1max, "].\n")

The interval of uncertainty for q1 is [ 4.760744 ,  5.450373 ].

cat("The current incumbent solution has q1 = ", incumbent$q1, " and objective value ", incumbent$obj, ".\n")

The current incumbent solution has q1 =  5.450373  and objective value  5.916804 .

Golden section search

To find the optimum, we perform golden section search on the interval from q1min to q1max. GS search inserts two symmetrically positioned points in the interval, giving us a total of four points including the endpoints. Since we are maximizing, if the maximum occurs at either end, we drop the opposite endpoint (leaving three points) and insert a new point in the symmetric position to the remaining interior point. If the maximum occurs at an interior point, we drop the endpoint that is further away from the maximum and again add a symmetric point.

# Set up the four points in the initial interval and get their objective values.
gs <- 0.618 # approximate golden section factor
gs1 <- 1 - gs
interval <- c(q1min, gs * q1min + gs1 * q1max, gs * q1max + gs1 * q1min, q1max)
obj <- c(findQ2(list(q1=q1min))$obj, 0, 0, findQ2(list(q1=q1max))$obj)
temp <- findQ2(list(q1=interval[2]))
obj[2] <- temp$obj
checkIncumbent(temp)
temp <- findQ2(list(q1=interval[3]))
obj[3] <- temp$obj
checkIncumbent(temp)
# While the interval width exceeds the gap limit, apply GS search.
while(interval[4] - interval[1] > gaptol) {
  # Find the index of the point with greatest objective value.
  i <- which.max(obj)
  # Replace the opposite endpoint.
  if (i >= 3) {
    # Replace the lower endpoint.
    interval[1] <- interval[2]
    obj[1] <- obj[2]
    # The higher interior point becomes the lower interior point.
    interval[2] <- interval[3]
    obj[2] <- obj[3]
    # Insert a new upper interior point.
    interval[3] <- gs1 * interval[1] + gs * interval[4]
    temp <- findQ2(list(q1 = interval[3]))
    checkIncumbent(temp)
    obj[3] <- temp$obj
  } else {
    # Replace the upper endpoint.
    interval[4] <- interval[3]
    obj[4] <- obj[3]
    # The lower interior point becomes the higher interior point.
    interval[3] <- interval[2]
    obj[3] <- obj[2]
    # Insert a new lower interior point.
    interval[2] <- gs * interval[1] + gs1 * interval[4]
    temp <- findQ2(list(q1 = interval[2]))
    checkIncumbent(temp)
    obj[2] <- temp$obj
  }
}
# Clean up.
rm(gs, gs1, temp)

Display the incumbent solution.

cat("Incumbent solution:\nq1 = ", incumbent$q1, "\nq2 = ", incumbent$q2, "\nx = ", incumbent$x$value, "\nObjective value = ", incumbent$obj, ".\n")

Incumbent solution:
q1 =  5.450373 
q2 =  0.4664311 
x =  1 0.1334608 -8.212684e-07 
Objective value =  5.916804 .

# Clean up x and see how badly constraints are violated.
v <- incumbent
v$x <- incumbent$x$value %>% pmin(., 1) %>% pmax(., 0)
cat("Adjusted x vector = ", v$x, "\n")

Adjusted x vector =  1 0.1334608 0 

cat("Common LHS of all three constraints = ", lhs(v), "\n")

Common LHS of all three constraints =  126.6189 

cat("RHS of (1) = ", rhs1(v), "\n")

RHS of (1) =  126.6186 

cat("RHS of (2) = ", rhs2(v), "\n")

RHS of (2) =  126.6188 

cat("RHS of (3) = ", rhs3(v), "\n")

RHS of (3) =  126.6186