Gestión óptima en microrredes con soporte fotovoltaico e hidrógeno verde

Revista Iberoamericana de Autom´atica e Inform´atica Industrial 22 (2025) 135-145 www.revista-riai.org Gesti´on ´optima en microrredes con soporte fotovoltaico e hidr´ogeno verde Rub´en Moliner-Heredia, Carlos Vivas∗, Francisco R. Rubio Dpto. de Ingenier´ıa de Sistemas y Autom´atica, Escuela T´ecnica Superior de Ingenier´ıa (Universidad de Sevilla) Camino de los Descubrimientos, s/n, 41092, Sevilla. Resumen La gesti´ on de la energ´ ıa es esencial para un correcto control de una microrred. Es importante que dicha gesti´ on tenga en cuenta la optimizaci´ on de la vida ´ util de los componentes de la microrred. Uno de los posibles comportamientos que se desean evitar es el apagado y encendido continuo de determinados componentes, tales como las celdas de combustible y los electrolizadores. En este art´ ıculo se ha propuesto un algoritmo usando control predictivo basado en modelo (MPC) utilizando diversas restricciones de tiempos m´ ınimos de activaci´ on de dichos componentes para evitar efectos perjudiciales sobre los equipos de hidr´ ogeno verde. Adem´ as, se ha propuesto un m´ etodo de transmisi´ on entre iteraciones del algoritmo para que la aplicaci´ on de dichas restricciones sea compatible con errores en las predicciones del algoritmo y de perturbaciones en el sistema. Para validar el algoritmo propuesto, se ha desarrollado un modelo simplificado de la microrred, y se han realizado simulaciones y comparaciones con varios algoritmos utilizando Matlab. Palabras clave: Control de recursos de Energ´ ıa Renovable, Control Predictivo basado en Modelo, Redes el´ ectricas inteligentes, Sistemas de gesti´ on y distribuci´ on de la Energ´ ıa, Control inteligente de sistemas el´ ectricos de potencia Optimal management in microgrids with photovoltaic and green hydrogen support Abstract Energy management is essential in order to develop a correct control of a microgrid. It is important that said management takes into account the optimization of the useful life of the microgrid components. One of the possible behaviors to be avoided is the continuous switching on and offof certain components, such as the fuel cell and the electrolyzer. In this paper, an algorithm using Model Predictive Control (MPC) using several minimum activation time constraints for these components in order to avoid these detrimental effects on the green H2 equipment has been proposed. Besides, a transmission method between iterations of the algorithm in order to allow compatibility of the application of the aforementioned constraints with the algorithm prediction errors and disturbances in the system has also been proposed. To validate the proposed algorithm, a simplified model of the microgrid has been proposed, and several simulations and comparisons with several algorithms using Matlab have been performed. Keywords: Control of renewable energy resources, Model predictive control, Smart grids, Energy and Distribution Management Systems, Intelligent control of power systems 1. Introducci´ on Durante la ´ ultima d´ ecada, los sistemas el´ ectricos han evolucionado significativamente, centr´ andose en la b´ usqueda de fuentes renovables y respetuosas con el medio ambiente. Esto se debe a una combinaci´ on de factores tecnol´ ogicos, tales como las mejoras en el rendimiento de los componentes de los elementos de captaci´ on de energ´ ıa renovable; factores sociales, dado que la sociedad es cada vez m´ as consciente de los problemas del medio ambiente (Bordons et al., 2020), as´ ı como de ∗Autor para correspondencia: [email protected] Attribution-NonCommercial-ShareAlike 4.0 International (CC BY-NC-SA 4.0) To cite this article: Moliner-Heredia, R., Vivas, C., Rubio, F. R., 2025. Optimal management in microgrids with photovoltaic and green hydrogen support. Revista Iberoamericana de Automática e Informática Industrial 22, 135-145. https://doi.org/10.4995/riai.2024.21922 factores geopol´ ıticos, tales como los requisitos de descarbonizaci´ on, los esfuerzos por lograr una autonom´ ıa energ´ etica y el aumento de los precios de la energ´ ıa. Todos estos factores han promovido el desarrollo de sistemas el´ ectricos de consumo y generaci´ on descentralizada y autosuficiente. Estos nuevos sistemas suelen tomar la forma de microrredes, conectadas a niveles de baja potencia. Una microrred se puede definir como un conjunto de generadores de energ´ ıa, de cargas controlables o incontrolables y de sistemas de almacenamiento de energ´ ıa, que puede funcionar conectado a la red o de manera aislada (Kumar and Saravanan, 2017). El control de las microrredes se puede definir como una estructura jer´ arquica separada en distintos niveles de control (Mahmoud et al., 2017): uno primario, independiente, que act´ ua sobre los componentes electr´ onicos en ´ ordenes temporales de reducida magnitud para garantizar los valores correctos de voltaje y frecuencia de la microrred; un secundario, que gestiona cualquier sobrecompensaci´ on causada por la gesti´ on del nivel primario; y el terciario, encargado de gestionar y optimizar flujos de potencia entre posibles microrredes interconectadas y la red el´ ectrica general. En la pr´ actica se articula en torno al concepto de sistema de gesti´ on de la energ´ ıa (EMS) que proporciona servicios de regulaci´ on de flujos de potencia, gesti´ on eficiente de la energ´ ıa, balances de voltaje/frecuencia, etc . No obstante, en la literatura se pueden encontrar otros m´ etodos de distribuci´ on de la gesti´ on de las distintas tareas que debe realizar el control de una microrred (Sen and Kumar, 2018). Los EMS tienen el objetivo de monitorizar los datos recibidos desde la microrred, predecir la evoluci´ on de la generaci´ on y la demanda, calcular el procedimiento de intercambio de energ´ ıa entre fuentes de generaci´ on, sistemas de almacenamiento y la red el´ ectrica que d´ e lugar a ´ optimos energ´ eticos o econ´ omicos, y gestionar la organizaci´ on de la aplicaci´ on de dichos pasos. No obstante, durante el proceso de optimizaci´ on hay que lidiar con la incertidumbre inherente a la predicci´ on del tiempo, del consumo, de los precios de la energ´ ıa proporcionada por la red, o a la aparici´ on de perturbaciones, cortes o fallos del sistema (Cabrera-Tobar et al., 2022). Para minimizar el impacto de estas incertidumbres, en la literatura se han propuesto diversas t´ ecnicas de control. Estas se han clasificado seg´ un su m´ etodo en dos grupos Ali et al. (2021): basadas en m´ etodos de inteligencia artificial, como redes neuronales (Lopez-Garcia et al., 2020) o m´ etodos de l´ ogica borrosa Chen et al. (2012), y basadas en m´ etodos cl´ asicos, como m´ etodos estoc´ asticos (Hu et al., 2016) o control predictivo basado en modelo (MPC) (Marquez et al., 2021). Concretamente, el m´ etodo MPC puede gestionar las incertidumbres y perturbaciones gracias a su mecanismo retroactivo, adem´ as de que permite incluir restricciones f´ ısicas y el comportamiento futuro del sistema. Tambi´ en permite aplicar restricciones temporales del comportamiento del sistema dentro de su algoritmo, de forma que se puedan evitar encendidos y apagados continuos de los componentes de la microrred (la cual es una soluci´ on ´ optima local), modo de operaci´ on ´ este que los deteriora en general y reduce su vida ´ util. En este art´ ıculo se ha desarrollado un algoritmo para la gesti´ on de la energ´ ıa de una microrred utilizando control de horizonte deslizante MPC, al cual se han a˜ nadido una serie de restricciones adicionales que restringen de manera laxa el tiempo m´ ınimo de activaci´ on de la celda de combustible y el electrolizador, con el objetivo de evitar la activaci´ on y desactivaci´ on continua de dichos componentes, evitando as´ ı su desgaste prematuro. Para validar el algoritmo se han realizado diversos ensayos de simulaci´ on en Matlab. La estructura del art´ ıculo se muestra a continuaci´ on: primero, en la Secci´ on 2 se describe el objetivo del art´ ıculo junto con los requisitos solicitados al algoritmo. En la Secci´ on 3 se describe la microrred y su modelo simulado. En la Secci´ on 4 se describe el algoritmo de control propuesto y se detalla la obtenci´ on de algunas restricciones. En la Secci´ on 5 se explica el caso de estudio, y en la Secci´ on 6 se muestran los resultados y su an´ alisis. Por ´ ultimo, en la Secci´ on 7 se detallan las conclusiones. 2. Objetivo El objetivo de este art´ ıculo es desarrollar un algoritmo de control de una microrred en su nivel terciario, el sistema de gesti´ on de la energ´ ıa. Los requisitos que se desea que satisfaga el algoritmo de control de la microrred son los siguientes: 1. Gesti´ on de la demanda de potencia. El algoritmo propuesto debe siempre satisfacer las necesidades de la demanda. 2. Optimizaci´ on econ´ omica. El algoritmo propuesto debe reducir la cantidad de potencia comprada desde la red el´ ectrica. 3. Restricciones de tiempo m´ ınimo de activaci´ on. Dado que la activaci´ on y desactivaci´ on reiterada de los componentes de la microrred reduce su vida ´ util, el algoritmo debe impedir que se produzca este comportamiento. 4. Eficiencia general. El algoritmo debe evitar producir recirculaciones e intercambios energ´ eticos innecesarios entre sus componentes. 5. Robustez. El algoritmo debe ser capaz de ser robusto frente a errores en la predicci´ on de la generaci´ on y la demanda, as´ ı como de desconexiones repentinas de componentes de la microrred. 3. Descripci´ on del modelo 3.1. Componentes de la microrred La microrred propuesta consta de los siguientes componentes: una planta de paneles fotovoltaicos, que act´ uan como la unidad generadora de potencia renovable para la microrred; unas cargas cuya activaci´ on depende de la actividad de los usuarios, cuyas necesidades de consumo de potencia deben satisfacerse siempre; una unidad de bater´ ıas, que proporcionan un almacenaje de potencia r´ apido y a corto plazo; una conexi´ on con la red el´ ectrica, dado que la microrred debe ser capaz de alternar su funcionamiento entre su modo aislado y conectado cuando le sea indicado; y una unidad de hidr´ ogeno verde, que permite un almacenamiento de potencia a largo plazo, pero con el coste de tener un rendimiento reducido. La unidad de almacenamiento de hidr´ ogeno verde consta de un electrolizador, que separa el hidr´ ogeno del agua utilizando electricidad, un dep´ osito de hidr´ ogeno verde, donde el hidr´ ogeno se almacena en forma de hidruros met´ alicos, y una 136 Moliner-Heredia, R. et al. / Revista Iberoamericana de Automática e Informática Industrial 22 (2025) 125-135 Figura 1: Estructura de la microrred. celda de combustible, que produce electricidad al combinar el hidr´ ogeno con ox´ ıgeno para formar agua. En la Figura 1 se puede observar la estructura de esta microrred. 3.2. Modelo de simulaci´on de la microrred El modelo de la microrred propuesta est´ a conformado por un sistema discreto de ecuaciones lineales que representan los equilibrios de potencia de la microrred y el estado de carga de las bater´ ıas y el dep´ osito de hidr´ ogeno verde. Primero, para cualquier instante de tiempo discreto de la simulaci´ on tm, el equilibrio de potencia en la microrred se representa de la siguiente manera: PGRIDm(tm)=PEZm(tm)+PBAT,CHm(tm)+PDEMm(tm) −PGENm(tm)+PFCm(tm)+PBAT,DS m(tm),(1) donde PGRIDmes la potencia intercambiada con la red el´ ectrica, PBAT,CHmyPBAT,DS mrepresentan la potencia cargada y descargada de la unidad de bater´ ıas, PEZmes la potencia consumida por el electrolizador, PFCmla potencia generada por la celda de combustible, y PGENmyPDEMmla potencia generada y demandada, respectivamente. En este caso, PEZm,PBAT,CHm,PDEMm, PGENm,PFCm,PBAT,DS m∈R+, y por otro lado, PGRIDm(tm)∈R. Si el valor de PGRIDmes positivo, se est´ a vertiendo potencia desde la red hacia la microrred, y viceversa si es negativo. La evoluci´ on temporal del estado de carga de las bater´ ıas (S OCm) se puede modelizar utilizando la siguiente ecuaci´ on: S OCm(tm+1) =S OCm(tm)−KBAT,DS mPBAT,DS m(tm) +KBAT,CHmPBAT,CHm(tm),(2) donde S OCm(tm)∈R+|0≤S OCm(tm)≤100 y KBATm,CH,KBATm,DS ∈R+son los coeficientes de eficiencia que relacionan el intercambio de potencia con S OCm. Ambos coeficientes no var´ ıan con el paso del tiempo, es decir, no se considera el desgaste de la capacidad de las bater´ ıas. Por otro lado, la evoluci´ on temporal del estado del almacenamiento del dep´ osito de hidr´ ogeno verde (MHL) se modeliza utilizando la siguiente ecuaci´ on: MHLm(tm+1) =MHLm(tm)−KFCmPFCm(tm) +KEZmPEZm(tm),(3) donde MHLm(tm)∈R+|0≤MHLm(tm)≤100, y KEZm,KFCm∈R+son los coeficientes de eficiencia que relacionan la conversi´ on entre la potencia generada y consumida con la cantidad de hidr´ ogeno verde generada. Ambos coeficientes no var´ ıan con el paso del tiempo, es decir, no se considera el desgaste de la celda de combustible ni del electrolizador. Adem´ as, para el modelo se definen dos ecuaciones adicionales: PNETm(tm)=PGENm(tm)−PDEMm(tm),(4) donde PNETm∈Res la potencia neta que entra al modelo, y engloba a ambas variables no controlables. Por otro lado, se define: PBATm(tm)=PBAT,DS m(tm)−PBAT,CHm(tm),(5) donde PBATm∈Res la potencia total que procede o se almacena en la bater´ ıa. Nota: Dado que el modelo se define como un conjunto de ecuaciones lineales, y solamente los niveles de carga del dep´ osito y las bater´ ıas est´ an modelizados como variantes en el tiempo, se asume que el resto de componentes de la microrred presentan una activaci´ on y desactivaci´ on s´ ubitas. Dicha premisa puede ser una aproximaci´ on correcta si el tiempo de muestreo de la simulaci´ on es bastante m´ as alto que el tiempo de estabilizaci´ on durante la activaci´ on de los componentes reales del sistema. 4. Soluci´ on propuesta Para satisfacer los requisitos expuestos en la Secci´ on 2 para realizar el control y la gesti´ on de la energ´ ıa en la microrred, se ha propuesto un algoritmo basado en el m´ etodo de control predictivo por modelo (MPC). Este algoritmo calcula la l´ ınea de actuaci´ on ´ optima durante una ventana de predicci´ on, y aplica en el modelo simulado la acci´ on de control calculada para el primer instante de tiempo. Este c´ alculo se repite peri´ odicamente, y la acci´ on de control se actualiza conforme a los nuevos c´ alculos. Moliner-Heredia, R. et al. / Revista Iberoamericana de Automática e Informática Industrial 22 (2025) 125-135 137 4.1. Algoritmo MPC propuesto En este subapartado se muestran las ecuaciones del algoritmo MPC propuesto. El algoritmo se ejecuta con una determinada frecuencia fopt |tk−tk−1=1/fopt. Cada vez que se ejecuta, resuelve un problema de optimizaci´ on, cuya extensi´ on en el tiempo est´ a contenida en una ventana de predicci´ on deslizante (PW ={t0, ..., tk, ..., tf}). Nota: Estos instantes de tiempo tkno tienen por qu´ e coincidir con los instantes tmque usa el modelo simulado. No obstante, la frecuencia de actualizaci´ on del modelo simulado debe ser siempre mayor o igual a la frecuencia de optimizaci´ on fopt. Funci´on objetivo La funci´ on objetivo de este problema de optimizaci´ on consiste en minimizar la compra de energ´ ıa desde la red el´ ectrica (requisito 2) a lo largo de toda la ventana de predicci´ on: m´ ın tf X tk=t0 PGRID,PUR(tk),(6) donde PGRID,PUR ∈R+es una variable de decisi´ on del problema de optimizaci´ on que representa la potencia comprada desde la red. Para la resoluci´ on del problema de optimizaci´ on se ha separado la potencia de la red en potencia comprada (procedente de la red) y potencia vendida (hacia la red). La relaci´ on entre ambas potencias se detalla m´ as adelante en el art´ ıculo. La funci´ on objetivo est´ a sujeta a una serie de restricciones, que se muestran a continuaci´ on. En estos casos, se utiliza un doble s´ ımbolo de igual (==) para indicar que se trata de una restricci´ on de igualdad. Restricciones de dep´ositos y bater´ıas: S OC(tk+1) == S OC(tk)−KBAT,DS PBAT,DS (tk) +KBAT,CH PBAT,CH (tk),(7a) MHL(tk)== MHL(tk−1) −KEZ PEZ (tk)+KFC PFC(tk),(7b) S OCmin ≤S OC(tk)≤S OCMAX,(7c) MHLmin ≤MHL(tk)≤MHLMAX,(7d) donde (7a) y (7b) siguen la misma estructura que (2) y (3), respectivamente. Por tanto, S OC ∈R+es una variable de decisi´ on que representa al estado de carga de las bater´ ıas, y MHL ∈R+ es una variable de decisi´ on que representa el estado del dep´ osito de hidr´ ogeno, ambas expresadas en porcentaje sobre su respectivo total de capacidad. Las variables de decisi´ on PBAT,DS , PBAT,CH ,PEZ ,PFC ∈R+representan a la potencia descargada por las bater´ ıas hacia la microrred, la potencia de carga de las bater´ ıas desde la microrred, la potencia almacenada en forma de hidr´ ogeno por el electrolizador y la recuperada por la celda de combustible. Adem´ as, KBAT,DS ,KBAT,CH ,KEZ ,KFC ∈R+son los coeficientes de rendimiento estimados para el problema de optimizaci´ on, para la descarga y carga de la bater´ ıa, as´ ı como para el electrolizador y la celda de combustible. Por ´ ultimo, los l´ ımites de la bater´ ıa y el dep´ osito se expresan con S OCmin, S OCMAX,MHLmin,MHLMAX ∈R+≤100, donde evidentemente S OCmin ≤S OCMAX yMHLmin ≤MHLMAX . Restricciones de equilibrios de potencias: PNET,PRED(tk)== PGEN,PRED(tk)−PDEM,PRED(tk),(8a) PNET,PRED(tk)+PFC(tk)+PBAT,DS (tk)+PGRID,PUR(tk)== PEZ (tk)+PBAT,CH (tk)+PGRID,S AL(tk),(8b) donde PNET,PRED ∈Res una variable auxiliar que engloba las predicciones de la potencia generada predicha PGEN,PRED y la potencia consumida predicha PDEM,PRED, ambas fuente de incertidumbres en este problema de optimizaci´ on. Por otro lado, en (8b) se define el equilibrio de potencias para garantizar el requisito 1. Aqu´ ı se muestra la variable de decisi´ on PGRID,S AL ∈R+, que representa la cantidad de potencia vendida hacia la red; igualmente, el resto de variables de decisi´ on reales (potencias) forman parte de la presente restricci´ on. Restricciones de carga y descarga de las bater´ıas: δBAT,DS (tk)+δBAT,CH (tk)== 1,(9a) PBAT,DS (tk)−PMAX BAT,DS δBAT,DS (tk)≤0,(9b) PBAT,DS (tk)+ϵ1−δBAT,DS (tk)≥0,(9c) PBAT,CH (tk)−PMAX BAT,DS δBAT,CH (tk)≤0,(9d) PBAT,CH (tk)+ϵ1−δBAT,CH (tk)≥0,(9e) donde δBAT,DS ,δBAT,CH ∈Bson variables de decisi´ on booleanas que permiten coordinar y relacionar las variables de decisi´ on PBAT,DS yPBAT,CH . Estas restricciones se han desarrollando aplicando Mixed Linear Dynamics (MLD) para evitar restricciones con productos de variables reales con booleanas, y su desarrollo se muestra m´ as adelante en el art´ ıculo. Por otro lado, PMAX BAT,DS ∈R+es el l´ ımite superior de PBAT,DS . Restricciones de la red: δGRID,S AL(tk)+δGRID,PUR(tk)== 1,(10a) PGRID,S AL(tk)−PMAX GRID,S ALδGRID,S AL(tk)≤0,(10b) PGRID,S AL(tk)+ϵ1−δGRID,S AL(tk)≥0,(10c) δGRID,PUR(tk)−PMAX GRID,PURδGRID,PUR(tk)≤0,(10d) PGRID,PUR(tk)+ϵ1−δGRID,PUR(tk)≥0,(10e) donde δGRID,S AL,δGRID,PUR ∈Bson variables de decisi´ on booleanas que coordinan y relacionan las variables de decisi´ on reales PGRID,S AL yPGRID,PUR. Estas restricciones se han desarrollado usando MLD. Por otro lado, PMAX GRID,S AL,PGRID,PUR ∈R+ son los l´ ımites superiores de PGRID,S AL yPGRID,PUR, respectivamente. 138 Moliner-Heredia, R. et al. / Revista Iberoamericana de Automática e Informática Industrial 22 (2025) 125-135 Restricciones de la unidad de hidr´ogeno verde: PEZ (tk)≤PMAX EZ ,(11a) Pmin EZ δEZ (tk)≤PEZ (tk),(11b) PEZ (tk)+ϵ≤Pmin EZ +δEZ (tk)−Pmin EZ +PMAX EZ +ϵ,(11c) δEZ (tk)ϵ+Pmin EZ ≤PEZ (tk),(11d) PEZ (tk)+δEZ (tk)Pmin EZ ≤PMAX EZ δEZ (tk)+ϵ′,(11e) PFC(tk)≤PMAX FC ,(12a) Pmin FC δFC (tk)≤PFC(tk),(12b) PFC(tk)+ϵ≤Pmin FC +δFC (tk)−Pmin FC +PMAX FC +ϵ,(12c) δFC(tk)ϵ+Pmin FC ≤PFC (tk),(12d) PFC(tk)+δFC (tk)Pmin FC ≤PMAX FC δFC (tk)+ϵ′,(12e) donde de manera similar a las restricciones anteriores, δEZ y δFC ∈Bson variables de decisi´ on booleanas que relacionan las variables de decisi´ on reales PFC yPEZ, y las restricciones se han obtenido aplicando MLD. Por otro lado, PMAX EZ ,Pmin EZ ,PMAX FC ,Pmin FC ∈R+, donde Pmin EZ ,Pmin FC ≥0 son los l´ ımites m´ aximo y m´ ınimo en los que puede encontrarse la potencia consumida por el electrizador o procedente de la celda de combustible cuando est´ an encendidos. Cuando est´ an apagados, PFC yPEZ son 0, valor que puede ser menor que Pmin EZ yPmin FC . La conversi´ on de este umbral inferior de activaci´ on se especifica m´ as adelante. Por ´ ultimo, ϵyϵ′∈R+|ϵ′> ϵ son t´ erminos con una magnitud muy reducida, aplicada para evitar errores num´ ericos. Restricciones para evitar recirculaciones de potencia: δEZ (tk)+δFC(tk)≤1,(13a) δFC(tk)+δBAT,CH(tk)≤1,(13b) δEZ (tk)+δBAT,DS (tk)≤1.(13c) La ecuaci´ on (13a) evita que la unidad de hidr´ ogeno est´ e absorbiendo y vertiendo potencia al mismo tiempo a la microrred, pues debido al relativamente reducido rendimiento de conversi´ on entre potencia el´ ectrica e hidr´ ogeno, simplemente se estar´ ıa desperdiciando energ´ ıa. Es lo que considerar´ ıamos una recirculaci´ on espuria de potencia entre componentes de la microrred sin prop´ osito pr´ actico. Por otro lado, para evitar recirculaciones similares entre la bater´ ıa y el sistema de hidr´ ogeno, con las consiguientes p´ erdidas energ´ eticas, se implementan (13b) y (13c). Las restricciones (13a)-(13c) expresan por tanto la condiciones necesarias para evitar recirculaciones e intercambios energ´ eticos innecesarias de acuerdo al requisito 4. Restricciones de m´ınimo tiempo de activaci´on: tk+MATEZ −1 X m=tk δEZ (m)≥MATEZ [δEZ (tk)−δEZ (tk−1)],(14a) para tk∈ {t0, ..., tf−MATEZ +1}, tk+MATFC −1 X m=tk δFC(m)≥MATFC [δFC(tk)−δFC(tk−1)],(14b) para tk∈ {t0, ..., tf−MATFC +1}, tf X m=tk δEZ (m)−[δEZ (tk)−δEZ (tk−1)]≥0,(14c) para tk∈ {tf−MATEZ +2, ..., tf}, tf X m=tk δFC(m)−[δFC (tk)−δFC(tk−1)]≥0,(14d) para tk∈ {tf−MATFC +2, ..., tf}. Aqu´ ı, MATEZ ,MATFC ∈Nson los tiempos m´ ınimos de activaci´ on del electrolizador y la celda de combustible, respectivamente. Con estas restricciones se limita la cantidad m´ ınima de muestras del optimizador en las cuales el electrolizador y la celda de combustible van a tener que mantenerse activas en el caso de que se requiera su activaci´ on en alg´ un momento. Por tanto, a la hora de decidir la frecuencia de activaci´ on del algoritmo de optimizaci´ on fopt |tk−tk−1=1/fopt, se tienen que tener en cuenta los tiempos m´ ınimos de activaci´ on, dado que debe cumplirse que MAT =n(tk−tk−1), donde n∈N. Como resumen, en la Tabla 1 se muestra una recopilaci´ on de las variables de decisi´ on. Tabla 1: Variables de decisi´ on del problema de optimizaci´ on. Variables reales (∈R+) Variables booleanas (∈B) PGRID,S AL δGRID,S AL PGRID,PUR δGRID,PUR PBAT,CH δBAT,CH PBAT,DS δBAT,DS PEZ δEZ PFC δFC S OC - MHL - 4.2. Los estados de activaci´on y MLD Definici´on de los estados de activaci´on Se define como estado de activaci´ on δ∈Ba toda aquella variable de decisi´ on booleana vinculada a otra variable de decisi´ on real, cuyo objetivo es determinar si dicha variable es distinta de un valor determinado (normalmente 0), para diferenciar entre dos subvariables y determinar cu´ al es la activa, o para coordinar dos variables de decisi´ on cuya activaci´ on es incompatible. En este caso, un estado de activaci´ on est´ a en estado activado si est´ a a 1, y desactivado si est´ a a 0. Para el problema de optimizaci´ on que se ha presentado, se necesita determinar los estados de activaci´ on para las siguientes variables del algoritmo: Moliner-Heredia, R. et al. / Revista Iberoamericana de Automática e Informática Industrial 22 (2025) 125-135 139 Potencia de la red el´ ectrica. Se necesita distinguir entre la potencia comprada procedente de la red (PGRID,PUR) y aquella vendida hacia la red (PGRID,S AL), dado que solo la primera se utiliza en la funci´ on objetivo. Por tanto, hay que definir δGRID,PUR yδGRID,S AL. En este caso, PGRID(tk)=PGRID,PUR(tk)δGRID,PUR(tk) +PGRID,S AL(tk)δGRID,S AL(tk),(15) dado que o bien se est´ a comprando o vendiendo potencia, tal y como se expresa en (10a). Por tanto, hay que considerar la siguiente equivalencia:        [PGRID(tk)≤0]⇔δGRID,S AL(tk)=1, [PGRID(tk)>0]⇔δGRID,PUR(tk)=1.(16) Potencia de las bater´ ıas. La eficiencia de la bater´ ıa depende de si se carga (PBAT,CH ) o se descarga (PBAT,DS ), por lo que se debe definir δBAT,CH yδBAT,DS . En este caso, PBAT (tk)=PBAT,CH(tk)δBAT,CH(tk) +PBAT,DS (tk)δBAT,DS (tk),(17) dado que solo se puede estar o cargando o descargando las bater´ ıas, tal y como se expresa en (9a). La equivalencia queda de la siguiente forma:        [PBAT (tk)≤0]⇔δBAT,CH(tk)=1, [PBAT (tk)>0]⇔δBAT,DS (tk)=1.(18) Elementos de la unidad de hidr´ ogeno verde. Los tiempos m´ ınimos de funcionamiento del electrolizador (PEZ) y de la celda de combustible (PFC ) deben ser restringidos para evitar activaciones y desactivaciones continuas (requisito 3), as´ ı como para evitar recirculaciones energ´ eticas innecesarios (requisito 4). Por tanto, se debe definir δEZ yδFC. Adem´ as, tanto PEZ como PFC tienen un umbral inferior de activaci´ on, por lo que debe tenerse en cuenta a la hora de definir sus estados de activaci´ on:        [0≤PEZ ≤ϵ′]⇔[δEZ =0], hPmin EZ ≤PEZ ≤PMAX EZ i⇔[δEZ =1],(19a)        [0≤PFC ≤ϵ′]⇔[δFC =0], hPmin FC ≤PFC ≤PMAX FC i⇔[δFC =1].(19b) Aplicaci´on del MLD El problema que surge al utilizar los estados de activaci´ on es que las ecuaciones resultantes son funci´ on de un producto de variables de decisi´ on reales y booleanas. Para evitar problemas con los solvers, se han convertido las restricciones l´ ogicas que definen dichas ecuaciones en restricciones de desigualdad mediante el uso de Mixed Logical Dynamics (MLDs), tal y como se explica en (Bemporad and Morari, 1999). Las conversiones que han sido utilizadas se muestran a continuaci´ on: Caso si y solo si (iff). En este caso, la relaci´ on ha⊤x≤0i⇐⇒ [δ=1](20) se sustituye por el siguiente sistema de inecuaciones:        ϵ≤a⊤x−(m−ϵ)δ≤ ∞, −∞ ≤ a⊤x+Mδ≤M.(21) Caso producto mixto de variable decisi´ on real y booleana. En este caso, la relaci´ on z=δa⊤x(22) se sustituye por el siguiente sistema de inecuaciones:                    −∞ ≤ z−Mδ≤0, 0≤z−mδ≤ ∞, −∞ ≤ z−a⊤x+m(1 −δ)≤0, 0≤z−a⊤x+M(1 −δ)≤ ∞. (23) En estas definiciones, x∈Ryδ∈Bson las variables de decisi´ on, a∈Run vector de par´ ametros, y Mymlos valores m´ aximo y m´ ınimo que puede alcanzar a⊤x, teniendo en cuenta los l´ ımites que pueda tener ya definidos xpreviamente. Para el desarrollo de este algoritmo, se han utilizado estas conversiones para transformar las relaciones (16), (18) y (19) en las desigualdades (9) - (12), realizando las operaciones oportunas para tener en cuenta todos los casos considerados. 4.3. Restricciones de m´ınimo tiempo de activaci´on Restricciones dentro del problema de optimizaci´on Para evitar una indeseada activaci´ on y desactivaci´ on repetida de la celda de combustible y el electrolizador, dado que esto reduce su vida ´ util (requisito 3), se han incluido varias restricciones que limitan el m´ ınimo tiempo de activaci´ on de dichos componentes de la microrred. En (Wittmann et al., 2018), se presentan varias restricciones que pueden forzar una determinada cantidad de estados de activaci´ on a 1 si se detecta un flanco de subida. Este m´ ınimo n´ umero de estados activos se denomina m´ ınimo tiempo de activaci´ on (MAT). Estas restricciones se aplican para cada muestra dentro de la ventana de predicci´ on; sin embargo, dependiendo de en qu´ e secci´ on de la ventana de predicci´ on se encuentre la muestra, se aplica una variaci´ on u otra de estas restricciones. La secci´ on principal incluye toda la ventana de predicci´ on excepto las ´ ultimas MAT muestras. Aqu´ ı, en el caso de que haya un flanco de subida, los siguientes MAT estados de activaci´ on se fuerzan a 1, es decir, se activan: tk+MAT−1 X m=tk δ(m)≥MAT [δ(tk)−δ(tk−1)],(24) para tk∈ {t0, ..., tf−MAT +1}. En la secci´ on final, si se detecta un flanco de subida, todos los estados de activaci´ on restantes de la ventana de predicci´ on se fuerzan a 1: tf X m=tk δ(m)−[δ(tk)−δ(tk−1)]≥0,(25) para tk∈ {tf−MAT +2, ..., tf}. 140 Moliner-Heredia, R. et al. / Revista Iberoamericana de Automática e Informática Industrial 22 (2025) 125-135 Transmisi´on de la restricci´on entre ventanas de predicci´on En el apartado anterior se ha definido c´ omo asegurar el tiempo m´ ınimo de activaci´ on dentro del problema de optimizaci´ on. El siguiente paso es transmitir esa informaci´ on entre iteraciones del algoritmo MPC, de manera que en una determinada cantidad de siguientes iteraciones se mantenga el forzado de los estados de activaci´ on. En este caso, al algoritmo se le transmiten los ´ ultimos MAT estados de activaci´ on reales (δR∈B1×MAT ) de un determinado componente del sistema. Los estados de activaci´ on reales pueden diferenciarse de los calculados y decididos por el optimizador por diversos motivos, tales como fallos en el funcionamiento del componente que lo obligan a apagarse temporalmente, o en casos como que se haya ordenado a la celda de combustible o al electrolizador activarse pero el dep´ osito est´ e excesivamente vac´ ıo o lleno como para que sea viable la activaci´ on de estos componentes, y dicha orden equivocada se haya debido a un error en la predicci´ on del consumo o de la generaci´ on. El siguiente paso es la detecci´ on de la cantidad de muestras que lleva activo un componente, es decir, detectar la cantidad de 1s que δRtiene en el extremo final de su vector. Una funci´ on cuenta los ´ ultimos 1s consecutivos (n1), y dependiendo del resultado, el algoritmo procede de tres maneras distintas: El ´ ultimo valor de δRes 0. No hay ninguna activaci´ on que transmitir, por lo que el algoritmo procede directamente con (24). Todos los valores de δRson 1s (n1=MAT). Esto significa que el componente ha estado ya activo durante el tiempo m´ ınimo de activaci´ on, por lo que el algoritmo procede directamente con (24). El ´ ultimo valor de δRno es 0, se detecta una cantidad de 1s consecutivos, pero δRno es un vector de 1s (n1<MAT). En este caso, el algoritmo de optimizaci´ on fuerza primero los siguientes MAT −n1 estados de activaci´ on del optimizador a 1. Posteriormente, se contin´ ua con (24) hasta llegar al final de la secci´ on principal. De esta manera, en cada iteraci´ on del optimizador, n1 ir´ a reduciendo su valor hasta llegar al segundo caso. Por tanto, para evitar problemas con los intervalos de la secci´ on principal, se recomienda que la ventana de predicci´ on PW sea mayor que 2 ·MAT. Este m´ etodo de propagaci´ on de la informaci´ on del m´ ınimo tiempo de activaci´ on permite que en el caso de fallos de los componentes u ´ ordenes err´ oneas, el algoritmo siga funcionando, pues en la siguiente iteraci´ on recibir´ a un 0 como ´ ultimo valor del estado de activaci´ on real y seguir´ a realizando los c´ alculos sin detectar incompatibilidades entre el estado de activaci´ on real y el ordenado. Por tanto, con esto se hace robusto frente a errores en la predicci´ on y fallos s´ ubitos del sistema (requisito 5). Al mismo tiempo, en condiciones de funcionamiento normales estas restricciones evitar´ an el encendido y apagado continuo de los componentes de la microrred (requisito 3). 5. Casos de estudio Para validar el algoritmo propuesto, se han realizado los siguientes experimentos: Tabla 2: Par´ ametros de simulaci´ on usados para los experimentos. Par´ ametro Valor Unidades S OCini 50 ( %) S OCmin 10 ( %) S OCMAX 90 ( %) MHLini 60 ( %) MHLmin 10 ( %) MHLMAX 90 ( %) KBAT,CH 2,4·10−6kW−1 KBAT,DS 2,45 ·10−6kW−1 KEZ 7,5·10−7kW−1 KFC 1,125 ·10−6kW−1 KBAT,CHm8·10−7kW−1 KBAT,DS m8,167 ·10−7kW−1 KEZm2,5·10−7kW−1 KFCm3,75 ·10−7kW−1 Pmin EZ 1.0 kW PMAX EZ 3.0 kW Pmin FC 0.5 kW PMAX FC 3.0 kW PMAX BAT,DS 3.0 kW PMAX BAT,CH 3.0 kW PMAX GRID,PUR 50.0 kW PMAX GRID,S AL 50.0 kW Experimento 1: Comparaci´on con un caso simple. En este experimento se ha realizado una simulaci´on del algoritmo propuesto aplic´andolo al modelo descrito en la subsecci´on 3.2. Se realiza una simulaci´on de 2 d´ıas (48 horas), en el que el tiempo de muestreo de la simulaci´on es de 5 minutos, y la aplicaci´on del algoritmo MPC se realiza cada 15 minutos. La ventana de predicci´on del algoritmo (PW) es de 12 horas, y los tiempos m´ınimos de activaci´on del electrolizador y del combustible son de 1 y 2 horas, respectivamente. Se asume que el algoritmo tiene un conocimiento perfecto del modelo, por tanto, los par´ametros usados por el algoritmo y el modelo simulado son los mismos. ´Estos par´ametros y los valores de inicializaci´on del estado de las bater´ıas (S OCini) y del dep´osito de hidr´ogeno verde (MHLini) se muestran en la Tabla 2. La evoluci´on temporal de la generaci´on fotovoltaica y de consumo se observan en las figuras de la Secci´on 6. En este caso, se asume una predicci´on perfecta de la potencia neta. Usando los mismos datos de partida, en el primer experimento se compara este algoritmo con uno m´as simple, usado en equipos dom´esticos de gesti´on de bater´ıas. ´Este consiste en que cuando PNET > 0, entonces se cargan las bater´ıas, y cuando ´estas est´an llenas, se activa el electrolizador. Si se recibe m´as energ´ıa, entonces se vierte a la red. De manera similar con PNET < 0, que entonces se descargan primero las bater´ıas, y tras eso, se activa la celda de combustible. En ´ultimo caso, se compra potencia desde la red. Las condiciones de tiempo de activaci´on del electrolizador y de la celda de combustible ser´an iguales para ambos experimentos. Moliner-Heredia, R. et al. / Revista Iberoamericana de Automática e Informática Industrial 22 (2025) 125-135 141 Se repiten las simulaciones utilizando distintas evoluciones temporales (casos) de potencia neta de partida y se comparan ambos algoritmos utilizando como ´ındice el sumatorio de la cantidad de potencia comprada desde la red (PPGRID,PUR). Experimento 2: Comparaci´on con un algoritmo sin restricciones de m´ınimo tiempo de activaci´on (validaci´on del requisito 3). En este experimento se probara´ un al-goritmo MPC similar al propuesto en el Experimento 1, pero sin utilizar las restricciones de m´ınimo tiempo de activaci´on de la celda de combustible y del electrolizador, es decir, sin las restricciones (14). Se realizara´ una comparaci´on del comportamiento del algoritmo MPC com-pleto con el algoritmo MPC sin restricciones de m´ınimo tiempo de activaci´on. Los parámetros utilizados son los que se muestran en la Tabla 2. Experimento 3: Comparaci´on con un algoritmo sin restricciones de recirculaciones innecesarias (validaci´on del requisito 4). En este experimento se realizara´ una simulaci´on de un algoritmo MPC similar al del Experimento 1, pero en este caso sin utilizar las restricciones de recirculaciones innecesarias, mostradas en (13), que impiden que el dep´osito de hidr´ogeno verde se llene y vac´ıe de manera simult´anea o que haya intercambios innecesarios de energ´ıa entre la unidad de hidr´ogeno verde y las baterías. Posteriormente se realizara´ una comparaci´on del comportamiento del algoritmo MPC completo propuesto con este ´ultimo algoritmo MPC sin las restricciones de bucle. En este experimento, los parámetros utilizados son tambi´en los de la Tabla 2. Experimento 4: Comparaci´on teniendo en cuenta la incertidumbre del modelo (validaci´on del requisito 5). En este experimento se realizara´ una simulaci´on del algoritmo MPC del Experimento 1, pero se le a˜nadira´ incertidumbre al modelo aplicando una variaci´on en los par´ametros del sistema que se le introducen al optimizador, de manera que no coincidan con los del modelo. En este caso, KBAT,CH y KEZ aumentan un 10 %, y KBAT,DS y KFC se reducen un 10 %. El objetivo de este experimento sera´ comparar el efecto de la incertidumbre en el modelo con respecto al caso con conocimiento perfecto del modelo, centr´andose en ver cu´anto afecta una predicci´on imperfecta al comportamiento de la predicci´on y de la funci´on de coste. 6. Resultados y an´ alisis Los resultados de los experimentos se muestran a continuaci´ on: Resultados del Experimento 1. Para comparar la validez del algoritmo MPC propuesto compar´ andolo con el algoritmo simple, se han realizado simulaciones para 11 casos distintos de 2 d´ ıas (48 horas) de duraci´ on, cada uno con su propia curva de potencia neta. Los resultados de la funci´ on objetivo para cada algoritmo y caso se muestran en la Tabla 3. Se puede observar como la potencia comprada a la red por el algoritmo propuesto es menor en todos los casos. Por otro lado, para mostrar con detalle el comportamiento de la potencia intercambiada por los componentes de la microrred al utilizar ambos algoritmos se muestran respectivas gr´ aficas correspondientes al Caso 1. En las gr´ aficas, la potencia es positiva si ´ esta se vierte al sistema de la microrred, y negativa cuando se extrae potencia. En la Figura 2 se puede observar el comportamiento del sistema tras aplicar el algoritmo MPC propuesto. En este caso, se puede observar c´ omo la bater´ ıa act´ ua principalmente vertiendo potencia a la microrred cuando la potencia neta es negativa (es decir, que la demanda supera a la generaci´ on), puntualmente apoyada por la celda de combustible (hacia las 42 horas aproximadamente en la gr´ afica). Este apoyo es menor debido al mayor rendimiento de la bater´ ıa sobre la unidad de hidr´ ogeno, as´ ı como a las restricciones de tiempo m´ ınimo de activaci´ on, que impiden que se est´ e activando y desactivando constantemente. Por otro lado, se puede observar que en los tramos en los que la potencia neta es positiva (por lo que la generaci´ on supera a la demanda) tanto la bater´ ıa como el electrolizador absorben energ´ ıa para recargar, siempre dentro de sus l´ ımites, la bater´ ıa y el dep´ osito de hidr´ ogeno (aproximadamente a las 10 horas en la gr´ afica). Adem´ as, parte de la potencia sobrante se vierte (vende) hacia la red. El comportamiento del algoritmo frente al pico de potencia de las 10 horas con respecto al de las 37 horas difiere debido a que la funci´ on objetivo minimiza la potencia comprada respecto a la red, y no concreta ni da prioridad a la venta de potencia ni a la carga de unos elementos sobre otros. Por tanto, ambos comportamientos son v´ alidos seg´ un las premisas consideradas. En la Figura 3 se puede observar el comportamiento del sistema tras aplicar el algoritmo simple explicado en la secci´ on anterior. En este caso, se puede observar como la bater´ ıa funciona como compensaci´ on de potencia cuando la potencia neta es negativa, y viceversa, se carga cuando la potencia neta es positiva. Sin embargo, al funcionar de manera reactiva a los cambios del proceso, aparece escalonada, y tiene tendencia a comprar potencia de la red como compensaci´ on. El electrolizador y la celda de combustible no se activan, tal y como estaba definido en su algoritmo, porque no funcionan hasta que la bater´ ıa no sature. En el caso de que saturase, tanto el electrolizador como la celda de combustible se podr´ ıan activar y desactivar de manera continua, pues no tienen ninguna restricci´ on de tiempo que los limite. Resultados del Experimento 2. Para comparar la aplicaci´ on de las restricciones de m´ ınimo tiempo de activaci´ on del algoritmo MPC se han realizado 11 simulaciones de 48 horas cada una, como en el experimento anterior, usando las mismas curvas de potencia neta para cada caso, del algoritmo sin las restricciones de tiempo m´ ınimo de activaci´ on. Tal y como se puede observar en la Tabla 3, el sumatorio de potencia media comprada de la red PPGRID,PUR del algoritmo MPC sin m´ ınimo tiempo de activaci´ on es de magnitud similar al del algoritmo MPC propuesto. Esto es comprensible, dado que se trata de algoritmos en t´ erminos similares. 142 Moliner-Heredia, R. et al. / Revista Iberoamericana de Automática e Informática Industrial 22 (2025) 125-135 0510 15 20 25 30 35 40 45 −5 0 5 Tiempo (h) Potencia (kW) PBAT PEZ PFC PGRID PNET Figura 2: Comportamiento del algoritmo MPC propuesto en el Caso 1. 0510 15 20 25 30 35 40 45 −5 0 5 Tiempo (h) Potencia (kW) PBAT PEZ PFC PGRID PNET Figura 3: Comportamiento del algoritmo simple (Experimento 1) en el Caso 1. 0510 15 20 25 30 35 40 45 −5 0 5 Tiempo (h) Potencia (kW) PBAT PEZ PFC PGRID PNET Figura 4: Comportamiento del algoritmo MPC sin restricciones de m´ ınimo tiempo de activaci´ on (Experimento 2) en el Caso 1. Moliner-Heredia, R. et al. / Revista Iberoamericana de Automática e Informática Industrial 22 (2025) 125-135 143