Identificación de modelos dinámicos de sistemas de refrigeración con eyectores de geometría variable

Equation Chapter 1 Section 1 Proyecto Fin de Grado Ingeniería de las Tecnologías Industriales Identificación de modelos dinámicos de sistemas de refrigeración con eyectores de geometría variable Dpto. Ingeniería Energética Escuela Técnica Superior de Ingeniería Universidad de Sevilla Autor: Sofía Carmona Cuadrado Tutor: José Manuel Salmerón Lissén Co-tutor: Daniel Limón Marruedo Sevilla, 2021 iii Proyecto Fin de Grado Ingeniería de las Tecnologías Industriales Identificación de modelos dinámicos de sistemas de refrigeración con eyectores de geometría variable Autor: Sofía Carmona Cuadrado Tutor: José Manuel Salmerón Lissén Profesor titular Co-tutor: Daniel Limón Marruedo Catedrático Dpto. Ingeniería Energética Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2021 v Proyecto Fin de Carrera: Identificación de modelos dinámicos de sistemas de refrigeración con eyectores de geometría variable Autor: Sofía Carmona Cuadrado Tutor: José Manuel Salmerón Lissén Co-tutor: Daniel Limón Marruedo El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: Sevilla, 2021 El Secretario del Tribunal vii A mi familia A mis amigos A mis profesores ix Agradecimientos A mi familia, por estar siempre a mi lado y apoyarme constantemente. Por todas las oportunidades que me han brindado a lo largo de mi vida para llegar donde estoy y hacerme ser la persona que soy hoy. A mis tutores Daniel Limón y José Manuel Salmerón, que se han volcado conmigo desde el primer día a pesar de los obstáculos que se han presentado durante este proceso. También a mis amigos y compañeros de universidad, por su cariño, ánimos y consejos. Por acompañarnos mutuamente en este camino y por todos los momentos de estudio y desconexión vividos junto a ellos. Finalmente, agradecer a todos los profesores que han formado parte de mi etapa de aprendizaje dentro de la universidad, por aportarme las bases necesarias para afrontar cualquier reto. Muchísimas gracias a todos. Sofía Carmona Cuadrado Sevilla, 2021 4.2.3 Modelos de procesos 63 4.2.4 Modelos de funciones de transferencia 64 4.2.5 Modelos no lineales 66 4.3 Modelos de Caja Gris 70 4.3.1 Sistema inicial de segundo orden 70 4.3.2 Sistema inicial de orden configurable 71 5 Aplicación a un sistema de tres eyectores 73 5.1 El ciclo MERS 73 5.2 Evaluación del rendimiento estacional 74 5.3 Análisis de viabilidad del MERS 74 6 Conclusiones 77 Bibliografía 81 xvii Índice de Tablas Tabla 1. Calendario de eliminación gradual del consumo de sustancias HCFC 1 Tabla 2. Características de los eyectores según su fase y flujo 9 Tabla 3. Evolución de la eficiencia de un eyector 11 Tabla 4. Características de los principales refrigerantes. 13 Tabla 5. Ceros, polos y eficiencia tras estimar modelos de función de transferencia 65 Tabla 6. Regresores y eficiencia tras estimar modelos ARX no lineales 67 Tabla 7. Ceros, polos, eficiencia y unidades tras estimar modelos de Hammerstein-Wiener 68 xix Índice de Figuras Figura 1. Diagrama esquemático del TCRE [17]. 2 Figura 2. Comparación de la respuesta dinámica para un PID fijo y uno autoadaptativo. 4 Figura 3. Ciclo de compresión mecánica. (a) componentes del ciclo y (b) diagrama p-h. 6 Figura 4. Principales partes del eyector 6 Figura 5. Tipos de boquilla. 7 Figura 6. Modo de funcionamiento de un eyector subsónico (a) presión primaria fija y (b) contrapresión fija 8 Figura 7. Modo de funcionamiento de un eyector supersónico (a) presión primaria fija y (b) contrapresión fija 8 Figura 8. Tipos de cámaras de mezcla (a) CPM y (b) CAM 9 Figura 9. Boquilla variable 10 Figura 10. Parámetros geométricos de un eyector 11 Figura 11. Diagramas T-s de tres tipos de refrigerantes: (a) fluido húmedo, (b) fluido seco, y (c) fluido isentrópico 12 Figura 12. Diferentes configuraciones de ciclos de refrigeración con eyetores 14 Figura 13. Ciclo de refrigeración con eyector 14 Figura 14. Proceso de identificación. 16 Figura 15. Diagrama de bloques identificación en bucle abierto. 17 Figura 16. Diagrama de bloques identificación en bucle cerrado. 17 Figura 17. Diagrama de bloques de una estructura ARX no lineal 20 Figura 18. Diagrama de bloques de una estructura Hammerstein-Wiener 21 Figura 19. Curva de puntos experimentales (a) tomados a mano y (b) filtrados 24 Figura 20. Modelo de caja negra y sistema de primer orden. Resultados obtenidos en el primer escalón en (a) sistema controlado y (b) respuesta ante escalón unitario 27 Figura 21. Modelo de caja negra y sistema de primer orden. Resultados obtenidos en el segundo escalón en (a) sistema controlado y (b) respuesta ante escalón unitario 28 Figura 22. Modelo de caja negra y sistema de primer orden. Resultados obtenidos en el tercer escalón en (a) sistema controlado y (b) respuesta ante escalón unitario 29 Figura 23. Modelo de caja negra y sistema de primer orden. Resultados obtenidos en el cuarto escalón en (a) sistema controlado y (b) respuesta ante escalón unitario 29 Figura 24. Modelo de caja negra y sistema de primer orden. Resultados obtenidos en sistema completo en (a) sistema controlado y (b) respuesta ante escalón unitario 30 Figura 25. Modelo de caja negra y sistema de segundo orden. Resultados obtenidos en el primer escalón en (a) sistema controlado y (b) respuesta ante escalón unitario 32 Figura 26. Modelo de caja negra y sistema de segundo orden. Resultados obtenidos en el segundo escalón en (a) sistema controlado y (b) respuesta ante escalón unitario 32 Figura 27. Modelo de caja negra y sistema de segundo orden. Resultados obtenidos en el tercer escalón en (a) sistema controlado y (b) respuesta ante escalón unitario 33 Figura 28. Modelo de caja negra y sistema de segundo orden. Aplicación del par obtenido en el tercer escalón al sistema completo 34 Figura 29. Modelo de caja negra y sistema de segundo orden. Resultados obtenidos en el cuarto escalón en (a) sistema controlado y (b) respuesta ante escalón unitario 34 Figura 30. Modelo de caja negra y sistema de segundo orden. Aplicación del par obtenido en el cuarto escalón al sistema completo 35 Figura 31. Modelo de caja negra y sistema de segundo orden. Resultados obtenidos en el sistema completo en (a) sistema controlado y (b) respuesta ante escalón unitario 36 Figura 32. Modelo de caja negra y sistema de orden configurable. Resultados obtenidos en el primer escalón en el (a) sistema controlado y (b) respuesta ante escalón unitario 37 Figura 33. Modelo de caja negra y sistema de orden configurable. Aplicación del par obtenido en el primer escalón al sistema completo 37 Figura 34. Modelo de caja negra y sistema de orden configurable. Resultados obtenidos en el segundo escalón en el (a) sistema controlado y (b) respuesta ante escalón unitario 38 Figura 35. Modelo de caja negra y sistema de orden configurable. Aplicación del par obtenido en el segundo escalón al sistema completo 38 Figura 36. Modelo de caja negra y sistema de orden configurable. Resultados obtenidos en el tercer escalón en el (a) sistema controlado y (b) respuesta ante escalón unitario 39 Figura 37. Modelo de caja negra y sistema de orden configurable. Aplicación del par obtenido en el tercer escalón al sistema completo 40 Figura 38. Modelo de caja negra y sistema de orden configurable. Resultados obtenidos en el cuarto escalón en el (a) sistema controlado y (b) respuesta ante escalón unitario 40 Figura 39. Modelo de caja negra y sistema de orden configurable. Aplicación del par obtenido en el cuarto escalón al sistema completo 41 Figura 40. Modelo de caja negra y sistema de orden configurable. Resultados obtenidos en el sistema completo en (a) sistema controlado y (b) respuesta ante escalón unitario 42 Figura 41. Modelo de caja negra y sistema de orden configurable. Resultados obtenidos en el sistema multimodelo 43 Figura 42. Modelo de caja gris y sistema de primer orden. Resultados obtenidos en el primer escalón en (a) sistema controlado y (b) respuesta ante escalón unitario 44 Figura 43. (a) Representación del Lugar de las Raíces. (b) Respuesta ante escalón unitario 45 Figura 44. Modelo de caja gris y sistema de 1 cero y 2 polos. Resultados obtenidos en el primer escalón en (a) sistema controlado y (b) respuesta ante escalón unitario 47 Figura 45. Modelo de caja gris y sistema de 1 cero y 2 polos. Aplicación del par obtenido en el primer escalón al sistema completo 47 Figura 46. Modelo de caja gris y sistema de 1 cero y 2 polos. Resultados obtenidos en el segundo escalón en (a) sistema controlado y (b) respuesta ante escalón unitario 48 Figura 47. Modelo de caja gris y sistema de 1 cero y 2 polos. Aplicación del par obtenido en el segundo escalón al sistema completo 48 Figura 48. Modelo de caja gris y sistema de 1 cero y 2 polos. Resultados obtenidos en el tercer escalón en (a) sistema controlado y (b) respuesta ante escalón unitario 49 Figura 49. Modelo de caja gris y sistema de 1 cero y 2 polos. Aplicación del par obtenido en el tercer escalón al sistema completo 49 xxi Figura 50. Modelo de caja gris y sistema de 1 cero y 2 polos. Resultados obtenidos en el cuarto escalón en (a) sistema controlado y (b) respuesta ante escalón unitario 50 Figura 51. Modelo de caja gris y sistema de 1 cero y 2 polos. Aplicación del par obtenido en el cuarto escalón al sistema completo 50 Figura 52. Modelo de caja gris y sistema de 1 cero y 2 polos. Resultados obtenidos en el sistema completo en (a) sistema controlado y (b) respuesta ante escalón unitario 51 Figura 53. Modelo de caja gris y sistema de 1 cero y 2 polos. Resultados obtenidos en el sistema multimodelo 52 Figura 54. Modelo de caja gris y sistema de 1 cero y 3 polos. Resultados obtenidos en el segundo escalón en (a) sistema controlado y (b) respuesta ante escalón unitario 53 Figura 55. Modelo de caja gris y sistema de 1 cero y 3 polos. Aplicación del par obtenido en el primer escalón al sistema completo 53 Figura 56. Modelo de caja gris y sistema de 1 cero y 3 polos. Resultados obtenidos en el segundo escalón en (a) sistema controlado y (b) respuesta ante escalón unitario 54 Figura 57. Modelo de caja gris y sistema de 1 cero y 3 polos. Aplicación del par obtenido en el segundo escalón al sistema completo 54 Figura 58. Modelo de caja gris y sistema de 1 cero y 3 polos. Resultados obtenidos en el tercer escalón en (a) sistema controlado y (b) respuesta ante escalón unitario 55 Figura 59. Modelo de caja gris y sistema de 1 cero y 3 polos. Aplicación del par obtenido en el tercer escalón al sistema completo 55 Figura 60. Modelo de caja gris y sistema de 1 cero y 3 polos. Resultados obtenidos en el cuarto escalón en (a) sistema controlado y (b) respuesta ante escalón unitario 56 Figura 61. Modelo de caja gris y sistema de 1 cero y 3 polos. Aplicación del par obtenido en el tercer escalón al sistema completo 56 Figura 62. Modelo de caja gris y sistema de 1 cero y 3 polos. Resultados obtenidos en el sistema completo en (a) sistema controlado y (b) respuesta ante escalón unitario 57 Figura 63. Modelo de caja gris y sistema de 1 cero y 3 polos. Resultados obtenidos en el sistema multimodelo 58 Figura 64. Curva resultante de la interpolación. 60 Figura 65. Interfaz de la System Identification App 61 Figura 66. Selección de estructura ARX (a) panel de selección de pares a estimar y (b) pares recomendados 61 Figura 67. Comparación de las estructuras de modelos ARX con la deseada 62 Figura 68. Selección de estructura de espacios de estado (a) panel de selección del orden a estimar y (b) orden recomendado 62 Figura 69. Comparación de las estructuras de modelos de espacios de estado con la deseada 63 Figura 70. Panel de selección de la estructura del modelo de procesos 63 Figura 71. Comparación de las estructuras de modelos de procesos con la deseada 64 Figura 72. Panel de selección de la estructura del modelo de funciones de transferencia 64 Figura 73. Curvas resultantes tras estimar a partir de modelos de funciones de transferencia. 65 Figura 74. Panel de configuración de los modelos ARX no lineales (a) selección de regresores y (b) selección de propiedades 66 Figura 75. Curvas resultantes tras estimar a partir de modelos ARX no lineales. 67 Figura 76. Panel de configuración de los modelos H-W. Selección de bloques no lineales y lineal 68 Figura 77. Curvas resultantes tras estimar a partir de modelos Hammerstein-Wiener. 69 Figura 78. Identificación de caja gris utilizando el Toolbox. Sistema inicial de segundo orden 71 Figura 79. Identificación de caja gris utilizando el Toolbox. Sistema inicial de orden configurable 71 Figura 80. Diagrama del ciclo MERS 73 Figura 81. Simulación estacional del rendimiento del sistema multieyector [95] 74 Figura 82. Estados estables para diferentes áreas de la garganta de la boquilla [17]. 75 Figura 83. Histórico anual de temperaturas y cambios de eyector 75 Figura 84. Histórico de la semana del 26-03 de temperaturas y cambios de eyector 76 Figura 85. Histórico de la semana del 17-07 de temperaturas y cambios de eyector 76 xxiii Notación Ar Area ratio (relación del área) C controlador cc Condensation critical point (punto crítico de condensación) e Perturbaciones/ruido G Función de transferencia del modelo dinámico H Función que modela el ruido k Propiedades del ruido kp Ganancia proporcional M Número de Mach na Polos nb Ceros nk Retraso o Inicial r Referencia externa Rc Relación de compresión S Función de sensibilidad de salida Si Función de sensibilidad de entrada t Tiempo Td Tiempo derivativo Ti Tiempo integral u Entrada x Vector de variables de estado y Salida Z Set de datos μ Relación de arrastre ω Velocidad del compresor SIGLAS Y ABREVIATURAS ARMAX Autoregressive moving-average exogenous model (modelo autorregresivo exógeno de media móvil) ARX Autoregressive exogenous model (modelo autorregresivo exógeno) BA Bucle abierto BC Bucle cerrado BJ Box-Jenkins model (modelo de Box-Jenkins) CAM Contant area mixing chamber (cámara de mezcla de área constante) CC Cooling capacity (capacidad de enfriamiento) CFC Clorofluorocarburos COP Coefficient of performance (coeficiente de rendimiento) CPM Constant pressure mixing chamber (cámara de mezcla a presión constante) CRMC Constant rate of momentum change method (método de la tasa constante de cambio del momento) EES Engineering equation solver (solucionador de ecuaciones de ingeniería) ERS Ejector refrigeration system (sistema de refrigeración por eyección) GEI Gases de efecto invernadero GWP Global warming potential (potencial de calentamiento global) HCFC Hidrocloroflurorocarburos HFC Hidrofluorocarburos IFAC International federation of automatic control (federación internacional de control automático) MERS Multiple ejector refrigeration system (sistema de refrigeración multieyector) MIMO Multiple input – multiple output (entrada multiple – salida multiple) ODP Ozone depletion potential (potencial de agotamiento del ozono) OE Output error (error de salida) PEM Prediction error minimization (minimización de la predicción del error) RA Relación de arrastre SISO Single input – single output (una única entrada – una única salida) SITB System identification toolbox (herramienta de identificación de sistemas) SYSID Symposium of system identification (simposio de identificación de sistemas) TCRE Transcritical CO2 refrigeration system with an ejector (sistema de refrigeración de CO2 transcrítico con eyector) TERS Transcritical ejector refrigeration system (sistema de refrigeración transcrítico con eyector) 1 INTRODUCCIÓN a producción de frío ha sido una necesidad para el ser humano desde la antigüedad. En 1805 Oliver Evans diseñó la primera máquina de refrigeración por compresión utilizando vapor en lugar de líquido e identificando sus principales componentes. En cuanto al primer aparato de aire acondicionado moderno, Willis Carrier desarrolló el concepto de climatización al estudiar cómo evitar el aumento de la humedad relativa del aire enfriado en 1902. Actualmente, los sistemas de refrigeración desempeñan un papel importante utilizándose para la conservación de alimentos, para la climatización de edificios residenciales, comerciales o vehículos y en sectores como el industrial o el sanitario. El desarrollo técnico y el crecimiento económico de la mayoría de países han favorecido a la incorporación de estos sistemas y se estima que a lo largo del SXXI, la demanda global de refrigeración y climatización aumente en un 72% [1]. Los ciclos tradicionales se accionan con electricidad o calor, por lo que el consumo de energía eléctrica y combustibles fósiles seguirán aumentando considerablemente. Las estimaciones del Instituto Internacional de Frío de París (IIF/IIR) calcula que el 15% de toda la electricidad producida en el mundo se emplea en procesos de refrigeración y aires acondicionados, pudiendo llegar a crecer hasta el 17% en un periodo de cincuenta años. Como consecuencia, la liberación de sustancias nocivas para el medio ambiente también se ha visto disparada, produciendo grandes problemas medioambientales: El agotamiento de la capa de ozono y el calentamiento global son los mayores efectos secundarios que se han manifestado. Las emisiones de gases de efecto invernadero pueden dividirse en emisiones directas e indirectas. Las emisiones directas se deben a las fugas de refrigerante y las indirectas a la emisión de CO2 al consumir energía obtenida través de la combustión de combustibles fósiles. Algunos refrigerantes empleados en los ciclos de refrigeración como los clorofluorocarburos (CFC), los hidroclorofluorocarburos (HCFC) y los hidrofluorocarburos (HFC) son causantes del agotamiento de la capa de ozono y del calentamiento global. En 1974, Rowland y Molina presentaron su descubrimiento de que el cloro y el bromo liberados por los productos químicos sintéticos migran a la estratosfera y destruyen la capa de ozono estratosférico que actúa como escudo contra la dañina radiación solar ultravioleta [2]. Tras varios años de negociaciones, en 1987 se firmó el Protocolo de Montreal para regular la producción y el comercio de sustancias que dañan y agotan la capa de ozono, incluyendo los refrigerantes mencionados anteriormente. El protocolo establece una eliminación gradual del uso de los HCFCs (Tabla 1) y una reducción de los HFCs. Si todos los países firmantes del acuerdo cumplen los objetivos propuestos, la capa de ozono deberá haberse recuperado para 2050. 1996 2004 2010 2015 2020 2030 Congelación del consumo Reducción del 35% Reducción del 65% Reducción del 90% Reducción del 99.5% Eliminación total del consumo Tabla 1. Calendario de eliminación gradual del consumo de sustancias HCFC El calentamiento global es el segundo efecto nocivo producido por refrigerantes. En el caso de que no se disminuyera la emisión de gases de efecto invernadero, en la primera mitad del próximo siglo podría producirse un aumento de la temperatura media global mayor que cualquier otro en la historia [3]. Por ello, en 1997 se firmó el Protocolo de Kioto que pone en funcionamiento la convención Marco de las Naciones Unidas sobre el Cambio Climático, comprometiendo a los 187 países industrializados firmantes a limitar y reducir las emisiones de gases de efecto invernadero (GEI) conforme a las metas acordadas individualmente. La Convención pide a los países la adopción de políticas y medidas de mitigación e informes periódicos acerca de su cumplimiento. La necesidad de encontrar una solución a largo plazo ha derivado en varias ideas principales: recurrir a fuentes de energías renovables, utilizar refrigerantes naturales como el CO2 para reducir las emisiones contaminantes y, en materia de mejora del COP del sistema, se propone reutilizar el calor desprendido durante la cogeneración a través de ciclos de absorción o eyectores. El eyector, que es el corazón del sistema de refrigeración por inyección, fue inventado por Sir Charles Parsons en 1901 para extraer aire del condensador de una máquina de vapor. En 1910, Maurice Leblanc utilizó un eyector L Fundamentos y estado del arte 8 8 Figura 6. Modo de funcionamiento de un eyector subsónico (a) presión primaria fija y (b) contrapresión fija La tobera convergente-divergente hace que el flujo que trasncurre a través del eyector lo haga a velocidades supersónicas, por encima de la velocidad del sonido (número de Mach, M>1), generando una gran diferencia de presión. De nuevo, puede trabajar en tres modos diferentes (Figura 7): • Modo crítico o de doble bloqueo: La RA se mantiene constante debido al bloqueo que subren los flujos primario y secundario. • Modo subcrítico o de bloqueo único: La RA es lineal con la contrapresión y únicamente se encuentra bloqueado el fujo primario. • Modo de mal funcionamiento o flujo de retorno: El flujo secundario se invierte provocando la disfuncionalidad del eyector. Es el eyector utilizado en los ERS. Figura 7. Modo de funcionamiento de un eyector supersónico (a) presión primaria fija y (b) contrapresión fija 2.1.2.2.2 Ca mara de mezcla Los eyectores también pueden clasificarse según el tipo de cámara de mezcla [8], [30], distinguiendo entre cámaras de mezcla a presión constante (CPM) y cámaras de mezcla de área constante (CAM). También se han realizado estudios [31], [32] tratando de combinar las dos tecnologías mencionadas, obteniendo como resultado un eyector con una cámara de mezcla a velocidad constante y de momento variable (CRMC). En los eyectores CPM, Figura 8.a, la salida de la boquilla y donde comienza el proceso de mezcla de los fluidos tienen lugar en la cámara de succión para más adelante atravesar una zona de área constante. Estos eyectores obtienen menores contrapresiones y como consecuencia, un mejor rendimiento. La mezcla de los fluidos en los eyectores CAM Figura 8.b comienza, a diferencia de los CPM, en la cámara de mezcla pasando en todo momento por una sección de área constante. Los eyectores CAM, aunque tengan un menor desempeño, pueden soportar mayores flujos de masa. 9 Figura 8. Tipos de cámaras de mezcla (a) CPM y (b) CAM El eyector CRMC modifica la geometría del difusor de modo que se elimina el choque termodinámico dentro de él en las condiciones de funcionamiento del punto de diseño. Esto se consigue permitiendo que el momento del flujo varíe a una tasa constante a medida que cruza el difusor, derivando en un aumento gradual de la presión estática desde la entrada hasta la salida. Como consecuencia, se evita la pérdida de presión total asociada al proceso de choque que se puede encontrar en difusores convencionales. 2.1.2.2.3 Fase del fluido El estado de la materia del fluido de trabajo tiene un impacto significativo en la dinámica del mismo [33], especialmente en las ondas de choque que pueden tener lugar en el interior de los eyectores. Nº FASES FLUJO PRIMARIO FLUJO SECUNDARIO FLUJO DE SALIDA NOMBRE DEL EYECTOR CARACTERÍSTICAS 1 Vapor Vapor Vapor Eyector de chorro de vapor Posible flujo bifásico Posibles ondas de choque 1 Líquido Líquido Líquido Eyector de chorro líquido Sin ondas de choque Flujo monofásico 2 Vapor Líquido Líquido Eyector de condensación Flujo bifásico con condensación del flujo primario Fuertes ondas de choque 2 Líquido Vapor Bifásico Eyector bifásico Flujo bifásico Posibles ondas de choque Tabla 2. Características de los eyectores según su fase y flujo El flujo dentro del eyector puede ser monofásico (vapor-vapor o líquido-líquido) o bifásico dependiendo de las condiciones de los flujos primario y secundario. Los eyectores monofásicos han sido muy estudiados en la literatura y referencias como [34]–[37]. Un eyector bifásico se puede a su vez clasificar en función de la naturaleza del flujo: En el eyector de condensación el flujo primario se condensa en el eyector y en el bifásico puro, el flujo a la salida es bifásico. Hoy en día, se ha avanzado mucho en el estudio de los eyectores bifásicos, por lo que su conocimiento y modelado se ha ampliado [29], [38]–[43]. En la Tabla 2 se pueden observar separadamente las características tanto de los eyectores monofásicos como de los bifásicos. Fundamentos y estado del arte 10 10 2.1.2.2.4 Eyector fijo o variable Un estudio de optimización [44] concluyó que, para un punto de operación concreto, solo existe una única geometría optima del eyector. Así, se sugirió la utilizacion de un eyector con geometría variable que se adaptara a diversos puntos de funcionamiento [6]. Los eyectores con geometría variable consisten en insertar un husillo antes de la boquilla primaria (Figura 9) y que cambie su posición dentro de la garganta, alterando su área de sección transversal y dando lugar al control del caudal másico primario. Otros autores han comprobado satisfactoriamente el uso de esta tecnología, recomendando su utilización al verse incrementado el rendimiento del eyector [14], [15], [17], [45]. Figura 9. Boquilla variable 2.1.2.3 Parámetros y su influencia en el rendimiento Algunos de los parámetros más importantes para describir el rendimiento de los eyectores en los ciclos de refrigeración son los siguientes. • Relación de arrastre, 𝜇: Es la relación entre el caudal másico del flujo secundario, 𝑚󰇗𝑠, y el caudal másico del flujo primario, 𝑚󰇗𝑝. Evalúa la eficiencia del ciclo de refrigeración. 𝜇 = 𝑚󰇗𝑠 𝑚󰇗𝑝 (5) • Relación de compresión, 𝑅𝐶: Es la presión estática a la salida del difusor, 𝑝𝑐, dividida por la presión estática del flujo secundario, 𝑝𝑒. 𝑅𝐶=𝑝𝑐 𝑝𝑒 (6) • Coeficiente de rendimiento, 𝐶𝑂𝑃: Representa la relación entre la energía térmica de evaporación, 𝑄𝑒 (efecto de refrigeración), y el total de energía entrante en el ciclo (𝑄𝑔+𝐿𝑝). 𝐶𝑂𝑃 = 𝑄𝑒 𝑄𝑔+𝐿𝑝 (7) • Capacidad de refrigeracion, 𝐶𝐶: Es 𝐶𝐶 =𝑚󰇗𝑒(ℎ𝑒,𝑜𝑢𝑡 −ℎ𝑒,𝑖𝑛) (8) • Rendimiento del eyector, 𝜂𝑒𝑦𝑒𝑐𝑡𝑜𝑟: Hay muchas maneras de definir la eficiencia del eyector. Esta expresión [14] lo define como la relación etre la energía de compresión real recuperada y la energía teórica disponible en la corriente motriz. 𝜂𝑒𝑦𝑒𝑐𝑡𝑜𝑟 =(𝑚󰇗𝑔+𝑚󰇗𝑒)(ℎ𝑐,𝑖𝑛 −ℎ𝑒,𝑜𝑢𝑡) 𝑚󰇗𝑔(ℎ𝑔,𝑜𝑢𝑡 −ℎ𝑒,𝑜𝑢𝑡) (9) Boquilla primaria Husillo 11 La eficiencia en los diferentes componentes de un eyector subsónico (boquilla primaria convergente) variando parámetros geométricos y condiciones de operación del eyector ha sido motivo de estudio [38]. Las partes analizadas son la boquilla primaria, la cámara de mezclas y el difusor. Los parámetros elegidos son la relación de presión entre la entrada del flujo primario con la del secundario; la relación del diámetro de salida del flujo primario con el de la cámara de mezcla; y la distancia existente entre la salida de la boquilla primaria y la entrada a la cámara de mezcla (NXP). Llegando a las conclusiones expuestas en la Tabla 3. EVOLUCIÓN DE LA EFICIENCIA BOQUILLA PRIMARIA CÁMARA DE MEZCLA DIFUSOR AUMENTO DE LA RELACIÓN DE PRESIÓN Decrece Crece hasta llegar a un punto en que se mantiene, pudiendo llegar a decrecer si aumenta demasiado Decrece, aunque a menor nivel que las anteriores AUMENTO DE LA RELACIÓN DEL DIÁMETRO Crece Decrece Crece, aunque a menor nivel que las anteriores AUMENTO DEL NXP Despreciable Despreciable Despreciable Tabla 3. Evolución de la eficiencia de un eyector Otras investigaciones [43] han llegado a los mismos razonamientos incluyendo en sus análisis la boquilla secundaria y obteniendo que la succión aumenta a medida que se eleva la relación de presión y la eficiencia disminuye cuando se incrementa la relación del diámetro. También se menciona que la succión se puede ver alterada debido a la temperatura exterior del aire. Finalmente, se expone que el diseño geométrico es uno de los pasos más importantes para el funcionamiento correcto del eyector [46] y se intentan optimizar los parámetros geométricos del eyector mostrados en la Figura 10 Figura 10. Parámetros geométricos de un eyectorcomo la contrapresión y la relación de arrastre para así lograr un mayor rendimiento del eyector en condiciones de funcionamiento determinadas [34]. Un análisis de sensitividad concluye que X7 y X1 son los parámetros que más influyen en el rendimiento del eyector, el resto de las variables tienen un efecto casi imperceptible al lado de estas. Consecuentemente, la eficiencia del eyector se puede modificar de forma sustancial variando el valor de X7 y X1. Figura 10. Parámetros geométricos de un eyector 2.1.3 El fluido de trabajo Los eyectores son muy sensibles al fluido de trabajo en términos de operación y rendimiento [37]. Es por ello, que la selección del refrigerante apropiado es el paso más importante en el proceso de diseño de un sistema de refrigeración con eyectores. En el pasado, el criterio principal de selección era la maximización del rendimiento; con el paso del tiempo, otros factores como la seguridad o el coste también son considerados. A la hora de seleccionar el fluido de trabajo se deben tener en cuenta los siguientes requisitos: • Propiedades termofísicas del fluido favorables: El fluido debe tener un alto calor latente de Fundamentos y estado del arte 12 12 vaporización, un rango de temperatura del generador amplio, una alta temperatura crítica y la presión del fluido a la temperatura del generador no debe ser muy elevada. Los fluidos deberán tener una masa molecular alta, ya que se obtienen mayores eficiencias en el eyector; aunque no demasiado ya que, a mayor masa molecular, menor es el tamaño del eyector. También es recomendable que posean una viscosidad reducida y conductividad térmica alta. • Preservación del medio ambiente: El potencial de agotamiento de la capa de ozono (ODP) deberá ser nulo y el potencial de calentamiento global (GWP) muy bajo. • Cuestiones de seguridad del fluido: El fluido será quimicamente estable, no tóxico, no explosivo y no corrosivo. Se deberá consultar la norma ASHRAE (Norma de seguridad para sistemas de refrigeración y designación y clasificación de refrigerantes): ANSI/ASHRAE 15 Y 34-2019 [47]. • Otras características a tener en cuenta serán el coste del fluido y su disponibilidad en el mercado. La clasificación de los fluidos de trabajo se puede realizar de dos maneras: en función de la pendiente de la línea de vapor saturado en el diagrama T-s o a partir de su composición química molecular. Refiriéndose a la pendiente de la línea de vapor saturado en el diagrama T-s [48], los refrigerantes de agrupan entre fluidos secos, húmedos o isentrópicos (Figura 11). En el caso de los fluidos húmedos, pueden condensar y formarse gotas en el interior del eyector, lo que puede afectar seriamente al proceso dinámico del gas en el eyector y al rendimiento del mismo [49]. Para solucionar este problema, se sugiere recalentar el flujo primario para garantizar que la expansión en la boquilla tenga lugar en la región recalentada [50]. Para los fluidos secos e isentrópicos, en cambio, no se produce un cambio de fase durante el proceso de expansión en la boquilla y el recalentamiento no es necesario. Figura 11. Diagramas T-s de tres tipos de refrigerantes: (a) fluido húmedo, (b) fluido seco, y (c) fluido isentrópico Basándose en su composición química, los fluidos de trabajo se pueden clasificar en: • Halocarbonos: Que a su vez se clasifican en clorofluorocarbonos (CFC), hidroclorofluorocarbonos (HCFC), hidrofluorocarbonos (HFC) y hidrofluoroolefina (HFO). • Hidrocarburos (HC): Está constituido por compuestos formados por H y C como el metano, etano, propano, ciclopropano, butano, isobutano… • Refrigerantes compuestos que consisten en mezclas de halorcarbonos. • Otros refrigerantes como el agua, el aire, el amoniaco o el CO2. En la siguiente tabla se pueden ver las características de los principales refrigerantes (Tabla 4). La elección final del refrigerante viene dada por un compromiso entre la eficiencia en rangos de funcionamiento seleccionados y el impacto medioambiental. En este trabajo, el refrigerante elegido es el CO2. 13 TIPO NOMBRE REFRIGERANTE VAPOR WET/DRY MASA MOLECULAR (Kg/Kmol) PUNTO EBULL A 1atm (ºC) CALOR LATENTE (KJ/Kg) GWP ODP Halocarbonos CFC R11 WET 137.4 23.7 186.2 4750 1 R12 WET 120.9 -29.8 147.8 10900 1 R113 DRY 187.4 47.6 155.9 6130 0.85 R114 DRY 10.9 3.8 133.7 9180 0.58 HCFC R22 WET 86.5 -40.8 196.8 1790 0.05 R123 DRY 152.9 27.9 1775 77 0.01 R141B DRY 116.9 32.1 233.1 717 0.12 R142B DRY 100.5 -9.2 212 2220 0.06 R500 WET 99.3 -33.6 - 8100 0.61 R502 WET 111.6 -45.3 - 4600 0.31 HFC R134A WET 102 -26.1 190.9 1370 0 R152A WET 66.1 -24 295.8 133 0 R245FA DRY 134.1 15.1 199 1050 0 HFO R1234YF WET 114.1 -30 180.25 4 0 Hidrocarburos HC R290 WET 44.1 -42.1 360.3 20 0 R600 DRY 58.1 -0.5 376.1 20 0 R600A DRY 58.1 -11.8 344.6 20 0 Compuestos R407A WET 90.1 -45.3 235.57 2107 0 R407B WET 102.9 -46.8 275 2804 0 R410A DRY 72.6 -51.58 276 2088 0 Otros AMONIACO R717 WET 17 -33.3 1226.1 0 0 AGUA R718 WET 18 100 2477.2 0 0 CO2 R744 WET 44 -78.5 197.7 1 0 Tabla 4. Características de los principales refrigerantes. Fundamentos y estado del arte 14 14 2.1.4 El ciclo de refrigeración con eyector El eyector, como se ha explicado anteriormente, tiene la función de aumentar la presión de un fluido primario que sale del condensador a través de uno secundario que ha circulado por el evaporador y que produce el efecto frigorífico. En él, ambos fluidos se mezclan y el fluido mixto se devuelve al compresor. Existen diversas configuraciones de ciclos de refrigeración con eyector como se pueden observar en el esquema de la Figura 12. Figura 12. Diferentes configuraciones de ciclos de refrigeración con eyetores A los ERS se les pueden añadir modificaciones, derivando en los diferentes tipos de ciclo mencionados en la figura anterior. El ciclo de refrigeración con eyector más básico es el SERS, que contiene un único eyector y es el utilizado en la investigación realizada en este trabajo, sustituyendo un eyector fijo por uno de geomería variable (Figura 13). Figura 13. Ciclo de refrigeración con eyector 15 2.2 Identificación de sistemas El conocimiento exhaustivo del modelo dinámico de un proceso industrial permite comprender el comportamiento del sistema y sus partes críticas. Es necesario disponer de un modelo matemático preciso para alcanzar un estado de mejora continua e innovación, logrando así aumentar la eficiencia del proceso y la optimización del diagnóstico y detección de fallos. En el momento que no se tiene un modelo dinámico completamente definido, es requerido aplicar la identificación con el fin de construir un modelo matemático del sistema basado en las entradas y salidas del mismo. El concepto de modelo dinámico se consigue al relacionar matemáticamente las variables de salida y entrada del sistema. Los valores de las señales de salida dependerán tanto de los valores instantáneos de las señales de entrada como del comportamiento anterior del sistema. La descripción del modelo dinámico puede realizarse mediante ecuaciones diferenciales, funciones de transferencia, ecuaciones en el espacio de estados y modelos de polo-cero-ganancia; y en cuanto a la representación temporal, se recurre tanto al tiempo continuo como al discreto. Básicamente se pueden distinguir cuatro áreas problemáticas en la teoría del sistema [25]: modelado, análisis, estimación y control. • Modelado: Se debe encontrar un modelo matemático que describa adecuadamente la situación física. • Análisis: Se analiza el comportamiento de salida del sistema mediante simulación. • Estimación: Después de haber obtenido una estructura de modelo adecuada estable, identificable y observable, se estiman las variables desconocidas a partir de un conjunto de datos dado de variables de entrada-salida. • Control: Se definen dos tipos de estrategias de control, controles en bucle abierto y en bucle cerrado. En el control en bucle abierto, la salida no interfiere en el funcionamiento del sistema. En cambio, en el control en bucle cerrado, se produce realimentación, es decir, utiliza observaciones reales de la salida para el cálculo de la entrada. En consecuencia, la identificación del sistema se considera como un modelado aproximado para una aplicación específica sobre la base de los datos observados y del conocimiento previo del sistema. 2.2.1 El proceso de identificación A continuación, se describe el procedimiento de identificación como un proceso iterativo con varias fases, cuyo objetivo es alcanzar el modelo matemático adecuado (véase la Figura 14) en base a datos experimentales. El primer paso en el proceso de identificación es el diseño del experimento, donde se pone de manifiesto toda la dinámica del sistema teniendo en cuenta la entrada, perturbaciones y salida para luego realizar la simulación. Una vez simulado, se adquieren los datos experimentales sin manipular, que contienen imperfecciones que deben ser tratadas para conseguir una identificación más limpia. Seguidamente, en el tratamiento de datos se aplican filtros a los datos obtenidos tras la simulación para eliminar el ruido, perturbaciones, datos duplicados, discontinuidades o mediciones incorrectas. Una vez se han procesado los datos, entra en juego la parte más complicada del proceso, la elección de la estructura que depende de dos factores, el conocimiento previo del sistema y el tipo de datos que se quieren conseguir, teniendo esto en cuenta se puede enfocar la adquisición del modelo de diversas maneras que se explicarán más adelante. El siguiente paso viene dado por la selección del método de estimación de parámetros, donde se opta por un método de estimación y un criterio de bondad para expresar lo adecuadamente que se ajusta el modelo entre los datos observados y esperados, en este trabajo se utilizará el método de minimización del error. Finalmente, se procede a validar si el modelo identificado es o no lo suficientemente bueno para su uso previsto, en el caso de que no se considere apropiado, se realizará una nueva iteración y el proceso comenzará de nuevo en el paso que se estime necesario. Fundamentos y estado del arte 16 16 Figura 14. Proceso de identificación. 2.2.2 Elección de la estructura A la hora de construir un modelo, la elección de la estructura juega un papel crucial en el proceso de identificación. Para ello, es vital el conocimiento del proceso y del sistema a identificar [23]. El camino que se debe recorrer para definir la estructura viene determinado por tres pasos: 1. Elección del tipo del set de modelos: Selección entre modelos lineales y no lineales y modelos de caja blanca, negra o gris. 2. Elección del tamaño del set de modelos: Se plantea la selección del grado de los polinomios y los tipos de variables a incluir en el modelo. 3. Elección de la parametrización del modelo: Cuando un set de modelos ha sido seleccionado, todavía resta parametrizarlo, esto es, encontrar una estructura particular que se adecue al set de datos disponible y a la aplicación buscada. 2.2.3 Validación del modelo Una vez elegida la estructura, el proceso de identificación nos proporciona un modelo concreto para dicha configuración. Este modelo puede ser el mejor disponible, ¿pero es lo suficientemente bueno para el fin deseado? Comprobar si un modelo es el adecuado se conoce como validación del modelo. Para encontrar la estructura apropiada, se deberán probar diferentes modelos, comparar sus resultados y pasar un proceso de “control de calidad” antes de ser entregado el modelo final al usuario. Durante la validación se verifica la conveniencia del ajuste del modelo para ver si es o no suficiente y, en el caso de que el resultado no sea satisfactorio, el proceso de identificación vuelve a comenzar probando una nueva configuración. En el proceso de validación no importa cómo se ha estimado o construido el modelo, lo más importante es determinar su adecuación a las medidas iniciales sin haberlas reproducido y lo “lejos” que está el modelo de la 17 descripción real del mismo [51]. Entre las formas de evaluación existentes, se aplicarán comparaciones basadas en conjuntos de datos del sistema; no se pueden utilizar los mismos datos para la estimación y la validación. Los datos de estimación se componen de aquellos utilizados para determinar el modelo, mientras que la validación utilizará un conjunto de datos disponible que no se haya empleado para construir el modelo que se quiere evaluar (datos de validación). 2.2.4 Identificación en bucle abierto y cerrado Una vez comienza el proceso de identificación, es necesario reconocer ante qué sistema de control se encuentra el usuario y si este posee o no realimentación. Un sistema en bucle abierto es aquel en el que únicamente actúa el proceso sobre la señal de entrada y como resultado se obtiene una señal de salida independiente a la de entrada, aunque basada en ella (Figura 15). No existe retroalimentación hacia el controlador. Es la técnica de identificación utilizada por defecto a la hora de adquirir un modelo dinámico. Figura 15. Diagrama de bloques identificación en bucle abierto. En cambio, se está ante un sistema de control en bucle cerrado cuando la acción de control se encuentra en función de la señal de salida, es decir, existe realimentación como se puede ver en la Figura 16. En este caso, la señal de salida, además de ser salida del sistema, es utilizada como una de las entradas de él, aportándole información útil. La principal diferencia entre ambos (BA y BC) es el origen de los datos y no los métodos utilizados para adquirirlos. Figura 16. Diagrama de bloques identificación en bucle cerrado. Las primeras investigaciones en cuanto a identificación en bucle cerrado [24], [52]–[55] surgen cuando resulta imposible la identificación en BA por motivos de seguridad, económicos, o si el sistema contiene mecanismos de retroalimentación inherentes. Como puede ser observado en la Figura 2 perteneciente al artículo sobre el que se va a trabajar, el sistema dinámico que se necesita identificar se encuentra bajo realimentación, es decir, en bucle cerrado. Para la resolución de este problema, se recurre a varios enfoques: • Aproximación directa: La retroalimentación es ignorada y el sistema se identifica en BA, utilizando las mediciones obtenidas de la entrada y la salida. Este enfoque consite en aplicar directamente un método de predicción del error como si no existiera el feedback. Se aplica a sistemas con mecanismos de retroalimentación arbitrarios (desconocidos). • Aproximación indirecta: Se trata de identificar la función de transferencia en BC y determinar sus parámetros en BA valiéndose del conocimiento del controlador y en la señal de referencia. Se aplica a sistemas con retroalimentación no lineal. • Aproximación conjunta de entrada-salida: En este enfoque se consideran la entrada y la salida conjuntamente como la salida del sistema impulsado por una señal de entrada adicional. Se determinan los parámetros en BA a partir de una estimación de este sistema aumentado. En este trabajo, se utilizará el enfoque indirecto debido al conocimiento de datos de entrada y del controlador. Identificación en bucle cerrado mediante un algoritmo propio 24 24 Una vez determinadas las coordenadas de todos los puntos, se realiza una conversión de ellas, inicialmente en píxeles ya que se trata de una imagen, a los valores por defecto que aparecen en la Figura 2. La muestra de datos experimentales está compuesta por 281 puntos. 3.2 Acondicionamiento de datos Tras recopilar los puntos del proceso, es probable que el algoritmo de identificación no pueda utilizarlos, por lo que se debe realizar un filtrado de datos para suavizar posibles errores cometidos a la hora de la obtención de la muestra. Como podemos observar en la Figura 19.a, la muestra tiene varias discontinuidades debido a la extracción realizada a mano: • Los datos no están uniformemente muestreados en el tiempo. • Existe retroceso temporal entre algún que otro punto y el posterior. • Aparición de puntos duplicados. • Outliers ocasionales: Observaciones infrecuentes o puntos que no parecen seguir la distribución característica del resto de los datos. El procedimiento seguido consiste en la realización de un promediado de todos los puntos. Primero de todo, se recorre la muestra en busca de puntos que decrecen temporalmente y que por lo tanto son inaceptables, y se igualan al punto posterior. A continuación, se aplica un promediado tomando el punto en cuestión, el anterior y el siguiente. Como resultado, en la Figura 19.b se observa una suavización de la curva, sin cambios bruscos entre puntos. Figura 19. Curva de puntos experimentales (a) tomados a mano y (b) filtrados 3.3 Aproximación indirecta Analizando los datos obtenidos y la información disponible tanto en la Figura 19.b como en el artículo sobre el que se trabaja [17], se llega a la conclusión de que se debe llevar a cabo una identificación en bucle cerrado ya que se proporciona la estructura de realimentación. Los parámetros del controlador son los siguientes: 𝐾𝑝=−13 𝑚𝑚2·𝑀𝑃𝑎−1 𝑇𝑖=15 𝑠 𝑇𝑑=0.3 𝑠 La identificación en BC, como se ha explicado en el Capítulo 2, se puede aproximar de tres maneras: directa, indirecta y entrada-salida conjunta. Observando los datos obtenidos, se llega a la conclusión que el enfoque más adecuado será el indirecto debido a que únicamente están disponibles los datos correspondientes a la señal de 25 referencia externa r y la salida y. En primer lugar, se estimará el sistema en bucle cerrado a partir de estos datos para posteriormente inferir la función de transferencia del modelo (en bucle abierto) haciendo uso del controlador conocido. Se tiene un sistema de la forma, 𝑦(𝑡)=𝐺𝑜(𝑞)𝑢(𝑡)+𝑣(𝑡)=𝐺𝑜(𝑞)𝑢(𝑡)+𝐻𝑜(𝑞)𝑒(𝑡) (15) 𝑢(𝑡)=𝑟(𝑡)−𝐶(𝑞)𝑦(𝑡) (16) donde 𝐺𝑜 es la función de transferencia del modelo dinámico, 𝐻𝑜 modelará el ruido utilizando un filtro inversamente estable, e corresponde al ruido blanco. La entrada u se compone de la señal de referencia y de un controlador invariante en el tiempo representado por C. El modelo de la planta a identificar se describe como 𝑦(𝑡)=𝐺(𝑞,𝜃)𝑢(𝑡)+𝐻(𝑞,𝜃)𝑒(𝑡) (17) y 𝜃 son los parámetros a calcular. Se asume que el BC está bien definido y por tanto su ecuación resultante es: 𝑦(𝑡)=𝐺𝑜(𝑞)𝑆𝑜(𝑞)𝑟(𝑡)+𝑆𝑜(𝑞)𝑣(𝑡) (18) 𝑢(𝑡)=𝑆𝑜(𝑞)𝑟(𝑡)−𝐶(𝑞)𝑆𝑜(𝑞)𝑣(𝑡) (19) 𝑆𝑜 𝑖 y 𝑆𝑜 son las funciones de sensibilidad de entrada y salida, respectivamente, que miden qué tan sensible es una señal de perturbación. 𝑆𝑜(𝑞)=1 1+𝐺𝑜(𝑞)𝐶(𝑞) (20) Finalmente, el sistema en BC queda definido del siguiente modo reescribiendo la ecuación (15): 𝑦(𝑡)=𝐺𝐵𝐶(𝑞)𝑟(𝑡)+𝑣𝐵𝐶(𝑡)=𝐺𝑜(𝑞) 1+𝐺𝑜(𝑞)𝐶(𝑞)𝑟(𝑡)+ 1 1+𝐶(𝑞)𝐺𝑜(𝑞)𝑣(𝑡) (21) La aproximación indirecta estimará 𝐺𝐵𝐶 para luego adquirir la función de transferencia en bucle abierto 𝐺𝑚𝑜𝑑𝑒𝑙𝑜 a partir de 𝐺𝐵𝐶 =𝐺𝑚𝑜𝑑𝑒𝑙𝑜𝐶 1+𝐺𝑚𝑜𝑑𝑒𝑙𝑜𝐶 (22) Una ventaja del enfoque indirecto es que se le puede aplicar cualquier método de identificación a (21) para la obtención de 𝐺𝐵𝐶, ya que se aborda desde un planteamriento en BA. Se optará por la minimización del error de predicción, que consiste en la búsqueda de un modelo donde el error entre los datos de aprendizaje y los del propio modelo sea el mínimo posible. En el método de predicción del error, se calcularán los parámetros de 𝐺𝐵𝐶(𝑞,𝜃) de modo que el error resultante de la diferencia entre la salida simulada y la medida sea mínimo: 𝑒(𝑡)=𝑦(𝑡)−𝑦(𝑡)=𝑦(𝑡)−𝐺𝐵𝐶(𝑞)𝑟(𝑡) (23) Este procedimiento además permite parametrizaciones arbitrarias 𝐺𝐵𝐶(𝑞,𝜃), de modo que es natural que los parámetros se relacionen con las propiedades del sistema en BA, por lo tanto (22) quedaría de la siguiente manera: 𝐺𝐵𝐶(𝑞,𝜃) = 𝐺𝑚𝑜𝑑𝑒𝑙𝑜(𝑞,𝜃)𝐶(𝑞) 1+𝐺𝑚𝑜𝑑𝑒𝑙𝑜(𝑞,𝜃)𝐶(𝑞) (24) Identificación en bucle cerrado mediante un algoritmo propio 26 26 El objetivo final será alcanzar la función de transferencia de la planta (relación entrada-salida) realizando una integración numérica. Se toma la ecuación (22) y se obtiene el modelo de la planta a través del comando feedback en MatLab. 𝐺𝑚𝑜𝑑𝑒𝑙𝑜(𝑞,𝜃)= 𝐺𝐵𝐶(𝑞,𝜃)𝐶(𝑞) 1−𝐺𝐵𝐶(𝑞,𝜃)𝐶(𝑞) (25) C corresponde al controlador PID. 𝐶 =𝐾𝑐(1+𝑇𝑑𝑠+ 1 𝑇𝑖𝑠) (26) Se tratará de inferir 𝐺𝑚𝑜𝑑𝑒𝑙𝑜 a través de modelos de caja negra y más adelante de caja gris, tratando así de obtener los mejores resultados posibles. 3.4 Modelos de caja negra Es sabido que en este tipo de modelos la dinámica se define a partir de datos de entrada y salida conocidos. Primeramente, se obtendrá la función en BC del sistema controlado a partir de la salida y la señal de referencia minimizando el error, para más adelante inferir cómo deber ser la función de transferencia de la planta con los parámetros obtenidos. Se examinarán ahora los intentos de identificación en BC siguiendo un método de caja negra para deducir la función de transferencia del modelo de planta comparando diversas estructuras. 3.4.1 Sistema de primer orden En cuanto a la elección del orden del sistema, se desea que el resultado sea causal, es decir, que dependa exclusivamente de los valores presentes y/o pasados de la entrada; no depende de valores futuros al no poderse obtener una salida antes de aplicar la entrada. Con miras a la causalidad, la planta no debe tener más ceros que polos. Como se conoce que la identificación se aplica a un sistema controlado por un PID y que este añade dos ceros, se podría concluir que la función de tranferencia contendría al menos dos ceros y como mínimo dos polos en su estructura y el sistema inicial debería ser de segundo orden. Sin embargo, es posible aplicar un sistema de primer orden dado que se produce cancelación de polos. Este fenómeno consiste en la cancelación de la dinámica de la planta con un cero del controlador produciéndose una cancelación polo-cero; uno de los polos de la función de transferencia identificada es igual al cero del controlador, por lo que se anulan, permitiendo así utilizar una estructura de primer orden para la identificación de la planta. Un sistema de primer orden viene dado por la siguiente ecuación: 𝐺𝐵𝐶 =𝐾 1+𝜏𝑠 (27) K corresponde a la ganancia estática y τ es la constante de tiempo del sistema. A continuación, se mostrarán los resultados obtenidos a través de de los ensayos que se han realizado a un sistema de primer orden, analizando también la validación del modelo. 3.4.1.1 Aproximación en un único escalón En referencia a la Figura 2, se aprecia la presencia de cuatro escalones correspondientes a distintos puntos de funcionamiento de la planta. Por tal motivo, en este primer planteamiento se evaluará cada escalón por separado y se sacará un modelo válido para todo el sistema centrándose únicamente en el intervalo de puntos seleccionado. A la hora de la valoración de la bondad del modelo, se evaluará 𝐺𝐵𝐶 en todo el ensayo. 27 Los datos experimentales no están uniformemente distribuidos en el tiempo, por lo que no se puede resolver la identificación a través del método de los mínimos cuadrados. Para analizar cada escalón por separado, se ha creado una función llamada mystep. Esta función transforma el modelo a identificar al espacio de estados, para más adelante, con el comando ode45 que resuelve ecuaciones diferenciales, minimizar el error y así inferir la función de transferencia. Escalón 1: El intervalo de puntos que conforman los datos de aprendizaje de este primer escalón es (1:48). La función de transferencia en BC a inferir es (27) y los parámetros a identificar serán K y τ. Par óptimo obtenido: 𝑃𝑎𝑟1=[𝐾,𝜏]=[1.0877,3.0542] Error de 𝑃𝑎𝑟1 respecto al escalón 1: 𝐸𝐶𝑀𝑒𝑠𝑐1 =0.3996 Función de transferencia en BC: 𝐺𝐵𝐶1=𝐾 1+𝜏𝑠 =1.0877 3.054𝑠+1 Una vez obtenida la funcion de transferencia del modelo controlado, se aprecia en la Figura 20.a que el orden a inferir deberá ser mayor, ya que la salida estimada no oscila como la deseada. Aunque el 𝐸𝐶𝑀𝑒𝑠𝑐1 sea pequeño, es insuficiente. Función de transferencia para el modelo completo evaluando el escalón 1: 𝐺𝑚𝑜𝑑𝑒𝑙𝑜1 =𝐺𝐵𝐶1𝐶 1−𝐺𝐵𝐶1𝐶=−1.0877𝑠 0.1191𝑠3+0.3936𝑠2+0.01508𝑠−0.0007596 Figura 20. Modelo de caja negra y sistema de primer orden. Resultados obtenidos en el primer escalón en (a) sistema controlado y (b) respuesta ante escalón unitario La Figura 20.b muestra la respuesta del sistema dinámico frente a una entrada escalón unitario. Analizando la gráfica, se observa que la respuesta transitoria no decae, por lo que se llega a la conclusión de que el sistema es inestable y por lo tanto no válido para la planta. Al ser inestable el sistema, el poceso de validación no es realizado. Escalón 2: El intervalo de puntos correspondiente a los datos de aprendizaje de este segundo escalón es (49:107). Par óptimo obtenido para el escalón 2: 𝑃𝑎𝑟2=[1.1059,2.9696] Error de 𝑃𝑎𝑟2 con respecto al escalón 2: Identificación en bucle cerrado mediante un algoritmo propio 28 28 𝐸𝐶𝑀𝑒𝑠𝑐1 =0.7745 Función de transferencia en BC: 𝐺𝐵𝐶2=1.1059 2.9696𝑠+1 De nuevo, se confirma en la Figura 21.a que la salida estimada no oscila como la deseada, por lo que se debería aumentar el orden del sistema. Función de transferencia para el modelo completo evaluando el escalón 2: 𝐺𝑚𝑜𝑑𝑒𝑙𝑜2 =−1.106𝑠 0.1148𝑠3+0.3819𝑠2+0.01197𝑠−0.0009178 Figura 21. Modelo de caja negra y sistema de primer orden. Resultados obtenidos en el segundo escalón en (a) sistema controlado y (b) respuesta ante escalón unitario Se muestra en la Figura 21.b que el sistema vuelve a ser inestable, por lo que tampoco se validará. Escalón 3: El intervalo de puntos correspondiente a los datos de aprendizaje del tercer escalón es (108:177). Par óptimo obtenido para el escalón 3: 𝑃𝑎𝑟3=[1.0851,2.2065] Error de 𝑃𝑎𝑟3 con respecto al escalón 3: 𝐸𝐶𝑀𝑒𝑠𝑐3 =1.0657 Función de transferencia en BC: 𝐺𝐵𝐶3=1.0851 2.2065𝑠+1 La Figura 22.a muestra que la salida no oscila, aunque el error es pequeño. Para obtener mayor precisión sería conveniente aumentar el orden del sistema, como se ha deducido en los escalones anteriores. Función de transferencia para el modelo completo evaluando el escalón 3: 𝐺𝑚𝑜𝑑𝑒𝑙𝑜3 =−1.0851𝑠 0.08605𝑠3+0.2835𝑠2+0.008064𝑠+0.0007373 29 Figura 22. Modelo de caja negra y sistema de primer orden. Resultados obtenidos en el tercer escalón en (a) sistema controlado y (b) respuesta ante escalón unitario La Figura 22.b enseña una respuesta inestable del sistema frente a una entrada esalon unitario. Escalón 4: El intervalo de puntos correspondiente a los datos de aprendizaje de este cuarto escalón es (178:281). Par óptimo obtenido: 𝑃𝑎𝑟4=[1.0572,1.3181] Error de 𝑃𝑎𝑟4 con respecto al escalón 4: 𝐸𝐶𝑀𝑒𝑠𝑐4 =1.632 Función de transferencia en BC: 𝐺𝐵𝐶4=1.0572 1.3181𝑠+1 La Figura 23.a muestra que la salida continúa sin oscilar. Función de transferencia para el modelo completo evaluando el escalón 4: 𝐺𝑚𝑜𝑑𝑒𝑙𝑜4 =−1.0572𝑠 0.05141𝑠3+0.1691𝑠2+0.003992𝑠−0.0004954 Figura 23. Modelo de caja negra y sistema de primer orden. Resultados obtenidos en el cuarto escalón en (a) sistema controlado y (b) respuesta ante escalón unitario Identificación en bucle cerrado mediante un algoritmo propio 30 30 Se muestra de nuevo en la Figura 23.b la respuesta del sistema; analizándola, se deduce que el sistema es inestable. 3.4.1.2 Aproximación del modelo completo Como se ha deducido del apartado anterior, el análisis de los escalones de manera aislada acumula un error considerable que hace el modelo poco conveniente. Se plantea entonces el problema de buscar un modelo dinámico apto que se centre en el sistema completo en vez de en intervalos pertenecientes a un escalón aislado. Otra función capaz de evaluar los datos no uniformemente muestreados en el tiempo ha sido creada para inferir la función de transferencia en el modelo completo. En este caso, mySim evaluará todos los escalones en vez de uno determinado al tener como argumento de entrada un vector que serpara los datos en diferentes escalones. Su funcionamiento será similar al de la función utilizada para cada escalón. Par óptimo obtenido: 𝑃𝑎𝑟=[1.0587,1.0863] Error del par con respecto al sistema: 𝐸𝐶𝑀 =4.0563 Función de transferencia en BC: 𝐺𝐵𝐶 =1.0587 1.0863𝑠+1 La Figura 24.a muestra una vez más que la salida aunque se ajuste a los saltos entre escalones es lineal, no oscila. Función de transferencia para el modelo completo: 𝐺𝑚𝑜𝑑𝑒𝑙𝑜 =−1.059𝑠 0.04237𝑠3+0.1389𝑠2+0.001784−0.0005087 La respuesta del sistema dinámico frente a una entrada escalón unitario es representada en la Figura 24.b. Analizándola, se interpreta que nuevamente el sistema resulta inestable y por lo tanto no válido para la planta. Figura 24. Modelo de caja negra y sistema de primer orden. Resultados obtenidos en sistema completo en (a) sistema controlado y (b) respuesta ante escalón unitario Como se ha observado en los sistemas de primer orden, todos ellos resultan ser inestables ante una respuesta escalón unitario. Tras no ser aceptable el modelo propuesto, se propone un aumento del orden del sistema, que pasará a ser de grado 2. 31 3.4.2 Sistema de segundo orden Es necesario aumentar el orden del sistema con el fin de adquirir un modelo dinámico más exacto y estable. Inicialmente en el sistema de primer orden, al variar los parámetros, únicamente se modificaba la velocidad de respuesta del sistema; en un sistema de segundo orden, la alteración de los parámetros cambia la forma total de respuesta. La mayoría de los sistemas industriales se comportan como uno de segundo orden, al exhibir una amplia gama de respuestas analizables. La función de transferencia para un sistema de segundo orden está diseñada en función de la ganancia K, la frecuencia natural 𝑤𝑛 y el coeficiente de amortiguamiento 𝜁. 𝐺𝐵𝐶 =𝐾𝑤𝑛 2 𝑠2+ 2 𝜁𝑤𝑛𝑠+𝑤𝑛 2 (28) Los parámetros de (28) están ligados al comportamiento físico de la respuesta y a la situación de los polos en el plano “S”. Únicamente conociendo el valor de 𝜁, podemos determinar el tipo de sistema y la forma de su respuesta: • Subamortiguado: (0<𝜁 <1) • Críticamente amortiguado: (𝜁 =1) • Sobreamortiguado: (𝜁 <1) • Oscilatorio: (𝜁 =0) • Inestable: (𝜁 <0) Se procederán a realizar los mismos experimentos que en el punto anterior, modificando el orden del sistema. 3.4.2.1 Aproximación en un único escalón Al presentarse cuatro escalones diferenciados, se planteará la evaluación del sistema mediante intervalos que comprendan un escalón completo, adaptando más adelante el modelo a la planta completa. Escalón 1: El intervalo de puntos correspondiente a los datos de aprendizaje de este primer escalón es (1:48). La función de transferencia en BC a inferir es (28) y los parámetros a identificar serán K, 𝜁 y 𝑤𝑛. Se utiliza el mismo algoritmo desarrollado para el sistema de primer orden, cambiando exclusivamente el número de parámetros y la estructura de la función de transferencia en BC a calcular. Par óptimo obtenido: 𝑃𝑎𝑟1=[𝐾,𝜁 ,𝑤𝑛]=[1.0292,0.3125,0.4020] Error de 𝑃𝑎𝑟1 respecto al escalón 1: 𝐸𝐶𝑀𝑒𝑠𝑐1 =0.0158 Función de transferencia en BC: 𝐺𝐵𝐶1=𝐾𝑤𝑛 2 𝑠2+ 2 𝜁𝑤𝑛𝑠+𝑤𝑛 2=0.1663 𝑠2+0.2513𝑠+0.1616 Una vez obtenida la funcion de transferencia del modelo controlado, se aprecia en la Figura 25.a que el orden podría ser adecuado, ya que la salida estimada oscila y se asemeja a la deseada. Además, el error se ve reducido considerablemente frente a la estructura de primer orden. Función de transferencia para el modelo completo evaluando el escalón 1: 𝐺𝑚𝑜𝑑𝑒𝑙𝑜1 =𝐺𝐵𝐶1𝐶 1−𝐺𝐵𝐶1𝐶=−0.1663𝑠 0.039𝑠4+0.1398𝑠3+0.04115𝑠2+0.001564𝑠−4.092·10−5 Identificación en bucle cerrado mediante un algoritmo propio 32 32 Figura 25. Modelo de caja negra y sistema de segundo orden. Resultados obtenidos en el primer escalón en (a) sistema controlado y (b) respuesta ante escalón unitario La Figura 25.b muestra la respuesta del sistema dinámico frente a una entrada escalón unitario. Analizando la gráfica, se observa que la respuesta transitoria no decae, por lo que se llega a la conclusión de que el sistema es inestable y por lo tanto no válido para la planta. La inestabilidad puede venir dada por el hecho de que se están forzando los polos en BC. De nuevo, al resultar inestable, no se validará el modelo resultante. Escalón 2: El intervalo de puntos correspondiente a los datos de aprendizaje del segundo escalón es (49:107). Par óptimo obtenido para el escalón 2: 𝑃𝑎𝑟2=[1.0143,0.2735,0.3989] Error de 𝑃𝑎𝑟2 respecto al escalón 2: 𝐸𝐶𝑀𝑒𝑠𝑐2 =0.0622 Función de transferencia en BC: 𝐺𝐵𝐶2=0.1614 𝑠2+0.2182𝑠+0.1591 Figura 26. Modelo de caja negra y sistema de segundo orden. Resultados obtenidos en el segundo escalón en (a) sistema controlado 33 y (b) respuesta ante escalón unitario La salida en la Figura 26.a se asemeja a la deseada, ajustándose poco a poco las oscilaciones. Función de transferencia para el modelo completo evaluando el escalón 2: 𝐺𝑚𝑜𝑑𝑒𝑙𝑜2 =−0.1614𝑠 0.039𝑠4+0.1385𝑠3+0.03694𝑠2+0.001595𝑠−1.972·10−5 Se observa en la Figura 26.b que el sistema es inestable. Escalón 3: El intervalo de puntos correspondiente a los datos de aprendizaje del segundo escalón es (108:177). Par óptimo obtenido para el escalón 3: 𝑃𝑎𝑟3=[0.9869,0.2105,0.4911] Error de 𝑃𝑎𝑟3 respecto al escalón 3: 𝐸𝐶𝑀𝑒𝑠𝑐3 =0.0641 Función de transferencia en BC: 𝐺𝐵𝐶3=0.238 𝑠2+0.2067𝑠+0.2411 Como se ilustra en la Figura 27.a, la salida estimada es muy similar a la que se quiere obtener. Función de transferencia para el modelo completo evaluando el escalón 3: 𝐺𝑚𝑜𝑑𝑒𝑙𝑜3 =−0.238𝑠 0.039𝑠4+0.1381𝑠3+0.03566𝑠2+0.002202𝑠+2.738·10−5 Figura 27. Modelo de caja negra y sistema de segundo orden. Resultados obtenidos en el tercer escalón en (a) sistema controlado y (b) respuesta ante escalón unitario Como se puede apreciar en la Figura 27.b, el sistema es estable, subarmortiguado y podría ser válido. A continuación, se realizará la valoración de la bondad del modelo, donde se evaluará 𝐺𝐵𝐶3 en todo el ensayo y se obtendrá el error para la planta completa. Para ello, se hará uso de la función mySim. Error de 𝑃𝑎𝑟3 con respecto a la planta: 𝐸𝐶𝑀3=4.0602 La Figura 28 revela que la aplicación del par obtenido para el escalón 3 aplicado al sistema completo no se adecua al resultado esperado. Las oscilaciones no se ajustan a cada escalón, únicamente al tercero, cosa que podría resultar esperada ya que se ha utilizado el modelo de ese escalón para todo el sistema. Identificación en bucle cerrado mediante un algoritmo propio 40 40 Figura 37. Modelo de caja negra y sistema de orden configurable. Aplicación del par obtenido en el tercer escalón al sistema completo Escalón 4: El intervalo de puntos correspondiente a los datos de aprendizaje del cuarto escalón es (178:281) y el par óptimo obtenido es el siguiente: 𝑃𝑎𝑟4=[0.4380,0.4044,0.7568,2.2440,0.6399] Error de 𝑃𝑎𝑟4 respecto al escalón 4: 𝐸𝐶𝑀𝑒𝑠𝑐4 =0.1193 Función de transferencia en BC: 𝐺𝐵𝐶4=0.0438𝑠2+0.4044𝑠+0.7568 𝑠3+2.244𝑠2+0.6399𝑠+0.7568 En este último escalón, el error es mayor que en los anteriores. La Figura 38.a demuestra que, de todos modos, la estimación se ajusta aceptablemente a lo esperado. Función de transferencia para el modelo completo evaluando el escalón 4: 𝐺𝑚𝑜𝑑𝑒𝑙𝑜4 =−0.438𝑠3−0.4044𝑠2−0.7568𝑠 0.039𝑠5+0.2004𝑠4+0.2526𝑠3+0.04627𝑠2+0.002041𝑠 Figura 38. Modelo de caja negra y sistema de orden configurable. Resultados obtenidos en el cuarto escalón en el (a) sistema 41 controlado y (b) respuesta ante escalón unitario La Figura 38.b muestra la respuesta del sistema estable ante una entrada escalón unitario y para la validación se vuelve a evaluar 𝐺𝐵𝐶4 en la planta completa adquiriendo el siguiente error en el 𝑃𝑎𝑟4 con respecto a la planta: 𝐸𝐶𝑀4=3.7429 El ECM en cada escalón se mantiene dentro de unos valores similares. Como se puede ver en la Figura 39, el sistema controlado 𝐺𝐵𝐶4 se ajusta casi perfectamente al último escalón, acrecentándose el error en los demás escalones. Se aprecia que en los escalones 1, 2 y 3 hay presencia de sobreoscilaciones. A continuación, se modificará la estructura, tratando de inferir el modelo para el sistema completo en vez de para un único escalón. Figura 39. Modelo de caja negra y sistema de orden configurable. Aplicación del par obtenido en el cuarto escalón al sistema completo 3.4.3.2 Aproximación del modelo completo El proceso de identificación tomando únicamente escalones aislados para luego ampliar el estudio al sistema completo no ha proporcionado resultados adecuados. Se realiza un nuevo planteamiento que busca hallar un modelo válido para el conjunto de toda la planta. Par óptimo obtenido: 𝑃𝑎𝑟=[1.3953,1.3062,2.2823,7.9803,2.4311] Error del par con respecto al sistema: 𝐸𝐶𝑀 =2.3393 Función de transferencia en BC: 𝐺𝐵𝐶 =1.3953𝑠2+1.3062𝑠+2.2823 𝑠3+7.9803𝑠2+2.431𝑠+2.2823 Se obtiene una respuesta controlada del modelo razonable, aunque la oscilación en todos los escalones es la misma en vez de variar como se puede ver en la Figura 40.a. La diferencia del tiempo de establecimiento ante cada cambio de escalon puede ser debido a que realmente no se está utilizando un mismo modelo para todo el sistema, sino uno especifico para cada punto de funcionamiento. Función de transferencia para el modelo completo: 𝐺𝑚𝑜𝑑𝑒𝑙𝑜 =−1.3953𝑠3−1.3062𝑠2−2.2812𝑠 0.039𝑠5+0.3868𝑠4+0.9086𝑠3+0.2033𝑠2+0.009749𝑠 Identificación en bucle cerrado mediante un algoritmo propio 42 42 Figura 40. Modelo de caja negra y sistema de orden configurable. Resultados obtenidos en el sistema completo en (a) sistema controlado y (b) respuesta ante escalón unitario En esta ocasión, al aplicar una entrada escalón unitario a la estructura de orden superior al sistema completo, se consigue la estabilidad (Figura 40.b) que se incumplía en modelos de grado inferior. Igualmente, el ajuste en cada escalón no es exacto y se concluye con que este método de identificación no es el mas idóneo ya que se intuye la no linealidad presente en el modelo. Se pasará a buscar una función de transferencia que ajuste con más exactitud la salida combinando los modelos de cada escalón. 3.4.3.3 Aproximación multimodelo En las simulaciones realizadas escalón a escalón, a medida que varía la referencia, la respuesta en BC también lo hace. Resulta razonable afirmar que el comportamiento del modelo es no lineal, es decir, que cambia en función del punto de funcionamiento en que se encuentre el sistema. Se ha comprobado que el modelado de sistemas es una tarea compleja al adoptarse modelos lineales aproximados que no representan adecuadamente el sistema real y sus no linealidades. Es por ello, que se ha buscado otro enfoque para la identificación del sistema no lineal. Inicialmente, se propuso realizar una variación de la ganancia para ver si se podrían solventar las desigualdades observadas en los diferentes escalones del sistema identificado, llegando a la conclusión de que no justifica la diferencia entre las sobreoscilaciones y la rapidez de las mismas. Al no explicar la ganancia variable este fenómeno, se ha deducido que el sistema está sometido a una dinámica cambiante, en otras palabras, existe un modelo dinámico distinto en cada punto de funcionamiento, como resultado, se estaría frente a un multimodelo. Existe una amplia literatura en cuanto a sistemas multimodelo se refiere [58], [77]–[82]. Diversos autores han profundizado en esta técnica de identificación, que resulta favorable en cuanto al modelado de sistemas no lineales. Un sistema multimodelo está basado en la descomposición del problema global en un conjunto de modelos locales que representarán la dinámica de una región específica: un modelo para cada escalón. Cada submodelo describirá el comportamiento del sistema en las diferentes áreas de funcionamiento y el conjunto de todas ellas se aproximará de una manera más exacta al sistema real. De hecho, la solución de un sistema múltiple combinado se compone de un número de sistemas dinámicos sencillos que se fusionan [77], donde cada modelo local reproduce la dinámica del sistema en una región especifica del espacio. La mayoría de los procesos industriales son complejos y no lineales; la ventaja de utilizar un modelo múltiple es la posibilidad de que éste trabaje a lo largo de una trayectoria de funcionamiento fija que se compone de varios puntos de funcionamiento predeterminados. Cuando el proceso está funcionando en un determinado punto de operación, se utiliza el submodelo en cuyo rango está contenido el punto de funcionamiento en cuestión; cuando el punto de funcionamiento cambia de rango, automáticamente se pasa a utilizar el modelo local que lo incluya en su región. En definitiva, la descripción del sistema global será la integración de la familia de modelos locales que lo componen. Anteriormente no ha sido posible la utilización del multimodelo porque en los sistemas de primer y segundo orden existían tramos donde su respuesta transitoria resultaba ser inestable. En esta ocasión, todos los tramos 43 identificados cumplen con la condición de estabilidad y se procederá a valorar la eficiencia de la estructura. Los cuatro modelos locales obtenidos en cada escalón se fusionan y se construye un sistema no lineal, donde en cada situación, un modelo diferente es utilizado: 𝐺𝑚𝑜𝑑𝑒𝑙𝑜1 =0.04156𝑠3−0.21586𝑠2−0.04597 0.039𝑠5+0.1574𝑠4+0.1013𝑠3+0.0107𝑠2+0.0003081𝑠,𝑟𝑒𝑓=−0.5 𝐺𝑚𝑜𝑑𝑒𝑙𝑜2 =0.131𝑠3+0.1481𝑠2−0.6467𝑠 0.039𝑠5+0.2778𝑠4+0.544𝑠3+0.1754𝑠2+0.009502𝑠,𝑟𝑒𝑓 =0 𝐺𝑚𝑜𝑑𝑒𝑙𝑜3 =0.01738𝑠3+0.4167𝑠2−1.743𝑠 0.039𝑠5+0.3876𝑠4+0.9468𝑠3+0.3217𝑠2+0.01763𝑠,𝑟𝑒𝑓 =0.5 𝐺𝑚𝑜𝑑𝑒𝑙𝑜4 =−0.438𝑠3−0.4044𝑠2−0.7568𝑠 0.039𝑠5+0.2004𝑠4+0.2526𝑠3+0.04627𝑠2+0.002041𝑠,𝑟𝑒𝑓 =1 El ECM medio de todos los modelos es el siguiente: 𝐸𝐶𝑀𝑚𝑒𝑑𝑖𝑜 = 0.0476 Como se demuestra en la Figura 41, todos los escalones se ajustan de manera correcta ya que se ha utilizado un sistema local para cada uno de ellos. El 𝐸𝐶𝑀𝑚𝑒𝑑𝑖𝑜 es muy pequeño y la respuesta ante una entrada escalón unitario en cada tramo es estable. Figura 41. Modelo de caja negra y sistema de orden configurable. Resultados obtenidos en el sistema multimodelo En el caso de que la referencia no sean los valores exactos a los de la gráfica, se realizará una ponderación entre un modelo y otro. Por ejemplo, imagine que se tiene una referencia de 0.2, este modelo se encontraría entre los escalones 2 y 3. El modelo dinámico asociado a esta referencia sería el siguiente: 𝐺𝑚𝑜𝑑𝑒𝑙𝑜0.2 =0.6∙𝐺𝑚𝑜𝑑𝑒𝑙𝑜2 +0.4∙𝐺𝑚𝑜𝑑𝑒𝑙𝑜3 Exceptuando este sistema multimodelo, los modelos de caja negra no han proporcionado resultados óptimos. Los errores en las estructuas de órdenes pequeños eran grandes y los sistemas resultaban inestables. La estructura de orden configurable ha sido la que ha arrojado los menores ECMs y la mayor estabilidad, aunque el tiempo de computación del algoritmo resultaba mayor que en los otros sistemas. 3.5 Modelos de caja gris Al identificar a través de modelos de caja negra, se han obtenido resultados inestables y con errores elevados que no se ajustaban a la curva deseada, estas disconformidades podrían aparecer al inferir la función de transferencia controlada, forzando sus polos en BC. En consecuencia, se ha optado por repetir los ensayos anteriores utilizando un modelo de caja gris, donde se identifica directamente la función de transferencia de la Identificación en bucle cerrado mediante un algoritmo propio 44 44 planta a partir de datos de entrada y salida conocidos. A continuación, se expondrán los resultados obtenidos al realizar diversos experimentos aplicando un modelo de caja gris al proceso de identificación. 3.5.1 Sistemas de primer y segundo orden Como se ha comentado anteriormente, en los modelos de caja negra la función de transferencia del modelo se obtenía algebraicamente una vez encontrado un par óptimo para la función de transferencia controlada. En cambio, en los modelos de caja gris, la realimentación y el controlador son introducidos en la función que minimiza el ECM de la función de transferencia en BA, es decir, la identificación en este caso devuelve los parámetros del modelo de planta; no se deben realizar operaciones a posteriori para hallar la función de transferencia del modelo. En los modelos de caja gris no se podrán reproducir las estructuras de primer y segundo orden utilizadas en los de caja negra. Esto se debe a que, al estar infiriendo la función realimentada, los ceros y polos de un sistema de primer o segundo orden no son suficientes para constituir un sistema válido. Para probar esta afimación, se han introducido estructuras de primer y segundo orden en el algoritmo y los resultados obtenidos se analizarán seguidamente. En un sistema de primer orden, los valores óptimos de la función de transferencia del modelo obtenidos a través del primer escalon presentan una magnitud demasiado elevada. 𝑃𝑎𝑟1orden =[1.3568,0.5105]·1016 Como muestra la Figura 42.b, ocurre lo mismo con la duración del régimen transitorio ante una entrada escalón unitario. Figura 42. Modelo de caja gris y sistema de primer orden. Resultados obtenidos en el primer escalón en (a) sistema controlado y (b) respuesta ante escalón unitario Con respecto al sistema de segundo orden, el tiempo de computación del algoritmo creado es indefinido, ni siquiera es capaz de devolver ningún valor y, por lo tanto, se prueba que tampoco es una estructura válida. Como consecuencia, surge la duda sobre qué estructura se debe utilizar para inferir el modelo de la planta. Para encontrar un modelo adecuado, se ha decidido utilizar el método del Lugar de las Raíces para determinar de manera manual estructuras con un número de ceros y polos determinados que proporcionen buenos resultados. 3.5.2 El Lugar de las Raíces La respuesta dinámica de un sistema de control en BC depende de la posición de los polos de su función de transferencia, que a su vez dependen de la posición de los polos y ceros en BA. El Lugar de las Raíces es un método gráfico que permite conocer la evolución dinámica de un sistema cuando se produce una modificación 45 en el valor de uno o varios parámetros y ajustar el valor de dicho parámetro para obtener el comportamiento deseado del sistema. Para la construcción del Lugar de las Raíces, se deben marcar en el plano la situación de los polos y ceros de la función de transferencia. A partir de aquí, una serie de reglas son aplicadas. Las reglas están basadas en la relación existente entre los polos y los ceros de la funcion de transferencia y las raíces de la ecuación característica y son las siguientes: 1. Número de ramas: El número de ramas es igual al número de polos en BA. Cada rama comienza en un polo y termina en un cero. Las ramas que no pueden llegar a ningún cero terminan en infinito. 2. Lugar de las Raíces sobre el eje real del plano: Un punto perteneciente al eje real pertenecerá al Lugar de las Raíces si la suma del número de polos y ceros situados a su derecha es impar. 3. Simetría: El Lugar de las Raíces debe ser simétrico respecto al eje real; es decir, cada polo complejo tiene su conjugado. 4. Asíntotas e intersección de las mismas: Cuando el número de polos en BA es mayor que el número de ceros, las ramas que no terminan en ningún cero tienden asintóticamente al infinito. Todas las asíntotas cortarán al eje real en el centroide. 5. Ángulos de salida y llegada de las ramas: Son los ángulos con los que parte la rama desde cada polo o con los que llegan a cada cero. 6. Puntos de dispersión o confluencia de ramas: Si una serie de ramas convergen a un punto, de ese mismo punto deben salir tantas ramas como entraron. Como el lugar de las raíces es simétrico, los puntos de dispersión y confluencia deben aparecer en el eje real. 7. Intersección con el eje imaginario: Los puntos de corte con el eje imaginario corresponden con el valor de ganancia K que hacen el sistema en BA críticamente estable. Esto se obtiene mediante el criterio de estabilidad de Routh. 8. Suma de raíces: Si la ecuación característica se deja en forma de polonomio mónico, la suma de las raíces es constante e igual al coeficiente del término 𝑠𝑛−1, cambiado de signo. 9. Determinación de la ganancia K : En todos aquellos casos en los que sea necesaria la obtención del valor de la ganancia K para la cual algún punto del plano pertenece al Lugar de las Raíces se realizará aplicando la condición del módulo. Figura 43. (a) Representación del Lugar de las Raíces. (b) Respuesta ante escalón unitario En este trabajo, el método del lugar de las raíces será aplicado a través del Control System Toolbox de MatLab. Esta herramienta proporciona algoritmos y aplicaciones para analizar, diseñar y ajustar sistemas de control Identificación en bucle cerrado mediante un algoritmo propio 46 46 lineales de forma metódica. Se ha utilizado la herramienta de diseño SISO configurada con el método del Lugar de las Raíces para modelar de manera interactiva una estructura potencialmente estable a partir de unas condiciones de diseño. El cumplimento de las especificaciones se hará en base al trazado del Lugar de las Raíces y las regiones definidas por las condiciones de diseño, desplazando los polos en BC a través del ajuste de la ganancia K de la ecuación característica. Estas condiciones se han obtenido a partir de la Figura 2 y son: • Sobreocilación <40%, que determina el ángulo de las rectas que parten del orígen. Esta condición entrega el conjunto de puntos donde se pueden ubicar los polos en BC ajustando correctamente la ganancia del sistema. • Frecuencia natural >𝜋 5 ⁄≈0.6 , que es el radio del semicírculo cuyo centro se posiciona en el origen. Los ceros y polos podrán ubicarse dentro de este semicírculo. Las líneas negras de la Figura 43.a corresponden a las restricciones de diseño, los círculos y las cruces a los ceros y los polos respectivamente, y, las líneas azules al Lugar de las Raíces. En la Figura 43.b se muestra la gráfica de respuesta del sistema a la entrada escalón unitario. 3.5.3 Sistema de 1 cero y 2 polos Añadiendo 1 cero y 2 polos el sistema parece resultar estable, por lo que se repetirán los experimentos con esta estructura. 𝐺𝑚𝑜𝑑𝑒𝑙𝑜 =𝑝𝑎𝑟(1)∗(𝑠+𝑝𝑎𝑟(2)) (𝑠+𝑝𝑎𝑟(3))∗(𝑠+𝑝𝑎𝑟(4)) (30) 3.5.3.1 Aproximación en un único escalón Se evaluará cada escalón por separado para finalizar con la validación del modelo obtenido en el ensayo completo. En el caso de que la respuesta del sistema ante una entrada en escalón sea inestable, dicho proceso de validación será omitido y se pasará a evaluar el siguiente escalón. Escalón 1: El intervalo de puntos correspondiente a los datos de aprendizaje de este primer escalón es (1:48). La función de transferencia del modelo será directamente estimada. Par óptimo obtenido: 𝑃𝑎𝑟1=[0.0137,86.6505,0.2642,0.0160] Error de 𝑃𝑎𝑟1 respecto al escalón 1: 𝐸𝐶𝑀𝑒𝑠𝑐1 =0.008 Función de transferencia del modelo evaluando el primer escalón: 𝐺𝑚𝑜𝑑𝑒𝑙𝑜1 =0.0137𝑠+1.187 𝑠2+0.2802+0.004227 El ECM obtenido en este primer escalón es muy pequeño. Evaluando la Figura 44, se deduce que el par podría ser una buena opción ya que se ajusta muy adecuadamente a la curva y el tiempo hasta alcanzar el régimen permanente ante una respuesta escalón unitario no es demasiado elevado. Debido a esto, se procederá a validar este modelo en el sistema completo. En la valoración de la bondad del modelo se evaluará 𝐺𝑚𝑜𝑑𝑒𝑙𝑜1 en todo el ensayo y se obtendrá el error para la planta completa. Error de 𝑃𝑎𝑟1 con respecto a la planta: 𝐸𝐶𝑀1=6.6598 47 Figura 44. Modelo de caja gris y sistema de 1 cero y 2 polos. Resultados obtenidos en el primer escalón en (a) sistema controlado y (b) respuesta ante escalón unitario La Figura 45 revela que la aplicación del par obtenido para el escalón 1 al sistema completo se ajusta muy bien en el primer escalón, pero no ol hace tanto a medida que se avanza en la figura. El error también es algo más elevado de lo esperado. Figura 45. Modelo de caja gris y sistema de 1 cero y 2 polos. Aplicación del par obtenido en el primer escalón al sistema completo Escalón 2: El intervalo de puntos correspondiente a los datos de aprendizaje del segundo escalón es (49:107). Par óptimo obtenido: 𝑃𝑎𝑟2=[0.009,126.3854,0.1803,0.0574] Error de 𝑃𝑎𝑟2 respecto al escalón 2: 𝐸𝐶𝑀𝑒𝑠𝑐2 =0.0963 Función de transferencia del modelo evaluando el segundo escalón: 𝐺𝑚𝑜𝑑𝑒𝑙𝑜2 =0.009𝑠+1.137 𝑠2+0.2377𝑠+0.0135 Al representar la curva obtenida frente a la deseada en la Figura 46.a, se observa que no se ajusta del todo a las oscilaciones y presenta un mayor error con respecto al escalón anterior. De todas formas, como la respuesta ante escalón unitario (Figura 46.b) es estable y el tiempo de establecimiento reducido, se sigue considerando un modelo adecuado. El modelo para este segundo escalón se validará en el sistema completo. Identificación en bucle cerrado mediante un algoritmo propio 48 48 Figura 46. Modelo de caja gris y sistema de 1 cero y 2 polos. Resultados obtenidos en el segundo escalón en (a) sistema controlado y (b) respuesta ante escalón unitario Error de 𝑃𝑎𝑟2 con respecto a la planta: 𝐸𝐶𝑀2=6.586 La aplicación del par óptimo calculado para el escalón 2 al sistema completo muestra un buen ajuste en los dos primeros escalones (Figura 47), pero a partir del tercero, las oscilaciones son insuficientes. El error con respecto a la planta sigue siendo elevado. Figura 47. Modelo de caja gris y sistema de 1 cero y 2 polos. Aplicación del par obtenido en el segundo escalón al sistema completo Escalón 3: El intervalo de puntos correspondiente a los datos de aprendizaje del tercer escalón es (108:177). Par óptimo obtenido: 𝑃𝑎𝑟3=[0.0479,35.4105,0.1002,0.1163] Error de 𝑃𝑎𝑟3 respecto al escalón 3: 𝐸𝐶𝑀𝑒𝑠𝑐3 =0.1266 Función de transferencia del modelo evaluando el tercer escalón: 𝐺𝑚𝑜𝑑𝑒𝑙𝑜3 =0.0479𝑠+1.696 𝑠2+0.2165+0.01165 49 Figura 48. Modelo de caja gris y sistema de 1 cero y 2 polos. Resultados obtenidos en el tercer escalón en (a) sistema controlado y (b) respuesta ante escalón unitario En la Figura 48 se aprecia que de nuevo la curva contiene las oscilaciones necesarias, pero no se ajusta exactamente a la deseada. También se puede ver que el sistema es estable, por lo que de nuevo se validará este escalón para el sistema completo. Error de 𝑃𝑎𝑟3 con respecto a la planta: 𝐸𝐶𝑀3=3.6706 En este caso, el ECM ha disminuido. Se puede observar en la Figura 49 que en el escalón donde mejor se ajusta la gráfica es el tercero. Esto tiene sentido, pues se ha tratado de adaptar el par obtenido en este mismo escalón para el sistema completo. Figura 49. Modelo de caja gris y sistema de 1 cero y 2 polos. Aplicación del par obtenido en el tercer escalón al sistema completo Escalón 4: El intervalo de puntos correspondiente a los datos de aprendizaje del cuarto escalón es (178:281). Par óptimo obtenido: 𝑃𝑎𝑟4=[0.5255,5.0905,0.0435,0.0435] Error de 𝑃𝑎𝑟4 respecto al escalón 4: 𝐸𝐶𝑀𝑒𝑠𝑐4 =0.1499 Función de transferencia del modelo evaluando el cuarto escalón: Identificación en bucle cerrado mediante un algoritmo propio 56 56 Una vez más, la respuesta ante escalón es estable y rápida, y el ajuste de la curva adecuado como se puede apreciar en la Figura 60. Figura 60. Modelo de caja gris y sistema de 1 cero y 3 polos. Resultados obtenidos en el cuarto escalón en (a) sistema controlado y (b) respuesta ante escalón unitario Error de 𝑃𝑎𝑟4 con respecto a la planta: 𝐸𝐶𝑀4=3.8835 En la Figura 61 queda representada la aplicación de 𝑃𝑎𝑟4 al sistema completo. Como era de esperar, el mejor ajuste se consigue en el escalón número 4, descendiendo la calidad del mismo en los escalones inferiores. Figura 61. Modelo de caja gris y sistema de 1 cero y 3 polos. Aplicación del par obtenido en el tercer escalón al sistema completo En cada uno de los escalones del sistema con 1 cero y 3 polos, el error mejora con respecto al ensayado con 1 cero y 2 polos y el tiempo de respuesta ante escalón unitario es similar. En cuanto a la validación del sistema completo con los modelos obtenidos en cada escalón, se vuelve a demostrar que no se ajusta de manera adecuada. 3.5.4.2 Aproximación modelo completo Se tratará de inferir una estructura de 1 cero y 3 polos para el sistema completo. Par óptimo obtenido: 57 𝑃𝑎𝑟=[3.1617,0.7072,0.3611,−0.0382,0.8699] Error del par con respecto al sistema: 𝐸𝐶𝑀 =2.3139 Función de transferencia del modelo completo: 𝐺𝑚𝑜𝑑𝑒𝑙𝑜 =3.162𝑠+2.236 𝑠3+1.193𝑠2+0.2671𝑠−0.01199 Figura 62. Modelo de caja gris y sistema de 1 cero y 3 polos. Resultados obtenidos en el sistema completo en (a) sistema controlado y (b) respuesta ante escalón unitario De nuevo, al identificar la función de transferencia del modelo, la respuesta ante escalón unitario resulta inestable (Figura 62.b). Al comparar la salida obtenida y la deseada, se advierte en la Figura 62.a que no se ajustan en los últimos escalones al igual que ha ocurrido en todas las pruebas donde se intenta inferir el modelo de la planta en el sistema completo. Estos dos motivos, prueban que el modelo de la planta es definitivamente no lineal y no es coherente seguir probando esta aproximación. 3.5.4.3 Aproximación multimodelo Finalmente, se volverán a tomar los modelos obtenidos para cada escalón como modelos locales particularizados en función de la referencia de salida: 𝐺𝑚𝑜𝑑𝑒𝑙𝑜1 =0.9708𝑠+1.369 𝑠3+1.326𝑠2+0.3497𝑠+0.00304,𝑟𝑒𝑓 =−0.5 𝐺𝑚𝑜𝑑𝑒𝑙𝑜2 =0.9205𝑠+1.106 𝑠3+1.119𝑠2+0.2534𝑠+0.008872,𝑟𝑒𝑓 =0 𝐺𝑚𝑜𝑑𝑒𝑙𝑜3 =0.4938𝑠+2.133 𝑠3+1.211𝑠2+0.3905𝑠+0.01434,𝑟𝑒𝑓 =0.5 𝐺𝑚𝑜𝑑𝑒𝑙𝑜4 =2.486𝑠+1.08 𝑠3+0.5182𝑠2+008931𝑠+0.005119,𝑟𝑒𝑓=1 El ECM medio de todos los modelos es el siguiente: 𝐸𝐶𝑀𝑚𝑒𝑑𝑖𝑜 = 0.0668 La representación de los sistemas locales integrados con la curva deseada en la Figura 63 se ajustan de manera adecuada. Además, el error ha disminuido con respecto al sistema de 1 cero y 2 polos. Identificación en bucle cerrado mediante un algoritmo propio 58 58 Figura 63. Modelo de caja gris y sistema de 1 cero y 3 polos. Resultados obtenidos en el sistema multimodelo Si se comparan las tres estructuras multimodelo evaluadas (orden configurable en modelo de caja negra, 1 cero y 2 polos y 1 cero y 3 polos en modelos de caja gris), la que presenta un 𝐸𝐶𝑀𝑚𝑒𝑑𝑖𝑜 más reducido es la primera de ellas. El motivo de ello es que el sistema de orden configurable tiene mayor orden que las estructuras analizadas en los modelos de caja gris. En este tipo de modelos, no se han añadido más ceros y polos porque el tiempo de computación aumentaba considerablemente cada vez que se modificaba la estructura. Igualmente, si se hubieran elegido modelos con un orden muy elevado, el ECM sería menor, aunque no en una magnitud tan relevante como la tardanza en el cálculo de los pares adecuados. Por estos motivos, en el siguiente capítulo se utilizarán las herramientas que proporciona MatLab para la identificación de sistemas, con el objetivo principal de disminuir el tiempo de computación y de mejorar el ajuste de los resultados obtenidos. 59 4 IDENTIFICACIÓN EN BUCLE CERRADO MEDIANTE EL TOOLBOX DE MATLAB A computación del algoritmo desarrollado de manera manual ha resultado ser muy lenta y los resultados no se han ajustado exactamente a los deseados. Por esta razón se ha decidido dejar el primer método de identificación a un lado. Siguiendo con la búsqueda de un modelo dinámico adecuado, en este capítulo se utilizará el System Identification Toolbox que ofrece MatLab. Esta herramienta proporciona funciones de MatLab, bloques de Simulink y una aplicación para crear modelos matemáticos de sistemas dinámicos basados en datos medidos de entrada-salida. Permite crear y utilizar modelos de sistemas dinámicos que no se pueden modelar fácilmente a partir de especificaciones o principios básicos. Puede utilizar datos de entrada y salida del dominio del tiempo y el dominio de la frecuencia para identificar funciones de transferencia de tiempo continuo y tiempo discreto, modelos de procesos y modelos de espacio de estados. La caja de herramientas proporciona técnicas de identificación tales como máxima verosimilitud, minimización de errores de predicción (PEM) e identificación de sistemas subespaciales. Para representar la dinámica de un sistema no lineal, puede estimar los modelos de Hammerstein-Wiener y los modelos ARX no lineales con no linealidades de redes sigmoideas, redes Wavelet y particiones en árbol. El Toolbox también realiza la identificación de sistemas de caja gris para estimar los parámetros de un modelo definido por el usuario. Puede utilizar el modelo identificado para la predicción de la respuesta del sistema y el modelado de planta en Simulink. También soporta el modelado de datos de series temporales y la predicción de series temporales. El toolbox de identificación de matlab (SITB) se publicó por primera vez en 1988. Desde entonces ha estado en continuo desarrollo hasta alcanzar su versión actual, la 9.14 que salió con la versión 2021a de matlab. Se ha informado de varios hitos en los Simposios SYSID de la IFAC [83]–[92]. Los cambios producidos a lo largo de los años han conducido a una mejor apreciación de la identificación de sistemas en la comunidad general de sistemas y control. En el manual del usuario [93] se pueden encontrar las diversas aplicaciones y modos de identificación disponibles en el Toolbox de MatLab. El software System Identification Toolbox estima los parámetros del modelo minimizando el error entre la salida del modelo y la respuesta medida. Para determinar la función de transferencia, la caja de herramientas utiliza la diferencia entre la salida del modelo y la salida medida. El criterio de minimización es una norma ponderada del error de predicción, y los algoritmos de estimación ajustan los parámetros de la estructura del modelo de forma que la norma de este error sea lo más pequeña posible. A su vez, se puede configurar el algoritmo de estimación mediante la modificación del criterio de minimización para centrar la estimación en un rango deseado y especificando las opciones de optimización para los algoritmos de estimación iterativos. Mediante esta herramienta, se tratará de inferir el modelo dinámico en bucle cerrado utilizando de nuevo modelos de caja negra y caja gris, y se compararán con los resultados obtenidos previamente a través del algoritmo diseñado de manera manual. L Solo se progresa cuando se piensa que se puede hacer algo más. - Guillermo Marconi - Identificación en bucle cerrado mediante el Toolbox de MatLab 60 60 4.1 Remuestreo uniforme de los datos experimentales La herramienta de identificación trabaja con datos uniformemente muestreados ya sea en el dominio del tiempo o el dominio de la frecuencia. Al haber tomado los datos de manera manual a partir de una imagen, no están distribuidos de manera equidistante en el tiempo. Es debido a este obstáculo que se necesita realizar una interpolación de los datos experimentales para que queden distribuidos uniformemente en el tiempo y así poder utilizar satisfactoriamente la herramienta de identificación de sistemas que proporciona MatLab. Para la interpolación se utiliza el comando resample que remuestrea datos uniformes o no uniformes a una nueva tasa fija especificando el método de interpolación junto con cualquiera de los argumentos de las sintaxis anteriores de este grupo. El método elegido es el de interpolación cúbica a trozos que preserva la forma de la curva. En la Figura 64 se puede apreciar que la curva uniformemente muestreada resultado de la interpolación es exactamente igual a la de datos experimentales no uniformes en el tiempo. Figura 64. Curva resultante de la interpolación. 4.2 Modelos de Caja Negra Para los modelos de caja negra se ha utilizado la aplicación de identificación de sistemas que proporciona MatLab. Esta app realiza la estimación de modelos lineales y no lineales de manera interactiva importando datos de entrada-salida medidos del dominio del tiempo y el dominio de la frecuencia (Figura 65). Estas estructuras de modelos varían en complejidad dependiendo de la flexibilidad que se necesite para tener en cuenta la dinámica y el ruido de su sistema. Se puede elegir una estructura, configurar el orden del modelo y calcular sus parámetros para que se ajusten a los datos de respuesta medidos. Los datos pueden preprocesarse realizando operaciones tales como la eliminación de tendencias, filtrado, remuestreo y reconstrucción de datos ausentes. A la hora de estimar y validar el sistema, se comparan los modelos identificados, analizando sus propiedades, calculando sus límites y validando con respecto al conjunto de datos de prueba. Tras la validación, si las estructuras de modelo simples no producen buenos resultados, se pueden seleccionar estructuras de modelo más complejas: especificando un orden de modelo superior para la misma estructura de modelo lineal, modelando explícitamente el ruido o utilizando una estructura de modelo no lineal. Finalmente, la estructura seleccionada será la más simple que proporcione el mejor ajuste a los datos medidos. A la hora de inferir el modelo dinámico, se ha hecho uso de diversas estructuras para tratar de buscar la que más se acerque al resultado deseado. 61 Figura 65. Interfaz de la System Identification App 4.2.1 Modelos ARX La primera estructura elegida es la ARX. El orden de la misma será elegido por la aplicación entre un rango (Figura 66.a). Una vez estimamos el modelo, la aplicación recomienda los pares óptimos (Figura 66.b). Los tres pares recomendados son: • Verde: ARX690, 6 polos y 9 ceros • Azul: ARX14130, 14 polos y 13 ceros • Rojo: ARX1515, 15 polos y 15 ceros Figura 66. Selección de estructura ARX (a) panel de selección de pares a estimar y (b) pares recomendados Probando las diferentes configuraciones, se han encontrado otras que también presentan un máximo de Identificación en bucle cerrado mediante el Toolbox de MatLab 62 62 precisión: • ARX12130, 12 polos y 13 ceros • ARX13130, 13 polos y 13 ceros Las gráficas de estas estructuras con respecto a la deseada están representadas en la Figura 67 y se percibe que la curva se ajusta en todos los escalones, logrando una precisión máxima del 89.95%. Visualmente, se ajusta mejor a los escalones intermedios que al primero y al último. Figura 67. Comparación de las estructuras de modelos ARX con la deseada 4.2.2 Modelos de espacios de estado Se pasa a estimar la estructura a partir de modelos de espacios de estado. De nuevo se le indica a la aplicación la búsqueda del mejor modelo en un rango de 15 (Figura 68.a) y sugiere que es el de orden 7 (Figura 68.b). Figura 68. Selección de estructura de espacios de estado (a) panel de selección del orden a estimar y (b) orden recomendado 63 Como se puede ver en la Figura 69, el sistema de orden 7 adquiere un ajuste del 82.67%; sin embargo, probando otros órdenes, los sistemas de orden 11, 13 y 15 superan en precisión al que nos sugiere la aplicación, llegando a un ajuste máximo del 90.03%. La gráfica se ajusta adecuadamente a los cuatro escalones, aunque se observan mayores diferencias en el cuarto. Figura 69. Comparación de las estructuras de modelos de espacios de estado con la deseada 4.2.3 Modelos de procesos El siguiente tipo de modelo a estimar es el de procesos. La estructura se puede elegir configurando el número de polos que pueden ser reales o que el sistema pueda resultar subamortiguado (contiene polos complejos); también es posible añadir un cero, un retraso y/o un integrador (Figura 70). Figura 70. Panel de selección de la estructura del modelo de procesos Tras probar diversas configuraciones, las que han presentado mejores resultados han sido: • P2U: 2 polos (sistema subamortiguado). Identificación en bucle cerrado mediante el Toolbox de MatLab 64 64 • P2DZU: 2 polos (sistema subamortiguado), 1 cero y un retraso. • P3ZU: 3 polos (sistema subamortiguado) y 1 cero. Una vez más los modelos obtenidos delineados en la Figura 71 se ajustan en mayor o menor medida a todos los escalones y la precisión máxima obtenida es similar a los ensayos anteriores (89.18%). Figura 71. Comparación de las estructuras de modelos de procesos con la deseada 4.2.4 Modelos de funciones de transferencia Las estructuras estimadas en esta simulación serán funciones de transferencia con polos y ceros elegidos por el usuario (Figura 72). No se han detectado resultados prometedores hasta que el modelo no cuenta con al menos dos polos, donde se comienza a ver un incremento notable de la eficiencia del modelo estimado. Se analizarán estructuras de dos polos en adelante. Figura 72. Panel de selección de la estructura del modelo de funciones de transferencia 65 Se han estimado modelos de entre dos y diez polos. Las configuraciones de cada número de polos con mayor porcentaje de eficiencia están dispuestas en la Tabla 5. Analizando la tabla, se percibe que a medida que se incrementa el número de ceros y polos la precisión aumenta, pero varía muy poco. Se llegan a niveles de casi el 90% de exactitud con la función deseada. A medida que se aumenta el orden de la función de transferencia se presenta un dilema: a mayor grado, mayor complejidad de la función, cuestión que no interesa a nivel computacional y algebraico. La elección del modelo final deberá encontrar un equilibrio entre un orden reducido y un mínimo de eficiencia. FUNCIÓN DE TRANSFERENCIA Nº POLOS Nº CEROS AJUSTE TF_2_0 2 0 88.86% TF_3_2 3 2 89.18% TF_4_4 4 4 89.52% TF_5_1 5 1 89.39% TF_6_2 6 2 89.64% TF_7_2 7 2 89.61% TF_8_4 8 4 89.72% TF_9_6 9 6 89.70% TF_10_6 10 6 89.71% Tabla 5. Ceros, polos y eficiencia tras estimar modelos de función de transferencia En la Figura 73 se pueden distinguir las funciones de transferencia estimadas para cada par de polos y ceros óptimo. El máximo ajuste se consigue en la configuracióne de 8 polos y 4 ceros; sin embargo, si se quiere trabajar con una estructura con un orden menor, una estructura adecuada sería la de 6 polos y 2 ceros o la de 4 polos y 4 ceros. Figura 73. Curvas resultantes tras estimar a partir de modelos de funciones de transferencia. Identificación en bucle cerrado mediante el Toolbox de MatLab 72 72 a los obtenidos utilizando el algoritmo diseñado en el Capítulo 3. Este método a través de la herramienta de MatLab no ha resultado igual de satisfactorio que el de caja negra. Esto es debido a la condición de no linealidad que mantiene el sistema a identificar. Dos de los métodos de identificación de caja negra disponibles en el System Identification Toolbox son para sistemas no lineales y, como consecuencia, el ajuste en esos casos es muy exacto. 73 5 APLICACIÓN A UN SISTEMA DE TRES EYECTORES e tiene un ciclo de refrigeración multieyector. Gracias a lo obtenido en los capítulos anteriores es posible determinar el tiempo que tarda en estabilizarse el ciclo cada vez que se produce un cambio de un eyector a otro. Este capítulo analizará la viabilidad de este MERS a partir de un histórico de temperaturas anual de un municipio sevillano, donde se analizarán los cambios de un eyector a otro en función de la temperatura exterior, tratando de determinar si en alguna situación se produciría una variación que no debiera. 5.1 El ciclo MERS Las investigaciones sugieren que, para obtener el COP óptimo, la geometría del eyector deberá ser diferente para cada condición de funcionamiento [94]–[102]. A parte de utilizar eyectores de geometría variable frente a un sistema de un eyector fijo, también se puede adoptar un sistema multieyector. Figura 80. Diagrama del ciclo MERS Un MERS utiliza eyectores conectados paralelamente entre sí, habilitando o deshabilitando cada uno de ellos para modificar la capacidad y el rendimiento del sistema. La literatura trata sistemas de tres y cuatro eyectores simples con geometrías optimizadas para aumentar el rendimiento en condiciones de funcionamiento flexibles, S La creatividad es la inteligencia divirtiéndose. - Albert Einstein - Aplicación a un sistema de tres eyectores 74 74 siendo un sistema de tres eyectores el que presenta una mayor eficiencia [103]. El diagrama de los componentes del MERS utilizado es el representado en la Figura 80. 5.2 Evaluación del rendimiento estacional La geometría de cada uno de los tres eyectores perteneciente al sistema ha sido optimizada para el aumento de la eficiencia frente a condiciones de funcionamiento variables. Para cada generación y temperatura exterior, se optimiza el MERS basándose en el método de la dirección conjugada utilizando el sofware EES [95]. Una vez se tienen los tres eyectores con las geometrías óptimas, se simula el rendimiento de cada eyector frente a diferentes temperaturas exteriores y el rendimiento estacional del MERS en un clima cálido considerando el perfil anual de un municipio sevillano. Se observa en la Figura 81 que cada eyector presenta una curva de rendimiento que se ajusta a un rango determinado de temperaturas donde el punto de corte entre curva y curva deberá ser a la temperatura a la que se deberá cambiar de un eyector a otro para lograr un mayor COP del sistema. Figura 81. Simulación estacional del rendimiento del sistema multieyector [95] La temperatura exterior a la que se activará el sistema es de 25ºC. Como se observa, se requieren tres ratios de área diferentes para cada módulo optimizado para maximizar el COP medio, requiriendo mayores valores de Ar para temperaturas de evaporación más bajas. De este modo, cuando el rendimiento de un eyector cae por debajo del rendimiento del eyector consecutivo correspondiente, el primer módulo se apaga mientras el segundo se activa, maximizando el rendimiento en todo el rango de funcionamiento. En concreto, el Ar optimizado para el rango más bajo de temperaturas exteriores es de 11,01, que corresponde a un COPcc crítico de 0,37. A continuación, el módulo cambia al Ar de 9,56 por encima de 31ºC, consiguiendo un COPcc de 0,36. Y, por último, un Ar de 8,3 para temperaturas exteriores superiores a 35ºC, relacionado con un COPcc de 0,34. Como resultado, el COP estacional obtenido mediante el uso de tres módulos optimizados del eyector es de 0,3 [95]. En definitiva, a mayor relación de área, mayor COP alcanzable, aunque menor es la temperatura crítica de condensación. 5.3 Análisis de viabilidad del MERS Gracias al histórico de temperaturas exteriores, se pueden ver todos los cambios realizados por el sistema entre un eyector y otro. Sin embargo, surge la cuestión de si es factible o no realizar la permutación entre un eyector u otro ya que en el proceso se pasa por un régimen transitorio que desajusta ligeramente el COP. Como se ha 75 observado durante el proceso de identificación en los capítulos anteriores, la salida equivale a la presión del enfriador de gas y la entrada corresponde a la variable manipulable que es la geometría del eyector. Cada vez que la presión del enfriador de gas varía, el eyector modifica su geometría para alcanzar el máximo COP posible. El artículo de Yang He [17] proporciona un gráfico que relaciona la variación de la geometría del eyector con la presión del enfriador de gas y el COP del sistema (Figura 82). Figura 82. Estados estables para diferentes áreas de la garganta de la boquilla [17]. Se ha elegido un modelo dinámico identificado en el capítulo 3 teniendo en cuenta el ECM total del sistema. El modelo seleccionado ha sido el sistema de orden configurable multimodelo presente en el apartado 3.4.3.3 y se asumirá esta misma dinámica para el MERS. En la Figura 2 se observa el tiempo que tarda en estabilizarse el sistema ante cambios de presión. El régimen transitoio no supera los 50s en ninguna ocasión, por lo que, en el MERS, cada vez que se realice una conmutación de eyectores se asumirá el mismo tiempo de establecimiento. Figura 83. Histórico anual de temperaturas y cambios de eyector En la Figura 83 se muestra el histórico anual de temperaturas y los consecuentes cambios que realizan los eyectores. En los primeros meses se puede apreciar que, al no superar los 25ºC la temperatura a lo largo del día, el sistema de refrigeración permanece apagado. A medida que avanzan los primeros meses y las temperaturas 0 0.5 1 1.5 2 2.5 3 3.5 0 5 10 15 20 25 30 35 40 45 50 02000 4000 6000 8000 10000 Eyector Temperatura (°C) Horas Temperatura Eyector Aplicación a un sistema de tres eyectores 76 76 aumentan, el sistema comienza a encencerse y utilizar únicamente un eyector. Progresivamente, a medida que se acerca a las estaciones calurosas, se van encendiendo los eyectores 2 y 3 y empiezan a utilizarse todos ellos en un mismo día. De nuevo, cuando las temperaturas vuelcen a disminuir, se producen menos permutaciones entre los eyectores. Ahora se analizarán una serie de semanas correspondientes a diferentes épocas del año para poder ver mejor estas transiciones entre eyectores. El primer encendido del año del eyector 1 ocurre en la hora 2054, día 85, que correspondería a la semana del 26 de marzo (Figura 84). La puesta en marcha es viable ya que de una hora a otra no hay una variación brusca de temperatura y, por lo tanto, gracias al bajo tiempo de establecimiento, no supondría una pérdida de rendimiento. Figura 84. Histórico de la semana del 26-03 de temperaturas y cambios de eyector Ahora se fija la vista en una semana calurosa, la del 17 de julio representada en la Figura 85. En estos 7 días, las temperaturas superan los 40ºC, es decir, el eyector 3 estará bastante activo, intercalando su funcionamiento con el 1 y el 2. Igualmente, en ciertos momentos de la semana, al descender la temperatura por debajo de los 25ºC, el sistema estará en fase de reposo. De nuevo, no se producen oscilaciones grandes de temperaturas y los cambios entre eyectores se producen en amplios rangos de tiempo, por lo que resulta viable la transición entre eyectores. Figura 85. Histórico de la semana del 17-07 de temperaturas y cambios de eyector Finalmente se concluye que, con los datos proporcionados, la implementación de un MERS es viable ya que el tiempo de establecimiento es pequeño. Sin embargo, esto es debido a que las variaciones de temperatura se presentan en el gráfico de manera horaria. Si el histórico de temperaturas estuviera representado minuto a minuto, podría evaluarse de manera más precisa el efecto del transitorio en los cambios de eyector. 0 1 2 0 5 10 15 20 25 30 2039 2059 2079 2099 2119 2139 2159 2179 2199 Eyector Temperatura (°C) Horas Temperatura Eyector 0 1 2 3 4 10 15 20 25 30 35 40 45 4740 4760 4780 4800 4820 4840 4860 4880 4900 4920 4940 Eyector Temperatura (°C) Horas Temperatura Eyector 77 6 CONCLUSIONES a producción de frío ha sido una necesidad para el ser humano desde la antigüedad. Actualmente, los sistemas de refrigeración desempeñan un papel muy importante tanto en la industria como en la vida cotidiana. El consumo de energía eléctrica y combustibles fósiles ha aumentado considerablemente y como consecuencia, se han manifestado efectos secundarios nocivos como el agotamiento de la capa de ozono y el calentamiento global. En 1987 se firmó el Protocolo de Montreal para regular la producción y el comercio de sustancias que dañan y agotan la capa de ozono y en 1997 se firmó el Protocolo de Kioto, comprometiendo a los países firmantes a limitar y reducir las emisiones de GEI. La investigación sobre la refrigeración con eyectores ha tenido un fuerte auge a partir de la preocupación por el medio ambiente, el consumo de la energía y el desarrollo sostenible, ya que se tratan de equipos económicos, fiables y mejoran el rendimiento del ciclo en el que están integrados. Su función es la de aumentar la presión de un fluido secundario a través de uno primario. El eyector está formado por cuatro partes: la boquilla primaria, la boquilla secundaria, la cámara de mezcla y el difusor. A su vez, los eyectores se pueden clasificar según su tipo de boquilla principal, que puede ser convergente o convergente-divergente; su tipo de cámara de mezcla (CAM o CPM); la fase en la que se encuentre el fluido, monofásico o bifásico; y si el eyector tiene un área fija o variable. También se han estudiado los parámetros que influyen en el rendimiento del eyector, obteniendo que la succión aumenta a medida que se eleva la relación de presión y se puede ver alterada por la temperatura exterior del aire; la eficiencia disminuye cuando se incrementa la relación del diámetro. El diseño geométrico es uno de los pasos más importantes para el funcionamiento correcto del eyector y se trata de optimizar para conseguir la mayor eficiencia posible. Los ERS pueden modificar su estructura para adaptarse a la aplicación a la que van a ser sometidos. Finalmente, el fluido de trabajo utilizado es de vital importancia en términos de operación y rendimiento del ciclo y se pueden clasificar en función de la pendiente de la línea de vapor saturado en el diagrama T-s (húmedo, seco o isentrópico) o a partir de su composición química molecular (halocarbonos, hidrocarburos, compuestos y otros). El refrigerante elegido en este trabajo ha sido el CO2 debido a sus propiedades poco contaminantes; sin embargo, el rendimiento del TCRE varía mucho con condiciones de funcionamiento variables, especialmente a diferentes presiones del enfriador de gas. Por lo tanto, un eyector con estructura variable ampliaría sus condiciones efectivas de funcionamiento pudiendo alcanzar de manera consistente un rendimiento óptimo y estable ajustando el área de la garganta de tobera bajo condiciones variables. Dado que la presión del enfriador de gas es un parámetro clave, se ha desarrollado un controlador en cuasi cascada que ajusta el área de la boquilla del eyector del ciclo TCRE a una posición óptima para variaciones de esta presión. La base para el diseño del sistema de control es la característica dinámica del sistema. El objetivo es tratar de mantener un COP óptimo variando el área de la garganta de tobera frente a cambios en la presión del enfriador de gas. Para ello, cada vez que surge una variación en la presión, se estima una nueva área de la garganta de tobera utilizando un algoritmo PID (4) para reducir el error entre la presión real del enfriador de gas. Como no se tiene conocimiento de la dinámica, se debe identificar. La identificación sistemas trata de estimar modelos matemáticos de sistemas dinámicos a partir de observaciones y mediciones realizadas en el proceso basándose en entradas, salidas y perturbaciones. Dos aproximaciones han sentado las bases de esta disciplina: el método de predicción del error y el enfoque no paramétrico del (sub)espacio de estados. El conocimiento exhaustivo del modelo dinámico de un proceso industrial permite comprender el comportamiento del sistema y sus partes críticas. Es necesario aplicar la identificación con el fin de construir un modelo matemático del sistema basado en las entradas y salidas del mismo. El procedimiento de identificación es un proceso iterativo con varias fases, cuyo objetivo es alcanzar el modelo matemático adecuado en base a datos experimentales. Los pasos de este proceso son: 1. Diseño del del experimento. 2. Tratamiento de los datos experimentales. L Conclusiones 78 78 3. Elección de la estructura del modelo: • Tipo de modelo: modelos de caja blanca, caja negra, caja gris y modelos difusos. • Tipo de estructura: no paramétrica (análisis de correlación y espectral), paramétrica lineal (autorregresivos, OE, BJ, espacio de estados, modelos de función de transferencia) y paramétrica no lineal (autorregresivos, Hammerstein-Wiener, series de volterra y redes neuronales). 4. Selección del método de estimación de los parámetros. 5. Validación del modelo. Una vez comienza el proceso de identificación, es necesario reconocer ante qué sistema de control se encuentra el usuario y si este posee o no realimentación. Un sistema en bucle abierto es aquel en el que únicamente actúa el proceso sobre la señal de entrada. En cambio, en un sistema de control en bucle cerrado, la acción de control se encuentra en función de la señal de salida y se puede abordar mediante una aproximación directa, indirecta o conjunta. Como únicamente se poseen como datos el PID del controlador y una gráfica que hace referencia a la presión del enfriador de gas (Figura 2), lograr la correlación de esta presión con el área de la boquilla variable es esencial para lograr el COP óptimo del sistema para cada condición de funcionamiento. La identificación se deberá realizar en bucle cerrado debido a la realimentación a través de diversos métodos. El primer paso del proceso de identificación es el diseño del experimento, donde se han tomado datos a mano a través de MatLab de una gráfica que representaba la evolución de la presión del sistema frente al tiempo. Las coordenadas han sido devueltas como píxeles, por lo que se ha realizado una conversión de píxeles a los valores “reales” del tiempo y la presión. Tras recopilar los puntos, se realiza el tratamiento de datos para suavizar posibles errores cometidos a la hora de la obtención de la muestra. Este procedimiento ha consistido en la realización de un promediado de todos los puntos para obtener una suavización de la curva sin cambios bruscos. El resultado de la toma de datos experimentales y su filtrado puede observarse en la Figura 19. Dado que se proporcionan los datos del controlador, a la hora de realizar la identificación del sistema en bucle cerrado, se ha aplicado la aproximación indirecta. Se han elegido tanto un modelo de caja negra como uno de caja gris para inferir la función de transferencia del modelo dinámico. Inicialmente se ha aplicado un modelo de caja negra variando el orden del sistema. El algoritmo creado para el modelo de caja negra optimiza los parámetros de la función en bucle cerrado, para luego inferir la función de transferencia de la planta de manera algebraica. Se ha comenzado identificando un modelo de primer orden que no ha arrojado resultados satisfactorios, al no realizar ninguna oscilación la curva, poseer un error muy alto y resultar inestables las respuestas ante un escalón unitario. Se ha cambiado la estructura a una de segundo orden, donde los cálculos han mejorado, aunque no tanto como para resultar aceptable el modelo; en esta ocasión, la respuesta escalón a escalón oscila, sigue siendo inestable frente a una entrada escalón unitario y el modelo optimizado para la planta completa no se toma como válido (carece de oscilación y resulta inestable). Finalmente, se aplica un sistema de orden configurable añadiendo dos ceros y un polo al sistema de segundo orden. Las respuestas de este último sistema escalón a escalón se ajustan muy bien a las curvas deseadas y son estables frente a una entrada escalón unitario; sin embargo, cuando se valida el modelo de cada escalón en la planta completa, el error aumenta y únicamente uno o dos escalones se ajustan aceptablemente a la curva. Es por ello que se confirma la posible no lienalidad del sistema y se propone utilizar una aproximación multimodelo donde se utiliza un modelo local u otro dependiendo del rango de presiones en el que se encuentre el sistema aportando como resultados valores muy prometedores a la vez que un ECM medio pequeño. Ahora se pasa a realizar la identificación mediante un modelo de caja gris, donde se optimizan directamente los parámetros de la función de transferencia del modelo. En este caso, las estructuras de primer y segundo orden no han podido ser utilizadas ya que, al estar infiriendo la función realimentada, los ceros y polos de estos sistemas no son suficientes para constituir un sistema válido; es ahora cuando se ha recurrido al método del Lugar de las Raíces. Determinando la frecuencia natural, la sobreoscilación y añadiendo ceros y polos, es posible averiguar qué estructuras resultarán estables. Se han realizado los mismos experimentos que en el modelo de caja negra, pero esta vez para dos estructuras: una de ellas con dos polos y un cero y la otra con tres polos y un cero. En ambas, las respuestas escalón a escalón son precisas y estables, aunque al evaluar la planta completa el error 79 aumenta considerablemente. De nuevo se aplica un sistema no lineal multimodelo particularizando un modelo distinto para cada tramo de presiones del enfriador de gas. En los modelos de caja gris, el tiempo de computación ha resultado similar, aunque algo mayor que el de los modelos de caja negra. Además, el error en los sistemas más precisos de cada modelo es parecido. Como el tiempo de computación ha resultado elevado y los sistemas pueden mejorar en ajuste, se ha optado por utilizar el System Identificatión Toolbox que proporciona MatLab y así inferir una vez más el sistema dinámico de la planta. Para ello, primero se ha tenido que realizar una interpolación de los datos experimentales para que estuvieran uniformemente muestreados en el tiempo. Se han probado diversas estructuras lineales de caja negra: ARX, modelos en el espacio de estados, modelos de procesos y modelos de función de transferencia, que han rozado ajustes del 90%. Estas estructuras han proporcionado resultados análogos a los obtenidos a partir del algoritmo diseñado en el tercer capítulo, aunque en un tiempo de computación mucho mas reducido, permitiendo analizar un mayor número de estructuras del modelo. • Modelos ARX: 13 polos y 13 ceros, 89.95%. • Modelos del espacio de estados: Sistema de orden 15, 90.03%. • Modelos de procesos: 3 polos (sistema subamortiguado) y 1 cero, 89.18%. • Modelos de función de transferencia: 8 polos y 4 ceros, 89.72%. También se han probado modelos no lineales ARX (98.21% precisión) y de Hammerstein-Wiener (97.02% precisión), cuyo ajuste ha sido el más elevado, confirmando una vez más que el sistema a identificar es no lineal. A través del Toolbox de Matlab, también se ha propuesto inferir la planta a través de un modelo de caja gris, pero al no poder incluir no linealidades, los resultados no han sido prometedores y han salido semejantes a los obtenidos en el algoritmo diseñado del Capítulo 3. Finalmente, se ha querido probar en un MERS si el tiempo de establecimiento que surge al cambiar de un eyector a otro es lo suficientemente pequeño como para que resulte viable la transición. El MERS trieyector pretende, a partir de la temperatura exterior, elegir el eyector que garantice un mayor COP dependiendo del rango de temperaturas en que se encuentre. Para ello, se ha tomado un histórico de temperaturas anual de un municipio sevillano. Como modelo dinámico de referencia, se ha escogido el sistema multimodelo obtenido en los experimentos de caja negra realizados para un sistema de orden configurable (apartado 3.4.3.3). El régimen transitorio en el sistema multimodelo nunca supera los 50s, por ello, al tener el registro en formato horario, 50s resultaría un tiempo insignificante frente a una hora y el MERS sería completamente viable. Conclusiones 80 80 81 Bibliografía [1] M. Isaac and D. P. van Vuuren, “Modeling global residential sector energy demand for heating and air conditioning in the context of climate change,” Energy Policy, vol. 37, no. 2, pp. 507–521, Feb. 2009, doi: 10.1016/J.ENPOL.2008.09.051. [2] M. J. Molina and F. S. Rowland, “Stratospheric sink for chlorofluoromethanes: chlorine atom-catalysed destruction of ozone,” Nat. 1974 2495460, vol. 249, no. 5460, pp. 810–812, 1974, doi: 10.1038/249810a0. [3] J. W. J. Bouma, “Global warming and heat pumps,” Heat Pumps Energy Effic. Environ. Prog., pp. 33– 42, 1993, doi: 10.1016/B978-0-444-81534-7.50011-5. [4] W. B. Gosney, “Principles of refrigeration,” p. 666, 1982, Accessed: Oct. 29, 2021. [Online]. Available: https://books.google.com/books/about/Principles_of_Refrigeration.html?hl=es&id=_WV5QgAACAA J. [5] W. F. Stoecker, “Steam-jet refrigeration,” Refrig. Air Cond., pp. 194–205, 1958. [6] C. Li, Y. Li, W. Cai, Y. Hu, H. Chen, and J. Yan, “Analysis on performance characteristics of ejector with variable area-ratio for multi-evaporator refrigeration system based on experimental data,” Appl. Therm. Eng., vol. 68, no. 1–2, pp. 125–132, Jul. 2014, doi: 10.1016/j.applthermaleng.2014.04.031. [7] H. Vidal and S. Colle, “Simulation and economic optimization of a solar assisted combined ejectorvapor compression cycle for cooling applications,” Appl. Therm. Eng., vol. 30, no. 5, pp. 478–486, Apr. 2010, doi: 10.1016/J.APPLTHERMALENG.2009.10.008. [8] G. Besagni, R. Mereu, and F. Inzoli, “Ejector refrigeration: A comprehensive review,” Renewable and Sustainable Energy Reviews, vol. 53. Elsevier Ltd, pp. 373–407, Jan. 01, 2016, doi: 10.1016/j.rser.2015.08.059. [9] J. qiang Deng, P. xue Jiang, T. Lu, and W. Lu, “Particular characteristics of transcritical CO2 refrigeration cycle with an ejector,” Appl. Therm. Eng., vol. 27, no. 2–3, pp. 381–388, Feb. 2007, doi: 10.1016/j.applthermaleng.2006.07.016. [10] D. Li and E. A. Groll, “Transcritical CO2 refrigeration cycle with ejector-expansion device,” Int. J. Refrig., vol. 28, no. 5, pp. 766–773, Aug. 2005, doi: 10.1016/j.ijrefrig.2004.10.008. [11] S. Fangtian and M. Yitai, “Thermodynamic analysis of transcritical CO2 refrigeration cycle with an ejector,” Appl. Therm. Eng., vol. 31, no. 6–7, pp. 1184–1189, May 2011, doi: 10.1016/J.APPLTHERMALENG.2010.12.018. [12] Y. He, J. Deng, and Z. Zhang, “Thermodynamic study on a new transcritical CO2 ejector expansion refrigeration system with two-stage evaporation and vapor feedback,” http://dx.doi.org/10.1080/10789669.2014.929422, vol. 20, no. 6, pp. 655–664, Aug. 2014, doi: 10.1080/10789669.2014.929422. [13] C. Lucas and J. Koehler, “Experimental investigation of the COP improvement of a refrigeration cycle by use of an ejector,” in International Journal of Refrigeration, Sep. 2012, vol. 35, no. 6, pp. 1595–1603, doi: 10.1016/j.ijrefrig.2012.05.010. [14] F. Liu, Y. Li, and E. A. Groll, “Performance enhancement of CO 2 air conditioner with a controllable ejector,” in International Journal of Refrigeration, Sep. 2012, vol. 35, no. 6, pp. 1604–1616, doi: