I am using R and the desolve package to implement a model with events occuring at both specific times and events occuring when certain condtions are true (i.e. root events).
As some background code here demonstrates using root events which trigger when a condition is met, whilst code here demonstrates using multiple non-root events.
As would be expected if a root event is used the time step will not be a trigger and vice versa if a time trigger is used the root will not trigger.
I considered using a root that returns 0 when the model time matches the event time. I have implemented this below and it appears to function however my understanding is that when building models we should not assume that the model time will occur (indeed I understand this is the entire reason behind using events). Therefore I am unsure if this is good practice or if there is a better way to do this.
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)
Extension - co-occuring events
I expected if this is the solution I may go on to have issues where multiple events occur at the same time but from experimenting it appears that providing the same paramters are not modified by multiple co-occuring events this is not an issue. I presume the correct way to handle these instances would be document the order of priority and, ideally, print/log a warning that it has occured (clearly for time triggered events a check could be made before the model is run but root triggers will only be discovered during the run).