r - 您如何在 deSolve (R) 中实现系统动力学风格的流水线延迟？

Question

我正在尝试使用 R 中的 deSolve 对管道延迟进行建模。我有一个具有恒定输入 (work_arrival) 的库存 (worktodo)，我想要一个管道延迟执行 (work_rate)，其中库存以相同的到达率下降3 步延迟。目前，我能够初始化管道延迟，但它似乎在延迟后重置（开启 3 步，关闭 3 步，...）。它应该保持与 work_arrival 匹配。有任何想法吗？

####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

通过将工作到达更改为库存（状态变量），您可以延迟访问它。包 (deSolve) 似乎通过在执行计算时仅在其历史记录中保留状态变量来优化速度。

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