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 evitara 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 vinculara 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 algun 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 contien 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 selecciono 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. En la modelación aplicada, utilizar directamente este valor para un solo refrigerador implicaría asumir que el edificio clínico 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 las fallas de transmisión se distribuyen de manera uniforme sobre estas infraestructuras mediante el principio de razón insuficiente de Laplace, la probabilidad \(p\) de que un evento 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 absorbe la clínica local se comporta como una variable aleatoria binomial \(N_{local} \sim \mathcal{B}(N_{total}, p)\), cuya esperanza matemática es:
\[ \mathbb{E}[N_{local}] = \frac{516 \text{ cortes}}{8 \text{ subestaciones}} = 64.5 \approx 64 \text{ cortes al año} \tag{3.2}\]
Al ajustar al entero más cercano para mantener la naturaleza discreta de los eventos, 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 se presentan exactamente 5.33 cortes idénticos cada mes). La variabilidad mensual depende directamente de los factores climáticos estacionales de la región. Para no asignar los apagones 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 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 aislar 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 mantener la conservación de la suma total en 64), obtenemos la distribución estacional empírica de fallas 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 estas 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 instantánea de eventos \(\lambda(t)\) varía en función del tiempo para absorber la estacionalidad climática.
El modelo se construye mediante la suma de una tasa base constante y dos funciones de perturbación gaussianas no normalizadas. La elección de núcleos gaussianos obedece a dos razones analíticas fundamentales:
- Transición Teórica Suave: Los fenómenos climáticos (como una ola de calor o la época de huracanes) no inician ni terminan de forma abrupta. Una gaussiana permite modelar la transición gradual de las estaciones sin introducir discontinuidades de salto que volverían inestable la resolución numérica de la EDE.
- Correlación Física Estacional: Permite alinear y concentrar el mayor riesgo de que la EDE salte al estado de falla eléctrica (apagón) exactamente en los días donde el clima ambiental exterior es más adverso.
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)\]
Para calibrar estos coeficientes con los datos empíricos de Tuxtla Gutiérrez (2024), se realiza el siguiente desglose analítico:
Tasa Base (\(\lambda_{base} = 0.0667\)): Representa el riesgo mínimo de fondo debido a fallas fortuitas o mantenimiento de la red, independiente del clima. Se calcula desde los meses más estables (febrero y noviembre, con 2 cortes al mes). Al dividir este valor, obtenemos la tasa diaria: \[\lambda_{base} = \frac{2 \text{ cortes}}{30 \text{ días}} \approx 0.0667 \text{ cortes/día}\]
Medias Temporales (\(\mu_1 = 152\) y \(\mu_2 = 244\)): Ubican el centro máximo de los picos de riesgo. El primer pico se centra a principios de junio (el 1 de junio equivale al día 152 del año), alineándose con el inicio del verano más caluroso. El segundo se centra a inicios de septiembre (día 244), el momento históricamente más activo de la temporada de huracanes.
Amplitudes (\(A_1 = 0.2333\) y \(A_2 = 0.2333\)): Capturan la altura máxima de la campana en los días críticos. Como junio y septiembre registraron 9 cortes mensuales, la tasa en su pico es de \(9/30 = 0.30\) cortes/día. Dado que la tasa base ya aporta \(0.0667\), la magnitud neta que debe sumar la gaussiana es la diferencia: \(0.30 - 0.0667 = 0.2333\).
Varianzas Estacionales (\(\sigma_1 = 40\) y \(\sigma_2 = 30\)): Controlan el ancho temporal de las afectaciones. Para la ola de calor, se asigna una desviación estándar de \(40\) días (varianza \(2\sigma_1^2 = 3200\)) para garantizar que la campana sea lo suficientemente ancha como para cubrir desde finales de abril (cuando inician las altas temperaturas) hasta mediados de julio. Para las lluvias, una desviación estándar de \(30\) días (varianza \(2\sigma_2^2 = 1800\)) modela una perturbación que cubre por completo agosto y septiembre, extendiéndose hacia octubre.
Sustituyendo cada coeficiente calibrado, la función de intensidad operativa queda definida formalmente como:
\[\lambda(t) = 0.0667 + 0.2333 \exp\left(-\frac{(t - 152)^2}{3200}\right) + 0.2333 \exp\left(-\frac{(t - 244)^2}{1800}\right)\]
3.6 Validación del Modelo (Cálculo Integral)
Por los fundamentos de los Procesos de Poisson No Homogéneos, el número esperado de eventos 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\]
Para validar la consistencia analítica de nuestra función \(\lambda(t)\), particionamos el año de 365 días en 12 subintervalos (correspondientes a los días acumulados de cada mes) y resolvemos numéricamente sus respectivas integrales, contrastando 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 Integral \(\mathbb{E}[N]\) |
|---|---|---|---|
| Enero | \([0, 31]\) | 7 | 2.10 |
| Febrero | \([31, 59]\) | 2 | 2.06 |
| Marzo | \([59, 90]\) | 3 | 3.26 |
| Abril | \([90, 120]\) | 3 | 5.56 |
| Mayo | \([120, 151]\) | 6 | 8.56 |
| Junio | \([151, 181]\) | 9 | 8.75 |
| Julio | \([181, 212]\) | 7 | 8.16 |
| Agosto | \([212, 243]\) | 8 | 9.40 |
| Septiembre | \([243, 273]\) | 9 | 8.29 |
| Octubre | \([273, 304]\) | 4 | 4.63 |
| Noviembre | \([304, 334]\) | 2 | 2.38 |
| Diciembre | \([334, 365]\) | 4 | 2.10 |
La suma acumulada de las 12 integrales definidas a lo largo de todo el dominio anual arrojó un valor total de 65.25 eventos esperados. Al comparar este resultado con la meta de 64 cortes empíricos, se observa que el modelo añade, en promedio, una falla anual extra.
Esta desviación hacia arriba es una decisión metodológica deliberada para blindar la simulación de Montecarlo: al ensanchar las campanas de riesgo para cubrir desde abril hasta octubre de manera ininterrumpida, el modelo genera un escenario ligeramente más conservador y severo que la simple estadística aislada. Esto garantiza que las proyecciones de estrés térmico sobre la red de frío no subestimen el peligro real al que estarán sometidos los biológicos de la OMS durante las estaciones críticas en Tuxtla Gutiérrez.
3.7 Distribución Log-Normal de las Duraciones de las Fallas
Una vez modelada la frecuencia con la que ocurren las interrupciones eléctricas, el siguiente componente crítico para el simulador termodinámico consiste en determinar cuánto tiempo dura cada apagón (\(D\)). La duración de una falla es una variable aleatoria continua estrictamente positiva con una marcada asimetría hacia la derecha (la mayoría de los cortes se resuelven en pocas horas, pero existen eventos atípicos muy prolongados). Por este motivo, la teoría analítica establece que el comportamiento de las duraciones se modela mediante una distribución Log-normal:
\[ D \sim \text{LogNormal}(\mu, \sigma^2) \]
Para calibrar los parámetros de escala \(\mu\) (media de los logaritmos) y de forma \(\sigma\) (desviación estándar de los logaritmos) para cada régimen estacional, es necesario procesar los indicadores empíricos de la CFE bajo un esquema de desagregación temporal.
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 general de una interrupción en Tuxtla Gutiérrez fue de \(2.39\) horas por corte Tabla 3.1. Sabiendo que el factor de escala espacial redujo el riesgo del edificio a una esperanza de \(64\) cortes locales al año, podemos deducir la cantidad total de tiempo que el refrigerador experimentará desabasto eléctrico a lo largo de todo el año:
\[ \text{Horas Totales del Año} = 64 \text{ cortes/año} \times 2.39 \text{ horas/corte} = 152.96 \text{ horas/año} \]
Para distribuir estas \(152.96\) horas anuales entre los 12 meses respetando la estacionalidad de la red, utilizamos el indicador SAIDI acumulado del Informe Anual de la CFE 2024 Tabla 3.2. Primero, calculamos la diferencia marginal (\(\Delta \text{SAIDI}_i\)) para aislar el impacto de cada mes. Posteriormente, mediante una regla de tres simple, equiparamos el SAIDI total acumulado al final del año (\(10.812\) horas) con nuestras \(152.96\) horas de afectación local. La asignación de horas de falla para cada mes viene dada por:
\[ \text{Horas del Mes}_i = \left( \frac{\Delta \text{SAIDI}_i}{10.812} \right) \times 152.96 \]
Finalmente, para conocer la duración promedio esperada de un corte individual en un mes específico (\(E_i\)), se divide la cantidad de horas asignadas de ese mes entre el número de cortes calculados previamente Tabla 3.3 (\(E_i = \text{Horas del Mes}_i / \text{Cortes}_i\)):
| Mes | SAIDI Acumulado | \(\Delta\) SAIDI | Horas del Mes Calculadas | Cortes del Mes | Duración Media \(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 Régimen
Para alimentar la conmutación de regímenes de la EDE, agrupamos los meses en las tres temporadas climáticas analizadas. Sin embargo, realizar un promedio aritmético simple entre las duraciones mensuales introduciría un sesgo metodológico, ya que trataría por igual a un mes con 9 cortes que a uno con solo 2 cortes.
Para resolver esto, calculamos la media empírica (\(\mathbb{E}[D_{reg}]\)) y la varianza empírica (\(\mathbb{V}[D_{reg}]\)) de cada régimen utilizando estadística descriptiva ponderada, donde el peso de cada mes está determinado estrictamente por su número de cortes (\(\omega_i = \text{Cortes}_i\)). Analíticamente, es decir:
\[ \mathbb{E}[D_{reg}] = \frac{\sum \omega_i E_i}{\sum \omega_i} = \frac{\text{Total de Horas del Régimen}}{\text{Total de Cortes del Régimen}} \]
\[ \mathbb{V}[D_{reg}] = \frac{\sum \omega_i (E_i - \mathbb{E}[D_{reg}])^2}{\sum \omega_i} \]
Agrupando los datos bajo este criterio de ponderación, se obtienen las siguientes métricas muestrales por temporada:
Temporada de Calor (Mayo, Junio y Julio): Acumula \(22\) cortes y \(53.25\) horas de interrupción, resultando en una duración media de \(\mathbb{E}[D] = 2.420\) horas y una varianza de \(\mathbb{V}[D] = 0.373\).
Temporada de Lluvias (Agosto, Septiembre y Octubre): Suma \(21\) cortes y \(35.45\) horas de interrupción, generando una duración media de \(\mathbb{E}[D] = 1.688\) horas y una varianza de \(\mathbb{V}[D] = 0.343\).
Temporada Base e Invierno (Resto de meses): Concentra \(23\) cortes y \(64.25\) horas, arrojando una duración media de \(\mathbb{E}[D] = 2.793\) horas con una varianza más dispersa de \(\mathbb{V}[D] = 4.650\) debido a la anomalía de abril.
3.7.3 Transformación Analítica a Parámetros Log-Normales
Una vez que conocemos la media (\(\mathbb{E}[D]\)) y la varianza (\(\mathbb{V}[D]\)) del sistema, debemos mapearlas hacia los parámetros analíticos internos \(\mu\) y \(\sigma\) que requiere la función de distribución de probabilidad Log-normal.
La esperanza y la varianza de una variable aleatoria con distribución Log-normal están dadas por las siguientes dos identidades:
\[ \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) = (\mathbb{E}[D])^2 \left[\exp(\sigma^2) - 1\right] \]
Al resolver este sistema de ecuaciones no lineales de forma algebraica para despejar los parámetros objetivos en función de los momentos estadísticos conocidos, se deducen las fórmulas generales de calibración:
\[ \sigma = \sqrt{\ln\left(1 + \frac{\mathbb{V}[D]}{(\mathbb{E}[D])^2}\right)} \]
\[ \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) \]
Al sustituir las medias y varianzas de cada temporada en estas relaciones, 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 | Meses que lo integran | Media Empírica \(\mathbb{E}[D]\) | Varianza Empírica \(\mathbb{V}[D]\) | Parámetro \(\mu\) | Parámetro \(\sigma\) |
|---|---|---|---|---|---|
| Base (Invierno) | Ene, Feb, Mar, Abr, Nov, Dic | 2.793 | 4.650 | 0.7933 | 0.6838 |
| Calor (Ola de Verano) | May, Jun, Jul | 2.420 | 0.373 | 0.8529 | 0.2485 |
| Lluvias (Ciclones) | Ago, Sep, Oct | 1.688 | 0.343 | 0.4667 | 0.3371 |
3.8 Calibración de la Ecuación Diferencial Estocástica (EDE)
Enlazando la dinámica de la red eléctrica con la termodinámica, recordamos que la temperatura interna \(T(t)\) del refrigerador está gobernada por un proceso de Ornstein-Uhlenbeck con saltos y conmutación de régimen. Para que la transición hacia el modelo computacional sea transparente, desglosamos la ecuación abstracta en su forma operativa mediante un sistema de ecuaciones acopladas, el cual depende estrictamente del estado del suministro eléctrico \(\alpha(t)\):
\[ dT(t) = \begin{cases} \kappa_{on}(\theta_{set} - T(t))dt + \sigma_{on}dW(t) + \gamma dN_{puerta}(t) & \text{si } \alpha(t) = 1 \text{ (Con energía)} \\ \kappa_{off}(T_{amb}(t) - T(t))dt + \sigma_{off}dW(t) + \gamma dN_{puerta}(t) & \text{si } \alpha(t) = 2 \text{ (Corte de luz)} \end{cases} \tag{3.8}\]
Para la modelación estocástica, los parámetros de inercia térmica y tiempo de autonomía (holdover time) se calibraron utilizando los datos técnicos del refrigerador institucional Haier HBC-150, el cual cuenta con certificación de protección Grado A contra congelamiento en el catálogo oficial PQS de la OMS (Organización Mundial de la Salud 2024).
Para reflejar el comportamiento de este equipo bajo los rigurosos estándares normativos, calibramos cada coeficiente de la Ecuación 3.8 con base en la física del enfriamiento:
3.8.1 Régimen 1: Operación Normal (\(\alpha(t) = 1\))
Ocurre cuando hay energía eléctrica y el compresor funciona empujando la temperatura hacia su objetivo.
- Temperatura de Consigna (\(\theta_{set}\)): Se establece en \(4^\circ \text{C}\), el centro ideal del rango normativo de seguridad (2°C a 8°C).
- Velocidad de Enfriamiento (\(\kappa_{on}\)): Representa la capacidad de potencia del motor para recuperar la temperatura. Se calibra en \(\kappa_{on} = 0.1188\), una velocidad rápida propia de la compresión activa de acuerdo con (Organización Mundial de la Salud 2024).
- Ruido del Compresor (\(\sigma_{on}\)): El motor se prende y se apaga constantemente para mantener los 4°C. Acorde con un intervalo de confianza del 99.7% (aplicando la regla de las 3 sigmas de la distribución normal), queremos que las fluctuaciones máximas no excedan los \(\pm 2^\circ \text{C}\) para evitar llegar a los 2°C (riesgo de congelación) o a los 6°C. Sabiendo que \(3 \times SD = 2^\circ \text{C}\), la desviación estándar estacionaria debe ser \(SD \approx 0.667^\circ \text{C}\). Despejando la volatilidad de la varianza a largo plazo del proceso Ornstein-Uhlenbeck \(\left(SD = \frac{\sigma_{on}}{\sqrt{2\kappa_{on}}}\right)\): \[0.667 = \frac{\sigma_{on}}{\sqrt{2(0.1188)}} \implies \sigma_{on} \approx 0.325\]
3.8.2 Régimen 2: Falla Eléctrica (\(\alpha(t) = 2\))
Se activa cuando el Proceso de Poisson dice que hay un apagón y el enfriamiento activo se detiene.
- 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}\)): De acuerdo con las especificaciones del fabricante validadas por la OMS (Organización Mundial de la Salud 2024), el equipo Haier HBC-150 garantiza un Holdover Time (tiempo de autonomía) de 60 horas con 50 minutos (60.83 horas) a una temperatura ambiente extrema de 43°C antes de rebasar los 8°C. Usando la solución analítica determinista de la Ley de Enfriamiento de Newton desde los 4°C de consigna hasta el límite de 8°C: \[8 = 43 + (4 - 43) e^{-\kappa_{off} (60.83)}\] Al resolver la ecuación \(\left(\ln\left(\frac{35}{39}\right) = -60.83 \kappa_{off}\right)\), obtenemos un valor sumamente bajo de \(\kappa_{off} \approx 0.00178 \text{ hr}^{-1}\). Esta dramática diferencia entre un parámetro rápido (\(\kappa_{on}\)) y uno muy lento (\(\kappa_{off}\)) demuestra la efectividad de los revestimientos de hielo y justifica el uso de EDEs con conmutación de régimen.
- Volatilidad Externa (\(\sigma_{off}\)): Sin el compresor encendido, las variaciones aleatorias obedecen únicamente a micro-filtraciones de aire por las gomas de la puerta. Por ello, establecemos un ruido significativamente menor: \(\sigma_{off} = 0.05\).
3.8.3 Perturbaciones Discretas: Apertura de Puertas (\(N_{puerta}(t)\))
De acuerdo con las normativas internacionales de manejo de red de frío, el refrigerador debe abrirse idealmente un máximo de dos veces al día: una apertura rutinaria matutina para extraer los frascos necesarios para la jornada y una vespertina para resguardar los sobrantes.
- Tasa de Aperturas (\(\lambda_{puerta}\)): Dado que la apertura matutina es un evento operativo fijo, la verdadera estocasticidad recae en la tarde, ya que el ingreso de biológicos depende de si hubo sobrantes en los termos clínicos. Por ello, la perturbación aleatoria se restringe a este momento y se modela computacionalmente con un Proceso de Poisson Homogéneo con una tasa de \(\lambda = 1\) evento/día.
- Impacto de la Apertura (\(\gamma\)): Al tratarse de un equipo tipo cofre (Icelined Refrigerator), la pérdida térmica es mucho menor que en equipos verticales, ya que el aire frío, siendo más denso, tiende a permanecer en el fondo del gabinete. Sin embargo, basándonos en estudios de convección para una apertura prolongada típica de 5 minutos (tiempo estimado para organizar y contabilizar viales), se introduce una carga térmica inesperada. Asignamos un incremento de \(\gamma = +1^\circ \text{C}\) por evento. Esta es una perturbación instantánea que rompe la continuidad de la trayectoria y reduce críticamente el “colchón térmico” disponible en caso de ocurrir un apagón durante la noche.
3.9 Cinética de Degradación Térmica (Ecuación de Arrhenius)
La Ecuación Diferencial Estocástica calibrada anteriormente nos permite simular con alta precisión la temperatura interna del refrigerador \(T(t)\). Sin embargo, para evaluar el impacto 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 calor. Esta velocidad de degradación \(k(T)\) se modela mediante la Ecuación de Arrhenius (Chen y Kristensen 2009):
\[ k(T(t)) = A \exp\left(-\frac{E_a}{R \cdot T_K(t)}\right) \tag{3.9}\]
Donde la temperatura \(T_K(t) = T(t) + 273.15\) se expresa en grados Kelvin y los parámetros cinéticos para un biológico termosensible genérico se establecen como:
- Constante Universal de los Gases (\(R\)): \(8.314 \, \text{J/(mol} \cdot \text{K)}\).
- Energía de Activación (\(E_a\)): Se calibra en \(85,000 \, \text{J/mol}\) (\(85 \text{ kJ/mol}\)). Este valor representa la barrera energética que debe superarse para que las proteínas del antígeno comiencen a romperse irreversiblemente.
- Factor Pre-exponencial (\(A\)): Se ajusta a un valor de \(1.2 \times 10^{14}\). Este factor de frecuencia asegura que la tasa de degradación sea prácticamente nula cuando el equipo opera correctamente en el rango normativo de 2°C a 8°C, pero escale logarítmicamente durante una excursión térmica severa.
3.9.1 El Modelo de Daño Térmico Acumulado (\(DA\))
En la práctica clínica normativa, las variaciones térmicas 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 crítico ocurre exclusivamente cuando la red de frío se rompe.
Para modelar este protocolo operativo dentro de nuestra simulación, el Daño Acumulado (\(DA\)) se calcula como la integral de la tasa de Arrhenius, condicionada por una función indicadora \(\mathbb{I}_{\{T(t) > 8^\circ\text{C}\}}\) que “enciende” el reloj de degradación únicamente cuando la temperatura cruza el límite superior de seguridad:
\[ DA(t) = \int_{0}^{t} \mathbb{I}_{\{T(s) > 8\}} \cdot k(T(s)) \, ds \]
Esta formulación integral permite que el simulador de Montecarlo contabilice el estrés térmico exacto de cada trayectoria. Si durante un apagón prolongado la EDE empuja a \(T(t)\) muy por encima de los 8°C, el integrando crecerá rápidamente. Si al final de un año simulado el \(DA(t)\) rebasa un umbral de descarte crítico predefinido, el lote de vacunas se considerará como pérdida total.