[R-sig-dyn-mod] Accessing internal time-varying and fixed parameters in root function and event function in deSolve

[Dr.Martin Johnson]

Hi all,

I'm using R+desolve to model seaweed growth and wish to implement a function in which seaweed is harvested given a threshold value of an internal variable of the model (i.e. not a state variable so not accessed via the 'y' of (t, y, parms)). I would also like the threshold value to be controlled by a variable that I can pass to the model in parms.  Some (simplified) example code

Seaweed_model <- function(t, y, parms,...) {

  with(as.list(c(y, parms)), {
    #(lots of internal variables hidden here for brevity)

   I_top<-PAR(t)*exp(-K_d*(z-h_MA))                                                          # calculate incident irradiance at the top of the farm
    I_av <-(I_top/(K_d*z+N_f*a_cs))*(1-exp(-(K_d*z+N_f*a_cs)))          # calclulate the average irradiance throughout the height of the farm
    g_E <- I_av/((I_s)+I_av))                                                                          # light limitation scaling function


#(state variables not relevant to this problem omitted)

    dN_f        <- mu_g_EQT*N_s-(d_m*N_f)                                    # change in fixed nitrogen

    list(c(dNH4, dNO3,dN_s,dN_f,dD),
         g_E       = g_E,
        )
  })
}

Then I would like to trigger an event function when g_E reaches a threshold value

rootfunc  <- function(t,y,parms,...){
  return(g_E - threshold_value)  #how do I access these values from the ode solver call?
}

And the event

eventfunc <- function(t, y, parms,...){
  c(y[1:2],y[3]*harvest_fraction,y[4]*harvest_fraction,y[5]) #reduces biomass state variables fractionally
}

The harvest parameters are added in to the parms passed to the ode call

parms_harvest<-c(
  harvest_threshold=0.2,
  harvest_fraction=0.75
)

lim_harvest <- ode(times = times, func = Seaweed_model,  y = y0, parms = c(parms_porphyra,parms_farm,parms_harvest), events=list(func=eventfunc, root=TRUE),rootfun=rootfunc)

Am I trying to achieve the impossible? All of the examples for desolve events do a root evaluation on a state variable of the model (e.g. y['N_s'] is accessible in the example above) rather than internal variables of the model. g_E as written in the above rootfunc clearly will not work but I have no idea how to access it if it's even possible.  Substituting  a state variable of the model for g_E (not useful scientifically but good for testing purposes) I get error: object 'threshold_value' not found so it doesn't seem like the root function is able to access parms even.

I'm sure I've missed something fundamental here. Thanks so much in advance for you help!

Martin Johnson




	[[alternative HTML version deleted]]