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)) {
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 .