Search code examples
rmathematical-optimization

constrained nonlinear minimization with many variables


Here is a minimization problem I've meant to solve, but no matter what form or package I try it with, it never resolves itself.

The Problem is a transportation problem with a quadratic objective function. It is formulated as follows:

Minimize f(x), with f(x) being x' * C * x, subject to the equality constraints UI * x - ci = 0.

where C is a diagonal matrix of constants, UI is matrix with the values 0, 1, -1 in order to set up the constraints.

I'll provide an example that I have tried with two functions so far, nloptr from its likewise called package and constrOptim.

Here's an example for nloptr:

require(nloptr)
objective <- function(x) {return( list( "objective" = t(x) %*% C %*% x,
                                    "gradient" =  2* C %*% x )) }
constraints <- function(x) {return( list( "constraints" =   ui %*% x - ci,
                                        "jacobian" = ui))}

C <- diag(c(10,15,14,5,6,10,8))
ci <- c(20, -30, -10, -20, 40))
ui <- rbind( c(1,1,1,0,0,0,0),
             c(-1,0,0,1,0,0,0),
             c(0,-1,0,-1,1,1,0),
             c(0,0,-1,0,-1,0,1),
             c(0,0,0,0,0,-1,-1))
opts <- list("alorithm" = "NLOPT_GN_ISRES")

res <- nloptr( x0=x0, eval_f=objective, eval_g_eq = constraints, opts=opts)

When trying to solve this Problem with constrOptim, I face the problem that I have to provide starting values that are within the feasible region. However, I will ultimately have a lot of equations and don't really know how to set these starting points.

Here's the same example with constrOptim:

C <- diag(c(10,15,14,5,6,10,8))
ci <- c(20, -30, -10, -20, 40)
ui <- rbind(    c(1,1,1,0,0,0,0),
            c(-1,0,0,1,0,0,0),
            c(0,-1,0,-1,1,1,0),
            c(0,0,-1,0,-1,0,1),
            c(0,0,0,0,0,-1,-1))
start <- c(10,10,10,0,0,0,0)

objective <- function(x) { t(x) %*% C %*% x }
gradient <- function(x) { 2 * C %*% x }
constrOptim(start, objective, gradient, ui = ui, ci = ci)

Solution

  • Try this:

    co <- coef(lm.fit(ui, ci))
    co[is.na(co)] <- 0
    res <- nloptr( x0=co, eval_f=objective, eval_g_eq = constraints, 
      opts=list(algorithm = "NLOPT_LD_SLSQP"))
    

    giving:

    > res
    
    Call:
    nloptr(x0 = co, eval_f = objective, eval_g_eq = constraints, 
        opts = list(algorithm = "NLOPT_LD_SLSQP"))
    
    
    Minimization using NLopt version 2.4.0 
    
    NLopt solver status: 4 ( NLOPT_XTOL_REACHED: Optimization stopped because 
    xtol_rel or xtol_abs (above) was reached. )
    
    Number of Iterations....: 22 
    Termination conditions:  relative x-tolerance = 1e-04 (DEFAULT) 
    Number of inequality constraints:  0 
    Number of equality constraints:    5 
    Optimal value of objective function:  37378.6963822218 
    Optimal value of controls: 28.62408 -29.80155 21.17747 -1.375917 -17.54977 -23.6277 -16.3723