[R] optimizing with non-linear constraints

[Rainer M Krug]

Hi

I have posted the following question on stackoverflow [1] but haven't
received a response yet. Has somebody here any idea, how to formulate
non-linear constraints using e.g. the packages alabama or nlopr? The
question is pasted below.

I am really stuck with this.

Thanks,

Rainer

,----
| I am successful in using constrOptim() for linear constraints (thanks to
| a comment in my earlier question), but I now have a problem with
| non-linear constraints.
|
| Consider the following example:
|
| I want to fit a function FUN(X1, X2, ... X12)
|
| X1 = dep.a
| X2 = dep.b
| X3 = dep.c
| ... for all z0. , na. and zjoint.
|
| and have the following constraint
| (link to formulas:
| http://stackoverflow.com/questions/32845169/how-to-optimize-with-non-linear-constraints-in-r
| )
| |
| With some re-formating and h > 0 and LAI > 0 I get for thew first one
|
| -X1 <= LAI^X2 / X3 < 1 - X1
|
| where LAI is simply a number which constant for each fit.
|
| I saw the packages alabama and nloptr but I don't manage to get my head
| around how I can specify these types of constraints?
|
| Also a web search did not yield any helpful information for me.
|
| Could somebody provide some info how I can convert these constraints
| into arguments for the fitting functions (e.g. hin, heq in the alabama
| package)?
`----


Footnotes:
[1] http://stackoverflow.com/questions/32845169/how-to-optimize-with-non-linear-constraints-in-r

--
Rainer M. Krug, PhD (Conservation Ecology, SUN), MSc (Conservation Biology, UCT), Dipl. Phys. (Germany)

Centre of Excellence for Invasion Biology
Stellenbosch University
South Africa

Tel : +33 - (0)9 53 10 27 44
Cell: +33 - (0)6 85 62 59 98
Fax : +33 - (0)9 58 10 27 44

Fax (D): +49 - (0)3 21 21 25 22 44

email: Rai...@krugs.de

Skype: RMkrug

PGP: 0x0F52F982

[Ravi Varadhan]

Hi Rainer,
It is very simple to specify the constraints (linear or nonlinear) in "alabama" . They are specified in a function called `hin', where the constraints are written such that they are positive. Your two nonlinear constraints would be written as follows:

hin <- function(x, LAI) {
h <- rep(NA, 2)
h[1] <- LAI^x[2] / x[3] + x[1]
h[2] <- 1 - x[1] - LAI^x[2] / x[3]
h
}

Please take a look at the help page. If it is still not clear, you can contact me offline.
Best,
Ravi

Ravi Varadhan, Ph.D. (Biostatistics), Ph.D. (Environmental Engg)
Associate Professor, Department of Oncology
Division of Biostatistics & Bionformatics
Sidney Kimmel Comprehensive Cancer Center
Johns Hopkins University
550 N. Broadway, Suite 1111-E
Baltimore, MD 21205
410-502-2619


[[alternative HTML version deleted]]

______________________________________________
R-h...@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[Rainer M Krug]

Ravi Varadhan <ravi.v...@jhu.edu> writes:

> Hi Rainer,
> It is very simple to specify the constraints (linear or nonlinear) in
> "alabama" . They are specified in a function called `hin', where the
> constraints are written such that they are positive.

OK - I somehow missed the part that, when the values x are valid, i.e. in
the range as defined by the conditions, the result of hin(x) that they
are all positive.

> Your two nonlinear constraints would be written as follows:
>
> hin <- function(x, LAI) {
> h <- rep(NA, 2)
> h[1] <- LAI^x[2] / x[3] + x[1]
> h[2] <- 1 - x[1] - LAI^x[2] / x[3]
> h
> }

Makes perfect sense.

>
> Please take a look at the help page. If it is still not clear, you can contact me offline.

Yup - I did. But I somehow missed the fact stated above.

I am using constrOptim() and constrOptim.nl() for a paper
and am compiling a separate document which explains how to get the
constraints for the two functions step by step - I will make it
available as a blog post and a pdf.

I might have further questions concerning the different fitting
functions and which ones are the most appropriate in my case.

Thanks a lot,

Rainer

> Best,
> Ravi
>
> Ravi Varadhan, Ph.D. (Biostatistics), Ph.D. (Environmental Engg)
> Associate Professor, Department of Oncology
> Division of Biostatistics & Bionformatics
> Sidney Kimmel Comprehensive Cancer Center
> Johns Hopkins University
> 550 N. Broadway, Suite 1111-E
> Baltimore, MD 21205
> 410-502-2619
>
>
> [[alternative HTML version deleted]]
>

--
Rainer M. Krug
email: Rainer<at>krugs<dot>de
PGP: 0x0F52F982

[Ravi Varadhan]

I would recommend that you use auglag() rather than constrOptim.nl() in the package "alabama." It is a better algorithm, and it does not require feasible starting values.
Best,
Ravi

[Rainer M Krug]

Envoyé de mon iPhone

> Le 1 oct. 2015 à 15:17, Ravi Varadhan <ravi.v...@jhu.edu> a écrit :
>
> I would recommend that you use auglag() rather than constrOptim.nl() in the package "alabama." It is a better algorithm, and it does not require feasible starting values.

Thanks - that was one question I wanted to ask later.

I will do so,

Rainer

[Rainer M Krug]

Hi Ravi,

I would like come back to your offer. I have a problem which possibly is
caused by a bug or by something I don't understand:

My function to be minimised is executed even when an element in hin() is
negative.

My hin looks as follow:

--8<---------------cut here---------------start------------->8---
hinMahat <- function(x, hauteur, na, zjoint, y, LAI, ...) {
if (x[1] < 0) {
cat(names(list(...)), "\n")
cat(..., "\n")
cat(x, "|", hauteur, LAI, y, "\n")
}

h <- rep(NA, 8)
if (!missing(na)) {
x <- c(na, x )
}
if (!missing(y)) {
x <- c(x, y)
}
if (!missing(zjoint)) {
x <- c(x[1], zjoint, x[2])
}

##
dep <- hauteur * (0.05 + LAI^0.02 / 2) + (x[3] - 1)/20
h[1] <- dep
h[2] <- hauteur - dep
## if (h[2]==0) {
## h[2] <- -1
## }
##
z0 <- hauteur * (0.23 + LAI^0.25 / 10) + (x[3] - 1)/67
h[3] <- z0
## if (h[3]==0) {
## h[3] <- -1
## }
h[4] <- hauteur - z0
##
h[5] <- x[1]
##
h[6] <- x[2]
h[7] <- hauteur - x[2]
##
h[8] <- hauteur - dep - z0
if (any(h<=0)) {
cat(h, "\n")
cat("\n")
}
return(h)
}
--8<---------------cut here---------------end--------------->8---

the x contains up to three elements: c(na=, zjoint=, y=) and I fit these
three, unless one or two are specified explicitely.

The values going into hin are:

,----
| ... (z u ua za z0sol )
| 3 11 17 23 29 37 0.315 0.422 0.458 0.556 1.567 1.747 1.747 37 0.001
|
| x(na, zjoint): -8.875735 24.51316
| hauteur: 28
| na: 8.1
| y: 3
|
| the resulting hin() is:
| 16.09815 11.90185 11.19352 16.80648 -8.875735 24.51316 3.486843 0.708335
`----


Which is negative in element 5 as x[2]=na is negative.

So I would expect that the function fn is not evaluated. But it is, and
raises an error:

,----
| Error in wpLELMahat(z = z, ua = ua, na = ifelse(missing(na), par[1], na), :
| na has to be larger or equal than zero!
`----

Is this a misunderstanding on my part, or is it an error in the function
auglag?


Below is the function which is doing the minimisation.

If I replace auglag() with constrOptim.nl(), the optimisation is working
as expected.

So I think this is a bug in auglag?

Let me know if you need further information.

Cheers,

Rainer

--8<---------------cut here---------------start------------->8---
fitAuglag.wpLEL.mahat.single <- function(
z,
u,
LAI,
initial = c(na=9, zjoint=0.2*2, y=3),
na, zjoint, y,
h = 28,
za = 37,
z0sol = 0.001,
hin,
...
) {
if (missing(hin)) {
hin <- hinMahat
}

wpLELMin <- function(par, na, zjoint, y, z, u, ua, hauteur, za, z0sol, LAI) {
result <- NA
try({
p <- wpLELMahat(
z = z,
ua = ua,
na = ifelse(missing(na), par[1], na),
zjoint = ifelse(missing(zjoint), par[2], zjoint),
h = hauteur,
za = za,
z0sol = z0sol,
LAI = LAI,
y = ifelse(missing(y), par[3], y)
)
result <- sum( ( (p$u - u)^2 ) / length(u) )
},
silent = FALSE
)
## cat("From wpLELMin", par, "\n")
return( result )
}

ua <- u[length(u)]
result <- list()
result$method <- "fitAuglag.wpLEL.mahat.single"
result$initial <- initial
result$dot <- list(...)
result$z <- z
result$u <- u

result$fit <- auglag(
par = initial,
fn = wpLELMin,
hin = hin,
na = na,
zjoint = zjoint,
y = y,
##
z = z,
u = u,
ua = ua,
hauteur = h,
za = za,
z0sol = z0sol,
LAI = LAI,
...
)
result$wp <- wpLELMahat(
z = z,
ua = ua,
na = ifelse ( missing(na), result$fit$par["na"], na),
zjoint = ifelse ( missing(zjoint), result$fit$par["zjoint"], zjoint),
h = h,
za = za,
z0sol = z0sol,
LAI = LAI,
y = ifelse ( missing(y), result$fit$par["y"], y)
)

class(result) <- c(class(result), "wpLELFit")
return(result)
}
#+end_src--8<---------------cut here---------------end--------------->8---

[Rainer M Krug]

Please ignore - list members - accidentally CCd.

Rainer

Rainer M Krug <Rai...@krugs.de> writes:

> Hi Ravi,
>
> I would like come back to your offer. I have a problem which possibly is
> caused by a bug or by something I don't understand:
>
> My function to be minimised is executed even when an element in hin() is
> negative.
>
> My hin looks as follow:
>

[Ravi Varadhan]

Dear Rainer,
This is NOT a bug in auglag. I already mentioned that auglag() can work with infeasible starting values, which also implies that the function must be evaluable at infeasible values. A simple solution to your problem would be to fix up your objective function such that it evaluates to `Inf' or some large value, when the parameter values are not in the constrained domain. constrOptim.nl() is a barrier method so it forces the initial value and the subsequent iterates to be feasible.
Best,
Ravi
________________________________________
From: Rainer M Krug <Rai...@krugs.de>
Sent: Tuesday, October 6, 2015 9:20 AM
To: Ravi Varadhan
Cc: 'r-h...@r-project.org'
Subject: Bug in auglag?

[Rainer M Krug]

Thanks a lot Ravi for the clarification and it makes sense - bug in my
understanding acknowledged.

I'll change the objective function as suggested.

Thanks,

Rainer

Ravi Varadhan <ravi.v...@jhu.edu> writes:

> Dear Rainer,
> This is NOT a bug in auglag. I already mentioned that auglag() can
> work with infeasible starting values, which also implies that the
> function must be evaluable at infeasible values. A simple solution to
> your problem would be to fix up your objective function such that it
> evaluates to `Inf' or some large value, when the parameter values are
> not in the constrained domain. constrOptim.nl() is a barrier method
> so it forces the initial value and the subsequent iterates to be
> feasible.
> Best,
> Ravi
> ________________________________________
> From: Rainer M Krug <Rai...@krugs.de>
> Sent: Tuesday, October 6, 2015 9:20 AM
> To: Ravi Varadhan
> Cc: 'r-h...@r-project.org'
> Subject: Bug in auglag?
>
> Hi Ravi,
>
> I would like come back to your offer. I have a problem which possibly is
> caused by a bug or by something I don't understand:
>
> My function to be minimised is executed even when an element in hin() is
> negative.
>
> My hin looks as follow:
>