Operación óptima de redes inteligentes con gran penetración de generación renovable
Abstract
En la actualidad, el sector energético está experimentando importantes cambios. La liberalización de los mercados de energía, los intentos de dar respuesta a los problemas asociados al cambio climático y los avances tecnológicos son solo parte de los factores que lo hacen evolucionar. Pueden distinguirse tendencias comunes en el sector: la reducción del uso de combustibles fósiles y la descabornización, la electrificación de los consumos, el aumento de producción de energías renovables, la aplicación de nuevas tecnologías como la inteligencia artificial, la creciente sofisticación del perfil de los consumidores y el desarrollo de sistemas almacenamiento y de vehículos eléctricos son solo algunos ejemplos. El sector energético es reconocido globalmente como un elemento esencial en la búsqueda de un futuro más sostenible. Una forma que este tiene de dar respuesta a los desafíos que se plantean actualmente es el desarrollo y la implementación de redes de distribución de electricidad inteligentes, las conocidas como smart grids. El conjunto de tecnologías que las forman puede ayudar a optimizar los recursos reduciendo los costes asociados, a la vez que se consigue integrar eficazmente los nuevos sistemas energéticos tales como la generación descentralizada o los sistemas de almacenamiento. En el presente proyecto se ha realizado un programa informático basado en Python que recrea el comportamiento de una red de distribución de media tensión. Se ha creado un modelo de una red de referencia en Europa que incluye diversas cargas y generadores, sistemas de producción eléctrica renovables y sistemas de almacenamiento de energía, entre otros elementos. Igualmente, se ha configurado un problema de optimización matemática MINLP que actúa como sistema de gestión de energía para resolver el problema de Asignación de Unidades y de Despacho Económico en la red. El objetivo de este programa es hallar para un periodo de 24 horas los valores óptimos de las variables de la red de forma que se reduzcan los costes y se mantengan las condiciones operacionales deseadas. Se incluye un modelo estocástico para gestionar las incertidumbres asociadas a los precios del mercado eléctrico y a la generación renovable. El modelo se ha programado en Pyomo y se utiliza el solver BONMIN (al que se accede a través del servidor NEOS) para resolver el problema de optimización. Los diversos casos estudiados presentan las funcionalidades del programa. El modelo permite las desconexión de las cargas industriales para adaptar la demanda e incluye la posibilidad de establecer niveles de reserva de potencia y de limitar la potencia importada por la red. Un caso en el que existe una gran producción de energías renovables es igualmente estudiado. En todos estos casos el programa ha hallado soluciones óptimas y las variables siguen el comportamiento esperado. La solución del problema estocástico asociado se ha realizado mediante una aproximación del problema relajado NLP. A pesar de ello, la solución es similar a la del problema determinista equivalente. El modelo podría ampliarse en el futuro para incluir un enlace entre las dos subredes que componen la red o simular la presencia de vehículos eléctricos, entre otras posibilidades.
Full text
Proyecto Fin de Carrera Ingeniería de Telecomunicación Formato de Publicación de la Escuela Técnica Superior de Ingeniería Autor: F. Javier Payán Somet Tutor: Juan José Murillo Fuentes Dep. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2013 Proyecto Fin de Máster Máster en Ingeniería Industrial Operación Óptima de Redes Inteligentes con gran penetración de Generación Renovable Autor: Alfonso Aycart Lazo Tutor: Alejandro Marano Marcolini Dpto. Ingeniería Eléctrica Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2019
Proyecto Fin de Máster Máster en Ingeniería Industrial Operación Óptima de Redes Inteligentes con gran penetración de Generación Renovable Autor: Alfonso Aycart Lazo Tutor: Alejandro Marano Marcolini Profesor Contratado Doctor Dpto. Ingeniería Eléctrica Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2019
Proyecto Fin de Máster: Operación Óptima de Redes Inteligentes con gran penetración de Generación Renovable Autor: Alfonso Aycart Lazo Tutor: Alejandro Marano Marcolini El tribunal nombrado para juzgar el trabajo arriba indicado, compuesto por los siguientes profesores: Presidente: Vocal/es: Secretario: acuerdan otorgarle la calificación de: El Secretario del Tribunal Fecha:
Agradecimientos Me gustaría comenzar agradeciéndole a Alejandro Marano, mi tutor en este Trabajo de Fin de Máster, su disponibilidad total y su ayuda para conseguir llevar a buen puerto este proyecto. Ha sabido guiarme incluso a distancia y le estoy muy agradecido por ello. Igualmente, me gustaría agradecer a mis supervisores en Energías De Portugal por ver siempre con buenos ojos que realizase este proyecto al mismo tiempo que desarrollaba mi actividad profesional allí y por las facilidades que me han puesto a la hora de volver a Sevilla para presentarlo. Se lo agradezco particularmente a Mariana Gonçalves. Por supuesto, mis compañeros y amigos del Trainee Program de EDP merecen una mención especial por hacerme más llevaderos estos últimos meses de trabajo con su presencia. Por último, y como no podría ser de otra manera, se lo agradezco enormemente a mis padres. Me han apoyado a lo largo de mis años de universidad y sin ellos nunca habría podido tener tantas experiencias de estudio en el extranjero. Si no fuese por su constante insistencia, seguramente este proyecto seguiría sin estar terminado. Gracias. Alfonso Aycart Lazo Estudiante de Máster de Ingeniería Industrial Sevilla, 2019 I
Resumen En la actualidad, el sector energético está experimentando importantes cambios. La liberalización de los mercados de energía, los intentos de dar respuesta a los problemas asociados al cambio climático y los avances tecnológicos son solo parte de los factores que lo hacen evolucionar. Pueden distinguirse tendencias comunes en el sector: la reducción del uso de combustibles fósiles y la descabornización, la electrificación de los consumos, el aumento de producción de energías renovables, la aplicación de nuevas tecnologías como la inteligencia artificial, la creciente sofisticación del perfil de los consumidores y el desarrollo de sistemas almacenamiento y de vehículos eléctricos son solo algunos ejemplos. El sector energético es reconocido globalmente como un elemento esencial en la búsqueda de un futuro más sostenible. Una forma que este tiene de dar respuesta a los desafíos que se plantean actualmente es el desarrollo y la implementación de redes de distribución de electricidad inteligentes, las conocidas como smart grids. El conjunto de tecnologías que las forman puede ayudar a optimizar los recursos reduciendo los costes asociados, a la vez que se consigue integrar eficazmente los nuevos sistemas energéticos tales como la generación descentralizada o los sistemas de almacenamiento. En el presente proyecto se ha realizado un programa informático basado en Python que recrea el comportamiento de una red de distribución de media tensión. Se ha creado un modelo de una red de referencia en Europa que incluye diversas cargas y generadores, sistemas de producción eléctrica renovables y sistemas de almacenamiento de energía, entre otros elementos. Igualmente, se ha configurado un problema de optimización matemática MINLP que actúa como sistema de gestión de energía para resolver el problema de Asignación de Unidades y de Despacho Económico en la red. El objetivo de este programa es hallar para un periodo de 24 horas los valores óptimos de las variables de la red de forma que se reduzcan los costes y se mantengan las condiciones operacionales deseadas. Se incluye un modelo estocástico para gestionar las incertidumbres asociadas a los precios del mercado eléctrico y a la generación renovable. El modelo se ha programado en Pyomo y se utiliza el solver BONMIN (al que se accede a través del servidor NEOS) para resolver el problema de optimización. Los diversos casos estudiados presentan las funcionalidades del programa. El modelo permite las desconexión de las cargas industriales para adaptar la demanda e incluye la posibilidad de establecer niveles de reserva de potencia y de limitar la potencia importada por la red. Un caso en el que existe una gran producción de energías renovables es igualmente estudiado. En todos estos casos el programa ha hallado soluciones óptimas y las variables siguen el comportamiento esperado. La solución del problema estocástico asociado se ha realizado mediante una aproximación del problema relajado NLP. A pesar de ello, la solución es similar a la del problema determinista equivalente. El modelo podría ampliarse en el futuro para incluir un enlace entre las dos subredes que componen la red o simular la presencia de vehículos eléctricos, entre otras posibilidades. III
XÍndice Parámetros 27 Variables 28 Función objetivo 28 Restricciones 28 Restricciones auxiliares 30 Formulación completa 32 Problema estocástico 32 3.2.3 Recapitulación 34 3.3 Estructura del Programa 34 3.4 Importación de Datos 36 3.4.1 Datos deterministas 36 Datos de la red 37 Costes 37 Generación renovable 37 Baterías 38 3.5 Estructura y Datos del Problema Estocástico 39 3.5.1 Escenarios de precios 39 3.5.2 Escenarios de generación eólica 39 3.5.3 Escenarios de generación fotovoltaica 41 3.5.4 Estructura del problema estocástico 41 4 Resultados del Modelo 43 4.1 Caso Base 43 4.2 Limitación de Potencia Importada 45 4.3 Limitación de Potencia Importada y Coste de Desconexión Reducido 47 4.4 Alta Generación Renovable 48 4.5 Restricción de Reservas 49 4.6 Resultados Estocásticos 51 5 Conclusiones y Líneas Futuras 53 Apéndice A Modelo Determinista 55 Apéndice B Estructura del problema estocástico 83 Índice de Figuras 95 Índice de Tablas 97 Índice de Códigos 99 Bibliografía 101
Notación ◦CGrados Celsius CO2Dióxido de carbono L Longitud R Resistencia eléctrica X Reactancia Xsh Reactancia shunt C Capacidad eléctrica n Prefijo que indica un factor de 10−9 k Prefijo que indica un factor de 103 M Prefijo que indica un factor de 106 G Prefijo que indica un factor de 109 T Prefijo que indica un factor de 1012 F Faradios h Horas Hz Hercios m Metro ΩOhmios p.u. Sistema por unidad V Voltio VA Voltiamperio W Vatio Wh Vatio hora s.a. Sujeto a RCuerpo de los números reales CCuerpo de los números complejos ZCuerpo de los números enteros sen Función seno cos Función coseno ≤Menor o igual ≥Mayor o igual ∆Incremento ∈Pertenece a ∀Para todo ∑n i=1Sumatorio para todo i desde 1 hasta n cos ϕFactor de potencia cTVector de costes traspuesto XI
XII Notación tÍndice de las horas iÍndice de los nodos gÍndice de los generadores pv Índice de los generadores fotovoltaicos wÍndices de los aerogeneradores bs Índices de las baterías ni Índices de las cargas residenciales ind Índices de las cargas industriales CA gCoste fijo del generador g CB gCoste proporcional a la potencia del generador g Cind Coste de desconexión de la industria ind Pmax Potencia activa máxima Qmax Potencia reactiva máxima Smax Potencia aparente máxima pg,tPotencia activa del generador g en la hora t qg,tPotencia reactiva del generador g en la hora t peind,t Potencia activa máxima posible de la industria ind para la hora t qeind,t Potencia reactiva máxima posible de la industria ind para la hora t mptPotencia total importada por la red en la hora t ηc,bs Rendimiento de carga de la batería bs ηd,bs Rendimiento de descarga de la batería bs socbs,tEstado de carga de la batería bs en la hora t Gi j Susceptancia de la línea conectando los nodos i y j Bi j Conductancia de la línea conectando los nodos i y j urtReserva de subida de potencia a mantener en la hora t drtReserva de bajada de potencia a mantener en la hora t φt Porcentaje máximo de potencia reactiva importada respecto a la potencia activa importada πtPrecio de la electricidad importada en la hora t ug,tEstado del generador g en la hora t lsind,tPorcentaje de conexión de la industria Vi,t∠θi,t Fasor del voltaje en el nodo i para la hora t; amplitud del voltaje V y ánglo θ p∈Ωp Escenario de precios p dentro del conjunto de escenarios de precios µpProbabilidad de ocurrencia del escenario p s∈Ωs Escenario de generación solar s dentro del conjunto de escenarios generación solar µsProbabilidad de ocurrencia del escenario s v∈Ωv Escenario de generación eólica v dentro del conjunto de escenarios generación eólica µvProbabilidad de ocurrencia del escenario v
Acrónimos AMI Infraestructura de Medida Avanzada (Advanced Metering Infrastructure) BB Algortimo de Ramificación y Acotamiento(Branch-and- Bound) BS Sistema de Baterías CCS Captura y Almacenamiento de Carbono (Carbon Capture and Storage) CIGRE Consejo Internacional de Grandes Redes Eléctricas CHP Cogeneración (Combined Heat and Power) CHPD Cogeneración basada en Diesel CHPFC Cogeneración basada en Pila de Combustible DC Corriente Continua DSM Gestión de la Demanda (Demand-Side Management) DSO Operador de la Red de Distribución (Distribution System Operator) ED Problema de Despacho Económico (Economic Dispatch problem) EDP Energías De Portugal EE Eficiencia Energética EMS Sistema de Gestión de Energía (Energy Management System) ETSI Escuela Técnica Superior de Ingeniería FC Pila de combustible (Fuel Cell) IPCC Grupo Intergubernamental de Expertos sobre el Cambio Climático (Intergovernmental Panel on Climate Change) LP Programación Lineal MINLP Programación Entera Mixta No Lineal (Mixed-Integer Non-Linear Programming) MIP Programación Entera Mixta Lineal (Mixed-Integer Programming) NLP Programación No Lineal (Non-Linear Programming) OA Algoritmo de Aproximación Exterior (Outer- Approximation algorithm) OCDE Organización para la Cooperación y el Desarrollo Económicos XIII
XIV Acrónimos PCC Punto de Acoplamiento Común (Point of Common Coupling) PIB Producto Interior Bruto PV Panel Fotovoltaico QP Programación Cuadrática (Quadratic Programming) SG Redes de distribución inteligentes (Smart Grids) TOU Basado en el tiempo de utilización (Time-Of-Use) UC Problema de Asignación de Unidades (Unit Commitment problem) UE Unión Europea
1 Introducción En los últimos años, conceptos como cambio cimático, desarrollo sostenible y ecologismo han ganado relevancia en los debates sociales y políticos, tanto a nivel local como internacional. De igual manera, la situación demográfica actual se ve marcada por el crecimiento de la población mundial y el consecuente aumento del consumo eléctrico, en 2017 la demanda eléctrica global creció un 3,1% (frente al 2,1% de la demanda energética). El aumento del consumo se produce especialmente en los países en vías de desarrollo como China e India, representando estos dos países el 70% de dicho crecimiento [ 18 ]. Esta tendencia se prevé que continúe, con el consumo de China prácticamente triplicando al de la Unión Europea en 2040 de mantenerse las condiciones actuales (ver figura 1.2) [21]. 1950 1960 1970 1980 1990 2000 2010 2020 2030 2040 2050 2060 2070 2080 2090 2100 2345678 9 10 11 12 13 14 15 16 17 Observado Predicción media Rango 90% predicciones Población (miles de millones) Figura 1.1 Predicciones de crecimiento de la población mundial [13]. En este contexto, se suceden numerosos estudios que advierten sobre las consecuencias que puede tener el cambio climático en nuestro planeta e instan a la puesta en marcha de medidas concretas para paliar sus efectos adversos [ 5 ] (ver figura 1.3). Todo ello ha desembocado en una apuesta a nivel internacional por la sostenibilidad y la lucha contra el calentamiento global, como se aprecia en el gran peso que han ganado estos aspectos en los Objetivos de Desarrollo Sostenible [ 12 ][ 46 ] y en la celebración de numerosas cumbres globales en las que se intenta con un éxito moderado establecer acuerdos y protocolos [ 2 ]. Del mismo modo, esta inercia reguladora ocasiona que cada vez más medidas ecológicas sean asimiladas por las sociedades occidentales: la economía circular y 1
2Capítulo 1. Introducción Figura 1.2 Predicciones de crecimiento de la demanda de energía primaria [21]. el reciclaje, el impulso de pequeñas comunidades sostenibles o incluso el cambio en los hábitos de alimentación por razones medioambientales han cobrado relevancia en los últimos años [ 8 ][ 22 ]. Han llamado particularmente la atención las manifestaciones de estudiantes contra el cambio climático que tuvieron lugar el día 15 de mayo de 2019 internacionalmente y que denotan una importante preocupación por parte de estos temas de una gran parte de la sociedad [3]. Sin embargo, estas iniciativas se demuestran insuficientes a día de hoy para cumplir los objetivos marcados con el fin de limitar el ascenso de la temperatura del planeta por encima de 2 ◦ C (o 1,5 ◦C como especifica el último informe del IPCC [5]). Se plantea por lo tanto necesario el desarrollo de los sistemas eléctricos con la inclusión de generación renovable, la descentralización y el auto-consumo, así como la integración en la red de nuevos actores como los vehículos eléctricos. Para optimizar estos recursos será indispensable el uso de la tecnología. Igualmente, la proliferación de sensores y de herramientas como la inteligencia artificial y el análisis de datos masivos supondrán un gran paso en este sentido [ 54 ]. Es en este espacio donde las redes de distribución eléctrica inteligentes se vuelven más relevantes y necesarias. 1.1 Transición Energética Las formas de energía convencional son generalmente finitas, contaminantes (causando incluso problemas de salud [ 16 ], que son aún más notorios en países como China) o generan residuos, además de ocasionar el conocido efecto invernadero que provoca un aumento de la temperatura del planeta. Con el objetivo de reducir estos efectos del cambio climático, numerosos países están desarrollando hojas de ruta a tal efecto. Esta adaptación de los modelos económicos, de consumo y de generación de energía es lo que se denomina como transición energética. Las tendencias que se observan en el panorama energético global actual, y de una forma más acentuada en los países de la OCDE, y que se prevén ganen mayor importancia en el futuro son las siguientes: • Aumento de la electrificación, ocasionado por el aumento de la tecnología tanto en los hogares como en la industria y por la descarbonización de la economía con el objetivo de reducir las emisiones de gases de efecto invernadero. • Inversiones en la generación de energías renovables y acelerado desarrollo tecnológico de las mismas, especialmente de la energía eólica, solar y de las técnicas de almacenamiento de energía.
1.1 Transición Energética 3 Figura 1.3 Riesgos e impactos del cambio climático [36]. •Mayor interés por medidas de eficiencia energética y reducción del consumo. •Digitalización y optimización de los procesos de producción y gestión de la energía. •Aumento de la movilidad eléctrica, especialmente en zonas urbanas. • Sofisticación de los consumidores finales de energía, pivotando hacia un perfil de prosumidor. Estas características están alineadas con lo que se consideran los tres pilares de la transición energética, los cuales ya se reconocían al principio de la década [7]: • Reducción del consumo energético: es necesario actuar en los tres principales sectores consumidores de energía (residencial, industrial y transportes). Deben promoverse medidas de eficiencia energética (equipos más eficientes en la industria y el sector residencial), sistemas de recuperación de energía, el aislamiento térmico de edificios y el cambio hacia la movilidad eléctrica (ver figura 1.4). La promoción de modelos de producción y consumo basados en la economía circular también se muestra crucial a efectos de la reducción y revalorización de los residuos. • Electrificación de la economía, con especial atención al desarrollo de energías renovables y sistemas de almacenamiento de energía. También se plantea necesario invertir en el desarrollo de tecnologías relacionadas con la captura y el almacenamiento de CO2 (CCS), con el fin de
4Capítulo 1. Introducción Figura 1.4 Evolución necesaria del vehículo eléctrico para alcanzar el objetivo de 2 ◦ C según [ 17 ]. superar los desafíos técnicos que existen actualmente y reducir sus costes [ 48 ]. Sin embargo, recientes estudios alertan del limitado alcance que las tecnologías CCS pueden tener en el objetivo de limitar el aumento de la temperatura del planeta [20]. • Desarrollo de las infraestructuras eléctricas y de las interconexiones entre países. A estos efectos, la promoción de las redes de distribución eléctricas inteligentes permitirá optimizar el uso de los recursos energéticos e integrar al consumidor en la cadena. Siguiendo este razonamiento, la Comisión Europea ha establecido una hoja de ruta con la que pretende definir visiones comunes para estados e industria y fijar objetivos específicos en cuanto a mitigación del cambio climático. Así, para el año 2030, se espera conseguir un aumento de la eficiencia energética en la UE del 23,5%, que el 32% de la energía final consumida proceda de fuentes renovables y que las emisiones de CO2 se reduzcan en un 40% respecto a los niveles de 1990 (actualmente se proyecta que este número aumente hasta el 60% en 2050) [ 19 ]. Dichos objetivos deben ser compartidos por los estados miembros, los cuales tienen que implementar medidas concretas para alcanzarlos. Igualmente, se han llevado a cabo acciones concretas a nivel europeo, como el establecimiento de un mercado de carbono internacional que ocasiona un aumento del precio del CO2que hace más atractiva la inversión en tecnologías limpias. Todas estas medidas de desarrollo sostenible ocasionan importantes cambios en la industria y especialmente en el sector energético. Una de las consecuencias relevantes es el desacople que se produce entre el PIB y el consumo de energía en los países de la OCDE. 1.2 Energías Renovables Según lo explicado anteriormente, la generación de energía renovable se contorna como un elemento clave para la transición energética mundial. Estas fuentes de energía son virtualmente inagotables al provenir de medios naturales con capacidad de regeneración, y generalmente tienen un impacto limitado en el medio ambiente. Asimismo, debido a su carácter autóctono, no generan dependencia de potencias extranjeras como suele ocurrir con los combustibles fósiles. Estas fuentes de energía pueden clasificarse según su origen: •Hidráulica : consiste en aprovechar la energía cinética y potencial de las corrientes y saltos de agua para generar electricidad mediante el funcionamiento de turbinas. Tiene un rendimiento energético muy elevado, el precio de la energía es reducido, su vida útil es muy larga y es una fuente de energía limpia. Igualmente, los embalses permiten el almacenamiento de agua para abastecimiento humano y se utilizan como método de regulación del caudal de los ríos. También es posible transformar las grandes centrales hidráulicas en reversibles mediante la
1.2 Energías Renovables 5 implementación de sistemas de bombeo, los cuales pueden ser usados como almacenamiento energético. Sin embargo, esta tecnología puede causar importantes daños ambientales a los ecosistemas, especialmente debido a la construcción de grandes embalses que inundan importantes extensiones de terreno. •Eólica : se aprovecha la energía cinética del viento para mover turbinas y generar electricidad. Estas tecnologías pueden ser instaladas en zonas muy variadas con escaso impacto medioambiental. En los últimos años su coste de instalación se ha reducido al mismo tiempo que aumentaban los rendimientos, lo que las ha llevado a ser competitivas con respecto a las fuentes convencionales [32]. •Solar : la fuente de energía en este caso es la radiación electromagnética del Sol. Esta radiación se captura mediante diversas tecnologías, entre las que se incluyen los paneles fotovoltaicos (energía solar fotovoltaica) y los colectores solar (energía solar térmica). Según el medio utilizado, el resultado puede ser electricidad o energía térmica. La producción a partir de esta fuente es muy predecible, y los precios se han reducido enormemente en la última década debido a los avances tecnológicos y las economías de escala, por lo que han llegado a alcanzar la paridad de red en algunas regiones [ 32 ]. Las expectativas de evolución de esta tecnología han sido por lo tanto constantemente revisadas al alza. 0 500 1000 1500 2000 2500 3000 2015 2020 2025 2030 Capacidad instalada (GW) Año WEO 2006 WEO 2010 WEO 2018 x37 x28 x15 Figura 1.5 Actualización de las previsiones de capacidad instalada de energía solar [21]. •Geotérmica : consiste en extraer el calor natural del interior de la tierra. Existen reservorios de variada temperatura y profundidad, por lo que estas tecnologías se clasifican entre alta, media, baja y muy baja temperatura. Producen mínimos residuos y son una medida de gran ahorra energético. Por otro lado, pueden causar contaminación de aguas y tan solo está disponible en determinados lugares, aunque en aquellos en los que los recursos son de fácil acceso su coste es muy reducido. •Marina : puede considerarse como un tipo de energía hidráulica, aunque la anterior suele estar asociada con el aprovechamiento de los recursos fluviales y esta hace referencia exclusivamente a los mares y océanos. Hay variadas tecnologías aunque en general se encuentran en fase de desarrollo tecnológico (excepto la energía mareomotriz, que técnicamente no difiere mucho de la energía hidráulica convencional). Puede obtenerse energía de las olas, de las mareas, de las corrientes oceánicas, de la diferencia térmica entre las masas de agua a distintas profundidades y del gradiente de salinidad. Estas fuentes de energía muestras unos recursos potenciales muy importantes, siendo la energía térmica oceánica la de mayores recursos con 44,000 TWh/año (aproximadamente el doble del consumo energético mundial por año) [ 38 ].
12 Capítulo 1. Introducción (generalmente conectadas directamente a las redes de distribución a nivel local) se plantea igualmente un cambio en el comportamiento de las redes que es interesante analizar. Por otro lado, el lenguaje de programación Python es un lenguaje muy versátil, de código abierto, con una sintaxis sencilla e intuitiva y permite modelizar sistemas grandes y complejos gracias a sus características de lenguaje orientado a objetos. Concretamente, la modelización se hará con Pyomo (Python Optimization Modeling Objects), una colección de paquetes de Python utilizado para la formulación de problemas de optimización matemática [30]. Se utilizará un red de referencia creada por el Grupo de Trabajo C6.04.02 de CIGRE. El objetivo de este grupo de trabajo fue establecer un red de referencia común para el estudio de la integración de sistemas de generación distribuida [ 40 ]. La red de distribución de media tensión deriva de una red rural alemana que alimenta un pequeño núcleo urbano y las áreas rurales de alrededor. Está conformada por un total de 14 nodos repartidos en dos subredes con un voltaje de 20 kV que son alimentadas por sendos transformadores de 110/20 kV. La red cuenta con fuentes de generación distribuida (tanto renovable como convencional), sistemas de almacenamiento de energía y diversas cargas residenciales e industriales. 1.5.1 Alcance y objetivos del proyecto En el ámbito de este proyecto debe crearse un modelo de la red y de un sistema de gestión de la energía y analizar los resultados derivados de dicho modelo. La primera parte consiste en modelizar todos los atributos de la red en Pyomo para representarla de la forma más veraz posible. Se creará igualmente un problema de optimización matemática (siendo en este caso de tipo mixto y no lineal: MINLP; ver 2.2) que permita minimizar los costes asociados a la red en función de una serie de restricciones. Estas incluyen, entre otras, el control del voltaje en los nodos y de la potencia aparente las líneas, el balance de los flujos de potencia y el mantenimiento de reservas destinadas al control de frecuencia. Este tipo de problema se llama convencionalmente de Asignación de Unidades (Unit Commitment) si se utilizan escenarios estocásticos basados en predicciones para definir el orden en el que las centrales se pondrán en funcionamiento y de Despacho Económico (Economic Dispatch) si los datos son deterministas y se pretende conocer la potencia que cada unidad deberá producir para suplir la demanda. Se puede encontrar una explicación más detallada de los dos tipos de problema en el apartado 3.2.1. Posteriormente se analizarán los datos extraídos del modelo. Se estudiarán tanto los resultados de un modelo determinista como de un modelo estocástico. Para este último se crearán diferentes escenarios de precios y de generación renovable a partir de datos históricos. Se compararan ambos modelos y se analizarán elementos como el uso de baterías, la desconexión de las cargas no esenciales o la influencia de restricciones auxiliares. En los siguientes apartados se explicará la metodología seguida durante el proyecto, incluida una descripción de la herramienta Pyomo, y se presentará en mayor detalle la red estudiada, así como el problema de optimización a resolver y los datos deterministas y estocásticos utilizados. Posteriormente, se incluye un análisis de los resultados y una conclusión del proyecto con posibles ampliaciones a considerar. En los anexos puede encontrarse el código completo creado en Python.
2 Metodología En este apartado se ofrece una descripción del proceso seguido y las herramientas utilizadas durante la realización del proyecto, especialmente Pyomo, así como una explicación de los problemas de optimización matemática. 2.1 Proceso de Investigación y Análisis El proceso seguido para realizar este proyecto comienza con la documentación y la investigación sobre el sujeto en cuestión: las redes inteligentes, su modelización y la integración de renovables en las mismas. Tras este paso, se ha procedido a modelar la red y el problema de optimización, el cual ha sido resuelto con los datos obtenidos. Por último se han analizado los resultados para verificar el correcto funcionamiento del modelo desarrollado. 2.1.1 Documentación La primera etapa del proyecto ha estado claramente marcada por la investigación y la búsqueda de documentos bibliográficos sobre las redes de distribución eléctrica, las Smart Grids, los problemas de asignación de unidades (UC) y de despacho económico (ED) y la modelización y resolución de problemas de optimización matemática. Igualmente, se ha buscado información relativa al cambio climático, al contexto energético mundial, a la proliferación de energías renovables, a los desafíos asociados a ellas y los vectores de innovación que comienzan a surgir en el mercado. Para comprender el concepto de red inteligente y la tecnología asociada se ha recurrido a artículos mencionados anteriormente como [ 26 ] y [ 27 ]. Del mismo modo, se han estudiado artículos que hacen uso de la red de distribución de CIGRE para diversos estudios, especialmente [ 45 ], aunque también [ 23 ] donde se explican los distintos niveles de control presentes en una smart grid. Ha sido también de gran ayuda el artículo [ 35 ] donde se sientan las bases del problema de optimización matemática a resolver, así como el libro [ 47 ] donde se pueden encontrar ejemplos de problemas similares modelizados con el software GAMS. Igualmente ha sido necesario el estudio de la documentación relativa a Pyomo para la programación del código en Python [30] [10]. 2.1.2 Modelización La parte del proyecto dedicada a la modelización ha sido la más importante y laboriosa. Esta puede dividirse en dos bloques: la modelización de la red de referencia y la creación y resolución del problema de optimización. Esta última parte representa la creación de un sistema de gestión de energía en la red (Energy Management System). Inicialmente se creó un modelo reducido de la red y del problema de optimización. Se utilizaron datos deterministas e instrumentales para una única hora y un número reducido de nodos y se limitó 13
14 Capítulo 2. Metodología el número de restricciones del problema. Progresivamente se fueron añadiendo diferentes cargas y generadores para estudiar si su comportamiento era el esperado hasta terminar con un modelo determinista completo en el que ya se tenían en cuenta la totalidad de las restricciones, incluidas las auxiliares. Por último, se creó un modelo estocástico partiendo del determinista. En él se incluyen varios escenarios, cada uno con una probabilidad asociada, según el precio del mercado eléctrico y la generación fotovoltaica y eólica. 2.1.3 Obtención de datos La obtención de los datos necesarios ha sido un proceso que se ha llevado en paralelo a la creación del modelo y al análisis de los resultados preliminares, puesto que ha sido necesario modificarlos continuamente para comprobar la validez de los mismos. Los datos relativos a la red se han extraído del artículo donde se detallan las características de la red de referencia de media tensión en Europa [ 40 ]. La información sobre el precio del mercado eléctrico para todo un año puede encontrarse en la página web del operador del mercado ibérico [ 9 ] y los datos de viento y radiación solar se han tomado de bases de datos abiertas [6]. Para la creación de los escenarios estocásticos relativos al precio y la generación renovable se han utilizado los datos de todos los días laborales de un año completo. Se ha recurrido a la agrupación de los mismos en cuatro grupos diferenciados según sus similitudes en valor y forma de la curva mediante un programa basado en lenguaje R. La probabilidad de cada escenario se ha determinado según el número de días perteneciente a cada grupo. 2.1.4 Análisis de resultados En esta sección se estudiarán los datos producidos mediante la resolución del problema de optimización realizado anteriormente. Se hará una comparación entre los resultados de los modelos determinista y estocástico durante un periodo de 24 horas. Los indicadores a analizar son los siguientes: •Coste total de la red •Voltaje en los nodos •Potencia importada de la red y potencia producida por los generadores y baterías •Estado de los generadores en los distintos momentos del día (encendido o apagado) •Consumo de las cargas no esenciales Se estudiará el efecto de diversos elementos, como la inclusión de restricciones auxiliares (mantener ciertos niveles de reservas de energía o conseguir un voltaje determinado en el Punto de Acoplamiento Común o PCC), la influencia de la desconexión de cargas no esenciales (la industria en este caso) o un aumento considerable de la generación renovable, siempre verificando que se respeten las restricciones esenciales de voltaje y balances de potencia en los nodos. Por otro lado, se presentará el resultado del modelo estocástico que tiene en cuenta la incertidumbres asociadas a la generación renovable y al precio del mercado eléctrico. 2.2 Problemas de Optimización Matemática Una gran parte del proyecto consiste en la definición y resolución de un problema de optimización matemática. La optimización matemática es el conjunto de técnicas que permiten encontrar el valor óptimo de una función de acuerdo con un criterio establecido y de forma que los valores cumplan con unas restricciones definidas. Para un problema de minimización: Dada una función f:A→R, buscar un elemento x0en A tal que f(x0)≤f(x)para todo x en A.
2.2 Problemas de Optimización Matemática 15 Un problema de optimización consta por lo tanto de una función a optimizar (función objetivo) formada por variables (son los valores que deben ser optimizados) y por parámetros (son predefinidos y se mantienen constantes), y por una serie de ecuaciones e inecuaciones que también combinan variables y parámetros y que componen las restricciones. Según el grado de las ecuaciones o inecuaciones que componen las restricciones y la función objetivo, los problemas pueden ser lineales, cuadráticos o no lineares. Si las variables de la función objetivo solo pueden tomar valores enteros, la programación será de tipo entera, y si se incluyen tanto variables enteras como continuas, el tipo de problema será de programación mixta. Se puede encontrar una descripción de los principales problemas de optimización en la tabla 2.1. Tabla 2.1 Problemas de optimización. Tipo de Problema Descripción Matemática LP (Programación Lineal) min cTx s.a.Ax =b x≥0 x∈Rn,c∈Rn,A∈Rm×n,b∈Rm QP (Programación Cuadrática) min cTx+1 2xTQx s.a.Ax =b x≥0 x∈Rn,c∈Rn,A∈Rm×n Q∈Rn×n,b∈Rm NLP (Programación No Lineal) min f (x) s.a.g(x) = 0 h(x)≤0 xL≤x≤xU x∈Rn f:Rn→R g,h:Rn→Rm MIP (Programación lineal Entera Mixta) min cTx+dTy s.a.Ax =b x,y≥0 x∈Zn,y∈Rl,c∈Rn,d∈Rl A∈Rm×n,Q∈Rm×l,b∈Rm MINLP (Programación No Lineal Entera Mixta) min f (x,y) s.a.h(x,y) = 0 g(x,y)≥0 xL≤x≤xU x∈Rn,y∈Zm f:Rn→R g,h:Rn→Rm Tal y como se explica en la sección 3.2.1, el problema matemático correspondiente al presente proyecto es de tipo MINLP. Además de esta clasificación, existen otros criterios para clasificar los problemas. A continuación se presentan algunos de los más comunes:
16 Capítulo 2. Metodología • Programación convexa: son los problemas en los que el conjunto de restricciones es convexo. • Programación estocástica: en estos problemas al menos una de las variables incluidas en la función objetivo o las restricciones son aleatorias. • Programación heurística: se llega a una solución mediante la realización de suposiciones, por lo que estas son aproximadas. Suele aplicarse en problemas demasiado extensos o complejos y los métodos de resolución suelen ser iterativos. • Programación dinámica: consiste en dividir el problema principal en subproblemas con el fin de acelerar la resolución. 2.2.1 Algoritmos de resolución de problemas MINLP Existen diversos algoritmos para resolver problemas mixtos y no lineares, tanto exactos como heurísticos. La naturaleza del problema y sus dimensiones determinarán qué método es más útil en cada caso. Los principales algoritmos están basados en la relajación de las restricciones para hallar la solución del problema con variables enteras relajadas o restricciones lineares. Uno de los métodos más extendidos es el de Ramificación y Acotamiento (Branch-and-Bound). Consiste en crear bifurcaciones de los valores de las variables binarias (algoritmo de árbol) para dividir el problema en subproblemas más abordables y posteriormente ir descartando aquellas ramas que tengan un coste mayor que el de la mejor solución encontrada hasta el momento. En cada nodo de cálculo, el programa tan solo tiene que resolver un problema NLP (o linear si las ecuaciones han sido linealizadas), lo que simplifica el tiempo de computación, y evaluar las ramificaciones que pueden surgir del nodo para descartar las no factibles o que resultan en una solución peor. Conforme se avanza en la resolución del problema el algoritmo reduce el número de ramas activas. Otro algoritmo es la Aproximación Exterior (Outer-Approximatio). Este es un método iterativo en el que inicialmente se establece un valor de las variables binarias, se resuelve el problema continuo NLP y se utilizan los valores obtenidos como límite superior de la solución (si el problema NLP con las variables binarias fijas es factible, su solución corresponde a una solución del MINLP y por lo tanto puede usarse como límite del mismo). Posteriormente se crea un problema maestro de tipo MIP a través de la solución obtenida en el NLP y de la relajación del problema principal, se resuelve y se usa la solución como límite inferior. Se comparan ambos límites con respecto a una tolerancia predefinida y, si no se ha hallado la solución óptima, se repite el proceso: se resuelve el NLP con los valores de las variables binarias obtenidos en la solución del problema maestro para hallar un nuevo límite superior que se incorpora al problema maestro. El algoritmo LP/NLP basado en Ramificación y Acotamiento aplica el principio del método OA anterior a una búsqueda en árbol. En este caso, el problema maestro se define de forma dinámica para reducir el número de ramas activas progresivamente. Inicialmente se resuelven subproblemas LP en los nodos para obtener límites inferiores y luego, si tienen solución factible, se resuelve el NLP, cuyas soluciones se añaden al problema maestro como límites superiores para acotar las soluciones de los problemas LP en los nodos abiertos [ 41 ]. Este método evita la sucesiva resolución del problema maestro MIP, pero puede aumentar el número de NLP a resolver. 2.3 Pyomo Tanto para la modelización de la red como para la definición y resolución del problema matemático, Pyomo es la herramienta principal que se ha utilizado a lo largo del proyecto [ 29 ] [ 30 ]. Es una colección de paquetes de Python utilizado para la modelización de problemas de optimización (Python Optimization Modeling Objects) que permite crear grandes problemas con una sintaxis sencilla, incluir todos los paquetes adicionales de Python y también importar y exportar datos en los formatos habituales.
2.3 Pyomo 17 Problema original Valores iniciales de las variables enteras (Y0) Resolver prob. NLP fijandoYk Resolver prob. maestro MILP Añandir las restricciones al prob. maestro y obtener nuevas Y ZUZD< Tolerancia? Solución óptima: ZU (límite superior) Solución óptima: ZD (límite inferior) No Sí Solución óptima encontrada Problema original Relajar variables enteras Resolver NLP en un nodo de la lista Eliminar nodo ¿Óptimo? ¿Mejor solución que la actual? No No Sí Actualizar la solución actual ¿Lista de nodos vacía? Sí Sí No Solución óptima encontrada Ramificación y Acotamiento (BB) Aproximación Exterior (OA) Figura 2.1 Algoritmos para resolver problemas de optimización MINLP (fuente: elaboración propia). Además de Pyomo, se han utilizado otros programas para diversas acciones, especialmente Microsoft Excel para la introducción de datos a Pyomo y el entorno y lenguaje de programación R [ 42 ] para agrupar ciertos datos en escenarios estocásticos. Igualmente se ha utilizado la biblioteca Matplotlib para generar los gráficos directamente con el programa en Python [31] 2.3.1 Elementos de Pyomo Python es un lenguaje orientado a objetos. Esto quiere decir que usa objetos para manipular datos según su funcionalidad asociada. Los objetos se caracterizan por tres atributos principales: •Estado: conjunto de datos asignados al objeto. •Comportamiento: métodos y funcionalidades con las que puede operar el objeto. •Identidad: identificador del objeto que lo diferencia de otros. En el caso de Pyomo, se utilizan diversos objetos que son los índices (sets), los parámetros, las variables, las restricciones (constraints) y los objetivos [10].
18 Capítulo 2. Metodología Figura 2.2 Resumen de Pyomo. •Parámetros: representan los datos que tienen que ser proporcionados al modelo. • Variables: son los objetos cuyo valor es desconocido inicialmente y que finalmente representan la solución del problema. •Sets: índices con los que se indexan las variables y parámetros. • Constraints: restricciones que tienen que cumplir las variables. Se definen mediante ecuaciones formadas por los parámetros y las variables. • Objetivos: expresiones o ecuaciones definidas mediante los parámetros y las variables y que tienen que ser minimizadas o maximizadas en la medida de lo posible. Por lo tanto, los parámetros y las variables son indexados por los sets y sirven para representar los elementos de la red que se quiere modelizar. Las variables serán optimizadas mediante un solver para que cumplan todos las restricciones especificadas mientras que se intenta minimizar o maximizar el valor final de los objetivos en la medida de lo posible. Minimizar n ∑ j=1 Ci jXj s.a. n ∑ j=1 Ai jXj≥Bi∀i=1...m Xj≥0∀j=1...n En este problema matemático de programación lineal, los parámetros A,ByCestán indexados por los sets iyj. La variable Xtiene que cumplir las restricciones de forma que se minimice la función objetivo. Pyomo tiene la capacidad de definir dos tipos de modelo diferentes: abstractos y concretos. Los modelos abstractos son aquellos en los que se utilizan símbolos para representar los parámetros, cuyos valores no se especifican hasta el momento de calcular la solución. De esta forma se permite inicializar los valores mediante diversos métodos y reutilizar los modelos con mayor facilidad. En los
2.3 Pyomo 19 modelos concretos los valores de los parámetros son definidos directamente durante la construcción del modelo. Los modelos abstractos pueden asociarse con una plantilla vacía que se rellena con valores posteriormente mientras que los concretos sería una plantilla con los datos ya incluidos. Para este proyecto se ha optado por un modelo abstracto, de forma que los datos serán importados desde un archivo Excel, facilitando así la modificación de los mismos. 2.3.2 Solvers Para resolver los problemas de optimización, Pyomo utiliza softwares matemáticos ajenos. A estos programas se les denomina solvers y su función es recibir la descripción del problema (el modelo) y realizar los cálculos necesarios para obtener la solución óptima. El tiempo de cálculo dependerá del tamaño, de la complejidad del problema y de la capacidad de cada solver, por lo que una evaluación del rendimiento de Pyomo en comparación con otros programas no es factible, las diferencias radicarán principalmente en el solver utilizado por cada uno. Pyomo soporta una gran variedad de solvers y cuenta con interfaces específicas para algunos de ellos (ver figura 2.2) [ 10 ]. Por otro lado, también permite trabajar con cualquier solver que tenga una interfaz AMPL, así como resolver problemas matemáticos a través del servidor NEOS (Network-Enabled Optimization System), una aplicación basada en internet, de libre acceso y que cuenta con una gran variedad de solvers matemáticos, tanto comerciales como de código abierto [ 25 ]. A través de estos programas, o mediante la combinación de ellos, pueden resolverse prácticamente todo tipo de problemas de optimización. Sin embargo, la oferta de solvers de código abierto capaces de resolver problemas de optimización MINLP es limitada [ 34 ] [ 53 ]. Linealizar las ecuaciones es una opción plausible, pero a fin de conseguir la mayor precisión en los resultados se evitará esta opción a pesar de aumentar el tiempo de cálculo. Para aquellos problemas con restricciones no convexas, pueden aplicarse algoritmos heurísticos para encontrar soluciones aproximadas. Algunos de los principales solvers que pueden ser utilizados para resolver problemas MINLP son los siguientes: •BARON (Branch-And-Reduced Optimization Navigator): es un programa orientado a resolver problemas con restricciones no convexas de una forma exacta. Es uno de los solvers más avanzados del mercado pero no está disponible en código abierto. Utiliza un método de Ramificación y Acotación optimizado [49]. •BONMIN (Basic Open-Source Nonlinear Mixed INteger programmin): como su nombre indica, este solver es de código abierto desarrollado dentro del proyecto COIN-OR [ 1 ] y permite resolver el tipo de problemas tratados en este proyecto. Permite encontrar soluciones exactas para problemas convexos mediante diferentes algoritmos (BB, OA, LP/NLP o Híbridos) y aplica un método heurístico (BB con posibilidad de generar diferentes puntos iniciales de cálculo aleatoriamente) para encontrar óptimos locales si el problema es no convexo [24]. •COUENNE (Convex Over and Under ENvelopes for Nonlinear Estimation): este programa también es de código abierto y ha sido igualmente desarrollado dentro de COIN-OR. Es utilizado para calcular soluciones exactas de problemas no convexos mediante linealización y algoritmos de árbol (BB) [14]. •DICOPT (DIscrete and Continuous OPTimizer): en este caso, el solver no es de código abierto y solo puede aplicarse a problemas en los que los términos no lineares sean estrictamente continuos. Utiliza algoritmos de aproximación exterior [4]. El solver escogido para este proyecto es BONMIN debido a su carácter de código abierto y a su capacidad para resolver problemas MINLP tanto convexos como no convexos. Este será utilizado a través del servidor NEOS para facilitar su uso y evitar los problemas asociados a la instalación, para la cual se requieren programas adicionales y la compilación mediante Cygwin para sistemas operativos Windows.
20 Capítulo 2. Metodología 2.3.3 Extensiones: PySP para modelización estocástica Pyomo incluye una serie de extensiones para resolver problemas de mayor complejidad o de naturaleza distinta. Entre estas se incluye una extensión dedicada a la optimización dinámica, una para crear redes de problemas y otra para resolver problemas mediante optimización multinivel. La extensión que resulta más relevante para el proyecto es PySP, que está orientada a la modelización y resolución de problemas estocásticos [53]. Para utilizar este programa, se crea un modelo determinista que va a servir como referencia para el estocástico. Posteriormente, se generan los escenarios, usualmente en forma de árbol de nodos en el que se combinan todas las variables estocásticas; cada nodo tendrá una probabilidad asignada respecto al nivel superior. En el nivel más bajo del árbol aparecerán, por lo tanto, todos los escenarios posibles. Por último, se introducen los datos mediante un fichero de texto .dat. En el árbol de escenarios cada nivel representa una etapa de decisión. En ella se calculará un valor para las variables asignadas a la misma a partir de los valores de las que ya han sido calculadas en los niveles inferiores. En la sección 3.5 se hace una explicación más detallada de la estructura del problema estocástico. Esta extensión incluye dos métodos para resolver el problema: resolver directamente la conocida como "forma extensa" del problema, indexando las variables a los escenarios estocásticos e incluyendo probabilidades en la función objetivo (ver 3.2.2), o aplicar el algoritmo Progressive Hedging, que lo resuelve mediante una estrategia de descomposición. Según Watson y Woodruff en [ 52 ] y [ 53 ], este último método es más útil para problemas extensos y puede ser utilizado con efectividad como método heurístico si el problema incluye variables binarias, como es el caso del presente proyecto. El algoritmo Progressive Hedging es un método iterativo en el que se calculan valores medios de las variables para cada nivel según los valores que tomen las variables que se encuentran en los niveles inferiores. Para cada iteración, se compara la solución calculada con la anterior y se detiene el algoritmo cuando la diferencia sea menor que una tolerancia previamente definida [53]. Al contrario que en con el modelo determinista, esta extensión de Pyomo no permite utilizar la interfaz remota con el servidor NEOS, por lo que en este proyecto no se ha tenido acceso a un solver de código abierto capaz de resolver el problema MINLP en el caso estocástico. Para este caso, se ha optado por relajar las variables binarias y resolver el problema NLP resultante con el solver Ipopt y aproximar posteriormente el valor de las variables binarias a 0 ó 1.
3 Caso de Estudio En este apartado va a presentarse en detalle el objeto de estudio del proyecto. Se explicarán las características de la red escogida, se presentará el problema matemático y su formulación completa y se describirá el funcionamiento y la estructura del programa desarrollado en Python. 3.1 Red de Distribución de Referencia La red objeto de estudio de este proyecto es una red de referencia desarrollada por un grupo de trabajo del CIGRE (Consejo Internacional de Grandes Redes Elétricas) [ 40 ]. Está enmarcada dentro de un proyecto cuyo objetivo era crear un sistema de referencia internacional para facilitar la realización de tests, análisis y validación de los nuevos métodos y técnicas que surgen para la integración de la generación distribuida en las redes eléctricas, así como para el estudio de las tecnologías de smart grids. La necesidad de crear dicho sistema de referencia surge de la proliferación de las unidades de distribución descentralizadas, las cuales suelen ser conectadas a las redes de media y baja tensión, causando significativos cambios estructurales en los actuales sistemas de potencia. En el marco de dicho proyecto, se desarrolló una metodología para facilitar la modelización de redes y se crearon redes de referencia para alta, media y baja tensión (tanto para América del Norte donde las redes funcionan a 60 Hz como para Europa y redes de 50 Hz). Estas redes pueden ser utilizadas para diversos fines: estudio del impacto de la generación distribuida en el flujo de potencia de líneas de distribución o en las curvas de voltaje, estudio de la estabilidad de la tensión, estudio de protecciones de las redes o estudio de sistemas de gestión de energía (como es el caso del presente estudio) entre otros. En este proyecto se modelizará la red de media tensión de referencia para Europa. Esta red ya ha sido utilizada en estudios como [23]. La red de referencia modelizada tiene una potencia nominal de 20 kV y funciona a una frecuencia de 50 Hz. Es una representación de una red rural de media tensión alemana que alimenta una pequeña ciudad y sus alrededores. Cuenta con dos subredes alimentadas por sendos transformadores de 20 MVA y 110/20 kV de una subestación cercana [ 45 ]. La red modelizada consta de catorce nodos divididos entre las dos subredes. Ambas ramas se conectan por un interruptor que puede utilizarse para simular las posibles configuración de la red o estudiar el desempeño de diferentes elementos de conexión, como es el caso de [ 23 ], donde se evaluaron los beneficios de implementar conversores AC-DC para regular la transmisión de potencia. En este proyecto se considerará abierto. Hay un único nodo de transición en el que ninguna carga o generador está conectado, los otros trece tienen algún elemento asociado. Puede verse una representación de la red en la figura 3.1. 3.1.1 Líneas de la red Los conductores que conectan los nodos están caracterizados por los parámetros presentados en la tabla 3.1. Las líneas que conectan los nodos de la subred 1 (nodos 1 a 11) son subterráneos y tienen 21
28 Capítulo 3. Caso de Estudio • Las líneas de la red precisan de parámetros para definirlas. Estos derivan de los presentados en la tabla 3.1 y son la susceptancia, la conductancia y la potencia aparente máxima en las líneas: Gi j,Bi j ySmax i j . •Otros parámetros asociados a la red son los límites de voltaje en los nodos Vmin yVmax. • Para representar algunas restricciones auxiliares se necesitan de parámetros como las reservas para el control primario de frecuencia (de subida urt y de bajada drt ) y el porcentaje máximo de potencia reactiva importada permitido en la red φt∈[0,1]. Por otro lado, existen también parámetros que tienen un carácter estocástico pero que van a presentarse como deterministas en este caso general. Son los relativos a la producción de energías renovables y a los precios del mercado eléctrico: •Los precios horarios de la electricidad importada por la red son representados con πt. • La potencia producida por los generadores fotovoltaicos es ppv,t y la de los aerogeneradores pw,t. Variables Las variables son los elementos cuyo valor el programa calculará para encontrar la solución óptima. Son los siguientes: • Las potencias producidas por los generadores son variables ( pg,t , qg,t ), así como su estado (encendido o apagado), que se representa por la variable binaria ug,t. •La potencia reactiva de los generadores renovables también es variable: qpv,tyqw,t. • Como ya se ha mencionado, las demandas de la industria pueden ser reducidas o directamente desconectadas para mantener las condiciones en la red. Esto las convierte en variables ( pind,t , qind,t ) y hace necesario crear una variable que represente el porcentaje de conexión de estas cargas: lsind,t∈[0,1](siendo 0 totalmente desconectada). • El uso de las baterías está definido por la potencias de carga y descarga, potencia aparente y el estado de carga en cada hora: pcbs,t,pdbs,t,qbs,tysocbs,t. • Las variables de estado de la red, el voltaje en los nodos, se representa por el fasor Vi,t∠θi,t y la potencia aparente en las líneas es si j,t . El ángulo de los nodos conectados a los transformadores se establece como referencia a 0◦. Función objetivo La función objetivo representa los costes de la red y la finalidad del programa es minimizarla. Incluye los costes de importar energía a la red, los costes de generación de los generadores no renovables y los costes de desconexión de las cargas industriales para un periodo de 24 horas: ∑ t∈T πt·mpt+∑ g∈G ug,tCA g+CB g·pg,t+∑ ind∈I Cind 1−lsind,t·peind,t!(3.1) En esta ecuación, mpt representa la energía importada por la red. Al estar dividido el problema en periodos de una hora, todas las potencias son equivalentes a la energía consumida durante esa hora. Restricciones Las restricciones del problema constituyen los límites operacionales y de transmisión de los elementos de la red, así como algunos requisitos adicionales que pueden ser utilizados por el DSO para asegurar el correcto funcionamiento o anticiparse a posibles problemas [47] [35].
3.2 Problema de Optimización 29 • Balances de potencia en los nodos: la potencia producida en los generadores e importada por la red debe ser igual a la potencia consumida por las cargas y las pérdidas en las líneas (representadas por los componentes de la matriz de admitancia Gi j y Bi j ). Este balance debe cumplirse tanto para la potencia activa como para la reactiva: ∑ g∈⊂i pg,t+∑ pv∈⊂i ppv,t+∑ w∈⊂i pw,t+∑ bs∈⊂ipdbs,t−pcbs,t−∑ ni∈⊂i pni,t−∑ ind∈⊂i pind,t =∑ j∈N Vi,t·Vj,tGi j ·cos ϕi j,t+Bi j ·sin ϕi j,t∀i,∀t(3.2) ∑ g∈⊂i qg,t+∑ pv∈⊂i qpv,t+∑ w∈⊂i qw,t+∑ bs∈⊂i qbs,t−∑ ni∈⊂i qni,t−∑ ind∈⊂i qind,t =∑ j∈N Vi,t·Vj,tGi j ·sin ϕi j,t−Bi j ·cos ϕi j,t∀i,∀t(3.3) Los ángulos ϕi j,t son calculados como la diferencia entre el desfase de los nodos iyj: ϕi j,t=θi,t−θj,t. Las ecuaciones del balance de potencia se obtienen a partir del cálculo de la potencia compleja ( SB ) y la corriente ( IB ) en una barra en función de la matriz de admitancia ( YB ) y del vector de voltajes de barra (VB) [33]: SB=BB·I∗ B(3.4) Si=Pi+jQi(3.5) IB=YB·VB(3.6) • Límites operacionales de la red: el voltaje en los nodos y la corriente que circula por las líneas tienen limites operacionales. Vmin ≤Vi,t≤Vmax ∀i,∀t(3.7) si j,t≤Smax i j ∀i,∀t(3.8) El voltaje en los nodos de acoplamiento común (1 y 12) será establecido por el transformador de la subestación que los alimenta y se le permitirá un margen de variación más limitado que en el resto de los nodos, un 1%. La potencia aparente en las líneas se calcula a partir de la intensidad que circula por las líneas: |Ii j,t|2=Vi,t·cos θi,t−Vj,t·cos θj,t+Vi,t·sin θi,t−Vj,t·sin θj,tG2 i j +B2 i j(3.9) sabiendo que Ii j =Si j√3 V • Límites de los generadores: los generadores no renovables tienen límites de potencia activa y reactiva. Aquellos que alimentan la red a través de un inversor (todos excepto el generador diésel), pueden producir o consumir potencia reactiva dentro de unos límites establecidos por la potencia aparente máxima de los inversores. Es necesario incluir las variables binarias de los generadores (encendido/apagado) en estas restricciones para que la potencia máxima sea cero cuando su estado es apagado.
30 Capítulo 3. Caso de Estudio ug,t·Pmin g≤pg,t≤ug,t·Pmax g∀g,∀t(3.10) ug,t·Qmin g≤qg,t≤ug,t·Qmax g∀g,∀t(3.11) qg,t2≤Smax g2−pg,t2∀g6=CHPD,∀t(3.12) • Potencia reactiva en los generadores renovables: estos generadores deben cumplir la misma restricción que los anteriores en cuanto a potencia reactiva se refiere. qpv,t2≤Smax pv 2−ppv,t2∀pv,∀t(3.13) qw,t2≤(Smax w)2−pw,t2∀w,∀t(3.14) • Límites y estado de carga de las baterías: las baterías deben respetar los límites de potencia máxima (tanto producida como consumida), los límites del inversor y los límites de carga. Se añade igualmente una ecuación para hacer el balance de la energía almacenada (estado de carga). pdbs,t≤Pmax bs ∀bs,∀t(3.15) pcbs,t≤Pmax bs ∀bs,∀t(3.16) pdbs,t−pcbs,t2+qbs,t2≤Smax bs 2∀bs,∀t(3.17) SOCmin bs ≤socbs,t≤SOCmax bs ∀bs,∀t(3.18) socbs,t=socbs,t−1+ηc·pcbs,t−pdbs,t/ηd∀bs,∀t(3.19) No es necesario incluir una restricción para especificar que una batería no puede cargarse y descargarse para una misma hora t ya que el programa evitará estas situaciones de manera natural al optimizar el problema gracias a la presencia de rendimientos de carga y descarga. • Las cargas industriales pueden ser desconectadas durante algunas horas para disminuir la demanda de la red cuando los precios del mercado eléctrico son demasiado elevados o los generadores no tienen capacidad suficiente para alimentar a todos los consumidores. La posibilidad de desconectarlas hace necesario definir restricciones para calcular la demanda en función del porcentaje de desconexión (variable lsind,t ) y de la potencia máxima de las cargas para la hora t(peind,t). pind,t=lsind,t·peind,t∀ind,∀t(3.20) qind,t=lsind,t·qeind,t∀ind,∀t(3.21) Restricciones auxiliares Además de las condiciones anteriores, que deben cumplirse para asegurar el correcto funcionamiento de la red, pueden identificarse algunas restricciones adicionales que puede ser exigidas por el DSO [ 35 ]. Estas restricciones auxiliares pueden ser contratadas por el operador o exigidas en tiempo real para mantener el balance en la red.
3.2 Problema de Optimización 31 • Potencia reactiva máxima en la red: puede exigirse mantener un porcentaje máximo de potencia reactiva respecto a la activa en la potencia importada por la red. mqt≤φt·mpt∀t(3.22) domde mqt y mpt son respectivamente las potencias reactivas y activas importadas por la red en la hora t. • Reservas para el control primario de frecuencia: esta restricción está relacionada con el primer nivel de control de la smart grid. Son las llamadas reservas rotativas. Cuando sea necesario, los controladores del primer nivel gestionarán los generadores a nivel local aumentando o reduciendo su potencia con el fin de mantener la frecuencia de la red en 50 Hz. Sin embargo, es el control centralizado el que debe garantizar que dichos generadores tienen capacidad suficiente para adaptar su producción. Esto se consigue estableciendo niveles de reserva de potencia activa, tanto de subida ( urt ) como de bajada ( drt ). Una reserva de potencia de subida implica que los generadores o baterías tienen que ser capaces de inyectar dicha potencia a la red. Una reserva de potencia de bajada significa que los generadores deben poder reducir su potencia o las baterías aumentar su demanda en la medida establecida por la reserva: ∑ g∈G ug,tPmax g−pg,t+∑ bs∈B urbs,t≥urt∀t(3.23) ∑ g∈G ug,tpg,t−Pmin g+∑ bs∈B drbs,t≥drt∀t(3.24) donde urbs,t es el incremento de potencia que puede proporcionar la batería bs en la hora t y drbs,tel margen de reducción de potencia o aumento de demanda que tiene. En estas ecuaciones se considera que únicamente los generadores y los sistemas de baterías pueden ser utilizados para control de frecuencia, aunque también podría considerarse añadir la reducción de consumo de las cargas no esenciales. Esto último se corresponde con las reservas rotativas mencionadas en la sección 1.4. El aumento/reducción de potencia que una batería podría proporcionar es calculado según su potencia máxima y la potencia que estén produciendo o consumiendo durante la hora t. urbs,t=Pmax bs −pdbs,t+pcbs,t(3.25) drbs,t=Pmax bs −pcbs,t+pdbs,t(3.26) Si una batería que está descargándose necesita aumentar la potencia inyectada en la red, solo podrá hacerlo con el margen disponible hasta la potencia máxima. Por otro lado, si la batería está cargándose, la batería tiene un margen mayor de incremento de la potencia inyectada, ya que podría dejar de consumir la potencia de carga y pasar a descargarse con la potencia máxima. La figura 3.3 ilustra cómo se calcula la capacidad que tienen un generador produciendo una potencia pt y una batería consumiendo pc para aumentar su producción y asegurar una reserva de subida. Evidentemente, las baterías deben tener suficiente energía almacenada (o margen de carga) para hacer frente a estas variaciones. Sin embargo, hay que tener en cuenta que el control de frecuencia se realizaría únicamente durante escasos minutos en una hora t (se van a considerar 10 minutos). Por lo tanto, la energía almacenada/capacidad libre debe ser mayor que la potencia a producir/consumir dividida entre 6 (equivalente a carga/descargar durante 10 minutos).
32 Capítulo 3. Caso de Estudio Pmax pt Pmax pc Generador Batería urg Pmax urbs p = 0 Figura 3.3 Incremento de potencia disponible para control primario de frecuencia (reserva de subida). Por ejemplo, para las reservas de subida el razonamiento es el siguiente: si la capacidad de aumento de la potencia inyectada, calculada con la ecuación 3.25, dividida entre 6 es mayor que la energía disponible ( urbs,t 6>socbs,t−SOCmin bs ), significa que la batería no tiene suficiente energía para proporcionar la potencia necesaria durante diez minutos, por lo que la capacidad de las baterías de cara al control de frecuencia pasa a ser la energía almacenada urbs,t=socbs,t−SOCmin g. Para implementar esta comprobación en el código de Pyomo es necesario iterar en las soluciones. Primero ha de resolverse el problema de optimización sin tener en cuenta las limitaciones de energía en las baterías. En esta solución se comprueba si alguna de las baterías no tendría suficiente energía o margen de carga para asegurar la reserva durante 10 minutos. Si es el caso, se incluyen restricciones adicionales para tener en cuenta los niveles de carga en las baterías y horas específicas que incumplan la condición. Posteriormente, se comprueba si han surgido nuevos casos problemáticos y se continúa iterando si fuese necesario. Formulación completa En la tabla 3.4 puede encontrarse la formulación del problema completo para el problema determinista. Problema estocástico El problema estocástico tiene la peculiaridad de contar con parámetros de naturaleza aleatoria. Para abordar esta cuestión, se crean diferentes escenarios con los valores más característicos que pueden adoptar los parámetros y se le atribuye una probabilidad de ocurrencia a cada uno. Los valores de las variables en la solución de cada escenario son, por lo tanto, diferentes. Esto implica que para la formulación del problema estocástico es necesario crear un set adicional por cada parámetro estocástico presente en el problema: • Escenarios de precio: p∈Ωp con una probabilidad asociada a cada uno µp que resultan en los precios horarios del mercado eléctrico πp t. • Escenarios de generación fotovoltaica: s∈Ωs con una probabilidad asociada a cada uno µs que resultan en las potencias generadas s pv,t. • Escenarios de generación eólica: v∈Ωv con una probabilidad asociada a cada uno µv que resultan en las potencias generadas pv w,t.
3.2 Problema de Optimización 33 Tabla 3.4 Formulación del problema optimización determinista. Minimize ∑ t∈T πt·mpt+∑ g∈G ug,tCA g+CB g·pg,t+∑ ind∈I Cind 1−lsind,t·peind,t! Sujeto a: ∑ g∈⊂i pg,t+∑ pv∈⊂i ppv,t+∑ w∈⊂i pw,t+∑ bs∈⊂ipdbs,t−pcbs,t−∑ ni∈⊂i pni,t −∑ ind∈⊂i pind,t=∑ j∈N Vi,t·Vj,tGi j ·cos ϕi j,t+Bi j ·sin ϕi j,t∀i,∀t ∑ g∈⊂i qg,t+∑ pv∈⊂i qpv,t+∑ w∈⊂i qw,t+∑ bs∈⊂i qbs,t−∑ ni∈⊂i qni,t−∑ ind∈⊂i qind,t =∑ j∈N Vi,t·Vj,tGi j ·sin ϕi j,t−Bi j ·cos ϕi j,t∀i,∀t Vmin ≤Vi,t≤Vmax ∀i,∀t si j,t≤Smax i j ∀i,∀t ug,t·Pmin g≤pg,t≤ug,t·Pmax g∀g,∀t ug,t·Qmin g≤qg,t≤ug,t·Qmax g∀g,∀t qg,t2≤Smax g2−pg,t2∀g6=CHPD,∀t qpv,t2≤Smax pv 2−ppv,t2∀pv,∀t qw,t2≤(Smax w)2−pw,t2∀w,∀t pdbs,t≤Pmax bs ∀bs,∀t pcbs,t≤Pmax bs ∀bs,∀t pdbs,t−pcbs,t2+qbs,t2≤Smax bs 2∀bs,∀t SOCmin bs ≤socbs,t≤SOCmax bs ∀bs,∀t socbs,t=socbs,t−1+ηc·pcbs,t−pdbs,t/ηd∀bs,∀t pind,t=lsind,t·peind,t∀ind,∀t qind,t=lsind,t·qeind,t∀ind,∀t mqt≤φt·mpt∀t ∑ g∈G ug,tPmax g−pg,t+∑ bs∈B urbs,t≥urt∀t ∑ g∈G ug,tpg,t−Pmin g+∑ bs∈B drbs,t≥drt∀t La principal diferencia con la formulación presentada anteriormente radica en el hecho de que tanto los parámetros estocásticos como las variables, estarán indexados por estos sets. La función objetivo también se ve alterada y se formula de la siguiente forma: ∑ p∈Ωp µp·∑ s∈Ωs µs·∑ v∈Ωv µv·∑ t∈T πp t·mpp,s,v t+∑ g∈G ug,tCA g+CB g·pp,s,v g,t+∑ ind∈I Cind h1−lsp,s,v ind,ti·peind,t! (3.27)
34 Capítulo 3. Caso de Estudio 3.2.3 Recapitulación En este apartado se ha presentado el problema de optimización que será modelizado con Pyomo. Se ha formulado el problema determinista completo, presentado los índices, los parámetros, las variables, la función objetivo y las restricciones operacionales y auxiliares, y se ha introducido la formulación del problema estocástico. Este problema representa un sistema gestión de energía de smart grid cuya finalidad es determinar cuáles son los valores óptimos de generación y consumo de las cargas no esenciales para mantener las condiciones necesarias en la red y se minimizar los costes asociados. Con el fin de resolver el problema de Asignación de Unidades y tratar la naturaleza aleatoria de algunos parámetros, el modelo estocástico puede ser utilizado para determinar qué generadores se pondrán en funcionamiento en cada hora. Una vez hayan desaparecido las incertidumbres, se resolvería el problema de Despacho Económico con el modelo determinista fijando el estado de cada generador según la solución del problema anterior. De esta forma se ajustaría con exactitud la generación y la demanda en la red. 3.3 Estructura del Programa En esta sección se presentará brevemente la estructura del programa creado y se darán ejemplos representativos en los que se aprecie la sintáxis propia de Pyomo. El código completo puede encontrarse en el Apéndice A. En primer lugar hay un bloque dedicado a la declaración de sets, parámetros y variables. Este bloque corresponde a la modelización de la red de distribución en sí. Siempre debe crearse primero un modelo (en este caso abstracto, ya que los valores se inicializarán más adelante) que incluirá a los demás elementos. En el siguiente código se muestra a título informativo la declaración del modelo, de un set Tque va desde 1 hasta el número de horas definido, de un set Ique sirve para indexar respecto a los nodos y de una variable (el voltaje en los nodos V_i) que puede tomar valores reales no negativos, cuyo valor se inicializa a 20 kV y que está indexada mediante los sets anteriores: Código 3.1 Declaración del modelo, sets y variables en Pyomo. # Creación del modelo m = AbstractModel() # Indices i y t de los nodos m.T = RangeSet(1, m.h) m.I = RangeSet(1,m.n) # Voltaje en los nodos m.V_i = Var(m.I, m.T, within=NonNegativeReals, initialize=20) Algunos parámetros o variables son declarados mediante funciones, como es el caso de la potencia producida por los paneles fotovoltaicos (p_pv) que se forma para cada hora ta partir de un factor de producción f_pv multiplicado por la potencia máxima del sistema fotovoltaico Pmax_pv. La siguiente función actúa como sendos bucles for, iterando sucesivamente respecto al número de horas definido y a los generadores fotovoltaicos: Código 3.2 Declaración de parámetros mediante función en Pyomo. # Declaración de la potencia PV
3.3 Estructura del Programa 35 def p_pv_init(m, i, pv, t): return m.f_pv[t]*m.Pmax_pv[i,pv] m.p_pv = Param(m.PV, m.T, within=NonNegativeReals, initialize=p_pv_init) Un segundo bloque está formado por las restricciones y la función objetivo. Representan el problema de optimización propiamente hablando, es decir, el programa de gestión de la energía eléctrica en la red. La declaración de la función objetivo y las restricciones sigue una lógica parecida a la anterior, realizándose a través de funciones. Es importante aclarar que la importación de energía desde los transformadores ha sido modelizada como generadores sin límite y que no están sujetos a las mismas restricciones que estos. Código 3.3 Declaración de función objetivo y restricciones en Pyomo. # Declaración función objetivo def obj_expression(m): return (sum(m.price[t]*(m.p_g[1,"Red",t] + m.p_g[12,"Red",t]) + sum(m.u_g[i,g,t]*(m.CosteB_g[i,g]*m.p_g[i,g,t]+m.CosteA_g[i,g]) for (i,g) in m.G if g != "Red") + sum(m.C_ind[i,ind]*(m.pe_ind[i,ind,t]*(1 - m.ls[i,ind,t])) for (i,ind) in m.IND) for t in m.T)) m.OBJ = Objective(rule=obj_expression) # Limites de voltaje en nodos def Vbounds_const_rule(m, i, t): return m.Vmin <= m.V_i[i,t] <= m.Vmax m.Vbounds_const = Constraint(m.I, m.T, rule=Vbounds_const_rule) El siguiente bloque corresponde a la introducción de los datos. Como se detallará en el apartado 3.4, los parámetros son importados a partir de un archivo Excel para facilitar su manipulación. Esto se consigue a través de un elemento DataPortal. El programa busca en el archivo indicado el valor asignado a un rango de Excel y lo asocia al parámetro o set especificado. En lugar de un único valor, puede ser que el rango incluya múltiples valores, como es el siguiente caso en el cual se inicializa el set ARCS mediante una matriz que incluye información sobre qué nodos están directamente conectados por líneas eléctricas. El hecho de inicializar los valores supone crear un modelo concreto de Pyomo a partir del abstracto. Para ello se crea un instancia (instance) con la que se trabajará a partir de ese momento, puesto que contiene los elementos con sus valores ya asignados: Código 3.4 Asignación de valores a los parámetros de Pyomo a partir de Excel. data = DataPortal() # Importa los parámetros de entrada data.load(filename=’Datos.xlsx’, range=’horas’, param=’h’) data.load(filename=’Datos.xlsx’, range=’nodos’, param=’n’) # Carga los índices de la red data.load(filename=’No_escen.xlsx’, range=’arcos’, format=’set_array’, set=m.ARCS)
36 Capítulo 3. Caso de Estudio instance = m.create_instance(data) Por último, se incluyen las líneas de código dedicadas a la resolución del problema matemático. Como ya se ha indicado, en este caso se resolverá a través de la interfaz NEOS y mediante el solver BONMIN. También se crea un archivo .pickle que guarda la solución encontrada para su posterior uso. Pueden incluirse igualmente líneas de código dedicadas a la creación de un archivo PDF con gráficos y su presentación por pantalla: Código 3.5 Resolución del problema matemático y creación de gráficas en Pyomo. # Para resolver remotamente a través de la interfaz NEOS solver_manager = SolverManagerFactory(’neos’) results = solver_manager.solve(instance, opt=’baron’) # Guarda la solucion en results instance.solutions.store_to(results) # Archiva results: solucion & metadatos de ejecucion file = open("results.pickle", "wb") pickle.dump(results, file) file.close() # Gráfica del perfil de las cargas residenciales fig4 = plt.figure() fig4, f4_ax = plt.subplots() carga = [] for t in horas: carga.append(value(instance.f_house[t])) f4_ax.plot(horas, carga) f4_ax.set_ylim([0,1]) f4_ax.set(xlabel = ’Horas’, ylabel=’Potencia activa [p.u.]’, title=’Perfil de las cargas residenciales’) f4_ax.grid() fig4.savefig("cargas.pdf", bbox_inches=’tight’) plt.show() 3.4 Importación de Datos En esta sección se detallarán datos más específicos que son necesarios para el funcionamiento de la red. Algunos surgen de investigaciones o valores reales habitualmente utilizados y otros son suposiciones. Para facilitar la manipulación de los datos, estos son importados por el programa desde un archivo Excel. 3.4.1 Datos deterministas Los datos deterministas son aquellos que se utilizan como base de la simulación y muchos de ellos también estarán presentes en la simulación del caso estocástico.
3.4 Importación de Datos 37 Datos de la red Además de los datos ya indicados en la sección 3.1, a continuación se muestran ciertos parámetros necesarios para definir el problema determinista: Tabla 3.5 Parámetros de la red [45]. Parámetro Valor Horas 24 Nodos 14 Vbase 20 kV Desviación de voltaje máxima permitida 5% Smax 10 MVA φ30% donde φ es el porcentaje máximo de potencia reactiva importada con respecto a la potencia activa y que equivale a un factor de potencia (cos ϕ) aproximado de 0,95. Costes Los costes de los generadores y de desconexión de la industria han sido supuestos y tienen un carácter instrumental, ya que serán modificados para comprobar el correcto funcionamiento del programa en diferentes casos. A los generadores se les ha asignado un coste fijo instrumental para evitar que sean encendidos cuando no sean necesarios. En este caso, la industria tiene un coste de desconexión igualmente instrumental para evitar que sean desconectadas en condiciones normales de precio del mercado eléctrico. Tabla 3.6 Costes de la red. Elemento Coste fijo (CA) [€/MWh] Coste variable (CB) [€/MWh] CHPD 20 150 FC 10 100 CHPFC 10 100 Industria 200 - El coste de importación de la electricidad ha sido calculado a partir de los datos del mercado eléctrico español para el año 2018 en periodos de 24 horas [ 9 ]. Se han tratado los datos para eliminar los fines de semana y se ha hecho una media de los días restantes. Estos precios han sido incrementados en un 30% para determinar cuál será el coste de importar energía a la red en el modelo determinista (ver figura 3.4). Generación renovable Para hallar el perfil de generación de los aerogeneradores, se ha seguido un procedimiento similar al anterior a partir de los datos presentados en [ 6 ] para el año 2018 en Oporto. Los datos de
44 Capítulo 4. Resultados del Modelo En la figura 4.3a puede apreciarse que toda la energía consumida por la red es importada a través de la subestación o producida por los generadores renovables para evitar poner en marcha los generadores convencionales. La energía importada por la red sigue el perfil del consumo total (residencial e industrial): existe un pico prolongado durante la mañana y un pico mayor por la noche debido al aumento de las cargas residenciales (que son las que tienen mayor peso en esta red). La figura 4.2 sirve de comprobación de que las industrias no se han desconectado al tener estas un coste mayor al de importación de energía. 0 5 10 15 20 25 Horas (h) 0 20 40 60 80 100 Consumo (%) Consumo industrial Figura 4.2 Consumo industrial en el caso base. En la subred 1, para disminuir las pérdidas de carga en las líneas, el voltaje en el nodo de entrada (PCC) se mantiene por encima del voltaje de los nodos intermedios para disminuir las pérdidas de carga y varía entre sus límites superior e inferior (20 ±1 % kV). El voltaje en los nodos varía según la carga y la generación en los nodos; el aumento de potencia en un nodo ocasiona un incremento del voltaje y el aumento de la carga resulta en una caída. Este comportamiento va a mantenerse para todos los casos estudiados excepto cuando existe una importante generación en los nodos intermedios. Por otro lado, la subred 2 muestra una lógica diferente en el voltaje de sus nodos. Al ser una red muy capacitiva y que está formada únicamente por cargas, las tensiones en los nodos intermedios son más elevadas que en el nodo de PCC. Esta inversión se debe a que el carácter capacitivo tiene un efecto similar al de un generador en los nodos intermedios si la red tiene cargas reducidas (esto puede apreciarse también en la subred 1 para la hora 1, ya que las cargas en la red en ese momento son mínimas). El hecho de que el voltaje del nodo 12 se mantenga en su valor mínimo es también consecuencia de las capacidades presentes en la red. Normalmente, la tensión se mantendría en su valor máximo para minimizar las pérdidas de carga pero las elevadas capacidades ocasionan fugas de corriente a tierra, por lo que un voltaje mayor ocasionaría una pérdida mayor de la corriente capacitiva y un mayor coste. En la figura 4.3b se muestra el resultado que se obtendría si se eliminase el carácter capacitivo de las líneas. Las tensiones en los nodos PCC podrían mantenerse en su nivel máximo y ser mayores que las de los nodos intermedios para minimizar pérdidas. Esto ocurre para ambas redes, lo que implica que las horas en que la tensión del nodo 1 está en su valor mínimo es consecuencia del carácter capacitivo de las líneas. Las baterías funcionan siguiendo ambas un mismo patrón, por lo que se representa únicamente la suma de sus potencias y estados de carga. Almacenan energía importada la red cuando el precio es menor, notablemente durante la noche, para descargarse cuando repunta, que suelen corresponderse con los momentos de mayor demanda eléctrica, y reducir el coste de funcionamiento.
4.2 Limitación de Potencia Importada 45 19.7 19.8 19.9 20.0 20.1 20.2 Voltaje (kV) Voltaje en los nodos V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12 V13 V14 5 10 15 20 Energia (MWh) Energia importada por la red Subred nodo 12 Subred nodo 1 0 5 10 15 20 25 Horas (h) 0.0 0.5 1.0 1.5 2.0 2.5 Potencia (MW) Potencia y energia de las baterias (suma) SOC (MWh) Carga Descarga (a) Caso base. 19.85 19.90 19.95 20.00 20.05 20.10 20.15 20.20 Voltaje (kV) Voltaje en los nodos V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12 V13 V14 5 10 15 20 Energia (MWh) Energia importada por la red Subred nodo 12 Subred nodo 1 0 5 10 15 20 25 Horas (h) 0.0 0.5 1.0 1.5 2.0 2.5 Potencia (MW) Potencia y energia de las baterias (suma) SOC (MWh) Carga Descarga (b) Caso base sin capacidades en las líneas. Figura 4.3 Voltajes, potencias y energía en el caso base. 4.2 Limitación de Potencia Importada En este caso de estudio se limita la potencia que puede importar la subred 1 a 16,5 MW para estudiar el comportamiento de los generadores. Los valores del resto de parámetros se mantienen iguales a los del caso base. El coste total de la red es de 44620,22 € con una energía total importada por la red de 571,609 MWh . En comparación con el caso base, el coste ha aumentado mientras que la energía importada se ha reducido debido a que parte de la demanda se ha suministrado a través de los generadores o desconectando las industrias, lo que tiene un coste mayor que importar energía. Durante los periodos de mayor consumo, la potencia importada alcanza su nivel máximo permitido y es necesario encender los generadores para responder a la demanda. El generador que se utiliza es el de mayor capacidad dentro de los de menor coste: el sistema de cogeneración basado en pila de combustible (CHPFC). Igualmente, es necesario que las baterías descarguen potencia a la red durante esos momentos. Para ello, las baterías se cargan durante el día, al contrario de lo que ocurría en el caso base que únicamente lo hacían cuando el precio de electricidad era menor. Esto ocurre en las horas 11, 15, 16 y 17 con el fin de que las baterías tengan suficiente energía en los momentos críticos. Las industrias se desconectan en algunos momentos críticos, especialmente en la hora 19, ya que los generadores por sí solos no serían capaces de responder a la demanda. Para entender por qué no se utilizan todos los generadores a su máxima potencia en lugar de desconectar las industrias en menor medida hay que considerar sus costes fijos: a pesar de que el coste por MWh de desconectar la industria es mayor que el de los generadores, producir energía con ellos tendría un coste mayor para la red. Esto se demuestra con el cálculo realizado a continuación. En la hora 12, las industrias producen al 83% de su capacidad máxima. Dado que en esta hora se desconecta el 99% de la
46 Capítulo 4. Resultados del Modelo 19.7 19.8 19.9 20.0 20.1 20.2 Voltaje (kV) Voltaje en los nodos V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12 V13 V14 6 9 12 15 Energia (MWh) Energia importada por la red Subred nodo 12 Subred nodo 1 0.0 0.5 1.0 1.5 2.0 2.5 Potencia (MW) Potencia y energia de las baterias (suma) SOC (MWh) Carga Descarga 0 5 10 15 20 25 Horas (h) 0.0 0.1 0.2 0.3 0.4 Potencia (MW) Potencia de los generadores CHPFC 9 (a) Importación limitada. 19.7 19.8 19.9 20.0 20.1 20.2 Voltaje (kV) Voltaje en los nodos V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12 V13 V14 6 9 12 15 Energia (MWh) Energia importada por la red Subred nodo 12 Subred nodo 1 0.0 0.5 1.0 1.5 2.0 2.5 Potencia (MW) Potencia y energia de las baterias (suma) SOC (MWh) Carga Descarga 0 5 10 15 20 25 Horas (h) 0.0 0.2 0.4 0.6 Potencia (MW) Potencia de los generadores CHPFC 9 CHPD 9 (b) Importación limitada con coste CHPD reducido. Figura 4.4 Voltajes, potencias y energía con importación limitada. 0 5 10 15 20 25 Horas (h) Desconexion (%) 100 48 99 100 24 100 100 25 100 Desconexion cargas industriales Ind 3 Ind 1 Ind 7 Ind 10 Ind 9 (a) Importación limitada. 0 5 10 15 20 25 Horas (h) Desconexion (%) 100 34 100 24 100 100 Desconexion cargas industriales Ind 3 Ind 1 Ind 7 Ind 10 Ind 9 (b) Importación limitada y coste CHPD reducido. Figura 4.5 Desconexión de la industria con importación limitada. industria en el nodo 10 y 35% de la del 9, el descenso de consumo es de 0,223 MWh. ∆Eind =0,83 ·(0,57375 ·0,35 +0,068 ·0,99) = 0,223MWh Esto supone un coste de 44,51 € desconectando las industrias. Si esa energía se produjese con el generador diésel, el coste sería de 53,45 €: CCHPD =150 ·0,223 +20 =53,45MWh
4.3 Limitación de Potencia Importada y Coste de Desconexión Reducido 47 En la figura 4.4b y 4.5b se muestran los resultados obtenidos si el generador diésel tuviese un coste fijo de 11 € /h en lugar de 20 € /h. La totalidad de la potencia sería producida por el generador CHPD en la hora 11 y se reduciría la desconexión de la industria en la hora 19. 4.3 Limitación de Potencia Importada y Coste de Desconexión Reducido Con el fin de estudiar el correcto funcionamiento de la desconexión de la industria, se han mantenido las restricciones y los valores de los parámetros del caso anterior (red limitada a 16,5 MW, costes de los generadores del caso base y reservas a cero). Sin embargo, en este caso se ha reducido el coste de desconexión de la industria a 100 €/MWh, menos que el coste de los generadores. La energía total importada por la red es exactamente la misma a la del caso anterior, aunque los costes ahora son más reducidos gracias al menor coste de desconexión 44413,39 € Puede apreciarse en la figura 4.7 que vuelven a desconectarse las industrias durante los periodos de máximo consumo de energía, en este caso en mayor medida que en el anterior (existe desconexión también en la hora 10). La potencia que antes era suministrada por los generadores se consigue ahora mediante la desconexión de las industrias. El comportamiento de las baterías sigue la misma lógica que en el caso anterior, ya que, desde su punto de vista, el funcionamiento de la red es exactamente el mismo (en lugar de producir más con generadores, se reduce la demanda con la industria). 19.7 19.8 19.9 20.0 20.1 20.2 Voltaje (kV) Voltaje en los nodos V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12 V13 V14 6 9 12 15 Energia (MWh) Energia importada por la red Subred nodo 12 Subred nodo 1 0 5 10 15 20 25 Horas (h) 0.0 0.5 1.0 1.5 2.0 2.5 Potencia (MW) Potencia y energia de las baterias (suma) SOC (MWh) Carga Descarga Figura 4.6 Voltajes, potencias y energía con coste de desconexión reducido.
48 Capítulo 4. Resultados del Modelo 0 5 10 15 20 25 Horas (h) Desconexion (%) 100 58 99 99 100 100 100 100 14 69 100 Desconexion cargas industriales Ind 3 Ind 1 Ind 7 Ind 10 Ind 9 Figura 4.7 Desconexión de las industrias con coste de desconexión reducido. 4.4 Alta Generación Renovable Para este caso se ha aumentado en gran medida la capacidad máxima de los generadores renovables, de forma que en algunos momentos del día la red produzca más de lo que consume, por lo que pasa a exportar energía. Se han multiplicado por 10 las capacidades de los sistemas fotovoltaicos y por 15 la del aerogenerador. Las potencias producida por los generadores renovables en este caso son las siguientes: 0.00 0.15 0.30 0.45 0.60 Potencia (MW) Potencia de los generadores fotovoltaicos 0 5 10 15 20 25 Horas (h) 0.0 2.5 5.0 7.5 10.0 Potencia (MW) Potencia de los aerogeneradores Figura 4.8 Producción de los generadores renovables aumentada. El coste total de la red es de 30062,33 € (-32,2% respecto al caso base) con una energía total importada por la red de 385,35 MWh (reducción del 32,8%). Puede verse en la figura 4.9 que la curva de energía importada por la subred 1 se ha desplazado verticalmente hacia abajo con respecto a la del caso base, ya que el perfil de generación de los generadores renovables y de las cargas se ha mantenido. De hecho, durante los periodos de menor consumo (horas 1 a 5 y hora 24) la subred 1 pasa de importar a exportar energía, lo que permite vender la energía durante esas horas. Es precisamente en las horas de menor coste que parte de esa energía es utilizada para cargar las baterías. El perfil de carga de estas es muy similar al caso base,
4.5 Restricción de Reservas 49 19.8 20.0 20.2 20.4 20.6 20.8 21.0 Voltaje (kV) Voltaje en los nodos V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12 V13 V14 0 5 10 15 Energia (MWh) Energia importada por la red Subred nodo 12 Subred nodo 1 0 5 10 15 20 25 Horas (h) 0.0 0.5 1.0 1.5 2.0 2.5 Potencia (MW) Potencia y energia de las baterias (suma) SOC (MWh) Carga Descarga Figura 4.9 Voltajes, potencias y energía con alta generación de renovables. descargándose en las horas de mayor coste. El consumo industrial se mantiene en todo momento a su nivel máximo. Al contrario que los otros casos, el voltaje de los primeros nodos de la subred 1 siempre es menor que en los nodos intermedios ya que la producción renovable tiene lugar en estos últimos. Disminuye aún en mayor medida cuando existe una mayor generación en los nodos intermedios, es decir, cuando la producción de energía fotovoltaica se suma a la eólica en las horas centrales del día. La tensión en los nodos de la subred 2 no varía ya que esta red no cuenta con generadores renovables. 4.5 Restricción de Reservas Este último caso es el más complicado pues supone añadir las restricciones de reserva de potencia para control de frecuencia primario y realizar sucesivas iteraciones resolviendo el problema matemático para poder considerar la energía disponible en las baterías (ver apartado 3.2.2). Se van a mantener los valores de los parámetros del caso base y se va a establecer una reserva de 620 kW tanto para subida como parar bajada de potencia. El coste total de la red es de 44442,98 € con una energía total importada por la red de 574,05 MWh . Estos valores no divergen prácticamente nada de los del caso base (+0,2% para el coste) porque los costes derivados de mantener algunos generadores encendidos es muy limitado. Tal y como se aprecia en las figuras 4.10 y ?? , las baterías son las que proveen prácticamente en su totalidad las reservas de subida y bajada, debido a que tienen un menor coste que los generadores. Se limita la carga a alrededor de 200 kW de forma que siempre exista un margen de aumento de
50 Capítulo 4. Resultados del Modelo 19.7 19.8 19.9 20.0 20.1 20.2 Voltaje (kV) Voltaje en los nodos V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12 V13 V14 5 10 15 20 Energia (MWh) Energia importada por la red Subred nodo 12 Subred nodo 1 0.0 0.5 1.0 1.5 Potencia (MW) Potencia y energia de las baterias (suma) SOC (MWh) Carga Descarga 0 5 10 15 20 25 Horas (h) 0.00 0.08 0.16 0.24 Potencia (MW) Potencia de los generadores CHPFC 9 CHPD 9 Figura 4.10 Voltajes, potencias y energía con restricciones de reserva. 620 kW hasta la potencia máxima de carga conjunta (800 kW) para asegurar la reserva de bajada (si fuese necesaria, las baterías pueden aumentar el consumo de la red en ese margen). Ocurre lo mismo para la reserva de subida, en los momentos en los que las baterías descargan potencia a la red esta es limitada a un valor cercano a 200 kW. En la primera hora es necesario encender el generador diésel para hacer frente a la reserva de subida: pese a que la batería del nodo 5 tiene capacidad para aumentar su potencia de descarga en la red en 696 kW (600 kW de potencia máxima y 96 kW que consume de la red para cargarse), tan solo cuenta con 180 kWh almacenados, por lo que nunca sería capaz de producir la potencia necesaria durante 10 minutos (solo podría descargar 90 kWh debido a que debe mantener un mínimo estado de carga de 90 kWh). Lo mismo ocurre en la última hora para la batería del nodo 10: su estado de carga está al mínimo, por lo que es necesario encender el generador CHPFC para asegurar la reserva de subida. En este caso se necesita menos margen (la batería del nodo 10 tiene menor potencia que la del 5) y este se puede proveer con un menor coste gracias al CHPFC. En ambos casos los generadores están produciendo a su mínima potencia. Con respecto al voltaje en los nodos, se mantiene la lógica explicada anteriormente. Al comparar estos valores con los del caso base, se aprecia que la curva de voltaje en los nodos intermedios sigue
4.6 Resultados Estocásticos 51 de una forma más exacta la curva de consumo (o de importación de energía), ya que en este caso la carga o descarga de las baterías es mucho más limitada y no altera mucho el perfil de voltaje. 4.6 Resultados Estocásticos El objetivo del modelo estocástico es resolver el problema de Asignación de Unidades. Se pretende determinar con antelación qué generadores será necesario poner en marcha teniendo en cuenta la incertidumbre asociada a los precios de los mercados eléctricos y a las fuentes de energías renovables. Se ha establecido un límite a la potencia importada por la red de 16,5 MW, igual que en el caso 4.2, para poder estudiar la puesta en marcha de los generadores y se ha mantenido el resto de parámetros con los valores del caso base. Los parámetros estocásticos están representados por los escenarios especificados en la sección 3.5. Como se ha explicado en el apartado 3.5.4, no se ha tenido acceso a ningún solver con capacidad para resolver el problema MINLP, ya que el entorno PySP no permite acceder remotamente a la interfaz NEOS. Por lo tanto, se ha resuelto el problema relajando las variables binarias de estado de los generadores y posteriormente se ha aproximado el valor de estas variables a 0 ó 1 (apagado/encendido). Los valores de estas variables son calculados a partir de los resultados de cada escenario haciendo una media ponderada con las probabilidades asociadas a cada uno. El resultado se muestra en la figura 4.11. El coste máximo para el caso estocástico se ha dado en el escenario 134 (escenario 1 de precios, 3 de generación fotovoltaica y 4 de generación eólica). Este escenario tiene el mayor coste promedio de electricidad (85,9 € /MWh) y las menores producciones medias renovables (0,042 y 0,088 para PV y eólica respectivamente). Esto supone un coste total de 51941,5 € . El coste mínimo de 22499,44 € se da en el escenario 241, que representa las condiciones opuestas con el menor precio medio y la mayor generación renovable. El coste medio ponderado con la probabilidad de ocurrencia de cada escenario es de 44626,49 €. Hay que considerar, tal y como se explica en [ 52 ], que la aproximación realizada puede implicar importantes desviaciones en cuanto a la optimalidad de la solución. La potencia que se estaba considerando para los generadores con un valor de ug,t que haya sido aproximado a 0 deberá ser producida por otros medios. En este caso, al estar la red limitada, tendría que aumentarse la producción de los generadores encendidos, pero si estos ya hubieran alcanzado su límite, sería necesario encender otros, por lo que cambiaría la solución a la que se ha llegado. Por lo tanto, se plantea necesario utilizar un solver capaz de resolver problemas MINLP para poder conseguir una solución más exacta en el problema estocástico. A pesar de ello, la solución se asemeja en gran medida a la presentada en el caso determinista equivalente: el coste difiere en menos de 6 € y el estado de los generadores es muy similar (ver figura 4.4a). Por lo tanto, este modelo podría ser utilizado para anticipar los generadores que deberán ser puestos en marcha a pesar de que los costes finales difieran en mayor o menor medida.
52 Capítulo 4. Resultados del Modelo Off On Estado de los generadores CHPFC 9 Off On FC 10 Off On FC 5 0 5 10 15 20 25 Horas (h) Off On CHPD Figura 4.11 Estado de los generadores en el caso estocástico.
5 Conclusiones y Líneas Futuras Este proyecto ha tenido como resultado final la modelización de una red de media tensión de referencia en Europa y la creación de un programa de optimización matemática en Pyomo. Se entregan junto a esta memoria los archivos necesarios para la utilización del modelo determinista (archivo .py con el modelo, archivo Excel para introducción de datos y un archivo .py utilizado para crear las gráficas a partir de los resultados del modelo) y los archivos del modelo estocástico (archivo .py con el modelo de referencia, archivos .dat que definen el árbol de escenarios y los archivos Excel con los valores de los parámetros estocásticos en cada escenario). Todos los archivos cuentan con comentarios para facilitar su comprensión y se incluye también un archivo de texto con instrucciones para el uso y modificación de los modelos. El modelo creado sirve para determinar los valores óptimos de generación e importación de energía en la red para un periodo de 24 horas de forma que se mantengan las condiciones operacionales deseadas y se reduzcan los costes. Se han presentado distintos casos para estudiar las funcionalidades del modelo con resultados positivos. Es posible establecer restricciones de reserva de potencia para el control de frecuencia primario y se ha estudiado el comportamiento de los generadores y de desconexión de la industria mediante casos en los que la red tiene una capacidad de importación limitada y diversos costes. Se ha analizado igualmente un caso con alta generación de energías renovables en el que la red pasa a exportar energía por algunas horas. En todos estos casos, el modelo ha encontrado una solución óptima del problema y los valores hallados coinciden con el comportamiento esperado. En el caso estocástico, pese a que se ha tenido que aproximar la solución del problema relajado, los resultados son muy similares a los del caso determinista equivalente, por lo que podría considerarse este modelo para el problema de Asignación de Unidades. A continuación se especifican algunos factores o funcionalidades adicionales que podrían incluirse en el modelo como pasos a seguir del proyecto: • Conectar las subredes a través de un enlace entre los nodos 8 y 14 para estudiar configuraciones de la red o el comportamiento de diferentes enlaces. • Incluir la presencia de vehículos eléctricos que representen el comportamiento de cargas adicionales pero también tengan funcionalidades para actuar como baterías (vehicle-to-grid. • Aumentar la versatilidad de desconexión de las cargas para estudiar en mayor medida la gestión de la demanda (demand-side management) y sus efectos en la red. • Incluir la desconexión de las cargas reservas de potencia para el control primario de frecuencia. • Considerar la alimentación de los sistemas basados en pilas de combustible para modelizar su comportamiento de una manera más fidedigna. •Incluir costes de puesta en marca y aumento o bajada de potencia en los generadores. 53
60 Appendix A. Modelo Determinista - sum(m.q_house[i,ni,t] for (j,ni) in m.NI if j==i) - sum(m.q_ind[i,ind,t] for (j,ind) in m.IND if j==i) == m.V_i[i,t]*m.V_i[i,t]*(m.G_ij[i,i]*sin(0) - m.B_ij[i,i]* cos(0)) + sum(m.V_i[i,t]*m.V_i[j,t] *(m.G_ij[i,j]*sin(m.theta_i[i,t]-m.theta_i[j,t]) - m.B_ij[i,j]*cos(m.theta_i[i,t]-m.theta_i[j,t])) for j in m. I if (i,j) in m.ARCS)) m.reactflux_const = Constraint(m.I, m.T, rule=reactflux_const_rule) ## Calculo de potencia aparente en las lineas (aumenta el tiempo de resolucion del problema) #def Sij_const_rule(m, i, j, t): # A = (sqrt(3)/m.V)**2 # # return ((m.s_ij[i,j,t]**2)*A == (m.G_ij[i,j]**2 + m.B_ij[i,j]**2) # *((m.V_i[i,t]*cos(m.theta_i[i,t]) - m.V_i[j,t]*cos(m.theta_i [j,t]))**2 # +(m.V_i[i,t]*sin(m.theta_i[i,t]) - m.V_i[j,t]*sin(m.theta_i[ j,t]))**2)) # #m.Sij_const = Constraint(m.ARCS, m.T, rule=Sij_const_rule) # LIMITES EN LA RED # Limites de voltaje en nodos def Vbounds_const_rule(m, i, t): return m.Vmin <= m.V_i[i,t] <= m.Vmax m.Vbounds_const = Constraint(m.I, m.T, rule=Vbounds_const_rule) ## Limite de aparente en arcos #def Sijbounds_const_rule(m, i, j, t): # return m.s_ij[i,j,t] <= m.Smax_ij[t] # #m.Sijbounds_const = Constraint(m.ARCS, m.T, rule=Sijbounds_const_rule) # Mantener voltaje en nodos PCC def Vpcc1_const_rule(m, t): return m.Vmin_pcc <= m.V_i[1,t] <= m.Vmax_pcc m.Vpcc1_const = Constraint(m.T, rule=Vpcc1_const_rule) def Vpcc12_const_rule(m, t): return m.Vmin_pcc <= m.V_i[12,t] <= m.Vmax_pcc m.Vpcc12_const = Constraint(m.T, rule=Vpcc12_const_rule)
61 # Fijar el angulo en nodos PCC a cero def Theta1_const_rule(m, t): return m.theta_i[1,t] == 0 m.Theta1_const = Constraint(m.T, rule=Theta1_const_rule) def Theta12_const_rule(m, t): return m.theta_i[12,t] == 0 m.Theta12_const = Constraint(m.T, rule=Theta12_const_rule) # DESCONEXION DE LA INDUSTRIA # Desconexion de potencia activa def actls_const_rule(m, i, ind, t): return m.p_ind[i,ind,t] == m.ls[i,ind,t]*m.pe_ind[i,ind,t] m.actls_const = Constraint(m.IND, m.T, rule=actls_const_rule) # Calculo de potencia reactiva def reactls_const_rule(m, i, ind, t): return m.q_ind[i,ind,t] == m.p_ind[i,ind,t]*tan(m.phi_ind[i,ind]) m.reactls_const = Constraint(m.IND, m.T, rule=reactls_const_rule) # GENERADORES # Limite min de activa en generadores def Pgenlowbound_const_rule(m, i, g, t): return m.u_g[i,g,t]*m.Pmin_g[i,g] <= m.p_g[i,g,t] m.Pgenlowbound_const = Constraint(m.G, m.T, rule=Pgenlowbound_const_rule ) #Limite max de activa en generadores def Pgenupbound_const_rule(m, i, g, t): return m.p_g[i,g,t] <= m.u_g[i,g,t]*m.Pmax_g[i,g] m.Pgenupbound_const = Constraint(m.G, m.T, rule=Pgenupbound_const_rule) # Limite min de reactiva en generadores def Qgenlowbound_const_rule(m, i, g, t): return m.u_g[i,g,t]*m.Qmin_g[i,g] <= m.q_g[i,g,t]
62 Appendix A. Modelo Determinista m.Qgenlowbound_const = Constraint(m.G, m.T, rule=Qgenlowbound_const_rule ) # Limite max de reactiva en generadores def Qgenupbound_const_rule(m, i, g, t): if (g != "FC" and g != "CHPFC"): return m.q_g[i,g,t] <= m.u_g[i,g,t]*m.Qmax_g[i,g] else: return Constraint.Skip m.Qgenupbound_const = Constraint(m.G, m.T, rule=Qgenupbound_const_rule) # Restricción de reactiva segun potencia aparente para generadores con inversor def FCreact_const_rule(m, i, g, t): if (g == "FC" or g == "CHPFC"): return m.q_g[i,g,t]**2 <= (m.Smax_g[i,g]**2 - m.p_g[i,g,t]**2)*m. u_g[i,g,t] else: return Constraint.Skip m.FCreact_const = Constraint(m.G, m.T, rule=FCreact_const_rule) # RENOVABLES # Reactiva en PV def PVreact_const_rule(m, i, pv, t): return m.q_pv[i,pv,t]**2 <= m.Smax_pv[i,pv]**2 - m.p_pv[i,pv,t]**2 m.PVreact_const = Constraint(m.PV, m.T, rule=PVreact_const_rule) # Reactiva en eolica def Wreact_const_rule(m, i, w, t): return m.q_wind[i,w,t]**2 <= m.Smax_w[i,w]**2 - m.p_wind[i,w,t]**2 m.Wreact_const = Constraint(m.WIND, m.T, rule=Wreact_const_rule) # BATERIAS # Balance en las baterias def SOC_const_rule(m, i, bs, t): if (t == 1): soc_0 = m.soc_init[i,bs] else: soc_0 = m.soc[i,bs,t-1]
63 return (m.soc[i,bs,t] == soc_0 + m.pc_bs[i,bs,t]*m.rend_c[i,bs] - m.pd_bs[i,bs,t]/m.rend_d[i,bs]) m.SOC_const = Constraint(m.BS, m.T, rule=SOC_const_rule) # Limites de aparente en baterias def Sbsbounds_const_rule(m, i, bs, t): return ((m.pd_bs[i,bs,t] - m.pc_bs[i,bs,t])**2 + m.q_bs[i,bs,t]**2 <= m.Smax_bs[i,bs]**2) m.Sbsbounds_const = Constraint(m.BS, m.T, rule=Sbsbounds_const_rule) # Limites de las baterias def SOCbounds_const_rule(m, i, bs, t): return (m.soc_min[i,bs] <= m.soc[i,bs,t] <= m.soc_max[i,bs]) m.SOCbounds_const_rule = Constraint(m.BS, m.T, rule=SOCbounds_const_rule ) def PmaxDbs_const_rule(m, i, bs, t): return (m.pd_bs[i,bs,t] <= m.Pmax_bs[i,bs]) m.PmaxDbs_const = Constraint(m.BS, m.T, rule=PmaxDbs_const_rule) def PmaxCbs_const_rule(m, i, bs, t): return (m.pc_bs[i,bs,t] <= m.Pmax_bs[i,bs]) m.PmaxCbs_const = Constraint(m.BS, m.T, rule=PmaxCbs_const_rule) #====================================================================== # RESTRICCIONES AUXILIARES #====================================================================== # Porcentaje máximo de reactiva def Qin_const_rule(m, t): return ((m.q_g[1,’Red’,t] + m.q_g[12,’Red’,t]) <= m.f_Q[t]*(m.p_g[1,’Red’,t] + m.p_g[12,’Red’,t])) m.Qin_const = Constraint(m.T, rule=Qin_const_rule) # Reservas control de frecuencia (inicialmente calculadas sin tener en cuenta los niveles de carga de las baterías) def UpReserve_const_rule(m, t):
64 Appendix A. Modelo Determinista suma_bs = 0 for (i,bs) in m.BS: Ur_bs = (m.Pmax_bs[i,bs] - m.pd_bs[i,bs,t] + m.pc_bs[i,bs,t]) suma_bs += Ur_bs return (sum(m.u_g[i,g,t]*(m.Pmax_g[i,g] - m.p_g[i,g,t]) for (i,g) in m.G if g!=’Red’) + suma_bs >= m.ur[t]) m.UpReserve_const = Constraint(m.T, rule=UpReserve_const_rule) def DownReserve_const_rule(m, t): suma_bs = 0 for (i,bs) in m.BS: Dr_bs = (m.Pmax_bs[i,bs] + m.pd_bs[i,bs,t] - m.pc_bs[i,bs,t]) suma_bs += Dr_bs return (sum(m.u_g[i,g,t]*(m.p_g[i,g,t] - m.Pmin_g[i,g]) for (i,g) in m.G if g!=’Red’) + suma_bs >= m.dr[t]) m.DownReserve_const = Constraint(m.T, rule=DownReserve_const_rule) #====================================================================== # INTRODUCCION DE LOS DATOS #====================================================================== data = DataPortal() # Carga los parametros de entrada data.load(filename=’Datos.xlsx’, range=’horas’, param=’h’) data.load(filename=’Datos.xlsx’, range=’nodos’, param=’n’) # Carga los indices de la red data.load(filename=’Datos.xlsx’, range=’arcos’, format=’set_array’, set= m.ARCS) data.load(filename=’Datos.xlsx’, range=’G_i’, set=m.G) data.load(filename=’Datos.xlsx’, range=’PV_i’, set=m.PV) data.load(filename=’Datos.xlsx’, range=’W_i’, set=m.WIND) data.load(filename=’Datos.xlsx’, range=’NI_i’, set=m.NI) data.load(filename=’Datos.xlsx’, range=’IND_i’, set=m.IND) data.load(filename=’Datos.xlsx’, range=’BS_i’, set=m.BS) # Carga los parametros estocasticos data.load(filename=’Datos.xlsx’, range=’estoc’, select=(’Horas’,’ Precio_red’), param=m.price)
65 data.load(filename=’Datos.xlsx’, range=’estoc’, select=(’Horas’,’f_pv’), param=m.f_pv) data.load(filename=’Datos.xlsx’, range=’estoc’, select=(’Horas’,’f_w’), param=m.f_w) # Carga los parametros de los generadores data.load(filename=’Datos.xlsx’, range=’generadores’, select=(’Nodo’,’ Generador’,’CosteA’), param=m.CosteA_g) data.load(filename=’Datos.xlsx’, range=’generadores’, select=(’Nodo’,’ Generador’,’CosteB’), param=m.CosteB_g) data.load(filename=’Datos.xlsx’, range=’generadores’, select=(’Nodo’,’ Generador’,’Pmin’), param=m.Pmin_g) data.load(filename=’Datos.xlsx’, range=’generadores’, select=(’Nodo’,’ Generador’,’Pmax’), param=m.Pmax_g) data.load(filename=’Datos.xlsx’, range=’generadores’, select=(’Nodo’,’ Generador’,’Qmin’), param=m.Qmin_g) data.load(filename=’Datos.xlsx’, range=’generadores’, select=(’Nodo’,’ Generador’,’Qmax’), param=m.Qmax_g) data.load(filename=’Datos.xlsx’, range=’generadores’, select=(’Nodo’,’ Generador’,’Smax’), param=m.Smax_g) # Carga los parametros de las renovables data.load(filename=’Datos.xlsx’, range=’PV’, select=(’Nodo’,’PV’,’Pmax’) , param=m.Pmax_pv) data.load(filename=’Datos.xlsx’, range=’PV’, select=(’Nodo’,’PV’,’Smax’) , param=m.Smax_pv) data.load(filename=’Datos.xlsx’, range=’wind’,select=(’Nodo’,’Aerogen’,’ Pmax’), param=m.Pmax_w) data.load(filename=’Datos.xlsx’, range=’wind’,select=(’Nodo’,’Aerogen’,’ Smax’), param=m.Smax_w) # Carga los parametros de las cargas esenciales data.load(filename=’Datos.xlsx’, range=’estoc’, select=(’Horas’,’f_house ’), param=m.f_house) data.load(filename=’Datos.xlsx’, range=’cargas_ni’, select=(’Nodo’,’ Carga’,’Pmax’), param=m.Pmax_ni) data.load(filename=’Datos.xlsx’, range=’cargas_ni’, select=(’Nodo’,’ Carga’,’Phi’), param=m.phi_ni) # Carga los parametros de las cargas no esenciale data.load(filename=’Datos.xlsx’, range=’estoc’, select=(’Horas’,’f_ind’) , param=m.f_ind) data.load(filename=’Datos.xlsx’, range=’cargas_ind’, select=(’Nodo’,’ Carga’,’Pmax’), param=m.Pmax_ind) data.load(filename=’Datos.xlsx’, range=’cargas_ind’, select=(’Nodo’,’ Carga’,’Phi’), param=m.phi_ind) data.load(filename=’Datos.xlsx’, range=’cargas_ind’, select=(’Nodo’,’ Carga’,’Coste’), param=m.C_ind) # Carga los parametros de las baterias
66 Appendix A. Modelo Determinista data.load(filename=’Datos.xlsx’, range=’baterias’, select=(’Nodo’,’ Bateria’,’SOCinit’), param=m.soc_init) data.load(filename=’Datos.xlsx’, range=’baterias’, select=(’Nodo’,’ Bateria’,’SOCmin’), param=m.soc_min) data.load(filename=’Datos.xlsx’, range=’baterias’, select=(’Nodo’,’ Bateria’,’SOCmax’), param=m.soc_max) data.load(filename=’Datos.xlsx’, range=’baterias’, select=(’Nodo’,’ Bateria’,’Rcarga’), param=m.rend_c) data.load(filename=’Datos.xlsx’, range=’baterias’, select=(’Nodo’,’ Bateria’,’Rdescarga’), param=m.rend_d) data.load(filename=’Datos.xlsx’, range=’baterias’, select=(’Nodo’,’ Bateria’,’Smax’), param=m.Smax_bs) data.load(filename=’Datos.xlsx’, range=’baterias’, select=(’Nodo’,’ Bateria’,’Pmax’), param=m.Pmax_bs) # Carga los parametros de la red data.load(filename=’Datos.xlsx’, range=’Gij’, format=’array’, param=m. G_ij) data.load(filename=’Datos.xlsx’, range=’Bij’, format=’array’, param=m. B_ij) data.load(filename=’Datos.xlsx’, range=’Auxiliar’, select=(’Horas’,’Smax ’), param=m.Smax_ij) data.load(filename=’Datos.xlsx’, range=’Vbase’, param=’V’) # Carga los parametros de los constraints auxiliares data.load(filename=’Datos.xlsx’, range=’Auxiliar’, select=(’Horas’,’f_Q ’), param=m.f_Q) data.load(filename=’Datos.xlsx’, range=’Auxiliar’, select=(’Horas’,’Ur’) , param=m.ur) data.load(filename=’Datos.xlsx’, range=’Auxiliar’, select=(’Horas’,’Dr’) , param=m.dr) # Se crea el modelo concreto con los datos introducidos instance = m.create_instance(data) ## Linea de codigo opcional para crear un archivo GAMS con el modelo concreto #instance.write(’determ.gms’) #====================================================================== # RESOLUCION #====================================================================== # Para resolver remotamente a través de la interfaz NEOS solver_manager = SolverManagerFactory(’neos’) results = solver_manager.solve(instance, opt=’bonmin’)
67 results.write() #============================================================================== # Iteraciones para considerar el nivel de carga de las baterias en las reservas #============================================================================== # Comprobamos si hay restricciones de reserva y si es el caso se comprueba si hay casos conflictivos iterar = False comprob = False for t in instance.T: if value(instance.ur[t]) > 0 or value(instance.dr[t]) > 0: comprob = True break if comprob == True: for t in instance.T: for (i,bs) in instance.BS: Ur_bs = (instance.Pmax_bs[i,bs] - instance.pd_bs[i,bs,t] + instance.pc_bs[i,bs,t]) if value(Ur_bs)/6 >= value(instance.soc[i,bs,t] - instance. soc_min[i,bs]): iterar = True break Dr_bs = (instance.Pmax_bs[i,bs] + instance.pd_bs[i,bs,t] - instance.pc_bs[i,bs,t]) if value(Dr_bs)/6 >= value(instance.soc_max[i,bs] - instance. soc[i,bs,t]): iterar = True break if iterar == True: break # Se itera añadiendo restricciones en las baterias y horas que no tengan el nivel # de carga adecuado en la solucion anterior while (iterar == True): # Listas para la comprobacion posterior, guardan la batería y la hora para las que se añade restriccion instance.modif_up = [] instance.modif_down = []
68 Appendix A. Modelo Determinista def UpReserveAdd_const_rule(instance, t): suma_bs = 0 changed = False for (i,bs) in instance.BS: Ur_bs = (instance.Pmax_bs[i,bs] - instance.pd_bs[i,bs,t] + instance.pc_bs[i,bs,t]) if value(Ur_bs)/6 >= value(instance.soc[i,bs,t] - instance. soc_min[i,bs]): Ur_bs = instance.soc[i,bs,t] - instance.soc_min[i,bs] instance.modif_up.append(str(i)+’ ’+str(t)) changed = True suma_bs += Ur_bs # Si no tiene nivel de carga suficiente, se añade la restricción para tenerlo en cuenta if changed == True: return (sum(instance.u_g[i,g,t]*(instance.Pmax_g[i,g] - instance.p_g[i,g,t]) for (i,g) in instance.G if g!=’Red’) + suma_bs >= instance.ur[t]) else: return Constraint.Skip instance.UpReserveAdd_const = Constraint(instance.T, rule= UpReserveAdd_const_rule) def DownReserveAdd_const_rule(instance, t): suma_bs = 0 changed = False for (i,bs) in instance.BS: Dr_bs = (instance.Pmax_bs[i,bs] + instance.pd_bs[i,bs,t] - instance.pc_bs[i,bs,t]) if value(Dr_bs)/6 >= value(instance.soc_max[i,bs] - instance. soc[i,bs,t]): Dr_bs = instance.soc_max[i,bs] - instance.soc[i,bs,t] instance.modif_down.append(str(i)+’ ’+str(t)) changed = True suma_bs += Dr_bs # Si no tiene margen de carga suficiente, se añade la restricció n para tenerlo en cuenta if changed == True: return (sum(instance.u_g[i,g,t]*(instance.p_g[i,g,t] - instance.Pmin_g[i,g]) for (i,g) in instance.G if g!=’Red’) + suma_bs >= instance.dr[t]) else: return Constraint.Skip
69 instance.DownReserveAdd_const = Constraint(instance.T, rule= DownReserveAdd_const_rule) # Vuelve a resolverse el problema de optimizacion solver_manager = SolverManagerFactory(’neos’) results = solver_manager.solve(instance, opt=’bonmin’) results.write() # Comprobar si en la nueva solución hay algun caso conflictivo que no haya sido ya resuelto iterar = False comp = ’’ for t in instance.T: for (i,bs) in instance.BS: Ur_bs = (instance.Pmax_bs[i,bs] - instance.pd_bs[i,bs,t] + instance.pc_bs[i,bs,t]) if value(Ur_bs)/6 >= value(instance.soc[i,bs,t] - instance. soc_min[i,bs]): comp = str(i) + ’ ’ + str(t) if comp not in instance.modif_up: iterar = True break Dr_bs = (instance.Pmax_bs[i,bs] + instance.pd_bs[i,bs,t] - instance.pc_bs[i,bs,t]) if value(Dr_bs)/6 >= value(instance.soc_max[i,bs] - instance. soc[i,bs,t]): comp = str(i) + ’ ’ + str(t) if comp not in instance.modif_down: iterar = True break if iterar == True: break #====================================================================== # GRAFICA 1: Voltajes y Potencias #====================================================================== # Vector con las horas y los colores horas = np.arange(1,int(value(instance.h))+1) colores = [’b’,’g’,’r’,’c’,’m’,’y’,’k’,’#FFC0CB’,’#00FF00’,’#FFA500 ’,’#008080’,’#FF00FF’,’#4B0082’,’#D2691E’,’b’,’g’,’r’,’c’,’m’,’y’,’k ’,’#FFC0CB’,’#00FF00’,’#FFA500’,’#0080000’,’#FF00FF’,’#4B0082’,’# D2691E’] fig1 = plt.figure()
76 Appendix A. Modelo Determinista # Se guarda la altura máxima para los ejes if (h[t-1] + ls[i][t-1]) > Dmax: Dmax = h[t-1] + ls[i][t-1] f3_ax.bar(horas, ls[i], bottom=h, width=bar_width, align=’center ’, alpha=opacity, color=colores[i], edgecolor=None, label=’ Ind %i’%descon[i]) f3_ax.yaxis.set_major_locator(plt.MaxNLocator(5)) f3_ax.set_ylim([0, Dmax + 5]) f3_ax.set_yticklabels([]) f3_ax.set(xlabel=’Horas (h)’, ylabel=’Desconexion (%)’, title=’ Desconexion cargas industriales’) f3_ax.grid() f3_lgd = f3_ax.legend(bbox_to_anchor=(1.01, 1), loc=’upper left’, borderaxespad=0.) fig3.savefig("industria.pdf", bbox_inches=’tight’) plt.show() #====================================================================== # GRAFICA 4: Perfiles de las cargas #====================================================================== fig4 = plt.figure() params = {’figure.figsize’: (5, 5)} plt.rcParams.update(params) fig4, (f4_ax1, f4_ax2) = plt.subplots(2, sharex=True) # Perfil residencial res = [] for t in horas: res.append(value(instance.f_house[t])) f4_ax1.plot(horas, res) f4_ax1.set_ylim([0,1]) f4_ax1.ticklabel_format(useOffset=False) f4_ax1.set(ylabel=’Potencia activa [p.u.]’, title=’Perfil de las cargas residenciales’) f4_ax1.grid()
77 # Perfil industrial ind = [] for t in horas: ind.append(value(instance.f_ind[t])) f4_ax2.plot(horas, ind) f4_ax2.set_ylim([0,1]) f4_ax2.set(xlabel=’Horas (h)’, ylabel=’Potencia activa [p.u.]’, title=’ Perfil de las cargas industriales’) f4_ax2.grid() fig4.subplots_adjust(hspace=0.3) plt.setp([a.get_xticklabels() for a in fig1.axes[:-1]], visible=False) fig4.savefig("cargas.pdf", bbox_inches=’tight’) plt.show() #====================================================================== # GRAFICA 5: Precios de importación #====================================================================== fig5 = plt.figure() params = {’figure.figsize’: (5, 3)} plt.rcParams.update(params) fig5, f5_ax = plt.subplots() p = [] Pmax = 0 Pmin = 1000 #Crea el vector con los valores for t in horas: p.append(value(instance.price[t])) if p[-1] > Pmax: Pmax = p[-1] if p[-1] < Pmin: Pmin = p[-1] f5_ax.plot(horas, p) f5_ax.set_ylim([Pmin - 5, Pmax + 5])
78 Appendix A. Modelo Determinista f5_ax.set(xlabel = ’Horas (h)’, ylabel=’Coste (%s/MWh)’%(u"\N{euro sign }"), title=’Coste de importacion de energia’) f5_ax.grid() fig5.savefig("precios.pdf", bbox_inches=’tight’) plt.show() #====================================================================== # GRAFICA 6: Renovables #====================================================================== fig6 = plt.figure() params = {’figure.figsize’: (5, 5)} plt.rcParams.update(params) fig6, (f6_ax1, f6_ax2) = plt.subplots(2, sharex=True) # Perfil PV pv = [] #Crea el vector con los valores for t in horas: pv.append(value(instance.f_pv[t])) f6_ax1.plot(horas, pv) f6_ax1.yaxis.set_major_locator(plt.MaxNLocator(5)) f6_ax1.set_ylim([0,0.4]) f6_ax1.ticklabel_format(useOffset=False) f6_ax1.set(ylabel=’Potencia activa [p.u.]’, title=’Perfil de generacion fotovoltaica’) f6_ax1.grid() w = [] for t in horas: w.append(value(instance.f_w[t])) f6_ax2.plot(horas, w) f6_ax2.yaxis.set_major_locator(plt.MaxNLocator(5)) f6_ax2.set_ylim([0.2,0.6]) f6_ax2.set(xlabel=’Horas (h)’, ylabel=’Potencia activa [p.u.]’, title=’ Perfil de generacion eolica’) f6_ax2.grid()
79 fig6.subplots_adjust(hspace=0.3) plt.setp([a.get_xticklabels() for a in fig6.axes[:-1]], visible=False) fig6.savefig("renovables.pdf", bbox_inches=’tight’) plt.show() #====================================================================== # GRAFICA 7: Potencia generadores renovables #====================================================================== fig7 = plt.figure() params = {’legend.fontsize’: 10, ’legend.handlelength’: 3, ’figure.figsize’: (5,5)} plt.rcParams.update(params) fig7, (f7_ax1, f7_ax2) = plt.subplots(2, sharex=True) aux_pv = [] aux_wind = [] Pmax_pv = 0 Pmax_wind = 0 bar_width = 1 for t in horas: suma_pv = 0 suma_wind = 0 for (i,pv) in instance.PV: suma_pv += value(instance.p_pv[i,pv,t]) for (i,w) in instance.WIND: suma_wind += value(instance.p_wind[i,w,t]) aux_pv.append(suma_pv) aux_wind.append(suma_wind) if suma_pv > Pmax_pv: Pmax_pv = suma_pv if suma_wind > Pmax_wind: Pmax_wind = suma_wind f7_ax1.bar(horas, aux_pv, width=bar_width, align=’center’, alpha=opacity , edgecolor=None) f7_ax1.set_ylim([0, Pmax_pv+0.1])
80 Appendix A. Modelo Determinista f7_ax1.set(ylabel=’Potencia (MW)’, title=’Potencia de los generadores fotovoltaicos’) f7_ax1.grid() f7_ax1.ticklabel_format(useOffset=False) f7_ax1.yaxis.set_major_locator(plt.MaxNLocator(5)) f7_ax2.bar(horas, aux_wind, width=bar_width, align=’center’, alpha= opacity, edgecolor=None) f7_ax2.set_ylim([0, Pmax_wind+0.2]) f7_ax2.set(xlabel=’Horas (h)’, ylabel=’Potencia (MW)’, title=’Potencia de los aerogeneradores’) f7_ax2.grid() f7_ax2.yaxis.set_major_locator(plt.MaxNLocator(5)) fig7.subplots_adjust(hspace=0.3) plt.setp([a.get_xticklabels() for a in fig7.axes[:-1]], visible=False) fig7.savefig("renovables_potencias.pdf", bbox_inches=’tight’) plt.show() #====================================================================== # MOSTRAR RESULTADOS POR PANTALLA #====================================================================== results.write() instance.OBJ.display() pg = 0 for (i,g) in instance.G: if g == ’Red’: for t in horas: pg += (value(instance.p_g[i,g,t])) print(’La energía total importada es %f MWh’%pg) #====================================================================== # GUARDAR LOS RESULTADOS #====================================================================== # Guarda la solucion en results instance.solutions.store_to(results) # Archiva los resultados: el archivo es usado por Solucion.py para crear las graficas # sin necesidad de volver a resolver el problema de optimizacion file = open("results.pickle", "wb")
81 pickle.dump(results, file) file.close() # Escribe un archivo CSV con los resultados with open(’results.csv’, mode=’w’) as Result: results_writer = csv.writer(Result, delimiter=’;’, quotechar=’"’, quoting=csv.QUOTE_NONNUMERIC) results_writer.writerow([’OBJ’,value(instance.OBJ)]) # Variables for v in instance.component_objects(Var, active=True): if v.name == ’p_g’ or v.name == ’u_g’: for (nodo,x) in instance.G: for t in instance.T: results_writer.writerow([v[nodo,x,t],value(v[nodo,x,t ])]) elif v.name == ’pd_bs’ or v.name == ’pc_bs’ or v.name == ’soc’: for (nodo,x) in instance.BS: for t in instance.T: results_writer.writerow([v[nodo,x,t],value(v[nodo,x,t ])]) elif v.name == ’p_ind’ or v.name == ’ls’: for (nodo,x) in instance.IND: for t in instance.T: results_writer.writerow([v[nodo,x,t],value(v[nodo,x,t ])]) elif v.name == ’s_ij’: for (i,j) in instance.ARCS: for t in instance.T: results_writer.writerow([v[i,j,t],value(v[i,j,t])]) elif v.name == ’V_i’: for nodo in instance.I: for t in instance.T: results_writer.writerow([v[nodo,t],value(v[nodo,t])]) elif v.name == ’theta_i’: for nodo in instance.I: for t in instance.T: results_writer.writerow([v[nodo,t],value(v[nodo,t])]) # Parametros for v in instance.component_objects(Param, active=True): if v.name == ’p_house’: for (nodo,x) in instance.NI:
82 Apéndice A. Modelo Determinista for t in instance.T: results_writer.writerow([v.name,nodo,x,t,v[nodo,x,t]]) elif v.name == ’p_pv’: for (nodo,x) in instance.PV: for t in instance.T: results_writer.writerow([v.name,nodo,x,t,v[nodo,x,t]]) elif v.name == ’p_wind’: for (nodo,x) in instance.WIND: for t in instance.T: results_writer.writerow([v.name,nodo,x,t,v[nodo,x,t]]) elif v.name == ’price’: for t in instance.T: results_writer.writerow([v.name,t,v[t]])
Apéndice B Estructura del problema estocástico La estructura del modelo estocástico se especifica mediante un archivo de text denominado ScenarioStructurre.dat: Código B.1 Árbol de escenarios del problema estocástico. set Stages := FirstStage SecondStage ; set Nodes := RootNode # Precios1 (4 nodos PV y 4 nodos Wind) Node111 Node112 Node113 Node114 Node121 Node122 Node123 Node124 Node131 Node132 Node133 Node134 Node141 Node142 Node143 Node144 # Precios 2 Node211 Node212 Node213 Node214 83
84 Appendix B. Estructura del problema estocástico Node221 Node222 Node223 Node224 Node231 Node232 Node233 Node234 Node241 Node242 Node243 Node244 # Precios 3 Node311 Node312 Node313 Node314 Node321 Node322 Node323 Node324 Node331 Node332 Node333 Node334 Node341 Node342 Node343 Node344 # Precios 4 Node411 Node412 Node413 Node414 Node421 Node422 Node423 Node424 Node431 Node432 Node433
85 Node434 Node441 Node442 Node443 Node444 ; param NodeStage := RootNode FirstStage Node111 SecondStage Node112 SecondStage Node113 SecondStage Node114 SecondStage Node121 SecondStage Node122 SecondStage Node123 SecondStage Node124 SecondStage Node131 SecondStage Node132 SecondStage Node133 SecondStage Node134 SecondStage Node141 SecondStage Node142 SecondStage Node143 SecondStage Node144 SecondStage # Precios 2 Node211 SecondStage Node212 SecondStage Node213 SecondStage Node214 SecondStage Node221 SecondStage Node222 SecondStage Node223 SecondStage Node224 SecondStage Node231 SecondStage Node232 SecondStage Node233 SecondStage Node234 SecondStage Node241 SecondStage Node242 SecondStage Node243 SecondStage Node244 SecondStage # Precios 3
92 Appendix B. Estructura del problema estocástico Scenario131 Node131 Scenario132 Node132 Scenario133 Node133 Scenario134 Node134 Scenario141 Node141 Scenario142 Node142 Scenario143 Node143 Scenario144 Node144 # Precios 2 Scenario211 Node211 Scenario212 Node212 Scenario213 Node213 Scenario214 Node214 Scenario221 Node221 Scenario222 Node222 Scenario223 Node223 Scenario224 Node224 Scenario231 Node231 Scenario232 Node232 Scenario233 Node233 Scenario234 Node234 Scenario241 Node241 Scenario242 Node242 Scenario243 Node243 Scenario244 Node244 # Precios 3 Scenario311 Node311 Scenario312 Node312 Scenario313 Node313 Scenario314 Node314 Scenario321 Node321 Scenario322 Node322 Scenario323 Node323 Scenario324 Node324 Scenario331 Node331 Scenario332 Node332 Scenario333 Node333 Scenario334 Node334 Scenario341 Node341 Scenario342 Node342 Scenario343 Node343
93 Scenario344 Node344 # Precios 4 Scenario411 Node411 Scenario412 Node412 Scenario413 Node413 Scenario414 Node414 Scenario421 Node421 Scenario422 Node422 Scenario423 Node423 Scenario424 Node424 Scenario431 Node431 Scenario432 Node432 Scenario433 Node433 Scenario434 Node434 Scenario441 Node441 Scenario442 Node442 Scenario443 Node443 Scenario444 Node444 ; set StageVariables[FirstStage] := u_g[*,*,*] ; set StageVariables[SecondStage] := p_g[*,*,*] ls[*,*,*] soc[*,*,*] V_i[*,*] theta_i[*,*] ; param StageCost := FirstStage OBJ SecondStage OBJ ; param ScenarioBasedData := False ;
Índice de Figuras 1.1 Predicciones de crecimiento de la población mundial [13] 1 1.2 Predicciones de crecimiento de la demanda de energía primaria [21] 2 1.3 Riesgos e impactos del cambio climático [36] 3 1.4 Evolución necesaria del vehículo eléctrico para alcanzar el objetivo de 2 ◦C según [17] 4 1.5 Actualización de las previsiones de capacidad instalada de energía solar [21] 5 1.6 Evolución de la potencia renovable instalada a nivel mundial [17] 6 1.7 Esquema general de una red de distribución eléctrica inteligente [15] 7 1.8 Representación de los elementos de una red de distribución eléctrica inteligente (fuente: elaboración propia) 8 1.9 Categorías de DSM [39] 10 1.10 Impacto de diferentes tecnologías en el consumo eléctrico diario [39] 11 2.1 Algoritmos para resolver problemas de optimización MINLP (fuente: elaboración propia) 17 2.2 Resumen de Pyomo 18 3.1 Red de referencia [40] 22 3.2 Perfiles de las cargas de la red [45] 24 3.3 Incremento de potencia disponible para control primario de frecuencia (reserva de subida) 32 3.4 Coste medio de importación de energía a la red 38 3.5 Perfiles medios de generación renovable 38 3.6 Creación de escenarios de precios del mercado eléctrico 40 3.7 Escenarios de costes de la energía importada por la red 40 3.8 Creación de escenarios de generación de energía eólica 40 3.9 Escenarios de generación eólica 41 3.10 Escenarios de generación fotovoltaica 41 4.1 Producción de los generadores renovables 43 4.2 Consumo industrial en el caso base 44 4.3 Voltajes, potencias y energía en el caso base 45 4.4 Voltajes, potencias y energía con importación limitada 46 4.5 Desconexión de la industria con importación limitada 46 4.6 Voltajes, potencias y energía con coste de desconexión reducido 47 4.7 Desconexión de las industrias con coste de desconexión reducido 48 4.8 Producción de los generadores renovables aumentada 48 4.9 Voltajes, potencias y energía con alta generación de renovables 49 4.10 Voltajes, potencias y energía con restricciones de reserva 50 4.11 Estado de los generadores en el caso estocástico 52 95
Índice de Tablas 2.1 Problemas de optimización 15 3.1 Parámetros de las líneas de la red [45] 23 3.2 Parámetros de las cargas de la red [45] 24 3.3 Parámetros de los generadores de la red [45] 25 3.4 Formulación del problema optimización determinista 33 3.5 Parámetros de la red [45] 37 3.6 Costes de la red 37 3.7 Parámetros de las baterías 39 97
Índice de Códigos 3.1 Declaración del modelo, sets y variables en Pyomo 34 3.2 Declaración de parámetros mediante función en Pyomo 34 3.3 Declaración de función objetivo y restricciones en Pyomo 35 3.4 Asignación de valores a los parámetros de Pyomo a partir de Excel 35 3.5 Resolución del problema matemático y creación de gráficas en Pyomo 36 3.6 Ejecución del problema estocástico a través del terminal de comandos de Windows 42 A.1 Modelo determinista en Pyomo 55 B.1 Árbol de escenarios del problema estocástico 83 99
Bibliografía [1] COIN-OR: Computational Infrastructure for Operations Research – Open-source Software for the Operations Research Community. [2] Conference of the Parties (COP) | UNFCCC. [3] El grito de los jóvenes contra el cambio climático se convierte en global | Sociedad | EL PAÍS. [4] GAMS/DICOPT. [5] Global Warming of 1.5 ºC — IPCC Special Report. [6] IEEE PES Intelligent Systems Subcommitte - Open Data Sets. [7] La transición energética. Una Nueva Cultura de la Energía. [8] Lundi Vert. [9] Operador del Mercado Ibérico de Energía (OMIE). [10] Pyomo Documentation — Pyomo 5.1 documentation. [11] Sowee : Contrats énergie, Station connectée, Pilotage chauffage. [12] Sustainable Development Goals - United Nations Sustainable Development. [13] World Population Prospects - Population Division - United Nations. [14] Couenne, an exact solver for nonconvex MINLPs, (2006). [15] Electricity Grid Modernization: Progress Being Made on Cybersecurity Guidelines, but Key Challenges Remain to be Addressed, Tech. Report GAO-11-117, U.S. Government Accountability Office, 1 2011. [16] OMS | Los efectos sobre la salud, 2014. [17] Chapter 2 of Perspectives for the Energy Transition - Investment Needs for a Low-Carbon Energy System, Tech. report, OECD/IEA, 2017. [18] IEA Global Energy and CO2 Status Report, Tech. report, 2017. [19] EC Communication - A Clean Planet for all A European strategic long-term vision for a prosperous, modern, competitive and climate neutral economy, Tech. report, European Comission, 2018. 101