Как реализовать конвейерную задержку в стиле системной динамики в deSolve (R)?

Question

Я пытаюсь смоделировать задержку конвейера с помощью deSolve в R. У меня есть один запас (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()