Hi, I've been able to track the problem down.
Apparently I was misled by the documentation of the R package, which states
on page 5 that the argument 'hin' to the auglag function
defines the inequalty constraints, hin(x) >= 0.
However, this is not true. It should say
defines the inequality constraints, hin(x) <= 0.
As a result my inequality constraints in combination with the bound
constraints resulted in a feasible region that was the empty set and the
whole thing blew up.
Attached is a self-contained example in R, that can be run with
"Rscript libnlopt-debug.r".
The example solves the problem
min f(x; a, b): R^2 -> R
where a and b are two constant vectors. At the optimum, x1=mean(a) and
x2=mean(b) with both x1, x2 >= 0. Bound constraints are set to [0, 10].
It can be seen that if you define
hin(x) = (x1, x2)
the optimization fails with generic failure after 399 iterations with x1 and
x2 approaching zero, whereas if you define
hin(x) = (-x1, -x2)
the algorithm converges to the correct solution (x1=2.56, x2=2.64) in 29
iterations.
Regards.
Ernest
library(nloptr)
set.seed(123)
A=rpois(25, 2)
B=rpois(25, 3)
x0 = c(1, 1)
f <- function(x) {
-sum(A * log(x[1]) - x[1] + B * log(x[2]) - x[2])
}
df <- function(x) {
c(-sum(A / x[1] - 1), -sum(B / x[2] - 1))
}
hin <- function(x) {
-x # x gives generic failure
}
nloptr(x0,
eval_f=f,
eval_grad_f=df,
eval_g_ineq=hin,
eval_jac_g_ineq=function(x) nl.jacobian(x, hin),
lb=c(1e-15, 1e-15), ub=c(10, 10),
opts=list(
print_level=3,
algorithm='NLOPT_LD_AUGLAG',
xtol_rel=1e-10,
xtol_abs=1e-10,
maxeval=-1,
maxtime=60,
local_opts=list(
algorithm='NLOPT_LD_SLSQP',
xtol_rel=1e-10)))
cat('analytical solution:', mean(A), mean(B), '\n')
_______________________________________________
NLopt-discuss mailing list
[email protected]
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/nlopt-discuss
