1

我正在使用 R+desolve 对海藻生长进行建模，并希望实现一个函数，在给定模型内部变量的阈值的情况下收获海藻（即不是状态变量，因此不能通过 (t , y, 参数))。我还希望阈值由我可以在 parms 中传递给模型的变量控制。一些（简化的）示例代码

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,

        )
  })
}

然后我想在g_E达到阈值时触发一个事件函数

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

和事件

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

收获参数被添加到传递给 ode 调用的参数中

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)

desolve 事件的所有示例都对模型的状态变量（例如，在上面的示例中可以访问 y['N_s']）而不是模型的内部变量进行根评估。上面 rootfunc 中写的 g_E 显然不起作用，但我不知道如果可能的话如何访问它。用模型的状态变量替换 g_E（在科学上没有用，但对测试有好处）我得到错误：找不到对象'threshold_value'，所以看起来根函数甚至无法访问参数。

我确定我在这里错过了一些基本的东西。非常感谢您的帮助！

4

0 回答 0