SEP TECNOLÓGICO NACIONAL MÉXICO

INSTITUTO TECNOLÓGICO DE TIJUANA

DEPARTAMENTO DE INGENIERÍA ELÉCTRICA Y ELECTRÓNICA

POSGRADO EN CIENCIAS DE LA INGENIERÍA

INGENIERÍA BIOMÉDICA

ASIGNATURA:
MODELADO DE SISTEMAS FISIOLÓGICOS

BioMath - When engineering meets biology

BioMath
ELABORADO POR:
DR. PAUL ANTONIO VALLE TRUJILLO
[Link]@[Link]

TIJUANA, B.C., MÉXICO

 i

Contenido

Unidad Página

1. Información general 1
1.1. Competencias previas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2. Unidades . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.3. Porcentajes de evaluación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.4. Calendarización . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.5. Tabla de voltajes y corrientes para cada componente . . . . . . . . . . . . . . . . . . . . . 2
1.6. Tabla de transformadas de Laplace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.7. Referencias . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2. Sistemas de control retroalimentados 5
2.1. Elementos de sistemas de control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.2. Función de transferencia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.3. Diagramas de bloques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.3.1. Diagrama de bloques de un sistema en lazo cerrado . . . . . . . . . . . . . . . . . . 9
2.3.2. Función de transferencia en lazo cerrado . . . . . . . . . . . . . . . . . . . . . . . . 10
2.3.3. Reglas de simpli…cación de algebra de bloques . . . . . . . . . . . . . . . . . . . . . 11
2.4. Respuesta de un sistema ante distintas señales de entrada . . . . . . . . . . . . . . . . . . 13
2.4.1. De…nición de escalón unitario . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.4.2. De…nición de rampa unitaria . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.4.3. De…nición de impulso unitario . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.4.4. De…nición de función sinusoidal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.4.5. Características de la respuesta al escalón . . . . . . . . . . . . . . . . . . . . . . . . 16
2.5. Análisis del error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.6. El concepto de estabilidad . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.7. Criterios de estabilidad . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3. Controladores 23
3.1. Controlador PID . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.2. Controlador P . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.3. Controlador I . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.4. Controlador D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.5. Controlador PI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.6. Controlador PD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
3.7. Restador . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.8. Sumador . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
4. Modelado matemático 34
4.1. Función de transferencia, error y estabilidad del sistema de mecánica pulmonar . . . . . . 35
4.2. Modelo matemático de ecuaciones integro-diferenciales . . . . . . . . . . . . . . . . . . . . 36
4.3. Respuesta del sistema en lazo abierto . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
4.3.1. Multisim: lazo abierto . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
4.3.2. Laboratorio de Electrónica: lazo abierto . . . . . . . . . . . . . . . . . . . . . . . . 44
4.4. Diseño del controlador . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
4.4.1. Simulink: lazo abierto . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
4.4.2. Simulink: lazo cerrado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.4.3. Multisim: Lazo cerrado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
 ii

Contenido (Continuación)

Capítulo Página

4.5. Sistema de control en Python . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

5. Prácticas 77
5.1. Práctica: Diseño de controladores para sistemas de segundo orden . . . . . . . . . . . . . . 77
5.2. Práctica: Comparación entre los controladores P, I, PD, PI y PID . . . . . . . . . . . . . . 79
5.3. Práctica: Sistema musculoesquelético . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
5.4. Práctica: Sistema respiratorio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
5.5. Práctica: Sistema cardiovascular . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
5.6. Práctica: Sistema digestivo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
5.7. Práctica: Sistemas …siológicos en el espacio de estados: EDOs . . . . . . . . . . . . . . . . 91
5.8. Práctica: Regeneración de glóbulos rojos . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
6. Proyecto …nal 96
6.1. Objetivo del proyecto …nal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
6.2. Descripción del sistema . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
6.3. Modelo matemático . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
6.4. Simulaciones numéricas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
6.5. Implicaciones …siológicas (Conclusiones) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
 1

Capítulo 1

Información general

La competencia a desarrollar en esta asignatura es emplear el control clásico y las técnicas de

modelizado de sistemas para aplicarlas en procesos …siológicos y biológicos.

1.1. Competencias previas

1. Aplicar leyes de Kirchho¤, teorema de superposición y la transformada de Laplace.

2. Resolver ecuaciones diferenciales e integro-diferenciales y sistemas de ecuaciones.

3. Utilizar los ampli…cadores operacionales.

4. Comprender el funcionamiento de los sistemas …siológicos.

1.2. Unidades

1. Sistemas de control retroalimentados.

Competencia especí…ca a desarrollar: Comprender los conceptos básicos de la Ingeniería de Control,

su respuesta en el tiempo y criterios de estabilidad.

2. Modelado matemático.

Competencia especí…ca a desarrollar: Modelizar sistemas …siológicos y biológicos mediante

ecuaciones diferenciales, integro-diferenciales y función de transferencia para el análisis de su
respuesta.

3. Controladores.

Competencia especí…ca a desarrollar: Diseña controladores para sistemas mediante técnicas de

control clásico.

4. Proyecto …nal.

Competencia especí…ca a desarrollar: Modelizar un sistema que describa un proceso …siológicos o

biológicos y diseñar controlador aplicable a dicho sistema.
 2

1.3. Porcentajes de evaluación

Criterio Porcentaje
Prácticas 40 %
Mini-Exámenes 10 %
Examen 20 %
Proyecto: Póster 10 %
Proyecto: Artículo 20 %

1.4. Calendarización

Actividad Semanas
1: Prácticas 11

2: Examen 1

3: Elaboración del proyecto …nal 3

1.5. Tabla de voltajes y corrientes para cada componente

Componente Símbolo Voltaje Corriente

v (t)
Resistor R Ri (t)
R
Z
1 dv (t)
Capacitor C i (t) dt C
C dt
Z
di (t) 1
Inductor L L v (t) dt
dt L
 3

1.6. Tabla de transformadas de Laplace

Dominio del tiempo (t) Dominio de ‘s’

1: af (t) + bg (t) aF (s) + bG (s)

df (t)
2: f 0 (t) = sF (s) f (0)
dt
Z t
1
3: f (t) dt F (s)
0 s

4: Impulso unitario: (t) 1

1
5: Escalón unitario: 1
s

1
6: Rampa: t
s2

at 1
7: e
s+a

tk at 1
8: e
k! (s + a)k+1

tk 1
9:
k! sk+1

!
10: sin !t
s2 + !2

s
11: cos !t
s2 + !2
 4

1.7. Referencias

1. Khoo MC. Physiological Control Systems Analysis, Simulation, and Estimation. Wiley, 2018.

2. Ogata K. Ingeniería de Control Moderno. Prentice Hall, 2003.

3. Nise, NS. Control systems engineering. John Wiley & Sons, 2020.

4. Kuo BC. Sistemas de Control Automático. Prentice Hall, 1996.

5. Dorf RC and Bishop RH. Modern Control Systems. Prentice Hall, 2011.

6. Díaz M. Sistemas Lineales I. Instituto Tecnológico de Durango.

7. Milo R and Phillips R. Cell Biology by the Numbers. Garland Science, 2016.

8. Gar…nkel A, Shevtsov J and Guo Y. Modeling Life. The Mathematics of Biological Systems. Springer,
2017.

9. Britton NF. Essential Mathematical Biology. Springer, 2005.

10. Unbehauen H, editor. Control Systems, Robotics, and Automation -Volume III: System Analysis
and Control: Classical Approaches-III. EOLSS Publications; 2009.

11. Inverting Operational Ampli…er. Electronics Tutorials. [Link]. Last access

09/June/2022.

12. University of Michigan. Control tutorials for Matlab & Simulink. [Link]. Last access
09/June/2022.

13. Educational Figures Created with [Link]. Last access 09/June/2022.

14. Kind, T., Faes, T. J., Lankhaar, J. W., Vonk-Noordegraaf, A., and Verhaegen, M. Estimation
of three-and four-element windkessel parameters using subspace model identi…cation. IEEE
Transactions on Biomedical Engineering, 2010, 57(7), 1531–1538.

15. Tetschke M, Lilienthal P, Pottgiesser T, Fischer T, Schalk E, Sager S. Mathematical Modeling of

RBC Count Dynamics after Blood Loss. Processes. 2018; 6(9):157. doi: 10.3390/pr6090157.
 5

Capítulo 2

Sistemas de control retroalimentados

Un sistema de control es una interconexión de componentes cuyo objetivo es proporcionar una

respuesta deseada en un proceso especí…co de la denominada planta, o ‘el sistema a controlar’.

2.1. Elementos de sistemas de control

A continuación se presentan ocho conceptos básicos sobre los sistemas de control:

1. Sistema: Cualquier objeto físico, conjunto de elementos o estructuras compuestas que se van a
controlar.

2. Función de transferencia: Es una expresión que relaciona la salida y la entrada en un sistema.

3. Polos: Las raíces del denominador de la función de transferencia de un sistema.

4. Ceros: Las raíces del numerador de la función de transferencia de un sistema.

5. Proceso: Cualquier operación que se va a controlar.

6. Sistema lineal: Es un sistema en el que se puede aplicar el principio de superposición; el cual

establece que la respuesta producida por la aplicación simultánea de dos funciones excitadoras
distintas es la suma de las dos respuestas individuales.

7. Sistema no lineal: Es un sistema al que no se le puede aplicar el principio de superposición, no se

puede calcular la respuesta determinando una a la vez y sumando los resultados.

8. Perturbación: Es una señal que tiende a afectar negativamente el valor de la salida de un sistema.

Los elementos de un sistema de control son los siguientes:

1. La señal de entrada suele representar el valor deseado al cual se quiere llevar la respuesta del sistema,
también se conoce como referencia o respuesta deseada.

2. El sistema o planta es el conjunto de elementos o estructuras compuestas que se van a controlar.

3. El error o variable controlada es la diferencia entre la señal de entrada y la salida del sistema.
 6

4. El controlador es el dispositivo que, con base en el error, aplica la señal de control necesaria al
sistema para obtener el valor deseado en la salida.

5. La señal de control es la respuesta del controlador para eliminar el error entre la salida del sistema
y la señal de entrada.

6. La señal de salida es la respuesta que tiene el sistema a una señal de entrada, con o sin la acción de
un controlador.

Sistema de control en lazo abierto

Son sistemas de control en los que la salida no tiene efecto sobre la acción de control. En estos sistemas
la salida ni se mide ni se realimenta para compararla con la entrada.

Entrada Controlador Sistema Salida

Figura 2.1: Sistema de control en lazo abierto.

Sistema de control en lazo cerrado

Son sistemas de control en los que la señal de salida tiene efecto directo sobre la acción de control.

Error Se~
nal
Entrada + Controlador Sistema Salida
! de control

Retroalimentaci4
on

Figura 2.2: Sistema de control en lazo cerrado.

El controlador es la parte más importante en un sistema de control debido a que es el responsable de

su rendimiento. Un controlador es un dispositivo o un algoritmo cuyo propósito es mantener el valor de
la variable controlada en el valor deseado para el sistema.

El objetivo del controlador es eliminar el error entre la señal de entrada y la señal de salida del sistema.
 7

La con…guración del sistema de control en lazo cerrado se conoce como ‘compensación en serie o
cascada’ y es la más comúnmente utilizada, mientras que el controlador más utilizado en este tipo
de con…guración es el denominado controlador Proporcional-Integral-Derivativo (PID), el cual calcula
continuamente el valor del error entre la salida y la entrada de un sistema dinámico. El controlador PID
es un dispositivo que considera el presente, pasado y futuro del error estimado.

2.2. Función de transferencia

La función de transferencia de un sistema lineal invariante en el tiempo está de…nida como la relación
de la transformada de Laplace de la salida a la transformada de Laplace de la entrada bajo la suposición
de que todas las condiciones iniciales son cero. La función de transferencia es una expresión que relaciona
la salida y la entrada en términos de los parámetros del sistema.
Sea el sistema lineal invariante en el tiempo de…nido por la siguiente ecuación diferencial

dn y dn 1 y dy dm x dm 1 x dx
a0 + a1 + ::: + an 1 + an y = b0 m + b1 m 1 + ::: + bm 1 + bm x;
dtn dtn 1 dt dt dt dt

para
n m;

y en donde y es la salida del sistema y x es la entrada. La función de transferencia está dada por

L fy (t)g Y (s) b0 sm + b1 sm 1 + ::: + bm 1 s + bm

Función de transferencia : G (s) = = = ;
L fx (t)g x(0)=0 X (s) a0 sn + a1 sn 1 + ::: + an 1 s + an

donde los polos del sistema son las raíces del denominador de su función de transferencia,es decir, su
ecuación característica X (s) = 0 y los ceros son las raíces del numerador Y (s) = 0.
A partir del concepto de función de transferencia, es posible representar la dinámica de un sistema
mediante ecuaciones algebraicas en s. Si la potencia más alta de s en el denominador de la función de
transferencia es igual a n, el sistema se denomina sistema de n ésimo orden.

Pasos para obtener la función de transferencia

1. Plantear la ecuación integro-diferencial del sistema.

2. Obtener la transformada de Laplace suponiendo condiciones iniciales iguales a cero.

3. Encontrar la relación de la salida Y (s) con respecto a la entrada X (s).

 8

Características de la función de transferencia

La aplicación del concepto de función de transferencia está limitada a los sistemas descritos mediante
ecuaciones diferenciales lineales invariantes en el tiempo. La función de transferencia se usa extensamente
en el análisis y diseño de dichos sistemas.
A continuación, se presentan algunas características importantes relacionadas con la función de
transferencia:

1. La función de transferencia de un sistema es un modelo matemático porque es un método operacional

para expresar la ecuación diferencial que relaciona la variable de salida con la variable de entrada.

2. La función de transferencia es una propiedad de un sistema, independiente de la magnitud y

naturaleza de la entrada o función de excitación.

3. La función de transferencia incluye las unidades necesarias para relacionar la entrada con la salida;
sin embargo, no proporciona información acerca de la estructura física del sistema. Las funciones de
transferencia de muchos sistemas físicamente diferentes pueden ser idénticas.

4. Si se conoce la función de transferencia de un sistema, se estudia la salida o respuesta para varias

formas de entrada, con la intención de comprender la naturaleza del sistema.

5. Si se desconoce la función de transferencia de un sistema, puede establecerse experimentalmente

introduciendo entradas conocidas y estudiando la salida del sistema. Una vez establecida una función
de transferencia, esta proporciona una descripción completa de las características dinámicas del
sistema, a diferencia de su descripción física.

2.3. Diagramas de bloques

Un diagrama de bloques de un sistema es una representación grá…ca de las funciones que lleva a cabo
cada componente y el ‡ujo de las señales. El diagrama muestra las relaciones existentes entre los diversos
componentes. A diferencia de una representación matemática, un diagrama de bloques tiene la ventaja
de indicar el ‡ujo de las señales del sistema.
En un diagrama de bloques se enlazan todas las variables del sistema mediante bloques funcionales.
Un bloque es un símbolo que representa una operación matemática que se realiza sobre la señal para
producir la salida. Las funciones de transferencia de los componentes por lo general se introducen en los
bloques correspondientes los cuales se conectan mediante ‡echas para indicar la dirección del ‡ujo de la
señal. La señal solo puede pasar en la dirección de las ‡echas, por tanto, un diagrama de bloques de un
sistema de control muestra explícitamente una propiedad unilateral.
 9

La Figura 2.3 muestra un elemento del diagrama de bloques. La punta de ‡echa que señala el bloque
indica la entrada, y la punta de ‡echa que se aleja del bloque representa la salida.

R(s) G(s) C(s)

Figura 2.3: R (s) : Entrada, G (s) : Función de transferencia, C (s) : Salida.

Las dimensiones de la señal de salida del bloque son las dimensiones de la señal de entrada multiplicadas
por las dimensiones de la función de transferencia en el bloque.
El diagrama de bloques permite evaluar la contribución de cada componente al desempeño general
del sistema. En general, la operación funcional del sistema se aprecia con más facilidad si se examina el
diagrama de bloques que si se revisa el sistema físico mismo. Un diagrama de bloques contiene información
relacionada con el comportamiento dinámico, pero no incluye información de la construcción física del
sistema. En consecuencia, muchos sistemas diferentes y no relacionados pueden representarse mediante el
mismo diagrama de bloques.
En un diagrama de bloques la principal fuente de energía no se muestra explícitamente y el diagrama de
bloques de un sistema determinado no es único. Es posible dibujar varios diagramas de bloques diferentes
para un sistema, dependiendo del punto de vista del análisis.

2.3.1. Diagrama de bloques de un sistema en lazo cerrado

La Figura 2.4 muestra un ejemplo de un diagrama de bloques de un sistema en lazo cerrado. La salida
C(s) se realimenta al punto de suma, en donde se compara con la entrada de referencia R(s).
La naturaleza en lazo cerrado del sistema se indica con claridad en la …gura. La salida del bloque, C(s)
en este caso, se obtiene al multiplicar la función de transferencia G(s) por la entrada al bloque, E(s).

Sumador
E(s) Ramif icaci4
on
R(s) + G(s) C(s)
!

Retroalimentaci4
on

Figura 2.4: Diagrama de bloques de un sistema en lazo cerrado.

 10

Cualquier sistema de control lineal puede representarse mediante un diagrama de bloques formado
por puntos de suma, bloques y puntos de rami…cación.

2.3.2. Función de transferencia en lazo cerrado

Para el sistema que aparece en la Figura 2.5, la salida C(s) y la entrada R(s) se relacionan de la
siguiente manera:

C (s) = G (s) E (s) ;

B (s) = H (s) C (s) ;

E (s) = R (s) B (s) :

E(s)
R(s) + G(s) C(s)
!

B(s)

H(s)

Figura 2.5: Análisis de un sistema de lazo cerrado.

Sustituyendo E (s) y B (s) en C (s) se tiene:

C (s) = G (s) [R (s) H (s) C (s)] ;

C (s) = G (s) R (s) G (s) H (s) C (s) ;

C (s) [1 + G (s) H (s)] = G (s) R (s) ;

C (s) G (s)
= :
R (s) 1 + G (s) H (s)
 11

2.3.3. Reglas de simpli…cación de algebra de bloques

Es importante señalar que los bloques pueden conectarse en serie, sólo si la entrada de un bloque
no se ve afectada por el bloque siguiente. Si hay efectos de carga entre los componentes, es necesario
combinarlos en un bloque único.
Cualquier cantidad de bloques en cascada que representen componentes sin carga puede sustituirse con
un solo bloque, cuya función de transferencia sea simplemente el producto de las funciones de transferencia
individuales.

Figura 2.6: Reglas del álgebra de los diagramas de bloques.

 12

La simpli…cación de un diagrama de bloques mediante reordenamientos y sustituciones reduce de

manera considerable la labor necesaria para el análisis matemático subsecuente. Sin embargo, conforme
se simpli…ca el diagrama de bloques, las funciones de transferencia de los bloques nuevos se vuelven más
complejas debido a que se generan polos y ceros nuevos.

Figura 2.7: Reglas del álgebra de los diagramas de bloques.

 13

2.4. Respuesta de un sistema ante distintas señales de entrada

En esta sección se de…nen cuatro señales básicas para determinar la dinámica de un sistema, es decir,
el escalón, el impulso, la rampa y una función sinusoidal.

2.4.1. De…nición de escalón unitario

Es una función discontinua cuyo valor es 0 para cualquier argumento negativo y 1 para cualquier
argumento positivo y el cero. Se de…ne de la siguiente forma
8
< 0; x<0
us (x) = ;
: 1; x 0

y la grá…ca se ilustra en la Figura 2.8.

0.8

0.6
us (x)

0.4

0.2

-10 -8 -6 -4 -2 0 2 4 6 8 10
x

Figura 2.8: Grá…ca del escalón unitario.

La transformada de Laplace de la función escalón está dada por

1
L [us (x)] = :
s

2.4.2. De…nición de rampa unitaria

Es una función elemental real de un solo argumento, continua y diferenciable en todo su dominio
excepto en un punto (el inicio de la rampa). Se de…ne de la siguiente forma
8
< 0; x<0
ur (x) = ;
: x; x 0
 14

y la grá…ca se ilustra en la siguiente Figura 2.9.

10

6
ur (x)

-10 -8 -6 -4 -2 0 2 4 6 8 10
x

Figura 2.9: Grá…ca de la rampa unitaria.

La transformada de Laplace de la rampa está dada por

1
L [ur (x)] = :
s2

2.4.3. De…nición de impulso unitario

Los sistemas del mundo real, como los mecánicos, físicos, eléctricos y …siológicos se pueden someter a
una fuerza externa de gran magnitud, que solamente actúa durante un período de tiempo muy corto. La
función denominada como impulso unitario se puede utilizar como un modelo matemático para representar
dicha fuerza o perturbación externa.
La función impulso, Delta de Dirac o función es una función generalizada que toma el valor cero
en todos los números reales excepto en el cero donde es in…nita, su integral sobre la recta real es 1. La
función se de…ne de la siguiente forma
8
< +1; x=0
(x) = ;
: 0; x 6= 0

y además satisface la siguiente identidad

Z +1
(x) dx = 1:
1

La función no es una función en el sentido tradicional porque ninguna función de…nida en los números
reales tiene estas propiedades. Al considerar lo anterior el impulso unitario se de…ne de la siguiente forma
8
< 1; x=0
(x) = ;
: 0; x 6= 0
 15

y la grá…ca se ilustra en la siguiente Figura 2.10.

0.8

0.6
/(x)

0.4

0.2

-10 -8 -6 -4 -2 0 2 4 6 8 10
x

Figura 2.10: Grá…ca del impulso unitario.

2.4.4. De…nición de función sinusoidal

Una señal sinusoidal (oscilación sinusoidal o señal sinusoidal) es aquella de…nida de la siguiente forma

'
u (t) = A sin (!t + ') = A sin ! t ; ! > 0;
!

las características de la función sinusoidal son las siguientes.

1. El coe…ciente A > 0 se denomina amplitud de la función.

2. El período T (medido en s) de la función sinusoidal está de…nido por el cociente

2
T = :
!

3. El coe…ciente ! representa un estiramiento o compresión horizontal, con este se modi…ca el período

!
de la función, también se conoce como frecuencia angular (medida en rad=s). El factor se conoce
2
como la frecuencia (f medida en Hz) de la función, es decir

!
f= ;
2

y de…ne el número de oscilaciones o ciclos completos por unidad de tiempo t.

4. El primer ciclo de la función se de…ne en un intervalo

' 2 '
t2 ; + :
! ! !
 16

'
5. La relación de indica el desplazamiento de fase, es decir, el desplazamiento en el eje horizontal.
!
La grá…ca para la función sinusoidal básica se ilustra en la siguiente Figura 2.11, con una amplitud
A, ! = 1 y ' = 0.

A/2
u(x)

-A/2

-A
0 :/2 : 3 :/2 2: 3: 4:
x

Figura 2.11: Grá…ca de la función sinusoidal.

2.4.5. Características de la respuesta al escalón

Las características de la respuesta al escalón se describen en la Figura 2.12. El valor …nal de la respuesta
es el nivel en estado estacionario, si la entrada del sistema es un escalón unitario el valor …nal deberá ser
1 asumiendo que el error en estado estacionario es 0. Si el valor pico de la respuesta al escalón es mayor
al valor …nal, esta diferencia se denomina sobreimpulso. Frecuentemente, el sobreimpulso se expresa en
términos de porcentajes mediante la siguiente fórmula
Respuesta pico Valor …nal
Porcentaje de sobreimpulso = 100 %:
Valor …nal
El tiempo que toma a la respuesta alcanzar el valor pico se denomina tiempo pico. Existen otros dos
tipos de medidas para determinar la rapidez de la respuesta del sistema, la primera es el tiempo de subida

Tr = t90 % t10 % ;

donde t90 % es el tiempo en que la respuesta alcanza el 90 % de su valor …nal y t10 % es el tiempo en que
la respuesta alcanza el 10 % de su valor …nal. La otra medida es el tiempo de establecimiento (Ts ), este se
de…ne como el tiempo que le toma a la respuesta establecerse entre los límites % del valor …nal. Los
valores mínimo y máximo de este intervalo se de…nen como los límites de tolerancia dentro de los cuales la
respuesta al escalón permanecerá para cualquier tiempo mayor al tiempo de establecimiento. Los valores
de generalmente se encuentran en el rango de 1 % a 5 %.
 17

1.4
Respuesta pico (Peak response)

Sobreimpulso >1 V (Overshoot > 1 V)

1.2 Límites ±5%
(Limits ±5%)

Valor
Respuesta al escalón (Step final - Estado estacionario
Response)
(Final value - Steady state)
0.8

0.6 Tiempo de subida

(Rise Time)

0.4

0.2

Tiempo de establecimiento
Tiempo
(Settling (Time) (seconds)
time)
0
Tiempo pico (Peak time) Tiempo [segundos]

Figura 2.12: Características de la respuesta al escalón.

2.5. Análisis del error

La respuesta temporal de un sistema de control consiste en dos partes: la respuesta transitoria y la

estacionaria, tal como se ilustra en la Figura 2.13. La respuesta transitoria se de…ne como aquella que va
desde el estado inicial al estado …nal. La respuesta estacionaria o en estado estacionario se de…ne como la
forma en la que la salida del sistema se comporta cuando t ! 1.
El error en estado estacionario se de…ne como la diferencia entre la entrada y la salida de un sistema
cuando el límite en el tiempo tiende a in…nito. Este análisis solo es útil para sistemas estables, por este
motivo se debe determinar antes la estabilidad del sistema.
 18

Vs (t) [V ]

T ransitorio Estado estacionario

t [s]

Figura 2.13: Características de la respuesta al escalón: Transitorio y estado estacionario.

El error en estado estacionario se puede calcular para sistemas en lazo abierto o lazo cerrado con
retroalimentación unitaria sin importar el tipo de sistema o su entrada, por ejemplo, para un sistema en
lazo abierto de…nido por el siguiente diagrama

R(s) G(s) C(s)

Figura 2.14: Lazo abierto.

la salida está dada por

C (s) = R (s) G (s) ;

el error es la entrada menos la salida, es decir

E (s) = R (s) C (s) ;

y al realizar las sustituciones correspondientes se obtiene que el error del sistema en lazo abierto está dado
por

E (s) = R (s) [1 G (s)] :

Entonces, al suponer que el sistema es estable, el error en estado estacionario se calcula mediante el
 19

siguiente límite
e (t ! 1) = l m sE (s) = l m sR (s) [1 G (s)] :
s!0 s!0

Ahora, para un sistema en lazo cerrado de…nido por el siguiente diagrama

E(s)
R(s) + G(s) C(s)
!

Figura 2.15: Lazo cerrado.

se determina su función de transferencia mediante álgebra de bloques y se obtiene

C (s) G (s)
= ;
R (s) 1 + G (s)
esta fórmula se puede deducir al realizar las siguientes igualdades

E (s) = R (s) (1) C (s) ;

C (s) = E (s) G (s) ;

entonces al sustituir E (s) en C (s), se llega a lo siguiente

C (s) = [R (s) C (s)] G (s) ;

C (s) [1 + G (s)] = G (s) R (s) ;

por lo tanto, se comprueba la función de transferencia obtenida por el algebra de bloques. Debido a que
el error es la entrada menos la salida, se tiene la siguiente fórmula

E (s) = R (s) C (s) ;

por lo tanto, al realizar las sustituciones correspondientes se obtiene lo siguiente

G (s) R (s) G (s)

E (s) = R (s) = 1 R (s) ;
1 + G (s) 1 + G (s)
…nalmente, el error del sistema de control en lazo cerrado está dado por

R (s)
E (s) = :
1 + G (s)
 20

Entonces, al suponer que el sistema es estable, es decir, todos los polos tienen parte real negativa, el error
en estado estacionario se calcula mediante el siguiente límite

sR (s)
e (t ! 1) = l m sE (s) = l m :
s!0 s!0 1 + G (s)

2.6. El concepto de estabilidad

Un sistema dinámico estable es aquel que presenta una respuesta acotada a una entrada o perturbación
acotada. Las respuestas de un sistema se clasi…can de la siguiente manera:

1. Estable y sobreamortiguado.

2. Estable y subamortiguado.

3. Marginalmente estable, oscilatorio o no amortiguado.

4. Inestable.

La estabilidad de un sistema se determina al calcular los polos. Recordando que los polos de un sistema
son las raíces del denominador de su función de transferencia

Y (s) b0 sm + b1 sm 1 + ::: + bm 1 s + bm
= ;
X (s) a0 sn + a1 sn 1 + ::: + an 1 s + an
es decir, su ecuación característica X (s) = 0 y los ceros son las raíces del numerador Y (s) = 0. Una
condición necesaria y su…ciente para que un sistema en lazo abierto o con retroalimentación en lazo cerrado
sea estable es que todos los polos de la función de transferencia tengan parte real negativa, por lo tanto, el
sistema es inestable si alguno de ellos tiene parte real positiva y la respuesta será desacotada para cualquier
tipo de entrada. Si la ecuación característica tiene raíces simples en el eje imaginario y todas las demás
raíces en el semiplano izquierdo, entonces la salida en estado estacionario serán oscilaciones sostenidas
para una entrada acotada. Sin embargo, si la entrada es una función sinusoidal con una frecuencia igual
a la magnitud de las raíces en el eje imaginario, entonces la salida será descotada y este tipo de sistemas
se de…nen como marginalmente estables.

2.7. Criterios de estabilidad

Una forma de ilustrar las raíces en el plano es mediante el denominado lugar de las raíces (del inglés
Root Locus), es decir, la ubicación geométrica de los polos y los ceros en el plano complejo s de la función
de transferencia del sistema como se ilustra un ejemplo en la Figura 2.16.
 21

Root Locus
3
0.89 0.81 0.7 0.56 0.38 0.2
Imaginary Axis (seconds-1 )

2
0.95
1 0.988
5 4 3 2 1
0
0.988
-1
0.95
-2
0.89 0.81 0.7 0.56 0.38 0.2
-3
-5 -4 -3 -2 -1 0 1 2
Real Axis (seconds -1 )
Figura 2.16: Ejemplo de grá…ca del lugar de las raíces en el plano complejo, donde el símbolo ‘ ’indica
la ubicación de los ceros y símbolo ‘ ’la de los polos del sistema. El eje x pertenece a los números reales
y el eje y a los números imaginarios.

El método de estabilidad de Routh–Hurwitz proporciona una respuesta al concepto de estabilidad al

considerar la ecuación característica de un sistema, la cual al ser representada en su forma Laplaciana se
escribe como se muestra a continuación

X (s) = an sn + an 1s
n 1
+ + a1 s + a0 = 0;

y como se mencionó anteriormente, para determinar la estabilidad del sistema es necesario calcular las
raíces de la ecuación característica. El criterio Routh–Hurwitz es su…ciente y necesario para establecer
estabilidad en sistemas lineales. Para el caso particular de sistemas de segundo orden, el requisito que se
debe de cumplir es que todos los coe…cientes del polinomio sean del mismo signo, es decir, todos positivos
o todos negativos. Por ejemplo, dado el sistema

X (s) = a2 s2 + a1 s + a0 ;

el arreglo de Routh está dado por

sn an an 2 s2 a2 a0
sn 1 an 1 an 3 = s1 a1 0 ;
sn 2 b1 b2 s0 b1 0

donde
an 1 an 2 an an 3
b1 = = an 2 = a0 ;
an 1

de acuerdo al criterio de Routh-Hurtwitz, el número de raíces de X (s) con parte real positiva es igual
al número de cambios de signo en la primera columna del arreglo de Routh. Con base en lo anterior se
 22

establece lo siguiente: Si todos los coe…cientes del polinomio característico de un sistema de segundo orden
son positivos o negativos, entonces todas sus raíces son negativas y el sistema es estable.
Para el caso de un sistema con una ecuación característica de tercer orden, es decir,

X (s) = a3 s3 + a2 s2 + a1 s + a0 ;

el arreglo de Routh está dado por

sn an an 2 s2 a3 a1
sn 1 an 1 an 3 s2 a2 a0
= ;
sn 2 b1 b2 s1 a1 0
sn 3 c1 c2 s0 b1 0

donde

an 1 an 2i an an (2i+1)
bi = ;
an 1
b1 an (2i+1) an 1 bi+1
ci = ;
b1

por lo tanto

an 1 an 2i an an (2i+1) a2 a1 a3 a0
b1 = = ;
an 1 a2
b2 = 0;

b1 a0 a2 b2
c1 = = a0 ;
b1
c2 = 0;

entonces, para que un sistema de tercer orden sea estable, es necesario y su…ciente que los coe…cientes
sean positivos y que a2 a1 a3 a0 > 0, sin embargo, si se da el caso en que a2 a1 = a3 a0 , entonces el sistema
es marginalmente estable dado que un par de raíces se ubican en el eje imaginario del plano–s.
 23

Capítulo 3

Controladores

En este capítulo se muestra el análisis para obtener las funciones de transferencia correspondientes
a las diferentes combinaciones de los controladores P, I y D, así como la respuesta del restador y del
sumador.

3.1. Controlador PID

El controlador Proporcional–Integral–Derivativo (PID) calcula continuamente el valor del error entre

la salida y la entrada de un sistema dinámico. La con…guración del controlador PID se muestra en la
Figura 3.1.

Figura 3.1: Diagrama eléctrico del controlador Proporcional-Integral-Derivativo inversor.

Los tres elementos del controlador PID se describen a continuación:

P: Proporcional (kP ) al error en el instante t, representa el error presente. La ganancia proporcional

kP tiene el efecto de aumentar proporcionalmente la señal de control para el mismo nivel de error.
El controlador proporcional hace que el sistema en lazo cerrado reaccione más rápidamente, pero
también aumenta el sobreimpulso. Además, tiende a reducir, pero no a eliminar, el error en estado
estacionario.

I: Proporcional (kI ) a la integral del error en el instante t, representa la acumulación del error pasado.
 24

La ganancia integral kI ayuda a reducir el error en estado estacionario. Si hay un error persistente
y constante, el integrador va aumentando la señal de control y reduce el error. Sin embargo, el
inconveniente radica en que puede ocasionar que la respuesta sea más lenta y en algunos casos
oscilatoria debido a los cambios de signo en la señal del error.

D: Proporcional (kD ) a la derivada del error en el instante t, representa una predicción el error futuro. La
ganancia derivativa kD puede ocasionar que la señal de control se vuelva grande inclusive cuando la
magnitud del error es relativamente pequeña. Aunque esta ganancia tiende a agregar amortiguación
a la respuesta del sistema, lo que disminuye el sobreimpulso, no tiene ningún efecto sobre el error
en estado estacionario.

Por lo tanto, el controlador PID se entiende como un controlador que considera el presente, pasado y
futuro del error estimado.
Antes de calcular la función de transferencia del controlador PID, es necesario indicar que una
característica fundamental de un ampli…cador operacional ideal es que tiene una impedancia de entrada
in…nita y por la tanto una corriente de entrada nula, debido a esto la diferencia de potencial entre las
terminales v+ y v es cero, y el nodo V0 (t) se denomina ‘tierra virtual’. Esto permite establecer lo
siguiente:

V0 (t) = I0 (t) = 0;

Ie (t) = Ir (t) ;

entonces al analizar el circuito se tiene

Ve (t) V0 (t) d
Ie (t) = +[Ve (t) V0 (t)] Ce = Ir (t) ;
Re dt
Z
1
Vs (t) = Rr Ir (t) + Ir (t) dt ;
Cr
y al aplicar la transformada de Laplace con condiciones iniciales iguales a cero se obtiene
Ve (s)
Ir (s) = + Ve (s) Ce s;
Re
Ir (s)
Vs (s) = Rr Ir (s) + ;
Cr s
con lo cual se puede sustituir Ir (s) en Vs (s)
Rr Cr s + 1 1 + Re Ce s
Vs (s) = Ve (s) ;
Cr s Re
por lo tanto, la función de transferencia del controlador PID está dada por
Vs (s) Rr Cr s + 1 1 + Re Ce s Re Rr Ce Cr s2 + (Re Ce + Rr Cr ) s + 1
= = :
Ve (s) Cr s Re Re Cr s
 25

El signo negativo en la función de transferencia se debe a la con…guración del controlador con el

ampli…cador operacional, por lo cual es necesario compensar este signo negativo con el restador o un
inversor 1 . Al separar y simpli…car los términos del numerador con el denominador se obtienen la ganancia
proporcional (kP ) y los coe…cientes de integración (kI ) y derivación (kD ) como se muestra a continuación.

(Re Ce + Rr Cr ) s Ce Rr Ce
P : = + = + kP ;
Re Cr s Cr Re Cr
1 kI
I : = ;
Re Cr s s
Re Rr Ce Cr s2
D : = Rr Ce s = kD s:
Re Cr s

Para eliminar el impacto que pueda ocasionar la relación Ce =Cr de los capacitores se puede construir
el controlador PID como se ilustra en el diagrama de bloques de la Figura 3.2.

kP

R +
e(t) kI + u(t)
+

kD
d
dt

Figura 3.2: Diagrama del controlador Proporcional-Integral-Derivativo en paralelo.

La respuesta u (t) de esta con…guración del controlador PID es una suma ponderada de la señal e (t)
que representa el error entre la señal de referencia y la salida del sistema, la integral de e (t) y la derivada
e (t) como se indica a continuación
Z t
de (t)
u (t) = kP e (t) + kI e (t) dt + kD ;
0 dt

es decir
kD s2 + kP s + kI
u (t) = :
s
En las siguientes secciones se muestra solamente el diagrama eléctrico de las otras con…guraciones y
el cálculo de su respectiva función de transferencia.
1
Para realizar las simulaciones se descarta este signo negativo.
 26

3.2. Controlador P

En la Figura 3.3 se ilustra el diagrama eléctrico del controlador P.

Figura 3.3: Diagrama eléctrico del controlador Proporcional inversor.

A continuación, se presenta el análisis para obtener la función de transferencia del controlador P.

ie (t) = i0 (t) + ir (t) ;

i0 (t) = V0 (t) = 0;

ie (t) = ir (t) ;

Ve (t) V0 (t) Ve (t)

ie (t) = = ;
Re Re

V0 (t) Vs (t) Vs (t)

ir (t) = = :
Rr Rr

Función de transferencia:
Vs (t) Rr
= = kP :
Ve (t) Re
 27

3.3. Controlador I

En la siguiente Figura 3.4 se ilustra el diagrama eléctrico del controlador I.

Figura 3.4: Diagrama eléctrico del controlador Integral inversor.

A continuación, se presenta el análisis para obtener la función de transferencia del controlador I.

ie (t) = i0 (t) + ir (t)

i0 (t) = V0 (t) = 0

ie (t) = ir (t)

Ie (s) = Ir (s) ;

Ve (t) = ie (t) Re ! Ve (s) = Ie (s) Re ;

Z
1 Ir (s)
Vs (t) = ir (t) dt ! Vs (s) = :
Cr Cr s

Función de transferencia:
Vs (s) 1 kI
= = :
Ve (s) Re Cr s s
 28

3.4. Controlador D

En la siguiente Figura 3.5 se ilustra el diagrama eléctrico del controlador D.

Figura 3.5: Diagrama eléctrico del controlador Derivativo inversor.

A continuación, se presenta el análisis para obtener la función de transferencia del controlador D.

ie (t) = i0 (t) + ir (t) ;

i0 (t) = V0 (t) = 0;

ie (t) = ir (t) ;

Ie (s) = Ir (s) ;

d [Ve (t) V0 (t)] dVe (t) Ie (s)

ie (t) = Ce = Ce ! Ve (s) = ;
dt dt Ce s

V0 (t) Vs (t) Vs (t)

ir (t) = = ! Vs (s) = Ir (s) Rr :
Rr Rr

Función de transferencia:
Vs (s)
= Rr Ce s = kD s:
Ve (s)
 29

3.5. Controlador PI

En la siguiente Figura 3.6 se ilustra el diagrama eléctrico del controlador PI.

Figura 3.6: Diagrama eléctrico del controlador Proporcional-Integral inversor.

A continuación, se presenta el análisis para obtener la función de transferencia del controlador PI.

ie (t) = ir (t) ;

i0 (t) = V0 (t) = 0;

Ie (s) = Ir (s) ;

Ve (t) = ie (t) Re ! Ve (s) = Ie (s) Re ;

Z
1 Cr Rr s + 1
Vs (t) = ir (t) Rr ir (t) dt ! Vs (s) = Ir (s) :
Cr Cr s

Función de transferencia:

Vs (s) Rr Cr s + 1 Rr 1 kI
= = + = kP + :
Ve (s) Re Cr s Re Re Cr s s
 30

3.6. Controlador PD

En la siguiente Figura 3.7 se ilustra el diagrama eléctrico del controlador PD.

Figura 3.7: Diagrama eléctrico del controlador Proporcional-Derivativo inversor.

A continuación, se presenta el análisis para obtener la función de transferencia del controlador PD.

ie (t) = i0 (t) + ir (t) ;

i0 (t) = V0 (t) = 0;

ie (t) = ir (t) ;

Ie (s) = Ir (s) ;

Ve (t) V0 (t) d [Ve (t) V0 (t)] Ve (t) dVe (t) Re Ie (s)

ie (t) = + Ce = + Ce ! Ve (s) = ;
Re dt Re dt Re Ce s + 1

V0 (t) Vs (t) Vs (t)

ir (t) = = ! Vs (t) = Rr Ir (s) :
Rr Rr

Función de transferencia:

Vs (s) Ce Re Rr s + Rr Rr
= = + Rr Ce s = (kP + kD s) :
Ve (s) Re Re
 31

3.7. Restador

Con esta con…guración del ampli…cador se utilizan las entradas invertida y no invertida con una
ganancia de 1, para producir una salida igual a la diferencia entre las entradas. En casos particulares se
pueden elegir también otros valores de resistencias para ampli…car la salida. La con…guración se muestra
en la Figura 3.8.

Figura 3.8: Diagrama eléctrico del Restador.

A continuación, se presenta el análisis del circuito del restador inversor. Las corrientes están dadas
por

I0 (t) = 0;

I1 (t) = I3 (t) + I0 (t) = I3 (t) ;

I2 (t) = I4 (t) + I0 (t) = I4 (t) ;

y expresadas en función de los voltajes se tiene lo siguiente

V1 (t) V0 (t)
I1 (t) = ;
R
V2 (t) V0 (t)
I2 (t) = ;
R
V0 (t) Vs (t)
I3 (t) = ;
R
V0 (t)
I4 (t) = ;
R
por lo tanto, se tienen las siguientes igualdades
V1 (t) V0 (t) V0 (t) Vs (t)
= ;
R R

V2 (t) V0 (t) V0 (t)

= ;
R R
 32

despejando V0 (t) en ambas ecuaciones se llega, respectivamente a

2 V1 (t) + Vs (t) V1 (t) + Vs (t)

V0 (t) = ! V0 (t) = ;
R R 2

2 V2 (t) V2 (t)
V0 (t) = ! V0 = ;
R R 2

con base en lo anterior se escribe lo siguiente

V1 (t) + Vs (t) V2 (t)

= ;
2 2
entonces, la función de transferencia está dada por

Vs (t) = V2 (t) V1 (t) ;

y si se utiliza como restador inversor, la respuesta es la siguiente:

Vs (t) = [V1 (t) V2 (t)] :

 33

3.8. Sumador

En la Figura 3.9 se ilustra el diagrama eléctrico del sumador inversor.

Figura 3.9: Diagrama eléctrico del sumador inversor.

A continuación, se presenta el análisis para obtener la función de transferencia del sumador.

i1 (t) + i2 (t) = i3 (t) + i0 (t) ;

V1 (t) V0 (t)
i1 (t) = ;
R
V2 (t) V0 (t)
i2 (t) = ;
R
V0 Vs (t)
i3 (t) = ;
R

i0 (t) = V0 (t) = 0;
V1 (t) V2 (t) Vs (t)
+ = ;
R R R
Función de transferencia:
Vs (t) = [V1 (t) + V2 (t)] :
 34

Capítulo 4

Modelado matemático

El circuito RLC de la Figura 4.1 representa un sistema de segundo orden que modeliza, de manera
simpli…cada, la mecánica pulmonar.

L R Pao (t)
Pao (t) PA (t) Q(t)
L
Q(t) R
C

P0 C PA (t)

Figura 4.1: Modelo eléctrico de la mecánica pulmonar. Sistema pulmonar elaborado con [Link].

El resistor R representa una combinación de la resistencia al ‡ujo de aire [Q (t)] en las vías respiratorias,
el tejido pulmonar y la pared torácica. El inductor L representa la inertancia1 al ‡ujo de aire en las vías
respiratorias. El capacitor C representa la compliancia2 combinada de las vías respiratorias, el tejido
pulmonar y la pared torácica. Los componentes R y C representan respectivamente las propiedades
mecánicas resistivas y de almacenamiento del sistema respiratorio. Para las simulaciones numéricas se
deben utilizar los siguientes valores para cada elemento del circuito RLC.

R = 10 k ;

L = 1 H;

C = 220 F:

El objetivo del modelo eléctrico es predecir la respuesta dinámica de la presión alveolar [PA (t)] a
diferentes formas de onda de presión [Pao (t)] aplicadas en la apertura de las vías respiratorias, es decir
PA (t)
:
Pao (t)
1
La inertancia (del inglés inertance) es la magnitud necesaria en la diferencia de presión requerida en un ‡uido para causar
un cambio en la tasa de ‡ujo con respecto al tiempo.
2
La compliancia (del inglés compliance) pulmonar es una medida de la habilidad que tiene un pulmón para estirarse y
expandirse.
 35

Referencia

Khoo MC. Physiological Control Systems Analysis, Simulation, and Estimation. Wiley, 2018. Section 4, page 93.

4.1. Función de transferencia, error y estabilidad del sistema de

mecánica pulmonar

Ecuaciones que describen la dinámica del sistema al analizar por mallas

Z
dQ (t) 1
Pao (t) = L + Q (t) R + Q (t) dt;
dt C
Z
1
PA (t) = Q (t) dt;
C

ahora, al aplicar la transformada de Laplace se obtiene lo siguiente

Q (s)
Pao (s) = LsQ (s) + Q (s) R + ;
Cs

al realizar la suma algebraica de cada término y factorizar se obtiene

LCs2 + RCs + 1
Pao (s) = Q (s) ;
Cs

la salida del sistema se de…ne de la siguiente manera

Q (s)
PA (s) = ;
Cs

por lo tanto, la función de transferencia está dado por

Q (s)
PA (s) Cs
= ;
Pao (s) LCs2 + RCs + 1
Q (s)
Cs
y al simpli…car la expresión se llega al siguiente resultado

PA (s) 1
= :
Pao (s) LCs2 + RCs + 1

El error [E (s)] en estado estacionario se calcula aplicando el siguiente límite

PA (s) 1 1
E (s) = l m sPao (s) 1 =lm s 1 = 0;
s!0 Pao (s) s!0 s LCs2 + RCs + 1

este resultado implica que para el escenario donde se aplica el escalón unitario como entrada al sistema
se cumple lo siguiente
l m PA (t) = Pao (t) = 1 V;
t!1

sin embargo, no se obtiene información sobre el tiempo de establecimiento, lo cual incide sobre el error
cuando se aplican las otras entradas variantes en el tiempo como el impulso, la rampa y la función
 36

sinusoidal. El objetivo del controlador PID (o sus variantes) es eliminar el error ante cualquier señal de
entrada al sistema en un intervalo de tiempo adecuado para el problema que se pretende resolver.
En cuanto a la estabilidad del sistema en lazo abierto, de acuerdo con el criterio de Routh-Hurtwitz,
el sistema de mecánica pulmonar es estable dado que todos los coe…cientes del denominador en la función
de transferencia son del mismo signo, y por lo tanto, la parte real de los ceros es negativa (Re < 0), esto
implica que la salida del sistema será acotada para una señal de entrada (o perturbación) acotada.

4.2. Modelo matemático de ecuaciones integro-diferenciales

El modelo matemático de ecuaciones integro-diferenciales se formula despejando las variables

dependientes (corrientes o voltajes dependiendo del problema que se analiza) en las ecuaciones principales
del sistema, en este caso las variables están dadas por el ‡ujo de aire en las vías respiratorias [Q (t)], y la
presión alveolar [PA (t)], por lo tanto, al realizar el despeje se obtiene lo siguiente:
Z
dQ (t) 1 1
Q (t) = Pao (t) L Q (t) dt ;
dt C R

Z
1
PA (t) = Q (t) dt:
C

Cabe destacar que la variable Q (t) solamente se puede despejar de los términos donde no se esté
derivando o integrando. Se observa también que el modelo es lineal y de primer orden, es decir, la potencia
más alta en las derivadas es de primer grado, y los términos Pao (t), Q (t) y PA (t) dependen solamente
de la variable independiente t.

4.3. Respuesta del sistema en lazo abierto

En esta sección se presenta la respuesta del sistema de mecánica pulmonar a las señales de escalón,
impulso, rampa y sinusoidal en lazo abierto.

4.3.1. Multisim: lazo abierto

Se debe construir el circuito eléctrico en Multisim como se muestra en la Figura 4.2 con los siguientes
elementos.

1. Un inductor (Inductor ) de 1 H, del menú Place basic.

2. Un resistor (Resistor ) de 10 k , del menú Place basic.

3. Un capacitor (Capacitor ) de 220 F , del menú Place basic.

 37

Figura 4.2: Sistema de mecánica pulmonar en Multisim.

4. Un delay (Delay) de 1 s, del menú Place Source:

5. Un osciloscopio de cuatro canales (4 Channel Oscilloscope-XSC ).

El osciloscopio se debe con…gurar a 1 s=Div y 500 mV =Div, la simulación se debe realizar en un

intervalo de tiempo de 0 t 10 segundos con un período de muestreo de 1 ms.

6. Tres interruptores (Switch), del menú Place basic.

7. Tres tierras (Ground ), del menú Place Source.

8. Cuatro generadores de funciones (Function Generator-XFG).

Para obtener cada una de las señales se debe con…gurar el generador de funciones como se indica a
continuación.

1. Escalón: Seleccionar la forma de onda (Waveform) cuadrada y con…gurar el dispositivo como se

indica a continuación (Signal options):

F requency : 100 mHz

Duty cycle : 99 %
Amplitude : 500 mVp
Of f set : 500 mV
 38

Nota: Set rise / fall time: 10 ms:

2. Impulso: Seleccionar la forma de onda (Waveform) cuadrada y con…gurar el dispositivo como se

indica a continuación (Signal options):

F requency : 100 mHz

Duty cycle : 10 %
Amplitude : 500 mVp
Of f set : 500 mV

Nota: Set rise / fall time: 10 ms:

3. Rampa: Seleccionar la forma de onda (Waveform) dientes de sierra y con…gurar el dispositivo como
se indica a continuación (Signal options):

F requency : 100 mHz

Duty cycle : 99 %
Amplitude : 500 mVp
Of f set : 500 mV

4. Sinusoidal: Seleccionar la forma de onda (Waveform) sinusoidal y con…gurar el dispositivo como se

indica a continuación (Signal options):

F requency : 250 mHz

Duty cycle : 50 %
Amplitude : 1 Vp
Of f set : 0 V

Para realizar las simulaciones es necesario con…gurar los ajustes, para esto se debe ir al menú:
Simulate ! Interactive simulation settings, y se debe con…gurar lo indicado en la Figura 4.3.
Los parámetros indicados para ajustar son los siguientes.

1. Initial conditions: Set to zero.

2. Start time (TSTART): 0 s:

3. End time (TSTOP): 10 s:

4. Maximum time step (TMAX): 0:001 s = 1 ms = 1E 3:

5. Initial time step (TSTEP): 0:001 s = 1 ms = 1E 3:

 39

Figura 4.3: Ajustes de simulación.

Es importante mencionar que en algunas ocasiones es necesario modi…car la con…guración de los ajustes
para cada simulación. Ahora, se muestran las respuestas del sistema de mecánica pulmonar en Multisim
a cada una de las entradas descritas anteriormente, es decir, al escalón en la Figura 4.4, al impulso en la
Figura 4.5, a la rampa en la Figura 4.6, y a la función sinusoidal en la Figura 4.7.
 40

Figura 4.4: Respuesta del sistema de mecánica pulmonar al escalón unitario en Multisim.
 41

Figura 4.5: Respuesta del sistema de mecánica pulmonar al impulso unitario en Multisim.
 42

Figura 4.6: Respuesta del sistema de mecánica pulmonar a la rampa en Multisim.

 43

Figura 4.7: Respuesta del sistema de mecánica pulmonar a la función sinusoidal en Multisim.
 44

4.3.2. Laboratorio de Electrónica: lazo abierto

Comparar las respuestas de la emulación en Multisim del circuito RLC de la mecánica pulmonar
simpli…cada con los resultados obtenidos en el Laboratorio de Electrónica.
Para el generador de funciones RIGOL DG1022 utilice las siguientes con…guraciones para generar
cada una de las cuatro señales utilizadas en Multisim:

1. Escalón: Seleccionar la forma de onda pulso (pulse) y con…gurar el dispositivo como se indica a
continuación:

Señal: P ulse
F requency : 100 mHz
Amplitude : 1 Vpp
Of f set : 500 mVDC
Duty cycle : 99 %

2. Impulso: Seleccionar la forma de onda pulso (pulse) y con…gurar el dispositivo como se indica a
continuación:

Señal: P ulse
F requency : 100 mHz
Amplitude : 1 Vpp
Of f set : 500 mVDC
Duty cycle : 10 %

3. Rampa: Seleccionar la forma de onda rampa (ramp) y con…gurar el dispositivo como se indica a
continuación:

Señal: Ramp
F requency : 100 mHz
Amplitude : 1 Vpp
Of f set : 500 mVDC
Symmetry : 99 %

4. Sinusoidal: Seleccionar la forma de onda sinusoidal (sine) y con…gurar el dispositivo como se indica
a continuación:

Señal: Sine
F requency : 250 mHz
Amplitude : 2 Vpp
Of f set : 0 VDC
 45

4.4. Diseño del controlador

En esta sección se presenta el diseño de un controlador para el sistema de mecánica pulmonar, para
esto se utiliza la herramienta Tune de Simulink para la sintonización de las ganancias en un controlador
PID.

4.4.1. Simulink: lazo abierto

El lienzo principal que se debe construir en Simulink para obtener la respuesta en lazo abierto se
muestra en la Figura 4.8.

Figura 4.8: Diagrama de bloques para obtener la respuesta del sistema de mecánica pulmonar en lazo
abierto en Simulink.

El conjunto de bloques que se necesitan para construir el diagrama y su con…guración se indican a

continuación.

1. El bloque Step para generar el escalón unitario con los siguientes ajustes en sus parámetros:

a) Step time: 0

b) Initial value: 0

c) Final value: 1

d ) Sample time: 0:001

 46

2. El bloque Pulse Generator para generar el impulso unitario con los siguientes ajustes en sus
parámetros:

a) Amplitude: 1

b) Period (secs): 10

c) Pulse Width ( % of period): 10

d ) Phase delay (secs): 1

3. El bloque Ramp para generar la rampa con los siguientes ajustes en sus parámetros:

a) Slope: 1=10

b) Start time: 0

c) Initial output: 0

4. El bloque Sine Wave para generar la señal sinusoidal con los siguientes ajustes en sus parámetros:

a) Sine type: Time based

b) Time (t): Use simulation time

c) Amplitude: 1

d ) Bias: 0

e) Frequency (rad= sec): 1:5708

f ) Phase (rad): 0

g) Sample time: 0:001

5. El bloque Subsystem se utiliza para diseñar el modelo de ecuaciones integro-diferenciales.

6. El bloque Scope para visualizar la respuesta del sistema mediante series en el tiempo en el
osciloscopio.

7. El bloque Mux para multiplexar las señales y suministrarlas al osciloscopio.

8. El bloque Manual Switch para cambiar entre cada una de las cuatro señales de entrada.
 47

En la Figura 4.9 se muestra el diagrama de bloques correspondiente al sistema de mecánica pulmonar

para obtener la respuesta en lazo abierto al escalón, al impulso, a la rampa y a la función sinusoidal. Para
esto se necesita el modelo matemático de ecuaciones integro-diferenciales que describe al sistema, es decir,
Z
dQ (t) 1 1
Q (t) = Pao (t) L Q (t) dt ;
dt C R

Z
1
PA (t) = Q (t) dt:
C

Figura 4.9: Diagrama de bloques para representar el modelo matemático de ecuaciones integro-
diferenciales.

El conjunto de bloques que se necesitan para construir el modelo matemático de ecuaciones integro-
diferenciales se indican a continuación.

1. El bloque Gain se utiliza para representar los coe…cientes, es decir, los valores de la resistencia (R),
la inductancia (L) y la capacitancia (C).
Z
2. El bloque Integrator para integrar la corriente Q (t).

d
3. El bloque Derivative para derivar la corriente Q (t).
dt
4. El bloque Add para realizar la suma algebraica de los términos en el modelo de ecuaciones integro-
diferenciales.

La con…guración de los parámetros para realizar la simulación se presenta en la Figura 4.10.

 48

Figura 4.10: Con…guración de los parámetros para la simulación.

Los parámetros indicados para ajustar son los siguientes.

1. Start time: 0:0 s:

2. Stop time: 10:0 s:

3. Type: Variable-step.

4. Solver: ode45, ode23, ode15s u ode23s (según lo requiera el problema, sin embargo, se recomienda
iniciar con el ode45 ).

5. Max step size: 1e 3:

6. Relative tolerance: 1e 3:

La respuesta obtenida del sistema de mecánica pulmonar al escalón se ilustra en la Figura 4.11, al
impulso en la Figura 4.12, a la rampa en la Figura 4.13, y a la función sinusoidal en la Figura 4.14.
 49

Pao(t)
PA(t)

0.8

0.6
V(t) [volts]

0.4

0.2

0 1 2 3 4 5 6 7 8 9 10
Time

Figura 4.11: Respuesta del sistema de mecánica pulmonar al escalón unitario en Simulink.
 50

Pao(t)
PA(t)

0.8

0.6
V(t) [volts]

0.4

0.2

0 1 2 3 4 5 6 7 8 9 10
Time

Figura 4.12: Respuesta del sistema de mecánica pulmonar al impulso unitario en Simulink.
 51

Pao(t)
PA(t)

0.8

0.6
V(t) [volts]

0.4

0.2

0 1 2 3 4 5 6 7 8 9 10
Time

Figura 4.13: Respuesta del sistema de mecánica pulmonar a la rampa en Simulink.

 52

Pao(t)
PA(t)

0.5
V(t) [volts]

-0.5

-1

0 1 2 3 4 5 6 7 8 9 10
Time

Figura 4.14: Respuesta del sistema de mecánica pulmonar a la función sinusoidal en Simulink.
 53

4.4.2. Simulink: lazo cerrado

Primero, se debe complementar el lienzo mostrado en la Figura 4.8 con el sistema de control en lazo
cerrado como se indica en la Figura 4.15.

Figura 4.15: Diagrama de bloques en Simulink para diseñar el controlador y obtener la respuesta en lazo
cerrado del sistema de mecánica pulmonar.

Los nuevos bloques utilizados en el diagrama son los siguientes.

1. El bloque Sum para calcular el error entre la señal de entrada y la respuesta de sistema.

2. El bloque PID Controller para diseñar y aplicar el controlador al sistema de mecánica pulmonar.

El bloque del controlador PID implementa un controlador PID (PID, PI, PD, solo P o solo I) donde
la salida del bloque es una suma ponderada del error (la diferencia entre la señal de entrada y la señal de
salida), la integral del error y la derivada del error como se indica a continuación
Z t
de (t)
u (t) = kP e (t) + kI e (t) dt + kD ;
0 dt

donde e (t) es el error y los coe…cientes kP , kI y kD son las ganancias del controlador. El bloque admite
varios tipos y estructuras de control. Las principales opciones con…gurables en el bloque se indican a
continuación.

1. Tipo de controlador (Controller: PID, PI, PD, solo P o solo I).

2. Dominio de tiempo (Time domain: Continuous-time).

3. Condiciones iniciales (Initial conditions: internal ).

 54

A medida que cambian estas opciones, la estructura interna del bloque también lo hace al activar
diferentes subsistemas dentro de él. Para realizar la sintonización de las ganancias se deben con…gurar los
parámetros del bloque como se indica en la Figura 4.16.

Figura 4.16: Ventana de ajustes para el bloque del controlador PID.

Se observa que todas las ganancias se ajustan a 1, este valor se toma arbitrariamente debido a que es
necesario que estos campos contengan un número real para realizar la sintonización de las ganancias con
la herramienta Tune.
 55

El diseño del controlador tiene como objetivo una sintonización de las ganancias kP , kI
y kD que permitan obtener una respuesta del sistema en lazo cerrado con un tiempo de
establecimiento Ts 100 ms y un porcentaje de sobreimpulso 10 %, es decir, un voltaje o
respuesta pico 1:1 V .

Para alcanzar el objetivo anterior con los valores de resistencia, inductancia y capacitancia utilizados
en el sistema de mecánica pulmonar, se debe ajustar el tiempo de la respuesta [Response Time (seconds)]
a 0:014 s y el comportamiento transitorio (Transient Behavior ) a 0:9 como se muestra en la Figura 4.17.
El parámetro del comportamiento transitorio puede tomar valores entre 0 y 0:9, un comportamiento
más ‘agresivo’o más ‘robusto’respectivamente. Un comportamiento transitorio más robusto disminuye el
sobreimpulso en la respuesta, mientras que un comportamiento transitorio más agresivo tiende a producir
oscilaciones sostenidas en la respuesta.
Los valores sintonizados para las ganancias son los siguientes (Parámetros del controlador):

Controller Parameters
P [kP ] 300:0188
I [kI ] 7575:5772
D [kD ] 0:39159
N 520:9749

aquí es importante mencionar que si la ganancia derivativa kD < 0, se toma kD = 0, es decir,

la sintonización de las ganancias proporciona como resultado un controlador PI. Adicionalmente, el
rendimiento y robustez del controlador se indican a continuación:

Performance and Robustness

Rise time 0:0127 s
Settling time 0:0968 s
Overshoot 10 %
Peak 1:1 V
Closed-loop stability Stable

Ahora, al actualizar el bloque del controlador PID en el botón Update Block de la ventana de ajustes
se observa el resultado de la Figura 4.18.
La respuesta obtenida del sistema de mecánica pulmonar en lazo cerrado con el controlador PID al
escalón se ilustra en la Figura 4.19, al impulso en la Figura 4.20, a la rampa en la Figura 4.21, y a la
función sinusoidal en la Figura 4.22.
 56

Figura 4.17: Sintonización de las ganancias para el controlador PID y ajustes en el tiempo de respuesta
y comportamiento transitorio de la respuesta del sistema de mecánica pulmonar en lazo cerrado.
 57

Figura 4.18: Bloque del controlador con los valores ajustados para las ganancias kP , kI y kD .
 58

1.2 Pao(t)
PA(t)
PID

0.8
V(t) [volts]

0.6

0.4

0.2

0 1 2 3 4 5 6 7 8 9 10
Time

Figura 4.19: Respuesta del sistema de mecánica pulmonar en lazo cerrado con el controlador PID al
escalón unitario en Simulink.
 59

1.2 Pao(t)
PA(t)
PID

0.8

0.6
V(t) [volts]

0.4

0.2

-0.2

0 1 2 3 4 5 6 7 8 9 10
Time

Figura 4.20: Respuesta del sistema de mecánica pulmonar en lazo cerrado con el controlador PID al
impulso unitario en Simulink.
 60

Pao(t)
PA(t)
PID
1

0.8

0.6
V(t) [volts]

0.4

0.2

0 1 2 3 4 5 6 7 8 9 10
Time

Figura 4.21: Respuesta del sistema de mecánica pulmonar en lazo cerrado con el controlador PID a la
rampa en Simulink.
 61

Pao(t)
PA(t)
PID

0.5
V(t) [volts]

-0.5

-1

0 1 2 3 4 5 6 7 8 9 10
Time

Figura 4.22: Respuesta del sistema de mecánica pulmonar en lazo cerrado con el controlador PID a la
función sinusoidal en Simulink.
 62

4.4.3. Multisim: Lazo cerrado

Ahora, con base en las ganancias kP , kI y kD obtenidas con la herramienta Tune de Simulink, se
procede a determinar los valores de los componentes del controlador PID, cuya función de transferencia
está dada por la siguiente expresión

Vs (s) Re Rr Ce Cr s2 + (Re Ce + Rr Cr ) s + 1 Ce Rr 1 Ce kI
= = + + + Ce Rr s = + kP + + kD s;
Ve (s) Re Cr s Cr Re Re Cr s Cr s

al separar con respecto a los términos del numerador se obtiene el siguiente resultado

Vs (s) Ce Rr 1 Ce kI
= + + + Ce Rr s = + kP + + kD s;
Ve (s) Cr Re Re Cr s Cr s

y, a partir de lo anterior se identi…can los componentes para cada ganancia

Vs (s) Ce kI
= + kP + + kD s;
Ve (s) Cr s

por lo tanto

Rr
kP = = 300:0188;
Re

1
kI = = 7575:5772;
Re Cr

kD = Ce Rr = 0:39159:

Con base en estos resultados se propone un valor comercial para la capacitancia del capacitor Cr

6
Cr = 1 F = 1 10 F;

y se calcula el valor de la resistencia del resistor Re como se indica a continuación

1
Re = 6)
= 132 ;
(7575:5772) (1 10
entonces, el valor de resistencia del resistor Rr es de

Rr = (300:0188) (132) = 39602 ;

y el valor de la capacitancia del capacitor Ce deberá tener un valor de

0:39159 6
Ce = = 9: 888 1 10 ;
39602
 63

sin embargo, debido a que estos valores no son comerciales, se toman los más cercanos mostrados en
Multisim, es decir,

Re = 133 ;

Rr = 40 k ;

Ce = 10 F:

Los nuevos elementos utilizados se describen a continuación.

1. Dos ampli…cadores operacionales 741 para construir el restador y el controlador PID del menú Place
Analog.

2. Cuatro fuentes de voltaje (DC Power ), del menú Place Source, conectadas en serie y con…guradas
a 15 V .

3. Dos resistores para el controlador PID (Re y Rr ) y cuatro de 100 k para el restador.

4. Dos capacitores para el controlador PID (Ce y Cr ).

5. Un resistor (10 k ), un inductor (1 H) y un capacitor (220 F ) para comparar las respuestas en

lazo abierto y lazo cerrado.

6. Siete tierras.

El circuito del sistema de mecánica pulmonar, en lazo abierto y en lazo cerrado con el controlador
PID, se ilustra en la Figura 4.23, mientras que la respuesta obtenida del sistema en lazo cerrado al escalón
se ilustra en la Figura 4.24, al impulso en la Figura 4.25, a la rampa en la Figura 4.26, y a la función
sinusoidal en la Figura 4.27.
 64

Figura 4.23: Diagrama eléctrico en lazo cerrado del sistema de mecánica pulmonar con el controlador PID.
 65

Figura 4.24: Respuesta del sistema de mecánica pulmonar en lazo cerrado con el controlador PID al
escalón unitario en Multisim.
 66

Figura 4.25: Respuesta del sistema de mecánica pulmonar en lazo cerrado con el controlador PID al
impulso unitario en Multisim.
 67

Figura 4.26: Respuesta del sistema de mecánica pulmonar en lazo cerrado con el controlador PID a la
rampa en Multisim.
 68

Figura 4.27: Respuesta del sistema de mecánica pulmonar en lazo cerrado con el controlador PID a la
función sinusoidal en Multisim.
 69

4.5. Sistema de control en Python

En esta sección se muestra el código necesario para obtener la respuesta en lazo abierto y en lazo
cerrado con el controlador PID del sistema de mecánica pulmonar al escalón, al impulso, a la rampa y a
la función sinusoidal. Para esto se necesita la función de transferencia del sistema, la cual está dada por
la siguiente expresión
PA (s) 1
= = sys;
Pao (s) LCs2 + RCs + 1
donde

R = 10 k ;

L = 1 H;

C = 220 F;

y la función de transferencia del controlador PID

Vs (s) Re Rr Ce Cr s2 + (Re Ce + Rr Cr ) s + 1
= = P ID;
Ve (s) Re Cr s

donde

Re = 133 ;

Rr = 40 k ;

Ce = 10 F;

Cr = 1 F;

por lo tanto, la respuesta en lazo cerrado se obtiene mediante algebra de bloques como se indica a
continuación
P ID sys
sysP ID = :
1 + P ID sys
La respuesta obtenida del sistema de mecánica pulmonar en lazo abierto y en lazo cerrado con el
controlador PID se ilustra al escalón en la Figura 4.28, al impulso en la Figura 4.29, a la rampa en la
Figura 4.30, y a la función sinusoidal en la Figura 4.31.
 70

!pip install control #Modulo para implementar sistemas de control

!pip install slycot #Modulo con rutinas y solucionadores para sistemas de control

import numpy #Libreria especializada en el calculo numerico y el analisis de datos

import [Link] as plt #Libreria para la generacion de graficas
import control

x0,t0,tF,dt = 0,0,10,1E-3
N=round((tF-t0)/dt)+1 #Numero total de iteraciones
t=[Link](t0,tF,N) #Arreglo del tiempo de 0:dt:10 segundos
u1=[Link](N) #Escalon unitario
u2=[Link](N); u2[round(1/dt):round(2/dt)]=1 #Impulso
u3=([Link](t0,tF,N))/tF # Rampa
u4=[Link](1.5708*t) #Funcion sinosoidal, 1.5708 rad/s = 250 mHz

# Elementos del circuito RLC

R=10E3 #10 kOhms
L=1E-6 #1 uH
C=220E-6 #220uF
# Funcion de transferencia
num=[1]
den=[C*L,C*R,1]
sys=[Link](num,den)
print(sys)

# Sistema en lazo cerrado con controlador

Rr=40E3 #40 kOhms
Re=133 #133 Ohms
Cr=1E-6 #1 uF
Ce=10E-6 #10 uF
#Funcion de transferencia del controlador
numPID=[Rr*Re*Cr*Ce,Re*Ce+Rr*Cr,1]
denPID=[Re*Cr,0]
PID=[Link](numPID,denPID)

# Funcion de transferencia del sistema de control

X=[Link](sys,PID)
sysPID=[Link](X,1,sign=-1)
print(sysPID)
 71

# Respuesta al escalon unitario

fig1=[Link]()
[Link](t,u1, ’-’, color=[0.5,0.05,0.05], label=’$Pao(t)$’) #Entrada
ts,Vs=control.forced_response(sys,t,u1,x0)
[Link](t,Vs, ’-’, color=[0,0.25,0.4], label=’$PA(t)$’) #Salida (lazo abierto)
ts,pid=control.forced_response(sysPID,t,u1,x0)
[Link](t,pid, ’:’, linewidth=3, color=[0.3,0.5,0.2], label=’$PID$’) #Controlador (lazo cerrado)
[Link](True)
[Link](-0.5, 10)
[Link](-0.1, 1.2)
[Link](’$t$ $[segundos]$’)
[Link](’$V(t)$ $[volts]$’)
[Link](’Respuesta al escalon’)
[Link](loc=’lower right’)
[Link]()
fig1.set_size_inches(4,6)
[Link](’python_escalon.png’, dpi=600)
[Link](’python_escalon.pdf’)

# Respuesta al impulso unitario

fig2=[Link]()
[Link](t,u2, ’-’, color=[0.5,0.05,0.05], label=’$Pao(t)$’) #Entrada
ts,Vs=control.forced_response(sys,t,u2,x0)
[Link](t,Vs, ’-’, color=[0,0.25,0.4], label=’$PA(t)$’) #Salida (lazo abierto)
ts,pid=control.forced_response(sysPID,t,u2,x0)
[Link](t,pid, ’:’, linewidth=3, color=[0.3,0.5,0.2], label=’$PID$’) #Controlador (lazo cerrado)
[Link](True)
[Link](0, 10)
[Link](-0.2, 1.2)
[Link](’$t$ $[segundos]$’)
[Link](’$V(t)$ $[volts]$’)
[Link](’Respuesta al impulso’)
[Link](loc=’upper right’)
[Link]()
fig2.set_size_inches(4,6)
[Link](’python_impulso.png’, dpi=600)
[Link](’python_impulso.pdf’)
 72

# Respuesta a la rampa
fig3=[Link]()
[Link](t,u3, ’-’, color=[0.5,0.05,0.05], label=’$Pao(t)$’) #Entrada
ts,Vs=control.forced_response(sys,t,u3,x0)
[Link](t,Vs, ’-’, color=[0,0.25,0.4], label=’$PA(t)$’) #Salida (lazo abierto)
ts,pid=control.forced_response(sysPID,t,u3,x0)
[Link](t,pid, ’:’, linewidth=3, color=[0.3,0.5,0.2], label=’$PID$’) #Controlador (lazo cerrado)
[Link](True)
[Link](-0.5, 10)
[Link](-0.1, 1.2)
[Link](’$t$ $[segundos]$’)
[Link](’$V(t)$ $[volts]$’)
[Link](’Respuesta a la rampa’)
[Link](loc=’lower right’)
[Link]()
fig3.set_size_inches(4,6)
[Link](’python_rampa.png’, dpi=600)
[Link](’python_rampa.pdf’)

# Respuesta a la funcion sinosoidal

fig4=[Link]()
[Link](t,u4, ’-’, color=[0.5,0.05,0.05], label=’$Pao(t)$’) #Entrada
ts,Vs=control.forced_response(sys,t,u4,x0)
[Link](t,Vs, ’-’, color=[0,0.25,0.4], label=’$PA(t)$’) #Salida (lazo abierto)
ts,pid=control.forced_response(sysPID,t,u4,x0)
[Link](t,pid, ’:’, linewidth=3, color=[0.3,0.5,0.2], label=’$PID$’) #Controlador (lazo cerrado)
[Link](True)
[Link](-0.5, 10)
[Link](-1.2, 1.2)
[Link](’$t$ $[segundos]$’)
[Link](’$V(t)$ $[volts]$’)
[Link](’Respuesta a la funcion sinosoidal’)
[Link](loc=’lower right’)
[Link]()
fig4.set_size_inches(4,6)
[Link](’python_sinusoidal.png’, dpi=600)
[Link](’python_sinusoidal.pdf’)
 73

Respuesta al escalón
1.2

1.0

0.8

0.6
V(t) [volts]

0.4

0.2

0.0 Pao(t)
PA(t)
PID
0 2 4 6 8 10
t [segundos]

Figura 4.28: Respuesta del sistema de mecánica pulmonar, en lazo abierto y lazo cerrado, al escalón
unitario en Python.
 74

Respuesta al impulso
1.2
Pao(t)
PA(t)
PID
1.0

0.8

0.6
V(t) [volts]

0.4

0.2

0.0

0.2
0 2 4 6 8 10
t [segundos]

Figura 4.29: Respuesta del sistema de mecánica pulmonar, en lazo abierto y lazo cerrado, al impulso
unitario en Python.
 75

Respuesta a la rampa
1.2

1.0

0.8

0.6
V(t) [volts]

0.4

0.2

0.0 Pao(t)
PA(t)
PID
0 2 4 6 8 10
t [segundos]

Figura 4.30: Respuesta del sistema de mecánica pulmonar, en lazo abierto y lazo cerrado, a la rampa
unitario en Python.
 76

Respuesta a la función sinosoidal

1.0

0.5
V(t) [volts]

0.0

0.5

1.0 Pao(t)
PA(t)
PID
0 2 4 6 8 10
t [segundos]

Figura 4.31: Respuesta del sistema de mecánica pulmonar, en lazo abierto y lazo cerrado, a la función
sinusoidal en Python.
 77

Capítulo 5

Prácticas

5.1. Práctica: Diseño de controladores para sistemas de segundo orden

Objetivo. Diseñar un controlador que permita eliminar el error en la respuesta del circuito RLC de la
Figura 5.1. El diagrama de bloques del lienzo principal se ilustra en la Figura 5.2.

R R
Ve (t) Vs (t)

i1 (t) i2 (t)
L R

R C

Figura 5.1: Diagrama eléctrico.

donde

L = _____ H;

C = _____ F;

R = _____ :

Actividades

1. Calcular analíticamente la función de transferencia del sistema.

2. Establecer el modelo de ecuaciones integro-diferenciales.

3. Determinar el error en estado estacionario y la estabilidad del sistema en lazo abierto.

 78

Figura 5.2: Diagrama de bloques.

4. Diseñar el controlador con ayuda de Simulink. Utilizar el bloque ‘PID Controller ’y la herramienta
‘Tune’para diseñar un controlador y sintonizar los valores óptimos para cada una de las ganancias
kP , kI y kD .

5. Determinar la respuesta al escalón, impulso, rampa unitaria y a la función sinusoidal

[u (t) = sin !t j ! = 250 mHz = =2 rad=s] en el intervalo t 2 [0; 10] (segundos), en Python,
Simulink y Multisim del circuito RLC como se indica a continuación:

a) Lazo abierto.

b) Lazo cerrado con el controlador diseñado.

6. Construir los circuitos en protoboard para comprobar los resultados obtenidos en las simulaciones.

7. Discutir los resultados obtenidos y elaborar el reporte de la práctica.

 79

5.2. Práctica: Comparación entre los controladores P, I, PD, PI y PID

Objetivo. Analizar la respuesta en el tiempo y estabilidad de un circuito RLC de segundo de las Figura
5.3 orden en lazo abierto y en lazo cerrado cuando se aplican los controladores P, I, PD, PI, y PID. El
diagrama de bloques del lienzo principal se ilustra en la Figura 5.4. donde

L C
Ve (t) Vs (t)

i1 (t) i2 (t)

R1 R2

Figura 5.3: Diagrama eléctrico del circuito RLC de segundo orden.

L = _____ H;

C = _____ F;

R1 = _____ ;

R2 = _____ :

Actividades

1. Calcular analíticamente la función de transferencia del sistema de mecánica pulmonar.

2. Establecer el modelo de ecuaciones integro-diferenciales.

3. Determinar el error en estado estacionario y la estabilidad del sistema en lazo abierto.

4. Utilizar las siguientes ganancias en los controladores:

a) Ganancia proporcional, kP = 10.

b) Ganancia integral, kI = 1 106 .

 80

Figura 5.4: Diagrama de bloques cuando se aplican los controladores P, I, PD, PI, y PID.

c) Ganancia derivativa, kD = 1 10 5.

5. Determinar la respuesta al escalón, impulso, rampa unitaria y a la función sinusoidal

[u (t) = sin !t j ! = 250 mHz = =2 rad=s] en el intervalo t 2 [0; 10] (segundos), en Python,
Simulink y Multisim del circuito RLC como se indica a continuación:

a) Lazo abierto.

b) Lazo cerrado con los controladores P, I, PD, PI y PID.

6. Construir los circuitos en protoboard para comprobar los resultados obtenidos en las simulaciones.

7. Discutir los resultados obtenidos y elaborar el reporte de la práctica.

 81

5.3. Práctica: Sistema musculoesquelético

Michael Khoo (2018), modeliza un compartimento del sistema musculoesquelético mediante el

diagrama mecánico de la Figura 5.5.

F0
Cs

x1 (t)

F (t) x1 (t) ! x(t)

x(t)

Cp

Figura 5.5: Diagrama mecánico del sistema musculoesquelético.

Descripción del sistema. F0 representa la fuerza desarrollada por el elemento contráctil activo del
músculo, mientras que F (t) es la fuerza real que resulta después de tener en cuenta las propiedades
mecánicas del músculo, y se asume que F0 = F (t), donde 0 < < 1. R representa la amortiguación
viscosa inherente al tejido, mientras que Cp (elemento elástico paralelo) y Cs (elemento elástico en
serie) re‡ejan las propiedades de almacenamiento elástico del sarcolema y los tendones musculares,
respectivamente.
La con…guración paralela se realiza para considerar las restricciones mecánicas impuestas a los
componentes del modelo. Si el resorte Cp se estira en una longitud incremental x (t), toda la combinación
en serie de R y Cs también se extenderá en la misma longitud. Además, la suma de la fuerza transmitida
a través de las dos ramas de la con…guración paralela debe ser igual a F (t). Aunque la suma de las
extensiones de Cs y R tendrá que ser igual a x (t), las contribuciones individuales de longitud de Cs y
R no necesitan ser iguales. Por lo tanto, si se asume que Cs se estira una longitud x1 (t), entonces la
extensión en la combinación paralela de R y F0 será x (t) x1 (t). La velocidad con la que se extiende
el amortiguador representado por R se obtiene al derivar x (t) x1 (t) con respecto al tiempo, es decir,
d [x (t) x1 (t)] =dt.
 82

Reference

Khoo MC. Physiological Control Systems Analysis, Simulation, and Estimation. Wiley, 2018. Section 2, page 26.

Actividades

1. Convertir el modelo mecánico a un modelo eléctrico, utilizar la analogía de fuerza-voltaje.

2. Utilizar los siguientes valores para los componentes del circuito: F (t) = 1 V , = 0:5, R = 100 ,
Cs = 10 F y Cp = 100 F .

3. Calcular analíticamente la función de transferencia del sistema musculoesquelético, utilice el

principio de superposición.

4. Determinar el error en estado estacionario y la estabilidad del sistema en lazo abierto.

5. Diseñar el controlador con ayuda de Simulink. Utilizar el bloque ‘PID Controller ’y la herramienta
‘Tune’para sintonizar los valores óptimos para cada una de las ganancias kP , kI y kD . Construir el
diagrama de bloques como se indica en el diagrama 5.6.

Figura 5.6: Diagrama de bloques.

6. Determinar la respuesta a la función escalón unitario en el intervalo t 2 [0; 10] (segundos), en Python,
Simulink y Multisim en lazo abierto y en lazo cerrado con el controlador.

7. Construir los circuitos en protoboard para comprobar los resultados obtenidos en las simulaciones.

8. Discutir los resultados obtenidos y elaborar el reporte de la práctica.

 83

5.4. Práctica: Sistema respiratorio

Michael Khoo (2018), modeliza el sistema respiratorio mediante un circuito RLC de tercer orden. El
circuito eléctrico que representa esta dinámica se muestra en la Figura 5.7.

RC LC Paw (t) RP
Pao (t) PA (t)

Q(t) QA (t)
CL
CS Q(t) ! QA (t)
Ppl (t)

CW

P0

Figura 5.7: Sistema respiratorio modelizado mediante un circuito RLC de tercer orden.

Descripción del sistema. Las vías respiratorias se dividen en dos categorías: las vías respiratorias
más grandes o centrales y las vías respiratorias más pequeñas o periféricas, con resistencias mecánicas al
‡uido iguales a RC y RP , respectivamente. El efecto de la inercia al ‡ujo de gas en las vías respiratorias
centrales está dado por LC . El aire que ingresa a los alvéolos también produce una expansión de la cavidad
de la pared torácica con el mismo volumen. Esto está representado por la conexión de las conformidades
pulmonares (CL ) y de la pared torácica (CW ) en serie. Sin embargo, una pequeña fracción del volumen
de aire que ingresa al sistema respiratorio se desvía de los alvéolos como resultado de la distensibilidad
de las vías respiratorias centrales y la compresibilidad del gas. Este volumen derivado es muy pequeño en
circunstancias normales a frecuencias respiratorias regulares, pero se vuelve progresivamente sustancial si
una enfermedad conduce a una obstrucción de las vías respiratorias periféricas (es decir, un aumento de
RP ) o rigidez de los pulmones o de la pared torácica (es decir, disminución de CL o CW ). Este efecto se
considera mediante el elemento de derivación, CS , en paralelo con CL y CW .
Las presiones desarrolladas en los diferentes puntos de este modelo pulmonar son Pao (t) en la apertura
de la vía aérea pulmonar, Paw (t) en las vías respiratorias centrales, PA (t) en los alvéolos y Ppl (t) en el
espacio pleural (entre el parénquima pulmonar y la pared torácica), estas presiones se re…eren a P0 (t), la
presión ambiental, que se puede establecer en cero. El caudal volumétrico de aire que ingresa al sistema
respiratorio está dado por Q (t) y el ‡ujo entregado a los alvéolos por QA (t). Por lo tanto, el ‡ujo derivado
 84

de los alvéolos debe ser Q (t) QA (t).

Ahora, para un individuo sano, los valores de los parámetros para cada elemento son los siguientes:

1
RC = 1 cmH2 O s L ;

LC = 0:01 cmH2 O s2 L 1
;
1
CS = 0:005 L (cmH2 O) ;
1
RP = 0:5 cmH2 O s L ;
1
CL = 0:2 L (cmH2 O) ;
1
CW = 0:2 L (cmH2 O) ;

donde 1 cmH2 O se de…ne como un centímetro de agua, equivalente a 98:0665 P a (pascales), o 0:7356
mmHg (milímetros de mercurio), es decir, es una unidad de presión.
Mientras tanto, para un paciente con en…sema, se proponen los siguientes valores asumiendo una
distensibilidad pulmonar y una resistencia de las vías respiratorias periféricas:

1
RP = 7:5 (cmH2 O) s L ;
1
CL = 0:4 L (cmH2 O) ;

y los valores numéricos de los demás parámetros son los mismos.

Reference

Khoo MC. Physiological Control Systems Analysis, Simulation, and Estimation. Wiley, 2018. Section 2, page 23.

Actividades

1. Calcular analíticamente la función de transferencia del sistema pulmonar.

2. Establecer el modelo de ecuaciones integro-diferenciales.

3. Determinar el error en estado estacionario y la estabilidad del sistema en lazo abierto.

4. Diseñar el controlador con ayuda de Simulink. Utilizar el bloque ‘PID Controller ’y la herramienta
‘Tune’para sintonizar los valores óptimos para cada una de las ganancias kP , kI y kD . Construir el
diagrama de bloques como se indica en el diagrama 5.8.

5. Ilustrar el cambio del ‡ujo de aire y el volumen tidal en respuesta a las siguientes formas de onda
de presión sinusoidal en la apertura de la vía aérea [Pao (t)] :

a) 15 respiraciones por minuto con una amplitud (A) de 2:5 cmH2 O, es decir, respiración normal.
 85

Figura 5.8: Diagrama de bloques para ajustar las ganancias del controlador con base en la respuesta del
individuo sano.

b) 30 respiraciones por minuto con una amplitud (A) de 1:5 cmH2 O, es decir, respiración elevada
o taquipnea.

6. Determinar la respuesta a la función sinusoidal [u (t) = A sin !t] en el intervalo t 2 [0; 30] (segundos),
en Python, Simulink y Multisim en lazo abierto y en lazo cerrado con el controlador.

7. Discutir los resultados obtenidos y elaborar el reporte de la práctica.

 86

5.5. Práctica: Sistema cardiovascular

En el estudio de las enfermedades cardiovasculares se puede obtener mucha información a partir de las
características hemodinámicas del sistema vascular, como la resistencia periférica total, la distensibilidad
arterial total y la impedancia característica de la aorta proximal. Un modelo utilizado para describir las
características vasculares es el denominado modelo del efecto Windkessel. Otto Frank sentó las bases de
este modelo a principios del siglo XX y postuló que la aorta podría representarse mediante un manómetro
con una sección distensible concentrada y una resistencia periférica también concentrada. A partir de lo
anterior, la analogía eléctrica se realiza mediante una combinación de un capacitor y una resistencia en
paralelo. En investigaciones posteriores, este modelo clásico de dos elementos se amplió con un tercer
elemento para tener en cuenta la impedancia de la parte proximal del lecho vascular pulmonar y a un
cuarto elemento para tener en cuenta la inercia de la sangre.

Descripción del sistema. El modelo de Windkessel de cuatro elementos contiene dos elementos
dinámicos. Por lo tanto, se necesitan dos estados para describir la dinámica. El vector de estados se
conforma por las variables FL (t) denotando el ‡ujo a través de la inercia arterial total, y la variable Pp (t)
representando la presión sobre la distensibilidad arterial. Entonces, asumiendo Pa (t) como la presión
arterial de entrada, y en consecuencia a Fa (t) como el ‡ujo hacia la aorta o arteria pulmonar, el modelo
de cuatro elementos se observa en la Figura 5.9.

FZ (t)
Pa (t) Pp (t)
Fa (t) L

FL (t) FC (t) FR (t)

C R

Figura 5.9: Modelo Windkessel de cuatro elementos del sistema cardiovascular.

 87

Los parámetros son Z, C, R y L, que representan respectivamente la impedancia característica del

lecho vascular pulmonar (aorta y arteria pulmonar), la distensibilidad aérea total, la resistencia periférica
y la inertancia arterial. Este modelo tiene muchas ventajas importantes, por ejemplo:

1. Su sencillez, unos pocos elementos interconectados son su…cientes para reproducir la dinámica
principal del sistema cardiovascular.

2. Existe una clara analogía entre los elementos eléctricos y los componentes hidráulicos implicados en
el efecto Windkessel. En consecuencia, se relacionan fácilmente con el signi…cado hemodinámico y
el acoplamiento ventrículo-arterial.

Reference

Kind, T., Faes, T. J., Lankhaar, J. W., Vonk-Noordegraaf, A., and Verhaegen, M. Estimation of three-and four-

element windkessel parameters using subspace model identi…cation. IEEE Transactions on Biomedical Engineering,

2010, 57(7), 1531–1538.

Actividades

1. Calcular analíticamente la función de transferencia del sistema pulmonar.

2. Establecer el modelo de ecuaciones integro-diferenciales.

3. Determinar el error en estado estacionario y la estabilidad del sistema en lazo abierto.

4. Diseñar el controlador con ayuda de Simulink. Utilizar el bloque ‘PID Controller ’y la herramienta
‘Tune’para sintonizar los valores óptimos para cada una de las ganancias kP , kI y kD . Construir el
diagrama de bloques como se indica en el diagrama 5.10.

5. Ilustrar el cambio de la presión sobre la distensibilidad arterial [Pp (t)] en respuesta a la

presión arterial de entrada Pa (t) = 0:8 + sin !t, calcule ! para 60 latidos por minuto
[!(rad=s) = 2 f (Hz)]. Utilice los siguientes valores para representar diferentes estados en el
paciente.
Parámetro Hipotenso Normotenso Hipertenso Unidades
Z 0:02 0:033 0:05 mmHg s ml 1

1
C 2 1:5 0:7 ml (mmHg)
R 0:6 0:95 1:4 mmHg s ml 1

L 0:005 0:01 0:02 mmHg s2 ml 1

 88

Figura 5.10: Diagrama de bloques para ajustar las ganancias del controlador con base en la respuesta del
individuo sano.

6. Determinar la respuesta a la función sinusoidal en el intervalo t 2 [0; 10] (segundos), en Python,

Simulink y Multisim en lazo abierto y en lazo cerrado con el controlador.

7. Elaborar el diagrama biológico del sistema con [Link].

8. Discutir los resultados obtenidos y elaborar el reporte de la práctica.

 89

5.6. Práctica: Sistema digestivo

Utilice el circuito RLC mostrado en la Figura 5.11 para describir la dinámica sobre la ingesta
de alimento desde el esófago hacia el estómago y los intestinos delgado y grueso. Los resistores
pueden representar la resistencia …siológica de las capas mucosas y los capacitores las capacidades de
almacenamiento del estómago e intestinos.

Ve (t)
i0 (t)
.R
Ve (t)

i1 (t) C R i2 (t)

-1 R
i3 (t) i4 (t)
,1 C

Vs (t)

-2 R
i5 (t) i6 (t)
,2 C

i0 (t) L
Vs (t)

Figura 5.11: Analogía del sistema digestivo mediante un diagrama eléctrico. Dibujo del sistema digestivo
elaborado con [Link].

Actividades

1. Describir, con base en suposiciones …siológicas, el aparato digestivo utilizando el circuito RLC.

2. Calcular analíticamente la función de transferencia del sistema digestivo.

3. Establecer el modelo de ecuaciones integro-diferenciales.

 90

4. Determinar el error en estado estacionario y la estabilidad del sistema en lazo abierto.

5. Diseñar el controlador con ayuda de Simulink. Utilizar el bloque ‘PID Controller ’y la herramienta
‘Tune’para sintonizar los valores óptimos para cada una de las ganancias kP , kI y kD . Construir el
diagrama de bloques como se indica en el diagrama 5.12.

6. Elegir una entrada, justi…cándola …siológicamente, para obtener la respuesta en el intervalo t 2 [0; 10]
(segundos), en Python, Simulink y Multisim en lazo abierto y en lazo cerrado con el controlador.

7. Construir los circuitos en protoboard para comprobar los resultados obtenidos en las simulaciones.

8. Discutir los resultados obtenidos y elaborar el reporte de la práctica.

Figura 5.12: Diagrama de bloques.

Instrucciones generales
La Práctica se va a entregar en formato de artículo de investigación para que se tome de ejemplo para
el Proyecto Final. La plantilla de LATEX se puede descargar del siguiente vínculo:

Descargar plantilla.
 91

5.7. Práctica: Sistemas …siológicos en el espacio de estados: EDOs

El sistema simpli…cado de la mecánica pulmonar es representado por la Figura 4.1, la cual se muestra
a continuación donde el resistor R representa una combinación de la resistencia al ‡ujo de aire [Q (t)]

L R Pao (t)
Pao (t) PA (t) Q(t)
L
Q(t) R
C

P0 C PA (t)

Figura 5.13: Modelo eléctrico de la mecánica pulmonar. Sistema pulmonar elaborado con [Link].

en las vías respiratorias, el tejido pulmonar y la pared torácica. El inductor L representa la inertancia
al ‡ujo de aire en las vías respiratorias. El capacitor C representa la compliancia combinada de las vías
respiratorias, el tejido pulmonar y la pared torácica. Los componentes R y C representan respectivamente
las propiedades mecánicas resistivas y de almacenamiento del sistema respiratorio.

Referencia

Khoo MC. Physiological Control Systems Analysis, Simulation, and Estimation. Wiley, 2018. Section 4, page 93.

El modelo de ecuaciones integro-diferenciales del sistema se muestra a continuación:

Z
dQ (t) 1 1
Q (t) = Pao (t) L Q (t) dt ;
dt C R

Z
1
PA (t) = Q (t) dt:
C

mientras que la función de transferencia está dada por la siguiente expresión:

PA (s) 1
= 2
;
Pao (s) LCs + RCs + 1

donde

R = 10 k ;

L = 1 H;

C = 220 F:
 92

A partir de estos resultados, realice las siguientes actividades.

Actividades

1. Establecer el modelo de Ecuaciones Diferenciales Ordinarias (EDOs) de primer orden con los
siguientes dos procedimientos:

a) Utilizar la función de transferencia para obtener el sistema de EDOs.

b) Aplicar cambio de variable.

2. Para el circuito RLC de segundo orden de la Figura 5.14, establecer el modelo de EDOs con el
procedimiento de preferencia.

L
Ve (t)
VL (t)

iL (t) iR (t) R C VC (t)

Figura 5.14: Circuito RLC.

3. Determinar la respuesta de ambos sistemas al escalón, impulso, rampa unitaria y a la función

sinusoidal [u (t) = sin !t j ! = 250 mHz = =2 rad=s] en el intervalo t 2 [0; 10] (segundos), en
Python aplicando la función [Link].

4. Discutir los resultados obtenidos y elaborar el reporte de la práctica.

 93

5.8. Práctica: Regeneración de glóbulos rojos

Los glóbulos rojos son células discoides bicóncavas responsables del suministro de oxígeno a los tejidos
desde los pulmones, y del transporte de dióxido de carbono desde los tejidos hacia a los pulmones. El
transporte de oxígeno se realiza por el complejo proteico contenido en la hemoglobina. Los humanos
adultos sanos tienen un total de 2 1013 a 3 1013 eritrocitos en un cualquier instante de tiempo los
hombres tienen alrededor de 5 a 6 millones de eritrocitos por microlitro de sangre y las mujeres tienen
alrededor de 4 a 5 millones de eritrocitos por microlitro de sangre.

Célula precursora
Producción hematopoyética

Riñón Proeritroblastos

Eritropoyetina Eritrocitos

Reducción
Oxigenación tisular

Disminución Reducción
de eritrocitos

Figura 5.15: Eritropoyesis. Elaborado con [Link].

La eritropoyesis (proceso de producción de glóbulos rojos, vea la Figura 5.14) comienza en la médula
ósea, donde las células madre se transforman a través de varias etapas hasta ser descargadas en el torrente
sanguíneo como reticulocitos sanguíneos, donde maduran rápidamente hasta convertirse en eritrocitos
completamente desarrollados, proceso que tiene una duración de aproximadamente 20 días. El tiempo de
vida media de los glóbulos rojos en el torrente sanguíneo es de unos 120 días, después de este tiempo son
destruidos por los fagocitos. Este proceso es necesario ya que los glóbulos rojos ya no pueden dividirse
ni regenerarse. El objetivo de la eritropoyesis es controlar el número de glóbulos rojos de modo que el
suministro de oxígeno a los tejidos coincida con la demanda en el cuerpo.
El siguiente modelo matemático de tres EDOs de primer orden es un modelo mecanicista de
 94

compartimento para la eritropoyesis después de la pérdida de sangre en personas sanas, fenómeno que se
puede modelizar como un proceso dinámico no lineal.
(Base x3 ) x1
x_ 1 = (X0 k1 x1 ) + ; (5.1)
Base
x_ 2 = (k1 x1 k2 x2 ) ; (5.2)

x_ 3 = (k2 x2 x3 ) ; (5.3)

donde x1 (t) representa el compartimento de células precursoras eritroides altamente proliferantes con
respecto a la eritropoyetina, x2 (t) describe el compartimento no proliferante de células precursoras
eritroides con respecto a la eritropoyetina, y x3 (t) el compartimento para los eritrocitos maduros y
reticulocitos sanguíneos.
Con respecto a los parámetros, X0 re‡eja la cantidad absoluta de células que se destinan al linaje
eritroide y que maduran en el primer compartimento de células precursoras eritroides. Las tasas de
transición y las tasas de mortalidad entre los compartimentos están dadas por k1 , k2 y , estas tasas
son independientes de la hormona eritropoyetina. La compensación de la pérdida de sangre se describe
mediante un término de retroalimentación de los eritrocitos a las células en proliferación basado en la
pérdida fraccional de eritrocitos. Con base en lo anterior, se introducen los parámetros y , que se utilizan
para la descripción de las características individuales de la eritropoyesis. Se asume que cada individuo
tiene un recuento medio de eritrocitos indicado por el parámetro Base . A continuación, se presenta la
tabla de valores y unidades de los parámetros.

Parámetro Valor Unidades

Base 885:42 g
X0 7:3785 g=d{a
k1 0:125 d{as 1

k2 0:1667 d{as 1

0:00833 d{as 1

1:28
0:56 d{as 1

donde las unidades para xi (t) están dadas por g=d{a;es decir, gramos de las células por día. Si se considera
la variable de control como una transfusión sanguínea [u (t)], el modelo se reescribe como se muestra a
continuación
(Base x3 ) (1 u (t)) x1
x_ 1 = (X0 k1 x1 ) + ; (5.4)
Base
x_ 2 = (k1 x1 k2 x2 ) ; (5.5)

x_ 3 = (k2 x2 x3 ) + u (t) x3 : (5.6)

 95

Referencia

Tetschke M, Lilienthal P, Pottgiesser T, Fischer T, Schalk E, Sager S. Mathematical Modeling of RBC Count

Dynamics after Blood Loss. Processes. 2018; 6(9):157. [Link]

Actividades

1. Resolver el sistema de EDOs de primer orden de eritropoyesis (5.1)–(5.3) en Python aplicando la

función [Link].

2. Resolver el sistema de EDOs de primer orden de eritropoyesis (5.1)–(5.3) en Python aplicando un

método numérico de diferencias …nitas.

3. Utilizar las condiciones iniciales x1 (0) = 50, x2 (0) = 30 y x3 (0) = 650.

4. Considerar el caso de transfusión sangínea al individuo mediante la siguiente función:

8
< 1; 0 t 1;
u (t) =
: 0; 1 < t < 1;

y resolver el sistema en Python aplicando el método numérico de diferencias …nitas.

5. Discutir los resultados obtenidos y elaborar el reporte de la práctica.

 96

Capítulo 6

Proyecto …nal

6.1. Objetivo del proyecto …nal

Diseñar un circuito RLC que modelice la dinámica de un sistema …siológico y un

controlador que elimine el error entre la entrada y la salida del sistema.

La plantilla de LATEX la pueden descargar del siguiente vínculo:

Descargar plantilla: Proyecto …nal

Consideraciones para la evaluación del proyecto:

1. Elaborado en español con una buena ortografía y redacción tendrá una cali…cación máxima de 20 %
de la cali…cación …nal.

2. Elaborado en inglés con una buena ortografía y redacción tendrá una cali…cación máxima de 30 %
de la cali…cación …nal.

3. La cali…cación se asignará con base en el desarrollo de cada sección en el artículo.

4. La plantilla proporcionada es la sugerida, sin embargo, pueden seleccionar una diferente de las
disponibles en el software de preferencia.

El contenido del artículo se describe a continuación, sin embargo, se debe utilizar la Práctica: Mecánica
pulmonar, como la base para desarrollar cada sección del proyecto.

6.2. Descripción del sistema

En esta sección se debe describir el sistema propuesto con base en la literatura consultada. A
continuación, se presentan ejemplos de referencias que deben incluir en esta sección:

1. Página de internet: Matlab. Dynamic system, [Link], 2019, Último acceso:

25/Marzo/2020.
 97

2. Artículo de investigación: Valle Paul A., Coria Luis N., Plata, Corina., Salazar Yolocuauhtli.
CAR-T Cell Therapy for the Treatment of ALL: Eradication Conditions and In Silico
Experimentation. Hemato 2.3 (2021): 441-462, doi: 10.3390/hemato2030028.

3. Libro: Gar…nkel, A., Shevtsov, J. and Guo, Y., Modeling life: the mathematics of biological systems,
Springer, 2017, doi: 10.1007/978-3-319-59731-7.

El archivo ‘[Link]’contiene el código necesario para mostrar estas referencias.

Se puede utilizar cualquier proceso …siológico del cuerpo humano para realizar la propuesta del sistema.
El modelizado inicial del sistema propuesto se debe presentar mediante un circuito eléctrico con las
siguientes características:

1. Debe contener al menos dos mallas o un nodo con al menos cuatro elementos.

2. La función de transferencia debe ser al menos de segundo orden. Se pueden utilizar dos capacitores,
dos inductores o un capacitor y un inductor.

3. La señal de entrada del sistema depende explícitamente del proceso que se modeliza. Si el sistema
lo requiere, se pueden utilizar distintos tipos de señales o distintas frecuencias de la misma señal
para representar el proceso …siológico.

4. Elaborar el diagrama biológico del sistema con [Link].

En la Figura 6.1 se ilustra el ejemplo del modelo eléctrico de la mecánica pulmonar indicando cada
uno de los elementos, voltajes y corrientes involucradas en la dinámica del circuito. Se deben dejar claros
los puntos que representan la entrada y la salida en el sistema.
Los valores y unidades del circuito se deben escribir en un arreglo como se indica a continuación:

1
RC = 1 cmH2 O s L ;

LC = 0:01 cmH2 O s2 L 1
;
1
CS = 0:005 L (cmH2 O) ;
1
RP = 0:5 cmH2 O s L ;
1
CL = 0:2 L (cmH2 O) ;
1
CW = 0:2 L (cmH2 O) ;
 98

RC LC Paw (t) RP
Pao (t) PA (t)

Q(t) QA (t)
CL
CS Q(t) ! QA (t)
Ppl (t)

CW

P0

Figura 6.1: Modelo eléctrico de la mecánica pulmonar.

Adicionalmente, se deben identi…car los elementos cuyo valor cambia cuando se modeliza el sistema
…siológico con una enfermedad, por ejemplo, los elementos

1
RP = 7:5 (cmH2 O) s L ;
1
CL = 0:4 L (cmH2 O) ;

toman estos valores cuando el sistema presenta en…sema pulmonar.

Se recomienda que las Figura sea en formato PDF o PNG, con una buena resolución ( 600 DP I) y
las respectivas etiquetas legibles o descripción detallada.

6.3. Modelo matemático

En esta sección se debe presentar el desarrollo su…ciente para obtener el modelo de ecuaciones integro-
diferenciales y la función de transferencia del sistema. Se deben identi…car las ecuaciones principales del
sistema:
Z
dQ (t) 1
Pao (t) = Q (t) RC + LC + [Q (t) QA (t)] dt; (6.1)
dt CS
Z Z
1 1 1
[Q (t) QA (t)] dt = QA (t) RP + + QA (t) dt; (6.2)
CS CL CW
Z
1 1
PA (t) = + QA (t) dt; (6.3)
CL CW
a partir de las cuales se obtiene el modelo de ecuaciones integro-diferenciales:

1
Q (t) = [:::] ; (6.4)
RC
 99

1
QA (t) = [:::] ; (6.5)
RP
y la función de transferencia al aplicar la transformada de Laplace a las Ecuaciones (6.1)–(6.3), como se
indica a continuación:

Pao (s) = :::; (6.6)

PA (s) = :::; (6.7)

::: = :::; (6.8)

al realizar las operaciones algebraicas correspondientes con las Ecuaciones (6.6)–(6.8) se obtiene como
resultado la siguiente función de transferencia:

PA (s) 0
= 3 2
:
PaO (s) 3s + 2s + 1s + 0

Realizar una breve discusión sobre la utilidad del modelo de ecuaciones integro-diferenciales y la
función de transferencia para analizar sistema propuesto. Seguir el procedimiento desarrollado en la
Práctica 8.

6.4. Simulaciones numéricas

En esta sección se deben presentar y describir el desarrollo de las simulaciones numéricas. El único
software obligatorio para realizar las simulaciones es Simulink por el diseño del controlador. Primero, se
debe presentar el diagrama de bloques principal como se indica en la Figura 6.8.

Figura 6.2: Diagrama de bloques principal.

Con el siguiente orden se deben presentar las Figuras correspondientes a las demás simulaciones:
 100

1. Subsistema con el diagrama de bloques correspondiente al modelo matemático de ecuaciones integro-

diferenciales.

2. Simulación comparando la señal de entrada, la respuesta del individuo saludable y la del individuo
enfermo. Esto se debe repetir para el caso en que se utilicen diferentes señales de entrada o diferentes
frecuencias en la señal.

3. Evidencia del diseño del controlador PID. Se debe realizar la sintonización de las ganancias kP , kI
y kD con la herramienta Tune.

4. Simulación comparando la señal de entrada, la respuesta del individuo saludable, la del individuo
enfermo y la del individuo enfermo con el controlador diseñado. Esto se debe repetir para el caso
en que se utilicen diferentes señales de entrada o diferentes frecuencias en la señal.

Se deben indicar los valores de los parámetros del bloque de la señal de entrada, sin embargo, no es
necesario una Figura de esto. Recomendaciones generales para la elaboración de las Figuras:

1. Todas las Figuras deben tener una descripción breve pero su…ciente para comprender lo que se está
ilustrando.

2. Acomodar las Figuras de tal manera que no quede solamente una por página.

3. El fondo de las Figuras debe de ser en color blanco.

4. La información de los ejes x (tiempo) y y (nivel de voltaje), se debe observar de manera clara.

5. Dependiendo de cada sistema, 1 segundo puede representar 1 hora, 1 d{a, etc.

6. Se recomienda que todas las Figuras sean en formato PDF o PNG para conservar la calidad.

6.5. Implicaciones …siológicas (Conclusiones)

En esta sección se deben discutir las implicaciones …siológicas considerando los siguientes puntos:

1. La utilidad del modelizado matemático para describir un sistema …siológico.

2. La comparación entre el sistema con los valores del individuo sano y los valores del individuo enfermo.

3. El tipo de controlador diseñado para llevar la respuesta del individuo enfermo a la del individuo
sano.
 101

4. La analogía del controlador como una terapia real para tratar la enfermedad.

5. La importancia de poder realizar simulaciones numéricas para explorar diferentes estrategias y evitar
hacer pruebas en pacientes reales.