3 Modelo Matemático para Optimizar la Operación de un Sistema Hidroeléctrico.
El objetivo de este capítulo es presentar una formulación rigurosa del problema de optimización de la operación de un sistema hidroeléctrico como un Proceso de Decisión de Markov (MDP), junto con una descripción detallada del algoritmo de Fitted Q-Iteration (FQI) utilizado para aproximar la política óptima. Se enfatizará la conexión entre los elementos matemáticos formales y su implementación numérica, así como la validación operativa mediante simulación en lazo cerrado.
3.1 Formulación como Proceso de Decisión de Markov (MDP)
Siguiendo la notación de Puterman (1994), se define el MDP como la colección de objetos \[\{(T, \mathcal{S}, \mathcal{A}, \mathbb{P}, \mathcal{R})\},\] donde cada componente se especifica a continuación.
3.1.1 Estructura Temporal \((T)\)
Se considera un horizonte temporal de operación compuesto por \(T = 50\) años históricos (periodo 1959–2008), en donde cada uno de ellos se divide en \(q = 24\) quincenas. Con el fin de reducir la dimensionalidad temporal y capturar la estacionalidad, las quincenas se agrupan en \(M = 6\) etapas hidrológicas mediante la función de agregación \(\mu: \{0,\dots,23\} \to \{0,\dots,5\}\) definida por \[\mu(q) = \begin{cases} 0, & 0 \leq q <10\\ 1, & 10 \leq q <14\\ 2, & 14 \leq q < 16\\ 3, & 16 \leq q < 18\\ 4, & 18 \leq q < 20\\ 5, & 20 \leq q \leq 23\\ \end{cases} \tag{3.1}\] de modo que \(m = \mu(q)\) es la etapa hidrológica correspondiente a la quincena \(q\). Esta partición reproduce exactamente la agrupación implementada en el código (etapas_quinc), verificada como correcta.
Dos resoluciones temporales distintas. Es indispensable distinguir, desde ahora, dos usos diferentes de esta estructura temporal: (i) para construir el conjunto de aprendizaje del algoritmo de FQI (Sec. Sección 3.3.1), la afluencia se agrega a nivel de etapa (\(M=6\) decisiones por año) y no se recorre la serie quincena por quincena; (ii) para desplegar y evaluar la política ya aprendida (Sec. Sección 3.4), la decisión sí se reconsulta en cada una de las 24 quincenas del año, usando en cada una la afluencia real de esa quincena, aunque los parámetros vigentes (capacidad máxima, curva guía, grilla de acción) sean los de la etapa \(m=\mu(q)\) que la contiene.
3.1.2 Espacio de estados \((\mathcal{S})\)
El sistema está compuesto por dos presas: La Angostura \((i = 1)\) y Malpaso \((i =2)\). El volumen almacenado en la presa \(i\) en la quincena \(q\) del año \(t\) se representa por \(V_{i,t,q} \in \mathbb{R}_{+} \cup \{0\}\) (millones de metros cúbicos, Mm³).
A diferencia del enfoque de Programación Dinámica Estocástica (PDE), en el que el almacenamiento se discretiza en celdas y el estado es un índice entero, en FQI el almacenamiento se mantiene como variable continua. El parámetro \(\Delta = 600\) Mm³ (el cual permite subdividir la capacidad de almacenamiento de cada presa) y los enteros \(N_1 = 27\), \(N_2 = 16\) (capacidades discretas máximas de cada presa) —heredados de la discretización de la PDE precedente— no se usan aquí para discretizar el estado, sino únicamente para fijar las capacidades máximas continuas de cada presa, \[V_{1,\max} := N_1 \Delta = 16{,}200 \ \text{Mm}^3, \qquad V_{2,\max} := N_2 \Delta = 9{,}600\ \text{Mm}^3, \tag{3.2}\] y para convertir a volumen la curva guía, definida originalmente en unidades de índice (Sec. Sección 3.1.5).
El espacio de estados del MDP es entonces \[ \begin{aligned} \mathcal{S} &= [0,V_{1,\max}]\times[0,V_{2,\max}] \times\{0,\dots,5\} \\ &= \left\{s=(s_1,s_2,m): s_i\in[0,V_{i,\max}],\ i=1,2,\ m\in\{0,\dots,5\}\right\}. \end{aligned} \tag{3.3}\]
Grilla gruesa de referencia. Aunque el estado en sí no se discretiza, tanto la construcción del conjunto de aprendizaje (Sec. Sección 3.3.1) como la aproximación del operador de maximización (Sec. Sección 3.3.3) requieren evaluar funciones en un conjunto finito y manejable de puntos de almacenamiento. Se define para ello una grilla de referencia de \(N_s=9\) puntos equiespaciados por presa, \[ \mathcal S_i^{\text{grid}} := \Big\{0,\ \tfrac{V_{i,\max}}{8},\ \tfrac{2V_{i,\max}}{8},\ \dots,\ V_{i,\max}\Big\}, \qquad i=1,2, \tag{3.4}\] con paso \(V_{1,\max}/8=2{,}025\) Mm³ y \(V_{2,\max}/8=1{,}200\) Mm³ respectivamente. Esta grilla no es una discretización del MDP: es un dispositivo de muestreo y de cómputo; la política final se evalúa siempre en el valor continuo real del almacenamiento.
3.1.3 Espacio de acciones \((\mathcal{A})\)
Para cada etapa \(m\) se define una unidad de extracción de agua \(u_m = (100,200,300,300,400,250)\) Mm³, y un conjunto de índices factibles \(\mathcal K_{i,m}\) heredado de la grilla de decisión de la PDE: \[ \begin{aligned} \mathcal{K}_{0} &= \{0,\dots,20\}, & \mathcal{K}_{1} &= \{0,\dots,15\}, & \mathcal{K}_{2} &= \{0,\dots,7\}, \\ \mathcal{K}_{3} &= \{0,\dots,9\}, & \mathcal{K}_{4} &= \mathcal{K}_{5} = \{0,\dots,12\}. \end{aligned} \] (idéntico para ambas presas en cada etapa). Estos índices no enumeran directamente el conjunto de acciones factibles del FQI: se usan únicamente para obtener la cota superior del turbinado factible en la etapa, \[ A_{\max}(m) := u_m \cdot \max(\mathcal K_{m}). \tag{3.5}\] Por ejemplo, para \(m=0\): \(u_0=100\), \(\max(\mathcal K_0)=20\), de modo que \(A_{\max}(0)=2{,}000\) Mm³ (así, un par de índices como \((k_1,k_2)=(3,5)\) de la PDE correspondería a turbinar \(300\) y \(500\) Mm³ respectivamente —dentro de la cota—, pero el FQI no restringe la acción a estos múltiplos enteros de \(u_m\)).
El espacio de acciones factible para el MDP continuo es el intervalo \[ \mathcal{A}_m := [0,A_{\max}(m)]\times[0,A_{\max}(m)] = \{a=(a_1,a_2)\} \tag{3.6}\] con \(a_i\) el volumen turbinado en la presa \(i\), continuo. Al igual que con el estado, para fines de muestreo (Sec. Sección 3.3.1) y de maximización aproximada (Sec. Sección 3.3.3) este intervalo se discretiza en una grilla gruesa de \(N_a=9\) puntos equiespaciados, \[ \mathcal A_m^{\text{grid}} := \Big\{0,\ \tfrac{A_{\max}(m)}{8},\ \dots,\ A_{\max}(m)\Big\}, \qquad |\mathcal A_m^{\text{grid}}|=9. \tag{3.7}\]
| Etapa \(m\) | \(u_m\) (Mm³) | \(\max(\mathcal K_m)\) | \(A_{\max}(m)\) (Mm³) |
|---|---|---|---|
| 0 | 100 | 20 | 2,000 |
| 1 | 200 | 15 | 3,000 |
| 2 | 300 | 7 | 2,100 |
| 3 | 300 | 9 | 2,700 |
| 4 | 400 | 12 | 4,800 |
| 5 | 250 | 12 | 3,000 |
Los valores de \(u_m\) y \(\mathcal K_{m}\) se seleccionaron, en el diseño original de la PDE, para que \(A_{\max}(m)\) no exceda la capacidad instalada de las plantas ni los caudales históricos máximos registrados; esa cota se conserva íntegramente en el modelo continuo de FQI.
3.1.4 Dinámica de transición \((\mathbb{P})\)
Balance de masa a nivel de etapa (usado para el conjunto de aprendizaje). Para construir las transiciones de entrenamiento (Sec. Sección 3.3.1), la afluencia se agrega a nivel de etapa: sea \[ \xi_{i,t,m} := \sum_{q\,:\,\mu(q)=m} W_{i,t,q} \tag{3.8}\] la afluencia histórica del año \(t\) agregada sobre las quincenas de la etapa \(m\) (para la presa \(i\)), donde \(W_{i,t,q}\) es el afluente neto observado en la quincena \(q\) del año \(t\) (dato histórico). El balance de masa teórico (sin restricción de capacidad) al aplicar la acción \(a_i\in\mathcal A_m^{\text{grid}}\) desde un estado de la grilla \(s_i\in\mathcal S_i^{\text{grid}}\) es \[ \tilde s_i := s_i + \xi_{i,t,m} - a_i, \qquad i=1,2. \tag{3.9}\]
Balance de masa a nivel de quincena (usado en el despliegue de la política). Al simular la política ya aprendida sobre la trayectoria histórica completa (Sec. Sección 3.4), la decisión se toma en cada quincena \(q\), usando el afluente histórico de esa misma quincena: \[ V_{i,t,q+1} = V_{i,t,q} + W_{i,t,q} - a_i, \qquad q=0,\dots,22, \tag{3.10}\] y, al cerrar el año (\(q=23\)), la transición conecta con el primer registro del año siguiente, \[ V_{i,t+1,0} = V_{i,t,23} + W_{i,t+1,0} - a_i. \tag{3.11}\] En ambos casos \(a_i\) es el volumen turbinado decidido para la etapa \(m=\mu(q)\) vigente en esa quincena (Sec. Sección 3.1.3); \(a_i\) no cambia de valor dentro de la misma quincena, pero sí puede cambiar de una quincena a otra dentro de la misma etapa, porque la política se reconsulta con el estado continuo corriente en cada quincena.
Proyección al conjunto factible. En ambos regímenes, el volumen teórico \(\tilde s_i\) (o \(V_{i,t,q+1}\)) puede exceder \(V_{i,\max}\) o resultar negativo. El volumen físicamente realizable se obtiene proyectando sobre el intervalo factible, \[ s_i' := \operatorname{proj}_{[0,V_{i,\max}]}(\tilde s_i) = \min\big(\max(\tilde s_i,0),\,V_{i,\max}\big), \tag{3.12}\] donde \(\operatorname{proj}_{[0,V_{i,\max}]}\) denota al operador que da el punto del intervalo \([0,V_{i,\max}]\) más cercano a \(\tilde s_i\), y las magnitudes de exceso y faltante se registran como derrame y déficit, respectivamente, \[ \text{derrame}_i := (\tilde s_i - V_{i,\max})^+, \qquad \text{déficit}_i := (-\tilde s_i)^+, \tag{3.13}\] con \(x^+:=\max(x,0)\).
La etapa evoluciona de manera determinista y cíclica, \(m'=(m+1)\bmod M\), y el estado siguiente es \(s'=(s_1',s_2',m')\).
Desde la perspectiva de Puterman (1994), la dinámica del MDP se resume en el núcleo de transición \(\mathbb P(\cdot\,|\,s,a)\), que en este modelo no se especifica en forma analítica cerrada: toda la incertidumbre proviene de la afluencia histórica \(\xi_{i,t,m}\) (o \(W_{i,t,q}\)), tratada como una variable aleatoria cuya ley se aproxima por la distribución empírica de las \(T=50\) realizaciones históricas disponibles para cada etapa —no por una tabla de probabilidades de transición sobre estados discretos, puesto que \(\mathcal S\) es ahora un continuo—.
3.1.5 Función de recompensa \((\mathcal{R})\)
3.1.5.1 Modelo hidráulico-energético linealizado
La energía generada depende del volumen turbinado y de una aproximación lineal de la carga hidráulica. Para cada presa \(i\) se define la función de carga continua \[ \begin{aligned} H_i(s_i,s_i') &= H_{i,\max}\cdot \frac{s_i+s_i'}{2\,V_{i,\max}}, \\[4pt] H_{1,\max} &= 120\ \text{m (Angostura)}, \qquad H_{2,\max} = 100\ \text{m (Malpaso)}. \end{aligned} \tag{3.14}\] Es decir, la carga varía linealmente entre \(0\) m (presa vacía) y la carga máxima propia de cada presa cuando ésta alcanza su capacidad plena, y se aproxima la carga efectiva durante la etapa/quincena mediante el promedio de la carga al inicio (\(s_i\)) y al final (\(s_i'\)) del intervalo.
Sea \(\eta = 0.9\) la eficiencia de conversión hidráulica-eléctrica y \(g = 9.81\) m/s² la aceleración de la gravedad. La energía generada por la presa \(i\) al turbinar \(a_i\) (Mm³) bajo la carga \(H_i(s_i,s_i')\) es \[ E_i(s,a,s') = \frac{\eta\, g\, H_i(s_i,s_i')\,(a_i\times 10^6)}{3600\times 10^6} \quad [\text{GWh}], \tag{3.15}\] donde \(10^6\) convierte Mm³ a m³ y \(3600\times10^6\) convierte julios-equivalentes a GWh (con la densidad del agua \(\rho=1000\) kg/m³ absorbida algebraicamente en la constante \(3600\)). La energía total generada por el sistema en la etapa/quincena es \[ E = E_1+E_2. \tag{3.16}\]
3.1.5.2 Estructura de penalizaciones
La recompensa incorpora penalizaciones lineales, con coeficientes calibrados para actuar como restricciones blandas de tipo big-M (su magnitud domina ampliamente el valor energético del agua, ver ejemplo numérico más abajo). Se definen, por presa \(i\), las siguientes tres penalizaciones:
Derrame (excedente sobre la capacidad máxima): \[ \Pi_i^{\text{derr}} = c_{\text{derr}}\cdot(\tilde s_i - V_{i,\max})^+, \qquad c_{\text{derr}} := \frac{C_{\text{derr}}}{\Delta}. \tag{3.17}\]
Déficit (faltante respecto de un balance no negativo). Proporcional a la magnitud del déficit, no una constante fija. \[ \Pi_i^{\text{def}} = c_{\text{def}}\cdot(-\tilde s_i)^+, \qquad c_{\text{def}} := \frac{C_{\text{def}}}{\Delta}. \tag{3.18}\]
Violación de la curva guía. Se penaliza cuando el almacenamiento proyectado excede la curva guía (un techo operativo de resguardo). \[ \Pi_i^{\text{CG}} = c_{\text{CG}}\cdot\big(s_i' - \text{CG}_i(m)\big)^+, \qquad c_{\text{CG}} := \frac{C_{\text{CG}}}{\Delta}. \tag{3.19}\]
con \(C_{\text{derr}}=C_{\text{def}}=10{,}000\), \(C_{\text{CG}}=5{,}000\), de modo que \(c_{\text{derr}}=c_{\text{def}}\approx16.67\) y \(c_{\text{CG}}\approx8.33\) por Mm³. A modo de referencia de escala: turbinar 1 Mm³ bajo una carga representativa de 100 m genera, por Ecuación 3.15, apenas \(\approx0.245\) GWh —es decir, \(c_{\text{derr}}\) es unas \(68\) veces mayor que el valor energético marginal del mismo volumen, confirmando su papel de penalización dominante y no de costo económico real.
\(\text{CG}_i(m)\) es la curva guía de la presa \(i\) en la etapa \(m\), definida originalmente en unidades de índice, \(\text{CG}_1^{\text{idx}}=(25,22,24,25,26,27)\), \(\text{CG}_2^{\text{idx}}=(16,15,14,14,15,16)\), y convertida a volumen mediante \[ \text{CG}_i(m) := \Delta\cdot \text{CG}_i^{\text{idx}}(m). \tag{3.20}\]
La función de recompensa total, sumada sobre ambas presas, es \[ \mathcal{R}(s,a,s') = \sum_{i=1}^{2}\Big[E_i - \Pi_i^{\text{derr}} - \Pi_i^{\text{def}} - \Pi_i^{\text{CG}}\Big]. \tag{3.21}\]
En el presente trabajo no se modela una distribución analítica explícita de la afluencia; en su lugar se construye un conjunto de aprendizaje empírico, descrito con precisión en la Sec. Sección 3.3.1 —y que, como ahí se detalla, no consiste en recorrer trayectorias históricas de decisiones reales, sino en evaluar la dinámica y la recompensa sobre una malla exhaustiva de pares estado-acción contra la totalidad de las realizaciones históricas de afluencia. Este enfoque se alinea con un esquema de Batch Reinforcement Learning “parcialmente libre de modelo” (Castelletti et al. (2010); Ernst et al. (2005)).
3.2 Ecuación de optimalidad de Bellman
Se adopta la notación \(Q^*\) para la función de valor óptima, por ser función del par estado-acción \((s,a)\) —lo que la literatura de aprendizaje por refuerzo llama función de valor-acción— y por consistencia con el nombre de la variable en el código (Q).
Conforme a la teoría clásica de MDPs descontados (Puterman (1994), Teorema 6.2.3), existe una única función de valor-acción óptima \(Q^*\) que satisface la ecuación de optimalidad de Bellman. La razón intuitiva es el factor de descuento \(\gamma=0.95<1\): la influencia de una decisión de hoy sobre lo que ocurre dentro de, digamos, 50 etapas está multiplicada por \(\gamma^{50}\approx0.08\), y sigue decayendo geométricamente cuanto más lejos se mire hacia el futuro. Ese decaimiento garantiza que, sin importar con qué aproximación inicial se empiece —incluso \(Q^{(0)}\equiv0\), como hace el algoritmo (Sec. Sección 3.3.1 y siguientes)—, al aplicar la ecuación de Bellman una y otra vez las sucesivas aproximaciones se acercan cada vez más entre sí, hasta converger a una única \(Q^*\). Esto es válido tanto si \(\mathcal S\) es finito (el caso clásico de Puterman) como si tiene una componente continua, como en este modelo (Sec. Sección 3.1.2). (Formalmente, esta propiedad se conoce como que el operador de Bellman es una contracción, y su justificación rigurosa se apoya en el teorema del punto fijo de Banach; no se desarrolla aquí por no ser el objeto central de este capítulo.)
Dado que la incertidumbre proviene de la afluencia histórica \(\xi=(\xi_1,\xi_2)\), cuya ley se aproxima por su distribución empírica \(\widehat P_m\) sobre las \(T=50\) realizaciones históricas de la etapa \(m\) (Sec. Sección 3.3.1), la ecuación de optimalidad se escribe como una esperanza sobre dicha distribución empírica —y no como una suma sobre un conjunto discreto de estados siguientes, formalismo que ya no aplica al ser \(\mathcal S\) un continuo—: \[ Q^{*}(s,a) = \mathbb{E}_{\xi\sim\widehat P_m}\Big[\,\mathcal R(s,a,s') + \gamma\max_{a'\in\mathcal A_{m'}} Q^{*}(s',a')\,\Big], \quad s' = f(s,a,\xi), \tag{3.22}\] con \(f\) la función de transición y proyección de las ecuaciones Ecuación 3.9–Ecuación 3.12, y \(m'=(m+1)\bmod M\).
Dado que no se dispone de \(\widehat P_m\) en forma cerrada sino de sus \(T\) realizaciones muestrales, el algoritmo de FQI (Sec. siguiente) trabaja con la versión por muestra de esta ecuación —el objetivo empírico de Bellman, Ecuación 3.25—, evaluada sobre cada una de las transiciones del conjunto de aprendizaje, y no con la esperanza poblacional (Ecuación 3.22) directamente.
3.3 Aproximación numérica mediante Fitted Q-Iteration (FQI)
La Programación Dinámica Estocástica está afectada por la maldición de la dimensionalidad: el costo computacional crece exponencialmente con el número de variables de estado, y requiere además un modelo de transición explícito (Castelletti et al. (2010)). Por ello, el presente trabajo aproxima \(Q^*\) mediante Fitted Q-Iteration, un algoritmo de aprendizaje por refuerzo batch que ajusta iterativamente un aproximador funcional supervisado a partir de un conjunto de transiciones muestreadas.
3.3.1 Conjunto de aprendizaje “parcialmente libre de modelo”
Concretamente, el conjunto de aprendizaje \(\mathcal F=\{(x_j, r_j, s'_{1,j}, s'_{2,j}, m'_j)\}_{j=1}^N\) se construye enumerando el producto cartesiano de: la grilla de estado Ecuación 3.4, la grilla de acción Ecuación 3.7, las \(M=6\) etapas, y las \(T=50\) realizaciones históricas de afluencia por etapa Ecuación 3.8. El vector de características de cada muestra es \[ x_j = (s_{1,j}, s_{2,j}, m_j, a_{1,j}, a_{2,j}) \in \mathcal S_1^{\text{grid}}\times\mathcal S_2^{\text{grid}}\times\{0,\dots,5\}\times\mathcal A^{\text{grid}}_{m_j}\times\mathcal A^{\text{grid}}_{m_j}, \tag{3.23}\] y, para cada combinación de estado-acción-etapa, se generan \(T=50\) muestras —una por cada año histórico \(t\)— aplicando Ecuación 3.9, Ecuación 3.12 y Ecuación 3.21 con la afluencia real \(\Xi_{i,t,m_j}\) de ese año. El tamaño total del conjunto de aprendizaje es \[ N = M\times|\mathcal S_1^{\text{grid}}|\times|\mathcal S_2^{\text{grid}}|\times|\mathcal A_1^{\text{grid}}|\times|\mathcal A_2^{\text{grid}}|\times T = 6\times9\times9\times9\times9\times50 = 1{,}968{,}300. \tag{3.24}\]
Esta receta de muestreo —barrer exhaustivamente estado y acción sobre una grilla, pero conservar todas las realizaciones históricas de la perturbación en cada punto, sin colapsarlas a una tabla de probabilidades ni ajustarles una distribución paramétrica— es la que Castelletti et al. (2010) denominan “parcialmente libre de modelo” (partial model-free): la dinámica del balance de masa sí se modela analíticamente, pero la ley de la afluencia se deja no paramétrica.
3.3.2 Iteración de Bellman empírica
Se inicializa \(Q^{(0)}(s,a) \equiv 0\) para todo \((s,a)\). Para cada transición \((x_j,r_j,s'_{1,j},s'_{2,j},m'_j)\in\mathcal F\) del conjunto de aprendizaje, se define el objetivo de Bellman de la iteración \(n\) por \[ y^{(n)}_j := r_j + \gamma \max_{a'\in\mathcal A_{m_j'}} Q^{(n-1)}(s'_{1,j},s'_{2,j},m_j',a'), \tag{3.25}\] que estima el valor-acción óptimo de \(x_j\) combinando la recompensa observada \(r_j\) y el mejor valor alcanzable desde el estado siguiente, según la aproximación de la iteración anterior. El cómputo del máximo sobre \(a'\), dado que \((s'_{1,j},s'_{2,j})\) es en general un punto continuo que no coincide con la grilla Ecuación 3.4, se describe en la Sec. Sección 3.3.3.
3.3.3 Aproximación del operador de maximización sobre estados continuos
**Esta subsección resuelve un punto que el texto anterior dejaba implícito sin justificación: cómo evaluar \(\max_{a'}Q^{(n-1)}(s',a')\) cuando \(s'\) es continuo. Evaluarlo por fuerza bruta para cada una de las \(N\approx 2\times10^6\) muestras sería prohibitivo, así que el algoritmo precomputa, para cada etapa \(m\), una superficie de valor máximo sobre la grilla de estado y la interpola:
Evaluación sobre grilla. Para cada etapa \(m\), se evalúa \(Q^{(n-1)}\) en todo el producto \(\mathcal S_1^{\text{grid}}\times\mathcal S_2^{\text{grid}}\times\{m\}\times\mathcal A_1^{\text{grid}}(m)\times\mathcal A_2^{\text{grid}}(m)\), y se maximiza sobre las dos coordenadas de acción: \[ \widehat V^{(n-1)}(s_1^{(p)},s_2^{(q)},m) := \max_{a_1,a_2\in\mathcal A_m^{\text{grid}}} Q^{(n-1)}(s_1^{(p)},s_2^{(q)},m,a_1,a_2), \qquad p,q=1,\dots,9, \tag{3.26}\] obteniéndose una matriz \(9\times9\) por etapa.
Interpolación bilineal. Sobre esta matriz se construye un interpolador continuo \(\widetilde V^{(n-1)}(\cdot,\cdot\,;m):\mathbb R^2\to\mathbb R\) (interpolación multilineal sobre grilla regular), con extrapolación lineal fuera de \([0,V_{1,\max}]\times[0,V_{2,\max}]\).
Sustitución en el objetivo de Bellman. El máximo de Ecuación 3.25 se aproxima evaluando el interpolador de la etapa siguiente en el punto continuo genuinamente alcanzado: \[ \max_{a'\in\mathcal A_{m_j'}} Q^{(n-1)}(s'_{1,j},s'_{2,j},m_j',a') \ \approx\ \widetilde V^{(n-1)}\big(s'_{1,j},s'_{2,j};\,m_j'\big). \tag{3.27}\]
Esta aproximación reduce el costo de \(O(N\times|\mathcal A^{\text{grid}}|^2)\) evaluaciones del regresor a \(O\big(M\times|\mathcal S^{\text{grid}}|^2\times|\mathcal A^{\text{grid}}|^2 + N\big)\), al reutilizar la misma superficie interpolada para las \(N/M\) muestras de cada etapa.
3.3.4 Aproximación funcional por regresión no paramétrica (ExtraTrees)
Una vez calculados los objetivos de Bellman \(y^{(n)}\) (Ecuación 3.25, con el máximo aproximado según Ecuación 3.27), se ajusta un aproximador \(Q^{(n)}\) que minimiza el error cuadrático entre las predicciones \(Q^{(n)}(x_j)\) y los objetivos \(y^{(n)}_j\) sobre todas las transiciones del conjunto de aprendizaje. Se utiliza, como propone Ernst et al. (2005), un ensamble de ExtraTrees (árboles extremadamente aleatorizados), adecuado para aproximar funciones de valor en espacios continuos.
Hiperparámetros (corregidos respecto de la versión previa):
| Hiperparámetro | Símbolo | Valor implementado | |
|---|---|---|---|
| Número de árboles | \(M_{\text{tree}}\) | 300 | |
| Profundidad máxima | — | sin límite (max_depth=None) |
|
| Muestras mínimas por hoja | \(n_{\min}\) | 50 (\(=T\)) | |
| Variables candidatas por corte | \(K\) | \(n=5\) (max_features=1.0) |
El valor \(n_{\min}=T=50\) no es arbitrario: al exigir que cada hoja del árbol contenga al menos tantas muestras como años históricos de afluencia, se favorece que cada hoja agrupe una fracción representativa de las \(T\) realizaciones de \(\xi\) asociadas a un mismo \((s,a)\) aproximado, de modo que el promedio dentro de la hoja aproxime, por la ley (empírica) de los grandes números, la esperanza condicional \(\mathbb E_{\xi\sim\widehat P_m}[\cdot]\) que aparece en Ecuación 3.22. Cada árbol particiona el espacio estado-acción en regiones (hojas) mediante cortes con umbral aleatorio (no exhaustivo, a diferencia de CART/Random Forest) sobre las \(K=n=5\) variables de entrada, y predice, para un punto de consulta, el promedio del objetivo de regresión entre las muestras de entrenamiento de su misma hoja. La predicción del ensamble es el promedio de las \(M_{\text{tree}}=300\) predicciones individuales.
Nota sobre la Figura de árbol de decisión. La versión previa incluía una figura con la leyenda “Árbol de decisión con profundidad máxima de 14 y un mínimo de 5 muestras por hoja”. Dado que el código real no limita la profundidad (
max_depth=None) y usa \(n_{\min}=50\), esa leyenda —y, en la práctica, la figura misma, que mostraría un árbol truncado no representativo— debe corregirse o sustituirse por una ilustración esquemática que no afirme una profundidad máxima inexistente.
El algoritmo repite este procedimiento iterativamente.
3.3.5 Criterio de paro basado en desempeño simulado
Un criterio de convergencia del tipo punto-fijo —detener cuando \(\|Q^{(n)}-Q^{(n-1)}\|\) cae bajo una tolerancia— no es aplicable cuando el aproximador es un ensamble de Extra-Trees: la aleatorización propia de la construcción de cada árbol (elección aleatoria de umbrales de corte) hace que \(Q^{(n)}\) y \(Q^{(n-1)}\) difieran por ruido de ajuste incluso si ya se hubiese alcanzado, en esencia, la función de valor-acción óptima; esta distancia no converge a cero. Así lo argumentan explícitamente Castelletti et al. (2010), y así lo implementa el código, que en su lugar monitorea el desempeño operativo simulado de la política inducida en cada iteración.
Simulación de la política golosa. Cada \(\Delta h\) iteraciones (en el código, \(\Delta h=4\)), se simula la operación del sistema bajo la política golosa de \(Q^{(n)}\) sobre la trayectoria histórica completa, en resolución de quincena (Sec. Sección 3.4), obteniendo energías \(e^{(i)}_{t,q}\), derrames \(\delta^{(i)}_{t,q}\) y déficits \(\epsilon^{(i)}_{t,q}\) por presa.
Funcional de desempeño. Se define \[ \widehat J_n := \sum_{t,q}\big(e^{(1)}_{t,q}+e^{(2)}_{t,q}\big) - c_{\text{derr}}\sum_{t,q}\big(\delta^{(1)}_{t,q}+\delta^{(2)}_{t,q}\big) - c_{\text{def}}\sum_{t,q}\big(\epsilon^{(1)}_{t,q}+\epsilon^{(2)}_{t,q}\big), \tag{3.28}\] una medida de desempeño agregada, no descontada y evaluada con afluencias reales quincenales —distinta, por tanto, del objetivo de ajuste Ecuación 3.25, que opera a nivel de una transición individual y con descuento \(\gamma\)—.
Regla de paro. Se detiene la iteración la primera vez que se observa una caída de \(\widehat J_n\) respecto del checkpoint inmediatamente anterior (una “reversión”): \[ n^\star := \min\{\, n-\Delta h \ :\ \widehat J_{n} < \widehat J_{n-\Delta h} \,\}, \tag{3.29}\] interpretando la reversión como evidencia de que, a partir de ese punto, predominan las fluctuaciones aleatorias del ajuste sobre el aprendizaje genuino (Castelletti et al. (2010)). Si no se observa ninguna reversión dentro del número máximo de iteraciones permitido, se usa como respaldo el mejor checkpoint observado, \(n^\star:=\arg\max_n \widehat J_n\). La aproximación final es \(Q^{*}:=Q^{(n^\star)}\), recuperada de una caché de modelos serializados guardada en cada punto de verificación (necesaria porque, por la aleatorización del ensamble, \(Q^{(n^\star)}\) no podría reproducirse exactamente reentrenando desde cero).
3.4 Política óptima y simulación operativa
Una vez obtenida \(Q^{*}\), se deriva la política óptima (golosa) \(\pi^{*}:\mathcal S\to\mathcal A\), \[ \pi^{*}(s) := \operatorname*{arg\,max}_{a\in\mathcal A_m^{\text{grid}}} Q^{*}(s,a), \tag{3.30}\] donde la maximización se realiza, en la práctica, mediante búsqueda exhaustiva sobre la grilla de acción de la etapa vigente Ecuación 3.7 (no sobre el conjunto continuo \(\mathcal A_m\)), evaluando \(Q^{*}\) en el estado continuo real —no discretizado— del sistema.
La política se valida mediante simulación en lazo cerrado sobre la serie histórica completa, siguiendo la práctica estándar de validación de políticas de control aprendidas documentada en Bertsekas (2019). El algoritmo, corregido para ser fiel al código, opera así:
Inicialización: se fijan los volúmenes iniciales \(V_{1,0}=10{,}000\) Mm³, \(V_{2,0}=6{,}000\) Mm³ (valores del código; no se especificaban antes).
Mapeo estacional: para cada quincena \(q\) (recorriendo \(t=1,\dots,T\) y, dentro de cada año, \(q=0,\dots,23\)), se determina \(m=\mu(q)\) y la grilla de acción \(\mathcal A_m^{\text{grid}}\) vigente.
Discretización: \(s_i=\mathcal D(V_i,N_i)\)— paso eliminado: el estado usado para consultar \(Q^*\) es el volumen continuo real \((V_{1,t,q}, V_{2,t,q})\), sin redondeo a índice alguno.Selección de acción, por búsqueda exhaustiva sobre la grilla de acción de 9×9 combinaciones de la etapa vigente: \[ (a_1^*,a_2^*) = \operatorname*{arg\,max}_{(a_1,a_2)\,\in\,\mathcal A_1^{\text{grid}}(m)\times\mathcal A_2^{\text{grid}}(m)} Q^{*}(V_{1,t,q},V_{2,t,q},m,a_1,a_2). \]
Balance hídrico, con el afluente de la misma quincena \(q\) (corrección respecto de la versión previa, que usaba \(q+1\)): \[ V_{i,t,q+1} = V_{i,t,q} + W_{i,t,q} - a_i^*, \qquad i=1,2. \]
Proyección física, incluyendo tanto el derrame como el déficit (la versión previa omitía el déficit): \[ V_i^{\text{new}} = \operatorname{proj}_{[0,V_{i,\max}]}(V_{i,t,q+1}), \qquad \text{derrame}_i = (V_{i,t,q+1}-V_{i,\max})^+, \qquad \text{déficit}_i = (-V_{i,t,q+1})^+. \]
Cálculo energético, con la carga continua y las cargas máximas distintas por presa (corrección de Ecuación 3.14): \[ H_i = H_{i,\max}\cdot\frac{V_{i,t,q}+V_i^{\text{new}}}{2\,V_{i,\max}}, \qquad E_i = \frac{\eta g\, H_i\,(a_i^*\times10^6)}{3600\times10^6}. \]
Registro y actualización: se almacenan energía, derrame, déficit y almacenamiento; el estado se actualiza a \(V_i^{\text{new}}\) para la siguiente quincena.
3.4.1 Métricas de desempeño y validación estadística
Las tablas siguientes (Tabla 3.3, Tabla 3.4) y la conclusión que las acompañaba en la versión previa reportan un derrame total del sistema de \(\approx 34{,}053\) Mm³ y una lectura crítica sobre la mala gestión de la política. Esta cifra es incompatible con la propia validación documentada en el encabezado del script corregido, que reporta “0 Mm³ de derrame y 0 Mm³ de déficit en ambas presas” al ejecutar exactamente el algoritmo aquí descrito. Es decir, las cifras de la versión previa corresponden, con alta probabilidad, a una ejecución de una versión anterior y no fiel del modelo (la que originó la auditoría documentada en el encabezado de FQI_corregido.py), no a la versión actual. No se sustituyen aquí por cifras nuevas para no fabricar resultados: se recomienda volver a ejecutar el script corregido sobre los datos reales y regenerar ambas tablas, así como reescribir el párrafo de interpretación en función de los resultados que se obtengan.
simulacion_FQI_corregida.csv).
| Métrica | Angostura | Malpaso | Sistema Total |
|---|---|---|---|
| Energía promedio quincenal [GWh] | (pendiente) | (pendiente) | (pendiente) |
| Energía anual estimada [GWh] | (pendiente) | (pendiente) | (pendiente) |
| Contribución al sistema [%] | (pendiente) | (pendiente) | 100.0 |
| Métrica | Angostura | Malpaso | Sistema Total |
|---|---|---|---|
| Derrame total [Mm³] | (pendiente) | (pendiente) | (pendiente) |
| Déficit total [Mm³] | (pendiente) | (pendiente) | (pendiente) |
Una vez regeneradas las tablas a partir de una corrida real del script, la interpretación debe contrastar específicamente los totales de derrame y déficit contra la validación reportada en el propio código (0 Mm³ en ambos rubros) antes de emitir un juicio sobre la calidad de la política aprendida; de mantenerse una discrepancia con lo documentado en el script, ésta debería investigarse como posible diferencia entre el entorno de ejecución (datos, semilla aleatoria, versión de librerías) y no asumirse automáticamente como una falla de la política.
3.5 Conclusiones del capítulo
- Síntesis de contribuciones matemáticas y numéricas.
- Reafirmación de la coherencia entre formulación teórica, implementación algorítmica y validación operativa —pendiente de completar una vez regeneradas las métricas de la Sec. Sección 3.4.1—.
- Declaración de cierre alineada con los objetivos de la tesis.