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

Do I see it correctly that your model returns the new state and not a derivative? And, is it a discrete time model and not continuous? Then solver `lsodar` (and all the others supported by `dede`) are not the right ones. The `iteration` solver would do, but it does not support delays. The good news is, that such a solver for discrete time steps should not be too difficult to implement in R directly. — tpetzoldt, Mar 05 '20 at 23:06
I completely agree that iteration would be the ideal situation. However, I'm trying to keep the method and package because it is easy to manipulate and well known for those people converting SD models over to R. I think I found a way to keep make this happen, but I had to make my rate a stock. The package (deSolve) and dede only seem to keep stocks in their history memory for lag access. This seems to greatly improve the speed. I benchmarked deSolve against several standard iterators in R and can't beat it. My speed also slows dramatically when I try to incorporate lags. — JD Caddell, Apr 05 '20 at 02:21

score 2 · Accepted Answer · answered Apr 05 '20 at 02:26

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

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

1 Answers1