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