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

Suppose there is a muskrat plague in a particular area. At first, there were 100 muskrats. The autonomous increase in the number of muskrats per muskrat per year amounts to an average of 20 muskrats per muskrat per year. Suppose that each year, 10 licences are granted to set muskrat traps. The licences are only valid for one year and each person holding a licence may set 10 traps. Assume the number of muskrats caught per trap is proportional to the number of muskrats and a catch rate per trap which is close to 0.2, say, on average 0.2, minimally 0.195, and maximally 0.205.

  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

La variable de estado (la población de muskrats) puede crecer sin limite, estabilizarse o desaparecer dependiendo de la tasa de captura. Existe un punto en donde la población se mantiene pero con la minima variación el sistema cambia y por lo tanto cambia la población de muskrats.

  1. Modela para un período de 10 años y determina cuál debería ser la tasa de captura por trampa deseada para evitar el crecimiento exponencial de las ratas almizcleras.
library("deSolve")
# Variables exógenas o parámetros
parameters <- c(autonomous.increase = 0.2, # muskrats/muskrat/year
                licences.granted = 10, # licences/year
                traps.per.licence = 10, # traps/licence
                catch.rate.per.trap = 0.2) # muskrats/trap/year
# Condiciones iniciales de las variables de estado (aquí se escriben sin "d")
InitialConditions <- c(muskrats = 100)
# Tiempo de modelación
times <- seq(0, # tiempo inicial
             10, # fin de la simulación
             0.1) # intervalo de tiempo
# Método de integración
intg.method <- c("rk4")
# Función del caso
muskrat.plague <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    # Variables endógenas
    traps.set <- licences.granted * traps.per.licence # traps/year
    catch.rate <- catch.rate.per.trap * muskrats # muskrats/trap/year 
    # Variables de flujo
    increase.rate <- autonomous.increase * muskrats # muskrats/year
    catch.rate.total <- catch.rate * traps.set # muskrats/year
    # Variables de estado (siempre inician con un "d")
    dmuskrats <- increase.rate - catch.rate.total
    list(c(dmuskrats)) # 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
muskrat.simulation <- ode(y = InitialConditions, 
                           times = times, 
                           func = muskrat.plague, 
                           parms = parameters, 
                           method = intg.method)
# Gráfica
plot(muskrat.simulation, col = c("blue"))