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

[Thomas Petzoldt]

Hi Nicola,

the error indicates that the history buffer is exhausted. You can 
increase the buffer by setting mxhist:

y_out_2 <- dede(y = yinit, times = times_sd, func = supplydrive,
                 parms = parms_sd, method = "lsoda",
                 control = list(mxhist = 1e6))

This is found on the ` ?dede help page.

Unfortunately, this solves only part of the problem in your case, 
because then the solver seems to get stuck due to numerical problems. 
Some of the states grow exponentially to extreme values.

To see this, you can put a print or cat in the model function:

cat("time=", t, ", slag=", slag, ", llag=", llag, ",
     taulag=", taulag, "\n")

... or you can try a brute-force simulation with another solver, e.g. 
"vode" and a modified absolute tolerance (either a large value or zero). 
The relative tolerance remains unchanged.

times_sd <- seq(0, 500, by = 1)
y_out_2 <- dede(y = yinit, times = times_sd, func = supplydrive,
                 parms = parms_sd, method = "vode",
                 control = list(mxhist = 1e6), atol = 0)

plot(y_out_2)

The output contains then at least part of the solution and the plot 
shows where the it ran into trouble. Setting the time step to a larger 
value (e.g. by=1) speeds up simulation, because a too small time step 
just limits internal step size adaptation.

Thomas



On 10.03.2022 at 19:06 Nicola Gambaro wrote:
> 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
> _______________________________________________
> R-sig-dynamic-models mailing list
> R-sig-dynamic-models using r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models

-- 
Dr. Thomas Petzoldt
https://tu-dresden.de/Members/thomas.petzoldt