SCP Notas Curso
SCP Notas Curso
SCP Notas Curso
9 786074 771268
Simulación y Control de
Procesos
Por
Departamento de Energía
Enero 2009
A mi esposa, padres, y
hermanos.
Índice general
1. Introducción 16
1.1. Terminología del control de procesos . . . . . . . . . . . . . . . . . . . . . . . . . 18
1.2. Diagramas de bloques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
1.3. Elementos y configuraciones básicas de un sistema de control . . . . . . . . . . . . 22
1.3.1. Configuración retroalimentada . . . . . . . . . . . . . . . . . . . . . . . . . 22
1.3.2. Configuración antealimentada . . . . . . . . . . . . . . . . . . . . . . . . . 24
1.4. Metodología del diseño de control . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
1.4.1. Objetivos y especificaciones de desempeño . . . . . . . . . . . . . . . . . . 26
1.4.2. Sensores y actuadores . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
1.4.3. Modelado matemático dinámico . . . . . . . . . . . . . . . . . . . . . . . . 32
1.4.4. Controlador . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
1.4.5. Simulación dinámica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
1.4.6. Implementación del controlador . . . . . . . . . . . . . . . . . . . . . . . . 35
1.5. Desarrollo histórico del control de procesos . . . . . . . . . . . . . . . . . . . . . . 36
1.5.1. La revolución industrial . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
1.5.2. Control clásico PID . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
1.5.3. Control moderno . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
1.5.4. Aplicaciones industriales . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
1.5.5. Control post-moderno . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
1.6. Diagrama de flujo del diseño de control . . . . . . . . . . . . . . . . . . . . . . . . 41
1.7. Aplicaciones de control de procesos . . . . . . . . . . . . . . . . . . . . . . . . . . 43
1.8. Referencias básicas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
2
2. Modelado Matemático Dinámico de Procesos 48
2.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
2.2. Formulación de modelos matemáticos . . . . . . . . . . . . . . . . . . . . . . . . . 49
2.2.1. Balances de materia y energía . . . . . . . . . . . . . . . . . . . . . . . . . 50
2.3. Formulación de balances . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
2.4. Formulación matemática de los balances . . . . . . . . . . . . . . . . . . . . . . . 52
2.4.1. Acumulación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
2.4.2. Entrada-salida por convección . . . . . . . . . . . . . . . . . . . . . . . . . 53
2.4.3. Entrada-salida por difusión o dispersión . . . . . . . . . . . . . . . . . . . 54
2.4.4. Entrada-salida por transporte interfase . . . . . . . . . . . . . . . . . . . . 55
2.4.5. Producción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
2.5. Modelos matemáticos dinámicos de procesos . . . . . . . . . . . . . . . . . . . . . 56
2.5.1. Tanques agitados isotérmicos en serie . . . . . . . . . . . . . . . . . . . . . 56
2.5.2. Reactor lote biológico de tanque agitado . . . . . . . . . . . . . . . . . . . 57
2.5.3. Reactor semi-lote de tanque agitado . . . . . . . . . . . . . . . . . . . . . . 60
2.5.4. Reactor continuo de tanque agitado . . . . . . . . . . . . . . . . . . . . . . 61
2.5.5. Reactores continuos de tanque agitado en serie . . . . . . . . . . . . . . . . 63
2.5.6. Intercambiador de calor . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
2.5.7. Reactor tubular flujo tapón . . . . . . . . . . . . . . . . . . . . . . . . . . 67
2.5.8. Reactor tubular con dispersión axial . . . . . . . . . . . . . . . . . . . . . 69
2.5.9. Columna de destilación binaria ideal continúa . . . . . . . . . . . . . . . . 71
2.6. Referencias básicas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
3. Simulación de Procesos 75
3.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
3.2. Simuladores de procesos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
3.2.1. ChemCad1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
3.2.2. Aspen Plus2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
3.3. Programas orientados a la simulación . . . . . . . . . . . . . . . . . . . . . . . . . 80
3
3.3.1. Matlab-Simulink . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
3.4. Lenguajes de programación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
3.4.1. Métodos numéricos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
3.5. Simulación dinámica de procesos con Matlab . . . . . . . . . . . . . . . . . . . . . 87
3.5.1. Sistemas de ODEs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
3.5.2. Sistemas de PDEs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
3.6. Simulación dinámica de procesos . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
3.6.1. Reactor lote . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
3.6.2. Reactor continuo de tanque agitado . . . . . . . . . . . . . . . . . . . . . . 102
3.6.3. Reactor tubular de flujo tapón . . . . . . . . . . . . . . . . . . . . . . . . . 104
3.6.4. Columna de destilación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
3.7. Referencias básicas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
4
4.7.2. Diagramas de bifurcación . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132
4.8. Comportamiento dinámico no-lineal en procesos químicos . . . . . . . . . . . . . . 132
4.8.1. Múltiples estados estacionarios . . . . . . . . . . . . . . . . . . . . . . . . . 132
4.8.2. Oscilaciones simples y complejas . . . . . . . . . . . . . . . . . . . . . . . . 137
4.8.3. Sensibilidad paramétrica y puntos calientes . . . . . . . . . . . . . . . . . . 139
4.8.4. Formación de patrones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140
4.9. Transformada de Laplace en control . . . . . . . . . . . . . . . . . . . . . . . . . . 141
4.9.1. Transformada de Laplace . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143
4.9.2. Transformada de Laplace de un modelo de primer orden con forzamiento . 145
4.9.3. Transformada de Laplace de un modelo de segundo orden con forzamiento 146
4.9.4. Transformada de Laplace con Matlab . . . . . . . . . . . . . . . . . . . . . 146
4.10. Función de transferencia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147
4.10.1. Polos y ceros . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148
4.10.2. Función de transferencia con Matlab . . . . . . . . . . . . . . . . . . . . . 149
4.10.3. Función de transferencia de un CSTR . . . . . . . . . . . . . . . . . . . . . 150
4.11. Respuesta de sistemas lineales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151
4.11.1. Respuesta de sistemas de primer orden . . . . . . . . . . . . . . . . . . . . 152
4.11.2. Respuesta de sistemas de segundo orden . . . . . . . . . . . . . . . . . . . 154
4.11.3. Sistemas de orden . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156
4.12. Identificación empírica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 158
4.12.1. El método de curva de reacción . . . . . . . . . . . . . . . . . . . . . . . . 158
4.12.2. Identificación empírica de una columna de destilación . . . . . . . . . . . . 163
4.12.3. Identificación empírica de un reactor tubular . . . . . . . . . . . . . . . . . 165
4.12.4. Identificación escalón con Matlab . . . . . . . . . . . . . . . . . . . . . . . 165
4.13. Referencias básicas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 167
5
5.2. Implementación numérica del control PID . . . . . . . . . . . . . . . . . . . . . . 174
5.2.1. Control P y PI en Matlab . . . . . . . . . . . . . . . . . . . . . . . . . . . 174
5.2.2. Control PID en Simulink . . . . . . . . . . . . . . . . . . . . . . . . . . . . 176
5.2.3. Controlador PI en un CSTR con Matlab . . . . . . . . . . . . . . . . . . . 179
5.2.4. Controlador PID en un CSTR con Simulink . . . . . . . . . . . . . . . . . 183
5.3. Función de transferencia a lazo cerrado . . . . . . . . . . . . . . . . . . . . . . . . 185
5.3.1. Función de transferencia a lazo cerrado con control P y PI . . . . . . . . . 188
5.4. Localización de raíces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 189
5.4.1. Plano complejo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 189
5.4.2. Polos a lazo cerrado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 190
5.4.3. Diagramas polares . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 190
5.4.4. Criterio de estabilidad de Routh . . . . . . . . . . . . . . . . . . . . . . . . 192
5.4.5. Diagrama de localización de raíces . . . . . . . . . . . . . . . . . . . . . . . 194
5.5. Análisis de la respuesta en frecuencia . . . . . . . . . . . . . . . . . . . . . . . . . 195
5.5.1. Diagramas de Bode . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 199
5.5.2. Diagrama de Nyquist . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 199
5.5.3. Criterios de estabilidad de Nyquist y Bode . . . . . . . . . . . . . . . . . . 201
5.5.4. Incertidumbre y robustez: margen de fase y margen de ganancia . . . . . . 202
5.5.5. Análisis de frecuencia con Matlab . . . . . . . . . . . . . . . . . . . . . . . 204
5.6. Referencias básicas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 208
6
6.4.1. Forma canónica controlable . . . . . . . . . . . . . . . . . . . . . . . . . . 220
6.4.2. Forma canónica observable . . . . . . . . . . . . . . . . . . . . . . . . . . . 221
6.5. Retroalimentado de estados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 221
6.6. Ubicación de polos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 223
6.6.1. Inspección directa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 223
6.6.2. Formula de Ackermann . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 224
6.6.3. Colocación de polos de un CSTR . . . . . . . . . . . . . . . . . . . . . . . 225
6.6.4. Colocación de polos con Matlab . . . . . . . . . . . . . . . . . . . . . . . . 227
6.7. Observadores de estados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 228
6.8. Referencias básicas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 230
7. Epílogo 231
Bibliografía 233
7
Índice de figuras
8
4-1. Estabilidad e inestabilidad de un punto de equilibrio. . . . . . . . . . . . . . . . . 113
4-2. Estabilidad asintótica, local y global. . . . . . . . . . . . . . . . . . . . . . . . . . 124
4-3. Comportamiento dinámico de sistemas de segundo orden en el plano de fase. . . . 131
4-4. Diagrama de Van Herdeen. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134
4-5. Plano de fase del CSTR. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
4-6. Diagrama de bifurcación del CSTR. . . . . . . . . . . . . . . . . . . . . . . . . . . 137
4-7. Oscilaciones periódicas en un CSTR con reacción autocatalítica. . . . . . . . . . . 139
4-8. Oscilaciones caóticas en un CSTR con reacción autocatalítica. . . . . . . . . . . . 140
4-9. Patrón espacio-temporal en un reactor tubular con dispersión axial. . . . . . . . . 142
4-10. Respuesta de un sistema de primer orden. . . . . . . . . . . . . . . . . . . . . . . 153
4-11. Respuesta de un sistema de segundo orden. . . . . . . . . . . . . . . . . . . . . . . 155
4-12. Parámetros de un modelo de segundo orden sub-amortiguado. . . . . . . . . . . . 156
4-13. Respuesta de un sistema de orden . . . . . . . . . . . . . . . . . . . . . . . . . 157
4-14. Ajuste de un modelo de primer orden. . . . . . . . . . . . . . . . . . . . . . . . . . 159
4-15. Ajuste de un modelo de segundo orden sub-amortiguado. . . . . . . . . . . . . . . 161
4-16. Identificación empirica de una columna de destilación. . . . . . . . . . . . . . . . . 163
4-17. Identificación empírica de un reactor tubular biológico. . . . . . . . . . . . . . . . 166
4-18. Simulación a un cambio escalón con Matlab. . . . . . . . . . . . . . . . . . . . . . 166
9
5-11. Estabilidad en el plano imaginario. . . . . . . . . . . . . . . . . . . . . . . . . . . 191
5-12. Diagrama de localización de raíces de un CSTR. . . . . . . . . . . . . . . . . . . . 196
5-13. Diagrama de Bode de una función de transferencia de segundo orden. . . . . . . . 200
5-14. Diagrama de Nyquist de una función de transferencia de segundo orden. . . . . . . 200
5-15. Diagrama de Nyquist para un CSTR a lazo cerrado con = 20 . . . . . . . . . 205
5-16. Diagrama de Nyquist para un CSTR a lazo cerrado con = 225 . . . . . . . . 206
5-17. Diagrama de Nyquist para un CSTR a lazo cerrado con = 30 . . . . . . . . . 206
5-18. Margen de ganancia y margen de fase con Matlab. . . . . . . . . . . . . . . . . . . 207
5-19. Diagrama de Nyquist de un sistema a lazo cerrado marginalmente estable. . . . . 207
10
Lista de Acrónimos
5
Debido a su aceptación universal los acrónimos se derivan de sus siglas en ingles.
11
Prefacio
2. Los objetivos académicos del curso son: (i) introducir conceptos de simulación y control de
procesos, (ii) formular y simular numéricamente modelos matemáticos dinámicos de proce-
sos, (iii) introducir herramientas de análisis dinámico y (iv) diseñar, analizar e implementar
a nivel numérico controladores con técnicas de control clásico y moderno.
La finalidad de estas notas es cubrir los objetivos académicos del curso de simulación y control
de procesos con la profundidad necesaria para cumplir el periodo trimestral del curso. En estas
notas se incluyen herramientas de cómputo del software Matlab, con el toolbox de Simulink. Por
un lado, el uso de paquetes de software libre y comercial, así como la generalización del uso de
equipos de cómputo en cursos de licenciatura, permite a los estudiantes simplificar algunos cálculos
que se hacían en años recientes en forma manual o gráfica. Por otro lado, el uso de software
con herramientas de control de procesos, permiten asimilar y afianzar en forma más rápida los
conceptos y metodologías que se imparten en el curso de simulación y control de procesos.
12
En estas notas se ha tratado de minimizar los desarrollos matemáticos largos y las pruebas
formales de algunos conceptos de control que se presentan, sin embargo, al ser un curso que se
apoya fuertemente en el uso de modelos matemáticos, análisis dinámico y métodos de simulación,
existe un nivel de matemático mínimo que los estudiantes deben de conocer y manejar en forma
adecuada. La teoría de control se ha desarrollado por ingenieros químicos, mecánicos y eléctricos y
los principios de control siempre se expresan en forma matemática y son potencialmente aplicables
a cualquier situación concreta. Así, los conceptos y métodos de simulación y control de procesos
que se abordan en estas notas se aplican si el proceso es químico, biológico, eléctrico o mecánico. Es
importante enfatizar que el éxito en la aplicación de los principios de control automático depende
en gran medida del conocimiento básico que se tenga en el campo de aplicación potencial, ya sea
de ingeniería, física, medicina, economía, ciencias sociales, etc.
El contenido que aquí se incluye se ha proporcionado durante al menos 1 año en el curso de
simulación y control de procesos de la UAM-A, a lo largo del cual se ha mejorado con la retroali-
mentación y respuesta de los estudiantes. En el primer Capítulo se presenta un panorama general
al control de procesos y algunas aplicaciones seleccionadas. En el Capítulo 2 se presentan nociones
de modelado matemático dinámico de procesos y se ilustran con ejemplos clásicos de ingeniería
química. En el Capítulo 3 se presentan herramientas de simulación y control de procesos, haciendo
énfasis en el uso de Matlab. En el Capítulo 4 se presentan conceptos y herramientas de análisis de
sistemas lineales y no-lineales. En el Capítulo 5 se presenta el diseño, implementación numérica y
análisis del control PID. En el Capítulo 6 se presentan nociones de técnicas de control moderno.
En el Capítulo 7 se presentan comentarios finales sobre el contenido que se cubre en estas notas.
Después del Capítulo 7 se presentan las referencias que se usan a lo largo del texto, las cuales
en general no se citan dentro del texto principal por claridad en la presentación, así como otras
referencias que el autor considera de especial interés para profundizar en conceptos y herramientas
de simulación y control de procesos. Al final de cada capítulo se proporcionan además referencias
a textos básicos que se recomiendan para que el lector complemente la información que se pro-
porciona. En la sección de apéndices se presenta una breve introducción a Matlab y Simulink. El
objetivo del autor es que el material que se proporciona sea auto-contenido, interesante, entendible
y útil.
13
Contacto y Agradecimientos
Estas notas se prepararon por el profesor Héctor Puebla del Departamento de Energía de
la UAM-A. Los comentarios y preguntas sobre este documento se pueden dirigir a la dirección
electrónica hpuebla@correo.azc.uam.mx. Actualizaciones y material relacionado a las notas que se
presentan se pueden consultar en la página:
: _
14
de 50 artículos y congresos internacionales arbitrados (Physics Letters A, Industrial Engineering
Chemistry Reserach, Chemical Engineering Science, Systems and Control Letters, IEEE Tran.
Circuits and Systems, IFAC e IEEE Conferences, etc.). El Dr. Puebla es revisor regular de revistas
internacionales (Chemical Engineering Journal, Physics Letters, Bulletin of Mathematical Biology,
etc.) y de tesis de posgrado (UAM, UNAM, IPN, etc.). Sus intereses de investigación incluyen
aplicaciones de control, análisis y control de sistemas biológicos y aplicaciones en ingeniería del
caos. De Septiembre del 2002 a Junio del 2006 el Dr. Puebla se desempeñó como investigador
científico en el Programa de Investigación en Matemáticas Aplicadas y Computación del Instituto
Mexicano del Petróleo. Desde Junio del 2006 el Dr. Puebla se encuentra laborando como profesor-
investigador en el Departamento de Energía de la Universidad Autónoma Metropolitana unidad
Azcapotzalco. Desde el 2004 pertenece al Sistema Nacional de Investigadores y desde el 2005 es
miembro regular del European Society for Mathematical and Theoretical Biology.
15
Capítulo 1
Introducción
Arranque del proceso: Durante el inicio de las operaciones de un proceso es necesario llevar
las unidades del proceso desde una condición inicial determinada a la condición de operación
deseable dictada por los cálculos del diseño de equipos para la producción deseada.
Paro del proceso: En muchos procesos industriales se llevan a cabo uno o dos paros del
proceso anualmente, ya sea por mantenimiento del proceso o por el descenso o aumento
16
temporal de la producción de los productos que ofrecen. El paro de un proceso es una
operación dinámica.
Interacciones complejas entre variables del proceso. El ejemplo más simple de interacción
entre variables es la relación entre la temperatura y la concentración en reactores químicos
no-isotérmicos. Otras interacciones entre variables de un proceso se deben a reacciones auto
catalíticas y sistemas con reciclo de masa y energía.
Variables no medibles. La medición de las variables clave de un proceso puede ser un proble-
ma complicado. En muchas situaciones no es práctico o posible medir directamente variables
del proceso tales como concentraciones, parámetros mecánicos, etc.
Restricciones sobre variables del proceso. Las restricciones sobre variables del proceso están
determinadas principalmente por limitaciones físicas y mecánicas del proceso y factores de
seguridad. Las limitaciones físicas se deben a valores máximos y mínimos que pueden tomar
algunas variables del proceso, tales como concentraciones y flujos negativos. Las limitaciones
mecánicas se deben a la capacidad de equipos mecánicos de implementar en forma rápida
parámetros muy altos o bajos, tales como flujos elevados, temperaturas muy bajas, etc.
Finalmente, debido a que tanto el factor humano como el capital invertido en un proceso son
de máxima prioridad, en muchos casos es necesario imponer valores máximos permisibles
de operación en un proceso, tales como valores límites en temperaturas, concentraciones,
presiones, etc.
17
arriba mencionados. En términos muy generales, el controlar un proceso significa influenciar su
comportamiento de modo que se obtenga una meta deseada.
El control de procesos está presente en una gran variedad de ámbitos, por ejemplo, procesos in-
dustriales, procesos de transporte, procesos de generación y transmisión de energía, mecatrónica,
electrónica, medicina, economía, entre otros. El control de procesos en el ámbito industrial in-
volucra el control de variables tales como temperatura, presión, flujo, etc. Desde una perspectiva
de planeación y administración, involucra maximizar beneficios de operación, garantizar calidad
de producto, incrementar seguridad, reducir el impacto al ambiente y mejorar condiciones de
operación.
El control de procesos es fundamental para lograr:
2. Minimización de desperdicios.
Señales: Una señal es una función del tiempo, que se define típicamente para tiempos mayores
o iguales a cero.
18
Sistema no lineal: Un sistema que contiene elementos no lineales y no satisface el principio
de superposición.
Estado estacionario: Son las condiciones de operación en las cuales no existe cambio respecto
al tiempo. Matemáticamente este corresponde a tener todas las derivadas del tiempo igual
a cero.
Modelo dinámico: Un modelo dinámico es un modelo que describe la evolución a través del
tiempo.
Perturbaciones: Son variables que fluctúan y causan que las variables de un proceso se
muevan del valor de operación deseado.
Variable manipulada o entrada de control: Una variable que se puede elegir de manera que
induzca al sistema a un comportamiento deseado.
Variable controlada: Es la variable que debe ser mantenida, o controlada, en algún valor
deseado.
19
Control retroalimentado: El control utiliza la medición directa de la variable (o variables) a
controlar para ajustar el valor (o valores) de la variable manipulada. Así, la información de
la variable controlada se retroalimenta a la variable manipulada.
Planta: Es conveniente referirse al sistema que se compone por el proceso y el actuador como
la planta.
Lazo abierto: El sistema controlado no tiene conexión entre la salida del sistema y su entrada.
Líneas con flecha: Indican señales o flujo de información. La flecha indica la dirección de la
señal o el flujo de información.
Bloques: Cada bloque denota un sistema que en general se describe con un modelo matemá—
tico.
20
Sistema
Señal con dirección directa
Escalamiento de señal
Función rampa
Comparación de señales
2-a Un sistema con entrada y salida La entrada se procesa en el sistema y genera una
señal de salida . La señal de salida se puede expresar como = ∗ =
2-c Dos sistemas, 1 y 2 , colocados en paralelo. La señal entra en cada uno de los sistemas,
a su vez, la salida de los dos sistemas converge a un punto de sumatoria el cual genera la
salida . En este caso = 1 + 2 .
21
se retroalimenta a la entrada del proceso e influye en el procesamiento subsiguiente que la
produce.
2. Actuador: Los actuadores se usan para manipular alguna variable de entrada que afecte la
variable controlada.
5. Controlador: La función del controlador es, dada una medición y un punto de referencia,
aplicar una señal de salida del controlador adecuada, que logre que el desempeño deseado
del sistema se cumpla a pesar de perturbaciones en el sistema.
Existen dos configuraciones básicas que se usan en el diseño de esquemas de control: anteali-
mentada y retroalimentada.
22
u y (a)
G
(b)
y1 u y2
G1 G2
(c)
G1
+ y
u
+
G2
y (d)
u +
G2
-
G1
23
ajustar la entrada de control para mantener la variable controlada en el valor deseado.
Simplicidad,
La principal desventaja es que las perturbaciones se corrigen hasta que se detectan a la salida
del proceso. La Figura (1-3-a) muestra una configuración retroalimentada básica.
los cambios de la entrada de control se efectúan antes de que las perturbaciones afecten la
señal de salida.
Dos fuertes desventajas son el requerimiento de modelos de perturbaciones, los cuales general-
mente no están disponibles y que este tipo de configuración no tiene capacidad de posicionamiento
y seguimiento de la variable controlada. La Figura (1-3-b) muestra la configuración antealimentada
básica.
24
G5
perturbación
(a)
yref G1 G2 G3
u y
G4
sensor
G5
perturbación
G4 (b)
sensor
yref (-)
G1 G2 G3
u y
25
En la práctica, es común que se combinen las dos configuraciones para el control de un proceso.
La combinación de ambas configuraciones proporciona las ventajas de las dos configuraciones con
la principal desventaja de que se requiere el diseño y operación de dos controladores.
3. Modelar, analizar y simular el comportamiento del proceso real a controlar por un conjunto
de ecuaciones diferenciales.
Objetivos de control
Desde un punto de vista de ingeniería de proceso los objetivos de control pueden ser:
26
1. Regulación: Se entiende como llevar un estado del sistema controlado a algún punto de
operación o equilibrio deseado. Para un proceso que se encuentra operando en estado esta-
cionario la regulación implica que las perturbaciones o ruidos que actúan sobre el sistema
deben tener poco efecto sobre la variable regulada.
3. Seguimiento: Se refiere a llevar un estado del sistema controlado a una trayectoria deseada.
Esto es, algunas variables en el sistema deben de responder en forma particular a señales de
referencia que dependen del tiempo.
Especificaciones de desempeño
Uno de los principales problemas del diseño de un sistema de control es cumplir con ciertos
parámetros de desempeño con el sistema a lazo cerrado. Por un lado, si el sistema con control
reacciona con demasiado impulso pueden surgir problemas de sobre-disparos en las variables con-
troladas y en la entrada de control que se aplica, lo cual puede ser muy riesgoso. Por otro lado, si
el sistema a lazo cerrado muestra una reacción baja la variable controlada puede tardar mucho en
llegar a la referencia deseada, lo cual puede causar problemas de eficiencia en la productividad.
Aun más, es posible que con un esquema de control retroalimentado mal aplicado se produzcan
inestabilidades en un sistema que a lazo abierto era estable, surjan oscilaciones en una respues-
ta previamente monótona, o bien se genere una alta sensibilidad a ruido de medición. Así, una
tarea principal de un ingeniero de control es diseñar un sistema de control que cumpla ciertas
especificaciones de desempeño a lazo cerrado.
Las especificaciones de desempeño describen cómo debe comportarse el sistema con control
o a lazo cerrado, por ejemplo: atenuación de perturbaciones, tiempo de regulación, restricciones
en valores máximos y mínimos de la variable controlada, comportamiento transitorio, etc. En
27
secciones posteriores se revisarán especificaciones de desempeño en un contexto formal usando la
respuesta de un sistema en el dominio del tiempo y con herramientas de análisis en frecuencia.
Sensores
El ingeniero de control debe decidir cuáles señales en el sistema serán medidas y con qué sensor.
En un proceso industrial, por ejemplo el ingeniero de control debe decidir cuáles son las variables
críticas a medir. Además de elegir el tipo de transductor y el hardware para la adquisición de datos,
etc. Las señales que se miden comúnmente en procesos industriales son temperaturas, presiones,
flujos, niveles, velocidad, peso y posición.
Las propiedades deseables de los sensores son:
Confiabilidad.
Exactitud.
Sensibilidad.
Inmunidad a ruido.
No intrusividad.
Los sensores más comunes se usan para determinar estados de movimiento (posición, velocidad,
aceleración y tensión mecánica), de temperatura (temperatura y flujo de calor), de movimiento
de fluido (presión, flujo y nivel) y estados electromagnéticos (voltaje, corriente, carga, luz, flujo
magnético). En gran medida, el desempeño final de un sistema de control está determinado por la
calidad de los elementos de medición que se seleccionen. A continuación se describen los elementos
de medición más comunes en procesos químicos, recordando que la medición de la variable de
interés para fines de control se hace generalmente mediante la combinación de sensor y transmisor.
28
donde la temperatura se mide y están enlazados a un voltímetro en alguna posición remota. La
segunda forma de termómetro utiliza la variación de la resistencia con la temperatura. El tercer
tipo de termómetro se denomina pirómetro y mide la radiación emitida de una superficie caliente.
Medición de flujo Existen varios tipos de flujo. El flujo másico se refiere a la masa de fluido
que pasa en un punto dado por unidad de tiempo (e.g. kg/min). El flujo volumétrico se refiere al
volumen de fluido por unidad de tiempo(e.g. l/seg). En el caso de los gases compresibles, el flujo
volumétrico se referencia a ciertas condiciones de temperatura y presión. Otro tipo de flujo es la
velocidad de flujo, la cual se define como la velocidad de un fluido (e.g. m/seg). El método más
común de medición de flujo genera una caída de presión al pasar a través de una restricción en
una tubería. Cuando existe un gran volumen de flujo se utiliza un medidor de flujo de turbina. La
rotación de las hojas de la turbina se miden por un detector de posición para producir una señal
que es proporcional al flujo. El medidor de flujo de vórtice genera una señal lineal proporcional al
flujo y mide el flujo al detectar pequeños vórtices que se generan aguas abajo de una obstrucción
de un blafe en la trayectoria de flujo. Otro método para medir flujo usa una emisión ultrasónica
para medir el flujo, donde un cambio de frecuencia que se genera por un desplazamiento Doppler
(cambio de la frecuencia con la velocidad) se correlaciona con la medida del flujo. Este método
tiene la ventaja de que no presenta intrusión en la tubería.
Medición de velocidad La forma más común para medir la velocidad es un tacómetro de co-
rriente directa (DC), el cual es simplemente un generador DC con un voltaje de salida proporcional
a la velocidad de rotación.
Medición de peso Existen dos métodos básicos para calcular el peso de un objeto: peso por
tensión y el pesado por balance de fuerzas. En el primero, el objeto a ser pesado distorsiona una
29
estructura de soporte y esta distorsión se mide para dar una indicación del peso. En el segundo
tipo, el peso del objeto se balancea por alguna fuerza (eléctrica, neumática o hidráulica) que se
mide.
Medición de nivel El método más simple usa el hecho de que la indicación de presión de un
líquido es directamente proporcional a la cabeza de líquido encima de el. Otros métodos usan
flotadores (cuya posición se puede medir) y técnicas que usan fuentes radiactivas de bajo nivel.
Actuadores
Una vez ubicados los sensores para informar el estado de un proceso, el siguiente paso es
determinar la forma de actuar sobre el sistema para hacerlo ir del estado actual al estado deseado.
Es el ingeniero de proceso quien al decidir las variables a controlar debe seleccionar el tipo y la
colocación de los actuadores. La selección del elemento final de control se apoya principalmente
en manuales técnicos y especializados que el lector puede consultar con los proveedores de este
tipo de dispositivos y su discusión está fuera de los objetivos de este texto.
Algunas propiedades deseables de los actuadores son:
Potencia.
Velocidad de respuesta.
Precisión de respuesta.
Costos.
Los actuadores se pueden operar en forma electro-mecánica, neumática (los cuales se controlan
por aire a presión), o hidráulicos (los cuales se controlan por la presión de un líquido como aceite
o agua). Los actuadores más comunes en la industria son motores, válvulas, bombas, compresores
y resistencias.
30
trabajo se obtiene de combustibles). A nivel industrial los tipos de motores más usados son de
inducción, continuos, sin escobetilla y paso a paso.
Válvulas Una válvula es un dispositivo que regula el flujo de líquidos, gases o sólidos al abrir,
cerrar completa o parcialmente el paso del fluido. Las válvulas se pueden operar en forma neumáti-
ca, por medio de motores, electro-neumáticamente, etc. Las válvulas de control más comunes son
las que se operan en forma neumática. Las válvulas más conocidas son la válvula de bola, válvula
check o de no retorno (que permite que el fluido pase en una sola dirección), válvula de mari-
posa (común cuando existen tuberías de gran diámetro y con un elevado flujo), válvula solenoide,
válvula de diafragma (que se usa principalmente en la industria farmacéutica), válvula de choque
(la cual controla el flujo a cierto coeficiente de flujo determinado por la apertura de la válvula),
válvula de globo, válvula de alivio (para liberar excesos de presión) y la válvula de reducción de
presión.
31
1.4.3. Modelado matemático dinámico
El uso de modelos en ingeniería química inicio desde 1950, pero el uso de modelos dinámicos,
a diferencia de los modelos en estados estacionario que se usan para análisis y diseño de plan-
tas, es más reciente. En gran medida, la llegada de herramientas más poderosas de cómputo ha
permitido un fuerte desarrollo de modelos dinámicos. Los modelos dinámicos son indispensables
cuando se desea hacer una evaluación de seguridad en una planta, al considerar arranques, paros
y operaciones anormales del proceso, así como en estudios de control y optimización de procesos.
La aplicación del modelado y la simulación permite:
Ayudar en experimentación.
Optimizar procesos.
En control de procesos, los modelos matemáticos dinámicos de procesos son útiles para:
1. Entender el proceso: Para poder lograr un diseño de control exitoso, es necesario tener una
comprensión del proceso a controlar. Esta comprensión usualmente se captura en un modelo
matemático.
2. Diseñar y evaluar sistemas de control: La teoría de control no trata con sistemas físicos reales
sino con la descripción matemática de los mismos a través de modelos matemáticos. Teniendo
un modelo, es posible predecir el impacto de distintos diseños posibles sin comprometer al
sistema real.
Dos enfoques para la construcción de modelos matemáticos dinámicos para fines de control de
procesos son:
Experimental. Se basa en considerar al sistema como una caja negra. En este enfoque se
propone una determinada estructura de modelo, a la que se varían los parámetros, ya sea
por prueba y error, o por algún algoritmo, hasta que el comportamiento dinámico del modelo
se ajusta al observado en la planta mediante ensayos.
32
Analítico. Se basa en el uso de leyes físicas (conservación de masa, energía y momento). El
modelo se obtiene a partir de las leyes fenomenológicas básicas que determinan las relaciones
entre todas las señales del sistema.
Una de las características más importantes del modelado es la necesidad de re-evaluar el modelo
físico y las ecuaciones que representan a un modelo físico con la finalidad de lograr coincidencia,
entre la predicción del modelo y el comportamiento real del proceso modelado.
1.4.4. Controlador
33
1.4.5. Simulación dinámica
Existen el modelo y los métodos pero los procedimientos son arduos y laboriosos que es más
sencillo y menos costoso la simulación.
Se requiere observar un sistema de evolución muy lenta, reduciendo la escala del tiempo.
34
Realizar estudios paramétricos, correr la simulación con diferentes valores de los paráme-
tros.
2. Simulación en tiempo real del sistema con el procesador del control operando a lazo abierto.
3. Simulación en tiempo real del procesador de control conectado al sistema a ser controlado.
35
2. Comunicación: La interconexión de sensores y actuadores requieren el uso de sistemas de
comunicación. El diseño de sistemas de comunicación y sus protocolos asociados es un as-
pecto cada vez más importante de la ingeniería de control moderna. La comunicación entre
dispositivos se puede hacer en forma analógica, ya sea en forma eléctrica o neumática, o por
medio de algún sistema basado en microprocesador.
36
1.5.1. La revolución industrial
En 1922 N. Minorsky introdujo un controlador compuesto por tres términos, conocido pos-
teriormente como el control proporcional-integral-derivativo o PID, para el direccionamiento de
barcos. En 1927, H.S. Black y colaboradores de los laboratorios de la Bell Telephone en los Esta-
dos Unidos desarrollaron la teoría de amplificadores retroalimentados. El problema que se abordó
estaba asociado a la comunicación masiva a grandes distancias, en las cuales se requería amplificar
periódicamente la señal de voz y minimizar efectos de ruido que también se podían amplificar.
Black demostró la utilidad del uso de esquemas de control retroalimentado para reducir la distor-
sión en los amplificadores repetidores.
La característica principal de la etapa del control clásico es que el análisis y diseño de sistemas
de control se baso en el domino de la frecuencia y no en el dominio del tiempo como se hizo en
la etapa de la revolución industrial. El análisis en el dominio de la frecuencia se fundamenta en
los desarrollos matemáticos que se generaron entre 1700 y 1900 por P.-S. de Laplace, J. Fourier y
37
A.L. Cauchy y las contribuciones a inicios de 1900 por H. Nyquist y H.W. Bode. En particular,
Nyquist en 1932, formulo un criterio de estabilidad que se basa en la gráfica polar de una función
compleja. Por otro lado, Bode en 1938 investigo la estabilidad de un sistema con control a través
de la representación de gráficas de respuestas de frecuencia de magnitud y fase de una función
compleja. Todos estos desarrollos y herramientas formaron la base para la elaboración de sistemas
de control que se aplicaron durante la segunda guerra mundial, tales como baterías anti-misiles y
sistemas de control de fuego.
La mayoría de los estudios teóricos y prácticos hasta la segunda guerra mundial fueron sistemas
de una entrada-una salida (SISO), es decir, involucraban el retroalimentado de una sola señal y
la corrección de una sola variable o estado. Después de la segunda guerra mundial, con el avance
de la teoría de control se comenzaron a desarrollar técnicas poderosas que permitían abordar
problemas con más variables a controlar y más de una entrada de control (conocidos como sistemas
multivariables), sistemas variantes en tiempo, así como muchos problemas no lineales.
Con la llegada de la época espacial, el diseño de control, principalmente en los Estados Unidos,
regresó a las herramientas de análisis en el dominio del tiempo que se desarrollaron a finales de
1800, abandonando en esta época las herramientas y técnicas en el dominio de la frecuencia del
control clásico. Los principales motivos de este cambio de propuesta en el diseño de control se
debe a que los problemas de control que surgieron en la época espacial tenían múltiples variables
a controlar y múltiples entradas a manipular, así como interacciones complejas entre las variables
del proceso. La teoría de control clásico o en el dominio de la frecuencia era inapropiada para
manejar esta clase de problemas, ya que los problemas que mejor se podían abordar eran sistemas
SISO.
La mayoría de los métodos de diseño de control en esta época se basaron en resolver algunos
problemas de control en forma optima (control óptimo) que pueden expresarse en la forma de
una ley retroalimentada y en el desarrollo de métodos de cómputo eficientes para resolverlos.
En particular, el control moderno se basa en gran medida en las contribuciones de R. Bellman,
L.S. Pontryagin y R. Kalman. Bellman en 1957 aplicó técnicas de programación dinámica en el
diseño de esquemas de control óptimo. En 1958, Pontryagin desarrolló el principio del máximo,
38
el cual resuelve problemas de control óptimo utilizando cálculo de variaciones. En 1960 Kalman
y colaboradores presentaron dos conceptos fundamentales en la teoría de control moderna: la
contolabilidad y observabilidad de un sistema. Kalman proporcionó además las ecuaciones de
diseño de un controlador conocido como regulador cuadrático lineal (LQR) y la teoría de filtros
óptimos.
39
del entendimiento del proceso, fuertes interacciones entre las variables, no-linealidades del proceso,
pocos modelos matemáticos precisos de los procesos industriales.
Desde 1980 dos técnicas de control han despertado gran interés a la aplicación de sistemas
industriales: (i) el control de modo predictivo (MPC) y (ii) el control robusto.
1. Control MPC: El MPC fue desarrollado en la industria a mediados de los 70s. Con la
capacidad de predicción, un controlador puede adaptarse para prevenir restricciones en las
variables de operación o en las entradas de control, en lugar de reaccionar después de que
haya pasado, como lo hace el PID. El control predictivo esta muy relacionado al control
óptimo. El control predictivo es una metodología de control que trata con dinámica compleja,
interacciones y perturbaciones medibles y es de fácil comprensión al operador de planta.
Además, este método considera restricciones del proceso y de operación.
2. Control robusto: La teoría del control robusto se desarrolló para tomar en cuenta la in-
certidumbre del modelo y analizar el efecto sobre la estabilidad y desempeño del sistema
de control. Esto implica que el controlador se debe de diseñar de modo que las especifica-
ciones de desempeño a lazo cerrado se logren a pesar de condiciones de operación variantes
o diferencias entre el proceso y su modelo.
Otras técnicas importantes de control que surgieron entre los años de 1975 y el 2000, que se
han usado con frecuencia en el control de procesos son:
Control adaptable: El control adaptable, como su nombre lo indica, trata de ajustarse a las
condiciones cambiantes de un proceso, las cuales se reflejan en cambios en los parámetros
del proceso. El control adaptable está provisto con algoritmos que ajustan los parámetros
inciertos con base a información de salida del proceso y un modelo matemático.
40
Control de modo deslizante: El control de modo deslizante forza a los estados del sistema
hacia una región, conocida como superficie de deslizamiento y trata de mantenerlos en ella.
Cualquier movimiento en la superficie deslizante es robusto con respecto a perturbaciones y
errores de modelado. El diseño del control de modo deslizante usualmente consiste de dos
pasos. Primero se diseña la superficie deslizante, de modo que el movimiento del sistema
cumple especificaciones de diseño. Posteriormente, se diseña una función de control para
forzar los estados del sistema hacia la superficie deslizante, lo cual usualmente se logra con
funciones de tipo interrupción.
Control basado en pasividad: Los sistemas pasivos se caracterizan por el hecho de que un
incremento neto en la energía almacenada, en cualquier instante de tiempo, es siempre menor
o igual a la energía entregada al sistema durante el mismo periodo. El concepto de pasividad
es un caso particular de la disipatividad de un sistema, que se define como la diferencia
entre la energía suministrada y el incremento en la energía almacenada en un sistema. El
control basado en pasividad se basa en encontrar una función de almacenamiento de energía,
de modo que el sistema a lazo cerrado sea pasivo y como una consecuencia, se obtiene la
estabilización del sistema.
Modelado matemático: Es necesario identificar los fenómenos más dominantes del proceso,
así como estimar los parámetros del modelo. Si el modelo matemático reproduce las obser-
vaciones experimentales se procede al siguiente paso, de otra forma se tiene que revisar el
modelo matemático.
Diseño de control: Una vez que la respuesta del modelo matemático es aceptable el siguiente
paso es definir la estrategia de control, para lo cual existen una gran variedad de algoritmos
41
Inicio
Identificar
Definir el
fenómenos y
Objetivo de
variables clave
Control y
alternativas
especificaciones
de desempeño
No
Identificar Definir y diseñar
Modelar y simular Si
fenómenos y la estrategia de
dinámicamente
variables clave control
El proceso
¿es aceptable
la respuesta
del modelo?
Modificar
la estrategia de
¿se cumple control
el objetivo y Modificar
especificaciones la estrategia de
de desempeño? control
Si No
Fin
42
de control, dos de las cuales se presentan en estas notas: el control clásico PID y el control
moderno con técnicas estado-espacio.
Un ejemplo practico y cotidiano muy útil para ilustrar los elementos principales de un sistema
de control es la toma de un baño. En este sistema se tienen dos corrientes de entrada de agua, una
caliente y otra fría. Las corrientes se mezclan para producir una corriente de salida a cierto flujo
y temperatura. Las variables manipulables o entradas de control son los flujos de las corrientes de
agua fría y caliente, las perturbaciones al proceso son las temperaturas de ambas corrientes y las
variables a regular pueden ser la temperatura o el flujo de la corriente que se mezcla. Para regular
el flujo y la temperatura de la corriente de salida se utilizan como actuadores las válvulas de las
corrientes de entrada. A través de las terminales nerviosas de la piel, las cuales envían una señal
al cerebro, es posible inferir el nivel de la temperatura de la corriente de salida. Los flujos de las
43
corrientes de entrada se pueden medir con sensores de flujo o en el caso más simple a través de la
sensación de presión de la corriente sobre el cuerpo.
Otro ejemplo típico de textos de control de procesos consiste en el control de nivel de líquido
en un tanque continuo agitado isotérmico. En este caso, el tanque continuo agitado isotérmico
recibe una corriente de entrada a ciertas condiciones de flujo y densidad. La corriente de salida del
tanque tiene una densidad, que en general es idéntica a la que existe dentro del tanque debido a la
agitación en el tanque. Otra consideración que es valida en muchas situaciones es que la densidad
de entrada es idéntica a la densidad de salida. Si el flujo de alimentación es idéntico al flujo de salida
el nivel de líquido en el tanque permanecerá constante. Sin embargo, los cambios en la demanda
del flujo de salida del tanque agitado provocan que exista un cambio en el nivel del líquido en el
tanque. Así, un objetivo de control es regular el nivel de líquido a un valor deseado. Para medir el
nivel en el tanque el dispositivo más sencillo y económico es un flotador. Un controlador determina
la entrada de control, que es el flujo de entrada al tanque agitado, el cual se puede medir con un
sensor de flujo. La implementación de la acción calculada por el controlador se lleva a cabo con
una válvula que regula la entrada de flujo al tanque agitado.
Los reactores continuos de tanque agitado isotérmico son comunes en aplicaciones que involu-
cran reacciones biológicas. La operación isotérmica permite que los microorganismos presentes
en este tipo de reactores no experimenten cambios bruscos de temperatura, los cuales afecten la
actividad de los mismos. Para el caso de un rector biológico para tratamiento de aguas residuales
la alimentación al CSTR isotérmico consiste en un flujo con contaminantes, denominados común-
mente como substratos. Los substratos se degradan por la presencia o adición de microorganismos
en el reactor. La regulación de la concentración de algún substrato o contaminante a un valor
deseado se considera como el principal objetivo de control. La entrada de control es el flujo de
alimentación al reactor. La principal perturbación es el cambio en la concentración de alimentación
de substrato al CSTR. La medición de substrato se puede hacer a través de toma de muestras y
análisis en laboratorio cada cierto tiempo. El actuador es una válvula para regular el flujo a la
44
entrada del CSTR isotérmico.
Se considera ahora un CSTR no isotérmico, el cual se alimenta por una corriente rica en
reactante , a ciertas condiciones de concentración y flujo de alimentación. Dentro del CSTR se
lleva a cabo una reacción exotérmica → → . El reactante se convierte al producto , pero
a altas temperaturas el producto experimenta una reacción adicional que conduce al producto
. El CSTR se enfría por medio de una chaqueta de enfriamiento. El objetivo de control es regular
la concentración de salida del CSTR a un valor de referencia deseado. Debido a la relación directa
entre la temperatura del reactor y la concentración, este objetivo de control se puede reformular
como la regulación de la temperatura de salida del CSTR a un valor de referencia deseado y así
aprovechar que la medición de temperatura es más rápida y económica. La variable manipulable
es el flujo de líquido de enfriamiento a la chaqueta de enfriamiento. Posibles perturbaciones al
proceso son cambios en la concentración y temperatura de la corriente de alimentación. Para
medir la temperatura dentro del CSTR se puede usar un sensor de temperatura. Para medir los
flujos de entrada y salida se usan sensores de flujo. El actuador es una válvula que modifica la
entrada de flujo a la chaqueta de enfriamiento.
Intercambiador de calor
45
Columna de destilación
Se considera una columna de destilación binaria, donde una mezcla binaria se alimenta con un
flujo variable y se tienen corrientes de salida en el fondo y en el domo de la columna. La corriente
en el domo de la columna pasa por un condensador y posteriormente a un tambor de destilado,
donde una parte se retorna al domo de la columna para su enriquecimiento adicional y otra parte
sale hacia otras unidades del proceso con una determinada composición de destilado. Parte de la
corriente que sale en el fondo de la columna se vaporiza y alimenta en el fondo de la columna. El
objetivo de control más simple es mantener la composición de destilado a un valor deseado cuando
la columna está sujeta a cambios en el flujo y la composición de alimentación. Para cumplir con
el objetivo de control la variable manipulable es la velocidad de flujo de reciclo a la columna. Las
composiciones de la columna se miden con sensores de composición y los flujos con sensores de
flujo.
Cadenas de suministro
Dosificación de drogas
46
El ejemplo más conocido del uso de esquemas de control es el tratamiento de la Diabetes tipo 1,
la cual se caracteriza por una deficiencia de la producción de insulina. El objetivo de control se
puede plantear como la regulación del nivel de glucosa en el plasma a través de la administración
de insulina externa. Las perturbaciones en este caso están asociadas al consumo de alimentos que
modifica los niveles de glucosa. La medición de glucosa en la sangre se realiza a través de sensores
de glucosa y la insulina externa se manipula a través de una bomba de dosificación de insulina.
Smith C.A. y Corripio, A.B. (2001). Control Automático de Procesos, Teoría y Práctica.
Limusa.
47
Capítulo 2
2.1. Introducción
Existen varios tipos de modelos disponibles y la naturaleza del sistema a ser descrito indicará
que elección de modelo proporcionará la representación más confiable de la realidad. Los modelos
se pueden clasificar como:
Dinámicos o estáticos: Los modelos dinámicos son los que varían con el tiempo y se usan
principalmente para fines de optimización dinámica y control de procesos. Los modelos
estáticos o en estado estacionario no dependen del tiempo y se usan principalmente en
cuestiones de diseño.
Distribuidos o agrupados: Los modelos distribuidos dependen de dos o más variables (ge-
neralmente tiempo y posiciones) y se representan comúnmente por ecuaciones diferenciales
parciales. Los modelos agrupados dependen de una sola variable (e.g., tiempo en modelos
dinámicos) y se representan por ecuaciones diferenciales ordinarias.
Estocásticos o determinísticos: Los modelos estocásticos son aquellos que contienen ele-
mentos aleatorios y el comportamiento futuro total no se puede predecir, pero si es posible
determinar las propiedades promedio. Los modelos determinísticos son sistemas en los cuales
su comportamiento futuro se determina por su pasado (sistemas causales), el pasado se re-
48
sume usualmente en el estado, es decir, se pueden describir a través de modelos matemáticos
derivados de principios de conservación.
No-lineales o lineales: Los modelos no-lineales son aquellos que contienen elementos no-
lineales y que no cumplen el principio de superposición. Los modelos lineales cumplen el
principio de superposición.
Continuos o discretos: Los modelos continuos son los que contienen variables y funciones
continuas y se describen comúnmente por ecuaciones diferenciales ordinarias o parciales.
Por otro lado, los modelos discretos contienen variables que están discretizadas en intervalos
definidos. En ingeniería de control es común trabajar con modelos discretos, debido a que la
implementación de un diseño de control se hace a través de dispositivos digitales, que hacen
uso de variables discretas, en ciertos intervalos de tiempo (muestreo).
5. Uso del modelo para el fin especifico, por ejemplo, para diseño, control, o cualquier otro
propósito.
49
2.2.1. Balances de materia y energía
Los principios básicos del modelado matemático de procesos son la aplicación de las ecuaciones
de conservación de masa y energía.
Balance de componentes
Para sistemas que tienen más de un componente o especie química el balance de masa se aplica
a cada componente. En este caso, cuando ocurre una transformación química, este cambio se debe
tomar en cuenta. Así, en el caso de masa que se produce por reacción el balance de componentes
se escribe como:
Balance de energía
Los balances de energía son necesarios siempre y cuando ocurran cambios importantes de
temperatura. Por ejemplo, cuando el calor de reacción provoca un cambio en la temperatura del
reactor. Los balances de energía son más complejos que los balances de masa, debido a que son
varios los procesos que pueden provocar cambios en la temperatura. Los balances de energía son
un enunciado de la primera ley de la termodinámica. El balance de energía se puede escribir como:
+(producción de energía)
50
2.3. Formulación de balances
Los pasos para establecer formalmente los balances de energía y masa son:
1. Elección de una región de balance: Se debe de elegir una región de balance donde las variables
del sistema sean constantes o cambien poco dentro del sistema. En este punto es muy útil
dibujar un esquema de la región de balance. La región de balance puede variar en forma
significante, dependiendo del área particular de interés del modelo. Por ejemplo en un reactor
tanque agitado la región de balance puede ser la mezcla reactante del tanque, debido a las
condiciones homogéneas dentro del sistema en cualquier instante de tiempo. En el caso de
un reactor tubular, la región de balance es un elemento diferencial de volumen dentro del
reactor.
2. Identificar corrientes que cruzan la frontera del sistema: Las corrientes que cruzan la frontera
de un sistema se pueden deber a convección, difusión, o procesos de transporte de energía y
masa en interfases.
3. Escribir los balances de masa en forma textual: Los balances de masa y energía se puede
abreviar en este paso como:
Aquí se pueden despreciar algunos términos del balance. Por ejemplo, si el sistema es dinámi-
co o en estado estacionario, si tiene entradas y salidas y si existe producción en el sistema.
4. Expresar cada término del balance formulado en el punto anterior en forma matemática con
variables medibles.
51
2.4. Formulación matemática de los balances
2.4.1. Acumulación
µ ¶
(acumulación de masa de i) = (2.5)
donde está dado en o y el tiempo , está dado en mı́n o La ecuación (2.5) en
términos de variables medibles se puede expresar como:
= (2.6)
= (2.8)
= (2.9)
52
() ( )
(2.10)
donde es el elemento diferencial de volumen.
El término de entrada-salida por convección esta asociado a la masa o energía que entra y
sale del sistema por transporte convectivo. Para el caso de la masa de componentes el término se
escribe como:
0 0 − (2.11)
donde (3 ) es el flujo volumétrico. Para la masa total del sistema se escribe como:
0 0 − (2.12)
0 0 0 0 − (2.13)
Donde los subíndices 0 y se refieren a las condiciones de entrada y salida del sistema respecti-
vamente. Así, el primer término de las ecuaciones anteriores es la entrada y el segundo término
representa la salida.
Para los sistemas distribuidos los términos de convección de masa y energía están dados por:
(0 − ) (2.14)
donde es la velocidad de flujo () que se expresa en unidades de distancia por tiempo,
es el área transversal de transferencia y los subíndices 0 y se refieren a la entrada y salida del
elemento diferencial Al tomar el límite a cero del incremento diferencial se obtienen las expresiones
53
conocidas de los términos convectivos:
( ) ( )
− − (2.16)
Debido a las variaciones de flujo y de velocidad local, efectos de pared y de turbulencia, algunas
moléculas se desplazarán más rápido que otras, de modo que es posible la presencia de dispersión
en el fluido. Para muchos procesos en ingeniería la dispersión es unidireccional, es decir sólo es
importante en una dirección. La contribución por difusión de flujo en procesos de ingeniería se
expresa usualmente por la ley de Fick:
= − (2.17)
donde es el coeficiente de difusión que se expresa en unidades de área por unidad de tiempo
(e.g., 2 ). Por otro lado, la difusión de energía se expresa de acuerdo a la ley de Fourier:
= − (2.18)
El signo negativo en las ecuaciones anteriores se refiere a que la difusión ocurre en la dirección de
concentración o temperatura decreciente.
Los procesos difusivos ocurren solo en sistemas distribuidos, de modo que al formular el balance
entrada-salida por difusión en un elemento diferencial y tomar el límite se llega a las expresión
de difusión para la masa y la energía:
2 2
(2.19)
2 2
Para el transporte difusivo que involucra sólidos porosos es común sustituir el coeficiente de
difusión por un coeficiente de dispersión equivalente, el cual incluye las características de los sólidos
porosos.
54
2.4.4. Entrada-salida por transporte interfase
El transporte de masa interfase es también una posible entrada o salida de masa o energía del
sistema. Comúnmente, el término de transporte en la interfase tiene la estructura:
(área de la interfase)(gradiente)
2.4.5. Producción
55
donde es la velocidad de reacción del componente por volumen y es positivo para productos y
negativo para reactantes Para expresar el término de producción de energía por reacción sólo se
multiplica el término de producción de masa (2.23) por el calor de reacción, ∆ es decir:
∆ = ∆ (2.24)
Un ejemplo clásico para ilustrar diseños de control es el caso de tanques conectados en serie.
Aquí se considera el caso más simple de dos tanques agitados isotérmicos conectados en serie
como se muestra en la Figura (2-1). En el tanque 1 entra una corriente a una velocidad de flujo
volumétrico 10 y una densidad 1 El tanque 1 tiene una altura de líquido 1 La corriente de
salida del tanque 1 1 entra al tanque 2 con una velocidad de flujo volumétrico 20 = 1 con
una densidad 2 El tanque 2 tiene una altura de líquido 2 y el flujo volumétrico de salida es 2 .
La densidad del sistema se considera constante por lo que 1 = 2 = Los balances de masa se
escriben como:
(1 1 )
= 10 − 1 (2.25)
(2 2 )
= 20 − 2
56
proporcional al nivel y a su sección transversal ):
√ √ √
= = (2.26)
1 10 1 p
= − 1 (2.27)
1 1
2 1 p 2 p
= 1 − 2
2 2
Para completar la descripción del modelo, además de asignar valores a los parámetros del
modelo 1 es necesario asignar condiciones iniciales a los estados 1 y 2 es decir:
1 (0) = 1
2 (0) = 2
Para el caso de operación en estado estacionario (, = 0), el modelo anterior se reduce a:
p
10 = 1 1 (2.28)
p p
1 1 = 2 2 = 10
lo cual establece que en estado estacionario, a densidad y área transversal constante, el flujo de
entrada al tanque 1 debe ser igual al flujo de salida del tanque 2.
El segundo ejemplo que se considera es el de un reactor lote biológico de tanque agitado como
el que se muestra en la Figura (2-2). Los reactores lotes se utilizan comúnmente en procesos que
requieren volúmenes pequeños de procesamiento o producción a pequeña escala de alto valor, tales
como los productos farmacéuticos, poliméricos, así como en muchas reacciones biológicas.
57
Corriente de
entrada tanque 1
Nivel en el
tanque 1, h1
Corriente de
salida tanque 1
Corriente de
entrada tanque 2
Nivel en el
tanque 1, h2
Corriente de
salida tanque 2
58
Tiodioglicol +
microorganismos
+ productos
= () − (0) = (2.29)
1
= − () (0) =
59
donde es la constante de velocidad de muerte de los microorganismos (1), es el factor
de producción aparente y () es la cinética enzimática, la cual está descrita por un modelo de
inhibición de substrato del tipo Andrews (Andrews, 1968):
máx
() = (2.30)
1 + +
= ( ) − (2.31)
= () − −
1 1
= − ( ) − () − + ( − )
+
=
60
Substrato
Microorganismos
+ substrato +
penicilina
1máx
( ) = (2.32)
+
2máx
() =
+ (1 + )
61
Reactivo A
Entrada del
liquido de
Enfriamiento Chaqueta de
enfriamiento
5 6 + 2 → 5 7
5 7 + 2 → 5 8 ()2 (2.33)
Considerando densidad constante los balances de masa para los componentes y son:
= (0 − ) − 1 − 3 2 (2.34)
= − + 1 − 2
62
es el volumen del reactor y son las constante cinéticas de reacción que están dadas por
expresiones del tipo Arrenihus:
= 0 exp(− ) = 1 2 3 (2.35)
1 1
= (0 − ) − ∆1 1 − ∆2 2
1
− ∆3 3 2 + ( − ) (2.36)
1
= ( − ) +
donde ∆ es la entalpía de reacción, es la densidad de la mezcla reactante, y son las
capacidades caloríficas de la mezcla reactante y del fluido de enfriamiento respectivamente,
es el coeficiente de transferencia de calor global, es el área de transferencia de calor entre el
reactor y el medio de enfriamiento, 0 es la temperatura de entrada, es la masa del fluido de
enfriamiento y es el calor que se remueve por el medio de enfriamiento. Las condiciones iniciales
son (0) = (0) = (0) = y (0) =
63
Corriente de
entrada reactor 1
Reactor 1
Reactor 2
Reactor 3
Corriente de
salida reactor 3
64
para una reacción → en cada uno de los reactores se pueden escribir como:
1
1 = (0 − 1 ) − 1 (2.37)
2
2 = (1 − 2 ) − 2
3
3 = (2 − 3 ) − 3
donde es la constante de velocidad de reacción de primer orden y son los volúmenes de los
reactores ( = 1 2 3).
Generalizando las ecuaciones anteriores se obtiene:
= (−1 − ) − (2.38)
que describe los balances de masa para reactores en serie. En este caso el reactor recibe como
alimentación la descarga del reactor anterior, es decir − 1 y descarga y alimenta al reactor + 1
Los intercambiadores de calor se usan ampliamente en todo tipo de industrias con la finalidad
de aumentar o disminuir temperatura a corrientes de fluidos. Con las necesidades de integración de
energía de los procesos industriales, los intercambiadores de calor se usan además para intercambiar
temperaturas entre corrientes del proceso.
La descripción dinámica de intercambiadores de calor se puede hacer de dos formas: (i) modelos
de parámetros agrupados (i.e. descritos por ODEs) y (ii) modelos de parámetros distribuidos (i.e.
descritos por PDEs). Debido a que las variaciones de los estados, i.e. temperaturas, ocurren no solo
en el tiempo sino además a lo largo de la posición espacial, los modelos de parámetros distribuidos
se ajustan mejor a la naturaleza de los intercambiadores de calor.
El uso de PDEs para describir a modelos de parámetros distribuidos complica el análisis,
la simulación y el diseño de esquemas de control, debido a que la teoría para PDEs es más
compleja de aplicar que la que existe para ODEs. De esta forma, en muchos casos se prefiere
usar la aproximación de modelos de parámetros agrupados para la descripción de la dinámica de
intercambiadores de calor.
65
Correinte del fluido
de servicio, v2
T1=f(t,z)
1 1
= −1 + (2 − 1 ) (2.39)
1 1
2 2
= −2 − (2 − 1 )
2 2
donde es el coeficiente de transferencia de calor global, 1 y son las capacidades caloríficas
y las masas totales de las corrientes ( = 1 2) respectivamente.
Para completar la descripción del modelo es necesario establecer las condiciones frontera y
condiciones iniciales. Las condiciones frontera más comunes son condiciones frontera de Dirichlet:
66
de Newmann:
( )
=0 (2.41)
y de Danckwerts:
( )
( ) − = 0 (2.42)
En este caso las condiciones frontera para los términos convectivos están dadas por:
Los reactores tubulares son de los más importantes para las industrias biotecnológicas, de
refinación del petróleo y poliméricas. Un reactor tubular es cualquier recipiente cilíndrico a través
del cual los reactantes fluyen y ocurre una reacción química o bioquímica en su interior. Es decir,
es una operación continua en la cual hay un movimiento constante de uno o de todos los reactivos
en una dirección espacial promedio y en el que no se hace intento alguno por inducir el mezclado
entre los elementos del fluido en ningún punto a lo largo de la dirección de flujo.
El uso de reactores tubulares biológicos de cama empacada (i.e. enzimas o microorganismos
inmovilizados) para algunas reacciones bioquímicas permite mayor estabilidad térmica, uso de un
mayor volumen celular, más eficiencia en la extracción de productos y velocidades de flujo mayores
respecto a otro tipo de reactores. Algunos ejemplos son la remoción de contaminantes de efluentes
industriales, la producción de antibióticos y la producción de etanol.
Se considera un reactor biológico tubular isotérmico de cama empacada como el que se muestra
en la Figura (2-7), donde se llevan a cabo dos reacciones, una de crecimiento autocatalítico de
microorganismos y la otra de desactivación de microorganismos:
→ → (2.45)
67
Productos +
Substrato microorganismos
desactivados
z=0 z=L
Cama de microorganismos
= ( ) −
= − − 1 ( ) (2.46)
= − +
máx
( ) = (2.47)
+
68
y condiciones frontera tipo Dirichlet:
El mecanismo por el cual el calor o la masa se dispersan en una cama empacada a través de la
cual un fluido está fluyendo ha sido objeto de una gran atención. Se ha supuesto que la dispersión
de calor y masa en la dirección axial y radial es un proceso difusivo sobre-impuesto al flujo
convectivo, de esta forma el efecto global se describe por un conjunto de ecuaciones diferenciales
parciales parabólicas. Los reactores tubulares en los cuales es importante el término de dispersión
son principalmente los reactores catalíticos de cama empacada. Algunos ejemplos del uso de este
tipo de reactores son la síntesis de amoníaco, la producción de ácido sulfúrico y muchos procesos
de refinación del petróleo.
Los modelos seudo-homogéneos son los más simples para el modelado del reactor de cama
fija catalítico. La consideración básica que se hace es que el reactor puede describirse como una
entidad que consiste solamente de una fase simple. En el modelo se supone que toda la superficie
catalítica está totalmente expuesta a las condiciones del medio y que tiene las mismas condiciones
a las del fluido y puede entonces ser completamente descrito por las variables del medio (tempe-
ratura, concentración, presión). Esto es estrictamente válido cuando el factor de efectividad (i.e., la
relación de la velocidad de reacción con resistencias difusionales sobre la velocidad de reacción a las
condiciones del medio) es igual a la unidad ( = 10) para todas las reacciones y/o componentes.
Esto permite que las velocidades intrínsecas de reacción se usen en el modelo de reacción sin
ninguna modificación.
El sistema que se considera se muestra en la Figura (2-8) y es un reactor tubular no isotérmico
en el cual se lleva a cabo una reacción exotérmica de primer orden irreversible → En términos
de los parámetros adimensionales de la Tabla 1, los balances de materia y energía están dados
por: ½ ¾
1 1 2 1 1 2
= − + (1 − 1 ) exp (2.50)
1 2 1 + 2
69
Fluido de enfriamiento
x2w, Temperatura de enfriamiento
Reactivo A Productos
z=0 z=L
Cama de Catalizador
Cuadro 2.1: Parámetros adimensionales del reactor tubular con dispersión axial
½ ¾
2 1 2 2 1 2 2
= 2
− + (1 − 1 ) exp (2.51)
2 1 + 2
+ (2 − 2 )
70
y de no deslizamiento o de Newmann a la salida del reactor:
¯
1 ( ) ¯¯
= 0 (2.53)
¯=
¯
2 ( ) ¯¯
= 0
¯=
El último ejemplo que se presenta es el caso de una columna de destilación binaria ideal
continua. La destilación es una de las operaciones unitarias más comunes en la industria del
petróleo. En la columna de destilación una mezcla que entra a la columna se separa en forma de
vapor (el componente más volátil o con menor punto de ebullición) y se eleva a lo largo de la
columna hasta llegar al domo. Parte de la alimentación (el componente menos volátil) desciende
a lo largo de la columna hasta alcanzar el fondo. El vapor en el domo de la columna, que es más
rico en el componente más volátil, se envía a un condensador. Para aumentar la calidad de los
componentes que se separan, es común que parte del vapor que se condensa en el domo de la
columna se recicle a la columna, así como parte del líquido se recicla a la columna en forma de
vapor después de pasar por un re-hervidor.
Se considera el modelo clásico presentado por Luyben de una columna continua de destilación
binaria ideal como se muestra en la Figura (2-9). La columna tiene platos teóricos ideales. Una
mezcla binaria con flujo molar y composición de alimentación se alimenta en el plato .
En el fondo de la columna (etapa + 1) se tienen moles con fracción mol de líquido La
corriente del líquido del fondo se separa en una parte que envía a un re-hervidor que alimenta a
la columna con vapor a una velocidad de flujo y composición y otra parte que sale con una
velocidad de flujo y composición El vapor en el domo de la columna (etapa 1) pasa por un
condensador total y el condensado se envía a un tambor de reflujo, el cual tiene una retención de
líquido con composición de líquido Parte del condensado se envía al plato 1 de la columna
como un reflujo a una velocidad 0 y una composición El producto del tambor de condensado
se remueve con una composición y una velocidad La retención de líquido en cada plato
de la columna es con composición de líquido y composición de vapor correspondiente
Algunas simplificaciones para formular el modelo matemático son: (i) presión constante, (ii)
el sistema es ideal, con el equilibrio descrito por una volatilidad relativa constante, (iii) no existen
71
Corriente
liquida de
entrada
Corriente de
Corriente de destilado
alimentación
Corriente de vapor
de entrada
Corriente
liquida de
salida Fondos
72
pérdidas de calor, (iv) los platos de la columna tienen 100 % de eficiencia y se comportan como
platos teóricos y (v) el reflujo y la alimentación están saturados ( = 1 donde es la calidad
de la alimentación). Las ecuaciones del modelo de la columna de describen a continuación.
En el tambor de reflujo se tiene:
= 1 − ( + ) (2.54)
en la sección de enriquecimiento ( = 1 2 ):
= (−1 − ) + (+1 − ) (2.55)
=
= ( + 1)
en el plato de alimentación:
= + −1 − 0 + +1 − 0 (2.56)
0 = + = +
0 = − (1 − ) =
= 0 (−1 − ) + 0 (+1 − ) (2.57)
en el re-hervidor:
= 0 − − 0 (2.58)
= 0 − 0
73
El equilibrio entre las fases líquida y gas para cada plato ( = 1 2 ) se expresa como:
= (2.59)
1 + ( − 1)
Luyben W.L. (1991). Process. Modelling, Simulation and Control for Chemical Engineers.
McGraw Hill.
74
Capítulo 3
Simulación de Procesos
3.1. Introducción
La simulación de procesos es la representación de un proceso por medio de un modelo matemáti-
co y su solución numérica para obtener información sobre el comportamiento del proceso.
La simulación de procesos se utiliza en:
Pruebas de planta piloto: El uso de simuladores puede ayudar a obtener buenas estimaciones
de las condiciones de operación del proceso real a partir de datos obtenidos en plantas piloto.
Diseño de equipos: El simulador de procesos proporciona los datos requeridos para el diseño
detallado de diferentes equipos.
Simulación de plantas existentes: La simulación de una planta existente puede ser útil si
existe la necesidad de cambiar las condiciones de operación.
La simulación de procesos también puede ser útil para encontrar la mejor estrategia para
elevar o disminuir la producción, aumentar la eficiencia de la operación a través de la integración
75
de energía, reconfigurar la planta para diferentes materias primas o requerimientos de productos
con diferentes composiciones.
La simulación de procesos se puede formular en 5 pasos:
1. Identificar el problema. En este paso se debe de definir las características del problema a
ser estudiado. Si es un problema de diseño y optimización aquí se define que se trata de
una simulación en estado estacionario. Por otro lado, si se trata de un problema de control
de procesos se tiene que realizar una simulación dinámica. En general, los modelos que se
usan para diseño de equipos son muy detallados y los modelos que se usan para control de
procesos pueden ser simples para no hacer muy complejo el problema del diseño de control.
3. Seleccionar la plataforma de simulación: Para los ingenieros químicos, existen tres platafor-
mas principales para realizar la simulación de un proceso,
76
Los programas orientados a simulación más conocidos son el Mathematica, Matlab,
Berkley-Madonna y Polymath.
5. Analizar los resultados obtenidos. Es común que al comenzar a usar alguna plataforma de
simulación se generen errores en la alimentación de las ecuaciones del modelo y la asignación
de parámetros, si esto ocurre, es probable que no se generen resultados de la simulación.
Por otro lado, si los resultados que se obtienen no tienen sentido físico, en este punto se
pueden identificar posibles errores e inconsistencias en los modelos y datos alimentados. La
evaluación formal de los resultados que se obtienen de una simulación se puede completar
con las herramientas descritas en el siguiente capítulo.
77
Puede usarse para modelar procesos típicos a los encontrados en las industrias químicas y
petroleras.
Dos de los simuladores de procesos de mayor uso son el Aspen Plus y el ChemCad.
3.2.1. ChemCad1
CC-Steady State: El cual incluye una base de datos de componentes químicos, métodos
termodinámicos y operaciones unitarias que permiten realizar la simulación en estado esta-
cionario de procesos químicos.
78
3.2.2. Aspen Plus2
Aspen Plus puede ser muy útil para diseñar mejores plantas e incrementar la producción en
plantas existentes, se pueden cambiar especificaciones, tales como la configuración del diagrama
de flujo, condiciones de operación y composición de la alimentación, para correr nuevos casos y
analizar alternativas y para analizar los resultados, se puede generar gráficas, reportes y archivos
extendidos.
Además de las características generales de los simuladores comerciales, con Aspen Plus se
pueden ejecutar un amplio rango de tareas adicionales. Entre ellas: (i) Realizar análisis de sensi-
bilidad y casos de estudio. (ii) Generar gráficas y tabular los resultados. (iii) Cálculo y regresión
de propiedades físicas. (iv) Evaluar los costos de la planta. (v) Optimizar procesos. (vi) Obtener
parámetros de diseño de equipo.
Aspen Plus además contiene datos, propiedades, modelos de operaciones unitarias, reportes
y otras características y capacidades desarrolladas para aplicaciones especificas a varios tipos de
industrias. Por ejemplo, para aplicaciones de simulación dinámica y control de procesos se tienen
los módulos de Aspen Dynamics y Aspen HYSYS Dynamics.
Aspen HYSYS. Este módulo esta dirigido para la simulación dinámica de procesos de la
industria del petróleo, procesamiento de gas y refinación. Aspen HYSYS tiene una base de
datos muy completa para el cálculo preciso de propiedades físicas, propiedades de transporte
y comportamiento de fase para las industrias del petróleo y refinación. Los modelos de las
operaciones unitarias son escalables y permiten múltiples niveles de detalle en geometría de
equipos para aumentar el nivel de fidelidad de las aplicaciones. Aspen HYSYS incluye varias
herramientas para el monitoreo y control de procesos, tales como diferentes controladores y
79
operaciones lógicas.
80
3.3.1. Matlab-Simulink
Matlab fue creado originalmente para resolver problemas matriciales, sin embargo, actualmente
es posible resolver prácticamente cualquier problema numérico de ciencias e ingeniería a través de
paquetes específicos (toolboxes), así como visualizar datos en una gran diversidad de presentaciones
gráficas. Matlab tiene además capacidades para escribir programas en un lenguaje de programación
sencillo, que comparte muchas características con otros lenguajes de programación, tales como
Visual Fortran, Visual Basic y C.
Para resolver sistemas de ecuaciones diferenciales que describen a un proceso en Matlab se
tienen dos opciones: (i) solución con librerías de métodos numéricos para resolver ODEs. (ii)
solución con la programación de los métodos numéricos.
Uno de los paquetes de herramientas de Matlab que es de gran utilidad para realizar cálculos
de simulación y control de procesos es el Simulink, el cual es un programa dirigido principalmente
a la simulación de sistemas dinámicos, ya que cuenta con características adicionales para sistemas
dinámicos que complementan las herramientas de objetivo general de Matlab. En Simulink el
usuario puede alimentar las ecuaciones diferenciales directamente, o bien a través del uso de iconos
gráficos para determinadas funciones (derivada, integración, sumatoria y funciones determinadas,
tales como cambios escalón y pulso, entre otras) y enlazarlas en un diagrama de flujo en forma
lógica.
81
Fortran se dispone del paquete IMSL que es una biblioteca que contiene más de 1000 subrutinas
de análisis numérico escritas en Fortran, entre las que se encuentran solución a sistemas lineales,
sistemas de valores y vectores propios, interpolación y aproximación, integración y diferenciación,
ecuaciones diferenciales, ecuaciones no lineales, optimización. La librería IMSL cuenta además con
funciones especiales de amplio uso en cálculos de ciencias e ingeniería.
En esta Sección se presentan los fundamentos de tres métodos numéricos que se usan común-
mente en el análisis y la simulación de procesos: (i) El método de Newton-Raphson, (ii) el método
de Euler y (iii) el método de Runge-Kutta. El método de Newton-Raphson se usa para encontrar
las raíces o puntos de equilibrio de un modelo. Los métodos de Euler y de Runge-Kutta son méto-
dos para la solución numérica de ecuaciones diferenciales ordinarias o problemas de valor inicial.
Se presenta además el método de discretización por diferencias finitas, el cual forma la base de
los tres métodos numéricos mencionados, así como para resolver en forma numérica sistemas de
ecuaciones diferenciales parciales.
Diferencias Finitas
El método de diferencias finitas trata con puntos específicos dentro un intervalo o dominio,
denominados puntos malla. El dominio está dividido en intervalos equidistantes, aunque tal con-
sideración de intervalos iguales no es necesaria. Si se supone que se tiene una función continua
() al usar una expansión de Taylor se pueden deducir fórmulas en diferencias para la primer
y segunda derivada involucrando solamente los valores en −1 y +1 donde indica el punto
malla de la discretización.
Desarrollando la serie de Taylor para (+1 ) y (−1 ) con un incremento de discretización
∆ se obtiene:
82
derivada:
(+1 ) − ( ) ( ) ∆ 2 ( )
= + + (3.3)
∆ 2! 2
( ) − (−1 ) ( ) ∆ 2 ( )
= − + (3.4)
∆ 2! 2
cada expresión es correcta hasta el orden de ∆ También se pueden sumar las ecuaciones (3.3) y
(3.4), reordenar y dividir por ∆ con lo que se obtiene:
La cual es correcta hasta el orden de ∆2 Si se restan las ecuaciones (3.4) y (3.3), reordenan y
dividen por ∆2 se obtiene una expresión para la segunda derivada:
Despreciando los términos de alto orden en las expresiones (3.5) y (3.6) se obtiene:
las cuales son las formas conocidas como primer y segunda derivada centradas, las cuales utilizan
información de puntos de discretización antes y después del punto de evaluación.
El método de diferencias finitas tiene la ventaja de su fácil formulación, sin embargo es necesario
en algunos casos un elevado número de puntos malla para obtener una buena precisión.
Método de Newton-Raphson
83
y
f(x) df(x)/dx
x
x* x2 x1 x0
( )
(+1 ) ≈ ( ) + (+1 − ) (3.9)
el proceso se puede repetir hasta que converge a un punto fijo o raíz, en donde se cumple que
(+1 ) = 0 El algoritmo que se obtiene es:
∙ ¸−1
( )
+1 = − ( ) (3.10)
La idea del método es comenzar con una punto de inicio cercano a la raíz, se reemplaza la
función por la recta tangente en ese valor, se iguala a cero y se despeja y se aplican las iteraciones
que se deseen. La Figura (3-1) muestra la aproximación numérica del método de Newton-Raphson
a la raíz ∗ de una función = ().
Método de Euler
El método de Euler es el método más simple para resolver problemas de valor inicial de la
forma general:
= ( ) (0 ) = 0 (3.11)
84
dxn+2/dt
dxn+1/dt
y
dxn/dt
f(x)
t
tn tn+1 tn+2
en un intervalo de [0 ] dada una condición inicial 0 El método en términos generales consiste
en tres pasos: (i) discretizar la derivada del tiempo , (ii) representar () en los puntos de
discretización, ( ) y (iii) aproximar la derivada en cada punto de discretización usando
los valores de () actuales y anteriores, ( ). La Figura (3-2) muestra las aproximaciones
numéricas que se hacen con el método de Euler en los puntos de discretización a una función ()
Una forma de obtener el algoritmo es a partir del desarrollo de Taylor de +1 en términos del
incremento de tiempo de discretización ∆ :
∆ 2 2
+1 ≈ + ∆ + + (3.12)
2! 2
+1 ≈ + ∆ para = 0 1 2 − 1 (3.13)
donde:
= ( ) (3.14)
La idea es iniciar el algoritmo (3.13) con el valor inicial que se conoce 0 , con lo cual se ob-
tienen valores actualizados +1 . El incremento de tiempo de la discretización se debe seleccionar
85
cuidadosamente y comúnmente depende del compromiso entre la precisión deseada y el costo com-
putacional. En problemas prácticos el método de Euler no es muy útil porque para obtener una
precisión razonable se necesita un paso ∆ muy pequeño.
Método de Runge-Kutta
X
+1 ≈ + ∆ (3.15)
=1
donde:
= ( + ∆ + ∆) = 1 2 (3.16)
y , , son constantes propias del esquema numérico. Los esquemas Runge-Kutta pueden ser
explícitos o implícitos dependiendo de las constantes del esquema.
Para el método de Runge-Kutta de cuarto orden, el cual es el algoritmo más usado para resolver
numéricamente un sistema de ecuaciones diferenciales ordinarias se tiene:
∆
+1 = + (1 + 22 + 23 + 4 ) (3.17)
6
donde:
1 = ( )
∆ ∆
2 = ( + 1 + )
2 2
∆ ∆
3 = ( + 2 + )
2 2
4 = ( + ∆3 + ∆)
En el algoritmo (3.17) el valor +1 está determinado por el valor presente más el producto
86
del tamaño del intervalo de tiempo ∆ por un promedio ponderado de pendientes.
Para resolver sistemas de ecuaciones diferenciales no-lineales en Matlab se tienen dos opciones:
(i) Solución con librerías de métodos numéricos para resolver ODEs. (ii) Solución con la progra-
mación de los métodos numéricos.
Para resolver sistemas de ecuaciones diferenciales ordinarias es conveniente crear dos archivos
desde el editor de Matlab, conocidos como archivos m por la extensión que utilizan. En el primero de
ellos se definen los parámetros y las ecuaciones a resolver. En el segundo se escriben la instrucción
que llama a la librería del método numérico a utilizar, las condiciones iniciales, el tiempo de
integración y se le proporciona formato a los resultados de la solución numérica.
Matlab cuenta con varios integradores de ecuaciones diferenciales ordinarias, los de mayor uso
son el 45 y el 23, los cuales se basan en métodos de Runge-Kutta por parejas. El 45
se recomienda para problemas en los cuales no existe rigidez numérica y se desea una precisión
moderada con bajo costo computacional. El 23 es adecuado para problemas que presentan
problemas de rigidez numérica elevada.
Para ilustrar la solución de sistemas de ODEs en Matlab considere el ejemplo modelo de Lotka-
Volterra, el cual se origino en los trabajos independientes de Lotka y Volterra en los años 20. El
matemático italiano Volterra estudio los registros de las pesquerías del Mar Adriático Superior y
observó que cuando la pesca había disminuido drásticamente, la proporción de los depredadores
87
había aumentado. Este hecho lo llevo a estudiar ese problema de una manera más general, lo
cual derivó en la primer teoría sistematizada de la dinámica de poblaciones. El modelo de Lotka-
Volterra se puede describir por las siguientes reacciones hipotéticas:
→1 +
+ →2 + (3.18)
→3
= 1 − 2 (3.19)
= 2 − 3 2
al tomar parámetros adimensionales se obtienen las ecuaciones conocidas del modelo de Lotka-
Volterra:
1
= 1 1 − 2 1 2 (3.20)
2
= 3 1 2 − 4 2
donde los parámetros son positivos. Para completar la descripción del modelo se definen las
condiciones iniciales como:
2 (0) = 2
88
function = ( );
% asignación de parámetros
1 = 10; 2 = 20; 3 = 20; 4 = 10;
% crear un arreglo matricial de renglones por 1 columna
% con la estructura dx=zeros(número de ecuaciones,1);
= (2 1);
% ecuaciones del modelo de lotka-volterra
(1) = 1 ∗ (1) − 2 ∗ (1) ∗ (2);
(2) = 3 ∗ (1) ∗ (2) − 4 ∗ (2);
figure(1);
( (: 1) (: 2)), (’Presas, x_1, Predadores x_2’), (’Tiempo, t’)
La solución que ofrece Matlab con los resolvedores que tiene programados en algunos casos
presentan problemas de convergencia. Aún más, al usar los métodos de Matlab no es sencillo
almacenar y presentar resultados de variables que no son los estados del sistema. Una alternativa
para tales casos es usar las capacidades de programación de Matlab y programar rutinas de métodos
numéricos de integración de ODEs. Uno de los métodos más eficientes para tal efecto es el Runge-
Kutta de cuarto orden.
Para el ejemplo del modelo de Lotka-Volterra, la estructura de los archivos m es la siguiente:
89
Archivo de ecuaciones: ()
function = ( )
% asignación de parámetros
= 1 :
() = ∗ ;
1 = ∗ (@ ());
= ();
= 1 : 2
90
() = () + (1() + 2 ∗ 2() + 2 ∗ 3() + 4())6;
1() = (1);
2() = (2);
figure(1)
Para resolver sistemas de ecuaciones diferenciales parciales en Matlab se tienen dos opciones:
(i) Solución con librerías de métodos numéricos para resolver PDEs. (ii) Solución con librerías
de métodos numéricos para resolver ODEs con la discretización de derivadas con el método de
diferencias finitas.
91
2.5
1.5
0.5
0
0 5 10 15 20 25 30 35 40 45 50
Tiempo, t
La solución con librerías de Matlab se puede usar para resolver un sistema de PDEs de la
forma: µ ¶ µ µ ¶¶ µ ¶
−
= + (3.22)
para:
0 (3.23)
≤ ≤
dos condiciones que deben de cumplirse son: (i) los elementos de la diagonal de la matriz de la
ecuación (3.22) deben ser todos cero o todos positivos y (ii) los elementos de la matriz de la
92
ecuación (3.25) deben ser todos cero o todos tomar valores diferentes de cero.
La sintaxis de la instrucción para resolver sistemas de PDEs es la siguiente:
es una sub-función m donde se definen las condiciones iniciales y tiene la sintaxis:
0 = ()
es una sub-función m donde se definen las condiciones frontera y tiene la sintaxis:
[ ] = ( ) Los subíndices y se refieren a las
condiciones a la entrada y salida del sistema.
Se considera el ejemplo del reactor tubular con dispersión axial, descrito en la Sección 2.5, que
se modela por las ecuaciones (2.50)-(2.53). En términos de la estructura general de solución con
librerías de Matlab. Si se toma 1 = 1 y 2 = 2 las ecuaciones (2.50)-(2.53) se pueden escribir
como:
⎡ ⎤ ⎡ ⎤ ⎡ ⎤
1
1 1 (1 )
⎣ ⎦ ⎣ ⎦ = ⎣ 1 ⎦+
1 2 1
2
(2 )
⎡ ⎤
−1 +
⎢ ⎥
⎢ ⎥
⎢ (1 − 1 ) exp( 1+22 ) ⎥
⎢ ⎥ (3.26a)
⎢ 1 ⎥
⎢ − 2 + (1 − 1 ) ⎥
⎣ ⎦
exp( 1+22 ) + (2 − 2 )
⎡ ⎤ ⎡ ⎤⎡ ⎤ ⎡ ⎤
1 − 10 − 11 1 0
⎣ ⎦+⎣ ⎦⎣ ⎦ = ⎣ ⎦ (3.26b)
2 − 20 − 12 2 0
⎡ ⎤ ⎡ ⎤⎡ ⎤ ⎡ ⎤
0 1 1 0
⎣ ⎦ + ⎣ ⎦⎣ ⎦ = ⎣ ⎦ (3.26c)
0 1 2 0
93
La estructura de los archivos m para los valores de los parámetros 1 = 50 2 = 50
= 1 = 0158 = 110 = 00 = 200 2 = 00 10 = 00 20 = 00 y condiciones
iniciales cero es:
% archivo general
= 0;
= [0 : 005 : 1];
= [0 : 05 : 10];
% re-asignar variables
= 110; = 00; = 200; 2 = 00; 10 = 00; 20 = 00;
= [1; 1];
94
= [1 ∗ (1); 2 ∗ (2)];
= [1; 2];
0 = ()
0 = [00; 00];
10 = 00; 20 = 00; 1 = 50; 2 = 50; 1 = 1 1; 2 = 1 2;
= [−1; −2];
= [0; 0];
= [1; 1];
% finalización de programa
Para utilizar las librerías de Matlab de solución de ODEs en la solución de PDE se utilizan
esquemas de diferencias finitas. Para el ejemplo del reactor tubular con dispersión axial (2.50)-
(2.53), la discretización de las segundas derivadas por esquemas en diferencias finitas centradas y
95
Concentracion, x 1
0.5
0
10
0.8 1
5 0.6
0.2 0.4
0 0
Tiempo, t Posicion, x
20
Temperatura, x2
10
0
10
0.8 1
5 0.6
0.2 0.4
0 0
Tiempo, t Posicion, x
Figura 3-4: Simulación dinámica de un reactor con dispersión axial. Perfiles espacio-temporales.
96
de las primeras derivadas por esquemas de diferencias finitas retrasadas de las variables 1 y 2
están dadas por:
¯ ¯
11 − 10 ¯¯ 21 − 20 ¯¯
¯ = 1 (10 − 11 ) ¯ = 2 (20 − 21 ) (3.31)
∆ =0 ∆ =0
¯ ¯
1 +1 − 1 ¯¯ 2+1 − 2 ¯¯
¯ = 0 ¯ =0 (3.32)
∆ = ∆ =
= ( )
97
% definir los parámetros del modelo
= 1;
= ( ∗ 2 1);
= 2 : − 1
98
11 = (() − ( − 1));
1 = ();
= 50;
= 1 :
99
0() = 0; 0( + ) = 0;
figure(1)
figure(2)
Se considera el modelo del reactor lote biológico presentado en la Sección 2.5.2 que se describe
por las ecuaciones (2.29-2.30). Los parámetros del modelo son máx = 0227 = 604 = 998
100
1.5
Concentracion
1
x=L
0.5
x = L/2
x=0
0
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
Tiempo, t
15
Temperatura
10
0
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
Tiempo, t
= 78 × 10−3 = 028 = 01 y 2 = 38 Los archivos de la simulación dinámica
son:
= ( )
= ;
= ;
= (2 1);
101
1
Microorganismos, x1
0.8
0.6
0.4
0.2
0
0 10 20 30 40 50 60 70 80 90 100
Tiempo, t
3
Substrato, x2
0
0 10 20 30 40 50 60 70 80 90 100
Tiempo, t
La Figura (3-6) muestra la evolución dinámica de los microorganismos y del substrato dentro
del reactor. Como se puede observar, la degradación biológica del substrato (tiodiglicol) se alcanza
en un tiempo aproximado de 80 unidades de tiempo y debido a los efectos combinados de la
cinética tipo Andrews y la desactivación de microorganismos, la cultura de microorganismos (A.
xylosoxydans SH91) alcanza un crecimiento máximo alrededor de 50 unidades de tiempo.
Se considera el modelo del reactor continuo de tanque agitado con cinética de Van de Vusse
presentado en la Sección 2.5.4 que se describe por las ecuaciones (2.34-2.36). Los parámetros para la
simulación son = = 8256 1 = 51000 3 = 1049 10 = 1287 × 1012 20 = 9043 × 106
102
20 = 9043 × 106 ∆1 = 1 = 42 ∆2 = 2 = −110 ∆3 = 3 = −4185 =
= 86668 = = 308285 1 = = 01 1 = = 3522 × 10−4 1 =
97583 2 = 85600 3 = 85600 = −6239 1 = 10000 2 = 5000 3 = 900
y 4 = 1000
Los archivos de la simulación dinámica son:
= ( )
= 8256; 1 = 51000; 3 = 1049; 10 = 128712; 20 = 90436; 30 = 90439;
1 = 42; 2 = −110; 3 = −4185; = 86668; = 308285; = 01;
= 3522 − 4; 1 = 97583; 2 = 85600; 3 = 85600; = −6239;
= (4 1);
(2 1 1); ( (: 1) (: 2)) (0Concentraciones, x_1, x_20) (0Tiempo0)
(2 1 2); ( (: 3) (: 4)) (0Temperaturas, x_3, x_40) (0Tiempo0)
103
3000
Concentraciones x1, x2
x1
x2
2000
1000
0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Tiempo, t
120
Temperaturas, x3, x4
100
x3
80 x4
60
40
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Tiempo, t
Se considera el modelo del reactor continuo biológico de flujo tapón que se presentó en la
Sección 2.5.7 que se describe por las ecuaciones (2.46-2.49). Las PDE del modelo se discretización
en la posición espacial por el método de diferencias finitas, lo cual conduce a un conjunto de ODEs
dadas por:
máx
= − (3.34)
+
− −1
=− − 1 máx (3.35)
∆ +
104
Para obtener resultados de simulación con una precisión satisfactoria con diferencias finitas, podría
ser necesario un elevado número de puntos malla. Los valores numéricos de los parámetros que se
usan en la simulación son los reportados por Dochain et al., los cuales corresponden a un reactor
planta piloto anaerobio para tratamiento de aguas residuales: = 0022 1 = 04 = 005
−1 = 0002 3 máx = 035 −1 = 04 = 02489 = 3734 A
través de simulaciones numéricas se encontró que 100 nodos o puntos malla son suficientes para
reproducir los resultados obtenidos por Dochain et al. al utilizar colocación ortogonal.
Los archivos de la simulación dinámica son:
= ( )
= 1;
= ( ∗ 2 1);
% nodo 1
% nodo 2 - N
= 2 :
105
Archivo de solución de ecuaciones: ()
= 100;
= 1 :
0() = 41758;
0( + ) = 02784;
function dx = destil(t,x)
106
20
Microorganismos, X
15
x=L
10 x = L/2
0
0 50 100 150
Tiempo, t
4
Substrato, S
0
0 50 100 150
Tiempo, t
0.8
Composicion plato 1
0.6
R = 2.25
0.4
R = 3.25
R = 1.25
0.2
-0.2
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
Tiempo, t
-5
x 10
10
8
Composicion plato 20
0
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
Tiempo, t
107
= 200; = 400; = 50; 1 = 1575; = 225; = 246;
= 1 + (1 − ) ∗ ;
= ( + 1);
= ∗ ;
1 = + ∗ ;
= 1 − 1;
% relación de equilibrio
= 1 :
= ( 1);
= 2 : − 1
= + 1;
= : − 1
108
() = (1 ∗ ( − 1) − 1 ∗ () − ∗ ());
(2 1 1); ( (: 1)) (0Composicion plato 10) (0Tiempo, t0)
(2 1 2); ( (: 20)) (0Composicion plato 200) (0Tiempo, t0)
La Figura (3-9) muestra la evolución dinámica de la composición en el domo y en el fondo de
la columna para diferentes relaciones de reciclo. Se puede observar que al aumentar la relación de
reciclo se incrementa la composición de destilado que se obtiene. Esto se debe a que se logra un
mayor enriquecimiento del vapor que se condensa en el domo de la columna.
Bequette B.W. (2003). Process Control: Modelling, Design and Simulation. Upper Saddle
River, Prentice Hall.
Farlow, S.J. (1993). Partial Differential Equations for Scientists and Engineers. Dover Pu-
blications, Inc. New York.
Luyben W.L. (1991). Process. Modelling, Simulation and Control for Chemical Engineers.
McGraw Hill.
109
Capítulo 4
Análisis de Procesos
4.1. Introducción
El análisis de procesos tiene el objetivo de caracterizar las soluciones analíticas o numéricas de
los modelos matemáticos que describen un proceso.
Para modelos lineales, algunas herramientas para caracterizar las soluciones son:
Transformada de Laplace.
Trasformada de Fourier.
Para modelos no-lineales las herramientas son más complejas y se basan principalmente en el
uso de métodos numéricos:
Soluciones numéricas.
Planos de fase.
Linealización.
Diagramas de bifurcación.
Exponentes de Lyapunov.
Exponentes fractales.
110
Las soluciones analíticas de un modelo matemático proporcionan la máxima cantidad de infor-
mación acerca del proceso que se modela, sin embargo, solo son posibles para modelos matemáticos
lineales y pocos modelos no-lineales. Para modelos matemáticos no-lineales es posible obtener solu-
ciones numéricas para algún conjunto de condiciones, lo cual conduce a una menor información
del proceso.
La caracterización de las soluciones de un modelo que describe un proceso es útil en simulación
y control de procesos por dos razones principales:
1. Entender la operación del proceso: Para llevar a cabo una operación eficiente de un proceso
es primordial conocer los comportamientos posibles del proceso que se está analizando.
2. Diseñar en forma exitosa un sistema de control: Para llevar a cabo un diseño de control exi-
toso es indispensable entender la operación del sistema a ser controlado, así como determinar
el comportamiento dinámico del proceso a lazo cerrado.
•
= ( ) (4.1)
111
Las ecuaciones de diseño de sistemas de control se pueden representar de la forma:
•
= ( ) (4.2)
= ( )
•
= () (4.3)
Una solución () de la ecuación (4.2) usualmente corresponde a una curva en el espacio de
estados (el espacio vectorial que se forma de los estados del sistema) conforme va desde
cero a infinito. Esta curva se refiere generalmente como trayectoria del estado o trayectoria
del sistema.
112
Figura 4-1: Estabilidad e inestabilidad de un punto de equilibrio.
•
= () (4.4)
113
donde () es una matriz de × Un sistema es lineal si satisface dos condiciones:
Los sistemas invariantes en tiempo también se les llama autónomos y los sistemas variantes
en tiempo como no-autónomos.
114
4.4.1. Variables de desviación
En ingeniería de control, una operación normal del sistema puede ocurrir alrededor de un
punto de operación, de modo que es conveniente introducir variables de desviación respecto al
punto de operación. Así, si la condición de operación normal corresponde al punto (∗ ∗ ) las
variables de desviación son, ̄ = − ∗ y ̄ = − ∗ . El punto de operación normal (∗ ∗ ) puede
variar dependiendo de la aplicación específica. En el contexto de diseño u optimización de procesos
pueden ser las condiciones de operación promedio de operación o el valor de diseño nominal del
proceso. En el contexto de control de procesos se relaciona a los valores de operación deseados o
puntos de referencia. El uso de variables de desviación proporcionan las bases para la linealización
de modelos.
1. Se considera una función no-lineal, con entrada () y salida (). La relación entre () y
() se obtiene mediante:
= () (4.7)
¯ ¯
2 ¯2
¯ 1
() ≈ (∗ ) + ¯ ( − ∗ ) + ¯ ( − ∗ )2 + (4.8)
¯∗ 2! 2 ¯∗
donde ∗ = (∗ ).
115
4. La ecuación (4.9) se puede escribir en términos de las variables de desviación ̄ = − ∗ y
̄ = − ∗ como: ¯
¯¯
̄ = ̄ (4.10)
¯∗
lo cual indica que ̄ es proporcional a ̄.
La ecuación (4.10) es un modelo matemático lineal para la función no lineal (4.2) cerca del
punto de operación en términos de variables de desviación (̄ ̄)
La expansión en series de Taylor se puede generalizar para un sistema de funciones no-lineales
de la forma (4.2), por medio de la siguiente aproximación:
⎡ ⎤
1 1 1
⎢ 1 2 ⎥
⎢ 2 2 2⎥
⎢ 1 ⎥
=⎢
⎢
2 ⎥
⎥ (4.11)
⎢ ⎥
⎣ ⎦
1
∗ ∗
donde la matriz se conoce como el Jacobiano del sistema no-lineal (4.2). La aplicación del
Jacobiano al sistema (4.2) conduce a un sistema lineal de la forma:
•
= (4.12)
•
= ( ) (4.13)
= ( )
116
conduce a las matrices y :
⎡ ⎤ ⎡ ⎤
1 1 1 1 1 1
⎢ 1 2 ⎥
⎢ 1 2 ⎥
⎢ 2 2 2 ⎥ ⎢ 2 2 2⎥
⎢ 1 ⎥ ⎢ 1 ⎥
= ⎢
⎢
2 ⎥
⎥ =⎢
⎢
2 ⎥
⎥
⎢ ⎥ ⎢ ⎥
⎣ ⎦ ⎣ ⎦
1
1
⎡ ⎤ ⎡ ⎤
1 1 1 1 1 1
⎢ 1 2 ⎥ ⎢ 1 2 ⎥
⎢ 2 2 2 ⎥ ⎢ 2 2 2⎥
⎢ 1 ⎥ ⎢ 1 ⎥
= ⎢
⎢
2 ⎥
⎥ =⎢
⎢
2 ⎥
⎥
⎢ ⎥ ⎢ ⎥
⎣ ⎦ ⎣ ⎦
1
1
de modo que un modelo lineal para fines de control está dado por:
•
= + (4.14)
= +
En muchos casos, solo se manipula una entrada de control y la medición es comúnmente una
función lineal de los estados, además de que no se ve afectada por la entrada de control. Así, es
común que la matriz de la entrada de control sea un vector columna, la matriz de mediciones
sea un vector renglón y que la matriz de entrada de control afectando las mediciones sea cero.
Es conveniente clarificar dos puntos:
Los modelos lineales sirven para predecir el comportamiento local (o en una zona donde
es valida la linealización), de un sistema no-lineal. Esto es, los modelos linealizados no son
útiles para estudiar el comportamiento lejos del punto de linealización.
La dinámica que despliega un sistema no lineal es mucho más rica que la de un sistema lineal,
debido a la presencia de fenómenos no lineales tales como múltiples puntos de equilibrio,
ciclos límite, caos, etc.
117
4.4.4. Puntos de equilibrio y linealización con Matlab
Con Matlab es posible calcular los puntos de equilibrio de un sistema no-lineal (o las raíces de
un sistema de ecuaciones algebraicas no-lineales) y realizar la linealización de un sistema no-lineal.
Puntos de equilibrio
En el punto de equilibrio las derivadas respecto al tiempo son cero, de forma que las ecuaciones
de un modelo no-lineal se pueden expresar de la forma general:
2 (1 2 ) = 0
=
(1 2 ) = 0
donde () son funciones no-lineales de los estados Por ejemplo, para el modelo de Lotka-
Volterra descrito por la ecuación (3.20) se tiene:
1 (1 2 ) = 1 1 − 2 1 2 = 0 (4.16)
2 (1 2 ) = 3 1 2 − 4 2 = 0
y es fácil mostrar que el modelo de Lotka-Volterra tiene dos puntos de equilibrio: (i) la solución
trivial = (0 0) y (ii) = (4 1 2 )
Para encontrar la solución (o puntos de equilibrio), para un sistema de la forma (4.15) en
Matlab se utilizan las funciones solve y fsolve.
solve:
fsolve:
118
1. Escribir en un archivo las funciones no lineales a resolver en la forma,
= nombre()
solve:
1. Archivo ():
= ()
1 = 1; 2 = 2; 3 = 2; 4 = 1;
= [1 1 − 2 1 2 ; 3 1 2 − 4 2 ];
0 = [1; 1];
La solución que produce Matlab en el primer caso es 1 = (0 4 ) y 2 = (0 1 2 ) y en el
segundo caso es = (05 05)
119
Linealización
1 2
[1 (1 2 ); 2 (1 2 );
=
; (1 2 )];
= [1 2 ];
1 2 3 4 1 2
= [1 1 − 2 1 2 ; 3 1 2 − 4 2 ; ]
120
El resultado que se obtiene es:
⎡ ⎤ ⎡ ⎤
1 1
1 − 2 2 −2 1
=⎣ 1 2 ⎦ =⎣ ⎦ (4.18)
2 2
1 2
3 2 3 1 − 4
∗1 ∗2 ∗1 ∗2
Se considera el modelo del CSTR con cinética de Van de Vusse descrito por las ecuaciones
(2.34)-(2.36). Las funciones no-lineales de un modelo para fines de control, ( ) están dadas
por:
1 3 2
1 ( ) = (10 − 1 ) − 01 exp(− )1 − 03 exp(− ) (4.19)
3 3 1
1 2
2 ( ) = − 2 + 01 exp(− )1 − 02 exp(− )2
3 3
1 1 1 2
3 ( ) = (30 − 3 ) − ∆1 01 exp(− )1 − ∆2 02 exp(− )2 −
3 3
1 3 2
∆3 03 exp(− )1 + (4 − 3 )
3
1
4 ( ) = (3 − 4 ) +
donde = ( ) son los estados del CSTR y = es la variable manipulable o entrada
de control. En CSTRs no-isotérmicos es común elegir el medio de enfriamiento como la variable a
manipular, con la finalidad de regular la temperatura dentro del reactor, lo cual en forma indirecta
permite regular la concentración de salida en el reactor.
Con los parámetros de la simulación dinámica de la Sección 3.6.2, el punto de equilibrio que
se obtiene a través de Matlab es, ∗ = (1810 10591 1106 1106).
La linealización del CSTR está dada por:
⎡ ⎤ ⎡ ⎤
1 1 1 1 1
⎢ 1 2 3 4 ⎥ ⎢ ⎥
⎢ 2 2 2 2 ⎥ ⎢ 2 ⎥
⎢ ⎥ ⎢ ⎥
=⎢
⎢
1 2 3 4 ⎥
⎥ =⎢
⎢
⎥
⎥
3 3 3 3 3
⎢ ⎥ ⎢ ⎥
⎣ 1 2 3 4 ⎦ ⎣ ⎦
4 4 4 4 4
1 2 3 4 ∗1 ∗2 ∗3 ∗4 ∗1 ∗2 ∗3 ∗4
121
lo cual conduce a:
⎡ ⎤
1 1
0 0
⎢ 1 3 ⎥
⎢ 2 ⎥
⎢ 01 exp(−
1
) − − 02 exp(− 2
) 0 ⎥
=⎢
⎢
3 3 3 ⎥
⎥ (4.20)
3 3
⎢ − 1 ∆2 02 exp(−
2
) ⎥
⎣ 1 3 3 ⎦
0 0
−
∗1 ∗2 ∗3 ∗4
⎡ ⎤
0
⎢ ⎥
⎢ ⎥
⎢ 0 ⎥
=⎢
⎢
⎥
⎥ (4.21)
⎢ 0 ⎥
⎣ ⎦
01
∗1 ∗2 ∗3 ∗4
donde:
substituyendo los parámetros de la simulación dinámica del CSTR y eligiendo el punto de linea-
lización como el punto de equilibrio, se obtiene el modelo lineal:
⎡ ⎤
−8256 0 ∼0 0
⎢ ⎥
⎢ ⎥
⎢ ∼0 −8 256 ∼0 0 ⎥
=⎢
⎢
⎥
⎥ (4.23)
⎢ ∼0 ∼0 −39 085 308285 ⎥
⎣ ⎦
0 0 86668 −86668
122
donde los términos menores a 1 × 10−20 se han tomado como cero.
4.5. Estabilidad
La estabilidad es un concepto clave en la operación y control de un proceso. La razón principal
es que los resultados de una operación inestable son impredecibles, con posibles daños a equipo
de proceso e incluso afectando el bienestar de los trabajadores. Para sistemas con control es
importante evaluar la estabilidad por las siguientes razones:
Los sistemas de control deben de ser estables bajo todas las condiciones de operación posibles.
2. Estabilidad local. Las trayectorias que inician dentro de cierta región del espacio de estados
se mantienen dentro de esta región para todo tiempo futuro.
3. Estabilidad global. Las trayectorias que inician en cualquier región del espacio de estados se
mantienen dentro del espacio de estados para todo tiempo futuro.
La Figura (4-2) ilustra los tres conceptos de estabilidad arriba descritos. Las flechas continuas
son trayectorias estables, las flechas punteadas indican trayectoria inestables, la curva cerrada con
una línea punteada es una región local del sistema y las curvas cerradas es el espacio de estados
donde las variables del sistema existen.
Para sistemas lineales la estabilidad se puede determinar de los valores propios de la matriz
del sistema (4.12). En particular, se tienen los siguientes resultados:
123
Figura 4-2: Estabilidad asintótica, local y global.
se dice que es estable, con respecto al punto de equilibrio ∗ si todos los valores propios de
la matriz tienen parte real no positiva. El sistema es estable asintóticamente, si todos los
valores propios de tienen parte real negativa. En el último caso, la matriz se denomina
Hurwitz.
∗ = −−1 (4.26)
•
si k() − ∗ k → 0 conforme → ∞ Si la matriz es Hurwitz, el sistema lineal = +
es estable asintóticamente.
124
Si la matriz es Hurwitz, existen dos números positivos y Λ tal que:
Es decir, las trayectorias del sistema lineal (4.27) permanecen acotadas conforme → ∞.
La propuesta más útil para estudiar la estabilidad de un sistema no-lineal se debe al trabajo de
Lyapunov, el cual incluye dos métodos para el análisis de la estabilidad: el método de linealización
y el método directo.
Si el sistema linealizado es estrictamente estable (i.e., si todos los valores propios de tienen
parte real negativa), entonces el punto de equilibrio del sistema no-lineal es asintóticamente
estable.
125
Si el sistema linealizado es inestable (i.e., si al menos un valor propio de tiene parte real
positiva), entonces el punto de equilibrio del sistema no-lineal es inestable.
•
= ( ) (4.30)
= ( )
El sistema (4.30) es estable si la respuesta del sistema esta acotada en magnitud cuando el
sistema se sujeta a una entrada de control acotada. Un sistema que exhibe una respuesta
no acotada a una entrada de control acotada es inestable. Una función de entrada acotada es
una función del tiempo que siempre se mantiene dentro de ciertos límites conforme transcurre
el tiempo. Por ejemplo, la función sinusoidal y la función escalón son entradas acotadas.
Un sistema con tal propiedad se dice que es de entrada acotada-salida acotada o BIBO. Aun
más, se pude probar que un sistema que es estable asintóticamente es además BIBO estable.
•
= (4.31)
•
= + (4.32)
126
la estabilidad se puede obtener con Matlab del cálculo de valores propios de la matriz :
⎡ ⎤
1
⎢ 11 12 ⎥
⎢ ⎥
⎢ 21 22 2 ⎥
=⎢
⎢
⎥
⎥ (4.33)
⎢ ⎥
⎣ ⎦
1 2
de la siguiente forma:
= [11 12 1 ; 21 22 2 ; ; 1 2 ; ]
3. Si todos los valores propios de la matriz son negativos, el sistema (4.31) o (4.32) es estable,
de otra forma es inestable.
1. Componentes de la matriz :
= [−8256 0 0 0; 0 − 8256 0 0;
2. Instrucción: ()
3. Se obtiene como resultado = (−8256 −119779 −5974 − 8256) por lo tanto el modelo
linealizado del CSTR, en el punto de equilibrio ∗ es estable.
127
4.6. Clasificación de puntos de equilibrio
Para clasificar los puntos de equilibrio de un sistema no-lineal, es conveniente estudiar los
puntos de equilibrio de un sistema de segundo orden de la forma:
•
1 = 1 (1 2 ) (4.34)
•
2 = 2 (1 2 )
al linealizar el modelo (4.34) en el punto de equilibrio (∗1 ∗2 ) se obtiene el modelo lineal de
segundo orden: ⎡ ⎤
1 1
•
=⎣ 1 2 ⎦ = (4.35)
2 2
1 2 ∗1 ∗2
Los puntos de equilibrio del sistema no-lineal (4.34), se pueden clasificar por el comportamiento
dinámico de las trayectorias del sistema lineal de segundo orden (4.35), el cual esta dictado por
los valores propios de la matriz
La clasificación de los puntos de equilibrio de de acuerdo a los valores propios, 1 2 es:
Nodo estable, si ambos valores propios son reales y negativos, 1 2 0 Por ejemplo el
sistema:
•
1 = −21 − 32 (4.36)
•
2 = −1 − 22
Nodo inestable, si ambos valores propios son reales y positivos, 1 2 0 Por ejemplo
el sistema:
•
1 = 21 − 32 (4.37)
•
2 = −1 + 22
128
Silla, si ambos valores propios son reales y un valor propio es negativo y otro positivo,
1 0 2 . Por ejemplo el sistema:
•
1 = 2 (4.38)
•
2 = 1 + 2
Foco estable, si los valores propios son complejos, = ± y la parte real es negativa,
1 2 0 Por ejemplo el sistema:
•
1 = −1 − 2 (4.39)
•
2 = 1 − 2
Foco inestable, si los valores propios son complejos y la parte real es positiva, 1 2 0
Por ejemplo el sistema:
•
1 = 1 − 2 (4.40)
•
2 = 1 + 2
Centro, si los valores propios son complejos y la parte real es cero, 1 = 2 = 0 Por
ejemplo el sistema:
•
1 = −2 (4.41)
•
2 = 1
129
4.7. Plano de fase y diagrama de bifurcación
Con la finalidad de visualizar las propiedades dinámicas de los puntos de equilibrio y de un
sistema, es útil introducir el plano de fase y el diagrama de bifurcación.
La solución de ecuaciones diferenciales puede ser descrita, en forma algebraica como una función
del tiempo, o en forma gráfica por una curva, que indica la evolución en el tiempo de dos estados del
sistema a partir de sus condiciones iniciales. A excepción de algunos casos especiales, las soluciones
algebraicas no se pueden derivar en términos de funciones simples, de modo que la descripción
gráfica es una mejor opción. Estas gráficas se conocen como planos de fase.
Los planos de fase se dibujan en el espacio fase, el cual es un sistema coordenado con un eje
para cada variable (para sistemas de dimensión mayor a tres se grafican en combinaciones de
solo dos estados). A partir de un conjunto de condiciones iniciales, las soluciones evolucionan a
lo largo de alguna curva en el espacio fase. La dinámica que ocurre, antes de que las trayectorias
se establezcan en un estado estacionario, se denomina la dinámica transitoria de la trayectoria.
El plano de fase se puede dividir en zonas en las cuales se despliegan diferentes propiedades. Las
fronteras de estas zonas se denominan separatrices.
El concepto de estabilidad en el plano de fase se deduce del comportamiento de las trayectorias
y como convergen a un punto de equilibrio. Por ejemplo, si la trayectoria se establece en algún
punto fijo se tiene un comportamiento estable. En el comportamiento inestable las trayectorias
divergen al acercarse al punto de equilibrio.
La Figura (4-3) muestra las trayectorias dinámicas de los sistemas (4.36)-(4.41) para diferentes
condiciones iniciales. En los casos de sistemas estables (foco y nodo estables), las trayectorias
convergen desde las diferentes condiciones iniciales hacia el punto de equilibrio del sistema. Para
sistemas inestables (foco inestable, nodo inestable y silla), las trayectorias divergen desde las
condiciones iniciales hacia valores al infinito. El caso del centro, es un sistema que se mantiene
oscilando alrededor de una trayectoria que se genera desde las condiciones iniciales establecidas.
130
1 3
0.5 1
0 -1
-2
-0.5 -3
-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 -4 -3 -2 -1 0 1 2 3 4
1 5
0.4 2
0.2 1
0 0
-0.2 -1
-0.4 -2
-0.6 -3
-0.8 -4
-1 -5
-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 -5 -4 -3 -2 -1 0 1 2 3 4 5
1.5
5
1
3
2
0.5
1
dx1/dt = - x2
0 dx2/dt = x1 0
-1
-0.5
-2
dx1/dt = x2
-3 dx2/dt = x1 + x2
-1
-4
-1.5
-1.5 -1 -0.5 0 0.5 1 1.5 -5
-5 -4 -3 -2 -1 0 1 2 3 4 5
131
4.7.2. Diagramas de bifurcación
Cuando un sistema, lineal o no lineal, involucra uno o varios parámetros, es posible que al
variar el valor de alguno de los parámetros cambie la naturaleza cualitativa de los puntos de equi-
librio y la forma en que convergen las trayectorias del sistema. Por ejemplo, conforme cambia un
parámetro, los estados estacionarios pueden llegar a ser inestables y conducir a oscilaciones esta-
bles, o un sistema con un estado estacionario estable puede ser reemplazado por sistemas con varios
estados estacionarios estables. El punto donde se produce un cambio en la naturaleza cualitativa
de la solución, al cambiar un parámetro del sistema, se conoce como punto de bifurcación. La
construcción de una gráfica que muestra el cambio de la solución a largo plazo de alguna variable,
al cambiar algún parámetro del sistema, se conoce como diagrama de bifurcación.
132
abordado en la literatura y esta asociada a reacciones exotérmicas. El caso más común de múltiples
estados estacionarios, que se refiere como bi-estabilidad, es el de dos estados estacionarios estables
separados por un estado inestable.
El ejemplo clásico, para introducir el problema de múltiples estados estacionarios en procesos
químicos, es un CSTR no-isotérmico, donde se lleva a cabo la reacción exotérmica →
propuesto originalmente por Aris y Amudson en 1958. El modelo matemático del CSTR de Aris
y Amudson se describe por las ecuaciones:
= (0 − ) − 0 exp(− ) (4.42)
∆
= (0 − ) + 0 exp(− ) + ( − )
1
= (10 − 1 ) − 1 exp(−2 ) (4.43)
2
= (20 − 2 ) + 1 exp(−2 ) + (2 − 2 )
Para construir el diagrama de Van Heerden se consideran las ecuaciones del CSTR en estado
estacionario:
133
300
250
Calor generado/removido
200
150
removido
100 generado
50
0
350 375 400 425 450 475 500
Temperatura, x2
10 exp(−2 )
(20 − 2 ) + (2 − 2 ) = − (4.46)
+ exp(−2 )
esta ecuación se puede interpretar físicamente de acuerdo a Van Heerden (1958). El lado izquierdo
de la ecuación es la velocidad total de remoción de calor y el lado derecho de la ecuación es
la velocidad de generación de calor, las cuales en estado estacionario deben de ser iguales. Los
estados estacionarios se determinan al graficar las velocidades del calor generado y calor removido
como una función de la temperatura del reactor. La intersección de las curvas de calor generado y
calor removido representan los puntos de operación. La Figura (4-4) muestra el diagrama de Van
Herdeen para el CSTR.
La característica principal en reactores, que da lugar a la multiplicidad de los estados esta-
cionarios, es la forma sigmoidal de la curva de generación de calor contra la temperatura y debido
a que la transferencia de calor depende linealmente de la temperatura, las curvas de generación
y remoción de calor pueden dar lugar a más de una intersección. Para reacciones exotérmicas,
134
conforme la reacción procede, los reactantes se agotan, lo cual tiende a causar un descenso en la
velocidad de reacción. Al mismo tiempo, el aumento de calor incrementa la temperatura y causa
que la velocidad de reacción se incremente, a través de la dependencia de Arrhenius de la veloci-
dad de reacción sobre la temperatura. Estos dos efectos de la velocidad de reacción, asociados con
la conversión de los reactantes, conducen a una dependencia no monotónica de la velocidad de
reacción sobre la concentración de reactantes, conduciendo a la posibilidad de múltiples estados
estacionarios.
Al analizar la Figura (4-4) se observa que para reacciones irreversibles se pueden presentar tres
casos:
2. El calor desprendido por la reacción es más que el necesario, el fluido se calienta y la con-
versión será prácticamente completa.
3. Un caso intermedio con tres soluciones (dos estables, una inestable) para las ecuaciones de
balance de materia y energía.
Un punto de operación es estable si la remoción de calor es mayor que calor generado por
la reacción. El reactor retornara inherentemente a su temperatura original si ocurre una pertur-
bación. Un punto inestable se presenta cuando un ligero incremento en la temperatura del reactor,
resulta en una velocidad de generación de calor mayor que la capacidad de remoción de calor. Así,
la temperatura en el reactor aumentara a un punto de operación estable de alta temperatura y
conversión. De igual forma, una ligera disminución en la temperatura reduce la velocidad de reac-
ción, tal que la chaqueta enfría el reactor hasta un punto de operación estable de baja temperatura
y conversión. Se observa que si la pendiente de la curva de producción de calor es mayor que la de
la curva de remoción de calor, entonces el estado estacionario es inestable. Si la pendiente de la
curva de calor generado es menor de la curva de remoción de calor, entonces el estado estacionario
es estable.
135
650
600
550
x2
500
450
400
350
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
x1
La Figura (4-5) muestra el plano de fase del CSTR para diferentes condiciones iniciales. Depen-
diendo de la condición inicial, el CSTR converge a estados estacionarios estables, que corresponden
a conversiones y temperaturas bajas y altas. El punto de operación inestable no se alcanza para
cualquier condición inicial, debido a que cualquier perturbación lleva el sistema a un punto de
operación estable.
Se puede observar que, de acuerdo a la forma de convergencia de las trayectorias a los puntos
de equilibrio, el punto de equilibrio de baja conversión es del tipo nodo estable y el punto de
equilibrio de alta conversión es del tipo foco estable.
136
480
460
440
Temperatura, x2
420
400
380
360
340
320
310 320 330 340 350 360 370 380
Temperaturadeenfriamiento, x2w
Para sistemas no-lineales que presentan estados estacionarios inestables el sistema puede evolu-
cionar a la aparición de oscilaciones sostenidas alrededor del estado estacionario inestable. En el
espacio fase definido por las variables de estado, las oscilaciones sostenidas corresponden a la
137
evolución hacia un ciclo límite (centro en sistemas lineales). Tales oscilaciones se caracterizan por
su amplitud y su periodo. Otro tipo de comportamiento oscilatorio que puede surgir es el caos,
el cual ha sido estudiado intensivamente en sistemas físicos, químicos y biológicos. En el espacio
fase, las oscilaciones caóticas corresponden a la evolución hacia un denominado atractor extraño.
Estas oscilaciones irregulares se caracterizan por su sensibilidad a las condiciones iniciales, lo cual
toma en cuenta la naturaleza impredecible de la dinámica caótica.
Para ilustrar las oscilaciones en procesos químicos se considera un sistema de reacciones auto-
catalíticas que se llevan a cabo en forma isotérmica en un CSTR:
+ 2 →1 3 = −1 2
→2 = −2 (4.47)
+ 2 →3 3 = −3 2
= (0 − ) − 1 2
= (0 − ) − 2 + 1 2 + 3 2 (4.48)
= (0 − ) − 3 2
1
= 1 − 1 − 1 1 22
2
= 1 − (1 + 2 )3 + 1 1 1 22 + 2 3 3 22 (4.49)
3
= 1 − 3 − 3 3 22
Las Figuras (4-7) y (4-8) muestran comportamientos oscilatorios periódicos y caóticos, así como el
plano de fase correspondiente, que se generan con el modelo (4.49) y los parámetros 1 = 180000
2 = 800 3 = 4000 1 = 15 con dos valores de la relación de especies de a en la
138
0.08
0.07
0.06
0.05
0.04
x1
0.03
0.02
0.01
0
0 0.5 1 1.5 2 2.5 3
Tiempo, t
0.1
0.08
0.06
x1
0.04
0.02
0
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
x3
Los reactores tubulares de flujo tapón no exhiben fenómenos de múltiples estados estaciona-
rios, oscilaciones y formación de patrones, debido en gran medida a la ausencia de un efecto de
retroalimentación en el sistema. Sin embargo, los reactores tubulares de flujo tapón presentan
dos efectos dinámicos interesantes: (i) la sensibilidad paramétrica y (ii) la presencia de puntos
calientes.
139
0.08
0.06
0.04
x1
0.02
0
0 0.5 1 1.5 2 2.5 3
0.1
0.08
0.06
x1
0.04
0.02
0
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
x3
140
a azul, asociados con la reducción-oxidación de la especies reactantes. La clase de reacciones que
producen este tipo de patrones de coloración se conocen como reacciones de Belousov-Zhabotinskii.
Turing formuló el primer mecanismo de formación de patrones en 1952, en el cual propuso
que el acoplamiento de la difusión con la cinética no-lineal es suficiente para producir soluciones
estacionarias que conducen a no homogeneidades espacio-temporales. El mecanismo propuesto por
Turing fue capaz de explicar y reproducir la formación de patrones en química y biología.
Las reacciones químicas no-lineales y la difusión no son los únicos mecanismos que conducen
a la formación de patrones espacio-temporales en un sistema químico. Un grupo de investigación
en Bélgica, conocido como el grupo de Bruselas, realizó estudios relacionados a la formación de
patrones al estudiar el origen químico de tales inestabilidades y concluyeron que el acoplamiento
entre cinéticas no-lineales y un mecanismo retroalimentado, son las causas de inestabilidades
estacionarias. El mecanismo retroalimentado puede ser térmico, como en el caso de reacciones
químicas, o en el caso de reacciones autocatalíticas un retroalimentado químico. En años recientes
también se ha reportado la formación de patrones debido al acoplamiento entre los fenómenos de
convección, difusión y reacción. Algunos reactores que han reportado la presencia de formación
de patrones son los reactores tubulares con reciclo, autotérmicos, operados con inversión de flujo
periódica y forzados.
El reactor tubular con dispersión axial que se presentó en la Sección 2.5.8 presenta la formación
de un patrón espacio-temporal oscilatorio para el siguiente conjunto de parámetros: 1 = 50
2 = 50 = 1 = 0185 = 140 = 30 = 200 2 = 00 10 = 00 20 = 00
La Figura (4-9) muestra la dependencia en el espacio y en el tiempo de las variables del reactor,
la concentración y la temperatura. La presencia de oscilaciones se muestra en los perfiles respecto
al tiempo.
141
1
Concentracion, x 1
0.5
0
10 9 8 7 1
6 0.8 0.9
5 4 0.6 0.7
3 0.4 0.5
2 1 0.2 0.3
0 0 0.1
Tiempo,t Posicion,x
10
Temperatura, x2
0
10 9 8 7 1
6 0.8 0.9
5 4 0.6 0.7
3 0.4 0.5
2 1 0.2 0.3
0 0 0.1
Posicion,x
Tiempot
142
4.9.1. Transformada de Laplace
Función escalón: La función escalón ideal, de magnitud es una función cuyo valor per-
manece cero hasta que este toma el valor de instantáneamente en el tiempo de inicio
establecido arbitrariamente en = 0
⎧
⎨ 0 0
() = (4.50)
⎩ ≥0
() = (4.52)
Función pulso rectangular ideal: Esta es una función que permanece cero hasta que toma
el valor de instantáneamente en el tiempo de inicio establecido arbitrariamente en = 0
143
mantiene este nuevo valor por unidades de tiempo y retorna a su valor inicial.
⎧
⎪
⎪ 0 0
⎪
⎨
() = 0 (4.53)
⎪
⎪
⎪
⎩ 0
La función pulso rectangular está dada también por, () = [() − ( − )] La trans-
formada de Laplace está dada por:
() = (1 − exp(−)) (4.54)
Función impulso ideal: La función impulso ideal de magnitud es una función que se
representa matemáticamente como:
() = (4.57)
Por otro lado, se tiene que si () representa la función escalón unitario, entonces:
() = () (4.58)
Función rampa ideal: Es una función cuyo valor inicia en = 0 y se incrementa linealmente
144
con una pendiente constante la cual se designa aquí como :
⎧
⎨ 0 0
() = (4.59)
⎩ 0
() = (4.60)
2
() = (4.62)
2 + 2
Se considera una ecuación de primer orden sujeta a un término externo o de forzamiento ():
•
() + () = () (4.63)
()( + ) = 0 + ()
0
() = + () (4.64)
+ +
145
con la condición inicial cero, la ecuación transformada es:
() = () (4.65)
+
•• •
() + 2() + 2 () = () (4.66)
() = () (4.68)
2 + 2 + 2
Para obtener la transformada directa e inversa de Laplace con Matlab se procede como sigue:
= ()
= ()
146
para obtener la solución evaluada en algún valor de o
Se consideran dos funciones simples: (i) () = exp() + 3 sin() y (ii) () = 1( − )
2. Escribir la instrucción: (i) = (); (ii) = ( 1) Lo que conduce a los
resultados:
•
= + (4.69)
= +
(0 ) = 0 (4.70)
donde () es el estado, () es la entrada, () es la salida, (0 ) es la condición inicial. y
son matrices de dimensiones adecuadas y coeficientes constantes.
Al tomar la transformada de Laplace para el sistema (4.69) con condiciones iniciales cero 0 = 0
se tiene:
147
que se puede re-escribir de la forma:
= ()()
la ecuación (4.73) se conoce como la función de transferencia del sistema (4.69), debido a que
transfiere una entrada () a una salida () y se puede considerar como un modelo entrada-
salida del sistema (4.69). La función de transferencia () permite una caracterización completa
entrada-salida de un sistema.
La función de transferencia () se puede escribir como () = ( − )−1 debido a que
el término de entrada de control comúnmente no está presente en la medición de un sistema. La
función de transferencia también se puede escribir como la relación de dos polinomios en términos
de la variable compleja de la transformada de Laplace:
()
() = ( − )−1 = (4.74)
()
−1
+ −1 + + 1 +
= −1
(4.75)
+ −1 + + 1 +
Los polinomios del numerador y denominador de la ecuación (4.75) se pueden factorizar en tér-
minos de sus raíces de acuerdo a:
donde 1 a son las raíces del numerador y 1 a son las raíces del denominador. De la
expresión (4.75) se observa que el orden del denominador es y el orden del numerador es
El orden relativo del sistema es la diferencia entre el orden del denominador y el orden del
148
numerador Para sistemas causales siempre y se conocen como sistemas estrictamente
propios. Si ≥ el sistema es propio y si el sistema es impropio.
Los ceros de un sistema lineal son las raíces del polinomio del numerador de la función de
transferencia, i.e., las raíces 1 a En tanto que los polos del sistema lineal son las raíces del
polinomio del denominador de la función de transferencia, i.e., las raíces 1 a Los ceros de la
función de transferencia no tienen efecto sobre la estabilidad del sistema. Por otra parte, los polos
determinan la estabilidad del sistema. Si todos los polos son negativos el sistema es estable, de
otra forma el sistema es inestable.
Las funciones de transferencia que no tienen polos ni ceros positivos o cero se conocen como
funciones de transferencia de fase mínima, en tanto que las que tienen polos y/o ceros negativos
son funciones de transferencia de fase no mínima. Los sistemas con funciones de transferencia
de fase mínima se denominan sistemas de fase mínima, en tanto que aquellos con funciones de
transferencia con fase no mínima se denominan sistemas de fase no mínima. Los sistemas de fase
no mínima son lentos en su respuesta, debido a un comportamiento defectuoso al inicio de la
respuesta.
Con Matlab es posible calcular: (i) la función de transferencia a partir de los polos y ceros de
un sistema, (ii) la función de transferencia a partir de un sistema lineal de la forma (4.69) y (iii)
obtener un sistema lineal de la forma (4.69) a partir de la función de transferencia.
1. Escribir las raíces del numerador (ceros) y del denominador (polos) de la función de
transferencia de la forma:
149
1. Escribir los componentes de las matrices y (la cual generalmente es cero):
= [11 12 1 ; 21 22 2 ; ; 1 2 ; ]
= [11 12 1 ; 21 22 2 ; ; 1 2 ; ]
= [11 12 1 ; 21 22 2 ; ; 1 2 ; ]
= [11 12 1 ; 21 22 2 ; ; 1 2 ; ]
= ( )
El modelo lineal del CSTR de Aris y Amudson (4.42) en el punto de equilibrio ∗ = ∗1 = 05
y ∗ = ∗2 = 400 considerando como entrada de control la temperatura de la chaqueta de
enfriamiento, = y como variable medible la temperatura del reactor, = = 2 está
dado por:
⎡ ⎤ ⎡ ⎤
• −2 −003125 0
= ⎣ ⎦ +⎣ ⎦ (4.77)
200 425 1
05400
= [0 1]
150
1. Escribir los componentes de las matrices y :
= [0; 1]
= [0 1]
= [0]
+2
= (4.78)
2 − 225 − 225
La función de transferencia (4.78) se puede escribir en términos de los ceros y polos del sistema
como:
( + 2)
= (4.79)
( + 075)( − 30)
de (4.79) se puede observar que el CSTR tiene un cero en −2 un polo negativo en −075 y un
polo positivo en 30 por lo que se puede concluir que el CSTR en el punto (05 400) es de fase no
mínima e inestable.
151
Perturbación impulso: El impulso es una idealización de una perturbación repentina de corta
duración. Este tipo de modelo puede representar perturbaciones de carga instantánea, así
como errores de medición.
Perturbación escalón: El modelo escalón se usa para representar una perturbación de carga
sostenida o una desviación en una medición.
Perturbación rampa: La rampa es una señal que es cero para un tiempo negativo y se
incrementa linealmente para el tiempo positivo. La rampa es útil para representar errores
de medición que se incrementan y perturbaciones que inician repentinamente sin parar. En
la práctica, las perturbaciones están acotadas, sin embargo la rampa es una idealización útil
en algunos casos.
Perturbación sinusoidal: La onda seno es un prototipo para una perturbación periódica. Con
una elección adecuada de la frecuencia, es posible representar perturbaciones de carga que
presentan fluctuaciones de baja frecuencia, así como ruido de medición de alta frecuencia.
Se considera un sistema lineal de primer orden forzado por una perturbación escalón, que por
conveniencia, se escribe de la forma:
•
() + −1
() = () (4.80)
donde es una constante de tiempo del sistema y es un parámetro que se conoce como la
ganancia en estado estacionario del proceso. Un ejemplo de un proceso de ingeniería es el de
un sistema de un tanque continuo agitado isotérmico, en el cual entra una corriente con una
concentración 0 = 1 opera en estado estacionario y en cierto instante de tiempo cambia a un
valor 0 = 2 El modelo matemático del tanque agitado es:
•
() = − () + () (4.81)
La Figura (4-10) muestra la respuesta entrada-salida para los valores = 4 mı́n, = 36
1 = 10 y diferentes valores del cambio escalón 2 en = 1000 mı́n
152
2
1.8
1.6
1.4
1.2 u= 1.5
u = 1.0 u= 2.0
y = C(t)
1 u= 0.5
u= 0.0
0.8
0.6
0.4
0.2
El sistema lineal de primer orden responde de una manera suave y con crecimiento exponencial
antes de llegar al nuevo valor en estado estacionario. Para perturbaciones positivas el nuevo valor en
estado estacionario es mayor y para perturbaciones negativas el nuevo estado estacionario es menor,
sin embargo, la magnitud de las perturbaciones son simétricas respecto al estado estacionario
inicial. De la Figura (4-10) se puede ver que el tiempo en el cual se alcanza el nuevo estado
estacionario es prácticamente el mismo sin importar el orden de cambio de la perturbación externa.
La solución en el domino del tiempo permite encontrar un parámetro que caracteriza la dinámi-
ca de sistemas de primer orden, la constante de tiempo del proceso, el cual no se modifica al
cambiar la magnitud del valor de la perturbación escalón. Por ejemplo, la solución al modelo (4.81)
está dada por:
() = (1 − exp(− )) (4.82)
donde en este caso = = y corresponde a un tiempo característico del proceso, que en
este caso es el tiempo de residencia promedio. Cuando = la solución del sistema (4.82) se
puede escribir como:
() = (1 − exp(−1)) = 0632 (4.83)
153
de modo que el tiempo característico es el tiempo en el cual se alcanza el valor del 632 %
del nuevo valor en estado estacionario cuando el sistema se somete a una perturbación escalón.
Un sistema forzado de segundo orden se describe por un ecuación diferencial de segundo orden
con el término de forzamiento :
2 −1
+ 2 + −2
= (4.84)
2
()
() = = 2 2 (4.85)
() + 2 + 1
la ecuación característica del sistema (4.85) está dada por la ecuación de segundo orden en el
denominador de (4.85):
2 2 + 2 + 1 (4.86)
() = =
( 2 + 2 + 1) ( − 1 )( − 2 )
de modo que la respuesta del sistema de segundo orden depende de la ubicación de las raíces 1
y 2 . Así, se pueden distinguir tres casos: (i) respuesta sobre-amortiguada cuando 1 (raíces
154
0.7
0.6
0.5
0.5
1.0
1.5
0.4 0.2
0.8
y
0.3
0.2
0.1
0 10 20 30 40 50 60 70 80 90 100
Tiempo
reales y distintas), (ii) respuesta críticamente amortiguada cuando = 1 (raíces iguales), (iii)
respuesta sub-amortiguada cuando 1 (raíces conjugadas).
De los tres casos, la respuesta sub-amortiguada, 1 es la que se presenta con mayor
frecuencia en los sistemas de control. La Figura (4-11) muestra la respuesta entrada-salida para
valores de = 05, = 436 y diferentes valores del factor de amortiguamiento
Para valores de menores de 10 (sistema sub-amortiguado), el sistema responde con sobredis-
paros y oscilaciones antes de llegar a un nuevo valor en estado estacionario. Para valores de = 10
(sistema críticamente amortiguado), la respuesta es suave y de tipo exponencial. Para valores de
mayores a 10 (sistema sobre-amortiguado), la respuesta sigue una forma sigmoidal antes de llegar
a un nuevo valor en estado estacionario. Se puede observar en la Figura (4-11) que la respuesta
de sistemas de segundo orden sobre-amortiguados y críticamente amortiguados es muy similar a
la respuesta de sistemas de primer orden. Para sistemas de segundo orden sub-amortiguados, los
parámetros que caracterizan la respuesta dinámica son tres:
2. El sobredisparo, el cual se define como la diferencia entre el valor más alto que se alcanza
155
0.7
0.6
b
0.5 sobredisparo
0.4
y
0.3
0.2
tiempo de asentamiento
0.1
tiempo de elevación
0
3. El tiempo de elevación, que es el valor en el cual se alcanza por primera vez el nuevo
valor en estado estacionario.
La Figura (4-12) muestra los tres parámetros para el sistema (4.84) con = 02
•
1 () = − 1 () + () (4.88)
•
2 () = − 2 () + 1
•
3 () = − 3 () + 2
•
4 () = − 4 () + 3
156
1
0.9
0.8
n=5
0.7
n=4
0.6 n=3
n=2
y = Ci
0.5 n=1
0.4
0.3
0.2
0.1
la relación entrada-salida se obtiene al derivar la ecuación del último tanque hasta que aparezca
la entrada aplicada al sistema de tanques. Después de una serie de manipulaciones algebraicas se
obtiene:
•••• ••• •• •
4 () + 4 + ( )2 (3 + 2 + 1 ) = () (4.89)
que es un sistema de cuarto orden.
La Figura (4-13) muestra la respuesta entrada-salida para un sistema de 5 tanques en serie
de igual volumen y valores = 4 mı́n, = 36 1 = 10 y 2 = 10 en = 500
mı́n La respuesta entrada-salida para cada uno de los tanques sigue una forma sigmoidal. En
el último tanque la forma sigmoidal es más pronunciada que la respuesta entrada-salida de los
tanques anteriores. En muchos casos, la respuesta entrada-salida de un sistema de orden se
puede aproximar con la respuesta entrada-salida de sistemas de primer y segundo orden.
157
4.12. Identificación empírica
Los modelos matemáticos que resultan de emplear los principios fundamentales, por lo general
tienen la gran limitación que da lugar a modelos muy complejos. En general, la construcción de
modelos matemáticos de la mayoría de procesos reales requiere un gran esfuerzo de ingeniería para
formular las ecuaciones, determinar todos los valores de los parámetros y resolver las ecuaciones,
usualmente a través de métodos numéricos. Este esfuerzo se justifica cuando son necesarias predic-
ciones muy precisas de la respuesta dinámica sobre un amplio rango de condiciones de operación.
Para fines de control, un modelo simple puede simplificar el diseño, simulación, o análisis de los
esquemas de control. Un método alternativo para construir modelos matemáticos diseñados especí-
ficamente para el control de procesos es la identificación empírica. Los modelos que se desarrollan
en este método proporcionan la relación dinámica entre variables de salida y entrada.
El modelado o identificación empírica consiste en desarrollar modelos, a partir de datos de
respuesta, a una perturbación de entrada. Los modelos se determinan realizando cambios pequeños
en las variables de entrada alrededor de las condiciones de operación nominales. La respuesta
dinámica resultante se usa para determinar el modelo. Este procedimiento es esencialmente una
linealización experimental del proceso, que es válida para alguna región alrededor de las condiciones
nominales. Los dos métodos de identificación más conocidos son: (i) el método de curva de reacción
y (ii) los métodos en base a principios estadísticos.
3. Obtener los datos de la respuesta de salida hasta que el proceso logre un nuevo estado
estacionario.
158
2
u = 2.0
1.8
1.6
y 2 - y 1 = 2.0 - 1.0 = 1.0
1.4
1.2
u = 1.0
1
y
0.8
0.6
0.4
0.2 t = 10.0
Un modelo de primer orden en el domino del tiempo y de Laplace se puede escribir como:
()
+ = = (4.90)
() + 1
∆
= = (4.91)
∆
159
La Figura (4-14) muestra los datos que se requieren para el ejemplo del tanque agitado (4.81). En
este caso:
2−1
= = 1 = 10 (4.92)
2−1
y el modelo empírico de primer orden es:
() 1
10 + = 1 = (4.93)
() 10 + 1
2 ()
2 2
+ 2 + = = 2 2 (4.94)
() + 2 + 1
en este caso la ganancia del proceso , la constante de tiempo del proceso y el factor de
amortiguamiento , definen el comportamiento dinámico del proceso de segundo orden.
De las respuestas de sistemas de segundo orden a una entrada escalos que se obtuvieron en la
Sección 4.11, se puede ver que en los casos de ≥ 1 las características dinámicas se pueden ajustar
a un modelo de primer orden. Para el caso sub-amortiguado, 1 el ajuste de la respuesta
de un modelo de segundo orden sub-amortiguado requiere los siguientes datos: (i) El incremento
de la perturbación escalón que se aplica ∆ = 2 − 1 (ii) El incremento del valor en estado
estacionario ∆ = 2 − 1 (iii) La amplitud del sobredisparo, 1 (iv) La amplitud de otro pico
antes de que el sistema se establezca en el nuevo valor en estado estacionario, . (v) El tiempo
entre dos picos sucesivos, . Los cálculos del modelo empírico de segundo orden sub-amortiguado
son:
∆
=
∆
1
−1
log( 1 )
= q
(4.95)
4 + [ −1 log( 1 )]2
2 1
q
= 1 − 2
2
La Figura (4-15) muestra los datos para el cálculo de los parámetros para el ejemplo (4.84)
160
0.7
0.6
b1 = 0.23 b2 = 0 .0 5
0.5
u2 = 1.0
0.4
y
0.3
0.1
u1 = 0.0
0
045 − 00
= = 045
10 − 00
1
2−1
log( 023
005
)
= q = 0236
1
4 2 + [ 2−1 log( 023
005
)]2
125 p
= 1 − 02362 = 19332
2
2
1932 + 2(0236)(193) + = 045 (4.96)
2
o en el dominio de Laplace:
() 045
= 2 2
(4.97)
() 193 + 2(0236)(193) + 1
161
Modelo de primer orden con tiempo muerto
()
= exp(− ) (4.98)
() + 1
donde () denota la salida, () la entrada, la ganancia en estado estacionario, la constante
de tiempo del proceso y el tiempo muerto. El único dato adicional al modelo de primer orden es
el tiempo muerto, el cual se obtiene de manera directa de la respuesta entrada-salida del sistema,
como el tiempo que tarda el sistema en responder con un cambio en la variable de salida.
Algunos sistemas muestran una respuesta inversa al aplicarse una perturbación de entrada, la
cual se debe a efectos competitivos entres diferentes dinámicas del proceso. La respuesta inversa
se puede considerar como un tiempo muerto, al tomar el tiempo en el cual el sistema regresa a
una respuesta directa y dividir este valor entre dos, es decir:
= (4.99)
2
162
0.9
0.88 X: 8.599
X: 5.434 Y: 0.8668
Composicion plato 1
Y: 0.8601
0.86 X: 4.246
Y: 0.8487
0.84
X: 5.441
Y: 0.8343
0.82 X: 8.636
Y: 0.826
0.8
1 2 3 4 5 6 7 8 9 10
Tiempo, t
-5
x 10
4.1
Composicion plato 20
4 X: 8.826
Y: 3.935e-005
X: 5.228
X: 4.231
Y: 3.886e-005
3.9 Y: 3.87e-005
3.8
X: 5.285
Y: 3.822e-005 X: 8.847
Y: 3.791e-005
3.7
3.6
0 1 2 3 4 5 6 7 8 9 10
Tiempo, t
163
gativas son:
8668 − 8487
+1 = = 8 044 4 × 10−2 +
1 = 0434 (4.100)
225 ∗ 11 − 225
826 − 8487
−1 = = 0100 89 −1 = 0441
225 ∗ 09 − 225
8 044 4 × 10−2
()+
1 = (4.102)
0434 + 1
0100 89
()−
1 =
0441 + 1
2 888 9 × 10−6
()+
20 = exp(−02) (4.103)
0228 + 1
3 511 1 × 10−6
()−
20 = exp(−02)
0285 + 1
En este caso, a pesar de que el sistema responde diferente a perturbaciones positivas y negativas,
la diferencia no es significante y se puede establecer un modelo de respuesta escalón con valores
promedio de ambas perturbaciones. De esta forma se obtienen los modelos:
009
()1 = (4.104)
0438 + 1
32 × 10−6
()20 = exp(−02)
025 + 1
164
4.12.3. Identificación empírica de un reactor tubular
Los reactores tubulares son otra clase de procesos químicos para los cuales es común traba-
jar con modelos lineales entrada-salida, que se obtienen por el método de curva de reacción. Se
considera el modelo del reactor biológico presentado en la sección 2.5.7 cuando se somete a una
perturbación escalón de ± 10 % en el flujo de alimentación, = después de que el sistema
alcanza el estado estacionario. La respuesta de concentración de substrato a la salida del reactor
se muestra en la Figura (4-17).
Los cálculos del modelo lineal son:
03393 − 02599
+ = = 3970 +
= 52 (4.105)
0002 ∗ 11 − 0002
01881 − 02599
− = = 3590 −
= 57
0002 ∗ 09 − 0002
3780
() = (4.106)
55 + 1
que es un modelo lineal que establece la relación entrada-salida del reactor tubular biológico de
flujo tapón.
A partir de un modelo lineal de la forma (4.69) con Matlab se puede obtener la respuesta a un
escalón por medio del siguiente procedimiento.
165
0.4
X: 245
0.35
Y: 0.3393
X: 156.2
Y: 0.3215
0.3
Substrato, S
X: 128.9
Y: 0.2599
0.25
0.15
0 50 100 150 200 250
Tiempo, t
Step Response
1
0.9
0.8
0.7
0.6
Amplitude
0.5
0.4
0.3
0.2
0.1
0
0 10 20 30 40 50 60
Time (sec)
166
Por ejemplo, si se considera el sistema (4.88), el modelo en el espacio de estados es:
⎡ ⎤ ⎡ ⎤
−436 0 0 0 436
⎢ ⎥ ⎢ ⎥
⎢ ⎥ ⎢ ⎥
• ⎢ 436 −436 0 0 ⎥ ⎢ 0 ⎥
= ⎢
⎢
⎥ + ⎢
⎥ ⎢
⎥
⎥ (4.107)
⎢ 0 436 −436 0 ⎥ ⎢ 0 ⎥
⎣ ⎦ ⎣ ⎦
0 0 436 −436 0
= [1 0 0 0]
= [1 0 0 0]
= [0]
La Figura (4.27) muestra la gráfica que se genera con las instrucciones anteriores. Se puede ver
de la Figura (4.27) que la respuesta a un cambio escalón unitario es equivalente a los resultados
que se muestran en la Figura (4-13) que se realizaron a través de la programación y simulación
numérica del modelo (4.88).
167
Smith C.A. y Corripio, A.B. (2001). Control Automático de Procesos, Teoría y Práctica.
Limusa.
Bequette B.W. (2003). Process Control: Modelling, Design and Simulation. Upper Saddle
River, Prentice Hall.
Luyben W.L. (1991). Process. Modelling, Simulation and Control for Chemical Engineers.
McGraw Hill.
Froment, G.F., Bischoff, K.B. (1990). Chemical Reactor Analysis and Design. John Wiley
and Sons, Inc.
Elnashaie, S.S.E.H., Elshishini, S.S. (1993). Modelling, Simulation and Optimization of In-
dustrial Fixed Bed Catalytic Reactors. Topics in Chemical Engineering. Edited by R. Hughes.
Gordon and Breach Science Publishers. OPA.
Slotine, J.-J. E., Li, W. (1991) Applied nonlinear control. Prentice Hall.
168
Capítulo 5
5.1.1. Introducción
Existen muchos controladores diferentes, pero el controlador retroalimentado que se usa más
ampliamente es el algoritmo proporcional integral y derivativo (PID). Su popularidad se debe a
las siguientes razones:
2. Simplicidad de su estructura.
4. No es costosa su implementación.
169
Existen muchas variaciones, pero la definición del algoritmo PID es:
Z
() = () + () + [()] (5.1)
0
= (5.2)
=
Los tres modos principales del control son el modo proporcional ( ), el modo integral
() y el modo derivativo ().
170
Sin embargo, la respuesta del modo integral es lenta y por lo tanto se usa comúnmente
combinada con otros modos.
En la práctica industrial es común combinar los tres modos. La selección de los mejores modos
del controlador depende además fuertemente de las características del proceso. En general, los
métodos de simulación se usan para probar diferentes estrategias y modos de control.
Algunas recomendaciones de combinación de los modos del control PID son:
Control P: En muchos casos, se puede considerar que el mejor control es el más simple, que
en el control PID corresponde al control puramente proporcional. El control proporcional
se puede recomendar cuando es aceptable tener una desviación en estado estacionario. En
control de procesos el control proporcional se usa con frecuencia para el control de nivel y
de presión.
Control PID: La combinación de los tres modos permite que la variable controlada llegue a
la referencia deseada en forma más rápida y con menos sobredisparos que el control PI. La
171
desventaja del control PID es su sensibilidad al ruido de medición, de modo que el control
PID se recomienda para procesos con respuesta lenta y mediciones poco susceptibles al ruido
de medición.
5.1.3. Sintonizado
El controlador PID se puede ajustar al cambiar los valores de la ganancia y las constantes
de tiempo de integración y tiempo derivativo El objetivo del sintonizado del controlador
es elegir los parámetros del controlador en forma adecuada para obtener un desempeño deseado.
Esto usualmente significa que las variables controladas deben de retornar al valor deseado a pesar
de perturbaciones en el sistema y seguir cambios en la señal de referencia en forma aceptable.
El controlador se puede ajustar por prueba y error, ya sea a través de experimentación o por
simulaciones numéricas. El método de prueba y error consiste en los siguientes pasos:
3. Agregar la acción integral con valores altos de Reducir los valores de en factores de 2
hasta que la respuesta sea oscilatoria, donde = 0 .
172
Controlador
P 1( )
PI 09( ) 33
PID 12( ) 20 2
Cuadro 5.1: Parámetros del sintonizado de Ziegler Nichols
Controlador
P
(1 + 3 )
PI
(09 + 12 ) 30+3
9+20
PID
(4
3
+
12
) 32+6
13+8
4
12+2
Método de Ziegler-Nichols
Este método se basa en la identificación empírica a lazo abierto por el método de curva de
reacción. Los datos que se requieren de la curva de reacción son el incremento de la variable
controlada y la entrada de control, ∆ y ∆ Los parámetros del método de Ziegler-Nichols son la
pendiente a través del punto de inflexión, normalizada por la magnitud ∆ es decir = ∆
y por su intersección con el eje (tiempo), denotado como . Las relaciones de sintonizado
completas se resumen en la Tabla 5.1.
Método de Cohen-Coon
El método de Cohen y Coon se basa también en datos del método de identificación empírica
a lazo abierto. Los valores teóricos de los parámetros del controlador se resumen en la Tabla 5.2.
Los parámetros y están dados por = ∆∆ y = ∆
La ganancia ultima, 0 es la ganancia que lleva el sistema a oscilaciones sostenidas solo
con acción proporcional. La frecuencia de las oscilaciones en la ganancia ultima se denotan como
0 donde 10 se denomina periodo ultimo. El periodo ultimo se determina experimentalmente
al incrementar desde valores pequeños hasta el valor en el cual surgen las oscilaciones. Los
parámetros del controlador se calculan de 0 y 0 de acuerdo a las reglas dadas en la Tabla 5.3.
173
Controlador
P 050
PI 0450 1120
PID 060 120 180
Cuadro 5.3: Parámetros del sintonizado del método de ganancia ultima
Control proporcional:
= ();
= ;
= − ;
= ;
174
( = )
= ∗ ;
= ;
Control proporcional-integral.
= ;
() = ;
() = ; =
= ( + 1 1);
( + 1) = ;
175
en el paso 1 y 2.
( = )
= ∗ + ∗ ( + 1);
= ;
5. Para graficar los valores que toma la entrada de control, agregar al final del archivo
que contiene las ecuaciones del modelo la instrucción,
esta instrucción genera una curva punteada en una de las figuras que se generan en el
archivo que contienen la instrucción de integración de ecuaciones.
6. Agregar una condición inicial con valor de cero en el archivo que contiene la instruc-
ción de integración.
En el control PID se tiene que tener cuidado al implementar numéricamente la acción deriva-
tiva. Esto debido a que la derivada del error de regulación es sensible a ruido de medición en la
variable controlada y puede conducir a cambios importantes y frecuentes de la variable manipula-
da. Aun más, la acción derivativa ideal no es realizable debido a que es impropia. Una alternativa
es usar filtros de primer orden para implementar una versión filtrada de la acción derivativa:
≈
+
176
simout
To Workspace
Clock
simout1
1
PID To Workspace1
s+1
PID Controller Transfer Fcn Step
simout2
To Workspace2
Subtract
simout3
To Workspace3
177
Paso (a)
simout
Clock To Workspace
Paso (g)
simout2 simout3 Paso (h)
To Workspace2 To Workspace3
Paso (c)
Paso (b)
1
PID simout1
Step s+1
Subtract PID Controller Transfer Fcn To Workspace1
Paso (d) Paso (e)
Paso (f)
d) Conectar la salida del módulo de comparación a la entrada del módulo del controlador
PID (Substract +
− → PID Controller).
e) Conectar la salida del controlador PID a la entrada del módulo de la función de trans-
ferencia (PID Controller → Transfer Fcn).
178
h) Conectar un módulo de datos de entrada al espacio de trabajo a la señal de salida del
controlador PID (PID Controller → To Workspace 3 ).
d) Cambiar en el módulo Step el tiempo de entrada (step time) y los valores iniciales
(initial value) y finales (final value) del cambio de escalón, que corresponden a las
referencia deseadas,
g) Cambiar en el módulo Transfer Fcn los valores de los polinomios del numerador y
denominador de la función de transferencia (ver Figura (5-3)).
h) Cambiar en el módulo PID Controller los valores de las ganancias proporcional, integral
y derivativa (ver Figura (5-3)).
Se considera el modelo del reactor continuo de tanque agitado de Aris y Amudson presentado
en la Sección 4.8.1 que se describe por las ecuaciones (2.34-2.36). El objetivo de control es controlar
la temperatura del reactor, = = 2 en el punto inestable del CSTR, 2 = 4000 La variable
manipulable es la temperatura de la chaqueta de enfriamiento, = = 2 La entrada de control
se activa en el tiempo de 100 minutos.
Los archivos de la simulación dinámica y la implementación del control PI son:
179
Figura 5-3: Asignación de parámetros de módulos en Simulink.
180
400
Entrada de control, u
350
300
Kc = 0.5
250 Kc = 5.0
Kc = 25.0
200
0 10 20 30 40 50 60 70 80 90 100
Tiempo, t
1200
1000
Temperatura, x2
800
600
400
200
0
0 10 20 30 40 50 60 70 80 90 100
Tiempo, t
= ( )
= 250; = 10; = 100000; = 10; = 2000; = 10; = 3500;
= 3500;
% control PI
= (2);
= 40000;
= − ;
= 250;
= 05;
( = 50)
181
= ∗ + ∗ (3);
= ;
= (3 1);
(3) = ;
182
650
600 Ki = 0.1
550 Ki = 0.25
Temperatura, x2 Ki = 0.5
500
450
400
350
0 50 100 150
Tiempo, t
1400
1200
Entrada de control, u
1000
800
600
400
200
0
0 50 100 150
Tiempo, t
control para un valor de la ganancia proporcional de = 25 y tres de la ganancia integral. Los
valores de las ganancias proporcionales e integrales se han asignado por un método de prueba y
error simple. Se observa que en los tres casos, el controlador lleva la temperatura del reactor a la
referencia deseada. Para valores moderados de la ganancia integral y proporcional se obtiene el
menor esfuerzo de control y el mejor desempeño a lazo cerrado. Por una parte valores muy bajos
de la ganancia integral vuelven muy lenta la respuesta del controlador y por otro lado, valores
muy grandes de la ganancia integral hacen que el esfuerzo inicial de control sea muy grande.
−023 −023
() = = 2
(5.4)
( + 04)( + 03) 012 + 07 + 1
183
t
Clo ck
T iem po
E n tra d a d e c o n tro l
-0 .2 3
y
0 .1 2 s2 + 0 .7 s+ 1
P e rtu rb a c i o n S alida m edida
Fu n cio n d e
e sc a l o n
tra n sfe re n c i a
la cual tiene dos polos ubicados en −04 y −03 y puede representar modelos entrada-salida con
respuesta de segundo orden, tales como columnas de destilación, intercambiadores de calor y
reactores químicos. La Figura (5-6) muestra el diagrama de flujo en Simulink para simular una
perturbación escalón, que se aplica en 5 unidades de tiempo con una magnitud de 0.1. La Figura
(5-7) muestra la respuesta escalón de la función de transferencia de segundo orden (5.4).
En la Figura (5-8) se muestra el diagrama de flujo de la implementación de un controlador
PID en Simulink para la función de transferencia (5.4). El objetivo de control es regular la salida
de la función de transferencia (5.4) a una referencia de 0 en las primeras 5 unidades de tiempo y
posteriormente regular a una referencia de 0.1. La Figura (5-9) muestra el desempeño a lazo cerrado
del control proporcional, proporcional-integral y proporcional-integral-derivativo. Los parámetros
del controlador se obtienen a través del método de sintonizado de Ziegler-Nichols. La Figura (5-8)
muestra los datos necesarios para el sintonizado de Ziegler-Nichols, de donde se obtiene:
−0012 + 00075
= = −00225
56 − 54
−00225
= = −0225
01
el sintonizado de Ziegler-Nichols se obtiene de los cálculos del Cuadro 5.1. Por ejemplo, para el
184
0
t = 0.1
m = (y 2 - y 1)/(x 2 - x1)
-0.005 = (-0.12 + 0.075)/(5.6 - 5.4))
= -0.0225
-0.015
-0.02
-0.025
5 5.2 5.4 5.6 5.8 6 6.2 6.4 6.6 6.8 7 7.2 7.4 7.6 7.8 8
Tiempo, t
Figura 5-7: Parámetros de sintonizado del método de Ziegler-Nichols para una función de trans-
ferencia de segundo orden.
12 12
= = = −5333
(01)(−0225)
−5333
= = = = −26665 (5.5)
2 (2)(01)
01
= = = −5333 = −266
2 2
Como se puede ver en la Figura (5-9) el controlador P presenta una desviación en estado esta-
cionario importante y los controladores PI y PID permiten regular la variables de salida al valor de
referencia deseado. Debido a la acción derivativa el controlador proporcional responde más rápido,
que el control PI, para llevar al sistema a la referencia deseada.
185
t
Clock
Tiempo
u
yref Entrada de control
set point
-0.23
PID y
Perturbacion 0.12s2 +0.7s+1
Suma PID Controller Salida medida
escalon Funcion de
transferencia
0.16
0.14
0.12
0.1 Kp = -44.44
Kp = -40.00, Ki = -121.21
0.08 Kp = -53.33, Ki = -266.66, Kd = -2.66
y ref
0.06
0.04
0.02
0 2 4 6 8 10 12
Figura 5-9: Desempeño del control PID para un sistema de segundo orden.
186
(a )
y ref G c(s ) G v(s ) G p(s)
u y
c o n tro la d o r a c tu a d o r p ro c e so
G m (s)
sensor
(b )
y ref e (s) y (s)
G (s)
y H (s)
H (s)
Figura 5-10: Diagramas de bloques de un esquema de control retroalimentado: (a) Completo. (b)
Simplificado.
(5-10-b). () es el producto de las funciones de transferencia de la ruta directa, () =
() () () que recibe como entrada la señal de referencia () () es el producto de la
ruta retroalimentada, que recibe como entrada la señal de salida (). La función de transferencia
que relaciona () con () se denomina la función de transferencia a lazo cerrado.
A partir del diagrama de bloques (5-10-b) se pueden establecer las siguientes relaciones:
187
al sustituir las ecuaciones (5.6) y (5.7) en (5.8) se obtiene:
()
() = () (5.9)
[1 + ()()]
que indica que la función de transferencia a lazo cerrado se puede expresar como la relación de la
función de transferencia de la ruta directa entre 1 más el producto de las funciones de transferencia
de la ruta directa y la ruta retroalimentada.
() = (5.10)
0 + 1
y se consideran dos variaciones del control PID: (i) acción proporcional pura y (ii) acción proporcional-
integral.
se considera además que los elementos de medición y actuación tienen velocidades de respuesta
muy rápidas, en comparación a la dinámica del sistema (5.10), de modo que las funciones de
transferencia respectivas se pueden despreciar.
Después de algunas manipulaciones algebraicas, las funciones de transferencia a lazo cerrado
188
para los controladores proporcional y proporcional-integral están dadas por:
() () = () (5.13)
0 + 1 +
+
() () = () (5.14)
0 2 + (1 + ) +
el valor de () en estado estacionario se puede obtener al usar el teorema del valor final:
para un cambio escalón unitario en la señal de referencia, () = 1 la respuesta en estado
estacionario es:
() () = (5.15)
1 +
() () = =1 (5.16)
de la ecuación (5.15) se puede observar que para tener un error en estado estacionario cero,
() = 1 se requieren valores muy altos de la ganancia proporcional del controlador . Por otra
parte, de la ecuación (5.16) se observa que la acción integral elimina la desviación en estado
estacionario.
Los números complejos se pueden representar en forma gráfica en forma análoga a la re-
presentación gráfica de números reales. Esta representación gráfica de los números complejos se
denomina el plano complejo. En el plano complejo la coordenada de un punto o número complejo
= Re() + Im() son sus partes real e imaginaria. Si = + entonces se localiza a
unidades a la derecha del origen y unidades por encima del origen (si y/o son negativos,
entonces se localiza a la izquierda y/o por abajo del origen). El eje horizontal o eje se denomina
el eje real y el eje vertical o eje se denomina eje imaginario.
189
5.4.2. Polos a lazo cerrado
donde son los polos a lazo cerrado, los cuales son además las raíces de la ecuación característica
y son los ceros a lazo cerrado.
La posición de los polos de la función de transferencia en el plano imaginario determina la
naturaleza del comportamiento dinámico del sistema a lazo cerrado. La Figura (5-11) muestra
las propiedades de estabilidad de un sistema dependiendo la ubicación de sus polos en el plano
complejo.
De la Figura (5-11) se puede observar lo siguiente:
Si los polos se ubican en el semiplano izquierdo, es decir valores negativos de la parte real,
el sistema a lazo cerrado es estable.
Si los polos se ubican en el semiplano derecho, es decir valores positivos de la parte real, el
sistema a lazo cerrado es inestable.
Para polos con parte real más alejada del origen el sistema reacciona más rápido hacia la
estabilidad (polos negativos) o inestabilidad (polos positivos).
Para polos con parte imaginaria más alejada del origen la respuesta estable (polos negativos)
o inestable (polos positivos) se vuelve más oscilatoria.
190
Im
11
1
Foco Foco 1.5
1
x 10
0.5
estable inestable 0.5
0
0
-0.5
-0.5
1.5
Centro -1
-1 -1.5
0 2 4 6 8 10 1 0 5 10 15 20 25 30
0.5
Nodo -0.5
Nodo
-1
- estable -1.5
inestable +
0 5 10 15 20 25 30
0 13
Re
1 x 10
6
0.8
4
Silla
0.6
0.4 2
0.2 20
x 10
0
3 0
-0.2 2 -2
-0.4
-0.6
1 -4
-0.8 0
0 2 4 6 8 10 -6
0 2 4 6 8 10
-1
-2
-3
0 5 10 15 20 25 30
191
El módulo o valor absoluto de un número complejo es la distancia desde el origen en el plano
complejo. Asi,
p
|| = 2 + 2 (5.18)
en donde los coeficientes son cantidades reales, i.e., 6= 0 Es decir, se elimina cualquier
raíz cero.
3. Si todos los coeficientes son positivos, se ordenan los coeficientes del polinomio en renglones
192
y columnas de acuerdo con el patrón o arreglo siguiente:
0 2 4 6
−1 1 3 5 7
−2 1 2 3 4
−3 1 2 3 4
−4 1 2 3 4 (5.21)
2 1 2
1 1
0 1
1 2 − 0 3 1 4 − 0 5 1 6 − 0 7
1 = 2 = 1 =
1 1 1
La evaluación de las continúa hasta que todos los resultados son cero. Se sigue el mismo
patrón de multiplicación en forma cruzada los coeficientes de los dos renglones anteriores al
evaluar las , las , las , etc. Es decir:
1 3 − 1 2 1 5 − 1 3 1 7 − 1 4
1 = 2 = 3 =
1 2 3
4. El criterio de estabilidad de Routh plantea que el número de raíces de la ecuación (5.20) con
partes reales positivas es igual al número de cambios de signo de los coeficientes de la primera
columna del arreglo (5.21). Debe señalarse que no es necesario conocer los valores exactos de
los términos de la primera columna, ya que sólo se necesitan los signos. La condición necesaria
y suficiente para que todas las raíces de la ecuación (5.20) se encuentren en el semiplano
izquierdo del plano es que todos los coeficientes de la ecuación (5.20) sean positivos y que
todos los términos de la primera columna del arreglo (5.21) tengan signo positivo.
193
de control lineal, sobre todo porque no sugiere cómo mejorar la estabilidad relativa ni cómo
estabilizar un sistema inestable. Sin embargo, es posible determinar los efectos de cambiar uno o
dos parámetros de un sistema si se examinan los valores que producen inestabilidad.
2 1 2 − 225
1 − 225 0 (5.24)
0 2 − 2 25
de modo que para que el CSTR sea estable ( todas las raíces de la ecuación característica se
encuentran en el semiplano izquierdo del plano complejo) todos los términos de la primer columna
de (5.24) deben de ser positivos. Por lo tanto:
225 (5.25)
W. R. Evans diseñó un método sencillo para encontrar las raíces de una ecuación característica.
Este método se denomina método del lugar geométrico de las raíces y en él se grafican las raíces de
la ecuación característica para todos los valores de un parámetro del sistema. El lugar geométrico de
las raíces o diagrama de localización de raíces es la gráfica de las raíces de la ecuación característica
de un sistema, con o sin control, en el plano complejo, en función de un parámetro variable (una
ganancia, un cero del controlador, etc). Mediante el método del lugar geométrico de las raíces, el
ingeniero de control puede predecir los efectos que tiene el introducir un controlador con modo
194
proporcional en un sistema a lazo abierto.
Para la función de transferencia del CSTR (4.78) a lazo abierto, la Figura (5-12) muestra el
diagrama de localización de raíces.
La Figura (5-12) muestra el efecto de la ganancia de un controlador proporcional sobre los polos
del sistema a lazo cerrado. Para una = 00 los polos del sistema se ubican en 30 y en −075 y
se tiene un sistema inestable. Para = 225 los polos del sistema se ubican en ∼ 00 ± 15 y se
tiene un sistema críticamente estable, lo cual concuerda con el criterio de Routh presentado en la
sección anterior. Para = 112 los polos del sistema se ubican en −45± ∼ 00 Para ∼ ∞
los polos se ubican en −20± ∼ 00.
()
() = = (5.26)
() + 1
195
Root Locus
2.5
0.96 0.92 System: h0.85 0.7 0.45
Gain: 10.1 System: h
Pole: -3.93 + 1.59i Gain: 2.25
0.978
2 Damping: 0.927 Pole: -0.00622 + 1.49i
Overshoot (%): 0.0422 Damping: 0.00417
Frequency (rad/sec): 4.24 Overshoot (%): 98.7
Frequency (rad/sec): 1.49
1.5
0.991
System: h System: h
0.998 Gain: 11.2 Gain: 0
0.5 Pole: -4.5 - 0.000201i Pole: -0.75
Damping: 1 Damping: 1
Overshoot (%): 0 Overshoot (%): 0
Imaginary Axis
0.991
-1.5 System: h
Gain: 10.4 System: h
Pole: -4.08 - 1.38i Gain: 2.25
Damping: 0.947 Pole: -0.00806 - 1.5i
-2 Overshoot (%): 0.0094 Damping: 0.00539
0.978
Frequency (rad/sec): 4.31 Overshoot (%): 98.3
Frequency (rad/sec): 1.5
0.96 0.92 0.85 0.7 0.45
-2.5
-10 -8 -6 -4 -2 0 2 4
Real Axis
196
con transformada de Laplace:
1
() = (5.28)
2 + 2
donde 1 y son la amplitud y frecuencia de la señal sinusoidal, se obtiene:
1
() = (5.29)
( + 1)(2 + 2 )
1 1
() = exp(− ) + √ sin( + ) (5.30)
1 + 2 2 1 + 2 2
donde:
= tan−1 (−) = − tan−1 () (5.31)
se puede observar en la ecuación (5.30) que a largo plazo (i.e., → ∞) el término exponencial se
desvanece, exp(− ) → 0 de modo que se obtiene:
1
() = √ sin( + ) (5.32)
1 + 2 2
1
2 = √ (5.33)
1 + 2 2
y con un ángulo + donde es un ángulo que esta definido por la ecuación (5.31). Con base a
la respuesta en frecuencia de un sistema de primer orden es conveniente introducir las siguientes
relaciones:
197
Ángulo de fase: Es la cantidad en radianes o grados en la cual la señal de salida se atrasa
(con signo negativo) o adelanta (con signo positivo) respecto a la señal de entrada.
= √ (5.36)
1 + 2 2
1
= √ (5.37)
1 + 2 2
= tan−1 (− ) (5.38)
y se puede observar que las relaciones anteriores son una función de la frecuencia de entrada
El análisis de respuesta a la frecuencia consiste esencialmente en determinar la relación de
amplitud, la relación de magnitud y el ángulo de fase de un sistema que se somete a diferentes
frecuencias, con la finalidad de utilizar esta información en el análisis y síntesis de sistemas con
control.
El análisis de respuesta en frecuencia se puede realizar a través de dos métodos: El método ex-
perimental y por medio de la transformación de la función de transferencia del proceso a controlar
cuando se aplica una señal sinusoidal. En el primero el sistema se somete a entradas sinusoidales
con diferentes frecuencias a través de un generador de frecuencia variable y se registra la respuesta
del sistema. En el segundo método se utilizan las propiedades matemáticas de los números com-
plejos, de modo que con la finalidad de obtener la relación de magnitud y el ángulo de fase se
substituye por la variable compleja en la función de transferencia del proceso a controlar y
se calcula la magnitud y el argumento del número complejo que resulta. La magnitud es igual a
la relación de amplitud y el argumento es igual al ángulo de fase.
Para ilustrar el segundo método se considera nuevamente el sistema de primer orden (5.26), al
substituir por se obtiene:
() = (5.39)
+ 1
la relación de amplitud es:
| |
= |()| = =√ (5.40)
| + 1| + 1
198
el ángulo de fase está dado por el ángulo de la expresión del número complejo:
= tan−1 ( )
= = 20 log10
199
Bode Diagram
20
Magnitude (dB)
-20
-40
-60
-80
0
-45
Phase (deg)
-90
-135
-180
-2 -1 0 1 2
10 10 10 10 10
Frequency (rad/sec)
Nyquist Diagram
1.5
2 dB 0 dB -2 dB
-4 dB
1
4 dB
-6 dB
6 dB
0.5
10 dB -10 dB
Imaginary Axis
20 dB -20 dB
0
-0.5
-1
-1.5
-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1
Real A xis
200
(0 0) es la relación de amplitud y se define como:
p
= |( 1 )| = Re (1 )2 + Im ( 1 )2 (5.43)
por otro lado, el ángulo que se forma entre el punto con frecuencia 1 y el eje real es el ángulo de
fase:
Im ( 1 )
= tan−1 (5.44)
Re (1 )
La Figura (5-14) muestra el diagramas de Nyquist de la función de transferencia (5.42).
El análisis de la respuesta en frecuencia es una herramienta útil para diseñar sistemas de control
retroalimentados. En particular, el análisis de la respuesta en frecuencia se utiliza en sistemas de
control para:
Estudiar las características de estabilidad de un sistema con control al usar los diagramas
de Bode y/o de Nyquist de la función de transferencia a lazo abierto, que está dada por:
Si el sistema es estable a lazo abierto, para que el lazo cerrado sea internamente estable es
necesario y suficiente que el diagrama de Nyquist de () () no encierre al punto (−1 0)
Si el sistema es inestable a lazo abierto, con polos en el semiplano derecho, entonces para
que el lazo cerrado sea internamente estable es necesario y suficiente que el diagrama de
Nyquist de () () encierre veces en sentido anti-horario al punto (−1 0)
201
Si el diagrama de Nyquist de () () pasa por el punto (−1 0) el sistema a lazo cerrado
tiene polos exactamente sobre el eje imaginario. Esta situación se conoce como condición de
estabilidad crítica.
si la ecuación característica tiene alguna raíz positiva el sistema a lazo cerrado es inestable, lo cual
se cumple si el diagrama de Nyquist de la función de transferencia (5.46) encierra el punto (−1 0)
conforme la frecuencia toma valores desde 0 a ∞
El criterio de estabilidad de Nyquist se puede extrapolar al diagrama de Bode, donde el punto
(−1 0) corresponde a un ángulo de fase de −180◦ y una magnitud de 1 o equivalentemente un
módulo de 0 . Los resultados básicos del criterio de estabilidad de Bode son:
El término incertidumbre se refiere a las diferencias o errores entre los modelos y la realidad.
Las incertidumbres de un modelo pueden ser paramétrica o estructural.
202
operación conduciendo a cambios en los parámetros del modelo y conocimiento imperfecto
de los parámetros del modelo.
Los modelos de procesos que se usan comúnmente en el control de procesos son modelos lineales
invariantes en tiempo, los cuales representan la dinámica verdadera del proceso sólo en forma
aproximada, de modo que es importante poder cuantificar el impacto que tendrá la incertidumbre
de un modelo en el diseño de un sistema de control.
Los diagramas de Bode y de Nyquist permiten cuantificar la robustez de un lazo de control.
La robustez se refiere a la capacidad de un controlador a tolerar incertidumbres de un modelo
matemático de un proceso. Entonces, la estabilidad robusta se refiere a la capacidad que tiene un
controlador para mantener la estabilidad del lazo al aplicarse a la planta verdadera.
Dos criterios que se usan para determinar cuantitativamente las capacidades de estabilidad
robusta de un controlador son el margen de ganancia y el margen de fase.
El margen de ganancia permite calcular la ganancia adicional que llevaría el lazo cerrado a
la condición de estabilidad crítica. El margen de ganancia se define como:
1
= (5.47)
|()|∗
El margen de fase cuantifica el retardo de fase puro que puede agregarse antes de alcanzar
la condición de estabilidad crítica. El margen de fase se define como:
= ∗ − (−180◦ ) (5.48)
203
Así, el margen de fase y el margen de ganancia, los cuales se pueden determinar fácilmente
por medio de Matlab, se pueden ver como factores de seguridad que se deben de considerar en el
sintonizado de controladores retroalimentados para tomar en cuenta incertidumbres en constantes
de tiempo, ganancias en estado estacionario y tiempos muertos.
El comando de Matlab calcula las magnitudes y los ángulos de fase de la respuesta en
frecuencia de un sistema en tiempo continuo, lineal e invariante con el tiempo. La sintaxis es la
siguiente.
Para el ejemplo (5.42) el diagrama de Bode se obtiene con Matlab por medio de las siguientes
instrucciones:
Con Matlab se puede realizar el diagrama de Nyquist con el comando La sintaxis es
la siguiente.
Las Figuras (5-15)-(5-17) muestran los diagramas de Nyquist, que se obtienen a través de
Matlab, para la funciones de transferencia del CSTR a lazo cerrado:
µ ¶
+2
() () = ( ()) (5.49)
2 − 225 − 225
204
1
4 dB 2 dB 0 dB -2 dB
-4 dB
0.8
6 dB
0.6 -6 dB
0.4 10 dB -10 dB
0.2
20 dB -20 dB
-0.2
-0.4
-0.6
-0.8
-1
-1.8 -1.6 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0
Figura 5-15: Diagrama de Nyquist para un CSTR a lazo cerrado con = 20
con un controlador proporcional y tres valores de la ganancia proporcional, = [2 225 30].
Para este ejemplo, la función de transferencia a lazo abierto tiene un polo positivo en −30 y de
acuerdo al criterio de Routh que se aplicó en la Sección 4.4.1, para valores de 225 el lazo
cerrado es estable. Para = 20 al aplicar el criterio de estabilidad de Nyquist se tiene que el
sistema a lazo cerrado es inestable, debido a que el punto (−1 0) se encierra una vez en sentido
de las manecillas del reloj. Para = 225 se tiene que el sistema a críticamente estable, debido
a que diagrama de Nyquist pasa por el punto (−1 0). Finalmente, para = 30 se tiene que
el sistema es estable, debido a que diagrama de Nyquist encierra el punto (−1 0) en el sentido
anti-horario una vez y la función de transferencia a lazo abierto tiene un polo en el semiplano
derecho.
205
1
4 dB 2 dB 0 dB -2 dB
-4 dB
0.8
6 dB
0.6 -6 dB
0.4 10 dB -10 dB
0.2
20 dB -20 dB
-0.2
-0.4
-0.6
-0.8
-1
-2 -1.8 -1.6 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0
Figura 5-16: Diagrama de Nyquist para un CSTR a lazo cerrado con = 225
1
4 dB 2 dB 0 dB -2 dB
-4 dB
0.8
6 dB
0.6 -6 dB
0.4 10 dB -10 dB
0.2
20 dB -20 dB
-0.2
-0.4
-0.6
-0.8
-1
-3 -2.5 -2 -1.5 -1 -0.5 0
Figura 5-17: Diagrama de Nyquist para un CSTR a lazo cerrado con = 30
206
Bode Diagram
Gm = 13.3 dB (at 5.48 rad/sec) , Pm = 101 deg (at 1.85 rad/sec)
50
Magnitude (dB)
-50
-100
-150
0
-90
Phase (deg)
-180
-270
-1 0 1 2 3
10 10 10 10 10
Frequency (rad/sec)
Nyquist Diagram
5
1
Imaginary Axis
-1
-2
-3
-4
-5
-2 -1 0 1 2 3 4 5 6
Real A xis
207
1. Escribir la instrucción: [ ] = ()
5
() = (5.50)
3 + 92 + 30 + 40
que es estable a lazo abierto, se determina el margen de fase y el margen de ganancia con
cuando se aplica una acción proporcional con = 100 es decir:
µ ¶
5
() () = ( ()) (5.51)
3 + 92 + 30 + 40
La Figura (5-18) muestra los resultados. El margen de fase es 1010◦ a una frecuencia de 185
sec y el margen de ganancia en es 133 y en magnitud 46 a una frecuencia de 548
sec Los resultados se pueden verificar al multiplicar la ganancia del controlador por el
margen de ganancia en magnitud 46019 lo cual conduce a un valor de ganancia proporcional
de = 46019 que es el valor con el cual se obtiene un sistema a lazo cerrado marginalmente
inestable, como se muestra en el diagrama de Nyquist (5-19).
Smith C.A. y Corripio, A.B. (2001). Control Automático de Procesos, Teoría y Práctica.
Limusa.
208
Baez-Lopez, D. (2006). Matlab con Aplicaciones a la Ingeniería, Fisica y Finanzas. Al-
faomega, Mexico, DF.
Bequette B.W. (2003). Process Control: Modelling, Design and Simulation. Upper Saddle
River, Prentice Hall.
209
Capítulo 6
210
matricial de un sistema para fines de control en variables de estado está dado por las ecuaciones:
•
= + (6.1)
= +
cualquier sistema dinámico lineal e invariante en el tiempo se puede expresar de la forma (6.1).
El modelo estado-espacio (6.1) consiste de dos ecuaciones:
1. Las ecuaciones de estado, que describen la evolución de los estados como una función de
la entrada y de las variables de estado y corresponden a un conjunto de ecuaciones
diferenciales ordinarias.
2. Las ecuaciones de salida, las cuales relacionan el valor de las señales de salida a las señales
de los estados y las entradas y corresponden a ecuaciones algebraicas.
Los modelos estado-espacio son una forma natural de los modelos de sistemas de ingeniería
de procesos. La razón es que los modelos de ingeniería de procesos se derivan de los principios
de conservación y los modelos resultantes se pueden transformar en la forma estado-espacio en la
mayoría de los casos de interés.
Para ilustrar como un sistema de orden se puede escribir como un modelo estado-espacio de
la forma (6.1) se considera un sistema de cuarto orden:
•••• ••• •• •
() + 1 () + 2 () + 3 () + 4 () = () (6.2)
para pasar a un sistema de ecuaciones diferenciales de primer orden se hacen los siguientes
cambios de variables:
1 =
• •
1 = = 2
• ••
2 = = 3 (6.3)
• •••
3 = = 4
• •••• ••• •• •
4 = = −1 () − 2 () − 3 () − 4 () + ()
= −1 4 − 2 3 − 3 2 − 4 1 +
211
el conjunto de ecuaciones de primer orden (6.3) se pueden escribir en forma matricial como:
⎡ ⎤ ⎡ ⎤
0 1 0 0 0
⎢ ⎥ ⎢ ⎥
⎢ ⎥ ⎢ ⎥
• ⎢ 0 0 1 0 ⎥ ⎢ 0 ⎥
= ⎢
⎢
⎥ + ⎢ ⎥
⎥ ⎢ ⎥ (6.4)
⎢ 0 0 0 1 ⎥ ⎢ 0 ⎥
⎣ ⎦ ⎣ ⎦
−4 −3 −2 −1
h i
= 1 0 0 0
6.3.1. Controlabilidad
De una manera informal se puede decir que la controlabilidad es una característica del sistema
que indica la capacidad que tiene una señal de entrada para influir en la evolución del estado.
Formalmente, se dice que un sistema es controlable en el tiempo 0 si se puede llevar de cualquier
estado inicial (0 ) a cualquier otro estado, mediante un vector de control en un intervalo de
tiempo finito.
212
Se considera el siguiente sistema:
̇ = + (6.5)
= +
6.3.2. Observabilidad
̇ = (6.7)
=
213
descrito mediante las ecuaciones (6.7) es completamente observable sí y sólo sí la matriz de ×:
⎡ ⎤
⎢ ⎥
⎢ ⎥
⎢ ⎥
= ⎢
⎢ ..
⎥
⎥ (6.8)
⎢ . ⎥
⎣ ⎦
−1
Cálculo de la controlabilidad.
Cálculo de la observabilidad:
214
6.3.4. Controlabilidad de un sistema de tanques agitados en serie
Se considera el ejemplo de dos tanques agitados en serie que se ilustra en la Figura (6-1). La
salida del tanque 1 1 entra al tanque 2 que tiene como salida la corriente 2 . Se desea evaluar
la capacidad de una entrada de control para controlar el nivel de los tanques 1 y 2 Se examinan
dos ubicaciones de la entrada de control: (i) la entrada de control es un flujo que entra al segundo
tanque y (ii) la entrada de control es el flujo de alimentación al primer tanque.
1 1 p
= − 1 (6.9)
1
2 1 p 2 p
= 1 − 2 +
2 2
1
= −1 1 (6.10)
2
= 1 1 − 2 2 +
= [ ]
donde: ⎡ ⎤⎡ ⎤ ⎡ ⎤
−1 0 0 0
= ⎣ ⎦⎣ ⎦=⎣ ⎦
1 −2 1 −2
215
(ii) entrada de flujo al
tanque 1
Nivel en el
tanque 1, h1
Nivel en el
tanque 1, h2
Corriente de
salida tanque 2
216
de modo que: ⎡ ⎤
0 0
= ⎣ ⎦ (6.12)
1 −2
1 1 p
= − 1 (6.13)
1
2 1 p 2 p
= 1 − 2
2 2
1
= − 1 1 (6.14)
2
= 1 1 − 2 2
en este caso: ⎡ ⎤⎡ ⎤ ⎡ ⎤
−1 0 1 −1
= ⎣ ⎦⎣ ⎦=⎣ ⎦
1 −2 0 1
de modo que: ⎡ ⎤
0 −1
= ⎣ ⎦ (6.16)
1 1
Del ejemplo se puede concluir que cuando la entrada de control se ubica en el primer tanque
el sistema es controlable, lo que implica que es posible llevar el estado del sistema, 1 y 2 a otro
estado deseado por medio de una entrada de control Cuando la entrada de control se ubica
217
en el segundo tanque no se cumple la propiedad de controlabilidad, lo cual se puede verificar
intuitivamente, ya que la ubicación de la entrada de control en el segundo tanque no tiene ningún
efecto en el estado del primer tanque y no es posible llevarlo a algún estado deseado con la
manipulación de la entrada de control.
Para el caso (ii) del ejemplo anterior se desea evaluar la observabilidad para dos casos: (i)
se dispone solo de la medición de la salida del tanque 1, = 1 1 y (ii) se dispone solo de la
medición de la salida del tanque 2, = 2 2 .
1
= − 1 1
2
= 1 1 − 2 2 (6.17)
= 1 1
donde:
⎡ ⎤
h i −1 0
= 1 0 ⎣ ⎦
1 −2
h i
= −12 0
218
de modo que: ⎡ ⎤
1 0
= ⎣ ⎦ (6.19)
−12 0
1
= − 1 1
2
= 1 1 − 2 2 (6.20)
= 2 2
en este caso:
⎡ ⎤
h i −1 0
= 0 2 ⎣ ⎦
1 −2
h i
= 1 2 −22
de modo que: ⎡ ⎤
0 2
= ⎣ ⎦ (6.22)
1 2 −22
219
6.4. Formas canónicas
La representación en variables de estado de un sistema no es única. De las diferentes formas de
la representación en variables de estado, algunas formas son más convenientes que otras para algún
fin específico. A estas formas particulares que ponen en evidencia características específicas del
sistema o que facilitan la aplicación de algún procedimiento se les conoce como formas canónicas.
Las formas canónicas se obtienen a partir de una transformación de estado o transformación de
similaridad, con la elección de una nueva variable de estado de la forma:
= (6.23)
donde es una matriz no singular. Para estados iniciales cero y entradas iguales, la respuesta de
las formas canónicas y cualquier representación de variables de estado de la forma (6.1) coincide.
••• •• • • ••
() + 1 () + 2 () + 3 () = 3 + 2 () + 3 () (6.26)
220
la forma canónica controlable es:
⎡ ⎤ ⎡ ⎤
0 1 0 0
⎢ ⎥ ⎢ ⎥
⎢ ⎥ ⎢ ⎥
̇ = ⎢ 0 0 1 ⎥ + ⎢ 0 ⎥ (6.27)
⎣ ⎦ ⎣ ⎦
−3 −2 −1 1
h i
= 3 2 1
Para el sistema general de orden (6.24), la forma canónica observable está dada por:
⎡ ⎤ ⎡ ⎤
0 −
⎢ ⎥ ⎢ ⎥
⎢ ⎥ ⎢ ⎥
⎢ −−1 ⎥ ⎢ ⎥
̇ = ⎢ ⎥ + ⎢ −1 ⎥ (6.28)
⎢ ⎥ ⎢ ⎥
⎢ ⎥ ⎢ ⎥
⎣ ⎦ ⎣ ⎦
−1 −1 1
h i
= 0 0 1
221
forma:
= − (6.30)
donde es una entrada externa escalar (o referencia), es una ganancia escalar para la pre-
compensación de la entrada externa y es un vector de ganancias retroalimentadas.
Sin perdida de generalidad si se considera igual a cero (lo cual ocurre comúnmente), el
sistema a lazo cerrado que se forma de aplicar (6.30) en (6.1) está dado por:
•
= ( − ) + () (6.31)
=
la matriz
= − (6.32)
| − | = 0 (6.33)
para el sistema a lazo cerrado (6.31) la estabilidad se determina de las raíces de la ecuación
característica que se obtienen de la expresión:
de modo que es posible modificar las propiedades de estabilidad del sistema a lazo abierto (6.1)
a través del retroalimentado de estados (6.30) por medio de una elección adecuada del vector de
ganancias .
El retroalimentado de estados de la forma (6.30) tiene las siguientes propiedades,
222
4. Modifica los polos del sistema (6.1).
2. Elegir polos deseados a lazo cerrado. Donde es el orden del sistema a lazo abierto (6.1).
3. Calcular el polinomio característico a lazo cerrado a partir de los polos a lazo cerrado
deseados:
( − 1 )( − 2 )( − ) = + 1 −1 + + (6.35)
223
5. Igualar la ecuación (6.35) con la ecuación característica a lazo cerrado (6.36). Al desarrollar
el producto e igualar coeficientes se obtendrá un sistema de ecuaciones con incógnitas.
6. Resolver el sistema de ecuaciones con incógnitas para obtener los valores del vector de
ganancias .
Para sistemas de la forma (6.1), de una entrada y una salida, Ackermann derivo una formula
para calcular la matriz de ganancias:
que permiten que una ley de control de retroalimentados de estados de la forma (6.30) ubique los
polos de un sistema a lazo cerrado en cualquier posición deseada.
La formula de Ackermann se desarrolla en 5 pasos:
3. A partir de los polos que se desean ubicar con el retroalimentado de estados (6.30) construir
un polinomio de grado de la forma:
224
donde −1 es la inversa de la matriz de controlabilidad.
que tiene polos a lazo abierto en −075 y 30 es decir es un sistema inestable a lazo abierto. Se
desean colocar los polos a lazo cerrado en −30 y −20 con un retroalimentado de estados de la
forma (6.30).
Inspección directa
1. Controlabilidad:
= [ ]
⎡ ⎤⎡ ⎤
−2 −003125 0
= ⎣ ⎦⎣ ⎦
200 425 1
⎡ ⎤
−00315
= ⎣ ⎦
425
⎡ ⎤
0 −003125
= ⎣ ⎦ (6.42)
1 425
225
4. Ecuación característica a lazo cerrado:
| − − | = 0
¯ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ¯
¯ h i¯
¯ 1 0 −2 −003125 0 ¯
¯ ⎣ ⎦−⎣ ⎦−⎣ ⎦ 1 2 ¯ = 0
¯ ¯
¯ 0 1 200 425 1 ¯
¯⎡ ⎤ ⎡ ⎤ ⎡ ⎤¯
¯ ¯
¯ 0 −2 −003125 0 0 ¯
¯⎣ ⎦−⎣ ⎦−⎣ ⎦¯ = 0
¯ ¯
¯ 0 200 425 1 2 ¯
¯⎡ ⎤¯
¯ ¯
¯ +2 003125 ¯
¯⎣ ⎦¯ = 0
¯ ¯
¯ −200 − 1 − 425 − 2 ¯
2 + (−2 − 225) + (−22 + 0031 251 − 225) = 0 (6.44)
6. Resolver el sistema:
5 = −225 − 2
Formula de Ackermann
226
3. Polinomio de grado :
() = 2 + 5 + 6 (6.48)
() = 2 + 5 + 6 (6.50)
⎡ ⎤2 ⎡ ⎤ ⎡ ⎤
−2 −003125 −2 −003125 1 0
= ⎣ ⎦ + 5⎣ ⎦ + 6⎣ ⎦
200 425 200 425 0 1
⎡ ⎤
−63 −02
= ⎣ ⎦
1450 39
5. Vector de ganancias :
227
1. Matrices y : = [−2 − 003125; 200 425]; = [0 1];
•
e = e
+ + ( − e
) (6.52)
e es un estimado del estado del sistema (6.1) y es una matriz de observación. Se puede
donde
e más
ver que el observador (6.52) es una copia del sistema (6.1) en términos del estado estimado
una corrección dada por un gradiente ( − e
)
De (6.52) y (6.1) se obtiene:
• •
e) = ( −
( − e) − ( − e
) (6.53)
228
como = :
• •
e) = ( −
( − e) − ( − e
) (6.54)
e) − ( −
= ( − e)
•
= ( − ) (6.57)
de modo que con una elección adecuada de la matriz de observación es posible que todos los
valores propios de − tengan parte real negativa y por lo tanto:
e) = 0
lı́m ( −
→∞
conforme → ∞
El retroalimentado de estados que se forma con los estados estimados está dado por:
= − e
(6.60)
229
6.8. Referencias básicas
Dos referencias importantes que el lector puede consultar para más información sobre los temas
que se presentan en este capítulo son:
230
Capítulo 7
Epílogo
Entender las bases de las técnicas de control modernas por técnicas estado-espacio.
231
Usar el Matlab y Simulink para cálculos de simulación y control de procesos.
En la parte final de cada capítulo se proporcionaron referencias de libros de textos básicos que
el lector puede consultar para ampliar la información que se presenta aquí. Entre los temas que
no se abordan y merecen atención para alumnos interesados en el tema de simulación y control de
procesos se encuentran:
232
Bibliografía
[1] Alhumaizi, K., Aris R. (1994). Chaos in a simple two-phase reactor. Chaos Solitons Fractals
4, 1985-2014.
[2] Alvarez, J., Alvarez, J. (1988). Analysis and control of fermentation processes by optimal
and geometric methods. Proceedings of the 1988 American Control Conference, June 1988,
Atlanta, GA, 1112-1117.
[3] Alvarez-Ramirez, J., Puebla, H. (2001). On classical PI control of chemical reactors. Chem.
Eng. Sci. 56, 2111.
[4] Alvarez-Ramirez, J., Puebla, H., Espinosa, G. (2001). A cascade control strategy for a space
nuclear reactor system. Ann. Nuclear Energy 28, 93.
[5] Alonso, A., Yidstie, E. (1996). Process systems, passivity and the second law of thermody-
namics. Comp. Chem. Eng. 20, S1119-S1124.
[6] Andrews, J. F. (1968). A mathematical model for the continuous culture of microorganisms
utilizing inhibitory substrate. Biotechnol. Bioeng. 10, 707-723.
[7] Ang, S., Braatz, R.D. (2004). Experimental projects for the process control laboratory. Chem.
Eng. Educ. 36, 182-185.
[8] Aris, R. (1969). Elementary Chemical Reactor Analysis. Prentice Hall, Inc., Englewood Cliffs,
N.J.
[9] Aris, R. (1994). Mathematical Modelling Techniques. Dover Publication Inc., New York.
[10] Aris R., Amundson, N.R. (1958). An analysis of chemical reactor stability and control-I, II,
III. Chem. Eng. Sci. 7, 121-155.
233
[11] Astrom, K.J., Wittenmark, B. (1990). Computer-Controlled Systems: Theory and Design.
Prentice Hall, Englewood Cliffs, New Jersey.
[13] Baez-Lopez, D. (2006). Matlab con Aplicaciones a la Ingeniería, Fisica y Finanzas. Al-
faomega, Mexico, DF.
[14] Bailey, J.E., Ollis, D.F. (1987). Biochemical Engineering Fundamentals. McGraw Hill, Inc.,
New York.
[15] Bastin, G., Dochain,D. (1990). On L ine Estimation and Adaptive Control of Bioreactors.
Elsevier, New York.
[16] Bequette B.W. (1997). An undergraduate course in process dynamics. Comp. Chem. Eng. 21
suppl., S261-S266.
[17] Bequette, B.W. (1991). Nonlinear control of chemical processes: A review. Ind. Eng. Chem.
Res. 30, 1391-1413.
[18] Bequette B.W. (2003). Process Control: Modelling, Design and Simulation. Upper Saddle
River, Prentice Hall.
[19] Bequette B.W. (2005). A laptop-based studio course for process control. IEEE Cont. Sys.
Magazine. 25, 45-49.
[20] Burns, R.S. (2001). Advanced Control Engineering. Butterworth-Heinemann, MA, USA.
[21] Christofides, P.D. (2001). Control of nonlinear distributed process systems: Recent develop-
ments and challenges. AIChE J 47, 514-518.
[22] del-Muro Cuellar, B., Velasco-Villa, M., Puebla, H., Alvarez-Ramirez, J. (2005). Model ap-
proximation for dead-time recycling systems. Ind. Eng. Chem. Res. 44), 4336-4343.
[23] Denbigh, K.G., Turner, J.C.R. (1972). Chemical Reactor Theory: An Introducion. Cambridge
University Press., New York.
234
[24] Dochain, D., Babary, J.P., Tali-Maamar, N. (1992). Modeling and adaptive control of non-
linear distributed parameter bioreactors via orthogonal collocation. Automatica 28, 873-883.
[25] Dorato, P.(1987). A historical review of robust control. IEEE Cont. Syst. Magazine April,
44-47.
[26] Doran, M.P. (1995). Bioprocess Engineering Principles. Academic Press Limited, San Diego,
C.A.
[27] Edgar T.F., Ogunnaike B.A., Muske K.R. (2006). A global view of graduate process control
education. Comp. Chem. Eng. 30, 1763-1774.
[28] Elnashaie, S.S.E.H., Elshishini, S.S. (1993). Modelling, Simulation and Optimization of In-
dustrial Fixed Bed Catalytic Reactors. Topics in Chemical Engineering. Edited by R. Hughes.
Gordon and Breach Science Publishers. OPA (Amsterdam).
[29] Farlow, S.J. (1993). Partial Differential Equations for Scientists and Engineers. Dover Pub-
lications, Inc. New York.
[30] Ferino, I., Rombi, E. (1999). Oscillating reactions. Catalysis Today 52, 291-305
[31] Fisher, D.G. (1991). Process control: an overview and personal perspective. Can. J. Chem.
Eng. 69, 5-26.
[32] Froment, G.F., Bischoff, K.B. (1990). Chemical Reactor Analysis and Design. John Wiley
and Sons, Inc.
[33] Garyhan, P., Elnashaie S. (2004). Static/dynamic bifurcation and chaotic behavior of an
ethanol fermentor. Ind. Eng. Chem. Res. 43, 1260.
[34] Glass, L., Mackey, M.C. (1988). From Clocks to Chaos: The Rhythms of Life. Princeton
University Press.
[35] Heafner, J.W. (1996). Modeling Biological Systems. Principles and Applications. Chapman &
Hall, New York, NY.
[36] Ingham, J. (2007). Chemical Engineering Dynamics: An Introduction to Modeling and Com-
puter Simulation. Wiley-VCH.
235
[37] Jensen, K.F., Ray, W.H. (1982). The bifurcation behavior of tubular reactors. Chem. Eng.
Sci. 37, 199-222.
[38] Kailath, T. (1980). Linear Systems. Prentice Hall, Englewood Cliffs, New Jersey.
[40] Khalil, H.K. (1996). Nonlinear Systems. Prentice-Hall. Upper Saddle River, N.J.
[41] Klatt, K.-U., Engell, S. (1998). Gain scheduling trajectory control of a continuous stirred
tank reactor. Comp. Chem. Eng. 22, 491-502.
[42] Lee, T-S., Pham, M.Q., Weigand, W.A., Harvey, S.P., Bentley, W.E. (1996). Bioreactor strate-
gies for the treatment of growth-inhibitory waste: an analysis of thiodiglycol degradation, the
main hydrolysis product of sulfur mustard. Biotechnol. Prog. 12, 533-539.
[43] Lee, J.S., Chang, K.S. (1996). Applications of chaos and fractals in process systems engineer-
ing. J. Proc. Cont. 6, 71-87.
[44] Luyben, W.L. (1990). Process Modeling, Simulation and Control for Chemical Engineers.
McGraw-Hill, New York, NY.
[45] Morari, M., Zafiriou, E. (1989). Robust Process Control. Prentice-Hall, New York.
[46] Morari, M., Gentilini, A. (2001). Challenges and opportunities in process control: biomedical
processes. AIChE J. 47, 2140—2143.
[48] Murray, R.M. (2003). Control in an Information Rich World: Report of the Panel on
Future Directions in Control, Dynamics, and Systems. SIAM 2003 [Online]. Available:
http://www.cds.caltech.edu/~murray/cdspanel.
[49] Ogata, K. (1995). Discrete-Time Control Systems. Prentice-Hall, Upper Saddle River, NJ.
[50] Ogata, K. (1997). Modern Control Engineering. Prentice-Hall, Upper Saddle River, NJ.
[51] Ortega, R., Loria, A., Nicklasson, P., Sira-Ramiírez, H. (1998). Passivity Based Control of
Euler L agrange Systems. Springer-Verlag, London.
236
[52] Perlmutter, D.D. (1972). Stability of chemical reactors. Prentice-Hall, Englewood Cliffs, N.J.
[53] Puebla, H., Alvarez-Ramirez, J. (2005). A cascade feedback control approach for hypnosis.
Annals Biomed. Eng. 33, 1449-1463.
[54] Riascos, C.A.M., Pinto, J.M. (2004). Optimal control of bioreactors: a simultaneous approach
for complex systems. Chem. Eng. J. 99 23—34.
[55] Sastry, S. (1999). Nonlinear Systems: Analysis, Stability, and Control. Springer.
[56] Shelden R.A., Dunn I.J. (2001). Dynamic simulation: modeling processes, the environment,
the world. Chem. Eng. Progress December, 44-48.
[57] Shinskey, F.G. (1988). Process Control Systems. McGraw-Hill, New York.r-Verlag, 1999.
[58] Shinskey F.G. (2002). Process control: as taught vs. as practiced. Ind. Eng. Chem. Res. 41,
3745-3750.
[59] Skogestad S., Morari M. (1988). Understanding the dynamic behavior of distillation columns.
Ind. Eng. Chem. Res. 27, 1848-1862.
[60] Slotine, J.-J. E., Li, W. (1991) Applied nonlinear control. Prentice Hall, London.
[61] Stephanopoulos, G. (1984). Chemical Process Control: An Introduction to Theory and Prac-
tice. Prenctice-Hall, NJ, USA.
[62] Subramanian, S., Balakotaiah, V. (1996). Classification of steady-state and dynamic behavior
of distributed reactor models. Chem. Eng. Sci. 51, 401-421.
[63] Szederkenyi, G., Kristensen, N.R., Hangos, K.M., Bay Jorgensen, S.B. (2002). Nonlinear
analysis and control of a continuous fermentation process. Comp. Chem. Eng. 26, 659—670.
[64] Torres-Robles, R., Castro-Arellano, J.J. (2003). Análisis y Simulación de Procesos de Refi-
nación del Petróleo. Alfaomega, México.
[65] Uppal, A., Ray, W.H., Poore, A.B. (1974). On the dynamic behaviour of continuous stirred
tank reactors. Chem. Eng. Sci. 29, 967.
237
[66] Utz, T., Hagenmeyer, V., Mahn, B. (2007). Comparative evaluation of nonlinear model predic-
tive and flatness-based two-degree-of-freedom control design in view of industrial application.
J. Process Cont.17, 129—141.
[67] Van Heerden, C. (1953). Autothermic process. properties and reactor design. Ind. Eng. Chem.
45, 1242-1253.
[68] Vidyasagar, M. (1978). Nonlinear Systems Analysis. Prentice Hall, Englewood Cliffs, NJ.
[69] Yin, J., Jensen, M.K. (2003). Analytic model for transient heat exchanger response. Int. J.
Heat Mass Transfer 46, 3255—3264
238
Apéndice A
Introducción a Matlab1
Para iniciar se hace un doble clic sobre el icono de Matlab en la pantalla principal de la PC,
lo cual despliega una pantalla similar a la que se muestra en la Figura (A-1).
En la pantalla inicial se puede visualizar el directorio de trabajo en la parte superior derecha.
El espacio de trabajo, donde se introducen las funciones en instrucciones para trabajar se ubica
adelante del símbolo con el cursor parpadeando. La ayuda en línea se puede llamar con el
comando help, con la tecla F1, o bien con el icono de ayuda de la barra de herramientas.
Algunas herramientas básicas de Matlab son las siguientes:
239
Figura A-1: Pantalla inicial de Matlab.
240
2. Operaciones con matrices: Además de las operaciones básicas de matrices, Matlab cuenta
con los siguientes comandos para operaciones con matrices:
Para las operaciones aritméticas con vectores y matrices, además de las operaciones arit-
méticas que aplica a escalares, Matlab acepta la multiplicación, división y exponenciación
de elemento por elemento,
donde los elementos de este vector son los coeficientes del polinomio en orden descendente.
En esta notación, se deben de incluir los coeficientes de términos cuyo coeficiente es cero.
Para operaciones aritméticas básicas Matlab reconoce los siguientes símbolos:
polyint() integración
241
4. Archivos : Matlab reconoce dos tipos de archivos que se conocen como archivos m debido a
que usan la extensión m: (i) Archivos de escritura: Un archivo de escritura consiste archivos
que almacenan instrucciones, expresiones o funciones sin que se admitan argumentos de
entrada o salida. (ii) Archivos de función. Un archivo de función es una función definida por
el usuario que usa Matlab a través de argumentos de entrada para proporcionar argumentos
de salida. La utilidad de los archivos de función es almacenar funciones u operaciones que
sean repetitivas y poder utilizarlas o modificarlas cuando sea necesario. Para crear un archivo
m del tipo función se siguen los siguientes pasos:
= ()
En las líneas siguientes se escriben las ecuaciones o instrucciones a resolver con Matlab
con el formato general,
= ()
Grabar el archivo de función. El nombre del archivo de función debe ser el mismo que
el nombre de la función más la extensión .
Para profundizar en el uso de Matlab para cálculos de ingeniería se recomienda al lector con-
sultar algunas de las muchas referencias que existen sobre el tema. La ayuda que ofrece en línea
Matlab es también muy recomendable. Para detalles sobre comandos específicos el lector puede
utilizar la función de Matlab y consultar el material relacionado al comando deseado. Por ejemplo:
help plot, help subplot, help if. Sin embargo, muchos lectores identificaran en forma
obvia el uso de la mayoría de los comandos de Matlab, debido a que su programación es muy simple
e intuitiva. A continuación se presentan la sintaxis y el uso de comandos de uso más frecuente en
este texto.
242
Funciones básicas de escritorio de Matlab
. help ... despliega en pantalla la información sobre el comando o función deseada.
. Termina la sesión actual de Matlab y sale del programa. Al usar el comando quit la
historia de los comandos que se usaron previamente en Matlab se mantiene en la memoria
de Matlab para la siguiente sesión.
@ Crea una función para su manejo en archivos de Matlab de modo que se pueden pasar
entre archivos como si fueran variables.
Tres puntos. Indican un cambio de línea dentro de una instrucción o archivo en Matlab.
Útil para escribir en varias líneas instrucciones muy grandes y poder visualizarlas en una
sola pantalla.
243
Operaciones con matrices en Matlab.
. Crea una matriz en la cual todos los elementos son iguales a cero.
. Una lazo que ejecuta una serie de acciones hasta concluir un contador y continúa
la ejecución del programa con la instrucción .
. Un lazo que ejecuta una serie de acciones si se cumple algun condicionante.
Figuras en Matlab
. Divide la ventana de la figura actual en varios panales. Se especifica por el número
de paneles P y su localización en formato renglón R y columna C, es decir, subplot(P,R,C).
. El comando plot(x,y) genera una figura para los conjuntos de datos de igual tamaño
en y .
Especificaciones de color: ‘y’ amarillo, ‘r’ rojo, ‘g’ verde, ‘b’ azul, ‘k’ negro, etc.
244
Especificaciones de línea: ‘−’ línea continua, ‘- -’ línea discontinua, ‘:’ línea punteada, ‘+’
datos marcados con el símbolo +, ‘’ datos marcados con un cuadro, ‘’ datos marcados con
un circulo, etc.
245
Apéndice B
Introducción a Simulink1
246
Figura B-1: Pantalla inicial de Simulink.
247
Figura B-2: Ventana de contrucción de un modelo en Simulink.
Cada aplicación contiene conjuntos de bloques y funciones definidas que se despliegan al hacer
doble clic sobre cada aplicación. Los bloques y funciones se despliegan a la derecha de la pantalla
(B-1). Algunos de los bloques y funciones definidas de mayor uso son los siguientes:
1. Integrador (Integrator), localizado en Simulink → Commonly Used Blocks. Sirve para inte-
grar una señal o función continua.
248
2. Scope (Scope), localizado en Simulink → Commonly Used Blocks. Sirve para visualizar en
forma gráfica resultados de la simulacion del modelo.
6. Reloj (clock), localizado en Simulink → Sources. Sirve para generar un vector del tiempo
de simulación del modelo.
7. Leer datos desde el espacio de trabajo de Matlab (From Workspace), localizado en Simulink
→ Sources. Sirve para leer datos de entrada al modelo desde el espacio de trabajo de Matlab.
8. Escalón (Step), localizado en Simulink → Sources. Sirve para generar una señal escalón.
9. Función (Fcn), localizado en Simulink → User-Defined Functions. Sirve para introducir una
función definida por el usuario.
11. Controlador PID (PID Controller), localizado en Simulink Extras → Aditional Linear. Sirve
para introducir un controlador PID.
249
3. Hacer las conexiones de los bloques del modelo al seleccionar y mantener con el botón
izquierdo del ratón de la computadora la salida o entrada de cada bloque y arrastrar hasta
otra salida, entrada o línea.
4. Para terminar de definir el modelo se introducen los cambios necesarios en los parámetros de
cada bloque. Esto se hace a través de dar un doble clic sobre cada bloque a ser modificado.
250
Apéndice C
1. Formar grupos de 4 integrantes: Los grupos promueven la interacción entre ellos y com-
plementar sus debilidades y fortalezas de los conocimientos del curso. El trabajo en grupo
también permite que los estudiantes desarrollen trabajos con un mayor grado de complejidad
respecto al desarrollo de proyectos individuales.
2. Búsqueda o asignación de un caso de estudio: Debido a que se permite que los estudiantes
seleccionen un caso de estudio de una amplia variedad de problemas, ellos se encuentran
más motivados para entender y desarrollar el caso de estudio de su elección. Los casos de
estudio que se han desarrollado incluyen reactores químicos y biológicos, intercambiadores
de calor, rutas metabólicas y celulares. Es muy importante verificar que el caso de estudio
presente el modelo dinámico y que los datos de los parámetros del modelo estén completos.
251
mana del curso. En la semana 11 se entrega un reporte final del proyecto en archivo PDF
siguiendo las guías de autores de revistas internacionales (titulo, autores, resumen, palabras
clave, secciones, referencias, letra a 12 puntos, márgenes amplios, doble espacio, referencias
estandarizadas, lista de figuras, figuras, etc.).
2. Modelado matemático: En esta sección se presenta el modelo matemático del caso de estudio.
Se debe de justificar la derivación del modelo matemático e identificar los términos que
contiene. Aquí también se obtienen los puntos de equilibrio del modelo matemático.
4. Análisis en el dominio del tiempo: En esta sección se realiza la linealización del modelo
matemático, la evaluación de la estabilidad a lazo abierto, el plano de fase y un estudio de
sensibilidad paramétrica.
5. Diseño y análisis del control clásico PID: En esta sección se obtiene la función de transferencia
del sistema a lazo abierto y cerrado, se obtiene los diagramas de localización de raíces, Bode
y de Nyquist y se diseña e implementa un control tipo PID.
6. Diseño y análisis del control con técnicas en el espacio de estados: En esta sección se evalúan
las propiedades de controlabilidad y observabilidad y se deriva una ley de control por retro—
alimentado de estados para el caso de estudio.
7. Conclusiones: En esta sección se presentan los comentarios finales de los resultados obtenidos
con el caso de estudio y una reflexión sobre el procedimiento y la utilidad del proyecto.
252