3 Estudio de Caso: Calibración Termodinámica, Eléctrica y Cinética de Degradación
3.1 Motivación y Contexto del Estudio
El presente capítulo consolida el marco teórico desarrollado previamente, trasladándolo de una formulación matemática abstracta a un escenario de aplicación altamente crítico: la red de frío para la preservación de vacunas en Tuxtla Gutiérrez, Chiapas. La elección de esta región geográfica como caso de estudio se debe a la combinación de condiciones climatológicas de calor extremo y a la vulnerabilidad de la infraestructura eléctrica local ante fenómenos meteorológicos. Bajo este entorno, los refrigeradores institucionales que resguardan los biológicos de alto valor inmunológico operan bajo condiciones de estrés continuo.
Para garantizar la validez y la utilidad del modelo computacional, se evitará la parametrización arbitraria de la Ecuación Diferencial Estocástica (EDE). Un modelo probabilístico con tasas de intensidad o volatilidades carentes de fundamento físico produciría trayectorias desconectadas de la realidad operativa. Por consiguiente, el objetivo central de este capítulo es la calibración rigurosa de las variables del sistema.
El proceso de estimación y asignación de parámetros se desarrollará en las siguientes tres fases analíticas:
- Caracterización Estocástica del Suministro Eléctrico: Mediante el análisis de datos oficiales de la Comisión Federal de Electricidad (CFE), se modelará la frecuencia y la duración de las interrupciones del servicio eléctrico. Estas estimaciones permitirán estructurar la conmutación de la cadena de Markov y los Procesos de Poisson No Homogéneos acoplados a la ecuación.
- Calibración Termodinámica del Equipo Operativo: A partir de las especificaciones técnicas y requerimientos normados por la Organización Mundial de la Salud (OMS), se deducirán analíticamente los coeficientes de reversión a la media (\(\kappa\)) y las volatilidades térmicas (\(\sigma\)) que dictan el enfriamiento activo y la inercia pasiva del refrigerador de estudio.
- Cinética de Degradación Biológica: Finalmente, se vinculará el comportamiento térmico dada por la EDE con la Ecuación de Arrhenius, estableciendo el marco cuantitativo necesario para evaluar la pérdida de viabilidad farmacológica que sufren las vacunas durante algún fallo.
3.2 Datos Base e Histórico de Interrupciones Eléctricas
El objetivo de esta fase es modelar la tasa de ocurrencia y la escala temporal de las interrupciones del suministro eléctrico en el área geográfica de Tuxtla Gutiérrez, Chiapas. Para evitar supuestos teóricos alejados de la realidad operativa de la red local, los datos de partida no se asignaron de manera arbitraria, sino que se obtuvieron a través de los registros oficiales de la Comisión Federal de Electricidad (CFE). Mediante una solicitud de información pública de transparencia dirigida a la CFE, se obtuvo el registro histórico del promedio mensual de interrupciones y su duración para el municipio de Tuxtla Gutiérrez.
A partir de dicha solicitud, se recabó el comportamiento histórico del municipio para el periodo comprendido entre los años 2019 y 2025. La Tabla 3.1 contiene dicha información, presentando el número de interrupciones promedio por mes y la duración promedio en horas por cada corte de energía.
| Año | Número de interrupciones promedio por mes | Duración Promedio de las interrupciones (Horas) |
|---|---|---|
| 2019 | 25 | 2.02 |
| 2020 | 42 | 1.73 |
| 2021 | 47 | 2.03 |
| 2022 | 40 | 2.59 |
| 2023 | 49 | 2.23 |
| 2024 | 43 | 2.39 |
| 2025 | 38 | 2.02 |
El análisis de la Tabla 3.1 muestra un incremento significativo en la frecuencia de fallas a partir del año 2020, estabilizándose en una franja superior a los 40 cortes mensuales. Este cambio se debe al crecimiento de la carga térmica residencial e industrial y a la saturación de los transformadores de distribución en la mancha urbana.
Para la calibración de la EDE en este trabajo de investigación, se seleccionó el año 2024 como el escenario de línea base. Dicho año presenta una media mensual de 43 interrupciones por mes con una duración promedio por corte de 2.39 horas. Al proyectar esta tasa mensual a una tasa anual, se tiene que el número total de interrupciones al año esperadas a nivel municipal (\(N_{total}\)) es:
\[ N_{total} = 43 \text{ cortes/mes} \times 12 \text{ meses} = 516 \text{ cortes/año} \tag{3.1}\]
3.3 Corrección de Escala Espacial (Factor \(\kappa_{espacial}\))
3.3.1 Justificación de la Reducción de Escala
El valor de \(516\) cortes al año obtenido en la Ecuación 3.1 es un dato agregado de todo el municipio, así que utilizar directamente este valor para un solo refrigerador implicaría asumir que el centro de almacenamiento de vacunas sufre más de un apagón diario, lo cual no es físicamente real. Este error ocurre al ignorar que un edificio está conectado a un solo nodo de la red, no a todos simultáneamente.
Para solucionar esto, introducimos un factor de escala espacial. De acuerdo con el Programa para el Desarrollo del Sistema Eléctrico Nacional (PRODESEN) y los diagramas de infraestructura del Centro Nacional de Control de Energía (CENACE) (Secretaría de Energía 2024), la mancha urbana de Tuxtla Gutiérrez está alimentada de manera sectorizada por exactamente 8 subestaciones principales de distribución (Tuxtla Céntrica, Oriente, Poniente, Terán, Tuxtla Norte, Tuxtla Dos, Mactumatzá y Real del Bosque).
3.3.2 La Esperanza Matemática Local
Asumiendo que los cortes de energía se distribuyen de manera uniforme sobre estas subestaciones, la probabilidad \(p\) de que un corte afecte a la subestación de interés (Tuxtla Norte, encargada de alimentar la zona de hospitales y laboratorios de la Colonia Maya) es \(p = 1/8\).
Por consiguiente, el número de cortes que tiene el centro de almacenamiento de vacunas se comporta como una variable aleatoria binomial \(N_{local} \sim \mathcal{B}(N_{total}, p)\), cuya esperanza matemática es:
\[ \mathbb{E}[N_{local}] = N_{total} \times p = \frac{516 \text{ cortes}}{8 \text{ subestaciones}} = 64.5 \approx 64 \text{ cortes al año} \tag{3.2}\]
Al redondear la esperanza a un número entero, establecemos una base realista de 64 interrupciones esperadas al año para nuestra simulación.
3.4 Distribución Temporal de las Fallas
Las interrupciones no ocurren de manera homogénea a lo largo del año (es decir, no ocurren exactamente 5.33 cortes cada mes). La variabilidad mensual depende directamente de los factores climáticos estacionales de la región. Para no asignar los apagones de cada mes de manera aleatoria, utilizamos el comportamiento documentado en el Informe Anual de la CFE 2024 (Comisión Federal de Electricidad 2024), específicamente en su página 82, donde se registra la evolución de los indicadores de continuidad del servicio a lo largo de los meses.
Para comprender la estructura de estos datos oficiales, la CFE evalúa la continuidad del suministro a través de dos indicadores internacionales de confiabilidad:
3.4.1 1. Índice de la Frecuencia de Interrupción Promedio, en número (SAIFI)
El indicador SAIFI (System Average Interruption Frequency Index) contabiliza el número promedio de interrupciones sostenidas que experimenta un usuario durante un periodo determinado (generalmente un año o un mes), es decir si se tiene una red eléctrica con un total de \(N_T\) usuarios atendidos. Si en un intervalo de observación ocurren \(M\) eventos de interrupción, donde cada evento \(i \in \{1, 2, \dots, M\}\) afecta a una cantidad de \(N_i\) usuarios, el índice SAIFI se expresa formalmente como:
\[ \displaystyle \text{SAIFI} = \frac{\sum_{i=1}^{M} N_i}{N_T} \tag{3.3}\]
3.4.2 2. Índice de Duración Promedio de Interrupción, en minutos (SAIDI)
El indicador SAIDI (System Average Interruption Duration Index) mide la duración total acumulada de las interrupciones sostenidas que experimenta el usuario promedio dentro del periodo de evaluación. Si cada evento de falla \(i\) tiene una duración temporal de \(r_i\) horas, el índice SAIDI se define analíticamente mediante:
\[ \displaystyle \text{SAIDI} = \frac{\sum_{i=1}^{M} N_i r_i}{N_T} \tag{3.4}\]
A partir de estas definiciones, la Tabla 3.2 presenta la evolución mensual acumulada de los indicadores SAIDI y SAIFI reportada en la página 82 del informe anual de la CFE.
| Mes | SAIDI Acumulado (Horas) | SAIFI Acumulado |
|---|---|---|
| Enero | 0.998 | 0.013 |
| Febrero | 1.110 | 0.017 |
| Marzo | 1.285 | 0.022 |
| Abril | 2.793 | 0.027 |
| Mayo | 4.120 | 0.037 |
| Junio | 5.774 | 0.053 |
| Julio | 6.557 | 0.065 |
| Agosto | 7.905 | 0.080 |
| Septiembre | 8.592 | 0.096 |
| Octubre | 9.063 | 0.103 |
| Noviembre | 9.314 | 0.106 |
| Diciembre | 10.812 | 0.113 |
Dado que el informe de la CFE publica el indicador SAIFI de forma acumulada, es necesario calcular la diferencia marginal (\(\Delta\)) respecto al mes anterior para obtener el peso estadístico de cada mes de forma individual:
\[ \Delta \text{SAIFI}_i = \text{SAIFI}_i - \text{SAIFI}_{i-1} \tag{3.5}\]
Para distribuir los 64 cortes anuales calculados para la subestación local respetando esta estacionalidad, aplicamos una regla de tres simple basada en la proporcionalidad marginal. Suponiendo que el valor total del SAIFI al final del año (\(0.113\)) representa el 100% de las fallas (equivalente a los 64 cortes), entonces el incremento mensual generado en cada periodo (\(\Delta \text{SAIFI}_i\)) representará una fracción directamente proporcional de cortes. Expresado de forma general para cualquier mes \(i\), la asignación viene dada por:
\[ \text{Cortes}_i = \left( \frac{\Delta \text{SAIFI}_i}{\text{SAIFI}_{total}} \right) \times \mathbb{E}[N_{local}] \tag{3.6}\]
Sustituyendo los valores de nuestro escenario (\(\text{SAIFI}_{total} = 0.113\) y \(\mathbb{E}[N_{local}] = 64\)), la ecuación del modelo es:
\[ \text{Cortes}_i = \left( \frac{\Delta \text{SAIFI}_i}{0.113} \right) \times 64 \tag{3.7}\]
Al hacer este cálculo y redondeando al entero más cercano (realizando un ajuste para conservar la suma total en 64 cortes anuales), obtenemos la distribución estacional de las fallas eléctricas presentada en la Tabla 3.3.
| Mes | \(\Delta\) SAIFI (Mensual) | Cortes Exactos Calculados | Cortes Asignados (Redondeo) |
|---|---|---|---|
| Enero | 0.013 | 7.36 | 7 |
| Febrero | 0.004 | 2.26 | 2 |
| Marzo | 0.005 | 2.83 | 3 |
| Abril | 0.005 | 2.83 | 3 |
| Mayo | 0.010 | 5.66 | 6 |
| Junio | 0.016 | 9.06 | 9 |
| Julio | 0.012 | 6.79 | 7 |
| Agosto | 0.015 | 8.49 | 8 |
| Septiembre | 0.016 | 9.06 | 9 |
| Octubre | 0.007 | 3.96 | 4 |
| Noviembre | 0.003 | 1.70 | 2 |
| Diciembre | 0.007 | 3.96 | 4 |
El análisis de la Tabla 3.3 revela con claridad los dos momentos estacionales de mayor vulnerabilidad y frecuencia de cortes en el suministro eléctrico de la región. El pico más severo se registra en Junio (9 cortes), un fenómeno ocasionado por las altas temperaturas extremas registradas en Tuxtla Gutiérrez durante la primavera-verano, las cuales provocaron una saturación en la red de distribución eléctrica debido al uso masivo y simultáneo de sistemas de aire acondicionado. El segundo periodo crítico ocurre durante Agosto y Septiembre (8 y 9 cortes respectivamente), coincidiendo directamente con el máximo histórico de la temporada de lluvias, tormentas tropicales y huracanes en el Sureste mexicano, eventos climáticos que suelen derribar líneas de alta tensión y afectar transformadores locales. En contraste, la época invernal se consolida como el periodo de máxima estabilidad en el servicio.
3.5 Modelado del Proceso de Poisson No Homogéneo (PPNH)
Para integrar las interrupciones eléctricas dentro de la Ecuación Diferencial Estocástica (EDE) en tiempo continuo \(t\), donde \(t \in [0, 365]\) representa los días del año, utilizamos un Proceso de Poisson No Homogéneo, cuya intensidad \(\lambda(t)\) varía en función del tiempo para tener en cuenta la estacionalidad climática de la región.
El modelo se construye mediante la suma de una tasa base de fondo y dos funciones de perturbación gaussianas no normalizadas. La elección de núcleos gaussianos permite modelar la transición gradual de las estaciones meteorológicas sin introducir discontinuidades de salto, concentrando el riesgo exactamente en los días donde el clima somete a la red a mayor estrés.
La función analítica propuesta para la intensidad diaria de fallas es:
\[ \lambda(t) = \lambda_{base} + A_1 \exp\left(-\frac{(t - \mu_1)^2}{2\sigma_1^2}\right) + A_2 \exp\left(-\frac{(t - \mu_2)^2}{2\sigma_2^2}\right) \tag{3.8}\]
3.5.1 Calibración mediante la Regla Empírica (68-95-99.7)
Para deducir la dispersión exacta de las campanas (\(\sigma_1\) y \(\sigma_2\)), recurrimos a la regla empírica de la distribución normal, la cual establece que el 95.4% de los datos se concentra a una distancia de \(\pm 2\sigma\) respecto a la media. Igualando esta amplitud (un ancho total de \(4\sigma\)) a la duración física en días de cada temporada climática en Chiapas, obtenemos la siguiente calibración analítica:
- Tasa Base (\(\lambda_{base} = 0.0667\)): Representa el riesgo mínimo constante por fallas accidentales o mantenimientos de la red, independiente del clima. Se calcula a partir de los meses más estables (febrero y noviembre, con 2 cortes al mes), resultando en \(2 \text{ cortes} / 30 \text{ días} \approx 0.0667\) cortes/día.
- Ola de Calor (\(\mu_1 = 152\) y \(\sigma_1 = 45\)): El primer pico térmico se centra a principios de junio (día 152). La temporada de calor abarca históricamente 6 meses (de marzo a agosto, \(\approx 180\) días). Aplicando la regla empírica, \(4\sigma_1 = 180 \implies \sigma_1 = 45\). Elevando al cuadrado, la varianza estacional queda en \(2\sigma_1^2 = 4050\).
- Temporada de Lluvias (\(\mu_2 = 244\) y \(\sigma_2 = 30\)): El segundo pico se centra a inicios de septiembre (día 244), el momento más severo de tormentas. Esta temporada abarca 4 meses (julio a octubre, \(\approx 120\) días). Aplicando la regla empírica, \(4\sigma_2 = 120 \implies \sigma_2 = 30\). Su varianza estacional queda en \(2\sigma_2^2 = 1800\).
3.5.2 Sistema de Amplitudes Cruzadas (\(A_1\) y \(A_2\))
Para que los picos máximos en \(t=152\) y \(t=244\) alcancen exactamente una intensidad de \(0.30\) cortes/día (equivalente a los 9 cortes mensuales reportados empíricamente en junio y septiembre), no basta con restar la tasa base de \(0.0667\). Al estar las campanas cercanas en el tiempo, la “cola” derecha del pico de calor contribuye a la intensidad durante la temporada de lluvias, y viceversa.
Para obtener las amplitudes exactas (\(A_1\) y \(A_2\)), se resuelve el siguiente sistema de ecuaciones evaluado en los centros \(\mu_1\) y \(\mu_2\):
\[ \begin{cases} \lambda(152) = 0.0667 + A_1 (1) + A_2 \exp\left(-\frac{(152 - 244)^2}{1800}\right) = 0.30 \\ \lambda(244) = 0.0667 + A_1 \exp\left(-\frac{(244 - 152)^2}{4050}\right) + A_2 (1) = 0.30 \end{cases} \]
Evaluando los términos exponenciales: \[ \begin{cases} A_1 + 0.009 A_2 = 0.2333 \\ 0.123 A_1 + A_2 = 0.2333 \end{cases} \]
Al resolver este sistema, obtenemos los valores de las amplitudes: \(A_1 = 0.2314\) y \(A_2 = 0.2048\). Sustituyendo todos los coeficientes, la función de intensidad es:
\[ \lambda(t) = 0.0667 + 0.2314 \exp\left(-\frac{(t - 152)^2}{4050}\right) + 0.2048 \exp\left(-\frac{(t - 244)^2}{1800}\right) \tag{3.9}\]
Como se aprecia en la Figura 3.1, el área bajo la curva de la función de intensidad \(\lambda(t)\) se divide en cuatro regiones climáticas fundamentales para reflejar el comportamiento térmico de Tuxtla Gutiérrez:
- Estabilidad Invernal (Área Gris, Días \(0 \le t < 62\) y \(304 < t \le 366\)): Abarca desde el 1 de enero hasta principios de marzo, y del 1 de noviembre al 31 de diciembre. Representa el periodo de menor demanda térmica y mayor estabilidad operacional de la red eléctrica, donde la intensidad responde únicamente a la tasa base de fallas accidentales o mantenimientos de la red (\(\lambda_{base} = 0.0667\)).
- Temporada de Calor y Estiaje (Área Naranja, Días \(62 \le t < 196\)): Comprende del 2 de marzo al 15 de julio, alcanzando su pico máximo en \(t = 152\) (1 de junio). En este periodo, las temperaturas extremas provocan la saturación de los transformadores por el uso intensivo de aire acondicionado.
- La Canícula (Área Púrpura, Días \(196 \le t < 227\)): Abarca aproximadamente 30 días, desde el 15 de julio hasta el 15 de agosto. Corresponde a la zona de transición e intersección entre ambas campanas gaussianas. Este fenómeno meteorológico se caracteriza por una disminución temporal de las precipitaciones y un repunte del calor extremo a mitad del verano, lo que somete al sistema eléctrico a un estrés continuo antes de la llegada de las tormentas severas.
- Temporada de Tormentas y Ciclones (Área Azul, Días \(227 \le t \le 304\)): Comprende del 15 de agosto al 31 de octubre, con su punto crítico centrado en \(t = 244\) (1 de septiembre). Representa la máxima incidencia de precipitación y eventos meteorológicos extremos capaces de derribar líneas de alta tensión.
3.6 Validación del Modelo y Análisis de la Anomalía de Enero
Por los fundamentos de los Procesos de Poisson No Homogéneos, el número esperado de cortes en cualquier intervalo \((t_1, t_2]\) está determinado por la integral definida de su función de intensidad:
\[ \mathbb{E}[N(t_1, t_2)] = \int_{t_1}^{t_2} \lambda(t) dt \tag{3.10}\]
Para validar la consistencia analítica de nuestra función \(\lambda(t)\), particionamos el año de 366 días (considerando que el año 2024 fue bisiesto, por lo que febrero aporta 29 días al dominio) en 12 subintervalos correspondientes a los días acumulados de cada mes. Al resolver numéricamente sus respectivas integrales mediante cuadratura gaussiana, contrastamos el área teórica bajo la curva contra los objetivos empíricos de la CFE:
| Mes | Intervalo de Integración \([t_{i-1}, t_i]\) | Objetivo Empírico (CFE) | Esperanza calculada \(\mathbb{E}[N]\) |
|---|---|---|---|
| Enero | \([0, 31]\) | 7 | 2.15 |
| Febrero | \([31, 60]\) | 2 | 2.37 |
| Marzo | \([60, 91]\) | 3 | 3.82 |
| Abril | \([91, 121]\) | 3 | 6.12 |
| Mayo | \([121, 152]\) | 6 | 8.73 |
| Junio | \([152, 182]\) | 9 | 8.74 |
| Julio | \([182, 213]\) | 7 | 8.39 |
| Agosto | \([213, 244]\) | 8 | 9.20 |
| Septiembre | \([244, 274]\) | 9 | 7.70 |
| Octubre | \([274, 305]\) | 4 | 4.27 |
| Noviembre | \([305, 335]\) | 2 | 2.31 |
| Diciembre | \([335, 366]\) | 4 | 2.09 |
La suma acumulada de las 12 integrales arroja un total de 65.90 cortes esperados a lo largo del año. Este resultado demuestra una congruencia matemática robusta frente a los 64 cortes teóricos calculados inicialmente. La diferencia de \(+2\) cortes anuales funciona como un margen de seguridad metodológico, garantizando que el modelo numérico no subestime el estrés real sobre la red de frío.
Finalmente, es importante destacar la inexactitud observada exclusivamente en el mes de Enero (7 cortes empíricos frente a los \(2.15\) esperados por la integral). Dado que enero es un periodo climatológicamente estable (sin huracanes ni olas de calor), esta desviación es originada muy probablemente por mantenimientos extraordinarios o contingencias locales impredecibles ocurridas específicamente a inicios del año 2024.
3.7 Distribución Log-Normal de las Duraciones de las Fallas eléctricas
Una vez modelada la frecuencia con la que ocurren las interrupciones eléctricas mediante el Proceso de Poisson, el siguiente componente crítico es determinar la duración exacta de cada apagón (\(D\)). En teoría de probabilidad aplicada a la ingeniería de confiabilidad, el tiempo de reparación o restablecimiento de un sistema no sigue a una distribución normal estándar.
La duración de una falla eléctrica es una variable aleatoria continua que presenta las siguientes dos características principales:
- Dominio estrictamente positivo: El tiempo de interrupción no puede tomar valores negativos (\(D > 0\)).
- Asimetría positiva (Cola pesada hacia la derecha): La inmensa mayoría de las fallas locales se resuelven en pocas horas (e.g., restablecimiento de fusibles o reconexiones de subestación), pero existe una probabilidad baja de eventos atípicos severos cuya reparación demora días enteros (e.g., explosión de transformadores o colapso de torres por huracanes).
Debido a que una distribución Normal clásica asignaría probabilidad a tiempos negativos y subestimaría el riesgo de cortes prolongados, la literatura establece que las duraciones de las interrupciones eléctricas se modelan mediante una Distribución Log-normal:
\[ D \sim \text{LogNormal}(\mu, \sigma^2) \tag{3.11}\]
Para calibrar el parámetro \(\mu\) (media de los logaritmos) y el parámetro \(\sigma\) (desviación estándar de los logaritmos) para cada temporada, es necesario utilizar los indicadores de la CFE.
3.7.1 Desglose Mensual de Horas de Interrupción
De acuerdo con el registro histórico municipal consolidado para el año 2024, la duración promedio de un corte en Tuxtla Gutiérrez fue de \(2.39\) horas por corte Tabla 3.1. Sabiendo que el factor de escala espacial redujo el riesgo de la unidad de salud a un promedio de \(64\) cortes al año, podemos calcular la cantidad total de tiempo que el refrigerador clínico experimentará desabasto eléctrico a lo largo de todo el año:
\[ \text{Horas Totales Anuales} = 64 \text{ cortes/año} \times 2.39 \text{ horas/corte} = 152.96 \text{ horas} \]
Para distribuir estas \(152.96\) horas a lo largo de los 12 meses del año, respetando la intensidad de daño a la infraestructura propia de cada época del año, utilizamos el indicador de duración acumulada SAIDI del Informe Anual de la CFE 2024 Tabla 3.2.
El primer paso consiste en calcular la diferencia marginal (\(\Delta \text{SAIDI}_i = \text{SAIDI}_i - \text{SAIDI}_{i-1}\)) para obtener el peso de cada mes individual. Posteriormente, aplicamos una regla de 3 simple: asumiendo que el valor total del SAIDI acumulado total (\(10.812\) horas) representa el 100% de la afectación anual (equivalente a nuestras \(152.96\) horas proyectadas), entonces el incremento de cualquier mes \(i\) representará una fracción directamente proporcional de horas de corte para el centro de almacenamiento. Analíticamente, la ecuación de asignación es:
\[ \text{Horas del Mes}_i = \left( \frac{\Delta \text{SAIDI}_i}{\text{SAIDI}_{total}} \right) \times \text{Horas Totales Anuales} \]
Sustituyendo los valores de nuestro escenario (\(\text{SAIDI}_{total} = 10.812\) y \(\text{Horas Totales} = 152.96\)), el modelo de asignación mensual queda como:
\[ \text{Horas del Mes}_i = \left( \frac{\Delta \text{SAIDI}_i}{10.812} \right) \times 152.96 \]
Finalmente, para conocer la duración media esperada de un corte individual en un mes específico (\(E_i\)), se divide el número de horas asignadas a dicho mes entre la cantidad de cortes de ese mes calculados previamente en la Tabla 3.3. Es decir, \(E_i = \text{Horas del Mes}_i / \text{Cortes}_i\). Al hacer estos cálculos se obtienen los siguientes resultados:
| Mes | SAIDI Acumulado | \(\Delta\) SAIDI | Duración Acumulada de las Interrupciones en el Mes | Cortes del Mes | Duración Promedio \(E_i\) (Horas) |
|---|---|---|---|---|---|
| Enero | 0.998 | 0.998 | 14.12 | 7 | 2.017 |
| Febrero | 1.110 | 0.112 | 1.58 | 2 | 0.790 |
| Marzo | 1.285 | 0.175 | 2.48 | 3 | 0.827 |
| Abril | 2.793 | 1.508 | 21.33 | 3 | 7.110 |
| Mayo | 4.120 | 1.327 | 18.77 | 6 | 3.128 |
| Junio | 5.774 | 1.654 | 23.40 | 9 | 2.600 |
| Julio | 6.557 | 0.783 | 11.08 | 7 | 1.583 |
| Agosto | 7.905 | 1.348 | 19.07 | 8 | 2.384 |
| Septiembre | 8.592 | 0.687 | 9.72 | 9 | 1.080 |
| Octubre | 9.063 | 0.471 | 6.66 | 4 | 1.665 |
| Noviembre | 9.314 | 0.251 | 3.55 | 2 | 1.775 |
| Diciembre | 10.812 | 1.498 | 21.19 | 4 | 5.298 |
3.7.2 Cálculo de Medias y Varianzas por Temporada Climática
Para integrar con precisión la duración de los cortes dentro de los regímenes estacionales de la EDE, es necesario agrupar los datos mensuales de la Tabla 3.5 en las tres temporadas climáticas previamente delimitadas en el análisis del Proceso de Poisson. La conformación dichos grupos se debe a las siguientes razones termodinámicas:
- Temporada de Calor (Abril a Julio): Agrupa el periodo de máximo estiaje y temperaturas extremas, culminando en la canícula. Durante este cuatrimestre, la saturación térmica de los equipos de distribución eléctrica es máxima.
- Temporada de Lluvias (Agosto a Octubre): Concentra el periodo crítico de huracanes y tormentas, caracterizado por daños mecánicos sobre la red eléctrica.
- Estabilidad Invernal (Enero a Marzo, Noviembre y Diciembre): Meses de baja exigencia térmica ambiental, donde los apagones ocurren mayoritariamente a factores de mantenimiento.
Realizar un promedio aritmético simple de las duraciones \(E_i\) entre los meses de cada temporada introduciría un sesgo grave, ya que asignaría el mismo peso a un mes con 9 cortes que a uno con solo 2 cortes. Para garantizar la rigurosidad estadística, se calcula la media empírica (\(\mathbb{E}[D]\)) y la varianza empírica (\(\mathbb{V}[D]\)) mediante un modelo de estadística descriptiva ponderada, donde el peso de cada mes (\(\omega_i\)) está dado por su frecuencia de cortes (\(\omega_i = \text{Cortes}_i\)).
Asi, las formulaciones para la esperanza y la varianza ponderada están dadas por:
\[ \mathbb{E}[D] = \frac{\sum_{i=1}^{n} \omega_i E_i}{\sum_{i=1}^{n} \omega_i} = \frac{\text{Horas Acumuladas de la Temporada}}{\text{Cortes Totales de la Temporada}} \tag{3.12}\]
\[ \mathbb{V}[D] = \frac{\sum_{i=1}^{n} \omega_i (E_i - \mathbb{E}[D])^2}{\sum_{i=1}^{n} \omega_i} \tag{3.13}\]
Al sustituir los valores de los meses correspondientes en las Ecuación 3.12 y Ecuación 3.13, se obtienen los siguientes resultados para cada temporada:
- Temporada de Calor (Abr, May, Jun, Jul): Acumula un total de \(25\) cortes y \(74.58\) horas de interrupción, arrojando una duración media alta de \(\mathbb{E}[D] = 2.983\) horas y una varianza considerable de \(\mathbb{V}[D] = 2.651\), esto se debe a las anomalías de colapso de red registradas durante las olas de calor primaverales.
- Temporada de Lluvias (Ago, Sep, Oct): Acumula \(21\) cortes pero únicamente \(35.45\) horas de interrupción. Esto nos da una media notablemente más baja de \(\mathbb{E}[D] = 1.688\) horas y una varianza muy estable de \(\mathbb{V}[D] = 0.343\), indicando que las fallas por lluvia son frecuentes pero suelen ser restablecidas con rapidez.
- Estabilidad Invernal (Ene, Feb, Mar, Nov, Dic): Acumula \(18\) cortes y \(42.92\) horas de interrupción. Su duración media es de \(\mathbb{E}[D] = 2.384\) horas con una varianza alta de \(\mathbb{V}[D] = 2.667\), provocada por mantenimientos atípicos observados en los bordes del año (especialmente en enero y diciembre).
3.7.3 Transformación Algebraica a Parámetros Log-Normales
Una vez conocidas las medias (\(\mathbb{E}[D]\)) y varianzas (\(\mathbb{V}[D]\)), estas deben transformarse hacia los parámetros \(\mu\) y \(\sigma\) que conforman la función de distribución Log-normal.
La teoría de la probabilidad nos dice que la esperanza y la varianza de una variable con distribución Log-normal dependen de sus parámetros a través del siguiente sistema de ecuaciones no lineales:
\[ \mathbb{E}[D] = \exp\left(\mu + \frac{\sigma^2}{2}\right) \] \[ \mathbb{V}[D] = \left[\exp(\sigma^2) - 1\right] \exp(2\mu + \sigma^2) \]
Sustituyendo la primera ecuación en la segunda, la varianza puede reescribirse como \(\mathbb{V}[D] = (\mathbb{E}[D])^2 \left[\exp(\sigma^2) - 1\right]\). A partir de esta simplificación, se despeja el sistema para obtener las fórmulas inversas de calibración:
\[ \sigma = \sqrt{\ln\left(1 + \frac{\mathbb{V}[D]}{(\mathbb{E}[D])^2}\right)} \tag{3.14}\]
\[ \mu = \ln(\mathbb{E}[D]) - \frac{\sigma^2}{2} = \ln\left( \frac{(\mathbb{E}[D])^2}{\sqrt{\mathbb{V}[D] + (\mathbb{E}[D])^2}} \right) \tag{3.15}\]
Al sustituir las medias y varianzas empíricas de cada temporada en las ecuaciones anteriores, obtenemos los parámetros exactos de la distribución Log-normal que modelará la duración de las fallas eléctricas en cada temporada:
| Temporada Climática | Meses que integran lo Integran | Media Ponderada \(\mathbb{E}[D]\) | Varianza Ponderada \(\mathbb{V}[D]\) | Parámetro \(\mu\) | Parámetro \(\sigma\) |
|---|---|---|---|---|---|
| Estabilidad (Invierno) | Ene, Feb, Mar, Nov, Dic | 2.384 | 2.667 | 0.6766 | 0.6202 |
| Calor (Estiaje y Canícula) | Abr, May, Jun, Jul | 2.983 | 2.651 | 0.9627 | 0.5106 |
| Lluvias (Tormentas) | Ago, Sep, Oct | 1.688 | 0.343 | 0.4668 | 0.3372 |
3.7.4 Visualización de las Densidades de Probabilidad (PDF)
Para comprender físicamente el comportamiento de los apagones en cada época del año, es fundamental analizar la Función de Densidad de Probabilidad (PDF, por sus siglas en inglés) de las tres distribuciones calibradas anteriormente. Matemáticamente, la función de densidad para una variable aleatoria \(D\) con distribución Log-normal se define como:
\[ f(x; \mu, \sigma) = \frac{1}{x \sigma \sqrt{2\pi}} \exp\left( - \frac{(\ln x - \mu)^2}{2\sigma^2} \right), \quad \text{para } x > 0 \tag{3.16}\]
Donde \(x\) representa la duración del corte eléctrico en horas, y el área bajo la curva en un intervalo \([a, b]\) representa la probabilidad exacta de que un apagón dure entre \(a\) y \(b\) horas.
A continuación, se grafican las densidades de probabilidad para cada temporada climática, utilizando los parámetros \(\mu\) y \(\sigma\) deducidos en la Tabla 3.6. Para facilitar la comparación visual del riesgo, todas las gráficas se evalúan en una ventana de observación de 0 a 12 horas.
3.8 Planteamiento del Modelo
En la literatura de ingeniería térmica y control, el comportamiento de la temperatura de una cámara frigorífica suele modelarse utilizando el enfoque de Ecuaciones Diferenciales Estocásticas de “caja gris” (Grey-box SDEs). Al representar el circuito térmico equivalente del refrigerador y añadir un término de ruido estocástico para absorber la incertidumbre de las fluctuaciones no medidas, la ecuación diferencial resultante tiene la forma de un proceso de Ornstein-Uhlenbeck. Este enfoque ha sido validado exitosamente en la estimación predictiva de unidades de refrigeración en trabajos como el de (Costanzo et al. 2013).
Sin embargo, para el problema logístico de la red de frío en Tuxtla Gutiérrez, el modelo clásico resulta limitado. Si bien el proceso de Ornstein-Uhlenbeck describe a la perfección la reversión a la media (la acción estabilizadora del termostato) y el ruido térmico interno de fondo, es insuficiente para modelar alteraciones drásticas en las condiciones de frontera.
Por este motivo, en la presente investigación se añadió a la ecuación diferencial estocástica clásica una arquitectura híbrida que incluye conmutación de régimen (regime-switching) y procesos de salto discreto:
En primer lugar, el estado del suministro eléctrico actúa como un interruptor determinado por un Proceso de Poisson No Homogéneo. Esta conmutación permite vincular el clima exterior con el riesgo operativo de la red; cuando el proceso de Poisson registra un apagón, el modelo simula un cambio radical en la termodinámica del equipo, pasando de un régimen de convección forzada (enfriamiento activo) a un régimen de convección natural pasiva (calentamiento inercial).
En segundo lugar, se incorpora un Proceso de Poisson Compuesto, denotado como \(J(t)\), para simular las perturbaciones del sistema ocasionadas por la operación clínica diaria. Matemáticamente, este término de salto se define como \(J(t) = \sum_{i=1}^{N_{puerta}(t)} \gamma_i\), donde la magnitud del choque térmico \(\gamma_i\) posee una distribución de probabilidad preestablecida (en este modelo, una distribución normal truncada) y \(N_{puerta}(t)\) es un Proceso de Poisson No Homogéneo, el cual es estrictamente independiente del movimiento browniano \(W(t)\). Cada apertura de la puerta del refrigerador por parte del personal médico inyecta estos choques instantáneos al equipo. Modelar estos eventos mediante la adición de \(dJ(t)\) a la ecuación diferencial permite simular cómo el calor externo entra de golpe y desplaza momentáneamente la temperatura interna antes de que la inercia térmica logre estabilizarla nuevamente.
3.9 Calibración de la Ecuación Diferencial Estocástica (EDE)
Para conectar el comportamiento aleatorio de la red eléctrica con la termodinámica del sistema, recordamos que la temperatura interna \(T(t)\) del equipo de refrigeración se describe mediante un proceso de Ornstein-Uhlenbeck, y que el modelo a usar incorpora perturbaciones por saltos y conmutación de régimen. Para facilitar la transición hacia la simulación computacional, desglosamos la ecuación diferencial en su forma operativa. Esta se presenta como un sistema de ecuaciones acopladas, cuyos términos de deriva y difusión dependen estrictamente del estado del suministro eléctrico \(\alpha(t)\):
\[ dT(t) = \begin{cases} \kappa_{on}(\theta_{set} - T(t))dt + \sigma_{on}dW(t) + dJ(t) & \text{si } \alpha(t) = 1 \text{ (Con energía)} \\ \kappa_{off}(T_{amb}(t) - T(t))dt + \sigma_{off}dW(t) & \text{si } \alpha(t) = 2 \text{ (Corte de luz)} \end{cases}, \quad 0 \leq t \leq T \tag{3.17}\]
Como se observa en el segundo escenario de la Ecuación 3.17, el término de saltos discretos (\(+ dJ(t)\)) se elimina por completo durante un corte de luz, esto debido a los protocolos operativos normativos: ante una falla eléctrica, el personal médico tiene estrictamente prohibido abrir el refrigerador. Al mantener la puerta sellada, se protege la inercia térmica provista por las botellas de agua y los paquetes refrigerantes, garantizando que la dinámica de la temperatura dependa exclusivamente de la pérdida pasiva hacia el ambiente exterior.
Cada término de la ecuación representa lo siguiente:
- \(T(t)\): La temperatura interna del equipo en el tiempo \(t\).
- \(\theta_{set}\): La temperatura ideal u objetivo configurada en el termostato.
- \(T_{amb}(t)\): La temperatura ambiente exterior de la ciudad de Tuxtla Gutiérrez.
- \(\kappa_{on}\): La velocidad de enfriamiento del equipo, la cual refleja la potencia del compresor cuando hay suministro eléctrico.
- \(\kappa_{off}\): La velocidad de pérdida térmica hacia el exterior durante un apagón, la cual depende de la calidad del aislamiento del equipo.
- \(\sigma_{on}\): La volatilidad térmica originada por la turbulencia y mezcla del aire interno cuando los ventiladores están encendidos.
- \(\sigma_{off}\): La volatilidad térmica presente durante un corte de luz, cuando el fluido interno transita a un estado de convección natural estática.
- \(dW(t)\): El incremento del movimiento browniano estándar, el cual captura las fluctuaciones aleatorias continuas.
- \(dJ(t)\): El incremento del Proceso de Poisson Compuesto que inyecta los choques térmicos al sistema.
- \(\gamma_i\): La magnitud aleatoria del choque térmico instantáneo provocado por la entrada de calor exterior en la \(i\)-ésima apertura del equipo.
- \(N_{puerta}(t)\): El Proceso de Poisson que modela el número de veces que el personal médico abre la puerta del refrigerador.
Para garantizar el rigor de la simulación, los coeficientes térmicos se dedujeron utilizando las especificaciones técnicas del refrigerador institucional Haier HBC-150. De acuerdo con el catálogo oficial PQS de la Organización Mundial de la Salud (OMS) (Organización Mundial de la Salud 2024, 89-90), este equipo cuenta con certificación de protección Grado A y está diseñado específicamente para mantener la red de frío en condiciones climáticas severas.
Con el propósito de aproximar de manera precisa el comportamiento de este equipo bajo los rigurosos estándares operativos de la OMS, cada coeficiente de la Ecuación 3.17 se calculó a partir de las leyes físicas de transferencia de calor y la mecánica de fluidos.
A continuación, se detalla el cálculo de cada parámetro de la Ecuación 3.17 basándose en las pruebas estandarizadas de este equipo.
3.9.1 Régimen 1: Operación Normal (\(\alpha(t) = 1\))
Este estado ocurre cuando hay suministro eléctrico. El compresor y los ventiladores del equipo funcionan activamente para llevar y mantener la temperatura interna hacia su nivel ideal.
Temperatura de Consigna (\(\theta_{set}\)): Se establece en \(5^\circ \text{C}\). Esta elección no es arbitraria; representa exactamente el punto medio del rango normativo de seguridad de la OMS (2°C a 8°C). Al centrar el termostato en 5°C, el modelo adquiere un margen simétrico de seguridad de \(3^\circ \text{C}\) hacia arriba (evitando la degradación por calor) y \(3^\circ \text{C}\) hacia abajo (evitando la destrucción de la vacuna por congelamiento).
Velocidad de Enfriamiento (\(\kappa_{on}\)): Este parámetro mide la rapidez con la que el motor abate el calor. Sus unidades son \(\text{hr}^{-1}\) (tasas de cambio por hora). Para calcularlo, recurrimos a la Ley de Enfriamiento de Newton (Vollmer 2009), la cual establece que la temperatura en el tiempo \(t\) sigue una trayectoria exponencial: \(T(t) = \theta_{set} + (T_0 - \theta_{set}) e^{-\kappa t}\). El protocolo de la OMS dicta un Pull-down time (tiempo de abatimiento) de 22 horas para que el equipo logre bajar desde un calor ambiental extremo de \(43^\circ \text{C}\) hasta entrar a la zona segura de \(8^\circ \text{C}\). Sustituyendo estos valores y nuestro objetivo de \(5^\circ \text{C}\), obtenemos: \[8 = 5 + (43 - 5) e^{-\kappa_{on} (22)}\] Al despejar la ecuación \(\left(\ln\left(\frac{3}{38}\right) = -22 \kappa_{on}\right)\), el resultado es \(\kappa_{on} \approx 0.1154 \text{ hr}^{-1}\).
Volatilidad y Ruido del Compresor (\(\sigma_{on}\)): El motor del refrigerador no enfría de manera perfecta; se enciende y se apaga constantemente, generando turbulencias. Dado que la solución analítica de un Proceso de Ornstein-Uhlenbeck sigue una distribución de probabilidad Normal (Gaussiana), podemos usar la “regla empírica” de la estadística para acotar este ruido. Para asegurar con un 99.7% de confianza que estas fluctuaciones naturales no rompan nuestro margen de seguridad de \(\pm 3^\circ \text{C}\), requerimos que tres desviaciones estándar abarquen exactamente esos tres grados (\(3 \times SD = 3^\circ \text{C}\)), lo que implica una desviación estándar de \(SD = 1^\circ \text{C}\).
La varianza de un proceso de Ornstein-Uhlenbeck clásico depende del tiempo \(t\). Sin embargo, dado que el refrigerador opera de manera continua durante bastante tiempo, la inercia inicial se desvanece y el sistema alcanza un estado estacionario (el límite cuando \(t \to \infty\)). Sustituyendo en la fórmula asintótica de la desviación estándar \(\left(SD = \frac{\sigma_{on}}{\sqrt{2\kappa_{on}}}\right)\): \[1 = \frac{\sigma_{on}}{\sqrt{2(0.1154)}} \implies \sigma_{on} = \sqrt{0.2308} \approx 0.4804\]
3.9.2 Régimen 2: Falla Eléctrica (\(\alpha(t) = 2\))
Este estado se activa cuando el Proceso de Poisson registra un apagón. El motor se detiene, y la dinámica del sistema pasa a depender exclusivamente del aislamiento térmico del refrigerador frente al clima exterior.
Temperatura Exterior (\(T_{amb}(t)\)): La temperatura ambiente actúa como el atractor del sistema. Dado que en Tuxtla Gutiérrez el calor oscila drásticamente a lo largo del día, esta variable se modela computacionalmente como una función senoidal que alcanza su pico máximo a las 15:00 hrs. Para nuestra calibración analítica del peor escenario (estiaje), evaluamos este atractor en un máximo absoluto de \(43^\circ \text{C}\). En este escenario, el término \((T_{amb} - T(t))\) se vuelve enorme y genera un calentamiento fuertemente acelerado.
Velocidad de Pérdida Térmica (\(\kappa_{off}\)): Representa la lentitud con la que el calor exterior logra penetrar el equipo. Se mide en \(\text{hr}^{-1}\). En el manual de la OMS (Organización Mundial de la Salud 2024, 89-90), el equipo Haier HBC-150 reporta un Hold-over time (tiempo de autonomía) de 60 horas con 50 minutos (\(60.833\) horas). Físicamente, esto significa que gracias a la enorme masa de hielo y agua en su interior, el equipo puede resistir rodeado de un calor de 43°C durante casi 61 horas antes de que su temperatura suba desde los \(5^\circ \text{C}\) ideales hasta el umbral de fallo normativo (\(10^\circ \text{C}\) en las pruebas de autonomía de la OMS). Usando nuevamente la Ley de Enfriamiento de Newton (Vollmer 2009) hacia el atractor externo: \[10 = 43 + (5 - 43) e^{-\kappa_{off} (60.833)}\] Al despejar la ecuación \(\left(\ln\left(\frac{33}{38}\right) = -60.833 \kappa_{off}\right)\), el resultado es \(\kappa_{off} \approx 0.00232 \text{ hr}^{-1}\).
La enorme diferencia de magnitud entre \(\kappa_{on}\) (rápido) y \(\kappa_{off}\) (muy lento) comprueba matemáticamente la alta efectividad térmica de los paquetes refrigerantes.
Volatilidad Externa (\(\sigma_{off}\)): Al cortarse la electricidad y apagarse los ventiladores, el aire interior deja de mezclarse a la fuerza (convección forzada) y pasa a un estado de reposo (convección natural estática). Según la mecánica de fluidos, esta ausencia de turbulencia reduce las variaciones aleatorias en un orden de magnitud. Por ello, la volatilidad térmica sin energía se calcula como la décima parte de la volatilidad activa: \(\sigma_{off} = \dfrac{\sigma_{on}}{10} \approx \frac{0.4804}{10} \approx 0.048\).
3.9.3 Perturbaciones Discretas: Apertura de Puertas (\(N_{puerta}(t)\))
De acuerdo con las normativas internacionales de manejo de la red de frío, la manipulación de los equipos de refrigeración debe limitarse estrictamente a los días hábiles (de lunes a viernes). Durante estas jornadas, el refrigerador debe abrirse idealmente un máximo de dos veces: una apertura rutinaria matutina (típicamente a las 8:00 hrs) para extraer los biológicos necesarios para el turno, y una apertura vespertina destinada a resguardar los sobrantes. Por el contrario, durante los fines de semana y días festivos, la unidad debe permanecer completamente sellada para maximizar la preservación de la inercia térmica, anulando por completo la inyección de calor exterior generada por el factor humano.
Dado que la apertura de las 8:00 hrs es un evento operativo determinista y predecible, la aleatoriedad del sistema recae en el resto de la jornada, condicionada por la variabilidad de la demanda clínica. Para aproximar el modelo matemático a la realidad, aislar el factor humano y evaluar sistemáticamente la resiliencia termodinámica del equipo frente a fallas eléctricas, el Proceso de Poisson se segmenta en tres escenarios de simulación. Cada escenario persigue un objetivo analítico particular, definiendo tasas de intensidad (\(\lambda\)) dependientes del tiempo:
- Escenario A (Control Absoluto): Diseñado como línea base para evaluar la eficiencia térmica pura del equipo, respondiendo a la pregunta de cuánto tiempo el refrigerador es capaz de mantener la cadena de frío sin intervención humana. Simula días de inactividad clínica, como fines de semana o días festivos, con lo cual la intensidad del proceso de Poisson es nula en todo momento (\(\lambda = 0\)). Las variaciones de la temperatura interna dependen exclusivamente de los cortes eléctricos y la temperatura ambiente exterior.
- Escenario B (Apego Normativo): Su objetivo es establecer el margen de seguridad del equipo operando bajo un entorno clínico ideal. Respeta la única apertura extra permitida por la norma para resguardar biológicos sobrantes. Esta perturbación se modela como un Proceso de Poisson con una intensidad de \(1\) evento/día. Operativamente, esta tarea se realiza al final de la jornada de consulta. Al distribuir esta probabilidad en un horario vespertino de 2 horas (de 15:00 a 17:00 hrs), la intensidad local se calcula como: \[\lambda_{tarde} = \frac{1 \text{ evento}}{2 \text{ horas}} = 0.5 \text{ eventos/hora}\] Fuera de este horario, la intensidad del proceso es estrictamente cero.
- Escenario C (Estrés Operativo y Contingencia): Busca cuantificar la vulnerabilidad del sistema ante situaciones de crisis, analizando cómo el exceso de aperturas agota el “colchón térmico” de la unidad dejándola indefensa si ocurre un apagón. Modela la saturación clínica donde la demanda sobrepasa la capacidad de los termos auxiliares de las enfermeras, obligando a realizar “viajes de resurtido” al refrigerador principal. Con base en los cercos epidemiológicos registrados en Chiapas durante 2025, el volumen de aplicación superó las 415,000 dosis mensuales, triplicando el consumo regular (aprox. 130,000 dosis). Al escalar operativamente este incremento, la necesidad de extracción matutina se triplica (1 apertura fija a las 8:00 hrs + 2 aperturas de resurtido). El retorno vespertino se mantiene normativo (1 apertura si es que sobran vacunas para resguardar nuevamente). Esto divide el proceso en dos ventanas activas:
- Ventana Matutina (Resurtido): Se esperan 2 eventos distribuidos de 09:00 a 14:00 hrs. \[\lambda_{mañana} = \frac{2 \text{ eventos}}{5 \text{ horas}} = 0.4 \text{ eventos/hora}\]
- Ventana Vespertina (Sobrantes): Se mantiene en \(0.5 \text{ eventos/hora}\) de 15:00 a 17:00 hrs.
3.9.4 Impacto de la Apertura y Salto Térmico (\(\gamma_i\))
Al tratarse de un equipo tipo cofre (Icelined Refrigerator), la pérdida térmica por convección es físicamente menor que en equipos verticales, ya que el aire frío es más denso y tiende a permanecer en el fondo. Sin embargo, el ingreso de aire caliente es inevitable y su magnitud no es una constante; depende del tiempo que la puerta permanezca abierta y de la correcta disposición del balasto térmico (botellas de agua).
Para modelar rigurosamente este choque térmico, nos basamos en los estudios experimentales del National Institute of Standards and Technology (NIST) (Chojnacky et al. 2010). En sus pruebas de uso clínico, el NIST demostró que:
- En el mejor escenario (aperturas rápidas con botellas de agua como barrera), el incremento de temperatura en la vacuna es mínimo, limitándose a \(0.16^\circ\text{C}\).
- En un escenario de descuido o falta de balasto térmico, el ingreso de calor genera un incremento adicional de \(1.2^\circ\text{C}\), resultando en un choque máximo esperado de \(1.36^\circ\text{C}\) (\(0.16 + 1.2\)).
Dado que en la práctica clínica intervienen múltiples factores humanos aleatorios e independientes, el Teorema del Límite Central justifica modelar la magnitud del salto \(\gamma_i\) como una variable aleatoria con distribución Normal (\(\mathcal{N}(\mu, \sigma^2)\)).
La media esperada (\(\mu\)) se calcula como el punto medio entre los límites operativos del NIST: \[\mu = \frac{0.16 + 1.36}{2} = 0.76^\circ\text{C}\]
Para determinar la varianza, aplicamos la regla empírica (regla de las tres sigmas), asumiendo que el 99.7% de los eventos clínicos caerán dentro de los límites físicos observados. Igualando la distancia desde la media hasta el extremo superior con \(3\sigma\): \[3\sigma = 1.36 - 0.76 = 0.60 \implies \sigma = 0.20^\circ\text{C}\]
Restricción Termodinámica (Filtro de Truncamiento): Las colas de una distribución Normal estándar se extienden hasta menos infinito. En el contexto físico, un valor de \(\gamma_i < 0\) implicaría que al abrir la puerta el refrigerador absorbe frío del exterior cálido, lo cual viola la Segunda Ley de la Termodinámica. Asimismo, el estudio del NIST establece que el ingreso mínimo absoluto de calor es de \(0.16^\circ\text{C}\).
Para garantizar la viabilidad física del proceso estocástico durante las simulaciones de Montecarlo, la variable aleatoria \(\gamma_i\) se define como una distribución truncada inferiormente, aplicando un filtro de maximización: \[\gamma_i = \max(0.16, X) \quad \text{donde } X \sim \mathcal{N}(0.76, 0.20^2) \tag{3.18}\] De esta forma, la Ecuación Diferencial Estocástica absorbe la variabilidad operativa del personal médico sin riesgo de arrojar perturbaciones termodinámicamente imposibles.
3.10 Cinética de Degradación Térmica (Ecuación de Arrhenius)
La Ecuación Diferencial Estocástica calibrada anteriormente permite simular con alta precisión la trayectoria térmica interna del refrigerador, \(T(t)\). Sin embargo, para evaluar el impacto logístico y en la salud pública, es necesario traducir estas variaciones térmicas en pérdida de viabilidad del biológico.
La degradación de las vacunas (particularmente las vivas atenuadas, como el biológico contra el Sarampión o la Poliomielitis) es un proceso de desnaturalización proteica que se acelera exponencialmente con el incremento de temperatura. Desde una perspectiva fisicoquímica, la concentración de ingrediente activo decae siguiendo una cinética de primer orden, cuya velocidad de reacción \(k(T)\) se modela mediante la Ecuación de Arrhenius (Sinko 2011):
\[ k(T(t)) = A \exp\left(-\frac{E_a}{R \cdot T_K(t)}\right) \tag{3.19}\]
Cada parámetro de la ecuación representa lo siguiente:
- \(k(T(t))\): La constante de velocidad de degradación térmica en el instante \(t\), expresada en \(\text{s}^{-1}\).
- \(T_K(t)\): La temperatura interna del fluido convertida a escala absoluta (grados Kelvin), calculada como \(T(t) + 273.15\).
- \(R\): La Constante Universal de los Gases, con un valor fijo de \(8.314 \, \text{J/(mol} \cdot \text{K)}\).
- \(E_a\): La Energía de Activación, calibrada en \(85,000 \, \text{J/mol}\) (\(85 \text{ kJ/mol}\)). Representa la barrera energética mínima que las moléculas del antígeno deben superar para iniciar el proceso irreversible de ruptura estructural.
- \(A\): El Factor Pre-exponencial (o factor de frecuencia), ajustado a un valor empírico de \(1.2 \times 10^{13} \, \text{s}^{-1}\). Este parámetro asegura que la tasa de degradación sea prácticamente nula cuando el equipo opera correctamente dentro de la red de frío (2°C a 8°C), pero que escale de forma logarítmica severa ante una excursión térmica crítica.
3.10.1 Modelo de Daño Térmico Acumulado y Merma de Inventario
En la práctica clínica normativa, las variaciones térmicas que oscilan dentro del rango seguro (2°C a 8°C) se consideran tolerables, y su impacto en la caducidad a corto plazo es ignorado. El daño logístico ocurre exclusivamente cuando la red de frío se rompe. Para incorporar este protocolo operativo al modelo matemático, la degradación se condiciona mediante una función indicadora \(\mathbb{I}_{\{T(t) > 8^\circ\text{C}\}}\), la cual activa el proceso de pérdida térmica únicamente cuando la temperatura cruza el límite superior de seguridad.
Bajo este esquema, la cantidad de vacunas viables en el inventario, \(I(t)\), se rige por una Ecuación Diferencial Ordinaria (EDO) de primer orden, donde la tasa de pérdida es proporcional al inventario actual (Sinko 2011):
\[ \frac{dI(t)}{dt} = -k(T(t)) \cdot I(t) \cdot \mathbb{I}_{\{T(t) > 8\}} \]
Al resolver esta ecuación diferencial mediante variables separables e integrando desde el inicio de la simulación (\(t=0\)) hasta un tiempo \(t\), se obtiene la cantidad exacta de inventario sobreviviente:
\[ \int_{I(0)}^{I(t)} \frac{dI}{I} = - \int_{0}^{t} k(T(s)) \cdot \mathbb{I}_{\{T(s) > 8\}} \, ds \implies I(t) = I(0) \exp\left( - \int_{0}^{t} k(T(s)) \cdot \mathbb{I}_{\{T(s) > 8\}} \, ds \right) \]
El término integral dentro de la función exponencial define formalmente el Daño Acumulado (\(DA\)):
\[ DA(t) = \int_{0}^{t} k(T(s)) \cdot \mathbb{I}_{\{T(s) > 8\}} \, ds \]
Esta formulación integral permite contabilizar analíticamente el estrés térmico exacto de cada trayectoria. Computacionalmente, para actualizar el estado del sistema en pasos discretos de tiempo \(\Delta t\) dentro de la simulación de Montecarlo, la solución analítica de la EDO estipula que la fracción de inventario perdido (merma instantánea) durante ese diferencial de tiempo se calcula extrayendo el complemento del factor de supervivencia:
\[ \text{Merma}(t) = I(t) \left( 1 - e^{-k(T(t)) \Delta t} \right) \]
Esta discretización asegura que el modelo numérico refleje con precisión la cinética fisicoquímica. En cada iteración, el algoritmo descuenta las dosis destruidas por la exposición al calor durante los apagones, permitiendo la activación automatizada de los protocolos logísticos de reabastecimiento si el inventario cae por debajo de los umbrales críticos de operación.
3.10.2 Exclusión del Riesgo por Congelamiento (\(T < 2^\circ\text{C}\))
Es importante destacar que el modelo de daño térmico propuesto se restringe exclusivamente a cuando la temperatura interna del equipo es mayor que 8°C (\(T(t) > 8^\circ\text{C}\)). Aunque el descenso de la temperatura por debajo de los \(2^\circ\text{C}\) (y especialmente hacia el punto de congelación a \(0^\circ\text{C}\)) representa un riesgo logístico crítico que destruye la viabilidad del biológico debido a la cristalización del agua, este fenómeno se excluye de la simulación por tres motivos teóricos y prácticos fundamentales:
- Incompatibilidad Cinética: El daño por congelamiento obedece a un cambio de fase físico y estructural, no a una reacción química progresiva. Por lo tanto, la Ecuación de Arrhenius carece de validez termodinámica para modelar este deterioro continuo, imposibilitando su incorporación bajo la misma Ecuación Diferencial Ordinaria.
- Deriva Térmica del Entorno: La estructura matemática de la Ecuación Diferencial Estocástica formulada no posee fuerzas de deriva que empujen el sistema hacia el sobreenfriamiento. Con suministro eléctrico, el equipo gravita hacia el termostato (\(\theta_{set} = 5^\circ\text{C}\)). Durante un apagón, la temperatura interna tiende hacia el ambiente cálido exterior (\(T_{amb}(t) \ge 5^\circ\text{C}\)), mientras que las perturbaciones por apertura de puertas (\(dJ(t)\)) únicamente inyectan energía térmica positiva.
- Arquitectura del Equipo: La simulación asume un equipo tipo cofre precalificado (Icelined Refrigerator). El diseño de estas unidades incorpora un revestimiento estabilizador (balasto de agua/hielo) cuya función termodinámica es precisamente absorber los picos de frío directo del compresor, protegiendo al biológico de caídas accidentales de temperatura.
Por consiguiente, bajo las condiciones operativas y climáticas simuladas, la probabilidad de sobreenfriamiento se considera nula. El análisis de vulnerabilidad del inventario queda acotado al estrés por calor derivado de la interrupción del suministro eléctrico.
3.11 Dinámica de Inventario y Demanda Clínica
Para evaluar la capacidad de la red de frío, el modelo termodinámico debe interactuar con las variaciones operativas del inventario. La simulación incorpora parámetros logísticos basados en la teoría de gestión de inventarios y en el calendario de salud pública mexicano.
3.11.1 Capacidad Logística y Reglas de Reabastecimiento
El modelo asume el uso de un equipo precalificado tipo cofre (Icelined Refrigerator), específicamente el modelo Haier HBC-150. Con un volumen útil de almacenamiento de aproximadamente 122 litros, la capacidad máxima del sistema se establece en \(5,000\) dosis. Esta cifra contempla el volumen físico de los empaques secundarios y el espacio necesario para garantizar la correcta convección del aire interno.
Para gestionar el riesgo de desabasto provocado por la merma térmica, se define un Punto de Reorden de \(1,000\) dosis. Este valor representa un stock de seguridad estratégico del 20%, el cual garantiza un colchón operativo para que la clínica continúe sus labores de inmunización mientras se activa el protocolo de reabastecimiento jurisdiccional.
3.11.2 Comportamiento Estacional de la Demanda
Ante la variabilidad de las aplicaciones diarias, el consumo se modela mediante una demanda sintética representativa, calibrada para reflejar la afluencia de un Centro de Salud Urbano de primer nivel. Esta demanda base no es constante, sino que obedece al comportamiento estacional del sistema de salud en México, tal como se detalla en la siguiente tabla:
| Mes | Demanda Base Diaria (Dosis) | Justificación Epidemiológica |
|---|---|---|
| Enero | 65 | Consumo regular de esquema básico. |
| Febrero | 110 | 1ª Jornada Nacional de Salud Pública. |
| Marzo | 87 | Regulación post-campaña. |
| Abril | 65 | Consumo regular de esquema básico. |
| Mayo | 110 | 2ª Jornada Nacional de Salud Pública. |
| Junio - Agosto | 65 | Periodo de consumo valle (vacaciones/lluvias). |
| Septiembre | 76 | Transición a temporada otoñal. |
| Octubre | 163 | Inicio de campaña invernal (Influenza/Neumococo). |
| Noviembre | 130 | Continuación de campaña invernal. |
| Diciembre | 87 | Desaceleración por periodo vacacional de invierno. |
3.11.3 Acoplamiento Estocástico Termodinámico-Logístico
La innovación principal de este bloque logístico radica en la dependencia entre el estrés térmico y el estrés de inventario. La demanda diaria real en el modelo, \(D(t)\), no es un escalar estático, sino una variable aleatoria que se acopla dinámicamente a las perturbaciones del proceso de Poisson \(N_{puerta}(t)\) calibrado en la sección anterior.
La regla de consumo se define a partir de los tres escenarios de simulación:
Escenario A (Control Absoluto): Al simular días de inactividad clínica (fines de semana o días festivos), el personal no abre el equipo (\(\lambda = 0\)). Consecuentemente, no hay extracción de biológicos: \[D_A(t) = 0\]
Escenario B (Apego Normativo): Simula una jornada clínica estándar con una única apertura matutina de extracción. El consumo diario es exactamente igual a la demanda base dictada por el mes correspondiente: \[D_B(t) = D_{base}\]
Escenario C (Estrés Operativo): Durante contingencias o cercos epidemiológicos, la demanda sobrepasa la capacidad de los termos auxiliares, forzando viajes de resurtido al equipo principal. Dado que el proceso general \(N_{puerta}(t)\) engloba todas las aperturas del día, definimos a \(N_{resurtido}(t)\) como el subconjunto de eventos ocurridos exclusivamente durante la ventana de máxima afluencia matutina (09:00 a 14:00 hrs).
Matemáticamente, por cada evento de resurtido generado, se extrae un bloque adicional de vacunas equivalente a la demanda base. Así, la demanda diaria total se convierte en una variable aleatoria dependiente del proceso de Poisson, definida por: \[D_C(t) = D_{base} + \left( N_{resurtido}(t) \times D_{base} \right)\]
Este acoplamiento garantiza que las trayectorias de Montecarlo con mayor cantidad de aperturas de puerta y por ende, mayor ingreso de calor y daño térmico correspondan exactamente a los días con mayor rotación de inventario, reflejando con precisión la saturación logística de una crisis sanitaria.