Modelado y control de unaCentral Hidroeléctrica
Sara Trujillo Naharro
Trabajo Fin de Máster
Máster en Automática, Robótica y Telemática
Tutor:
Eduardo Fernández Camacho
Departamento de Ingeniería de Sistemas y Automática
Escuela Técnica Superior de Ingeniería
Universidad de Sevilla
Diciembre de 2012
Índice
1. Introducción 11.1. Motivación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2. Resumen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
2. Energía hidráulica 32.1. Aprovechamiento de la energía hidráulica . . . . . . . . . . . . . . . . . . . . . . . . . . 32.2. Central hidroeléctrica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.2.1. Características de una central hidroeléctrica . . . . . . . . . . . . . . . . . . . . . 42.2.2. Tipos de centrales hidroeléctricas . . . . . . . . . . . . . . . . . . . . . . . . . . 42.2.3. Partes principales de una central hidroeléctrica . . . . . . . . . . . . . . . . . . . 5
3. Modelado de una central hidroeléctrica 63.1. Descripción del sistema . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63.2. Subsistemas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
3.2.1. Modelo del tramo de río . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63.2.2. Modelo del lago . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93.2.3. Modelo del conducto . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93.2.4. Modelo de la turbina . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103.2.5. Modelo de la bomba . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103.2.6. Modelo de conductos equipados con una turbina y una bomba . . . . . . . . . . . 10
4. Implementación del sistema en Simulink 124.1. Modelo de un tramo de río . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124.2. Modelo de una presa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134.3. Modelo de los lagos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144.4. Modelo del conducto U1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144.5. Modelo de una turbina . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154.6. Modelo de un conducto equipado con turbina y bomba . . . . . . . . . . . . . . . . . . . 164.7. Modelo del sistema completo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164.8. Parámetros del sistema . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
5. Control predictivo basado en modelo 245.1. Modelo de predicción en el espacio de estados . . . . . . . . . . . . . . . . . . . . . . . . 255.2. Función objetivo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
6. Diseño del controlador 286.1. Punto de funcionamiento y restricciones . . . . . . . . . . . . . . . . . . . . . . . . . . . 286.2. Linealización y modelo de predicción . . . . . . . . . . . . . . . . . . . . . . . . . . . . 286.3. Escenario de control y función objetivo . . . . . . . . . . . . . . . . . . . . . . . . . . . 306.4. Implementación del controlador en Simulink . . . . . . . . . . . . . . . . . . . . . . . . . 31
7. Simulación y resultados 327.1. Ponderación de la función objetivo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 327.2. Referencia de potencia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 327.3. Perturbaciones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
7.3.1. Caudal de entrada en Tramo R1 . . . . . . . . . . . . . . . . . . . . . . . . . . . 337.3.2. Caudal de entrada en Tramo R3 . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
7.4. Seguimiento de la potencia de referencia . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
8. Conclusiones 42
II
Índice de tablas1. Parámetros relacionados con el Subsistema 1 . . . . . . . . . . . . . . . . . . . . . . . . 192. Parámetros relacionados con el Subsistema 2 . . . . . . . . . . . . . . . . . . . . . . . . 203. Parámetros relacionados con el Tramo R1 del río . . . . . . . . . . . . . . . . . . . . . . . 204. Parámetros relacionados con el Tramo R2 del río . . . . . . . . . . . . . . . . . . . . . . . 215. Parámetros relacionados con el Tramo R3 del río . . . . . . . . . . . . . . . . . . . . . . . 216. Parámetros relacionados con el Tramo R4 del río . . . . . . . . . . . . . . . . . . . . . . . 227. Parámetros relacionados con el Tramo R5 del río . . . . . . . . . . . . . . . . . . . . . . . 228. Parámetros relacionados con el Tramo R6 del río . . . . . . . . . . . . . . . . . . . . . . . 239. Datos Generales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
III
Índice de figuras1. Estructura de generación en tiempo real . . . . . . . . . . . . . . . . . . . . . . . . . . . 12. Ejemplo de central hidroeléctrica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63. Visión de conjunto del sistema . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74. Sección de un tramo del río . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85. Discretización de un tramo de río . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86. Implementación en Simulink de la librería de Tramos . . . . . . . . . . . . . . . . . . . . 127. Implementación en Simulink de la librería de Presas . . . . . . . . . . . . . . . . . . . . . 138. Implementación en Simulink del Lago 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . 149. Implementación en Simulink del Conducto U1 . . . . . . . . . . . . . . . . . . . . . . . . 1510. Implementación en Simulink de la Turbina 1 . . . . . . . . . . . . . . . . . . . . . . . . . 1511. Implementación en Simulink de un conducto C1 equipado con turbina y bomba . . . . . . 1612. Esquema en Simulink de la central hidroeléctrica completa . . . . . . . . . . . . . . . . . 1713. Esquema en Simulink de tramo y presa . . . . . . . . . . . . . . . . . . . . . . . . . . . 1714. Esquema en Simulink del Subsistema 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . 1815. Esquema en Simulink del Subsistema 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . 1816. Estrategia de Control Predictivo basado en Modelo (MPC) . . . . . . . . . . . . . . . . . 2517. Implementación en Simulink del sistema de control para la central hidroeléctrica . . . . . 3118. Precio de la electricidad γ durante 24 horas . . . . . . . . . . . . . . . . . . . . . . . . . 3219. Referencia de potencia pre f a seguir en 24 horas . . . . . . . . . . . . . . . . . . . . . . . 3220. Variación del caudal de entrada en el Tramo R1 durante 24 horas . . . . . . . . . . . . . . 3321. Variación del caudal lateral en el Tramo R3 durante 24 horas . . . . . . . . . . . . . . . . 3322. Potencia total generada por la central durante 24 horas . . . . . . . . . . . . . . . . . . . 3423. Detalle de la potencia total generada por la central . . . . . . . . . . . . . . . . . . . . . . 3424. Nivel del Lago L1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3525. Nivel del Lago L2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3526. Nivel del Lago L3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3527. Nivel al final del tramo R1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3628. Nivel al final del tramo R2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3629. Nivel al final del tramo R3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3630. Nivel al final del tramo R4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3731. Nivel al final del tramo R5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3732. Nivel al final del tramo R6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3733. Caudal Turbina T1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3834. Caudal Turbina/Bomba C1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3835. Caudal Turbina T2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3836. Caudal Turbina/Bomba C2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3937. Caudal Turbina D1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3938. Caudal Turbina D2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3939. Caudal Turbina D3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4040. Caudal Turbina D4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4041. Caudal Turbina D5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4042. Caudal Turbina D6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4143. Nivel de las secciones de todos los tramos después de 24 horas . . . . . . . . . . . . . . . 41
IV
A Luis,
por darme tu ánimo, apoyo y alegría para afrontar
éste y todos los proyectos que empredemos juntos
1. Introducción
1.1. MotivaciónLa energía eléctrica no se puede almacenar en grandes cantidades y, por ello, es necesario que en cada
momento se genere la cantidad precisa que se necesite. Encender un interruptor o conectar un electrodo-méstico implica que en ese mismo instante una central eléctrica debe producir la electricidad necesaria.
El operador del mercado, responsable de este equilibrio, trabaja por un uso sostenible de la energía conel objetivo de contribuir al aplanamiento de la curva de la demanda mediante un cambio en los comporta-mientos a la hora de consumir energía, lo que supone una mayor eficiencia en el conjunto del sistema y unamejor integración de las energías renovables. Para lograr este cambio destacan las medidas de eficienciay ahorro energético, la discriminación horaria, la gestión automática de cargas o el servicio de gestión dedemanda de interrumpibilidad.
La gestión de la demanda es la planificación e implementación de distintas medidas destinadas a influiren el modo de consumir energía para que se modifique el perfil de consumo diario. Estas medidas contribu-yen a la reducción de las emisiones de CO2, a la mejor integración de las energías renovables en el sistemaeléctrico y a una mayor eficiencia energética del sistema en su conjunto.
Las medidas de gestión de la demanda se clasifican en cuatro grandes grupos en función del tipo deimpacto que tienen sobre la curva de la demanda: reducción del consumo, desplazamiento del consumo delas horas punta a las valle, llenado de horas valles y reducción del consumo en las horas punta.
Diariamente, como puede verse en la Figura 1, la demanda de energía eléctrica se cubre con la genera-ción de centrales basadas en distintas tecnologías, como son:
Hidráulica.
Nuclear.
Fuel/Gas.
Carbón.
Ciclo combinado.
Eólica.
Resto de generadores de Régimen Especial.
Figura 1: Estructura de generación en tiempo real
1
La energía hidráulica tiene la cualidad de ser renovable, pues no agota la fuente primaria al explotarla,y es limpia, ya que no produce en su explotación sustancias contaminantes de ningún tipo. Sin embargo, elimpacto medioambiental de las grandes presas, por la severa alteración del paisaje e, incluso, la inducción deun microclima diferenciado en su emplazamiento, ha desmerecido la bondad ecológica de este concepto enlos últimos años. Recientemente se están realizando centrales minihidroeléctricas, mucho más respetuosascon el ambiente y que se benefician de los progresos tecnológicos, logrando un rendimiento y una viabilidadeconómica razonables.
La motivación principal para realizar este trabajo fue la de estudiar un sistema basado en energías reno-vables que permitiera adecuarse a la demanda haciendo uso del concepto de almacenamiento. El objetivofundamental ha sido el de realizar el modelado y control de una central hidroeléctrica de bombeo, por lacapacidad de ésta para contribuir a que la oferta se adecue al perfil de demanda diario, proporcionandola flexibilidad necesaria al sistema para que, como decíamos al principio, en cada momento se genere lacantidad precisa que se necesite.
En el ajuste continuo de la producción a la demanda es necesario disponer de centrales cuya poten-cia pueda ser fácilmente regulable, con una gran flexibilidad de operación. Las centrales hidroeléctricaspresentan estas características jugando un papel muy importante en el conjunto del parque de centrales degeneración de energía eléctrica de cualquier país. Son instalaciones con una alta velocidad de respuesta antelos cambios de demanda, lo que quiere decir que en unos minutos pasan de estar paradas a dar la potencianominal. Esto no ocurre con las centrales de combustible fósil o nuclear, que necesitan horas, dependiendode las condiciones en las que se produzca el arranque de las mismas. Por todo esto, las centrales hidroeléc-tricas se convierten en las instalaciones más adecuadas para cubrir las puntas de demanda, así como paracubrir las bajas imprevistas de otras centrales.
1.2. ResumenEl agua proveniente de la evaporación de los océanos, además de servir para otros fines, tales como
riego, limpieza, enfriamiento, consumo, etc., que lo convierten en un liquido vital para los seres humanos,se utiliza también para accionar máquinas giratorias llamadas turbinas, que a su vez mueven generadoresque transforman la energía mecánica en energía eléctrica.
Las plantas hidroeléctricas aprovechan los caudales y caídas de agua. El sol calienta las masas de agua,de su evaporación se forman nubes y eventualmente lluvia que fluye a través de caudalosos ríos. El aguaen estos ríos tiene una enorme cantidad de energía mecánica potencial, y para aprovechar esta energía seseleccionan cauces de ríos que tienen algunas características importantes que incluyen amplio caudal deagua y diferencias importantes de altura en corta distancia.
El objetivo fundamental de este trabajo ha sido el de realizar el modelado y control de una central hi-droeléctrica de bombeo. En la sección 2, se introduce el concepto de energía hidráulica, se mencionan lostipos de centrales hidroeléctricas que existen atendiendo a distintos criterios y se definen las partes queconforman una central. En la sección 3, se aborda la obtención de un modelo matemático de una centralhidroeléctrica. En primer lugar, se describe el sistema a modelar. A continuación, se desarrollan los mode-los de cada uno de los subsistemas y las partes de que consta la central: lagos, tramos de río, conductos,turbinas, bombas, etc. En la sección 4, se describe la implementación en Matlab-Simulink de la totalidad delmodelo desarrollado en la sección anterior, incluyendo el código que describe cada parte del modelo usadoen los diferentes bloques utilizados. En la sección 5, se introduce la estrategia de control a implementar queconsiste en el control predictivo basado en modelo (MPC). En la sección 6, se aborda el diseño del con-trolador para la central hidroeléctrica modelada. En primer lugar, se determina el punto de funcionamientodel sistema y se analizan las restricciones del mismo. En segundo lugar, se obtiene un modelo linealizadoen torno al punto de funcionamiento y se desarrolla el modelo de predicción para el controlador. En tercerlugar, se plantea el escenario de control (seguimiento de una trayectoria de referencia para la potencia totalgenerada) describiendo la función objetivo para el problema de control óptimo a resolver. Por último, semuestra la implementación del controlador en Simulink. En la sección 7, se presentan los resultados obte-nidos de las simulaciones realizadas con el controlador diseñado actuando sobre el modelo no lineal de lacentral. Para finalizar, en la sección 8 se mencionan las conclusiones del trabajo.
2
2. Energía hidráulicaSe denomina energía hidráulica a aquella que se obtiene del aprovechamiento de las energías cinética y
potencial de la corriente del agua, saltos de agua o mareas. Se puede transformar a muy diferentes escalas,existiendo desde hace siglos pequeñas explotaciones en las que la corriente de un río mueve un rotor depalas y genera un movimiento aplicado, por ejemplo, en molinos rurales. Sin embargo, la utilización mássignificativa la constituyen las centrales hidroeléctricas de presas.
Cuando el Sol calienta la Tierra, además de generar corrientes de aire, hace que el agua del mar seevapore, ascienda por el aire y se mueva hacia las regiones montañosas, para luego caer en forma de lluvia.Esta agua se puede colectar y retener mediante presas. Parte del agua almacenada se deja salir para quemueva los álabes de una turbina engranada con un generador de energía eléctrica.
2.1. Aprovechamiento de la energía hidráulicaComo ya se ha dicho, antiguamente se aprovechaba la energía del agua utilizándose ruedas hidráulicas
para moler trigo, sin embargo, la posibilidad de emplear esclavos y animales de carga retrasó su aplica-ción generalizada hasta el siglo XII. La energía hidroeléctrica debe su mayor desarrollo al ingeniero civilbritánico John Smeaton, que construyó por vez primera grandes ruedas hidráulicas de hierro colado.
La hidroelectricidad tuvo mucha importancia durante la Revolución Industrial, impulsando a las indus-trias textiles y del cuero y los talleres de construcción de máquinas a principios del siglo XIX. Aunque lasmáquinas de vapor ya estaban perfeccionadas, el carbón era escaso y la madera poco satisfactoria comocombustible. La energía hidráulica ayudó al crecimiento de las nuevas ciudades industriales que se crearonen Europa y América, hasta la construcción de canales a mediados del siglo XIX. Las presas y los canaleseran necesarios para la instalación de ruedas hidráulicas sucesivas cuando el desnivel era mayor de cincometros, no siendo posible la construcción de grandes presas de contención todavía. El bajo caudal de aguadurante el verano y el otoño, unido a las heladas en invierno, obligaron a sustituir las ruedas hidráulicas pormáquinas de vapor en cuanto se pudo disponer de carbón.
Las formas más frecuentemente utilizadas para explotar la energía hidráulica son:
Desvío del cauce de agua
El principio fundamental de esta forma de aprovechamiento hidráulico de los ríos se basa en el hechode que la velocidad del flujo de éstos es básicamente constante a lo largo de su cauce, el cual siempre esdescendente. Este hecho revela que la energía potencial no es íntegramente convertida en cinética comosucede en el caso de una masa en caída libre, la cual se acelera, sino que ésta es invertida en las llamadaspérdidas, es decir, la energía potencial se "pierde" en vencer las fuerzas de fricción con el suelo, en eltransporte de partículas, en formar remolinos, etc. Entonces, esta energía potencial podría ser aprovechadasi se pudieran evitar las llamadas pérdidas y hacer pasar al agua a través de una turbina. El conjunto de obrasque permiten el aprovechamiento de la energía anteriormente mencionada reciben el nombre de centralhidroeléctrica o hidráulica.
Interceptación de la corriente de agua
Este método consiste en la construcción de una presa de agua que retenga el cauce de agua causando unaumento del nivel del río en su parte anterior a la presa de agua, el cual podría eventualmente convertirseen un embalse. El dique establece una corriente de agua no uniforme y modifica la forma de la superficiede agua libre del río antes y después de éste, que toman forma de las llamadas curvas de remanso. Elestablecimiento de las curvas de remanso determinan un nuevo salto geodésico aprovechable de agua.
2.2. Central hidroeléctricaEn una central hidroeléctrica se utiliza energía hidráulica para la generación de energía eléctrica. En
general, estas centrales aprovechan la energía potencial que posee la masa de agua de un cauce naturalen virtud de un desnivel, también conocido como salto geodésico. El agua en su caída entre dos niveles
3
del cauce se hace pasar por una turbina hidráulica la cual transmite la energía a un generador donde setransforma en energía eléctrica.
2.2.1. Características de una central hidroeléctrica
Las dos características principales de una central hidroeléctrica, desde el punto de vista de su capacidadde generación de electricidad, son:
La potencia, que está en función del desnivel existente entre el nivel medio del embalse y el nivelmedio de las aguas debajo de la central, y del caudal máximo turbinable, además de las característicasde las turbinas y de los generadores usados en la transformación.
La energía garantizada en un lapso de tiempo determinado, generalmente un año, que está en funcióndel volumen útil del embalse y de la potencia instalada. La potencia de una central puede variar desdeunos pocos MW, como en el caso de las minicentrales hidroeléctricas, hasta miles de MW como esel caso de las grandes centrales de Paraguay, Brasil o China, por citar algunas.
Las centrales hidroeléctricas y las centrales térmicas (que usan combustibles fósiles) producen la energíaeléctrica de una manera muy similar. En ambos casos la fuente de energía es usada para impulsar una turbinaque hace girar un generador eléctrico, que es el que produce la electricidad. Una central térmica usa calorpara, a partir de agua, producir el vapor que acciona las paletas de la turbina, en contraste con la plantahidroeléctrica, la cual usa la fuerza del agua directamente para accionar la turbina.
2.2.2. Tipos de centrales hidroeléctricas
Según su concepción arquitectónica
Centrales al aire libre, al pie de la presa, o relativamente alejadas de ésta. Están conectadas por mediode una tubería en presión.
Centrales en caverna, generalmente conectadas al embalse por medio de túneles, tuberías en presión,o por la combinación de ambas.
Según su régimen de flujo
Centrales de agua fluyente: También denominadas centrales de filo de agua o de pasada, utilizanparte del flujo de un río para generar energía eléctrica. Operan en forma continua porque no tienencapacidad para almacenar agua, no disponen de embalse. Turbinan el agua disponible en el momento,limitadamente a la capacidad instalada. En estos casos, las turbinas pueden ser de eje vertical, cuandoel río tiene una pendiente fuerte, u horizontal cuando la pendiente del río es baja.
Centrales de embalse: Es el tipo más frecuente de central hidroeléctrica. Utilizan un embalse parareservar agua e ir graduando el agua que pasa por la turbina. Es posible generar energía durante todoel año si se dispone de reservas suficientes. Requieren una inversión mayor.
Centrales de regulación: almacenamiento del agua que fluye del río capaz de cubrir horas de consumo.
Centrales de bombeo o reversibles: Una central hidroeléctrica reversible es una central hidroeléctricaque además de poder transformar la energía potencial del agua en electricidad, tiene la capacidad dehacerlo a la inversa, es decir, aumentar la energía potencial del agua (por ejemplo subiéndola a unembalse) consumiendo para ello energía eléctrica. De esta manera puede utilizarse como un métodode almacenamiento de energía (una especie de batería gigante). Están concebidas para satisfacerla demanda energética en horas pico y almacenar energía en horas valle. Aunque lo habitual es queestas centrales turbinen/bombeen el agua entre dos embalses a distinta altura, existe un caso particularllamado centrales de bombeo puro donde el embalse superior se sustituye por un gran depósito cuyaúnica aportación de agua es la que se bombea del embalse inferior.
4
Según su altura de caída del agua
Centrales de alta presión: son las centrales de más de 200 m de caída del agua, siendo usadas lasturbinas tipo Pelton.
Centrales de media presión: son las centrales con caída del agua de 20 a 200 m, siendo dominante eluso de turbinas tipo Francis, aunque también se puedan usar tipo Kaplan.
Centrales de baja presión: son centrales con desniveles de agua de menos de 20 m, siendo usadas lasturbinas tipo Kaplan.
Centrales de muy baja presión: son centrales correspondientes con nuevas tecnologías, pues llega unmomento en el cuál las turbinas tipo Kaplan no son aptas para tan poco desnivel. Suelen situarse pordebajo de los 4 m.
2.2.3. Partes principales de una central hidroeléctrica
La Figura 2 muestra un ejemplo de central hidroeléctrica. En ella se ilustran las diferentes partes:
Embalse: se denomina embalse a la acumulación de agua producida por una obstrucción en el lechode un río o arroyo que cierra parcial o totalmente su cauce. La obstrucción del cauce puede ocurrirpor causas naturales como, por ejemplo, el derrumbe de una ladera en un tramo estrecho del río oarroyo, la acumulación de placas de hielo o las construcciones hechas por los castores, y por obrasconstruidas por el hombre para tal fin, como son las presas. Los embalses generados al construir unapresa pueden tener la finalidad de:
• Regular el caudal de un río o arroyo, almacenando el agua de los períodos húmedos para utili-zarlos durante los períodos más secos para el riego, para el abastecimiento de agua potable, parapermitir la navegación o para diluir contaminantes. Cuando un embalse tiene más de un fin, sele llama de usos múltiples.
• Contener los caudales extremos de las avenidas o crecidas.
• Crear espacios para esparcimiento y deportes acuáticos.
• Crear una diferencia de nivel para generar energía eléctrica, mediante una central hidroeléctrica.
Presa: la presa se encarga de mantener el agua en un lugar alto para garantizar que tenga fuerzasuficiente para mover las turbinas. Se construye a través de un río, arroyo o canal, para almacenar elagua a fin de derivarla o regular su curso fuera del cauce.
Tubería forzada: es la tubería que lleva el agua desde la presa hasta las turbinas.
Turbina: la turbina se encarga de hacer girar el generador cuando recibe la fuerza del agua. Es lamáquina destinada a transformar en movimiento giratorio la fuerza viva del fluido.
Generador, transformador y líneas eléctricas: el generador es el encargado de producir la electricidad.La electricidad viaja desde el generador hasta unos transformadores, donde se eleva la tensión parapoder transportar la electricidad hasta los centros de consumo.
5
Figura 2: Ejemplo de central hidroeléctrica
3. Modelado de una central hidroeléctrica
3.1. Descripción del sistemaEl sistema considerado es una central hidroeléctrica de bombeo o revesible compuesta por varios sub-
sistemas interconectados entre sí, como se ilustra en la Figura 3.La central considerada en este trabajo consta de tres lagos (L1, L2 y L3) y un río que está dividido en
seis tramos (R1, R2, R3, R4, R5 y R6) cada uno de los cuales termina en una presa equipada con una turbinapara producción de potencia (D1, D2, D3, D4, D5 y D6). Los lagos y los tramos del río están conectados porun conducto (U1), conductos equipados con una turbina (T1 y T2) y conductos equipados con una bomba yuna turbina (C1 y C2). El río es alimentado por un caudal de entrada (qin) y un caudal lateral (qtributary).
Para simplificar el modelado del sistema se hacen las siguientes suposiciones:
Los conductos están conectados al fondo de los lagos y al fondo del lecho del río.
La secciones transversales de los lagos y de los tramos del río son rectangulares.
El ancho del tramo varía linealmente a lo largo del tramo.
La pendiente del lecho del río es constante a lo largo de cada tramo.
3.2. Subsistemas
3.2.1. Modelo del tramo de río
Para el modelado de un fluido incompresible en un cauce unidireccional se hará uso de las ecuacionesgenerales de Saint-Venant. Estas son ecuaciones no-lineales de primer orden que describen la conserva-ción de la masa y la conservación de la cantidad de movimiento, en función de dos variables: la seccióntransversal del agua s(t,z) y el flujo q(t,z). A su vez, ambas se pueden escribir en función de la principalcoordenada espacial del río z, y del tiempo t. El sistema de ecuaciones de Saint-Venant se puede expresar
6
Figura 3: Visión de conjunto del sistema
de la siguiente forma (aunque también existen otras):
∂q(t,z)∂ z
+∂ s(t,z)
∂ t= 0 (3.1)
1g
∂
∂ t
(q(t,z)s(t,z)
)+
12g
∂
∂ z
(q2(t,z)s2(t,z)
)+
∂h(t,z)∂ z
+ I f (t,z)− I0(z) = 0 (3.2)
donde s(t,z) es el área de la sección transversal, q(t,z) es la descarga a través de la sección s, h(t,z) laprofundidad del agua, I f (t,z) la pendiente equivalente de fricción (tasa a la que se pierde energía debida ala resistencia del canal), I0(z) pendiente del canal y g la aceleración de la gravedad.
La tasa equivalente de fricción se define como:
I f (t,z) =q(t,z)2 · p(t,z)4/3
k2str · s(t,z)10/3 (3.3)
donde kstr es el coeficiente de Gauckler-Manning-Strickler1 y p(t,z) el perímetro mojado de la seccióntransversal. Asumiendo que la sección transversal del río es rectangular, como se muestra en la Figura 4, sepueden escribir las siguientes relaciones:
s(t,z) = w(z)h(t,z) (3.4)p(t,z) = 2h(t,z)+w(z) (3.5)
donde w(z) es el ancho del río.Por otro lado, los tramos de la central hidráulica considerada no están aislados, sino que (aparte de estar
conectados entre sí) también reciben flujos laterales. Por tanto, en cada sección es necesario considerar laposible contribución de un caudal lateral, que se supone entra en el tramo perpendicularmente a la direcciónprincipal del flujo. El caudal lateral, afecta sólo a la ecuación de conservación de la masa, la cual se modificade la siguiente manera:
∂q(t,z)∂ z
+∂ s(t,z)
∂ t= ql(z) (3.6)
donde ql(z) es la descarga lateral por unidad de longitud.
1Este coeficiente cambia de acuerdo al tipo de lecho del rió. En el modelo desarrollado se considera constante a lo largo del río.
7
w(0)
h(0,
t)
Z
s(z,t)h(
Z,t)
w(Z)
p(z,t)
q(z,t)
Figura 4: Sección de un tramo del río
Discretización espacial
Una forma simple de implementar un sistema de ecuaciones en derivadas parciales es discretizarlo envarios sistemas de ecuaciones diferenciales ordinarias, sustituyendo las derivadas espaciales por sus corres-pondientes diferencias finitas. Una ecuación diferencial ordinaria es una relación que contiene funciones deuna sola variable independiente, y una o más de sus derivadas con respecto a esa variable, mientras que, lasecuaciones del sistema en derivadas parciales involucran derivadas parciales de varias variables.
Por tanto, es aconsejable discretizar el tramo en N secciones a lo largo de la dirección del flujo, cadauna de longitud dz = Z
N , siendo Z la longitud total del tramo. Para describir correctamente la física delsistema, los pasos de discretización de las variables q y h se han solapado, dado que cada altura depende delas descargas previas y hacia adelante, y viceversa. En la Figura 5 se muestra el esquema de discretizacióndescrito. La aproximación dada por las ecuaciones en diferencias finitas es:
Figura 5: Discretización de un tramo de río
8
∂hi(t)∂ z
≈ (hi+1(t)−hi(t))dz
∂qi(t)∂ z
≈ (qi(t)−qi−1(t))dz
(3.7)
donde hi es la altura al comienzo de la sección i y qi es la descarga en el medio de la sección i. Tomando qiny qout como los caudales de entrada y salida del tramo respectivamente y combinando las ecuaciones (3.1)a (3.7), el sistema completo de ecuaciones diferenciales ordinarias resultante que debe ser implementadopara simular el comportamiento de cada tramo es (la dependencia del tiempo de las variables se omite):
∂h1∂ t =− 1
w1· q1−qin−ql1
dz/2
∂q1∂ t = q1
w1h1· ql1
dz/2 −2q1
w1h1· q1−qin
dz/2 +
[1
w1
(q1h1
)2−gw1h1
]· h2−h1
dz +
+gw1h1
(I0−
[q2
1(w1+2h1)4/3
k2str(w1h1)
10/3
])
i = 2, · · · ,N
∂hi∂ t =− 1
wi· qi−qi−1−qli
dz/2
∂qi∂ t = qi
wihi· qli
dz −2qiwihi· qi−qi−1
dz +
[1wi
(qihi
)2−gwihi
]· hi+1−hi
dz +
+gwihi
(I0−
[q2
i (wi+2hi)4/3
k2str(wihi)10/3
])∂hN+1
∂ t=− 1
wN+1· qout −qN
dz/2(3.8)
donde wi representa el ancho del río al comienzo de la sección i, wN+1 y hN+1 representan el ancho del ríoy la altura al final del tramo respectivamente y qli es el caudal lateral total de la sección i. La pendiente dellecho del río I0 se supone constante. Ya que el ancho de los tramos del río se supone varía linealmente, apartir de los valores de w1 y wN+1 se calcula el ancho wi como:
wi = w1 +(i−1)(wN+1−w1)
N(3.9)
3.2.2. Modelo del lago
Tomando de nuevo qin(t) y qout(t) como los caudales de entrada y salida del lago bajo consideración,respectivamente, el volumen de agua contenida en los lagos varía de acuerdo a la siguiente ecuación:
∂v(t)∂ t
= qin(t)−qout(t) (3.10)
Dado que se supone que la sección transversal del lago es rectangular, la ecuación (3.10) se puedeexpresar equivalentemente como:
∂h(t)∂ t
=qin(t)−qout(t)
S(3.11)
donde h(t) es el nivel de agua en el lago y S es la superficie del lago.
3.2.3. Modelo del conducto
El flujo dentro del conducto puede ser modelado usando la ley de Bernouilli. Asumiendo que la sec-ción del conducto es mucho más pequeña que la superficie del lago, el flujo del lago L1 al lago L2 puedeexpresarse como:
qU1(t) = SU1sign(hL2(t)−hL1(t)+hU1)√
2g |hL2(t)−hL1(t)+hU1 | (3.12)
donde hL1(t) y hL2(t) son los niveles de agua de lagos L1 y L2, hU1 es la diferencia de altura en el conducto,SU1 es la es la sección del conducto y g es la aceleración de la gravedad.
9
Tomando x = hL2(t)−hL1(t)+hU1 , la ecuación (3.12) se puede escribir como:
qU1(t) = SU1
√2gsign(x)
√|x| (3.13)
Dado que la función sign(x)√|x| no es derivable en x = 0, y por tanto, no apropiada para muchas
técnicas de control, se hace la siguiente aproximación:
sign(x)√|x|≈ x
(x2 + ε4)1/4 (3.14)
Nótese que para ε = 0 las dos funciones son equivalentes, mientras que tomando un valor suficiente-mente pequeño de ε se consigue una buena aproximación2.
3.2.4. Modelo de la turbina
Para cada turbina se asume que se puede controlar directamente la descarga de la misma. La potenciaproducida viene dada por la siguiente ecuación:
pt(t) = ktqt(t)∆ht(t) (3.15)
donde kt es el coeficiente de la turbina, qt(t) es la descarga de la turbina y ∆ht(t) es la diferencia de alturas.
3.2.5. Modelo de la bomba
Las bombas pueden ser modeladas de forma similar a las turbinas. La potencia absorbida por una bombaviene dada por:
pp(t) = kpqp(t)∆hp(t) (3.16)
donde kp es el coeficiente de la bomba, qp(t) es la descarga de la bomba y ∆hp(t) es la diferencia de alturas.
3.2.6. Modelo de conductos equipados con una turbina y una bomba
Los conductos C1 y C2 están equipados con una bomba y una turbina y, por tanto, pueden usarse lasecuaciones (3.15) y (3.16) para expresar la cantidad de potencia generada o absorbida. Sin embargo, lasturbinas y las bombas en el conducto no pueden funcionar a la vez, y así debe ser expresado en el problemade control. Dependiendo de la formulación del problema de control y del método usado para resolver elproblema, se pueden usar diferentes modelos. A continuación, se muestran diferentes posibilidades para elmodelado de estos conductos de flujo bidireccional. Se asume que el flujo puede ser determinado por elcontrolador.
Modelo discontinuo
Se denota por qC1(t) al flujo a través del conducto C1. Se asume que:
qC1(t)≥ 0 cuando C1 funciona como una turbina.
qC1(t)< 0 cuando C1 funciona como una bomba.
Usando esta convención no es necesario expresar explícitamente que C1 puede funcionar como una turbinao una bomba, alternativamente. La potencia producida puede ser expresada como:
pC1(t) = kC1 (qC1(t))qC1(t)∆hC1(t) (3.17)
2 1ε
corresponde a la derivada de la aproximación en x = 0.
10
donde ∆hC1(t) es la altura del conducto, que depende del nivel de agua en el lago L1 y en el tramo R1 y:
kC1 (qC1(t)) ={
ktC1 cuando qC1(t)≥ 0kpC1 cuando qC1(t)< 0 (3.18)
donde ktC1 es el coeficiente de la turbina y kpC1 es el coeficiente de la bomba. El flujo en C1 está limitado:
qC1(t) ∈ [−qC1 p,max,−qC1 p,min]∪ [qC1t,min,qC1t,max] (3.19)
donde los valores qC1 p,max, qC1 p,min, qC1 p,max y qC1 p,min son positivos (los subíndices son t para turbina y ppara bomba).
La ecuación (3.18) y la restricción expresada en (3.19) hacen al modelo del conducto C1 discontinuo, ypor tanto, no apropiado para muchas técnicas de control.
Modelo suavizado
La ecuación (3.18) puede ser escrita como:
kC1 (qC1(t)) =12((1+ sign(qC1(t)))ktC1 +(1− sign(qC1(t)))kpC1) (3.20)
y simplificada usando la siguiente aproximación:
sign(x)≈x
(x2 + ε2)1/2 (3.21)
La restricción (3.19) puede ser simplificada imponiendo:
qC1(t) ∈ [−qC1 p,max,qC1t,max] (3.22)
Modelo del flujo doble
Introduciendo dos variables independientes para expresar el flujo en C1 se puede obtener otro modelosimplificado de conductos equipados con turbina y bomba:
qC1p(t): flujo cuando C1 funciona como una bomba.
qC1t (t): flujo cuando C1 funciona como una turbina.
Asumiendo estas nuevas variables, ambas positivas, se puede escribir:
qC1(t) = qC1t (t)−qC1p(t) (3.23)
y
pC1(t) =(
ktC1qC1t (t)− kpC1qC1p(t))
∆hC1(t) (3.24)
Las restricciones expresadas en (3.19) pueden ser escritas en términos de qC1p(t) y qC1t (t):
qC1p(t) ∈ [qC1 p,min,qC1 p,max]
qC1t (t) ∈ [qC1t,min,qC1t,max] (3.25)
11
4. Implementación del sistema en SimulinkEl modelo de la central hidráulica descrito en la sección 3 se ha implementado en el entorno Simulink
de Matlab. Éste ha sido un paso fundamental para el desarrollo del trabajo dado que la obtención de unmodelo preciso era necesaria para diseñar una estrategia de control fiable que satisfaga los objetivos delproyecto y cumpla con las restricciones físicas impuestas en las variables del sistema.
4.1. Modelo de un tramo de ríoPara implementar en Simulink el sistema de 2N + 1 ecuaciones diferenciales ordinarias descrito en
(3.8) en la sección 3.2.1, se ha utilizado una aproximación computacionalmente eficiente que hace uso defunciones de Matlab con vectores, lo que permite calcular en un solo paso los vectores completos de alturasy caudales en cada sección del tramo. El esquema en Simulink en el que se implementa el comportamientode cada tramo se muestra en la Figura 6. Se ha creado una librería para los tramos que se puede reutilizarcambiando únicamente los parámetros de la máscara del bloque.
Figura 6: Implementación en Simulink de la librería de Tramos
El bloque Matlab Function que calcula las derivadas de las variables de estado del tramo implementadoen la librería, contiene el siguiente código:
12
function [dH,dQ,nl] = reachR(u,H,Q,w,p)
%#codegen
N = p(1);
g = p(2);
kstr = p(3);
I0 = p(4);
dz = p(5);
nl = p(6);
qin = u(1);
qout = u(2);
ql = zeros(N,1);
ql(nl) = u(3);
dH = zeros(N+1,1);
dQ = zeros(N,1);
dH(1) = -(1/w(1))*((Q(1)-qin-ql(1))/(dz/2));
dQ(1) = ((Q(1))/(w(1)*H(1)))*(ql(1)/(dz/2)) - ...
((2*Q(1))/(w(1)*H(1)))*((Q(1)-qin)/(dz/2)) + ...
(((1/w(1))*((Q(1)/H(1))^2)) - ...
g*w(1)*H(1))*((H(2)-H(1))/dz) + ...
g*w(1)*H(1)*(I0-(((Q(1)^2)*((w(1)+2*H(1))^(4/3)))/...
((kstr^2)*((w(1)*H(1))^(10/3)))));
dH(2:N) = -(1./w(2:N)).*((Q(2:N)-Q(1:N-1)-ql(2:N))./dz);
dQ(2:N) = ((Q(2:N))./(w(2:N).*H(2:N))).*(ql(2:N)./dz) - ...
((2.*Q(2:N))./(w(2:N).*H(2:N))).*((Q(2:N)-Q(1:N-1))./dz) + ...
(((1./w(2:N)).*((Q(2:N)./H(2:N)).^2)) - ...
g.*w(2:N).*H(2:N)).*((H(3:N+1)-H(2:N))./dz) + ...
g.*w(2:N).*H(2:N).*(I0-(((Q(2:N).^2).* ...
((w(2:N)+2.*H(2:N)).^(4/3)))./ ...
((kstr^2).*((w(2:N).*H(2:N)).^(10/3)))));
dH(N+1) = -(1/w(N+1))*((qout-Q(N))/(dz/2));
4.2. Modelo de una presaEl esquema en Simulink en el que se implementa el comportamiento de la turbina de una presa se mues-
tra en la Figura 7. Se ha creado una librería para las presas que se puede reutilizar cambiando únicamentelos parámetros de la máscara del bloque.
Figura 7: Implementación en Simulink de la librería de Presas
13
El bloque Matlab Function que calcula la potencia suministrada por la turbina de la presa, dada por laexpresión (3.15), contiene el siguiente código:
function [qD,pD] = damD(Hout,Hin,u,p)
ktD = p;
hDin = Hin(end);
hDout = Hout(1);
qD = u;
dhD = hDin - hDout;
pD = ktD*qD*dhD;
4.3. Modelo de los lagosEl esquema en Simulink en el que se implementa el comportamiento de los lagos se muestra en la Figura
8, para el caso concreto del Lago 1.
Figura 8: Implementación en Simulink del Lago 1
El bloque Matlab Function que calcula la derivada de la altura del lago dada por la ecuación (3.11),contiene el siguiente código:
function dhL1 = lake1(p,u)
%#codegen
SL1 = p(1);
qinL1 = p(2);
qU1 = u(1);
qT1 = u(2);
qC1p = u(3);
qC1t = u(4);
qin = qinL1 + qU1 + qC1p;
qout = qT1 + qC1t;
dhL1 = (qin - qout)/SL1;
Los Lagos 2 y 3 son completamente análogos al anterior.
4.4. Modelo del conducto U1
El esquema en Simulink en el que se implementa el comportamiento del caudal que atraviesa el con-ducto U1 se muestra en la Figura 9.
El bloque Matlab Function que calcula el caudal a través del conducto entre el Lago 1 y el Lago 2presentada en (3.12), contiene el siguiente código:
14
Figura 9: Implementación en Simulink del Conducto U1
function qU1 = ductU1(p,u)
%#codegen
SU1 = p(1);
g = p(2);
epsilon = p(3);
hU1 = p(4);
hL1 = u(1);
hL2 = u(2);
x = hL2 - hL1 + hU1;
qU1 = SU1*sqrt(2*g)*(x/((x^2+epsilon^4)^(1/4)));
4.5. Modelo de una turbinaEl esquema en Simulink en el que se implementa el comportamiento de una turbina se muestra en la
Figura 10, para el caso concreto de la Turbina 1.
Figura 10: Implementación en Simulink de la Turbina 1
El bloque Matlab Function que calcula la potencia suministrada por la turbina, dada por la expresión(3.15), contiene el siguiente código:
15
function pT1 = turbineT1(p,u)
%#codegen
hT1 = p(1);
ktT1 = p(2);
hL1 = u(1);
hR2T1 = u(2);
qT1 = u(3);
dhT1 = hT1 + hL1 - hR2T1;
pT1 = ktT1*qT1*dhT1;
4.6. Modelo de un conducto equipado con turbina y bombaEl esquema en Simulink en el que se implementa el comportamiento de la turbina/bomba se muestra en
la Figura 11, para el caso concreto del conducto C1.
Figura 11: Implementación en Simulink de un conducto C1 equipado con turbina y bomba
Se ha elegido para la implementación el modelo de flujo doble descrito en la sección 3.2.6. El bloqueMatlab Function que calcula la potencia suministrada por la turbina o consumida por la bomba, dada por laexpresión (3.24), contiene el siguiente código:
function pC1 = t_pumpC1(p,u)
%#codegen
hC1 = p(1);
ktC1 = p(2);
kpC1 = p(3);
hL1 = u(1);
hR1C1 = u(2);
qC1t = u(3);
qC1p = u(4);
dhC1 = hC1 + hL1 - hR1C1;
pC1 = (ktC1*qC1t-kpC1*qC1p)*dhC1;
4.7. Modelo del sistema completoLa central hidráulica completa considerada en este proyecto se muestra en el esquema de Simulink de
la Figura 12.
16
Figura 12: Esquema en Simulink de la central hidroeléctrica completa
El sistema consta de seis tramos y seis presas y de dos subsistemas adicionales compuestos por lagos,conductos, conductos equipados con turbina y conductos equipados con turbina y bomba. Las entradas delsistema completo son los caudales que circulan por cada uno de los elementos del sistema y las salidas sonlas potencias de cada uno de los equipos que producen o consumen energía en el conjunto.
Los elementos básicos del sistema son los elementos que se han denominado tramo (Reach) y presa(Dam), presentados en la Figura 13.
Figura 13: Esquema en Simulink de tramo y presa
Las entradas a un tramo son los caudales que lo atraviesan (caudal de entrada en la sección inicial,caudal de salida en la sección final y caudal lateral (que entra en una determinada sección del tramo) y lassalidas son los vectores de alturas y de caudales en las secciones del tramo y el índice correspondiente a lasección que de entrada del caudal lateral, en caso de que lo hubiera.
La presa o dique, está a su vez equipada con una turbina que produce electricidad en función de ladiferencia de alturas entra la sección final del tramo anterior y la sección inicial del tramo siguiente, asícomo del caudal que circula por el dique. Sus entradas son los vectores de alturas de los tramos anteriory posterior y el caudal que circula por el dique, y sus salidas son el caudal y la potencia producida en laturbina.
Subsistema 1
Los elementos que conforman el subsistema uno pueden verse en la Figura 14.
17
Figura 14: Esquema en Simulink del Subsistema 1
Subsistema 2
Los elementos que conforman el subsistema dos pueden verse en la Figura 15.
Figura 15: Esquema en Simulink del Subsistema 2
4.8. Parámetros del sistemaLos parámetros seleccionados para la implementación del modelo se muestran en las Tablas 1 a 9.
18
Tabla 1: Parámetros relacionados con el Subsistema 1
Símbolo Descripción Unidades Valor
SL1 Superficie Lago 1 m2 10000
SL2 Superficie Lago 2 m2 5000
hL1min Nivel mínimo Lago 1 m 10,5
hL1max Nivel máximo Lago 1 m 13,5
hL2min Nivel mínimo Lago 2 m 5,5
hL2max Nivel máximo Lago 2 m 8,5
qinL1Caudal entrada Lago 1 m3
s 5
qinL2Caudal entrada Lago 2 m3
s 5
hU1 Diferencia de alturas Conducto U1 m 5
SU1 Sección transversal Conducto U1 m2 6
ktT1 Coeficiente Turbina T1J
m4 8000
ktC1 Coeficiente Turbina C1J
m4 8000
kpC1 Coeficiente Bomba C1J
m4 14000
hT1 Diferencia de alturas Conducto T1 m 223
hC1 Diferencia de alturas Conducto C1 m 200
qT1min Caudal mínimo Turbina T 1m3
s 0
qT1max Caudal máximo Turbina T1m3
s 20
qC1t,min Caudal mínimo Conducto C1 en modo turbina m3
s 0
qC1t,max Caudal máximo Conducto C1 en modo turbina m3
s 10
qC1 p,min Caudal mínimo Conducto C1 en modo bomba m3
s 0
qC1 p,max Caudal máximo Conducto C1 en modo bomba m3
s 5
19
Tabla 2: Parámetros relacionados con el Subsistema 2
Símbolo Descripción Unidades Valor
SL3 Superficie Lago 3 m2 10000
hL3min Nivel mínimo Lago 3 m 6
hL3max Nivel máximo Lago 3 m 9
qinL3Caudal entrada Lago 3 m3
s 10
ktT2 Coeficiente Turbina T2J
m4 8000
ktC2 Coeficiente Turbina C2J
m4 8000
kpC2 Coeficiente Bomba C2J
m4 14000
hT2 Diferencia de alturas Conducto T2 m 233
hC2 Diferencia de alturas Conducto C2 m 250
qT2min Caudal mínimo Turbina T2m3
s 0
qT2max Caudal máximo Turbina T2m3
s 20
qC2t,min Caudal mínimo Conducto C2 en modo turbina m3
s 0
qC2t,max Caudal máximo Conducto C2 en modo turbina m3
s 10
qC2 p,min Caudal mínimo Conducto C2 en modo bomba m3
s 0
qC2 p,max Caudal máximo Conducto C2 en modo bomba m3
s 5
Tabla 3: Parámetros relacionados con el Tramo R1 del río
Símbolo Descripción Unidades Valor
I01 Pendiente del lecho Tramo R1 - 0,0025
L1 Longitud Tramo R1 m 10000
w1,1 Ancho al inicio Tramo R1 m 30
wN+1,1 Ancho al final Tramo R1 m 50
N1 Nº de secciones discretización Tramo R1 - 20
ktD1 Coeficiente turbina Dique D1J
m4 8000
hR1min Nivel mínimo Dique D1 m 14,5
hR1max Nivel máximo Dique D1 m 17,5
qD1min Caudal mínimo Turbina Dique D1m3
s 5
qD1max Caudal máximo Turbina Dique D1m3
s 300
qin Caudal de entrada Tramo R1m3
s 200
LC1
Distancia desde el comienzo del Tramo R1
hasta la entrada lateral del Conducto C1m 5000
20
Tabla 4: Parámetros relacionados con el Tramo R2 del río
Símbolo Descripción Unidades Valor
I02 Pendiente del lecho Tramo R2 - 0,0015
L2 Longitud Tramo R2 m 8000
w1,2 Ancho al inicio Tramo R2 m 40
wN+1,2 Ancho al final Tramo R2 m 45
N2 Nº de secciones discretización Tramo R2 - 20
ktD2 Coeficiente turbina Dique D2J
m4 8000
hR2min Nivel mínimo Dique D2 m 16,5
hR2max Nivel máximo Dique D2 m 19,5
qD2min Caudal mínimo Turbina Dique D2m3
s 5
qD2max Caudal máximo Turbina Dique D2m3
s 300
LT1
Distancia desde el comienzo del Tramo R2
hasta la entrada lateral del Conducto T1m 4000
Tabla 5: Parámetros relacionados con el Tramo R3 del río
Símbolo Descripción Unidades Valor
I03 Pendiente del lecho Tramo R3 - 0,002
L3 Longitud Tramo R3 m 6000
w1,3 Ancho al inicio Tramo R3 m 40
wN+1,3 Ancho al final Tramo R3 m 55
N3 Nº de secciones discretización Tramo R3 - 20
ktD3 Coeficiente turbina Dique D3J
m4 8000
hR3min Nivel mínimo Dique D3 m 21,5
hR3max Nivel máximo Dique D3 m 24,5
qD3min Caudal mínimo Turbina Dique D3m3
s 5
qD3max Caudal máximo Turbina Dique D3m3
s 300
LtributaryDistancia desde el comienzo del Tramo R3
hasta la entrada del caudal lateralm 3000
qtributary Caudal lateral m3
s 30
21
Tabla 6: Parámetros relacionados con el Tramo R4 del río
Símbolo Descripción Unidades Valor
I04 Pendiente del lecho Tramo R4 - 0,0015
L4 Longitud Tramo R4 m 8000
w1,4 Ancho al inicio Tramo R4 m 55
wN+1,4 Ancho al final Tramo R4 m 45
N4 Nº de secciones discretización Tramo R4 - 20
ktD4 Coeficiente turbina Dique D4J
m4 8000
hR4min Nivel mínimo Dique D4 m 17,5
hR4max Nivel máximo Dique D4 m 20,5
qD4min Caudal mínimo Turbina Dique D4m3
s 5
qD4max Caudal máximo Turbina Dique D4m3
s 300
LC2
Distancia desde el comienzo del Tramo R4
hasta la entrada lateral del Conducto C2m 4000
Tabla 7: Parámetros relacionados con el Tramo R5 del río
Símbolo Descripción Unidades Valor
I05 Pendiente del lecho Tramo R5 - 0,0025
L5 Longitud Tramo R5 m 6000
w1,5 Ancho al inicio Tramo R5 m 50
wN+1,5 Ancho al final Tramo R5 m 60
N5 Nº de secciones discretización Tramo R5 - 20
ktD5 Coeficiente turbina Dique D5J
m4 8000
hR5min Nivel mínimo Dique D5 m 13,5
hR5max Nivel máximo Dique D5 m 16,5
qD5min Caudal mínimo Turbina Dique D5m3
s 5
qD5max Caudal máximo Turbina Dique D5m3
s 300
LT2
Distancia desde el comienzo del Tramo R5
hasta la entrada lateral del Conducto T2m 3000
22
Tabla 8: Parámetros relacionados con el Tramo R6 del río
Símbolo Descripción Unidades Valor
I06 Pendiente del lecho Tramo R6 - 0,002
L6 Longitud Tramo R6 m 8000
w1,6 Ancho al inicio Tramo R6 m 60
wN+1,6 Ancho al final Tramo R6 m 80
N6 Nº de secciones discretización Tramo R6 - 20
ktD6 Coeficiente turbina Dique D6J
m4 8000
hR6min Nivel mínimo Dique D6 m 11,5
hR6max Nivel máximo Dique D6 m 14,5
qD6min Caudal mínimo Turbina Dique D6m3
s 10
qD6max Caudal máximo Turbina Dique D6m3
s 300
hD6out Nivel después del Dique D6 m 2
Tabla 9: Datos Generales
Símbolo Descripción Unidades Valor
kstr Coeficiente de Gauckler-Manning-Strickler m1/3
s 30
g aceleración de la gravedad m 9,81
23
5. Control predictivo basado en modeloEl control predictivo basado en modelo (MPC) tiene como objetivo resolver de forma efectiva problemas
de control y automatización de procesos que se caractericen por presentar un comportamiento dinámicocomplicado, multivariable, y/o inestable. La estrategia de control utiliza el modelo matemático del procesoa controlar para predecir el comportamiento futuro de dicho sistema, y en base a este comportamientofuturo calcular la señal de control. El control predictivo integra disciplinas como el control óptimo, controlestocástico, control de procesos con retardo de tiempo, control multivariable y control con restricciones.
El MPC se enmarca dentro de los controladores óptimos, es decir, aquellos en los que las actuacionesresponden a la optimización de un criterio. El criterio a optimizar, o función de coste, está relacionado conel comportamiento futuro del sistema, que se predice gracias a un modelo dinámico del mismo, denomina-do modelo de predicción. El intervalo de tiempo futuro que se considera en la optimización se denominahorizonte de predicción. Dado que el comportamiento futuro del sistema depende de las actuaciones quese aplican a lo largo del horizonte de predicción, son éstas las variables de decisión respecto a las quese optimiza el sistema. La aplicación de estas actuaciones sobre el sistema conduce a un control en bucleabierto. La posible discrepancia entre el comportamiento predictivo y el comportamiento real del sistemacrean la necesidad de imponer cierta robustez al sistema incorporando realimentación del mismo. Esta reali-mentación se consigue gracias a la técnica del horizonte deslizante que consiste en aplicar las actuacionesobtenidas durante un periodo de tiempo, tras el cual se muestrea el estado del sistema y se resuelve unnuevo problema de optimización. De esta manera, el horizonte de predicción se va deslizando a lo largo deltiempo.
Una de las propiedades más atractivas del MPC es su formulación abierta, que permite la incorporaciónde distintos tipos de modelos de predicción, sean lineales o no lineales, monovariables o multivariables, y laconsideración de restricciones sobre las señales del sistema. Esto hace que sea una estrategia muy utilizadaen diversas áreas del control. El MPC es una de las pocas técnicas que permiten controlar sistemas conrestricciones incorporando éstas en el propio diseño del controlador.
En resumen, las principales ventajas de MPC son:
Formulación en el dominio del tiempo, lo cual le permite ser una técnica flexible e intuitiva.
Permite tratar con sistemas lineales y no lineales, monovariables y multivariables, utilizando la mismaformulación para los algoritmos del controlador.
La ley de control responde a criterios de optimización.
Permite la incorporación de restricciones en el diseño del controlador.
Como desventajas, cabe mencionar:
Requiere el conocimiento de un modelo dinámico del sistema suficientemente preciso.
Requiere un algoritmo de optimización, por lo que solo se puede implementar por medio de uncomputador.
Requiere un alto coste computacional, lo que hace difícil su aplicación a sistemas rápidos.
La Figura 16 resume la estrategia del MPC. En cada instante t y haciendo uso del modelo del proceso3
se predicen las futuras salidas para un determinado horizonte de predicción Ny. Estas salidas predichas,y(t +k | t)4 para k = 1, · · · ,Ny dependen de los valores conocidos hasta el instante t de las entradas y salidaspasadas y de las señales de control futuras u(t +k | t) para k = 0, · · · ,Ny−1 que se pretenden mandar al sis-tema y que son las que se quieren calcular. El conjunto de señales de control futuras se calcula optimizandoun determinado criterio en el que se pretende mantener el proceso lo más próximo posible a la trayectoriade referencia yr(t + k). Este criterio suele tomar la forma de una función cuadrática de los errores entre lasalida predicha y la trayectoria de referencia futura, incluyendo en muchos casos el esfuerzo de control. Siel criterio es cuadrático, el modelo lineal y no existen restricciones se puede obtener una solución explícita,
3El MPC puede usar como modelo de predicción cualquier descripción, como por ejemplo la respuesta impulsional, la respuesta
ante escalón, función de transferencia o el espacio de estados.4valor de la variable y en el instante t + k calculado en el instante t.
24
en otro caso se debe usar un método iterativo de optimización. La señal de control u(t | t) es enviada al pro-ceso mientras que las siguientes señales de control calculadas son desechadas, puesto que en el siguienteinstante de muestreo ya se conoce y(t +1) y se repite el proceso con este nuevo valor y todas las secuenciasson actualizadas. Se calcula por tanto u(t+1 | t+1) (que en principio sera diferente al u(t+1 | t) al disponerde nueva información), haciendo uso del concepto de horizonte deslizante.
Figura 16: Estrategia de Control Predictivo basado en Modelo (MPC)
5.1. Modelo de predicción en el espacio de estadosSea un sistema dado por su descripción en el espacio de estados en tiempo discreto:
x(k+1) = Ax(k)+Bu(k)+v(k)y(k) = Cx(k)+n(k) (5.1)
donde x(k) es el vector de estado de dimensión nx, u(k) es el vector de entradas manipulables de dimensiónnu, y(k) es el vector de salidas de dimensión ny, v(k) es un vector de variables aleatorias que representan elruido en el proceso y n(k) es un vector de variables aleatorias que representan el ruido en las medidas.
Se define el vector de incrementos de la señal de control como ∆u(k) = u(k)−u(k− 1). Teniendo encuenta esto, es conveniente expresar (5.1) en función del incremento:[
x(k+1)u(k)
]=
[A B0 Inu×nu
][x(k)
u(k−1)
]+
[B
Inu×nu
]∆u(k)+
[Inx×nx
0
]v(k)
y(k) =[
C 0][ x(k)
u(k−1)
]+n(k) (5.2)
Haciendo el siguiente cambio de variable:
x(k) =[
x(k)u(k−1)
](5.3)
se tiene que:
x(k+1) = Mx(k)+N∆u(k)+Hv(k)y(k) = Qx(k)+n(k) (5.4)
donde los valores de las matrices M, N y Q se obtienen comparando (5.2) con (5.4).En adelante, se supone que los vectores v(k) y n(k) son ruido blanco, por lo tanto:
E [v(k)] = 0∀kE [n(k)] = 0∀k
E[v(k) ·v(k)T ]= Rv
E[n(k) ·n(k)T ]= Rn (5.5)
25
donde las matrices de covarianza Rv y Rn son diagonales y se suponen constantes ∀k.La salida en el instante k+ j se puede calcular recursivamente de la siguiente manera:
y(k+1) = Qx(k+1)+n(k+1) = QMx(k)+QN∆u(k)+QHv(k)+n(k+1)y(k+2) = Qx(k+2)+n(k+2) = QMx(k+1)+QN∆u(k+1)+QHv(k+1)+n(k+2)
= QM2x(k)+QMN∆u(k)+QN∆u(k+1)+QMv(k)+QHv(k+1)+n(k+2)...
...
y(k+ j) = QM jx(k)+j−1
∑i=0
QM j−i−1N∆u(k+ i)+j−1
∑i=0
QM j−i−1Hv(k+ i)+n(k+ j) (5.6)
De esta forma, la predicción de la salida y en el instante k+ j es:
y(k+ j |k) = E [y(k+ j)]
= QM jE [x(k)]+j−1
∑i=0
QM j−i−1N∆u(k+ i)+j−1
∑i=0
QM j−i−1HE [v(k+ i)]+E [n(k+ j)] (5.7)
Teniendo en cuenta (5.5), finalmente resulta:
y(k+ j |k) = QM j ˆx(k)+j−1
∑i=0
QM j−i−1N∆u(k+ i) (5.8)
donde ˆx(k) es la estimación del estado del sistema en el instante k. Puesto en forma matricial se tiene:y(k+1 |k)y(k+2 |k)
...y(k+Ny |k)
=
QMQM2
...QMNy
ˆx(k)+
QN 0 · · · 0
QMN QN · · · 0...
.... . .
...QMNy−1N QMNy−2N · · · QN
∆u(k)∆u(k+1)
...∆u(k+Ny−1)
(5.9)
De forma compacta:
y = F ˆx(k)+G4u = f+G4u (5.10)
Esto significa que la evolución futura del sistema se puede descomponer en dos: la respuesta libre f,debida exclusivamente al estado del sistema en el instante actual manteniendo la entrada constante y larespuesta forzada G4u, debida exclusivamente a las acciones de control.
La matriz G es triangular inferior por bloques y puede obtenerse de forma sencilla calculando cadaelemento no nulo como:
Gi j = QMi− jN (5.11)
5.2. Función objetivoHabitualmente, la función objetivo a minimizar es:
J (∆u(k)) =Ny
∑j=1
((y(k+ j |k)−yr(k+ j))T Qy
j (y(k+ j |k)−yr(k+ j))+
+∆u(k+ j−1)T Ruj∆u(k+ j−1)
)(5.12)
26
donde Q j y R j son matrices diagonales de ponderación definidas positivas cuyos elementos penalizanlos errores cuadráticos entre la predicción de la salida y la referencia yr y los esfuerzos de control ∆urespectivamente a lo largo del horizonte de predicción.
Sustituyendo (5.10) en (5.12), la función objetivo se puede escribir como:
J (∆u) =124uT H4u+b4u+ f0 (5.13)
siendo
H = 2(GT QG+Ru)
b = 2(f−yr)T QyG
f0 = (f−yr)Qy(f−yr)T (5.14)
Las matrices Qyy Ru son matrices diagonales por bloques, cuyos elementos diagonales son las matrices
Qyj y Ru
j . Si no se consideran restricciones en las variables, se puede obtener una solución analítica para laseñal de control óptima, igualando a cero el gradiente de J:
∇J (4u) = 0⇒∆u∗ =−H−1bT = (GT QyG+Ru
)−1GT Qy(yr− f) =−KMPC · (f−yr) (5.15)
Si se consideran restricciones en las variables, como límites en la amplitud de las señales de controlu(k) ∈ [u,u]∀k o en la salida y(k) ∈ [y,y]∀k, estas pueden describirse, respectivamente, para el horizontede predicción Ny de la siguiente forma:
1uu≤ T4u+1uu(k−1)≤ 1uu1yy≤ f+G4u≤ 1yy (5.16)
donde 1u es una matriz de dimensión (nu×Ny)× nu formada por Ny matrices identidad Inu×nu , 1y es unamatriz de dimensión (ny×Ny)× ny formada por Ny matrices identidad Iny×ny y T es una matriz triangularinferior por bloques cuyos bloques no nulos son matrices matrices identidad Inu×nu . De forma compacta, sepuede escribir:
R4u≤ c (5.17)
donde:
R =
T−TG−G
c =
1uu−1uu(k−1)−1uu+1uu(k−1)
1yy− f−1yy+ f
(5.18)
De esta forma, el problema de optimización que se debe resolver es un problema de optimización cua-drática (QP):
min∆u(k) J (∆u(k)) =124uT H4u+b4u+ f0
s.a. R4u≤ c (5.19)
27
6. Diseño del controlador
6.1. Punto de funcionamiento y restriccionesPara diseñar el controlador MPC hace falta linealizar el modelo de la central entorno a un punto de
funcionamiento para conseguir un sistema equivalente a (5.4). Se utilizará como punto de funcionamientoel punto de equilibrio del sistema. El conjunto de ecuaciones definido por (3.11) y (3.8) definen un modelono-lineal de la forma x = f (x,u), donde:
x =[hL1 hL2 hL3 QR1 HR1 · · · QR6 HR6
]TQR j =
[q1, j · · · qN j , j
]THR j =
[h1, j · · · hN j+1, j
]Tu =
[qT1 qC1 qT2 qC2 qD1 qD2 qD3 qD4 qD5 qD6
]T (6.1)
son los vectores de estado y entrada, respectivamente. Para encontrar el punto de equilibrio (x0,u0) seresuelve mediante métodos numéricos el sistema de ecuaciones dado por:
x = f (x,u) = 0 (6.2)
Como puede verse en las Tablas 1 a 8, el modelo de centra hidroeléctrica considerado tiene restriccionesen los niveles de los lagos y al final de cada tramo, así como en los caudales de las diferentes turbinas ybombas del sistema:
hLimin ≤ hLi ≤ hLimax i = 1,2,3hR jmin ≤ hN j+1, j ≤ hR jmax j = 1, · · · ,6
qTkmin ≤ qTk ≤ qTkmax k = 1,2qDlmin ≤ qDl ≤ qDlmax l = 1, · · · ,6
−qCm p,max ≤ qCm ≤ qCmt,max m = 1,2 (6.3)
donde se ha considerado el modelo suavizado descrito en la sección 3.2.6 para los conductos Cm. Para poderplantear las restricciones en la forma (5.16), se define el vector de salida del sistema como:
y =[pTotal hLi hLi hLi hN1+1,1 hN2+1,2 hN3+1,3 hN4+1,4 hN5+1,5 hN6+1,6
]T (6.4)
donde pTotal es la potencia total en MW generada por la centra, dada por:
pTotal =(
10e−6)(
pT1 + pC1 + pT2 + pC2 + pD1 + pD2 + pD3 + pD4 + pD5 + pD6
)(6.5)
6.2. Linealización y modelo de predicciónCon el fin de obtener una descripción del sistema en el espacio de estados como la planteada en (5.4),
se linealiza el modelo en torno al punto de funcionamiento (x0,u0) calculando las matrices Jacobianas:
Ai j =∂ fi(x,u)
∂x j
∣∣∣∣x=x0,u=u0
Bi j =∂ fi(x,u)
∂u j
∣∣∣∣x=x0,u=u0
(6.6)
de tal forma que:
δ x(t) = Aδx(t)+Bδu(t) (6.7)
siendo δx(t) = x(t)− x0 y δu(t) = u(t)−u0. En el vector de salida elegido en (6.4) aparecen las poten-cias de las diferentes turbinas/bombas del sistema. Dichas potencias son función tanto de alturas como decaudales, por lo tanto se tiene que y = h(x,u). Dado que en la descripción en el espacio de estados buscadala salida solo depende del estado, es decir y(k) = Cx(k), es necesario ampliar el vector de estado con los
28
caudales de entrada. Para ello, se añaden al sistema actuadores con una dinámica lineal de primer orden, deforma que:
δu(t) =1
τs+1δu′(t) (6.8)
donde u′ es el nuevo vector de entrada y τ = 30 es la constante de tiempo de la función de transferencia. Laecuación (6.8) puede escribirse:
δ u(t) =1τ
(δu′(t)−δu(t)
)(6.9)
Usando (6.7) y (6.9) se tiene:[δ x(t)δ u(t)
]=
[A B0 − 1
τInu×nu
][δx(t)
δu
]+
[0
1τ
Inu×nu
]δu′(t) (6.10)
Haciendo el siguiente cambio de variable:
δz(t) =[
δx(t)δu
](6.11)
se tiene que:
δ z(t) = A′δz(t)+B′δu′(t) (6.12)
donde los valores de las matrices A′ y B′ se obtienen comparando (6.12) con (6.10). Ahora puede obtenerse:
Ci j =∂hi(x,u)
∂δ z j
∣∣∣∣x=x0,u=u0
(6.13)
lo que proporciona la relación δy(t) = Cδz(t), con δy(t) = y(t)−y0 = y(t)−h(x0,u0), para completar ladescripción en el espacio de estados. En resumen:
δ z(t) = A′δz(t)+B′δu′(t)δy(t) = Cδz(t) (6.14)
El siguiente paso es discretizar el sistema lineal obtenido. Para ello se utiliza el método de mantenedoresde orden cero de las entradas y un tiempo de muestreo de Ts = 10s, resultando:
δz(k+1) = A∗δz(k)+B∗δu′(k)δy(k+1) = C∗δz(k) (6.15)
Ahora se puede obtener el modelo de predicción correspondiente a la central hidroeléctrica aplicando a(6.15) los pasos descritos por (5.2) y (5.9). Dado que se han utilizado actuadores con una dinámica linealde primer orden, las restricciones en los caudales de entrada u se mantienen para las nuevas entradas u′,sin necesidad de hacer ninguna modificación adicional. De esta forma, los vectores con los límites para lasrestricciones son:
u =[qT1max qC1max qT2max qC1max qD1max qD2max qD3max qD4max qD5max qD6max
]T −u0
u =[qT1min −qC1max qT2min −qC2max qD1min qD2min qD3min qD4min qD5min qD6min
]T −u0
y =[2 · p0
Total hL1max hL2max hL3max hR1max hR2max hR3max hR4max hR5max hR6max]T −y0
y =[0 hL1min hL2min hL3min hR1min hR2min hR3min hR4min hR5min hR6min
]−y0 (6.16)
29
6.3. Escenario de control y función objetivoEl escenario de control que se plantea es el seguimiento de una trayectoria de referencia para la potencia
total generada por la central hidroeléctrica durante 24 horas. Se supone que la referencia cambia cada 30minutos, manteniéndose constante en dichos intervalos. El problema de control optimo a resolver es:
min J =47
∑k=0
γk
ˆ tk+1
tk
∣∣∣∣∣pre fk −
8
∑i=1
pi(t)
∣∣∣∣∣dt
s.a. x = f (x,u) (6.17)
donde tk = 1800k, γ es un vector de ponderación con el precio de la electricidad para cada intervalo de 30minutos, pre f
k es la referencia de potencia a seguir y pi(t) son las potencias de las diferentes turbinas/bombasdel sistema. Sabiendo que x = f (x,u) se usa para obtener el modelo de predicción de y, teniendo en cuenta
el concepto de horizonte deslizante Ny, que8∑
i=1pi(t) = pTotal y dada la elección del vector de salida y
planteada en (6.4), el problema (6.17) puede escribirse de la siguiente forma:
min∆u(k) J (∆u(k)) =Ny
∑j=1
(y(k+ j |k)−yre f (k+ j)
)TΓ j(y(k+ j |k)−yre f (k+ j)
)(6.18)
donde
y(k+ j |k) =[pTotal hLi hLi hLi hN1+1,1 hN2+1,2 hN3+1,3 hN4+1,4 hN5+1,5 hN6+1,6
]T(k+ j |k)
(6.19)
es el vector de salida predicho para el instante k+ j,
yre f (k+ j) =[pre f (k+ j) 0 0 0 0 0 0 0 0 0
]T (6.20)
es el vector de referencia en el instante k+ j y Γ j es una matriz diagonal, cuya diagonal es el vector
dΓ j =[γ(k+ j) 0 0 0 0 0 0 0 0 0
](6.21)
siendo γ(k+ j) el precio de la electricidad en el instante k+ j. Sustituyendo (5.9) en (6.18), se tiene:
J (∆u(k)) =Ny
∑j=1
((F j ˆx(k)+G j4u−yre f (k+ j)
)TΓ j(F j ˆx(k)+G j4u−yre f (k+ j)
)+
+∆u(k+ j−1)T Ruj∆u(k+ j−1)
)(6.22)
donde F j y G j son la fila j de las matrices F y G respectivamente. Se añade el término de ponderación delas señales de control con el objetivo de evitar cambios muy bruscos en los caudales de entrada, donde lamatriz Ru
j es una matriz identidad. Operando sobre (6.22) resulta:
J (∆u(k)) =124uT H4u+b4u+ f0 (6.23)
siendo
H = 2(GTΓG+Ru
)
b = 2(f−yre f )TΓG
f0 = (f−yre f )Γ(f−yre f )T (6.24)
Incorporando las restricciones al problema de la forma (5.16)-(5.18) utilizando (6.16), se concluye queel problema de optimización (6.17) a resolver para el seguimiento de trayectorias de referencia de la po-tencia generada por la central hidráulica es completamente análogo al problema de optimización cuadráticaQP planteado en (5.19):
min∆u(k) J (∆u(k)) =124uT H4u+b4u+ f0
s.a. R4u≤ c (6.25)
30
6.4. Implementación del controlador en SimulinkLa implementación del controlador puede verse en la Figura 17, donde se ha utilizado un bloque S-
Function que resuelve el problema de optimización mediante el algoritmo Interior Point Convex de lafunción quadprog de Matlab. Cabe destacar que el controlador obtenido mediante la linealización del sis-tema en torno al punto de equilibrio se aplica sobre el modelo no lineal de la central hidroeléctrica. Paraconseguir un correcto funcionamiento, se ha planteado un esquema en el que se recalcula la linealización encada paso k de muestreo en torno a los valores actuales del estado x y la entrada u, siempre que la derivadade la potencia total generada en el instante k−1 no esté por encima 1e−3. Es decir, si la salida no varia deun instante a otro, se conserva la última linealización realizada.
Figura 17: Implementación en Simulink del sistema de control para la central hidroeléctrica
Se supone que el vector de estados se puede medir completamente. Para simular ruido en las medidas,en cada instante k se añade al vector de estado ruido blanco de media x(k) y desviación típica 0,01. Se haelegido el tiempo de muestreo Ts = 10s para reducir el tiempo de simulación. El horizonte de predicciónutilizado es Ny = 60s. A priori puede parecer un horizonte muy pequeño. La razón fundamental para elegireste valor, es que dada la dimensión del sistema (249 estados y 10 entradas), si los tamaños de las ma-trices del modelo de predicción son demasiado grandes, pueden generar problemas de memoria y ademásralentizar considerablemente el optimizador.
Para comprobar el desempeño del controlador frente a perturbaciones se han modificado los caudales qiny qtributray de manera que estos varían cada 30 minutos en un intervalo de±5% en torno al valor definido enlas Tablas 3 y 5 respectivamente. Además, a cada caudal en el instante k, se le añade también ruido blancode media q(k) y desviación típica 0,1.
31
7. Simulación y resultados
7.1. Ponderación de la función objetivoLa Figura 18 muestra los valores del precio de la electricidad γ durante 24 horas utilizados en la simu-
lación.
0 5 10 15 2062
63
64
65
66
67
68
69
70Precio de la electricidad para 24 horas
γ [E
UR
/ W
h]
t [h]
Figura 18: Precio de la electricidad γ durante 24 horas
7.2. Referencia de potenciaLa Figura 19 muestra la referencia de potencia pre f a seguir en 24 horas. Esta referencia varía en torno
a la potencia nominal5 de la central y01 = 1749,5MW en un intervalo de ±1%.
0 5 10 15 201730
1735
1740
1745
1750
1755
1760
1765
1770Referencia de potencia
pre
f [M
W]
t [h]
Figura 19: Referencia de potencia pre f a seguir en 24 horas
5potencia generada en el equilibrio
32
7.3. Perturbaciones
7.3.1. Caudal de entrada en Tramo R1
La Figura 20 muestra la variación del caudal de entrada en el Tramo R1 durante 24 horas, antes de añadirruido blanco.
0 5 10 15 20190
192
194
196
198
200
202
204
206
208
210Caudal de entrada Tramo R1
qin
[m
3 / s
]
t [h]
Figura 20: Variación del caudal de entrada en el Tramo R1 durante 24 horas
7.3.2. Caudal de entrada en Tramo R3
La Figura 21 muestra la variación del caudal lateral en el Tramo R3 durante 24 horas, antes de añadirruido blanco.
0 5 10 15 2028.5
29
29.5
30
30.5
31
31.5Caudal lateral Tramo R3
qtr
ibuta
ry [m
3 / s
]
t [h]
Figura 21: Variación del caudal lateral en el Tramo R3 durante 24 horas
33
7.4. Seguimiento de la potencia de referenciaLa Figura 22 muestra la potencia total generada por la central durante 24 horas. Se observa como el
controlador permite seguir la referencia impuesta pese a las perturbaciones y los ruidos en las medidas. LaFigura 23 muestra el detalle de la potencia durante un intervalo de cambio en la referencia.
0 5 10 15 201730
1735
1740
1745
1750
1755
1760
1765
1770Potencia total generada
pT
ota
l [M
W]
t [h]
Figura 22: Potencia total generada por la central durante 24 horas
7.8 8 8.2 8.41738
1739
1740
1741
1742
1743
1744
1745
1746Potencia total generada (detalle)
pT
ota
l [M
W]
t [h]
Referencia
Potencia generada
Figura 23: Detalle de la potencia total generada por la central
34
Las Figuras 24 a 26 muestran los niveles de los Lagos.
0 5 10 15 2010
10.5
11
11.5
12
12.5
13
13.5
14
Nivel del Lago L1
hL1 [
m]
t [h]
Figura 24: Nivel del Lago L1
0 5 10 15 205
5.5
6
6.5
7
7.5
8
8.5
9
Nivel del Lago L2
hL1 [
m]
t [h]
Figura 25: Nivel del Lago L2
0 5 10 15 205.5
6
6.5
7
7.5
8
8.5
9
9.5
Nivel del Lago L3
hL3 [
m]
t [h]
Figura 26: Nivel del Lago L3
35
Las Figuras 27 a 32 muestran los niveles al final de los Tramos R1 a R6.
0 5 10 15 2014
14.5
15
15.5
16
16.5
17
17.5
18
Nivel al final del Tramo R1
hN
+1,1
[m
]
t [h]
Figura 27: Nivel al final del tramo R1
0 5 10 15 2016
16.5
17
17.5
18
18.5
19
19.5
20
Nivel al final del Tramo R2
hN
+1,2
[m
]
t [h]
Figura 28: Nivel al final del tramo R2
0 5 10 15 2021
21.5
22
22.5
23
23.5
24
24.5
25
Nivel al final del Tramo R3
hN
+1,3
[m
]
t [h]
Figura 29: Nivel al final del tramo R3
36
0 5 10 15 2017
17.5
18
18.5
19
19.5
20
20.5
21
Nivel al final del Tramo R4
hN
+1,4
[m
]
t [h]
Figura 30: Nivel al final del tramo R4
0 5 10 15 2013
13.5
14
14.5
15
15.5
16
16.5
17
Nivel al final del Tramo R5
hN
+1,5
[m
]
t [h]
Figura 31: Nivel al final del tramo R5
0 5 10 15 2011
11.5
12
12.5
13
13.5
14
14.5
15
Nivel al final del Tramo R6
hN
+1,6
[m
]
t [h]
Figura 32: Nivel al final del tramo R6
37
Las Figuras 33 a 42 muestran los caudales de entrada del sistema.
0 5 10 15 20
9.9
9.95
10
10.05
10.1
10.15
10.2
10.25
10.3
Caudal Turbina T1
qT
1 [
m3 /
s]
t [h]
Figura 33: Caudal Turbina T1
0 5 10 15 20−0.15
−0.1
−0.05
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
Caudal Turbina/Bomba C1
qC
1 [
m3 /
s]
t [h]
Bomba
Turbina
Figura 34: Caudal Turbina/Bomba C1
0 5 10 15 20
9.9
9.95
10
10.05
10.1
10.15
10.2
10.25
10.3
Caudal Turbina T2
qT
2 [
m3 /
s]
t [h]
Figura 35: Caudal Turbina T2
38
0 5 10 15 20−0.2
−0.1
0
0.1
0.2
0.3
0.4
0.5
0.6
Caudal Turbina/Bomba C2
qC
2 [
m3 /
s]
t [h]
Bomba
Turbina
Figura 36: Caudal Turbina/Bomba C2
0 5 10 15 20199.99
199.992
199.994
199.996
199.998
200
200.002
200.004
200.006
Caudal Turbina D1
qD
1 [
m3 /
s]
t [h]
Figura 37: Caudal Turbina D1
0 5 10 15 20209.99
209.995
210
210.005
Caudal Turbina D2
qD
2 [
m3 /
s]
t [h]
Figura 38: Caudal Turbina D2
39
0 5 10 15 20239.985
239.99
239.995
240
240.005
240.01
240.015
240.02
Caudal Turbina D3
qD
3 [
m3 /
s]
t [h]
Figura 39: Caudal Turbina D3
0 5 10 15 20239.985
239.99
239.995
240
240.005
240.01
240.015
240.02
Caudal Turbina D4
qD
4 [
m3 /
s]
t [h]
Figura 40: Caudal Turbina D4
0 5 10 15 20249.985
249.99
249.995
250
250.005
250.01
250.015
250.02
Caudal Turbina D5
qD
5 [
m3 /
s]
t [h]
Figura 41: Caudal Turbina D5
40
0 5 10 15 20249.99
249.995
250
250.005
250.01
250.015
250.02
250.025
Caudal Turbina D6
qD
6 [
m3 /
s]
t [h]
Figura 42: Caudal Turbina D6
La Figura 43 muestra el nivel de las secciones de todos los tramos después de 24 horas.
1 2 3 4 5 60
5
10
15
20
25Niveles de las secciones en todos los tramos depués de 24 horas
h [m
]
Tramos R
Figura 43: Nivel de las secciones de todos los tramos después de 24 horas
41
8. ConclusionesObservando los resultados de simulación del controlador diseñado actuando sobre el modelo no lineal
de la central, se concluye que la estrategia de control predictivo basado en modelo es satisfactoria. Lapotencia total generada por la central consigue seguir la trayectoria de referencia utilizada para un períodode 24 horas, como se muestra en las Figuras 22 y 23. Cabe destacar los siguientes aspectos:
Aunque se ha utilizado una estrategia de control lineal, el esquema de relinealización planteado per-mite conseguir el objetivo, a pesar de que el controlador actúa sobre un modelo no lineal.
Se han introducido en la simulación perturbaciones al sistema en forma de variación aleatoria delcaudal de entrada al tramo R1 y del caudal lateral del tramo R3, como se muestra en las Figuras 20 y21. Asimismo, se ha introducido ruido blanco en la medida del vector de estados. Incluso bajo estascircunstancias, el controlador es capaz de conseguir su objetivo.
El problema de control óptimo planteado permite obtener señales de control que hacen que se cum-plan todas las restricciones del sistema, tanto en las entradas (caudales de turbinas y bombas) comoen los estados (alturas de los lagos y al final de cada tramo del río).
La forma de la trayectoria de referencia utilizada se asemeja a una curva de demanda tipo en la cual hayperíodos de tiempo donde la demanda se sitúa por debajo de la potencia nominal y otros donde se sitúa porencima. En los resultados de simulación se observa que el sistema se comporta como cabe esperar:
En períodos donde la referencia está por debajo de la potencia nominal, los tramos C1 y C2 funcionancomo bombas (Figuras 34 y 36) produciéndose el llenado de los lagos (Figuras 24 a 26).
En períodos donde la referencia está por encima de la potencia nominal, los tramos C1 y C2 funcionancomo turbinas (Figuras 34 y 36) produciéndose el vaciado de los lagos (Figuras 24 a 26).
Este comportamiento coincide con el concepto de almacenamiento mencionado anteriormente, en el cualen períodos con poca demanda este tipo de centrales transforman energía mecánica (potencia consumidapor las bombas durante un cierto intervalo de tiempo) en energía potencial (agua almacenada en los lagos)para luego liberarla en períodos con alta demanda.
42
Bibliografía[1] Hierarchical and distributed model predictive control of large-scale systems, European FP7 STREP
project HD-MPC, Sept. 1 2008 - Dic. 31 2011, Online: http://www.ict-hd-mpc.eu/
[2] Savorgnan, C. Diehl, M. Control benchmark of a Hydro Power Valley, Version 1, Jul. 8 2011, Online:http://www.ict-hd-mpc.eu/benchmarks/HD-MPC_hydropower_valley_public_benchmark.zip
[3] Camacho, E. F. Bordons, C. Model Predictive Control, Advanced textbooks in control and signal pro-cesing, Springer-Verlang London Limited, 2004
[4] Petrone, F. Model Predictive Control of a Hydro Power Valley, PhD Tesis, Politecnico di Milano, 2009-2010
43
Top Related