R_using nlsLM() with constraints

Question

I am using nlsLM {minpack.lm} to find the values of parameters a and b of function myfun which give the best fit for the data set, mydata.

mydata=data.frame(x=c(0,5,9,13,17,20),y = c(0,11,20,29,38,45))

myfun=function(a,b,r,t){
  prd=a*b*(1-exp(-b*r*t))
  return(prd)
}

and using nlsLM

myfit=nlsLM(y~myfun(a,b,r=2,t=x),data=mydata,start=list(a=2000,b=0.05),
                  lower = c(1000,0), upper = c(3000,1))

It works. But now I would like to introduce a constraint which is a*b<1000. I had a look at the option available in nlsLM to set constraint via nls.lm.control. But it's not much of help. can somebody help me here or suggest a different method to to this?

Answer 1

To my knowledge, in nlsLM of minpack.lm package lower and upper parameter bounds are the only available constrains.
One possible solution to your problem is based on the use of the package nloptr, an R interface to the free open-source library NLopt for nonlinear optimization.
In the code below I use cobyla, an algorithm for derivative-free optimization with nonlinear inequality and equality constraints:

mydata <- data.frame(x=c(0,5,9,13,17,20), y=c(0,11,20,29,38,45))

myfun=function(a,b,r,t){
  prd=a*b*(1-exp(-b*r*t))
  return(prd)
}

library(nloptr)

sse <- function(ab,r,x,y){
   sum((y - myfun(ab[1],ab[2],r,x))^2)
}

nonlincons <- function(ab,r,x,y) {
   h <- numeric(1)
   h[1] <-  1000 - ab[1]*ab[2]
   return(h)
}

x0 <- c(2000,0.05)
optpar <- cobyla(x0=x0, fn=sse, lower = c(1000,0), upper = c(3000,1), 
                 hin=nonlincons, r=2, x=mydata$x, y=mydata$y, 
                 control = list(xtol_rel = 1e-12, maxeval = 20000))
(optpar <- optpar$par)

sum((mydata$y-myfun(optpar[1],optpar[2],r=2,mydata$x))^2)

The values of the optimal parameters are:

[1] 3.000000e+03 2.288303e-02

and the sum of squared errors is:

[1] 38.02078

Hope it can help you.

Answer 2

It seems that nlsLM doesn't accept constraints, one possible method is to use constrOptim function instead, but the drawback is that you have to turn the problem into the form that constrOptim accepts. By the way, I'm currently working on an R package that gives a nice interface for optim and constrOptim and make necessary transformation automatically, but it is not ready to use yet. Here is how to use constrOptim for this problem.

The first step is that constrOptim only accepts linear constraints, a * b < 1000 => log(a) + log(b) < log(1000), so you need to reparametrize the problem using log(a) and log(b), the constraints becomes: -log(a) - log(b) > -log(1000), log(a) >= log(1000), -log(a) >= -log(3000), -log(b) >= 0, ui and ci are matrix form for these constraints. And myfun1 is the reparametrize version for your function, and then you can use constrOptim to obtain your solution for log(a) and log(b).

mydata=data.frame(x=c(0,5,9,13,17,20),y = c(0,11,20,29,38,45))

myfun=function(a,b,r,t){
    prd=a*b*(1-exp(-b*r*t))
    return(prd)
}

myfun1 <- function(logab, data){
    sum((data$y - myfun(exp(logab[1]), exp(logab[2]), r = 2, t = data$x)) ^ 2)
}

ui <- rbind(c(-1,-1), c(1,0), c(-1,0), c(0,-1))
ci <- c(-log(1000), log(1000), -log(3000), 0)
init <- log(c(2000, 0.05))
r <- constrOptim(theta = init, f = myfun1, grad = NULL, ui = ui, ci = ci, data = mydata)