Basado en el Problema 6.2 del libro Small System Dynamics Models for Big Issues: Triple Jump towards Real-World Complexity de Erik Pruyt.

Supongamos que existe una plaga de cucarachas en una zona determinada de la ciudad. Para iniciar nuestra modelación consideraremos que hay 20,000 cucarachas. Existe un aumento autónomo en el número de cucarachas que asciende a 5 cucarachas por cucaracha por semana debido a su alta tasa de reproducción. Sin embargo, hay una tasa de decrecimiento autónomo relacionado con la muerte natural o por otro medio no relacionado a fumigación de las cucarachas, la cual es de 1 cucaracha por cucaracha por semana.

Para controlar esta plaga es necesario realizar fumigaciones periódicas con el fin de reducir el número de cucarachas. Cada año se conceden 200 licencias para realizar fumigaciones en la zona afectada, es decir 200 personas y/o empresas tienen la responsabilidad de fumigar. Las licencias son válidas solo durante un año, y cada titular de una licencia debe realizar 5 fumigaciones al año. El número de cucarachas eliminadas por fumigación depende del número de cucarachas que aún no han muerto ni han sido eliminadas por fumigaciones anteriores, así como de la eficiencia de las fumigaciones, la cual consideraremos que oscila entre el 15% y el 25%.

  1. En Vensim, dibuje el diagrama causal del problema (identifica los ciclos de retroalimentación existentes)

  1. En Vensim, dibuje el diagrama de flujo del problema

  1. Explica cuál será el comportamiento esperado de la variable de estado basándote en el comportamiento que ves en ambos diagramas

  2. Modela para un período de 5 años y determina cuál debería ser la eficiencia de las fumigaciones deseada para evitar el crecimiento exponencial de las cucarachas

library("deSolve")
library("ggplot2")

# Variables exógenas o parámetros
parametros <- c(tasa_decrec = 1, # Cucarachas que mueren mensualmente
                tasa_reproduccion = 5, # Cucarachas que nace mensualmente
                minima_fumigacion = 0.15, # Tasa de efectividad mínima de fumigación
                ef_fumigacion_max = 0.25, # Tasa de efectividad máxima de fumigación
                fum_licencia = 5, # Número de fumigaciones por año
                licencias = 200, # Número de licencias por año
                eficiencia_ideal = 0.20
                ) 

# Condiciones iniciales de las variables de estado (aquí se escriben sin "d")
InitialConditions <- c(cucarachas = 20000)

# Tiempo de modelación
times <- seq(0, # tiempo inicial
             12*5, # fin de la simulación
             1) # intervalo de tiempo

# Método de integración
intg.method <- c("rk4")

# Función del caso (se llena de abajo hacia arriba)
plagas <- function(t, state, parametros) {
  with(as.list(c(state, parametros)), {
    
    # Variables endógenas o auxiliares
    total_fumigaciones = (licencias*fum_licencia/52)
    exito_fumigacion = cucarachas * eficiencia_ideal
    
    # Variables de flujo
    nacimientos_cucarachas = cucarachas * tasa_reproduccion
    natural_cucarachas = cucarachas * tasa_decrec
    cucarachas_fumigacion = exito_fumigacion * total_fumigaciones
    
    # Variables de estado (siempre inician con un "d")
    dcucarachas = nacimientos_cucarachas - natural_cucarachas - cucarachas_fumigacion 
     
    list(c(dcucarachas)) # En la lista deben ir las variables que quieres graficar, dentro de la c() las variables de estado (las que inician con "d") y fuera de la c(), pero dentro del paréntesis de la lista, las demás variables que se deseen plotear.
  })
}

# Simulación del modelo
out <- ode(y = InitialConditions,
           times = times,
           func = plagas, # recuerda actualizar el nombre de tu función para cada caso que resuelvas
           parms = parametros,
           method = intg.method)

# Gráficos de resultados del modelo
plot(out, main = "Población de cucarachas", ylab = "Cucarachas", xlab = "Semanas")

df = as.data.frame(out)
ggplot(df, mapping = aes (x = times, y = cucarachas)) +
         geom_line (color = "blue") +
         geom_point (color = "green") +
         labs(title = "Población de cucarachas", x = "Semanas", y = "Cucarachas") +
         theme_minimal()