From Zero To Hero: Funcionalidades Prácticas de para El Análisis y Control de Procesos
From Zero To Hero: Funcionalidades Prácticas de para El Análisis y Control de Procesos
funcionalidades prácticas
de Matlab/Simulink para el
análisis y control de procesos
Comité editorial
Alberto Amaya
Patricia Noguera
Fabio Pavas
Veronique Bellanger
Fredy Chaparro
Jairo Peña
Pedro Benjumea
Edición
Editorial Universidad Nacional de Colombia
direditorial@unal.edu.co
www.editorial.unal.edu.co
Equipo editorial
Coordinación editorial: José Francisco Rengifo
Corrección de estilo: Marcela Garzón Gualteros
Pauta gráfica en LaTeX: Juan Carlos Villamil y Óscar Andrés Prado-Rubio
Diagramación: Óscar Andrés Prado-Rubio
Salvo cuando se especifica lo contrario, las figuras y tablas del presente volumen son
propiedad del autor.
Prefacio 9
5
Contenido
6
Contenido
Referencias 249
7
Prefacio
Sin duda, el control de procesos es un área interdisciplinaria y puede ser considerada como
una de parte vital de la Ingeniería Química. El control de procesos se encuentra en constante
evolución guiada por una cantidad creciente de contribuciones de investigación, la reducción
del costo de dispositivos digitales y un incremento del poder computacional. En consecuencia,
es necesario que los ingenieros químicos hagan frente a los desarrollos de automatización para
proponer un mejor diseño y operación de plantas modernas.
Con esto en mente, los cursos de control de procesos convencionales, dentro del programa de
pregrado de Ingeniería Química, se centran en temas como modelado y análisis dinámico de
procesos, instrumentación y teoría clásica de control de procesos. A pesar de la relevancia del
control de procesos, después de un par de años en frente de la asignatura, para mí es claro
que este no es uno de los cursos más populares del programa. Esto se debe principalmente a
que el análisis dinámico no se aborda normalmente en otros cursos, la terminología de control
fue creada por los desarrolladores de automatización (e.g., ingenieros electrónicos, mecánicos
y de sistemas), la relativa complejidad de la formulación matemática (los estudiantes no usan
con frecuencia el dominio de Laplace) y falta de habilidades de simulación.
Esas razones hacen que los estudiantes se vean saturados de nueva información durante el
curso y tengan dificultades para seguir las temáticas de forma exhaustiva. Por tanto, no mu-
chos se sienten alentados a seguir proyectos, estudios de posgrado y, finalmente, aspirar a una
carrera asociada al control de procesos. Se han hecho esfuerzos sustanciales para abordar los
problemas antes mencionados: por ejemplo, dar más relevancia a la comprensión de sistemas
dinámicos dentro del programa de Ingeniería Química, emplear literatura más amigable para
la enseñanza de la teoría de control de procesos y el refinamiento de temas para pensar en el
curso de Control Automático de Procesos. Sin embargo, todavía hay espacio para mejoras en
cuanto a cómo explotar herramientas computacionales como Matlab y Simulink en el curso de
Control Automático de Procesos. Ahora, una de las dificultades para esto es que la información
sobre implementaciones de software está dispersa (por ejemplo, documentación de soporte de
Matlab, pocos libros, páginas web) y es complejo desarrollar adecuadamente las habilidades
de simulación para análisis de sistemas dinámicos y la implementación de arquitecturas de
control.
9
en el concepto y el análisis en lugar de centrarse en las implementaciones de problemas en el
software. El contenido del libro ha sido seleccionado siguiendo las necesidades evidenciadas
en libros de texto de control de procesos más recientes. Por tanto, los estudiantes pueden
hacer un paralelo directo de las implementaciones con la teoría mediante un modelo ilustra-
tivo relevante para la Ingeniería Química, que los alienta a ir más allá de los casos de estudio
propuestos. El caso seleccionado es la reacción de Van de Vusse cuando se lleva a cabo un
reactor contínuo de mezcla completa (CSTR). Como tal, esta obra no fue pensado para que
los estudiantes adquieran los conceptos de cursos Control Automático de Procesos, pues la
literatura para esto está bien fundamentada; este libro muestra directamente con ejemplos
cómo resolver los problemas matemáticos asociados al control automático de procesos, con
especial atención en la herramienta de Simulink. De esta manera, aquellos estudiantes que
tienen falencias en programación son alentados a sobrepasar sus dificultades y a familiarizarse
con métodos sistemáticos de programación para control de procesos.
En el presente documento, todas las figuras fueron realizadas por el autor modificando de
capturas de pantalla de la interfaz del software Matlab/Simulink o gráficas de resultados del
mismo software, a excepción de aquellas donde se señale lo contrario.
10
Capítulo 1
Resumen
En esta sección se introducirán los ambientes de Matlab para versiones posteriores al 2018
que pueden ser utilizados para propósitos de control automático de procesos. Dependien-
do de la aplicación, un ambiente es más conveniente que otro. De manera convencional,
se utiliza el ambiente de programación en código útil directamente desde la ventana de
comandos o como archivo M-file. Teniendo en cuenta el nivel de programación esperado
en Matlab por parte de los estudiantes, no se profundiza. Seguidamente, se introduce una
herramienta simbólica y se presenta donde es útil. Por otro lado, se desglosa una herra-
mienta de Matlab diseñada específicamente para facilitar la implementación y evaluación
de estructuras control automático llamada Simulink. Por último, se muestra cómo los
ambientes de programación se pueden conectar y realizar programación jerárquica. Cabe
resaltar que, como tales, las herramientas numéricas y simbólicas utilizadas pertenecen a
Matlab y no a un ambiente determinado.
Por otro lado, no es tan común el uso del motor simbólico de Matlab. Dentro del curso de
control, algunas herramientas simbólicas complementan los esfuerzos numéricos en aras de
evitar obtener soluciones particulares a problemas que hacen más lentas las rutinas de solu-
ción. Adicionalmente, el motor simbólico de Matlab puede ser utilizado para validar cálculos
11
Capítulo 1. Introducción a las herramientas de Matlab
Desde este menú es posible iniciar un modelo desde cero o acceder a los ejemplos que trae
Matlab; también se selecciona el modelo en blanco para iniciar la interfaz de trabajo. La ven-
tana emergente se muestra en la figura 1.3a. Para complementar el trabajo, las herramientas
necesarias se encuentran en Simulink library browser, que se accede usando el botón marcado
en rojo y se muestra en la figura 1.3b. A modo de información, sin intención de empezar a
trabajar en Simulink, se introduce la herramienta. Por defecto, la librería de Simulink está
dividida en subcarpetas; en cada subcarpeta o sublibrería se encuentran los bloques predefi-
nidos que se usan habitualmente en la simulación dinámica de un sistema.
Las sublibrerías de Simulink son: Commonly used blocks, Continuous, Dashboard, Disconti-
nuities, Discrete, Logic and Bit Operations, Look-Up Tables, Math Operations, Model Verifi-
cation, Model-Wide Utilities, Ports and Subsystems, Signal Attributes, Signal Routing, Sinks,
12
Capítulo 1. Introducción a las herramientas de Matlab
Sources, User-Defined Functions y Additional Math and Discrete. Simulink mediante los blo-
ques permite usar diferentes representaciones de los modelos: ecuaciones diferenciales, funcio-
nes de transferencia, espacio de estados, representación cero-polo, entre otros. Aparecen otras
librerías para llevar a cabo aplicaciones más allá del alcance del curso de control de procesos
y, por ende, no se explorarán en este libro.
Dependiendo del tipo de simulación por realizar, se requerirán diferentes tipos de bloques y
configuraciones. Los bloques pueden añadirse a un modelo con clic sostenido, arrastrándolos
hasta la ventana de trabajo. Alternativamente, con el cursor se selecciona el bloque que se va
a usar y se digita Ctrl+I, o se da un clic derecho en el bloque y un clic izquierdo a “Add block
to new model”. Para conectar los bloques, se acerca el cursor a una de las entradas o salidas
hasta que este cambie en forma de cruz; se da clic izquierdo y de manera sostenida se arrastra
en dirección del tope contrario del otro bloque. La conexión es correcta cuando la línea deja
de ser roja y se vuelve continua. La punta de flecha de la línea indica la salida de la señal,
mientras el otro extremo es la entrada.
Las especificaciones mínimas requeridas de cada bloque se relacionan con la operación que
realizan dentro de la simulación del sistema. Para observar los parámetros necesarios del blo-
que se debe señalar el bloque con el cursor y dar doble clic.
13
Capítulo 1. Introducción a las herramientas de Matlab
Las simulaciones para utilizar en control de procesos se pueden dividir en cinco componentes:
• Fuentes (Sources): se encuentran las señales de entrada para modelos más importantes
o herramientas para construirlas.
• Manejo de señales (Signal routing): son las herramientas para procesar información (bus
de datos) entre bloques.
14
Capítulo 1. Introducción a las herramientas de Matlab
En cada uno de los bloques se encuentran los ítems que se deben configurar. Dependiendo de
la señal, se solicita especificar diferentes parámetros y, para ver los parámetros, solo hay que
hacer doble clic en cada bloque. Ejemplos de las propiedades se pueden ver en la figura 1.5.
Los parámetros son:
3. Señal rampa: la función rampa, a diferencia de la señal paso, no tiene un valor final cons-
tante y su incremento es sostenido en el tiempo. Por tal razón, es necesario especificar
la pendiente, el tiempo de arranque de la función y su valor inicial.
5. Señal aleatoria: señal de números aleatorios caracterizados por una media y varianza.
Nota: en algunas de las configuraciones se solicita un tiempo de muestreo. Para sistemas con-
tinuos, se deja el valor de cero, que corresponde a que la señal genera valores en cada tiempo
de integración, el cual es definido por el método numérico. Para sistemas discretos, este deberá
ser definido según el tiempo de muestreo establecido para toda la simulación.
Por último, se presenta una herramienta muy útil para construir señales compuestas de ma-
nera sencilla, que será útil tanto para el análisis de sistemas como para la evaluación de
controladores. Esto puede hacerse con el bloque constructor de señales (Signal Builder). La
interfaz original para manipulación se presenta en la figura 1.6a, donde por defecto una señal
15
Capítulo 1. Introducción a las herramientas de Matlab
pulso rectangular es mostrada. Los ejes por defecto van de 0 a 10. Para incrementar la varia-
ble independiente, se debe ir al menú Axes en la parte superior; posteriormente, se selecciona
Change Time Range y se especifican los límites de tiempo deseados. Después de hacer este
procedimiento, se puede modificar la señal en tiempos mayores a diez unidades de tiempo
(Nota: las unidades de tiempo serán definidas por el modelo matemático que se implemente).
Es posible modificar la señal de dos maneras, después de seleccionar un nodo: se puede sim-
plemente mover el nodo con el mouse y, por otro lado, cuando se selecciona un nodo se puede
ingresar nuevas coordenadas deseadas (Nota: las coordenadas de un nodo en el tiempo (t) no
pueden ser mayores a un nodo con un tiempo superior).
Un par de nodos también puede moverse paralelamente. Con el cursor nos movemos hasta
el segmento de línea que deseamos mover; el cursor cambia a dos líneas paralelas las cuales
pueden moverse exclusivamente de arriba abajo con el clic izquierdo del mouse sostenido. Una
señal modificada se muestra en la figura 1.6b.
Adicionalmente, la forma rectangular de la señal que se tiene por defecto puede modificarse
añadiendo nuevos nodos. Para esto se oprime shift+clic principal sobre el punto donde se desee
añadir el nodo. Finalmente, para agregar otra señal se debe seleccionar Signal del menú prin-
cipal y posteriormente New. Se pueden seleccionar diferentes tipos de señales (triangualares,
constantes, rectangulares) para crear sobre la misma base. Las señales creadas automática-
mente generan una salida adicional en el bloque Signal Builder para ser conectadas a otros
bloques.
16
Capítulo 1. Introducción a las herramientas de Matlab
1.2.2. Sinks
En la sección anterior se presentaron las posibilidades de configuración de bloques para seña-
les, pero no se han usado. Para poder visualizar los resultados se deben emplear los bloques de
sumidero, con los cuales se pueden observar en forma gráfica o numérica los resultados de la
simulación. Los bloques de la librería Sinks solo tienen puertos de entrada porque son los blo-
ques receptores de señales. Las opciones para manejo de resultados se muestran en la figura 1.7.
Entre las opciones para realizar gráficas, se encuentran los bloques Scope y XY Graph. El
bloque Scope grafica variables en función del tiempo y el bloque XY Graph requiere las seña-
les de los ejes de representación de las variables “y” y “x”. Por su parte, los bloques Display y
To Workspace reproducen la información de salida en forma numérica. El primero muestra la
información de manera digital sobre el mismo bloque, el cual permite la selección del formato
17
Capítulo 1. Introducción a las herramientas de Matlab
numérico para el despliegue de la información de salida (Short, long, entre otros). El bloque
de To Workspace envía la información al espacio de trabajo de Matlab utilizando un formato
de estructuras o matrices según se indique en su configuración. Este bloque es vital para la
conectividad entre Simulink y Matlab como se hablará más adelante.
Ejemplo 1.1.
Se desea construir una entrada paso que empiece en el tiempo 5 unidades, con un valor ini-
cial de 1 unidad hasta un valor final de 11 unidades (amplitud del paso, 10 unidades). Para
visualizar el resultado, se desea generar una gráfica en el mismo ambiente de Simulink.
Solución
Para realizar esta implementación, se describirá paso a paso la estrategia pues es la primera.
Un resumen de cómo debe lucir el progreso del trabajo se muestra en la figura 1.8.
18
Capítulo 1. Introducción a las herramientas de Matlab
• Hay que conectar los bloques, notar nuevamente que el bloque paso solo tiene una salida
y el bloque de graficar solo tiene una entrada. Para conectarlos simplemente arrastre
la flecha desde el bloque paso hasta el Scope. En las versiones más modernas, cuando
se colocan los bloques alineados, el mismo Simulink muestra en sombra que estos dos
bloques pueden conectarse. Ahora, si se siente profesional, puede usar un truco: señale
el primer bloque y oprima la tecla ctrl, sin soltarla, señale el bloque de Scope, conexión
rápida.
• ¡Sorpresa, no pasó nada! Bueno, para ver la respuesta se debe hacer doble clic sobre el
bloque Scope.
• FIN.
Algunos comentarios acerca de la solución: por defecto, la simulación se llevó a cabo para 10
unidades de tiempo (tiempo de simulación por defecto). Esto puede modificarse en el menú
principal de la ventana de trabajo donde aparece el número 10 en un recuadro llamado Si-
mulation stop time. Adicionalmente, Simulink usa un formato poco conveniente para mostrar
las figuras con un fondo negro, línea amarilla y ejes sin formato. Para hacer seguimiento a
simulaciones no es tan relevante; sin embargo, es imperativo cambiar ese formato.
19
Capítulo 1. Introducción a las herramientas de Matlab
Figura 1.8. Simulación de la señal paso, que muestra la configuración del bloque y la gráfica obtenida
De manera análoga, se puede visualizar cualquiera de las señales dispuestas en el menú Sour-
ces, por ejemplo, Ramp, Sine wave, Random number, Signal builder, etc.
Debido a que el formato de las figuras por defecto es muy malo para impresión, por el fondo
negro, se mostrará qué se puede hacer para mejorarlo. Para esto, se usa un botón similar a un
engranaje llamado Configuration properties, que se encuentra disponible también en el menú
View. Algunos parámetros importantes por configurar para las gráficas son:
• Número de puertos de entrada: si se desea mostrar más de una figura por ventana, se
sugiere utilizar el botón de Layout para elegir una distribución diferente.
• Mostrar leyenda: indica qué color corresponde a cada señal. Esto es útil cuando se usa
más de una señal en la figura.
• y-label: adición en la última versión de Matlab para poder poner el nombre del eje y.
En versiones anteriores de Simulink, con el mismo botón se puede configurar el formato de es-
tilo. Para la versión en cuestión, la configuración de estilo se encuentra en el botón View/Style.
20
Capítulo 1. Introducción a las herramientas de Matlab
Ejemplo 1.2.
Se desea cambiar el formato de la figura del ejemplo anterior por uno más parecido al de las
figuras de Matlab: fondo blanco, línea de color azul y más gruesa, nombre en los ejes y leyenda.
Solución
21
Capítulo 1. Introducción a las herramientas de Matlab
• Mux: permite combinar señales en un bus de datos común, en la configuración del bloque
se especifican cuántas señales de entrada se quieren combinar.
• Manual Switch: permite escoger, entre dos señales, cuál atraviesa el bloque.
• Switch: permite utilizar un criterio para escoger, entre señales, cuál atraviesa el bloque.
Cuando se hacen operaciones matemáticas entre buses de datos, normalmente no cambian las
dimensiones del vector.
Ejemplo 1.3.
Se desea combinar en un solo bus de datos la señal paso diseñada anteriormente con una onda
sinusoidal de amplitud cinco y BIAS cinco para ser graficadas enseguida sobre la misma fi-
gura. Adicionalmente, se desea separar de nuevo las señales y graficarlas en la misma ventana
de Scope, pero en diferente figura.
Solución
Para realizar la mezcla de buses de datos y su posterior separación se usarán las herramientas
Mux y Demux, respectivamente. Teniendo en cuenta que son dos señales de entrada para el
Mux y dos de salida para Demux; estas condiciones deben ser establecidas en la configuración
de cada bloque. Utilizando los conceptos anteriores para ubicar bloques y unirlos se construye
el diagrama mostrado en la parte izquierda de la figura 1.10. Nota: algo nuevo es la división
del bus de datos que sale del Mux para ser enviado al Demux. Para hacer eso se realiza la
unión de la línea con el Demux usando el botón secundario del mouse.
22
Capítulo 1. Introducción a las herramientas de Matlab
Más adelante, será útil poder extraer de un bus de datos solo un par de señales. Por ejemplo,
en el bus de datos vienen quince señales, pero solo se desean graficar dos. Entonces, el bloque
Selector es interesante para realizar la tarea.
• Sum y Add: se utiliza para sumar o restar señales, el número de señales para trabajar
y su signo se fijan en la configuración del bloque. El bloque Add permite manejar un
mayor número de señales.
• Substract: aunque el bloque anterior permite sumar y restar, Simulink ofrece este bloque
especializado para restar señales.
Estos bloques están caracterizados por tener entradas y salidas. En la mayoría de ellos el
número de variables no cambia al realizarse la operación aritmética.
23
Capítulo 1. Introducción a las herramientas de Matlab
Ejemplo 1.4.
Operaciones matemáticas
Se desea construir una señal onda cuadrada utilizando señales paso (no el Signal bouilder) y
posteriormente amplificar esa señal. La onda cuadrada empieza en el tiempo t = 5, pasando
de un valor inicial de 1 a 11 unidades. Una unidad de tiempo más adelante, la señal debe
regresar al valor original.
Solución
Se utiliza la plantilla anterior con la señal paso implementada. Para realizar la señal onda
cuadrada se utiliza la resta de dos señales paso trasladadas en el tiempo. La primera señal
es la que se diseñó con anterioridad y la segunda es una señal que se empieza en el tiempo 6
con valor inicial de 0 y una amplitud de 10 unidades. Cuando se restan estas dos señales se
genera la onda cuadrada. Para amplificar la señal se usa el bloque Gain con un valor de 2, lo
que significa que la señal de entrada se amplifica por 2. La implementación y los resultados
se muestran en la figura 1.11.
1.2.5. Modelos
La implementación de modelos se tratará más adelante; por ahora, queda implícito que se
requieren bloques especiales para implementar modelos en Simulink, dependiendo de la na-
turaleza del modelo (continuo, discontinuo, lineal, no lineal, en el dominio del tiempo, en el
dominio de Laplace, en el dominio de la frecuencia, entre otras posibilidades).
24
Capítulo 1. Introducción a las herramientas de Matlab
25
Capítulo 1. Introducción a las herramientas de Matlab
26
Capítulo 2
Resumen
Para abordar la cotidianidad en Ingeniería Química, es relevante notar que las no linea-
lidades aparecen constantemente en nuestras formulaciones matemáticas, por ejemplo:
velocidades de reacción (e.g., el término preexponencial de Arrhenius, cinéticas de al-
to orden, cinéticas de microorganismos, cinéticas enzimáticas, cinéticas no elementales),
comportamiento del pH (curva de titración), relaciones de equilibrio de fases, modelos
de actividad (efecto en destilación, extracción líquido-líquido, cristalización), ecuaciones
de estado, equilibrio químico, caída de presión en lechos empacados (reactores de flujo
piston, filtros), entre otros. Teniendo en cuenta que el curso de Control Automático de
Procesos se fundamenta en la aplicación de aproximaciones basadas en modelos para di-
señar y evaluar estructuras de control, el primer paso es contar con un entendimiento del
sistema a través del modelo soportada por evidencia experimental. En este capítulo se
utilizarán herramientas tanto de Matlab como de Simulink para la solución de sistemas de
ecuaciones diferenciales ordinarias correspondientes a problemas de valor inicial (PVI).
Para la solución se utilizarán funciones anónimas, funciones separadas, interpreted fun-
ctions y s-functions. Todas estas implementaciones son equivalentes; sin embargo, una u
otra es más conveniente dependiendo del tamaño de la aplicación. Por último, se hacen
algunos comentarios relacionados con la configuración de la solución.
2.1. Introducción
El punto de partida para analizar el comportamiento de sistemas dinámicos es poder resol-
ver las ecuaciones que los definen. En Ingeniería Química, los balances de materia, energía
y momento, pueden llevar a diferentes tipos de ecuaciones diferenciales. Cuando se tienen
sistemas de parámetros concentrados (i.e., sistemas donde los estados no dependen de la lo-
calización, como en sistemas de tanque agitado), las ecuaciones de conservación dinámicas
llevan a sistemas de ecuaciones diferenciales ordinarias con variable independiente el tiempo.
Estas ecuaciones diferenciales ordinarias son problemas de valor inicial y pueden ser resueltos
a través de métodos de paso disponibles en Matlab. Por otro lado, cuando se tienen sistemas
de parámetros distribuidos (i.e., sistemas donde los estados dependen de la posición, como en
27
Capítulo 2. Simulación de sistemas no lineales continuos
En este capítulo se analizan inicialmente los solvers con los que cuenta Matlab y se pro-
fundiza en el problema que implica tener ecuaciones diferenciales rígidas y opciones para
manejarlas. Enseguida, se presentan diferentes formas para resolver ecuaciones diferenciales
ordinarias desde el ambiente de Matlab a través de funciones anónimas y la definición de fun-
ciones (function). Adicionalmente, se presenta cómo resolver las mismas ecuaciones desde
el ambiente de Simulink utilizando los bloques definidos por el usuario llamados Interpreted
function y s-function. Las implementaciones en Simulink son útiles para usar configuraciones
especiales de los métodos y explotar interconexiones con Matlab para obtener figuras con un
buen formato.
Matlab-Simulink cuenta con una serie de funciones que se utilizan para resolver problemas de
valor inicial tanto de ecuaciones diferenciales ordinarias (EDO o ODE por sus siglas en inglés)
como de ecuaciones diferenciales algebraicas (DAE por sus siglas en inglés). La selección del
método por utilizar depende de las características de las ecuaciones que se deben resolver. En
la tabla 2.1 se muestra un resumen de los métodos y cuándo usarlos.
La selección del solver es muy importante y requiere un conocimiento del tipo de ecuaciones
que se van a manejar. Por defecto, se recomienda utilizar el ode45 como primera opción; sin
embargo, este método de paso variable tiene problemas con sistemas rígidos y toma mucho
tiempo para resolverlos o incluso podría no converger.
El problema radica en el algoritmo que se utiliza para determinar el tamaño de paso del
método. Los métodos numéricos adaptativos utilizan un tamaño de paso variable para op-
timizar el número de pasos en la solución de ecuaciones diferenciales; si un paso constante
es usado, esto implica malgastar energía en regiones con cambios graduales. Cuando se usa
un tamaño de paso variable, se toman tamaños de paso grandes si la derivada de la función
tiende a ser constante y se reduce el tamaño cuando hay cambios abruptos de pendiente.
Para lidiar con los cambios de pendiente se utiliza criterios del error local de truncamiento
en cada paso. Entonces este estimado decide si se aumenta o disminuye el paso. Por ejemplo,
el método Step-Halving usa el error local usando un orden de truncamiento y dos tamaños
de paso, o el método Runge-Kutta Fehlberg compara los resultados del método Runge-Kutta
RK-5 y RK-4 para determinar el nuevo tamaño de paso (Dato curioso: la razón del nombre
del ode45 es porque utiliza el método Runge-Kutta de 4.o y 5.o orden para comparar el error
28
Capítulo 2. Simulación de sistemas no lineales continuos
Tabla 2.1. Funciones de Matlab para la solución de ecuaciones diferenciales ordinarias y ecuaciones
diferenciales algebraicas
de truncamiento; entonces no se debe referir al método como “ode cuarenta y cinco” sino “ode
cuatro-cinco”. Esto aplica para los otros métodos igualmente).
La rigidez de ecuaciones diferenciales es una propiedad casi cualitativa, aunque se han pro-
puesto maneras de cuantificarlo como el índice de rigidez o razón de rigidez y ocurre cuando
se tienen fenómenos en diferentes escalas de tiempo. Esto implica que algunos estados del
sistema reaccionan muy rápido ante perturbaciones y otros no, lo cual es evidente en las cons-
tantes que acompañan los términos en la ecuación diferencial. Cuando un método numérico
convencional se enfrenta al problema, los fenómenos rápidos dominan la selección del tamaño
de paso forzando al método a realizar un paso muy pequeño, incluso para regiones donde
la función es suave, haciendo la simulación muy ineficiente. Esto ocurre pues al integrar el
sistema de ecuaciones con un método de orden n y una tolerancia de error local de 10−n ,
el paso de integración del algoritmo debe hacerse más pequeño que el valor indicado por la
estimación del error local debido a las restricciones impuestas por la región de estabilidad
numérica (Cellier y Kofman, 2006; Migoni y Kofman, 2007).
Por ende, estos sistemas requieren la utilización de algoritmos de orden de integración menor
o algoritmos especializados para sistemas rígidos. Estos son más eficientes en la determinación
del tamaño de paso y, por tanto, los tiempos de simulación son menores. Aunque esta solución
es la recomendada y funciona casi siempre para los casos evaluados en el curso de control, aún
hay retos para sistemas más complejos. Un problema común de los algoritmos para sistemas
29
Capítulo 2. Simulación de sistemas no lineales continuos
Para crear una función anónima, se utiliza el símbolo “@”. Para la ecuación diferencial del
reactor se usa:
dCA F 2
= (CAf − CA ) − k1 CA − k3 CA (2.1)
dt V
El código en Matlab que generaría la primera ecuación diferencial y la ejecuta es:
%--------------------------------------------------------------------------
k1 = 5/6; % Constante de reacción para A-->B (min^-1)
k2 = 5/3; % Constante de reacción para B-->C (min^-1)
k3 = 1/6; % Constante de reacción para 2A-->D (mol/(l min))
D = 4/7; % Tasa de dilución (min^-1)
Caf = 10; % Concentración de A a la entrada (mol/l)
dCa_dt = @(t,C)D*(Caf-C(1)) - k1*C(1) -k3*C(1)^2; % Función anónima
30
Capítulo 2. Simulación de sistemas no lineales continuos
Donde:
Para constituir un sistema de ecuaciones, se puede crear cada ecuación de manera indepen-
diente y luego agruparlas en una nueva variable o simplemente crearlas como un vector desde
el principio.
La estructura para resolver sistemas de ecuaciones es la misma para cualquier solver por
utilizar. Como ejemplo se emplea el método ode15s, que resuelve sistemas de ecuaciones
diferenciales originarias rígidas y sistemas de ecuaciones diferenciales algebraicas. Se utiliza
la variable options para configurar las opciones de la función. La estructura simplificada con
opciones de tolerancias del método de manera explícita es:
%--------------------------------------------------------------------------
options = odeset('AbsTol',1e-6,'RelTol',1e-3);
[t,C] = ode15s(odefun,tspan,y0,options)
%--------------------------------------------------------------------------
Donde:
31
Capítulo 2. Simulación de sistemas no lineales continuos
Ejemplo 2.1.
Se requiere resolver las ecuaciones diferenciales ordinarias para es sistema reactivo de Van
de Vusse en un reactor continuo de mezcla completa (CSTR) utilizando funciones anónimas.
Los parámetros, variables de entrada y condiciones iniciales se muestran en la tabla 6.1.
Solución
Para esta primera solución, se va a explotar funciones anónimas en Matlab para incluir las
ecuaciones dentro del mismo programa. Utilizar funciones anónimas de la manera en que se
muestra es útil cuando se tienen pocas ecuaciones diferenciales. El programa se divide en
tres secciones: inicialización (i.e., parámetros, variables de entrada y condiciones iniciales, de
acuerdo con la tabla 6.1), solución del modelo (i.e., planteamiento de las ecuaciones, configu-
ración del solver y su ejecución) y presentación de resultados (i.e., gráficas con su respectivo
formato).
Como aspectos importantes en esta implementación, se tiene la definición del tiempo de si-
mulación como un vector especial. Se utiliza la función linspace con el fin de generar un
vector de tiempo de 100 elementos entre el valor inicial y el tiempo total de simulación. Esto
no quiere decir que el algoritmo no use un paso variable durante su ejecución; solamente que
retorna los valores de los estados en los instantes de tiempo especificados por la función usada.
Si se desea ver los tiempos usados por el solver, entonces se debe definir el vector de tiem-
po como tspan = [0 tsim], donde tsim corresponde al valor final del tiempo de integración.
32
Capítulo 2. Simulación de sistemas no lineales continuos
%--------------------------------------------------------------------------
% Solución del sistema Reactor VDV
% Sistema de ecuaciones en el mismo script-función anónima
clc
clear
close all
%% Solución
%% Resultados
figure(1)
subplot(211)
plot(t,C(:,1),'b','LineWidth',2); grid;
ylabel('C_A[mol/l]','FontWeight','bold','FontSize',14)
set(gca,'FontSize',12,'FontWeight','bold')
subplot(212);grid
plot(t,C(:,2),'b','LineWidth',2); grid;
ylabel('C_B[mol/l]','FontWeight','bold','FontSize',14)
xlabel('Tiempo [min]','FontWeight','bold','FontSize',14)
set(gca,'FontSize',12,'FontWeight','bold')
%--------------------------------------------------------------------------
33
Capítulo 2. Simulación de sistemas no lineales continuos
Figura 2.1. Simulación del reactor Van de Vusse usando los valores mostrados en tabla 6.1
Nótese que se recomienda un formato de figuras que hacen que luzcan bien cuando se exporten
a documentos de texto. Este formato incluye: formato de líneas, tamaño de fuente y formato
de fuente en negrita.
Como alternativa al empleo de funciones anónimas para representar las ecuaciones diferencia-
les, se puede plantear un ambiente diferente al programa principal. Este ambiente corresponde
a un programa diferente diseñado para proveer unas salidas cuando se alimentan detarminadas
entradas. Una función por lo general se ubica en un archivo diferente con nombre exactamente
igual al de la función. Como el objetivo de los programas por realizar es solo tener un M-file,
se utilizará una sintaxis especial que permite tener tanto el programa maestro como la función
que contiene las ecuaciones diferenciales en el mismo archivo.
34
Capítulo 2. Simulación de sistemas no lineales continuos
%--------------------------------------------------------------------------
function ['Salidas']=Nombre('Entradas')
Algoritmo % Conjunto de operaciones por realizar en la función
Salidas = ...
end
%--------------------------------------------------------------------------
Donde “Salidas” es un conjunto de variables separadas por coma que se calculan dentro de la
función, “Nombre” es la referencia para la función (i.e., nombre que es usado para ejecutarla)
y “Entradas” es un conjunto de variables/parámetros separados por comas que son necesarios
para hacer los cálculos dentro de la función. Algo particularmente importante de las funciones
es cómo manejan la información. Desde el punto de vista del workspace, para una función solo
existen las variables/parámetros que se especificaron como entradas, a pesar de que existan en
el workspace durante la ejecución del programa. Análogamente, una vez se termina la ejecución
de la función, para el programa principal solo existen las variables/parámetros que se especi-
ficaron como salidas de las función (Nota: aunque esta es una condición para la ejecución de
funciones, existen maneras para evitar esta condición como las funciones evalin y assignin).
Ejemplo 2.2.
Se requiere resolver las ecuaciones diferenciales ordinarias para la reacción de Van de Vusse
en un reactor CSTR utilizando funciones. Los parámetros, variables de entrada y condiciones
iniciales se muestran en la tabla 6.1.
Solución
Se usan funciones con el fin de especificar las ecuaciones diferenciales, pero esto se hace dentro
de un mismo programa. Utilizar funciones es útil cuando se tienen muchas ecuaciones diferen-
ciales, acompañadas con ecuaciones complementarias como cinéticas y muchos parámetros.
Para lograr el uso de funciones en el mismo archivo que el programa que las llama, el programa
maestro debe ser también una función.
35
Capítulo 2. Simulación de sistemas no lineales continuos
Hint: algo importante en esta implementación es que el nombre del archivo debe ser
igual al de la función principal, en este caso VDVReactor_2. Hay que tener cuidado
con el con el nombre que se selecciona, pues no puede ser igual a uno de una función
existente en Matlab.
%--------------------------------------------------------------------------
% Solución del sistema Reactor VDV
% Sistema de ecuaciones en el mismo script-función
function VDVReactor_2
clc
clear all
close all
%% Solución
%% Resultados
figure(1)
subplot(211)
plot(t,C(:,1),'b','LineWidth',2); grid;
ylabel('C_A[mol/l]','FontWeight','bold','FontSize',14)
set(gca,'FontSize',12,'FontWeight','bold')
subplot(212);grid
plot(t,C(:,2),'b','LineWidth',2); grid;
ylabel('C_B[mol/l]','FontWeight','bold','FontSize',14)
xlabel('Tiempo [min]','FontWeight','bold','FontSize',14)
set(gca,'FontSize',12,'FontWeight','bold')
end
dC_dt = zeros(length(C),1);
36
Capítulo 2. Simulación de sistemas no lineales continuos
end
%--------------------------------------------------------------------------
Algo importante por notar en la implementación, es que los parámetros cinéticos y las va-
riables de entrada se especificaron como entradas a la función ODE_equations. Por tanto,
se definen antes de ejecutar el método numérico. Esta implementación es conveniente, pues
permite tener solo un programa para realizar los cálculos; además es útil para análisis de
sensibilidad, estimación de parámetros, optimización, control de procesos, entre otras aplica-
ciones.
Al ejecutar el código propuesto se obtiene como resultado la figura 2.2, idéntica a la figura
2.1 del ejemplo anterior.
Figura 2.2. Simulación del reactor Van de Vusse usando los valores mostrados en tabla 6.1
Nota: esta implementación tiene una desventaja. Como el programa principal es una función
que no tiene entradas ni salidas, al terminar la ejecución del programa ninguno de los resulta-
dos queda grabado en el workspace. Para subsanar este inconveniente, se pueden guardar los
resultados en un archivo de datos .mat usando el comando save dentro del programa principal.
37
Capítulo 2. Simulación de sistemas no lineales continuos
En los ejemplos 2.1 y 2.2 se mostraron dos opciones de configuración para el solver de las ecua-
ciones diferenciales ordinarias usando la función odeset. Dentro de esta función se pueden
configurar 21 propiedades de los algoritmos, útiles en diferentes aplicaciones; sin embargo, se
mostró solo la configuración por defecto de dos de ellas llamadas AbsTol y RelTol. Teniendo
en cuenta la relevancia de sintonizar la exactitud de los métodos, estas tolerancias pueden
disminuirse. Casos en los que se recomienda modificar las tolerancias son, por ejemplo: cuan-
do en la solución se ve ruido numérico, cuando se tienen concentraciones negativas cercanas
a cero que violan balances de materia o cuando la solución no se ve suavizada.
Finalmente, hay muchas opciones para hacer implementaciones más genéricas útiles para sis-
temas más grandes. A este nivel se puede:
• Usar variables llamadas “estructuras” para transportar información entre funciones. Una
estructura puede almacenar muchos parámetros y es posible transportar de una función
a otra solo llamando el nombre la variable principal. Por ejemplo, en la variable llamada
par, podemos almacenar los parámetros D = 4/7, Caf = 10, Cao = 2 y Cbo = 1,117,
así:
%--------------------------------------------------------------------------
par.D = 4/7; % Tasa de dilución (min^-1)
par.Caf = 10; % Concentración de A a la entrada (mol/l)
par.Cao = 2; % Concentración inicial de A (mol/l)
par.Cbo = 1.117; % Concentración inicial de B (mol/l)
%--------------------------------------------------------------------------
Donde par puede ser la entrada a una función y dentro de ella el algoritmo reconoce
los parámetros independientes para usar en las ecuaciones, entonces tenemos los cuatro
parámetros incluidos. Esto difiere al usar un vector par = [4/7 10 2 1.117], pues en
este último caso tenemos que saber que el primer elemento del vector es D, el segundo
es Caf , el tercero es Cao y el cuarto es Cbo . Usando estructuras, le ponemos el nombre
de la variable y no importa el orden en el que las definimos, pues siempre existirán como
parámetros dentro de la estructura. Esto es particularmente importante cuando tene-
mos modelos muy grandes con una cantidad de parámetros elevada o cuando estamos
desarrollando un modelo en el que continuamente estamos modificando el programa,
por lo que memorizar la posición de parámetros en un vector se torna una pesadilla.
38
Capítulo 2. Simulación de sistemas no lineales continuos
El problema con esta configuración es que para cada tasa de dilusión evaluada los vec-
tores de Ca y Cb pueden tener diferente tamaño, lo cual puede entrar en conflicto con
la matriz original de tamaño An,m,l . Como alternativa aparecen las variables llamadas
células. Una célula es un contenedor análogo a una matriz, pero en cada posición de
esta se puede almacenar cualquier tipo de información sin importar su naturaleza (e.g.,
vector, matriz, letras, archivos de texto) y sin importar su dimensión. La sintaxis de las
células es muy parecida a las de matrices para referirse a las posiciones, solo cambiando
el tipo de paréntesis usado. Estas pueden usarse así:
%--------------------------------------------------------------------------
A = [1 2 3; 4 5 6]; % Matriz 2x3 de ejemplo
B = [7 8; 9 10]; % Matriz 2x2 de ejemplo
R_cell{1} = A;
R_cell{2} = B;
%--------------------------------------------------------------------------
Para resolver ecuaciones diferenciales ordinarias en Simulink se utilizan los mismos métodos
numéricos que en Matlab, pues el motor de solución es igual. Sin embargo, la forma en que
se concibe el programa es parcialmente diferente. Para realizar las implementaciones para la
solución de las EDO en Simulink como ambiente principal de programación, se explotan ar-
quitecturas de programación análogas a la mostrada en la figura 1.12b. Entonces, el programa
principal es en Simulink y este llama programas específicos en Matlab (M-file) para realizar
tareas particulares. En esta sección se utilizarán dos funciones diferentes que permiten realizar
la simulación de sistemas representados por ecuaciones diferenciales ordinarias. Adicionalmen-
te, se emplean las herramientas de conectividad entre ambientes de trabajo para poder cargar
información a Simulink y posterior a la simulación realizar las gráficas con formato en Matlab.
39
Capítulo 2. Simulación de sistemas no lineales continuos
40
Capítulo 2. Simulación de sistemas no lineales continuos
entrada y uno de salida, lo que implica que todas las entradas se deben tener como una
serie de tiempo y la salida será otra serie de tiempo al terminar la simulación. Para las
ecuaciones diferenciales, el bus de entrada debe contener las entradas al sistema y los
estados (que son retroalimentados), con el fin de que la salida sea un bus con los valores
de dCi /dt. Este vector será integrado con otro bloque especial llamado Integrator. Este
bloque utiliza una aproximación numérica para evaluar la integral de la señal de entrada
en el tiempo actual, utilizando el valor de la entrada en el tiempo actual y los estados
en el tiempo anterior. Para este bloque hay que especificar las condiciones iniciales del
sistema (MatWorks, 2022i). La configuración por defecto de los bloques utilizados se
muestran en la figura 2.5. Para el caso de estudio, la configuración será modificada.
Figura 2.5. Configuración por defecto de los bloques Interpreted function e Integrator
41
Capítulo 2. Simulación de sistemas no lineales continuos
3. Programa en Matlab para realizar las figuras: una vez se terminan los cálculos en Si-
mulink, los datos de solución de las ecuaciones diferenciales son enviados al workspace
utilizando un bloque especial de Simulink llamado To workspace. Aprovechando que
los datos de la simulación existen en el workspace, se construye un programa con un
M-file que grafique los resultados con un formato ajustado para que se vean bien al ser
exportadas para un reporte.
Estas ideas son empleadas en el ejemplo 2.3 para la solución del sistema de ecuaciones dife-
renciales ordinarias del reactor Van de Vusse.
Ejemplo 2.3.
Se desea evaluar el comportamiento de los estados del sistema reactivo Van de Vusse en un
CSTR cuando se aplican perturbaciones tanto en la tasa de dilución como en la concentración
de entrada de A. Los valores de parámetros, condiciones iniciales y variables de entrada son
tomados de la tabla 6.1. La simulación se realiza por 30 min usando una Interpreted function,
donde en el tiempo=10 min se perturba la tasa de dilución con un cambio escalón de amplitud
de 1/5 min−1 . En el tiempo=20 min, se perturba la concentración de A en el alimento con
un cambio escalón de amplitud -2 mol/l. Finalmente, se desea graficar tanto las señales de
entrada como de salida del modelo en Matlab de manera automática.
Solución
Para resolver este problema, se utiliza la estrategia descrita anteriormente. Primero se realiza
la implementación en Simulink, seguido de la construcción del modelo matemático en Matlab
para, por último, realizar de manera automática las gráficas usando otro programa de Matlab.
En primera instancia, se utilizan los bloques mostrados en la tabla 2.2 para construir el mo-
delo mostrado en la figura 2.6.
42
Capítulo 2. Simulación de sistemas no lineales continuos
Figura 2.6. Implementación en Simulink del sistema Van de Vusse usando Interpreted fuction
Tabla 2.2. Tipo de variable asociada a cada bloque y su ubicación en la librería de Simulink
Variable Ubicación Bloque
Tasa de dilución Sources Step
Concentración de entrada Sources Step
Modelo User-defined functions Interpreted function
Manejo de señales Signal Routing Mux
Resultados Simulink Sinks Scope
Resultados Matlab Sinks To workspace
Conectividad Menú Model properties
Fuente: Elaboración propia
Por otro lado, nótese que se usaron unos bloques de conectividad llamados To workspace.
Este bloque permite enviar al workspace la información que se le alimenta en diferente tipo
de formato. Esto se hace para que se pueda construir un programa en Matlab que permita
graficar los resultados con un formato análogo al utilizado en implentaciones previas.
Hint: como dato curioso, los bloques de Simulink permiten asociar imágenes para
mejorar la visualización. Esto se puede hacer en el menú Mask. En este caso se puso
un reactor decorativo, que será usado de ahora en adelante en las simulaciones.
43
Capítulo 2. Simulación de sistemas no lineales continuos
Ahora, se construye en Matlab la función con las ecuaciones diferenciales, de la misma manera
que se hizo en la sección 2.3. El programa es:
%--------------------------------------------------------------------------
Ca = u(3); % Concentración de A
Cb = u(4); % Concentración de B
dC_dt = zeros(length(2),1);
44
Capítulo 2. Simulación de sistemas no lineales continuos
end
%--------------------------------------------------------------------------
Nótese que en el programa el orden de las entradas, donde las primeras dos corresponden a
las entradas en el mismo orden que fueron puestas en Simulink, seguidas de los estados en el
mismo orden en que serán puestas las ecuaciones diferenciales.
Por último, se construye un programa en Matlab para realizar las gráficas deseadas utilizando
los valores exportados por los bloques To workspace. Es importante tener en cuenta el for-
mato de datos usado por el workspace es series de tiempo. Por tanto, dentro de las variables
exportadas están tanto el tiempo como los datos (i.e., entradas para “Inputs” y estados pa-
ra “yout_NL”. Las series de tiempo usan la nomenclatura de variables tipo estructura, que
guardan variables dentro de variables (como ya se mencionó). Para facilitar el manejo de las
variables, en la primera parte del programa se extraen las variables a un formato convencional
de vector. Todas la gráficas se ubican en una sola figura usando el comando subplot.
Ahora bien, se requiere un truco final para correr automáticamente el programa de las figu-
ras cuando se termine la simulación de Simulink; aunque esto no es necesario, incorporarlo
facilita mucho la simulación. Para esto, se abre el menú Model properties que se encuentra en
la pestaña de Modelling/Model settings. Se abrirá la ventana mostrada en la figura 2.8 y se
selecciona la opción Callbacks. Esta funcionalidad de Simulink permite correr programas de
Matlab en diferentes momentos de la simulación: antes de empezar, cuando termina, cuando
se corre, cuando se cierra, entre otras. La opción que se usa es la StopFcn, que indica qué hacer
cuando el programa de Simulink termina la simulación. En la ventana de la derecha se escribe
el nombre del programa de Matlab que hace las figuras. De esta manera, cuando Simulink
termina la ejecución del programa, automáticamente se correrá el programa de Matlab que
hace las figuras, que en este caso se llama “Plots_NL”.
%--------------------------------------------------------------------------
% Programa para hacer la figuras
close all
% Entradas
45
Capítulo 2. Simulación de sistemas no lineales continuos
% Resultados
figure(1)
subplot(221)
plot(tiempo_In,D,'b','LineWidth',2)
ylabel('D[1/min]','FontWeight','bold','FontSize',14)
set(gca,'FontSize',12,'FontWeight','bold')
grid
subplot(222)
hold on
plot(tiempo,Ca,'b','LineWidth',2)
grid
ylabel('C_A[mol/l]','FontWeight','bold','FontSize',14)
set(gca,'FontSize',12,'FontWeight','bold')
subplot(223)
plot(tiempo_In,Caf,'b','LineWidth',2)
ylabel('Caf[mol/l]','FontWeight','bold','FontSize',14)
xlabel('Tiempo [min]','FontWeight','bold','FontSize',14)
set(gca,'FontSize',12,'FontWeight','bold')
grid
subplot(224)
plot(tiempo,Cb,'b','LineWidth',2)
grid
ylabel('C_B[mol/l]','FontWeight','bold','FontSize',14)
xlabel('Tiempo [min]','FontWeight','bold','FontSize',14)
set(gca,'FontSize',12,'FontWeight','bold')
%--------------------------------------------------------------------------
46
Capítulo 2. Simulación de sistemas no lineales continuos
Para obtener los resultados, se ejecuta el programa de Simulink desde la ventana de trabajo
y automáticamente se deben generar las siguientes gráficas (figura 2.9):
Figura 2.9. Resultados de la simulación del sistema Van de Vusse ante dos perturbaciones
47
Capítulo 2. Simulación de sistemas no lineales continuos
Figura 2.10. Estructura de un modelo en Simulink usando s-function y las ecuaciones que lo
describen. [Fuente: elaboración propia con base en MatWorks (2022y)]
Para resolver dicho modelo híbrido (i.e., continuo, discreto o mezcla), Simulink ejecuta un
algoritmo jerárquico en forma de script de Matlab en etapas (MatWorks, 2022y), las cuales
se ilustran en la figura 2.11. El proceso comienza con la inicialización. En esta etapa se es-
tablecen: tamaño de las señales, tipo de datos, información acerca del tiempo de muestreo;
se evalúan los bloques de parámetros, se determina el orden de ejecución del programa y se
fija la memoria. Este proceso se realiza una sola vez, para posteriormente entrar al ciclo de
simulación donde cada interacción se refiere a “paso de simulación”. En el ciclo, Simulink
ejecuta cada parte del programa siguiendo el orden que se definió en la inicialización a través
de funciones que son ejecutadas usando una variable llamada Flag. Estas funciones calculan
los estados, las derivadas y las salidas en cada tiempo de muestreo.
Si el sistema contiene estados continuos, hay un ciclo interno llamado “integración” donde se
resuelven las ecuaciones diferenciales. Este ciclo se repite hasta que se alcanza la exactitud
deseada para los cálculos de los estados continuos.
Para usar una s-function, se realiza una implementación análoga a la del caso anterior cuan-
do se empleó la Interpreted function. Entonces, el programa consta de una implementación en
Simulink y una en M-file. Como una extensión del caso anterior, se plantea una arquitectura
específica que se muestra en la figura 2.12. La arquitectura de la implementación es muy
parecida al caso anterior, con la diferencia de que no se usa el bloque de Integrator y no hay
retroalimentación de los estados. Para la configuración del bloque s-function, solo requie-
re especificar el nombre del la s-function guardada como una función M-file. El bloque de
esta función se encuentra en la librería de Simulink en la subsección de User defined functions.
48
Capítulo 2. Simulación de sistemas no lineales continuos
Figura 2.11. Ciclo de simulación de una s-function en Simulink. [Fuente: elaboración propia con
base en MatWorks (2022y)]
49
Capítulo 2. Simulación de sistemas no lineales continuos
Figura 2.12. Arquitectura genérica para implementar la solución de ecuaciones diferenciales en Si-
mulink a través de la función s-function
50
Capítulo 2. Simulación de sistemas no lineales continuos
51
Capítulo 2. Simulación de sistemas no lineales continuos
%
% The following outlines the general structure of an S-function.
2. Orden de Flag: en esta sección se muestra la serie de funciones que se ejecutan depen-
diendo del valor de la variable flag.
switch flag,
%%%%%%%%%%%%%%%%%%
% Initialization %
%%%%%%%%%%%%%%%%%%
case 0,
[sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes;
%%%%%%%%%%%%%%%
% Derivatives %
%%%%%%%%%%%%%%%
case 1,
sys=mdlDerivatives(t,x,u);
%%%%%%%%%%
% Update %
%%%%%%%%%%
case 2,
sys=mdlUpdate(t,x,u);
%%%%%%%%%%%
% Outputs %
%%%%%%%%%%%
case 3,
sys=mdlOutputs(t,x,u);
%%%%%%%%%%%%%%%%%%%%%%%
% GetTimeOfNextVarHit %
%%%%%%%%%%%%%%%%%%%%%%%
case 4,
52
Capítulo 2. Simulación de sistemas no lineales continuos
sys=mdlGetTimeOfNextVarHit(t,x,u);
%%%%%%%%%%%%%
% Terminate %
%%%%%%%%%%%%%
case 9,
sys=mdlTerminate(t,x,u);
%%%%%%%%%%%%%%%%%%%%
% Unexpected flags %
%%%%%%%%%%%%%%%%%%%%
otherwise
DAStudio.error('Simulink:blocks:unhandledFlag', num2str(flag));
end
En esta parte del programa solo se llaman a subsecciones; se puede ver que en cada caso
no hay código para ejecutar acciones, solo se llaman las funciones que son definidas más
adelante. Esta parte no debe ser modificada, aunque para usuarios más avanzados, esta
parte del programa puede ser extendida para incluir nuevas variables. A continuación
se muestra el sentido de cada una de las funciones definidas por los Flags.
%=============================================================================
53
Capítulo 2. Simulación de sistemas no lineales continuos
% mdlInitializeSizes
% Return the sizes, initial conditions, and sample times for the S-function.
%=============================================================================
%
function [sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes
%
% call simsizes for a sizes structure, fill it in and convert it to a
% sizes array.
%
% Note that in this example, the values are hard coded. This is not a
% recommended practice as the characteristics of the block are typically
% defined by the S-function parameters.
%
sizes = simsizes;
sys = simsizes(sizes);
%
% initialize the initial conditions
%
x0 = []; % Vector con las condiciones iniciales
%
% str is always an empty matrix
%
str = [];
%
% initialize the array of sample times
%
ts = [0 0]; % Configuración de los tiempos de muestreo
% end mdlInitializeSizes
%=============================================================================
54
Capítulo 2. Simulación de sistemas no lineales continuos
almacenar las ecuaciones diferenciales y los parámetros requeridos, así como en la sección
2.4.1. Hay que tener en cuenta que las entradas vienen definidas en el vector de entradas
(u). Las ecuaciones diferenciales se escriben en términos de los estados (x).
%=============================================================================
% mdlDerivatives
% Return the derivatives for the continuous states.
%=============================================================================
%
function sys=mdlDerivatives(t,x,u)
% end mdlDerivatives
%
%=============================================================================
% mdlUpdate
% Handle discrete state updates, sample time hits, and major time step
% requirements.
%=============================================================================
%
%=============================================================================
% mdlUpdate
% Handle discrete state updates, sample time hits, and major time step
% requirements.
%=============================================================================
%
function sys=mdlUpdate(t,x,u)
% end mdlUpdate
%
%=============================================================================
% mdlOutputs
% Return the block outputs.
%=============================================================================
6. mdlOutputs: en esta función se especifican cuáles serán las salidas del sistema; pueden
ser todos los estados, un subset de los estados o variables determinadas a partir de los
estados.
55
Capítulo 2. Simulación de sistemas no lineales continuos
function sys=mdlOutputs(t,x,u)
% end mdlOutputs
%=============================================================================
%=============================================================================
% mdlGetTimeOfNextVarHit
% Return the time of the next hit for this block. Note that the result is
% absolute time. Note that this function is only used when you specify a
% variable discrete-time sample time [-2 0] in the sample time array in
% mdlInitializeSizes.
%=============================================================================
%
function sys=mdlGetTimeOfNextVarHit(t,x,u)
% end mdlGetTimeOfNextVarHit
%
%=============================================================================
% mdlTerminate
% Perform any end of simulation tasks.
%=============================================================================
%
function sys=mdlTerminate(t,x,u)
sys = [];
% end mdlTerminate
56
Capítulo 2. Simulación de sistemas no lineales continuos
Ejemplo 2.4.
Se desea evaluar el comportamiento de los estados del sistema reactivo Van de Vusse en un
CSTR cuando se aplican perturbaciones tanto en la tasa de dilución como en la concentración
de entrada de A. Los valores de parámetros, condiciones iniciales y variables de entrada son
tomados de la tabla 6.1. La simulación se realiza por 30 min usando una s-function, donde
en el tiempo=10 min se perturba la tasa de dilución con un cambio escalón de amplitud de 1/5
min−1 . En el tiempo=20 min, se perturba la concentración de A en el alimento con un cambio
escalón de amplitud -2 mol/l. Finalmente, se desea graficar tanto las señales de entrada como
de salida del modelo en Matlab de manera automática.
Solución
Para resolver este problema se utiliza el mismo programa de Simulink del ejemplo 2.3 sus-
tituyendo la Interpreted function con la s-function. En primera instancia, se muestra la
construcción del modelo matemático como una s-function para seguidamente realizar de
manera automática las gráficas usando otro programa de Matlab (i.e., usando un Callbacks).
La arquitectura de Simulink se muestra en la figura 2.13.
Figura 2.13. Implementación en Simulink del sistema Van de Vusse usando s-function
57
Capítulo 2. Simulación de sistemas no lineales continuos
%
% The general form of an MATLAB S-function syntax is:
% [SYS,X0,STR,TS,SIMSTATECOMPLIANCE] = SFUNC(T,X,U,FLAG,P1,...,Pn)
%
% What is returned by SFUNC at a given point in time, T, depends on the
% value of the FLAG, the current state vector, X, and the current
% input vector, U.
%
% FLAG RESULT DESCRIPTION
% ----- ------ --------------------------------------------
% 0 [SIZES,X0,STR,TS] Initialization, return system sizes in SYS,
% initial state in X0, state ordering strings
% in STR, and sample times in TS.
% 1 DX Return continuous state derivatives in SYS.
% 2 DS Update discrete states SYS = X(n+1)
% 3 Y Return outputs in SYS.
% 4 TNEXT Return next time hit for variable step sample
% time in SYS.
% 5 Reserved for future (root finding).
% 9 [] Termination, perform any cleanup SYS=[].
%
%
% The state vectors, X and X0 consists of continuous states followed
% by discrete states.
%
% Optional parameters, P1,...,Pn can be provided to the S-function and
% used during any FLAG operation.
%
% When SFUNC is called with FLAG = 0, the following information
% should be returned:
%
% SYS(1) =
Number of continuous states.
% SYS(2) =
Number of discrete states.
% SYS(3) =
Number of outputs.
% SYS(4) =
Number of inputs.
% Any of the first four elements in SYS can be specified
% as -1 indicating that they are dynamically sized. The
% actual length for all other flags will be equal to the
% length of the input, U.
% SYS(5) = Reserved for root finding. Must be zero.
% SYS(6) = Direct feedthrough flag (1=yes, 0=no). The s-function
% has direct feedthrough if U is used during the FLAG=3
% call. Setting this to 0 is akin to making a promise that
% U will not be used during FLAG=3. If you break the promise
% then unpredictable results will occur.
% SYS(7) = Number of sample times. This is the number of rows in TS.
%
%
% X0 = Initial state conditions or [] if no states.
%
% STR = State ordering strings which is generally specified as [].
%
% TS = An m-by-2 matrix containing the sample time
% (period, offset) information. Where m = number of sample
% times. The ordering of the sample times must be:
58
Capítulo 2. Simulación de sistemas no lineales continuos
%
% TS = [0 0, : Continuous sample time.
% 0 1, : Continuous, but fixed in minor step
% sample time.
% PERIOD OFFSET, : Discrete sample time where
% PERIOD > 0 & OFFSET < PERIOD.
% -2 0]; : Variable step discrete sample time
% where FLAG=4 is used to get time of
% next hit.
%
% There can be more than one sample time providing
% they are ordered such that they are monotonically
% increasing. Only the needed sample times should be
% specified in TS. When specifying more than one
% sample time, you must check for sample hits explicitly by
% seeing if
% abs(round((T-OFFSET)/PERIOD) - (T-OFFSET)/PERIOD)
% is within a specified tolerance, generally 1e-8. This
% tolerance is dependent upon your model's sampling times
% and simulation time.
%
% You can also specify that the sample time of the S-function
% is inherited from the driving block. For functions which
% change during minor steps, this is done by
% specifying SYS(7) = 1 and TS = [-1 0]. For functions which
% are held during minor steps, this is done by specifying
% SYS(7) = 1 and TS = [-1 1].
%
% SIMSTATECOMPLIANCE = Specifices how to handle this block when saving and
% restoring the complete simulation state of the
% model. The allowed values are: 'DefaultSimState',
% 'HasNoSimState' or 'DisallowSimState'. If this value
% is not speficified, then the block's compliance with
% simState feature is set to 'UknownSimState'.
%
% The following outlines the general structure of an S-function.
%
switch flag,
%%%%%%%%%%%%%%%%%%
% Initialization %
%%%%%%%%%%%%%%%%%%
case 0,
[sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes;
%%%%%%%%%%%%%%%
% Derivatives %
%%%%%%%%%%%%%%%
case 1,
sys=mdlDerivatives(t,x,u);
59
Capítulo 2. Simulación de sistemas no lineales continuos
%%%%%%%%%%
% Update %
%%%%%%%%%%
case 2,
sys=mdlUpdate(t,x,u);
%%%%%%%%%%%
% Outputs %
%%%%%%%%%%%
case 3,
sys=mdlOutputs(t,x,u);
%%%%%%%%%%%%%%%%%%%%%%%
% GetTimeOfNextVarHit %
%%%%%%%%%%%%%%%%%%%%%%%
case 4,
sys=mdlGetTimeOfNextVarHit(t,x,u);
%%%%%%%%%%%%%
% Terminate %
%%%%%%%%%%%%%
case 9,
sys=mdlTerminate(t,x,u);
%%%%%%%%%%%%%%%%%%%%
% Unexpected flags %
%%%%%%%%%%%%%%%%%%%%
otherwise
DAStudio.error('Simulink:blocks:unhandledFlag', num2str(flag));
end
% end sfuntmpl
%
%=============================================================================
% mdlInitializeSizes
% Return the sizes, initial conditions, and sample times for the S-function.
%=============================================================================
%
function [sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes
%
% call simsizes for a sizes structure, fill it in and convert it to a
% sizes array.
%
% Note that in this example, the values are hard coded. This is not a
% recommended practice as the characteristics of the block are typically
% defined by the S-function parameters.
%
sizes = simsizes;
60
Capítulo 2. Simulación de sistemas no lineales continuos
sys = simsizes(sizes);
%
% initialize the initial conditions
%
x0 = [2;1.117]; % Se ponen las condiciones iniciales [Ca Cb]
%
% str is always an empty matrix
%
str = [];
%
% initialize the array of sample times
%
ts = [0 0]; % Sistema continuo
% end mdlInitializeSizes
%
%=============================================================================
% mdlDerivatives
% Return the derivatives for the continuous states.
%=============================================================================
%
function sys=mdlDerivatives(t,x,u)
% Estados
Ca = x(1); % Estado 1
Cb = x(2); % Estado 2
% Entradas
61
Capítulo 2. Simulación de sistemas no lineales continuos
% Parámetros
% Ecuaciones diferenciales
% end mdlDerivatives
%
%=============================================================================
% mdlUpdate
% Handle discrete state updates, sample time hits, and major time step
% requirements.
%=============================================================================
%
function sys=mdlUpdate(t,x,u)
% end mdlUpdate
%
%=============================================================================
% mdlOutputs
% Return the block outputs.
%=============================================================================
%
function sys=mdlOutputs(t,x,u)
% end mdlOutputs
%
%=============================================================================
% mdlGetTimeOfNextVarHit
% Return the time of the next hit for this block. Note that the result is
% absolute time. Note that this function is only used when you specify a
% variable discrete-time sample time [-2 0] in the sample time array in
% mdlInitializeSizes.
%=============================================================================
%
function sys=mdlGetTimeOfNextVarHit(t,x,u)
62
Capítulo 2. Simulación de sistemas no lineales continuos
% end mdlGetTimeOfNextVarHit
%
%=============================================================================
% mdlTerminate
% Perform any end of simulation tasks.
%=============================================================================
%
function sys=mdlTerminate(t,x,u)
sys = [];
% end mdlTerminate
%=============================================================================
Figura 2.14. Resultados de la simulación del sistema Van de Vusse ante dos perturbaciones usando
una s-function
63
Capítulo 2. Simulación de sistemas no lineales continuos
En esta sección se mostró cómo implementar la s-function desde la plantilla que Matlab pro-
vee como ejemplo. Sin embargo, las versiones recientes de Matlab ofrecen una opción llamada
el s-function builder (figura 2.15). Esta aplicación se encuentra en la librería Simulin-
k/User defined function/S-Function builder. Allí, la ventana del aplicativo guía al usuario
para definir las características del modelo a resolver y generar la plantilla diseñada para las
necesidades del usuario. En principio, se debe llegar a M-file análogo al mostrado en el ejemplo.
64
Capítulo 2. Simulación de sistemas no lineales continuos
Figura 2.16. Configuración del método de solución de ecuaciones diferenciales ordinarias en Simulink
de ver aspectos numéricos de la solución; por ejemplo, para el caso del reactor Van de Vusse,
se puede analizar la manera en que evoluciona el tamaño de paso durante la evolución de la
solución como se muestra en la figura 2.17. El análisis de esta información es vital para una se-
lección sistemática de los parámetros del método numérico en casos de estudio más complejos.
65
Capítulo 2. Simulación de sistemas no lineales continuos
Figura 2.17. Estadísticas del desempeño del método numérico durante la solución del reactor Van de
Vusse en Simulink
66
Capítulo 3
Resumen
Por lo general, en la literatura de control de procesos clásica se inicia con manipulaciones
matemáticas para representar ecuaciones diferenciales ordinarias lineales en diferentes
dominios para posteriormente analizarlas. En este capítulo se muestran herramientas
útiles para realizar la transformación de ecuaciones al dominio de Laplace y determinar
la función de transferencia y representación cero-polo, obtener espacio de estados y poder
migrar de una representación a otra. Adicionalmente, se explora el motor simbólico de
Matlab para realizar operaciones entre el dominio del tiempo y de Laplace. Finalmente,
se muestra una herramienta de Matlab que permite obtener modelos lineales a partir de
datos de simulación o experimentales.
3.1. Introducción
La estructura de este capítulo se basa en la estrategia para la obtención de modelos lineales y su
posterior uso, ya sea para representar el sistema o para encontrar la solución de las ecuaciones.
Teniendo en cuenta que el punto de partida son ecuaciones diferenciales ordinarias, se usa el
siguiente proceso (Seborg et al., 2011):
67
Capítulo 3. Sistemas lineales continuos
Parte de la primera etapa fue realizada en el capítulo anterior. Para cada una de estas etapas
restantes, las herramientas disponibles para llevarlas a cabo serán revisadas en las siguientes
secciones.
F (x) = 0 (3.1)
68
Capítulo 3. Sistemas lineales continuos
Donde:
%---------------------------------------------------------------------------------
Fobjcal = @(X)Ecuaciones(X,parámetros); % Función anónima
options = optimoptions('fsolve','Display','iter'); % Opciones
Co = [C0]; % Condiciones iniciales
[x,fval,exitflag,output]= fsolve(Fobjcal,Co,options); % Solucionador
%---------------------------------------------------------------------------------
En esta sintaxis se utiliza una función anónima llamada “Ecuaciones”, que puede estar ubi-
cada ya sea en otro programa (M-file) o en el mismo. Esta función retorna en un vector el
valor de las ecuaciones diferenciales que deben ser igualadas a cero por el optimizador (i.e.,
F(X) = 0) y la variable independiente son los estados (X) iterados por el solver para encon-
trar el estado estable. Puesto que las ecuaciones diferenciales son función de las entradas al
sistema y los parámetros, se propone usar estos últimos como entradas a la función anónima
utilizando el vector “parámetros”; entonces pueden ser definidos por fuera de las ecuaciones
dándole a la función más versatilidad.
En la variable options se guardan las modificaciones al método numérico que deseen rea-
lizarse. Para esto se utiliza el comando optimoptions. Teniendo en cuenta que la función
optimoptions aplica para el toolbox de optimización de Matlab, en primera instancia se tiene
que especificar cuál función se desea modificar; en este caso aplica para función fsolve. Se-
guidamente, se realizan las modificaciones al método numérico según la conveniencia del caso.
En la sintaxis mostrada, se especifica que se desean ver las iteraciones de la optimización;
esta configuración es muy útil para valorar en una etapa temprana si hay tendencia a que el
problema diverja o si hay demoras en la simulación.
Otras opciones para manipular son: tolerancia de terminación de la función (TolFun), tole-
rancia de terminación en X (TolX), máximo número de iteraciones permitido (MaxIter) y
máximo número de evaluación de funciones permitido (MaxFunEvals). En total, se tienen 50
parámetros para modificar.
69
Capítulo 3. Sistemas lineales continuos
En la última línea del código, se muestra cómo se ejecuta el solver para encontrar el valor de
los estados en el estado estable. Como variables de salida se tienen:
La variable de salida exitflag es muy importante, pues informa qué criterio de convergencia
ha sido alcanzado; es relevante aclarar que no necesariamente es una solución viable. Esto es
útil para problemas de optimización no convexa, donde se pueden tener varias soluciones o
ninguna. Los valores que puede tomar exitflag son:
Ejemplo 3.1.
Se desea determinar el valor de las concentraciones en estado estable para el reactor Van de
Vusse (VDV) cuando se tiene unos valores de entrada de D̄ = 4/7 min−1 y C̄Af = 10 mol l−1
(tabla 6.1).
Solución
Para resolver este problema, se va a emplear como base la rutina implementada en la sección
2.3.3, donde se utilizó un solo programa con múltiples funciones, ya que se recicla gran parte
del código. Entonces, se tiene un programa principal que contiene los parámetros/variables
requeridos, la ejecución de la rutina de solución fsolve y la presentación de resultados.
La rutina de solución ejecuta a la función anónima que contiene las ecuaciones diferenciales
llamada ODE_equations. Teniendo en cuenta que la optimización requiere unos estimados
iniciales para la solución, se toman las condiciones iniciales del sistema como estimados. La
implementación es:
70
Capítulo 3. Sistemas lineales continuos
%--------------------------------------------------------------------------
% Solución del sistema reactor VDV en estado estable
% Sistema de ecuaciones en el mismo script-función
% Método de solución : fsolve
function SS_VDVReactor
clc
clear
close all
%% Solución
%% Resultados
disp(' ')
disp('El estado estable del sistema es:')
disp(' ')
disp(['C_ss = [' num2str(C_ss(1)) ' ' num2str(C_ss(2)) ']'])
disp(' ')
disp(['Criterio de salida (Flag): ' num2str(ExitFlag)])
disp(' ')
disp('Estadísticas de la solución: ')
disp(' ')
disp(output)
end
dC_dt = zeros(length(C),1);
dC_dt(1) = D*(Caf-Ca) - k1*Ca -k3*Ca^2; % ODE Ca
dC_dt(2) = -D*Cb + k1*Ca - k2*Cb; % ODE Cb
71
Capítulo 3. Sistemas lineales continuos
end
%--------------------------------------------------------------------------
Cabe destacar que la función ODE_equations es la misma empleada en la sección 2.3.3. Sin
embargo, hay que notar que para la función anónima Fobjcal se tiene solo el vector de con-
centraciones como variable, pues el tiempo no se requiere. El resto de entradas que tiene la
función ODE_equations son interpretadas como parámetros y se debieron definir con ante-
rioridad.
Nótese que la variable de tiempo, llamada tss, no tiene uso aparente; se define para hacer del
tiempo un parámetro. Esto no afecta la solución, ya que las ecuaciones diferenciales no tienen
el tiempo de manera explícita. Adicionalmente, se modificaron más opciones en la configura-
ción del algoritmo incluyendo: la gráfica de la condición first order optimality (convergencia),
la evolución de las iteraciones, la tolerancia para la convergencia de la función, la tolerancia
para la convergencia de C y el número máximo de iteraciones. Por último, se muestran los
resultados impresos en la pantalla (solución, criterio de convergencia y datos de la optimiza-
ción). Los resultados arrojados por el programa son:
Equation solved.
C_ss = [3 1.117]
Estadísticas de la solución:
iterations: 3
funcCount: 12
algorithm: ’trust-region-dogleg’
firstorderopt: 1.7375e-12
message: ’Equation solved. fsolve completed because ...
72
Capítulo 3. Sistemas lineales continuos
Para este caso, la solución se logró solo en tres iteraciones dando como resultado C̄A = 3 mol/l
y C̄B = 1,117 mol/l. El criterio de solución fue Flag = 1, lo que indica que la función fsolve
convergió a una raíz viable. Como bonus, se pidió a la función retornar el valor del first order
optimality o condición Karush-Kuhn-Tucker (KKT) (véase firstorderopt). Esta condición
es una medida que indica qué tan cerca está la solución del óptimo basada en una función
auxiliar en términos del Lagrangian (Nocedal y Wright, 1999). Esta se conoce como una
condición necesaria, pero no suficiente:
3.2.2.1. Comentarios
73
Capítulo 3. Sistemas lineales continuos
muy demandantes en sistemas grandes, debido al algoritmo usado para el cómputo del pre-
condicionador (Coleman y Li, 1996).
• No convergencia del optimizador: esto puede ocurrir incluso cuando las condiciones
iniciales para el método sean factibles. En este caso, se recomienda cambiar el set de
condiciones iniciales e intentar de nuevo.
• Durante la solución, el optimizador itera valores de los estados inviables, por ejemplo:
concentraciones negativas. En este caso, se recomienda migrar del fsolve a un método
de optimización que permita el uso de restricciones. Como característica, se pueden
imponer límites al espacio de búsqueda y así evitar buscar los estados en zonas no
factibles. Esto se mostrará en la siguiente sección.
Donde:
74
Capítulo 3. Sistemas lineales continuos
Para este problema en particular, solo es de interés restringir la optimización con límites
inferior y superior para los estados (i.e., lb y ub); por tanto, las demás variables no se
modifican. Ahora, en el caso anterior se escribió una función en Matlab dedicada a retornar
el valor de las derivadas en un punto determinado de los estados y las entradas. El objetivo
del algoritmo fsolve era encontrar los estados que llevaran a cada una de las ecuaciones
diferenciales lo más cerca a cero. En este caso, se puede ver en el sistema de ecuaciones (3.2)
que la función objetivo es una sola, por lo que se propone usar el concepto de la norma-2
(magnitud del vector o longitud euclidiana) para agrupar todas las funciones en una sola (se
pueden plantear otras funciones objetivo):
"k=1 #1/p
p
X
kvkp = |vk | , p=2 (3.3)
N
%--------------------------------------------------------------------------------
Fobjcal = @(X)Ecuaciones(X,parámetros); % Función anónima
options = optimoptions('fmincon','Display','iter'); % Opciones del solucionador
A = []; % Matriz de restricciones lineales en desigualdad
B = []; % Matriz de restricciones lineales en desigualdad
Aeq = []; % Matriz de restricciones lineales en igualdad
Beq = []; % Matriz de restricciones lineales en igualdad
LB = [0 0]; % Límite inferior para los estados estables
UB = [inf inf]; % Límite superior para los estados estables
Nonlcon = []; % Restricciones no lineales
xo = [x0]; % Vector con las condiciones iniciales para la optimización
[x,fval,exitflag,output]= fmincon(Fobjcal,xo,A,B,Aeq,Beq,LB,UB,Nonlcon,options);
%-------------------------------------------------------------------------------
Nótese que la estructura de la función anónima, las opciones y variables de salida son análogas
a las mostradas en la sección 3.2.2. Adicionalmente, se tienen las matrices y vectores de las
restricciones lineales y no lineales en igualdad y desigualdad; estas restricciones no son usadas
y, por ende, se dejan como matrices vacías. Las únicas variables modificadas son LB y UF,
donde se restringen los estados a tener valores positivos. Finalmente, la función fmincon se
ejecuta indicando las restricciones y opciones. Los valores que puede tomar exitflag son
(MatWorks, 2022e):
75
Capítulo 3. Sistemas lineales continuos
Ejemplo 3.2.
Se desea determinar el valor de las concentraciones en estado estable para el reactor Van de
Vusse, cuando se tiene unos valores de entrada de D = 4/7 min−1 y CAf = 10 mol l−1 (tabla
6.1). Como precaución, se deben incluir restricciones para garantizar que las concentraciones
solo tomen valores positivos durante la solución.
Solución
Para resolver este problema, se va a emplear como base la rutina implementada en la sec-
ción 3.2.2. Se tiene entonces casi el mismo programa principal que contiene los parámetros
y variables requeridos. Enseguida, se muestra la ejecución de la rutina de solución fmincon
y posteriormente la presentación de resultados. La rutina de solución ejecuta la función anó-
nima que contiene las ecuaciones diferenciales y la definición de la función objetivo llamada
igualmente ODE_equations. Teniendo en cuenta que la optimización requiere unos estimados
iniciales para la solución, se toman las condiciones iniciales del sistema como estimados. La
implementación es:
76
Capítulo 3. Sistemas lineales continuos
%--------------------------------------------------------------------------
% Solución del sistema reactor VDV en estado estable
% Sistema de ecuaciones en el mismo script-función
% Método de solución: fmincon
function SS_VDVReactor_Const
clc
clear
close all
%% Solución
[C_ss,Fval,ExitFlag,output] = fmincon(Fobjcal,Co,A,B,Aeq,Beq,...
LB,UB,Nonlcon,options);
%% Resultados
disp(' ')
disp('El estado estable del sistema es:')
disp(' ')
disp(['C_ss = [' num2str(C_ss(1)) ' ' num2str(C_ss(2)) ']'])
disp(' ')
disp(['Valor de la función objetivo: ' num2str(Fval)])
disp(' ')
disp(['Criterio de salida (Flag): ' num2str(ExitFlag)])
disp(' ')
disp('Estadísticas de la solución: ')
disp(' ')
disp(output)
77
Capítulo 3. Sistemas lineales continuos
end
dC_dt = zeros(length(C),1);
dC_dt(1) = D*(Caf-Ca) - k1*Ca -k3*Ca^2; % ODE Ca
dC_dt(2) = -D*Cb + k1*Ca - k2*Cb; % ODE Cb
end
%--------------------------------------------------------------------------
Nótese que la función ODE_equations es casi la misma empleada en el ejemplo anterior, con
la diferencia de que la salida no es un vector con las derivadas, sino la función objetivo defi-
nida como la magnitud vectorial de estas (i.e., norma-2). Además, la misma modificación de
la función anónima es usada, donde el tiempo se remueve como variable independiente y se
usa un valor arbitrario de tss = 100.
En este caso, se modifican opciones del algoritmo que incluyen: la gráfica del valor de la fun-
ción objetivo (convergencia), la evolución de las iteraciones, tolerancia para la convergencia
de la función, la tolerancia para la convergencia de C y el número máximo de iteraciones.
Los resultados se imprimen en la pantalla como solución, criterio de convergencia, valor de la
función objetivo y datos de la optimización.
First-order Norm of
Iter F-count f(x) Feasibility optimality step
0 3 2.388187e+00 0.000e+00 2.204e+00
1 7 6.519876e-01 0.000e+00 2.456e+00 1.098e+00
2 11 2.687852e-01 0.000e+00 1.414e+00 3.638e-01
3 17 2.327100e-01 0.000e+00 2.010e+00 1.174e-01
4 21 7.300987e-02 0.000e+00 2.394e+00 9.941e-02
5 26 4.859785e-02 0.000e+00 2.338e+00 3.115e-02
6 30 2.971362e-02 0.000e+00 2.275e+00 3.152e-02
7 33 2.692830e-02 0.000e+00 2.357e+00 2.062e-02
8 36 8.679474e-03 0.000e+00 2.287e+00 1.186e-02
9 40 6.418413e-03 0.000e+00 2.167e+00 6.450e-03
10 44 5.922775e-03 0.000e+00 2.015e+00 3.533e-03
11 47 2.152741e-03 0.000e+00 2.070e+00 3.194e-03
12 51 5.414107e-04 0.000e+00 2.496e+00 1.143e-03
13 61 3.684503e-04 0.000e+00 2.388e+00 3.349e-04
14 64 1.899786e-04 0.000e+00 2.306e+00 1.809e-04
15 72 1.637121e-04 0.000e+00 1.008e+00 1.685e-04
78
Capítulo 3. Sistemas lineales continuos
fmincon stopped because the size of the current step is less than
the selected value of the step size tolerance and constraints are
satisfied to within the default value of the constraint tolerance.
C_ss = [3 1.117]
Estadísticas de la solución:
iterations: 24
funcCount: 111
constrviolation: 0
stepsize: 9.1308e-07
algorithm: ’interior-point’
firstorderopt: 2.2078
cgiterations: 10
message: ’Local minimum possible. Constraints satisfied
Para este caso, la solución se alcanzó después de 24 iteraciones dando como resultado C̄A = 3
mol/l y C̄B = 1,117 mol/l (se trata de la misma solución anterior). El criterio de solución
fue Flag = 2, lo que indica que el paso en el vector C es menor a la tolerancia definida en
las opciones. Como bonus se pidió a la función retornar la evolución del valor de la función
objetivo; esto es relevante para validar si el mínimo es una raíz. El resultado para la evolución
de la solución se muestra en la figura 3.2.
79
Capítulo 3. Sistemas lineales continuos
Figura 3.2. Evolución del valor de la función objetivo durante la solución del problema
Reto: determinar la influencia del algoritmo y otro set de condiciones iniciales al mo-
mento de realizar la optimización.
80
Capítulo 3. Sistemas lineales continuos
%--------------------------------------------------------------------------
[xeq] = trim('sys',x0)
%--------------------------------------------------------------------------
Donde:
Ejemplo 3.3.
Se desea determinar el valor de las concentraciones en estado estable para el reactor Van de
Vusse, cuando se tiene unos valores de entrada de D = 4/7 min−1 y CAf = 10 mol l−1 (tabla
6.1).
Solución
Para resolver el problema, se recicla la implementación del modelo no lineal en Simulink usado
en el ejemplo 2.4 de la sección 2.4.2, removiendo los bloques de conectividad con Matlab (i.e.,
To workspace) y el Callback que permitía ejecutar el programa que realiza las figuras en
Matlab. Este programa es guardado con el nombre de VDV_NLsimulation. La idea es utilizar
este modelo y ejecutar la función trim desde un programa maestro escrito como un M-file.
Ahora, lo primero que se requiere es reconfigurar las entradas en el programa de Simulink, por
lo que se definen entonces como el valor dado en el ejemplo antes de realizar perturbaciones.
Es decir, las entradas paso quedan con un valor constante en todo el dominio del tiempo y
este corresponde a los valores alrededor de los cuales se desea conocer el punto de equilibrio.
La implementación se muestra en la figura 3.3.
81
Capítulo 3. Sistemas lineales continuos
Figura 3.3. Implementación modificada del reactor Van de Vusse en Simulink usando una S-fuction
Nótese de la implementación que los valores de las entradas se dejaron especificados como
variables (i.e., D_eq y Caf_eq); de esta manera, estas pueden ser actualizadas desde el pro-
grama maestro. Para determinar el punto de equilibrio, se implementa un código en Matlab
que cargue los valores de las entradas y las condiciones iniciales para la iteración realizada
por el algoritmo de optimización. La implementación del M-file es:
%--------------------------------------------------------------------------
% Determinación de puntos de equilibrio
% Modelo implementado como una S-function
% Método de solución: Sequential Quadratic Programming (SQP)
clc
clear
close all
82
Capítulo 3. Sistemas lineales continuos
%--------------------------------------------------------------------------
Ceq =
3.0000
1.1170
Ueq =
0.5714
10.0000
83
Capítulo 3. Sistemas lineales continuos
a cada estado y cada entrada alrededor del estado estable, dando como resultado la ecuación
(3.5) en variables de desviación.
dx1
= f1 (x1 , x2 , ..., xn , u1 , u2 , ..., ur ) (3.4)
dt
dx2
= f2 (x1 , x2 , ..., xn , u1 , u2 , ..., ur )
dt
..
.
dxn
= fn (x1 , x2 , ..., xn , u1 , u2 , ..., ur )
dt
n r
dx0i ∂fi ∂fi
x0j + u0k
X X
= (3.5)
dt j=1
∂xj x̄,ū
k=1
∂uk x̄,ū
Donde:
x : vector de estados
x’ : vector de estados en variables de desviación
u : vector de entradas
u’ : vector de entradas en variables de desviación
f : función no lineal
x̄ : estado estable para los estados
ū : estado estable para las entradas
n : número de estados
r : número de entradas
i : índice para la ecuación diferencial
j : índice para cada estado
k : índice para cada entrada
La ecuación derivada para la linealización de un modelo de ecuaciones diferenciales ordinarias
se puede volver más compacta si se usa el concepto de la matriz jacobiana para las derivadas
con respecto tanto a los estados como a las entradas. Entonces, el sistema de ecuaciones puede
representarse por:
dx0
= Ax0 + Bu0 (3.6)
dt
Donde x0 y u0 son los vectores de estados y entradas en variables de desviación, y las matrices
A y B vienen dadas por la representación matricial, así:
84
Capítulo 3. Sistemas lineales continuos
Ejemplo 3.4.
Se desea linealizar el sistema de ecuaciones para el reactor Van de Vusse alrededor del estado
estable C̄A y C̄B de manera analítica.
Solución
dCA 2
=D (CAf − CA ) − k1 CA − k3 CA
dt
dCB
=D (CBf − CA ) + k1 CA − k2 CB
dt
2 , (b)
Se puede notar que en la primera ecuación hay tres parámetros de no linealidad: (a) el CA
la multiplicación de una entrada D con una salida CA y (c) la multiplicación de dos entradas
CAf y D. En la segunda ecuación, hay una sola fuente de no linealidad en la multiplicación
de la entrada D con una salida CB (nótese que CBF = 0). Ya que las dos ecuaciones son no
lineales, estas deben ser linealizadas utilizando el método de Taylor, lo que permite obtener
las derivadas parciales de las ecuaciones con respecto a las entradas y los estados.
85
Capítulo 3. Sistemas lineales continuos
dC 0 A
= −D̄ − k1 − 2k3 C̄A C 0 A + C̄Af − C̄A D0 + D̄C 0 Af
(3.9)
dt
dC 0 B
= k1 C 0 A + −k2 − D̄ C 0 B − C̄B D0
dt
Es muy importante el orden al realizar la linealización, por lo que se recomienda derivar
primero con respecto a estados y después entradas, esto con el fin de facilitar la representación
matricial del sistema conocida como espacio de estados. Utilizando la forma genérica del
espacio de estados se tiene:
Ċ 0 A −D̄ − k1 − 2k3 C̄A 0 C0A C̄Af − C̄A D̄ D0
= + (3.10)
Ċ 0 B C0B C 0 Af
k1 −k2 − D̄ −C̄B 0
C0A 1 0 C0A 0 0 D0
= +
C0B 0 1 C0B 0 0 C 0 Af
En este caso, se supone que todos los estados son medidos y, por tanto, la matriz C es la
identidad.
%--------------------------------------------------------------------------
syms 'Lista de variables independientes (estados, entradas y parámetros)'
Eq = ['Ecuaciones simbólicas'] % Ecuaciones separadas por comas
A = jacobian(Eq,['Estados']);
B = jacobian(Eq,['Entradas']]);
%--------------------------------------------------------------------------
Donde:
86
Capítulo 3. Sistemas lineales continuos
y entradas en el estado estable previamente estimado y usar el comando eval. Para mostrar
cómo se ejecutan de manera práctica la obtención y la edición del formato de las derivadas
simbólicas, se tiene el ejemplo 3.5.
Ejemplo 3.5.
Se desea linealizar el sistema de ecuaciones para el reactor Van de Vusse alrededor del estado
estable determinado en la sección 3.2 y conocer su valor numérico.
Solución
Se recomienda, para este caso, realizar de manera manual las derivadas para luego validarlas
con los resultados del programa. Se parte de los parámetros mostrados en la tabla 6.1, los
valores de entrada de D̄ = 4/7 min−1 y C̄Af = 10 mol l−1 , y los estimados de los estados
en el estado estable C̄A = 3 mol l−1 y C̄B = 1,117 mol l−1 . Con las indicaciones mostradas,
se definen las variables en principio como simbólicas: entradas, estados y parámetros del
modelo; posteriormente, se definen las ecuaciones diferenciales a linealizar en forma vectorial
y se ejecuta el comando jacobian. Para determinar el valor numérico de las derivadas, se
definen los valores numéricos de los parámetros, las entradas y los estados en el estado estable.
Finalmente, se usa el comando eval para conocer las matrices. El código de Matlab es:
%--------------------------------------------------------------------------
% Determinación del jacobiano simbólico
% Sistema de ecuaciones en el mismo script-función
% Método de solución: jacobian
clc
clear
%% Derivación simbólica
syms D Caf Ca Cb k1 k2 k3
87
Capítulo 3. Sistemas lineales continuos
%% Resultados
disp('Anum = ')
disp(' ')
disp(Anum)
disp('Bnum = ')
disp(' ')
disp(Bnum)
%--------------------------------------------------------------------------
Nótese que las ecuaciones diferenciales fueron introducidas de la misma manera que se usó
anteriormente para su solución numérica, lo cual permite reciclar el código realizado. Por otro
lado, hay que resaltar que Matlab fue concebido para efectuar operaciones numéricas y su
fuerte no es mostrar las ecuaciones simbólicas en un formato reconocible. Cuando las ecua-
ciones son muy complejas, es difícil reconocer cuál es realmente la forma de las ecuaciones
entre tantos paréntesis y operaciones aritméticas. En la implementación realizada, se utilizó
el comando pretty, que tiene como propósito mostrar las ecuaciones un poco más cercanas
al tipo matemático. Para la implementación realizada es muy importante definir todas las
variables simbólicas con valores numéricos; de otro modo, las matrices permanecerán como
simbólicas.
Asym =
/ - D - k1 - 2 Ca k3, 0 \
| |
\ k1, - D - k2 /
Bsym =
/ Caf - Ca, D \
| |
\ -Cb, 0 /
Anum =
88
Capítulo 3. Sistemas lineales continuos
-2.4048 0
0.8333 -2.2381
Bnum =
7.0000 0.5714
-1.1170 0
Teniendo en mente que se desea linealizar un sistema de ecuaciones diferenciales, en este caso
f (x) corresponde a cada ecuación diferencial a linealizar, x a cada estado/entrada y ∆x es
una perturbación muy pequeña, que define la exactitud de la aproximación. Para determinar
entonces todas las derivadas, se aprovecha la representación vectorial que se tiene del sistema
de ecuaciones para así ahorrar esfuerzo computacional. El uso de la aproximación se muestra
en el ejemplo 3.6.
Ejemplo 3.6.
Se desea linealizar numéricamente el sistema de ecuaciones para el reactor Van de Vusse al-
rededor del estado estable determinado en la sección 3.2.
89
Capítulo 3. Sistemas lineales continuos
Solución
De manera análoga al caso anterior, se parte de los parámetros mostrados en la tabla 6.1, los
valores de entrada de D̄ = 4/7 min−1 y C̄Af = 10 mol l−1 , y los estimados de los estados en
el estado estable C̄A = 3 mol l−1 y C̄B = 1,117 mol l−1 . Para este caso, se desea aprovechar
la implementación de las ecuaciones diferenciales realizada en la sección 2.3.3; por tanto, el
resto de código se adapta a la función previamente elaborada. Adicionalmente, se aprovecha la
derivada simbólica determinada en el ejemplo 3.5 para validar los resultados de la estimación
numérica.
La implementación es:
%--------------------------------------------------------------------------
% Linealización del reactor Van de Vusse
% Sistema de ecuaciones en el mismo script-funciones separadas
% Método de solución: Analítico versus numérico
function Jacobian_VDVReactor
clc
clear
close all
%% Solución
90
Capítulo 3. Sistemas lineales continuos
%% Resultados
format shortG
ErrA = (Aana - Anum) % Error absoluto de cada término en A
ErrB = (Bana - Bnum) % Error absoluto de cada término en B
end
function [A,B]=Janalitic(t,C,k1,k2,k3,Caf,D)
a11=-D-k1-2*k3*Ca;
a12=0;
a21=k1;
a22=-D-k2;
A = [a11 a12; a21 a22]; % Jacobiano analítico en A
b11=Caf-Ca;
b12=D;
b21=-Cb;
b22=0;
B = [b11 b12; b21 b22]; % Jacobiano analítico en B
end
function [A,B]=Jbackwards(t,C,k1,k2,k3,Caf,D)
U = [D Caf];
dC= 1e-3; % Perturbación para el método de diferenciación
F =ODE_equations(t,C,k1,k2,k3,U);
for i=1:length(C)
C(i) = C(i)+dC;
A(:,i)=(ODE_equations(t,C,k1,k2,k3,U)-F)./dC;
C(i) = C(i)-dC;
end
for i=1:length(U)
U(i) = U(i)+dC;
B(:,i)=(ODE_equations(t,C,k1,k2,k3,U)-F)./dC;
U(i) = U(i)-dC;
end
end
91
Capítulo 3. Sistemas lineales continuos
dC_dt = zeros(length(C),1);
dC_dt(1) = D*(Caf-Ca) - k1*Ca -k3*Ca^2; % ODE Ca
dC_dt(2) = -D*Cb + k1*Ca - k2*Cb; % ODE Cb
end
%--------------------------------------------------------------------------
Para realizar la diferenciación numérica, se explora la función que contiene las ecuaciones
diferenciales de forma vectorial. De esta manera no se realiza aproximación de la derivada
ecuación por ecuación y estado por estado, sino estado por estado para todas las ecuaciones
simultáneamente; así las matrices A y B se llenaran una columna a la vez y no elemento por
elemento.
Para materializar la estrategia, primero se calcula el valor de las derivadas sin perturbación
(i.e., f (x)) ejecutando la función ODE_equations alimentando el valor de entradas y estados
en el estado estable. Ahora, para la determinación del valor de las ecuaciones diferenciales
perturbadas (i.e., f (x + ∆x)), se modifica el vector de estados/entradas (uno a la vez) con un
valor muy pequeño (e.g., 1×10−3 ); entre más pequeño sea este número, se mejora la exactitud
de la aproximación numérica. Nótese que la perturbación de cada estado y entrada se hace en
el ciclo; después de aproximar la derivada, la perturbación debe removerse del vector original.
Anum =
-2.4049 0
0.83333 -2.2381
Bnum =
7 0.57143
-1.117 0
ErrA =
0.00016667 0
-1.3023e-13 -3.0997e-13
ErrB =
3.2685e-13 5.7043e-13
-1.8741e-13 0
92
Capítulo 3. Sistemas lineales continuos
Como se observa en los resultados, es importante resaltar que la aproximación numérica solo
tiene diferencias apreciables en el término a11.
Reto 1: teniendo en cuenta que no siempre se tiene la solución analítica para comparar
resultados, es importante familiarizarse con el uso del paso y los tipos de aproximaciones
de la derivada. El reto consiste en implementar la derivada centrada y evidenciar la
exactitud de los métodos en función del tamaño de la perturbación realizada.
Reto 2: cuando el sistema es muy pequeño, vale la pena implementar las ecuaciones
simbólicas directamente para calcular los valores de las matrices como se mostró en
el ejemplo 3.5. Sin embargo, para modelos grandes no se ve bien. Hay un comando de
Matlab que permite convertir el modelo simbólico directamente en una función anónima
llamado matlabFunction. Explorar esta opción y validar los resultados.
93
Capítulo 3. Sistemas lineales continuos
%--------------------------------------------------------------------------
[A,B,C,D]=linmod('Modelo',x,u);
[A,B,C,D]=linmod2('Modelo',x,u);
%--------------------------------------------------------------------------
Donde:
Algo importante es que el orden de estados y entradas debe ser el mismo usado en la im-
plementación del modelo en Simulink y este se mantiene en el resultado (i.e., modelo de
espacio de estados obtenido). Finalmente, cuando no se alimentan estados ni entradas, por
defecto los algoritmos utilizan cero, lo que significa que el modelo se linealiza alrededor de cero.
Ejemplo 3.7.
Se desea linealizar el sistema de ecuaciones para el reactor Van de Vusse alrededor del estado
estable determinado en la sección 3.2, a partir de un modelo implementado en Simulink.
Solución
Siguiendo la misma estructura, se parte de los parámetros mostrados en la tabla 6.1, los valores
de entrada de D̄ = 4/7 min−1 y C̄Af = 10 mol l−1 , y los estimados de los estados en el estado
estable C̄A = 3 mol l−1 y C̄B = 1,117 mol l−1 . En este caso se tienen dos programas: la imple-
mentación del modelo en Simulink y un código en Matlab que llama al modelo para realizar
la linealización. Aquí se desea aprovechar la implementación de las ecuaciones diferenciales
realizada en la sección 3.2.4. En esta implementación se sustituyen las entradas del modelo
por puertos llamados Inport y las salidas por puertos llamados Outport; esto es vital para
los algoritmos linmod. La implementación en Simulink se llama VDV_NLsimulation_ports.
Seguidamente se construye un programa en Matlab donde se definen las variables de entrada
y los estados en el estado estable (debe recordarse que los parámetros del modelo en este
caso están dentro de la S-function). Entonces, se aplican las funciones para determinar las
matrices del espacio de estado (i.e., A, B, C y D) usando ambas funciones planteadas. En la
figura 3.4 se muestra la nueva implementación de Simulink.
94
Capítulo 3. Sistemas lineales continuos
Figura 3.4. Implementación modificada del reactor Van de Vusse en Simulink reemplazando entradas
y salidas por puertos Inport y Outport
%--------------------------------------------------------------------------
% Linealización del reactor Van de Vusse
% Sistema en Simulink
% Método de solución: numérico-linmod
clc
clear
close all
%% Solución
[A1,B1,C1,D1]=linmod('VDV_NLsimulation_ports',Css,Uss);
[A2,B2,C2,D2]=linmod2('VDV_NLsimulation_ports',Css,Uss);
%% Resultados
format shortG
disp('A= ')
disp(' ')
disp(A1)
disp('B = ')
disp(' ')
disp(B1)
ErrA = (A1 - A2) % Diferencias absolutas de cada término en A
95
Capítulo 3. Sistemas lineales continuos
%--------------------------------------------------------------------------
A=
-2.4048 0
0.83333 -2.2381
B =
7 0.57143
-1.117 0
ErrA =
3.1084e-08 0
-5.6697e-12 1.6411e-11
ErrB =
0 0
0 0
96
Capítulo 3. Sistemas lineales continuos
%--------------------------------------------------------------------------
syms t
Ft = 'Función';
Fs = laplace(Ft);
pretty(Fs)
%--------------------------------------------------------------------------
Donde:
La función laplace puede ser aplicada a vectores de funciones y también modificar la varia-
ble de Laplace “s” que se tiene por defecto. La aplicación del procedimiento se muestra en el
ejemplo 3.8.
Ejemplo 3.8.
Solución
Se plantea un programa para obtener todas las transformadas con una sola ejecución. Inicial-
97
Capítulo 3. Sistemas lineales continuos
mente, se definen todas las variables simbólicas por utilizar, incluyendo los parámetros de las
funciones y el tiempo. A continuación, se construye un vector con todas las funciones a las que
se desea aplicar la transformada. Es importante tener en cuenta que las funciones individuales
son simbólicas, por lo que el vector también lo será. Después se ejecuta la función laplace
para obtener todas las transformadas en un vector y, finalmente, se muestran los resultados.
%--------------------------------------------------------------------------
% Obtención de transformadas de Laplace
% Método de solución: Matlab simbólico
clc
clear
syms t a w
%--------------------------------------------------------------------------
Entre las funciones usadas, la última es especial: se trata de una función paso o escalón uni-
tario retrasado dos unidades de tiempo. Es decir, es una función paso convencional a la que
se le aplicó el teorema de traslación de dos unidades de tiempo. Esto se conoce también como
una función a la que se le aplicó un transport delay. Como tal, esta función es discontinua y
debe plantearse de manera especial en Matlab. Para esto se usa la función heaviside, que
no es propiamente la función paso sino una aproximación de la discontinuidad (figura 3.5).
Esta función puede ser usada como un interruptor que prende y apaga otras funciones.
n 1 t≥0
u(t) = 0,5 t = 0 (3.13)
0 t<0
98
Capítulo 3. Sistemas lineales continuos
Figura 3.5. Aproximación de la función paso utilizando la función heaviside(t), evaluada para
un vector de tiempo entre −1 y 1
/ 1 a s a w 720 exp(-2 s) \
| -----, -------, -------, - --------, --------- |
| a + s 2 2 2 2 7 s |
\ s + w s + w (a - s) /
Se puede validar en la literatura la veracidad de las transformadas. Cabe resaltar que el resul-
tado de la transformada de la última función corresponde a la aplicación de la transformada
a la función paso y al uso del teorema de traslación de Laplace.
n1 t ≥ a
u(t) = (3.14)
0 t<a
99
Capítulo 3. Sistemas lineales continuos
Para realizar el procedimiento, se utiliza la función residue de Matlab. Esta función en-
cuentra numeradores (Ri ), polos (Pi ) y término residual (K) del cociente de dos polinomios
B(s)/A(s). Para el caso en el que no hay raíces repetidas se tiene el siguiente modelo de
expansión en fracciones parciales:
%--------------------------------------------------------------------------
B = ['Coeficientes'];
A = ['Coeficientes'];
[r,p,k] = residue(B,A);
%--------------------------------------------------------------------------
Donde:
100
Capítulo 3. Sistemas lineales continuos
Ejemplo 3.9.
Solución
De manera análoga al ejercicio 3.8, se plantea un solo programa para realizar todas las trans-
formaciones. Inicialmente, se definen todos los numeradores y denominadores en forma de
vectores; luego, se aplica la función residue y, por último, se imprimen los resultados.
%--------------------------------------------------------------------------
% Expansión en fracciones parciales
% Los tres tipos de raíces son evaluados
clc
101
Capítulo 3. Sistemas lineales continuos
clear
%% Planteamiento de ecuaciones
num1 = 1; % Numerador 1
den1 = [1 6 11 6 0]; % Denominador 1
num2 = [1 1]; % Numerador 2
den2 = [1 4 4 0]; % Denominador 2
num3 = [1 1]; % Numerador 3
den3 = [1 4 5 0 0]; % Denominador 3
[r1,p1,k1] = residue(num1,den1);
[r2,p2,k2] = residue(num2,den2);
[r3,p3,k3] = residue(num3,den3);
%% Resultados
disp('Sistema 1: [r p]')
disp(' ')
disp([r1 p1])
disp('Sistema 2: [r p]')
disp(' ')
disp([r2 p2])
disp('Sistema 3: [r p]')
disp(' ')
disp([r3 p3])
%--------------------------------------------------------------------------
El único aspecto con el que se debe tener cuidado en esta implementación es cómo se definen
los vectores de los polinomios. Los coeficientes se listan de manera descendente en potencias
de s. Si el polinomio no cuenta con un término, su coeficiente es cero. Nótese que la longitud
del vector de coeficientes es el orden del polinomio más uno.
Sistema 1: [r p]
-0.1667 -3.0000
0.5000 -2.0000
-0.5000 -1.0000
0.1667 0
Sistema 2: [r p]
-0.2500 -2.0000
0.5000 -2.0000
0.2500 0
Sistema 3: [r p]
102
Capítulo 3. Sistemas lineales continuos
Los resultados de los coeficientes y denominadores se colocaron en paralelo para que se vea la
correspondencia, pues el orden se mantiene. Hay que tener cuidado, pues la función residue
retorna las raíces del polinomio del denominador. Cuando se hace la expansión el denomina-
dor sería s − P (i). Otro aspecto para tener en cuenta es cuando se tienen raíces repetidas
(sistema 2). Para este caso, se observa que hay dos polos con valor de −2, entonces el primer
coeficiente corresponde a la menor potencia y el segundo a la multiplicidad.
Como complemento al programa, se puede validar los polos del sistema mediante el empleo
de la función roots.
%--------------------------------------------------------------------------
syms s
Fs = 'Función';
Ft = ilaplace(Fs);
pretty(Ft)
%--------------------------------------------------------------------------
Donde:
Ejemplo 3.10.
103
Capítulo 3. Sistemas lineales continuos
1
Y (s) = (3.18)
s4 + 6s3 + 11s2 + 6s
s+1
Y (s) =
s3 + 4s2 + 4s
s+1
Y (s) =
s4 + 4s3 + 5s2
Solución
Utilizando como base el esquema del programa para determinar la transformada de Laplace, se
plantea un solo programa para realizar todas transformaciones al dominio del tiempo. Inicial-
mente, se define la variable de Laplace “s” como simbólica; entonces, se definen las funciones
en un vector para posteriormente aplicar la antitransformada y, finalmente, se muestran los
resultados.
%--------------------------------------------------------------------------
% Obtención de antitransformadas de Laplace
% Método de solución: Matlab simbólico
clc
clear
syms s
Fs = [1/(s^4+6*s^3+11*s^2+6*s), (s+1)/(s^3+4*s^2+4*s),(s+1)/(s^4+4*s^3+5*s^2)];
Ft = ilaplace(Fs);
%--------------------------------------------------------------------------
Para esta implementación, hay que tener mucho cuidado con la ubicación de los paréntesis.
Una opción es usar la función pretty con las funciones planteadas para validar que quedaron
bien implementadas. Los resultados de ejecutar el programa son:
104
Capítulo 3. Sistemas lineales continuos
t exp(-2 t) exp(-2 t) 1
----------- - --------- + -
2 4 4
Estos resultados se pueden validar de manera sencilla usando tablas para obtener las anti-
transformadas de Laplace.
Como se ha visto en lo que va corrido de este documento, los sistemas se pueden representar
de diferentes maneras: partiendo de las ecuaciones diferenciales ordinarias, trasnformarlas en
ecuaciones diferenciales ordinarias lineales (en variables de desviación) y, posteriormente, con-
vertirlas a su representación lineal en espacio de estados o funciones de transferencia. Estas
transformaciones se denominan “realizaciones” de un modelo lineal. En el control de procesos
clásico es muy importante poder pasar de una realización del modelo a otra, pues cada una
ofrece unas ventajas en cuanto a uso o análisis de sistemas. En las siguientes secciones se espe-
cifican cómo construirlas y sus propiedades (estos modelos o realizaciones pueden representar
sistemas continuos o discretos; para aplicaciones de control de procesos se restringe solo a los
modelos continuos).
En Matlab, cada una de las representaciones posee un formato especial, es decir, un tipo de
estructura para definirlas y guardar información relevante asociada a estas. Para identificarlas,
cada modelo se define como objeto. Un objeto se puede entender como un contenedor de datos
especializado que encapsula un modelo e información relevante de este en una estructura. En
este sentido, representar modelos en forma de objetos permite manipular sistemas lineales de
una manera compacta en vez de manejar vectores, matrices o células.
Es importante definir los modelos como objetos; esto permite hacer transformaciones entre
estructuras, que es muy útil más adelante. Adicionalmente, los formatos de los modelos se
preservan para sus implementaciones en Simulink.
105
Capítulo 3. Sistemas lineales continuos
dx
= Ax + Bu (3.19)
dt
y = Cx + Du
Donde x, u y y son vectores. Para un sistema con n estados, m entradas y p salidas, se tiene
que las dimensiones de las matrices son:
A : n x n
B : n x m
C : p x n
D : p x m
%--------------------------------------------------------------------------
sys_ss = ss(A,B,C,D);
set(sys_ss,'Property1',Value1,...,'PropertyN',ValueN);
%--------------------------------------------------------------------------
Donde:
106
Capítulo 3. Sistemas lineales continuos
Y (s) 0 P olinomio(s)0
G(s) = =0 (3.20)
U (s) P olinomio(s)0
Donde:
Para crear un objeto que represente una función de transferencia, se utiliza la siguiente sin-
taxis:
%--------------------------------------------------------------------------
sys_tf = ft([numerador],[denominador]);
set(sys_ft,'Property1',Value1,...,'PropertyN',ValueN);
%--------------------------------------------------------------------------
Donde:
La función de transferencia se crea definiendo los polinomios del numerador y del denomi-
nador como vectores en orden descendente de potencias de s. Para crear la representación
de funciones de transferencia de sistemas MIMO, se puede concatenar funciones de transfe-
rencia individuales o utilizar argumentos como células de Matlab. Para el caso de células, el
numerador y el denominador son células con tantas filas como salidas y tantas columnas como
entradas. Entonces los vectores en el numerador {i.j} y en el denominador {i,j} se especifican
para la entrada j y la salida i (si todos los denominadores son iguales, se puede definir co-
mo denominador común) (MatWorks, 2022w). Las propiedades que pueden ser definidas son
análogas al caso de espacio de estados.
3.5.1.3. Cero-polo
La representación de un modelo en cero-polo es un arreglo de una función de transferencia
donde los polinomios se expresan de manera factorizada en términos de su ganancia, ceros y
107
Capítulo 3. Sistemas lineales continuos
Donde:
Para crear una representación en objeto cero-polo, se deben especificar los ceros, los polos y
la ganancia del sistema usando la siguiente sintaxis:
%--------------------------------------------------------------------------
sys_zp = zpk(Z,P,K)
set(sys_zp,'Property1',Value1,...,'PropertyN',ValueN)
%--------------------------------------------------------------------------
Donde:
Para sistemas MIMO, se deben introducir los vectores como células. En este orden de ideas,
Z y P son vectores dentro de células con tantas filas como salidas y tantas columnas como
entradas. Los vectores Z{i, j} y P {i, j} especifican los ceros y los polos para la función de
transferencia que relacionan la entrada j con la salida i. Por otro lado, K debe ser una
matriz con tantas filas como salidas y tantas columnas como entradas. Entonces, K(i, j) es la
ganancia de la función de transferencia de la entrada j con respecto a la salida i (MatWorks,
2022z).
108
Capítulo 3. Sistemas lineales continuos
dx
= Ax + Bu (3.22)
dt
y = Cx + Du
%--------------------------------------------------------------------------
sys_tf = ss2tf(A,B,C,D,ui);
sys_zp = ss2zp(A,B,C,D,ui);
sys_ss = tf2ss(num,den);
sys_zp = tf2zp(num,den);
sys_tf = zp2tf(z,p,k);
sys_ss = zp2ss(z,p,k);
Sys_tf = tf(sys);
Sys_zp = zp(sys);
Sys_ss = ss(sys);
%--------------------------------------------------------------------------
Donde:
Las primeras seis funciones son útiles cuando se cuenta con los valores desglosados de matrices,
numeradores/denominadores o vectores de ceros/polos y matriz de ganancias. Las siguientes
tres funciones son útiles cuando se cuenta con el modelo ya definido como objeto, bien sea en
espacio de estados, función de transferencia o cero-polo. En el ejemplo 3.11 se muestra cómo
aprovechar transformaciones.
109
Capítulo 3. Sistemas lineales continuos
Ejemplo 3.11.
Partiendo de la linealización del modelo del reactor Van de Vusse, obtener la representación
en funciones de transferencia y cero-polo.
Solución
%--------------------------------------------------------------------------
% Transformaciones entre modelos lineales
% Sistema de ecuaciones en el mismo script-funciones separadas
% Método de solución: analítico versus numérico
clc
clear
close all
[A,B,C,D]=linmod('VDV_NLsimulation_ports',Css,Uss); % Linealización
Sys_ss=ss(A,B,C,D,'StateName',{'Ca';'Cb'},'InputName',{'D';'Caf'},...
'OutputName',{'Ca';'Cb'}) % Crea espacio de estados
110
Capítulo 3. Sistemas lineales continuos
for i=1:length(Uss)
[num,den]=ss2tf(A,B,C,D,i); % Transforma el modelo
for j=1:length(Css)
G{j,i}= tf(num(j,:),den); % Genera el objeto de FT
disp(' ')
disp(['Función de transferencia: G(' num2str(i) ',' num2str(j) ')'])
G{j,i}
end
end
%--------------------------------------------------------------------------
Sys_ss =
A =
Ca Cb
Ca -2.405 0
Cb 0.8333 -2.238
B =
D Caf
Ca 7 0.5714
Cb -1.117 0
C =
Ca Cb
Ca 1 0
Cb 0 1
D =
D Caf
Ca 0 0
Cb 0 0
111
Capítulo 3. Sistemas lineales continuos
7 s + 15.67
---------------------
s^2 + 4.643 s + 5.382
-1.117 s + 3.147
---------------------
s^2 + 4.643 s + 5.382
ans =
0.5714 s + 1.279
---------------------
s^2 + 4.643 s + 5.382
0.4762
---------------------
s^2 + 4.643 s + 5.382
Sys_tf =
7
Ca: ---------
s + 2.405
-1.117 s + 3.147
Cb: ---------------------
s^2 + 4.643 s + 5.382
0.5714
Ca: ---------
112
Capítulo 3. Sistemas lineales continuos
s + 2.405
0.4762
Cb: ---------------------
s^2 + 4.643 s + 5.382
Reto 2: utilizando el teorema del valor final, validar que la ganancia dada por la
representación cero-polo corresponde a la determinada desde la función de transferencia.
113
Capítulo 4
Resumen
En el análisis de sistemas lineales se utiliza una serie de herramientas para anticipar el
comportamiento del sistema y poder usar esa información para el diseño y el control de
procesos. Teniendo en cuenta la pluralidad de los usuarios de control, se han desarrollado
diversas aproximaciones en el dominio del tiempo, en el dominio de Laplace y en el dominio
de la frecuencia para clasificar, cualificar y cuantificar el comportamiento de sistemas
lineales. Desde el punto de vista del control clásico de procesos, se utilizan respuestas del
sistema en el dominio del tiempo, así como la inspección de funciones de transferencia
para analizar sistemas. En esta sección se revisarán herramientas útiles para simular y
analizar el sistema lineal del reactor Van de Vusse. Cabe resaltar que el análisis de sistemas
lineales es un trabajo de interpretación de parámetros obtenidos de las representaciones;
sin embargo, es importante mostrar las herramientas para la extracción de la información
y ayudas gráficas, que son el fundamento para el análisis por realizar. En este capítulo
también se muestra cómo realizar simulaciones de sistemas lineales tanto en Matlab como
en Simulink y validar la aproximación lineal del sistema. Además, aborda cómo realizar
procesos de reducción de orden de modelos usando varias técnicas. Es importante tener
presente que el análisis de sistemas en el dominio de la frecuencia (i.e., Bode, Nyquist y
Nichols) está más allá del alcance de este documento.
4.1. Introducción
En este capítulo se mostrará cómo aprovechar el modelo lineal para determinar tres caracte-
rísticas fundamentales de los sistemas: estabilidad, controlabilidad y observabilidad, para lo
cual se utilizan diferentes opciones gráficas para representar los sistemas. Posteriormente, se
presentarán implementaciones en código de Matlab, herramientas del toolbox de control de
procesos y opciones en Simulink para evaluar la respuesta de sistemas lineales ante entradas
convencionales. Las implementaciones en Simulink se pueden acoplar con las simulaciones
realizadas en el capítulo anterior para así poder corroborar el poder predictivo de las apro-
ximaciones lineales. Finalmente, se exploran opciones para reducir el tamaño de los modelos
lineales y validar sus respuestas.
115
Capítulo 4. Simulación y análisis de sistemas lineales
4.2. Estabilidad
Partiendo de cualquiera de las representaciones del modelo, lo primero que hay que verificar
es su estabilidad. Debido al fuerte componente de métodos numéricos usados para la solución
de modelos dinámicos, se deben separar dos conceptos de estabilidad, para poderse enfocar en
la estabilidad inherente del sistema sin desconocer potenciales problemas numéricos durante
la solución.
• Error de redondeo: tiene que ver con el número de cifras significativas usadas en los
cálculos. Puede ser un problema para integraciones con cientos o miles de operaciones que
ocasiona desviaciones de la solución. Este problema se minimiza usando computadoras
modernas, con gran capacidad de procesamiento, ya que manejan doble precisión.
• Propagación del error: es preciso tener presente que los errores de truncamiento y re-
dondeo se acumulan durante la solución, en algunos casos de manera exponencial u
oscilatoria, lo que ocasiona que la solución numérica se aleje drásticamente de la real. A
este error se le conoce como error global y normalmente es un orden de magnitud mayor
al error local.
116
Capítulo 4. Simulación y análisis de sistemas lineales
La estabilidad numérica es algo que se presenta cuando se resuelven las ecuaciones diferencia-
les de manera numérica, y se apoya en una buena selección del método de integración como
se mostró en la sección 2.2. Para el análisis de sistemas se recomienda prestar atención a
problemas numéricos del método de solución en una etapa temprana, con el fin de enfocar
esfuerzos en la identificación de la estabilidad inherente del sistema.
Ya sea desde el punto de vista experimental o de la simulación del modelo no lineal, la es-
tabilidad BIBO se puede verificar aplicando una perturbación acotada (e.g., entrada paso,
onda, rampa acotada) y evaluando la respuesta del sistema si es acotada. Sin embargo, es
importante poder determinar la estabilidad del sistema desde el propio modelo matemático.
Se tiene entonces que para que un estado sea estable en un tiempo definido, se cumple que la
derivada de la función con respecto al tiempo es cero en ese punto. No obstante, lo anterior no
indica que sucede en la vecindad de este estado estable y, por tanto, si este punto de equilibrio
es de naturaleza estable o inestable.
Para un sistema con un solo estado, cuando se evalúa la derivada de la ecuación diferencial
con respecto al estado y se generan tres situaciones: si el valor es negativo, el sistema es
estable; si el valor es positivo, el sistema es inestable, y si la derivada vale cero, el sistema
puede ser estable o inestable (i.e., crítica o marginalmente estable). Cuando se tiene un sis-
tema con múltiples estados (i.e., múltiples ecuaciones diferenciales) se extrapola este concepto.
Partiendo del modelo no lineal, se realiza primero la linealización. La matriz del jacobiano A
corresponde a la derivada de las ecuaciones diferenciales con respecto a los estados. Se tiene
entonces que el sistema será estable si los valores propios (eigenvalores) de esa matriz tienen
parte real negativa. Se puede demostrar que los valores propios de la matriz A corresponden
a las raíces del polinomio característico de la función de transferencia; entonces, sus valores se
pueden relacionar directamente con la dinámica del sistema tanto en la velocidad de respues-
ta como en la forma de esta. El anterior es un criterio de estabilidad local y, por tanto, solo
es válido alrededor del punto donde fue determinado. Por otra parte, algo que este método
no indica es el tamaño del vecindario donde la estabilidad es válida; para esto se requieren
criterios de estabilidad más robustos como Lyapunov.
117
Capítulo 4. Simulación y análisis de sistemas lineales
Para calcular los valores propios de la matriz A o los polos del sistema (dado como polinomio
característico o como objeto de Matlab), se puede usar una de las siguientes funciones (todas
equivalentes):
%--------------------------------------------------------------------------
Val_propios = eig(A); % Valores propios de una matriz
Raices = roots(den); % Raíces del polinomio característico (den)
Polos = pole(sys); % Polos del sistema sys
%--------------------------------------------------------------------------
Donde:
%--------------------------------------------------------------------------
Sys = tf(num,den); % Definición del sistema SISO
pzmap(Sys); % Diagrama
sgrid % grid-opcional
%--------------------------------------------------------------------------
Como detalle especial, para este tipo de plano se puede usar la opción de líneas de guía con el
comando sgrid, el cual dibuja las líneas de igual factor de amortiguamiento y frecuencia na-
tural. Esta herramienta fue construida para sistemas SISO (una entrada-una salida), entonces
cuando se tiene un sistema MIMO (múltiples entradas-múltiples salidas) en representación de
espacio de estados, la herramienta solo grafica los polos. Para ver los ceros también se deben
obtener primero las funciones de transferencia individuales para luego graficarlas. De manera
alternativa se puede usar la función iopzmap para ver cada diagrama cero-polo independien-
temente para cada combinación entrada-salida. Para mostrar la identificación de polos y su
representación gráfica, se usa el ejemplo 4.1.
Ejemplo 4.1.
118
Capítulo 4. Simulación y análisis de sistemas lineales
x
x˙1 −4 0,3 1,5 0
1
x2
x˙2 = 1,2
−2 1
+ 1 u (4.1)
x˙3 −0,5 2 −3,5 x3 0
x1
y1 1 0 0 0
x2
y2 = 0 1
0
+
u
0
y3 0 0 1 x3 0
Solución
Para realizar esta tarea, se empieza con la definición del sistema en espacio de estados basado
en las matrices dadas. Posteriormente, todo el sistema se transforma en funciones de transfe-
rencia, mientras los valores propios de la matriz A se calculan directamente. Por otra parte,
del objeto de funciones de transferencia se extrae el denominador y se calculan las raíces.
Para validar lo anterior, se puede usar la función pole. Todos los resultados con comparados.
Finalmente, se construye el diagrama cero-polo a partir de las funciones de transferencia. El
código que realiza la tarea es:
%--------------------------------------------------------------------------
% Programa para validar la estabilidad de un sistema
% Desde espacio de estados y funciones de transferencia
clc
clear
close all
%% Modelos
%% Resultados
119
Capítulo 4. Simulación y análisis de sistemas lineales
pzmap(Sys_tf(1),Sys_tf(2),Sys_tf(3))
legend('FT 1','FT 2','FT 3')
ylabel('Eje Imaginario','FontWeight','bold','FontSize',14)
xlabel('Eje Real','FontWeight','bold','FontSize',14)
set(gca,'FontSize',12,'FontWeight','bold')
%--------------------------------------------------------------------------
El resultado es:
Sys_tf =
2 s + 7.85
3: -------------------------------
s^3 + 9.5 s^2 + 27.39 s + 16.79
120
Capítulo 4. Simulación y análisis de sistemas lineales
Nótese que el formato por defecto de la figura de pzmap fue modificado manualmente usando
el editor de propiedades, para de esta manera resaltar los ceros y los polos.
ẋ = Ax + Bu (4.2)
y = Cx + Du
121
Capítulo 4. Simulación y análisis de sistemas lineales
h i
Φc = B AB A2 B ··· An−1 B (4.3)
%--------------------------------------------------------------------------
Co = ctrb(A,B); % Definiendo las matrices
Co = ctrb(SYS); % Alimentando el sistema como objeto
Cont = length(A) - rank(Co); % Determinación de si Co es full rank
%--------------------------------------------------------------------------
Para sistemas muy grandes, en ocasiones el cálculo del rango de la matriz de controlabili-
dad tiene problemas e indica que el sistema es incontrolable sin serlo. Esto puede suceder
en sistemas como columnas de destilación, donde la matriz A es muy grande y muchos de
sus elementos son cero. Para estos casos, se recomienda validar la controlabilidad del sistema
usando la matriz de controlabilidad en staircase form a través del comando ctrbf.
4.3.2. Observabilidad
Un sistema es observable si, para una secuencia de estados y entradas, el estado actual del
sistema se puede determinar en un tiempo finito usando solo las salidas medidas. Este principio
se utiliza para el diseño de observadores o softsensors. Partiendo de la representación en
espacio de estados, se puede construir la matriz de observabilidad:
122
Capítulo 4. Simulación y análisis de sistemas lineales
C
CA
CA
2
Φo = (4.4)
..
.
n−1
CA
%--------------------------------------------------------------------------
Ob = obsv(A,C); % Definiendo las matrices
Ob = obsv(SYS); % Alimentando el sistema como objeto
Obsv = length(A) - rank(Ob);
%--------------------------------------------------------------------------
123
Capítulo 4. Simulación y análisis de sistemas lineales
2. impulse(sys): esta función permite obtener la respuesta impulso del sistema o grafi-
carla. Para sistemas donde la matriz D no es cero, la respuesta infinita del sistema para
el tiempo cero es ignorada.
Ejemplo 4.2.
Determinar la respuesta dinámica del sistema lineal Van de Vusse ante perturbación paso de
amplitud 1/7 en cada entrada, perturbación impulso y perturbación paso compuesta, donde
la tasa de dilución tiene una entrada paso de amplitud D = 1/7 (min−1 ) en el tiempo cero que
se combina con una CAf = 0,5 (mon/min) en el tiempo 10 min (i.e., retrasada en el tiempo).
Solución
Partiendo del modelo lineal estimado analíticamente para el reactor Van de Vusse, se crea
inicialmente un objeto en espacio de estados para el sistema. Seguidamente, se procede a usar
directamente las funciones especificando las características de las señales. La primera entrada
será simulada con la función step, la segunda con la función impulse y la última con la
función lsim, donde la entrada u se construye a partir de vectores (a esto se le llama series
de tiempo). El código generado es:
%--------------------------------------------------------------------------
% Simulación de respuestas de sistemas lineales
% Implementación en Matlab
clc
clear
close all
% Estado estable
124
Capítulo 4. Simulación y análisis de sistemas lineales
a11=-D-k1-2*k3*Cao;
a12=0;
a21=k1;
a22=-D-k2;
b11=Cafo-Cao;
b12=D;
b21=-Cbo;
b22=0;
Sys_ss=ss(A,B,C,D,'StateName',{'Ca';'Cb'},'InputName',{'D';'Caf'},...
'OutputName',{'Ca';'Cb'}); % Crea espacio de estados
%% Respuesta a perturbaciones
options = stepDataOptions('InputOffset',0,'StepAmplitude',1/7);
step(Sys_ss,[],options); % Respuesta paso
figure
impulse(Sys_ss) % Respuesta impulso
time = linspace(0,20);
u(:,1)= 1/7*heaviside(time);
u(:,2)= 0.5*heaviside(time-10);
figure
lsim(Sys_ss,u,time) % Respuesta a la entrada compuesta
%--------------------------------------------------------------------------
En esta aplicación se modifican las características del modelo lineal para nombrar las entra-
das, estados y salidas, a través de la función ss. Con la función stepDataOptions se cambian
las propiedades de la entrada paso, definiendo que su valor inicial es cero (InputOffset) y
que la amplitud del paso es D = 1/7 min−1 (StepAmplitude). Finalmente, nótese que en la
construcción de la señal compuesta se empleó la función heaviside para generar el paso con
retraso en el tiempo. Los resultados del programa son (figuras 4.2, 4.3 y 4.4):
125
Capítulo 4. Simulación y análisis de sistemas lineales
Figura 4.2. Respuesta del sistema lineal Van de Vusse ante una perturbación escalón en cada entrada
Figura 4.3. Respuesta del sistema lineal Van de Vusse ante una perturbación impulso en cada entrada
126
Capítulo 4. Simulación y análisis de sistemas lineales
Figura 4.4. Respuesta del sistema lineal Van de Vusse ante una perturbación paso de cada entrada
a diferente tiempo (para esta función, las señales de entrada también se muestran en la figura por
defecto)
Es importante recordar que para las funciones de transferencia y espacio de estados, las
variables tanto de entrada como de salida están en variables de desviación.
Con el fin de llevar a cabo el análisis en el dominio de la frecuencia, Matlab ofrece una serie de
funciones para que una vez definido el sistema se generen las respuestas para analizar. Como
se mencionó, el análisis en el dominio de la frecuencia normalmente no es prioritario en el
control clásico de procesos químicos. Sin embargo, se muestran las funciones como información:
127
Capítulo 4. Simulación y análisis de sistemas lineales
de frecuencias w. El rango de frecuencias debe ser un número real en radianes por unidad
de tiempo (por defecto en segundos). Cuando se tienen sistemas identificados, la función
también calcula las matrices de convarianza 2 × 2 (covH(i,j,k,:,:)) para la entrada
i-ésima a la salida j-ésima a la frecuencia w(k).
Todas estas funciones permiten realizar el procedimiento para varios sistemas MIMO de
manera simultánea.
Para usar la herramienta es necesario alimentar el modelo que se va a analizar. Este modelo
debe estar disponible en cualquiera de las representaciones para modelos lineales mostradas
con anterioridad. El modelo puede estar en el workspace de Matlab o ser cargado por un
archivo; para esto se debe ir al menú File/import y cargar el modelo. Una vez hecho lo
anterior, la herramienta por defecto muestra la respuesta paso del sistema. Para realizar las
otras gráficas, ir a Edit/Plot configurations y seleccionar qué tipo de respuestas se de-
sean analizar y cómo distribuirlas en las ventanas.
128
Capítulo 4. Simulación y análisis de sistemas lineales
Ejemplo 4.3.
Determinar la respuesta dinámica del sistema lineal Van de Vusse ante una perturbación paso
compuesta, donde la tasa de dilución tiene una entrada paso de amplitud D = 1/7 (min−1 )
en el tiempo cero que se combina con una CAf = 0,5 (mon/min) en el tiempo 10 min (i.e.,
retrasada en el tiempo).
Solución
Partiendo del modelo lineal estimado para el reactor Van de Vusse, se crea un objeto en espa-
cio de estados para el sistema como se hizo en el ejercicio 4.2. Posteriormente, se crea la señal
de entrada deseada y el vector de tiempo en el workspace. En la herramienta de Linear Sys-
tema Analyzer se carga el modelo y la señal de entrada como se muestra en las figuras 4.6 y 4.7.
129
Capítulo 4. Simulación y análisis de sistemas lineales
Figura 4.6. Ventana de importación del modelo Van de Vusse en espacio de estados desde el workspace
Figura 4.7. Ventana de importación de señal de entrada y construcción del vector de tiempo
Para importar el modelo a la herramienta, se usó el código del ejercicio de simulación de mo-
delos lineales de la sección 4.4.1 (ejemplo 4.2). En ese mismo ejemplo se generó la matriz de
señales de entrada u; es importante recalcar que el vector de tiempo se construyó directamente
en la herramienta. Adicionalmente, la herramienta ofrece una forma de construir señales de
entrada usando el botón Design signal, disponible en la sección de importar señal. Cuando se
simula el sistema se tiene como respuesta (figura 4.8):
130
Capítulo 4. Simulación y análisis de sistemas lineales
Figura 4.8. Respuesta del sistema Van de Vusse lineal ante la señal de entrada construida
Se puede validar que la respuesta obtenida es igual a la que arroja el ejemplo 4.2 para la
perturbación en cuestión.
Los modelos lineales continuos que se generaron con anterioridad en el ambiente de Matlab
se encuentran en la librería de Simulink llamada Continuous. Para el modelo de espacio de
estados se utiliza el bloque State-Space, el cual se muestra en la figura 4.9.
Para el modelo de espacio de estados, se accede a la ventana de configuración del bloque con
131
Capítulo 4. Simulación y análisis de sistemas lineales
doble clic sobre este. Análogamente a como se construyó en Matlab, hay que definir las matri-
ces A,B,C y D. Es muy importante establecer las dimensiones de cada matriz, especialmente la
matriz D que por lo general es nula, lo que no quiere decir que D = 0; las dimensiones de esta
matriz se definen por el número de salidas y de entradas del sistema. Las matrices que definen
el espacio de estados se pueden escribir directamente en la ventana de configuración con la
misma sintaxis de Matlab o se les puede asingar un nombre para que se importen desde el
workspace. Es por esto que la conexión con Matlab es importante, pues un código en Matlab
las puede calcular para después importarlas en Simulink. Mientras las matrices existan en el
workspace, pueden ser leídas directamente por Simulink. Adicionalmente, es posible especificar
las condiciones iniciales y los nombres para los estados.
El espacio de estados fue diseñado para representar un sistema MIMO y, por tanto, tiene
múltiples entradas y salidas. Sin embargo, es importante notar que el bloque solo tiene una
entrada y una salida; esto implica que esos puertos son multivariable y se requieren los bloques
Mux para combinar entradas en un solo bus de datos y Demux para separar las salidas del
bloque.
132
Capítulo 4. Simulación y análisis de sistemas lineales
Por último, se tiene el bloque para el modelo en representación cero-polo (figura 4.11). Esta
representación es análoga a la función de transferencia, con la diferencia de que se deben
proveer los numeradores y el denominador como los ceros y polos y no los coeficientes de cada
polinomio.
En el ejemplo 4.4 se muestra cómo utilizar estos bloques para realizar la misma simulación
de perturbaciones en la tasa de dilución y concentración en la entrada del reactor Van de Vusse.
133
Capítulo 4. Simulación y análisis de sistemas lineales
Ejemplo 4.4.
Determinar la respuesta dinámica del sistema lineal Van de Vusse ante una perturbación en
la tasa de dilución de amplitud D = 1/7 (min−1 ) en el tiempo 5 min combinada con un paso
negativo de CAf = 2 (mol/l) en el tiempo 10 min.
Solución
Para esta implementación se utilizará como programa maestro un código en Matlab que lla-
mará a los modelos implementados en Simulink. La idea es simular los modelos lineales en
representación de espacio de estados y funciones de transferencia con las mismas entradas y
determinar la diferencia entre los dos. Cuando se tienen los resultados de la simulación, se usa
Matlab para graficar los resultados.
En primera instancia, se construye el modelo lineal en espacio de estados del reactor usando
la linealización analítica del sistema y se asignan nombres a las matrices para que puedan ser
reconocidas por el modelo en Simulink. Nótese que en este caso no se define el objeto de espa-
cio de estados usando ss. A continuación, se obtienen las funciones de transferencia a partir
del modelo en espacio de estados usando el comando de transformación ss2tf que se ejecuta
para cada entrada. Como resultado, se tienen los vectores de numeradores y denominadores
para cada función de transferencia que se guardan en formato de célula.
Por otro lado, se construye un programa en Simulink para evaluar la respuesta de los dos
modelos en paralelo. Para el modelo en espacio de estados, basta con introducir las matrices
A, B, C y D con sus nombres en la configuración del bloque de espacio de estados. Ahora,
para el modelo en funciones de transferencia, la situación es un poco más compleja pues se
tiene que interconectar cada función usando el principio de superposición de sistemas lineales
(esto se detalla más adelante). En Simulink se utilizan los bloques To workspace para enviar
las respuestas de los sistemas a Matlab y poder graficarlos usando el formato modificado para
tener figuras visualmente más agradables. Desde el programa maestro en Matlab se ejecuta
el programa de Simulink usando la función sim.
%--------------------------------------------------------------------------
% Programa para cargar los parámetros de los modelos lineales
% y simular respuesta del reactor de Van de Vusse
% Modelos en Simulink y programa maestro en Matlab
clc
clear
134
Capítulo 4. Simulación y análisis de sistemas lineales
close all
% Parámetros cinéticos
% Estado estable
a11=-D-k1-2*k3*Cao;
a12=0;
a21=k1;
a22=-D-k2;
b11=Cafo-Cao;
b12=D;
b21=-Cbo;
b22=0;
for i=1:length(B)
[num{i},den{i}]=ss2tf(A,B,C,D,i); % Transforma el modelo
end
%% Simulación y resultados
sim('VDV_Linear_simulation')
% Entradas
tiempo_L=CL_ss.Time;
Ca_ss=CL_ss.Data(:,1);% Serie de tiempo de modelo en EE
Cb_ss=CL_ss.Data(:,2);% Serie de tiempo de modelo en EE
135
Capítulo 4. Simulación y análisis de sistemas lineales
figure(1)
subplot(231)
plot(tiempo_In,D,'b','LineWidth',2)
ylabel('D [1/min]','FontWeight','bold','FontSize',14)
set(gca,'FontSize',12,'FontWeight','bold')
grid
subplot(232)
hold on
plot(tiempo_In,Ca_ss,'b','LineWidth',2)
plot(tiempo_In,Ca_tf,'--r','LineWidth',2)
legend('Espacio de estados','Función de transferencia')
grid
ylabel('C_A [mol/l]','FontWeight','bold','FontSize',14)
set(gca,'FontSize',12,'FontWeight','bold')
subplot(233)
plot(tiempo_In,Err_Ca,'b','LineWidth',2)
ylabel('Error C_A [mol/l]','FontWeight','bold','FontSize',14)
set(gca,'FontSize',12,'FontWeight','bold')
grid
subplot(234)
plot(tiempo_In,Caf,'b','LineWidth',2)
ylabel('Caf [mol/l]','FontWeight','bold','FontSize',14)
xlabel('Tiempo [min]','FontWeight','bold','FontSize',14)
set(gca,'FontSize',12,'FontWeight','bold')
grid
subplot(235)
hold on
plot(tiempo_L,Cb_ss,'b','LineWidth',2)
plot(tiempo_L,Cb_tf,'--r','LineWidth',2)
legend('Espacio de estados','Función de transferencia')
grid
ylabel('C_B [mol/l]','FontWeight','bold','FontSize',14)
xlabel('Tiempo [min]','FontWeight','bold','FontSize',14)
set(gca,'FontSize',12,'FontWeight','bold')
subplot(236)
plot(tiempo_In,Err_Cb,'b','LineWidth',2)
ylabel('Error C_B [mol/l]','FontWeight','bold','FontSize',14)
xlabel('Tiempo [min]','FontWeight','bold','FontSize',14)
set(gca,'FontSize',12,'FontWeight','bold')
grid
%--------------------------------------------------------------------------
136
Capítulo 4. Simulación y análisis de sistemas lineales
Figura 4.12. Implementación en Simulink del modelo lineal del reactor Van de Vusse en espacio de
estados y funciones de transferencia
Para crear un subsistema Linear model (TF), primero se construyen las cuatro funciones de
transferencia y sus interconexiones mostradas en la figura 4.12 dentro del programa principal
de Simulink. Para introducir los numeradores y denominadores en cada función de transfe-
rencia, es preciso tener en cuenta que en el programa de Matlab las funciones de transferencia
fueron guardadas en células llamadas num y den. Por tanto, en los bloques de función de trans-
ferencia de la figura 4.13, los vectores de numeradores y denominadores deben ser llamados
como células.
Entonces, se seleccionan las funciones de transferencia que se desea colocar dentro del sub-
sistema y se da clic derecho sobre estas. Se busca en la ventana emergente de opciones la
llamada Create subsystem from selection. De esta manera, las funciones de transferen-
cia mostradas en la figura 4.13 desaparecen de la ventana principal y se reemplazan con un
solo bloque. Nótese que la máscara tiene las entradas y salidas con puertos llamados Inport y
Outport, los cuales se generan automáticamente cuando se crea la máscara. En consecuencia,
en el programa principal de Simulink solo queda un bloque representando las funciones de
transferencia.
137
Capítulo 4. Simulación y análisis de sistemas lineales
Figura 4.13. Implementación de las funciones de transferencia dentro de la máscara llamada Linear
model (TF)
Cuando se ejecuta el programa de Matlab, se obtienen los siguientes resultados (figura 4.14):
Figura 4.14. Comparación de las repuestas del sistema lineal Van de Vusse con modelos en espacio
de estados y funciones de transferencia
138
Capítulo 4. Simulación y análisis de sistemas lineales
Ejemplo 4.5.
Determinar la respuesta dinámica del sistema no lineal y lineal para el reactor Van de Vusse
ante una perturbación en la tasa de dilución de amplitud D = 1/7 (min−1 ) en el tiempo 5
min combinada con un paso negativo de CAf = 2 (mol/l) en el tiempo 10 min.
Solución
Para esta implementación se utilizará Simulink como programa maestro. La idea es simular
ambos modelos con las mismas entradas y determinar la diferencia entre los dos, tomando
como solución verdadera la salida del modelo no lineal. Cuando se tienen los resultados de la
simulación, se usa Matlab para realizar las figuras utilizando un Callback. El modelo no lineal
se recicla de la implementación de la S-function y el modelo lineal se implementa usando
el bloque de espacio de estados. La implementación se muestra en la figura 4.15. Nótese que,
para el modelo lineal, se transforman las entradas a variables de desviación y las salidas a
valores nominales. Otro aspecto que se implementa en esta simulación es que los parámetros
requeridos se cargan desde el workspace usando un Callback pero con la opción InitFcn, para
que corra el programa de los parámetros antes de ejecutar la simulación.
139
Capítulo 4. Simulación y análisis de sistemas lineales
Figura 4.15. Implementación del sistema lineal y no lineal para el reactor Van de Vusse en Simulink
El programa de Matlab para cargar los parámetros del modelo de espacio de estados antes de
correr la simulación de Simulink es:
%--------------------------------------------------------------------------
% Programa para cargar los parámetros de la función de transferencia
% Reactor de Van de Vusse
% Parámetros cinéticos
% Estado estable
a11=-D-k1-2*k3*Cao;
a12=0;
a21=k1;
a22=-D-k2;
b11=Cafo-Cao;
b12=D;
b21=-Cbo;
b22=0;
140
Capítulo 4. Simulación y análisis de sistemas lineales
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% Programa para hacer las figuras
close all
tiempo_L=C_SS.Time;
Ca_ML=C_SS.Data(:,1);
Cb_ML=C_SS.Data(:,2);
% Entradas
% Error
%% Resultados
figure(1)
subplot(231)
plot(tiempo_In,D,'b','LineWidth',2)
ylabel('D[1/min]','FontWeight','bold','FontSize',14)
set(gca,'FontSize',12,'FontWeight','bold')
grid
subplot(232)
hold on
plot(tiempo,Ca,'b','LineWidth',2)
plot(tiempo_L,Ca_ML,'--r','LineWidth',2)
legend('Modelo no lineal','Modelo lineal')
grid
ylabel('C_A[mol/l]','FontWeight','bold','FontSize',14)
set(gca,'FontSize',12,'FontWeight','bold')
141
Capítulo 4. Simulación y análisis de sistemas lineales
subplot(233)
plot(tiempo,ErrCa,'b','LineWidth',2)
grid
ylabel('Error Ca [ %]','FontWeight','bold','FontSize',14)
set(gca,'FontSize',12,'FontWeight','bold')
subplot(234)
plot(tiempo_In,Caf,'b','LineWidth',2)
ylabel('Caf[mol/l]','FontWeight','bold','FontSize',14)
xlabel('Tiempo [min]','FontWeight','bold','FontSize',14)
set(gca,'FontSize',12,'FontWeight','bold')
grid
subplot(235)
hold on
plot(tiempo,Cb,'b','LineWidth',2)
plot(tiempo_L,Cb_ML,'--r','LineWidth',2)
legend('Modelo no lineal','Modelo lineal')
grid
ylabel('C_B[mol/l]','FontWeight','bold','FontSize',14)
xlabel('Tiempo [min]','FontWeight','bold','FontSize',14)
set(gca,'FontSize',12,'FontWeight','bold')
subplot(236)
hold on
plot(tiempo,ErrCb,'b','LineWidth',2)
grid
ylabel('Error Cb [ %]','FontWeight','bold','FontSize',14)
xlabel('Tiempo [min]','FontWeight','bold','FontSize',14)
set(gca,'FontSize',12,'FontWeight','bold')
%--------------------------------------------------------------------------
Al comparar las repuestas de ambos sistemas, se obtienen los siguientes resultados para las
entradas, estados y el error del modelo lineal (figura 4.16):
142
Capítulo 4. Simulación y análisis de sistemas lineales
Figura 4.16. Respuesta de los modelos lineal y no lineal para el reactor Van de Vusse ante diferentes
perturbaciones
• Tiempo de elevación (tr, Rise time): tiempo que se tarda la respuesta en alcanzar por
primera vez el nuevo valor de estado estable.
• Tiempo pico (tp, Time to first peak): tiempo que tarda la respuesta en alcanzar su
primer valor máximo.
• Tiempo de asentamiento (ts, Settling time): tiempo requerido para que la respuesta
alcance y permanezca en la banda de ±5 % de su valor final (±1 % también es utilizado).
• Sobreimpulso (OS, Overshoot: OS = a/b): razón entre la cantidad que excede la res-
puesta en el primer pico y el valor del estado estable.
143
Capítulo 4. Simulación y análisis de sistemas lineales
• Razón de asentamiento (DR, Decay ratio: DR = c/a): razón entre la cantidad que
excede la respuesta en el segundo pico y el exceso del primero (exceso del valor de
estado estable).
• Periodo de oscilación (P, Period of oscillation): tiempo entre dos picos (o valles) conse-
cutivos.
Figura 4.17. Respuesta típica de un sistema de segundo orden subamortiguado donde se indican los
parámetros de desempeño a identificar. [Fuente: adaptado de Seborg et al. (2011)]
Cuando los sistemas no se comportan como uno de segundo orden amortiguado, solo algunos
de los parámetros anteriores pueden ser determinados. Ahora, si bien en la literatura hay
expresiones analíticas para determinar los parámetros en función de las constantes de la
función de transferencia (i.e., τ y ζ), estas solo son válidas si el sistema es propiamente de
segundo orden o se ha aproximado a uno de esa naturaleza (Seborg et al., 2011). Por tanto,
aún es relevante poder extraer estos parámetros directamente de la curva de respuesta.
Para facilitar esta tarea, Simulink ha incorporado una serie de herramientas en la ventana de
visualización de resultados a través del bloque Scope, usando el panel Screen cursors. Las
herramientas son:
• Cursor Measurements panel: permite tomar distancias entre puntos de alguna de las
respuestas graficadas en la figura. Esta herramienta ofrece dos opciones de medición:
Screen Cursors y Waveform Cursors. En Screen Cursors, se grafican dos cursores (x, y)
como punteros y la herramienta toma medidas entre estas dos coordenadas. Para el modo
Waveform Cursors, unos cursores verticales aparecen sobre la figura, los cuales pueden
moverse y en el panel de respuesta se dan las distancias entre los cursores y también la
ubicación de estos. Como opciones permite también fijar la diferencia en frecuencias de
144
Capítulo 4. Simulación y análisis de sistemas lineales
ambos cursores o sus posiciones. Esta herramienta es muy útil para identificar tiempo
de elevación, tiempo de asentamiento, tiempo de primer pico, periodo, sobreimpulso y
razón de asentamiento (MatWorks, 2022c).
• Peak finder panel: permite encontrar la localización de los picos en la señal graficada.
Los picos se definen como máximos locales en la figura (MatWorks, 2022o).
Nótese que, dentro de cada panel de herramientas, se pueden configurar las opciones de cada
una de estas. En el ejemplo 4.6 se muestra la visualización de cada panel.
Ejemplo 4.6.
1
G(s) = (4.5)
s2 + 0,5s + 1
Solución
145
Capítulo 4. Simulación y análisis de sistemas lineales
Figura 4.18. Implementación del sistema de segundo orden para evaluación de resultados
Figura 4.19. Respuesta del sistema de segundo orden ante perturbación paso, uso de la herramienta
Cursor Measurements panel
En esta figura, de manera manual se ubicaron los cursores en el punto del primer pico (sobre-
impulso) y cuando la señal entró al umbral del ±5 % del valor de la señal final. En el panel
de resultados se muestran la coordenadas de ambos puntos (figura 4.19).
146
Capítulo 4. Simulación y análisis de sistemas lineales
Figura 4.20. Respuesta del sistema de segundo orden ante perturbación paso, uso de la herramienta
Signal statistics panel
Figura 4.21. Respuesta del sistema de segundo orden ante perturbación paso, uso de la herramienta
Bilevel measurements panel
Finalmente, los resultados del panel de Peak finder ayuda a validar los valores de los picos,
su localización y determinar la razón de asentamiento (figura 4.22).
147
Capítulo 4. Simulación y análisis de sistemas lineales
Figura 4.22. Respuesta del sistema de segundo orden ante perturbación paso, uso de la herramienta
Peak finder panel
Reto: validar los resultados obtenidos usando las herramientas de medida de Simulink
con las ecuaciones analíticas para su estimación a partir de los parámetros de la función
de transferencia.
• Mientras se realiza el diagrama cero-polo del sistema, se evidencia que hay ceros cercanos
a polos.
148
Capítulo 4. Simulación y análisis de sistemas lineales
• Se diseña un controlador de alto orden y se desea implementar uno de bajo orden como
un PID. Un controlador de alto orden se diseña con métodos como el Gaussian lineal-
cuadrático o técnicas como la síntesis H∞.
Los métodos implementados en el Control System toolbox se basan en simplificación del mo-
delo por eliminación de segmentos (funciones 1 y 2), aproximación a otros modelos (función
3 y 4) o selección de modos (función 5) (MatWorks, 2022m). Entonces, se pueden obtener
diferentes representaciones finales de los modelos reducidos según el método utilizado. Debido
a esto, es importante validar cada modelo reducido con el modelo original para valorar las
consecuencias de la reducción. Normalmente, se utiliza el diagrama de Bode o la respuesta a
una entrada paso como se realizó en la sección 4.4.4.
En las siguientes secciones se abordan detalles de cada comando incluyendo el método usado
para realizar la reducción y las configuraciones especiales. Se anticipa que muchos de estos
conceptos trascienden el alcance de un curso introductorio a control de procesos, pero es im-
portante mencionarlos en aras de diferenciar aspectos de los métodos.
149
Capítulo 4. Simulación y análisis de sistemas lineales
%--------------------------------------------------------------------------
[sysr,u] = minreal(sys,tol)
sys_ss = ss(sys,'minimal')
%--------------------------------------------------------------------------
Donde:
Nótese que la segunda función que aparece en el código anterior es una forma alternativa de
obtener la realización mínima de un sistema en espacio de estados. Esta se lleva a cabo desde
la misma definición del modelo como un objeto en espacio de estados.
%--------------------------------------------------------------------------
msys = sminreal(sys)
%--------------------------------------------------------------------------
Donde:
El modelo que se obtiene con esta función no es necesariamente mínimo y puede tener mayor
orden que el que se obtiene con la función minreal. Sin embargo, usando sminreal la es-
tructura de estados de la función original se mantiene, mientras que con la función minreal
cambia.
150
Capítulo 4. Simulación y análisis de sistemas lineales
número mínimo de estados que se requieren. Los valores singulares de Hankel son una me-
dida de la contribución de cada estado en el comportamiento del sistema o energía de estos.
Entonces, estados asociados a bajos valores singulares de Hankel pueden ser eliminados para
simplificar el modelo preservando las características más importantes en términos de estabili-
dad, frecuencia y respuesta en el tiempo (MatWorks, 2022g). La sintaxis para usar la función
es:
%--------------------------------------------------------------------------
rsys = balred(sys,orders,[],opts)
%--------------------------------------------------------------------------
Donde:
La gráfica de los valores singulares de Hankel del sistema lineal invariante en el tiempo (LTI)
se obtiene con la función hsv = hsvd, con los argumentos (sys,options). Se usa la función
hsvdOptions(Name,Value) para configurar las tolerancias del error absoluto y relativo de
la descomposición entre estados estables e inestables, el intervalo de frecuencia o de tiempo
para calcular los valores singulares de Hankel o el offset para los límites estables e inestables
cumpliendo (Gawronski y Juang, 1990):
Para opciones personalizables se puede usar la función hsvplot para ajustar ejes, escala y
estilos de la figura en forma de estructura (MatWorks, 2022h).
Durante la reducción del modelo, si el sistema tiene polos inestables, la función primero
descompone el sistema en subsistemas (estable e inestable) y la aproximación a un modelo
reducido se hace con la parte estable. Para realizar modificaciones a las opciones de la aproxi-
mación se usa la función balredOptions. Las modificaciones incluyen el offset y la tolerancia
para la descomposición de la parte estable e inestable y hacen énfasis en los intervalos de fre-
cuencia de interés. Como opción especial, la función permite usar el método de truncamiento
opts = balredOptions('StateElimMethod','Truncate').
151
Capítulo 4. Simulación y análisis de sistemas lineales
La función encuentra una representación reducida donde las matrices de Gram de controlabi-
lidad y observabilidad son diagonales e iguales a las del sistema original. La diagonal de dichas
matrices es el vector de valores singulares de Hankel; por ende, bajos valores en el vector indi-
can que dichos estados pueden ser removidos sin afectar significativamente las características
dinámicas del sistema original. La función balreal se complementa con la función modred
para eliminar los estados de baja energía.
Dado el caso en que se tengan polos inestables, la función inicialmente descompone el sistema
en las partes estable e inestable. La parte estable es separada, balanceada y se añade a la par-
te inestable para retornar el sistema reducido. Cuando se tienen polos inestables, los valores
singulares de Hankel se muestran como infinito (inf). Análogo a la función para obtener los
valores singulares de Hankel, se utiliza la función options = hsvdOptions(Name,Value)
para configurar opciones de tolerancia, offset e intervalos de frecuencia y tiempo para calcular
las matrices de Gram.
%--------------------------------------------------------------------------
[sysb,g] = balreal(sys);
%--------------------------------------------------------------------------
Donde:
Esta función permite determinar un modelo reducido tras eliminar estados particulares de un
modelo lineal; por tanto, debe ser alimentada con un vector que indique la posición o índices
de los estados por eliminar. Esta función se usa mancomunadamente con la función balreal,
donde primero se realiza la reducción balanceada para aislar los estados de baja energía y luego
se eliminan con la función modred teniendo en cuenta el vector de valores singulares de Hankel.
%--------------------------------------------------------------------------
[sysb,g] = balreal(sys);
elim = (g<1e-8);
rsys = modred(sysb,elim,'method')
%--------------------------------------------------------------------------
Donde:
152
Capítulo 4. Simulación y análisis de sistemas lineales
Para reducir el modelo se tienen dos métodos: MatchDC y Truncate. El primer método es
seleccionado por defecto y garantiza tener la misma ganancia; el segundo elimina los estados
indicados directamente. Se espera tener una mejor respuesta en el dominio de la frecuencia
utilizando el método de truncamiento, sin embargo, no hay garantía de que la ganancia se
conserve (MatWorks, 2022n).
Para mostrar el uso de las funciones que permiten reducir modelos se plantea el ejemplo 4.7.
Ejemplo 4.7.
Determinar la respuesta dinámica del sistema lineal Van de Vusse y aproximaciones reducidas
del modelo lineal ante una perturbación en la tasa de dilución de amplitud D = 1/7 (min−1 )
en el tiempo 5 min combinada con un paso negativo de CAf = 2 (mol/l) en el tiempo 10 min.
Solución
Análogo a la simulación para validar el modelo lineal, para esta implementación se utiliza un
programa en Matlab como programa maestro que ejecuta los modelos lineales implementa-
dos en Simulink. La idea es simular todos modelos con las mismas entradas y determinar la
diferencia entre ellos. Cuando se tienen los resultados de la simulación, se usa Matlab para
generar las figuras.
El modelo lineal de referencia es la linealización analítica del sistema Van de Vusse en espa-
cio de estados. Se le asignan nombres a las matrices para que puedan ser reconocidas por el
modelo en Simulink. En este caso, para obtener modelos reducidos se utiliza la cancelación
cero-polo desde funciones de transferencia y la eliminación de estados usando realización ba-
lanceada.
Para llevar a cabo la cancelación cero-polo es necesario validar si hay ceros cercanos a po-
los en el sistema. Se realiza la transformación entonces de espacio de estados a funciones de
transferencia y se usa la función pzmap para evidenciarlo. Posteriormente, se aplica la función
minreal a las funciones de transferencia para hacer la cancelación. En cuanto a la reducción
balanceada del sistema, se utiliza para esto la función balreal aplicada directamente al mo-
delo en espacio de estados. A dicho sistema balanceado se le determinan los valores singulares
de Hankel y se elimina el estado de menor energía usando la función modred. Aprovechando
la oportunidad, se realizan ajustes de formato a la figura de valores singulares de Hankel.
153
Capítulo 4. Simulación y análisis de sistemas lineales
%--------------------------------------------------------------------------
% Programa para realizar diferentes realizaciones de un modelo lineal
% Caso de estudio del reactor de Van de Vusse
clc
clear
close all
% Parámetros cinéticos
% Estado estable
a11=-D-k1-2*k3*Cao;
a12=0;
a21=k1;
a22=-D-k2;
b11=Cafo-Cao;
b12=D;
b21=-Cbo;
b22=0;
%% Reducción de modelos
for i=1:length(B)
[num{i},den{i}]=ss2tf(A,B,C,D,i); % Transforma el modelo
for j=1:length(C)
154
Capítulo 4. Simulación y análisis de sistemas lineales
Sys_tf{j,i}=tf(num{i}(j,:),den{i});
Sys_tfr{j,i}= minreal(Sys_tf{j,i},tol);
end
end
figure
subplot(221)
pzmap(Sys_tf{1,1})
ylabel('Eje imaginario','FontWeight','bold','FontSize',14)
xlabel('Eje real','FontWeight','bold','FontSize',14)
set(gca,'FontSize',12,'FontWeight','bold')
subplot(222)
pzmap(Sys_tf{1,2})
ylabel('Eje imaginario','FontWeight','bold','FontSize',14)
xlabel('Eje real','FontWeight','bold','FontSize',14)
set(gca,'FontSize',12,'FontWeight','bold')
subplot(223)
pzmap(Sys_tf{2,1})
ylabel('Eje imaginario','FontWeight','bold','FontSize',14)
xlabel('Eje real','FontWeight','bold','FontSize',14)
set(gca,'FontSize',12,'FontWeight','bold')
subplot(224)
pzmap(Sys_tf{2,2})
ylabel('Eje imaginario','FontWeight','bold','FontSize',14)
xlabel('Eje real','FontWeight','bold','FontSize',14)
set(gca,'FontSize',12,'FontWeight','bold')
figure
[Sysb,g] = balreal(Sys_ss); % Realización balanceada
P = hsvoptions; % Ajuste de formato
P.Title.FontSize = 14; % Ajuste de formato
P.Title.FontWeight='bold'; % Ajuste de formato
P.XLabel.FontSize= 12; % Ajuste de formato
P.XLabel.FontWeight='bold'; % Ajuste de formato
P.YLabel.FontSize= 12; % Ajuste de formato
P.YLabel.FontWeight='bold'; % Ajuste de formato
P.TickLabel.FontSize= 12; % Ajuste de formato
P.TickLabel.FontWeight='bold'; % Ajuste de formato
h=hsvplot(Sysb,P); % Valores singulares de Hankel (VSH)
elim = (g<0.3); % Identificación de VSH menores a 0.3
Rsys = modred(Sysb,elim,'MatchDC'); % Eliminación de estados
%% Simulación y resultados
sim('VDV_Linear_simulation_MR')
% Entradas
tiempo_L=CL_ss.Time;
155
Capítulo 4. Simulación y análisis de sistemas lineales
figure
subplot(231)
plot(tiempo_In,D,'b','LineWidth',2)
ylabel('D [1/min]','FontWeight','bold','FontSize',14)
set(gca,'FontSize',12,'FontWeight','bold')
grid
subplot(232)
hold on
plot(tiempo_In,Ca_ss,'b','LineWidth',2)
plot(tiempo_In,Ca_tf,'--r','LineWidth',2)
plot(tiempo_In,Ca_bal,'--k','LineWidth',2)
legend('Modelo lineal','ZP-cancelación','Balreal')
grid
ylabel('C_A [mol/l]','FontWeight','bold','FontSize',14)
set(gca,'FontSize',12,'FontWeight','bold')
subplot(233)
hold on
plot(tiempo_In,Err_Ca_tf,'r','LineWidth',2)
plot(tiempo_In,Err_Ca_bal,'k','LineWidth',2)
legend('ZP-cancelación','Balreal')
ylabel('% Error C_A [mol/l]','FontWeight','bold','FontSize',14)
set(gca,'FontSize',12,'FontWeight','bold')
grid
subplot(234)
plot(tiempo_In,Caf,'b','LineWidth',2)
ylabel('Caf [mol/l]','FontWeight','bold','FontSize',14)
xlabel('Tiempo [min]','FontWeight','bold','FontSize',14)
set(gca,'FontSize',12,'FontWeight','bold')
grid
subplot(235)
hold on
plot(tiempo_L,Cb_ss,'b','LineWidth',2)
plot(tiempo_L,Cb_tf,'--r','LineWidth',2)
plot(tiempo_In,Cb_bal,'--k','LineWidth',2)
legend('Modelo lineal','ZP-cancelación','Balreal')
grid
ylabel('C_B [mol/l]','FontWeight','bold','FontSize',14)
xlabel('Tiempo [min]','FontWeight','bold','FontSize',14)
set(gca,'FontSize',12,'FontWeight','bold')
subplot(236)
156
Capítulo 4. Simulación y análisis de sistemas lineales
hold on
plot(tiempo_In,Err_Cb_tf,'r','LineWidth',2)
plot(tiempo_In,Err_Cb_bal,'k','LineWidth',2)
legend('ZP-cancelation','Balreal')
ylabel('% Error C_B [mol/l]','FontWeight','bold','FontSize',14)
xlabel('Tiempo [min]','FontWeight','bold','FontSize',14)
set(gca,'FontSize',12,'FontWeight','bold')
grid
Figura 4.23. Implementación en Simulink del modelo lineal del reactor Van de Vusse en espacio de
estados y los modelos reducidos por diferentes métodos
7 s + 15.67
---------------------
s^2 + 4.643 s + 5.382
157
Capítulo 4. Simulación y análisis de sistemas lineales
0.5714 s + 1.279
---------------------
s^2 + 4.643 s + 5.382
-1.117 s + 3.147
---------------------
s^2 + 4.643 s + 5.382
0.4762
---------------------
s^2 + 4.643 s + 5.382
Realizando el diagrama cero-polo para este sistema se obtiene el siguiente resultado (figura
4.24):
158
Capítulo 4. Simulación y análisis de sistemas lineales
Figura 4.24. Diagrama cero-polo para las funciones de transferencia del sistema completo Van de
Vusse
Como se observa en la figura, para las funciones de transferencia que relacionan cada entrada
con la concentración de A hay un cero superpuesto con un polo y, por ende, se puede cancelar.
Aplicando la función minreal a las funciones de transferencia se obtienen las versiones redu-
cidas (los formatos de la respuesta de Matlab se ajustaron para mostrar la correspondencia
entrada-salida):
7
---------
s + 2.405
0.5714
---------
s + 2.405
Para la realización balanceada del espacio de estados se determinan los valores singulares de
Hankel, de lo que se obtiene la siguiente figura de energía (figura 4.25):
159
Capítulo 4. Simulación y análisis de sistemas lineales
Figura 4.25. Valores singulares de Hankel para el sistema balanceado del reactor Van de Vusse
En principio, solo los estados de muy baja energía (tol=1x10^-8) deben ser eliminados. En
aras de realizar una eliminación para este sistema, se relaja el límite recomendado a tol=0.3
y se puede valorar las consecuencias de esta elección. El sistema reducido por el método de
balanceo entrada-salida es:
A =
x1
x1 -2.237
B =
u1 u2
x1 -2.593 -0.2187
C =
x1
y1 -2.419
y2 -0.9576
D =
u1 u2
y1 0.1061 0.001044
y2 -0.5254 -0.005166
Como se puede ver, el sistema reducido solo tiene un estado. Cabe resaltar que ahora la
matriz C ya no es la identidad y la matriz D tiene valores diferentes a cero. Finalmente, la
comparación de los modelos se hace simulándolos en Simulink y se obtienen los siguientes
resultados (figura 4.26):
160
Capítulo 4. Simulación y análisis de sistemas lineales
Figura 4.26. Evaluación del desempeño de los modelos reducidos comparados con la representación
analítica del sistema Van de Vusse en espacio de estados
De estos resultados se evidencia que la diferencia entre el modelo no lineal y el obtenido por
la cancelación cero-polo es casi imperceptible; al validarlos internamente, estos son del orden
de magnitud de error numérico. Para la segunda reducción, hay una pérdida de desempeño
del modelo debido a que el estado eliminado tenía una contribución relevante en la salida, lo
cual no es recomendable.
161
Capítulo 5
Resumen
Para el diseño, la implementación y la evaluación de estructuras de control se utilizan
una serie de herramientas numéricas que permiten la simplificación del modelo lineal,
determinación del emparejamiento entre las variables por controlar y manipular, la im-
plementación de los controladores y su sintonización para, por último, investigar la esta-
bilidad del lazo cerrado. En este capítulo se revisarán herramientas útiles aplicadas tanto
a casos de estudio modelo como al reactor Van de Vusse. Inicialmente, se presentan tres
estrategias para la identificación de modelos lineales a partir de datos generados por el
modelo no lineal, útil para realizar una reducción de orden del sistema que represente
mejor la respuesta en una región determinada. Porteriormente, se presentan dos meto-
dologías para el diseño de estructuras de control que perminten no solamente realizar
el emparejamiendo sino anticipar problemas potenciales de la implementación de múlti-
ples lazos de control. Enseguida, se muestran estrategias de implementación para control
on-off, Proporcional-Integral-Derivativo (PID) definido por el usuario y cómo usar el blo-
que PID predeterminado de Simulink. Esta sección es particularmente interesante pues,
se introduce cómo realizar pruebas de seguimiento de trayectorias o pruebas de rechazo
a perturbaciones ante la presencia de ruido. Así mismo, indica cómo evaluar diferen-
tes estrategias de sintonización de controladores tanto definidas por el usuario como la
herramienta de Matlab para llevarlo a cabo. Por último, se implementan estategias para
determinar la estabilidad del lazo cerrado de control tanto para control proporcional como
para controladores PID.
5.1. Introducción
Este capítulo está estructurado en tres partes principales. En la primera, se exploran dife-
rentes métodos para determinar la representación lineal de un sistema a partir de datos; esto
se conoce como ídentificación de sistemas (en inglés: system identification o data-driven mo-
dels). Identificar modelos a partir de series de datos es conveniente cuando no tenemos un
buen entendimiento del sistema para plantear un modelo basado en principios físicos o cuando
queremos obtener modelos de orden reducido a partir de simulaciones del modelo no lineal
(de manera alternativa a lo mostrado en la sección 4.5). En la segunda parte del capítulo,
163
Capítulo 5. Diseño y análisis de estructuras de control
Bien sea para la obtención de un modelo para un sistema complejo o para la reducción de
orden de un modelo no lineal con muchos estados, la identificación de sistemas es una herra-
mienta complementaria para el análisis de sistemas y el uso de modelos lineales para control,
pues se basa en el entendimiento de funciones de transferencia de bajo orden. Si el sistema
inherentemente es de orden superior, entonces se pueden utilizar técnicas de identificación de
sistemas para encontrar una función de transferencia de bajo orden que lo represente en la
vecindad de un punto de operación.
El Curve fitting es una aplicación gráfica de Matlab construida usando el Graphical User
Interface (GUI) para hacerla más amigable. Se puede acceder a ella directamente desde el
menú de aplicaciones de Matlab o escribiendo en el Command window cftool. La ventana
de inicio se muestra en la figura 5.1.
164
Capítulo 5. Diseño y análisis de estructuras de control
La ventana de trabajo se caracteriza por diferentes secciones descritas de arriba hacia abajo:
• Botones para realizar el ajuste o detenerlo: cuando se tiene activada la opción de au-
toajuste, al momento de realizar las modificaciones a las condiciones, automáticamente
la optimización se lleva a cabo. De otro modo, se usa el botón “fit” para correr la
optimización.
165
Capítulo 5. Diseño y análisis de estructuras de control
• Gráfica del resultado: muestra los datos usados para la regresión y los valores predichos
por el modelo.
• Tabla de ajuste: presenta el valor de los estadísticos que determinan la calidad del modelo
en una tabla.
La herramienta se usa de manera muy sencilla: como punto de partida, es necesario contar con
los datos para ajustar el modelo, los cuales pueden provenir de experimentos o de la simulación
de un sistema complejo no lineal (en Matlab o Simulink). Se requieren dos series de datos: la
de tiempo y la del estado de interés; entonces, el primer paso es realizar una simulación del
modelo no lineal en respuesta a una perturbación conocida. La naturaleza de la perturbación
es muy importante para la estimación, pues define la estructura del modelo y la utilidad de
los parámetros. Algunas señales útiles para la estimación son: paso, impulso, sinusoidales,
aleatorias estructuradas. Por su parte, para la identificación de procesos químicos, la entrada
paso es comúnmente usada.
Al tener los datos en el workspace, estos pueden ser reconocidos en la ventana de trabajo del
Curve fitting; simplemente se seleccionan de la lista de opciones que se despliega en cada
una de las series requeridas. Cuando los datos se han seleccionado, se grafican inmediatamente
y, por defecto, la herramienta realiza una regresión lineal; entonces aparece el modelo en la
parte de resultados con los parámetros de calidad del ajuste. En este punto se puede seleccionar
el tipo de función de transferencia a la que se le realizará el ajuste: primer orden, primer orden
más tiempo muerto, segundo orden, etc. Teniendo en cuenta que se manejan series de tiempo,
la ecuación definida por el usuario para el ajuste debe ser la solución analítica en el tiempo
de la función de interés. Nótese que esta solución analítica depende del tipo de perturbación
usada para obtener los datos. Cuando se actualiza la función, se obtienen los nuevos paráme-
tros y se ven los resultados gráficamente. El uso de la herramienta se muestra en el ejemplo 5.1.
Ejemplo 5.1.
Se desea determinar unas funciones de transferencia de bajo orden para el sistema Van de
Vusse alrededor del punto de operación D̄ = 4/7 min−1 y C̄Af = 10 mol l−1 (tabla 6.1). Se
166
Capítulo 5. Diseño y análisis de estructuras de control
tiene interés particular en las funciones de transferencia con respecto a la entrada de tasa de
dilución.
Solución
Para este problema, se va a utilizar la implementación del modelo no lineal en Simulink para
generar las series de tiempo requeridas para la estimación del modelo. Se perturba el sistema
con una entrada paso en la tasa de dilución de D = 4/7 min−1 a D = 5/7 min−1 en el tiempo
cero. El código que genera las series de tiempo es:
%--------------------------------------------------------------------------
% Programa para realizar identificación de modelos en el SIT_app
% SIT_app: System Identification Toolbox
% Modelo no lineal en Simulink
% Programa maestro en Matlab
clc
clear
close all
sim('VDV_NLsimulation')
167
Capítulo 5. Diseño y análisis de estructuras de control
sentarlos. Se utilizan los estadísticos de estimación para valorar la calidad del ajuste. Como
opción, se puede incluir los intervalos de confianza del predictor usando el menú de tools,
sección prediction bounds.
La ventana de configuración para una de las estimaciones se muestra en la figura 5.2, inclu-
yendo los resultados para la concentración de A:
Figura 5.2. Configuración de la herramienta Curve fitting para la determinación de una función
de transferencia de primer orden entre CA y D, R-square=1
168
Capítulo 5. Diseño y análisis de estructuras de control
Figura 5.3. Resultados del ajuste con la herramienta Curve fitting para la determinación de una
función de transferencia de primer orden más tiempo muerto entre CB y D, R-square=0.98
Reto 1: para una descripción más apropiada, encontrar la función que da el siguiente
resultado (figura 5.4):
Figura 5.4. Resultado de otra función usando Curve fitting para la determinación de una
función de transferencia entre CB y D
169
Capítulo 5. Diseño y análisis de estructuras de control
5.2.1.1. Comentarios
Es importante tener en cuenta que esta herramienta, al usar optimización determinística para
un problema de optimización no-convexo, no garantiza la convergencia al óptimo global ni
que realmente se tenga un buen ajuste. Dado el caso de que se presente algún problema con
el ajuste del modelo, es necesario verificar que las opciones de ajuste que vienen por defecto
sean las adecuadas para el sistema en cuestión. La opciones se pueden modificar haciendo clic
en el botón Fit options, donde aparecerá una ventana como se muestra en la figura 5.5. Algu-
nos parámetros importantes para modificar son: el algoritmo de optimización, los estimados
iniciales de los parámetros y los límites superior e inferior de los parámetros de ajuste. Así,
por ejemplo, cuando la dinámica del sistema es tan lenta que las constantes de tiempo son
exageradamente grandes, es probable que aparezca un error en el momento de hacer el ajuste.
Para ayudar al método, entonces se pueden poner valores apropiados en los límites superior
e inferior que estén de acuerdo con la dinámica del sistema. Esto último dependerá de cada
caso en especial y requiere conocimiento del sistema.
Figura 5.5. Ventana de configuración para realizar el ajuste en la herramienta Curve fitting
170
Capítulo 5. Diseño y análisis de estructuras de control
Por último, la herramienta Curve fitting es solo una interfaz gráfica de Matlab que usa los
métodos de optimización con resticciones. Es así que, para ajustes más complejos, por ejemplo
el ajuste de parámetros para el modelo cuando se usan ecuaciones diferenciales, se recomienda
desarrollar código directamente tanto para el problema de optimización como para determinar
los índices de desempeño del modelo. De esta manera, se tiene acceso por ejemplo a optimiza-
dores estocásticos, al uso de índices estadísticos más robustos para valorar el poder predictivo
del modelo o a la realización del proceso de validación propiamente del modelo y no simple-
mente estimación de parámetros. Este tipo de problemas está más allá del alcance de este libro.
dx
= Ax(t) + Bu(t) + Ke(t) (5.2)
dt
y(t) = Cx(t) + Du(t) + e(t)
Como se observa, esta realización genérica desacopla las entradas de las perturbaciones. Den-
tro de las opciones de las funciones, se puede fijar que ciertas matrices no sean estimadas para
darle así una configuración especial al modelo. Las funciones se utilizan de manera análoga.
Como ejemplo se muestra la siguiente sintaxis para usar la función n4sid:
%--------------------------------------------------------------------------
sys = n4sid(data,n,opt,Name,Value)
%--------------------------------------------------------------------------
171
Capítulo 5. Diseño y análisis de estructuras de control
Donde:
Una vez se estima el modelo, es necesario confirmar la representación obtenida. Esto se realiza
en el ejemplo 5.2 para un sistema multivariable.
Ejemplo 5.2.
Se desea determinar una representación en espacio de estados para el reactor Van de Vusse al-
rededor del punto de operación D̄ = 4/7 min−1 y C̄Af = 10 mol l−1 (tabla 6.1). Para efectuar
la estimación, se perturba el sistema con una entrada paso positiva en cada entrada alrededor
del punto de operación con amplitud de D = 1/7 min−1 y CAf = 0,5 mol l−1 .
Solución
En la segunda parte, se debe perturbar el modelo no lineal para generar las series de tiempo
de respuesta. Para esto se utiliza la implementación de modelo no lineal en Simulink con
entradas paso. En este caso, se propone utilizar 10 min de simulación, con la perturbación en
la tasa de dilución en el tiempo cero y en la concentración de entrada transcurridos 5 min.
Esta configuración se implementa directamente en Simulink en los bloques de las entradas,
donde los valores de las entradas antes y después de la perturbación se especifican desde el
programa principal en Matlab. En esta simulación en particular, se cambia el método numé-
rico de solución de ecuaciones diferenciales en Simulink a uno discreto para generar series de
tiempo con tamaño de paso constante, el cual se especifica también desde Matlab. Se utilizan
nuevamente los bloques To workspace para que las entradas y salidas estén disponibles en el
workspace para el programa principal en Matlab. Se ejecuta entonces el programa de Simulink
llamado VDV_NLsimulation.
En la tercera parte parte del programa se realiza la estimación del modelo continuo en espacio
de estados con ayuda del comando n4sid. Hay que recordar que para realizar la estimación
172
Capítulo 5. Diseño y análisis de estructuras de control
Como cuarto paso, se simula la respuesta de sistema lineal analítico y el identificado para
valorar inicialmente la representación de los modelos con respecto a las series de tiempo del
modelo no lineal. Los resultados finalmente son graficados.
%--------------------------------------------------------------------------
% Programa para realizar identificación de modelos en espacio de estados
% Modelo no lineal en Simulink
% Programa maestro en Matlab
clc
clear
close all
% Parámetros cinéticos
173
Capítulo 5. Diseño y análisis de estructuras de control
% Estado estable
a11=-Dss-k1-2*k3*Cao;
a12=0;
a21=k1;
a22=-Dss-k2;
b11=Cafss-Cao;
b12=Dss;
b21=-Cbo;
b22=0;
syst=ss(At,Bt,Ct,Dt);
% Aplicación de la perturbación
sim('VDV_NLsimulation')
174
Capítulo 5. Diseño y análisis de estructuras de control
sys = n4sid(Data,2,'Ts',0,'DisturbanceModel','none');
sys =
Continuous-time identified state-space model:
dx/dt = A x(t) + B u(t) + K e(t)
y(t) = C x(t) + D u(t) + e(t)
A =
x1 x2
x1 -7.406 -29.36
x2 0.8686 2.582
B =
D Caf
x1 149.4 10.35
x2 -17.05 -0.8423
175
Capítulo 5. Diseño y análisis de estructuras de control
C =
x1 x2
Ca 0.1266 0.7091
Cb 0.01879 0.2302
D =
D Caf
Ca 0 0
Cb 0 0
K =
Ca Cb
x1 0 0
x2 0 0
Parameterization:
FREE form (all coefficients in A, B, C free).
Feedthrough: none
Disturbance component: none
Number of free coefficients: 12
Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Estimated using N4SID on time domain data "Data".
Fit to estimation data: [98.08;98.45]%
FPE: 2.233e-13, MSE: 3.385e-06
Como se puede ver, no solamente se arrojan los resultados de la estimación sino también esta-
dísticas de incertidumbre y calidad del ajuste a través de las funciones idssdata, getpvec,
getcov y los valores de fit to estimation data, error final de predicción de Akaike (Akaike’s
Final Prediction Error-FPE) y promedio del cuadrado de los errores (Mean-Squared Error-
MSE). Los resultados de la simulación son (figura 5.7):
176
Capítulo 5. Diseño y análisis de estructuras de control
Figura 5.7. Comparación entre los resultados de los modelos del sistema (modelo no lineal), lineal
(modelo lineal obtenido analíticamente) y el identificado (modelo lineal identificado) usando las mismas
señales de entrada
Es evidente que el modelo lineal analítico se desvía del modelo no lineal debido a que este
es válido en una vecindad muy cercana al punto de operación. Por otro lado, este error no
aparece en el modelo estimado a causa de la naturaleza del ajuste realizado.
177
Capítulo 5. Diseño y análisis de estructuras de control
Ejemplo 5.3.
Se desea determinar unas funciones de transferencia de bajo orden para el sistema Van de
Vusse alrededor del punto de operación D̄ = 4/7 min−1 y C̄Af = 10 mol l−1 (tabla 6.1). Se
178
Capítulo 5. Diseño y análisis de estructuras de control
tiene interés en las funciones de transferencia con respecto a la entrada de tasa de dilución.
Solución
Para este problema se va a utilizar la implementación del modelo no lineal que se tiene en
Simulink para generar las series de tiempo requeridas para la estimación del modelo. Con
este objetivo, se perturba el sistema con una entrada paso en la tasa de dilución de D = 4/7
min−1 a D = 5/7 min−1 en el tiempo cero. El código que genera las series de tiempo es el
mismo del ejemplo 5.1. Es preciso tener presente que se usa un método numérico discreto para
la solución de las ecuaciones diferenciales con un tiempo de muestreo Ts=0.01; además, que
a la herramienta de identificación de sistemas se le alimentan datos en variables de desviación.
Figura 5.9. Configuración de la herramienta System Identification app con las series de tiem-
po para realizar la estimación
Para ver las series de tiempo cargadas, se señala la opción time plot que se encuentra debajo
de los bloques donde se importaron los datos. Para que salga cada serie de tiempo en una
figura independiente, es preciso asegurarse de que solo una señal esté señalada en los bloques
de importación de datos. Las figuras generadas por la herramienta no tienen un buen formato,
pero pueden ser exportadas a través de la opción File - copy figure. Teniendo la figura como
179
Capítulo 5. Diseño y análisis de estructuras de control
Figura 5.10. Series de tiempo cargadas como entradas para la estimación de funciones de transferencia
en la herramienta System Identification app
Una vez se tengan los datos cargados para la estimación, se despliega la ventana Estimate,
donde se selecciona la opción Transfer function model; después de esto, se evidencia la siguiente
ventana emergente (figura 5.11):
180
Capítulo 5. Diseño y análisis de estructuras de control
En esta ventana se indica el número de polos y ceros para identificar, además si se quiere
identificar un modelo continuo o discreto; en este caso, se indicó 2 polos y 1 cero, y que el
modelo deseado es continuo. Adicionalmente, aparecen opciones para tiempo muerto (i.e.,
I/O delay) y para algunas condiciones para la estimación. Para este ejemplo, se dejaron las
opciones por defecto, pero se forzó que las condiciones iniciales fueran cero. Una vez con-
figurada la estimación se procede con el botón Estimate y es entonces cuando aparece una
ventana llamada Plant Identification Process. Esta ventana muestra información del modelo
para identificar y la serie de tiempo, además de datos acerca de la estimación como el método
y las iteraciones realizadas para estimar los parámetros del modelo. Finalmente, muestra los
resultados de la estimación indicando el grado de ajuste con el coeficiente de determinación,
el FPE y MSE. La ventana resultante para la estimación de la función de transferencia entre
CA y D se muestra en la figura 5.12:
181
Capítulo 5. Diseño y análisis de estructuras de control
La función de transferencia identificada puede evidenciarse dando doble clic sobre el resultado
llamado por defecto tf1. La siguiente ventana se despliega (figura 5.13):
182
Capítulo 5. Diseño y análisis de estructuras de control
183
Capítulo 5. Diseño y análisis de estructuras de control
Figura 5.14. Resultados de la predicción del modelo entre CA y D, con ajuste de r2 =100 %
Figura 5.15. Resultados de la predicción del modelo entre CB y D con un ajuste de r2 =99.97 %
De lo anterior se evidencia el alto grado de ajuste teniendo en cuenta que el modelo está
prácticamente superpuesto con los datos usados para la estimación.
184
Capítulo 5. Diseño y análisis de estructuras de control
T
−1
λij = (Kij ) K (5.4)
ij
T
−1
Λ=K⊗ K (5.5)
185
Capítulo 5. Diseño y análisis de estructuras de control
Donde:
Como tal, Matlab no ofrece una función que permita hacer el cálculo, razón por la que deben
realizarse las operaciones de manera individual. En el caso de no tener una matriz cuadrada
de ganancias, indicando que no se tiene igual número de variables controladas y para mani-
pular, en principio no se puede obtener la matriz inversa de ganancias; para este caso se usa
el comando de la matriz pseudoinversa de Moore-Penrose con el comando pinv.
Ejemplo 5.4.
Solución
Como punto de partida para determinar la matriz RGA, se toma el modelo lineal del sistema
en cuestión. En principio se puede utilizar cualquiera de las implementaciones mostradas en la
sección 3.3, un modelo simplificado usando reducción de modelos o un modelo identificado; en
este caso se utiliza la linealización analítica del sistema. Para realizar el emparejamiento hay
que definir cuántos de los estados se desea controlar. Teniendo en cuenta que se tienen solo
dos estados, se supone que se busca controlarlos a ambos usando los dos grados de libertad
disponibles. La viabilidad de estos objetivos de control fue corroborada al determinar la ma-
triz de controlabilidad del sistema. Posteriormente, se debe determinar la matriz de ganancias
para el sistema que será 2 x 2. Se proponen dos estrategias para esto:
186
Capítulo 5. Diseño y análisis de estructuras de control
Una vez se tiene la matriz de ganancias, se aplica la ecuación para determinar la matriz
RGA. Finalmente, ambos resultados son comparados desde el punto de vista de las matrices
de ganancias y RGA. Para cuantificar el error, se realiza un rearreglo de la matriz como vector.
%--------------------------------------------------------------------------
% Programa para realizar el diseño de estructura de control
% Método - arreglo de ganancias relativas
clc
clear
close all
% Parámetros cinéticos
% Estado estable
a11=-Dss-k1-2*k3*Cao;
a12=0;
a21=k1;
a22=-Dss-k2;
b11=Cafss-Cao;
b12=Dss;
b21=-Cbo;
b22=0;
syst=ss(A,B,C,D);
for i=1:length(A)
for j=1:length(B)
K1(i,j)=-C(j,:)*inv(A)*B(:,i)+D(j,i);
end
187
Capítulo 5. Diseño y análisis de estructuras de control
end
RGA1=K1.*(inv(K1))'
for j=1:length(B)
[num,den]=ss2tf(A,B,C,D,j);
for i=1:length(A)
K2n(j,i)=poly2sym(num(i,:))/poly2sym(den);
end
end
x = 0;
K2 = eval(K2n);
RGA2=K2.*(inv(K2))';
%% Resultados
figure
subplot(121)
bar(categorical({'K11','K12','K21','K22'}),...
reshape((K1-K2)./K1,[1,prod(size(K1))]));
xlabel('Ganancias','FontWeight','bold','FontSize',14)
ylabel('Error ( %)','FontWeight','bold','FontSize',14)
set(gca,'FontSize',12,'FontWeight','bold')
subplot(122)
bar(categorical({'RGA11','RGA12','RGA21','RGA22'}),...
reshape((RGA1-RGA2)./RGA1,[1,prod(size(RGA1))]));
xlabel('Arreglo de ganancias relativas','FontWeight','bold','FontSize',14)
ylabel('Error ( %)','FontWeight','bold','FontSize',14)
set(gca,'FontSize',12,'FontWeight','bold')
K1 =
2.9109 0.5848
0.2376 0.0885
K2 =
2.9109 0.5848
0.2376 0.0885
RGA1 =
2.1717 -1.1717
-1.1717 2.1717
RGA2 =
188
Capítulo 5. Diseño y análisis de estructuras de control
2.1717 -1.1717
-1.1717 2.1717
Figura 5.16. Errores relativos en la determinación de las ganancias de lazo abierto y ganancias
relativas para el sistema Van de Vusse
Se observa en las figuras que las diferencias entre ambas aproximaciones para la determina-
ción de las ganancias y RGA es prácticamente un error numérico del orden de magnitud de
la exactitud relativa del punto flotante.
Este método es muy interesante, pues la descomposición en valores singulares (SVD, por sus
siglas en inglés), trae información importante de las ganancias y la direccionalidad multiva-
riable, que ayuda a entender sistemas en los que todas las entradas afectan todas las salidas
en diferente medida (i.e., sistemas interaccionantes).
189
Capítulo 5. Diseño y análisis de estructuras de control
ȳ = K ū (5.7)
|K − λI| = 0 (5.8)
Si uno de los eigenvalores es cero, entonces K es una matriz singular y el sistema al que repre-
senta será difícil de controlar. Ahora, si por lo menos uno de los eigenvalores es muy pequeño
en comparación con los otros, entonces se requerirán grandes cambios en una o más variables
manipuladas para poder controlar el sistema. Si bien esta interpretación de los eigenvalores
es útil, se considera que no es muy robusta para caracterizar sistemas interaccionantes; es allí
donde se propone utilizar mejor la descomposición en valores singulares. Algunas caracterís-
ticas de la descomposición en valores singulares son:
• Los primeros r valores singulares son números positivos, donde r es el rango de la matriz
K T · K.
190
Capítulo 5. Diseño y análisis de estructuras de control
• Por lo general, los valores singulares diferentes de cero se ordenan de mayor a menor
(i.e., σ1 : mayor, σr : menor valor singular diferente de cero).
Para utilizar la información contenida en los valores singulares, se recurre al concepto del
número de condición (CN, por sus siglas en inglés), el cual es la razón entre el mayor y menor
valor singular de la matriz K. Teniendo en cuenta el orden en el que se arrojan los valores
singulares, el número de condición se define como
σ1
CN = (5.9)
σr
Esto es análogo al caso de los eigenvalores: cuando se tiene que K es singular, el CN será
infinito. Ahora, si el sistema es mal condicionado (ill-conditioned), el número de condición
será muy grande (e.g., CN >10), lo que indica que pequeños cambios en el proceso hacen muy
difícil de controlar la planta.
Al llevar a cabo una comparación con la matriz RGA, se encuentra que esta última es esca-
lada (normalizada) y, por tanto, los valores contenidos en la matriz no dependen de la escala
de entradas y salidas. Por el contrario, la descomposición en valores singulares sí depende de
la escala.
Ejemplo 5.5.
Se requiere determinar en número de condición para el sistema Van de Vusse alrededor del
punto de operación D̄ = 4/7 min−1 y C̄Af = 10 mol l−1 (tabla 6.1), donde se desea contro-
lar la concentración de A y B en el reactor usando las dos entradas disponibles (sistema 2×2).
Solución
191
Capítulo 5. Diseño y análisis de estructuras de control
(i.e., Kij = −C jrow ∗ A−1 ∗ B icol + Dijrow ). Es entonces cuando se aplica la descomposición en
valores singulares y se determina el número de condición.
%--------------------------------------------------------------------------
% Programa para realizar el diseño de estructura de control
% Método - Descomposición en valores singulares
clc
clear
close all
% Parámetros cinéticos
% Estado estable
a11=-Dss-k1-2*k3*Cao;
a12=0;
a21=k1;
a22=-Dss-k2;
b11=Cafss-Cao;
b12=Dss;
b21=-Cbo;
b22=0;
syst=ss(A,B,C,D);
for i=1:length(A)
for j=1:length(B)
K(i,j)=-C(j,:)*inv(A)*B(:,i)+D(j,i);
end
end
192
Capítulo 5. Diseño y análisis de estructuras de control
disp(' ')
disp('La matriz de valores singulares es:')
disp(' ')
disp(S)
disp(' ')
disp('El número de condición es:')
disp(' ')
disp(CN)
2.9796 0
0 0.0398
74.8596
Los resultados obtenidos en el número de condición son congruentes con el ejemplo 5.4, donde
la matriz RGA anticipaba las interacciones en el sistema y posibles problemas para controlar
los dos estados simultáneamente.
Para mostrar cómo descomponer un sistema más grande que el reactor Van de Vusse y luego
evaluar qué estructuras de control son pertinentes, se propone el ejemplo 5.6.
Ejemplo 5.6.
Para el sistema representado por la matriz de ganancias que se muestra a continuación, de-
terminar el número de condición para el sistema completo y para todas las descomposiciones
posibles de sistemas 2 × 2.
193
Capítulo 5. Diseño y análisis de estructuras de control
0,48 0,90 −0,006
K = 0,52 0,95 0,008
0,90 −0,95 0,020
Solución
%--------------------------------------------------------------------------
% Programa para realizar el diseño de estructura de control
% Método - Descomposición en valores singulares
clc
clear
close all
disp(' ')
disp('La matriz de RGA es:')
disp(' ')
disp(RGA)
disp(' ')
disp('La matriz de valores singulares es:')
disp(' ')
194
Capítulo 5. Diseño y análisis de estructuras de control
disp(diag(S))
disp(' ')
disp('El número de condición es:')
disp(' ')
disp(CN)
for i=1:length(K)
for j=1:length(K)
Kp=K;
Kp(i,:) = [];
Kp(:,j) = [];
[U,S,V] = svd(Kp); % S contiene los valores singulares
RGAnew = Kp.*(inv(Kp))'; % RGA sistema reducido
CN(i,j) = S(1,1)/S(end); % Número de condición
NewIndIp = index;
NewIndJp = index;
NewIndIp(i)= [];
NewIndJp(j)= [];
Strc{i,j} = ['y' num2str(NewIndIp(1)) '-' 'y' num2str(NewIndIp(2));
'u' num2str(NewIndJp(1)) '-' 'u' num2str(NewIndJp(2))];
RGA_M{i,j}= RGAnew;
end
end
%% Resultados
disp(' ')
disp('V. Controladas Var. Manipuladas CN ')
for i=1:length(K)
for j=1:length(K)
disp([Strc{i,j}(1,:) ' ' Strc{i,j}(2,:) ' '...
' ' num2str(CN(i,j))])
end
end
%--------------------------------------------------------------------------
195
Capítulo 5. Diseño y análisis de estructuras de control
1.6183
1.1434
0.0097
166.5212
Según la matriz de RGA del sistema completo, el emparejamiento recomendado para este es
muy complejo y puede proponerse: y1 − u1 , y2 − u3 y y3 − u2. Sin embargo, en los valores
obtenidos del RGA se nota que los valores de λ < 0,5 implican problemas. Ahora, el valor del
número de condición para el sistema completo (CN = 166.5212) confirma que el control del
sistema completo es muy difícil. Por tanto, se realizan las reducciones del sistema y se nota
que el número de condición para los subsistemas y2 − y3 con u1 − u2 , y y1 − y3 con u1 − u2
tienen valores relativamente bajos. Para estos sistemas se extraen manualmente los RGA, con
el fin de validar los emparejamientos:
RGA_M{1,3}=
0.3662 0.6338
0.6338 0.3662
RGA_M{2,3}=
0.3602 0.6398
0.6398 0.3602
De esta última validación, se concluye que las estructuras más prometedoras para ser evaluadas
serían:
196
Capítulo 5. Diseño y análisis de estructuras de control
• Estructura 1: y2 − u2 y y3 − u1
• Estructura 2: y1 − u2 y y3 − u1
Nota: es preciso tener cuidado con los índices de las estructuras diseñadas.
5.3.2.1. Comentarios
La implementación de este tipo de controladores es muy sencilla debido a que funciona como
un condicional o interruptor, en la que la salida del controlador se conmuta entre dos valores
así:
(
pmax si e ≥ 0
p(t) = (5.10)
pmin si e < 0
Donde los valores de pmax y pmin son los valores de la acción de control cuando la variable
para controlar está por encima o por debajo del valor del set-point (error = yset − y). Nótese
que los valores de p(t) máximo y mínimo pueden ser 0 % y 100 % cuando se trata de válvulas
de control, pero pueden ser otros valores nominales cuando se refiere a otras variables mani-
puladas.
Ahora, como se observa, este controlador on-off puro comienza a operar tan pronto se tiene
un error positivo o negativo en relación con el valor deseado de la variable controlada. Esto
es conveniente para sistemas lentos y que no poseen ruido. Teniendo en cuenta que estos
dos últimos factores pueden no cumplirse, la acción de los controladores on-off se relaja al
introducir una banda muerta o brecha diferencial que reduce la sensibilidad del controlador
al error:
197
Capítulo 5. Diseño y análisis de estructuras de control
(
pmax si e ≥ δ
p(t) = (5.11)
pmin si e < −δ
Es necesario indicar cuándo debe encenderse y apagarse la acción del bloque dependiendo de
la señal alimentada (i.e., Switch on point y Switch off point); así como qué valor debe tomar
la acción de control p(t) en estas dos situaciones, es decir, pmax y pmin . En el ejemplo 5.7 se
muestra cómo implementar el control on-off y realizar diferentes tipos de pruebas de rechazo
a perturbaciones y cambio de set-point incluyendo ruido.
198
Capítulo 5. Diseño y análisis de estructuras de control
Ejemplo 5.7.
Determinar la respuesta dinámica del sistema no lineal del reactor Van de Vusse utilizando
un controlador on-off de brecha diferencial, evaluando la influencia de la configuración del
controlador (i.e., tamaño de la brecha diferencial), diferentes perturbaciones (i.e., cambio de
set-point y perturbación en entrada) y la presencia de ruido. Como estructura de control,
se establece controlar la concentración de salida del componente B, utilizando como variable
manipulada la concentración de entrada del reactivo A (CAf ). Las condiciones iniciales del
reactor corresponden al estado estable óptimo (i.e., D̄ = 4/7 (min−1 ), C̄Af = 10 (mol/l),
C̄A = 3 (mol/l) y C̄B = 1,117 (mol/l)).
Solución
Para el programa principal de Matlab es necesario definir las condiciones iniciales de las varia-
bles de entrada del sistema, la tasa de dilución y la concentración de entrada de A, estas son
requeridas por el programa de Simulink. Además se define el tiempo de simulación, el cual fue
definido con simulaciones preliminares. Seguidamente, se establecen parámetros requeridos por
el controlador como el tamaño del la brecha diferencial (2δ), el set-point de la concentración
de B y unas variables binarias para activar las pruebas de seguimiento de set-point o activa-
ción de ruido en las simulaciones. En esta sección también hay que difinir las perturbaciones
por aplicar en la tasa de dilución y cambio de set-point de la concentración de B, y en qué
momento ocurren. Estas últimas variables se requieren para bloques de Simulink. Finalmente,
el código de Matlab ejecuta el programa de Simulink llamado VDV_NLsim_on_off_control
y grafica los resultados.
%--------------------------------------------------------------------------
% Programa para realizar control on-off del reactor VDV
% Modelo no lineal en Simulink
% Programa maestro en Matlab
clc
clear
close all
199
Capítulo 5. Diseño y análisis de estructuras de control
sim('VDV_NLsim_on_off_control')
%% Resultados
200
Capítulo 5. Diseño y análisis de estructuras de control
grid on
subplot(224)
plot(tiempo,Cb,'b','LineWidth',2);
hold on
plot(tiempo,Cbset,'r','LineWidth',2)
plot(tiempo,Cbset+H_Cb/2,'--r','LineWidth',2)
plot(tiempo,Cbset-H_Cb/2,'--r','LineWidth',2)
legend('Cb','Cb_{set}')
xlabel('Tiempo [min]','FontWeight','bold','FontSize',14)
ylabel('C_B [mol/l]','FontWeight','bold','FontSize',14)
set(gca,'FontSize',12,'FontWeight','bold')
grid on
%--------------------------------------------------------------------------
En Simulink se prepara el programa mostrado en la figura 5.18, que contiene el modelo del
reactor Van de Vusse implementado como una S-function, el controlador y las perturbacio-
nes.
Figura 5.18. Implentación en Simulink del control on-off de la concentración de B utilizando CAf
como variable manipulada
201
Capítulo 5. Diseño y análisis de estructuras de control
Figura 5.19. Configuración definida para la prueba de set-point tracking utilizando el bloque de
Signal builder
Para el control se usa el bloque Relay dentro de una estructura de control retroalimentado
y la configuración mostrada en la figura 5.20. Los parámetros requeridos se obtienen desde el
programa maestro en Matlab; además, es preciso tener cuidado con los signos.
202
Capítulo 5. Diseño y análisis de estructuras de control
Finalmente, el programa incluye un bloque Switch para que la variable medida pueda uti-
lizar ruido o no, dependiendo de la variable Noise_on, que se define en el programa de
Matlab. Si se desea realizar la simulación con ruido, se debe establecer esta señal con el blo-
que Random number, el cual crea una señal aleatoria definida por su media y su varianza. En
este caso se usa un valor medio de cero y una varianza definida por la variable Noise_var.
Utilizando el programa principal, se tiene una herramienta que permite valorar diferentes
configuraciones del sistema; inicialmente, el desempeño del controlador frente a un cambio
único de set-point y el rechazo de una perturbación de la tasa de dilución como se muestra
en la figura 5.21.
Figura 5.21. Respuesta del sistema de control on-off ante cambios simples de set-point y tasa de
dilución
El mismo caso anterior puede ser evaluado en presencia de ruido, como se muestra en la figura
5.22, donde se observa el gran efecto que tiene en el desempeño de la variable manipuada y
controlada.
203
Capítulo 5. Diseño y análisis de estructuras de control
Figura 5.22. Respuesta del sistema de control on-off ante cambios simples de set-point y tasa de
dilución cuando se aplica ruido
Otra opción es realizar la prueba de set-point tracking más el rechazo de perturbaciones como
se muestra en la figura 5.23.
Figura 5.23. Respuesta del sistema de control on-off para la prueba de set-point tracking más rechazo
a perturbación en la tasa de dilución
204
Capítulo 5. Diseño y análisis de estructuras de control
Reto 1: realizar un programa aún más generalizado para ejecutar más pruebas y
generar una sola figura que evidencie el impacto de cada variable modificada.
5.4.1. Comentarios
Si bien el control on-off es sencillo de implementar y evaluar, es evidente que resulta en va-
riaciones cíclicas de la variable controlada debido a la acción tipo onda cuadrada que tiene la
variable manipulada. Posiblemente, las oscilaciones en la variable controlada no son necesa-
riamente un problema, aunque sí puede existir un estrés excesivo sobre el actuador debido a
los cambios abruptos a los que es expuesto cuando se tienen dinámicas rápidas del sistema o
brechas diferenciales muy pequeñas. Esto debe ser evaluado en cada caso para validar que no
se trate de un problema en la implementación real.
Z t
1 de(t)
∗ ∗
p(t) = p̄ + Kc e(t) + e(t )dt + τd (5.12)
τi 0 dt
Donde:
205
Capítulo 5. Diseño y análisis de estructuras de control
P (s) 1
= Kc 1 + + τd s (5.13)
E(s) τi s
Donde:
Figura 5.24. Implementación de un control PID en Simulink usando la estructura en paralelo. Nota:
los bloques utilizados son Gain, Integrator y derivative
206
Capítulo 5. Diseño y análisis de estructuras de control
P (s) 1 τd s
= Kc 1 + + (5.14)
E(s) τi s ατd s + 1
Donde α es un valor entre 0.05 y 0.2, con 0.1 como valor típico. Esta implementación resuelve
el problema de la condición de realización de la parte derivativa del controlador, así como
también reduce la sensibilidad ante el ruido y ayuda con el fenómeno del derivative kick. La
implementación en Simulink del PID con filtro derivativo se muestra en la figura 5.25. No
obstante, se recomienta definitivamente eliminar el fenómeno del derivative kick, y aplicar la
acción derivativa sobre la variable medida y no sobre el error (Seborg et al., 2011).
Figura 5.25. Implementación de un control PID en Simulink con filtro derivativo usando la estructura
en paralelo. Nota: los bloques utilizados son Gain, Integrator, y derivative
Como adición a la implementación en Simulink, se puede utilizar una acción para evitar un
fenómeno llamado reset windup, que afecta la parte integral del control PID cuando el siste-
ma tiene errores persistentes. Por ejemplo: cuando se inicia un proceso (start-up), la variable
controlada está alejada del valor del set-point. Esto implica que hay un error grande que se
mantiene por un intervalo de tiempo y, como consecuencia, la integral del error puede adqui-
rir valores muy grandes provocando que el actuador se sature. La saturación del actuador se
manifiesta, por ejemplo, cuando una válvula se abre al 100 % o se cierra totalmente. Incluso
207
Capítulo 5. Diseño y análisis de estructuras de control
cuando la variable controlada llegue al set-point, la acción integral sigue saturando el actuador
lo que genera un gran Overshoot difícil de eliminar. Para removerlo naturalmente, la acción
integral debería tener un error en el sentido contrario por un largo periodo de tiempo (Seborg
et al., 2011).
Para evitar el fenómeno de reset windup, se puede incorporar a la implementación del PID
una acción que borra la memoria de la acción integral. Por este efecto, esa acción se denomina
antireset windup. La implementación del controlador PID con antireset windup se muestra en
el ejemplo 5.8.
Ejemplo 5.8.
Determinar la respuesta dinámica del sistema no lineal del reactor Van de Vusse utilizando
un controlador PID con acción antireset windup, para evaluar la influencia de la sintoniza-
ción del controlador, la reacción ante diferentes perturbaciones (i.e., cambio de set-point y
perturbación en entrada) y la presencia de ruido. Como estructura de control, se establece
controlar la concentración de salida del componente B, utilizando como variable manipulada
la concentración de entrada del reactivo A (CAf ). La condiciones iniciales del reactor corres-
ponden al estado estable óptimo (i.e., D̄ = 4/7 (min−1 ), C̄Af = 10 (mol/l), C̄A = 3 (mol/l) y
C̄B = 1,117 (mol/l)).
Solución
De manera análoga al caso del control on-off, se propone utilizar la estrategia de programación
con Matlab como programa maestro y usar Simulink para la implementación del reactor y el
controlador. Se recicla la estrategia de conmutación para usar las diferentes perturbaciones y
adicionar ruido a la simulación. Como complemento al caso anterior, en la parte inicial del
programa de Matlab, donde se hacen las confirguraciones, se define una nueva variable llama-
da Useantiwindup, que se utilizará en la implementación del controlador PID para activar
la acción de antireset windup.
Adicionalmente, se diseñó una función llamada PID_tuner (no mostrada, para dejarla como
reto), la cual tiene el método de sintonización IMC (Internal Model Control) para calcular
las constantes del controlador (i.e., Kc, τi y τd ), que corresponden a la salida de la función
(Seborg et al., 2011). Teniendo en cuenta que la función de sintonización utiliza el modelo
linealizado representado en funciones de transferencia, esta función requiere como entradas
el índice de la variable que se va a controlar, el índice de la variable manipulada y qué tan
rápido se desea la respuesta en lazo cerrado en comparación con el lazo abierto. Para este
caso la variable para controlar es CB con índice 2 (segunda ecuación diferencial); la variable
manipulada es CAf con índice 2, pues es la segunda entrada y se desea que la respuesta en
lazo cerrado sea 5 veces más rápida que en lazo abierto, así que τc = τ /5. En la función,
208
Capítulo 5. Diseño y análisis de estructuras de control
la sintonización del parámetro del filtro diferencial se hace manualmente; además, la función
permite introducir los parámetros del controlador sin usar el método IMC. Posteriomente,
se definen las perturbaciones, el código de Matlab ejecuta el programa de Simulink llamado
VDV_NLsim_PID_control_pedal y grafica los resultados.
%--------------------------------------------------------------------------
% Programa para realizar control PID del reactor VDV
% PID-antiresetwindup usando bloque contrucción propia
% Modelo no lineal en Simulink
% Programa maestro en Matlab
clc
clear
close all
% Definiciones de controlador
% Sintonización
sim('VDV_NLsim_PID_control_pedal')
209
Capítulo 5. Diseño y análisis de estructuras de control
%% Resultados
%--------------------------------------------------------------------------
Ahora, en Simulink se prepara el programa mostrado en la figura 5.26, que contiene el modelo
del reactor Van de Vusse implementado como una S-function, el controlador PID en un
subsistema y las perturbaciones usadas anteriormente en el controlador on-off.
210
Capítulo 5. Diseño y análisis de estructuras de control
Figura 5.26. Implentación en Simulink del control PID de la concentración de B utilizando CAf como
variable manipulada
Se observa que la implementación en Simulink es casi igual que en el caso anterior, al sustituir
el bloque de Relay por un bloque de subsistema llamado PID antireset windup, que se
construyó con la implementación del controlador. Nótese que a la salida del controlador se
adiciona el p̄, llamado bias; este valor corresponde al valor de la variable manipulada en el
estado estable inicial. Lo anterior es necesarion pues el controlador está en variables de des-
viación y hay que retornarlo a variables nominales.
Dentro del bloque PID antireset windup se incluye el controlador PID mostrado en la
figura 5.25, al cual se adiciona una sección para borrar la memoria de la acción integral para
el caso en que el actuador se sature, modificación denominada acción antireset windup. Esta
acción es necesaria para controladores usados en arranque de equipos o situaciones en las que
se tiene grandes cambios de set-point, donde la memoria de la acción integral hace que se
sature el actuador provocando overshoots significativos (Seborg et al., 2011).
211
Capítulo 5. Diseño y análisis de estructuras de control
Figura 5.27. Implementación de un controlador PID con acción antireset windup para el control de
la concentración CB usando CAf
En la implementación del PID, la acción PID antireset windup se logra evaluando la señal
de salida del controlador antes y después de un saturador. Si ambas señales son iguales, su
resta hace que la acción PID antireset windup sea cero y el controlador PID con filtro
diferencial no se afecta. Sin embrago, si la señal del controlador es mayor al máximo o menor
al mínimo establecido en el saturador, restarlas genera un valor diferente de cero. Este valor es
multiplicado por el mismo valor de la acción integral (aunque podría ser diferente ganancia)
y se adiciona a las acciones proporcional, integral y derivativa. Al realizar esto, se genera una
acción contraria a la acumulación de error por la parte integral, que al final viene a disminuir
la acumulación de error en la memoria del controlador. Como se mencionó, se adicionó un
Switch con la variable Useantiwindup para activar el efecto PID antireset windup.
212
Capítulo 5. Diseño y análisis de estructuras de control
Figura 5.28. Respuesta del sistema de control PID con filtro derivativo y antireset windup ante
cambios simples de set-point y tasa de dilución
El caso anterior también se puede evaluar en presencia de ruido, como se muestra en la figura
5.29.
Figura 5.29. Respuesta del sistema de control PID con filtro derivativo y antireset windup ante
cambios simples de set-point y tasa de dilución cuando se aplica ruido
213
Capítulo 5. Diseño y análisis de estructuras de control
Figura 5.30. Respuesta del sistema de control PID con filtro derivativo y antireset windup para
la prueba de set-point tracking más rechazo a perturbación en la tasa de dilución
Reto 1: construir una función generalizada que permita realizar la sintonización del
controlador usando diferentes métodos y así valorar su efecto en la respuesta.
Reto 3: el controlador PID puede ser implementado como una S-function. Intentar
implementarlo para valorar su desempeño. Esto es práctico cuando de tienen varios
lazos de control en el sistema.
214
Capítulo 5. Diseño y análisis de estructuras de control
Simulink permite implementar de varias formas el PID en paralelo, como se muestra en las
siguientes ecuaciones:
!
P (s) 1 N
=P +I +D (5.15)
E(s) s 1 + N 1s
!!
P (s) 1 N
=P 1+I +D 1 (5.16)
E(s) s 1+ s
Los aspectos por configurar en la ventana principal para este bloque son (MatWorks, 2023):
215
Capítulo 5. Diseño y análisis de estructuras de control
por los bloques del modelo. En control de procesos, usando un modelo en ecuaciones
diferenciales, el tiempo de muestreo será definido por el método numérico que resuelve las
ecuaciones.
e) Ventana Main: esta ventana se encuentra debajo de la fórmula del controlador. Aquí se
tienen los parámetros del controlador y la opción de usar el Autotuning para encontrarlos.
216
Capítulo 5. Diseño y análisis de estructuras de control
• Integrator and filter initial conditions: para esta parte se selecciona si se desea espe-
cificar las condiciones iniciales del integrador y el filtro diferencial de manera interna
o externa con bloques específicos. Las condiciones iniciales para cada uno de estos
parámetros define cuál es la salida del controlador en el tiempo cero. Si el controlador
se establece en variables de desviación, se recomienda usar estos valores en cero y
especificar el bias por fuera del bloque.
• External reset: esta opción permite establecer una condición que reinicializa el in-
tegrador y el filtro diferencial. Cuando se activa, en el bloque de PID aparece una
nueva entrada, donde se debe especificar la condición que activa la reinicialización de
los parámetros a las condiciones especificadas en Integrator Initial condition y Filter
Initial condition en los puertos I0 y D0. La señal de activación puede ser una señal
creciente, decreciente o de nivel; aquí se podría seleccionar una opción para ignorar
la reinicialización durante el proceso de linealización.
• Tracking mode: es una opción que permite a la salida del controlador seguir una señal
especificada por el puerto TR (tracking) que aparece en el bloque de PID al activarla.
En este caso, la diferencia entre la señal de referencia y la salida del controlador se
retroalimentan a la entrada del integrador. Esta opción no es de uso convencional,
pero ofrece aplicaciones para evitar windup en estructuras de control que tienen varios
lazos de control o cuando se desea tener transiciones suaves en sistemas que alternan
el uso entre dos controladores.
• Limit output: por defecto se tiene -infinito y +infinito como límites inferior y superior
de la señal de salida del controlador. En esta casilla, se pueden modificar indepen-
dientemente estos valores para evitar conflicto con el modelo o los actuadores. Níotese
que la implementación del controlador está en variables de desviación, es decir que la
señal de salida del controlador sí puede tomar valores negativos siempre y cuando no
viole condiones de implementación una vez el bias se adciona a la salida de este.
• Anti-windup: esta opción permite activar la acción de “borrado de memoria” de la
acción integral del controlador cuando se presenta el fenómeno de Reset windup, a
través de dos mecanismos. El primero se denomina Clamping (represión), en el que
se introduce un circuito de retroalimentación análogo al mostrado en la sección 5.5.1.
Para esto, se retroalimenta al integrador la diferencia entre la señal saturada y la
insaturada de salida del controlador. Este mecanismo se recomienda para plantas
con tiempos muertos pequeños, pues puede tener un desempeño pobre para plantas
con tiempos muertos grandes. El segundo mecanismo se denomina back-calculation,
y corresponde a la adición a la estructura anterior de una ganancia a la señal re-
troalimentada. Se recomienda utilizar el mismo tiempo integral para esta ganancia
217
Capítulo 5. Diseño y análisis de estructuras de control
√
(I = K/τi ) para controladores PI o I ∗ D para las plantas con tiempos muertos
relativamente altos.
h) Ventana Data types: las opciones por configurar en esta ventana están asociadas a la
herramienta complementaria llamada Open Fixed-Point Tool, a la cual puede accesarse
directamente desde Data types. Su uso está más allá del alcance de este libro.
i) Ventana State Attributes: esta es una ventana con parámetros útiles cuando se desea ge-
nerar código. Aquí se pueden asignar nombres y atributos para diferentes partes del con-
trolador. Por ejemplo, se puede asignar un nombre al estado asociado al integrador o al
filtro diferencial para controladores continuos.
A pesar de la gran cantidad de opciones que posee el bloque del controlador PID, su uso es
relativamente sencillo e intuitivo una vez se tienen claros los fundamentos del controlador y el
significado de sus parámetros. En el ejemplo 5.9, se establece un paralelo de la implementación
realizada en el caso de estudio anterior, donde no solo se evidencia la simple configuración del
bloque de PID, sino el potencial de la heramienta para la sintonización de los parámetros del
controlador.
Ejemplo 5.9.
Determinar la respuesta dinámica del sistema no lineal del reactor Van de Vusse utilizando
un controlador PID con acción antireset windup, usando el bloque predeterminado de Si-
mulink para control PID de un grado de libertad, con el fin de evaluar la herramienta de
autosintonización del controlador y su influencia en la respuesta del sistema ante diferentes
perturbaciones (i.e., cambio de set-point y perturbación en entrada) y la presencia de ruido.
Como estructura de control, se establece controlar la concentración de salida del componente
B, utilizando como variable manipulada la concentración de entrada del reactivo A (CAf ).
Las condiciones iniciales del reactor corresponden al estado estable óptimo (i.e., D̄ = 4/7
(min−1 ), C̄Af = 10 (mol/l), C̄A = 3 (mol/l) y C̄B = 1,117 (mol/l)).
Solución
Siguiendo la misma línea de los controladores hasta ahora implementados, se propone utilizar
la estrategia de programación con Matlab como programa maestro y usar Simulink para la
implementación del reactor y el controlador.
218
Capítulo 5. Diseño y análisis de estructuras de control
transformación por realizar para convertir los parámetros de un PID acoplado a uno desaco-
plado. Esta transformación se puede llevar a cabo directamente en el código de Matlab o en
el bloque de Simulink.
%--------------------------------------------------------------------------
% Programa para realizar control PID del reactor VDV
% PID usando bloque de Simulink
% Modelo no lineal en Simulink
% Programa maestro en Matlab
clc
clear
close all
% Definiciones de controlador
sim('VDV_NLsim_PID_control_box')
219
Capítulo 5. Diseño y análisis de estructuras de control
%% Resultados
%--------------------------------------------------------------------------
Ahora, en Simulink se prepara el programa mostrado en la figura 5.32, que contiene el modelo
del reactor Van de Vusse implementado como una S-function, el controlador PID usando el
bloque de Simulink y las perturbaciones empleadas anteriormente en los controladores.
220
Capítulo 5. Diseño y análisis de estructuras de control
Figura 5.32. Implentación en Simulink del control PID de la concentración de B empleando CAf
como variable manipulada, utilizando el bloque de PID con un grado de libertad
Se observa que la implementación en Simulink es casi igual que en el caso anterior, al sustituir
el bloque de PID antireset windup por un bloque con la leyenda PID(s) con un extraño
símbolo al final de la palabra (se explicará más adelante). Una vez más, nótese que a la salida
del controlador se le adiciona el bias p̄.
Dentro del bloque PID controller se realiza la configuración del controlador PID utilizando
la implementación en paralelo en tiempo continuo. Por tanto, la ecuación del controlador será:
!
P (s) 1 N
=P +I +D (5.17)
E(s) s 1 + N 1s
En la ventana Main se escoge la fuente interna y se introducen los parámetros del controlador.
Aquí se puede elegir entre: utilizar los parámetros obtenidos de la función PID_tuner, leer los
parámetros desde el workspace o, como en este caso, utilizar la herramienta de sintonización
de Simulink para obtenerlos. La implementación modificada se muestra en la figura 5.33.
221
Capítulo 5. Diseño y análisis de estructuras de control
Figura 5.33. Ventana inicial de configuración de un controlador PID para el control de la concen-
tración CB usando CAf (Nota: los parámetros del controlador que aparecen fueron obtenidos con la
herramienta de sintonización de Matlab)
222
Capítulo 5. Diseño y análisis de estructuras de control
En la parte superior de la ventana de Autotuning se ven opciones para la planta (no hay
que modificar nada aquí), el controlador, la herramienta de diseño y herramientas de la sin-
tonización. En este punto es importante detenerse un poco y explorar dichas opciones. En
las opciones del controlador, Simulink recomienda por defecto una sintonización balanceada
entre una prueba de seguimiento de set-point y una de rechazo a perturbaciones, debido a
que la sintonización tiene un carácter antagónico. Cuando se oprime el botón de opciones del
controlador se despliega una ventana en la que se puede escoger entre las opciones Balanced,
Reference tracking e Input disturbance rejection. Seleccionando cada una de ellas se evidencia
el impacto en la respuesta del sistema en la figura principal.
Figura 5.35. Ventana de sintonización de controladores PID para el caso del ejemplo
En la siguiente pestaña, llamada Design, se puede seleccionar entre dominio de tiempo y fre-
cuencia, donde se recomienda seleccionar tiempo. Adicionalmente, es posible adicionar una
serie de figuras a la ventana principal para el análisis del sistema. Las opciones disponibles
223
Capítulo 5. Diseño y análisis de estructuras de control
se muestran en la figura 5.36, las cuales incluyen: Plant, Open-loop, Reference tracking, Con-
troller effort, Input disturbance rejection, Output disturbace rejection, entre otras. Cuando se
selecciona una opción, esta figura aparece compartiendo espacio con la figura por defecto de
la respuesta en lazo cerrado ante una perturbación escalón del set-point.
Figura 5.36. Opciones de figuras que se pueden analizar para la respuesta en lazo cerrado del sistema
sintonizado
Por ejemplo, se desea evaluar tanto la respuesta del sistema con una prueba de cambio de
set-point y una prueba de rechazo a perturbaciones, utilizando la sintonización balanceada y
una sintonización adaptada para cambios de set-point. Para esto, se cambia en la ventana de
diseño la opción de sintonización Balanced a la opción Reference tracking, y en la ventana
de Add plot se adiciona una figura de Input disturbance rejection. Cuando se realiza esto,
automáticamente aparece una nueva figura en la ventana de la herramienta de sintonización
y las respectivas respuestas como se observa en la figura 5.37. Nótese que los parámetros
del controlador que se presentan en la parte inferior han cambiado en comparación con la
figura 5.35 y, para actualizar el bloque de PID principal, hay que oprimir el botón llamado
Update block ubicado en la parte superior derecha.
224
Capítulo 5. Diseño y análisis de estructuras de control
Figura 5.37. Evaluación de una sintonización alternativa a la obtenida por defecto por la herramienta
Autotuning de Simulink
Si se desea evidenciar los parámetros del controlador diseñado inicialmente con los nue-
vos resultados cuantitativos del desempeño de estos controladores (i.e., Rise time, Settling
time, Overshoot, Peak time, closed loop stability, entre otros), se puede oprimir la opción
Show parameters en la parte superior y obtener los resultados mostrados en la figura 5.38.
Figura 5.38. Resultados de la sintonización por defecto y modificada hacia pruebas de seguimiento
de set-point: parámetros de los controladores e índices de desempeño de los controladores
225
Capítulo 5. Diseño y análisis de estructuras de control
Por otro lado, la herramienta de sintonización permite un ajuste manual de las caracterís-
ticas del desempeño del sistema en lazo cerrado en relación con Response Time (seconds) y
Trasnsient Behavior. Para Response Time (seconds) se puede ajustar el controlador más lento
o más rápido y para Trasnsient Behavior si se require un controlador más agresivo o más
robusto. Es posible intuir que para cada configuración ambos extremos son antagónicos; por
ejemplo: si se desea evaluar la posibilidad de tener una respuesta en lazo cerrado más rápida,
se puede ajustar Response Time (seconds) ya sea en la barra o directamente modificando el
número que aparece al lado derecho de la barra. Como prueba se movió un 25 % el indicador
de Response time (seconds) hacia la derecha, lo que implicó que el parámetro cambiara de
0.8081 a 0.244. Los resultados se muestran en la figura 5.39.
Figura 5.39. Resultados de la sintonización por defecto y modificada hacia pruebas de seguimiento
de set-point: parámetros de los controladores e índices de desempeño de los controladores
Tanto las configuraciones por defecto de diseño del controlador como las manuales, abren una
gama de posibilidades para implementar controladores y evaluarlos sin tener que recurrir a
simulaciones individuales desde el programa principal. Esto permite identificar potenciales
candidatos de una manera más rápida para luego ser evaluados y comparados. Una vez se
está satisfecho con el resultado, se actualizan los parámetros del controlador y se cierra la
ventana de sintonización.
Para las otras opciones de configuración del bloque del controlador PID, solamente se modifica
en este caso Ouput saturation para evitar que la variable manipulada tome valores superiores
a Caf max e inferiores a cero como se muestra en la figura 5.40, y se activa la opción Antireset
windup. Al activar esta opción, en el bloque principal del PID aparecerá una marca especial
en el nombre como ya se mencionó.
226
Capítulo 5. Diseño y análisis de estructuras de control
Figura 5.40. Configuración del bloque de PID en la opción Output Saturation incluyendo límites
inferior y superior para la variable manipulada y la acción Antireset windup
Cuando se ejecuta el programa para el escenario donde se tiene un solo cambio de set-point y
una perturbación en la tasa de dilución se obtiene (figura 5.41):
Figura 5.41. Respuesta del sistema de control PID sintonizado con la herramienta Autotuning ante
cambios simples de set-point y tasa de dilución
227
Capítulo 5. Diseño y análisis de estructuras de control
El caso anterior puede ser evaluado en presencia de ruido, como se muestra en la figura 5.42.
Figura 5.42. Respuesta del sistema de control PID sintonizado con la herramienta Autotuning ante
cambios simples de set-point y tasa de dilución cuando se aplica ruido
Otra opción es realizar la prueba de set-point tracking más el rechazo de perturbaciones como
se muestra en la figura 5.43.
Figura 5.43. Respuesta del sistema de control PID sintonizado con la herramienta Autotuning para
la prueba de set-point tracking más rechazo a perturbación en la tasa de dilución
228
Capítulo 5. Diseño y análisis de estructuras de control
Reto 1: seleccionar las diferentes opciones de diseño del controlador para respuestas
mejoradas ante perturbaciones o pruebas de cambio en set-point para valorar el carácter
antagónico de las sintonizaciones.
El método de Root locus determina y grafica las trayectorias de los polos de uno o más sistemas
SISO en lazo cerrado en función de una ganancia (K) de un controlador proporcional con
retroalimentación negativa. En otras palabras, el método de Root locus evalúa el efecto en la
localización de los polos de un sistema al variar la ganancia retroalimentada. La localización
de los polos indican directamente la estabilidad en lazo cerrado (MatWorks, 2022q). Si se
tiene el siguiente sistema:
num(s)
sys(s) = (5.18)
den(s)
Los polos del sistema en lazo cerrado para un controlador proporcional vienen dados por las
raíces del polinomio:
El objetivo del método es entonces determinar las raíces del polinomio anterior en función
de la ganancia del controlador proporcional (K) cuando varía de cero a infinito. De manera
conveniente, el método ajusta el vector de ganancias para producir una respuesta suavizada.
Teniendo en cuenta que se utiliza el Zero-pole map, los polos aparecen el la gráfica como “x”.
Como alternativa a la figura de localización de polos, se puede utilizar el método para que se
entregen el vector de ganancias y las raíces, o para que el usuario especifique las ganancias de
229
Capítulo 5. Diseño y análisis de estructuras de control
interés.
%----------------------------------------------------------------------------
rlocus(sys); % Graficar el resultado de localización de polos
[r,k] = rlocus(sys); % Obtener polos y ganancias
%----------------------------------------------------------------------------
Es preciso tener cuidado en cuanto al modo en que se construye el sistema lineal cuando está
compuesto por diversas funciones de transferencia. Los casos de referencia se muestran en la
figura 5.44 (Nota: ver reto 3 en el ejemplo 5.10).
Figura 5.44. Ejemplos de representaciones de lazos cerrados con retroalimentación negativa para
ser utilizados por la función rlocus (Nota: observe cómo el producto de funciones de transferencia
representan el proceso, actuadores y sensores)
La figura que se obtiene usando la función rlocus es especial, ya que permite visualizar in-
formación directamente desde la figura dando clic sobre el punto de interés. Allí aparecerá la
ganancia (K), el polo, el factor de amortiguamiento (Damping), el sobreimpulso (Overshoot
( %)) y la frecuencia (rad/s). Además, el punto en cuestión puede moverse y se evidencia el
cambio de la ganancia.
Para ejemplificar el uso de esta función de Matlab se presentan dos casos de estudio: en el
primero se muestra el uso de la función rlocus y en el segundo se evalúa la estabilidad de un
controlador más complejo.
230
Capítulo 5. Diseño y análisis de estructuras de control
Ejemplo 5.10.
Se desea conocer los límites de la ganancia de un control proporcional (K) para tener esta-
bilidad en lazo cerrado de un sistema lineal compuesto por varias funciones de transferencia
representando: el proceso, la dinámica del actuador y la dinámica del sensor. Las funciones
de transferencia son, respectivamente:
1
Gp(s) = (5.20)
5s + 1
1
Gv(s) = (5.21)
2s + 1
1
Gm(s) = (5.22)
s+1
El diagrama de bloques del sistema se muestra en la figura 5.45. Nótese que la ganancia K
deberá ser sustituida por valores numéricos.
Figura 5.45. Representación del lazo cerrado de la estructura de control en cuestión para el análisis
de estabilidad de un control proporcional con ganancia (K)
Solución
231
Capítulo 5. Diseño y análisis de estructuras de control
hacen que el sistema sea estable). Para esto se evalúa la respuesta paso del sistema utilizando
la función step únicamente sobre las funciones en el lazo abierto, sin incluir el sensor ni el
controlador. Posteriormente, se utiliza la función rlocus para graficar directamente la loca-
lización de los polos del lazo cerrado. El paso clave para esta función es la manera en que
se introducen las funciones de transferencia; la ganancia del controlador proporcional no se
incluye y debe prestarse atención a cómo se ponen las funciones de transferencia, junto con el
sensor en este caso, que correspondería al lazo cerrado. Originalmente, el código se construye
hasta aquí y se ejecuta para deteminar si hay una ganancia de controlador proporcional que
haga inestable el lazo cerrado; cabe la posibilidad de que no exista y el sistema sea estable
para cualquier ganancia.
El código de Matlab con las secciones descritas, para las cuales se han introducido algunas
configuraciones especiales para las figuras, es:
%--------------------------------------------------------------------------
% Programa para determinar estabiliad en lazo cerrado
% Caso: sistema de varias funciones de transferencia
% Métodos: rlocus y pole/zero location
clc
clear all
close all
%% Etapa de evaluación
232
Capítulo 5. Diseño y análisis de estructuras de control
figure
rlocus(Km*Gp*Gv*Gm); % Prueba rlocus - lazo cerrado
ylabel('Eje imaginario','FontWeight','bold','FontSize',14)
xlabel('Eje real','FontWeight','bold','FontSize',14)
title('Root locus','FontWeight','bold','FontSize',14)
set(gca,'FontSize',12,'FontWeight','bold')
%% Etapa de validación
figure
pzplot(CL1,CL2,CL3)
legend('Kc=1','Kc=12.6','Kc=15')
ylabel('Eje imaginario','FontWeight','bold','FontSize',14)
xlabel('Eje real','FontWeight','bold','FontSize',14)
title('Zero-Pole Map','FontWeight','bold','FontSize',14)
set(gca,'FontSize',12,'FontWeight','bold')
figure
subplot(311)
step(CL1)
ylabel('Respuesta','FontWeight','bold','FontSize',14)
set(gca,'FontSize',12,'FontWeight','bold')
grid on
subplot(312)
step(CL2)
ylabel('Respuesta','FontWeight','bold','FontSize',14)
set(gca,'FontSize',12,'FontWeight','bold')
grid on
subplot(313)
step(CL3)
ylabel('Respuesta','FontWeight','bold','FontSize',14)
set(gca,'FontSize',12,'FontWeight','bold')
grid on
Utilizando la primera parte del programa principal, se obtiene la respuesta paso del sistema
233
Capítulo 5. Diseño y análisis de estructuras de control
en lazo abierto como se muestra en la figura 5.46. Como característica destaca que se puede
obtener información de los valores de respuesta y de tiempo haciendo clic sobre la figura. Por
otro lado, si se hace doble clic sobre la figura se despliega la ventana de editor de propiedades
donde se pueden modificar propiedades de estilo y características de la respuesta (figura 5.47).
Figura 5.47. Modificación de las propiedades de estilo de la figura de respuesta paso con el comando
step
Como producto del método de rlocus se obtiene la figura 5.48, donde se ha ubicado el cursor
234
Capítulo 5. Diseño y análisis de estructuras de control
Figura 5.48. Respuesta del método rlocus para determinación de trayectoria de polos en lazo
cerrado
Tomando como resultado que la ganancia última es K = 12.6, se introduce la segunda parte
del código con el fin de construir lazos cerrados para el sistema con unas ganancias de control
proporcional iguales a:
%--------------------------------------------------------------------------
Kc1 = 1; % Ganancia de controlador P (Estable)
Kc2 = 12.6; % Ganancia de controlador P (Críticamente estable)
Kc3 = 15; % Ganancia de controlador P (Inestable)
%--------------------------------------------------------------------------
Para luego obtener la representación de los polos y ceros en el diagrama complejo imaginario
mostrado en la figura 5.49. Se observa que para la ganacia de CL1, el sistema sería estable;
para CL2 con K = 12.6, los polos del sistema se ubican sobre el eje correspondiente a una
parte real de cero, lo que lo hace críticamente estable y, para el caso de CL3, el sistema debería
ser inestable. Se muestra en la figura que el factor de amortiguamiento para este sistema es
cero, lo que valida que es crítamente amortiguado.
235
Capítulo 5. Diseño y análisis de estructuras de control
Figura 5.49. Representación de los polos y ceros en el plano complejo imaginario de los tres sistemas
CL1, CL2, CL3, que incorporan diferentes ganancias de un controlador proporcional
Para estos casos, se valida la respuesta del sistema ante una perturbación escalón, la cual se
muestra en la figura 5.50.
Figura 5.50. Respuesta de los sistemas CL1, CL2 y CL3 en lazo cerrado ante una perturbación escalón
unitario
236
Capítulo 5. Diseño y análisis de estructuras de control
Reto 1: demostrar analíticamente con el método de Routh que la ganancia última del
sistema es K = 12.6.
Reto 2: demostrar analíticamente con el método de ganancia última que K = 12.6.
Como segundo caso, se tiene la determinación de la estabilidad en lazo cerrado para el sistema
de reacción Van de Vusse cuando se implementa un controlador PID. Esto es particularmente
importante, pues la acción integral disminuye la estabilidad de un sistema mientras la acción
derivativa otorga estabilidad al lazo de control.
Ejemplo 5.11.
Se desea analizar la estabilidad en lazo cerrado para el reactor Van de Vusse tanto para un
control proporcional como para un PID, teniendo en cuenta que el control integral disminuye
la estabilidad del sistema mientras que el control derivativo la incrementa. Como estructura
de control, se establece controlar la concentración de salida del componente B, utilizando
como variable manipulada la concentración de entrada del reactivo A (CAf ). La condiciones
iniciales del reactor corresponden al estado estable óptimo (i.e., D̄ = 4/7 (min−1 ), C̄Af = 10
(mol/l), C̄A = 3 (mol/l) y C̄B = 1,117 (mol/l)).
Solución
237
Capítulo 5. Diseño y análisis de estructuras de control
El código de Matlab con las secciones descritas, para las cuales se han introducido algunas
configuraciones especiales para las figuras, es:
%--------------------------------------------------------------------------
% Programa para determinar estabiliad en lazo cerrado
% Caso: reactor Van de Vusse
% Métodos: rlocus y pole/zero location
clc
clear all
close all
% Estado estable
a11=-D-k1-2*k3*Cao;
a12=0;
a21=k1;
a22=-D-k2;
b11=Cafo-Cao;
b12=D;
b21=-Cbo;
b22=0;
238
Capítulo 5. Diseño y análisis de estructuras de control
Sys_ss=ss(A,B,C,D,'StateName',{'Ca';'Cb'},'InputName',{'D';'Caf'},...
'OutputName',{'Ca';'Cb'}); % Crea espacio de estados
figure
step(Sys); % Prueba paso del sistema
ylabel('Respuesta','FontWeight','bold','FontSize',14)
xlabel('Tiempo','FontWeight','bold','FontSize',14)
title('Respuesta paso','FontWeight','bold','FontSize',14)
set(gca,'FontSize',12,'FontWeight','bold')
figure
rlocus(Sys); % Prueba rlocus lazo cerrado
ylabel('Eje imaginario','FontWeight','bold','FontSize',14)
xlabel('Eje real','FontWeight','bold','FontSize',14)
title('Root locus','FontWeight','bold','FontSize',14)
set(gca,'FontSize',12,'FontWeight','bold')
figure
pzplot(CL1); % Diagrama polo-zero
ylabel('Eje imaginario','FontWeight','bold','FontSize',14)
xlabel('Eje real','FontWeight','bold','FontSize',14)
title('Zero-Pole Map','FontWeight','bold','FontSize',14)
set(gca,'FontSize',12,'FontWeight','bold')
figure
step(CL1)
ylabel('Respuesta','FontWeight','bold','FontSize',14)
xlabel('Tiempo','FontWeight','bold','FontSize',14)
title('Respuesta paso','FontWeight','bold','FontSize',14)
set(gca,'FontSize',12,'FontWeight','bold')
239
Capítulo 5. Diseño y análisis de estructuras de control
Sys =
0.4762
---------------------
s^2 + 4.643 s + 5.382
Nótese que debido a las propiedades especiales asignadas al objeto de espacio de estados, la
entrada y la salida de interés se muestran en los resultados. El resultado de la perturbación
escalón a esta función de transferencia se presenta en la figura 5.51.
Figura 5.51. Valoración de la estabilidad del sistema en lazo abierto usando la prueba paso
Con la prueba de localización de polos por medio del rlocus se obtiene la figura 5.52. A
diferencia del ejemplo anterior, en este caso no se cruza el plano de parte real negativa a
positiva, los que indica que para este sistema no hay límite de estabilidad.
240
Capítulo 5. Diseño y análisis de estructuras de control
Figura 5.52. Trayectorias de la localización de los polos del lazo cerrado para un controlador propor-
cional y el sistema Van de Vusse
Para el controlador PID, el objeto del controlador y el lazo cerrado tienen las siguientes
representaciones:
C1 =
1 s
Kp + Ki * --- + Kd * --------
s Tf*s+1
CL1 =
241
Capítulo 5. Diseño y análisis de estructuras de control
La gráfica del los polos del sistema y respuesta escalón se muestra en las figuras 5.53 y 5.54.
Figura 5.53. Representación de los polos y ceros en el plano complejo imaginario del lazo cerrado del
controlador PID
Figura 5.54. Respuesta del sistema en lazo cerrado ante una perturbación escalón unitario
242
Capítulo 5. Diseño y análisis de estructuras de control
243
Capítulo 6
Caso de estudio
Resumen
En este capítulo se presenta el modelo básico utilizado para implementar los métodos y herra-
mientas abordados en este libro. El caso de estudio escogido es de gran pertinencia para Ingeniería
Química, con una baja complejidad en aras que el estudiante no tenga limitaciones para su im-
plementación. Sin embargo, este caso de estudio tiene los elementos básicos que se encuentran
en casos de estudio más retadores.
k k
A −−−−1−−→ B −−−−2−−→ C
k
2A −−−−3−−→ D
245
Capítulo 6. Caso de estudio
dCA F 2
= (CAf − CA ) − k1 CA − k3 CA (6.1)
dt V
dCB F
= (CBf − CA ) + k1 CA − k2 CB
dt V
Para el modelo se conocen los siguientes parámetros, variables de entrada y condiciones ini-
ciales:
246
Capítulo 6. Caso de estudio
247
Referencias
Referencias
249
Referencias
MatWorks (2022h). hsvplot: plot Hankel singular values and return plot handle. https:
//se.mathworks.com/help/control/ref/dynamicsystem.hsvplot.html.
MatWorks (2022j). linmod: extract continuous-time linear state-space model around operating
point. https://la.mathworks.com/help/simulink/slref/linmod.html.
MatWorks (2022k). linmod2: extract continuous-time linear state-space model around opera-
ting point. https://la.mathworks.com/help/simulink/slref/linmod2.html.
250
Referencias
Migoni, G. y Kofman, E. (2007). Linearly implicit discrete event methods for stiff ODEs.
Part I: theory [conference]. In RPIC 2007-XII Reunion de Trabajo en Procesamiento de la
Información y Control, Río Gallegos, Argentina.
Rosen, C., Vrecko, D., Gernaey, K., Pons, M.-N., y Jeppsson, U. (2006). Implementing ADM1
for plant-wide benchmark simulations in Matlab/Simulink. Water Science and Technology,
54(4), 11–19.
Seborg, D. E., Edgar, T. F., Mellichamp, E. A., y Doyle, F. J. (2011). Process dynamics and
control (third edition). John Wiley & Sons.
251
Índice alfabético
253
Índice alfabético
Funciones separadas 25, 32-33, 84, 88, Graphical User Interface (GUI) 126, 162,
108 175
eigenvalores 115, 188-189
error 26-27, 51-52, 58-59, 89, 91, 99, 114, H
133-134, 136, 139-141, 149, 154-155,
Hankel 148-151, 153, 158, 248
159, 164, 168, 174-175, 185-187, 195-
Heaviside 96-97, 123, 166
196, 203-206, 210
espacio de estados 11, 65, 83-84, 92, 99, 104, I
106-109, 111, 117, 122, 125, 128-130,
132, 134-138, 147-149, 151-152, 155, Identificación de sistemas 173
158-159, 169-170, 173, 175, 212, 238 Toolbox 175-182
representación 66, 92, 103, 106-107, 116, impulse 122-123, 126
119-120, 132, 137,169-170, 184, 235 Integrator 38-39, 41, 46, 204-205, 215
estabilidad Interior-point 74-75
BIBO 115, 235 Internal Model Control (IMC) 206-207, 236-
inherente 114-116 237
lazo cerrado 115, 161-162, 183, 206, 222, Interpreted function 26, 38-41, 46, 55
224, 227, 229-230, 233-236, 239-240
numérica 27, 114-115 J
estado estable 65, 69, 70, 74, 77, 79, 83, 85,
Jacobiano 88-89, 115
87-88, 90, 115, 141-143
Simbólico 84-85
solución 66-68, 81-82, 92, 144
Numérico 87, 91
exactitud 36, 46, 87-88, 90-91, 187
exitflag 67-69, 73, 75 K
F Karush-Kuhn-Tucker (KKT) 71
Filosofía de programación 10, 47 L
Matlab como ambiente maestro 24
Simulink como ambiente maestro 24, 38, Laplace 7, 22, 65-66, 94-99, 101-103, 105,
47 107, 113, 204
Flag 46, 48-51, 55-59, 67-71, 73, 75, 77 Transformada 66, 94-98, 105, 107
Filtro derivativo 203, 205, 211-212, 236 Transformada inversa 66, 98, 101-103
fmincon 72-75, 77-78, 247 Levenberg-Marquardt 67
Fortran 23, 45 linmod 91-93, 108, 248
freqresp 125 lsim 52, 58-59, 79, 81, 92-93, 108, 122-123,
fsolve 66-68, 70-73 126, 165, 170, 172-173, 197-198, 207,
Función de transferencia 105-111, 131, 142, 217
143, 146, 164-167, 179-181, 204, 212,
230, 235-236, 238 M
Representación 65, 132, 147, 162, 235
math operators 12, 21
G matrices 16, 36-37, 82-83, 85-87, 90-92, 103-
104, 107-108, 117, 120-121, 126, 130-
Gain 21-22, 183, 204-205, 249 132, 148-151, 169-171, 184-186, 188,
Ganacia del controlador 204, 227, 230 191-192
Gram 149-150 método Runge-Kutta 26, 114
254
Índice alfabético
255
Índice alfabético
T
tiempo de asentamiento (settling time) 141,
143-144, 223
tiempo derivativo 236
tiempo de elevación (rise time) 141, 223
tiempo integral 204, 215, 236
tiempo pico (time to first peak) 141
trim 78-79, 81, 249
trust-region 67-68, 70-71
trust-region-dogleg 67-68, 70-71
trust-region-reflective 74, 78
V
Valores singulares véase Hankel
Z
Zero-pole map 227, 231, 237
256
Óscar Andrés Prado-Rubio es ingeniero químico de
la Universidad Nacional de Colombia (2005), máster en
Automatización de Procesos Industriales de la Univer-
sidad de Ibagué (2007) y PhD en Ingeniería Química
y Bioquímica de la Universidad Técnica de Dinamarca
(2010). Cuenta con un postdoctorado en la Universidad
Técnica de Dinamarca (2012) y se desempeñó como es-
pecialista en membranas y como gerente de proyectos
en LiqTech International S.A. (2013), Dinamarca. Pos-
teriormente, se vinculó a la Universidad Nacional de Co-
lombia como profesor asociado y, recientemente, como
director de Investigación y Extensión de la Sede Mani-
zales (2018-2021). Prado-Rubio ha sido profesor visitante
en la Universidad de Newcastle (2016), Universidad de
Guanajuato (2017), Universidad Técnica de Dinamarca
(2018) y en el Instituto Nacional de las Ciencias Aplica-
das de Francia (2018). Actualmente, se desempeña como
Investigador Senior en la Universidad Técnica de Dina-
marca (2022-2024).