Diseño del sistema de control de un evaporador

i Equation Chapter 1 Section 1 Trabajo Fin de Grado Grado de Ingeniería en Tecnologías Industriales Diseño del Sistema de Control de un Evaporador Autor: Diego Sokolowski Barrón Tutor: Prof. Dr. Francisco Javier Gutiérrez Ortiz Departamento de Ingeniería Química y Ambiental Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2016 iii Trabajo Fin de Grado Grado de Ingeniería en Tecnologías Industriales Diseño del Sistema de Control de un Evaporador Autor: Diego Sokolowski Barrón Tutor: Dr. Francisco Javier Gutiérrez Ortiz Profesor Titular Dep. de Ingeniería Química y Ambiental Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2016 v Trabajo Fin de Grado: Diseño del Sistema de Control de un Evaporador Autor: Diego Sokolowski Barrón Tutor: Prof. Dr. Francisco Javier Gutiérrez Ortiz El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: El Secretario del Tribunal Fecha: vii Agradecimientos A mi familia y amigos por su constante apoyo durante toda la carrera en general y durante este proyecto en particular. A mi tutor Francisco Javier Gutiérrez Ortiz por su dedicación y tiempo prestado. Diego Sokolowski Barrón Sevilla, 2016 ix Resumen En este proyecto se desarrolla un modelo matemático de un evaporador de tubos verticales con recirculación forzada, y se diseña un sistema de control. En primer lugar, se definen las ecuaciones y variables que constituyen el proceso, especificando las salidas y entradas del sistema. Luego, se identifican los emparejamientos salida-entrada más idóneos para un control SISO, mediante la RGA. Posteriormente, se implementa el modelo no lineal en Simulink-MATLAB®, para poder observar las respuestas del sistema frente a cambios en sus variables de entrada. Gracias a la utilidad Ident Toolbox de MATLAB® se identifican modelos simplificados lineales, a partir de los cuales se aborda el diseño de los controladores por realimentación más adecuados. Finalmente, se implementan controladores avanzados anticipativos y en cascada, con el objetivo de reducir los efectos indeseables provocados por perturbaciones. Paralelamente a todo este proceso, se ha desarolllado una interfaz gráfica GUI en Simulink, que permite al usuario interactuar con el sistema controlado para observar los resultados directamente a tiempo real. 1 1 INTRODUCCIÓN 1.1 Objetivo del proyecto El objetivo de este proyecto es diseñar el sistema de control para un evaporador de tubos largos verticales de recirculación forzada a partir del desarrollo previo de un modelo matemático del sistema. Adicionalmente, se trabajará en un entorno gráfico de simulación que permita observar las respuestas del sistema en tiempo real al variar las variables de perturbación o los puntos de consigna del sistema controlado. Dicho entorno, que se implementará en MATLAB®-Simulink, permite el estudio y la visualización de los resultados del proceso de una manera más rápida e intuitiva. Dado que los evaporadores industriales disponen de múltiples entradas manipulables, perturbaciones, y salidas, la complejidad de la solución de control radica en la naturaleza multivariable de la planta a estudiar. Por tanto, el desarrollo de este documento se centra en tres bloques principales:  Recopilación de información y obtención del modelo dinámico.  Desarrollo de la solución de control del sistema multivariable.  Diseño del entorno gráfico. En comparación con otros estudios previos, este trabajo profundiza en mayor detalle en el modelo dinámico del sistema y ofrece una interfaz gráfica que dota al conjunto de un carácter didáctico. 1.2 El sistema evaporador Un evaporador es un intercambiador de calor en el que se produce una transferencia de energía térmica entre un fluido caliente en estado gaseoso (típicamente vapor de agua) y un fluido líquido a menor temperatura. El fluido gaseoso se enfría y se condensa, y el fluido líquido se calienta. Su nombre se debe al cambio de estado producido en el líquido, pues parte de éste se evapora al recibir el calor del otro fluido. Los evaporadores se encuentran en todo tipo de industrias. Uno de sus usos más frecuentes es como sistema de refrigeración, en equipos como cámaras frigoríficas, sistemas de aire acondicionado y neveras. También son comúnmente utilizados en procesos químicos de precipitación, cristalización y extracción de líquidos, entre otros. Dependiendo de su aplicación y carga térmica su diseño es diferente, variando su tamaño, capacidad y estructura. Existen diferentes configuraciones, entre las que destacan: evaporador de película descendente, evaporador de película ascendente, evaporador de tubos verticales horizontales o verticales con recirculación forzada, evaporador de placas, etc. Debido a su especial diseño, los evaporadores con tubos verticales y recirculación forzada son adecuados para concentrar disoluciones con tendencia a la cristalización, incrustación, o que tienen gran viscosidad. Se usan mayoritariamente en la industria alimenticia e industria farmacéutica y química, así como en procesos medioambientales como el tratamiento de aguas. Es el tipo de evaporador elegido en este proyecto. Dicho sistema evaporador consta de: Equipos  Intercambiador de calor de tubos verticales (evaporador).  Separador flash líquido-gas.  Condensador.  Tubería de recirculación entre el separador líquido-gas y el intercambiador de calor, así como una bomba de impulsión para la recirculación. 2 Corrientes  Una entrada de alimentación de disolución en el circuito de recirculación.  Una entrada de vapor en la cámara exterior a los tubos del intercambiador de calor.  Una corriente de entrada de agua refrigerada en el condensador.  Una salida de concentrado del circuito de recirculación.  Una salida de condensado en el intercambiador de calor.  Una salida de condensado procedente del flujo del vapor del separador que pasa a través del condensador . El proceso de funcionamiento es el siguiente: 1) La solución líquida que se quiere concentrar es alimentada por el tubo de entrada inferior al evaporador y fluye de abajo hacia arriba a través de los tubos impulsada por la bomba. 2) En los tubos del evaporador el líquido se caliente y parte del agua se evapora debido al vapor presente en la cámara que rodea a los tubos. De este modo, se genera en el interior de los tubos una mezcla líquido-vapor saturada. En concreto, al aumentar el caudal de alimentación de vapor (ya sea por la variación de la presión de alimentación proveniente de una caldera o por la variación en la apertura de la válvula), aumentará la presión en el interior de la cámara del intercambiador, lo que incrementará también su temperatura y, por tanto, el flujo de calor hacia el líquido contenido en los tubos. 3) A la salida de los tubos del evaporador, el vapor de agua se separa del líquido concentrado en un separador líquido-vapor. 4) El vapor de agua producido es conducido a través de un condensador. 5) El líquido concentrado del separador es recirculado de nuevo hacia el evaporador mediante una bomba. Aguas arriba de la bomba se encuentra una salida de concentrado, y posteriormente de nuevo la entrada del líquido a concentrar. El esquema del sistema evaporador con sus variables se muestra en la Figura 1.1 y Tabla 1.1 [6]. Figura 1.1 Esquema de funcionamiento Evaporado 3 Variables Descripción Punto nominal Máximo Mínimo Unidades F1 Caudal de alimentación 10.0 kg/min F2 Caudal de líquido salida evaporador 52.0 - - kg/min F2′ Caudal de producto 2.0 5 0 kg/min F3 Caudal de recirculación 50.0 - - kg/min F4 Caudal de vapor 8.0 - - kg/min F5 Caudal de condensado 8.0 - - kg/min X1 Concentración de alimentación 5.0 % (kg s. kg d.)* X2 Concentración en el separador 25.0 40 0 % (kg s. kg d.) X2′ Concentración de producto de salida 25.0 40 0 % (kg s. kg d.) T1 Temperatura de alimentación 40.0 OC T2 Temperatura de producto de salida 84.6 - - OC T4 Temperatura del vapor de agua 80.6 - - OC L2 Nivel del separador 1.0 2 0 m P0 Presión de vapor de la caldera 500 kPa P2 Presión en el evaporador 50.5 100 0 kPa F100 Caudal de vapor de alimentación 9.3 - - kg/min XF100 Apertura de válvula de caudal de vapor 0.5 1 0 - Fc Caudal de condensado de la carcasa 9.3 - - kg/min T100 Temperatura de vapor de alimentación 119.9 - - OC P100 Presión de vapor de alimentación 194.7 - - kPa Q100 Flujo de calor en el intercambiador 339.0 - - kW F200 Caudal de refrigeración 208.0 400 0 kg/min T200 Temperatura de entrada líquido de refrigeración 25.0 OC Q200 Flujo de calor en el condensador 307.9 - - kW T201 Temperatura de salida líquido de refrigeración 46.1 - - OC Tabla 1.1 Variables del sistema * kg s. (kg de soluto), kg d. (kg de disolución) 4 5 2 MODELO DINÁMICO Para obtener el modelo dinámico se usan principalmente las ecuaciones de conservación de energía y de masa, así como las ecuaciones de transferencia de calor y las ecuaciones de estado termodinámicas. Se toma como referencia base el módelo de Newell and Lee, al que se añaden diversas modificaciones en base al modelado de plantas químicas [7]. 2.1 Balance de masa y energía. Evaporador 2.1.1 Modelo del intercambiador propuesto por Newell & Lee El modelo propuesto por Newell and Lee [6] estipula las siguientes consideraciones: La temperatura del líquido es: T2=0.5616P2+0.3126X2+48.43 (2–1) La fórmula anterior proviene de la linealización de la curva de equilibrio temperatura-presión de líquido saturado del agua alrededor del punto nominal de funcionamiento, e incluye un término que aproxima el punto de ebullición debido a la presencia del soluto (aumento ebulloscópico). Dicho término se obtiene linealizando las gráficas de Duhring para un líquido en el que se encuentra disuelto NaCl, que es la sal usada en este modelo. La dinámica del balance de energía se considera muy rápida, por lo que: F4= (Q100−F1·Cp(T2−T1))/λ (2–2) Cp es la capacidad calorífica de la disolución y se asume constante e igual a 0.07 kW/K(kg/min) (4.2 KJ/kgK para agua saturada a 80 ºC) y 𝜆 es el calor latente de vaporización de la disolución y se asume constante igual a 38.5 kW/(kg/min) (2304.11 KJ/kg a 50.5 KPa ). 2.1.2 Modelo final del evaporador Las ecuaciones de Newell y Lee conducen a un modelo demasiado simplificado que considera una circunstancia inadmisible en este proyecto: La temperatura del líquido en el evaporador es la temperatura final T2 en el mismo instante en el que entra en el inicio de los tubos. Es decir, el modelo de Newell y Lee considera que el intercambiador de calor de tubos largos actúa como un tanque homogéneo simplificado de mezcla perfecta. Para solucionar esto, se propone un modelado diferente, en el que el evaporador se divide en “N” compartimentos secuenciales, contínuos e iguales en tamaño. Cada compartimento actúa como un tanque homogéneo, teniendo por tanto cada uno de ellos corrientes de entrada y de salida. Puede subdividirse cada corriente en una fracción líquida y otra vapor, pero para mayor claridad, se tratará como si hubiese dos corrientes (líquida y vapor) tanto en la entrada como en la salida de los compartimentos. La temperatura en el interior de cada compartimento es, por tanto, igual a la temperatura de las corrientes de salida de ese compartimento. La siguiente Figura aclara la idea expuesta: 6 Figura 2.1 Esquema subdivisión el evaporador La entrada del evaporador, F2, es la suma de la corriente de alimentación F1 y la corriente de recirculación F3. Planteando los balances de energía y de masa de cada subdivisión se obtiene el siguiente sistema de ecuaciones, para n={1,…,N}: M NCpLdT2,n dt =(F2,n−1T2,n−1−F2,nT2,n)CpL+(F4,n−1T2,n−1−F4,nT2,n)Cpv −(F4,n−F4,n−1)Hv(T2,n)+hintAint(Tw,n−T2,n)/N (2–3) M NdX2,n dt =F2,n−1X2,n−1−F2,nX2,n (2–4) T2,n=0.5616P2+0.3136X2,n+48.43 (2–5) F2,n=(F2,n−1+F4,n−1)−F4,n (2–6) MwCpwdTw,n dt =(Q100−hintAint(Tw,n−T2,n)) (2–7) Cuanto mayor sea el valor de “N”, mejor será el modelado del evaporador, pero mayor será la complejidad al aumentar el número de incógnitas y ecuaciones. En este proyecto se toma el valor N=3, pues se ha calculado previamente hasta N=5 y a partir de N=3 la mejoría es cada vez más pequeña, aumentando excesivamente la complejidad de los diagramas en Matlab®-Simulink. En cuanto a las constantes y variables introducidas en el anterior sistema de ecuaciones, a continuación se recogen sus descripciones y valores conocidos:  𝐌: Masa de líquido en el evaporador. Se asume constante e igual a 20 kg.  𝐂𝐩𝐋: Capacidad calorífica de la disolución líquida. Se asume constante e igual al Cp del agua líquida (0.08 kW/K(kg/min)), según las tablas de líquido-gas saturado del agua [4].  𝐂𝐩𝐯: Capacidad calorífica del vapor de agua. Se asume constante e igual a 0.025 kW/K(kg/min), 7 según las tablas de líquido-gas saturado del agua [4].  𝐓𝟐,𝐧−𝟏 y 𝐓𝟐,𝐧: Temperatura de las corrientes de entrada y de salida en cada compartimento, respectivamente. Son variables del sistema de ecuaciones y tienen unidades de K.  𝐅𝟐,𝐧−𝟏 y 𝐅𝟐,𝐧: Caudal másico (kg/min) del fluido líquido de entrada y de salida en cada compartimento, respectivamente. Son variables.  𝐅𝟒,𝐧−𝟏 y 𝐅𝟒,𝐧: Caudal másico (kg/min) del fluido gaseoso de entrada y de salida en cada compartimento, respectivamente. Son variables.  𝐍: número total de divisiones realizadas en el evaporador (3 en este caso), n={1,…,N}.  𝐇𝐯(𝐓𝟐,𝐧): Entalpía del vapor evaluada a la temperatura T2,n. Unidades: kW·min/kg. Teniendo en cuenta que Hv(40 C o)=Hv(n=1)=42.892 kWmin/kg y que Hv(84.6 C o)=Hv(n=3)= 44.178 kWmin/kg [4], se linealiza entre los dos puntos para obtener una expresión lineal de la entalpía en función de la porción del intercambiador que se esté analizando (de un total de N porciones): Hv(n)−Hv(0)=Hv′(n)·(n−n(0))→ Hv(n)=42.892+44.178−42.892 3−1 (n−1)=42.892+0.643(n−1) (2–8)  𝐗𝟐,𝐧−𝟏 y 𝐗𝟐,𝐧: Concentración de la disolución a la entrada y a la salida de la porción, respectivamente. Unidades medidas en % (gramos de soluto entre gramos de disolución).  𝐏𝟐: Presión de operación en el interior del evaporador (en kPa). La ecuación (2–7) proviene del estudio de la transferencia de calor entre el vapor de entrada en la cámara y el líquido que circula en el interior de los tubos, según se detalla en el apartado siguiente: 2.1.2.1 Transferencia de calor en los tubos interiores del intercambiador de calor El vapor que entra en el intercambiador proporciona un suministro de calor que calienta la pared exterior de los tubos interiores. Este calentamiento produce un cambio de temperatura en la pared exterior e interior de los tubos, que produce un aumento de la temperatura del fluido que circula por el interior de los tubos (ver Figura 2.2) [4]. El calor del vapor de entrada (Q100) se invierte, por tanto, en calentar la pared de los tubos (MwCpwdTw,n dt ) y el fluido interior (Q=UAi(T100−T2)=1 (KL)−1+(hintAint)−1(Tw1,n− T2,n)=KL(Tw1−Tw2)=hintAint(Tw2,n−T2,n)), como se ve reflejado en la euación (2–7). Figura 2.2 Sección longitudinal interior de una de las mitades simétricas del intercambiador de calor  𝐐𝟏𝟎𝟎: Ver Tabla 1.1. 8  𝐓𝐰,𝐧: Temperatura de la pared interior de los tubos del intercambiador en la porción “n”. Variable medida en K. Su valor nominal inicial se estima Tw,n 0= 88 C o. En Figura 2.2 se usa el subíndice Tw2,n para esta variable para distinguir adecuadamente la pared exterior de la interior. Desde este punto en adelante se utiliza Tw,n para dicha variable a fin de simplificar los índices.  𝐡𝐢𝐧𝐭𝐀𝐢𝐧𝐭: Producto de la constante de convección entre la pared interior de los tubos y la corriente interior a los mismos por el área interior de los tubos. Constante e igual a 100.11 kW/K. Dicho valor proviene de la ecuación evaluada en torno al punto nominal de (2–7): p.n.→Q100 0=hintAint(Tw,n 0−T2,n 0)→hintAint=F100 0·λs Tw,n 0−T2,n 0=9.3·36.6 88−84.6=100.11 kW/K (2–9)  𝐌𝐰: Masa de los tubos. Esta constante es un valor intrínseco del intercambiador y depende lógicamente del material y el tamaño del mismo. Como estimación se toma el valor Mw=200 kg.  𝐂𝐩𝐰: Capacidad calorífica de los tubos del intercambiador. Depende del material de los mismos. Se estima constante e igual a 0.02 kW/K(kg/min) [4]. De esta manera, quedan definidas todos las parámetros y variables del sistema de ecuaciones. Para analizar los grados de libertad más adelante, hay que tener en cuenta que el primer compartimento y el último en los que se divide el intercambiador poseen las particularidades indicadas en la Tabla 2.1. n=1 Aclaración n=3 Aclaración 𝐅𝟐,𝟎=𝐅𝟏+𝐅𝟑′ La entrada de la primera porción es el líquido de alimentación 𝐅𝟐,𝟑=𝐅𝟐 La salida de la última porción es igual a la salida total del intercambiador 𝐅𝟒,𝟎=𝟎 No hay vapor en la entrada del intercambiador 𝐅𝟒,𝟑=𝐅𝟒 La salida de la última porción es igual a la salida total del intercambiador 𝐓𝟐,𝟎=𝐂𝐩(𝐅𝟏𝐓𝟏+𝐅𝟑′𝐓𝟐,𝐍 𝟎) (𝐅𝟏+𝐅𝟑′)𝐂𝐩 Temperatura de entrada del líquido de alimentación 𝐓𝟐,𝟑=𝐓𝟐 Temperatura del líquido de salida 𝐗𝟐,𝟎=𝐗𝟏𝐅𝟏+𝐗𝟐,𝐍 𝟎𝐅𝟑′ 𝐅𝟏+𝐅𝟑′ Concentración del líquido de alimentación 𝐗𝟐,𝟑=𝐗𝟐 Concentración del líquido de salida Tabla 2.1 Particularidades primera y última porción El calor transferido al líquido de proceso también se puede expresar en función de la temperatura del vapor en el interior de la cámara: Q100,n=UA1(T100−T2,n)→Q100=∑Q100,n n (2–10) donde UA1 es el coeficiente global de transferencia de calor entre el vapor de la cámara y el líquido de recirculación de los tubos. Evaluando la ecuación anterior en el punto nominal de funcionamiento se puede obtener un valor aproximado de UA1: Q100 0=UA1(T100 0−T20)→UA1=339 119.9−84.6=9.6 kW/K (2–11) 9 La presión de vapor P100 es una variable que determina la temperatura del vapor en condiciones de saturación. Se puede obtener una ecuación que relaciona la temperatura del vapor con su presión aproximando la curva de equilibrio temperatura-presión de vapor saturado por su linealización en el punto nominal de funcionamiento: T100=0.1538P100+90.0 (2–12) 2.2 Balance de masa y energía. Separador. La Figura 2.3 muestra las variables que conforman la dinámica del separador, la bomba de impulsión y la salida de producto: Figura 2.3 Variables separador, bomba y producto El balance de masa realizado sobre el total del líquido de proceso del sistema se puede formular como sigue, efectuado en el separador: ρAdL2 dt=(F2+F3)−(F2′+F3′) (2–13) El producto ρA (ρ es la densidad del líquido y A el área transversal del separador) se asume como constante de valor 20 kg/m. Como se asume que la variación de masa con respecto al tiempo en el evaporador es nula, se obtiene: dM dt=(F3′+F1)−(F2+F3+F4)=0→F2+F3+F4=F3′+F1 (2–14) Sustituyendo la ecuación anterior en (2–13) se llega al siguiente resultado: ρAdL2 dt=F1−F4−F2′ (2–15) Desarrollando el balance de masa de la sal disuelta en el separador separador se obtiene: ρAd(L2X2′) dt =(F2+F3)X2−(F2′+F3′)X2′→ ρAL2dX2′ dt+ρAX2′dL2 dt=(F2+F3)X2−(F2′+F3′)X2′ (2–16) Finalmente, al sustituir resulta: 16 Ecuación lineal alrededor del punto de funcionamiento (dominio de Laplace) ec. no lineal Elemento del evaporador 𝛒𝐀·𝐬𝐋𝟐=𝐅𝟏−𝐅𝟒−𝐅𝟐′ (2–15) separador 𝐗𝟐′=(𝛒𝐀𝐋𝟐 𝟎 (𝐅𝟏𝟎+𝐅𝟑𝟎−𝐅𝟒𝟎)𝐬+𝟏)−𝟏𝐗𝟐 (2–17) separador 𝐂·𝐬𝐏𝟐=𝐅𝟒−𝐅𝟓 (2–18) separador 𝐓𝟒=𝟎.𝟓𝟎𝟕𝐏𝟐 (2–19) separador 𝐌 𝐍𝐂𝐩𝐋𝐬𝐓𝟐,𝐧−𝟏=(𝐓𝟐,𝐧−𝟏 𝟎𝐂𝐩𝐯+𝐇𝐯(𝐓𝟐,𝐧 𝟎))𝐅𝟒,𝐧−𝟏−(𝐓𝟐,𝐧 𝟎𝐂𝐩𝐯+𝐇𝐯(𝐓𝟐,𝐧 𝟎))𝐅𝟒,𝐧 +(𝐓𝟐,𝐧−𝟏 𝟎𝐂𝐩𝐋)𝐅𝟐,𝐧−𝟏−(𝐓𝟐,𝐧 𝟎𝐂𝐩𝐋)𝐅𝟐,𝐧+(𝐅𝟐,𝐧−𝟏 𝟎𝐂𝐩𝐋+𝐅𝟒,𝐧−𝟏 𝟎𝐂𝐩𝐯)𝐓𝟐,𝐧−𝟏 −(𝐅𝟐,𝐧 𝟎𝐂𝐩𝐋+𝐅𝟒,𝐧 𝟎𝐂𝐩𝐯−𝐡𝐢𝐧𝐭𝐀𝐢𝐧𝐭 𝐍)𝐓𝟐,𝐧+ (𝐡𝐢𝐧𝐭𝐀𝐢𝐧𝐭 𝐍)𝐓𝐰,𝐧 (2–3) Intercambiador 𝐌 𝐍·𝐬𝐗𝟐,𝐧=(𝐗𝟐,𝐧−𝟏 𝟎)𝐅𝟐,𝐧−𝟏+(𝐅𝟐,𝐧−𝟏 𝟎)𝐗𝟐,𝐧−𝟏−(𝐗𝟐,𝐧 𝟎)𝐅𝟐,𝐧−(𝐅𝟐,𝐧 𝟎)𝐗𝟐,𝐧 (2–4) Intercambiador 𝐓𝟐,𝐧=𝟎.𝟓𝟔𝟏𝟔𝐏𝟐+𝟎.𝟑𝟏𝟑𝟔𝐗𝟐,𝐧 (2–5) Intercambiador 𝐅𝟐,𝐧=𝐅𝟐,𝐧−𝟏+𝐅𝟒,𝐧−𝟏−𝐅𝟒,𝐧 (2–6) Intercambiador 𝐌𝐰𝐂𝐩𝐰·𝐬𝐓𝐰,𝐧=𝐐𝟏𝟎𝟎−𝐡𝐢𝐧𝐭𝐀𝐢𝐧𝐭𝐓𝐰,𝐧+𝐡𝐢𝐧𝐭𝐀𝐢𝐧𝐭𝐓𝟐,𝐧 (2–7) Intercambiador 𝐐𝟏𝟎𝟎=𝐔𝐀𝟏𝐓𝟏𝟎𝟎−𝐔𝐀𝟏𝐓𝟐 (2–10) Intercambiador 𝐓𝟏𝟎𝟎=𝟎.𝟏𝟓𝟑𝟖𝐏𝟏𝟎𝟎 (2–12) Intercambiador 𝐐𝟐𝟎𝟎=𝐔𝐀𝟐 𝟏+ 𝐔𝐀𝟐 𝟐𝐂𝐩𝐅𝟐𝟎𝟎 𝟎𝐓𝟒−𝐔𝐀𝟐 𝟏+ 𝐔𝐀𝟐 𝟐𝐂𝐩𝐅𝟐𝟎𝟎 𝟎𝐓𝟐𝟎𝟎+𝐔𝐀𝟐𝐐𝟐𝟎𝟎 𝟎 𝐅𝟐𝟎𝟎 𝟎(𝟐𝐂𝐩𝐅𝟐𝟎𝟎 𝟎+𝐔𝐀𝟐)𝐅𝟐𝟎𝟎 (2–22) Condensador 𝐓𝟐𝟎𝟏=𝐓𝟐𝟎𝟎+𝟏 𝐂𝐩𝐅𝟐𝟎𝟎 𝟎𝐐𝟐𝟎𝟎−𝐐𝟐𝟎𝟎 𝟎 𝐂𝐩𝐅𝟐𝟎𝟎 𝟎𝟐𝐅𝟐𝟎𝟎 (2–23) Condensador 𝐅𝟓=𝛌−𝟏𝐐𝟐𝟎𝟎 (2–24) Condensador 𝟏𝟖𝐩𝟏𝟎𝟎 𝟎 𝟎.𝟎𝟖𝟐𝐓𝟏𝟎𝟎 𝟎𝐂𝐩𝐯·𝐬𝐓𝟏𝟎𝟎=𝐔𝐀𝟏𝐓𝟐−𝐔𝐀𝟏𝐓𝟏𝟎𝟎+(𝐇𝐯(𝐓𝟏𝟎𝟎 𝟎)−𝐂𝐩𝐯𝐓𝟏𝟎𝟎 𝟎)𝐅𝟏𝟎𝟎 +(𝐇𝐯(𝐓𝟏𝟎𝟎)−𝐂𝐩𝐯𝐓𝟏𝟎𝟎)𝐅𝟏𝟎𝟎 𝟎−(𝐇𝐋(𝐓𝟏𝟎𝟎 𝟎)−𝐂𝐩𝐯𝐓𝟏𝟎𝟎 𝟎)𝐅𝐜 −(𝐇𝐋(𝐓𝟏𝟎𝟎)−𝐂𝐩𝐯𝐓𝟏𝟎𝟎)𝐅𝐜𝟎 (2–29) Cámara 17 𝐅𝟏𝟎𝟎=(𝟏𝟒.𝟐·𝟔𝟎𝐊𝐯𝐬√𝐏𝟎𝟎𝛒𝟎𝟎)𝐱𝐅𝟏𝟎𝟎+ ( 𝟏𝟒.𝟐·𝟔𝟎𝐊𝐯𝐬·𝐱𝐅𝟏𝟎𝟎 𝟎·𝛒𝟎𝟎 𝟐√𝐏𝟎𝟎𝛒𝟎𝟎 ) 𝐏𝟎 + ( 𝟏𝟒.𝟐·𝟔𝟎𝐊𝐯𝐬·𝐱𝐅𝟏𝟎𝟎 𝟎·𝐩𝟎𝟎 𝟐√𝐏𝟎𝟎𝛒𝟎𝟎 ) 𝛒𝟎 (2–31) Válvulas 𝐅𝐜=𝐊𝐯𝐬𝐜·𝟔𝟎·𝛒𝐚𝐠𝐮𝐚·𝐱𝐜 𝟐√𝐩𝟏𝟎𝟎 𝟎+𝛒𝐚𝐠𝐮𝐚𝐠𝐡 𝟏𝟎𝟓−𝟏𝐏𝟏𝟎𝟎 (2–33) Válvulas Tabla 3.1 Linealizaciones en torno al punto de funcionamiento Los puntos nominales de las variables, dependerán del evaporador que se esté usando. Para este proyecto, se usan los puntos nominales recogidos en la Tabla 1.1. Por otro lado, los puntos nominales de las variables del modelo del evaporador dividido en ‘N’ tanques homogéneos, se calculan con el modelo estático, anulando las derivadas de dichas ecuaciones y despejando los valores de equilibrio. El código utilizado para el cálculo de estos puntos de funcionamiento se pueden ver en Anexo A. n 𝐓𝟐,𝐧 ( 𝐂 𝐨) 𝐓𝐰,𝐧( 𝐂 𝐨) 𝐗𝟐,𝐧(%) 𝐅𝟒,𝐧(𝐤𝐠 𝐦𝐢𝐧) 𝐅𝟐,𝐧(𝐤𝐠 𝐦𝐢𝐧) 1 83.82 87.2 22.5 2.25 57.75 2 84.2 87.6 23.7 5.15 54.85 3 84.6 88 25 8 52 Tabla 3.2 Puntos de funcionamiento variables intermedias evaporador El sistema lineal proporciona una aproximación del modelo no lineal alrededor del punto nominal de funcionamiento. Estos resultados permiten obtener una idea aproximada de la relación entre las variables. De esta manera, se puede observar qué pasaría frente a cambios en ciertos puntos de funcionamiento (diferente comportamiento según el evaporador usado), y el diferente peso que tienen unas variables sobre otras. Por ejemplo, en la ecuación (2–23) linealizando en el punto nominal: 𝐓𝟐𝟎𝟏=1·𝐓𝟐𝟎𝟎+0.069·𝐐𝟐𝟎𝟎−0.1·𝐅𝟐𝟎𝟎 (3–5) De donde se puede concluir que la variable que más afecta a la temperatura de salida del refrigerante, es la temperatura de entrada del mismo, seguido del caudal de refrigerante, ante cambios unitarios en las variables con las unidades elegidas. Aunque con estas ecuaciones se pueden obtener las funciones de transferencia del sistema, el proceso es tedioso, por lo que más adelante (capítulo 5), se identificarán estas funciones G(s) por otro método mas rápido y comúnmente utilizado si se dispone del sistema físico en funcionamiento o de un simulador. 18 19 4 MODELO EN SIMULINK-MATLAB La implementación del modelo dinámico se ha realizado con la utilidad SIMULINK de MATLAB®. Esta herramienta permite crear un diagrama de bloques que integra las ecuaciones del modelo. En dicho diagrama, se pueden cambiar manualmente los valores de los puntos de equilibrio, entradas y perturbaciones para observar mediante un entorno gráfico el efecto de estos cambios en las salidas seleccionadas. Esta herramienta ofrece por tanto un método rápido e intuitivo para analizar y observar la dinámica del evaporador. Se ha desarrollado el modelo de Simulink que implementa el modelo no lineal del evaporador, denominado “modelo_nolineal_final.slx”. Para abrir el archivo, se realiza la siguiente secuencia: 1. Abrir la aplicación “Matlab”, escribir “Simulink” y pulsar el botón “enter” para iniciar la aplicación. El tiempo de carga puede ser de un par de minutos. Es necesario tener los archivos “parametros.m” y “modelo_nolineal_final.slx” en el directorio de trabajo de Matlab (“Current Folder”). 2. Ejecutar el archivo “parámetros.m” y pulsar sobre el botón “Run”. Este archivo carga todos los puntos nominales de funcionamiento, así como las constantes del sistema y otros parámetros. 3. Abrir el archivo “modelo_nolineal_final.slx”. Aquí se puede observar el diagrama completo ya creado del modelo dinámico. 4. Para ejecutar la simulación, sólo es necesario especificar el periodo de tiempo en la parte superior de la ventana y pulsar el botón “Run”. Al hacer esto, y transcurridos unos segundos, se habrá realizado la simulación para el tiempo especificado (en minutos, según las unidades utilizadas en el modelo) y se puede acceder a las ventanas “Scope” para observar la respuesta de las variables de salida. 5. Para realizar diferentes simulaciones, sólo es necesario cambiar los parámetros deseados, que son, típicamente, las varaibles manipulables y las perturbaciones, y volver a pulsar el botón “Run”. Por defecto, la configuración inicial del archivo dispone de todos los parámetros en sus puntos nominales de funcionamiento. En la Figura 4.1 se puede observar la pantalla principal del diagrama de bloques del modelo no lineal: Figura 4.1 Pantalla principal modelo no lineal Simulink  A la izquierda de la pantalla hay una columna de escalones. Estos escalones representan los puntos de funcionamiento de las 3 variables manipulables y las 5 perturbaciones del sistema. Las ganancias conectadas a continuación son el valor por el que se multiplica el punto de funcionamiento, dando como resultado las variaciones de estas variables, modificables por el usuario. 20  Hay tres bloques principales: “EVAPORADOR(Intercambiador de calor)”, “SEPARADOR” y “CONDENSADOR”. Pulsando sobre cada uno de ellos se pueden observar varios niveles de bloques que comprenden todas las ecuaciones desarrolladas en el capítulo 2 de este proyecto. Para ver dichos niveles se puede abrir directamente el archivo de Simulink adjunto al proyecto o acudiendo al ANEXO B donde se recogen también las imágenes de estos diagramas.  Diversos bloques tienen un color específico. Estos colores representan la diferente funcionalidad de los bloques, en concreto: o Azul: Variables manipulables del sistema, típicamente escalones. o Amarillo: Variables de perturbación del sistema, típicamente escalones. o Rojo: Salidas del sistema, típicamente bloques “Scope” que permiten observar su evolución temporal. o Verde: Variables intermedias. Para la simulación, el usuario sólo tiene que cambiar las ganancias conectadas a los escalones para modificar el valor de estas variables con respecto al punto de funcionamiento. Es necesario tener en cuenta los límites establecidos en la Tabla 1.1 A modo de ejemplo, para observar el efecto de una variación en xF100 sobre X2’, se introduce el valor 1.1 en “Gain” (incremento del 10% de xF100 con respecto a su punto de funcionamiento), obteniéndose el siguiente resultado al pulsar sobre el bloque Scope X2’: Figura 4.2 Ejemplo de respuesta de simulación En los capítulos siguientes se utiliza este modelo para las simulaciones, así como la modificación pertinente del diagrama para la implementación de control, como se indicará en los correspondientes apartados y que viene recogida también en el ANEXO B. A continuación se muestra y se comenta el efecto de cada variable manipulabre sobre cada variable de control: 21  Un cambio en la abertura de la válvula XF100 provoca un cambio de tipo “2o orden” en la concentración final de producto. Esto se debe al cambio de temperatura y, por tanto, al cambio de calor transferido que se produce al líquido de los tubos al introducir una mayor o menor cantidad de flujo vapor en la cámara del intercambiador.  Un cambio en la abertura de la válvula XF100 provoca un cambio en la cantidad de vapor producido en el evaporador, que afecta directamente a la presión P2. En este caso la dinámica es más lenta que la evolución de X2’, y es de primer orden.  El nivel del separador se ve directamente afectado por el aumento o disminución de cantidad de vapor producida por el cambio en la abertura de XF100, presentando una respuesta de tipo integrador. Figura 4.3 Respuestas variables de control frente a variación en XF100 Figura 4.4 Respuestas variables de control frente a variación en F200  Al producirse una variación en F200, X2’ reacciona como un sistema de primer orden. En la Iustración 4.4 se refleja la evolución producida en dicha variable de control. Al aumentar la refrigeración, el sistema evaporador es capaz de enfriar un mayor caudal de vapor, lo que permite una mayor evaporación y por tanto una mayor concentración de producto. De forma opuesta es la respuesta del sistema al disminuir la refrigeración.  La presión del evaporador dismiye al aumentar la refrigeración, y aumenta al bajar la corriente de refrigeración, presentando una respuesta que se ajusta a un sistema de primer orden. 22  L2 varía en relación a X2 con una respuesta de tipo integrador. Al cambiar la corriente refrigerante, el mayor o menor caudal de vapor producido afecta directamente a la cantidad de líquido en el separador. Figura 4.5 Respuestas variables de control frente a variación en F2  Frente a cambios en la corriente de producto, la variación de la concentración de producto y la variación de la presión del evaporador. Esto se debe a las consideraciones tomadas con respecto a F3 (ver 2.6 Consideraciones).  El nivel de líquido del separador es directamente afectado por el caudal de salida de producto. Así, al aumentar dicho caudal, disminuye la cantidad de líquido en el separador, y viceversa. La respuesta que se obtiene es de tipo integrador. Figura 4.6 Respuestas variables de control frente a variación en P0 23  Los cambios en la presión de vapor proveniente de una caldera afectan a las variables de control de manera similar a los cambios producidos por variaciones en la abertura de la válvula. Esto se debe a que el flujo de entrada a la cámara de evaporador F100 depende proporcionalmente tanto de XF100 como de P0. Figura 4.7 Respuestas variables de control frente a variación en F1  Los cambios en la corriente de alimentación son ilustrativos de la no linealidad del sistema. Así, un aumento de dicho caudal produce una caída de concentración de producto mucho más rápida que la subida de concentración provocada por una disminución de caudal de alimentación (ambos con respuesta de 2o orden). Dicho comportamiento no lineal también se ve reflejado en la variación de la presión del evaporador frente a cambios en f1.  La relación entre el nivel del separador y el caudal de alimentación es más directa. Aumentar el caudal de entrada implica, según el balance másico, un aumento en el nivel de líquido del separador, y viceversa. Este cambio responde a una función de tipo integrador. Figura 4.8 Respuestas variables de control frente a variación en X1 24  La concentración del líquido de alimentación afecta directamente a la concentración final de producto. X2’ evoluciona como un sistema de 2 orden frente a cambios en X1’.  La presión del evaporador responde de manera inversa ante los cambios en la concentración de entrada. Un aumento de X1 produce una disminución en P2 y viceversa. Esto se debe a que una mayor mayor concentración de sal en el líquido de los tubos provoca que frente al a la misma cantidad de calor transferido se produzca una menor cantidad de vapor de agua, lo que disminuye la presión del sistema.  El nivel del separador está relacionado, con la cantidad de vapor de agua producido y la concentración final de producto, presentando una respuesta de tipo integrador. Figura 4.9 Respuestas variables de control frente a variación en T1  Los cambios positivos en la temperatura del líquido de alimentación producen un aumento de la presión del evaporador, ya que el sistema pierde capacidad de refrigeración al aumentar la temperatura de la mezcla sin variar el caudal de entrada en el condensador. La respuesta de P2 frente a T1 es de tipo lineal.  Frente a cambios T1, X2’ reacciona con una respuesta de segundo orden proporcional a la entrada, y L2 reacciona con una respuesta de tipo integrador, inversamente proporcional con respecto a cambios en T1. 25 Figura 4.10 Respuestas variables de control frente a variación en T200  Los cambios en la temperatura del líquido refrigerante afectan directamente a la presión del evaporador. Cuanto mayor sea la temperatura de refrigeración, menos vapor se podrá condensar y por tanto mayor será la presión del evaporador (y viceversa). La respuesta de P2 frente a cambios en T200 es de primer orden.  El nivel del separador cambio de manera directamente proporcional frente a cambios en T200. La respuesta de L2 se asemeja a un sistema integrador con retardo.  La concentración de producto final es inversamente proporcional a los cambios en la temperatura del refrigerante. El sistema no puede condensar la cantidad de vapor del sistema funcionando en torno a su punto nominal, por lo que la concentración final de producto disminuye ante un aumento en la temperatura del refrigerante. La respuesta es de segundo orden. 32 Figura 5.7 Respuesta sistema nolineal-lineal. X2’ frente a escalón amplitud ± 10% en F200, y ajuste Figura 5.8 Respuesta sistema nolineal-lineal. P2 frente a escalón amplitud ± 10% en F200, y ajuste Figura 5.9 Respuesta sistema nolineal-lineal. L2 frente a escalón amplitud ± 10% en F200, y ajuste 33 Figura 5.10 Respuesta sistema nolineal-lineal. L2 frente a escalón amplitud ± 10% en F2, y ajuste En todas las simulaciones se repite el mismo patrón:  Color Azul. Las respuestas lineales correspondientes a G+(s) se ajustan mucho mejor a la respuesta no lineal frente a escalones positivos en la variable manipulable del 10% de amplitud. Sin embargo, frente a entradas de -10% de amplitud, la bondad del ajuste es peor que las otras 2 funciones de transferencia. Esto se debe a cuestiones de diseño, pues G+(s) se diseña en función de entradas de +10% de amplitud con respecto al punto de funcionamiento.  Color Verde. Las respuestas lineales correspondientes a G−(s) se ajustan mucho mejor a la respuesta no lineal frente a escalones en la variable manipulable del -10% de amplitud. Al contrario que en el caso anterior, es frente a escalones positivos en la variable manipulable donde las respuestas lineales de G−(s) tienen un peor ajuste. Se debe de nuevo al diseño de la función de transferencia.  Color Rojo. Las respuestas lineales correspondientes a G(s) ofrecen una bondad de ajuste intermedia entre G+(s) y G−(s). Utilizar esta función es la opción más adecuada frente a simulaciones en las que se cambien las variables manipulables tanto con escalones positivos como negativos con respecto a su punto de funcionamiento. Cuanto mayor sea la amplitud de entrada con respecto al punto de funcionamiento de la variable manipulable, peor será la bondad del ajuste del sistema lineal. Esto se debe a la no linealidad del sistema. En este proyecto se ha considerado en las simulaciones escalones de hasta el 20% con respecto al punto nominal de las variables manipulables en las simulaciones, pues en estos casos los valores alcanzados en régimen permanente de las variables controladas lineales tienen un error con respecto a las no lineales menor siempre al 25%. La pérdida de bondad de simulación conforme mayor sean estas amplitudes se traduce en una peor efectividad de los controladores diseñados a la hora de seguir referencias y rechazar perturbaciones. 5.2 Identificación Gd(s) Siguiendo el mismo procedimiento usado en el apartado anterior, se define la matriz Gd(s) : [X2′(s) P2(s) L2(s)]=Gd(s)· [ P0(𝑠) F1(s) X1(𝑠) T1(s) T200(s) ] (5–6) La matriz dada por la ecuación (5–7) representa las funciones de transferencias frente a perturbaciones de escalón positivo del 10% con respecto al punto nominal de funcionamiento. 34 Gd+(s)= [ 23.275 4.269s+1e−0.3s 19.07 34.84s+1e−0.02s −0.0671 s −6.729 2.675s+1e−0.3s −0.0263 1.287𝑠+1𝑒−0.3𝑠 0.0504 s 4.376 4.489s+1e−0.9s −0.579 45.0s+1e−0.11s 0.00233 se−s 0.21 2.67s+1e−0.8s 0.21 43.45s+1 −0.000841 s −0.37 42.34s+1e−0.3s 1.55 44.0𝑠+1 0.000883 se−2s ] (5–7) Utilizando escalones de -10% de amplitud con respecto al punto de funcionamiento, se obtiene la matriz (5–8): Gd−(s)= [ 13.133 2.785s+1e−0.3s 9.093 50.24s+1 −0.0716 s −16.65 5.706s+1e−0.3s 150(1+0.017𝑠) (1000𝑠+1) 0.0478 s e−0.02𝑠 4.343 4.604s+1e−0.9s −0.7 56.1s+1e−0.3s 0.00250 se−s 0.195 2.79s+1e−0.5s 0.198 41.92s+1 −0.000831 s −0.416 45.63s+1e−0.4s 1.53 43.5s+1e−0.024s 0.000911 se−2s ] (5–8) Finalmente se toma Gd(s) como media aritmética de Gd+(s) y Gd-(s) : Gd(s)= [ 18.204 3.527s+1e−0.3s 14.08 42.54s+1e−0.01s −0.0694 s −11.69 4.188s+1e−0.3s −0.0263 1.287𝑠+1𝑒−0.3𝑠 0.0491 s 4.4 4.54s+1e−0.9s −0.64 50.6s+1e−0.21s 0.00242 se−s 0.203 2.73s+1e−0.65s 0.204 42.69s+1 −0.000836 s −0.393 44.0s+1e−0.35s 1.54 43.8s+1e−0.012s 0.000897 se−2s ] (5–9) Las funciones de transferencia anteriores se han calculado mediante la utilidad System Ident Toolbox de Matlab. A continuación se muestran las comparaciones entre las respuestas de las variables de control del sistema no lineal y de las funciones lineales, así como la bondad de los ajustes, frente a variaciones en las variables de perturbación (a la izquierda frente a cambios del +10% y a la derecha frente a cambios del -10%). El color negro representa la respuesta de las variables de control del sistema no lineal frente a escalones en la perturbación de amplitud ± 10%. El color azul representa la respuesta de las variables de control del sistema lineal Gd+(s) (5–7) frente a escalones en la perturbación de amplitud ± 10%, el color verde representa la respuesta de las variables de control del sistema lineal Gd−(s) (5–4) frente a escalones en la perturbación de amplitud ± 10%, y el color rojo representa la respuesta de las variables de control del sistema lineal Gd(s) (5–9) frente a escalones en la perturbación de amplitud ± 10%. Figura 5.11 Respuesta sistema nolineal-lineal. X2’ frente a escalón amplitud ± 10% en P0, y ajuste 35 Figura 5.12 Respuesta sistema nolineal-lineal. P2 frente a escalón amplitud ± 10% en P0, y ajuste Figura 5.13 Respuesta sistema nolineal-lineal. L2 frente a escalón amplitud ± 10% en P0, y ajuste Figura 5.14 Respuesta sistema nolineal-lineal. X2’ frente a escalón amplitud ± 10% en F1, y ajuste 36 Figura 5.15 Respuesta sistema nolineal-lineal. P2 frente a escalón amplitud ± 10% en F1, y ajuste Figura 5.16 Respuesta sistema nolineal-lineal. L2 frente a escalón amplitud ± 10% en F1, y ajuste Figura 5.17 Respuesta sistema nolineal-lineal. X2’ frente a escalón amplitud ± 10% en X1, y ajuste 37 Figura 5.18 Respuesta sistema nolineal-lineal. P2 frente a escalón amplitud ± 10% en X1, y ajuste Figura 5.19 Respuesta sistema nolineal-lineal. L2 frente a escalón amplitud ± 10% en X1, y ajuste Figura 5.20 Respuesta sistema nolineal-lineal. X2’ frente a escalón amplitud ± 10% en T1, y ajuste 38 Figura 5.21 Respuesta sistema nolineal-lineal. P2 frente a escalón amplitud ± 10% en T1, y ajuste Figura 5.22 Respuesta sistema nolineal-lineal. L2 frente a escalón amplitud ± 10% en T1, y ajuste Figura 5.23 Respuesta sistema nolineal-lineal. X2’ frente a escalón amplitud ± 10% en T200, y ajuste 39 Figura 5.24 Respuesta sistema nolineal-lineal. P2 frente a escalón amplitud ± 10% en T200, y ajuste Figura 5.25 Respuesta sistema nolineal-lineal. L2 frente a escalón amplitud ± 10% en T200, y ajuste 5.3 Elección emparejamientos SISO Los emparejamientos de salida-entradas más favorables para un sistema SISO se pueden analizar con más profundidad usando la matriz DRGA (dynamic relative gain array) [10]. La matriz DRGA es similar a la RGA (matriz de ganancias relativas estáticas), con la particularidad de que para su cálculo se utilizan las funciones de transferencia completas, no solo las ganancias estáticas. El análisis de esta matriz permite interpretar qué emparejamientos son más favorables para realizar un control monovariable, información muy útil sobre las interacciones entre variables del sistema. Con estos datos se pueden escoger las parejas entrada-salida que más interactúan entre sí y menos con las demás. Se define DRGA: DRGA(G(s))=Λ(s)≜G(s)× (G(s)−1)T (5–10) Siendo G(s) la función de transferencia de la variable de control frente a la variable manipulable del sistema. ‘x’ representa el producto elemento por elemento de las dos matrices. Cada elemento ‘λij’ de la DRGA representa: λij=lim s→0(Gij(s))todos los lazos abiertos Gij(s))todos los lazos cerrados) (5–11) 40 Aunque la ecuación (5–10) es correcta, implica cálculo matricial cuyos elementos son funciones de Laplace. Operar este sistema matricial es complicado e induce a error incluso para programas informáticos. Sin embargo, una simplificación del cálculo fue introducida en el artículo publicado en 2010 por Wuhua Hu, WenJian Cai y Gaoxi Xiao. [12]. Dicha simplificación consiste en cambiar la matriz G(s) por la matriz de ganancias G(0), pero tanto los integradores y diferenciadores puros que hubiere se sustituyen por la variable “ε”, que tiene un valor pequeño (típicamente un par de ordenes inferior a la ganancia más baja de G(s)). El resultado es el mismo que la RGA, pero con la ventaja de que se pueden considerar los integradores y diferenciadores en el cálculo, aspecto que no se podía desarrollar anteriormente. Aplicando al caso anterior: G(0,ε)=[171 159.6 0 0.02 −0.07923 0 −0.6875 ε−5.84·10−5 ε−0.0500 ε] (5–12) Λ=RGA([ 171 159.6 0 0.02 −0.07923 0 −0.6875 𝜀−5.84·10−5 𝜀−0.0500 𝜀])𝜀→1·10−7= [171 159.6 0 0.02 −0.07923 0 −0.6875 𝜀−5.84·10−5 𝜀−0.0500 𝜀]X ( [171 159.6 0 0.02 −0.07923 0 −0.6875 𝜀−5.84·10−5 𝜀−0.0500 𝜀]−1 ) 𝑇= [171 159.6 0 0.02 −0.07923 0 −0.6875 𝜀−5.84·10−5 𝜀−0.0500 𝜀] 𝑋 [0.004733 0.001195 −0.06508 9.534 −10.21 −131.1 0 0 −20𝜀]= XF100 F200 F2′ [0.809 0.191 0 0.191 0.809 0 0 0 1] X2′ P2 L2 (5–13) Sea yi la salida del sistema y mj la entrada (conjunto proceso-válvula)[7] :  λij=1. Significa que la ganancia estática de yi respecto de mj es independiente de si los lazos están abiertos o cerrados. Se trata de un caso in interacción. En este caso, el emparejamiento salida-entrada es el ideal para el lazo de control.  λij=0. Este caso no indica si el correspondiente emparejamiento es factible o no. El control depende enteramente de los demás lazos.  0<λij<1. El emparejamiento yi-mj interactúa con los demás lazos. Esta interacción es mayor cuanto más se acerque el término a 0.5, y menor cuanto más se acerque a 1.  λij<0. El signo de la ganancia estática de yi respecto de mj depende de si los lazos están abiertos o cerrados. La ganancia cambia de signo cuando los demás lazos están cerrados.  λij>1. La ganancia estática de yi respecto de mj es menor con los restantes lazos cerrados que cuando están abiertos. 41  λij=∞. La ganancia estática de yi respecto de mj con los otros lazos cerrados es 0. Atendiendo a dichas propiedades de la DRGA, las parejas más adecuadas para el modelo equivalente SISO son las siguientes, puesto que son las que corresponden a los elementos de la matriz más cercanos a la unidad: X2′→XF100 P2→F200 L2→F2′ Por otro lado, cabe destacar que aunque se eligan estos lazos, cualquier variación de las variables XF100 y F200 afectarán inevitablemente a las salidas P2, X2’, lo que es significativo de la naturaleza multivariable del sistema. 5.4 Control SISO equivalente con controles por realimentación El sistema de control que se va a utilizar es la realimentación negativa. Mediante el uso de controladores que actúan sobre las entradas y que toman como referencia el error en la salida, se busca solventar los problemas de seguimiento de la referencia de salida y de rechazo de perturbaciones. El esquema típico de esta solución de control se puede observar en la Figura 5.26: Figura 5.26 Modelo de control realimentado El sistema realimentado consta de los siguientes elementos:  Yr(s): Referencia. Valor de la salida que se quiere obtener.  e(s): Error medido. Diferencia entre la referencia y la variable controlada medida.  Gc(s): Controlador. Elemento de control que recibe el error y proporciona la salida para el elemento final de control (en este caso una válvula) para corregir dicho error.  m(s). Entrada del elemento final de control.  Gv(s): Función de transferencia de la válvula.  u(s): Entrada del proceso (salida de la válvula).  Gp(s): Sistema. Función que modela el comportamiento del sistema físico a controlar.  Y(s): Salida del sistema. 48  Ts. Tiempo de subida. Es el tiempo que tarda la salida en alcanzar por primera vez el valor de la salida en régimen permanente.  Te5%. Tiempo de establecimiento al 5%. Es el tiempo que tarda la salida en alcanzar el régimen estacionario (con una desviación del 5% con respecto al régimen permanente). A continuación se muestra la respuesta de X2’ frente a una cadena de perturbaciones, manteniendo la referencia en el punto nominal. En concreto, en el minuto 20 se introduce un escalón de amplitud -10% en T1 y en el minuto 60 un escalón de amplitud +10% en la perturbación P0. Figura 5.31 Control realimentado X2’-XF100 frente a perturbaciones Como se puede observar, la perturbación P0 afecta mucho más a la sobreoscilación de X2’ que la perturbación T1. Esto cobra sentido físicamente pues una variación en la presión de la caldera provocará un cambio no deseado en el caudal de vapor de entrada que el controlador intentará corregir. Aunque en el régimen permanente la perturbación se termina rechazando al completo, la sobreoscilación provocada puede suponer un problema para los requerimientos de concentración de producto. Para ello, más adelante se introducirán otras técnicas para rechazar dichas perturbaciones con más efectividad. Para verificar la importancia del filtro derivativo, en la Figura 5.32 se ilustra el tipo de respuesta que se obtiene en la salida del sistema en caso de factor derivativo N=0.1 y N=20 frente a la misma referencia anterior, donde se puede observar la amplificación del ruido. La respuesta que se obtiene muestra el efecto negativo sobre la variable manipulable de un factor derivativo bajo, pues existen muchos picos de difícil seguimiento para la válvula real. Por otro lado, la salida X2’ tiene efectos de amplificación de ruido, indeseables en todo caso. 49 Figura 5.32 Respuesta de X2′% y xF100% frente a cambios en la referencia con diferentes filtros PID El código utilizado en Matlab para presentar las gráficas es el que aparece a continuación. El resto de códigos están recogidos en el Anexo A. Código 3.1 figure %a= introducir aquÌ la variable de Scope Simulink de la salida %b= introducir aquÌ la variable de Scope Simulink de la entrada ax1 = subplot(2,1,1); % top subplot ax2 = subplot(2,1,2); % bottom subplot plot(ax1,... a(:,1),a(:,3),'black--',... a(:,1),a(:,2),'r') grid(ax1,'on') title(ax1,'salida en % con controlador PID') ylabel(ax1,'X2¥(%)') xlabel(ax1,'time(minutos)') plot(ax2,... b(:,1),b(:,2),'b') grid(ax2,'on') title(ax2,'entrada en % con controlador PID') ylabel(ax2,'xF100 (%)') xlabel(ax2,'time(minutos)') 50 5.4.3 Control realimentado para el emparejamiento 𝐏𝟐−𝐅𝟐𝟎𝟎 GP%(s)=(P2%(s) F200%(s))=−0.3170 43.82s+1 (5–23) Aplicando las relaciones descritas en la Tabla 5.1, se calculan los parámetros de los controladores PI y PID (𝜏𝑐=𝜏/4): PI { 𝐤𝐜=1 −0.3170 43.82 43.82/4+0=−𝟏𝟐.𝟔𝟏𝟖 𝐭𝐢=min(43.82; 4·(4.823 4+0))=𝟒𝟑.𝟖𝟐 (5–24) PID { 𝐤𝐜=1 −0.31702·43.82+0 2·43.82 4+0=−𝟏𝟐.𝟔𝟏𝟖 𝐭𝐢=43.82+02=𝟒𝟑.𝟖𝟐 𝐭𝐝=43.82·0 2·43.82+0=𝟎 (5–25) Destaca que tanto el PI como el PID sean iguales según el diseño SIMC, debido a que la acción derivativa es nula. Ésto tiene una explicación: El sistema simplificado no tiene retardo, puesto que la respuesta no lineal original no presenta un crecimiento inicial suave (derivada ≠ 0). Esto indica por tanto que la respuesta no lineal presenta un comportamiento típico de primer orden, lo que hace innecesario el uso de un PID (segundo orden) para anular los polos. Con el polo único del PI se consigue este objetivo. Frente a esta situación, la elección es clara a favor del PI, pues la inclusión de la acción derivativa no sólo no mejoraría la respuesta, sino que incluiría problemas ya mencionados como la amplificación de ruido. El esquema de control es exactamente el mismo que el utilizado en la Figura 5.27, sustituyendo los correspondientes parámetros del controlador, y la pareja entrada-salida a controlar, así como las ganancias de ajuste de adimensionalización. La Figura 5.33 muestra la respuesta del sistema, con el controlador PI implementado, ante un cambio inicial del -10% en la referencia P2 y, posteriormente, a los 100 minutos, un cambio del +10% con respecto al punto nominal de funcionamiento. Este lazo es un claro ejemplo de un controlador de acción directa (véase la ganancia negativa). Frente a un aumento de la P2, el controlador reacciona proporcionando un mayor valor F200, y viceversa. 51 Figura 5.33 Respuesta de P2% y F200% frente a cambios en la referencia con PI La respuesta es notoriamente más lenta que en el emparejamiento X2-xF100, lo cual no es un problema puesto que P2 no es una variable candidata a sufrir exigentes seguimientos de la referencia. Esto se debe a la diferencia en las constantes de tiempo. Aún así, se puede conseguir una mayor rapidez en el seguimiento de la referencia aumentando Kc, aunque con la siguiente limitación:  Puesto que el ruido de salida se realimenta (la salida se compara con la referencia y posteriormente se introduce en el controlador), un valor alto de Kc significará un aumento de ruido significativo en la entrada, situación que se quiere evitar o reducir todo lo posible para que no se produzcan continuamente cambios bruscos en la válvula que manipula el caudal F200. La interacción con el lazo X2’-XF100 también infuye en la rapidez de la respuesta. Dicho problema de interacción no se puede solucionar completamente con las técnicas de control monovariable. Finalmente, se ajusta la sintonización del controlador controlador manualmente a los siguientes parámetros: Este modelo es un buen ejemplo también para observar cómo un aumento desmesurado de la ganancia del controlador afectaría negativamente a la entrada (pese a la mejora de la salida si la válvula física real 𝒌𝒄= −15 %𝑃2/%𝐹200 𝒕𝒊=𝟒𝟎 𝑚𝑖𝑛 𝑺.𝑶.~𝟎% 𝒕𝒆𝟓%= 180 𝑚𝑖𝑛 Figura 5.34 Respuesta de P2% y F200% frente a cambios en la referencia con PI final 52 pudiese soportar tales variaciones en la entrada), y cómo mejora la salida frente grandes escalones en la referencia con la ayuda del sistema “anti-windup” incorporado [7]. La Figura 5.35 muestra para Kc=-100 dicho efecto en la entrada y la salida para una referencia del 30% con respecto al punto nominal, así como una comparación con el efecto antiwindup desactivado: Figura 5.35 Respuesta de P2% y F200% con kc= -100, con y sin antiwindup Es preferible sin duda un valor de ganancia menor en este caso para que no haya tanto ruido en la entrada, pese a la pérdida de rapidez en el seguimiento de la referencia. Hay que tener en cuenta que P2 es la variable menos importante de las 3 variables controladas, y no se requiere, por tanto, una excesiva precisión en dicho seguimiento. A continuación se muestra la respuesta de L2’ frente a una cadena de perturbaciones, manteniendo la referencia en el punto nominal. En concreto, en el minuto 20 se introduce un escalón de amplitud +20% en T200. Figura 5.36 Control realimentado P2-F200 frente a perturbaciones 53 El rechazo a la perturbación T200 es poco eficaz. Aunque en régimen permanente el sistema consigue rechazar al completo la perturbación, la sobreoscilación es considerable y se mantiene durante largo tiempo. Más adelante se mejorará la respuesta de P2 frente a T200 con otras técnicas de control. 5.4.4 Control realimentado emparejamiento 𝐋𝟐−𝐅𝟐′ GP%(s)=(L2%(s) F2′%(s))=−0.125 s (5–26) Aplicando las relaciones descritas para un sistema integrador de la Tabla 5.1, se calculan los parámetros de los controladores PI (𝜏𝑐=1): PI{𝐤𝐜=1 −0.125·21=−16 𝐭𝐢=2·1=2 (5–27) Como en el caso anterior, al no haber retardo no es necesario incluir acción derivativa para el correcto funcionamiento del control. Este caso, puesto que es un sistema integrador, en realidad implementando un controlador que solo tuviese acción proporcional se conseguría un error en régimen permanente nulo frente a cambios en la referencia: Gbc=Gc·GP 1+Gc·GP=Kc·KP s 1+KcKP s=KcKP s s+KcKP s=1 1 KcKPs+1→r.p.,s=0→Gbc=1 (5–28) Sin embargo, el interés real de control en esta variable no es el seguimiento de una referencia, sino el rechazo de perturbaciones. Puesto que el evaporador presentará problemas graves de funcionamiento si se alcanzan los límites del separador, se mantendrá como referencia fija el nivel del evaporador en L2=1 metro. En cuanto a las perturbaciones, sí es necesario en este caso un término integrador para que en régimen permanente el error provocado por una perturbación sea igual a 0, y se implementa por tanto un PI con los parámetros obtenidos en la ecuación (5–27). Si hay cambios en una perturbación (F1 en este caso, comportamiento de tipo integrador) y no en la referencia, la función de transferencia en bucle cerrado queda como sigue: Gbc=Gd 1+Gc·GP=Gd 1+Tis+1 TisKc·KP 𝑠=Gd·Tis2 Tis2+(Tis+1)kck→r.p.,s=0→Gbc=0 (5–29) Lo cual indica que frente a un cambio en la perturbación F1, con un controlador PI la respuesta L2 tiene un error en régimen permanente=0. En la Figura 5.37 se representa la salida L2 y la variable manipulada F2’ frente a escalones del +10% y -10% en la perturbación F1, manteniendo una referencia L2 fija en el punto nominal de funcionamiento, para el sistema realimentado con el controlador PI y con el controlador P. 54 Figura 5.37 Respuesta de L2% y F2% frente a la perturbación F1 con controlador PI y P Queda clara por tanto la elección a favor del PI. La mejor estrategia para aumentar la rapidez de reacción ante la perturbación F1 es aumentar ligeramente la ganancia del controlador, teniendo en cuenta el hecho de que, como en el caso anterior, el ruido en la entrada se puede amplificar indeseablemente. Teniendo esto en cuenta, se declaran los siguientes parámetros, obteniendo la respuesta de la Figura 5.38 para el mismo escalón en la perturbación del caso anterior. 𝒌𝒄= −25 %𝐿2/%𝐹2 𝒕𝒊=𝟐 𝑚𝑖𝑛 Figura 5.38 Respuesta de L2% y F2% frente a la perturbación F1 con controlador PI final 55 5.5 Control avanzado para rechazo de perturbaciones: acción anticipativa y control en cascada El control descrito hasta ahora ofrece buenas soluciones frente a cambios en la referencia, así como error nulo en régimen permanente para perturbaciones. Sin embargo, el efecto sobreoscilatorio que pueden provocar dichas perturbaciones puede ser significativo, y es de interés reducirlo lo máximo posible. Para ello, se va a utilizar en este proyecto técnicas de control anticipativo y control en cascada [7]. Control Anticipativo (feed-forward) Los sistemas de control por realimentación requieren que exista error en la variable de proceso a controlar para ejercer la acción correctora. Esta forma de actuar implica un cierto retraso en la acción de control y, como consecuencia, una corrección no del todo eficiente frente a perturbaciones externas. La idea básica del control anticipativo es medir las perturbaciones y actuar sobre el proceso inmediatamente que se produzcan, sin tener que esperar a que afecten a la variable que se está controlando. Para ello se ha de disponer de un modelo de comportamiento del proceso frente a las perturbaciones (Gd(s), calculado anteriormente). La Figura 5.39 muestra el diagrama de bloques de la acción anticipativa aplicada a un sistema realimentado. Figura 5.39 Control por acción anticipativa o feedforward Se definen los siguientes elementos nuevos:  GFT (s): Función de transferencia del sensor-transmisor que mide la perturbación.  GF(s): Función de transferencia del control anticipativo. Idealmente, para que la variación en la variable de control sea nula frente a una perturbación, se debe cumplir: GD(s)+GV(s)·Gp(s)·GF(s)·GFT(s)=0 (5–30) Por lo tanto, la función de transferencia ideal del control anticipativo es: GF(s)=− GD(s) GV(s)Gp(s)GFT(s) (5–31) Nótese que en la práctica dicho control no perfecto debido a la imposibilidad de medición de todas las perturbaciones, la inexactitud de la función de transferencia del proceso, errores en las medidas, etc. 56 Control en Cascada En muchos procesos es muy sencillo detectar una perturbación antes de que tenga un efecto apreciable sobre la variable de proceso a controlar. Esta detección, que se lleva a cabo mediante la medida de alguna variable interna, permite actuar rápida e intensamente sobre el proceso, evitando desviaciones importantes en la variable a controlar. Esta idea es particularmente útil cuando existen perturbaciones que afectan a la propia variable de proceso que se manipula para controlarlo (perturbaciones a la entrada). [7] Para estos casos, se utiliza el esquema de control en cascada, representado en la Figura 5.40. El control en cascada introduce un control por realimentación secundario (esclavo), que tiene como objetivo corregir la acción del elemento final de control (una válvula en este caso). Este control se produce anterioridad al lazo primario (maestro), frente a una perturbación que afecta directamente a dicho elemento final de control. Figura 5.40 Control en cascada Donde GCM(s), GCE(s) y GTE(s) representan respectivamente las funciones de transferencia del controlador maestro, el controlador esclavo y el sensor transmisor usado para medir la entrada u(s). Todas las pruebas de simulación desde este punto en adelante se realizan sobre el sistema SISO completo con los 3 bucles de realimentación (ver modelo completo en el Anexo B). De esta manera, esta última parte de diseño de control contempla seguimiento de referencias y rechazo de perturbaciones para las 3 variables de control con sus lazos de realimentación cerrados. 5.5.1 Acción anticipativa para 𝐅𝟏 en el emparejamiento X2’-XF100 Atendiendo a la Figura 5.39, el control anticipativo ofrece la siguiente solución de control: GF(s)=− GD(s) GV(s)Gp(s)GFT(s)=− −11.67 4.165s+1e−0.3s ( %C kg/min) 171 2.475s+1 e−0.93s(%C 1)·100−0 15−7 ( % kg/min)=2.311s+0.934 712.2s+171 1 % (5–32) Nota: el resultado de la división contiene un retardo positivo, eliminado por ser irrealizable. Un retardo positivo implicaría una acción de control que actuase antes de que se produjese la propia peturbación, lo cual no es posible. los límites de las perturbaciones adoptados se muestran en la Tabla 6.1. 57 Teniendo en cuenta que el controlador ofrece un valor de m(s) porcentual, se añade un último término para adimensionalizar correctamente: GF%(s)=XF100 F1(1 %)·XF100 XF100(% 1)=2.311s+0.934 712.2s+171 ·100−0 1−0 =231.1s+93.4 712.2s+171% % (5–33) Al introducir la función GF en el modelo para X2’-xF100 en Simulink, se obtienen las siguientes respuestas para escalones de ± 10% con respecto al punto nominal en la perturbación F1 en los minutos 10 y 30, respectivamente: Figura 5.41 respuesta X2’% y xF100% frente a escalón en F1, con y sin control anticipativo Los resultados con el sistema controlado añadiendo la acción anticipativa de F1 son significativamente más satisfactorios. Como se puede observar, la sobreoscilación provocada en la variable de control por la variación en F1 es mucho menor en el sistema con acción anticipativa (valor de pico 64% frente a 74%). Esto se debe al control feed-forward, que modifica la salida que proporciona el controlador nada más se detecta la perturbación (ver minuto 10 y 30). En el modelo sin acción anticipativa, la salida del controlador no se ve afectada hasta que el control por realimentación evalúa la salida con la referencia, lo cual es un proceso mucho más lento. El control no es perfecto porque variar F1 implica variar X2, que implica variar XF100 y, por tanto, afecta también a L2 (sistema multivariable). 5.5.2 Control prealimentado para 𝐓𝟐𝟎𝟎 en el emparejamiento 𝐏𝟐−𝐅𝟐𝟎𝟎 La perturbación que afecta más claramente a la presión del sistema es la temperatura del agua de refrigeración a la entrada del condensador. Dicha temperatura influye directamente en la cantidad de vapor que se puede condensar y, por tanto, en la presión que existirá en el separador según la cantidad de vapor existente. El controlador anticipativo resulta como sigue: GF(s)=− GD(s) GV(s)Gp(s)GFT(s)=− 1.543 43.78s+1 (KPa C o) −0.0792 43.82s+1(KPa kg/min)· 100−0 37.5−12.5 (%C o)=16.91s+0.3858 3.469s+0.0792 kg/min % (5–34) Teniendo en cuenta que el controlador ofrece un valor de m(s) porcentual, se añade un último término para adimensionalizar correctamente: 64 5.6 Control SISO completo con todos los controladores El diagrama de bloques utilizado en este apartado se muestra en el anexo B (diagrama control SISO completo). Incluye todas las técnicas de control desarrolladas en este capítulo, actuando en conjunto. Están incluidos, por tanto:  3 lazos de realimentación para los emparejamientos XF100 - X2’, F200 - P2 y F2’ - L2.  control feed-forward para F1 - X2’, T200 - P2 y F1 – L2’.  Control en cascada para P0 – X2’. 5.6.1 Control SISO completo frente a perturbaciones A continuación se compara, por un lado, la respuesta de las variables de salida y variables manipulables del sistema SISO completo con todas las estrategias de control previamente desarrolladas implementadas, y por otro, la respuesta de dichas variables sin los controles feedforward y sin el control en cascada, frente a perturbaciones. La tabla siguiente muestra la secuencia de variación de las variables de perturbación: Tiempo(min) 20 40 60 80 100 120 140 160 180 200 P0(bar) +10 -5 F1(kg/min) -15 +5 X1(%) +10 T1(ºC) +20 +10 T200(ºC) +15 +5 Tabla 5.2 Cadena de cambios en perturbaciones La Figura 5.49 muestra la respuesta de todas las variables manipulables y de control del sistema, frente a los cambios mostrados en la Tabla 5.2. Resultados  Clara mejoría en las variables de control en cuanto a sobreoscilaciones. El control anticipativo y en cascada permiten actuar rápidamente frente a las perturbaciones y evitar picos en el transitorio. Este efecto es especialmente notable en la respuesta de X2’ y L2’; los picos del transitorio que se producen en el sistema sin el conrol FF y cascada se reducen hasta en un 92%.  En cuanto a las entradas, frente a las perturbaciones se producen cambios más rápidos en orden a rechazarlas en el sistema con control FF y cascada. Esto se debe a la reacción inmediata de la variable manipulable en cuanto se percibe la perturbación. Esta reacción produce un error más pequeño en la salida, pero se aumenta el esfuerzo de control, y por lo tanto el esfuerzo que tiene que hacer la válvula.  Los controladores diseñados se basan en funciones de transferencia lineales y se implementan en el sistema no lineal. Esta implementación crea lógicamente que la respuesta de las variables de control no sea exactamente la deseada, pues cuanto más alejadas del punto de funcionamiento nominal se encuentren los valores de las perturbaciones, peor será el control al dejar de ser representativas las funciones de transferencia del sistema no lineal. En concreto, cuando se superan en este proyecto valores del ± 20 % con respecto al punto de funcionamiento de las perturbaciones, el control empieza a ser cada vez menos efectivo y las variables de control tienen un mayor tiempo de establecimiento en 65 régimen permanente y unos picos de transitorio más acusados.  El diseño de cada controlador se realiza individualmente para cada pareja variable de control – variable manipulable, manteniendo el resto dev ariables manipulables y perturbaciones en su punto de funcionamiento. Gracias a la DRGA, los emparejamientos elegidos son los más óptimos, pero no dejan de ser un estudio SISO con limitaciones en las interacciones entre variables (aspecto mejorable con un sistema MIMO). Estas interacciones se pueden ver en la Figura 5.49. En concreto, en el minuto 100 de la variable de control X2’%, se observa como existen 2 picos de transitorio al introducirse una variación en F1, en vez de establecerse directamente en la referencia tras el primer pico. Esto se debe a la interacción con otras variables. Figura 5.49 respuesta salidas y entradas del sistema SISO completo frente a cadena de perturbaciones sin control anticipativo-cascada y con él 5.6.2 Control SISO completo frente a seguimiento de referencias A continuación se compara, por un lado, la respuesta de las variables de salida y variables manipulables del sistema SISO completo con todas las estrategias de control previamente desarrolladas implementadas, y por otro, la respuesta de dichas variables sin los controles feedforward y sin el control en cascada, frente a seguimiento de referencias. La tabla siguiente muestra la secuencia de variación de las referencias: Tiempo(min) 20 40 60 80 100 120 140 X2’ref (%) +10 -10 P2 ref (%) -20 0 L2 ref (%) -15 0 Tabla 5.3 Cadena de cambios en referencias 66 Respuesta de todas las variables de entrada y de control del sistema: Figura 5.50 respuesta salidas y entradas del sistema SISO completo frente a cadena de referencias sin control anticipativo-cascada y con él  La respuesta en este caso de las entradas y las salidas es prácticamente idéntica para ambos modelos (sin y con control avanzado integrado). De hecho, el único motivo por el cual se observa una ligera diferencia en X2’% es que el controlador PID diseñado para el modelo en cascada no es exactamente el mismo que el diseñado sin cascada. Esta similitud casi exacta radica en el hecho de que el diseño de controladores anticipativos y en cascada se basan en el rechazo de perturbaciones, pero no modifican la respuesta frente a cambios en las referencias.  Cabe destacar la lenta velocidad de cambio de P2%: la entrada se satura entre 0-400 kg/min por lo que una vez alcanzado dichos límites, como ocurre en este caso, la pendiente de la presión de salida es limitada en unos valores máximo y mínimo. La saturación es una no linealidad muy importante, pues en el diseño de controladores lineales no se tiene en cuenta, por ejemplo, que las válvulas estén limitadas a los valores 0%-100% (completamente cerrada – completamente abierta). Se pueden dar situaciones en el control lineal en las que el controlador ofrezca un valor de salida mayor o menor a esos límites, lo cual es imposible. Para solventar este problema en el sistema no lineal se añaden bloques de saturación, que crean una simulación más realista y limitan la velocidad de respuesta de las variables de control. Se pueden observar también como F2’ y XF100 alcanzan límites de saturación.  Son notables otras no linealidades en la respuesta. Bajar la presión del sistema es más complicado que subirla. Así, cuando se pide una referencia del -20% con respecto al punto de funcionamiento de P2, la 67 entrada F200 se satura en 400 kg/min y aun así la pendiente de bajada de P2 es poco pronunciada. Sin embargo, al volver a requerir un incremento de la referencia, la entrada F200 se satura en 0 kg/min y la presión del evaporador alcanza el valor de referencia en mucho menos tiempo.  Se comprueba la idoneidad de los emparajemanientos variable de control – variable manipulable adoptados tras el análisis de la DRGA. Sin embargo, se pueden apreciar efectos típicos del sistema multivariable. L2 presenta picos en los minuto 20 y 80 debido a los cambios de referencia en X2’ y por tanto en XF100. P2 también se ve afectada en el minuto 80 con un pequeño pico debido al cambio en XF100. Por su parte, X2’ también presenta picos en el transitorio debido a los cambios de F200 en el minuto 100 y de F2’ en el minuto 120. 5.7 Comparación con otros autores La mayoría de autores ([2],[3],[11]) que abordan el modelo del evaporador de tubos largos verticales con recirculación forzada toman como referencia el modelo dinámico de Newell and Lee [6]. Puesto que aquí se modifica este modelo, los resultados que se obtienen son diferentes a otros documentos, aunque sí se pueden distinguir las similitudes y diferencias siguientes:  En cuanto a las entradas y perturbaciones, el modelo de Newell and Lee considera que F3 es una perturbación (ver Capitulo 2), y que P100 es una entrada manipulable del sistema. Puesto que en realidad dichas variables no pueden modificarse directamente, en este proyecto F3 no es una perturbación.  Los emparejamientos escogidos para el control SISO sí coinciden con los escogidos por otros autores [2],[3],[11]. Debido a que en otras publicaciones una variable manipulable es P100, los valores de la matriz RGA son diferentes.  La solución de control realimentado que ofrecen consiste también en integrar controladores PI en los 3 bucles de realimentación. La ganancia de estos controladores es por lo general más alta, teniendo en cuenta que es calculado por el método de Ziegler and Nichols y que no se tiene en cuenta el ruido en las señales.  Con respecto a control prelimentado, Newell and Lee elige también rechazar F1 con control feedforward sobre el bucle de realimentación utilizado para controlar X2. Este proyecto mejora los siguientes aspectos con respecto a publicaciones anteriores:  El modelo dinámico ofrece resultados más acertados y coherentes con respecto al comportamiento real de las variables de control. En la Figura 5.51 [6] se muestran las respuestas de las variables de control X2 y P2 frente a escalones en P100 y F200, respectivamente, según el modelo de Newell and Lee. En dicho modelo, X2 evoluciona siguiendo un comportamiento subamortiguado, y en la respuesta de P2 existen cambios de pendiente abruptos, sin cambios en la variable manipulable. Estos comportamientos son singulares e impropios de un evaporador.  Los controladores se han diseñado por teorías de control modernas (SIMC-Skogestad).  Se desarrolla un modelo completo de simulación en Simulink y se incluyen sensores y ruido blanco en la implementación. Los resultados obtenidos como consecuencia de estas variaciones son mejores y más realistas con respecto a otros estudios previos sobre la materia en cuestión. 68 Figura 5.51 Respuestas X2 frente a +10 KPa en P100 y P2 frente a +10 kg/min en F200. Modelo Newell and Lee 69 6 INTERFAZ USUARIO EN TIEMPO REAL Partiendo como base del modelo de Simulink presentado en el Capítulo 4, se ha desarrollado una GUI (Graphical User Interface) que permite al usuario modificar los principales parámetros del sistema de control y observar los resultados en tiempo real sin necesidad de modificar nada directamente sobre el modelo de Simulink. Para ejecutar la GUI y hacerla funcionar el usuario debe seguir los siguientes pasos: 1. Abrir la aplicación MATLAB® e introducir “Simulink” en la línea de comandos. 2. Abrir el archivo “CONTROLSISO_DEFINITIVO_REALTIME.slx” y pulsar sobre el “bloque Scope” en rojo nombrado como “SALIDAS %” y el bloque azul “ENTRADAS%”. Esto permitirá observar la respuesta de las salidas y las entradas en la gráficas emergentes. 3. Compilar y ejecutar (botón “RUN”) de los archivos “parámetros.m” y “miGUI.m”. Al final de este se desplegará la GUI, que se muestra en la Figura siguiente: Figura 6.1 Pantalla GUI En esta pantalla se puede interactuar con diversos elementos. Lo primero que se debe hacer para ello es pulsar “Buscar” y elegir el archivo “CONTROLSISO_DEFINITIVO_REALTIME.slx”. Posteriormente, al pulsar sobre “Cargar .slx” la GUI quedará enlazada con el sistema de Simulink y se podrá actuar sobre él. En la pantalla se pueden observar 3 bloques en rojo y 5 bloques en amarillo. Corresponden a las referencias de las salidas y a las perturbaciones, respectivamente, en sus valores absolutos. Además, se puede observar una imagen que muestra el esquema del sistema evaporador. Cada bloque tiene una barra “slider”, en la que se puede deslizar el bloque interior de un lado a otro para cambiar el valor de la ganancia con respecto al punto nominal de funcionamiento para modificar el valor de la referencia o perturbación. El valor de dicha ganancia se muestra en el bloque rectangular blanco inferior, donde también se puede introducir un valor numérico concreto. Además, en cada extremo de la barra hay 2 bloques con 2 números, que indican el valor mínimo y máximo que puede tomar la ganancia y valor absoluto de la variable. A modo de ejemplo, el primer bloque rojo, (“Referencia X2’ ”) está limitado entre 0% y 40%, por lo que la 70 ganancia que podemos introducir está comprendida entre 0 y 1.6 (25*0=0, 25*1.6=40). El funcionamiento de la GUI consiste en trasladar los valores que se introducen en ella a las ganancias conectadas a los bloques “step” de las referencias y perturbaciones de Simulink (ver Anexo B), modificándolas a la par que lo hace el usuario. Por otro lado, aunque los límites de las variables de salida se conocen (ver Tabla 1.1), los límites de las perturbaciones son un dato desconocido puesto que éstas son externas y no dependen del sistema de estudio. En este proyecto se estiman los siguientes límites: Punto nominal Límite Inferior (Ganancia) Límite Inferior (Valor Absoluto) Límite Superior (Ganancia) Límite Superior (Valor Absoluto) 𝐏𝟎𝟎 = 5 bar 0.5 2.5 bar 1.5 7.5 bar 𝐅𝟏𝟎= 10 kg/min 0.7 7 kg/min 1.5 15 kg/min 𝐗𝟏𝟎 = 5 % 0 0 % 2 10 % 𝐓𝟏𝟎 = 40 OC 0.5 20 OC 1.5 60 OC 𝐓𝟐𝟎𝟎 𝟎 = 25 OC 0.5 12.5 OC 1.5 37.5 OC Tabla 6.1 Límites perturbaciones 4. Pulsando “Start Simulation” se inicia la simulación, en la que se pueden ir variando las ganancias y ver sus consecuencias directas. Un ejemplo de ello se muestra en la Figura 6.2. Pulsando en “Stop Simulation” se puede parar la simulación y reiniciarla de nuevo con “Start Simulation”. Figura 6.2 Simulación GUI El código utilizado para la creación de la GUI se recoge en el Anexo A. Por otro lado, para conseguir que la simulación funcione ralentizada simulando el efecto del tiempo real, se utiliza el bloque de Simulink “RealTime Pacer” [13]. Se ha creado otra GUI, denominada: “miGUI_nocontrol.m”, que se puede utilizar de la misma manera para la simulación del modelo no lineal no controlado: “modelo_nolineal_final_REALTIME.slx”. 71 7 CONCLUSIONES Y LÍNEAS FUTURAS 7.1 Conclusiones En este proyecto se ha trabajado en el desarrollo del modelo dinámico de un evaporador de recirculación forzada y de su sistema e control, así como en su implementación en Simulink para poder analizar los resultados expresados gráficamente. Tomando como base los resultados obtenidos a lo largo de los capítulos anteriores, se destacan las siguientes conclusiones:  El modelo dinámico ofrece un comportamiento admisible para la dinámica de las variables de salida. Sin embargo, el modelo sigue siendo simplificado y tiene limitaciones, aunque los resultados son satisfactorios bajo las consideraciones de funcionamiento tomadas.  Se han mejorado modelos anteriores y se han incluido cambios importantes. En concreto, en cuanto al modelo dinámico, se ha incluido el modelo de tanques homogéneos en serie para el evaporador, se ha dotado de dinámica a determinadas variables que otros autores no consideran, se ha incluido un purgador de condensado, se ha considerado la presión de vapor saliente de una caldera como perturbación, se ha añadido el modelo de la válvula del caudal de vapor de entrada y se ha modificado la localización de la bomba de impulsión para aproximar el modelo a un evaporador real. Por otro lado, en cuanto al control, se han añadido no linealidades en la simulación, ruido blanco, y se han adimensionalizado las variables de control y variables manipulables para una correcta implementación en controladores P, PI, PID. Además, se utilizan técnicas de diseño de parámetros de controladores más modernas (SIMC) y se hace un estudio de acciones anticipativas y control en cascada. Se han añadido sensores en el modelo de control.  La utilización de Simulink es de gran ayuda en la visualización de los resultados. Con el modelo que se ha creado en este proyecto se pueden modificar fácilmente las variables del sistema y observar las diferentes respuestas en cuestión de segundos, todo ello en una interfaz gráfica intuitiva. La división de diagramas en diferentes niveles permite observar con claridad la implementación de las diferentes ecuaciones. Así, el primer nivel muestra la configuración del evaporador en bloques y las entradas, perturbaciones y salidas principales, lo que convierte al archivo en una herramienta adecuada para la docencia y útil para usuarios que no dispongan de conocimientos avanzados del programa, ya que no se requieren apenas conocimientos del mismo para ejecutarlo y cambiar los principales parámetros. De manera adicional, el diagrama de Simulink sirve como base para futuras ampliaciones del proyecto.  Se ha comprobado la utilidad de la RDGA a la hora de hacer la elección de los mejores emparejamientos SISO. La elección surgida de este análisis ha permitido obtener resultados muy favorables en un sistema que no incluye un control multivariable.  El control realimentado del sistema no lineal con controladores sintonizados por el método SIMC ofrece buenos resultados en el seguimiento de referencias. No se producen errores graves como inestabilidad ni existen errores en régimen permanente ni grandes sobreoscilaciones. A su vez, se comprueba la utilidad de herramientas como “anti-windup” y el filtro derivativo.  La utilización de control prealimentado y control en cascada permite mejorar considerablemente el rechazo frente a perturbaciones. Integrados estos controladores con los de realimentación, se consigue un modelo de control SISO capaz de disminuir el efecto de las perturbaciones y de realizar un seguimiento efectivo de las referencias. 7.2 Líneas futuas 1. Para mejorar la simulación y los resultados de futuros estudios, se puede profundizar más en el modelo dinámico planteado. Hay diversas variables cuyas dinámicas se han despreciado en este 72 proyecto, y que pueden ser objeto de desarrollo, como, por ejemplo, la dinámica del condensador. 2. El siguiente paso lógico con respecto al sistema de control consiste en implementar técnicas de control multivariable para mejorar los problemas que surgen de la interacción entre variables. Para ello, existen diversos campos de de estudio como el regulador multivariable lineal (LMR), así como controles avanzados adaptativos y predictivos. 73 Apéndice A. Códigos Matlab A.1 Parámetros y puntos nominales de funcionamiento Código A.1 parametros.m %parámetros sistema pA= 20; M= 20; C= 4; Cp= 0.07; yr= 38.5; ys= 36.6; UA2= 6.84; pi= 3.1416; %puntos de equilibro (puntos nominales de funcionamiento) P0=5; xF100E=0.5; X1E=5; X2E=25; F1E=10; F2E=2; F3E=50; F4E=8; F5E=8; T1E=40; T2E=84.6; T3E=80.6; L2E=1; P2E=50.5; F100E=9.3; T100E=119.9; P100E=194.7; Q100E=339; F200E=208; T200E=25; T201E=46.1; Q200E=307.9; A.2 Cálculo de puntos nominales de funcionamiento, (evaporador N compartimentos) Código A.2 pnominalesevap5.m function [matsol]= pnominalesevap5 (F1, T1, X1, F100, P2, n) Cp=0.08; ys=36.6; Cpv= 0.025;... UA=100.11; syms T21 X21 F41 F21 TW1; eq1= F1*Cp*T1-F21*Cp*T21-F41*42.9+UA*(TW1-T21)/n; 80 % --- Executes on button press in pushbutton_quitarslx. function pushbutton_quitarslx_Callback(hObject, eventdata, handles) % hObject handle to pushbutton_quitarslx (see GCBO) % eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) % --- Executes on button press in pushbutton_run. function pushbutton_run_Callback(hObject, eventdata, handles) % hObject handle to pushbutton_run (see GCBO) % eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) mystring = get(hObject,'String'); status = get_param(bdroot,'simulationstatus'); if strcmp(mystring,'Start Simulation') % Check the status of the simulation and start it if it's stopped if strcmp(status,'stopped') set_param(bdroot,'simulationcommand','start') end % Update the string on the pushbutton set(handles.pushbutton_run,'String','Stop Simulation') elseif strcmp(mystring,'Stop Simulation') % Check the status of the simulation and stop it if it's running if strcmp(status,'running') set_param(bdroot, 'SimulationCommand', 'Stop') end % Update the string on the pushbutton set(handles.pushbutton_run,'String','Start Simulation') else warning('Unrecognized string for pushbutton_run') %#ok<WNTAG> end % Assign handles and the startstop object to the base workspace assignin('base','miGUI_handles',handles) assignin('base','run_hObject',handles.pushbutton_run) function edit_l2_Callback(hObject, eventdata, handles) % hObject handle to edit_l2 (see GCBO) % eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) % Hints: get(hObject,'String') returns contents of edit_l2 as text % str2double(get(hObject,'String')) returns contents of edit_l2 as a double value = get(hObject,'String'); %valores lÌmites 2-0 if (str2double(value)>2.0) 81 value=num2str(2.0); end if (str2double(value)< 0) value=num2str(0); end % Update the model's gain value set_param([bdroot '/GainC'],'Gain',value) set(handles.edit_l2,'String',value); slider_position =str2double(value); set(handles.slider_l2,'Value',slider_position); % Update simulation if the model is running status = get_param(bdroot,'simulationstatus'); if strcmp(status,'running') set_param(bdroot, 'SimulationCommand', 'Update') end guidata(hObject,handles); % --- Executes during object creation, after setting all properties. function edit_l2_CreateFcn(hObject, eventdata, handles) % hObject handle to edit_l2 (see GCBO) % eventdata reserved - to be defined in a future version of MATLAB % handles empty - handles not created until after all CreateFcns called % Hint: edit controls usually have a white background on Windows. % See ISPC and COMPUTER. if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) set(hObject,'BackgroundColor','white'); end % --- Executes on slider movement. function slider_l2_Callback(hObject, eventdata, handles) % hObject handle to slider_l2 (see GCBO) % eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) % Hints: get(hObject,'Value') returns position of slider % get(hObject,'Min') and get(hObject,'Max') to determine range of slider slider_position = get(hObject,'Value'); value = num2str(slider_position); % Update the model's gain value set_param([bdroot '/GainC'],'Gain',value) % Set the value of the gain edit box set(handles.edit_l2,'String',value); % Update simulation if the model is running status = get_param(bdroot,'simulationstatus'); if strcmp(status,'running') set_param(bdroot, 'SimulationCommand', 'Update') end guidata(hObject,handles); 82 % --- Executes during object creation, after setting all properties. function slider_l2_CreateFcn(hObject, eventdata, handles) % hObject handle to slider_l2 (see GCBO) % eventdata reserved - to be defined in a future version of MATLAB % handles empty - handles not created until after all CreateFcns called % Hint: slider controls usually have a light gray background. if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) set(hObject,'BackgroundColor',[.9 .9 .9]); end function edit_p2_Callback(hObject, eventdata, handles) % hObject handle to edit_p2 (see GCBO) % eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) % Hints: get(hObject,'String') returns contents of edit_p2 as text % str2double(get(hObject,'String')) returns contents of edit_p2 as a double value = get(hObject,'String'); %valores lÌmites 2-0 if (str2double(value)>2.0) value=num2str(2.0); end if (str2double(value)< 0) value=num2str(0); end % Update the model's gain value set_param([bdroot '/GainB'],'Gain',value) set(handles.edit_p2,'String',value); slider_position =str2double(value); set(handles.slider_p2,'Value',slider_position); % Update simulation if the model is running status = get_param(bdroot,'simulationstatus'); if strcmp(status,'running') set_param(bdroot, 'SimulationCommand', 'Update') end guidata(hObject,handles); % --- Executes during object creation, after setting all properties. function edit_p2_CreateFcn(hObject, eventdata, handles) % hObject handle to edit_p2 (see GCBO) % eventdata reserved - to be defined in a future version of MATLAB % handles empty - handles not created until after all CreateFcns called % Hint: edit controls usually have a white background on Windows. % See ISPC and COMPUTER. if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) 83 set(hObject,'BackgroundColor','white'); end % --- Executes on slider movement. function slider_p2_Callback(hObject, eventdata, handles) % hObject handle to slider_p2 (see GCBO) % eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) % Hints: get(hObject,'Value') returns position of slider % get(hObject,'Min') and get(hObject,'Max') to determine range of slider slider_position = get(hObject,'Value'); value = num2str(slider_position); % Update the model's gain value set_param([bdroot '/GainB'],'Gain',value) % Set the value of the gain edit box set(handles.edit_p2,'String',value); % Update simulation if the model is running status = get_param(bdroot,'simulationstatus'); if strcmp(status,'running') set_param(bdroot, 'SimulationCommand', 'Update') end guidata(hObject,handles); % --- Executes during object creation, after setting all properties. function slider_p2_CreateFcn(hObject, eventdata, handles) % hObject handle to slider_p2 (see GCBO) % eventdata reserved - to be defined in a future version of MATLAB % handles empty - handles not created until after all CreateFcns called % Hint: slider controls usually have a light gray background. if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) set(hObject,'BackgroundColor',[.9 .9 .9]); end function edit_t200_Callback(hObject, eventdata, handles) % hObject handle to edit_t200 (see GCBO) % eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) % Hints: get(hObject,'String') returns contents of edit_t200 as text % str2double(get(hObject,'String')) returns contents of edit_t200 as a double value = get(hObject,'String'); %valores lÌmites 1.5-0.5 if (str2double(value)>1.5) value=num2str(1.5); end if (str2double(value)< 0.5) value=num2str(0.5); end % Update the model's gain value 84 set_param([bdroot '/GainH'],'Gain',value) set(handles.edit_t200,'String',value); slider_position =str2double(value); set(handles.slider_t200,'Value',slider_position); % Update simulation if the model is running status = get_param(bdroot,'simulationstatus'); if strcmp(status,'running') set_param(bdroot, 'SimulationCommand', 'Update') end guidata(hObject,handles); % --- Executes during object creation, after setting all properties. function edit_t200_CreateFcn(hObject, eventdata, handles) % hObject handle to edit_t200 (see GCBO) % eventdata reserved - to be defined in a future version of MATLAB % handles empty - handles not created until after all CreateFcns called % Hint: edit controls usually have a white background on Windows. % See ISPC and COMPUTER. if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) set(hObject,'BackgroundColor','white'); end % --- Executes on slider movement. function slider_t200_Callback(hObject, eventdata, handles) % hObject handle to slider_t200 (see GCBO) % eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) % Hints: get(hObject,'Value') returns position of slider % get(hObject,'Min') and get(hObject,'Max') to determine range of slider slider_position = get(hObject,'Value'); value = num2str(slider_position); % Update the model's gain value set_param([bdroot '/GainH'],'Gain',value) % Set the value of the gain edit box set(handles.edit_t200,'String',value); % Update simulation if the model is running status = get_param(bdroot,'simulationstatus'); if strcmp(status,'running') set_param(bdroot, 'SimulationCommand', 'Update') end guidata(hObject,handles); % --- Executes during object creation, after setting all properties. function slider_t200_CreateFcn(hObject, eventdata, handles) % hObject handle to slider_t200 (see GCBO) % eventdata reserved - to be defined in a future version of MATLAB 85 % handles empty - handles not created until after all CreateFcns called % Hint: slider controls usually have a light gray background. if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) set(hObject,'BackgroundColor',[.9 .9 .9]); end function edit_t1_Callback(hObject, eventdata, handles) % hObject handle to edit_t1 (see GCBO) % eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) % Hints: get(hObject,'String') returns contents of edit_t1 as text % str2double(get(hObject,'String')) returns contents of edit_t1 as a double value = get(hObject,'String'); %valores lÌmites 1.5-0.5 if (str2double(value)>1.5) value=num2str(1.5); end if (str2double(value)< 0.5) value=num2str(0.5); end % Update the model's gain value set_param([bdroot '/GainG'],'Gain',value) set(handles.edit_t1,'String',value); slider_position =str2double(value); set(handles.slider_t1,'Value',slider_position); % Update simulation if the model is running status = get_param(bdroot,'simulationstatus'); if strcmp(status,'running') set_param(bdroot, 'SimulationCommand', 'Update') end guidata(hObject,handles); % --- Executes during object creation, after setting all properties. function edit_t1_CreateFcn(hObject, eventdata, handles) % hObject handle to edit_t1 (see GCBO) % eventdata reserved - to be defined in a future version of MATLAB % handles empty - handles not created until after all CreateFcns called % Hint: edit controls usually have a white background on Windows. % See ISPC and COMPUTER. if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) set(hObject,'BackgroundColor','white'); end % --- Executes on slider movement. function slider_t1_Callback(hObject, eventdata, handles) % hObject handle to slider_t1 (see GCBO) % eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) % Hints: get(hObject,'Value') returns position of slider 86 % get(hObject,'Min') and get(hObject,'Max') to determine range of slider slider_position = get(hObject,'Value'); value = num2str(slider_position); % Update the model's gain value set_param([bdroot '/GainG'],'Gain',value) % Set the value of the gain edit box set(handles.edit_t1,'String',value); % Update simulation if the model is running status = get_param(bdroot,'simulationstatus'); if strcmp(status,'running') set_param(bdroot, 'SimulationCommand', 'Update') end guidata(hObject,handles); % --- Executes during object creation, after setting all properties. function slider_t1_CreateFcn(hObject, eventdata, handles) % hObject handle to slider_t1 (see GCBO) % eventdata reserved - to be defined in a future version of MATLAB % handles empty - handles not created until after all CreateFcns called % Hint: slider controls usually have a light gray background. if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) set(hObject,'BackgroundColor',[.9 .9 .9]); end function edit_x1_Callback(hObject, eventdata, handles) % hObject handle to edit_x1 (see GCBO) % eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) % Hints: get(hObject,'String') returns contents of edit_x1 as text % str2double(get(hObject,'String')) returns contents of edit_x1 as a double value = get(hObject,'String'); %valores lÌmites 2-0 if (str2double(value)>2) value=num2str(2); end if (str2double(value)< 0) value=num2str(0); end % Update the model's gain value set_param([bdroot '/GainF'],'Gain',value) set(handles.edit_x1,'String',value); slider_position =str2double(value); set(handles.slider_x1,'Value',slider_position); 87 % Update simulation if the model is running status = get_param(bdroot,'simulationstatus'); if strcmp(status,'running') set_param(bdroot, 'SimulationCommand', 'Update') end guidata(hObject,handles); % --- Executes during object creation, after setting all properties. function edit_x1_CreateFcn(hObject, eventdata, handles) % hObject handle to edit_x1 (see GCBO) % eventdata reserved - to be defined in a future version of MATLAB % handles empty - handles not created until after all CreateFcns called % Hint: edit controls usually have a white background on Windows. % See ISPC and COMPUTER. if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) set(hObject,'BackgroundColor','white'); end % --- Executes on slider movement. function slider_x1_Callback(hObject, eventdata, handles) % hObject handle to slider_x1 (see GCBO) % eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) % Hints: get(hObject,'Value') returns position of slider % get(hObject,'Min') and get(hObject,'Max') to determine range of slider slider_position = get(hObject,'Value'); value = num2str(slider_position); % Update the model's gain value set_param([bdroot '/GainF'],'Gain',value) % Set the value of the gain edit box set(handles.edit_x1,'String',value); % Update simulation if the model is running status = get_param(bdroot,'simulationstatus'); if strcmp(status,'running') set_param(bdroot, 'SimulationCommand', 'Update') end guidata(hObject,handles); % --- Executes during object creation, after setting all properties. function slider_x1_CreateFcn(hObject, eventdata, handles) % hObject handle to slider_x1 (see GCBO) % eventdata reserved - to be defined in a future version of MATLAB % handles empty - handles not created until after all CreateFcns called % Hint: slider controls usually have a light gray background. if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) set(hObject,'BackgroundColor',[.9 .9 .9]); end function edit_f1_Callback(hObject, eventdata, handles) % hObject handle to edit_f1 (see GCBO) 88 % eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) % Hints: get(hObject,'String') returns contents of edit_f1 as text % str2double(get(hObject,'String')) returns contents of edit_f1 as a double value = get(hObject,'String'); %valores lÌmites 1.5-0.7 if (str2double(value)>1.5) value=num2str(1.5); end if (str2double(value)< 0.7) value=num2str(0.7); end % Update the model's gain value set_param([bdroot '/GainE'],'Gain',value) set(handles.edit_f1,'String',value); slider_position =str2double(value); set(handles.slider_f1,'Value',slider_position); % Update simulation if the model is running status = get_param(bdroot,'simulationstatus'); if strcmp(status,'running') set_param(bdroot, 'SimulationCommand', 'Update') end guidata(hObject,handles); % --- Executes during object creation, after setting all properties. function edit_f1_CreateFcn(hObject, eventdata, handles) % hObject handle to edit_f1 (see GCBO) % eventdata reserved - to be defined in a future version of MATLAB % handles empty - handles not created until after all CreateFcns called % Hint: edit controls usually have a white background on Windows. % See ISPC and COMPUTER. if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) set(hObject,'BackgroundColor','white'); end % --- Executes on slider movement. function slider_f1_Callback(hObject, eventdata, handles) % hObject handle to slider_f1 (see GCBO) % eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) % Hints: get(hObject,'Value') returns position of slider % get(hObject,'Min') and get(hObject,'Max') to determine range of slider slider_position = get(hObject,'Value'); value = num2str(slider_position); 89 % Update the model's gain value set_param([bdroot '/GainE'],'Gain',value) % Set the value of the gain edit box set(handles.edit_f1,'String',value); % Update simulation if the model is running status = get_param(bdroot,'simulationstatus'); if strcmp(status,'running') set_param(bdroot, 'SimulationCommand', 'Update') end guidata(hObject,handles); % --- Executes during object creation, after setting all properties. function slider_f1_CreateFcn(hObject, eventdata, handles) % hObject handle to slider_f1 (see GCBO) % eventdata reserved - to be defined in a future version of MATLAB % handles empty - handles not created until after all CreateFcns called % Hint: slider controls usually have a light gray background. if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) set(hObject,'BackgroundColor',[.9 .9 .9]); end function edit_p0_Callback(hObject, eventdata, handles) % hObject handle to edit_p0 (see GCBO) % eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) % Hints: get(hObject,'String') returns contents of edit_p0 as text % str2double(get(hObject,'String')) returns contents of edit_p0 as a double value = get(hObject,'String'); %valores lÌmites 1.5-0.5 if (str2double(value)>1.5) value=num2str(1.5); end if (str2double(value)< 0.5) value=num2str(0.5); end % Update the model's gain value set_param([bdroot '/GainD'],'Gain',value) set(handles.edit_p0,'String',value); slider_position =str2double(value); set(handles.slider_p0,'Value',slider_position); % Update simulation if the model is running status = get_param(bdroot,'simulationstatus'); if strcmp(status,'running') set_param(bdroot, 'SimulationCommand', 'Update') end guidata(hObject,handles); % --- Executes during object creation, after setting all properties. 96 Evaporador. Nivel 2 (Concentraciones y caudales líquidos) Intercambiador de Calor. Nivel 2 (Cámara y válvulas) Evaporador. Nivel 2 (Balance de masa, vapor de proceso) 97 Evaporador. Nivel 2 (Líquido de proceso-temperatura T2) Evaporador. Nivel 2 (Q100 -> Q) 98 Evaporador. Nivel 2 (Balance de energía, líquido de proceso F4) B.2 MODELO NO LINEAL. CONTROL SISO COMPLETO La siguiente imagen muestra el control completo del sistema. Para ello, se crea un subsistema de todo el “Nivel 0” visto al inicio de este apéndice, y se añaden los controles prealimentados, realimentados y en cascada. 99 Diagrama Control SISO Completo 100 101 Índice de figuras Figura 1.1 Esquema de funcionamiento Evaporado 2 Figura 2.1 Esquema subdivisión el evaporador 6 Figura 2.2 Sección longitudinal interior de una de las mitades simétricas del intercambiador de calor 7 Figura 2.3 Variables separador, bomba y producto 9 Figura 4.1 Pantalla principal modelo no lineal Simulink 19 Figura 4.2 Ejemplo de respuesta de simulación 20 Figura 4.3 Respuestas variables de control frente a variación en XF100 21 Figura 4.4 Respuestas variables de control frente a variación en F200 21 Figura 4.5 Respuestas variables de control frente a variación en F2 22 Figura 4.6 Respuestas variables de control frente a variación en P0 22 Figura 4.7 Respuestas variables de control frente a variación en F1 23 Figura 4.8 Respuestas variables de control frente a variación en X1 23 Figura 4.9 Respuestas variables de control frente a variación en T1 24 Figura 4.10 Respuestas variables de control frente a variación en T200 25 Figura 5.1 respuesta en scope-Simulink 28 Figura 5.2 Pantalla de estimación de modelo Ident 29 Figura 5.3 test ajuste respuesta 29 Figura 5.4 Respuesta sistema nolineal-lineal. X2’ frente a escalón amplitud ± 10% en XF100, y ajuste 31 Figura 5.5 Respuesta sistema nolineal-lineal. P2 frente a escalón amplitud ± 10% en XF100, y ajuste 31 Figura 5.6 Respuesta sistema nolineal-lineal. L2 frente a escalón amplitud ± 10% en XF100, y ajuste 31 Figura 5.7 Respuesta sistema nolineal-lineal. X2’ frente a escalón amplitud ± 10% en F200, y ajuste 32 Figura 5.8 Respuesta sistema nolineal-lineal. P2 frente a escalón amplitud ± 10% en F200, y ajuste 32 Figura 5.9 Respuesta sistema nolineal-lineal. L2 frente a escalón amplitud ± 10% en F200, y ajuste 32 Figura 5.10 Respuesta sistema nolineal-lineal. L2 frente a escalón amplitud ± 10% en F2, y ajuste 33 Figura 5.11 Respuesta sistema nolineal-lineal. X2’ frente a escalón amplitud ± 10% en P0, y ajuste 34 Figura 5.12 Respuesta sistema nolineal-lineal. P2 frente a escalón amplitud ± 10% en P0, y ajuste 35 Figura 5.13 Respuesta sistema nolineal-lineal. L2 frente a escalón amplitud ± 10% en P0, y ajuste 35 Figura 5.14 Respuesta sistema nolineal-lineal. X2’ frente a escalón amplitud ± 10% en F1, y ajuste 35 Figura 5.15 Respuesta sistema nolineal-lineal. P2 frente a escalón amplitud ± 10% en F1, y ajuste 36 Figura 5.16 Respuesta sistema nolineal-lineal. L2 frente a escalón amplitud ± 10% en F1, y ajuste 36 Figura 5.17 Respuesta sistema nolineal-lineal. X2’ frente a escalón amplitud ± 10% en X1, y ajuste 36 102 Figura 5.18 Respuesta sistema nolineal-lineal. P2 frente a escalón amplitud ± 10% en X1, y ajuste 37 Figura 5.19 Respuesta sistema nolineal-lineal. L2 frente a escalón amplitud ± 10% en X1, y ajuste 37 Figura 5.20 Respuesta sistema nolineal-lineal. X2’ frente a escalón amplitud ± 10% en T1, y ajuste 37 Figura 5.21 Respuesta sistema nolineal-lineal. P2 frente a escalón amplitud ± 10% en T1, y ajuste 38 Figura 5.22 Respuesta sistema nolineal-lineal. L2 frente a escalón amplitud ± 10% en T1, y ajuste 38 Figura 5.23 Respuesta sistema nolineal-lineal. X2’ frente a escalón amplitud ± 10% en T200, y ajuste 38 Figura 5.24 Respuesta sistema nolineal-lineal. P2 frente a escalón amplitud ± 10% en T200, y ajuste 39 Figura 5.25 Respuesta sistema nolineal-lineal. L2 frente a escalón amplitud ± 10% en T200, y ajuste 39 Figura 5.26 Modelo de control realimentado 41 Figura 5.27 Diagrama de bloques Simulink X2-xF100 45 Figura 5.28 Esquema Simulink PID 45 Figura 5.29 Respuesta de X2'% y xF100% frente a cambios en la referencia para PI y PID 46 Figura 5.30 Respuesta de X2'% y xF100% frente a cambios en la referencia para PID final 47 Figura 5.31 Control realimentado X2’-XF100 frente a perturbaciones 48 Figura 5.32 Respuesta de X2'% y xF100% frente a cambios en la referencia con diferentes filtros PID 49 Figura 5.33 Respuesta de P2% y F200% frente a cambios en la referencia con PI 51 Figura 5.35 Respuesta de P2% y F200% con kc= -100, con y sin antiwindup 52 Figura 5.36 Control realimentado P2-F200 frente a perturbaciones 52 Figura 5.37 Respuesta de L2% y F2% frente a la perturbación F1 con controlador PI y P 54 Figura 5.39 Control por acción anticipativa o feedforward 55 Figura 5.40 Control en cascada 56 Figura 5.41 respuesta X2’% y xF100% frente a escalón en F1, con y sin control anticipativo 57 Figura 5.42 Salida P2% y entrada F200% frente a escalón en T200, sin y con control anticipativo 58 Figura 5.43 Salida L2% y entrada F2’% frente a escalón en F1, sin y con acción anticipativa 59 Figura 5.44 Salida L2% y entrada F2’% frente a escalón en F1, con FF y resto de bucles abiertos 60 Figura 5.45 Salida L2% y entrada F2’% frente a escalón en F1, sin y con control anticipativo 2 60 Figura 5.46 Esquema de Simulink para el control en cascada del emparejamiento xF100-X2’ 61 Figura 5.47 : salida X2’% y entrada xF100% frente a perturbación P0 sin y con control en cascada 63 Figura 5.48 Control en cascada + acción anticipativa 63 Figura 5.49 respuesta salidas y entradas del sistema SISO completo frente a cadena de perturbaciones sin control anticipativo-cascada y con él 65 Figura 5.50 respuesta salidas y entradas del sistema SISO completo frente a cadena de referencias sin control anticipativo-cascada y con él 66 Figura 5.51 Respuestas X2 frente a +10 KPa en P100 y P2 frente a +10 kg/min en F200. Modelo Newell and Lee 68 Figura 6.1 Pantalla GUI 69 103 Figura 6.2 Simulación GUI 70 104 105 Bibliografía [1] Anibal Alberto Bizama Soto Blog. (2012). Recuperado el 13 de Julio de 2016, de http://anibalbizama.blogspot.com.es/2012/11/9-sistema-de-control-de-procesos.html [2] Cao, Y. (2010). CONSTRAINED SELF-OPTIMIZING CONTROL VIA DIFFERENTIATION. Cranfield, UK. [3] Dittmar, R. (2015). Decentralized SISO Active Disturbance Rejection Control of the Newell-Lee forced circulation evaporator. Heide, Alemania. [4] Fernández Benítez, J., & Corrochano Sánchez, C. (2012). Cuadernos de Transmisión de Calor. Madrid: Sección de Publicaciones ETSII-UPM. [5] imm. (s.f.). Recuperado el 13 de Julio de 2016, de http://www.iim.unsj.edu.ar/control/ [6] Newell, R., & Lee, P. (1989). Applied Process Control: A Case Study. Brisbane, Australia: Prentice Hall. [7] Ollero de Castro, P., & Fernández Camacho, E. (2012). Instrumentación y control de plantas químicas. Sevilla: Sintesis. [8] Samson. (2012). Application Notes for Valves Sizing. Sizing examples. Frankfurt Am main. [9] Sastrón, F. (2013). Apuntes Ingeniería de Control, ETSII-UPM. Madrid. [10] Skogestad, S., & Postlethwaite, I. (2001). Multivariable Feedback Control, Analysis and Design. John Wiley and Sons. [11] Switalski, S. (1999). Control system with the specific performance index for an evaporator. Gliwice, Polonia. [12] Wuhua Hu, Wen-Jian Cai, & Gaoxi Xiao. (2010). Relative Gain Array for MIMO Processes Containing Integrators and/or Differentiators. Singapore. [13] Gautam Vallabha. Real-Time pacer for Simulink. Mathworks. Recuperado el 10 de Julio de 2016, de www.mathworks.com