How do you implement a system dynamics style pipeline delay in deSolve (R)?

Question

I'm attempting to model a pipeline delay with deSolve in R. I have one stock (worktodo) that has a constant input (work_arrival) and I want a pipeline delay execution (work_rate) where the stock goes down by the same arrival rate with a 3 step delay. Currently, I am able to initialize the pipeline delay, but it seems to reset after the lag (on for 3 steps, off for 3 steps, ...). It should stay on to match the work_arrival. Any ideas?

####System Dyanmics Model - Pipeline Delay
library(deSolve)
library(tidyverse)

#model setup
finaltime  =  50
initialtime  =  0
timestep  =  1

#create a time vector
simtime <- seq(initialtime, finaltime, by= timestep)

#add auxs
auxs <- c(
   work_arrival = 50
)

#add stocks
stocks <- c(
   worktodo= 600 )



# This is the model function
model <- function(time, stocks, auxs){
  with(as.list(c(stocks, auxs)),{
#add aux calculations

   tlag <- 3
   if(time < tlag){
      work_rate = 0
   }
   else{
      ylag <- lagderiv(time - tlag)
      work_rate <- ylag
   }

   #if(time == 3) print(structure(ylag))


#add stock calculations

   worktodo  =  work_arrival - work_rate

#return data
return(list(c(

   worktodo),
   work_rate = work_rate,
   work_arrival = work_arrival))
  })
}

data <- data.frame(dede(y= stocks, times = simtime, func = model, parms = auxs, method = "lsodar"))

df <- data %>% 
   pivot_longer(-time, names_to = 'variable')


ggplot(df, aes(time, value, color = variable))+
   geom_line(size =1.25)+
   theme_minimal()

Answer 1

By changing the work arrival to a stock (state variable) you are able to access it as a lag. The package (deSolve) seems to optimize speed by keeping only state variables in its history as it performs calculations.

####System Dyanmics Model - Pipeline Delay
library(deSolve)
library(tidyverse)

#model setup
finaltime  =  50
initialtime  =  0
timestep  =  1

#create a time vector
simtime <- seq(initialtime, finaltime, by= timestep)

#add auxs
auxs <- c(
  work_arrival = 50
)

#add stocks
stocks <- c(
  worktodo= 600 ,
  work_arrival_stock = 50
  )



# This is the model function
model <- function(time, stocks, auxs){
  with(as.list(c(stocks, auxs)),{
    #add aux calculations
    #work_arrival_stock_depletion = work_arrival_stock
    tlag <- 3
    if(time < tlag){
      work_rate = 0
    }
    else{
      ylag <- lagvalue(time - tlag)[2] #[2] grabs the value of the second stock
      work_rate <- ylag
    }

    #if(time == 3) print(structure(ylag))


    #add stock calculations
    worktodo  =  work_arrival - work_rate
    work_arrival_stock = 0


    #return data
    return(list(c(
      worktodo,
      work_arrival_stock),
      work_rate = work_rate,
      work_arrival = work_arrival))
  })
}

data <- data.frame(dede(y= stocks, times = simtime, func = model, parms = auxs, method = "lsodar"))

df <- data %>% 
  pivot_longer(-time, names_to = 'variable')


ggplot(df, aes(time, value, color = variable))+
  geom_line(size =1.25)+
  theme_minimal()