Estrés en cuerpos de agua

library(deSolve)

InitialConditions <- c(
  niños = 1e6, #personas
  jovenes = 1e6 * 0.5, #personas
  adultos = 2e6, #personas
  tercera.edad = 1e6 * 0.5, #personas
  reservas.de.agua = 100e6 #m3
  )

times <- seq(2010,  #initial time
             2060,   #end time
             1  ) #time step, year

estres.hidrico <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {
    #Endogenous auxiliary variables
    consumo.total.niños <- niños * consumo.per.capita.niños
    consumo.total.jovenes <- jovenes * consumo.per.capita.jovenes
    consumo.total.adultos <- adultos * consumo.per.capita.adultos
    consumo.total.tercera.edad <- (tercera.edad * consumo.per.capita.tercera.edad) * (1 - ahorro.total.tercera.edad) 
      
    #Flow variables
    nacimientos <- (adultos * tasa.de.natalidad)
    pubertad <- niños / tiempo.promedio.niños
    madurez <- jovenes / tiempo.promedio.jovenes
    vida.adulta <- adultos / tiempo.promedio.adultos
    muertes <- tercera.edad / tiempo.promedio.tercera.edad
    consumo <- consumo.total.niños + consumo.total.jovenes + consumo.total.adultos + consumo.total.tercera.edad
      
    #Métricas de desempeño
    deficit.hidrico <- consumo / reservas.de.agua
    estres.hidrico.adulto.activo <- (consumo.total.niños + consumo.total.jovenes + consumo.total.tercera.edad) / adultos #ojo con los parentesis
    
    #State (stock) variables (d)
    dniños <- nacimientos - pubertad
    djovenes <- pubertad - madurez
    dadultos <- madurez - vida.adulta
    dtercera.edad <- vida.adulta - muertes
    dreservas.de.agua <- recarga - evaporacion - consumo # max para que no tome valores negativos en las reservas de agua
    
    list(c(dniños, djovenes, dadultos, dtercera.edad, dreservas.de.agua),
         estres.hidrico.adulto.activo = estres.hidrico.adulto.activo,
         deficit.hidrico = deficit.hidrico)
  })
}

parameters<-c( 
  tasa.de.natalidad = 25/1000, #25 por cada 1000 adultos
  tiempo.promedio.niños = 15, #años
  tiempo.promedio.jovenes = 10, #años
  tiempo.promedio.adultos = 40, #años
  tiempo.promedio.tercera.edad = 20, #años
  recarga = 5e6, #m3
  evaporacion = 10e6, #m3 
  consumo.per.capita.niños = 100 / 1000, #litros cúbicos
  consumo.per.capita.jovenes = 100 / 1000, #litros cúbicos
  consumo.per.capita.adultos = 150 / 1000, #litros cúbicos
  consumo.per.capita.tercera.edad = 150 / 1000, #litros cúbicos # convertido para el examen
  ahorro.total.tercera.edad = 20/100 #porcentaje del consumo
  )


intg.method<-c("rk4")

out <- ode(y = InitialConditions,
           times = times,
           func = estres.hidrico,
           parms = parameters,
           method =intg.method )

plot(out,
     col=c("blue"))