CODEWITHSUNDEEP
rconstraintsnls
rm167
rm167

Reputation: 1247

R_using nlsLM() with constraints

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?

Upvotes: 5

Views: 3742

Answers (2)

Consistency
Consistency

Reputation: 2922

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)

Upvotes: 2

Marco Sandri
Marco Sandri

Reputation: 24252

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.

Upvotes: 7

Related Questions