## (PASO 1)
#Instalar librerias
library("deSolve")
library(gridExtra)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:gridExtra':
## 
##     combine
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyr)
library(ggplot2)

Plantilla para modelación de dinámica de sistemas

Como primer paso haremos un listado del tipo de variables en el modelo:

Variables de Estado (Stock variables) Déficit Hídrico Industrial Déficit Hídrico Agrícola Déficit Hídrico Doméstico

Variables de flujo (flow variables) demanda.excedente.industrial demanda.excedente.agricola demanda.excedente.domestica control.consumo.industrial control.consumo.agricola control.consumo.domestico

Variables auxiliares endógenas (endogenous auxiliary variables)

Parámetros de simulación (variables en la frontera del sistema o exogenous auxiliary variables) Demanda Industrial Demanda Agrícola Demanda Doméstica Tecnificación y Modernización Industrial Tecnificación y Modernización Agrícola Tecnificación y Modernización Doméstico Tasa de déficit global

Esta clasificación es útil para organizar el espacio de trabajo en Rstudio.

Chunk con todos los elementos para modelar el caso.

library(deSolve)

# -------------------------------
# 1. Parámetro de personalización (Matrícula A00836033)
# -------------------------------
X <- 3

# -------------------------------
# 2. Condiciones iniciales
# -------------------------------
InitialConditions <- c(
  deficit.industrial = 20,
  deficit.agricola   = 50,
  deficit.domestica  = 15
)

# -------------------------------
# 3. Parámetros exógenos
# -------------------------------
parameters <- c(
  demanda.industrial       = X + 5,   # 8
  demanda.agricola         = X + 10,  # 13
  demanda.domestica        = X + 3,   # 6
  tasa.deficit             = 2,
  tecnificacion.industrial = 0.06,
  tecnificacion.agricola   = 0.09,
  tecnificacion.domestica  = 0.04
)

# -------------------------------
# 4. Horizonte temporal
# -------------------------------
times <- seq(0, 40, by = 1/3)

# -------------------------------
# 5. Función del modelo
# -------------------------------
modelo.hidrico <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    
    # ---------------------------
    # Variables de flujo de entrada
    # ---------------------------
    demanda.excedente.industrial <- demanda.industrial * tasa.deficit
    demanda.excedente.agricola   <- demanda.agricola * tasa.deficit
    demanda.excedente.domestica  <- demanda.domestica * tasa.deficit
    
    # ---------------------------
    # Variables de flujo de salida
    # ---------------------------
    control.consumo.industrial <- deficit.industrial * tecnificacion.industrial
    control.consumo.agricola   <- deficit.agricola * tecnificacion.agricola
    control.consumo.domestica  <- deficit.domestica * tecnificacion.domestica
    
    # ---------------------------
    # Ecuaciones diferenciales
    # ---------------------------
    dDeficit.industrial <- demanda.excedente.industrial - control.consumo.industrial
    dDeficit.agricola   <- demanda.excedente.agricola   - control.consumo.agricola
    dDeficit.domestica  <- demanda.excedente.domestica  - control.consumo.domestica
    
    # ---------------------------
    # Variable auxiliar endógena
    # ---------------------------
    deficit.hidrico.total <- deficit.industrial + deficit.agricola + deficit.domestica
    
    # ---------------------------
    # Salida del modelo
    # ---------------------------
    list(
      c(dDeficit.industrial, dDeficit.agricola, dDeficit.domestica),
      demanda.excedente.industrial = demanda.excedente.industrial,
      demanda.excedente.agricola   = demanda.excedente.agricola,
      demanda.excedente.domestica  = demanda.excedente.domestica,
      control.consumo.industrial   = control.consumo.industrial,
      control.consumo.agricola     = control.consumo.agricola,
      control.consumo.domestica    = control.consumo.domestica,
      deficit.hidrico.total        = deficit.hidrico.total
    )
  })
}

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

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

# Convertir a data frame
out_df <- as.data.frame(out)

# -------------------------------
# 8. Gráfica básica
# -------------------------------
plot(out, col = "darkblue")

out_df      <- as.data.frame(out)
out_df$anio <- out_df$time

library(ggplot2)
library(dplyr)
library(tidyr)

sectores_plot <- out_df %>%
  select(anio, deficit.industrial, deficit.agricola, deficit.domestica) %>%
  pivot_longer(
    cols = -anio,
    names_to = "sector",
    values_to = "deficit"
  ) %>%
  mutate(
    sector = recode(
      sector,
      "deficit.industrial" = "Industrial",
      "deficit.agricola"   = "Agrícola",
      "deficit.domestica"  = "Doméstico"
    )
  )

ultimos_puntos <- sectores_plot %>%
  group_by(sector) %>%
  filter(anio == max(anio))

ggplot(sectores_plot, aes(x = anio, y = deficit, group = sector)) +
  geom_line(aes(color = sector), linewidth = 1.5, alpha = 0.95) +
  geom_point(
    data = ultimos_puntos,
    aes(color = sector),
    size = 3
  ) +
  geom_text(
    data = ultimos_puntos,
    aes(label = round(deficit, 1), color = sector),
    hjust = -0.2,
    size = 3.8,
    show.legend = FALSE
  ) +
  geom_vline(
    xintercept = 12,
    linetype = "22",
    linewidth = 0.8,
    color = "gray10"
  ) +
  annotate(
    "label",
    x = 12,
    y = max(sectores_plot$deficit) * 0.18,
    label = "Año 12:\nfin del periodo de observación",
    size = 3.2,
    label.size = 0.2,
    fill = "white",
    color = "gray10"
  ) +
  scale_color_manual(
    values = c(
      "Industrial" = "darkblue",
      "Agrícola"   = "darkgreen",
      "Doméstico"  = "darkred"
    )
  ) +
  coord_cartesian(xlim = c(0, 42)) +
  labs(
    title = "Evolución del déficit hídrico por sector",
    subtitle = "Región ficticia La Seca | Simulación de 40 años con intervalos de 4 meses",
    x = "Tiempo (años)",
    y = "Déficit hídrico acumulado",
    color = "Sector",
    caption = "Los valores al final de cada curva muestran el déficit alcanzado al año 40."
  ) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(face = "bold", size = 16),
    plot.subtitle = element_text(size = 11, color = "gray30"),
    axis.title = element_text(face = "bold"),
    legend.position = "top",
    legend.title = element_text(face = "bold"),
    panel.grid.minor = element_blank(),
    panel.grid.major.x = element_blank(),
    plot.caption = element_text(color = "gray40", hjust = 0)
  )
## Warning in annotate("label", x = 12, y = max(sectores_plot$deficit) * 0.18, :
## Ignoring unknown parameters: `label.size`

Déficit Total

ultimo_total <- out_df %>%
  filter(anio == max(anio))

ggplot(out_df, aes(x = anio, y = deficit.hidrico.total)) +
  geom_area(alpha = 0.18, fill = "#6A4C93") +
  geom_line(linewidth = 1.6, color = "#4A266A") +
  geom_point(
    data = ultimo_total,
    size = 3.2,
    color = "#4A266A"
  ) +
  geom_text(
    data = ultimo_total,
    aes(label = round(deficit.hidrico.total, 1)),
    hjust = -0.2,
    size = 3.8,
    color = "#4A266A"
  ) +
  geom_vline(
    xintercept = 12,
    linetype = "22",
    linewidth = 0.8,
    color = "gray45"
  ) +
  annotate(
    "segment",
    x = 12, xend = 12,
    y = 0, yend = max(out_df$deficit.hidrico.total) * 0.35,
    color = "gray45",
    linewidth = 0.7
  ) +
  annotate(
    "label",
    x = 14.2,
    y = max(out_df$deficit.hidrico.total) * 0.28,
    label = "La política se observa\nformalmente hasta el año 12",
    size = 3.2,
    fill = "white",
    color = "gray25",
    label.size = 0.2
  ) +
  coord_cartesian(xlim = c(0, 42)) +
  labs(
    title = "Trayectoria del déficit hídrico total",
    subtitle = "Indicador agregado de presión hídrica regional",
    x = "Tiempo (años)",
    y = "Déficit hídrico total",
    caption = "El sombreado resalta la acumulación del déficit a lo largo del horizonte de simulación."
  ) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(face = "bold", size = 16),
    plot.subtitle = element_text(size = 11, color = "gray30"),
    axis.title = element_text(face = "bold"),
    panel.grid.minor = element_blank(),
    panel.grid.major.x = element_blank(),
    plot.caption = element_text(color = "gray40", hjust = 0)
  )
## Warning in annotate("label", x = 14.2, y = max(out_df$deficit.hidrico.total) *
## : Ignoring unknown parameters: `label.size`