Adding model-fitting constraints to a subset of parameters using FME::modFit or optim

I am attempting to fit a differential equation model by minimizing least squares, where three of the four parameters to be optimized over can't fall below zero. Other than that restriction, I would prefer not apply any further constraints.

Here is the data, model and helper functions:

library(deSolve)
library(FME)

# Data
y.data <- data.frame(P = c(57537.4049311277, 58610.5565091595, 59380.7326528789, 59831.1412677501, 59956.9401326381, 59770.9438648549, 59339.7636243282, 58743.8966682831, 58062.6867570216, 57372.6802888231, 56729.6957034719, 56118.5748350006, 55507.8890166606, 54867.5694355373, 54168.9851576162, 53385.1758149806, 52500.082193378, 51534.2686505716, 50516.106686327), 
SSL = c(5918.49121715316, 6223.93819671156, 6522.34271426757, 6817.47760344158, 7115.49425740418, 7423.50644959141, 7748.62985898112, 8097.51802623853, 8472.32658437055, 8872.9327092582, 9294.62369710465, 9730.87922007949, 10175.6848086032, 10623.7661461289, 11075.7535333951, 11538.2076285642, 12018.4500914707, 12518.681172246, 13039.7328357312),
S = c(61.8032786885246, 69.8469945355191, 71.1475409836066, 75.0163934426229, 65.6393442622951, 63.7049180327869, 60.0437158469945, 67.9344262295082, 66.5573770491803, 67.0327868852459, 57.6010928961749, 43.7158469945355, 57.224043715847, 65.1366120218579, 49.2131147540984, 30.7158469945355, 34.3715846994535, 12.3661202185792, 27.4972677595628))

# Initial state conditions
y0 <- c(P=y.data[1,]$P, SSL=y.data[1,]$SSL, S=y.data[1,]$S)

# Initial free parameter values
a <- 0.00075
b <- 0.004464286
rs <- 0.01
ks <- 100
p0 <- c(a=a,b=b,rs=rs,ks=ks)

# Series along which to integrate
t <- 0:18
l <- length(t)

# Differential equation model
LVfn_dc <- function(time, y, param) {
  P = y[1]
  SSL = y[2]
  S = y[3]
  with(as.list(param), {
    dPdt <- ((-0.4101*4)*time^3) + ((18.058*3)*time^2) - ((297.9*2)*time) + 1511.2
    dSSLdt <- ((7.2468*2)*time) + 265.67      
    dSdt <- (rs*S)*((ks-S-(a*P)-(b*SSL))/ks)
    return(list(c(dPdt,dSSLdt,dSdt)))
  })
}

# Function to solve differential equation
sde <- function(param, time, f, y) {
  sol <- deSolve::ode(y = y, 
                      times = time, 
                      func = f, 
                      parms = param,
                      method = "rk4",
                      rtol = 1e-15, atol = 1e-15, maxsteps = 1e9)
  return(sol)
}

# Least squares function for optimization using FME::modFit
LS_modFit <- function(param, time, falala, y, y.values){
  # call sde function to compute y.theta(t) after initializing
  y.theta.all <- c()
  y.theta <- c()
  call.sde <- sde(param,time,falala,y)
  y.theta.all <- call.sde[,-1] # disregard t column
  # only keep y.theta(t.i) values for each data point t.i
  for(i in 1:length(time)){
    if(time[i]%%1 == 0){
      y.theta <- rbind(y.theta,y.theta.all[i,])
    }
  }
  error <- sum((y.values - y.theta)^2)
  return(error)
}

I have tried using both FME::modFit and stats::optim, only bounding the lower limits of the constrained parameters, but this produces "non-finite finite-difference value" error messages.

# Optimizing with limited lower bounds
fit <- FME::modFit(f = LS_modFit, p = p0,
                     time=t, falala=LVfn_dc, y=y0, y.values=y.data,
                     upper=c(Inf,Inf,Inf,Inf), lower=c(0,0,-Inf,1),
                     method = "L-BFGS-B")

When I use non-infinite, but generous upper and lower bounds for all four parameters (that don't produce infinite values when run with the sde function), I still receive non-finite type errors.

# Optimizing with upper and lower bounds for all parameters
fit <- FME::modFit(f = LS_modFit, p = p0,
                     time=t, falala=LVfn_dc, y=y0, y.values=y.data,
                     upper=c(0.02,0.02,1,400), lower=c(0,0,-2,65),
                     method = "L-BFGS-B")

Any insight as to where I'm going wrong?

Solution

Not all optimizers support constraints natively, so FME uses a transformation for others. I see two different ways to get it work:

Use finite constraints, e.g. 100 or 1e6
use an optimizer that supports constraints natively, e.g. Portor bobyqa.

I modified also the following:

Instead of the loop, I used the pre-defined function modCost and added timeas an independent variable to y.data. Note that it is important to define a weighting method if the state variables have different scales
replaced rk4 by lsoda.
I recommend also to reduce rtoland rtol as it makes the simulation slow and does not help much, because the optimizers have their own tolerances.

library(deSolve)
library(FME)

# Data
y.data <- data.frame(
  time = 0:18,
  P = c(57537.4049311277, 58610.5565091595, 59380.7326528789, 59831.1412677501, 59956.9401326381, 59770.9438648549, 59339.7636243282, 58743.8966682831, 58062.6867570216, 57372.6802888231, 56729.6957034719, 56118.5748350006, 55507.8890166606, 54867.5694355373, 54168.9851576162, 53385.1758149806, 52500.082193378, 51534.2686505716, 50516.106686327),
  SSL = c(5918.49121715316, 6223.93819671156, 6522.34271426757, 6817.47760344158, 7115.49425740418, 7423.50644959141, 7748.62985898112, 8097.51802623853, 8472.32658437055, 8872.9327092582, 9294.62369710465, 9730.87922007949, 10175.6848086032, 10623.7661461289, 11075.7535333951, 11538.2076285642, 12018.4500914707, 12518.681172246, 13039.7328357312),
  S = c(61.8032786885246, 69.8469945355191, 71.1475409836066, 75.0163934426229, 65.6393442622951, 63.7049180327869, 60.0437158469945, 67.9344262295082, 66.5573770491803, 67.0327868852459, 57.6010928961749, 43.7158469945355, 57.224043715847, 65.1366120218579, 49.2131147540984, 30.7158469945355, 34.3715846994535, 12.3661202185792, 27.4972677595628)
)

# Initial state conditions
y0 <- c(P=y.data[1,]$P, SSL=y.data[1,]$SSL, S=y.data[1,]$S)

# Initial free parameter values
a <- 0.00075
b <- 0.004464286
rs <- 0.01
ks <- 100
p0 <- c(a=a,b=b,rs=rs,ks=ks)

# Series along which to integrate
t <- 0:18
l <- length(t)

# Differential equation model
LVfn_dc <- function(time, y, param) {
  P = y[1]
  SSL = y[2]
  S = y[3]
  with(as.list(param), {
    dPdt <- ((-0.4101*4)*time^3) + ((18.058*3)*time^2) - ((297.9*2)*time) + 1511.2
    dSSLdt <- ((7.2468*2)*time) + 265.67
    dSdt <- (rs*S)*((ks-S-(a*P)-(b*SSL))/ks)
    return(list(c(dPdt,dSSLdt,dSdt)))
  })
}

# Function to solve differential equation
sde <- function(param, time, f, y) {
  sol <- deSolve::ode(y = y,
                      times = time,
                      func = f,
                      parms = param,
                      method = "lsoda", # lsoda is more stable and faster than rk4
                      rtol = 1e-6, atol = 1e-6)
  return(sol)
}

# Least squares function for optimization using FME::modFit
LS_modFit <- function(param, time, falala, y, y.values){
  call.sde <- sde(param, time, falala, y)
  ## important! choose appropriate weighting method,
  ## depending on scale of variables
  ## options:  one of "none", "std", "mean"
  modCost(call.sde, y.values, weight="std")
}


# Optimizing with limited lower bounds
fit <- FME::modFit(f = LS_modFit, p = p0,
                     time=t, falala=LVfn_dc, y=y0, y.values=y.data,
                     ## infinite "constraints" work only
                     ## for a subset of solvers, e.g. Port and bobyqa
                     upper=c(Inf,Inf,Inf,Inf), lower=c(0,0,-Inf,1),
                     ## workaround for the other optimizers
                     ## recommended: don't make limits too big
                     #upper=c(1e6,1e6,1e6,1e6), lower=c(0,0,-1e6,1),
                     method = "Port")

guess <- sde(p0, t, LVfn_dc, y0)
fitted <- sde(fit$par, t, LVfn_dc, y0)

plot(guess, fitted, obs=y.data, mfrow=c(1, 3))
legend("bottomleft", col=1:2, lty=1:2, legend=c("guess", "fitted"))

solution of the ODE system with fitted parameters

^{Created on 2023-04-01 with reprex v2.0.2}