Добавление ограничений для подгонки модели к подмножеству параметров с помощью FME::modFit или optim

Вы можете рассмотреть возможность использования пакета DEoptim R. Я попробовал, и у меня не было ошибки с вашим кодом:

library(DEoptim)
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)
}


obj_DEoptim <- DEoptim(fn = LS_modFit, lower = c(0,0,-2,65), upper = c(0.02,0.02,1,400), 
                       time = t, falala = LVfn_dc, y = y0, y.values = y.data)

Не все оптимизаторы изначально поддерживают ограничения, поэтому FME использует преобразование для других. Я вижу два разных способа заставить его работать:

1. Используйте конечные ограничения, например. 100 или 1e6
2. используйте оптимизатор, который изначально поддерживает ограничения, например. Port или bobyqa.

Я также изменил следующее:

• Вместо цикла я использовал предопределенную функцию modCost и добавил timeв качестве независимой переменной к y.data. Обратите внимание, что важно определить метод взвешивания, если переменные состояния имеют разные масштабы.
• заменено rk4 на lsoda.
• Я рекомендую также уменьшить rtol и rtol, так как это делает симуляцию медленной и мало помогает, потому что у оптимизаторов есть свои допуски.

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"))

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