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).
In deSolve: Can't understand how to early stop ode solver with root functions I sketched the general idea of the root->action event mechanism. This means that desolve actively searches for a root if a sign change in one component of the root function is detected. As such the result of the root function should be a list of a continuous functions that change sign at the root. For the fixed times it should be a list of t-event_time(k).
To avoid code duplication, and thus for higher flexibility, in eventfunc simply use yroot=rootfunc(t,y,parms).