Diagramas Causales, Estructura y Comportamiento de Sistemas Dinámicos

Instrucciones y recomendaciones

  1. Cada estudiante debe entregar de manera individual su trabajo y es responsable de manera individual por la calidad de los productos entregados.
  2. La colaboración grupal es deseable durante la resolución de esta tarea, pero cada estudiante debe demostrar capacidad crítica individual en la resolución de los ejercicios descritos en este documento.
  3. La presentación de tu trabajo es muy importante. Sé conciso y claro en tu exposición, describe adecuadamente tus gráficos y empléalos para hacer más claros y contundentes tus argumentos. No escribas argumentos o ideas incompletas, concluye tus pensamientos de manera adecuada.
  4. Si empleas tablas y/o gráficos, asegúrate de incluir notas para cada uno de ellos. Cada elemento gráfico en tu exposición debe explicarse por sí mismo. Toma como referencia las anotaciones que Sterman hace en figuras y tablas en el libro de texto.
  5. Revisa cuidadosamente cada inciso, responde a la totalidad de las preguntas.
  6. No te limites a la información presentada en cada caso. Si consideras que los casos son incompletos busca de manera individual más información y referencia adecuadamente tus fuentes.
  7. Es altamente recomendable consultar el libro de texto para resolver los ejercicios descritos en esta tarea.

Problema 1: Preguntas de Opción Múltiple (20 puntos)

Revisa el libro de Pruyt (2013) y responde a las preguntas de opción múltiple listadas a continuación. Al entregar tu actividad sólo indica la opción seleccionada.

Capítulo 3:

Multiple Choice Question 2, página 69 (1 punto), D

Multiple Choice Question 3, página 70 (1 punto), C

Multiple Choice Question 4, página 70 (1 punto), D

Multiple Choice Question 5, página 71 (1 punto), C

Multiple Choice Question 6, página 71 (1 punto), B

Multiple Choice Question 8, página 72 (1 punto), A

Multiple Choice Question 9, página 72 (1 punto), D

Multiple Choice Question 12, página 74 (1 punto), A

Multiple Choice Question 13, página 74 (1 punto), A

Multiple Choice Question 14, página 75 (1 punto), A

Multiple Choice Question 17, página 76 (1 punto), B

Multiple Choice Question 19, página 78 (1 punto), D

Capítulo 7:

Multiple Choice Question 1, página 104 (1 punto), D

Multiple Choice Question 2, página 104 (1 punto), C

Multiple Choice Question 3, página 104 (1 punto), C

Multiple Choice Question 6, página 106 (1 punto), D

Multiple Choice Question 8, página 107 (1 punto), A

Multiple Choice Question 11, página 108 (1 punto), C

Multiple Choice Question 12, página 109 (1 punto), A

Multiple Choice Question 16, página 112 (1 punto), C

Problema 2: Pneumonic Plague (20 puntos)

Introducción:

On 3 August 2009, media reported about an outbreak of pneumonic plague in north-west China:

“A second man has died of pneumonic plague in a remote part of north-west China where a town of more than 10,000 people has been sealed off. […] Local officials in north-western China have told the BBC that the situation is under control, and that schools and offices are open as usual. But to prevent the plague [from] spreading, the authorities have sealed off Ziketan, which has some 10,000 residents. About 10 other people inside the town have so far contracted the disease, according to state media. No-one is being allowed [to] leave the area, and the authorities are trying to track down people who had contact with the men who died. […] According to the WHO, pneumonic plague is the most virulent and least common form of plague. It is caused by the same bacteria that occur in bubonic plague – the Black Death that killed an estimated 25 million people in Europe during the Middle Ages. But while bubonic plague is usually transmitted by flea bites and can be treated with antibiotics, [pneumonic plague, which attacks the lungs, can spread from person to person or from animals to people], is easier to contract and if untreated, has a very high case-fatality ratio.” (BBC 03/08/2009)

Descripción del caso:

You are asked to make a SD model of this outbreak. Use the following assumptions: The total population of Ziketan amounted initially to 10000 citizens. New infections make that citizens belonging to the susceptible population become part of the infected population, which initially consists of just 1 person. The number of infections equals the product of the infection ratio, the contact rate, the susceptible population, and the infected fraction. Initially, the normal contact rate amounts to 50 contacts per week and the infection ratio to a staggering 75% per contact. The infected fraction equals of course the infected population over the sum of all other subpopulations. If citizens from the infected population die, they enter the statistics of the deceased population, else they are quarantined to recover. The recovering could be modeled simplistically as (1 − fatality ratio) ∗ infected population / recovery time. Suppose for the sake of simplicity that the average recovery time and the average decease time are both 2 days. The fatality ratio depends on the antibiotics coverage of the population which –in this poor part of China– is 0% at first.

The fatality ratio is 90% at 0% antibiotics coverage of the population and 15% at 100% antibiotics coverage of the population. Assume for the sake of simplicity that those belonging to the recovering population do not pose any threat of infection, either because they are really quarantined or because they are not contagious any more.

Preguntas del caso:

2.1. Desarrolla el diagrama de flujo del caso (5 puntos, producto: diagrama de flujo)

knitr::include_graphics("C:\\Users\\Javier Cazares\\OneDrive\\Escritorio\\EGAP\\Materias\\Modelación de Sistemas\\Neumonic Plague.png")

2.2. Construye el modelo de dinámica de sistemas basado en R y simula el modelo por un periodo de 1 mes.(5 puntos, productos: hipótesis dinámica y gráfico que la respalde)

# Cargamos el paquete deSolve
library(deSolve)

# Definimos la función que describe las ecuaciones diferenciales del modelo
pneumonic_plague <- function(time, state, parameters) {
  with(as.list(c(state, parameters)), {
    # Calculamos la fracción infectada
    infected_fraction <- I / (S + I + R)
    
    # Calculamos los flujos
    infections <- contact_rate * infection_ratio * S * infected_fraction
    recoveries <- (1 - fatality_ratio) * I / recovery_time
    deaths <- fatality_ratio * I / death_time
    
    # Calculamos las tasas de cambio
    dS <- -infections
    dI <- infections - recoveries - deaths
    dR <- recoveries
    dD <- deaths
    
    # Retornamos las tasas de cambio
    return(list(c(dS = dS, dI = dI, dR = dR, dD = dD)))
  })
}

# Definimos los parámetros del modelo
parameters <- c(contact_rate = 50,  # contactos por semana
                infection_ratio = 0.75,  # probabilidad de infección por contacto
                recovery_time = 2,  # días
                death_time = 2,  # días
                fatality_ratio = 0.9)  # porcentaje

# Definimos las condiciones iniciales utilizando los nombres de las variables utilizados en las ecuaciones
initial_conditions <- c(S = 9999, I = 1, R = 0, D = 0)

# Definimos el tiempo de simulación
times <- seq(0, 30, by = 1)  # simulamos 30 días (1 mes)

# Simulamos el modelo
out <- ode(y = initial_conditions,
           times = times,
           func = pneumonic_plague,
           parms = parameters)

# Graficamos los resultados y añadimos la leyenda directamente
matplot(times, out[, -1], type = "l", col = 1:4, lty = 1, xlab = "Días", ylab = "Número de Personas", main = "Evolución del Brote de Peste Neumónica")
legend("topright", legend = c("Susceptibles", "Infectados", "Recuperados", "Fallecidos"), col = 1:4, lty = 1)

Hipótesis Dinámica:

2.3. Grafica y comenta la evolución de infections, deaths, recovering population y deceased population (5 puntos, producto: gráficos de evolución y comentarios del comportamiento)

# Simulamos el modelo
out <- ode(y = initial_conditions,
           times = times,
           func = pneumonic_plague,
           parms = parameters)

# Extraemos los datos de las columnas para cada grupo poblacional
time <- out[, "time"]
susceptible <- out[, "S"]
infected <- out[, "I"]
recovered <- out[, "R"]
deceased <- out[, "D"]

# Graficamos los resultados
plot(time, susceptible, type = "l", col = "blue", ylim = c(0, max(susceptible)), xlab = "Días", ylab = "Número de Personas", main = "Evolución del Brote de Peste Neumónica")
lines(time, infected, col = "red")
lines(time, recovered, col = "green")
lines(time, deceased, col = "black")

# Añadimos una leyenda para clarificar las líneas
legend("topright", legend = c("Susceptibles", "Infectados", "Recuperados", "Fallecidos"), col = c("blue", "red", "green", "black"), lty = 1)

Comentarios

Infections: Se espera un rápido aumento al inicio debido a la alta tasa de contacto y de infección, seguido por una disminución conforme la población susceptible se reduce y más infectados se recuperan o fallecen.

Deaths: Las muertes aumentarán conforme al pico de infecciones,reflejando a la alta tasa de mortalidad. Luego, disminuirán a medida que el número de casos activos disminuta.

Recovering Population: Aumentará después del pico de infecciones, reflejando el tiempo de recuperación de 2 días y la proporción de infectados que sobreviven.

Decesead Population: Seguirá una trayectoria parecida a las muertes, acumulándose con el tiempo hasta que el brote sea controlado.

2.4. Suponga ahora que la antibiotics coverage of the population es del 100%. Compara la evolución de de infections, deaths, recovering population y deceased population con las obtenidas en la pregunta anterior. (5 puntos, productos: código en R, gráficas y comentarios comparando ambos modelos)

# Definimos la función modificada con 100% de cobertura de antibióticos
pneumonic_plague_100_antibiotics <- function(time, state, parameters) {
  with(as.list(c(state, parameters)), {
    # Calculamos la fracción infectada
    infected_fraction <- I / (S + I + R + D)
    
    # Calculamos los flujos
    infections <- contact_rate * infection_ratio * S * infected_fraction
    recoveries <- (1 - fatality_ratio_100) * I / recovery_time
    deaths <- fatality_ratio_100 * I / death_time
    
    # Tasas de cambio
    dS <- -infections
    dI <- infections - recoveries - deaths
    dR <- recoveries
    dD <- deaths
    
    # Retornamos las tasas de cambio y los flujos
    return(list(c(dS, dI, dR, dD)))
  })
}

# Parámetros modificados para 100% cobertura de antibióticos
parameters_100 <- parameters
parameters_100["fatality_ratio_100"] <- 0.15  # Fatality ratio al 15%

# Simulamos el modelo con 100% cobertura
out_100 <- ode(y = initial_conditions,
               times = times,
               func = pneumonic_plague_100_antibiotics,
               parms = parameters_100)
# Extraer los datos del modelo
time <- out[, "time"]
susceptible <- out[, "S"]
infected_0 <- out[, "I"]
recovered_0 <- out[, "R"]
deceased_0 <- out[, "D"]

infected_100 <- out_100[, "I"]
recovered_100 <- out_100[, "R"]
deceased_100 <- out_100[, "D"]

# Graficar los resultados de ambos escenarios
matplot(time, cbind(susceptible, infected_0, infected_100, recovered_0, recovered_100, deceased_0, deceased_100), 
        type = "l", lty = c(1, 1, 2, 1, 2, 1, 2), col = c("blue", "red", "red", "green", "green", "black", "black"),
        xlab = "Días", ylab = "Número de Personas", main = "Comparación de Escenarios de Cobertura de Antibióticos")
legend("topright", legend = c("Susceptibles", "Infectados 0%", "Infectados 100%", "Recuperados 0%", "Recuperados 100%", "Fallecidos 0%", "Fallecidos 100%"), 
       col = c("blue", "red", "red", "green", "green", "black", "black"), lty = c(1, 1, 2, 1, 2, 1, 2))

Infections: Se espera que el número de infecciones inicialmente sea similar en ambos casos pero, dado que hay mayor cobertura de antibióticos los infectados podrán tomar menos tiempo en recuperarse pero más en morir, lo que podría compensar para tener un comportamiento similar a lo largo del tiempo.

Deaths:La cantidad de muertes se espera que sea considerablemente menor asumiendo un escenario con 100% de cobertura de antibióticos debido a la reducción drástica de la tasa de mortalidad.

Recovering Popultion: Una menor tasa de mortalidad (escenario de 100% de cobertura de antibióticos) conlleva a su vez un aumento en el número de recuperaciones.

Decesead Population: La población fallecida será mucho menor en el escenario de 100% de cobertura de anibióticos.

Problema 3: Administración de Cadena de Suministro (30 puntos)

Introducción:

After following this SD course, you decide to start working as a model-based business analyst. Your first job consists in helping out the manager of a luxury car company that produces hand-made sports cars. The company faces severe inventory, production and workforce fluctuations following sales fluctuations. The manager wants you to make a model to understand the interactions between production, inventory and workforce in order to be able to reduce these fluctuations.

Descripción del caso:

Modelando la estructura de una firma

Smart as you are, you decide to make, first of all, a model in equilibrium with constant sales of 100 cars per month. These sales empty the stock of inventory. The inventory initially consists of 300 cars and is replenished by a production inflow which equals the size of the workforce (i.e. the employees) times the productivity of an average worker. The average productivity is currently 1 car per person per month.

The target inventory currently equals the sales times an inventory coverage of 3 months. This target inventory is used to calculate the inventory correction which is used to calculate the target production. The inventory correction is calculated as: ( target inventory - inventory ) / time to correct inventory. The time to correct inventory is currently 2 months. The target production is then the sales plus the inventory correction.

The target production divided by the productivity of the average worker results in the target workforce. The discrepancy between the target workforce and the actual workforce currently drives the ‘hiring and firing’ strategy of the company as follows: net hire rate = ( target workforce - workforce ) / time to adjust workforce of 10 months. The workforce could be modeled as a stock variable regulated by the net hire rate. Since you will start out from equilibrium, you could take the target workforce as the initial value of the workforce.

Preguntas Iniciales del caso

3.1. Construye un modelo de dinámica de sistemas basado en el caso anterior, simula el modelo por un periodo de 100 meses. Asegúrate que el modelo este en equilibro (2 puntos, productos: modelo en R, envía tu modelo con la versión final de tu tarea, emplea el siguiente formato para nombrar tu modelo: TuNombre_cadena_suministro.R)

library(deSolve)

# Definición de la función del modelo que incluye la salida de Sales
model <- function(time, state, parameters) {
  with(as.list(c(state, parameters)), {
    
    # La producción es igual a la fuerza laboral multiplicada por la productividad
    Production <- Workforce * Productivity
    
    # Cálculos según las ecuaciones del modelo
    Target_Inventory <- Sales * Inventory_Coverage
    Inventory_Correction <- (Target_Inventory - Inventory) / Time_to_Correct_Inventory
    Target_Production <- Sales + Inventory_Correction
    Target_Workforce <- Target_Production / Productivity
    Net_Hire_Rate <- (Target_Workforce - Workforce) / Time_to_Adjust_Workforce
    
    # Tasas de cambio
    dInventory <- Production - Sales
    dWorkforce <- Net_Hire_Rate
    
    # Retorno de las tasas de cambio junto con Sales para visualización
    return(list(c(Inventory = dInventory, Workforce = dWorkforce), Sales = Sales))
  })
}

# Parámetros
parameters <- c(Productivity = 1, Inventory_Coverage = 3, Time_to_Correct_Inventory = 2,
                Time_to_Adjust_Workforce = 10, Sales = 100)

# Condiciones iniciales
state <- c(Inventory = 300, Workforce = 100)

# Tiempos para la simulación
times <- seq(0, 100, by = 1)

# Simulación del modelo
out <- ode(y = state, times = times, func = model, parms = parameters)

# Convertir resultados a un data frame para facilitar el manejo
results <- as.data.frame(out)

# Resultados de la simulación con Sales incluida en el gráfico
matplot(results$time, results[, c("Inventory", "Workforce", "Sales")], type = "l", col = c("blue", "red", "green"),
        xlab = "Meses", ylab = "Número de Personas", main = "Simulación de la Gestión de la Cadena de Suministro")
legend("topright", legend = c("Inventario", "Fuerza Laboral", "Ventas"), col = c("blue", "red", "green"), lty = 1)

Este modelo se mantiene en equilibrio a lo largo del tiempo.

Ventas:Las ventas se pueden ver con una línea constante a lo largo del tiempo, reflejando que se mantienen estables con un valor de 100 autos por mes ya que de esta manera configuramos el modelo.

Inventario: La linea del inventario se mantiene completamente plana y cosntante en 300 unidades a lo largo del tiempo ya que el sistema esta en equilibrio desde el inicio. Es decir, que la producción se ajusta perfectamente para compensar las ventas y mantener el inventario al nivel deseado. Esto sucede porque la producción y las ventas son iguales, ventas constantes = 100, producción = 100 al mes. La fuerza laboral es suficiente y constante, de manera que el inventario no muestra fluctuaciones ni disminuciones.

Fuerza Laboral: La fuerza laboral se mantiene constante en 100 personas, indicando que no existe la necesidad de ajustar la cantidad de empleados. Esta estabilidad en la fuerza laboral sugiere que el tamaño actual es el indicado para alcanzar las metas de producción requeridas para satisfacer las ventas y mantener el inventario al nivel deseado.

3.2. Ahora corre el modelo con especificando haciendo que la variable sales cambie de 100 a car per month a 150 cars per month. Muestra y describe el comportamiento del modelo para las siguientes variables: sales, inventory y workforce (3 puntos, productos: gráficos describiendo el comportamiento del sistema, comentarios que describan el comportamiento del sistema).

library(deSolve)

# Definición de la función del modelo que incluye Sales en la salida
model <- function(time, state, parameters) {
  with(as.list(c(state, parameters)), {
    
    # La producción es igual a la fuerza laboral multiplicada por la productividad
    Production <- Workforce * Productivity
    
    # Cálculos según las ecuaciones del modelo
    Target_Inventory <- Sales * Inventory_Coverage
    Inventory_Correction <- (Target_Inventory - Inventory) / Time_to_Correct_Inventory
    Target_Production <- Sales + Inventory_Correction
    Target_Workforce <- Target_Production / Productivity
    Net_Hire_Rate <- (Target_Workforce - Workforce) / Time_to_Adjust_Workforce
    
    # Tasas de cambio
    dInventory <- Production - Sales
    dWorkforce <- Net_Hire_Rate
    
    # Retorno de las tasas de cambio y Sales para visualización
    list(c(Inventory = dInventory, Workforce = dWorkforce), Sales = Sales)
  })
}

# Parámetros iniciales con ventas ajustadas desde el inicio
parameters <- c(Productivity = 1, Inventory_Coverage = 3, Time_to_Correct_Inventory = 2,
                Time_to_Adjust_Workforce = 10, Sales = 150)

# Condiciones iniciales ajustadas
state <- c(Inventory = 300, Workforce = 100)  # Ajuste según nuevas ventas

# Tiempos para la simulación
times <- seq(0, 100, by = 1)

# Simulación del modelo
out_cars_150 <- ode(y = state, times = times, func = model, parms = parameters)

# Convertir resultados a un data frame para facilitar el manejo
results_cars_150 <- as.data.frame(out_cars_150)

# Ajustando los márgenes para crear espacio para la leyenda fuera del gráfico
par(mar = c(5, 4, 4, 8) + 0.1)  # Aumenta el margen derecho para hacer espacio

# Creando el gráfico
matplot(results_cars_150$time, results_cars_150[, c("Inventory", "Workforce", "Sales")], type = "l", 
        col = c("blue", "red", "green"),
        xlab = "Meses", ylab = "Número de Unidades", main = "Simulación con Ventas Iniciales de 150")

# Colocando la leyenda fuera del área del gráfico, en el margen derecho extendido
legend("topright", inset = c(-0.3, 0), legend = c("Inventario", "Fuerza Laboral", "Ventas"), 
       col = c("blue", "red", "green"), lty = 1, xpd = TRUE)

Ventas: Las ventas se muestran como una línea recta sobre las 150 unidades puesto que esta cantidad se da como una constante en nuestro modelo.

Inventario: Los inventarios muestran una caída inicial rápida ya que al inicio no tiene suficiente para manejar el incremento inmediato en ventas (no sé inicio en equilibrio en este caso) y la producción (limitada inicialmente por la fuerza laboral disponible) no pudo satisfacer inmediatamente la nueva demanda. Posteriormente, los inventarios suben y se estabilizan, indicand oque la producción se ha ajustado para igualar y mantener el nuevo nivel de ventas. Esta estabilización es el resultado de los ajusted en la fuerza laboral que incrementaron la capacidad de producción.

Fuerza Laboral:La fuerza laboral muestra un incremento considerable al inicio en respuesta al déficit inicial en producción frenta a la nueva demanda. El aumento de la fuerza laboral permite a la empresa incrementar la producción para satisfacer las ventas más elevadas. Despues, la fuerza laboral parece estabilizarse a un nuevo nivel más alto que el inicial puesto que es necesario para mantener la producción adecuada al nivel aumentado de ventas. La estabilización y su mantenimiento indican que la empresa ha encontrado el balance necesario entre producción y demanda bajo las nuevas condiciones de venta.

3.3. ¿Cuánto tiempo toma para que el sistema vuelva al equilibrio dinámico? (5 puntos, productos: gráfico con anotaciones describiendo el tiempo que tarda el sistema en volver al equilibrio).

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
# Asegúrate de que results_cars_150 es el data frame correcto y está disponible
# Por ejemplo, si results_cars_150 es el resultado de una simulación, asegúrate de que ya esté convertido a data.frame
# Supongamos que results_cars_150 ya está correctamente formateado y nombrado

# Verificar las columnas disponibles en results_cars_150
if(!"Inventory" %in% names(results_cars_150) || !"Workforce" %in% names(results_cars_150)) {
  stop("Columnas 'Inventory' o 'Workforce' no encontradas en results_cars_150")
}

# Filtrar para encontrar el primer mes donde ambos, el inventario y la fuerza laboral, están cerca de los valores objetivo
stabilization_point <- results_cars_150 %>%
  filter(abs(Inventory - 450) < 5, abs(Workforce - 150) < 5) %>%
  summarize(EarliestStabilizationTime = min(time)) %>%
  pull(EarliestStabilizationTime) # Usamos pull para extraer el valor directamente

# Mostrar el mes de estabilización
stabilization_point
## [1] 55
print(stabilization_point)
## [1] 55
library(deSolve)
library(dplyr)

# Suponiendo que 'results_cars_150' contiene los datos de la simulación y está en formato de data.frame
# También asumimos que 'time', 'Inventory', y 'Workforce' están correctamente definidos en 'results_cars_150'

# Plot principal
plot(results_cars_150$time, results_cars_150$Inventory, type = "l", col = "blue",
     ylim = c(100, 600),  # Estableciendo límites para el eje Y desde 100 hasta 600
     xlab = "Meses", ylab = "Número de Unidades", main = "Tiempo de Equilibrio del Sistema",
     axes = TRUE, ann = TRUE)

# Añadir fuerza laboral al gráfico
lines(results_cars_150$time, results_cars_150$Workforce, col = "red")

# Añadir línea de ventas
lines(results_cars_150$time, rep(150, length(results_cars_150$time)), col = "green", lty = 2) 

# Leyenda para entender el gráfico
legend("topright", legend = c("Inventario", "Fuerza Laboral", "Ventas"),
       col = c("blue", "red", "green"), lty = 1, bty = "n")

# Si ya has determinado el punto de estabilización, puedes añadir una línea y texto
# Asegúrate de que 'stabilization_point' es un número y está definido
# Ejemplo: stabilization_point <- 40
stabilization_point <- 40 # Asumimos un valor para el punto de estabilización

abline(v = stabilization_point, col = "darkgrey", lwd = 2, lty = 2)
text(stabilization_point, 550, "Estabilización", pos = 4, cex = 0.8)

Toma 55 meses en que se empieza llegar al nuevo equilibrio dinámico.

A partir de dicho mes se puede observar como ambas variables reducen considerablemente su volatilidad y comienzan a ser más estables y cada variable empieza a osilar ligeramente arriba/abajo de su punto de equiulibrio, es decir, inventarios alrededor de 450 y fuerza laboral alrededor de 150.

3.4. Construye el diagrama de fase del sistema y describe el comportamiento esperado empleando este diagrama de fase (5 puntos, productos: diagrama de fase y explicación de comportamiento).

library(ggplot2)

# Suponiendo que los resultados de la simulación están en 'results'
# Actualizar el código para usar 'linewidth' en lugar de 'size'

ggplot(results_cars_150, aes(x = Inventory, y = Workforce)) +
  geom_path(linewidth = 1, color = "blue") +  # Usa 'linewidth' para especificar el grosor de la línea
  labs(x = "Inventario", y = "Fuerza Laboral", title = "Diagrama de Fase: Inventario vs. Fuerza Laboral") +
  theme_minimal() +
  annotate("point", x = results_cars_150$Inventory[1], y = results_cars_150$Workforce[1], colour = "red", size = 3) +
  annotate("point", x = results_cars_150$Inventory[nrow(results)], y = results_cars_150$Workforce[nrow(results_cars_150)], colour = "green", size = 3) +
  annotate("text", x = results_cars_150$Inventory[1], y = results_cars_150$Workforce[1] - 10, label = "Inicio", color = "red", hjust = 0.5, vjust = -1) +
  annotate("text", x = results_cars_150$Inventory[nrow(results_cars_150)], y = results_cars_150$Workforce[nrow(results_cars_150)] + 10, label = "Fin", color = "green", hjust = 0.5, vjust = 1)

El sistema comienza con un inventario relativamente alto en comparación con la fuerza laboral (x3 veces).

Sin embargo, esto fue cambiando a lo largo del tiempo. Apenas empezó la simulación el inventario fue disminuyendo mientras que aumentaba la fuerza laboral. Esto suecedió deibdo a que las ventas iniciales eran más elevadas de lo que la fuerza laboral podía producir inicialmente, provocando que el inventario cayera al principio y la fuerza laboral aumentará con el propósito de poder cubrir las necesidades de esta elevada demanda.

La trayectoria comienza a formar un bucle grande que gradualmente se estabiliza en un ciclo más pequeño, a medida que las oscilasiones se van reduciendo conforme el sistema se acerca a un estado de equilibrio dinámico. Las oscilasiones se presentan dados los ajustes en la fuerza laboal y en la producción en respuesta a las variaciones del inventario.

Cuando el ciclo llega al final de su trayectoria se puede observar que el sistema se aproxima a su equilibrio dinámico donde tanto el inventario como la fuerza laboral fluctúan alrededor de sus puntos estables. Este equilibrio refleja el ajuste a la fuerza laboral y los niveles de inventarios que se alinean con la demanda (ventas) y la capacidad de producción (oferta). En este caso, estos puntos oscilan cerca de 150 para la fuerza laboral y de 450 para los inventarios.

Este sistema muestra una respuesta de ajuste dinámico donde puede haber una reacción exagerada en término de reducción de inventario o aumento en la fuerza laboral, seguido de ajustes para alinear mejor la producción con la demanda y los inventarios deseados.

Después de algunos periodos de ajuste y oscilación, el sistema tiende a estabilizarse en un ciclo o punto de equilibrio donde las fluctuaciones son cada vez más mínimas y predecibles. Esto permite mejorar la eficiencia y efectividad de los recursos.

Modelando la estructura de la cadena de suministro

Now, suppose that the company you work for is part of a supply chain consisting of 3 manufacturing companies, all three with the (same) structure developed in the first part. In fact, you are working for the assembler who sells the final product to the final customers. This assembler is supplied by a tier one supplier, and the tier one supplier is supplied by the tier 2 supplier. The tier 2 supplier faces extremely undesirable oscillations in demand and production.

Preguntas adicionales del caso

3.5. Asume que la producción del ensamblador (i.e. assembler) constituye la demanda de la firma en el primer nivel de suministro (i.e. tier 1 supplier), y que la producción de la firma en la primera línea de suministro constituye la demanda de la firma en el segundo nivel de suministro (i.e. tier 2 supplier). Expande tu modelo original, asume que la estructura intérnate de estas firmas adicionales es idéntica a la estructura de la firma ensambladora (5 puntos, productos: modelo en R, envía tu modelo con la versión final de tu tarea, emplea el siguiente formato para nombrar tu modelo: TuNombre_cadena_suministro_completo.R).

library("deSolve")

model.cars.supply.chain <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
     
    # Endogenous auxiliary variables Assembler
    Target.Inventory <- Sales * Inventory.Coverage
    Inventory.Correction <- (Target.Inventory - Inventory) / Time.To.Correct.Inventory
    Target.Production <- Sales + Inventory.Correction
    Target.Workforce <- Target.Production / Productivity.Of.Average.Worker
  
    # Flow variables Assembler
    Production <- Workforce * Productivity.Of.Average.Worker
    Net.Hire.Rate <- (Target.Workforce - Workforce) / Time.To.Adjust.Workforce

    # State (stock) variables Assembler
    dInventory <- Production - Sales
    dWorkforce <- Net.Hire.Rate

    # Endogenous auxiliary variables Tier 1 supplier
    Target.Inventory1 <- Production * Inventory.Coverage
    Inventory.Correction1 <- (Target.Inventory1 - Inventory1) / Time.To.Correct.Inventory
    Target.Production1 <- Production + Inventory.Correction1
    Target.Workforce1 <- Target.Production1 / Productivity.Of.Average.Worker
    Sales1 <- Production

    # Flow variables Tier 1 supplier
    Production1 <- Workforce1 * Productivity.Of.Average.Worker
    Net.Hire.Rate1 <- (Target.Workforce1 - Workforce1) / Time.To.Adjust.Workforce

    # State (stock) variables Tier 1 supplier
    dInventory1 <- Production1 - Production
    dWorkforce1 <- Net.Hire.Rate1

    # Endogenous auxiliary variables Tier 2 supplier
    Target.Inventory2 <- Production1 * Inventory.Coverage
    Inventory.Correction2 <- (Target.Inventory2 - Inventory2) / Time.To.Correct.Inventory
    Target.Production2 <- Production1 + Inventory.Correction2
    Target.Workforce2 <- Target.Production2 / Productivity.Of.Average.Worker
    Sales2 <- Production1

    # Flow variables Tier 2 supplier
    Production2 <- Workforce2 * Productivity.Of.Average.Worker
    Net.Hire.Rate2 <- (Target.Workforce2 - Workforce2) / Time.To.Adjust.Workforce
  
    # State (stock) variables Tier 2 supplier
    dInventory2 <- Production2 - Production1
    dWorkforce2 <- Net.Hire.Rate2
      
    list(c(dInventory, 
           dInventory1, 
           dInventory2, 
           dWorkforce, 
           dWorkforce1, 
           dWorkforce2), 
         Sales = Sales, 
         Sales1 = Sales1, 
         Sales2 = Sales2)
  })
}

# Parameters
parameters <- c(Sales = 150, # Cars / Month
                Productivity.Of.Average.Worker = 1, # Car/person/month
                Inventory.Coverage = 3, # Months
                Time.To.Correct.Inventory = 2, # Months
                Time.To.Adjust.Workforce = 10 # Months
)

# Initial Conditions
InitialConditions <- c(Inventory = 300, # cars
                       Inventory1 = 300, # cars
                       Inventory2 = 300, # cars
                       Workforce = 100, # people
                       Workforce1 = 100, # people
                       Workforce2 = 100) # people
            
# Temporal Configuration
times <- seq(
  0,   # Initial time, Months
  200, # End time, Months
  1    # Time step, Months
)

# Método de Integración
intg.method <- c("rk4")
  
# Simulación
out <- ode(y = InitialConditions,
           times = times,
           func = model.cars.supply.chain,
           parms = parameters,
           method = intg.method)

# Gráfica individual de cada variable
variable_names <- colnames(out)[-1]
colors <- rainbow(length(variable_names))

par(mfrow = c(3, 3), mar = c(4, 4, 2, 1)) # Configurar la disposición de los gráficos en una cuadrícula de 3x3

for (i in 2:ncol(out)) {
  plot(out[, 1], out[, i], type = "l", col = colors[i-1], xlab = "Meses", ylab = variable_names[i-1], main = variable_names[i-1])
}

3.6. Describe la dinámica de comportamiento que resulta de integrar estas nuevas firmas al modelo para las siguientes variables: sales, inventory y workforce (5 puntos, productos: gráficos describiendo el comportamiento del sistema para las tres firmas, comentarios que describan el comportamiento del sistema).

library(ggplot2)
library(reshape2)

# Convertir los resultados a un data frame
results <- as.data.frame(out)
results$time <- out[, "time"]

# Crear datos largos para ggplot
results_long <- melt(results, id.vars = "time", variable.name = "variable", value.name = "value")

# Filtrar y crear los plots necesarios
# Gráfico de Inventarios
inventory_data <- results_long[grep("Inventory", results_long$variable),]
g_inventory <- ggplot(inventory_data, aes(x = time, y = value, color = variable)) +
  geom_line() +
  labs(title = "Dinámica de Inventario en la Cadena de Suministro", y = "Inventario", x = "Tiempo (meses)") +
  scale_color_manual(values = c("red", "orange", "green"), labels = c("Assembler", "Tier 1 Supplier", "Tier 2 Supplier")) +
  theme_minimal()

# Gráfico de Fuerza Laboral
workforce_data <- results_long[grep("Workforce", results_long$variable),]
g_workforce <- ggplot(workforce_data, aes(x = time, y = value, color = variable)) +
  geom_line() +
  labs(title = "Dinámica de la Fuerza Laboral en la Cadena de Suministro", y = "Fuerza Laboral", x = "Tiempo (meses)") +
  scale_color_manual(values = c("red", "orange", "green"), labels = c("Assembler", "Tier 1 Supplier", "Tier 2 Supplier")) +
  theme_minimal()

# Gráfico de Ventas
sales_data <- results_long[grep("Sales", results_long$variable),]
g_sales <- ggplot(sales_data, aes(x = time, y = value, color = variable)) +
  geom_line() +
  labs(title = "Dinámica de Ventas en la Cadena de Suministro", y = "Ventas", x = "Tiempo (meses)") +
  scale_color_manual(values = c("blue", "purple", "pink"), labels = c("Assembler", "Tier 1 Supplier", "Tier 2 Supplier")) +
  theme_minimal()

# Mostrar los gráficos
print(g_inventory)

print(g_workforce)

print(g_sales)

INVENTARIOS:

Tier 2 Supplier (verde): Este nivel muestra la mayor volatilidad en el inventario a lo largo del tiempo. Inicia con un aumento significativo, seguido de fluctuaciones pronunciadas que parecen estabilizarse parcialmente hacia el final del periodo simulado, aunque todavía con cierta variabilidad. Esta volatilidad podría indicar una respuesta directa a las fluctuaciones en la demanda proveniente del Tier 1 Supplier, que es su cliente directo.

Tier 1 Supplier (naranja): Este nivel muestra fluctuaciones menos pronunciadas en comparación con el Tier 2, pero aún significativas. Inicialmente, sigue una tendencia similar al Tier 2, lo que sugiere una dependencia directa de las fluctuaciones en su producción y las necesidades del ensamblador. A lo largo del tiempo, estas fluctuaciones disminuyen, indicando posiblemente una mejor adaptación o ajustes en la gestión de su inventario y producción.

Assembler (rojo): El nivel del ensamblador muestra las fluctuaciones más controladas y estables en su inventario. Inicia con ciertas oscilaciones que se suavizan considerablemente a medida que avanza el tiempo. Esto sugiere que el ensamblador puede estar mejor preparado para manejar variaciones en la demanda o que no depende tanto de variables externas.

FUERZA LABORAL:

Tier 2 Supplier (verde): Este nivel experimenta las fluctuaciones más extremas en la fuerza laboral a lo largo del tiempo, lo cual indica una alta sensibilidad a las condiciones de demanda y producción. Las grandes oscilaciones podrían ser indicativas de respuestas bruscas a los cambios en la producción del Tier 1 Supplier, lo cual puede llevar a ineficiencias y costos aumentados por la constante reajuste de personal.

Tier 1 Supplier (naranja): Aunque este nivel también muestra variabilidad significativa, es menos extrema que la del Tier 2 Supplier. Esto puede reflejar una mejor capacidad para amortiguar las fluctuaciones en la demanda proveniente del ensamblador. Sin embargo, sigue siendo susceptible a las variaciones y puede beneficiarse de una mejor planificación y estrategias de gestión de recursos humanos para minimizar los picos y valles.

Assembler (rojo): La fuerza laboral en el nivel del ensamblador es notablemente más estable que en los niveles de los proveedores, aunque muestra algunas fluctuaciones moderadas. Esto puede ser indicativo de que el ensamblador tiene mejores sistemas de gestión de la fuerza laboral o que está menos afectado por las variaciones inmediatas en la demanda en comparación con sus proveedores.

VENTAS

Assembler (azul): Las ventas del ensamblador son constantes y estables a lo largo del tiempo, manteniéndose en 150 unidades mensuales. Esto indica que la demanda de los productos finales del ensamblador es estable, lo cual es ideal para mantener una planificación de producción y fuerza laboral equilibrada.

Tier 1 Supplier (púrpura): Las ventas del Tier 1, que son esencialmente la producción del ensamblador, muestran pequeñas fluctuaciones que siguen de cerca las ventas del ensamblador, pero con una ligera volatilidad que puede ser el resultado de pequeños ajustes en la producción o en el inventario.

Tier 2 Supplier (rosa): Este nivel muestra una volatilidad extremadamente alta en las ventas, que son impulsadas por la producción del Tier 1 Supplier. Las grandes fluctuaciones pueden indicar que el Tier 2 es el más afectado por cualquier cambio en la demanda que se produce en niveles superiores de la cadena. Estas grandes oscilaciones pueden causar problemas significativos en la planificación de la producción y la gestión de la fuerza laboral.

3.7. Describe e implementa una política que ayude a disminuir la inestabilidad de la cadena de suministro. Compara en un mismo gráfico el comportamiento del sistema con y sin política (5 puntos, productos: gráficos describiendo el comportamiento del sistema para las tres firmas con y sin política, comentarios que describan el comportamiento del sistema).

library("deSolve")

model.cars.supply.chain <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    
    if (apply_policy) {
      # Políticas implementadas
      Time.To.Adjust.Workforce <- 5
      Time.To.Correct.Inventory <- 3
      Inventory.Coverage <- 2
    }

    # Cálculos del modelo para Assembler
    Target.Inventory <- Sales * Inventory.Coverage
    Inventory.Correction <- (Target.Inventory - Inventory) / Time.To.Correct.Inventory
    Target.Production <- Sales + Inventory.Correction
    Target.Workforce <- Target.Production / Productivity.Of.Average.Worker
    Production <- Workforce * Productivity.Of.Average.Worker
    Net.Hire.Rate <- (Target.Workforce - Workforce) / Time.To.Adjust.Workforce
    dInventory <- Production - Sales
    dWorkforce <- Net.Hire.Rate

    # Cálculos del modelo para Tier 1 Supplier
    Target.Inventory1 <- Production * Inventory.Coverage
    Inventory.Correction1 <- (Target.Inventory1 - Inventory1) / Time.To.Correct.Inventory
    Target.Production1 <- Production + Inventory.Correction1
    Target.Workforce1 <- Target.Production1 / Productivity.Of.Average.Worker
    Sales1 <- Production
    Production1 <- Workforce1 * Productivity.Of.Average.Worker
    Net.Hire.Rate1 <- (Target.Workforce1 - Workforce1) / Time.To.Adjust.Workforce
    dInventory1 <- Production1 - Production
    dWorkforce1 <- Net.Hire.Rate1

    # Cálculos del modelo para Tier 2 Supplier
    Target.Inventory2 <- Production1 * Inventory.Coverage
    Inventory.Correction2 <- (Target.Inventory2 - Inventory2) / Time.To.Correct.Inventory
    Target.Production2 <- Production1 + Inventory.Correction2
    Target.Workforce2 <- Target.Production2 / Productivity.Of.Average.Worker
    Sales2 <- Production1
    Production2 <- Workforce2 * Productivity.Of.Average.Worker
    Net.Hire.Rate2 <- (Target.Workforce2 - Workforce2) / Time.To.Adjust.Workforce
    dInventory2 <- Production2 - Production1
    dWorkforce2 <- Net.Hire.Rate2
    
    return(list(c(dInventory, dInventory1, dInventory2, dWorkforce, dWorkforce1, dWorkforce2),
                Sales = Sales, Sales1 = Sales1, Sales2 = Sales2))
  })
}

# Parameters and Initial Conditions as before

# Simulation without policy
parameters["apply_policy"] <- FALSE
out_no_policy <- ode(y = InitialConditions, times = times, func = model.cars.supply.chain, parms = parameters, method = intg.method)

# Simulation with policy
parameters["apply_policy"] <- TRUE
out_with_policy <- ode(y = InitialConditions, times = times, func = model.cars.supply.chain, parms = parameters, method = intg.method)

# Function to plot comparisons
plot_comparison <- function(data1, data2, variable) {
  plot(data1[, "time"], data1[, variable], type = "l", col = "red", ylim = range(c(data1[, variable], data2[, variable])), main = paste("Comparison of", variable))
  lines(data2[, "time"], data2[, variable], col = "blue")
  legend("topright", legend = c("Without Policy", "With Policy"), col = c("red", "blue"), lty = 1)
}

# Plot comparisons for each variable
plot_comparison(out_no_policy, out_with_policy, "Inventory")

plot_comparison(out_no_policy, out_with_policy, "Workforce")

plot_comparison(out_no_policy, out_with_policy, "Sales")

POLTICIAS IMPLEMENTADAS

  1. Reducción del Tiempo para Ajustar la Fuerza Laboral (Time.To.Adjust.Workforce):

De 10 a 5 meses: Este cambio significa que las empresas pueden responder más rápidamente a las necesidades de ajustar su fuerza laboral. Al poder contratar o despedir trabajadores en un período más corto, la empresa puede adaptarse más eficientemente a las demandas cambiantes del mercado o a las variaciones en la producción. Esto debería ayudar a minimizar los períodos de sobre o subempleo, reduciendo así los costos laborales innecesarios y mejorando la respuesta a las fluctuaciones de la demanda.

La implementación de las políticas resulta en una notable estabilización del inventario. Las fluctuaciones se suavizan considerablemente y el inventario se mantiene en un nivel más constante y predecible. Esto es indicativo de una respuesta más equilibrada y controlada a las variaciones en la demanda y producción, evitando los extremos de sobreproducción o escasez.

2.Aumento del Tiempo para Corregir el Inventario (Time.To.Correct.Inventory): De 2 a 3 meses: Al aumentar el tiempo que se toma para ajustar el inventario al nivel deseado, estás permitiendo una respuesta menos abrupta a los cambios en el inventario. Esto puede ayudar a evitar decisiones precipitadas que podrían llevar a una sobreproducción o escasez debido a reacciones exageradas a las fluctuaciones menores del mercado. Un enfoque más moderado y gradual para la gestión de inventarios puede contribuir a una mayor estabilidad en la cadena de suministro.

3.Reducción de la Cobertura de Inventario (Inventory.Coverage):

De 3 a 2 meses: Al disminuir la cantidad de inventario que se mantiene en relación con las ventas mensuales, la política está orientada a reducir el capital inmovilizado en inventarios no vendidos. Esto puede hacer que la empresa sea más ágil y menos susceptible a cambios en la demanda del consumidor. Una cobertura de inventario más baja también puede disminuir los costos de almacenamiento y reducir el riesgo de obsolescencia del producto.

Estos ajustes en conjunto están diseñados para hacer que la cadena de suministro sea más eficiente y menos propensa a los ciclos extremos de boom y caída, lo que debería resultar en una operación más suave y predecible, crucial para la planificación estratégica y financiera a largo plazo de la empresa.

RESULTADOS

Inventarios:La implementación de la política ha estabilizado notablemente los niveles de inventario. La gráfica muestra una línea casi plana, lo que indica que el inventario se mantiene constante y cerca del nivel objetivo a lo largo del tiempo. Esto sugiere que las medidas implementadas han sido efectivas para controlar las fluctuaciones de inventario.

Fuerza Laboral: Las políticas parecen haber reducido la amplitud de las fluctuaciones en la fuerza laboral, aunque todavía se observan algunos movimientos. Esto podría indicar que mientras la política ha ayudado, podría necesitar ajustes para optimizar aún más la estabilidad de la fuerza laboral.

Ventas: Las ventas permanecen constantes en ambas situaciones, lo que indica que las políticas implementadas no han tenido un impacto directo en las ventas. Sin embargo, la estabilidad en las ventas es un factor positivo, ya que permite a la cadena de suministro planificar con mayor certeza.

Problema 4: Planificación a largo plazo de nuevas ciudades (30 puntos)

Introducción:

The next few decades will be characterized by further urbanization. Part of it will –especially in Asia and Africa– be caused by so-called ‘New Towns’. The concept of ‘New Town’ is used for a particular type of urban development: a more or less independent new town relatively distant from existing cities. New urban developments within existing (older) urban areas (New Town in Town) or large urban extensions to the edges of the major cities is generally referred to as urban renewal, not as New Towns. The New Town concept is not new. Worldwide, there are hundreds of new towns (Brasilia, Chicago, San Francisco, Milton Keynes, many Eastern European, Russian, South American, Chinese, Korean and French cities) and many new ones are under construction (especially in Asia and Africa). Although the number of inhabitants of these cities varies from less than 50000 to more than a million, and although they are very different in design and function, there are interesting dynamic parallels and the same ‘mistakes’ –in spite of good examples like Detroit– are being made over and over again. One of the SD classics –Urban Dynamics by Jay W. Forrester (1969)– shows that SD is an appropriate method to help to understand New Town structures and dynamics and to test policies to solve undesirable evolutions. Suppose that you are hired as an adviser-modeler by a few rising new towns to test policies to solve the ‘root causes’ of undesirable dynamics. Since you are following this SD course, you decide to make and subsequently use a SD model. Use the information below, which is based on a modified version of George Richardson’s URBAN1 model.

Descripción del caso:

The population

The size of the population of a new town changes through immigration, births, emigration, and deaths. Suppose that the new town you are about to model already has 50000 inhabitants, and has a birth rate of 3%, and a death rate of 1.5%. And suppose that there is constant emigration with a normal emigration rate of 7% per year. Immigration could be modeled as the product of the current population, the normal immigration rate, the job availability multiplier for immigration, and the housing availability multiplier for immigration. This housing availability multiplier for immigration is a function of the households to houses ratio: suppose that if the households to houses ratio equals 1, the multiplier equals 1, that if the households to houses ratio equals 1.5, the multiplier equals 0.25, that if the households to houses ratio equals 2, the multiplier equals 0, that if the households to houses ratio equals 0, the multiplier equals 1.4, and that if the households to houses ratio equals 0.5, the multiplier equals 1.3. The number of households depends on the size of the population and the average size of households, in that part of the world still 4. Suppose that the normal immigration rate equals 10%.

Houses

The number of houses, initially equal to 14000, increases by means of construction of houses and decreases through demolition of houses. The average demolition rate of houses (without additional policies) equals 1.5% per year. The construction of houses could be modeled as a 3rd order delay (delayed with 2 years) of the product of the land availability multiplier for houses, the housing scarcity multiplier, the number of houses, and the construction rate of houses of 7% per year. The housing scarcity multiplier is a function of the households to houses ratio connecting following points (0, 0.2), (0.5, 0.3), (1, 1), (1.5, 1.7), and (2, 2). The land availability multiplier for houses is a function of the land fraction occupied: for a land fraction occupied of 0% the multiplier equals 0.4, for a land fraction occupied of 25% the multiplier equals 1, for a land fraction occupied of 50% the multiplier equals 1.5, for a land fraction occupied of 75% the multiplier equals 1, and for a land fraction occupied of 100% the multiplier equals 0. The land fraction occupied corresponds to the sum of the land use of all businesses and the land use of all houses, divided by the total area. Suppose that the useful total area of the new town is 5000 hectare, that the amount of land per house is 0.05 hectare, and that the amount of land per business (i.e. for each business structure) is 0.1 hectare.

Businesses and labor force

The number of businesses (i.e. business structures), initially 1000 business, increases through construction of business structures and decreases through demolition of business structures with an average demolition rate of business structures of 2.5% per year.

Construction of business structures could be modeled as the product of the land availability multiplier for business structures, the business labor force multiplier, the number of businesses, and the construction rate of business structures of 7% per year.

This land availability multiplier for business structures is a function of the land fraction occupied: the multiplier is of course 0 for 100% of the land fraction occupied, it equals 1 for 0% of the land fraction occupied, and it equals 1.5 for 50% of the land fraction occupied. The business labor force multiplier is a function of the labor force to jobs ratio connecting (0, 0.2), (0.5, 0.3), (1, 1), (1.5, 1.7), and (2, 2). The aforementioned job availability multiplier for immigration is also a function of the labor force to jobs ratio connecting following points (0, 2), (0.5, 1.75), (1, 1), (1.5, 0.25), and (2, 0.1). This labor force to jobs ratio depends of course on (i) the size of the labor force which equals the product of the population and the labor force to population ratio of 35%, and of (ii) the number of jobs which equals the number of businesses times the initial number of jobs per business structure which amounts to 18 FTEs. Model following two key performance indicators too: the unemployment ratio and the housing vacancy ratio. Both are per definition between 0 and 100%.

Preguntas del caso:

4.1. Construye un modelo de dinámica de sistemas basado en el caso anterior, simula el modelo por un período de 200 meses (5 puntos, productos: modelo en R)

# Declaración de Modelo
model <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {
  
# ENDOGENOUS AUXILIARY VARIABLES
    
# Population
    households <- population/average.size.households
    households.houses.ratio <- households/houses
    housing.availability.multiplier.immigration <- approx(c(1,1.5,2,0,0.5),
                                                          c(1,0.25,0,1.4,1.3),
                                                          xout=households.houses.ratio)$y
    
# Houses
    housing.scarcity.multiplier <- approx(c(0,0.5,1,1.5,2),
                                          c(0.2,0.3,1,1.7,2),
                                          xout=households.houses.ratio)$y
    land.use.all.businesses <- businesses*land.per.business
    land.use.all.houses <- houses*land.per.house
    land.fraction.occupied <- (land.use.all.businesses+land.use.all.houses)/total.area
    land.availability.multiplier.houses <- approx(c(0,0.25,0.5,0.75,1),
                                                  c(0.4,1,1.5,1,0),
                                                  xout=land.fraction.occupied)$y
    
# Businesses
    land.availability.multiplier.business.structures <- approx(c(1,0,0.5),
                                                               c(0,1,1.5),
                                                               xout=land.fraction.occupied)$y
    labor.force <- population*labor.force.to.population.ratio
    jobs <- businesses*initial.number.of.jobs.per.business.stucture
    labor.force.to.jobs.ratio <- labor.force/jobs
    business.labor.force.multiplier <- approx(c(0, 0.5, 1, 1.5, 2),
                                              c(0.2, 0.3, 1, 1.7, 2),
                                              xout=labor.force.to.jobs.ratio)$y
    job.availability.multiplier.immigration <- approx(c(0, 0.5, 1, 1.5, 2),
                                                      c(2, 1.75, 1, 0.25, 0.1),
                                                      xout=labor.force.to.jobs.ratio)$y
    
# FLOW VARIABLES
    
# Population
    emmigration <- population*normal.emmigration.rate
    births <- population*birth.rate
    deaths <- population*death.rate
    
# Houses
    demolition.of.houses <- houses*demolition.rate.of.houses
    construction.of.houses.actual <- (land.availability.multiplier.houses*housing.scarcity.multiplier*houses*construction.rate.of.houses)
    construction.of.houses.delayed <- construction.in.transit/average.time.of.delay
    
# Businesses
    demolition.of.business.structures <- businesses*demolition.rate.of.business.structures
    construction.of.business.structures <- land.availability.multiplier.business.structures*business.labor.force.multiplier*businesses*construction.rate.of.business.structures
    immigration <- population*normal.immigration.rate*job.availability.multiplier.immigration*housing.availability.multiplier.immigration
    
# Key Performance Indicators (KPI)
    unemployment.ratio <- ifelse(labor.force.to.jobs.ratio <= 1, 0, min(1.0, labor.force.to.jobs.ratio-1))
    housing.vacancy.ratio <- max(0, 1-households.houses.ratio)
    
# Stock variables
    
# Population
    dpopulation <- immigration+births-emmigration-deaths
    
# Houses
    dhouses <- construction.of.houses.actual-demolition.of.houses
    
# Businesses
    dbusinesses <- construction.of.business.structures-demolition.of.business.structures
    
# Material delays
    dconstruction.in.transit <- construction.of.houses.actual-construction.of.houses.delayed
    
    return(list(c(dpopulation, dhouses, dbusinesses, dconstruction.in.transit)))
  })
}
#Parameters

parameters <- c(average.time.of.delay=2,
                birth.rate=.03,
                death.rate=.015,
                normal.emmigration.rate=.07,
                normal.immigration.rate=.1,
                average.size.households=4,
                demolition.rate.of.houses=.015,
                construction.rate.of.houses=.07,
                total.area=5000,
                land.per.house=.05,
                land.per.business=.1,
                labor.force.to.population.ratio=.35,
                initial.number.of.jobs.per.business.stucture=18,
                demolition.rate.of.business.structures=.025,
                construction.rate.of.business.structures=.07)

#Condiciones Iniciales
InitialConditions <- c(population=50000,
                       houses=14000,
                       businesses=1000,
                       construction.in.transit=0)

# Configuración Temporal
times <- seq(
  0,   # Initial time, Months
  200, # End time, Months
  1    # Time step, Months
)

# Método de Integración
intg.method <- "rk4"

# Simulación del modelo
out <- ode(
  y = InitialConditions,
  times = times,
  func = model,
  parms = parameters,
  method = intg.method
)

plot(out)

4.2. Realiza gráficos de la evolución de las empresas, las viviendas y la población y de los efectos de éstos en la tasa de desempleo (unemployment ratio) y la tasa de viviendas vacías (housing vacancy ratio). ¿Qué problemas puedes derivar de estos gráficos? (5 puntos, productos: gráficos y comentarios del comportamiento)

# Declaración de Modelo
model<-function(t, state, parameters) {
  with(as.list(c(state,parameters)), {

#ENDOGENOUS AUXILIARY VARIABLES

#Population
  households<-population/average.size.households
  households.houses.ratio<-households/houses
  housing.availability.multiplier.immigration<-approx(c(1,1.5,2,0,0.5),
                                                        c(1,0.25,0,1.4,1.3),
                                                        xout=households.houses.ratio)$y
    
#Houses
  housing.scarcity.multiplier<-approx(c(0,0.5,1,1.5,2),
                                        c(0.2,0.3,1,1.7,2),
                                        xout=households.houses.ratio)$y
  land.use.all.businesses<-businesses*land.per.business
  land.use.all.houses<-houses*land.per.house
  land.fraction.occupied<-(land.use.all.businesses+land.use.all.houses)/total.area
  land.availability.multiplier.houses<-approx(c(0,0.25,0.5,0.75,1),
                                                c(0.4,1,1.5,1,0),
                                                xout=land.fraction.occupied)$y
    
#Businesses
  land.availability.multiplier.business.structures<-approx(c(1,0,0.5),
                                                             c(0,1,1.5),
                                                             xout=land.fraction.occupied)$y
  labor.force<-population*labor.force.to.population.ratio
  jobs<-businesses*initial.number.of.jobs.per.business.stucture
  labor.force.to.jobs.ratio<-labor.force/jobs
  business.labor.force.multiplier<-approx(c(0, 0.5, 1, 1.5, 2),
                                            c(0.2, 0.3, 1, 1.7, 2),
                                            xout=labor.force.to.jobs.ratio)$y
  job.availability.multiplier.immigration<-approx(c(0, 0.5, 1, 1.5, 2),
                                                    c(2, 1.75, 1, 0.25, 0.1),
                                                    xout=labor.force.to.jobs.ratio)$y

#FLOW VARIABLES

#Population
  emmigration<-population*normal.emmigration.rate
  births<-population*birth.rate
  deaths <-population*death.rate

#Houses
  demolition.of.houses<-houses*demolition.rate.of.houses
  construction.of.houses.actual<-(land.availability.multiplier.houses*housing.scarcity.multiplier*houses*construction.rate.of.houses)
  construction.of.houses.delayed<-construction.in.transit/average.time.of.delay

  #Businesses
  demolition.of.business.structures<-businesses*demolition.rate.of.business.structures
  construction.of.business.structures<-land.availability.multiplier.business.structures*business.labor.force.multiplier*businesses*construction.rate.of.business.structures
  immigration<-population*normal.immigration.rate*job.availability.multiplier.immigration*housing.availability.multiplier.immigration

#Key Performance Indicators (KPI)
  unemployment.ratio<-ifelse(labor.force.to.jobs.ratio<=1,0,min(1.0,labor.force.to.jobs.ratio-1))
  housing.vacancy.ratio<-max(0,1-households.houses.ratio)


#Stock variables

#Population
  dpopulation<-immigration+births-emmigration-deaths
    
#Houses
  dhouses<-construction.of.houses.actual-demolition.of.houses

#Businesses
  dbusinesses<-construction.of.business.structures-demolition.of.business.structures

#Material delays
  dconstruction.in.transit<-construction.of.houses.actual-construction.of.houses.delayed

list(c(dpopulation, dhouses, dbusinesses,dconstruction.in.transit),
       unemployment.ratio=unemployment.ratio,
       housing.vacancy.ratio=housing.vacancy.ratio)
  })
}

# Parameters
parameters<-c(average.time.of.delay=2,
              birth.rate=.03,
              death.rate=.015,
              normal.emmigration.rate=.07,
              normal.immigration.rate=.1,
              average.size.households=4,
              demolition.rate.of.houses=.015,
              construction.rate.of.houses=.07,
              total.area=5000,
              land.per.house=.05,
              land.per.business=.1,
              labor.force.to.population.ratio=.35,
              initial.number.of.jobs.per.business.stucture=18,
              demolition.rate.of.business.structures=.025,
              construction.rate.of.business.structures=.07)

InitialConditions <-c(population=50000,
                      houses=14000,
                      businesses=1000,
                      construction.in.transit=0)

#Configuración Temporal
times <- seq(
  0,   # Initial time, Months
  200, # End time, Months
  1    # Time step, Months
)

#Metodo de Integración
intg.method<-c("rk4")

#Simulación del modelo
out <- ode(
  y = InitialConditions,
  times = times,
  func = model,
  parms = parameters,
  method = intg.method
)

# Gráfico de resultados
plot(out, col = c("blue"))

  1. Población: La población muestra una disminución inicial significativa antes de estabilizarse. Esto podría indicar un alto índice de emigración o mortalidad inicial que luego se estabiliza debido a cambios en las políticas de inmigración o mejoras en las condiciones de vida. La caída inicial de la población puede afectar la base de consumidores y la fuerza laboral disponible, lo cual podría tener efectos negativos en la economía local.

  2. Viviendas: El número de viviendas también disminuye inicialmente y luego se nivela. Este comportamiento podría ser resultado de una tasa de demolición que supera la construcción nueva o de una falta inicial de inversión en construcción de viviendas.La disminución de viviendas puede exacerbar problemas de alojamiento y llevar a un aumento en el costo de vivienda debido a la escasez.

  3. Negocios: Similar a la población y las viviendas, el número de negocios disminuye y luego muestra fluctuaciones antes de estabilizarse. Esto puede reflejar retos económicos que llevan al cierre de negocios o a una inversión conservadora en nuevos negocios. La reducción en el número de negocios puede llevar a un aumento del desempleo y una reducción en la producción económica local.

  4. Construcción en Tránsito: Hay un pico inicial en la construcción en tránsito que rápidamente desciende a cero. Esto sugiere una respuesta inicial robusta a la necesidad de construcción que se completa rápidamente. El rápido descenso en la construcción podría indicar una falta de proyectos futuros, afectando el crecimiento económico a largo plazo.

  5. Tasa de Desempleo: El desempleo inicialmente es alto, pero disminuye rápidamente y se estabiliza en un nivel bajo. Este patrón es probablemente el resultado de la estabilización de la población y el ajuste en el número de negocios.Aunque la tasa de desempleo baja, el impacto inicial alto puede tener efectos duraderos en la comunidad.

  6. Tasa de Viviendas Vacías: La tasa de viviendas vacías muestra un aumento significativo en la mitad del período antes de empezar a disminuir, lo cual puede ser resultado de las fluctuaciones en la construcción de viviendas y la demografía cambiante. Un alto número de viviendas vacías puede llevar a la degradación de propiedades y barrios, reduciendo el valor del inmueble en la región.

El modelo sugiere que la comunidad está enfrentando desafíos significativos en términos de mantener una población estable, asegurar suficientes viviendas y negocios, y manejar el desempleo y las viviendas vacías. Las políticas efectivas deberían enfocarse en estimular la economía, apoyar la construcción de viviendas y negocios, y fomentar una fuerza laboral adecuada para las necesidades de los negocios. Esto podría incluir incentivos para nuevos negocios, programas para el desarrollo de vivienda asequible, y estrategias para aumentar la inmigración o reducir la emigración.

4.3. Forrester’s Urban Dynamics shows that additional demolition of empty bad quality housing in case of high housing vacancy ratios and rezoning could prevent urban decay to kick in or worsen. Thus add a variable additional demolition rate to demolish beyond normal demolishing above a 10% housing vacancy ratio. Let the additional demolition rate increase linearly from 0% per year for a 10% vacancy ratio to 5% per year for a 15% vacancy ratio, and linearly from 5% per year for a 15% vacancy ratio to 50% per year for a 100% housing vacancy ratio. Add this additional effect to the demolition of houses too (30 puntos, producto: modelo actualizado en R [10 puntos], gráficos [10 puntos] y conclusiones del comportamiento de las variables y el sistema [10 puntos])

library("deSolve")

# Declaración de Modelo
model_func <- function(t, state, params) {
  with(as.list(c(state, params)), {
    
    # Population
    households <- population / average.size.households
    households.houses.ratio <- households / houses
    housing.availability.multiplier.immigration <- approx(
      x = c(1, 1.5, 2, 0, 0.5),
      y = c(1, 0.25, 0, 1.4, 1.3),
      xout = households.houses.ratio
    )$y
    
    # Houses
    housing.scarcity.multiplier <- approx(
      x = c(0, 0.5, 1, 1.5, 2),
      y = c(0.2, 0.3, 1, 1.7, 2),
      xout = households.houses.ratio
    )$y
    land.use.all.businesses <- businesses * land.per.business
    land.use.all.houses <- houses * land.per.house
    land.fraction.occupied <- (land.use.all.businesses + land.use.all.houses) / total.area
    
    land.availability.multiplier.houses <- approx(
      x = c(0, 0.25, 0.5, 0.75, 1),
      y = c(0.4, 1, 1.5, 1, 0),
      xout = land.fraction.occupied
    )$y
    
    # Businesses
    land.availability.multiplier.business.structures <- approx(
      x = c(1, 0, 0.5),
      y = c(0, 1, 1.5),
      xout = land.fraction.occupied
    )$y
    
    labor.force <- population * labor.force.to.population.ratio
    jobs <- businesses * initial.number.of.jobs.per.business.structure
    labor.force.to.jobs.ratio <- labor.force / jobs
    
    business.labor.force.multiplier <- approx(
      x = c(0, 0.5, 1, 1.5, 2),
      y = c(0.2, 0.3, 1, 1.7, 2),
      xout = labor.force.to.jobs.ratio
    )$y
    
    job.availability.multiplier.immigration <- approx(
      x = c(0, 0.5, 1, 1.5, 2),
      y = c(2, 1.75, 1, 0.25, 0.1),
      xout = labor.force.to.jobs.ratio
    )$y
    
    # FLOW VARIABLES
    
    # Population
    emmigration <- population * normal.emmigration.rate
    births <- population * birth.rate
    deaths <- population * death.rate
    
    # Houses
    demolition.of.houses <- houses * demolition.rate.of.houses
    construction.of.houses.actual <- land.availability.multiplier.houses * housing.scarcity.multiplier * houses * construction.rate.of.houses
    construction.of.houses.delayed <- construction.in.transit / average.time.of.delay
    
    # Businesses
    demolition.of.business.structures <- businesses * demolition.rate.of.business.structures
    construction.of.business.structures <- land.availability.multiplier.business.structures * business.labor.force.multiplier * businesses * construction.rate.of.business.structures
    immigration <- population * normal.immigration.rate * job.availability.multiplier.immigration * housing.availability.multiplier.immigration
    
    # Key Performance Indicators (KPI)
    housing.vacancy.ratio <- max(0, 1 - households.houses.ratio)
    
    if (housing.vacancy.ratio > 0.1) {
      additional.demolition.rate <- approx(
        x = c(0.1, 0.15, 1),
        y = c(0, 0.05, 0.5),
        xout = housing.vacancy.ratio
      )$y
    } else {
      additional.demolition.rate <- 0
    }
    
    additional.demolition <- houses * additional.demolition.rate
    total.house.demolition <- demolition.of.houses + additional.demolition
    
    unemployment.ratio <- ifelse(labor.force.to.jobs.ratio <= 1, 0, min(1, labor.force.to.jobs.ratio - 1))
    
    # Stock variables
    # Population
    dpopulation <- immigration + births - emmigration - deaths
    
    # Houses
    dhouses <- construction.of.houses.actual - total.house.demolition
    
    # Businesses
    dbusinesses <- construction.of.business.structures - demolition.of.business.structures
    
    # Material delays
    dconstruction.in.transit <- construction.of.houses.actual - construction.of.houses.delayed
    
    list(c(dpopulation, dhouses, dbusinesses, dconstruction.in.transit), 
         unemployment.ratio = unemployment.ratio, 
         housing.vacancy.ratio = housing.vacancy.ratio, demolition.of.business.structures=demolition.of.business.structures)
  })
}

# Parámetros
params <- list(
  average.time.of.delay = 2,
  birth.rate = 0.03,
  death.rate = 0.015,
  normal.emmigration.rate = 0.07,
  normal.immigration.rate = 0.1,
  average.size.households = 4,
  demolition.rate.of.houses = 0.015,
  construction.rate.of.houses = 0.07,
  total.area = 5000,
  land.per.house = 0.05,
  land.per.business = 0.1,
  labor.force.to.population.ratio = 0.35,
  initial.number.of.jobs.per.business.structure = 18,
  demolition.rate.of.business.structures = 0.025,
  construction.rate.of.business.structures = 0.07
)

# Condiciones iniciales
init_state <- c(
  population = 50000,
  houses = 14000,
  businesses = 1000,
  construction.in.transit = 0
)

# Configuración Temporal
times <- seq(
  0,   # Initial time, Months
  200, # End time, Months
  1    # Time step, Months
)

# Metodo de Integración
intg.method <- "rk4"

# Simulación del modelo
out <- ode(
  y = init_state,
  times = times,
  func = model_func,
  parms = params,
  method = intg.method
)

# Gráfico de resultados
plot(out, col = c("darkblue"))

  1. Población: La población muestra un descenso inicial significativo antes de estabilizarse. Esto podría ser resultado de altas tasas de emigración o mortalidad al principio del período simulado, o posiblemente por un efecto de las variables de inmigración que no compensan las salidas.

  2. Viviendas: El número de viviendas inicialmente aumenta y luego se estabiliza. Esto sugiere que las políticas de construcción están respondiendo efectivamente a la necesidad inicial, pero después el crecimiento se ralentiza, posiblemente debido a la limitación de espacio disponible o a una disminución en la necesidad de nuevas viviendas.

  3. Empresas: Similar a las viviendas, el número de empresas inicialmente aumenta y luego experimenta un declive antes de estabilizarse. Este declive puede deberse a un aumento en la demolición o a un ambiente de negocios que se vuelve menos favorable, lo cual podría estar influenciado por la disminución de la población y por tanto una reducción en la demanda de servicios y productos.

  4. Construcción en tránsito: Hay un pico inicial que luego cae rápidamente, lo que indica que los proyectos de construcción que estaban en proceso se completan. La rápida disminución podría ser un indicativo de una gestión eficiente en la finalización de los proyectos de construcción.

  5. Tasa de desempleo: La tasa de desempleo se mantiene baja y relativamente estable, lo que sugiere que hay un buen equilibrio entre la fuerza laboral y los empleos disponibles. Esto es un indicador positivo de salud económica en la simulación.

  6. Tasa de viviendas vacías: Este indicador muestra una volatilidad considerable, con picos que pueden indicar períodos de sobrecapacidad donde el número de viviendas supera la demanda. La presencia de picos también podría indicar problemas temporales en el mercado de vivienda, como podría ser un retardo en la reacción del mercado a los cambios en la población o en la economía.

  7. Demolición de estructuras de negocios: La tasa de demolición muestra una tendencia decreciente, lo cual podría estar alineado con las políticas de reducción de construcción o con una estabilización en el sector de negocios después de un ajuste inicial.

El sistema modelado muestra una dinámica compleja influenciada por la interacción entre políticas de población, vivienda y negocios. Para optimizar la gestión urbana y prevenir problemas como el exceso de viviendas vacías o la inestabilidad en el sector empresarial, sería crucial implementar políticas adaptativas que puedan responder dinámicamente a las condiciones cambiantes del mercado y las necesidades de la población.