在 R desolve 中结合根事件函数和非根事件函数
Combining root and non-root event functions in R desolve
我正在使用 R 和 desolve 包来实现一个模型,该模型具有在特定时间发生的事件和在某些条件为真时发生的事件(即根事件)。
一些后台代码here demonstrates using root events which trigger when a condition is met, whilst code here演示了如何使用多个非根事件。
正如预期的那样,如果使用根事件,时间步将不会触发,反之亦然,如果使用时间触发器,根将不会触发。
我考虑过在模型时间与事件时间匹配时使用 returns 0 的根。我已经在下面实现了它并且它似乎起作用但是我的理解是在构建模型时我们不应该假设模型时间会发生(事实上我理解这是使用事件背后的全部原因)。因此,我不确定这是否是好的做法,或者是否有更好的方法。
library("desolve")
yini <- c(temp = 18, heating_on = 1, eventoccurs = 0)
temp <- function(t,y, parms) {
dy1 <- ifelse(y[2] == 1, 1.0, -0.5)
dy2 <- 0
dy3 <- 0
list(c(dy1, dy2, dy3))
}
event_times <- c(1, 5, 20)
rootfunc <- function(t, y, parms) {
time_in_event <- 1
if(t %in% event_times){ time_in_event = 0}
yroot <- c(y[1] - 18, y[1] - 20, time_in_event)
return(yroot)
}
eventfunc <- function(t, y, parms) {
time_in_event <- 1
if(t %in% event_times){ time_in_event = 0}
yroot <- c(y[1] - 18, y[1] - 20, time_in_event)
whichroot <- which(abs(yroot) < 1e-6) # specify tolerance
if(whichroot == 2) { y[2] <- 0 } else { y[2] <- 1 }
if(whichroot == 3) { y[3] <- 1 } else { y[3] <- 0 }
return(y)
}
times <- seq(from=0, to=20,by=0.1)
out <- lsode(times=times, y=yini, func = temp, parms = NULL,
rootfun = rootfunc, events = list(func=eventfunc, root = TRUE))
plot(out, lwd=2)
扩展 - 共同发生的事件
我预计如果这是解决方案,我可能会继续遇到多个事件同时发生的问题,但从实验来看,提供相同的参数似乎不会被多个同时发生的事件修改,这不是问题。我认为处理这些实例的正确方法是记录优先级顺序,理想情况下,print/log 警告它已经发生(显然对于时间触发事件,可以在模型 运行 但根触发器只会在 运行).
期间被发现
在中我勾勒出了root->action事件机制的大致思路。这意味着如果检测到根函数的一个分量发生符号变化,desolve
会主动搜索根。因此,根函数的结果应该是一个在根处改变符号的连续函数的列表。对于固定时间,它应该是 t-event_time(k)
.
的列表
为避免代码重复,从而获得更高的灵活性,在 eventfunc
中只需使用 yroot=rootfunc(t,y,parms)
.
我正在使用 R 和 desolve 包来实现一个模型,该模型具有在特定时间发生的事件和在某些条件为真时发生的事件(即根事件)。
一些后台代码here demonstrates using root events which trigger when a condition is met, whilst code here演示了如何使用多个非根事件。
正如预期的那样,如果使用根事件,时间步将不会触发,反之亦然,如果使用时间触发器,根将不会触发。
我考虑过在模型时间与事件时间匹配时使用 returns 0 的根。我已经在下面实现了它并且它似乎起作用但是我的理解是在构建模型时我们不应该假设模型时间会发生(事实上我理解这是使用事件背后的全部原因)。因此,我不确定这是否是好的做法,或者是否有更好的方法。
library("desolve")
yini <- c(temp = 18, heating_on = 1, eventoccurs = 0)
temp <- function(t,y, parms) {
dy1 <- ifelse(y[2] == 1, 1.0, -0.5)
dy2 <- 0
dy3 <- 0
list(c(dy1, dy2, dy3))
}
event_times <- c(1, 5, 20)
rootfunc <- function(t, y, parms) {
time_in_event <- 1
if(t %in% event_times){ time_in_event = 0}
yroot <- c(y[1] - 18, y[1] - 20, time_in_event)
return(yroot)
}
eventfunc <- function(t, y, parms) {
time_in_event <- 1
if(t %in% event_times){ time_in_event = 0}
yroot <- c(y[1] - 18, y[1] - 20, time_in_event)
whichroot <- which(abs(yroot) < 1e-6) # specify tolerance
if(whichroot == 2) { y[2] <- 0 } else { y[2] <- 1 }
if(whichroot == 3) { y[3] <- 1 } else { y[3] <- 0 }
return(y)
}
times <- seq(from=0, to=20,by=0.1)
out <- lsode(times=times, y=yini, func = temp, parms = NULL,
rootfun = rootfunc, events = list(func=eventfunc, root = TRUE))
plot(out, lwd=2)
扩展 - 共同发生的事件
我预计如果这是解决方案,我可能会继续遇到多个事件同时发生的问题,但从实验来看,提供相同的参数似乎不会被多个同时发生的事件修改,这不是问题。我认为处理这些实例的正确方法是记录优先级顺序,理想情况下,print/log 警告它已经发生(显然对于时间触发事件,可以在模型 运行 但根触发器只会在 运行).
在desolve
会主动搜索根。因此,根函数的结果应该是一个在根处改变符号的连续函数的列表。对于固定时间,它应该是 t-event_time(k)
.
为避免代码重复,从而获得更高的灵活性,在 eventfunc
中只需使用 yroot=rootfunc(t,y,parms)
.