Modelado y control de sistemas de climatización de edificios en TRNSYS

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 Trabajo Fin de Grado Grado en Ingeniería de Tecnologías Industriales Modelado y Control de Sistemas de Climatización de edificios en TRNSYS Autor: Alberto Oliveros Iglesias Tutor: Daniel Limón Marruedo Joaquín García Ordoñez Dpto. Ingeniería de Sistemas y Automática Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2020 Trabajo Fin de Grado Grado en Ingeniería de Tecnologías Industriales Modelado y Control de Sistemas de Climatización de edificios en TRNSYS Autor: Alberto Oliveros Iglesias Tutor: Daniel Limón Marruedo Catedrático de Universidad Joaquín García Ordoñez Dpto. Ingeniería de Sistemas y Automática Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2020 Trabajo Fin de Grado: Modelado y Control de Sistemas de Climatización de edificios en TRNSYS Autor: Alberto Oliveros Iglesias Tutor: Daniel Limón Marruedo Joaquín García Ordoñez 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: A mi familia A mis amigos A mis profesores Agradecimientos El objetivo de estas lineas es agredecer a todas aquellas personas me han apoyado y ayudado a conseguir este objetivo durante estos cuatro años. En primer lugar, me gustaría agredecer a mis padres todo el esfuerzo que han realizado para que yo haya podido cumplir este sueño, por toda la ayuda recibida de su parte y todo el apoyo brindado para la consecución de mis objetivos. Por otro lado, querría dedicar unas palabras a mi tutor Daniel Limón quien me ha dado la oportunidad de poder realizar el trabajo bajo su tutela y las facilidades que me ha propuesto en todo mommento para la realización de este proyecto. No quiero olvidarme tampoco del que ha sido uno de los grandes artífices de estre trabajo, Joaquín García, que desde el primer momento ha mostrado plena disposición para sacar este trabajo adelante, sin importar ni el día ni la hora, y siempre ha presentado interés en como se ha ido desarrollando el proyecto en cada momento. Por último y no menos importante, no quiero olvidarme de la que ha sido mi "familia" durante mi estancia en esta ciudad. A ellos queiero agradecerles todos esos momentos que hemos vivido durante nuestra etapa universitaria. En especial quiero hacer mención a dos personas, mi Hermanita y Pablo, que en este último año se han convertido en un pilar fundamental para mi. Alberto Oliveros Iglesias Escuela Técnica Superior de Ingeniería Sevilla, 2020 III XÍndice 7 Control MPC 61 7.1 Observador de estados 61 7.2 Control MPC 62 7.3 Aplicación al caso de estudio 63 7.4 Resultados obtenidos 68 8 Ampliación del edificio 73 8.1 Definición del sistema 73 8.2 Incorporación de equipos 77 8.3 Identificación del sistema 80 8.4 Control LQR 80 8.5 Control MPC 85 9 Conclusiones y Ampliaciones 91 Bibliografía 93 Índice de Tablas 2.1 Salidas del bloque del edificio (Type 56) 12 2.2 Entradas del bloque del edificio (Type 56) 13 3.1 Coeficientes cide la CAP de la enfriadora. 22 3.2 Coeficientes pi j de la Potencia del intercambiador . 22 3.3 Tabla resumen de las variables empleadas por cada equipo. 22 5.1 Resumen de parámetros a considerar. 40 5.2 Variables definidas en TRNBuil 41 5.3 Coeficientes utilizados para el cálculo de la densidad. 45 XI Índice de Figuras 1.1 Evolución del consumo energético por sectores en Espala (Fuente IEA) [17] 2 2.1 Configuración del edificio 6 2.2 Asignación de nombre a cada sala 7 2.3 Ventana para la modificación de características de una sala. 8 2.4 Características del muro exterior norte de la Sala 1 9 2.5 Numeración resultante tras la definición de los muros 10 2.6 Esquema conceptual de conexiones 11 2.7 Conexiones realizadas entre el edificio y el scope 12 2.8 Conexiones realizadas entre el tiempo y el edificio 14 2.9 Resultado de realizar todas las conexiones 14 2.10 Radiación incidente en las caras Norte y Sur 15 2.11 Radiación incidente en las caras Este y Oeste 15 2.12 Temperaturas de la sala 1 hasta la 4 considerando la influencia atmosférica 16 2.13 Temperaturas de la sala 5 hasta la 7 considerando la influencia atmosférica 16 3.1 Objeto de estudio en el Capítulo 3 19 3.2 Esquema del circuito refrigerador 20 3.3 Esquema del circuito representando las variables 23 3.4 Esquema del conexionado a realizar 27 3.5 Resultados obtenidos tras la modificación de las Aperturas (Salas 1 a 4) 28 3.6 Resultados obtenidos tras la modificación de las Aperturas (Salas 5 a 7) 28 3.7 Resultados obtenidos tras la modificación del fcp(Salas 1 a 4) 29 3.8 Resultados obtenidos tras la modificación del fcp (Salas 5 a 7) 30 4.1 Esquema entradas y salidas del sistema 32 4.2 Esquema obtención de datos en TRNSYS 32 4.3 Entradas del sistema de la Prueba 2 33 4.4 Salidas del sistema de la Prueba 2 34 4.5 Entradas del sistema de la Prueba 8 35 4.6 Salidas del sistema de la Prueba 8 35 4.7 Resultados obtenidos para decidir el número de estados a considerar 36 4.8 Resultados en la identificación de la Prueba 2 37 4.9 Resultados en la identificación de la Prueba 8 38 4.10 Resultados tras la aplicación de escalones unitarios en el espacio de estados 38 5.1 Cuadro resultante tras la implementación de personas en TRNSYS 40 5.2 Cuadro resultante tras la implementación de las luces en TRNSYS 41 5.3 Cuadro resultante tras la implementación de los ordenadores en TRNSYS 41 5.4 Esquema del conexionado a realizar para introducir el simulador 47 5.5 Situación de la Sala 2 durante la simulación 47 5.6 Situación de la Sala 4 durante la simulación 48 XIII XIV Índice de Figuras 5.7 Situación del Pasillo durante la simulación 48 6.1 Esquema del observador utilizado 52 6.2 Diagramas de bloques del control LQR 54 6.3 Esquema de implementación del controlador LQR en TRNSYS 57 6.4 Resultados obtenidos de la implementacion del LQR (Salas 1 hasta 4) 57 6.5 Resultados obtenidos de la implementacion del LQR (Salas 5 hasta 7) 58 6.6 Resultados obtenidos de la implementacion del LQR (Aperturas 1 hasta 4) 58 6.7 Resultados obtenidos de la implementacion del LQR (Aperturas 5 hasta 7 y fcp) 59 7.1 Bloque del objeto del controlador (ctrl) 67 7.2 Esquema de implementación del controlador MPC en TRNSYS 68 7.3 Resultados obtenidos de la implementacion del MPC (Salas 1 hasta 4) 69 7.4 Resultados obtenidos de la implementacion del MPC (Salas 5 hasta 7) 69 7.5 Resultados obtenidos de la implementacion del MPC (Aperturas 1 hasta 4) 70 7.6 Resultados obtenidos de la implementacion del MPC (Aperturas 5 hasta 7 y fcp) 70 7.7 Resultados obtenidos de la implementacion del MPC (Temperatrua del agua de salida de la enfriadora) 71 8.1 Configuración del edificio con tres plantas 74 8.2 Esquema conexionado para el modelo de trs plantas 74 8.3 Radiación incidente sobre las caras Norte y Sur del nuevo edificio. 75 8.4 Radiación incidente sobre las caras Este y Oeste del nuevo edificio. 75 8.5 Evolución de las temperaturas en el edificio de tres plantas (Salas 1 a 4) 76 8.6 Evolución de las temperaturas en el edificio de tres plantas (Salas 5 a 7) 76 8.7 Cambios que se producen por la apertura de válvulas en el nuevo modelo (Salas 1 a 4) 78 8.8 Cambios que se producen por la apertura de válvulas en el nuevo modelo (Salas 5 a 7) 78 8.9 Cambios que se producen por la variación de fcp en el nuevo modelo 79 8.10 Resultados de la identificación realizada (Entradas 1 a 10) 81 8.11 Resultados de la identificación realizada (Entradas 11 a 21) 82 8.12 Temperaturas y Referencias con un controlador LQR (Planta Baja) 83 8.13 Temperaturas y Referencias con un controlador LQR (Planta Media) 84 8.14 Temperaturas y Referencias con un controlador LQR (Planta Alta) 84 8.15 Comparativa entre los valores de fcp y de la temperatura a la salida del chiller (Control LQR) 85 8.16 Temperaturas y Referencias con un controlador MPC (Planta Baja) 87 8.17 Temperaturas y Referencias con un controlador MPC (Planta Media) 88 8.18 Temperaturas y Referencias con un controlador MPC (Planta Alta) 88 8.19 Comparativa entre los valores de fcp y de la temperatura a la salida del chiller (Control MPC) 89 1 Introducción y Estado del Arte Hoy en día se está llevando a cabo un gran crecimiento del uso de la energía disponible en el mundo, lo que supone grandes impactos medioambientales como pueden ser la destrucción de la capa de ozono o el calentamiento global, entre otros. Se puede estimar que la energía consumida tanto por edificios comerciales como residenciales ha sufrido un aumento entre el 20 % y 40 % [ 13 ], lo cual puede deberse en una primera idea al crecimiento de la población mundial que se viene sucediendo hasta ahora. Es por ello, que uno de los grandes retos del futuro erradica en la eficiencia de los edificios. Todos estos puntos van a ser tratados un poco más en detalle a continuación. 1.1 Motivación y desarrollo del proyecto Como se ha enunciado, el rápido crecimiento de la población mundial está provocando un crecimiento en la demanda energética lo cual produce algunos inconvenientes medioambientales como el de la destrucción de la capa de ozono, el calentamiento global o el cambio climático. Algunas predicciones muestran que este crecimiento se va a seguir produciendo debido entre otros factores al gran crecimiento económico que están experimentando algunas naciones. Uno de los hechos que lo muestran es el gran crecimiento que ha sufrido la demanda de energía primaria situado por encima de la tasa de crecimiento poblacional. Además, se tienen datos de que el PIB de las naciones ha aumentado, así como la renta per cápita de los ciudadanos. Con estos datos se puede demostrar que el crecimiento económico está ligado en gran parte a una mayor demanda energética, lo cual conlleva un aumento del CO2 per cápita que se produce, más concretamente un 5 % (según International Energy Agency) [17]. Este aumento de las mejoras en las condiciones económicas de los ciudadanos es la base en el incremento de la demanda energética en los edificios, la cual ha aumentado (alrededor de un 16 % ) hasta niveles equiparables con sectores como la industria y el transporte. También hay que tener en cuenta que la demanda energética que se produce depende de factores como la localización geográfica o el tamaño de las zonas residenciales. A modo de ejemplo en la Figura 1.1 se muestra la evolución que ha sufrido la demanda energética en cada uno de los sectores principales a lo largo de los años. En ella se puede observar que a pesar de la crisis económica que afectó a nuestro país, la demanda energética apenas sufrió descensos comparables con los que si que sufrió la industria y el trasnporte lo cual también se ve un poco influido por el clima que se desarrolla en España. Debido a este fenómeno que causa un aumento de las emisiones de CO2 que se generan en el entorno de la climatización de viviendas y edificios, la eficiencia energética de estos se ha convertido en uno de los grandes retos mundiales. Debido a esto, se ha producido un aumento en el uso de sistemas HVAC, del inglés Heating, Ventilation and Air Conditioning. Hay datos que indican que alrededor del 50 % de la energía consumida en las viviendas corresponde al acondicionamiento climático. Además se estima que el crecimiento de la energía demanda por este hecho se dispare en los próximos años. Hasta ahora, se ha relacionado el gran consumo energético con las viviendas. Sin embargo, en los edificios no residenciales como oficinas u hoteles, se produce una gran cantidad de demanda energética para poder llevar a cabo los servicios que prestan. Casi la mitad de la energía que se consume en las estancias no residenciales puede ser atribuible a oficinas y comercios, seguidas de hoteles, restaurantes, colegios y hospitales. De todas esta energía, se estima que la mitad es usada para el acondicionaniento o lo que es lo mismo, es consumida por equipos HVAC. 1 2Capítulo 1. Introducción y Estado del Arte Figura 1.1 Evolución del consumo energético por sectores en Espala (Fuente IEA) [17]. Según [ 13 ] se puede estimar que un 20 % de la energía consumida por un país es destinada a la climatización, de forma que se puede afirmar que una gran parte de energía demandada es destinada al uso de equipos HVAC y por tanto es importante lograr en intentar conseguir una mayor eficiencia en estos equipos, de forma que con menos energía demanda, sean capaces de poder aportar la potencia necesaria para acondicionar térmicamente espacios de diversa índole. Para ello, se pueden desarrollar sistemas de control para estos equipos con el objetivo de lograr una gestión eficiente, lo cual es una tarea compleja ya que genera la necesidad de disponer de modelos dinámicos que permitan emular de forma realista el comportamiento de los edificios. Debido a este fenómeno, en los últimos años se han producido grandes avances en este campo de estudio [ 1 ] y se han desarrollado numerosos entornos de simulación como TRNSYS, Dymola o Energy Plus, los cuales permiten modelar el comportamiento térmico de edificios e implementar en ellos distintos equipos HVAC. Entre los softwares existentes, se va a emplear TRSNSY, que es un paquete de software desarrollado por la universidad de Wisnconsin-Madison, el cual ha sido elegido por factores como: • Dicho software incorpora una sencilla interfaz gráfica que permite la creación sencilla de modelos y permite obtener el transitorio que sufren las temperaturas modificadas por la implementación de sistemas de climatización. • TRNSYS permite incorporar parámetros muchos más realistas a la simulación como son los acoplamientos por flujos de aire entre distintos espacios, mal funcionamiento de equipos, ganancias provocadas por personas o apertura de puertas entre otras. • Permite gran extensibilidad con programas externos como MATLAB ® , lo cual es francamente útil para poder implementar distintos tipos de controladores sobre el sistema. • Permite hacer uso de librerías DLL externas, las cuales permiten representar componentes que hayan sido creados y modelados por el usuario. Por todo esto, y algunas desventajas que presentan el resto de programas enunciados, TRNSYS ha sido el entorno elegido. En él se deben implementar los distintos equipos que modelan la máquina enfriadora, la cual debe ser diseñada teniendo en cuenta la eficiencia de la operación realizada. Uno de los problemas a los que debe enfrentarse es la complejidad de crear modelos matemáticos de cada uno de los equipos que componen la máquina enfriadora, pues los fabricantes no detallan mucha información acerca de esta. Para ello, durante este proyecto se van a ir analizando cada uno de los pasos que deben ser realizados para lograr una correcta implementación de un sistema inspirado en uno de los laboratorios que se encuentra en la Escuela Técnica Superior de Ingeniería de la Universidad de Sevilla. Cuando ya se disponga de dicho modelo, se van a aplicar sobre éste distintas técnicas de control que permitan una mejora de la eficiencia energética. 1.1 Motivación y desarrollo del proyecto 3 Mas concretamente se emplearán las técnicas LQR (Linear Quadratic Regulator) y MPC (Model Predictive Control), el cual se ha probado que puede llegar a ahorrar hasta un total de un 30 % si se compara con otras técnicas de control algo más clásicas [ 11 ]. Sin embargo, una de las desventajas de este tipo de técnicas de control es el gran coste computacional que conlleva, aunque se puede llevar hasta un segundo plano debido a que posee la capacidad de optimizar un balance entre el coste energético y las condiciones de confort. A modo de resumen, el objetivo de este trabajo es la implementación de distintas técnicas de control con el fin de observar la influencia que tiene en los parámetros que marcan la eficiencia de equipos que permiten la refrigeración de edificios. 2 Definición del sistema La ciencia puede divertirnos y fascinarnos, pero es la Ingeniería la que cambia el mundo. Isaac Asimov En este primer capítulo, se persigue la definición del sistema sobre el que se basa toda la realización de dicho trabajo. Es por ello que se trata de una parte fundamental, ya que de esta dependerá la calidad de los resultados obtenidos. En las próximas páginas, se podrá visualizar y caracterizar el sistema, y posteriormente, se ilustrará como se deben implementar correctamente todas las características necesarias para su construcción. Este capítulo se apoya fundamentalmente en la herramienta TRNSYS, un software comercial desarrollado por la universidad de Wisconsin-Madison (EE.UU) orientado, entre otras cosas, a la simulación energética de edificios, considerando aspectos realistas como pérdidas por flujo de aire, ganancias internas del sistema (personas, luces...), así como pérdidas propias del edificio [18]. Una vez se ha definido el entorno de trabajo, ya se está en disposición de comenzar el diseño del edificio sobre el cual se apoyará el trabajo restante. 2.1 Diseño del edificio El diseño del edificio que se va a realizar está basado en la construcción de los laboratorios ubicados en la Escuela Técnica Superior de Ingeniería de la Universidad de Sevilla, más concretamente, se puede asimilar a los que se encuentran en el Laboratorio 1, planta 2. Para el diseño, si lo asemejamos a los laboratorios, el esquema estará compuesto por un total de seis despachos, tres de ellos con ubicación norte y los otros tres con ubicación sur, además de un pasillo común para lograr la comunicación entre despachos. Para poder ilustrar de una forma más clara la configuración, se ha realizado el esquema visible en la Figura 2.1, el cual ha sido creado con la herramienta online de Google SketchUp [16]. Como se puede observar en la configuración propuesta, todas las estancias están equipadas de una ventana para poder aportar al edificio un intercambio de calor provocado por la radiación solar, y que podrá ser considerado como una fuente de perturbaciones en el modelo. Con respecto a las dimensiones de los despachos se han diseñado todos de la misma dimensión con un área de unos 16 m2 , siendo cuadrados de 4m por cada dimensión, la altura de los techos, se sitúa a unos 3.5m sobre el nivel del suelo. Todos ellos, como ya se ha enunciado anteriormente, disponen de una ventana localizada en las paredes exteriores (norte y sur) con una dimensión de 1.75 mde ancho por 1.20 mde alto. En relación al pasillo, como es lógico, tendrá una dimensión a lo largo de 12 m2 , mientras que su ancho es algo más pequeño que el de los despachos, siendo este de 3m . Al igual que los despachos, el techo está situado a una altura de 3.5m . Además, también dispone de una ventana situada en la cara oeste del edificio, con las mismas dimensiones que las instaladas en los despachos. Para finalizar con el diseño preliminar del edificio, se ha de saber que para facilitar comunicación entre los despachos y el pasillo se disponen de unas puertas de dimensión de 2.07 m de alto por 0.91 m de ancho. Al 5 12 Capítulo 2. Definición del sistema con el símbolo . Entonces si se hace click sobre este y luego uno se dirige a realizar la unión verá que cuando se sitúa encima de uno de los bloques previamente insertados, aparece una nube 8 puntos. A pesar de la existencia de tantos puntos de unión, lo cual se ha realizado para permitir mejores visualizaciones en el archivo, su significado es el mismo. Lo que sí es realmente importante es la dirección en la que se realiza dicho enlace, puesto que primero se debe pulsar sobre el bloque que queramos que nos aporte unas ciertas salidas sobre el segundo bloque que las tomará como entradas. Comenzamos, en primer lugar, con la conexión más sencilla, siendo esta la que se da entre el bloque que simula al edificio y el scope insertado. Una vez se han realizado la unión (desde el edificio hasta el scope), aparecerá una ventana emergente, donde se muestran por un lado las salidas del bloque del edificio y por el otro, las entradas al scope. Las entradas al scope son variables (número variables en cada eje y limites valores del eje Y) que pueden ser sometidas a modificación en el interior del propio bloque. Las variables de salida del bloque del edificio son las mostradas en la Tabla 2.1, de las cuales solo serán de interés las temperaturas de las habitaciones. Para realizar la conexión basta con unir los nombres de las temperaturas con las que han sido creadas automáticamente en el bloque de la representación. Tabla 2.1 Salidas del bloque del edificio (Type 56). Definición Notación Temperatura de cada una de las salas T1,T2...T7 Calor intercambiado por cada una de las salas QSens1,QSens2...QSens7 Para enlazar de una forma adecuada las variables y obtener una buena representación, se deben realizar las conexiones mostradas en la Figura 2.7. De esta forma, ya se tienen correctamente conectados ambos bloques y se puede continuar con la otra conexión. Figura 2.7 Conexiones realizadas entre el edificio y el scope. Para poder establecer todas las conexiones, ya solo es necesario realizar la correspondiente al bloque atmosférico y el edificio. Para ello, se vuelve a habilitar el icono de las conexiones y se hace la unión desde el bloque del tiempo hacia el edificio. Al igual que antes, aparece una ventana emergente donde entre las variables que se pueden visualizar, hay que destacar las siguientes: • La temperatura de bulbo seco (Dry bulb temperature), es la temperatura medida con un termómetro, el cual se encuentra protegido de fenómenos como la radiación solar y la humedad. Esta debe ser conectada tanto a la temperatura ambiente (TAMB), como a la temperatura de la superficie del suelo (TSGRD), ambas disponible entre las multiples entradas del edificio. • La temperatura efectiva del cielo (Effective sky temperature) es un concepto bastante abstracto, pero aquí basta con realizar la conexión entre esta y la entrada el edificio que representa lo mismo, denominada con TSKY • La humedad relativa del ambiente viene dada por la variable Percent relative humidity, que como se sabe representa la cantidad de agua o vapor que está presente en el aire. En el bloque del edificio 2.5 Simulación en TRNSYS 13 existe una variable que representa lo mismo, denominada RELHUMAMB, por lo que ambas deben ser conectadas. • Para la determinación de la radiación incidente en el edificio, la base de datos dispone de las distintas variaciones que se producen en el ángulo cenital (es el formado con la normal del punto donde nos encontramos) y el acimutal (es el que forma la proyección del Sol con la dirección sur). Estas dos variables están representadas por Solar zenith angle y Solar azimuth angle respectivamente, mientras que en el edificio están definidas como AZEN y AAZM. • Para el cálculo de la radiación total horizontal, se ha de especificar la cantidad de radiación incidente (Total horizontal radiation) y la Horizontal beam radiation, que la cantidad de radiación solar recibida antes de ser dispersada por la atmósfera. Además, también se distingue el parámetro que determina el ángulo de incidencia horizontal. Para ello, estas tres variables deben ser conectadas con IT_H_0_0, IB_H_0_0 yAI_H_0_0 respectivamente. •Al igual que sucede con la cara horizontal, el resto de orientaciones (Norte, Sur, Este y Oeste) deben ser caracterizadas con las correspondientes cantidades totales de radiación solar, así como con la beam radiation y el ángulo de incidencia sobre cada superficie. En este caso, las superficies aparecen numeradas desde la 1 hasta la 4 , por lo que hay que realizar la conexión sabiendo la relación entre los números y las superficies, del tal forma que: –La superficie 1se corresponde con la orientación norte (N_180_90). –La superficie 2con la cara este del edificio (S_0_90). –La superficie 3con la orientación oeste (E_270_90). –Por último, la 4con la orientada hacia el sur (W_90_90). • El último parámetro a tener en cuenta la reflectancia de la tierra (Ground reflectance), que es la fracción de la radiación solar que es reflejada por el suelo exterior al edificio. Esta variable, debe estar conectada a la que en el edificio aparece ilustrada como (GRDREF), que esta representa el mismo concepto. Una vez se han aclarado estos conceptos algo más tediosos, en la Tabla 2.2 se muestra un resumen de todas las variables de entradas al bloque que representa el edificio diseñado. Por su parte, en la Figura 2.8, se muestra el resultado de realizar todo el conexionado correspondiente. Debido a las limitaciones de espacio en la impresión de la pantalla las conexiones que no se muestran son las correspondientes a los ángulos de incidencia sobre las superficies. Tabla 2.2 Entradas del bloque del edificio (Type 56). Definición Notación Temperaturas exteriores TAMB, TSKY Humedad relativa del ambiente RELHUMAAMB Ángulos cenitales y acimutal AZEN AAZM Incidencias totales en cada orientaciones IT_N_180_90 a Beam radiation IB_N_180_90 a Ángulo de incidencia sobre cada superficie AI_N_180_90 a Reflectancia del suelo GRDRED aEsta variable se da para todas las orientaciones del edificio (norte, sur, este y oeste). Una vez se han realizado las conexiones correspondientes, se puede proceder a la simulación del archivo y comprobar como el edificio se ve afectado por su interacción con el ambiente. Como se expuso en el mapa conceptual del enlace entre los distintos bloques, el resultado de realizarlo en TRNSYS es el que se puede visualizar en la Figura 2.9. Antes de realizar las simulaciones, habría que manipular los parámetros de la simulación. Esto se realiza en la ventana resultante que se genera cuando se accede a Assembly/Control Cards, donde se deberá establecer un periodo de simulación, así como un tiempo de muestreo. Para el periodo de simulación, hay que saber que la base de datos del tiempo posee valores almacenados siendo el primer valor el del primer dia del año, por lo que para realizar una simulación en una época estival, se deberán introducir valores medios, situados estos alrededor de 5000 hr . Para ello, deben modificarse los campos Simulation start time y 14 Capítulo 2. Definición del sistema Figura 2.8 Conexiones realizadas entre el tiempo y el edificio. Simulation stop time. En relación al tiempo de muestreo, este puede establecerse en unas unidades distintas a las que se han especificado los límites de la simulación, en el caso de estudio un tiempo de muestreo de 1min es suficiente ya que los sistemas donde se analizan temperaturas de espacios son más lentos que muchos otros. Para realizar dicha modificación, se debe hacer en el campo Simulation time step. De aquí en adelante, se mostrarán los diversos comportamientos del sistema creado cuando este se ve influido por las condiciones atmosféricas que se dan en el exterior. Para ello, con el esquema realizado y habiendo definido los parámetros de simulación de forma adecuada y coherente, basta con ejecutarlo pulsando F8 . Si todo se ha realizado correctamente, se debe visualizar el comportamiento en una ventana auxiliar generada por TRNEXE. Figura 2.9 Resultado de realizar todas las conexiones. Para la presentación en esta memoria descriptiva, dichas gráficas resultantes, van a ser las realizadas con MATLAB ® , ya que las obtenidas con la herramienta de TRNSYS, no son de tan buena calidad como las obtenidas con la otra herramienta. Para ello, obviamente se debe realizar una comunicación entre ambos programas. Esta comunicación, en esta ocasión, se supone conocida, ya que será objeto de estudio en Capítulo 3. Estas simulaciones tendrán una duración de un día ( 24 hr ), por lo que se podrán observar todos los cambios producidos en cada una de las salas cuando se modifican los valores de las distintas incidencias de la radiación solar a lo largo del día. 2.5 Simulación en TRNSYS 15 Como se puede observar en la Figura 2.10 y Figura 2.11, durante las primeras horas del día, las radiaciones incidentes a primera horas de la mañana son mayores en la cara este del edificio, mientras que a última hora, este fenómeno se produce en la orientada hacia el oeste como es lógico. Sin embargo, en las caras norte y sur, las radiaciones son más duraderas en el tiempo pero presentan una menor incidencia. Figura 2.10 Radiación incidente en las caras Norte y Sur. Figura 2.11 Radiación incidente en las caras Este y Oeste. 16 Capítulo 2. Definición del sistema Con esta incidencia, las tempretauras de las habitaciones, representadas en la Figura 2.12 y Figura 2.13, resultan: Figura 2.12 Temperaturas de la sala 1 hasta la 4 considerando la influencia atmosférica. Figura 2.13 Temperaturas de la sala 5 hasta la 7 considerando la influencia atmosférica. Como se puede observar, las gráficas resultantes son coherentes, y habiendo una temperatura exterior que ronda los 32 ◦C en su pico máximo, las temperaturas de las salas ascienden a un valor comprendido entre los 27 ◦C y 28 ◦C . Además, aunque no se aprecie a simple vista, cuando la radicación incidente es máxima en la 2.5 Simulación en TRNSYS 17 cara oeste (última hora de la tarde), es en esa sala donde se produce un repunte en la temperatura alcanzando un valor de 27.57 ◦C, siendo este superior al del resto de los espacios en ese instante. Con esto se da por concluido la redacción del Capítulo 2, donde se ha explicado, el modelado del sistema y, se han expuesto los resultados obtenidos tras el modelado, observando que son coherentes, considerando así que el sistema es apto para proceder con los siguientes pasos. 3 Incorporación de Equipos Durante la redacción del correspondiente capítulo se pretende dar a conocer los distintos elementos que componen el sistema HVAC (Heating, Ventilation and Air Conditioning) que será empleado para acondicionar la temperatura de uno de los espacios que componen el edificio. Para ello, durante todo el capítulo se va a ir exponiendo uno a uno todos los elementos empleados y, una vez se disponga del esquema completo y el modelado de cada uno de ellos, se va a estudiar cómo debe implementarse en TRNSYS, usando para ello la herramienta MATLAB ® , por lo que también debe de explicarse cómo se comunican y que variables usan ambos programas. Por último, para poder comprobar si todo ello se ha realizado correctamente, se procede a la presentación de los resultados obtenidos y el correspondiente análisis de todos ellos. A modo de ejemplo, el objetivo de este capítulo puede expresarse en la Figura 3.1, donde todo lo que se encuentra bajo el rectángulo de color rojo, será lo que vaya a ser tratado. Figura 3.1 Objeto de estudio en el Capítulo 3. Para comenzar el desarrollo del capítulo, se partirá por la definición de cada uno de los componentes y cada uno de los modelos que se van a implementar para poder lograr la construcción del sistema que permita acondicionar las temperaturas de cada una de las salas [11] [1]. 3.1 Descripción del sistema y modelo de esquipos En primer lugar, antes de entrar en materia es conveniente saber qué es un sistema HVAC. Un sistema HVAC básicamente es un conjunto de equipos que deben cumplir tres funciones: la primera de ellas, hace referencia al calentamiento de una estancia para poder aumentar así su temperatura; el segundo objetivo es la ventilación de espacios, que consiste en la renovación de aire interior de un habitáculo mediante la aportación del aire exterior con el objetivo de asegurar su calidad; la última función es la de acondicionar la sala enfriándola para lograr así una disminución de su temperatura. Todos estos, son los pilares básicos sobre los que se apoyan los sistemas HVAC, pero durante la realización de este proyecto, el trabajo se centrará exclusivamente sobre uno de ellos, que se basará en la refrigeración. 19 20 Capítulo 3. Incorporación de Equipos Lo que se pretende es tener la posibilidad de mantener la temperatura de cada uno de los espacios, que se encuentran en un modelo diseñado, en unos niveles razonables cuando la temperatura exterior es superior a la deseada, es decir, cuando el edificio en cuestión se encuentra en una época estival. Dicho esto, ya se deja claro que el objetivo a conseguir es el enfriamiento mediante una extracción de energía, de las salas que componen un edificio. Para el desarrollo del resto del capítulo, el desarrollo explicativo se realizará para el caso particular que ha sido elaborado en el Capítulo 2, donde de forma resumida, se dispone de un edificio compuesto de seis salas (tres de ellas ubicación norte y las otras tres con orientación sur), todas ellas comunicadas entre sí mediante un pasillo. Volviendo al objeto del capítulo, como se ha dicho, el objetivo primordial es el enfriamiento de las habitaciones. Para ello se propone un sistema de refrigeración basado en un circuito cerrado de tuberías de agua. El funcionamiento de dicho circuito es el siguiente: el agua de la que se dispone en dicho circuito de hace pasar, mediante el impulso provocado por una bomba, por una enfriadora o chiller que permite disminuir su temperatura. Una vez se ha conseguido enfriar el agua, esta se distribuye uniformemente en un sistema de siete tuberías que desemboca en cada una de las salas. Para poder realizar el intercambio energético que permita la disminución de temperatura se hace uso de un fancoil o intercambiador de calor. La potencia intercambiada en cada espacio depende fundamentalmente del caudal de agua que circula por cada fancoil, dicho caudal se encuentra regulado por sendas válvulas de tres vías controlables que permite la modificación del caudal que circula por cada intercambiador. Por último, una vez el agua ha atravesado el intercambiador y su correspondiente temperatura ha aumentado, todos los caudales de agua son recogidos mediante un mezclador para lograr que el agua vuelva a ser enfriada por el chiller. Mediante este ciclo, se consigue extraer una determinada potencia a cada una de las salas y conseguir así disminuir la temperatura de cada una de ellas. A modo de resumen, lo anteriormente expuesto puede visualizarse en la Figura 3.2. En dicho esquema se ha simplificado el número de válvulas para poder ilustrar su composición en profundidad y poder detallar así donde se ubica cada uno de los componentes. Figura 3.2 Esquema del circuito refrigerador. A continuación, se va a ir exponiendo uno a uno cada uno de los elementos y las ecuaciones que modelan cada uno de ellos, dichas ecuaciones han sido elaboradas con datos de fabricantes comerciales para un determinado modelo [11]. En primer lugar, el elemento más básico que se encuentra es la válvula de tres vías, cuya función es básicamente regular la cantidad de caudal que se deriva al fancoil y por el canal de derivación. El parámetro que caracteriza su comportamiento se denomina apertura y puede tomar valores comprendidos entre cero y uno. El valor de dicho parámetro representa el porcentaje del caudal parcial (es el que se distribuye por cada módulo, es igual al que atraviesa cada intercambiador más el que es derivado) que pasa por el intercambiador, de forma que el resto es desviado por el canal de derivación manteniendo su temperatura constante. En las ecuaciones que describen el comportamiento, se denota como qps al que pasa por el intercambiador, qpd al 3.1 Descripción del sistema y modelo de esquipos 21 que es derivado y qp al caudal parcial que pasa por cada unidad. Con esto, las ecuaciones que modelan el comportamiento son qps =qpA qpd =qp(1−A)(3.1) donde, el parámetro A, representa la apertura de cada una de las válvulas. Hay que indicar que, durante todo el desarrollo que se va a realizar, las unidades en las que se expresan los caudales serán l/min . Por otro lado, se tiene que el caudal total que circula por esquema, representado por q, cumple q= 7 ∑ i=1 qpi (3.2) que, para el caso de estudio aplicado, se repartirá de forma equitativa por cada una de las salas en cuestión, por lo que el caudal parcial de cada uno de los módulos será una séptima parte del caudal total de partida que se ha establecido en el diseño del equipo. El siguiente elemento básico que se encuentra en el esquema es la mezcladora, cuya función es garantizar la mezcla de todos los conductos disponibles, los cuales incluyen los que han atravesado el fancoil tanto como los que han sido desviados. A la salida de la enfriadora el caudal es conocido ya, pues es la suma de todos los caudales parciales, denotado como q , sin embargo, el cálculo de la temperatura ( Twe ) es algo más tedioso, pues se tienen distintos caudales a distintas temperaturas. La temperatura resultante del agua se puede obtener realizando una media ponderada en base a los caudales, de cada una de las temperaturas que llegan a la mezcladora. De forma matemática, esto puede expresarse como Twe =∑7 i=1(Tpsqps +Tpeqpd ) q,(3.3) donde, Tps ( ◦C ) representa la salida de la temperatura del agua del intercambiador y Tpe ( ◦C ) es la temperatura del agua del canal de derivación, o lo que es lo mismo la temperatura del agua a la salida de la enfriadora, puesto que se supone que no se producen pérdidas energéticas en dicho canal. En relación a la enfriadora, esta se puede definir como una máquina frigorífica cuyo objetivo es el de enfriar un medio líquido. Lo que interesa con respecto a este equipo es, que para dada una capacidad nominal ( CAPN ) y un factor de carga parcial ( f cp ), se pueda determinar la temperatura ( Tws ) a la que es enfriado el caudal de agua que atraviesa dicho equipo. El caudal que entra está definido por la variable q ( l/s ) y por la temperatura Twe ( ◦C ) que, si no se consideran pérdidas en el conducto, esta coincide con la de la salida de la mezcladora. El factor de carga parcial que afecta a la enfriadora puede obtenerse mediante el cálculo que puede observarse en f cp =cpρq(Twe −Tws) CAP (3.4) donde, la cp es la capacidad calorífica del fluido, q(l/s) representa el caudal del trabajo de la máquina, Twe la temperatura a la entrada, Tws la de la salida, y la variable CAP, representa la capacidad de la enfriadora, la cual depende de su valor nominal ( CAPN ) y de los valores que toman la temperatura del aire exterior, Tae(◦C) , y de la temperatura del agua a la salida de la enfriadora, Tws . En forma de ecuación, todo ello puede expresarse como CAP =CAPN(c0+c1Tws +c2T2 ws +c3Tae +c4T2 ae +c5TwsTae),(3.5) en la cual pueden observarse algunos parámetros ci , los cuales han sido obtenidos por identificación por datos aportados por el fabricante. Los valores de estos coeficientes son los que pueden observarse en la Tabla 3.1. A modo de resumen, la enfriadora puede ser considerada como un bloque, al cual se le aporta un valor del factor de carga parcial ( f cp ) y una temperatura de entrada del agua ( Twe ), y este devuelve, en función de la temperatura ambiente (Tae), una temperatura de salida del agua (Tws). El último elemento que se va a emplear es el facncoil o intercambiador de calor. Este elemento está constituido por un conjunto de tubos por los que circula el fluido en cuestión, dicho fluido ve modificada su 28 Capítulo 3. Incorporación de Equipos Figura 3.5 Resultados obtenidos tras la modificación de las Aperturas (Salas 1 a 4). Figura 3.6 Resultados obtenidos tras la modificación de las Aperturas (Salas 5 a 7). Una vez se obtienen los resultados de las simulaciones que han sido realizadas, es conveniente destacar y obtener algunas conclusiones de los resultados obtenidos: • La primera conclusión que se puede obtener con una vista rápida sobre las gráficas, se encuentra en que se puede afirmar que, cuando la apertura de una válvula se modifica hacia un valor superior (se abre más la válvula), la temperatura de la sala o espacio disminuye, mientras que si esta se cierra se produce un aumento de la temperatura. Explicar este fenómeno de forma física no es complicado, 3.4 Resultados de las simulaciones 29 pues lo que realmente ocurre cuando se produce una modificación en alguna de las aperturas, es una variación del caudal que atraviesa el intercambiador, de forma que la potencia que absorbe el fluido será mayor cuanto mayor sea la apertura y viceversa. • Si uno se centra en lo que pasa en el caso concreto de la sala 5 (Figura 3.6), en el cual la apertura se mantiene constante, se puede observar que, a pesar de ello, la temperatura de la sala aumenta, aunque en menor medida que cuando se produce una modificación. Esto viene causado porque en la simulación se produce un cambio en algunas de las aperturas en el mismo instante de tiempo, modificando entonces el valor de la temperatura del agua que entra en el fancoil, ya que, aunque el fcp se mantiene constante (la capacidad de enfriamiento es la misma), la temperatura a la salida de los intercambiadores ya no es la misma, por lo que tampoco será la misma a la entrada de la enfriadora, generando así una temperatura distinta a la salida de este, causado todo ello por el principio de conservación de la energía, ya que si la máquina enfriadora genera la misma potencia frigorífica en todo momento, el que una sala se lleve más potencia implica que se la quita a la otras. • La última nota aclaratoria que sería necesario realizar, se puede centrar en las sobreoscilaciones que se producen en algunos casos cuando se aumenta la apertura de las válvulas, como ocurre en el caso de la sala 1 (Figura 3.5), donde se puede apreciar una pequeña sobreoscilación. Este fenómeno se produce debido a que en las simulaciones, TRNSYS autoajusta los valores mediante funciones de transferencia de primer orden cuando los parámetros no sufren muchas variaciones. Como se puede ver, las aperturas permanecen constantes durante casi toda la simulación excepto cuando se produce el cambio, por lo que tiende ajustar las curvas. Cuando se produce así el cambio, justo en el instante anterior, el sistema vuelve a ajustar, generando así una pequeña sobreoscilación que es corregida posteriormente. Para comprobar que ocurre cuando se modifica el valor del fcp, se realiza otra simulación en la que, manteniendo el valor de las aperturas constantes, se produce una única variación en el valor del fcp. Los resultados pueden observarse en las Figuras 3.5 y 3.6. Figura 3.7 Resultados obtenidos tras la modificación del fcp(Salas 1 a 4). En las simulaciones que se han realizado variando los valores que toma el fcp, se observa que, cuando dicho valor se aumenta, se puede ver que las temperaturas de las habitaciones disminuyen. Esto puede explicarse si se entiende cual es la función de la enfriadora, ya que cuando se produce un aumento en el fcp, la capacidad de la enfriadora aumenta, y por tanto la temperatura a la que sale el agua de la enfriadora disminuye, logrando 30 Capítulo 3. Incorporación de Equipos Figura 3.8 Resultados obtenidos tras la modificación del fcp (Salas 5 a 7). así, que la potencia absorbida por el agua tras su paso por el fancoil también aumente, generando por tanto una disminución en la temperatura de la sala. También se puede destacar que el enfriamiento producido en la sala 7, que es la que representa al pasillo, es inferior al del resto de salas, ya que el espacio que representa el pasillo es superior al de las salas, por lo que este necesita mayor potencia para lograr el mismo enfriamiento. 4 Identificación del sistema Un o de los pasos más importantes para poder lograr una correcta caracterización del sistema corresponde a la identificación del mismo, ya que los modelos térmicos del edificio son muy complejos por tener muchos factores a tener en cuenta y obtener ecuaciones físicas de ellos es muy difícil. Sin embargo, para las técnicas de control que van a ser empleadas es necesario disponer de estas ecuaciones para tener un conocimiento del sistema. Es por ello que se quiere realizar una identificación en espacio de estados que explique la dinámica del sistema [ 10 ], el cual puede definirse como un sistema de ecuaciones diferenciales que se apoyan en el uso de un vector de estados. El correspondiente modelo deberá presentar una forma como x(k+1) = Ax(k)+Bu(k)(4.1a) y(k) = Cx(k),(4.1b) donde x∈Rnx representa el vector de estados, u∈Rnu al conjunto de variables de control e y∈Rny representa las salidas del sistema. Para poder llevar a cabo dicha identificación, se va a hacer uso de un paquete de identificación de MATLAB ® . Para ello, se debe definir cuáles son las entradas y salidas del sistema en cuestión, preparar las pruebas necesarias para lograr la identificación, y por último, aplicar las herramientas necesarias para lograr el objetivo en cuestión. 4.1 Definición de entradas y salidas Antes de poder proceder con la identificación de un sistema cualquiera, es importante saber que variables del sistema van a ser consideradas como entradas y cuales como salidas, ya que esto define las dimensiones de cada uno de los vectores y matrices que caracterizan un espacio de estados. Para el caso de estudio de este proyecto, se parte de un sistema al cual se le quieren implementar distintos tipos de estrategias de control para poder obtener referencias deseadas en los valores de las temperaturas de cada una de las salas, siendo en este caso concreto un total de siete salas. Para poder aportar la energía capaz de producir un enfriamiento en dichas salas, se hace uso de un equipo frigorífico compuesto de una enfriadora, un intercambiador o fancoil en cada una de las salas, siete válvulas que modelan el flujo de fluido por cada uno de los intercambiadores, y una mezcladora que permite la unión de todos los flujos para poder impulsarlos de nuevo en el chiller o enfriadora. Cada uno de estos equipos poseen distintas variables que permiten modelar el comportamiento del mismo así como la potencia que se aporta a cada una de las salas. En este caso se toman como parámetros de control las válvulas que modelan el caudal que circula por cada fancoil, y el fcp que caracteriza a la enfriadora, el cual permite un mayor o menor enfriamiento del fluido que lo atraviesa. De esta forma, se puede definir un bloque llamado sistema que englobe tanto el modelo del edificio como el del conjunto de equipo que constituyen la máquina enfriadora. Dicho sistema, se puede modelar mediante las variables de entradas definidas, siendo estas un total de ocho, las cuales generan distintas salidas. Todo esto que acaba de ser enunciado, es lo que se puede extraer de la Figura 4.1. Con esto, ya se puede empezar a dimensionar algunos de los elementos que constituyen el modelo de espacio de estados. Más concretamente, el vector yestará compuesto de un total de 7 componentes que representan las temperaturas; el vector u es el relativo a las variables de control, y por tanto almacena las aperturas (de la componente uno a la siete) y el fcp (en la última componente) en un total de 8 componentes; 31 32 Capítulo 4. Identificación del sistema el vector de estados x todavía no se le puede determinar su dimensión exacta, pues tendrá que ser ajustado para obtener los mejores resultados posibles. Figura 4.1 Esquema entradas y salidas del sistema. 4.2 Procedimiento identificación Para poder realizar la identificación correspondiente al sistema, se debe hacer uso de un toolbox de identificación disponible en MATLAB ® . más concretamente el denominado como System Identification Toolbox. Antes de entrar a conocer debe usarse dicho paquete, hay que saber qué tipo de datos y la forma que deben adoptar estos para que la identificación que se produzca sea exitosa. El objetivo es poder aportar unos datos al toolbox en el que se pueda ver cómo afecta la modificación de cada una de las entradas sobre cada una de las salidas del sistema. Para ello, se pretende realizar tantas simulaciones como número total de entradas se dispongan. En cada una de esas simulaciones se debe ir generando un escalón sobre cada una de las entradas para comprobar así qué efecto produce esa entrada sobre cada una de las salidas. En este tipo de pruebas, las entradas que permanecen constantes lo hacen sobre un valor, siendo este el punto de funcionamiento. De esta forma la identificación obtenida es considerada como un sistema lineal en espacio de estados en torno a un punto de operación. Dicho esto, se puede comenzar a realizar cada una de las pruebas en cuestión, las cuales deben implementarse en TRNSYS mediante el esquema visible en la Figura 4.2. El modelado de cada una de las entradas en todas las pruebas se realiza mediante MATLAB ® , el cual genera unos valores de aperturas y fcp que son introducidos en el modelo del equipo enfriador y este genera cada una de las potencias cedidas en las habitaciones del edificio. Figura 4.2 Esquema obtención de datos en TRNSYS. Para poder realizar esta identificación se va a considerar una temperatura exterior del ambiente constante, de forma que cuando el tiempo atmosférico sea introducido, este sea considerado como una perturbación no conocida del sistema. Otros de los parámetros que deben ser modelados en el sistema son los puntos de operación elegidos para cada una de las salidas, los cuales han tomado un valor de 0.3 para las aperturas de las válvulas y 0.1 para el factor de carga parcial. El momento en el que se produce un cambio en las entradas del sistema será justo en la mitad de la simulación, momento en el que se produce un cambio en las entradas del sistemaes cuando el sistema se ha estabilizado correctamente antes de producir nuevas variaciones. Para ilustrar todo esto desde otra perspectiva, el código propuesto para modelar los distintos escalones y general variaciones sobre cada una de las salidas es: 4.2 Procedimiento identificación 33 %% Se diseñan las entradas para el modelado %% Se haran 8 pruebas, donde se modificaran las entradas del sistema Apertura = 0.3*ones(7,nTimeSteps+1); fcp = 0.1*ones(1,nTimeSteps+1); %% Establecemos 14 pruebas %% % DESCOMENTAR LO QUE PROCEDA para cada prueba % Apertura(1,:) = [0.3*ones(1,round((nTimeSteps+1)/2)) 0.6*ones(1,nTimeStepsround((nTimeSteps+1)/2-1))]; % Prueba 1 % Apertura(2,:) = [0.3*ones(1,round((nTimeSteps+1)/2)) 0.6*ones(1,nTimeStepsround((nTimeSteps+1)/2-1))]; % Prueba 2 % Igual hasta la prueba 8 Cuando se dispone de cada una de las entradas definidas, ya solo es necesario introducir el valor que toman cada una de ellas en un tiempo de muestreo en el modelo de la máquina que se ha implementado en el Capítulo 3 y almacenando cada una de las temperaturas de las salas del edificio, que son las salidas del sistema. Para poder almacenar los datos y disponer de ellos cuando se desee, se puede hacer uso del comando save de MATLAB ® , el cual permite almacenar algunas o todas las variables que están disponibles en el workspace, tal y como se muestra a continuación: %% Se almacenan los resultados de la simulación %% save(’prueba_1’,’time’,’Apertura’,’fcp’,’temperatura_habitaciones’); Esta línea de código debe ejecutarse cuando ya se ha acabado de realizar la simulación, ya que será cuando se disponga al completo de todos los datos obtenidos. Estas simulaciones se tienen que realizar tantas veces como número de entradas tenga el sistema, siendo para este caso concreto un número igual a ocho. De esta forma se deben disponer de un total de ocho resultados en los que se dispongan de todas las variables de entradas y sus correspondientes influencias en cada una de las salidas. Con estos resultados, ya se puede proceder a realizar la identificación de un espacio de estados del sistema completo. Figura 4.3 Entradas del sistema de la Prueba 2 . 34 Capítulo 4. Identificación del sistema Los datos de partida para realizar la identificación, como ya se ha dicho son un total de ocho pruebas realizadas, las cuales deben ser cargadas una a una en MATLAB ® . Los datos cargados, son los correspondientes a la simulación completa realizado con TRNSYS, por lo que en ellos hay algunos periodos de la simulación que no son de interés puesto que en los primeros tiempos se debe alcanzar un régimen permanente para los puntos de operación. Esto hecho genera la necesidad de recortar los valores de los vectores de entrada y salida del sistema. Por otro lado, para lograr una identificación en torno al punto de equilibrio, los valores de las entradas y salidas deben estar centrados en cero, de forma que habrá que restarle el punto de operación a las entradas y salidas del sistema. El punto de operación de las salidas puede ser considerado como el valor que toma la primera componente del vector de las temperaturas tras haber realizado la selección necesaria. A modo de ejemplo se van a mostrar dos de las pruebas tanto con las entradas y salidas tras haber realizado la selección y el centrado necesario. Figura 4.4 Salidas del sistema de la Prueba 2 . En las Figuras 4.3 y 4.4 se ilustran los resultados obtenidos tras la realización de la prueba 2, en la que el escalón de las entradas se produce en la apertura de la segunda válvula y se ve como el resto se mantienen completamente constante durante el resto del tiempo mostrado. Para comprobar que las modificaciones que se producen en las salidas son coherentes, se ve como la temperatura de la sala 2 baja algunos grados y el resto sufren un pequeño aumento, lo cual es lógico porque la temperatura a la que sale el agua del segundo intercambiador es algo superior a la de los instantes anteriores y como el fcp permanece constante, la capacidad de enfriamiento es las misma. Por lo que la potencia aportada a cada sala disminuye, produciendo así un aumento de temperatura. Por otro lado, en las Figuras 4.5 y 4.6 se muestran los resultados relativos a la última de las pruebas a realizar en la que se produce un escalón en el factor de carga parcial (fcp) de la enfriadora manteniendo todas las aperturas en el valor del punto de operación que se ha establecido. Con respecto al resultado obtenido, para poder entenderlo, hay que saber que cuando se produce un aumento en el valor del fcp lo que se produce es un aumento de la capacidad de enfriamiento del chiller, es decir, una disminución de la temperatura del agua cuando llega a la salida de éste. De esta forma aunque el caudal de agua que circula por cada intercambiador no varía (no se modifica ninguna de las aperturas), si lo hace la temperatura del fluido, logrando así que el intercambio energético que se produce entre la sala y el fancoil sea superior a los instantes anteriores, logrando así que la temperatura de cada una de las salas sufra un descenso siendo menor en el pasillo por las dimensiones de éste. 4.2 Procedimiento identificación 35 Figura 4.5 Entradas del sistema de la Prueba 8 . A la vista de los resultados obtenidos tras la selección correcta de los datos de interés y su correspondiente centrado en torno al punto de equilibrio, ya se está en disposición de hacer uso del toolbox de identificación disponible en MATLAB ® . Con la respectiva explicación de cómo hacer uso de él, se va a ir explicando cada una de las funciones de interés y como hacer un uso correcto de ellas. Dichas funciones deben ser aplicadas sobre los datos que se han obtenidos hasta ahora. Figura 4.6 Salidas del sistema de la Prueba 8 . 36 Capítulo 4. Identificación del sistema Para poder crear un modelo en espacio de estados, se deben encapsular todos los datos de las medidas del sistema que se quieren identificar, creando para ello un vector fila cuyas variables son de tipo IDdata, donde cada una de las componentes almacena las entradas y salidas de cada una de las pruebas realizadas, por lo que en caso en cuestión se creará uno de ocho componentes. Para ello basta con hacer uso de uno de los comandos de los que se disponen en el toolbox, más concretamente la función iddata, a la cual hay que pasarle un vector columna con todas las salidas obtenidas, otro vector columna con las diferentes entradas y, por último, hay que indicar que son datos discreto con un tiempo de muestreo igual a uno. De esta forma, la función debe presentar el esquema siguiente: IDdata{i}=iddata(Salidas,Entradas,1). De esta forma, se dispone de un vector formado por objetos de tipo iddata, los cuales deben ser integrados en un único objeto de tipo iddata, el cual albergue todos los datos de las pruebas realizadas. Para ello, se dispone de una función conocida con el nombre de merge, del inglés unir, la cual permite realizar dicha acción. Para ello solo es necesario pasar como argumentos cada una de las componentes del vector construido anteriormente IDmerged =merge(IDdat{1 : nexp}). Por último, ya solo falta obtener el espacio de estados resultante de la identificación, para lo cual el toolbox dispone de una función que hace dicha identificación. La función a la que se debe llamar es ssest que es capaz de estimar un espacio de estados tanto en discreto como continuo para un número de estados determinado a partir de los datos aportados. Uno de los parámetros de diseño es el número de estados a usar. Para poder estimar dicho número MATLAB ® dispone de un asistente que te asigna a priori la mejor opción llamando a la función ssmodel7=ssest(IDmerged,7 : 14,0Ts0,1), donde los argumentos de entrada son el conjunto de datos agrupados en un único objeto, los estados entre los que se desea buscar mejor opción y se debe indicar que el sistema a obtener se quiere en tiempo discreto con un tiempo de muestreo igual a uno. Si se prueba con unos valores que varíen entre uno y catorce estados, se obtienen los resultados de la Figura 4.7. En ella se puede observar que la mejor de las opciones propuesta es la de ocho estados, aunque esta habría que compararla con las versiones de siete y ocho, pues los resultados que se obtienen haciendo esto generan una primera idea, la cual no debe ser tomada como definitiva y se debe perfilar un poco más. Figura 4.7 Resultados obtenidos para decidir el número de estados a considerar. Tras realizar pruebas para los tres valores propuestos, se llega a la conclusión de que el mejor de los resultados que se obtienen es para un vector de estados formado por ocho componentes. A continuación, se van a analizar los resultados que se han obtenido. 4.3 Resultados obtenidos 37 4.3 Resultados obtenidos Como se acaba de decir ya se ha acordado que el espacio de estados que va a ser usado es el compuesto de un vector de ocho estados, pues este es el que mejor se comporta y modela el comportamiento real. Dicho modelo es un sistema de ecuaciones como x(k+1) = Ax(k)+Bu(k)(4.2a) y(k) = Cx(k),(4.2b) en el cual ya se pueden caracterizar cada una de las componentes: donde x∈R8 , u∈R8 e y∈R7 , de forma que las matrices A∈R8x8 , B∈R8x8 y C∈R7x8 , lo cual permite definir por completo el modelo matemático del sistema. Para poder afirmar que la identificación realizada se comporta de forma similar a como lo realizaría el sistema, se van a representar algunas comparaciones entre el sistema identificado y los datos reales. Además, una forma de observar el comportamiento del sistema en espacio de estados obtenidos se puede visualizar aplicando un escalón sobre cada una de las entradas y ver cómo afecta a cada una de las salidas. Para comenzar con el análisis en la Figura 4.8 se muestra la comparación entre los resultados obtenidos tras la identificación de la prueba 2, en la que el escalón se produce en la apertura de la segunda válvula. Con los resultados que se observan en la imagen, se puede afirmar que son buenos, ya que para todas las salidas el régimen permanente se ajusta por completo y en la cual se produce el escalón, el ajuste es prácticamente ideal. El transitorio de las otras salidas no se ajusta tan bien debido a que presenta un cero de fase no mínima, los cuales no pueden identificarse con los estados establecidos, de todas formas esto no causa ningun inconveniente ya que interesa la constante de tiempo y el régimen permanente (sin error). Figura 4.8 Resultados en la identificación de la Prueba 2. Por otro lado, en la Figura 4.9, se pueden ver los resultados obtenidos para la prueba 8 en comparación con los datos reales. En ellos se puede observar que el ajuste en régimen permanente es casi perfecto y durante el transitorio es ajuste que se produce casi que se ajusta al valor real. Con estos resultados, se puede afirmar que el modelo obtenido es bueno, ya que los resultados restantes que no se han mostrado en esta memoria por cuestión de espacio, son muy similares a los obtenidos en el caso de la prueba dos, ya que el comportamiento de las variaciones en las aperturas de las válvulas es muy similar en las habitaciones, así como su influencia sobre las demás salas. 44 Capítulo 5. Simulador de Perturbaciones for k = 1:numero_dias num_personas_sala1 = [num_personas_sala1 0*ones(1,llegada1(1)) num_personas _llegada1(1).*ones(1,(llegada2(1)-llegada1(1))) num_personas_llegada 2(1).*ones(1,(llegada3(1)-llegada2(1))) num_personas_llegada3(1).*ones (1,salida1(1)-llegada3(1)) num_personas_salida1(1).*ones(1,salida2(1)- salida1(1)) num_personas_salida2(1).*ones(1,salida3(1)-salida2(1)) num_ personas_salida3(1).*ones(1,hora_cierre-salida3(1)) 0*ones(1,dia-hora_ cierre)]; end Con ello, ya se puede dar por concluida la forma de implementar el número de personas. Se pasa ahora al estudio en profundidad de cómo establecer si las luces de cada una de las salas se encuentran encendidas o no. Como se expuso anteriormente, las luces de cada sala se encontrarán encendidas siempre y cuando haya una persona en el interior de la sala en cuestión. Para poder determinar esto, se debe hacer uso de los vectores anteriormente creados que indican el número de personas que llegan y salen en cada momento de cada una de las salas ( num_personas_llegadai y num_personas_salidai ). Cuando se disponga de dicha información, ya se puede construir el vector completo que indica cual es el estado de las luces en cada tiempo de muestreo. En código siguiente únicamente se detallará con la llegada 1 para reducir el texto introducido. for i = 1:7 if num_personas_llegada1(i)>0 luces_llegada1(i) = 1; end end luz_S1 = [luz_S1 0*ones(1,llegada1(1)) luces_llegada1(1).*ones(1,(llegada2(1)- llegada1(1))) luces_llegada2(1).*ones(1,(llegada3(1)-llegada2(1))) luces_ llegada3(1).*ones(1,salida1(1)-llegada3(1)) luces_salida1(1).*ones(1,salida 2(1)-salida1(1)) luces_salida2(1).*ones(1,salida3(1)-salida2(1)) luces_ salida3(1).*ones(1,hora_cierre-salida3(1)) 0*ones(1,dia-hora_cierre)]; Dicho esto, ya se tiene identificado por completo cual es el estado de cada una de las luces de las salas en cada momento. Para poder determinar el estado de los ordenadores en cada momento de la simulación, se parte de la hipótesis de que cuando una persona llega a la sala, esta puede encender o un ordenador de forma completamente aleatoria. Para ello, se generan vectores ( num_ordenadores_llegadai y num_ordenadores_salidai ) para cada una de las llegadas en los que se indica si se enciende o no un ordenador, siempre con la condición de que haya llegado una nueva persona a la sala. Adicionalmente, se han de saturar los valores respectivos, de forma que el número de ordenadores nunca pueda ser negativo en los distintos turnos de salida. Esto resulta: % Se crean los posibles encendidos de los ordenadores num_ordenadores_llegada1 = round((1-0).*rand(1,7)+0).*num_personas_llegada1; num_ordenadores_llegada2 = num_ordenadores_llegada1+round((1-0).*rand(1,7)+0) .*(num_personas_llegada2-num_personas_llegada1); % Se crean los posibles apagados de los ordenadores num_ordenadores_salida1 = num_ordenadores_llegada3 - (num_personas_llegada3 - num_personas_salida1); num_ordenadores_salida2 = num_ordenadores_salida1 - (num_personas_salida1 - num _personas_salida2); % Se saturan los valores para cada una de las salidas for i = 1:7 if num_ordenadores_salida1(i) < 0 num_ordenadores_salida1(i) = 0; end end Con ello, ya solo es necesario formar el vector completo del estado de los ordenadores en cada momento, para poder aportar al modelo del edificio los valores necesarios en cada tiempo de muestreo. 5.3 Implementación en MATLAB®45 ordenadores_sala1 = [ordenadores_sala1 0*ones(1,llegada1(1)) num_ordenadores_ llegada1(1).*ones(1,(llegada2(1)-llegada1(1))) num_ordenadores_llegada2(1) .*ones(1,(llegada3(1)-llegada2(1))) num_ordenadores_llegada3(1).*ones(1, salida1(1)-llegada3(1)) num_ordenadores_salida1(1).*ones(1,salida2(1)- salida1(1)) num_ordenadores_salida2(1).*ones(1,salida3(1)-salida2(1)) num_ ordenadores_salida3(1).*ones(1,hora_cierre-salida3(1)) 0*ones(1,dia-hora_ cierre)]; De esta forma, ya solo hace falta indicar cómo se va a realizar el modelado de las infiltraciones que se producen cuando abren o cierran las puertas. Para ello, tras la búsqueda de diversa información, se ha llegado a la conclusión de que la mejor opción para cuantificar el valor de dicho caudal de aire que pasa de una sala a otra se puede determinar mediante la fórmula de Pham and Oliver (1983) [ 4 ], en la cual permite el cálculo de la infiltración como I(kg/hr) = 3600 0.226 Aρi(gH)0.5ρi−ρo ρi0.52 1+(ρo/ρi)0.333 1.5 ,(5.1) donde I representa a las infiltraciones producidas, A es el área de la puerta por la que se produce dicha infiltración, g indica el valor de la gravedad en m/s2 , H hace referencia a la altura de la puerta, mientras que ρi y ρo son los valores de las densidades en cada uno de los espacios que se comunican. Como se puede observar, solo son necesarios los cálculos de cada una de las dos densidades. Como el cálculo de las infiltraciones depende de la densidad del aire de la sala en cada momento, dicho procedimiento se debe realizar junto con las determinaciones de la densidad en cada tiempo de muestreo. Para poder obtener las densidades, se recurre a la fórmula de la densidad revisada (CIPM-2007) [ 14 ] que se establece como ρa=pMA Z RT 1−xv1−Mv Ma,(5.2) siendo p(Pa) la presión a la que se encuentran las salas, T(K) la temperatura de cada una de las salas, Ma(g/mol) la masa molar del aire seco, Mv(g/mol) la masa molar del vapor de agua, R(J/molK) = 8.314475 y, por último, Z y xv que representan el factor de compresibilidad y la fracción molar de vapor de agua respectivamente, las cuales se calculan con xv=h(α+βp+γt2)exp(AT 2+BT +C+D/T) p(5.3a) Z=1−p T[a0+a1t+a2t2+(b0+b1t)xv+(c0+c1t)x2 v]+ p2 T2(d+ex2 v),(5.3b) donde h (tanto por uno) representa la humedad relativa de la sala en cuestión, T son las temperaturas en grafos K mientras que ten ◦Cy los valores de los coeficientes son los expresados en la Tabla 5.3 Tabla 5.3 Coeficientes utilizados para el cálculo de la densidad.. Coe f Valor Coe f Valor A1.2378847 ×10−5K−2a01.58123 ×10−6KPa−1 B−1.9121316 ×10−2K−1a1−2.9331 ×10−8Pa−1 C33.93711047 a21.1043 ×10−10K−1Pa−1 D−6.3431645 ×103K b05.707 ×10−6KPa−1 α1.00062 b1−2.051 ×10−8Pa−1 β3.14 ×10−8Pa−1c01.9898 ×10−4KPa−1 γ5.610 ×10−7K−2c1−2.376 ×10−6Pa−1 e−0.765 ×10−8K2Pa−2d1.83 ×10−11K2Pa−2 Con estas ecuaciones, ya se está en disposición de comenzar el cálculo en cada uno de los tiempos de muestreo, en los cuales se debe leer en primer lugar, cada una de las temperaturas de las salas y con ellas, 46 Capítulo 5. Simulador de Perturbaciones proceder al cálculo de la densidad. Este cálculo implementado, se puede expresar como: % Densidad para la Sala 1 Xv_S1 = (RH/100) * (alpha + beta*P_salas + gamma*T_S1^2) *(exp (A*Tk_S1^2 + B* Tk_S1 + C+ D/Tk_S1) / P_salas); Z = 1-(P_salas/Tk_S1)*(a0+a1*T_S1+a2*T_S1^2+(b0+b1*T_S1)*Xv_S1+(c0+c1*T_S1)*Xv_ S1^2)+(P_salas^2/Tk_S1^2)*(d+e*Xv_S1^2); densidad_S1 = (P_salas*Ma)/(Z*R*Tk_S1)*(1-Xv_S1*(1-Mv/Ma)); Hay que tener en cuenta que previamente se han tenido que inicializar todas las variables y coeficientes correspondientes para que el cálculo se pueda llevar a cabo. Para poder establecer cuando si y cuando no se va a producir una infiltración, se debe hacer uso de los que almacenan los datos de llegada y permiten, por tanto, saber cuándo se abre y cierra una puerta. Cuando este fenómeno se produce, el caudal de aire que pasa de un lado a otro depende de las densidades del aire en cada uno de los dos lados, de forma que se produce un flujo de aire de zonas de mayor densidad hacia zonas de menor valor. Para ello, en cada momento se debe comprobar cuál es el valor de la densidad en cada sala y establecer el sentido de flujo para poder otorgar el valor en el sentido correcto. Esto puede traducirse en código como lo siguiente: if ((todas_llegada1(1+7*w) == nStep && infiltraciones_llegada1(1+7*w) == 1) || (todas_llegada2(1+7*w) == nStep && infiltraciones_llegada2(1+7*w) == 1) || (todas_llegada3(1+7*w) == nStep && infiltraciones_llegada3(1+7*w) == 1) || (todas_salida1(1+7*w) == nStep && infiltraciones_salida1(1+7*w) == 1) || ( todas_salida2(1+7*w) == nStep && infiltraciones_salida2(1+7*w) == 1) || ( todas_salida3(1+7*w) == nStep && infiltraciones_salida3(1+7*w) == 1) || ( hora_cierre == nStep && infiltraciones_cierre(1+7*w) == 1)) if densidad_S1 < densidad_Pasillo infiltraciones_a_S1 = 0; infiltraciones_Pasillo_de_S1 = 3600*0.226*Area_puerta*sqrt(g*H)* densidad_Pasillo*sqrt((densidad_Pasillo-densidad_S1)/densidad_ Pasillo)*(2/(1+(densidad_S1/densidad_Pasillo)^0.333))^1.5 else infiltraciones_a_S1 = 3600*0.226*Area_puerta*sqrt(g*H)*densidad_S1*sqrt ((densidad_S1-densidad_Pasillo)/densidad_S1)*(2/(1+(densidad_Pasillo /densidad_S1)^0.333))^1.5 infiltraciones_Pasillo_de_S1 = 0; end end Hecho esto, ya se puede concluir con la explicación un poco más en detalle de cómo se ha elaborado el simulador, otorgando un valor a todas las variables de entrada en el bloque del edificio. Ya solo falta comprobar que todo funciona correctamente y establecer la conexión entre todos los bloques en un modelo de TRNSYS. 5.4 Resultados obtenidos Para poner a prueba todo lo que ha sido elaborado hasta ahora, se van a implementar todos los archivos que han sido elaborados hasta ahora en TRNSYS, para poder realizar las posteriores simulaciones. Para empezar en una plantilla de TRNSYS se deben realizar las distintas conexiones entre los bloques del tiempo atmosférico, el bloque del edificio y la llamada a MATLAB®, como se puede ver en la Figura 5.4 Las conexiones que se deben realizar entre todos los componentes son: • La conexión que se realiza entre le bloque del tiempo atmosférico y el edificio es la relativa a la incidencia de cada una de las variables exteriores sobre el edificio en cuestión, las cuales ya han sido tratadas de forma mucho más detallada con anterioridad. • La unión Simulador-Edificio, aporta al edificio todos los valores instantáneos de cada una de las variables que han sido detalladas con los distintos códigos de MATLAB ® y las cuales indican el 5.4 Resultados obtenidos 47 número de personas, ordenadores, estado de las luces y las infiltraciones producidas a través de las puertas. • El edificio únicamente debe aportar al simulador el valor que toma cada una de las temperaturas de las salas creadas en cada tiempo de simulación y que son necesarias para el cálculo de las densidades. Figura 5.4 Esquema del conexionado a realizar para introducir el simulador. Implementando este esquema de conexión ya se puede proceder a la simulación para observar cómo se comportan cada una de las salas. Dicha simulación corresponderá a un día completo con un tiempo de muestreo, y cuyo fin es observar cómo todos los parámetros se generan de forma aleatoria y cómo afectan a cada una de las temperaturas de cada sala. Para ello se va a estudiar que ocurre en tres de las salas durante la simulación completa. En la Figura 5.5 puede observarse qué ocurre en la Sala número dos durante un periodo igual a 2 días. en el primer día solo asiste una persona a la sala y permanecen en ella hasta la hora del cierre, además enciende el ordenador. Puede verse como la luz permanece encendida durante todo el periodo. Cuando dicha persona entra en la sala, se produce una infiltración hacia ella ya que la densidad de la sala es mayor a la del pasillo, mientras que cuando esta abandona la sala, no se produce ninguna infiltración hacia la sala y si se produce en el pasillo, ya que las densidades han cambiado. Se puede observar que cuando una persona llega a la sala y enciende el ordenador se produce un lógico aumento de la temperatura de dicha sala, la cual disminuye cuando ya no se encuentra nadie en su interior. En el segundo día ocurre algo similar, pero en este caso llega haber un total de 2 personas y dos ordenadores en la sala, por lo que se observa que se produce un mayor aumento de la temperatura durante ese periodo. Sin embargo, la infiltración provocada es menor, lo cual indica que las temperaturas de las salas y el pasillo en ese instante no son muy dispares. Figura 5.5 Situación de la Sala 2 durante la simulación. 48 Capítulo 5. Simulador de Perturbaciones Figura 5.6 Situación de la Sala 4 durante la simulación. En la Figura 5.6 se observa cómo evoluciona la temperatura de la sala 4 y como varían las distintas variables. Una de las cosas que hay que destacar durante el primer día es cómo afecta la llegada de 3 personas a la temperatura, se puede observar como aumenta drásticamente y como aumenta también la pendiente con la que crece la temperatura. También se observa cómo aumenta en mayor medida cuando se produce la llegada de una persona que enciende el ordenador. Para comprobar el correcto funcionamiento se puede observar el comportamiento de las luces, que permanecen encendidas siempre que haya una persona. Figura 5.7 Situación del Pasillo durante la simulación. 5.4 Resultados obtenidos 49 Por último, en la Figura 5.7 se puede observar cómo el tamaño de sala muestra gran influencia sobre las evoluciones de la temperaturas, ya que el pasillo casi triplica al tamaño de una sala y las temperaturas alcanzadas son algunos grados inferiores a la del resto de las salas. Esto puede explicar que cuando se produce la salida de las personas, las infiltraciones se produzcan hacia el pasillo, puesto que posee una menor temperatura, mientras que a primera hora de la mañana cuando se producen las primeras llegadas, es el pasillo el que tiene mayores temperaturas, debido que presenta un menor contacto con el exterior, y por ello las infiltraciones se producen en sentido contrario. Como se puede observar en las correspondientes simulaciones, se puede afirmar el correcto funcionamiento del simulador y el comportamiento completamente aleatorio que este presenta, de forma que se puede aplicar a cualquier simulación que se realice independientemente de la duración de ésta. 6 Control LQR En este capítulo, se pretende aplicar al sistema diseñado un controlador LQR para el seguimiento de referencias[ 5 ]. Un controlador LQR puede ser definido como un controlador en bucle cerrado para sistemas MIMO (Multiples Inputs Multiples Outputs). Para su diseño se parte del sistema identificado en el Capítulo 4. Dicho sistema era enunciado como x(k+1) = Ax(k)+Bu(k)(6.1a) y(k) = Cx(k),(6.1b) donde, como ya se dijo en su momento, x∈Rnx representa el vector de estados, u∈Rnu al conjunto de variables de actuación presentes en el sistema, por último, y∈Rny que hace referencia a las temperaturas de las habitaciones, las cuales se desean controlar. Durante este capítulo se llevará a cabo una explicación teórica del controlador LQR, así como de la versión del controlador LQR con efecto integral. Para ello, el primer lugar, se debe introducir un observador de estados, que, a partir de la salida del sistema y la actuación del controlador, este sea capaz de calcular el estado del sistema. Al final de capítulo, se exponen los resultados obtenidos de la aplicación de dicho controlador al sistema elaborado en capítulos anteriores. 6.1 Observador de estados Un observador de estados [ 7 ] se puede definir como un bloque que estima las variables de estado de un sistema, usando como base las variables de control y salida de un sistema. Para su explicación se parte del sistema expuesto en la Ecuación 6.1 que representa el sistema en cuestión. Se sabe que el modelo matemático del observador es igual que el modelo de la planta al que se le debe añadir un término de error que representa la diferencia entre la salida medida y la observada. Se debe saber que ˆx hace referencia durante todo el capítulo al vector de estados observado. Representando lo expuesto, se tiene ˆx(k+1) = Aˆx(k)+Bu(k) +L(y−Cˆx(k)) = (A−LC)ˆx(k)+Bu(k)+Ly(k).(6.2) En la Ecuación 6.2, el segundo y tercer término pueden agruparse de forma que construyan un vector u como composición de las variables de control y las salidas del sistema, tal que u=u(k) y(k).(6.3) La estructura que se ha implementado puede visualizarse gráficamente en la Figura 6.1. Realizando la composición del vector u , la ecuación del observador puede reescribirse de una forma más compacta como ˆx(k+1) = (A−LC)ˆx(k)+B Lu.(6.4) 51 52 Capítulo 6. Control LQR Figura 6.1 Esquema del observador utilizado. Con esto, ya se dispone de la ecuación del observador de estados (Ecuación 6.4), donde si se ha determinado correctamente el valor de L, se puede obtener una estimación del estado actual del sistema. La determinación de dicha ganancia en este caso se obtiene mediante la ayuda de MATLAB ® , haciendo uso del comando dlqr cuyas entradas son dos matrices que simulan el comportamiento del sistema, además de dos matrices que ponderas la intensidad con la que se modelan las variables. En este caso concreto, para el cálculo de dicha ganancia L, el comando a ejecutar sería L=dlqr(AT,CT,W,V)T .(6.5) De esta forma, se concluye con la explicación del primer paso necesario para la elaboración del esquema de control. La existencia del observador es algo fundamental debido a como se ha definido el sistema, ya que sin él no sería posible la estimación del vector de estados que lo caracteriza en cada momento. El siguiente paso, como es lógico es el cálculo de la variable de control. 6.2 Control LQR El controlador LQR del inglés Linear Quadratic Regulator, es un algoritmo de control óptimo, a partir del cual se obtiene una acción de control para lograr resolver el problema de control de un sistema, que en este caso será de tiempo discreto [ 6 ]. Para enunciar el problema de resolución del controlador, se parte de las ecuaciones del sistema en cuestión, representado como x(k+1) = Ax(k)+Bu(k)(6.6a) y(k) = Cx(k).(6.6b) El objetivo del control LQR, es el cálculo de una acción de control proporcional al estado del sistema. En forma de ecuación, esto puede verse como u(k) = −Kx(k),(6.7) donde, esa matriz de ganancias (K), se obtiene de la minimización de una función de coste, J, sujeta a la dinámica del sistema en cuestión. Esta función J, es elaborada ponderando el peso de los valores de tanto los estados como de las señales de actuación, todos ellos al cuadrado, de la siguiente manera J= ∞ ∑ k=0xT kQxk+uT kRuk,(6.8) 6.3 Control LQR con efecto integral 53 siendo las matrices Q >0 y R >0 sendos pesos que permiten variaciones en la forma de controlar. Para entender esto mejor se exponen dos ejemplos: • Si se toman unos valores de Q=1 y R=100 , se produce una acción de control mucho más suave, ya que se pondera mucho más el término de las variables de control, por lo que para que se minimice correctamente, ese término debe ser lo más pequeño posible. • Por el contrario, con unos valores de Q=100 y R=1 , la acción de control debe llevar el vector de estados a 0, por lo que como no se pondera tanto la acción de control, se producirá una acción más agresiva sobre las variables a controlar. Observando la importancia que toman ambos parámetros, cuando se proceda al diseño del controlador se debe buscar un compromiso entre la brusquedad en el controlador, que puede provocar sobreoscilaciones, y la velocidad de control, de forma que el sistema alcance el punto deseado en un tiempo razonable. De la Ecuación 6.8, se demuestra que la constante proporcional que minimiza la función de coste, cumple la siguiente ecuación: K= (R+BTSB)−1(BTSA).(6.9) Además se identifica la matriz S como solución de la ecuación de Riccati en tiempo discreto, la cual puede enunciarse como ATSA −S−(ATSB)(BTSB +R)−1(BTSA)+Q=0.(6.10) Por lo que, una vez se dispone de dicha matriz S, habiendo solucionado la ecuación de Riccati, es trivial el cálculo de la acción de control que se debe aplicar al sistema. Sin embargo, este tipo de controlador no contempla el problema de seguimiento o tracking, ya que de lo que se encarga es de la minimización de los estados del sistema. Por ello, se recurre al control LQR con efecto integral, el cual permite abordar el seguimiento de una referencia. 6.3 Control LQR con efecto integral Para abordar el problema del seguimiento de referencias con control LQR, lo que se conoce como control con efecto integral [ 2 ], se comienza planteando las referencias como si fuesen perturbaciones que afectan al sistema con el objetivo de eliminar el error, definido como la diferencia entre la referencia y la salida del sistema, cuando el sistema en cuestión se encuentra en régimen permanente. Además, se define el sistema en términos incrementales haciendo las siguientes declaraciones: ∆x(k) = x(k)−x(k−1) ∆u(k) = u(k)−u(k−1) ∆y(k) = y(k)−y(k−1).(6.11) Haciendo esta definición, el error puede ser definido como ∆e(k) = ∆r(k)−∆y(k) , Si además se trata la referencia o perturbación como un escalón donde en régimen permanente se cumple que r(k+1)−r(k) = 0 , resulta que el incremento del término del error puede expresarse como ∆e(k) = ∆r(k)−∆y(k) = −∆y(k).(6.12) Por otro lado, el sistema identificado del que se parte, por linealidad, puede expresarse también en términos incrementales, de forma que se tiene ∆x(k+1) = A∆x(k)+B∆u(k) ∆y(k) = C∆x(k).(6.13) Combinando las ecuaciones 6.12 y 6.13, el término del error puede escribirse en una como e(k) = e(k−1)−C∆x(k) , lo que permite reconstruir la primera ecuación del sistema, en una análoga pero ampliada 7 Control MPC En este capítulo, se abordará el diseño y sintonización de un controlador MPC (Model Predictive Control), que representa un controlador predictivo, para el sistema desarrollado en el Capítulo 4, el cual podía ser definido como x(k+1) = Ax(k)+Bu(k)(7.1a) y(k) = Cx(k),(7.1b) donde, como ya se dijo en su momento, x∈Rnx representa el vector de estados, u∈Rnu al conjunto de variables de actuación presentes en el sistema, por último, y∈Rny que hace referencia a las temperaturas de las habitaciones, las cuales se desean controlar. Para ello, en el capítulo presente, se enunciará en que consiste el control MPC [ 15 ] y el necesario uso de un observador de estados para obtener el estado actual del sistema en cada tiempo de muestreo, posteriormente se aplicarán dichos conocimientos al caso concreto de estudio y poder diseñar así una ecuación de optimización y, por último se analizarán los resultados obtenidos. 7.1 Observador de estados Al igual que ocurre con el controlador LQR, para el controlador MPC, es necesaria la existencia de un observador de estados [ 7 ]. Para el diseño de dicho observador se puede tomar como referencia la Sección 6.1, aunque se enuncie brevemente la estructura de dicho observador. El observador de estados permite calcular, como ya se ha enunciado en otras ocasiones, el estado actual en el que se encuentra el sistema, tomando como datos de partida el valor de las señales de salida del sistema y los valores de las variables de control, se calculan los estados en los que se encuentra el sistema en cada momento. Para dicho observador se emplea una estructura ampliada que permite realizar el cálculo de las perturbaciones existenetes en el sistema que serán de utilidad para la función de coste del controlador. Para ello, usa una estructura del tipo ˆx(k+1) = ˆx(k+1) d(k+1)=A0 0I−LC Iˆx(k)+B 0Lu,(7.2) donde el vector u es una composición de los vectores que agrupan variables de control y las salidas del sistema, es decir, presenta una estructura del tipo u=u(k) y(k).(7.3) Al igual, el cálculo de la matriz L, se puede realizar mediante la resolución del problema LQR que representa el sistema y, en un caso más práctico con MATLAB ® , con ayuda del comando dlqr (usado para sistemas en tiempo discreto) y del paso de los correctos argumentos, se pueden obtener sin ningún tipo de problema los valores correspondientes. 61 62 Capítulo 7. Control MPC 7.2 Control MPC El control MPC, del inglés Model Predictive Control, es una estrategia de control que presenta sus bases en el control óptimo, el cual se basa en tres puntos claves: el uso de un modelo para predecir la salida del sistema en una serie de tiempos de muestreo futuros; el cálculo de la mejor serie de secuencias futuras que minimicen una función de coste correspondiente; la aplicación en cada tiempo de muestreo de la primera señal de control, obviando el resto que han sido calculadas. Entre las ventajas que presentan este tipo de esquemas de control [ 3 ], se destaca que es aplicable a casi cualquier tipo de sistema, trata de forma sencilla los casos de sistemas multivariables, así como la implementación intrínseca de restricciones a tener en cuenta en la fase de minimización. Sin embargo, también presenta algunos inconvenientes como es la necesidad de requerir un modelo muy preciso para la predicción de la salida del sistema además de la carga computacional que conlleva dicho algoritmo. Como se ha dicho, la primera fase del control MPC consiste en la realización de una predicción de la salida y estado del sistema, cuyo objetivo es que con datos anteriores de los que se disponga, obtener y determinar el estado actual en el que se encuentra el sistema. Para ello, se debe disponer de un modelo que proporcione información del tipo dx dt =f(x,u,t)(7.4a) y=h(x,u,t)(7.4b) x(to) = xo,(7.4c) donde x∈Rnx representa los estados, u∈Rnu son las variables de control, y∈Rny las salidas del sistema. Además, se indica una condición inicial del estado del sistema, la cual en el instante cero suele definirse como nula. Como el caso de estudio es un sistema discreto, el modelo de predicción que se usa es el disponible para tiempos discretos o también conocido como modelo de espacio de estados, dicho modelo ya ha sido expuesto en numerosas ocasiones, pues también ha sido objeto de estudio en capítulos anteriores. El modelo es el siguiente x(k+1) = Ax(k)+Bu(k)(7.5a) y(k) = Cx(k),(7.5b) Con dichas ecuaciones, se dispone ya de un modelo de predicción bastante correcto para la realización del controlador. Como se dijo, la segunda fase del problema consiste en la minimización de una función de coste con el objetivo de obtener unas secuencias de actuaciones optimas. Para su diseño se debe analizar el problema en cuestión y elegir una que acorde con las características de cada uno. Para ello, se suele partir de una básica que satisfaga un Linear Quadratic Problem el cual permita la minimización de una función de coste básica a la que se le vayan añadiendo diversos términos. Esta función de coste básica puede verse como la medición de la desviación entre la diferencia entre los estados y variables de control con un valor cero. Por ello, dicha función puede representarse en su forma más básica como V(x(0),u) = N−1 ∑ k=0 [x(k)TQx(k)+ u(k)TRu(k)]+ x(N)TPx(N) s.a.=x(k+1) = Ax(k)+Bu(k),(7.6) donde N representa el número de tiempos de muestreos futuros a considerar, conocido como horizonte de predicción. El objetivo con dicha función es encontrar una secuencia de actuaciones futuras que minimicen lo máximo posible su valor. Estas secuencias de actuaciones, están representadas con u={u(0),u(1),u(2),...,u(N−1)}.(7.7) Como puede observarse, además del diseño de la función objetivo, se deben sintonizar los parámetros Q y R, donde en función de ambos se puede conseguir una respuesta de un tipo u otro. Cuando el valor de Q es bastante alto con respecto al de R, lo que se persigue es el conseguir una rápida minimización de los estados, es decir, que el sistema lleve los estados al punto de equilibrio lo antes posible. Si por el contrario, el valor de 7.3 Aplicación al caso de estudio 63 R es superior al de Q, se consiguen que, penalizando la rapidez del sistema, los valores de las actuaciones sean lo más pequeños posibles. Resumiendo esto, el problema de minimización puede enunciarse como m´ ın uV(x(0),u).(7.8) La última fase es la correspondiente a la asignación de las primeras variables aplicables al primer tiempo de muestreo obtenidas durante la realización del problema, ya que de todas las obtenidas en cada tiempo de muestreo, obviando por tanto las variables obtenidas para tiempos de muestreo posteriores. De esta forma, el valor a aplicar es u∗=u(1).(7.9) Con esto, se puede concluir una breve explicación acerca del control MPC. Además de lo expuesto anteriormente, hay que considerar que se a la función de minimización se le pueden añadir diversas restricciones, como pueden ser: límites en las variables de control, límites en las salidas del sistema y variaciones máximas y mínimas de las actuaciones entre otras. En el siguiente apartado, todo esto se estudiará un poco más en detalle para el caso concreto de estudio, que es el aplicado al edificio. 7.3 Aplicación al caso de estudio En este caso, los conceptos expuestos anteriormente, se llevarán a la aplicación del sistema concreto del que se dispone, el edificio compuesto de las siete salas. Una de las ventajas que se observan en este controlador es la posibilidad de incorporar restricciones, lo cual cambia un poco la caracterización del sistema utilizado hasta ahora. Se sabe que los modelos de las máquinas obtenidos presentan su rango de operación en un rango de la temperatura del agua que atraviesa la enfriadora entre 5◦Cy15◦C. Por esta razón y dadas las circunstancias de que al controlador MPC se puedan añadir restricciones, el nuevo sistema identificado tendrá una nueva salida a controlar que será la temperatura del agua en la salida de la enfriadora (Twsnext ). De esta forma, la estructura del nuevo sistema identificado es x(k+1) = Ax(k)+Bu(k)(7.10a) y(k) = Cx(k),(7.10b) donde x∈R8 representa el vector de estados, u∈R8 al conjunto de variables de control, por último, y∈R8a las temperaturas de las habitaciones y la temperatura del agua a la salida de la enfriadora. Por otro lado, la función de coste [ 11 ] [ 9 ] está compuesta por varios términos, donde todos ellos, representan un fin a conseguir en el cálculo de la secuencia óptima de acciones de control. En dicha función de coste se tienen que establecer una serie de restricciones entre las que caben destacar las siguientes: • La primera restricción a tener en cuenta es la dinámica que sigue el sistema, es decir el modelo de espacio de estados que ha sido elaborado. Con esta dinámica se calculan algunas variables que son necesaria en la ecuación de predicción. • Una de las restricciones se encuentran en los valores que pueden tomar las variables de control. En relación a las aperturas de las válvulas, estas pueden variar su apertura entre cero y uno, por lo que los valores calculados deben encontrarse entre ambos valores. Por otro lado, y de forma análoga ocurre con el factor de carga parcial, el cual debe estar comprendido en el mismo rango de valores que las aperturas. • De forma complementaria se establecen restricciones en las variables de salida. Las temperaturas de las salas están representadas por las primeras variables, mientras que la última de ellas representa la temperatura del agua de salida en la enfriadora, en la cual es realmente importante asegurar su cumplimiento ya que los equipos están diseñados para un rango de valores determinados. • Por último, en un estado estacionario ficticio, la dinámica del sistema debe seguir cumpliendose, sabiendo ademas que los valores de las ecuaciones se mantienen constante y no varían. Esto puede expresarse como que el estado siguiente debe ser igual que el anterior. 64 Capítulo 7. Control MPC La función de coste propuesta, junto con el conjunto de restricciones que debe satisfacer, es m´ ın x,u,xs,us N−1 ∑ k=0x(k)−xs u(k)−usTQ0 0Rx(k)−xs u(k)−us+[x(N)−xs]TP[x(N)−xs] +[ys +d−yre f ]TS[ys +d−yre f ]+[us−uop]TJ[us−uop](7.11a) s.a x(k+1) = Ax(k)+ Bu(k)(7.11b) y(k) = Cx(k)(7.11c) xs=Axs+Bus(7.11d) ys=Cxs(7.11e) x(0) = ˆx0(7.11f) u∈U(7.11g) us∈U(7.11h) y∈Y(7.11i) ys∈Y(7.11j) donde, us∈R8 representa a las variables de control en el estado estacionario, ys∈R8 a las variables de salida medidas en el estado estacionario, xs∈R8 son los estados estacionarios, Y representa el conjunto de restricciones aplicables a las salidas del sistema, U hacen referencia a las restricciones en las variables de control, x0∈R8 es el estado inicial del que parte el sistema en cada tiempo de muestreo, al igual que ocurre con los valores de la variable d∈R8 que ambos son calculados mediante el observador y, por último, yre f ∈R8 son el conjunto de referencias establecidas para cada una de las salidas. A continuación, se exponen de una forma más detallada cada uno de los términos que componen la función a minimizar: • El primero de ellos, representa el seguimiento de los estados que se deben lograr durante el cálculo, es decir, la diferencia entre el estado en cada iteración (x(k)) y el estado estacionario al que se debe llegar. Este término es acumulativo, o en otras palabras, debe calcularse y tantos tiempos de muestreo como valga el parámetro del horizonte de predicción (N). Por otro lado, la ponderación de este término está marcada por el valor de la matriz diagonal Q. • Otro de los términos que se debe calcular y sumar N veces, representa la diferencia entre el valor de la actuación en cada instante y su valor estacionario. Dicho término cobra mas o menos protagonismo en base al valor de la matriz diagonal R. • En relación a los términos que no deben acumularse N veces, el primero que encontramos, pondera el valor el estado justo en el tiempo de muestreo después del valor de N. Este término, se conoce como coste terminal y se encuentra ponderado por el valor de la matriz P, la cual es el resultado de la solución de la ecuación de Riccati en tiempo discreto. El coste terminal surge del control óptimo LQR y garantiza la estabilidad del controlador MPC. ATPA −P−(ATPB)(BTPB +R)−1(BTPA)+Q=0.(7.12) • Uno de los términos que asegura el cumplimiento de la referencia en el control, está ponderado por la matriz S, donde se debe prestar especial atención, puesto que ahora el sistema de estudio está compuesto de una nueva salida, la temperatura del agua, la cual no debe seguir fielmente una referencia si no que debe mantenerse dentro de unos límites establecidos. Para ello el peso de la última componente de la diagonal de la matriz S (conocida como peso de seguimiento), debe ser bastante inferior al de las otras para aportar un poco de holgura al no cumplimiento de la referencia. Volviendo al término en cuestión, este básicamente evalúa y pondera la diferencia en valor absoluto entre el valor de la referencia y el de la salida en estado estacionario, incluyendo en este las posibles perturbaciones. • Por último se incluye también un término económico en la función, el cual permite minimizar el factor de carga parcial lo máximo posible, tomando así valores más altos para las aperturas de las válvulas con respecto al fcp. Gracias a esto, se garantizará un menor consumo de energía eléctrica y por tanto un menor coste económico. Este término ni mucho menos tiene por objetivo el aportar una cifra o valor a dicho coste. 7.3 Aplicación al caso de estudio 65 A la hora de establecer la sintonización de cada uno de los valores de los parámetros hay que llegar a un compromiso entre rapidez, sobreoscilación y el cumplimiento de la referencia. Para ello, hay que hacer una modificación adecuada en los valores de las matrices de forma que cuando se aumentan los valores de cada una de ellas se obligue a realizar una mayor minimización en dicho término. En relación a la implementación de dicho bucle de control en MATLAB ® , el proceso puede dividirse en varios pasos: en primer lugar, se deben leer cada una de las salidas del sistema, una vez se hayan leído, con ayuda del observador y el controlador, se calcula la señal de actuación y, por último, se aplican dichos valores calculados al sistema en cuestión. Antes de realizar este bucle , se deben realizar una serie de inicializaciones que permiten poder realizar el proceso de forma exitosa. En primer lugar, para lograr la determinación de cada una de las matrices que componen la función de minimización se deben realizar dos acciones. Por un lado, se debe usar la función dlqr, con la cual se resuelve la ecuación de Riccati y aporta el valor de la matriz P. El resto de matrices pueden definirse manualmente sabiendo que son matrices diagonales. De forma que el código es: % Calculo de L observador W = 0.01.*eye(nx++ny); V = 0.001.*eye(ny); L = dlqr(A_gorro’,C_gorro’,W,V); L = L’; % Calc K, P Q = 1.*eye(nx); R = 0.1.*blkdiag(1,1,1,1,1,1,1,10); [K,P]=dlqr(A,B,Q,R); K=-K; % Definicion de S y J S = 100000.*blkdiag(1,1,1,1,1,1,1,0.01); J = 0.01.*blkdiag(0.1,0.1,0.1,0.1,0.1,0.1,0.1,1); % Termino para las restricciones blandas epsilon = 1000; Como se puede observar en el código, aparece un término epsilon (ε) , el cual es usado para la incorporación de una restricción blanda en las variables de salida medida, pues se producen algunos problemas cuando se introduce restricción en estas variables de forma que se facilita el proceso de iteración y convergencia del algoritmo de optimización introduciendo una pequeña holgura para lograr realizar la simulación. Este término, tiene un peso muy grande para evitar que la holgura calculada sea muy grande y por tanto las restricciones no se cumplan fielmente. Para la creación del controlador MPC, se debe recurrir a la librería Yalmip [ 19 ]. Para ello, primero se deben establecer los puntos de operación sobre los que se han linealizado el sistema de espacio de estados, ya que para la introducción de restricciones estos puntos hay que tenerlos en cuenta. Para la declaración de restricciones se una estructura del tipo params en la cual se deben almacenar los valores máximos y mínimos que pueden tomar las variables, si estos valores se introducen en términos absolutos (números físicos reales), es cuando hay que restarle los puntos de operación sobre los que se ha linealizado el sistema. Cuando ya se dispone de la estructura, se debe crear unas variables especiales con la función Polyhedron, las cuales serán usadas en la definición del controlador. En otras palabras, se tiene que: % Puntos de operacion params.Uop = [0.3.*ones(nu-1,1);0.1]; params.Yop = [21.4408;21.2162;21.4408;21.4304;21.2619;21.4980;23.1351;12.8903]; % Restricciones establecidas params.Ylb = [-50.*ones(ny-1,1);5]-params.Yop; params.Yub = [50.*ones(ny-1,1);15]-params.Yop; params.Ulb = 0.*ones(nu,1)-params.Uop; params.Uub = 1.*ones(nu,1)-params.Uop; 66 Capítulo 7. Control MPC % Creación de los polihedros FU = [diag(1./params.Uub);diag(1./params.Ulb)]; PU = Polyhedron(FU,ones(2*nu,1)); FY = [diag(1./params.Yub);diag(1./params.Ylb)]; PY = Polyhedron(FY,ones(2*ny,1)); Para la creación del objeto del controlador, se deben diferenciar principalmente dos tipos de variables, las llamadas constr, son las que representan a las restricciones, mientras las correspondientes a cost son las que se aprecian en la función a minimizar. Mediante esto, la definición puede realizarse como: % Restricciones iniciales y estacionarias constr = constr + [x{1} == x0]; constr = constr + [xs == A*xs + B*us]; constr = constr + [ys == C*xs]; for k = 1:N % Funciones costes acumulativas cost = cost + [x{k}-xs]’ * Q * [x{k}-xs]; cost = cost + [u{k}-us]’ * R * [u{k}-us]; cost = cost + [epsilon * sum(s{k}).^2]; % Restricciones dinamicas constr = constr + [x{k+1} == A*x{k} + B*u{k}]; constr = constr + [y{k} == C*x{k}]; % Restricciones constr = constr + [PY.A * y{k} <= PY.b + s{k}]; constr = constr + [PU.A * u{k} <= PU.b]; constr = constr + [PY.A * ys <= PY.b + s{k}]; constr = constr + [PU.A * us <= PU.b]; end % Resto de costes cost = cost + [(x{N+1}-xs)’*P*(x{N+1}-xs)]; cost = cost + [(ys-ref+d)’*S*(ys-ref+d)]; cost = cost + [(us-operacion)’*J*(us-operacion)]; % Creación del objeto MPC options = sdpsettings(’solver’,’quadprog’,’verbose’,1); ctrl = optimizer(constr,cost,options,vertcat(x0,d,ref),vertcat(u{1},cost)); Como se puede ver, este código lo que crea es, básicamente un bloque llamado ctrl, al cual se le aportan los valores del estado calculado por el observador, la referencia y el conjunto de estados calculados por el observador referentes a las salidas del sistema (sistema ampliado), el cual está representado por la variable d. Lo expuesto aquí, puede verse representado en la Figura 7.1. Una vez se han realizado estos cálculos y definiciones previas, se puede proceder al primer paso que se debe realizar en el bucle de control, la lectura de las salidas del sistema, es decir, las temperaturas de las salas y la temperatura del agua a la salida de la enfriadora. Para ello se hace uso de la variable llamada trnInputs, que es usada para la comunicación entre MATLAB ® y TRNSYS. Se debe prestar atención y calcular dichas variables en términos relativos al punto de operación, pues son los puntos donde han sido linealizado el sistema. Dicho esto, se tiene el siguiente fragmento de código. % Lectura de todas las salidas del sistema Tae = trnInputs(1); y_actual(1,nStep) = trnInputs(2); 7.3 Aplicación al caso de estudio 67 Figura 7.1 Bloque del objeto del controlador (ctrl). y_actual(2,nStep) = trnInputs(3); y_actual(3,nStep) = trnInputs(4); y_actual(4,nStep) = trnInputs(5); y_actual(5,nStep) = trnInputs(6); y_actual(6,nStep) = trnInputs(7); y_actual(7,nStep) = trnInputs(8); y_actual_control(:,nStep-1)=[y_actual(:,nStep);Twsnext(1,nStep-1)]-params.Yop; El segundo paso corresponde al cálculo de la acción de control, pero antes se debe proceder a la estimación del estado actual del sistema, obtenido gracias al observador de estados. Hay que tener en cuenta que se debe usar un sistema ampliado para poder lograr la obtención del parámetro d. Para ello se emplea un observador de la forma ˆx(k+1) = A0 0I−LC Iˆx(k)+B 0Lu(k) y(k).(7.13) La estimación calculada es la entrada al bloque de control, junto con la referencia que se quiere conseguir en cada tiempo de muestreo. Para comprobar que el control se produce de forma exitosa, se extrae del bloque de control también el valor del coste se tiene en cada tiempo de muestreo. Esto se puede ver como: % Se obtiene la estimacion de estados u_observador(:,nStep-1) = [u_actual(:,nStep-1); y_actual_control(:,nStep-1)]; x_observador(:,nStep) = (A_gorro-L*C_gorro)*x_observador(:,nStep-1) + [B_gorro L]*u_observador(:,nStep-1); % Se calcula los valores de la variable de control ctrlsol = ctrl([x_observador(1:nx+ny,nStep);yref(:,nStep)]); u_actual(:,nStep) = ctrlsol(1:nu); cost(:,nStep) = ctrlsol(nu+1); Por último, para la elaboración del bucle de control, solo es necesario aplicar las variables obtenidas a cada uno de los elementos físicos que representan cada una. Para ello, solo es necesaria la correcta asignación de cada una de las componentes sus elementos correspondientes. % Se asignan las variables de control Apertura(:,nStep) = u_actual(1:nu-1,nStep) + params.Uop(1:end-1); fcp(1,nStep) = u_actual(end,nStep) + params.Uop(end); El código anterior representa el último paso del bucle de control a implementar. Con esto, ya solo es necesario establecer y diseñar el esquema de simulación completo en TRNSYS y, asignar las correspondientes variables entre ellas. De esta forma que el esquema a implementar en dicho programa es que puede visualizarse en la Figura 7.2. Como se ha visto, el controlador genera unos valores de los parámetros de Apertura y fcp óptimos para realizar un control sobre el sistema en cuestión. Estos valores deben ser aplicados, en primer lugar, sobre los 68 Capítulo 7. Control MPC Figura 7.2 Esquema de implementación del controlador MPC en TRNSYS. equipos que generan o calculan la potencia necesaria que se debe aportar sobre cada una de las habitaciones. Para ello, las conexiones entre los bloques deben ser las siguientes: • Una conexión entre el bloque del tiempo y el del controlador MPC que permita la toma del dato de la temperatura exterior del ambiente, la cual es necesaria para caracterizar el funcionamiento de la enfriadora o chiller. • Por otro lado, la conexión realizada entre el edificio y el controlador que permite la aportación de cada una de las temperaturas de las habitaciones a dicho controlador. Estas son las entradas del controlador. Para obtener la temperatura del agua a la salida de la enfriadora, se realiza mediante la variable interna disponible en el archivo de MATLAB®. • De forma análoga, pero en sentido contrario, el bloque de control tiene que aportar cada una de las potencias calculadas a través de los datos de aperturas y fcp del controlador. Estas potencias son las encargadas de regular la temperatura de las habitaciones. • Por último, se tiene la conexión correspondiente entre el tiempo atmosférico y el edificio, la cual es necesaria para su caracterización y modelado de interacciones con el ambiente, es decir, el modelado de las pérdidas energéticas que se producen. Con todo esto, ya se puede proceder a realizar las correspondientes simulaciones las cuales van a verificar el correcto funcionamiento del controlador elaborado, además de mostras las evoluciones que presenten el resto de variables presentes en dicho controlador. 7.4 Resultados obtenidos Para corroborar que la implementación del controlador MPC explicado ha sido correcta, y comprobar así su funcionamiento, se realizará una simulación en TRNSYS con el objetivo de obtener una serie de resultados que puedan ser analizados con el objetivo de extraer una serie de conclusiones. Para dicha simulación se producirán diversos cambios en las temperaturas de las habitaciones para comprobar las respuestas que genera el sistema para llegar a cumplir dicha referencia que se ha establecido previamente. Las simulaciones realizadas tendrán una duración de 800 min , donde se producirán un total de cuatro cambios en las referencias de cada una de las salas, estos cambios se producen cada 200 min . Además, también se tiene como objetivo el comprobar que se satisfacen todas las restricciones que han establecido en la realización del control MPC. En estas simulaciones el control se aplica siempre sobre el sistema real y no sobre el linealizado. Los resultados de la simulación realizada donde se muestran las evoluciones de las temperaturas que se han obtenido, pueden ser observados en las Figuras 7.3 y 7.4. Como se observa, las referencias impuestas son seguidas y mantenidas correctamente y, el tiempo de establecimiento es extremadamente pequeño, por lo que el control es bastante aceptable. Se puede observar que las pequeñas perturbaciones provocadas por las dinámicas que representan las pérdidas por infiltraciones también son superadas con éxito. En las Figuras 7.5 y 7.6 se puede observar las distintas evoluciones que se producen en las variables de control. Si se observan dichas gráficas, se llega a la conclusión de que la rapidez que se puede observar en el seguimiento de las referencias se puede lograr gracias a los picos que se producen en las actuaciones. Cuando dichos picos se producen hacía un valor superior al que se encuentran (válvula se abre), se observa como la temperatura de la sala correspondiente disminuye, mientras que cuando se producen a un valor inferior (válvula se cierra), las temperaturas aumentan. 7.4 Resultados obtenidos 69 Figura 7.3 Resultados obtenidos de la implementacion del MPC (Salas 1 hasta 4). Figura 7.4 Resultados obtenidos de la implementacion del MPC (Salas 5 hasta 7). De forma algo similar, pero algo más compleja, puesto que influye sobre la temperatura de todas las salas, ocurre con la evolución del fcp. Sin embargo, lo que hay que destacar realmente, es como los valores de las aperturas de las válvulas siempre se encuentran en un valor superior al del fcp, lo que genera siempre un menor consumo que si se produjese en el sentido contrario. Este hecho provoca en la curva de la temperatura del agua algo menos de estabilidad, pues como el fcp siempre va a tender a ser más bajo que las aperturas, el agua a la salida de la enfriadora va a ver como su temperatura se verá incrementada, pero siempre manteniéndose 76 Capítulo 8. Ampliación del edificio Figura 8.5 Evolución de las temperaturas en el edificio de tres plantas (Salas 1 a 4). En los resultados que se obtienen se pueden extraer algunas ideas que verifican que la implementación llevada a cabo ha sido correcta. En primer lugar, se puede observar que las tendencias que siguen de las temperaturas son similares a las que se produce en la temperatura ambiente en un día completo, mostrando un aumento hasta las últimas horas de la tarde y una disminución desde entonces hasta las primeras horas de la mañana cuando se producen los primeros rayos de sol. Figura 8.6 Evolución de las temperaturas en el edificio de tres plantas (Salas 5 a 7). 8.2 Incorporación de equipos 77 Por otro lado, se puede observar como a última hora de la tarde las temperaturas de las salas situadas en la última planta sufren un notable aumento debido a la propiedad que tiene el aire caliente de dirigirse hacia las zonas superiores del edificio por tener una menor densidad cuando este aumenta su temperatura. De esta forma se puede observar el aumento de temperatura en la última planta, de forma que, el hecho de que no sean muy superiores a las de las otras plantas se puede explicar sabiendo que las condiciones iniciales de las que parte la simulación son idénticas en todas las salas, por lo que esa mayor diferencia de temperaturas entre una sala y otra, podrá ser visualizada si se simula durante bastante más tiempo. Con esto se concluye el diseño e implementación del edificio de estudio y se puede proceder a implementar la nueva máquina enfriadora que va a ser usada para lograr un acondicionamiento de temperatura en la sala. 8.2 Incorporación de equipos Para poder completar el modelado del sistema completo que debe ser analizado, se deben incluir los modelos de los equipos que componen la máquina enfriadora encargada de la refrigeración del edificio. Para ello se puede seguir como referencia el Capítulo 3, si bien habrá que tener en cuenta algunos detalles que serán matizados en este apartado. El equipo empleado para la refrigeración va a estar formado por los mismos equipos que cuando se trataba el caso de una única planta, lo cuales pueden enumerarse y definirse de una forma algo resumida como: • Un conjunto de válvulas de tres vías que permiten regular el caudal que es derivado a los distintos fancoils. El parámetro fundamental que modela el funcionamiento de cada una de las válvulas es la apertura, por lo que será una de las variables de control del conjunto de equipos completo (máquina enfriadora). Como el nuevo edificio está constituido por 21 salas, se deben considerar por tanto un total de 21 válvulas de tres vías. • Otros elementos clave de la máquina enfriadora son los distintos fancoils o intercambiadores de calor que permiten realizar el intercambio energético entre el fluido que los atraviesa y el aire de las salas, absorbiendo el fluido cierta cantidad de calor de las salas para poder disminuir así su temperatura. La máquina enfriadora estará formada por un total de 21 equipos como estos. • El elemento clave que garantiza el correcto funcionaniento es la enfriadora o chiller cuya función es disminuir la temperatura del agua del equipo completo para poder garantizar así que esta pueda absorber el calor de la sala. Es aquí donde se puede observar una de las modificaciones respecto al caso con una única sala, pues la potencia elegida para dicha máquina ya no es válida para el nuevo sistema puesto que las dimensiones se han triplicado y la potencia demanda por el nuevo edificio es muy superior a la anterior. Por ello se ha elegido un chiller con una capacidad nominal de 73.5kW . La variable que permite modelar su funcionamiento es el factor de carga parcial o fcp que determina la cantidad de potencia que se emplea para enfriar el agua que circula por su interior. Hay que destacar que las ecuaciones y los coeficientes que modelan el comportamiento de cada uno de los equipos se siguen manteniendo ya que el funcionamiento de los equipos no cambia, ya que únicamente se hace uso de más unidades para poder adaptar la máquina al nuevo edificio. El esquema que recoge como se deben conectar todos los equipos y el conjunto de variables que usa cada uno de ellos tampoco presenta ninguna diferencia con el del caso nuevo, dicho diagrama es el que se puede observar en la Figura 3.3. Con todo lo que se ha enunciado hasta ahora, el código en MATLAB ® que permite implementar los equipos de una forma correcta sigue el mismo esquema que los casos anteriores realizando una pequeña modificación en parámetros como son el caudal de agua total del equipo que ahora se puede cuantificar en 210 l/min , de forma que sigan circulando 10 l/min por cada válvula. Además, el número de veces que hay que llamar a cada una de las válvulas y fancoils pasan de ser 7 a 21, ya que se ha triplicado el número de estos dispositivos. Esto puede observarse en el siguiente fragmento de código: %% Tratamos en primer lugar a la valvula for i=1:21 [qsi(i,nStep) qbi(i,nStep)] = Valve(qin(i,1),Apertura(i,nStep)); end %% Ahora se intercambia calor entre la habiacion y el fancoil for i=1:21 78 Capítulo 8. Ampliación del edificio [Potencias(i,nStep) Tout(i,nStep)] = fanCoil(Twsnext(1,nStep-1),qsi(i,nStep )); end Para comprobar que se ha realizado correctamente, se va a realizar una simulación en la que se lleven a cabo distintos cambios en los valores que toman las aperturas de las válvulas y poder observar así la influencia que tiene cada una de las variables (Aperturas y fcp) sobre la temperatura de todas las salas del edificio. Figura 8.7 Cambios que se producen por la apertura de válvulas en el nuevo modelo (Salas 1 a 4). Figura 8.8 Cambios que se producen por la apertura de válvulas en el nuevo modelo (Salas 5 a 7). 8.2 Incorporación de equipos 79 En las Figuras 8.7 y 8.8 aparecen representada cada una de las evoluciones que siguen las temperaturas de las salas existentes. En relación a la nomenclatura de la leyenda hay que puntualizar el significado de PB (Planta Baja), PM (Planta Media) y PA (Planta Alta). De las figuras pueden extraerse algunas ideas que verifiquen que la implementación del modelo de la máquina enfriadora ha sido el correcto: • La primera conclusión que se puede obtener con una vista rápida sobre las gráficas, se encuentra en que se puede afirmar que, cuando la apertura de una válvula se modifica hacia un valor superior (se abre más la válvula), la temperatura de la sala o espacio disminuye, mientras que si esta se cierra se produce un aumento de la temperatura. Al igual que ocurría con el modelo de una planta, este fenómeno se puede explicar de forma física, por lo que se concluye que se ha realizado correctamente. • En las gráficas se puede observar cómo no todas ellas parten de la misma temperatura, lo cual es lógico pues el valor inicial de las aperturas no es el mismo y por ello estas han variado durante instantes previos a la serie de datos que han sido escogidos para realizar la representación, ya que lo representado no es más que un fragmento de la simulación realizada. La otra variable que influye en el comportamiento de la máquina enfriadora es el factor de carga parcial, que para observar su influencia sobre la temperatura de las salas se va a realizar una simulación en la que se aumente su valor (de un 0.2 a un 0.3) , de forma que la temperatura del agua a la salida del chiller disminuya considerablemente. Figura 8.9 Cambios que se producen por la variación de fcp en el nuevo modelo. Como se puede observar en los resultados de la Figura 8.9, cuando dicho valor se aumenta, se puede ver que las temperaturas de las habitaciones disminuyen, ya que la temperatura a la que sale el agua de la enfriadora disminuye, logrando así, que la potencia absorbida por el agua tras su paso por el fancoil también aumente, generando por tanto una disminución en la temperatura de la sala. También se puede destacar que el enfriamiento producido en las salas que representan los pasillos es inferior al del resto de salas, ya que el espacio que representa dichas salas es superior al del resto, por lo que este necesita mayor potencia para lograr el mismo enfriamiento. De esta forma se da por concluida la implementación de los equipos que permiten enfriar las salas, con lo que ya se dispone del modelo del nuevo edificio con los nuevos equipos, que constituyen el nuevo sistema completo. Ahora es necesario obtener un modelo matemático mediante identificación que permita modelar el sistema completo. 80 Capítulo 8. Ampliación del edificio 8.3 Identificación del sistema Para llevar a cabo la identificación, al igual que se hizo en el caso de una planta, se va a recurrir al toolbox de identificación que está disponible en MATLAB ® . Antes de plantear el modelo matemático y definir el conjunto de variables que lo conforman, es importante saber cuáles van a ser las entradas y las salidas del sistema en cuestión. Al igual que ocurre con el caso de una planta, las entradas del sistema siguen siendo las aperturas de las válvulas y el fcp, mientras que las salidas del sistema son cada una de las temperaturas de las distintas salas existentes. De esta forma se disponen de un total de 22 entradas (21 aperturas de válvulas y el fcp) y de 21 salidas (todas las temperaturas de las salas). De esta forma ya se puede acotar el caso de estudio y poder establecer el modelo de espacio de estados que modelará el comportamiento del sistema completo. Estas ecuaciones pueden establecerse como x(k+1) = Ax(k)+Bu(k)(8.1a) y(k) = Cx(k),(8.1b) en el cual ya se pueden caracterizar cada una de las componentes: donde x∈Rnx , u∈R22 e y∈R21 . Como se puede observar ya solo es necesario conocer el número de estados que definen el sistema. Para poder determinar el espacio de estados se deben definir el primer lugar el conjunto de pruebas que se deben realizar para poder aplicar el algoritmo de identificación empleado. El objetivo principal de las simulaciones a realizar es poder ilustrar cómo afecta la variación de cada una de las entradas sobre las salidas del sistema, generando para ello una modulación en forma de escalón sobre las entradas de forma que el cambio de la entrada sea instantáneo y de una amplitud determinada, siendo como es lógico el resto de las entradas constantes durante toda la simulación. Adicionalmente, una de las condiciones de las pruebas a realizar es que la temperatura ambiente va a ser considerada constante durante los ensayos, de forma que posteriormente sea considerada como una fuente de perturbaciones que afecta a las salas. Para realizar las mencionadas pruebas, en los escalones sobre los valores de las aperturas van a producirse variaciones de un 0.3, siendo el aumento en el valor del factor de carga parcial de un 0.1. De esta forma se deben realizar un total de 22 pruebas, una por cada entrada disponible, siendo estas muy similares a las que se realizaban en el Capítulo 4 (más concretamente en la Figura 4.3) Una vez que se hayan realizado todas las pruebas en cuestión se puede aplicar el código de identificación estudiado que permite determinar las matrices que conforman el modelo de espacio de estados. Tras realizar diversas pruebas con el número de estados a utilizar, se ha llegado a la conclusión de que los mejores resultados se obtienen para un total de 22 estados. Con esto, el parámetro nx ya queda definido de una forma adecuada. Si se aplica el código de identificación y se representan los resultados obtenidos, estos se pueden mostrar en las Figuras 8.10 y 8.11. En ellas se puede corroborar que el sistema ha sido indentificado de forma adecuada, ya que cuando se produce una modificación en las aperturas de las válvulas se puede observar cómo únicamente afecta a la sala donde se produce esa modificación. Por otro lado, cuando se produce una modificación en el valor del factor de carga parcial, se puede observar como este influye en el valor de todas las temperaturas de las habitaciones. Con esto ya se ha conseguido una identificación que refleja fielmente el comportamiento del sistema y puede hacerse uso de él para poder aplicar los controladores que han sido aplicados en el caso de una planta. 8.4 Control LQR Para comenzar la aplicación de estrategias de control con el objetivo de intentar lograr y mantener en una referencia deseada en las temperaturas de las habitaciones existentes en el edificio, se va a proponer el uso de un controlador LQR (Linear Quadratic Regulator). Para dicho controlador se parte del conocimiento del modelo matemático de espacio de estados que se ha elaborado en el apartado anterior, este modelo ha sido definido como x(k+1) = Ax(k)+Bu(k)(8.2a) y(k) = Cx(k),(8.2b) 8.4 Control LQR 81 Figura 8.10 Resultados de la identificación realizada (Entradas 1 a 10). 82 Capítulo 8. Ampliación del edificio Figura 8.11 Resultados de la identificación realizada (Entradas 11 a 21). 8.4 Control LQR 83 siendo x∈R22 , u∈R22 e y∈R21 . Para poder estimar los estados reales en los que se encuentra el sistema, se debe emplear un observador de estados que permita obtenerlos a partir de las variables de control y las salidas que estas generan en el sistema. La estructura del observador de estados que se emplea puede enunciarse matemáticamente de la forma ˆx(k+1) = (A−LC)ˆx(k)+B Lu,(8.3) de forma que u es un vector en forma de columna que almacene los valores de las actuaciones y las salidas (en ese orden) del instante k. De esta forma ya se puede estimar los distintos parámetros que en cierta forma definen el estado en el que se encuentra el sistema. La estrategia de control que se emplea está basada en una modificación del típico controlador LQR comúnmente conocido, denotado como control LQR con efecto integral. En dicho controlador, la ley de control que debe aplicarse en cada tiempo de muestreo se puede enunciar como ∆u(k) = Kz∆x(k) e(k−1)=Kx∆x(k)+Kee(k−1),(8.4) en la que las constantes proporcionales por las que se ven multiplicadas el error y la variación de los estados, resultan de aplicar la ecuación de Ricatti sobre un nuevo sistema ampliado el cual resulta ser ∆x(k+1) e(k)=A0 −C I ∆x(k) e(k−1)+B 0∆u(k) ∆y(k) = C∆x(k).(8.5) Mediante este modelado, ya se puede proceder a la implementación de los controladores en MATLAB ® como se ha definido en el Capítulo 6, realizando pequeñas modificaciones de las dimensiones en algunos de los parámetros que modelan en controlador, para poder adaptarlo ahora a las nuevas dimensiones del sistema ampliado. Figura 8.12 Temperaturas y Referencias con un controlador LQR (Planta Baja). Para comprobar el correcto funcionamiento del controlador LQR implementado, se va a realizar una simulación en la que se producirán un total de cuatro cambios de referencia cada 200 min , de forma que se pueda observar cómo evoluciona cada temperatura llegando a alcanzar el valor deseado. Los resultados que 84 Capítulo 8. Ampliación del edificio se obtienen son los que se pueden observar en las Figuras 8.12, 8.13 y 8.14. Figura 8.13 Temperaturas y Referencias con un controlador LQR (Planta Media). Con los resultados obtenidos, se puede afirmar que el controlador ha sido diseñado e implementado correctamente, de forma que otorga a las actuadores los valores que permiten aportar la potencia necesaria para lograr el mantenimiento de las salas a la temperatura deseada. Figura 8.14 Temperaturas y Referencias con un controlador LQR (Planta Alta). 8.5 Control MPC 85 Adicionalmente, para mostrar la coherencia de los resultados obtenidos se va a representar en la Figura 8.15 la comparación entre los valores que toman el factor de carga parcial (fcp) y la temperatura del agua a la salida de la enfriadora o chiller. En dicha imagen se puede observar que en esos instantes los valores del fcp son muy bajos, lo que no quiere decir que se haya usado una máquina demasiado potente, ya que al principio de la simulación (durante el transitorio que no se ha mostrado), cuando las temperaturas de las salas evolucionan con la temperatura del ambiente y el agua que circula por la máquina refrigeradora todavía no se encuentra acondicionada a la temperatura de trabajo, el esfuerzo que debe realizar el chiller es alto y demanda por tanto mayores valores del fcp que son cercanos a la unidad. Por otro lado, en la imagen se puede observar como la temperatura del agua se mantiene siempre entre los niveles correctos de funcionamiento, situados estos entre 5◦C y 15◦C , de forma que cuando la temperatura vaya a sobrepasar tanto el límite inferior como superior, se produce una variación positiva o negativa para evitar este fenómeno. Este hecho justifica la tendencia que siguen ambas curvas representadas y aporta consistencia también a los resultados obtenidos. Se debe destacar que los valores de las aperturas de cada una de las válvulas existentes no han sido mostrados para evitar así un elevado número de gráficas en este memoria meramente descriptiva, pero si se debe destacar que los valores obtenidos muestran unos valores que concuerdan con un aumento o disminución de sus valores, en función se quieran disminuir o aumentar las temperaturas respectivamente. Figura 8.15 Comparativa entre los valores de fcp y de la temperatura a la salida del chiller (Control LQR). Con estos resultados obtenidos se puede dar por finalizado este apartado relativo a la implementación y diseño de un controlador LQR con efecto integral y se puede ver la importancia de diseñarlo de forma adecuada para poder generar así un enfriamiento de las habitaciones generando el menor gasto energético posible. 8.5 Control MPC Como se ha estudiado en capítulos anteriores, se va a estudiar cómo se comporta y cómo se diseña un controlador predictivo MPC para el sistema en cuestión. En este nuevo objeto de estudio se permite incluir restricciones menos severas en el seguimiento de referencias por parte de las salidas, lo que permite poder establecer unos límites en los valores que toma la temperatura del agua a la salida de chiller, ya que este parámetro es muy crucial a la hora de caracterizar el comportamiento de la enfriadora o chiller. Los valores de la temperatura del agua que marcan el buen funcionamiento de los equipos varían entre 5◦Cy15◦C. Bibliografía [1] A. Cardona, D. Limón, J.F. Coronel, and L. Pérez-Lombard, Creación de casos de estudio para la gestión energética de sistemas HVAC en edificios basados en trnsys, 2019. [2] Limón D., Diapositivas de la asignatura Control de Procesos Industriales, 4 º GIERM, Escuela Técnica Superior de Ingeniería, Universidad de Sevilla. [3] Rodríguez D., Diapositivas de la asignatura Complementos de Control, 4 º GITI, Escuela Técnica Superior de Ingeniería, Universidad de Sevilla. [4] A.M. Foster, M.J. Swain, R. Barrett, and S.J. James, Experimental verification of analytical and CFD predictions of infiltration through cold store entrances. [5] E. Hendricks, O. Jannerup, and P. Sørensen, Linear systems control deterministic and stochastic methods, 1ºed., 2008. [6] Ogata K, Sistemas de control en tiempo discreto, 2ºed., 1996. [7] Ogata K., Ingeniería de control moderna, 5ºed., 2010. [8] S.A. Klein, W.A. Beckman, J.W. Mitchell, J.A. Duffie, N.A. Duffie, J.C. T.L. Freeman, Mitchell, J.E. Braun, B.L. Evans, J.P. Kummer, R.E. Urban, A. Fiksel, J.W. Thornton, N.J. Blair, P.M. Williams, D.E. Bradley, T.P. McDowell, M. Kummert, and D.A. Arias, Trnsys 16 a transient system simulation program. Volume 1, Getting Started., 5ºed., 2007. [9] Daniel Limón, Ignacio Alvarado, Teodoro Alamo, and Eduardo F Camacho, Mpc for tracking piecewise constant references for constrained linear systems, Automatica 44 (2008), no. 9, 2382–2387. [10] Lennart Ljung, System identification, Wiley Encyclopedia of Electrical and Electronics Engineering (1999), 1–19. [11] J.G. Ordonez, J.M. Nadales, D.Limón, Luis Perez-Lombard, Juan F. Coronel, and A. Cardona, Modelo de sistema de refrigeración de edificio orientado al control y gestión de la eficiencia energética, 2020. [12] José L. Vázquez Otero, Análisis de pérdidas energéticas por infiltración de aire, IES Universidade Laboral de Culleredo, 2016. [13] Luis Pérez-Lombard, José Ortiz, and Christine Pout, A review on buildings energy consumption information, Energy and buildings 40 (2008), no. 3, 394–398. [14] A Picard, R S Davis, M Glaser, and K Fujii, Revised formula for the density of moist air (CIPM-2007), 2008. [15] Rawlings and James B., Model predictive control: Theory and design / by james b. rawlings and david q. mayne., 2013. [16] Página web de Google SketchUp, https://www.sketchup.com/es/plans-and-pricing/sketchup-free. [17] Página web de la International Energy Agency, https://www.iea.org/, 2020. 93 94 Bibliografía [18] Página web del distrubidor de TRNSYS español, https://aiguasol.coop/es/software-energia/trnsys-18- simulacion-energetica/. [19] Página web Multi Parametric Toolbox (MPT3), https://yalmip.github.io/solver/mpt/. [20] University of Wisconsin-Madison © 2005 by the Solar Energy Laboratory, Trnsys 17 a transient system simulation program. Volume 5, Multizone Building modeling with Type56 and TRNBuild, 7 º ed., 2012. [21] University of Wisconsin Madison © 2005 by the Solar Energy Laboratory, Trnsys 17 a transient system simulation program. volume 3, standard component library overview, 8ºed., 2014. [22] University of Wisconsin-Madison © 2005 by the Solar Energy Laboratory, Trnsys 17 a transient system simulation program. Volume 4, Mathematical Reference, 10ºed., 2014.