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.