您如何在deSolve（R）中实现系统动力学样式的管线延迟？

Question

我正在尝试使用R中的deSolve对流水线延迟进行建模。我有一个库存（worktodo）的输入为常数（work_arrival），并且我希望流水线延迟执行（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()