[R-sig-dyn-mod] Problem with dede's lagvalue (deSolve)

[Nicola Gambaro]

Dear R dynamics community

I am trying to solve a dynamic model based on DDEs with deSolve’s dede solver. However, I am encountering problems understanding and solving errors related to the lagvalue function. Here’s the code:

----------------------------------------------------
library(deSolve)
supplydrive <- function(t, y, parms) {
  with(as.list(c(y, parms)),{
    #y[1] = demand
    #y[2] = price
    #y[3] = supply
    #y[4] = extraction rate
    #y[5] = production stock
    #y[6] = manufacturing stock
    #y[7] = in-use stock
    #y[8] = waste management stock
    
    #define delays
    #Extraction rate // s
    slag <- t - s
    if (slag <= 0) {
      e_slag <- y[4]
    } else {
      e_slag <- lagvalue(slag, 4)
    }
    #Stocks // l
    llag <- t - l
    if (llag <= 0) {
      s3_llag <- y[7]
      s2_llag <- y[6]
    } else {
      s3_llag <- lagvalue(llag, 7)
      s2_llag <- lagvalue(llag, 6)
    }
    #Price // tau
    taulag <- t - tau
    if (taulag <= 0) {
      pricelag <- y[2]
      dpricelag = 0
    } else {
      pricelag <- lagvalue(taulag, 2)
      dpricelag <- lagderiv(taulag, 2)
    }
    
    #waste and dissipation flows
    w1 = eta*y[4] + gamma*(1-eta)*e_slag#delay s<---
    w4 = theta*(alpha4*s3_llag + alpha7*s2_llag)#delay l
    d1 = lambda*y[5]
    d2 = lambda*y[6]
    d3 = lambda*y[7]
    d4 = lambda*y[8]
    totald = d1+d2+d3+d4
    
    #supply, demand, price
    dy1 = r*y[1]
    dy2 = c*y[2]*((y[1]-y[3])/y[1])
    dy3 = y[4] - (w1 + w4 + totald)
    
    #extraction rate
    if (dpricelag > 0)
      dy4 = y[4] * (dpricelag/pricelag) #delay tau
    else
      dy4 = z*dy1
    
    #stocks
    dy5 = y[4] - alpha1*y[5] + alpha2*y[6] + alpha5*y[8] - w1 - d1
    dy6 = alpha1*y[5] - alpha2*y[6] - alpha3*y[6] + alpha6*y[8] - alpha7*y[6] - d2
    dy7 = alpha3*y[6] - alpha4*y[7] - d3
    dy8 = alpha4*y[7] - alpha5*y[8] - alpha6*y[8] + alpha7*y[6] - w4 - d4
    
    list(c(dy1, dy2, dy3, dy4, dy5, dy6, dy7, dy8))
  })
}

yinit <- c(1, 1, 1, 1, 1, 1, 1, 1)
times_sd <- seq(0, 200, by = 0.1)
parms_sd <- c(alpha = c(0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1), 
                 tau = 120, r = 0.1, c = 0.8, s=2, l=1, 
                 eta=0.2, gamma=0.3, theta=0.5, 
                 lambda=0.01, z=0.2)
y_out_2 <- dede(y = yinit, times = times_sd, func = supplydrive, parms = parms_sd, method = "lsoda”)

—————————————————

This runs fine. However, if the times sequence length is increased to, say,

times_sd <- seq(0, 500, by = 0.1) ,

then I get the following error

Error in lagvalue(taulag, 2) : 
illegal input in lagvalue - lag, 355.98, too large, at time = 475.98

However, this shouldn’t be the case because the conditional statements in the function should make sure that the lagtime is not to a time before the beginning of the simulation (t=0). Is this a memory issue? How could I go about solving this issue?

Many thanks,
Nicola