PONTIFICIA UNIVERSIDAD CATÓLICA DEL PERÚ ESCUELA DE GRADUADOS MAESTRÍA EN INGENIERÍA DE CONTROL Y AUTOMATIZACIÓN “IDENTIFICACIÓN Y CONTROL MULTIVARIABLE DE UNA PLANTA PILOTO DE DESALINIZACIÓN POR OSMOSIS INVERSA” Tesis para optar por el grado de MAGISTER EN INGENIERÍA DE CONTROL Y AUTOMATIZACIÓN Autor: Jimmy Mendoza Montalvo Asesor: PhD Ing. Javier Sotomayor Moriano LIMA – PERÚ 2016 1 DEDICATORIA A mis padres Mario y Clara, y a mi tía Frida, quienes con su infinito amor y paciencia han sabido encaminar mi formación. A mis hermanos, Mario, Yaskha, John y Milton, en quienes confío y deseo todo lo mejor. A Esperanza, mi compañera de vida, quien me alienta siempre para continuar y culminar todos mis proyectos. 2 AGRADECIMIENTOS Un especial agradecimiento a mis profesores por todo el conocimiento brindado para hacer de mí, un excelente profesional. A mis colegas de la maestría quienes siempre me han apoyado en el aprendizaje de esta nueva disciplina de la ingeniería. Finalmente, a Concytec por haber hecho posible realizar mis estudios de posgrado con profunda dedicación. 3 RESUMEN El presente trabajo de tesis tiene por intención establecer los principios de la tecnología de desalinización de agua de mar, la cual constituye una solución prometedora a la escasez de agua en el mundo. Teniendo en cuenta la actual situación de nuestro país, en el que a pesar de ser uno de los primeros países con mayor disponibilidad de recursos hídricos, se hace evidente la necesidad de mejorar los sistemas de distribución de agua y más aún poner en marcha las nuevas tecnologías que se ofrecen para suplir la carencia de este elemento fundamental, implementando además adecuadas técnicas de control que permitan lograr una tasa de producción constante y una pureza aceptable. Una de las tecnologías que se ha impuesto hoy en día, debido a sus bajos costos de producción y mantenimiento respecto de las técnicas de destilación, es la desalinización por ósmosis inversa, que en términos generales es un medio con el que se logra obtener agua pura a partir del agua de mar mediante la aplicación de una presión relativamente baja a través de una membrana semipermeable. En el Perú, esta tecnología a gran escala es poco aprovechada, por lo que es necesario desarrollar una mayor investigación sobre la misma con el objeto de consolidar los conocimientos que nos permitan operar de manera eficiente las plantas de desalinización actuales y futuras. Es así que nace la propuesta de equipar a la universidad con una planta piloto a escala de desalinización por ósmosis inversa, y la cual se expone en la presente tesis desde su diseño hasta la implementación en el Laboratorio de Control Avanzado de la PUCP considerando todos los detalles técnicos y económicos. Esta planta permitirá a los alumnos generar investigación acerca del aumento de la productividad y la calidad del agua, así como la mejora en la eficiencia energética y la seguridad del funcionamiento de este tipo de plantas. Esta tesis presenta además, la metodología a seguir para la correcta identificación del sistema en estudio con el fin de obtener un modelo adecuado que describa el comportamiento dinámico del proceso de desalinización por ósmosis inversa. Por último, se diseña un sistema de control clásico multivariable con la correspondiente propuesta de implementación práctica en la planta piloto de desalinización. 4 AGRADECIMIENTOS El autor de este trabajo de tesis, Jimmy Mendoza Montalvo, agradece el apoyo del Programa Nacional de Innovación para la Competitividad y Productividad (Innóvate Perú) entidad que financió el proyecto 207-FINCyT-IA-2013, en el marco del cual se desarrolló la presente tesis: “Identificación y Control Multivariable de una Planta Piloto de Desalinización por Osmosis Inversa”. 5 ÍNDICE INTRODUCCIÓN ................................................................................................................ 1 CAPÍTULO 1. ESTUDIO DE LAS PLANTAS DESALINIZADORAS DE AGUA DE MAR. .......................................................................................................... 3 1.1. Estado del Arte ........................................................................................................ 3 1.2. Descripción de una Planta Desalinizadora de Agua de Mar por Ósmosis Inversa . 5 CAPÍTULO 2. MODELADO Y CONTROL DE PLANTAS DESALINIZADORAS POR OSMOSIS INVERSA. ................................................................... 10 2.1. Modelado de Plantas Desalinizadoras de Agua de Mar por Ósmosis Inversa ...... 10 2.1.1. Modelado Matemático y Modelado Mediante Identificación de Sistemas ... 10 2.1.2. Comparación de Modelos de Plantas Desalinizadoras de Agua de Mar por Ósmosis Inversa .......................................................................................................... 17 2.2. Control de Plantas Desalinizadoras de Agua de Mar por Osmosis Inversa .......... 20 2.2.1. Revisión de las Estrategias de Control .......................................................... 20 2.2.2. Control Clásico ............................................................................................. 21 2.2.3. Control Avanzado ......................................................................................... 23 2.3. Objetivo de la Tesis .............................................................................................. 29 CAPÍTULO 3. IDENTIFICACIÓN Y DISEÑO DEL SISTEMA DE CONTROL DE UNA PLANTA PILOTO DE DESALINIZACIÓN POR ÓSMOSIS INVERSA. ............................................................................................... 30 3.1. Introducción .......................................................................................................... 30 3.2. Descripción de la Planta Piloto de Desalinización por Ósmosis Inversa del Laboratorio de Control y Automatización de la PUCP.................................................... 33 3.2.1. Implementación de la Planta Piloto .............................................................. 33 3.2.2. Características Generales de la Planta Piloto. ............................................... 36 3.2.3. Características Técnicas de la Planta Piloto. ................................................. 37 3.2.4. Adquisición de Datos de la Planta Piloto. ..................................................... 50 3.3. Identificación de la Planta Piloto. ......................................................................... 51 3.3.1. Generalidades ................................................................................................ 51 3.3.2. Metodología de la Identificación .................................................................. 52 3.3.3. Métodos de Identificación ............................................................................. 54 3.3.4. Convergencia e Identificabilidad .................................................................. 55 3.3.5. Diseño del Experimento ................................................................................ 56 3.3.6. Identificación No Paramétrica....................................................................... 57 3.3.7. Identificación Paramétrica Fuera de Línea ................................................... 61 3.4. Diseño del Sistema de Control de la Planta Piloto ................................................ 77 3.4.1. Modelo Seleccionado para el Diseño del Controlador .................................. 77 3.4.2. Interacción entre Lazos ................................................................................. 78 3.4.3. Emparejamiento de Variables ....................................................................... 80 3.4.4. Desacoplamiento del Sistema ....................................................................... 81 3.4.5. Control Multivariable por Desacoplo utilizando PID ................................... 86 3.4.6. Análisis de Resultados .................................................................................. 88 CAPÍTULO 4. PROPUESTA DE IMPLEMENTACIÓN PRÁCTICA DEL SISTEMA DE CONTROL. .................................................................... 92 4.1. Introducción. ......................................................................................................... 92 4.2. Algoritmo de Control. ........................................................................................... 92 4.3. Hardware del Sistema de Control. ........................................................................ 95 4.4. Software del Sistema de Control ........................................................................... 96 CONCLUSIONES............................................................................................................... 99 RECOMENDACIONES ................................................................................................... 100 BIBLIOGRAFÍA............................................................................................................... 101 ANEXOS 107 LISTA DE FIGURAS Figura 1.1. Esquema del proceso de desalinización de agua de mar por ósmosis inversa de la planta de Llobregat, Barcelona, España. ................................................................................. 7 Figura 1.2. Membrana de ósmosis inversa de tipo enrollamiento en espiral. ......................... 8 Figura 2.1. Esquema de los subsistemas de una planta de RO. ............................................ 13 Figura 2.2. Esquema de simulación de la planta de RO. ....................................................... 13 Figura 2.3. Diagrama de simulación de la planta de desalinización. .................................... 18 Figura 2.4. Respuesta del flujo volumétrico de permeado al escalón de presión en la alimentación para varios modelos. ........................................................................................ 19 Figura 2.5. Respuesta de la conductividad eléctrica al escalón de presión y pH en la alimentación para varios modelos. ........................................................................................ 19 Figura 2.6. Técnicas de control utilizadas en plantas de desalinización. .............................. 20 Figura 2.7. Comparación de la respuesta de flujo. ................................................................ 22 Figura 2.8. Comparación de la respuesta de conductividad. ................................................. 22 Figura 2.9. Control adaptivo directo de plantas no lineales usando redes neuronales. ......... 23 Figura 2.10. Comportamiento del flujo de agua producto. ................................................... 24 Figura 2.11. Comportamiento de la salinidad del agua producto.......................................... 24 Figura 2.12. Respuestas temporales del sistema de control del bastidor de ósmosis inversa con controladores DMC vs PID. ........................................................................................... 26 Figura 2.13. Señales de control del sistema de control del bastidor de ósmosis inversa con controladores DMC vs PID. .................................................................................................. 26 Figura 2.14. Flujo de permeado y señal de control para el controlador PID y varios modelos de la planta. ........................................................................................................................... 27 Figura 2.15. Resultados para un cambio en el set-point del 10% de flujo de permeado utilizando un controlador PI. ................................................................................................. 28 Figura 2.16. Resultados para un cambio en el set-point del 10% de flujo de permeado utilizando un controlador MPC. ............................................................................................ 28 Figura 3.1. Diagrama P&ID de la planta piloto de desalinización. ....................................... 31 Figura 3.2. Planta piloto de desalinización por ósmosis inversa. .......................................... 36 Figura 3.3. Tanques de alimentación (izquierda), tanque de agua producto y agua de rechazo (centro), y tanque de ácido (izquierda). ................................................................................ 38 Figura 3.4. Bombas de baja presión (superior izquierda), alta presión (superior derecha), dosificación (inferior izquierda) y recirculación (inferior derecha). ..................................... 39 Figura 3.5. Membranas de ósmosis inversa. ......................................................................... 40 Figura 3.6. Filtro ultravioleta. ............................................................................................... 41 Figura 3.7. Filtro multimedia y filtro carbón ........................................................................ 42 Figura 3.8. Filtro cartucho de 5μm. ...................................................................................... 43 Figura 3.9. Controlador lógico programable. ........................................................................ 44 Figura 3.10. Transmisor de presión. ..................................................................................... 44 Figura 3.11. Transmisor de pH. ............................................................................................ 45 Figura 3.12. Transmisor de flujo. .......................................................................................... 46 Figura 3.13. Transmisor de conductividad. ........................................................................... 47 Figura 3.14. Transmisor de temperatura. .............................................................................. 48 Figura 3.15. Variador de frecuencia...................................................................................... 48 Figura 3.16. Válvula de control. ........................................................................................... 49 Figura 3.17. Transmisor de pH y conductividad. .................................................................. 50 Figura 3.18. Esquema de identificación. ............................................................................... 52 Figura 3.19. Flujo del proceso de identificación. .................................................................. 53 Figura 3.20. Respuesta del flujo y de la conductividad al escalón de presión. ..................... 59 Figura 3.21. Respuesta de la conductividad al escalón de pH. ............................................. 59 Figura 3.22. Respuesta del flujo volumétrico a la señal PRBS de presión. .......................... 62 Figura 3.23. Respuesta de la conductividad a la señal PRBS de presión. ............................. 63 Figura 3.24. Respuesta de la conductividad a la señal PRBS de pH. ................................... 63 Figura 3.25. Diagrama de la estructura ARX (C=1) y ARMAX (C≠1). ............................... 64 Figura 3.26. Valor FPE en la estimación del modelo para 𝑔11. ............................................ 71 Figura 3.27. Valor FPE en la estimación del modelo para 𝑔21. ............................................ 72 Figura 3.28. Valor FPE en la estimación del modelo para 𝑔22. ............................................ 72 Figura 3.29. Validación cruzada del modelo ARMAX [2 2 2] con los datos adquiridos para la estimación del modelo 𝑔11. ............................................................................................... 73 Figura 3.30. Validación cruzada del modelo ARMAX [2 2 2] con los datos adquiridos para la estimación del modelo 𝑔21. ............................................................................................... 73 Figura 3.31. Validación cruzada del modelo ARMAX [2 2 2] con los datos adquiridos para la estimación del modelo 𝑔22. ............................................................................................... 73 Figura 3.32. Polos y ceros del modelo estimado para 𝑔11. ................................................... 74 Figura 3.33. Polos y ceros del modelo estimado para 𝑔21. ................................................... 74 Figura 3.34. Polos y ceros del modelo estimado para 𝑔22. ................................................... 74 Figura 3.35. Respuesta al escalón unitario del modelo estimado para 𝑔11. .......................... 75 Figura 3.36. Respuesta al escalón unitario del modelo estimado para 𝑔21. .......................... 75 Figura 3.37. Respuesta al escalón unitario del modelo estimado para 𝑔22. .......................... 75 Figura 3.38. Modelo con estructura P-canónica.................................................................... 78 Figura 3.39. Control descentralizado de un sistema multivariable. ...................................... 78 Figura 3.40. Emparejamiento de variables. ........................................................................... 82 Figura 3.41. Sistema desacoplado. ........................................................................................ 82 Figura 3.42. Sistema con desacoplador simplificado. ........................................................... 83 Figura 3.43. Adicción de polo para forzar la realizabilidad del desacoplador. ..................... 84 Figura 3.44. Diagrama del sistema con desacoplo en Simulink. .......................................... 84 Figura 3.45. Respuesta del sistema para un escalón de presión de 40 psi. ........................... 85 Figura 3.46. Aproximación de modelos. ............................................................................... 86 Figura 3.47. Estructura del sistema de control propuesto. .................................................... 88 Figura 3.48. Respuesta del flujo del permeado. .................................................................... 88 Figura 3.49. Respuesta de la conductividad eléctrica del permeado. .................................... 89 Figura 3.50. Respuesta del flujo del permeado. .................................................................... 90 Figura 3.51. Respuesta de la conductividad eléctrica del permeado. .................................... 90 Figura 4.1. Diagrama de flujo del sistema de control. .......................................................... 93 Figura 4.2. Diagrama de flujo del controlador (subrutina). .................................................. 94 Figura 4.3. Hardware del sistema de control. ....................................................................... 95 Figura 4.4. Interfaz del software RSLogix 5000. .................................................................. 96 Figura 4.5. Interfaz del software FactoryTalk View Studio ME. .......................................... 97 Figura 4.6. Diagrama funcional de la planta piloto de desalinización. ................................. 97 Figura 4.7. Simulación de la planta desalinizadora de agua de mar por ósmosis inversa. .... 98 LISTA DE TABLAS Tabla 1.1. Energía y costos de producción de agua desalinizada por tecnología. .................. 4 Tabla 1.2. Tabla de aplicación de cada proceso de desalinización. ........................................ 5 Tabla 2.1. Diferentes modelos de plantas RO en la literatura. .............................................. 18 Tabla 3.1. Tabla comparativa de propuestas de instrumentación¹. ....................................... 32 Tabla 3.2. Líneas de la planta piloto. .................................................................................... 34 Tabla 3.3. Componentes principales del sistema de ósmosis inversa. .................................. 37 Tabla 3.4. Componentes principales del sistema de control. ................................................ 37 Tabla 3.5. Especificaciones técnicas de los tanques. ............................................................ 38 Tabla 3.6. Especificaciones técnicas de las bombas. ............................................................ 39 Tabla 3.7. Especificaciones técnicas de las membranas. ...................................................... 40 Tabla 3.8. Especificaciones técnicas del filtro ultravioleta. .................................................. 41 Tabla 3.9. Especificaciones técnicas del filtro multimedia y filtro carbón. .......................... 42 Tabla 3.10. Especificaciones técnicas del filtro cartucho 5μm. ............................................ 43 Tabla 3.11. Especificaciones técnicas del controlador. ......................................................... 43 Tabla 3.12. Especificaciones técnicas del transmisor de presión. ......................................... 44 Tabla 3.13. Especificaciones técnicas del sensor de pH ....................................................... 45 Tabla 3.14. Especificaciones técnicas del transmisor de flujo. ............................................. 46 Tabla 3.15. Especificaciones técnicas del sensor de conductividad. .................................... 47 Tabla 3.16. Especificaciones técnicas del transmisor de temperatura. ................................. 47 Tabla 3.17. Especificaciones técnicas del variador de frecuencia. ....................................... 48 Tabla 3.18. Especificaciones técnicas de la válvula de control. ........................................... 49 Tabla 3.19. Especificaciones técnicas del transmisor de pH y conductividad. ..................... 50 Tabla 3.20. Rango de linealidad del proceso de desalinización. ........................................... 58 Tabla 3.21. Características del proceso. ................................................................................ 60 Tabla 3.22. Determinación del periodo básico de la PRBS. ................................................. 62 Tabla 3.23. Órdenes de los polinomios A(q), B(q) y C(q) de los modelos candidatos. ........ 67 Tabla 3.24. Parámetros estimados del polinomio 𝐴(𝑞) para el modelo 𝑔11. ........................ 67 Tabla 3.25. Parámetros estimados del polinomio 𝐵(𝑞) para el modelo 𝑔11. ........................ 67 Tabla 3.26. Parámetros estimados del polinomio 𝐶(𝑞) para el modelo 𝑔11. ........................ 68 Tabla 3.27. Parámetros estimados del polinomio 𝐴(𝑞) para el modelo 𝑔21. ........................ 68 Tabla 3.28. Parámetros estimados del polinomio 𝐵(𝑞) para el modelo 𝑔21. ....................... 68 Tabla 3.29. Parámetros estimados del polinomio 𝐶(𝑞) para el modelo 𝑔21. ........................ 69 Tabla 3.30. Parámetros estimados del polinomio 𝐴(𝑞) para el modelo 𝑔22. ........................ 69 Tabla 3.31. Parámetros estimados del polinomio 𝐵(𝑞) para el modelo 𝑔22. ....................... 69 Tabla 3.32. Parámetros estimados del polinomio 𝐶(𝑞) para el modelo 𝑔22. ........................ 69 Tabla 3.33. Tabla de resultados para la estimación del modelo 𝑔11. .................................... 70 Tabla 3.34. Tabla de resultados para la estimación del modelo 𝑔21. .................................... 71 Tabla 3.35. Tabla de resultados para la estimación del modelo 𝑔22. .................................... 71 Tabla 3.36. Fórmulas de Skogestad para la sintonización del PID. ...................................... 87 Tabla 3.37. Valores de sintonización utilizados para el proceso estudiado. ......................... 87 Tabla 3.38. Resultados del sistema de control. ..................................................................... 91 INTRODUCCIÓN El agua es esencial para la salud, el bienestar, el desarrollo económico así como para la preservación del medio ambiente. Sin embargo, de acuerdo a la Organización de las Naciones Unidas, dos de cada diez personas en el mundo no tienen acceso a este recurso. La escasez del agua, su baja calidad y un saneamiento deficiente ocasionan millones de muertes cada año, afectando la seguridad alimenticia, las opciones de sustento y las oportunidades de educación. Del total de agua disponible en la tierra sólo un 2.5% se presenta en forma de agua dulce y es aprovechado en mayor proporción, mientras que la mayor cantidad de ésta se encuentra en los océanos. Es decir, nuestro principal recurso hídrico presenta una alta salinidad y por ello la desalinización del agua de mar es una importante alternativa. Sin embargo, los costos relacionados al consumo de energía de las plantas desalinizadoras restringen el crecimiento de su uso, en especial en países como el Perú, por lo cual, es indispensable el desarrollo de técnicas de control que permitan abaratar costos en los procesos de estas plantas. El presente trabajo de tesis expone el estudio de las planta de desalinización por ósmosis inversa, tecnología que ha ganado terreno en el ámbito de la desalinización debido a la constante mejora en la tecnología de membranas semipermeables para altos flujos a presiones relativamente bajas. En el primer capítulo se describe el estado del arte de los sistemas de desalinización de agua de mar, subrayando la importancia de la tecnología de osmosis inversa frente a la las demás. Se describe además el flujo del proceso de desalinización en una planta industrial con esta tecnología, esquematizando todas las etapas que intervienen en la misma. En el segundo capítulo se muestran los modelos más destacados en la literatura y se realiza una comparación entre las respuestas al escalón, estudiando su comportamiento para luego seleccionar un modelo que sirva para el diseño del controlador. Luego se mencionan las distintas estrategias de control aplicadas en los sistemas de desalinización por ósmosis inversa. 1 En el tercer capítulo se expone el diseño de una planta piloto, así como su implementación en el laboratorio de control y automatización de la PUCP mediante una evaluación técnica y económica. Se realiza la descripción de las características generales de la planta piloto, además de las características técnicas de los equipos y de la instrumentación que la componen. Luego, se realiza la identificación de la planta piloto utilizando la instrumentación mencionada y por último se diseña el controlador utilizando la estrategia de desacoplamiento de variables para la posterior aplicación de un controlador convencional. En el cuarto capítulo se elabora una propuesta de implementación práctica del sistema de control, describiendo el algoritmo del mismo considerando los modos de operación de la planta piloto. Se exhibe además el flujo de información en el hardware de control. Se desarrolla el código correspondiente utilizando el software del fabricante del controlador, y se elabora una interfaz gráfica para la interacción con el sistema de control. Se realizan las pruebas de simulación utilizando una rutina con el comportamiento dinámico de la planta. Finalmente, se realizan las conclusiones y recomendaciones, incorporándose además las referencias bibliográficas y los anexos que se utilizaron para el desarrollo del presente trabajo de tesis. 2 CAPÍTULO 1. ESTUDIO DE LAS PLANTAS DESALINIZADORAS DE AGUA DE MAR. 1.1. Estado del Arte Aproximadamente 97.5% del agua en nuestro planeta está situado en los océanos, mientras que, del 2.5% que representa el agua dulce, aproximadamente el 70% se encuentra en forma de glaciares y nieve, y el 30% es agua subterránea, agua de ríos y lagos, y humedad de aire (Voutchkov, 2013). Por tal motivo, la desalinización de agua de mar es la vía para aprovechar la principal fuente de agua, y la cual típicamente se realiza mediante dos tipos de tecnologías, térmica y de separación por membranas. La primera comprende las técnicas de destilación flash multietapa (MSF), destilación de múltiple efecto (MED) y la destilación con compresión de vapor (VC), en tanto que, la segunda contiene a las técnicas de electrodiálisis (ED) y ósmosis inversa (RO). La desalinización a escala industrial, tiene sus inicios en la primera parte del siglo XX; sin embargo, la expansión y propagación de la industria desalinizadora ocurre en el periodo de 1960 a 1980, donde la tecnología utilizada era la desalinización térmica, más precisamente, la técnica de destilación flash multietapa. Durante este periodo el estándar de la industria consistía en plantas desalinizadoras con 24 etapas y una capacidad de producción de 27 276 m³ por día (Alatiqi, et al., 1999). Esta técnica se basa en el principio de que al reducir abruptamente la presión de agua de mar por debajo del valor de su presión de equilibrio, ocurre una evaporación súbita o una ebullición explosiva de la misma. En 1980, se aplica la técnica de la destilación de múltiple efecto que utiliza el mismo principio que el proceso anterior, radicando su principal diferencia en la 3 forma en que se lleva a cabo la evaporación. En las plantas con destilación de múltiple efecto se utilizan varios evaporadores del tipo película delgada, con los cuales se logran mejores coeficientes de transferencia de calor que los que se pueden obtener en las plantas con destilación flash multietapa, donde se produce la evaporación súbita en forma directa (Ramillo, et al., 2003). En ese mismo año, se aplica la técnica de destilación por compresión de vapor, el cual provoca la condensación del agua en una superficie que permite la transferencia de calor del condensado hacia la salmuera ubicada en el lado opuesto de dicha superficie, resultando en la vaporización. Aquí el compresor es esencial para aumentar la presión en el lado de vapor y reducir la presión en el lado del agua de alimentación, o salmuera, para bajar su temperatura de ebullición (Cotruvo, et al., 2010). Durante el periodo de 1980 a 1999 la técnica de ósmosis inversa (RO) empieza a ganar terreno en el ámbito de la desalinización de agua debido a la constante mejora en la tecnología de membranas semipermeables para altos flujos, operación a presiones relativamente bajas y capaces de soportar las duras condiciones del agua de mar de alimentación y salmuera. Esta técnica ha evolucionado más rápidamente que cualquier otra tecnología de desalinización, principalmente debido a su competitivo consumo de energía y a los costos de producción de agua los cuales se aprecian en la Tabla 1.1 (Voutchkov, 2013). Tabla 1.1. Energía y costos de producción de agua desalinizada por tecnología. Energía / Costo MSF MED VC RO Energía utilizada (kWh/m³) 5.7 – 7.8 12.7 – 15.0 8.0 – 12.0 2.5 – 4.0 Costo de producción (USD/m³) 0.7 – 3.5 0.9 – 4.0 1.0 – 3.5 0.5 – 3.0 La ósmosis inversa es una técnica para extraer sólidos disueltos del agua, tales como las sales, utilizando una membrana semipermeable, la cual posee una permeabilidad alta respecto del agua y muy baja para las sales. A diferencia de las técnicas anteriores, este proceso no involucra ningún cambio de fase. En este proceso, el agua pasa a través de la membrana impulsada por una bomba que eleva su presión hasta un valor superior al de su presión osmótica natural. Para este 4 propósito se utiliza típicamente una bomba de alta presión de 54 a 82 bar (Ramillo, et al., 2003). Por último, en los sistemas de tratamiento basados en electrodiálisis, una corriente continua se hace pasar a través del agua, impulsando los iones a través de las membranas a los electrodos de carga opuesta. Estos sistemas poseen gran número, entre 300 y 600 pares, de membranas de intercambio catiónico y aniónico separados por divisores de flujo para evitar que se peguen entre sí y transmitir el flujo de agua desalinizada a través y fuera de las membranas. Esta última técnica puede tratar aguas con alta turbiedad y bioensuciamiento mejor que los sistemas de ósmosis inversa, sin embargo a un costo mayor. La Tabla 1.2, muestra el rango de salinidad de agua para el cual las técnica de destilación (MSF, MED o VC), ósmosis inversa y electrodiálisis pueden ser aplicadas de manera rentable en la desalinización de agua (Voutchkov, 2013). Tabla 1.2. Tabla de aplicación de cada proceso de desalinización. Rango de salinidad de agua para Tecnología una aplicación rentable (mg/L) Destilación 20 000 – 100 000 Osmosis Inversa 50 – 46 000 Electrodiálisis 200 – 3 000 Dado que nuestro territorio limita con océano pacífico, el cual posee una salinidad de aproximadamente 45 000 mg/L es recomendable el uso de la tecnología de ósmosis inversa. Por ello, se describe a continuación el proceso llevado a cabo en una planta desalinizadora de agua de mar en la que se utiliza dicha tecnología. 1.2. Descripción de una Planta Desalinizadora de Agua de Mar por Ósmosis Inversa En general, cualquier fuente natural de agua presenta sólidos en dos formas, estas son suspendidas y disueltas. Los sólidos en suspensión existen como forma de partículas insolubles, como por ejemplo escombros, organismos marinos, limo, coloides, etc., mientras que los sólidos disueltos se presentan en forma soluble, como el caso de los iones de minerales tales como cloruro, sodio, calcio, magnesio, etc. 5 En la actualidad, prácticamente todas las plantas de desalinización por ósmosis inversa incorporan dos etapas para la remoción secuencial de los sólidos en suspensión y sólidos disueltos de las fuentes de agua. Para una clara comprensión del proceso, podemos ilustrar el diagrama de la planta Llobregat situada en la ciudad de Barcelona, España, en la Figura 1.1. Sin embargo, la descripción dada corresponde a la mayoría de plantas desalinizadoras que utilizan esta tecnología. El proceso inicia con la captación de agua de mar a través de tomas a mar abierto utilizando sistemas de bombeo. Aquí se usan filtros en el ingreso del agua de mar hacia los conductos para evitar el ingreso de partículas de gran tamaño u organismos marinos, así como oxidantes para evitar su crecimiento en las tuberías de conducción. Una vez captada el agua comienza la etapa de pretratamiento, la cual consiste inicialmente en un procedimiento de remoción de arena, sedimentación y flotación por aire disuelto para minimizar el contenido de materiales gruesos y para proteger de la sobrecarga de sólidos a las instalaciones de filtración aguas abajo. El proceso continúa con la filtración en medios granulares, o filtración convencional, donde se filtra la fuente de agua a través de una o más capas de medios granulares (por ejemplo, carbón de antracita, arena de sílice, granate). Los filtros utilizados para el tratamiento del agua salina son típicamente unidades bicapa y se componen de una capa gruesa de carbón de antracita sobre una cama de arena fina (Cipollina, et al., 2009). Para mejorar el proceso de filtración se utiliza algunos componentes químicos como coagulantes y floculantes mientras que para suprimir la cristalización de partículas minerales en proceso de formación en la superficie de la membrana de ósmosis inversa se utilizan anti incrustantes. Sin embargo, existen casos en que el agua contiene altos niveles de compuestos orgánicos (concentración de carbono orgánico total mayor a 6 mg/L) y de sólidos en suspensión (turbidez promedio mensual excede los 20 NTU) y por ello se aplican sistemas de filtración de dos etapas. En esta configuración, la primera etapa de filtración está diseñada principalmente para eliminar los sólidos gruesos y orgánicos en forma suspendida. Los filtros de la segunda etapa están configurados para retener los sólidos finos y limos y para retirar una porción (20 a 40 por ciento) de los compuestos orgánicos solubles contenidos en el agua salina por biofiltración. 6 7 REMINERALIZACIÓN ÓSMOSIS INVERSA IMPULSIÓN A (SEGUNDO PASO) DEPÓSITO ABIERTO ÓSMOSIS INVERSA (PRIMER PASO) DEPÓSITO Y BOMBEO DE AGUA PRODUCTO FILTRACIÓN TANCADA (FILTROS BICAPA) TRATAMIENTO DE EFLUENTES BOMBEO EN FILTRACIÓN MICROFILTRACIÓN RECHAZO CERRADA (FILTROS CARTUCHO) EMISARIO EDAR LLOBREGAT OBRA DE LLEGADA EVACUACIÓN DE FANGOS FILTRACIÓN ABIERTA BOMBEO DE (FILTROS BICAPA) FLOTACIÓN AGUA DE MAR CAPTACIÓN DE AGUA DE MAR Figura 1.1. Esquema del proceso de desalinización de agua de mar por ósmosis inversa de la planta de Llobregat, Barcelona, España. El agua filtrada es transportada mediante bombas de baja presión para vencer la resistencia de los últimos elementos del pretratamiento que son los filtros con tamaño de pasajes en micrones o filtros cartucho. En esta tubería de conducción se adiciona soluciones ácidas y básicas para la regulación del pH, así como oxidantes para evitar el ensuciamiento de la membrana. Comienza así la etapa de tratamiento, con el incremento de la presión del agua de alimentación la cual debe superar la presión osmótica del agua de mar. La magnitud de esta presión depende de la ubicación geográfica y generalmente se encuentra entre 55 a 70 bares (Voutchkov, 2013). Para superar la presión que garantiza la separación de las sales de la fuente de agua salina se utilizan diversos tipos de bombas, entre los que se encuentran las bombas reciprocantes o de desplazamiento positivo, bombas centrífugas, bombas turbina vertical, bombas de caja partida y bombas multietápicas de sección anular. El agua salina con elevada presión es entregada al sistema de ósmosis inversa. Este sistema está conformado por bastidores que operan de forma independiente y en los cuales se montan los recipientes de ósmosis inversa. Dichos recipientes contienen membranas con configuraciones comerciales de tipo enrollamiento en espiral, fibra hueca, placas, y tubulares, donde cada una tiene sus propias características. Las configuraciones de fibra hueca y enrollamiento en espiral generalmente tienen características de operación más favorables de rendimiento en relación con el costo, y se utilizan con mayor frecuencia (Cotruvo, et al., 2010). La Figura 1.2 muestra el tipo de membrana utilizado en el presente trabajo de tesis, no obstante el principio es el mismo para las demás configuraciones. Figura 1.2. Membrana de ósmosis inversa de tipo enrollamiento en espiral. 8 La presión del agua de alimentación aplicada contrarresta la presión osmótica y supera las pérdidas de presión que se producen cuando el agua viaja a través de la membrana, manteniendo así el agua dulce en el lado de baja salinidad de la membrana, hasta que ésta agua sale del recipiente de membrana. Las sales contenidas en el lado de la alimentación de la membrana son retenidas y elevan su concentración de sales. Como resultado se obtienen dos corrientes, la primera de agua con baja salinidad (producto o permeado) a baja presión y la segunda, agua con elevada salinidad (rechazo, concentrado o salmuera) a alta presión. Es posible aprovechar la energía hidráulica del concentrado y suministrarla en forma de presión adicionándola a la alimentación, esto reduciría la energía requerida para el sistema de ósmosis inversa, por lo que, en operación la bomba de alta presión trabajaría como una bomba de mediana presión. Luego de ello, el rechazo a baja energía se reúne junto a parte de los efluentes del sistema de filtración, previamente tratados, para ser retornados al mar considerando el impacto medioambiental. Los residuos del tratamiento de estos efluentes son retirados de la planta desalinizadora para ser dispuestos como relleno sanitario. Por otro lado, el agua producto puede ser llevada hacia otra etapa de tratamiento por ósmosis inversa dependiendo del requerimiento de la calidad de agua, de otro lado seguirá su curso hacia la etapa de postratamiento. Es necesario resaltar que mientras las membranas semipermeables rechazan todos los sólidos suspendidos, no son una barrera absoluta para los sólidos disueltos (minerales y orgánicos). Algunos sólidos disueltos acompañarán el paso de agua dulce a través de la membrana, por lo que la tasa de agua y paso de sales son las dos características clave del rendimiento de las membranas de ósmosis inversa (Voutchkov, 2013). El permeado generado por los bastidores de ósmosis inversa es enviado a la etapa de postratamiento, donde es estabilizado mediante la adición de cal, o con el contacto con calcita, y la adición de dióxido de carbono para proporcionar un nivel adecuado de alcalinidad y dureza para la protección del sistema de distribución de agua producto contra la corrosión. Finalmente el agua es almacenada, acondicionada y desinfectada antes de su entrega a los usuarios finales. 9 CAPÍTULO 2. MODELADO Y CONTROL DE PLANTAS DESALINIZADORAS POR OSMOSIS INVERSA. 2.1. Modelado de Plantas Desalinizadoras de Agua de Mar por Ósmosis Inversa 2.1.1. Modelado Matemático y Modelado Mediante Identificación de Sistemas Existen, en la literatura, dos enfoques en el desarrollo de modelos para el proceso de ósmosis inversa. El primer enfoque es el modelado matemático, donde a través de leyes físicas y mediante ecuaciones de balance se llega a expresar una función matemática que describa el comportamiento dinámico de la planta. El segundo enfoque es la identificación de sistemas, donde se presenta el modelado como una caja negra en la cual se desarrollan funciones de transferencia o ecuaciones en el espacio de estados que relacionan las variables de entrada y las variables de salida utilizando para ello señales de entrada especiales, como la señal impulso, escalón o binarias pseudoaleatorias. El presente capítulo expone el trabajo conjunto realizado sobre el tema (Carrasco, et al., 2014). Uno de los primeros estudios fue realizado por Alatiqi (Alatiqi, et al., 1989), quien realizó el modelado de una planta piloto de desalinización por ósmosis inversa en Kuwait. Este procedimiento consistió en la identificación del comportamiento dinámico que gobierna el proceso de desalinización, donde debido a que éste presenta varias variables de entrada y salida se utilizó un enfoque multivariable. El modelo consideró como variables manipuladas la presión, que fue manejado mediante el estrangulamiento de la válvula de rechazo, y el pH del alimentado, 10 mediante la adición de una solución ácida. Se tomó como variables controladas el flujo volumétrico de permeado y, debido a que la calidad de agua desalinizada depende de la conductividad eléctrica de la misma, ésta fue la segunda variable de salida. Como resultado se obtuvo una expresión general para la representación del sistema dado por (2.1). 𝒀 = 𝑮𝑼 (2.1) 𝐹 𝑔11 𝑔12 𝑃 [ ] = [𝑔 𝑔 ] [ ] (2.2) 𝐶 21 22 𝑝𝐻 Donde 𝐹 es el flujo volumétrico de permeado, 𝐶 es la conductividad eléctrica del permeado, 𝑃 es la presión de la alimentación, y 𝑝𝐻 es el pH de la alimentación. Las ecuaciones finales en el dominio de Laplace se expresan en (2.3), (2.4), (2.5) y (2.6). 𝐹 0.002(0.056𝑠 + 1) = 𝑔 𝑃 11 = (2.3) (0.003𝑠2 + 1𝑠 + 1) 𝐹 = 𝑔 = 0 (2.4) 𝑝𝐻 12 𝐶 −0.51(0.35𝑠 + 1) = 𝑔21 = (2.5) 𝑃 (0.213𝑠2 + 7𝑠 + 1) 𝐶 −57(0.32𝑠 + 1) = 𝑔22 = (2.6) 𝑝𝐻 (0.6𝑠2 + 1.8𝑠 + 1) Estas expresiones se obtuvieron mediante un procedimiento el cual inició con el sistema en estado estacionario. Luego se aplicó una señal escalón a la primera variable de entrada sin variar la segunda y se observó la respuesta de las variables controladas, repitiendo el procedimiento con la variación de la segunda variable de entrada. Posteriormente se procesaron los datos resultantes y se aproximaron a funciones de transferencia de segundo grado con un cero en el numerador. Robertson (Robertson, et al., 1996) y Abbas (Abbas, 2006) utilizaron este modelo para el diseño de un controlador predictivo DMC, mientras que Assef (Assef, et al., 1997) utilizó la misma técnica de identificación de parámetros que planteaba Alatiqi para modelar una planta con cuatro membranas de celulosa, obteniendo también 11 funciones de transferencia de segundo orden, con la particularidad de presentar un retardo de 120 segundos en la variable de conductividad eléctrica debido a la presión y pH de la alimentación. Burden (Burden, et al., 2001) propuso una variación del modelo anterior considerando las mismas variables para una planta piloto con membranas de fibra hueca, y excitando el sistema mediante una señal escalón para la identificación de la misma. El modelo obtenido fue utilizado para el desarrollo del controlador CMPC. Riverol y Pilipovik (Riverol & Pilipovik, 2005) aplicaron la identificación de sistemas y llegaron a un modelo con desacoplamiento perfecto, considerando las mismas variables que Alatiqi. En general, los modelos han sido forzados a una respuesta con una determinada forma, esto es, una función de segundo grado expresado en (2.7). 𝐾𝑖(𝜏𝑖𝑠 + 1) 𝑔𝑖 (2.7) ( = 𝑠) 𝜏 2𝑠2 𝑖 + 2𝜁𝑖𝜏𝑖𝑠 + 1 En lo que respecta al modelado matemático de las plantas desalinizadoras, Gambier (Gambier, et al., 2007) hace hincapié en el hecho que si bien en la literatura actual encontramos modelos en tiempo-espacio y modelos generados por la identificación de parámetros, no hay mucha información al respecto acerca de modelos dinámicos de parámetros concentrados obtenidos por la aplicación de leyes físicas en una planta desalinizadora. La importancia de obtener estos modelos es la de implementar sistemas tolerantes a fallos y conocer sobre las respuestas transitorias y estables del sistema. Es por ello que la mayor contribución y el objetivo de ese trabajo es presentar un modelo de parámetros concentrados basados en leyes físicas. El proceso de modelado se dividió en varios pasos. El primero fue dividir la planta en subsistemas, considerando que usualmente una planta desalinizadora de agua de mar por osmosis inversa se divide en cinco subsistemas: generación de energía, pretratamiento del agua, bombeo y recuperación de energía, unidad de RO y pos- tratamiento (Figura 2.1). 12 sustancias agua de Unidad de Desalinización Unidad sustancias químicas mar químicas Eléctrica Subsistema de Subsistema de Subsistema de Subsistema de Bombeo con Subsistema de Generación Desalinización Pretratamiento recuperador de Postratamiento de Potencia (Unidades RO) energía suministro de energía eléctrica Figura 2.1. Esquema de los subsistemas de una planta de RO. Los sistemas de pretratamiento y postratamiento fueron modelados como tanques, por lo que el enfoque se centró en establecer el modelo para la membrana de ósmosis inversa. La unidad de ósmosis inversa a su vez se dividió en los subsistemas de salmuera, permeado y la membrana, y para la obtención del modelo se asumió lo siguiente:  La membrana responde a un sistema de filtrado ideal a través de difusión.  Los flujos dentro de la unidad de ósmosis inversa son laminares.  La unidad de ósmosis inversa es alimentada con agua. Para el subsistema de salmuera se describió el modelado con ecuaciones para la masa, flujo, momento, temperatura, energía y concentración de sal. Se realizó además el modelado para el subsistema de permeado de igual manera. (Figura 2.2). Flujo constante de Flujo de permeado agua de alimentación Unidad de Concentración del permeado Ósmosis Temperatura variable Inversa del agua de alimentación Flujo de salmuera Válvula de apertura variable Figura 2.2. Esquema de simulación de la planta de RO. Para el modelado de la membrana se utilizó ecuaciones especiales usando la membrana de fibra hueca, luego se simuló y se realizó experimentos en el sistema en lazo abierto y haciendo uso de los parámetros de la membrana. El primer experimento consistió en analizar el flujo de permeado y la concentración ante una determinada temperatura (18°C) y a diferentes porcentajes de apertura de la válvula (0-100%), la conclusión es que a menor porcentaje de apertura de la válvula 13 disminuye la sal concentrada en el flujo de permeado lo que significa un agua de mayor calidad y que necesita menor postratamiento. El segundo experimento consistió en mantener constante la válvula al 50% e incrementar la temperatura de 18°C a 23°C lo que resultó que haya más agua de permeado. Esto provoca un cambio en los coeficientes de la membrana y la concentración de permeado crece y luego decae. Las ecuaciones que se propusieron muestran la relación del flujo y la concentración respecto a la temperatura. Finalmente, en el tercer experimento, la apertura de la válvula y la temperatura se cambian dinámicamente. Cuando la apertura de la válvula cae de 100 a 50% el flujo de permeado se incrementa linealmente mientras que la concentración de sal decrece exponencialmente. Un similar efecto se observa cuando se incrementa la temperatura del agua que ingresa de 18°C a 24°C. Como conclusión en este trabajo, se obtiene que el modelo matemático propuesto se ajusta a la respuesta de una planta real al igual que los modelos propuestos en la literatura, con excepción de la necesidad de ajustar algunos parámetros propios de la planta real. Esto fue resultado importante, pues demostró la posibilidad de estimar más modelos matemáticos de igual relevancia. Otro ejemplo de modelado matemático es presentado por Al-haj Ali (Al-haj Ali, et al., 2008), en cuyo estudio se desarrolló un modelo dinámico sencillo para controlar la calidad del agua utilizando ósmosis inversa generando un modelo teórico donde cada membrana es descrita por dos ecuaciones diferenciales ordinarias utilizando para ello leyes físicas de conservación. El modelo consideró las suposiciones siguientes:  La solución contiene sólo un soluto y un solvente (solución binaria).  La permeabilidad del agua es constante e independiente de la presión.  El flujo dentro del tubo es turbulento.  La polarización de la concentración en la superficie de la membrana es descrita por el modelo de película de Nerst.  La presión osmótica puede ser representado por la ecuación de Van’t Hoff.  En el lado de permeado la mezcla es uniforme en la dirección radial.  Se desprecian los efectos de difusión axial y angular.  La caída de presión radial es despreciable. 14 Para describir el modelo se utilizaron dos ecuaciones implícitas las cuales se resuelven iterativamente para conseguir el flujo de permeado y la concentración en la superficie de la membrana. Luego se obtiene la concentración de sales de la salmuera así como la presión a la que se encuentra. Por último, se halla el flujo de salmuera y se calcula la concentración de sales del agua producto. Este modelo es probado en un módulo con condiciones de operación de temperatura ambiente de 30° C, presión de alimentación de 35 bar, concentración de la solución de 2000 ppm y velocidad de entrada de 6.8 m/s. La membrana utilizada tuvo un diámetro de 1.25 cm y una longitud de 60cm. Al aplicar el modelo en la experimentación de este módulo se obtuvo una buena correlación frente al mismo, lo cual indica que existe una dependencia lineal de la presión de salida con respecto a la velocidad de entrada. Asimismo se midió la concentración axial en diferentes posiciones de la membrana para contrastarse con el modelo. Se aplicó una entrada escalón de 0.5 a 0.55 m³/s en el flujo de alimentación resultando en una disminución de la presión de alimentación. En consecuencia, la concentración de sales en la superficie de la membrana disminuyó y la presión de salida (salmuera) aumentó. Sin embargo, esta experimentación debe modificarse para aplicaciones industriales dado que la alimentación sólo fue agua levemente salada. McFall (McFall, et al., 2008) desarrolló un modelo matemático detallado de la planta, para luego aplicar un controlador no lineal que incluye componentes feedback y feed-forward. El modelo de la planta está basado en los principios básicos de conservación: balance microscópico de energía cinética, balance local de materia y balance microscópico de materia en un volumen de control, y para ello se consideró lo siguiente:  El fluido es incompresible,  La masa y el volumen interno son constantes,  El agua en el módulo viaja en flujo pistón (no hay retromezcla ni difusión axial), 15  La presión osmótica puede relacionarse con el total de sólidos disueltos en la superficie de la membrana.  La fricción en las paredes de las membranas y tuberías es despreciable con respecto a las pérdidas hidráulicas a través de las válvulas y membranas. El balance de materia se realizó alrededor de las válvulas de control manipuladas, el balance local de materia en el punto de derivación del bypass y para el balance de envoltura se tomó un volumen de control diferencial en la dirección del flujo. La solución del modelo se realizó iterativamente resolviendo las ecuaciones diferenciales del balance de envoltura utilizando el método de disparo. La síntesis del controlador se realizó en base a una perturbación originada en la concentración de la alimentación. Teóricamente, el modelo puede ser mejorado considerando la variación de la presión en la dirección del flujo de alimentación dentro de los bastidores. Es posible usar ecuaciones diferenciales parciales para describir los perfiles de velocidad y concentración en cualquier punto, asimismo tomar en cuenta el gradiente de concentración en la dirección radial y axial dentro de las membranas y finalmente incluirlo en el modelo la dependencia de la presión osmótica con la temperatura. Por otra parte, un modelado específico fue presentado por Syafiie (Syafiie, et al., 2008), cuyo objetivo fue modelar el comportamiento dinámico de los componentes del subsistema de pretratamiento (filtro de arena y cartuchos de asentamiento). Los parámetros que se utilizaron en las ecuaciones del modelado del filtro de arena fueron el caudal de ingreso y de salida, el diferencial de masa del agua que pasa a través de ellos y de las partículas suspendidas. Se modeló además, la cantidad de bacteria viviente y componentes inorgánicos en la solución de agua salada, así como su relación y reacción con los componentes químicos adicionados. Por último, se simuló el sistema usando el software EcoSimPro. Así también se han obtenido modelos matemáticos con fines más específicos, como por ejemplo el desarrollado por Bartman (Bartman, et al., 2009), quien propuso un algoritmo de control predictivo basado en modelo para controlar la inversión del flujo en un sistema de desalinización por ósmosis inversa. Basado en el hecho de que es fundamental evitar la formación de incrustaciones en las paredes de las 16 membranas y en que los métodos convencionales, adición de agentes anti- incrustantes y flushing, presentan muchos inconvenientes y desventajas respecto al método de inversión de flujo, desarrolló un modelo para el control de la inversión de flujo, con lo cuales buscó determinar la trayectoria óptima entre las condiciones normales de operación y un nuevo estado estacionario utilizando el control predictivo basado en modelo. El desarrollo del modelo se basó en un balance global de materia y en balances locales de energía alrededor de las válvulas del sistema, y luego se tomó la concentración efectiva en las unidades de ósmosis inversa como un promedio ponderado de las concentraciones de la alimentación y rechazo. Se utilizó además, una ecuación algebraica para relacionar la presión del sistema con los demás parámetros y una ecuación para la presión osmótica basada en la concentración y temperatura. 2.1.2. Comparación de Modelos de Plantas Desalinizadoras de Agua de Mar por Ósmosis Inversa Los modelos más usados para plantas desalinizadoras de la literatura actual se resumen en el trabajo de Gambier (Gambier, et al., 2007), en el cual se indica que para hallar las funciones de transferencia mostradas en la Tabla 2.1 se utilizaron experimentos en plantas desalinizadoras reales teniendo como elementos de entrada la presión en la membrana y el pH de la alimentación así como elementos de salida, el flujo volumétrico y la conductividad eléctrica del permeado. En la Figura 2.3 se muestra lo expuesto. En la Tabla 2.1, la notación [1] hace referencia al modelo propuesto por Assef (Assef, et al., 1997), la notación [2] al modelo propuesto por Riverol (Riverol & Pilipovik, 2005) y la notación [3] hace referencia al modelo propuesto por Alatiqi (Alatiqi, et al., 1989). 17 Tabla 2.1. Diferentes modelos de plantas RO en la literatura. Funciones de Transferencia Entradas Ref. Válvula de control de presión Válvula de entrada de ácido −0.155(0.375𝑠 + 1) 𝑭 𝐺11(𝑆) = 𝐺 ( 12(𝑆) = 0 0.22𝑠 + 1)(2.51𝑠 + 1) 2.48𝑒−120𝑠 0.45𝑒−120𝑠 𝑪 𝐺21(𝑆) = 𝐺 = (114𝑠 + 1)(113𝑠 + 1) 22(𝑆) (104𝑠 + 1)(100𝑠 + 1) 4.74 𝑷 𝐺31(𝑆) = − 𝐺32(𝑆) = 0 (1.45𝑠 + 1) 0.077 𝒑𝑯 𝐺41(𝑆) = 0 𝐺42(𝑆) = − (21.2𝑠 + 1) 𝑷 𝒑𝑯 0.0045(0.104𝑠 + 1) 𝑭 𝐺11(𝑆) = 𝐺12(𝑆) = 0 0.012𝑠2 + 𝑠 + 1 −0.12𝑠 + 0.22 10(−3𝑠 + 1) 𝑪 𝐺21(𝑆) = 0.1𝑠2 𝐺 + 0.3𝑠 + 1 21(𝑆) = 𝑠2 + 5𝑠 + 1 0.002(0.056𝑠 + 1) 𝑭 𝐺11(𝑆) = 𝐺 = 0 0.003𝑠2 12(𝑆) + 0.1𝑠 + 1 −0.51(0.35𝑠 + 1) −57(0.32𝑠 + 1) 𝑪 𝐺21(𝑆) = 2 𝐺 ( ) = 0.213𝑠 + 0.7𝑠 + 1 22 𝑆 0.6𝑠2 + 1.8𝑠 + 1 Presión antes de la Flujo del permeado membrana pH de entrada Conductividad del permeado Figura 2.3. Diagrama de simulación de la planta de desalinización. Se utilizó el software Simulink de Matlab para contrastar los modelos representados por sus funciones de transferencia. Los resultados mostrados fueron obtenidos por la simulación de estos modelos ante entradas de tipo escalón que representan el porcentaje de apertura de las válvulas de presión y de pH respectivamente, dando como resultado las Figura 2.4 y Figura 2.5. 18 [3] [2] [1] Salidas Figura 2.4. Respuesta del flujo volumétrico de permeado al escalón de presión en la alimentación para varios modelos. Figura 2.5. Respuesta de la conductividad eléctrica al escalón de presión y pH en la alimentación para varios modelos. De acuerdo con estos resultados, se ha seleccionado el modelo propuesto por Alatiqi (Alatiqi, et al., 1989) dado en (2.2) a (2.6) para realizar la identificación del sistema mediante pruebas de simulación. Esto debido a que en los demás modelos presentados en la Tabla 2.1 se tiene una ganancia positiva en el modelo que relacionan la presión de la alimentación y la conductividad eléctrica del permeado, así como en el modelo que relaciona el pH de la alimentación con la conductividad eléctrica del permeado. En este último modelo, además, exhibe un comportamiento de fase no mínima. 19 2.2. Control de Plantas Desalinizadoras de Agua de Mar por Osmosis Inversa 2.2.1. Revisión de las Estrategias de Control Como una parte esencial de la industria de la desalinización por ósmosis inversa, se requiere del control de procesos para la operación en condiciones óptimas, así como para el incremento del ciclo de vida de la planta. En términos generales, el objetivo de los lazos de control en un proceso de ósmosis inversa es mantener una tasa de producción constante y a la vez lograr una pureza aceptable. Inicialmente, Mindler (Mindler & Epstein, 1986) propuso un sistema de control on-off el cual tiene la desventaja de que es necesario tener una gran unidad de almacenamiento para compensar las variaciones en la demanda y evitar el frecuente encendido y apagado del sistema de desalinización. Más tarde, en 1999, Alatiqi (Alatiqi et al., 1999) realizó una revisión de los controladores y de la instrumentación utilizada en las plantas de ósmosis inversa, encontrando que los controladores comúnmente utilizados son los de tipo PID, especialmente PI, y posteriormente propuso un sistema de control multilazo, el cual incluía un controlador de presión y dos controladores de pH. De la igual manera, en esta sección se abarca la mayoría de las técnicas de control de plantas desalinizadoras, los cuales se resumen en el diagrama de la Figura 2.6. Controladores de plantas de RO Controladores Controladores clásicos avanzados PID Control Control Control no Control Robusto predictivo lineal óptimo Zieglers Algortimos CMPC DMC Nichols inteligentes Algoritmo Enjambre inmuno- de genético párticulas Figura 2.6. Técnicas de control utilizadas en plantas de desalinización. 20 2.2.2. Control Clásico El control PID es el más utilizado en la industria pero presenta la limitación de ser aplicable únicamente a plantas tipo SISO. Como se observó anteriormente, el sistema de desalinización por osmosis inversa es de tipo multivariable, por lo que aplicar este tipo de control implica el desacoplamiento de las funciones de transferencia. A continuación se presentan diversos estudios sobre la aplicación de controladores PID y cuya sintonización se compara con el resultado de la sintonización por el método de Ziegler y Nichols para plantas desalinizadoras por ósmosis inversa. Kim (Kim, et al., 2008) presentó el diseño de controladores PID para una planta de osmosis inversa, realizándose una comparación entre la sintonización por el método de Zieglers-Nichols y un algoritmo inmuno-genético. El objetivo principal del sistema de control fue mantener el pH del agua producto constante y el modelo utilizado fue el modelo de Alatiqi (Alatiqi, et al., 1989). La propuesta estuvo basada en dos partes, un algoritmo genético que se inspira en el mecanismo de la selección natural, proceso biológico en el cual el individuo más fuerte tiene más probabilidad de ser ganador en un entorno de competencia, y un algoritmo inmune que simula el funcionamiento del cuerpo humano en cuanto a la creación de anticuerpos y antígenos, con una capacidad de aprendizaje, memoria, reconocimiento de patrones y recreación de anticuerpos relevantes de forma rápida. Este algoritmo posibilitó la autorregulación del sistema, y representó una opción factible y sencilla en el diseño del controlador pues únicamente se requirió de una función objetivo. Los resultados de simulación mostraron un tiempo de estabilización de 147 segundos para el flujo y de 25 segundos para la conductividad. Sin embargo estos resultados no fueron comprobados experimentalmente. Rathore (Rathore, et al., 2013) implementó un algoritmo de optimización por enjambre de partículas al sistema de control del proceso de desalinización por 21 ósmosis inversa con el objetivo de sintonizar los controladores PID. Tomó como modelo el obtenido por Alatiqui, y aplicó el método de desacoplamiento, para lograr que las variables manipuladas afecten sólo a una variable de salida. Una vez desacoplado el sistema aplicó el algoritmo de optimización por enjambre de partículas para sintonizar los parámetros 𝐾𝑝, 𝐾𝑖 y 𝐾𝑑 de cada controlador PID. Este algoritmo intenta mejorar la solución utilizando posibles candidatos, denominados “partículas”, basándose en una simple fórmula matemática con el cual se intenta conseguir el objetivo actualizando la velocidad de éstas y su posición en el espacio de búsqueda, de manera similar al desplazamiento de una bandada de aves o un grupo de peces. Por último, utilizó el método de Ziegler y Nichols para comparar sus resultados observándose en ellos una notable mejoría con respecto al tiempo de establecimiento y el sobreimpulso (Figura 2.7 y Figura 2.8). 1.4 1.2 1.0 0.8 0.6 0.4 0.2 0 0 1 2 3 4 5 6 7 8 9 10 Tiempo (seg) × 104 Figura 2.7. Comparación de la respuesta de flujo. 800 700 600 500 400 300 200 100 0 -100 0 1 2 3 4 5 6 7 8 9 10 Tiempo (seg) Figura 2.8. Comparación de la respuesta de conductividad. 22 Conductividad (µS/cm) Flujo (gpm) 2.2.3. Control Avanzado 2.2.3.1. Control Adaptativo Sobre este tipo de controlador, es posible encontrar el informe de Rao (Rao, et al., 1994), quien mencionó algunas técnicas de control inteligente para el mejoramiento de la automatización de procesos de desalinización, entre los cuales indica una combinación entre el control adaptivo y las redes neuronales utilizando el esquema mostrado en la Figura 2.9. Sin embargo no presentó ningún resultado con respecto al uso de esta técnica de control. Modelo de 𝑦𝑚 referencia 𝑟 𝑒 Suma 𝑢 𝑦 Red Planta 𝑝 neuronal No Lineal Muestreo Muestreo con retardo con retardo Figura 2.9. Control adaptivo directo de plantas no lineales usando redes neuronales. Por otro lado, se menciona el trabajo de Chaaben (Chaaben, et al., 2011), donde se explicó la importancia de la energía solar en zonas desiertas y áridas en las que la calidad de agua es mala, motivo por el cual se construyeron sistemas de desalinización por osmosis inversa aprovechando la energía fotovoltaica, sin embargo, este tipo de energía posee intermitencias elevadas que ocasionan la obturación de las membranas. Es por ello que se diseñó un Controlador Adaptivo por Modelo de Referencia de Estructura Variable (VSMRAC), cuyo algoritmo aseguró el valor mínimo del error de seguimiento entre la planta y el modelo de referencia, incluso ante la presencia de fuentes 23 intermitentes. El controlador diseñado no fue implementado, por lo que sólo se realizaron simulaciones los cuales se presentan en la Figura 2.10 y Figura 2.11. 2.8 Qs modelo Qs medido 2.6 2.4 2.2 2.0 1.8 1.6 1.4 0 2 4 6 8 10 12 14 16 18 20 Tiempo (seg) Figura 2.10. Comportamiento del flujo de agua producto. 0.40 0.35 Cs modelo Cs medido 0.30 0.25 0.20 0.15 0.10 0.05 0 2 4 6 8 10 12 14 16 18 20 Tiempo (seg) Figura 2.11. Comportamiento de la salinidad del agua producto. 2.2.3.2. Control Predictivo En el tema de control predictivo, Burden (Burden, et al., 2001) presentó una aplicación experimental de control avanzado y optimización en un módulo de ósmosis inversa de membrana de fibra hueca. El objetivo del estudio fue comparar el rendimiento de un controlador estándar proporcional-integral (PI) con el rendimiento de un controlador predictivo basado en restricciones 24 Salinidad del agua producto, Cs (g/L) Flujo de agua producto, Qs (lt/min) (CMPC). Este trabajo estuvo basado en la planta piloto utilizada por Assef (Assef, et al., 1997), sin embargo sustituyó las membranas arrolladas en espiral por un módulo de membrana de fibra hueca de 4". Se planteó, además, que la estrategia de control adecuada debía permitir la manipulación de la tasa de flujo de permeado manteniendo la calidad del producto con el fin de poder ajustar la producción de agua para satisfacer la demanda en una planta. El algoritmo del control predictivo empleado es una estrategia de control basada en el modelo que aplica control feedback y feedforward para procesos multivariables y utiliza una rutina de optimización que permite manejar restricciones; además, tanto si se trata del coste de las materias primas o las condiciones para una operación segura, se pueden considerar las limitaciones inherentes de todos los procesos. Por esta razón, este tipo de controladores ha crecido en popularidad en la industria de procesos, siendo aplicable a sistemas no cuadrados, es decir, para el caso en que un proceso tiene más salidas que entradas. Para el controlador PI, en el que se utilizó las constantes de ajuste de Ziegler y Nichols, la simulación produjo una respuesta aceptable a excepción de un desplazamiento de la conductividad por debajo del punto de consigna; sin embargo, de acuerdo a este estudio, esta estrategia de control no pudo controlar la planta piloto adecuadamente para un cambio de la referencia de la tasa de flujo de permeado debido a una no linealidad en el efecto del pH de alimentación como medio de control de la conductividad del permeado. A partir de estos resultados, se planteó el controlador CMPC cuyo propósito es operar la planta piloto mediante el mantenimiento de las cuatro variables de proceso dentro de rangos especificados. A diferencia de la estrategia de control PI, la estrategia CMPC pudo controlar todas las salidas del proceso y fue capaz de maximizar el flujo de permeado, mientras que mantuvo el resto de resultados de los procesos dentro de los límites, en particular el pH de alimentación. Robertson (Robertson, et al., 1996), utilizando también el modelo desarrollado por Alatiqi, propuso la implementación de controladores de matriz dinámica (DMC), los cuales fueron capaces de monitorear y controlar la 25 temperatura y pH de la alimentación, el flujo, la conductividad y la presión del permeado. El enfoque DMC permitió operar la planta para lograr varios flujos de permeado sin afectar la operación global ni la calidad del producto. De otra parte, Moncada (Moncada-Valerio, et al., 2012) enunció la necesidad de establecer un control óptimo para aumentar la productividad, reducir el consumo energético e incrementar la vida útil de la membrana de ósmosis inversa. Por ello, propuso un controlador con un modelo interno para predecir el comportamiento dinámico de la salida de la planta. Se realizó la simulación mediante Simulink para compararlo con un controlador PI tal como se muestra en la Figura 2.12 y la Figura 2.13. Figura 2.12. Respuestas temporales del sistema de control del bastidor de ósmosis inversa con controladores DMC vs PID. Figura 2.13. Señales de control del sistema de control del bastidor de ósmosis inversa con controladores DMC vs PID. 26 2.2.3.3. Control Robusto La mayoría de plantas desalinizadoras por ósmosis inversa se controlan mediante el uso de controladores clásicos; es decir, controladores PID o PI debido a la practicidad con la cual se manejan. Sin embargo, estas plantas cambian a menudo sus parámetros debido al ensuciamiento de la membrana semipermeable. Así, los parámetros iniciales que se tomaron en cuenta para el diseño del controlador cambian luego de horas de operación del sistema. Estas fluctuaciones se reflejan en los cambios de los puntos de operación y consecuentemente producen parámetros inciertos (Gambier & Badreddin, 2011). Para tratar con esta incertidumbre en el cambio de los parámetros se puede aplicar el enfoque del control robusto; sin embargo tanto a los fabricantes como a los operarios de las plantas no les agrada la idea de cambiar esta clásica estrategia de control debido a que se les resulta más complejo la comprensión de las técnicas avanzadas de control. Es por esta razón que existen pocos estudios realizados bajo este enfoque. Gambier (Gambier & Badreddin, 2011) desarrolla un controlador PID del bucle de control de caudal de permeado usando optimización paramétrica multi-objetivo, con lo cual logra menor sensibilidad a cambios de parámetros en la planta. Los resultados de simulación del modelo de una planta real, mostrados en la Figura 2.14, indican que el método propuesto produjo un desempeño satisfactorio para varias condiciones de operación. : Flujo de permeado en lt/h : Apertura de válvula Modelo Inferior Modelo Superior Modelo Aleatorio Figura 2.14. Flujo de permeado y señal de control para el controlador PID y varios modelos de la planta. 27 Por otro lado, Al-Haj (Al-haj Ali, et al., 2010) desarrolló un controlador MPC, el cual proporcionó buenas características de robustez, considerando una estrategia más realista y cercana a la práctica industrial. Los resultados se indican en la Figura 2.15 y la Figura 2.16. Tiempo de muestreo Tiempo de muestreo Tiempo de muestreo Tiempo de muestreo Figura 2.15. Resultados para un cambio en el set-point del 10% de flujo de permeado utilizando un controlador PI. Tiempo de muestreo Tiempo de muestreo Tiempo de muestreo Tiempo de muestreo Figura 2.16. Resultados para un cambio en el set-point del 10% de flujo de permeado utilizando un controlador MPC. 28 2.3. Objetivo de la Tesis Considerando que los procesos que presentan un mayor significado en el funcionamiento de las plantas desalinizadoras de agua de mar por ósmosis inversa son la variación del flujo volumétrico y la variación de la conductividad eléctrica, del agua producto, las cuales son afectadas por la influencia de cambios en las variables de presión y pH del agua de alimentación, esta tesis presenta el siguiente objetivo general: “Modelar una planta de desalinización por osmosis inversa y diseñar un sistema de control basado en control clásico multivariable”. Para cumplir con este objetivo, es necesario alcanzar los siguientes objetivos específicos: 1. Identificar mediante pruebas de simulación el modelo de una planta desalinizadora por ósmosis inversa. 2. Diseñar y simular el sistema de control mediante controladores clásicos multivariables. 3. Presentar una propuesta de implementación los controladores en la planta piloto de la PUCP. 4. Realizar tareas de análisis de resultados. 29 CAPÍTULO 3. IDENTIFICACIÓN Y DISEÑO DEL SISTEMA DE CONTROL DE UNA PLANTA PILOTO DE DESALINIZACIÓN POR ÓSMOSIS INVERSA. 3.1. Introducción Se implementó una planta piloto a escala realizando previamente una evaluación técnica y económica. Esto es, se realizó un diseño preliminar de una planta de desalinización valiéndonos de leyes de conservación para el posterior dimensionamiento de los componentes del sistema, considerando la ubicación, el presupuesto y la capacidad de producción requerida, la cual fue de 2 m³/día, suficiente para efectuar los ensayos de laboratorio necesarios en la actividad académica. Luego de ello, se solicitó los servicios de un proveedor para la ejecución del proyecto, lográndose definir la configuración del sistema propuesto. En la Figura 3.1 se muestra la configuración final de las líneas así como los equipos establecidos. Concluida la instalación del sistema de tuberías y equipos que constituyen la planta piloto se procedió a realizar la puesta en marcha. Esto nos permitió conseguir el objetivo del proceso de desalinización, ratificándose el mismo mediante un analizador portátil de conductividad eléctrica. Culminada esta etapa, se procedió a determinar la instrumentación necesaria para la implementación del sistema de control, la misma que se seleccionó de acuerdo a criterios técnicos y económicos. El número de sensores y actuadores estuvo acorde con la literatura, y en la cual se añadieron otros con el fin de proteger las membranas y evaluar las distintas estrategias de control. 30 31 Figura 3.1. Diagrama P&ID de la planta piloto de desalinización. Para la selección de los mismos se realizó un cuadro comparativo de propuestas de instrumentación el cual se indica en la Tabla 3.1. Tabla 3.1. Tabla comparativa de propuestas de instrumentación¹. Item Equipo Opción 1 Opción 2 Opción 3 Descripción Costo Descripción Costo Descripción Costo 1 Flujo DF05AR15 $2,428 8705T $3,670 volumétrico² Kobold Rosemount 2 Presión² CEN-3245 $645 PTP31 $422 2051TG $1,005 Kobold Endress Rosemount 3 pH² APM-1E-1 $3,377 CPF81D $2,156 3900-01-10 $1,864 Endress Rosemount 4 Conductividad - - CLS15D $3,287 228 $0³ eléctrica² Endress Rosemount 5 Temperatura² TDA-15H $836 - - 0068N11N00A $839 Kobold Rosemount 6 Válvula EV260B $1,290 - - Proporcional Danfoss 7 Variador de VLT 5004 $1,156 Sinamics V20 $1,081 PowerFlex 525 $485 Frecuencia Danfoss Siemens Allen Bradley 8 Fuente de Sitop Smart $621 1606XLE240E $215 24 VDC Siemens 9 Bomba de 522-A-N5-IEC $4,406 TPG800NHH $886 Dosificación Neptune Seko ¹ Información obtenida en febrero de 2015. ² Sensor – Transmisor. ³ Incluido con el sensor de pH. De acuerdo a esta información se procedió con la adquisición de los instrumentos, considerando los rangos de medición, condiciones de operación, precisión de las mediciones presentadas en sus respectivas fichas técnicas, y por último, se consideró el valor económico de cada instrumento. Una vez recibido los instrumentos, se definió la ubicación de los mismos en la planta piloto teniendo en cuenta las condiciones necesarias para una correcta medición en el caso de los sensores. Así también se definió la ubicación de los actuadores para una adecuada variación en las variables de entrada. 32 Posteriormente se coordinó con el proveedor de la planta piloto para la implementación de los instrumentos de acuerdo a lo establecido. Se incorporó además un controlador lógico programable CompactLogix y una pantalla táctil PanelView de Allen-Bradley los cuales fueron facilitados por la universidad. Para finalizar la implementación de la planta piloto, se procedió a instalar todos los elementos de seguridad e incorporar también dos contenedores metálicos para situar el controlador, el variador de velocidad, la fuente de alimentación de 24 VDC, la pantalla táctil, el cableado eléctrico y los pulsadores necesarios para la puesta en marcha del sistema. 3.2. Descripción de la Planta Piloto de Desalinización por Ósmosis Inversa del Laboratorio de Control y Automatización de la PUCP. 3.2.1. Implementación de la Planta Piloto La planta piloto de desalinización por ósmosis inversa situada en el Laboratorio de Control y Automatización de la PUCP, fue diseñada luego de un amplio estudio del estado del arte de este tipo de procesos, luego de lo cual se puso en funcionamiento utilizando componentes industriales, obteniéndose un módulo (planta piloto con instrumentación) que permita la realización de ensayos de laboratorio aplicados a la investigación. El módulo considera las etapas de pre- tratamiento y tratamiento de una planta desalinizadora de agua de mar, omitiéndose la etapa del pos-tratamiento, en la que se adiciona componentes minerales al agua producto para su distribución de acuerdo a normativas legales. Es decir, su construcción se basó en las etapas en las cuales se eliminan los componentes orgánicos y las sales presentes en el agua de alimentación, utilizando la tecnología de ósmosis inversa definida en el capítulo 1 del presente documento. Con fines prácticos, en el diseño de la planta se consideró su operación con agua salobre con una concentración de sales máxima aproximada de 18000 ppm, que difiere de la concentración del agua de mar. Sin embargo esta restricción no modifica significativamente la dinámica de este tipo de procesos, por lo que es válida la aplicación y comprobación en ésta, de la eficacia de las distintas técnicas 33 de la identificación de sistemas, de las diversas técnicas de control clásico y avanzado, así como su aplicación en el ámbito de la detección de fallos. En la Figura 3.1, se presenta el diagrama de tuberías e instrumentación (P&ID) de la planta piloto, donde se muestran las ubicaciones e interconexiones de los distintos componentes. Mediante este diagrama se describe el proceso, el cual inicia con el suministro de agua salobre desde el tanque de alimentación (TK-001) utilizando una bomba de baja presión (P-001). El tanque de alimentación (de 700 litros de capacidad) contiene agua salobre la cual previamente se preparó con 12 kg de sal y agua de la red para obtener una concentración de sales de 18000 ppm. La bomba de baja presión sirve de apoyo para que el flujo pueda vencer las caídas de presión en la línea 2 (ver Tabla 3.2 y Figura 3.1) originadas por los filtros hasta su llegada a la succión de la bomba de alta presión (P-002). Tabla 3.2. Líneas de la planta piloto. Línea Descripción Desde Hasta 1 3/4"-SS150-FDW-001 TK-001 P-001 2 3/4"-SS150-FDW-002 P-001 P-002 3 3/4"-SS300-FDW-003 P-002 U-001 4 3/4"-SS150-DW-004 U-001 TK-003 U-002 5 3/4"-SS300-BRR-005 U-001 TK-004 U-002 6 3/4"-SS300-FDW-006 3/4"-SS300-P-003 U-002 3/4"-SS300-BRR-005 7 1/4"-PEL-AS-007 TK-002 3/4"-SS150-P-002 8 3/4"-SS150-FDW-008 TK-003 TK-001 TK-004 9 3/4"-SS150-DR-009 TK-003 Exterior TK-004 F-002 F-003 10 3/4"-SS300-DW-010 US-02 3/4"-SS150-P-002 3/4"-SS300-P-006 11 3/4"-SS150-FDW-011 TK-001 Exterior 34 En este tramo se puede realizar la dosificación del ácido clorhídrico o sulfúrico (tanque TK-002) para la regulación del pH en el flujo de alimentación, mediante el uso de la bomba dosificadora US-02 a través de la línea 7. Una vez superada la etapa de pre-tratamiento, inicia la etapa de tratamiento con el bombeo a alta presión del agua de alimentación, utilizando la bomba P-002, por toda la línea 3 con dirección a las membranas de ósmosis inversa (U-001 y U-002). Aquí observamos la posibilidad de la operación del sistema en tres configuraciones distintas, gracias a la línea 6. Estas tres configuraciones son las siguientes: operación del sistema utilizando sólo la membrana U-001, operación del sistema utilizando las dos membranas en serie y por último, utilizando éstas en paralelo. El proceso objeto de estudio se sitúa en esta etapa, pues la separación de las sales del agua de alimentación es lograda gracias a la membrana de ósmosis inversa. El agua de rechazo o agua con elevada concentración de sales es conducida por la línea 5 y colectada en el tanque TK-004, mientras el agua producto (permeado) o agua con baja concentración de sales es dirigida por la línea 4 hacia el tanque TK-003. Las mediciones son tomadas de los transmisores de caudal (FIT-01) y conductividad (AE-02), las cuales registran las variables en el permeado para elaborar el algoritmo de control en el controlador multivariable UC-01. Este último envía las señales de control (4-20mA) a los actuadores, estos son la bomba dosificadora US-02 y el variador de velocidad US-01. Así mismo podemos tomar las medidas de pH, presión y temperatura del agua de alimentación antes del ingreso a la membrana utilizando los sensores AE-01, PIT-01 y TIT-03 respectivamente. El transmisor AIT-01 se encarga de recibir las señales de los sensores de pH y conductividad para enviarlo al controlador UC-01. Cabe mencionar, que es posible reemplazar el uso del variador de frecuencia conectado a la bomba de alta presión por la actuación sobre la válvula de control UV-01 para realizar la regulación de la presión en la alimentación. La diferencia radica en la reducción del consumo de energía eléctrica suministrado a la electrobomba P-002 al usar el variador de frecuencia. Terminada la experiencia en la planta piloto, la bomba P-003 posibilita la reutilización del agua a través de la recirculación de la mezcla de agua producto y 35 agua de rechazo hacia el tanque de alimentación. Para este caso, debe considerarse la adición de una solución básica como el hidróxido de sodio para restaurar el pH al valor inicial de la experiencia. Finalmente, es recomendable la limpieza mensual de las membranas habilitando la línea 10 y suministrando agua potable mediante la bomba P-003 y el correspondiente juego de válvulas para evitar el flujo a través de las líneas 2 y 8. 3.2.2. Características Generales de la Planta Piloto. Las características generales de la planta piloto son: Proceso : Desalinización por ósmosis inversa Volumen de alimentación : 0.7 m³ Caudal de alimentación : 0.45 m³/h Presión de alimentación : 190 m Concentración inicial : 18000 ppm Solución Ácida : Ácido clorhídrico Caudal de permeado : 2 m³/día Variables controladas : Caudal del permeado Conductividad del permeado Variables manipuladas : Presión de la alimentación pH de la alimentación Controlador : CompactLogix Visualización : Pantalla LCD 15” Uso : Investigación, docencia Figura 3.2. Planta piloto de desalinización por ósmosis inversa. 36 3.2.3. Características Técnicas de la Planta Piloto. Para el conocimiento de la calidad de los instrumentos utilizados en el sistema de control de la planta piloto se ofrece un extracto de las especificaciones técnicas de cada componente. La Figura 3.2 muestra la instalación completa en la cual se encuentran instalados los componentes en mención, mientras que la Tabla 3.3 y la Tabla 3.4 muestran los componentes principales del sistema de ósmosis inversa y del sistema de control respectivamente. Tabla 3.3. Componentes principales del sistema de ósmosis inversa. Ítem Descripción Código 1 Tanque de alimentación TK-001 2 Tanque de agua producto TK-003 3 Tanque de agua de rechazo TK-004 4 Tanque de ácido clorhídrico TK-003 5 Bomba de baja presión P-001 6 Bomba de alta presión P-002 7 Bomba de recirculación P-003 8 Membrana de ósmosis inversa U-001 U-002 9 Filtro ultravioleta F.001 10 Filtro multimedia F-002 11 Filtro carbón F-003 12 Filtro cartucho 5 μm F-004 Tabla 3.4. Componentes principales del sistema de control. Ítem Descripción Código 1 Un controlador lógico programable UC-01 2 Sensor c/ Transmisor de presión PIT-01 3 Sensor c/ Transmisor de pH AE-01 4 Sensor c/ Transmisor de caudal FIT-01 5 Sensor c/ Transmisor de conductividad AE-02 6 Sensor c/ Transmisor de temperatura TIT-03 7 Válvula de control UV-01 8 Variador de frecuencia US-01 9 Bomba dosificadora US-02 10 Pantalla LCD 15" 37 De acuerdo al listado de componentes, se describe a continuación cada uno de ellos. Los tanques TK-001, TK-003 y TK-004 almacenan tanto el agua de alimentación, con concentración de sales de 18000 ppm, como el agua producto y el agua de rechazo de manera independiente (ver Tabla 3.5 y Figura 3.3). El tanque TK-002, que contiene ácido clorhídrico (HCl), permite el cierre hermético de la solución considerando las precauciones necesarias para su manipulación. Tabla 3.5. Especificaciones técnicas de los tanques. Características Tanques TK-001 TK-002 TK-003 TK-004 Material¹ SS 316 PE SS 316 SS 316 Agua Agua Agua de Fluido HCl salobre producto rechazo Volumen 700 lt 80 lt 250 lt 250 lt ¹ SS 316 = acero inoxidable tipo 316, PE = polietileno. Figura 3.3. Tanques de alimentación (izquierda), tanque de agua producto y agua de rechazo (centro), y tanque de ácido (izquierda). La bomba de baja presión P-001 sirve de apoyo para vencer las pérdidas de presión causadas por los filtros, mientras que la bomba de alta presión P- 002 logra suministrar la presión requerida por las membranas de osmosis 38 inversa. Junto a la bomba de dosificación US-02 constituyen los elementos finales de control. La bomba de recirculación P-003 es un elemento secundario que favorece el reciclaje del agua utilizada (ver Figura 3.4). En la Tabla 3.6 se indican las especificaciones técnicas de estos equipos. Tabla 3.6. Especificaciones técnicas de las bombas. Características¹ Bombas P-001 P-002 P-003 US-02* Marca Epli Berkeley Salmson Seko Modelo MGPS10G Multi-H203 TGP800 Material² SS 316 SS 316 SS 316 PP Caudal máx. 6 m³/h 3.4 m³/h 6.5 m³/h 0.018 m³/h Presión máx.³ 61 m 210 m 51 m 124 m Potencia máx. 2.2 kW 1.64 kW 1.25 kW 0.0239 kW Velocidad 3420 RPM 3450 RPM 3500 RPM - Frecuencia 60 Hz 60 Hz 60 Hz - ¹ Información de acuerdo a catálogos correspondientes. ² SS 316 = acero inoxidable tipo 316. ³ Presión expresada en metros de columnas de agua. * En proceso de adquisición. Imagen de bomba manual actual. Figura 3.4. Bombas de baja presión (superior izquierda), alta presión (superior derecha), dosificación (inferior izquierda) y recirculación (inferior derecha). 39 Las membranas de ósmosis inversa U-001 y U-002 (Figura 3.5), con ayuda de la bomba de alta presión, facilitan la separación de gran porcentaje de sales contenidas en la alimentación. Sus especificaciones técnicas se indican en la Tabla 3.7. Tabla 3.7. Especificaciones técnicas de las membranas. Características¹ Membranas U-001 / U-002 Marca KOCH Modelo 4040-HR Presión de operación² 176 m Presión máxima² 422 m Flujo de permeado 8 m³/d Porcentaje de rechazo 99.55% pH permitido 4 - 11 Temperatura máxima de operación 45 °C Presión diferencial máxima² 42 m ¹ Información de acuerdo a catálogo correspondiente. ² Presión expresada en metros de columna de agua. Figura 3.5. Membranas de ósmosis inversa. El filtro ultravioleta F-001 (Figura 3.6) es una mejor alternativa al uso de cloro en la desinfección de virus y bacterias, pues su funcionamiento no implica la adición de sustancias químicas. Sus especificaciones técnicas se presentan en la Tabla 3.8. 40 Tabla 3.8. Especificaciones técnicas del filtro ultravioleta. Características¹ Filtro Ultravioleta F-001 Marca Sterilight Silver Modelo S12Q-PA Presión máxima² 88 m Temperatura máxima 40 °C Turbidez³ < 1 NTU Hierro < 0.3 mg/L Manganeso < 0.05 mg/L Dureza < 120 mg/L Taninos < 0.1 mg/L Material de carcasa SS 304 Potencia consumida 48 W ¹ Información de acuerdo a catálogo correspondiente. ² SS 304 = acero inoxidable tipo 304. ³ Unidades nefelométricas de turbidez. Figura 3.6. Filtro ultravioleta. El filtro multimedia F-002 (Figura 3.7) permite la remoción de los sólidos suspendidos mediante las capas internas las cuales frenan a las partículas de mayor tamaño en las capas superiores y a las de menor tamaño, en las inferiores. El filtro carbón F-003 elimina microorganismos, patógenos, amplia gama de químicos como combustibles, y ciertos metales como el plomo, cadmio o mercurio. Sus especificaciones técnicas se muestran en la Tabla 3.9. 41 Tabla 3.9. Especificaciones técnicas del filtro multimedia y filtro carbón. Características¹ Filtros Multimedia Filtro Carbón F-002 F-003 Marca General Electric General Electric Modelo M21 C21 Presión máxima 50 m 50 m Temperatura máxima 40 °C 40 °C Medios filtrantes Antracita Carbón Arena Activado Garnet Grava Material de carcasa Polietileno Polietileno Fibra de vidrio Fibra de vidrio ¹ Información de acuerdo a catálogo correspondiente. Figura 3.7. Filtro multimedia y filtro carbón El filtro cartucho de 5 μm F-004 (Figura 3.8) tiene el propósito de proteger las membranas de ósmosis inversa capturando partículas que hayan pasado por etapas previas del pretratamiento, evitando así daños o ensuciamientos prematuros en las mismas. Sus especificaciones técnicas se presentan en la Tabla 3.10. 42 Tabla 3.10. Especificaciones técnicas del filtro cartucho 5μm. Características¹ Filtro Cartucho 5 µm F-004 Marca KOKO Modelo KKFS10-21 Tipo Cartucho Material Polipropileno Presión máxima 85 m Temperatura máxima 38 °C ¹ Información de acuerdo a catálogo correspondiente. Figura 3.8. Filtro cartucho de 5μm. El controlador UC-01 (Figura 3.9) es el encargado de elaborar la ley de control para mantener la referencia de flujo volumétrico y conductividad eléctrica del agua producto. Sus especificaciones técnicas se presentan en la Tabla 3.11. Tabla 3.11. Especificaciones técnicas del controlador. Características¹ Controlador UC-01 Marca Allen-Bradley Modelo CompactLogix 5370 L3 Memoria de usuario 2 MB Tareas del controlador 32 Programas por tarea 100 Módulos expansores locales 16 Puntos de E/S de expansión locales 512 Tarjeta de memoria Flash SD 1 GB ¹ Información de acuerdo a catálogo correspondiente. 43 Figura 3.9. Controlador lógico programable. El transmisor de presión PIT-01 (Figura 3.10) entrega la información sobre la presión de la alimentación al controlador, sin embargo no se utiliza para la elaboración de la señal de control. Sirve para cuantificar el efecto de la variación del porcentaje de cierre de la válvula o la variación de la señal del variador de frecuencia. Sus especificaciones técnicas se presentan en la Tabla 3.12. Tabla 3.12. Especificaciones técnicas del transmisor de presión. Características¹ Transmisor de Presión PIT-01 Marca Rosemount Modelo 2051T Sensor Capacitivo Rango de medición 0 - 55.2 bar Precisión ± 0.075% Salida 4-20 mA Protocolo HART Material SS 316L Voltaje de alimentación 10.5 - 42.4 VDC ¹ Información de acuerdo a catálogo correspondiente. Figura 3.10. Transmisor de presión. 44 El sensor de pH AE-01 (Figura 3.11) entrega la información sobre la medición del pH de la alimentación y éste es enviado al controlador mediante el transmisor AT-01, sin embargo no se utiliza para la elaboración de la señal de control. Sirve para cuantificar el efecto de la variación del caudal volumétrico de la bomba dosificadora. Sus especificaciones técnicas se presentan en la Tabla 3.13. Tabla 3.13. Especificaciones técnicas del sensor de pH Características¹ Sensor de pH AE-01 Marca Rosemount Modelo 3900 Rango de medición 0-14 Precisión ± 0.01 Salida 4-20 mA Protocolo HART Material SS 316L Voltaje de alimentación 24 VDC ¹ Información de acuerdo a catálogo correspondiente. Figura 3.11. Transmisor de pH. El transmisor de flujo FIT-01 (Figura 3.12) se encarga de enviar la información del caudal volumétrico del permeado al controlador. Éste último elabora la señal de control utilizando además el valor de la conductividad. Sus especificaciones técnicas se presentan en la Tabla 3.14. 45 Tabla 3.14. Especificaciones técnicas del transmisor de flujo. Características¹ Transmisor de Flujo FIT-01 Marca Rosemount Modelo 8705 Sensor Magnético Rango de medición 0 - 8.5 m³/h Precisión ± 0.25% Presión máxima 19.6 bar Temperatura máxima 80 °C Salida 4-20 mA Protocolo HART Material SS 316L Voltaje de alimentación 24 VDC ¹ Información de acuerdo a catálogo correspondiente. Figura 3.12. Transmisor de flujo. El sensor de conductividad AE-02 (Figura 3.13) entrega la información de la conductividad eléctrica del permeado y éste es enviado al controlador mediante el transmisor AT-01 para la elaboración de la señal de control en conjunto con la señal de flujo volumétrico. Sus especificaciones técnicas se presentan en la Tabla 3.15. 46 Tabla 3.15. Especificaciones técnicas del sensor de conductividad. Características¹ Sensor de Conductividad AE-02 Marca Rosemount Modelo 228 Rango de medición hasta 2 S/cm Precisión ± 1% Salida 4-20 mA Protocolo HART Material SS 316L Voltaje de alimentación 24 VDC ¹ Información de acuerdo a catálogo correspondiente. Figura 3.13. Transmisor de conductividad. El transmisor de temperatura TIT-01 (Figura 3.14) proporciona la temperatura del agua de alimentación. Esta información es considerada en la operación de los equipos, pues de acuerdo a sus características estos presentan restricciones en la temperatura. Sus especificaciones técnicas se presentan en la Tabla 3.16. Tabla 3.16. Especificaciones técnicas del transmisor de temperatura. Características¹ Transmisor de Temperatura TIT-01 Marca Rosemount Modelo 68 Sensor PT 100 Rango de medición 0-100°C Precisión ± 0.03% Salida 4-20 mA Protocolo HART Material SS 316L Voltaje de alimentación 24 VDC ¹ Información de acuerdo a catálogo correspondiente. 47 . Figura 3.14. Transmisor de temperatura. El variador de frecuencia US-01 (Figura 3.15) es el primer actuador del sistema de control que se encarga de regular la velocidad de giro de la bomba P-002 con el fin de variar la presión suministrada hacia las membranas. Sus especificaciones técnicas se presentan en la Tabla 3.17. Figura 3.15. Variador de frecuencia. Tabla 3.17. Especificaciones técnicas del variador de frecuencia. Variador de Frecuencia Características¹ US-01 Marca Allen-Bradley Modelo PowerFlex 525 Potencia máxima 5 HP Factor de Potencia 0.98 Entrada -10 a +10 VDC 4-20 mA Precisión ±0.5% Temperatura ambiente máxima 50 °C ¹ Información de acuerdo a catálogo correspondiente. 48 La válvula de control UV-01 (Figura 3.16) se utiliza como una alternativa para el primer actuador. Se encarga de efectuar la acción de control de acuerdo al valor indicador por el controlador. Posee un actuador de tipo eléctrico con posicionador además de una entrada de 4 a 20 mA. Sus especificaciones técnicas se presentan en la Tabla 3.18. Tabla 3.18. Especificaciones técnicas de la válvula de control. Características¹ Válvula de Control UV-01 Marca Valbia Actuador eléctrico VB030-P Válvula 2500-3/4" Torque máximo 30 Nm Presión máxima 40 bar Entrada 4-20 mA 0-10 VDC Temperatura máxima 55°C Voltaje de alimentación 24 VDC ¹ Información de acuerdo a catálogo correspondiente. Figura 3.16. Válvula de control. El transmisor de pH y conductividad AIT-01 (Figura 3.17) se encarga de recibir las señales de los sensores AE-01 y AE-02, mostrar el valor de las mismas en un indicador y enviarlos hacia el controlador. Sus especificaciones técnicas se presentan en la Tabla 3.19. 49 Tabla 3.19. Especificaciones técnicas del transmisor de pH y conductividad. Transmisor de pH y Características¹ Conductividad AIT-01 Marca Rosemount Modelo 1056 Entradas 2 (4-20mA) Salidas 2 (4-20mA) Protocolo HART Voltaje de alimentación 24 VDC ¹ Información de acuerdo a catálogo correspondiente. Figura 3.17. Transmisor de pH y conductividad. 3.2.4. Adquisición de Datos de la Planta Piloto. Para la aplicación de las diversas técnicas de identificación de sistemas será necesario acondicionar la toma y envío de las señales utilizando una tarjeta de adquisición de datos DAQ – PCI 6229 National Instruments, proporcionada por el Laboratorio de Control y Automatización de la PUCP. Esto permitirá generar el tipo de señal que es necesario para lograr una correcta identificación de la dinámica del sistema. De otro lado, los algoritmos de control se podrán implementar en el controlador lógico programable UV-01 (Figura 3.9). 50 3.3. Identificación de la Planta Piloto. 3.3.1. Generalidades En la sección previa se detalló el funcionamiento de la planta piloto de desalinización por osmosis inversa comparando su estructura con una planta de desalinización de agua de mar mediante la misma tecnología. Por otro lado, en el capítulo anterior se subrayaron las variables más importantes consideradas en el sistema de control. Es así que en adelante se realizará la tarea de identificación del sistema mediante pruebas de simulación con el propósito de mostrar la metodología para la obtención de un modelo que represente el comportamiento dinámico del proceso de desalinización. Es posible aplicar el enfoque teórico o de primeros principios para desarrollar un modelo preciso y detallado basado en la teoría de solución-difusión y transferencia de materia (Jiang, et al., 2014), sin embargo este enfoque conduciría a modelos generalmente complejos y no lineales que deben ser sometidos a un proceso de simplificación y linealización además de presentar el inconveniente de requerir un conocimiento muy especializado sobre la tecnología del proceso (Behar & Martínez, 2010). Empleando técnicas de identificación de sistemas, se puede llevar a cabo un análisis experimental que permita el desarrollo de modelos matemáticos mediante la medición de las señales de entrada y salida del sistema, considerando que el modelo obtenido gobernará el comportamiento dinámico de estas señales y no describirá la estructura interna del sistema. Estos modelos son aproximaciones suficientes para muchas áreas de aplicación tales como el control automático (Isermann & Münchhof, 2010). Por ello, se emplearán las señales de entrada y salida para determinar el comportamiento temporal dentro de una clase de modelos matemáticos (Figura 3.18). El error (desviación) entre el proceso real o sistema y su modelo matemático deberá ser tan pequeño como sea posible (Isermann & Münchhof, 2010). 51 Señal de Señal de Proceso Salida Entrada u[k] (desconocido) y[k] Modelo Identificación Figura 3.18. Esquema de identificación. 3.3.2. Metodología de la Identificación El procedimiento de identificación posee una secuencia lógica resumida en la Figura 3.19 (Ljung, 1999), y cuyos pasos se describen como sigue: 1. Diseño del experimento, involucra la adquisición de datos, elección de variables a medir, tipo de señal de entrada, elección del periodo de muestreo y el tratamiento previo de la información. 2. Selección de la estructura del modelo, la cual es derivada a partir del conocimiento previo del proceso y de las perturbaciones. 3. Formulación del criterio, se establece para dar una medición de cuán bien se ajusta un modelo a los datos experimentales. 4. Cálculo del modelo, formulado como un problema de optimización, donde el mejor modelo es el que mejor se ajusta a los datos de acuerdo al criterio establecido. 5. Validación del modelo identificado, evaluar si éste puede ser considerado como una buena descripción del proceso. 52 Conocimiento a priori Diseño del Experimento Selección de la Datos Estructura del Modelo Formulación del Criterio Cálculo del Modelo Validación del Modelo No Revisar ¿Satisfactorio? Si Emplear Figura 3.19. Flujo del proceso de identificación. En la práctica, el procedimiento de la identificación de sistemas es iterativo y cuando se investiga un proceso del cual se tiene un conocimiento previo pobre, entonces es razonable empezar con un análisis de la respuesta transitoria o frecuencial para conseguir una estimación preliminar (burda) de la dinámica y de las perturbaciones (Åström & Wittenmark, 1997). Es muy probable que el primer modelo obtenido no supere el proceso de validación por lo que se debe revisar varios de los pasos del procedimiento de identificación. El modelo puede ser deficiente debido a varias razones tales como (Ljung, 1999):  El procedimiento numérico falló al encontrar el mejor modelo de acuerdo al criterio seleccionado.  El criterio no fue bien seleccionado.  El conjunto de modelos no fue el apropiado, pues no contienen una buena descripción del sistema.  El conjunto de datos no proporciona la suficiente información para la selección de buenos modelos. 53 Es importante señalar además, que la tarea de identificación debe asumir la condición de linealidad del proceso, para con ello evitar las dificultades asociadas a la teoría de control no lineal y porque la aproximación lineal resulta plausible en muchos casos (Behar & Martínez, 2010). Para la aplicación de lo mencionado anteriormente, se propuso realizar la tarea de la identificación del sistema multivariable como una serie de sistemas SISO, por lo que se debe encontrar un modelo matemático que represente la influencia de una variable de entrada, 𝑢𝑖, sobre una variable de salida, 𝑦𝑖. 3.3.3. Métodos de Identificación De acuerdo a la definición presentada, los diferentes métodos de identificación paramétrica pueden ser clasificados por los siguientes criterios (Isermann & Münchhof, 2010):  Clase de modelos matemáticos, los cuales pueden ser identificación no paramétrica, donde se utilizan modelos no paramétricos, que proporcionan una relación entre una cierta entrada y su correspondiente respuesta mediante una tabla o curva característica muestreada (tales como las respuestas al impulso, al escalón o respuestas en frecuencia presentadas en forma gráfica o tabulada); y por otro, la identificación paramétrica, donde se utilizan modelos paramétricos, esto es, ecuaciones que contienen explícitamente los parámetros del proceso.  Clase de señales de prueba empleadas, estas pueden ser señales determinísticas, aquí se encuentran los métodos tales como el análisis de Fourier y la estimación espectral de parámetros; estocásticas (aleatorias), se mencionan el análisis espectral; o seudoestocásticas (determinísticas con propiedades similares a las señales estocásticas), dentro de las que encontramos el análisis por correlación y la estimación paramétrica. 54  Cálculo del error entre el modelo y el proceso, estos pueden ser: error de entrada, error de salida y error de la ecuación generalizada. Típicamente se utiliza el error de salida, por ejemplo, cuando se usa la respuesta al impulso como modelo, y el error generalizado si se emplea las ecuaciones diferenciales, ecuaciones en diferencias o funciones de transferencia. Se incluyen además otra amplia distinción de los métodos los cuales son (Åström & Wittenmark, 1997):  Identificación fuera de línea, donde la información es almacenada para luego ser transferida a un ordenador donde se evaluarán y procesarán los datos.  Identificación en línea, en este caso la identificación se lleva a cabo paralelamente al experimento, esto es, el ordenador se encuentra acoplado con el proceso realizando la estimación con los datos que se encuentran disponibles en cada instante determinado. 3.3.4. Convergencia e Identificabilidad Las estimaciones que se determinaron en un horizonte de tiempo finito convergerán hacia los valores reales de las funciones de correlación y, por tanto: lim 𝐸{𝒆} = 0 (3.1) 𝑁→∞ Siempre que el modelo coincida en estructura y orden con el proceso (Isermann & Münchhof, 2010). El vector de parámetros 𝜽 de un modelo es identificable, si los valores estimados ?̂? convergen a los parámetros verdaderos 𝜃0 en el cuadrado medio (Isermann & Münchhof, 2010). Esto significa que: lim 𝐸{?̂?(𝑁)} = 𝜃0 (3.2) 𝑁→∞ 55 Para el criterio de mínimos cuadrados las condiciones necesarias se resumen en las siguientes: 1. El orden 𝑚 y el tiempo muerto 𝑑 son conocidos. 2. La señal de entrada 𝑢[𝑘] = 𝑈[𝑘] − 𝑈∞ es medida con exactitud y su valor 𝑈∞ es conocida con exactitud. 3. La señal de entrada 𝑢[𝑘] cumple con la condición de excitación persistente. 4. El proceso debe ser estable, controlable y observable. 5. La perturbación estocástica 𝑣[𝑘] la cual es superpuesta sobre la señal de salida 𝑦[𝑘] = 𝑌[𝑘] − 𝑌∞ debe ser estacionaria. El valor estacionario 𝑌∞ debe ser conocido con exactitud y debe corresponder a 𝑈∞. 6. El error 𝑒[𝑘] no debe estar correlacionado y 𝐸{𝑒[𝑘]} = 0. 3.3.5. Diseño del Experimento De acuerdo a la metodología propuesta, se procedió a realizar el estudio del experimento con el objeto de obtener un conocimiento a priori del proceso de desalinización en la planta piloto así como la instrumentación utilizada para el desarrollo de la investigación. La adquisición de los datos se puede realizar mediante la tarjeta DAQ – PCI 6229 National Instruments con su respectiva bornera para las conexiones de los sensores y actuadores. Para facilitar el procedimiento de identificación es posible utilizar el software Matlab con su módulo Simulink el cual permite la interconexión con la tarjeta DAQ. Los datos adquiridos por lo general no están listos para ser utilizados en el desarrollo del modelo, por lo que deben ser sometidos a un preprocesamiento utilizando un filtro pasa-bajos, antes de usar algún algoritmo de estimación. Las razones por las que se afecta la calidad de los datos, además del ruido, son la existencia de valores atípicos, que son datos los cuales no se ajustan en gran medida a otras partes de la información general, y por otro lado está la carencia de ciertos datos producto del mal funcionamiento del sensor o de las intermitencias del suministro eléctrico (Tangirala, 2014). 56 Para la selección del periodo de muestreo, se tomó en cuenta que valores mucho más pequeños que la constante de tiempo del sistema conducen a problemas numéricos debido a que los polos se aproximan al círculo unitario, además es posible que el modelo reproduzca con mayor exactitud el proceso en bandas de altas frecuencias. En caso contrario, un valor más grande que la constante de tiempo del sistema producirá datos con poca información acerca de la dinámica del mismo (Ljung, 1999). Por ello, se aplicó la recomendación de un valor similar al de la décima parte de la constante de tiempo principal del sistema (Kuo, 1997). 3.3.6. Identificación No Paramétrica Con el objeto de proveer una perspectiva acerca de las características (determinísticas) más importantes de la dinámica del proceso tales como el tiempo de retardo, ganancia y constante de tiempo (Tangirala, 2014), se realizaron las pruebas de simulación de la identificación no paramétrica. Determinación del Rango de Linealidad de la Planta Piloto Para cumplir con la condición de linealidad del proceso, es necesario determinar el rango donde el comportamiento presenta tal característica. Esta operación se realizará enviando señales al sistema en lazo abierto, para luego obtener una curva de valores de señal de entrada versus su respectivo valor de salida cuando ésta se encuentra en estado estacionario. Este procedimiento se resume en los siguientes pasos: 1. Se fija la entrada (variable manipulada 𝑢𝑖) a un valor constante. Se registra la salida del proceso cuando esta señal ha alcanzado el estado estable (𝑦𝑖). 2. Graficamos el punto (𝑢𝑖, 𝑦𝑖). 3. Repetimos el paso 1 registrando los valores de entrada en una tabla, considerando variaciones de 5%, donde el 100% corresponde a la apertura total del actuador. 4. En base a los puntos hallados se grafica la curva de 𝑢 versus 𝑦. 57 5. Se determina el rango donde el proceso presenta un comportamiento lineal. Un procedimiento similar fue aplicado por Alatiqui para determinar el rango de linealidad del proceso de desalinización, el cual se presenta en la Tabla 3.20. Por lo que, de esta misma manera debe aplicarse el procedimiento para la determinación del rango de linealidad del proceso en la planta piloto. Tabla 3.20. Rango de linealidad del proceso de desalinización. Variables Unidades Rango Lineal Flujo gpm 0.85-1.25 Presión psi 800-1000 Conductividad μS/cm 400-450 pH – 6-7 Determinación de la Respuesta al Escalón La respuesta al escalón es la señal más simple que puede utilizarse para la identificación y es aproximada por ejemplo mediante la apertura o cierre súbito de una válvula. Un escalón ideal es una señal cuyo tiempo de crecimiento es cero, lo que físicamente no puede ser realizada por la necesidad de una energía infinita; sin embargo si el tiempo de crecimiento inicial es varias veces más pequeño que la constante de tiempo más pequeña del sistema, el error que se introduce en la identificación puede considerarse despreciable (Behar & Martínez, 2010). Con esta premisa, se procedió a realizar una variación súbita de la presión dentro de todo el rango de linealidad y se capturó la respuesta del flujo volumétrico y de la conductividad eléctrica. Posteriormente se realizó una variación súbita del pH dentro de todo el rango de linealidad y se capturó la respuesta de la conductividad eléctrica. Las características dinámicas del proceso estudiado se presentan en la Figura 3.20 y la Figura 3.21. 58 Figura 3.20. Respuesta del flujo y de la conductividad al escalón de presión. Figura 3.21. Respuesta de la conductividad al escalón de pH. De las Figura 3.20 y la Figura 3.21, se determinan los parámetros mencionados. Estos son indicados en la Tabla 3.21. 59 Tabla 3.21. Características del proceso. Respuestas Características P – F P – C pH – F Retardo de tiempo 0 seg. 0 seg. 0 seg. Constante de tiempo principal 3 40 150 Tiempo de asentamiento 15 100 400 Ganancia 0.002 -0.51 -57 Orden del proceso 2 2 2 Con esta información, es posible afirmar que el comportamiento dinámico de cada proceso puede ser representado por un modelo de segundo orden (Smith & Corripio, 1991) con la estructura presentada en la ecuación (3.3). 𝐾 𝑔𝑖𝑗 = (3.3) (𝜏1𝑠 + 1)(𝜏2𝑠 + 1) Donde: 𝑔𝑖𝑗(𝑠) : Función de transferencia que relaciona la variación de la i-ésima variable entrada a la j-ésima variable de salida. 𝐾 : Ganancia estática del proceso. 𝜏1, 𝜏2 : Constantes de tiempo. Es decir, este modelo representará a cada proceso mostrado en la Figura 3.20 y en la Figura 3.21, así tenemos para 𝑖 = 1,2 y 𝑗 = 1,2 lo siguiente: 𝐹(𝑠) 𝑔11(𝑠) = (3.4) 𝑃(𝑠) 𝐶(𝑠) 𝑔21(𝑠) = (3.5) 𝑃(𝑠) 𝐶(𝑠) 𝑔22(𝑠) = (3.6) 𝑝𝐻(𝑠) 60 Sin embargo, como se mencionó anteriormente estos modelos sólo nos proporcionan una perspectiva de las características más importantes del proceso, por lo que será necesario determinar un modelo con un mayor grado de exactitud para obtener resultados satisfactorios en la implementación del controlador propuesto. 3.3.7. Identificación Paramétrica Fuera de Línea En contraste con lo anterior, la identificación paramétrica está basada en la minimización de ciertas señales de error por medio de métodos de regresión estadística y ha sido complementada con métodos especiales para sistemas dinámicos (Isermann & Münchhof, 2010). Aquí se estudian los métodos para estimar modelos paramétricos, que son el resultado de parametrizar las respuestas de los modelos de la planta y del ruido, en lugar de coeficientes de respuesta. Para ello se debe proporcionar la información como el retardo y el orden del sistema, lo que usualmente es obtenido de la identificación no paramétrica (Tangirala, 2014). En lo que a la condición experimental respecta, es importante que la señal de entrada sea persistentemente excitada. En términos generales, esto implica que todos los modos del sistema deben estar excitados durante el experimento de identificación (Söderström & Stoica, 1989). Para ello utilizaremos una señal binaria seudoaleatoria (PRBS), el cual es una señal periódica, determinística con propiedades similares al ruido blanco (Ljung, 1999). Determinación de la PRBS y la Correspondiente Respuesta del Sistema Para la generación de la PRBS se determinó su periodo básico aplicando las siguientes reglas (Chen & Gu, 2000): 2𝜋𝑇𝑠 𝑇𝑚𝑎𝑥 ≈ (3.7) 0.15 2𝜋𝑇𝑠 2𝜋𝑇𝑠 < 𝑇𝑚𝑖𝑛 < (3.8) 10 5 61 𝑇𝑃𝑅𝐵𝑆 ≈ 𝑇𝑚𝑖𝑛 + 𝑇𝑚𝑎𝑥 (3.9) donde, 𝑇𝑠 : Periodo de muestreo 𝑇𝑚𝑎𝑥: Ancho de pulso máximo. 𝑇𝑚𝑖𝑛: Ancho de pulso mínimo. La selección del tiempo de muestreo se llevó a cabo utilizando el criterio de la décima parte de la constante de tiempo principal del sistema. Con la información de la Tabla 3.21 se obtuvo el periodo básico de la PRBS correspondiente a (3.4), (3.5) y (3.6), este resultado se muestra en la Tabla 3.22. Tabla 3.22. Determinación del periodo básico de la PRBS. Planta 𝑻𝒔 𝑻𝒎𝒊𝒏 𝑻𝒎𝒂𝒙 𝑻𝑷𝑹𝑩𝑺 (seg) (seg) (seg) (seg) 𝑔11 0.3 1 13 14 𝑔21 4 6 168 174 𝑔22 15 19 629 648 Se aplicó, como regla general, una duración de la señal PRBS equivalente a 20 ciclos completos. La Figura 3.22, Figura 3.24 y Figura 3.23 muestran las señales PRBS anteriormente determinadas con las respectivas respuestas de las variables de salida. Figura 3.22. Respuesta del flujo volumétrico a la señal PRBS de presión. 62 Figura 3.23. Respuesta de la conductividad a la señal PRBS de presión. Figura 3.24. Respuesta de la conductividad a la señal PRBS de pH. Selección de la Estructura del Modelo Paramétrico La tarea más importante y más difícil del procedimiento de identificación es la selección de la estructura del modelo (Ljung, 1999). Esta estructura es derivada del conocimiento previo del proceso y de las perturbaciones. En algunos casos, el único conocimiento a priori es que el proceso puede ser descrito como un sistema lineal en un rango de operación en particular (Åström & Wittenmark, 1997). En la identificación de la planta piloto se consideró la estructura ARX (Auto-Regresiva Controlada) y la estructura ARMAX (Auto-Regresiva de Media Móvil Controlada), siendo esta última de mayor interés debido a que describe las propiedades de la perturbación añadiendo flexibilidad al modelo ARX mediante la descripción de la ecuación del error como un movimiento 63 de la media del ruido blanco (Ljung, 1999). La ecuación (3.10) y la Figura 3.25 representan las estructuras ARX (si 𝐶(𝑞) = 1) y ARMAX (si 𝐶(𝑞) ≠ 1), donde, a partir de la identificación no paramétrica se considera un proceso sin retardo (ver Tabla 3.21). 𝐴(𝑞)𝑦(𝑡) = 𝐵(𝑞)𝑢(𝑡) + 𝐶(𝑞)𝑒(𝑡) (3.10) donde: 𝐴(𝑞) = 1 + 𝑎 𝑞−1 + ⋯+ 𝑎 𝑞−𝑛𝑎 1 𝑛 𝑎 𝐵(𝑞) = 𝑏 𝑞−1 1 + ⋯+ 𝑏 𝑞−𝑛𝑏 𝑛 𝑏 𝐶(𝑞) = 𝑐 𝑞−1 1 + ⋯+ 𝑐 𝑞−𝑛𝑐 𝑛 𝑐 𝑪(𝒒) 𝒆 𝑨(𝒒) 𝑩(𝒒) 𝒖 𝒚 𝑨(𝒒) Figura 3.25. Diagrama de la estructura ARX (C=1) y ARMAX (C≠1). La ecuación (3.10) puede escribirse como un término que represente la dinámica del proceso y otro que represente la dinámica de la perturbación, como se indica en (3.11). 𝑦(𝑡) = 𝐺𝑝(𝑞, 𝜽)𝑢(𝑡) + 𝐺𝑣(𝑞, 𝜽)𝑒(𝑡) (3.11) Tal que: 𝐵(𝑞) 𝐶(𝑞) 𝐺𝑝(𝑞, 𝜽) = , 𝐺𝑝(𝑞, 𝜽) = (3.12) 𝐴(𝑞) 𝐴(𝑞) Los parámetros ajustables están dado por (3.13), considerando que en el caso de la estructura ARX se omitirán los valores de 𝑐1 hasta 𝑐𝑛 . 𝑐 𝜽 = [𝑎 … 𝑎 𝑏 …𝑏 𝑐 … 𝑐 ]𝑇1 𝑛𝑎 1 𝑛𝑏 1 𝑛𝑐 (3.13) El objetivo entonces es determinar los valores de 𝑛𝑎, 𝑛𝑏 y 𝑛𝑐 que logren ajustar el modelo al proceso estudiado utilizando las señales de entrada 𝑢(𝑘) y las señales de salida 𝑦(𝑘). 64 Formulación del Criterio Una vez que el modelo es seleccionado, lo que sigue es su estimación, esto esencialmente es un problema de optimización. Un criterio típico de estimación tiene una forma de minimización de una función de los errores de predicción (Tangirala, 2014). Este criterio, para sistemas discretos, es usualmente representado por (3.14) (Åström & Wittenmark, 1997). 𝑛 𝐽(𝜃) = ∑ 𝑓(𝜀(𝑡)) (3.14) 𝑖=0 Donde 𝐽(𝜃) : Función de costo de la estimación 𝑓(𝜀(𝑡)) : Función de ponderación del error 𝜀 : Error de entrada, error de salida o un error generalizado La función 𝑓 se basa generalmente en una medición de la distancia. En este trabajo se utiliza el método de mínimos cuadrados, los cuales minimizan la distancia euclideana entre los valores predichos y observados. Estimación de los Parámetros De acuerdo a lo anterior, se utilizó el método de mínimos cuadrados para la estimación de los parámetros del modelo. En el problema general de los mínimos cuadrados se supone que la variable calculada está dada por la expresión (3.15) (Åström & Wittenmark, 1997). ?̂? = 𝜃1𝜑1(𝑥) + 𝜃2𝜑2(𝑥) + … + 𝜃𝑛𝜑𝑛(𝑥) (3.15) Donde 𝜑1, 𝜑2, … , 𝜑𝑛 son funciones conocidas, 𝜃1, 𝜃2, … , 𝜃𝑛 son los parámetros desconocidos y los pares (𝑥𝑖, 𝑦𝑖) son obtenidos del experimento. El principio establece que los parámetros deben seleccionarse de tal manera que la función de costo (3.14), donde 𝑓 = 𝜀2 𝑖 y 𝜀𝑖 = 𝑦𝑖 − ?̂?𝑖, sea mínima. Para simplificar los cálculos, se introduce la siguiente notación: 65 𝝋 = [𝜑1 𝜑2 … 𝜑𝑛]𝑇 (3.16) 𝜽 = [𝜃1 𝜃2 … 𝜃 ]𝑇𝑛 (3.17) 𝒚 = [𝑦1 𝑦2 … 𝑦 𝑇 𝑛] (3.18) 𝜺 = [𝜀1 𝜀2 … 𝜀 𝑇 𝑛] (3.19) 𝚽 = [𝝋𝑻(𝑥 𝑻 1) 𝝋 (𝑥2) … 𝝋𝑻(𝑥 𝑇 𝑁)] (3.20) Formulando el problema en forma compacta, podemos escribir la función de costo con la expresión (3.21). 1 1 𝑱(𝜽) = 𝜺𝑇𝜺 = ‖𝜺‖2 (3.21) 2 2 Donde 𝜺 = 𝒚 − ?̂?. Por ello, la ecuación (3.15) se expresará como en (3.22). ?̂? = 𝚽𝜽 (3.22) Entonces, se debe determinar el vector de parámetros 𝜽 tal que ‖𝜺‖ sea mínima. La solución para este problema es única y se representa mediante la ecuación (3.23) si 𝚽𝐓𝚽 es no singular. ?̂? = (𝚽𝐓𝚽)−1𝚽𝐓𝒚 (3.23) Los algoritmos de identificación se encuentran implementados en el software Matlab, en los cuales se utiliza el criterio mencionado. La función 𝑎𝑟𝑥 implementa el método de estimación de mínimos cuadrados que utiliza la factorización QR de ecuaciones lineales sobredeterminadas, mientras que la función 𝑎𝑟𝑚𝑎𝑥 estima el modelo para una sistema de tres entradas y una salida utilizando el método de estimación iterativo ARMAX (Ljung, 2015). Para la aplicación de estos algoritmos, se realizó una lista de modelos candidatos con estructuras ARX y ARMAX y cuyos órdenes se muestran en la Tabla 3.23. En esta lista no se considera el tiempo de retardo debido a lo anteriormente expuesto. 66 Tabla 3.23. Órdenes de los polinomios A(q), B(q) y C(q) de los modelos candidatos. Estructura na nb nc ARX 2 2 0 ARX 3 3 0 ARMAX 1 1 1 ARMAX 2 2 2 ARMAX 3 3 3 ARMAX 4 4 4 Se utilizó estos modelos propuestos para la identificación de cada relación expresada en (3.4), (3.5) y (3.6), obteniéndose los parámetros de los polinómios 𝐴(𝑞), 𝐵(𝑞) y 𝐶(𝑞) para cada caso. Se consideró para ello aproximadamente un 70% de los datos presentados en la Figura 3.22, Figura 3.23 y Figura 3.24, con el objeto de utilizar la otra parte de los datos en el procedimiento de validación de los modelos. En la Tabla 3.24, Tabla 3.25 y Tabla 3.26 se muestran los resultados para la relación entre la presión de la alimentación y el flujo volumétrico del permeado establecida en (3.4). Tabla 3.24. Parámetros estimados del polinomio 𝑨(𝒒) para el modelo 𝒈𝟏𝟏. Estructura 𝑎1 𝑎2 𝑎3 𝑎4 ARX22 -0.4386 -0.4295 0 0 ARX33 -0.2683 -0.3005 -0.2619 0 ARMAX111 -0.9286 0 0 0 ARMAX222 -1.672 0.6973 0 0 ARMAX333 -1.6863 0.5967 0.1021 0 ARMAX444 -0.5949 -1.0607 0.4354 0.2499 Tabla 3.25. Parámetros estimados del polinomio 𝑩(𝒒) para el modelo 𝒈𝟏𝟏. Estructura 𝑏1 𝑏2 𝑏3 𝑏4 ARX22 -0.0000039 0.0002668 0 0 ARX33 -0.0000129 0.0001944 0.0001652 0 ARMAX111 0.0001676 0 0 0 ARMAX222 0.0001015 -0.0000504 0 0 ARMAX333 0.0000144 0.0001587 -0.0001482 0 ARMAX444 0.000009 0.0001871 0.0000366 -0.0001736 67 Tabla 3.26. Parámetros estimados del polinomio 𝑪(𝒒) para el modelo 𝒈𝟏𝟏. Estructura 𝑐1 𝑐2 𝑐3 𝑐4 ARMAX111 -0.7004 0 0 0 ARMAX222 -1.7328 0.7538 0 0 ARMAX333 -1.7452 0.6526 0.1016 0 ARMAX444 -0.6483 -1.0778 0.4866 0.2609 De igual manera, en la Tabla 3.27, Tabla 3.28 y Tabla 3.29 se muestran los resultados para la relación entre la presión de la alimentación y la conductividad eléctrica del permeado establecida en (3.5), y en la Tabla 3.30, Tabla 3.31 y Tabla 3.32, los resultados para la relación entre el pH de la alimentación y la conductividad eléctrica del permeado establecida en (3.6). Tabla 3.27. Parámetros estimados del polinomio 𝑨(𝒒) para el modelo 𝒈𝟐𝟏. Estructura 𝑎1 𝑎2 𝑎3 𝑎4 ARX22 -0.7329 -0.1583 0 0 ARX33 -0.5226 -0.3935 0.0537 0 ARMAX111 -0.9384 0 0 0 ARMAX222 -1.2946 0.3653 0 0 ARMAX333 -1.7568 0.7474 0.0278 0 ARMAX444 -0.9886 -0.649 0.6844 -0.0154 Tabla 3.28. Parámetros estimados del polinomio 𝑩(𝒒) para el modelo 𝒈𝟐𝟏. Estructura 𝑏1 𝑏2 𝑏3 𝑏4 ARX22 0.0030 -0.0720 0 0 ARX33 0.0040 -0.0594 -0.0269 0 ARMAX111 -0.0487 0 0 0 ARMAX222 -0.0070 -0.0337 0 0 ARMAX333 0.0007 -0.0582 0.0481 0 ARMAX444 0.0009 -0.0591 0.0077 0.0346 68 Tabla 3.29. Parámetros estimados del polinomio 𝑪(𝒒) para el modelo 𝒈𝟐𝟏. Estructura 𝑐1 𝑐2 𝑐3 𝑐4 ARMAX111 -0.1151 0 0 0 ARMAX222 -1.0356 0.3830 0 0 ARMAX333 -1.8169 0.8000 0.0322 0 ARMAX444 -1.0422 -0.6575 0.742 -0.0162 Tabla 3.30. Parámetros estimados del polinomio 𝑨(𝒒) para el modelo 𝒈𝟐𝟐. Estructura 𝑎1 𝑎2 𝑎3 𝑎4 ARX22 -0.5959 -0.2243 0 0 ARX33 -0.3245 -0.35 -0.0793 0 ARMAX111 -0.8888 0 0 0 ARMAX222 -1.0828 0.2065 0 0 ARMAX333 -0.2102 -0.7664 0.204 0 ARMAX444 0.3987 -0.4608 -0.7322 0.2137 Tabla 3.31. Parámetros estimados del polinomio 𝑩(𝒒) para el modelo 𝒈𝟐𝟐. Estructura 𝑏1 𝑏2 𝑏3 𝑏4 ARX22 0.1366 -10.531 0 0 ARX33 0.3075 -7.7292 -6.6973 0 ARMAX111 -6.7084 0 0 0 ARMAX222 0.1882 -7.2279 0 0 ARMAX333 0.3115 -7.5376 -5.7145 0 ARMAX444 0.3273 -7.3075 -10.3495 -6.543 Tabla 3.32. Parámetros estimados del polinomio 𝑪(𝒒) para el modelo 𝒈𝟐𝟐. Estructura 𝑐1 𝑐2 𝑐3 𝑐4 ARMAX111 -0.2778 0 0 0 ARMAX222 -1.1594 0.2457 0 0 ARMAX333 -0.2778 -0.8041 0.2408 0 ARMAX444 0.3466 -0.5264 -0.7753 0.2515 69 Validación y Selección del Modelo Una vez obtenido los modelos, es necesario revisarlos con el fin de revelar sus insuficiencias. Este procedimiento puede realizarse utilizando respuestas al escalón, al impulso, analizando los polos y zeros, y los errores de predicción (Åström & Wittenmark, 1997). En este caso se realizó el análisis utilizando el criterio del error final de predicción (Akaike, 1981) dado por (3.24), donde 𝑁 es el número de muestras y 𝑝 el número de parámetros del modelo. Así también, se usó el índice de performance (FIT) dado por (3.25), donde 𝑦 es la salida medida, ?̂? la salida estimada e ?̅? la media; además de una validación cruzada. Es importante recalcar que para la validación se utilizó la segunda parte de los datos, los mismos que no fueron utilizados en el procedimiento de estimación de los parámetros. 𝑁 + 𝑝 𝐹𝑃𝐸(𝑝) = 𝐽(𝜽) (3.24) 𝑁 − 𝑝 ‖𝑦 − ?̂?‖ 𝐹𝐼𝑇 = (1 − ) × 100% (3.25) ‖𝑦 − ?̅?‖ La Tabla 3.33 resume los resultados para la estimación del modelo que representa la expresión (3.4). De la misma manera, la Tabla 3.34 y la Tabla 3.35 muestran los resultados para la estimación del modelo que representa a las expresiones (3.5) y (3.6) respectivamente. En estas tablas se agrega también el valor de la función de costo. Tabla 3.33. Tabla de resultados para la estimación del modelo 𝒈𝟏𝟏. Ítem Estructura FIT(%) FPE 𝐽(𝜽) 1 ARX22 79.30 0.0007 0.0007 2 ARX33 80.41 0.0007 0.0007 3 ARMAX111 79.77 0.0007 0.0007 4 ARMAX222 82.70 0.0005 0.0005 5 ARMAX333 83.00 0.0005 0.0005 6 ARMAX444 83.04 0.0005 0.0005 70 Tabla 3.34. Tabla de resultados para la estimación del modelo 𝒈𝟐𝟏. Ítem Estructura FIT(%) FPE 𝐽(𝜽) 1 ARX22 90.11 9.8236 9.6775 2 ARX33 90.74 8.6806 8.4884 3 ARMAX111 86.24 19.0192 18.8055 4 ARMAX222 91.25 7.7700 7.5973 5 ARMAX333 92.95 5.1066 4.9380 6 ARMAX444 92.98 5.1323 4.9083 Tabla 3.35. Tabla de resultados para la estimación del modelo 𝒈𝟐𝟐. Ítem Estructura FIT(%) FPE 𝐽(𝜽) 1 ARX22 94.69 1.7311 1.7015 2 ARX33 95.21 1.4260 1.3898 3 ARMAX111 93.26 2.7915 2.7554 4 ARMAX222 95.92 1.0349 1.0085 5 ARMAX333 95.94 1.0392 0.9999 6 ARMAX444 96.00 1.0282 0.9768 Se representó mediante la Figura 3.26, Figura 3.27 y Figura 3.28, el valor de la función FPE pues es un buen criterio de selección debido a que además de considerar la función de costo también penaliza la complejidad del modelo dado por el número de parámetros 𝑝. Figura 3.26. Valor FPE en la estimación del modelo para 𝒈𝟏𝟏. 71 Figura 3.27. Valor FPE en la estimación del modelo para 𝒈𝟐𝟏. Figura 3.28. Valor FPE en la estimación del modelo para 𝒈𝟐𝟐. A partir de estos gráficos es posible apreciar que a partir de la estructura 4 (armax222) ya no existe variación significativa del valor de la función FPE. Por otro lado, la Tabla 3.33, Tabla 3.34 y Tabla 3.35 nos indican que los modelos armax presentan un mejor índice de performance (FIT), ello debido al considerar la estimación de la parte estocástica de manera separada a la dinámica del proceso. Cabe resaltar que la función de costo, en estas tablas, nos indican que no es necesario incrementar el orden del modelo, pues esto implicaría incrementar polos, los cuales no aportan dinámica significativa al modelo. En este sentido, se seleccionó el modelo con estructura armax222, como el modelo que representa, con mayor exactitud, el comportamiento dinámico del proceso de desalinización. En la Figura 3.29, Figura 3.30 y Figura 3.31 se muestra la validación cruzada del modelo seleccionado con los datos experimentales para estimar los parámetros que representen la relación (3.4), (3.5) y (3.6) respectivamente. 72 Figura 3.29. Validación cruzada del modelo ARMAX [2 2 2] con los datos adquiridos para la estimación del modelo 𝒈𝟏𝟏. Figura 3.30. Validación cruzada del modelo ARMAX [2 2 2] con los datos adquiridos para la estimación del modelo 𝒈𝟐𝟏. Figura 3.31. Validación cruzada del modelo ARMAX [2 2 2] con los datos adquiridos para la estimación del modelo 𝒈𝟐𝟐. Para finalizar el proceso de validación se muestra el gráfico de polos y ceros, así como las respuestas al escalón de los modelos seccionados en la Figura 3.32, Figura 3.33, Figura 3.34, Figura 3.35, Figura 3.36 y Figura 3.37. 73 Figura 3.32. Polos y ceros del modelo estimado para 𝒈𝟏𝟏. Figura 3.33. Polos y ceros del modelo estimado para 𝒈𝟐𝟏. Figura 3.34. Polos y ceros del modelo estimado para 𝒈𝟐𝟐. 74 Figura 3.35. Respuesta al escalón unitario del modelo estimado para 𝒈𝟏𝟏. Figura 3.36. Respuesta al escalón unitario del modelo estimado para 𝒈𝟐𝟏. Figura 3.37. Respuesta al escalón unitario del modelo estimado para 𝒈𝟐𝟐. 75 Los modelos encontrados quedan representados entonces por las expresiones (3.26), (3.27) y (3.28). 𝑏 𝑞−1 + 𝑏 𝑞−2 𝑐 𝑞−1 + 𝑐 𝑞−2 111 211 111 2 𝑓(𝑡) = 𝑝(𝑡) + 11 𝑒(𝑡) (3.26) 1 + 𝑎 𝑞−1 + 𝑎 𝑞−2 1 + 𝑎 𝑞−1 + 𝑎 𝑞−2 111 211 111 211 𝑏 𝑞−1 + 𝑏 𝑞−2 𝑐 𝑞−1 + 𝑐 𝑞−2 1 2 1 2 𝑐(𝑡) = 21 21 𝑝(𝑡) + 21 21 𝑒(𝑡) (3.27) 1 + 𝑎 𝑞−1 + 𝑎 𝑞−2 1 + 𝑎 𝑞−1 + 𝑎 𝑞−2 121 221 121 221 𝑏 𝑞−1 + 𝑏 𝑞−2 1 2 𝑐 𝑞−1 1 + 𝑐2 𝑞−2 𝑐(𝑡) = 22 22 𝑝ℎ(𝑡) + 22 22 𝑒(𝑡) (3.28) 1 + 𝑎 𝑞−1 + 𝑎 𝑞−2 1 2 1 + 𝑎 𝑞−1 1 + 𝑎 𝑞−2 22 22 22 222 Donde: 𝑎1 = −1.6720 , 𝑎 11 2 = 0.6973 11 𝑏1 = 0.0001015 , 𝑏2 = −0.0000504 11 11 𝑐1 = −1.7328 , 𝑐2 = 0.7538 11 11 𝑡𝑚𝑢𝑒𝑠𝑡𝑟𝑒𝑜 = 0.3 𝑠𝑒𝑔 11 𝑎1 = −1.2946 , 𝑎2 = 0.3653 21 21 𝑏1 = −0.0070 , 𝑏 21 2 = −0.0337 21 𝑐1 = −1.0356 , 𝑐2 = 0.3830 21 21 𝑡𝑚𝑢𝑒𝑠𝑡𝑟𝑒𝑜 = 4 𝑠𝑒𝑔 21 𝑎1 = −1.0828 , 𝑎2 = 0.2065 22 22 𝑏1 = 0.1882 , 𝑏2 = −7.2279 22 22 𝑐1 = −1.1594 , 𝑐 22 2 = 0.2457 22 𝑡𝑚𝑢𝑒𝑠𝑡𝑟𝑒𝑜 = 15 𝑠𝑒𝑔 22 76 3.4. Diseño del Sistema de Control de la Planta Piloto 3.4.1. Modelo Seleccionado para el Diseño del Controlador En esta sección se presentará el diseño de un controlador descentralizado que permita mantener en un punto de consigna el valor del flujo volumétrico y la conductividad eléctrica del permeado. Para ello se utilizará la técnica del desacoplamiento, que en términos generales tiene el objetivo de minimizar las interacciones que existen entre las variables en un proceso MIMO. Para el diseño del controlador multivariable se ha seleccionado el modelo propuesto por Alatiqi (Alatiqi, et al., 1989) dado en (2.2) a (2.6). Pues, tal como se mencionó anteriormente, en los demás modelos presentados en la Tabla 2.1 se tiene una ganancia positiva en el modelo que relacionan la presión de la alimentación y la conductividad eléctrica del permeado, así como en el modelo que relaciona el pH de la alimentación con la conductividad eléctrica del permeado, además este último modelo se exhibe también un comportamiento de fase no mínima (Gambier, et al., 2006). Arreglando las ecuaciones en forma matricial se obtiene (3.29). 𝒀 = 𝑮𝑼 (3.29) 𝐹 𝑔11 𝑔12 𝑃 [ ] = [𝑔 ] [ ] (3.30) 𝐶 21 𝑔22 𝑝𝐻 Donde: 0.002(0.056𝑠 + 1) 0 𝑔 𝑔 (0.003𝑠2 11 12 + 0.1𝑠 + 1) 𝐺 = [ 𝑔 ] = (3.31) 21 𝑔22 −0.51(0.35𝑠 + 1) −57(0.32𝑠 + 1) [(0.213𝑠2 + 0.7𝑠 + 1) (0.6𝑠2 + 1.8𝑠 + 1)] Cabe mencionar que, de acuerdo a la literatura, las unidades de tiempo de este modelo se encuentran en minutos, por lo que para fines prácticos se convirtió las unidades del modelo a segundos. Esta operación se realizó considerando que en la dinámica y el control del proceso la variable independiente es el tiempo y en consecuencia, la unidad de 𝑠 es el recíproco del tiempo (Smith & Corripio, 1991). 77 De esta manera se obtiene el modelo dado en (3.32) con su correspondiente diagrama de bloques mostrado en la Figura 3.38. 0.002(3.36𝑠 + 1) 0 (10.8𝑠2 + 6𝑠 + 1) 𝐺 = (3.32) −0.51(21𝑠 + 1) −57(19.2𝑠 + 1) [(766.8𝑠2 + 42𝑠 + 1) (2160𝑠2 + 108𝑠 + 1)] 𝒖𝟏 𝒚𝟏 𝒈𝟏𝟏 + 𝒈𝟏𝟐 𝒈𝟐𝟏 𝒖𝟐 𝒚𝟐 𝒈𝟐𝟐 + Figura 3.38. Modelo con estructura P-canónica. 3.4.2. Interacción entre Lazos En la práctica, un sistema de control multivariable puede tratarse lazo a lazo si la interacción entre estos es despreciable (Figura 3.39), sin embargo, existen situaciones cuando puede haber una interacción considerable entre los diferentes lazos de control, debido a esto puede ser difícil controlar tales sistemas lazo a lazo (Åström & Hägglund, 2009). 𝒈 … 𝒈 𝒓 + 𝒄𝟏 𝟎 𝒖 𝟏𝟏 𝟏𝒏 𝒚 [ ⋱ ] [ ⋮ ⋱ ⋮ ] 𝟎 𝒄 𝒈 𝒏 𝒏𝟏 ⋯ 𝒈𝒏𝒏 − Figura 3.39. Control descentralizado de un sistema multivariable. 78 Por lo tanto surge la necesidad de analizar tales interacciones evaluando algún método que garantice su minimización, para con ello garantizar un adecuado control de las variables del proceso de desalinización. Los dos métodos predominantes para la reducción de las interacciones entre variables son la separación de la dinámica y el desacoplamiento. La separación de la dinámica se logra sintonizando los controladores PID de modo que la constante de tiempo de un lazo cerrado difiera en un factor de diez veces la constante de tiempo del segundo lazo cerrado. El desacoplamiento se logra añadiendo una señal feedforward de la salida de un controlador hacia la salida del otro controlador (Vilanova & Visioli, 2012). Una manera sencilla de investigar el efecto de la interacción consiste en analizar cómo se ve influenciada la ganancia estática del proceso de un lazo por las ganancias en los otros lazos (Åström & Hägglund, 2009). Para ello, consideramos el sistema, dado por la Figura 3.38, con una cantidad de variables similar al proceso estudiado, y la expresamos mediante (3.33) y (3.34). 𝑌1(𝑠) = 𝑔11(𝑠)𝑈1(𝑠) + 𝑔12(𝑠)𝑈2(𝑠) (3.33) 𝑌2(𝑠) = 𝑔21(𝑠)𝑈1(𝑠) + 𝑔22(𝑠)𝑈2(𝑠) (3.34) Investigando cómo la ganancia estática en el primer lazo es influenciada por el controlador en el segundo lazo, Bristol asumió que el segundo lazo estaba en control perfecto, lo que quiere decir que la salida del segundo lazo es cero, es decir, ahora se tiene (3.35) y (3.36). 𝑌1(𝑠) = 𝑔11(𝑠)𝑈1(𝑠) + 𝑔12(𝑠)𝑈2(𝑠) (3.35) 0 = 𝑔21(𝑠)𝑈1(𝑠) + 𝑔22(𝑠)𝑈2(𝑠) (3.36) Eliminando la entrada 𝑈2(𝑠) de la ecuación (3.35) se obtiene (3.37). 𝑔11(𝑠)𝑔22(𝑠) − 𝑔12(𝑠)𝑔21(𝑠) 𝑌1(𝑠) = 𝑈 (𝑠) (3.37) 𝑔 (𝑠) 1 22 79 La razón de las ganancias estáticas del primer lazo cuando el segundo lazo está abierto y cuando el segundo lazo está cerrado está dado por (3.38). 𝑔11(0)𝑔22(0) 𝜆 = (3.38) 𝑔11(0)𝑔22(0) − 𝑔12(0)𝑔21(0) El parámetro 𝜆 se denomina índice de interacción de Bristol (o ganancia relativa) para sistemas TITO (Bristol, 1966) y se puede interpretar también como la interacción de señales a bajas frecuencias. Este parámetro puede ser generalizado para sistemas con muchas entradas y muchas salidas, la idea es comparar las ganancias estáticas de una salida cuando los otros lazos están abiertos con las ganancias cuando todas las otras salidas son cero (Åström & Hägglund, 2009). 3.4.3. Emparejamiento de Variables Para lograr el desacoplamiento del sistema multivariable, es preciso resolver el problema del emparejamiento de las variables de control con las variables manipuladas. Para ello definimos la ganancia relativa o índice de interacción de Bristol para el emparejamiento de variables 𝑦𝑖 − 𝑢𝑗 como la relación de dos ganancias representativas, siendo la primera la ganancia del proceso en un lazo aislado y, la segunda, la ganancia aparente del proceso en el mismo lazo cuando todos los lazos son cerrados (Qing-Guo, et al., 2008). Esto se expresa mediante (3.39). (𝜕𝑦𝑖⁄𝜕𝑢𝑗)𝑢𝑘≠𝑗𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡𝑒 𝜆𝑖𝑗 = = 𝑔 −1 𝑖𝑗[𝐺 (0)]𝑗𝑖 (3.39) (𝜕𝑦𝑖⁄𝜕𝑢𝑗)𝑦𝑙≠𝑖𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡𝑒 Así también definimos la matriz de ganancias relativas (RGA), Λ(𝐺) en (3.40). Λ(𝐺) = [𝜆𝑖𝑗] = 𝐺(0) ⨂ 𝐺−𝑇(0) (3.40) Donde ⨂ es el producto hardamart y 𝐺−𝑇(0) es la transpuesta de la inversa de 𝐺(0). Una vez calculada esta matriz se debe seguir este criterio (Tham, 1999): 80  𝜆𝑖,𝑗 ≈ 1: El par 𝑦𝑗 − 𝑢𝑖 es un buen candidato para un lazo independiente SISO.  𝜆𝑖,𝑗 < 0: La función de transferencia aparente tiene un signo diferente al bucle abierto uno, por lo que existe la posibilidad de un error de cálculo de la señal de retroalimentación, lo que lleva a la inestabilidad.  𝜆𝑖,𝑗 positivo (pequeño): La "ganancia" aparente del bucle cerrado es mucho más alta que la del bucle abierto. Esto puede causar la degradación del rendimiento o incluso la inestabilidad cuando se cierra el bucle. Se identifica una entrada que tiene poco efecto en una salida determinada en lazo abierto pero con efecto significativo en lazo cerrado debido al acoplamiento y la retroalimentación. Probablemente no es un buen candidato emparejamiento.  𝜆𝑖,𝑗 positivo (grande): La ganancia de bucle abierto es mayor que el del bucle cerrado, por lo que el control puede llegar a ser ineficaz cuando se cierra el resto de los bucles. Para el caso del proceso de desalinización dado por las expresiones (2.2) a (2.6) la matriz RGA resulta ser una matriz identidad de 2 × 2, por lo que aplicando el criterio se concluye que el emparejamiento de variables crearía el primer lazo de control dado por la variable controlada, flujo volumétrico de permeado, y la variable manipulada, presión de alimentación; asimismo tendríamos el segundo lazo de control con la variable controlada, conductividad eléctrica del permeado, y la variable manipulada, pH de la alimentación (ver Figura 3.40). Este es el mejor escenario, si desea controlar el sistema utilizando un controlador PID descentralizado. 3.4.4. Desacoplamiento del Sistema En los casos cuando el control multilazo (Figura 3.40) no alcanza efectivamente las especificaciones deseadas, una posible estrategia para abortar el control MIMO podría ser transformar la función de transferencia matricial en una matriz diagonal, esta estrategia se conoce como desacoplamiento (Albertos & Sala, 2004). La Figura 3.41 muestra la nueva estructura del sistema utilizando el desacoplador. 81 𝑭 − 𝑺𝑷 + 𝑷 𝑭 𝑪𝟏 𝑷𝒓𝒐𝒄𝒆𝒔𝒐 𝑪𝑺𝑷 + 𝒑𝑯 𝑪 𝑪𝟐 − Figura 3.40. Emparejamiento de variables. El desacoplamiento es una manera sencilla de tratar con la dificultad creada por las interacciones en los lazos. La idea consiste en diseñar un controlador que reduzca el efecto de la interacción. Idealmente, los cambios en un punto de consigna deberían afectar sólo a la salida del proceso correspondiente (Åström & Hägglund, 2009). 𝑷∗ 𝑷 𝑭 𝑫𝒆𝒔𝒂𝒄𝒐𝒑𝒍𝒂𝒅𝒐𝒓 𝑷𝒓𝒐𝒄𝒆𝒔𝒐 𝑫 𝑮 𝒑𝑯∗ 𝒑𝑯 𝑪 Figura 3.41. Sistema desacoplado. Existen dos maneras de lograr el desacoplamiento (Albertos & Sala, 2004): mediante la cancelación feedforward de términos de acoplamiento cruzados (enfoque en función de transferencia), y en base a la medición de estados vía una ley de retroalimentación (enfoque en espacio de estados). En este trabajo se realizó el desacoplamiento utilizando el enfoque en función de transferencia calculando un desacoplador simplificado (Luyben, 1970) calculado mediante la expresión (3.41). La estructura del sistema se muestra en la Figura 3.42. −𝑔12 −𝑔21 𝐷12 = ; 𝐷21 = (3.41) 𝑔11 𝑔22 82 𝑫𝒆𝒔𝒂𝒄𝒐𝒑𝒍𝒂𝒅𝒐𝒓 𝑷𝒓𝒐𝒄𝒆𝒔𝒐 ∗ 𝑷 𝑷 𝒈 𝑭 + 𝟏𝟏 + 𝑫𝟏𝟐 𝒈𝟏𝟐 𝑫𝟐𝟏 𝒈𝟐𝟏 ∗ 𝒑𝑯 𝒑𝑯 + 𝒈𝟐𝟐 + 𝑪 Figura 3.42. Sistema con desacoplador simplificado. Utilizando el modelo del proceso de desalinización dado por (3.36), se procedió a calcular el desacoplador mediante la expresión (3.41), donde debido a que el término 𝑔12 es cero sólo se obtuvo la expresión para 𝐷21. Este único desacoplador se expresa mediante (3.42). 0.51 (2160𝑠2 + 108𝑠 + 1) (21𝑠 + 1) 𝐷21 = − (3.42) 57 (766.8𝑠2 + 42𝑠 + 1) (19.2𝑠 + 1) Sin embargo, es importante que el desacoplador sea físicamente realizable (Corriou, 2010), esto implica que la red de desacoplo debe poseer elementos que sean propios, causales y estables. De acuerdo a (3.42) esta condición no se satisface debido a que la función de transferencia no es propia, por lo que fue necesario forzar la realizabilidad. Esto se logró insertando un polo simple estable con la multiplicidad adecuada (Garrido, et al., 2012), como se indica en (3.43). 1 (3.43) (𝜆𝑠 + 1)𝑟 Para no alterar la dinámica del desacoplador se fijó por conveniencia un polo de - 0.5, es decir un valor de 𝜆 = 2. En la Figura 3.43, se aprecia la lejanía del polo adicionado (en rojo) con respecto a los polos originales de la dinámica del desacoplador. 83 Figura 3.43. Adicción de polo para forzar la realizabilidad del desacoplador. De esta manera el desacoplador del sistema estará determinado por la expresión (3.44). 0.51 (2160𝑠2 + 108𝑠 + 1) (21𝑠 + 1) 1 𝐷21 = − (3.44) 57 (766.8𝑠2 + 42𝑠 + 1) (19.2𝑠 + 1) (2𝑠 + 1) Luego de resolver el problema de la realizabilidad del controlador, se procedió a simular el sistema desacoplado en Simulink tal como se indica en la Figura 3.44. Figura 3.44. Diagrama del sistema con desacoplo en Simulink. En la Figura 3.45 se muestra la respuesta de las variables de flujo volumétrico y conductividad eléctrica de permeado para un escalón en la presión en la alimentación de 800 a 840 psi manteniendo el valor del pH de la alimentación en un valor de 7. Las líneas a tramos indican un desacoplamiento ideal, mientras que la 84 línea continua describe el comportamiento dinámico del sistema utilizando el descoplador (3.44). Figura 3.45. Respuesta del sistema para un escalón de presión de 40 psi. Se pudo lograr entonces la reducción de la interacción que existe entre la presión de la alimentación y la conductividad eléctrica del permeado. No obstante, es necesario resaltar que esta técnica presenta algunas dificultades como (Skogestad & Postlethwaite, 2005): 1. Ser muy sensible a los errores de modelado e incertidumbres. 2. Puede no ser deseable su uso frente al rechazo de perturbaciones. 3. Si la planta posee polos y ceros en el semiplano derecho, entonces el desacoplamiento generalmente introducirá polos y ceros extra en el sistema en lazo cerrado. La mayoría de estos problemas se evitan si se utiliza un desacoplamiento parcial cuando la planta es triangular superior o inferior (Skogestad & Postlethwaite, 2005), como es el caso del proceso de desalinización. 85 3.4.5. Control Multivariable por Desacoplo utilizando PID Se ha demostró que es posible la reducción de las interacciones en el proceso de desalinización, por lo que luego se procedió a diseñar el controlador PID para cada lazo del sistema multivariable. Para realizar ello, se utilizó el método de Skogestad (Skogestad, 2003) o método SIMC. El método de sintonización propuesto utiliza modelos con estructuras específicas, y además presenta una serie de fórmulas para la aproximación de los modelos cuyas estructuras no se encuentren cubiertas por este método. Utilizando estas aproximaciones se logró reducir los modelos con el objeto de aplicar el mismo. Las ecuaciones (3.45) y (3.46) muestran los modelos simplificados correspondientes a las expresiones 𝑔11 y 𝑔22 de la ecuación (3.32). En la Figura 3.46 se verifica la aproximación de los modelos estudiados. 𝐾 0.002 𝑔∗ 11 = = (3.45) 𝑇𝑠 + 1 3.2𝑠 + 1 𝐾 −57 𝑔∗ 22 = = (3.46) (𝑇1𝑠 + 1)(𝑇2𝑠 + 1) (81.5𝑠 + 1)(7.3𝑠 + 1) Figura 3.46. Aproximación de modelos. 86 Determinados los modelos simplificados, se procede a utilizar el método SIMC, el cual para procesos sin retardo de tiempo se resume en la Tabla 3.36 (Haugen, 2004). Tabla 3.36. Fórmulas de Skogestad para la sintonización del PID. Proceso 𝑲𝒄 𝑻𝒊 𝑻𝒅 𝐾 1 𝑘 𝑇 0 𝑠 𝐾𝑇 1 𝐶 𝐶 𝐾 𝑇 min{𝑇, 𝑘1𝑇𝐾𝑇 𝐶} 0 𝑇𝑠 + 1 𝐶 𝐾 1 𝑘1𝑇𝐶 𝑇 (𝑇𝑠 + 1)𝑠 𝐾𝑇𝐶 𝐾 𝑇1 min{𝑇 , 𝑘 𝑇 (𝑇 𝑠 + 1)(𝑇 𝑠 + 1) 𝐾𝑇 1 1 𝐶} 𝑇2 1 2 𝐶 𝐾 1 4𝑇 4𝑇 𝑠2 4𝐾𝑇2 𝐶 𝐶 𝐶 El autor explica que el valor de 𝑘1 generalmente es 1, sin embargo puede incrementarse para conseguir un mejor atenuación de las perturbaciones. Además indica que el valor de 𝑇𝐶 debe mayor que cero, por lo que estos parámetros han sido variados para la sintonización de cada lazo. Los valores utilizados para cada caso se resumen en la Tabla 3.37. Tabla 3.37. Valores de sintonización utilizados para el proceso estudiado. Lazo Parámetros k1 𝑻𝑪 𝑲𝒄 𝑻𝒊 𝑻𝒅 PID 1 1 10 160 3.2 0 PID 2 3 30 -0.05 81.5 0 La estructura del controlador PID usado en el lazo de control, de acuerdo al método propuesto, es la estructura en serie o cascada la cual se indica en (3.47). 𝑇 𝑠+1 𝐶(𝑠) = 𝐾𝑐 ( 𝑖 ) (𝑇 𝑠 + 1) (3.47) 𝑇𝑖𝑠 𝑑 Estos controladores fueron implementados en el software Simulink junto al desacoplador dado por (3.44). La estructura del sistema de control propuesto para se muestra en la Figura 3.47. 87 Figura 3.47. Estructura del sistema de control propuesto. 3.4.6. Análisis de Resultados Utilizando la estructura establecida en la Figura 3.47, se simuló el sistema de control, considerando un cambio en la referencia del flujo volumétrico de permeado de 1 gpm a 1.2 gpm a los 100 segundos de iniciado el control. De manera simultánea se procedió a realizar un cambio en la referencia de la conductividad eléctrica del permeado, de un valor de 450 μS/cm a 410 μS/cm. Los resultados se muestran en la Figura 3.48 y Figura 3.49. Figura 3.48. Respuesta del flujo del permeado. 88 Figura 3.49. Respuesta de la conductividad eléctrica del permeado. De la Figura 3.48 y Figura 3.49, observamos que el sistema de control, compuesto por los controladores PID y el desacoplador establecido, logra alcanzar de manera eficaz los puntos de consigna de las variables de control. Se aprecia además en líneas rojas el comportamiento dinámico del sistema de control sin el uso del desacoplador. Se ha mostrado entonces, que las interacciones son mayores cuando se omite este componente. Por otro lado, se realizó una variación en la referencia de 1.2 a 1.12 gpm, simulando un cambio en la demanda del flujo de agua desalinizada. Posteriormente se aplicó una perturbación en el flujo volumétrico del permeado, simulando una fuga de agua producto de 0.05 gpm, verificándose, en la Figura 3.50 y Figura 3.51, la rápida respuesta del sistema de control ante esta perturbación de tipo escalón, manteniéndose el objetivo de control en la referencia deseada. 89 Perturbación Figura 3.50. Respuesta del flujo del permeado. Perturbación Figura 3.51. Respuesta de la conductividad eléctrica del permeado. 90 Los resultados fueron cuantificados utilizando el índice IAE (integral del error absoluto), el cual brinda información sobre la calidad y efectividad del controlador. Este índice es calculado mediante la expresión (3.48). 𝐼𝐴𝐸 = ∫|𝑒(𝑡)|𝑑𝑡 (3.48) Donde, para el caso de simulación por computador, se presenta la versión discreta del mismo. Este viene representado por la expresión (3.49). 𝑛 𝐼𝐴𝐸 = ∑|𝑒(𝑘)| (3.49) 𝑘=1 Así, se resumen los resultados en la Tabla 3.38: Tabla 3.38. Resultados del sistema de control. Flujo Conductividad PI y PI PI PI desacoplado desacoplado 𝑒𝑠𝑠 0 0 0 𝑚 [%] 0 34 0 𝑡𝑠𝑠 50 300 200 IAE sin perturbación 20 1.75 × 104 1.17 × 104 IAE con perturbación 33 2.79 × 104 1.34 × 104 Debido a que el flujo de permeado sólo depende de la presión aplicada en la alimentación, la respuesta no varía al aplicar el control con desacoplamiento. La Tabla 3.38 nos muestra que la conductividad eléctrica del permeado posee un menor tiempo de establecimiento (𝑡𝑠𝑠) y un sobreimpulso (𝑚) nulo al aplicar la estrategia de control con desacoplamiento. Para todos los casos el error en estado estacionario (𝑒𝑠𝑠) fue nulo. Así también, el menor valor del índice IAE nos indica que el sistema de control se desempeña de manera eficaz en este último caso. Se concluye entonces que el sistema de control con desacoplamiento, logra el objetivo establecido dentro del rango de linealidad del proceso multivariable. 91 CAPÍTULO 4. PROPUESTA DE IMPLEMENTACIÓN PRÁCTICA DEL SISTEMA DE CONTROL. 4.1. Introducción. En el presente capítulo se presenta la propuesta de implementación del sistema de control de la planta piloto de desalinización. Por consiguiente, se describe el algoritmo utilizado para lograr el objetivo de control del flujo volumétrico y conductividad eléctrica del agua producto, así también se describe el hardware necesario para la elaboración de la señal de control. Por último se desarrolla una interfaz gráfica que posibilitará el monitoreo y control de la planta piloto, la misma que será instalada en una pantalla táctil y operará a nivel cliente. 4.2. Algoritmo de Control. De acuerdo con el capítulo anterior, se presentó un controlador clásico por desacoplo como solución al problema del control de la planta piloto de desalinización por ósmosis inversa, por lo que en la Figura 4.1 se ha esquematizado el diagrama de flujo del sistema de control. Este sistema contempla los modos de operación automático y manual, sin embargo, de lo mencionado anteriormente, no se elaboró una configuración cliente-servidor debido a que el sistema es una planta piloto instalada en el laboratorio de la PUCP con fines académicos. Para el caso del control manual, la producción de agua desalinizada y la calidad de la misma dependerá únicamente de los valores ingresados por el operador para la actuación en la presión y el pH de la alimentación mediante la variación de la 92 frecuencia en bomba de alta presión y la regulación del flujo de ácido en la bomba dosificadora (P-002 y US-02 en la Figura 3.1). INICIO Modo Leer modo No ¿Automático? Si Almacenar las variables flujo y Presión y pH conductividad del de la permeado, presión y alimentación pH de la alimentación P, pH (F, C, P, pH) Set-point de flujo y Calcular los errores Enviar las señales conductividad e = F – F de control a cada 1 sp del permeado e2 = Csp - C actuador Fsp y Csp Elaborar las señales de Obtener el valor de Parámetros control: presión y las variables F y C del pH de la mediante de los controlador alimentación transmisores de (P, pH) flujo y conductividad del permeado Aplicar las señales de control mediante cada actuador Obtener el valor de las variables F y C mediante de los transmisores de flujo y conductividad del permeado FIN Figura 4.1. Diagrama de flujo del sistema de control. Para el caso del modo automático, se inicia con el almacenamiento de las variables manipuladas y controladas, luego, mediante el ingreso del punto de consigna en las variables controladas, Fsp y Csp, se calculan los errores e1 y e2. Estos errores ingresan en la subrutina constituida por el controlador, y cuyo diagrama de flujo se detalla en 93 la Figura 4.2. Esta subrutina inicia con la obtención de los errores mencionados para calcular tres señales: la primera es calculada por el controlador 𝑃𝐼1, dado por (3.47), y donde se utilizan los parámetros 𝐾𝑝, 𝐾𝑖 y 𝐾𝑑 de la Tabla 3.37, la segunda señal es elaborada por el desacoplador presentado en (3.42); y la tercera señal es elaborada por el controlador 𝑃𝐼2 el cual fue presentado por (3.47) y cuyos parámetros se establecieron en la Tabla 3.37. Las dos últimas señales se adicionan para elaborar la señal de control pH, mientras que la primera señal, P, se enviará directamente al actuador correspondiente. Todas las señales se implementarán en su versión discreta con un periodo de muestreo de 1 segundo. Esta subrutina culmina con el envío de las señales de control hacia la interfaz gráfica así como a los actuadores mencionados. INICIO Errores Obtener los e1 y e2 errores e1 e2 Mediante PI1(z) Mediante D21(z) Mediante PI2(z) Parámetros calcular la calcular la calcular la Parámetros Kp1, Ki1, acción de acción de acción de Kp2, Ki2, Kd1 control control control Kd2 P pH1 pH2 Sumar la acción de control y la señal del desacoplador. Obtener la señal Obtener la señal Presión de la pH de la de control de control alimentación alimentación P pH Enviar la señal de control P y pH FIN Figura 4.2. Diagrama de flujo del controlador (subrutina). Seguidamente, se aplican las señales de control, P y pH, mediante el variador de frecuencia y la bomba dosificadora, respectivamente. Por último, se almacenan todas las variables para la siguiente iteración del bucle de control. 94 4.3. Hardware del Sistema de Control. El algoritmo de control presentando en la sección anterior será implementado en un controlador CompactLogix 5370 L3, el cual requerirá de parte del hardware de la planta piloto mencionado en el capítulo 3 del presente documento. Estos equipos son presentados en el diagrama de la Figura 4.3, el cual resume el flujo de información a través de cada uno de ellos. El controlador recabará la información de la lectura de los sensores de flujo volumétrico y conductividad eléctrica del permeado, Rosemount 8705 y Rosemount 228. Para esta última variable será necesario enviar la señal al controlador mediante el transmisor Rosemount 1056, el cual se indica en la Figura 3.17. Posteriormente, el controlador elaborará las señales de control de acuerdo a lo establecido en el apartado anterior. Estas señales de control serán enviadas hacia los actuadores los cuales son el variador de frecuencia y la bomba de dosificación, PowerFlex 525 y Seko TGP, respectivamente. Los resultados de respuesta transitoria de las variables de control se mostrarán en el PanelView, el cual contendrá una interfaz gráfica con la capacidad de comunicación con el controlador para el envío y recepción de las señales. Bomba dosificadora Seko TGP (alimentación) Sensor de flujo Rosemount 8705 (permeado) Controlador CompactLogix 5370 L3 Interfaz gráfica en pantalla táctil PanelView 1500 Sensor de conductividad Rosemount 228 Variador de (permeado) frecuencia PowerFlex 525 (alimentación) Figura 4.3. Hardware del sistema de control. 95 4.4. Software del Sistema de Control Se realizaron las pruebas de implementación del algoritmo de control utilizando el software RSLogix 5000, cuya interfaz se aprecia en la Figura 4.4. En él se desarrolló una tarea periódica denominada Controlador, donde se ejecutó la rutina ControlPI_Desacoplo en un intervalo de tiempo de 1 segundo. Dicha rutina contiene el código correspondiente al algoritmo presentado en la Figura 4.1. Figura 4.4. Interfaz del software RSLogix 5000. Para verificar el correcto funcionamiento del controlador diseñado en el capítulo anterior, se incorporó otra rutina que describió el comportamiento dinámico de la planta desalinizadora, dada por (3.32), con un periodo de 0.2 segundos el cual fue menor al tiempo de actualización de la acción de control. Ambas rutinas fueron descargadas a un PLC virtual denominado RSLogix Emulate 5000 que permitió simular el comportamiento del sistema de control. Por otro lado, se desarrolló una interfaz gráfica para la interacción con el sistema de control, la misma que será instalada en el PanelView 1500. Para ello, se utilizó el software FactoryTalk View (Figura 4.5) para el diseño de las ventanas que facilitan la visualización de las variables. 96 Figura 4.5. Interfaz del software FactoryTalk View Studio ME. Se creó una ventana en la que se realizó el diagrama funcional de la planta piloto de desalinización por ósmosis inversa, mostrado en la Figura 4.6. Asimismo, se incluyó el selector correspondiente al modo de operación y un botón de apertura de la ventana del controlador (UC-01). Figura 4.6. Diagrama funcional de la planta piloto de desalinización. 97 De la misma manera se diseñó la ventana del controlador, al que se agregó gráficos de tendencia de las variables manipuladas y variables controladas. La Figura 4.7 muestra esta ventana, donde además se incluye las casillas correspondientes al ingreso de los parámetros del controlador. Toda la interfaz fue enlazada con las variables del PLC virtual para verificar la interacción entre las mismas. Figura 4.7. Simulación de la planta desalinizadora de agua de mar por ósmosis inversa. Como se aprecia en este gráfico, se ingresó un cambio en la referencia de 1 a 1.05 gpm, alcanzándose la referencia en un tiempo aproximado de 50 segundos. Luego se modificó la referencia de la conductividad eléctrica de 450 a 430 μS/cm, esto se logró en un tiempo de 2 minutos. Posteriormente se añadió un cambio en la referencia de flujo de 1.05 a 1 gpm, además de una perturbación en el flujo simulando un fuga de agua producto de 0.05 gpm, con resultados satisfactorios en el control de las variables de flujo y conductividad eléctrica del permeado. La presente propuesta se implementará en la planta piloto reemplazando el código de la rutina PlantaDesalinizadora por un algoritmo de lectura de las variables de control enviadas por los transmisores correspondientes. Por último, se incorporará una tarea adicional de escritura de las variables manipuladas para el correcto envío de los datos hacia los actuadores. 98 CONCLUSIONES  Se implementó una planta piloto de desalinización por osmosis inversa en el laboratorio de la PUCP realizando previamente una evaluación técnica y económica. Se implementó además la instrumentación de la misma para la aplicación de la identificación de sistemas y la estrategia de control.  Se estudió el comportamiento dinámico del proceso de desalinización mediante ósmosis inversa en una planta piloto, determinándose que la variación del flujo volumétrico del agua producto depende únicamente de la presión de la alimentación, mientras que la variación de la conductividad eléctrica del permeado depende tanto de esta última variable como también del pH de la alimentación.  Se modeló una planta de desalinización por ósmosis inversa mediante el uso de las técnicas de identificación de sistemas. Los modelos obtenidos presentan una estructura ARMAX de segundo orden sin retardo de tiempo, los cuales describen el comportamiento dinámico del proceso objeto de estudio con una aproximación promedio del 90%.  Se diseñó un sistema de control clásico multivariable basado en desacoplamiento para controlar de forma efectiva una planta de desalinización mediante ósmosis inversa. Los resultados muestran una reducción del tiempo de establecimiento en un 33% con respecto del control PI descentralizado en el control de la variable de conductividad eléctrica. Además, se verifica una disminución del índice IAE en un 52% para este mismo caso, lo cual indica una mayor atenuación de las perturbaciones utilizando la técnica de control PI con desacoplamiento.  Se desarrolló una propuesta de implementación práctica del controlador diseñado, incluyendo además, una interfaz gráfica para la interacción con el sistema de control de la planta piloto. 99 RECOMENDACIONES  Aplicar el sistema de control diseñado a la planta piloto de desalinización mediante el reemplazo de la rutina correspondiente al modelo del comportamiento dinámico de la planta.  Estudiar la variante del algoritmo de identificación recursiva para la aplicación en un sistema de control que contemple todo el rango de linealidad del proceso de desalinización por ósmosis inversa.  Comparar el desempeño del sistema de control diseñado en la planta piloto frente a las distintas técnicas de control mencionadas en la literatura, con el objeto de determinar las ventajas y desventajas de los mismos. 100 BIBLIOGRAFÍA Abbas, A., 2006. Model predictive control of a reverse osmosis desalination unit. Desalination, 194(1-3), p. 268–280. Abbas, A. & Al-Bastaki , N., 2005. Modeling of an RO water desalination unit using neural networks. Chemical Engineering Journal, November, 114(1-3), p. 139–143. Akaike, H., 1981. Modern development of statistical methods. Trends and Progress in System Identification, Volumen 7, pp. 169-184. Alatiqi, I. M., Ettouney, H. M. & El-Dessouky, H. T., 1999. Process control in water desalination industry: an overview. Desalination, 126(1-3), p. 15–32. Alatiqi, I. M., Ghabris, A. H. & Ebrahim, S., 1989. System identification and control of reverse osmosis desalination. System identification and control of reverse osmosis desalination, Volumen 75, pp. 119-140. Albertos, P. & Sala, A., 2004. Multivariable Control Systems: An Engineering Approach. 2004 ed. Heidelberg(Berlin): Springer. Al-haj Ali, M., Ajbar, A., Alhumaizi, K. & Ali, E., 2008. Controlling Water Quality using Reverse Osmosis: The Development of Simplified Dynamic Model. UKACC Control. Al-haj Ali, M., Ajbar, A., Ali, E. & Alhumaizi, K., 2010. Robust model-based control of a tubular reverse-osmosis desalination unit. Desalination, May, 255(1-3), p. 129–136. Assef, J. Z., Watters, J. C., Deshpande, P. B. & Alatiqi, I. M., 1997. Advanced control of a reverse osmosis desalination unit. Journal of Process Control, 7(4), p. 283–289. Åström, K. J. & Hägglund, T., 2009. Control PID avanzado. Madrid(Madrid): Pearson Educación. Åström, K. J. & Wittenmark, B., 1997. Computer-Controlled Systems: Theory and Design. 3rd ed. Upper Saddle River(New Jersey): Prentice Hall. 101 Bartman, A. R., Christofides, P. D. & Cohen, Y., 2009. Nonlinear Model-Based Control of an Experimental Reverse-Osmosis Water Desalination System. Industrial & Engineering Chemistry Research, April, 48(13), p. 6126–6136. Bartman, A. R., McFall, C. W., Christofides, P. D. & Cohen, Y., 2009. Model-predictive control of feed flow reversal in a reverse osmosis desalination process. Journal of Process Control, 19(3), p. 433–442. Behar, A. & Martínez, M., 2010. Identificación y Control Adaptativo. 1st ed. Madrid: Pearson. Bristol, E. H., 1966. On a new measure of interaction for multivariable process control. IEEE Transactions on Automatic Control, pp. 133-134. Burden, A. C., Deshpande, P. B. & Watters, J. C., 2001. Advanced process control of a B-9 Permasep® permeator desalination pilot plant. Desalination, 133(3), p. 271–283. Carrasco, N. y otros, 2014. Estado del arte sobre el modelado y control de plantas desalinizadoras por ósmosis inversa, Lima: s.n. Chaaben, A. B., Andoulsi, R., Sellami, A. & Mhiri, R., 2011. MIMO Modeling Approach for a Small Photovoltaic Reverse Osmosis Desalination System. Journal of Applied Fluid Mechanics, Volumen 4, pp. 35-41. Chen, J. & Gu, G., 2000. Control Oriented System Identification: An H∞ Approach. 1st ed. s.l.:John Wiley & Sons. Cipollina, A., Micale, G. & Rizzuti, L., 2009. Seawater Desalination: Conventional and Renewable Energy Processes (Green Energy and Technology). 2009 ed. Heidelberger(Berlin): Springer. Corriou, J.-P., 2010. Process Control: Theory and Applications. 1st ed. Heidelberg(Berlin): Springer. Cotruvo, J. y otros, 2010. Desalination Technology: Health and Environmental Impacts. 1st ed. Boca Raton(Florida): CRC Press. Eykhoff, P., 1974. System Identification Parameter and State Estimation. 1st ed. Chichester: Wiley. Gambier, A. & Badreddin, E., 2011. A robust control approach for a reverse osmosis desalination plant. Proceedings of the 12th International Conference on Environmental Science and Technology, pp. 8-10. 102 Gambier, A., Krasnik, A. & Badreddin, E., 2007. Dynamic Modeling of a Simple Reverse Osmosis Desalination Plant for Advanced Control Purposes. American Control Conference, 2007. ACC '07, pp. 4854 - 4859. Gambier, A., Wellenreuther, A. & Badreddin, E., 2006. Optimal control of a reverse osmosis desalination plant using multi-objective optimization. Computer Aided Control System Design, 2006 IEEE International Conference on Control Applications, 2006 IEEE International Symposium on Intelligent Control, 2006 IEEE, October.pp. 1368 - 1373. García-Rodríguez, L. & Gómez-Camacho, C., 2001. Perspectives of solar-assisted seawater distillation. Desalination, May, 136(1-3), p. 213–218. Garrido, J., Vázquez, F. & Morilla, F., 2012. Diseño de Sistemas de Control Multivariable por Desacoplo con Controladores PID. X Simposio CEA de Ingenería de Control, Marzo.pp. 54-71. Haugen, F., 2004. PID Control. Trondheim(Sør-Trøndelag): Tapir Academics Press. Isermann, R. & Münchhof, M., 2010. Identification of Dynamic Systems. 2011 ed. Heidelberg(Berlin): Springer. Jiang, A. et al., 2014. Mathematical Modeling and Simulation of SWRO Process Based on Simultaneous Method. Journal of Applied Mathematics, Volume 2014, pp. 1-11. Keesman, K., 2011. System Identification. 2011 ed. Heidelberger(Berlin): Springer. Kim, J. S. y otros, 2008. Auto tuning PID controller based on improved genetic algorithm for reverse osmosis plant. Proceedings of World Academy of Science Engineering and Technology, Volumen 30, pp. 1001-1007. King, M., 2011. Process Control: A Practical Approach. 1st ed. Chichester: Wiley. Kumar, B., Dewal, M. L. & Mukherjee, S., 2013. Control and monitoring of MSF-RO hybrid desalination process. Control, Automation, Robotics and Embedded Systems (CARE), 2013 International Conference on, December.pp. 1 - 5. Kuo, B. C., 1997. Sistemas de Control Digital. 2da ed. Ciudad de México(Distrito Federal): CECSA. Ljung, 1999. System Identification: Theory for the User. 2nd ed. Upper Saddle River(New Jersey): Prentice Hall. Ljung, L., 2015. System Identification Toolbox™ User's Guide, Natick: The MathWorks, Inc.. 103 Ljung, L. & Söderström, T., 1983. Theory and Practice of Recursive Identification. Cambridge(Massachusetts): The MIT Press. Luyben, W. L., 1970. Distillation decoupling. AIChE Journal, March, Volume 16(2), p. 198– 203. McFall, C. W., Bartman, A., Christofides, P. D. & Cohen, Y., 2008. Control and Monitoring of a High Recovery Reverse Osmosis Desalination Process. Industrial & Engineering Chemistry Research, 26 July, 47(17), p. 6698–6710. Mhaskar, P. y otros, 2008. Isolation and handling of actuator faults in nonlinear systems. Automatica, 44(1), p. 53–62. Mindler, A. B. & Epstein, A. C., 1986. Chapter 2.6 Measurements and control in reverse osmosis desalination. Desalination, August, Volumen 59, pp. 343-379. Moncada-Valerio, J., Rivas-Perez, R. & Sotomayor-Moriano, J., 2012. Control predictivo multivariable de un bastidor de ósmosis inversa de una planta desalinizadora de agua de mar. Memoria XV Congreso Latinoamericano de Control Automático, Octubre. Park, J.-M.y otros, 2009. PID Tuning Based on RCGA Using Ziegler-Nichols Method. Transactions of the Korean Society for Noise and Vibration Engineering, 19(5), pp. 475-481. Pintelon, R. & Schoukens, J., 2001. System Identification: A Frequency Domain Approach. New York(New York): IEEE Press. Qing-Guo, W., 2002. Decoupling Control. 2003 ed. Heidelberg(Berlin): Springer. Qing-Guo, W., Zhen, Y., Wen-Jian, C. & Chang-Chieh, H., 2008. PID Control for Multivariable Processes (Lecture Notes in Control and Information Sciences). 2008 ed. Heidelberg(Berlin): Springer. Rahardianto, A., Shih, W.-Y., Lee, R.-W. & Cohen, Y., 2006. Diagnostic characterization of gypsum scale formation and control in RO membrane desalination of brackish water. Journal of Membrane Science, August, 279(1-2), p. 655–668. Ramillo, L., Gómez de Soler, S. & Coppari, N., 2003. Tecnologías de Proceso para Desalinización de Aguas. Comisión Nacional de Energía Atómica (CNEA, Buenos Aires), pp. 22-27. Rao, G. P. y otros, 1994. Towards improved automation for desalination processes, Part II: Intelligent control. Desalination, 97(1-3), pp. 507-528. 104 Rathore, N. S., Kundariya, N. & Narain, A., 2013. PID Controller Tuning in Reverse Osmosis System based on Particle Swarm Optimization. International Journal of Scientific and Research Publications, June, 3(6), pp. 1-5. Riverol, C. & Pilipovik, V., 2005. Mathematical Modeling of Perfect Decoupled Control System and Its Application: A Reverse Osmosis Desalination Industrial-Scale Unit. Journal of Automated Methods and Management in Chemistry, Volumen 2005(2), p. 50–54. Robertson, M. W. y otros, 1996. Model based control for reverse osmosis desalination processes. Desalination, April, 104(2), p. 59–68. Skogestad, S., 2003. Simple analytic rules for model reduction and PID controller tuning. Journal of Process Control, June, 13(4), p. 291–309. Skogestad, S. & Postlethwaite, I., 2005. Multivariable Feedback Control: Analysis and Design. 2nd ed. Hoboken(New Jersey): Wiley-Interscience. Smith, C. A. & Corripio, A. B., 1991. Control Automático de Procesos. 1ra ed. Ciudad de México(Distrito Federal): Limusa. Söderström, T. & Stoica, P., 1989. System Identification. Cambridge: Prentice Hall. Syafiie, S., Tadeo, F., Palacin, L. & De Prada, C., 2008. Modelling for dynamic simulation of pretreatment in Reverse Osmosis plants. Industrial Engineering and Engineering Management, 2008. IEEM 2008. IEEE International Conference on, December.pp. 1663 - 1667. Tangirala, A. K., 2014. Principles of System Identification: Theory and Practice. Boca Raton(Florida): CRC Press. Tham, M. T., 1999. Multivariable control: an introduction to decoupling control. Department of Chemical and Process Engineering, University of Newcastle upon Tyne, pp. 1-19. Torky, O. M., Elamvazuthi, I. & Hanif, N. H., 2009. PC based SCADA system for reverse osmosis desalination plants. Research and Development (SCOReD), 2009 IEEE Student Conference on, November.pp. 442 - 445. Vilanova, R. & Visioli, A., 2012. PID Control in the Third Millennium: Lessons Learned and New Approaches. 2012 ed. Heidelberg(Berlin): Springer. Voutchkov, N., 2013. Desalination Engineering: Planning and Design. 1st ed. New York(New York): McGraw-Hill Education. 105 Wade, N. M., 2001. Distillation plant development and cost update. Desalination, May, 136(1-3), p. 3–12. Zhu, Y., 2001. Multivariable System Identification. Eindhoven: Elsevier Science. Zilouchian, A. & Jafar, M., 2001. Automation and process control of reverse osmosis plants using soft computing methodologies. Desalination, 135(1-3), p. 51–59. 106 ANEXOS Anexo 1. Programa para el pre-procesamiento de los datos. close all; clc; % Preprocesamiento de los datos de respuesta al % escalón de presión. y1 = F1(1:50); y1p = mean(y1); y2 = C1(1:50); y2p = mean(y2); p1 = P1 - min(P1); f1 = F1 - y1p; c1 = C1 - y2p; p1 = dtrend(p1); f1 = dtrend(f1); figure(1); subplot(3,1,1); plot(time,f1,'b'); grid on; subplot(3,1,2); plot(time,c1,'b'); grid on; subplot(3,1,3); stairs(time,p1,'r'); grid on; % Preprocesamiento de los datos de respuesta al % escalón de pH. y3 = C2(1:50); y3p = mean(y3); ph2 = pH2 - min(pH2); c2 = C2 - y3p; figure(2); subplot(2,1,1); plot(time,c2,'b'); grid on; subplot(2,1,2); plot(time,ph2,'b'); grid on; Anexo 2. Programa para la generación de la señal binaria pseudoaleatoria. %% PBRS Presión - Flujo clear all; clc; load prbsdata1; load prbsdata2; load prbsdata3; Test = 15; tau = 3; Ts = tau/10; Tmax = (2*pi*Ts)/0.15; Tmin = (2*pi*Ts)/5; Tmax = ceil(Tmax); Tmin = ceil(Tmin); Tprbs = Tmax + Tmin; t1 = 0:Ts:20*Tprbs; t1=t1'; u1 = idinput(length(t1),'PRBS',[0 1/Tprbs],[0 200]); prbs = [t1 u1]; figure(1); clf; stairs(t1,u1,'Linewidth',1.5); grid on; save prbsdata1 Ts t1 u1; %% PBRS Presión - Conductividad clear all; clc; load prbsdata1; load prbsdata2; load prbsdata3; Test = 100; tau = 40; Ts = tau/10; Tmax = (2*pi*Ts)/0.15; Tmin = (2*pi*Ts)/5; Tmax = ceil(Tmax); Tmin = ceil(Tmin); Tprbs = Tmax + Tmin; t2 = 0:Ts:20*Tprbs; t2=t2'; u2 = idinput(length(t2),'PRBS',[0 20/Tprbs],[0 200]); prbs = [t2 u2]; figure(1); clf; stairs(t2,u2,'Linewidth',1.5); grid on; save prbsdata2 Ts t2 u2; %% PBRS pH - Conductividad clear all; clc; load prbsdata1; load prbsdata2; load prbsdata3; Test = 400; tau = 150; Ts = tau/10; Tmax = (2*pi*Ts)/0.15; Tmin = (2*pi*Ts)/5; Tmax = ceil(Tmax); Tmin = ceil(Tmin); Tprbs = Tmax + Tmin; t3 = 0:Ts:20*Tprbs; t3=t3'; u3 = idinput(length(t3),'PRBS',[0 20/Tprbs],[0 1]); t3p = t3(135:length(t3)) - t3(135); u3p = u3(135:length(t3)); prbs = [t3p u3p]; figure(1); clf; stairs(t3p,u3p,'Linewidth',1.5); grid on; save prbsdata3 Ts t3 u3; Anexo 3. Programa para la identificación paramétrica fuera de línea. %% Identificación del modelo presión – flujo ----------------------- close all; clear all; clc; Ts = 0.3; planta = iddata(f1,p1,Ts); idplot(planta); lg1 = 614; lg2 = 615; lgf = length(time); ze = planta(1:lg1); % para identificación ze = dtrend(ze); zv = planta(lg2:lgf); % para validación zv = dtrend(zv); m1 = arx(ze,[2 2 0]); m2 = arx(ze,[3 3 0]); m3 = armax(ze,[1 1 1 0]); m4 = armax(ze,[2 2 2 0]); m5 = armax(ze,[3 3 3 0]); m6 = armax(ze,[4 4 4 0]); comp = compare(zv,m1,m2,m3,m4,m5,m6); figure(2); plot(time(lg2:lgf),zv.y,'b','Linewidth',1); hold on; plot(time(lg2:lgf),comp{4}.y,'r','Linewidth',1); grid off; ylabel('F (gpm)','Fontsize',10); xlabel('Tiempo (seg)','Fontsize',10); legend('Datos','Modelo'); axis([time(lg2) time(lgf) -0.3 0.4]); % Parámetros de Polinomio A A1 = [ m1.report.Parameters.ParVector(1:2,1)' 0 0 ]; A2 = [ m2.report.Parameters.ParVector(1:3,1)' 0 ]; A3 = [ m3.report.Parameters.ParVector(1,1)' 0 0 0 ]; A4 = [ m4.report.Parameters.ParVector(1:2,1)' 0 0 ]; A5 = [ m5.report.Parameters.ParVector(1:3,1)' 0 ]; A6 = m6.report.Parameters.ParVector(1:4,1)'; disp('Polinomio A'); disp([A1;A2;A3;A4;A5;A6]); % Parámetros de Polinomio B B1 = [ m1.report.Parameters.ParVector(3:4,1)' 0 0 ]; B2 = [ m2.report.Parameters.ParVector(4:6,1)' 0 ]; B3 = [ m3.report.Parameters.ParVector(2,1)' 0 0 0 ]; B4 = [ m4.report.Parameters.ParVector(3:4,1)' 0 0 ]; B5 = [ m5.report.Parameters.ParVector(4:6,1)' 0 ]; B6 = m6.report.Parameters.ParVector(5:8,1)'; disp('Polinomio B'); disp([B1;B2;B3;B4;B5;B6]); % Parámetros de Polinomio C C3 = [ m3.report.Parameters.ParVector(3,1)' 0 0 0 ]; C4 = [ m4.report.Parameters.ParVector(5:6,1)' 0 0 ]; C5 = [ m5.report.Parameters.ParVector(7:9,1)' 0 ]; C6 = m6.report.Parameters.ParVector(9:12,1)'; disp('Polinomio C'); disp([C3;C4;C5;C6]); % Resultados fit1 = m1.report.Fit.FitPercent; fpe1 = m1.report.Fit.FPE; lfn1 = m1.report.Fit.LossFcn; mse1 = m1.report.Fit.MSE; fit2 = m2.report.Fit.FitPercent; fpe2 = m2.report.Fit.FPE; lfn2 = m2.report.Fit.LossFcn; mse2 = m2.report.Fit.MSE; fit3 = m3.report.Fit.FitPercent; fpe3 = m3.report.Fit.FPE; lfn3 = m3.report.Fit.LossFcn; mse3 = m3.report.Fit.MSE; fit4 = m4.report.Fit.FitPercent; fpe4 = m4.report.Fit.FPE; lfn4 = m4.report.Fit.LossFcn; mse4 = m4.report.Fit.MSE; fit5 = m5.report.Fit.FitPercent; fpe5 = m5.report.Fit.FPE; lfn5 = m5.report.Fit.LossFcn; mse5 = m5.report.Fit.MSE; fit6 = m6.report.Fit.FitPercent; fpe6 = m6.report.Fit.FPE; lfn6 = m6.report.Fit.LossFcn; mse6 = m6.report.Fit.MSE; res = [ fit1 fpe1 lfn1 mse1 fit2 fpe2 lfn2 mse2 fit3 fpe3 lfn3 mse3 fit4 fpe4 lfn4 mse4 fit5 fpe5 lfn5 mse5 fit6 fpe6 lfn6 mse6 ]; disp(' FIT FPE LF MSE '); disp(res); % Gráfico de barras: fit = [ fpe1 fpe2 fpe3 fpe4 fpe5 fpe6 ]; est = [ 1 2 3 4 5 6 ]; figure(3); bar(est,fit); xlabel('Estructura (ítem)','Fontsize',10); ylabel('FPE','Fontsize',10); axis([0 7 0 8e-4]); % Polos y ceros pol11 = pole(m4); zer11 = zero(m4); figure(4); th = 0:pi/50:2*pi; r = 1; xunit = r * cos(th); yunit = r * sin(th); plot(pol11(1),0,'bx'); hold on; plot(pol11(2),0,'bx'); plot(zer11(2),0,'ro'); plot(xunit,yunit,'k'); axis([-2 2 -1.2 1.2]); % Respuesta escalón unitario [y1 x1] = step(m4); figure(5); stem(x1,y1); axis([0 x1(length(x1))-10.9 0 2.2e-3]); ylabel('C (\mus/cm)','Fontsize',10); xlabel('Tiempo (seg)','Fontsize',10); % Respuesta impulso unitario [y2 x2] = impulse(m4); figure(6); stem(x2,y2); axis([0 x1(length(x1))-12.9 0 5e-4]) ylabel('C (\mus/cm)','Fontsize',10); xlabel('Tiempo (seg)','Fontsize',10); %% Identificación del modelo presión – conductividad --------------- close all; clear all; clc; Ts = 4; planta = iddata(c1,p2,Ts); idplot(planta); lg1 = 530; lg2 = 531; lgf = length(time); ze = planta(1:lg1); % para identificación ze = dtrend(ze); zv = planta(lg2:lgf); % para validación zv = dtrend(zv); m1 = arx(ze,[2 2 0]); m2 = arx(ze,[3 3 0]); m3 = armax(ze,[1 1 1 0]); m4 = armax(ze,[2 2 2 0]); m5 = armax(ze,[3 3 3 0]); m6 = armax(ze,[4 4 4 0]); comp = compare(zv,m1,m2,m3,m4,m5,m6); figure(2); plot(time(lg2:lgf),zv.y,'b','Linewidth',1); hold on; plot(time(lg2:lgf),comp{4}.y,'r','Linewidth',1); grid off; ylabel('C (\mus/cm)','Fontsize',10); xlabel('Tiempo (seg)','Fontsize',10); legend('Datos','Modelo'); axis([time(lg2) time(lgf) -70 120]); % Parámetros de Polinomio A A1 = [ m1.report.Parameters.ParVector(1:2,1)' 0 0 ]; A2 = [ m2.report.Parameters.ParVector(1:3,1)' 0 ]; A3 = [ m3.report.Parameters.ParVector(1,1)' 0 0 0 ]; A4 = [ m4.report.Parameters.ParVector(1:2,1)' 0 0 ]; A5 = [ m5.report.Parameters.ParVector(1:3,1)' 0 ]; A6 = m6.report.Parameters.ParVector(1:4,1)'; disp('Polinomio A'); disp([A1;A2;A3;A4;A5;A6]); % Parámetros de Polinomio B B1 = [ m1.report.Parameters.ParVector(3:4,1)' 0 0 ]; B2 = [ m2.report.Parameters.ParVector(4:6,1)' 0 ]; B3 = [ m3.report.Parameters.ParVector(2,1)' 0 0 0 ]; B4 = [ m4.report.Parameters.ParVector(3:4,1)' 0 0 ]; B5 = [ m5.report.Parameters.ParVector(4:6,1)' 0 ]; B6 = m6.report.Parameters.ParVector(5:8,1)'; disp('Polinomio B'); disp([B1;B2;B3;B4;B5;B6]); % Parámetros de Polinomio C C3 = [ m3.report.Parameters.ParVector(3,1)' 0 0 0 ]; C4 = [ m4.report.Parameters.ParVector(5:6,1)' 0 0 ]; C5 = [ m5.report.Parameters.ParVector(7:9,1)' 0 ]; C6 = m6.report.Parameters.ParVector(9:12,1)'; disp('Polinomio C'); disp([C3;C4;C5;C6]); % Resultados fit1 = m1.report.Fit.FitPercent; fpe1 = m1.report.Fit.FPE; lfn1 = m1.report.Fit.LossFcn; mse1 = m1.report.Fit.MSE; fit2 = m2.report.Fit.FitPercent; fpe2 = m2.report.Fit.FPE; lfn2 = m2.report.Fit.LossFcn; mse2 = m2.report.Fit.MSE; fit3 = m3.report.Fit.FitPercent; fpe3 = m3.report.Fit.FPE; lfn3 = m3.report.Fit.LossFcn; mse3 = m3.report.Fit.MSE; fit4 = m4.report.Fit.FitPercent; fpe4 = m4.report.Fit.FPE; lfn4 = m4.report.Fit.LossFcn; mse4 = m4.report.Fit.MSE; fit5 = m5.report.Fit.FitPercent; fpe5 = m5.report.Fit.FPE; lfn5 = m5.report.Fit.LossFcn; mse5 = m5.report.Fit.MSE; fit6 = m6.report.Fit.FitPercent; fpe6 = m6.report.Fit.FPE; lfn6 = m6.report.Fit.LossFcn; mse6 = m6.report.Fit.MSE; res = [ fit1 fpe1 lfn1 mse1 fit2 fpe2 lfn2 mse2 fit3 fpe3 lfn3 mse3 fit4 fpe4 lfn4 mse4 fit5 fpe5 lfn5 mse5 fit6 fpe6 lfn6 mse6 ]; disp(' FIT FPE LF MSE '); disp(res); % Gráfico de barras: fit = [ fpe1 fpe2 fpe3 fpe4 fpe5 fpe6 ]; est = [ 1 2 3 4 5 6 ]; figure(3); bar(est,fit); xlabel('Estructura (ítem)','Fontsize',10); ylabel('FPE','Fontsize',10); axis([0 7 0 20]); % Polos y ceros pol11 = pole(m4); zer11 = zero(m4); figure(4); th = 0:pi/50:2*pi; r = 1; xunit = r * cos(th); yunit = r * sin(th); plot(pol11(1),0,'bx'); hold on; plot(pol11(2),0,'bx'); plot(zer11(2),0,'ro'); plot(xunit,yunit,'k'); axis([-5 1 -1.5 1.5]); % Respuesta escalón unitario [y1 x1] = step(m4); figure(5); stem(x1,y1); axis([0 x1(length(x1)) -0.6 0]) ylabel('C (\mus/cm)','Fontsize',10); xlabel('Tiempo (seg)','Fontsize',10); % Respuesta impulso unitario [y2 x2] = impulse(m4); figure(6); stem(x2,y2); axis([0 x1(length(x1)) -0.015 0]) ylabel('C (\mus/cm)','Fontsize',10); xlabel('Tiempo (seg)','Fontsize',10); %% Identificación del modelo pH – conductividad -------------------- close all; clear all; clc; Ts = 15; planta = iddata(c2,ph2,Ts); idplot(planta); lg1 = 460; lg2 = 461; lgf = length(time); ze = planta(1:lg1); % para identificación ze = dtrend(ze); zv = planta(lg2:lgf); % para validación zv = dtrend(zv); m1 = arx(ze,[2 2 0]); m2 = arx(ze,[3 3 0]); m3 = armax(ze,[1 1 1 0]); m4 = armax(ze,[2 2 2 0]); m5 = armax(ze,[3 3 3 0]); m6 = armax(ze,[4 4 4 0]); comp = compare(zv,m1,m2,m3,m4,m5,m6); figure(2); plot(time(lg2:lgf),zv.y,'b','Linewidth',1); hold on; plot(time(lg2:lgf),comp{4}.y,'r','Linewidth',1); grid off; ylabel('C (\mus/cm)','Fontsize',10); xlabel('Tiempo (seg)','Fontsize',10); legend('Datos','Modelo'); axis([time(lg2) time(lgf) -70 100]); % Parámetros de Polinomio A A1 = [ m1.report.Parameters.ParVector(1:2,1)' 0 0 ]; A2 = [ m2.report.Parameters.ParVector(1:3,1)' 0 ]; A3 = [ m3.report.Parameters.ParVector(1,1)' 0 0 0 ]; A4 = [ m4.report.Parameters.ParVector(1:2,1)' 0 0 ]; A5 = [ m5.report.Parameters.ParVector(1:3,1)' 0 ]; A6 = m6.report.Parameters.ParVector(1:4,1)'; disp('Polinomio A'); disp([A1;A2;A3;A4;A5;A6]); % Parámetros de Polinomio B B1 = [ m1.report.Parameters.ParVector(3:4,1)' 0 0 ]; B2 = [ m2.report.Parameters.ParVector(4:6,1)' 0 ]; B3 = [ m3.report.Parameters.ParVector(2,1)' 0 0 0 ]; B4 = [ m4.report.Parameters.ParVector(3:4,1)' 0 0 ]; B5 = [ m5.report.Parameters.ParVector(4:6,1)' 0 ]; B6 = m6.report.Parameters.ParVector(5:8,1)'; disp('Polinomio B'); disp([B1;B2;B3;B4;B5;B6]); % Parámetros de Polinomio C C3 = [ m3.report.Parameters.ParVector(3,1)' 0 0 0 ]; C4 = [ m4.report.Parameters.ParVector(5:6,1)' 0 0 ]; C5 = [ m5.report.Parameters.ParVector(7:9,1)' 0 ]; C6 = m6.report.Parameters.ParVector(9:12,1)'; disp('Polinomio C'); disp([C3;C4;C5;C6]); % Resultados fit1 = m1.report.Fit.FitPercent; fpe1 = m1.report.Fit.FPE; lfn1 = m1.report.Fit.LossFcn; mse1 = m1.report.Fit.MSE; fit2 = m2.report.Fit.FitPercent; fpe2 = m2.report.Fit.FPE; lfn2 = m2.report.Fit.LossFcn; mse2 = m2.report.Fit.MSE; fit3 = m3.report.Fit.FitPercent; fpe3 = m3.report.Fit.FPE; lfn3 = m3.report.Fit.LossFcn; mse3 = m3.report.Fit.MSE; fit4 = m4.report.Fit.FitPercent; fpe4 = m4.report.Fit.FPE; lfn4 = m4.report.Fit.LossFcn; mse4 = m4.report.Fit.MSE; fit5 = m5.report.Fit.FitPercent; fpe5 = m5.report.Fit.FPE; lfn5 = m5.report.Fit.LossFcn; mse5 = m5.report.Fit.MSE; fit6 = m6.report.Fit.FitPercent; fpe6 = m6.report.Fit.FPE; lfn6 = m6.report.Fit.LossFcn; mse6 = m6.report.Fit.MSE; res = [ fit1 fpe1 lfn1 mse1 fit2 fpe2 lfn2 mse2 fit3 fpe3 lfn3 mse3 fit4 fpe4 lfn4 mse4 fit5 fpe5 lfn5 mse5 fit6 fpe6 lfn6 mse6 ]; disp(' FIT FPE LF MSE '); disp(res); % Gráfico de barras: fit = [ fpe1 fpe2 fpe3 fpe4 fpe5 fpe6 ]; est = [ 1 2 3 4 5 6 ]; figure(3); bar(est,fit); xlabel('Estructura (ítem)','Fontsize',10); ylabel('FPE','Fontsize',10); axis([0 7 0 3]); % Polos y ceros pol11 = pole(m4); zer11 = zero(m4); figure(4); th = 0:pi/50:2*pi; r = 1; xunit = r * cos(th); yunit = r * sin(th); plot(pol11(1),0,'bx'); hold on; plot(pol11(2),0,'bx'); plot(zer11(2),0,'ro'); plot(xunit,yunit,'k'); axis([-1 40 -2 2]); % Respuesta escalón unitario [y1 x1] = step(m4); figure(5); stem(x1,y1); axis([0 x1(length(x1)) -60 0]) ylabel('C (\mus/cm)','Fontsize',10); xlabel('Tiempo (seg)','Fontsize',10); % Respuesta impulso unitario [y2 x2] = impulse(m4); figure(6); stem(x2,y2); axis([0 x1(length(x1)) -0.6 0]) ylabel('C (\mus/cm)','Fontsize',10); xlabel('Tiempo (seg)','Fontsize',10); Anexo 4. Programa para la identificación paramétrica en línea. Identificador recursivo de máxima probabilidad con factor de olvido variable, para la estimación de un modelo con estructura ARMAX de segundo orden. Modelo a estimar: 0.3𝑞−1 + 0.2𝑞−2 1 − 0.1𝑞−1 𝑦(𝑡) = 𝑢(𝑡) + 𝑒(𝑡) 1 − 0.4𝑞−1 − 0.5𝑞−2 1 − 0.4𝑞−1 − 0.5𝑞−2 Diagrama de bloques en simulink: Diagrama del identificador recursivo: Código del identificador recursivo function [yest,theta_next,... lambda_next,varianza_next] = ... Identificador(y,phi,lmax,lmin,var_max,dvar_max,... kappa,theta,lambda,varianza) % Condiciones iniciales rho = 1e6; I = eye(5); persistent P psi dvarianza; if isempty(P) P = rho*I; end if isempty(psi) psi = zeros(5,1); end if isempty(dvarianza) dvarianza = 0; end % Algoritmo yest = phi'*theta; error = y-yest; psi = phi - theta(5,1)*psi; L = (P*psi)/(lambda + psi'*(P*psi)); P = (1/lambda)*(P-L*psi'*P); if (varianza >= var_max) || (dvarianza >= dvar_max) No = 0.1; else No = 10000; end sigma = varianza*No; lambda_n = 1 - (1 - psi'*L)*(error^2/sigma); if lambda_n < lmin lambda_n = lmin; elseif lambda_n > lmax lambda_n = lmax; end varianza_next = kappa*varianza + (1 - kappa)*(error^2); dvarianza = abs(varianza_next - varianza); lambda_next = lambda_n; theta_next = theta + L*error; Resultados del identificador recursivo: Anexo 5. Programa para la simulación de la planta DAM con desacoplamiento. % Modelo con Desacoplador Discreto clear all; clc; s = tf('s'); figure(1); clf; figure(2); clf; %% Función de Transferencia % FT en min G11min = 0.002*(0.056*s + 1)/(0.003*s^2 + 0.1*s + 1); G21min = -0.51*(0.350*s + 1)/(0.213*s^2 + 0.7*s + 1); G22min = -57*(0.320*s + 1)/(0.600*s^2 + 1.8*s + 1); % FT en seg t = 60; G11seg = 0.002*(0.056*(t*s) + 1)/(0.003*(t*s)^2 + 0.1*(t*s) + 1); G21seg = -0.51*(0.350*(t*s) + 1)/(0.213*(t*s)^2 + 0.7*(t*s) + 1); G22seg = -57*(0.320*(t*s) + 1)/(0.600*(t*s)^2 + 1.8*(t*s) + 1); % Tiempo Discreto: Ts = 1; tfinal = 600; ttt = 0:Ts:tfinal; G11d = c2d(G11seg,Ts,'zoh'); G21d = c2d(G21seg,Ts,'zoh'); G22d = c2d(G22seg,Ts,'zoh'); num11 = G11d.num{1}; den11 = G11d.den{1}; num21 = G21d.num{1}; den21 = G21d.den{1}; num22 = G22d.num{1}; den22 = G22d.den{1}; b1_11 = num11(1,2); b2_11 = num11(1,3); a1_11 = den11(1,2); a2_11 = den11(1,3); b1_21 = num21(1,2); b2_21 = num21(1,3); a1_21 = den21(1,2); a2_21 = den21(1,3); b1_22 = num22(1,2); b2_22 = num22(1,3); a1_22 = den22(1,2); a2_22 = den22(1,3); % Desacoplador: G21_22 = G21seg/(G22seg*(2*s+1)); G21_22d = c2d(G21_22,Ts,'zoh'); num21_22 = G21_22d.num{1}; den21_22 = G21_22d.den{1}; b1_21_22 = num21_22(1,2); b2_21_22 = num21_22(1,3); b3_21_22 = num21_22(1,4); b4_21_22 = num21_22(1,5); a1_21_22 = den21_22(1,2); a2_21_22 = den21_22(1,3); a3_21_22 = den21_22(1,4); a4_21_22 = den21_22(1,5); %% Simulación de Planta Discreta esc = 200; retar_esc = 100; % retardo en el escalón u1 = esc*ones(length(ttt),1); for i = 1:(length(ttt)) if i= 1000 F(k,1) = F(k,1) + perturbacion_F; end C(k,1) = y2d(k,1) + y3d(k,1) + ruidoC(k,1); % Entradas er1 = spF(k,1) - F(k,1); er2 = spC(k,1) - C(k,1); uk1 = - a1_c1*uk1_1 + b0_c1*er1 + b1_c1*er1_1; uk2 = - a1_c2*uk2_1 + b0_c2*er2 + b1_c2*er2_1; if uk1 > 200 uk1 = 200; elseif uk1 < -200 uk1 = -200; end % Desacoplador yk21_22 = - a1_21_22*yk21_22_1 - a2_21_22*yk21_22_2 ... - a3_21_22*yk21_22_3 - a4_21_22*yk21_22_4 ... + b1_21_22*uk1_1 + b2_21_22*uk1_2 ... + b3_21_22*uk1_3 + b4_21_22*uk1_4; ydes(k,1) = -yk21_22; uk2r = uk2 - yk21_22; if uk2r > 0.5 uk2r = 0.5; elseif uk2r < -0.5 uk2r = -0.5; end P(k,1) = uk1; pH(k,1) = uk2r; % Respuestas yk11 = - a1_11*yk11_1 - a2_11*yk11_2 + b1_11*uk1_1 + b2_11*uk1_2; yk21 = - a1_21*yk21_1 - a2_21*yk21_2 + b1_21*uk1_1 + b2_21*uk1_2; yk22 = - a1_22*yk22_1 - a2_22*yk22_2 + b1_22*uk2r_1 + b2_22*uk2r_2; % Actualizando variables yk21_22_4 = yk21_22_3; yk21_22_3 = yk21_22_2; yk21_22_2 = yk21_22_1; yk21_22_1 = yk21_22; yk11_2 = yk11_1; yk11_1 = yk11; yk21_2 = yk21_1; yk21_1 = yk21; yk22_2 = yk22_1; yk22_1 = yk22; uk1_4 = uk1_3; uk1_3 = uk1_2; uk1_2 = uk1_1; uk1_1 = uk1; uk2r_2 = uk2r_1; uk2r_1 = uk2r; uk2_2 = uk2_1; uk2_1 = uk2; er1_1 = er1; er2_1 = er2; k = k + 1; end % Condiciones iniciales: Fi = 1; % GPM Ci = 450; % uS/cm Pi = 800; % PSI pHi = 6.5; %% Graficos figure(1); subplot(2,1,1); stairs(tiempo,F+Fi,'b','Linewidth',1.5); grid on; xlabel('Tiempo (seg)'); ylabel('Flujo (GPM)'); title('Respuesta del Flujo debido a la Presión'); subplot(2,1,2); stairs(tiempo,P+Pi,'b','Linewidth',1.5); grid on; xlabel('Tiempo (seg)'); ylabel('Presión (PSI)'); title('Señal de Control u_1 (Presión)'); figure(2); subplot(2,1,1); stairs(tiempo,C+Ci,'m','Linewidth',1.5); grid on; xlabel('Tiempo (seg)'); ylabel('Conductividad (\muS/cm)'); title('Respuesta de la Conductividad debido a la Presión y el pH'); subplot(2,1,2); stairs(tiempo,pH+pHi,'b','Linewidth',1.5); grid on; xlabel('Tiempo (seg)'); ylabel('pH'); title('Señal de Control u_2 (pH)'); figure(3); stairs(tiempo,ydes,'m','Linewidth',1.5); grid on; xlabel('Tiempo (seg)'); ylabel('Amplitud'); title('Señal del Desacoplador'); Anexo 8. Código en texto estructurado de la rutina PlantaDesalinizadora de la tarea Planta, implementada en el software RSLogix 5000. // Modelo de la Planta Desalinizadora // Modelo en ecuación en diferencias //Lectura de Datos: uk1p := DPresion; uk2rp := DpH; // Respuesta del Flujo debido a la Presión yk11 := -a1_11*yk11_1 - a2_11*yk11_2 + b1_11*uk1p_1 + b2_11*uk1p_2; yk11_2 := yk11_1; yk11_1 := yk11; Flujo := yk11 + 1 + Perturbacion_Flujo; // Respuesta de la Conductividad debido a la Presión yk21 := -a1_21*yk21_1 - a2_21*yk21_2 + b1_21*uk1p_1 + b2_21*uk1p_2; yk21_2 := yk21_1; yk21_1 := yk21; // Respuesta de la Conductividad debido al pH yk22 := -a1_22*yk22_1 - a2_22*yk22_2 + b1_22*uk2rp_1 +b2_22*uk2rp_2; yk22_2 := yk22_1; yk22_1 := yk22; Conductividad := yk21 + yk22 + 450; // Almacenando valores anteriores uk1p_3 := uk1p_2; uk1p_2 := uk1p_1; uk1p_1 := uk1p; uk2rp_2 := uk2rp_1; uk2rp_1 := uk2rp; Anexo 9. Código en texto estructurado de la rutina ControlPI_Desacoplo de la tarea Controlador, implementada en el software RSLogix 5000. // Modo Automático IF Modo THEN // Error er1 := Flujo_SP - Flujo; er2 := Conductividad_SP - Conductividad; // Parámetros del controlador PI_1: Ts := 1; a1_c1 := -1; b0_c1 := Kp_1 + 0.5*Ki_1*Ts; b1_c1 := 0.5*Ki_1*Ts - Kp_1; a1_c2 := -1; b0_c2 := Kp_2 + 0.5*Ki_2*Ts; b1_c2 := 0.5*Ki_2*Ts - Kp_2; uk1c := -a1_c1*uk1c_1 + b0_c1*er1 + b1_c1*er1_1; uk2c := -a1_c2*uk2c_1 + b0_c2*er2 + b1_c2*er2_1; // Modelo del Desacoplador en Ecuación en Diferencias yk21_22 := - a1_21_22*yk21_22_1 - a2_21_22*yk21_22_2 - a3_21_22*yk21_22_3 - a4_21_22*yk21_22_4 + b1_21_22*uk1c_1 + b2_21_22*uk1c_2 + b3_21_22*uk1c_3 + b4_21_22*uk1c_4; IF uk1c>100 THEN uk1c:= 100; END_IF; IF uk1c<-100 THEN uk1c:=-100; END_IF; uk2rc := uk2c - yk21_22/100; IF uk2rc>0.5 THEN uk2rc:= 0.5; END_IF; IF uk2rc<-0.5 THEN uk2rc:=-0.5; END_IF; yk21_22_4 := yk21_22_3; yk21_22_3 := yk21_22_2; yk21_22_2 := yk21_22_1; yk21_22_1 := yk21_22; er2_1 := er2; er1_1 := er1; uk1c_4 := uk1c_3; uk1c_3 := uk1c_2; uk1c_2 := uk1c_1; uk1c_1 := uk1c; uk2c_1 := uk2c; pH_desacoplador := yk21_22; DPresion := uk1c; Presion := DPresion + 900; DpH := uk2rc; pH := DpH + 6.5; pH_uc2 := uk2c; // Modo Manual ELSE DPresion := Presion - 900; DpH := pH - 6.5; END_IF;