using SymPy, Plots, Markdown, LaTeXStrings, Printf
# Definimos las variables de nuestro sistema
# Especificamos que la constante de tiempo es siempre positiva
@syms t::real A::real K_p::real tau_p::positive s(t, A, K_p, tau_p, s)4 Sistemas lineales de primer orden
4.1 Definición de sistema lineal de primer orden
Un sistema de primer orden es aquel cuya salida \(y(t)\) puede ser modelada por una ecuación diferencial de primer orden como:
\[ a_1 \frac{\mathrm{d}y (t)}{\mathrm{d}t} + a_0 y (t) = b f (t) \tag{4.1}\]
donde \(f(t)\) es la entrada al sistema.
Si \(a_0 \neq 0\):
\[\frac{a_1}{a_0} \frac{\mathrm{d}y (t)}{\mathrm{d}t} + y (t) = \frac{b}{a_0} f(t)\]
Si se define \(\frac{a_1}{a_0} = \tau_p\) y \(\frac{b}{a_0} = K_p\) y se sustituye en la ecuación anterior se obtiene:
\[ \tau_p \frac{\mathrm{d}y (t)}{\mathrm{d}t} + y (t) = K_p f (t) \tag{4.2}\]
donde:
\(\tau_p\) es la constante de tiempo del proceso
\(K_p\) es la ganancia del proceso
Si \(y(t)\) y \(f(t)\) están definidos mediante la utilización de variables de desviación alrededor del estado estacionario, las condiciones iniciales son \(y(0)=0\) y \(f(0)=0\).
Operando se encuentra la función de transferencia de un proceso de primer orden:
\[G (s) = \frac{K_p}{\tau_p s + 1}\]
Los sistemas de primer orden son los más frecuentes en los procesos de la industria alimentaria, por ello su estudio es de gran importancia. Estos sistemas se caracterizan por:
Su capacidad de almacenar materia, energía o cantidad de movimiento. Esta capacidad está directamente relacionada con la ganancia del proceso.
Una resistencia asociada con el caudal de materia, energía o cantidad de movimiento. Esta resistencia o inercia viene dada por la constante de tiempo.
En el caso particular de que \(a_0 = 0\):
\[G (s) = \frac{K_p}{\tau_p s}\]
Se trata de aquellos sistemas de primer orden denominados integradores puros y se hablará de ellos más adelante.
4.2 Respuesta a una entrada en escalón
Para un escalón de altura \(A\) y un sistema de primer orden la salida \(y(s)\) es:
\[y(s) = \frac{A}{s} \frac{K_p}{\tau_p s + 1} \tag{4.3}\]
En tiempo real, invirtiendo las transformadas de Laplace, se obtiene:
\[y(t) = A K_p \left( 1 - \mathrm{e}^{-\frac{t}{\tau_p}} \right) \tag{4.4}\]
Podemos reproducir este cálculo utilizando SymPy. Como es habitual, empezamos definiendo las variables necesarias y cargando las librerías que vamos a utilizar:
Definición de la función de transferencia del proceso, \(G(s)\):
G = K_p/(tau_p*s + 1)\(\frac{K_{p}}{s \tau_{p} + 1}\)
La función de entrada es un escalón de altura \(A\), \(f(s)\):
f = A/s\(\frac{A}{s}\)
Por lo tanto, la respuesta del proceso \(y(s)\) es:
y_s = G*f\(\frac{A K_{p}}{s \left(s \tau_{p} + 1\right)}\)
Realizando la transformada inversa de Laplace de \(y(s)\) obtenemos la respuesta del proceso en tiempo real, \(y(t)\):
y = sympy.inverse_laplace_transform(y_s, s, t)\(A K_{p} \left(\theta\left(t\right) - e^{- \frac{t}{\tau_{p}}} \theta\left(t\right)\right)\)
Representando la función en coordenadas adimensionales, \(\frac{y(t)}{A K_p}\) frente a \(\frac{t}{\tau_p}\), se obtiene la típica salida de un sistema de primer orden:
4.3 Propiedades de un sistema de primer orden
Cabe destacar las siguientes características de cualquier sistema de primer orden:
4.3.1 Autorregulación
El proceso alcanza un nuevo estado estacionario sin necesidad de un sistema de control, tal como se muestra en la Figura 4.1.
Se puede comprobar que la respuesta a tiempos largos se encuentra acotada o que tiende a un valor concreto.
La salida del proceso en el nuevo estado estacionario es:
\[\lim_{t \to \infty} y (t) = A K_p\]
Podemos comprobar
sympy.limit(y, t, oo)\(A K_{p}\)
Cuanto mayor es la ganancia menor debe ser la entrada del sistema (perturbación) para producir el mismo efecto final.
Puede ocurrir que la constante \(a_0\) de la Ecuación 4.1 sea nula. Este tipo de procesos se conocen como integradores puros, ya que la salida es la integral de la entrada con el tiempo. Estos procesos pueden ser difíciles de controlar debido a que no presentan autorregulación. Ejemplos comunes de este tipo de sistemas son los tanques con líquidos, depósitos de gases y sistemas de almacenamiento de materias primas y productos.
4.3.2 Velocidad de la respuesta
Para calcular la velocidad de la respuesta hay que encontrar la pendiente de la respuesta:
\[\frac{\mathrm{d}\left[ \frac{y (t)}{A K_p} \right]}{\mathrm{d}t} = \frac{1}{\tau_p} e^{- \frac{t}{\tau_p}}\]
Para \(t=0\):
\[\frac{\mathrm{d}\left[ \frac{y (t)}{A K_p} \right]}{\mathrm{d}t} = \frac{1}{\tau_p}\]
Cuanto mayor sea \(\tau_p\), menor será la pendiente inicial de la respuesta del sistema y mayor será el tiempo necesario en alcanzar el nuevo estado estacionario.
Evaluando la Ecuación 4.4 para diferentes tiempo se obtiene la siguiente tabla:
| Tiempo transcurrido | \(1 \tau_p\) | \(2 \tau_p\) | \(3 \tau_p\) | \(4 \tau_p\) |
|---|---|---|---|---|
| \(y(t)\) respecto a su valor estacionario | 0.632 | 0.865 | 0.950 | 0.982 |
El cálculo de estos valores es muy sencillo:
# Calculamos los valores de y(t) para 0, 1, 2, 3 y 4
println("τ y(t)/A")
for τ ∈ 1:4
output = y(A=>1, K_p=>1, tau_p=>1, t=>τ)
@printf "%d %.3f" τ output
println()
endτ y(t)/A
1 0.632
2 0.865
3 0.950
4 0.982Transcurrido cuatro veces la constante de tiempo del proceso se puede asegurar que ha llegado el sistema al nuevo estado estacionario.
4.4 Respuesta a una función impulso
Al introducir un impulso de área \(A\) se obtiene la siguiente respuesta:
\[y (s) = \frac{K_p}{\tau_p s + 1} A\]
que en tiempo real es:
\[y (t) = \frac{K_p A}{\tau_p} \mathrm{e}^{- \frac{t}{\tau_p}}\]
Como en el caso de una entrada en escalón, podemos reproducir con SymPy el cálculo de la respuesta de un proceso de primer orden a una entrada en escalón:
using SymPy, Plots, LaTeXStrings
@syms t::real K_p::real A::real tau_p::positive s
f = A
G = K_p/(tau_p*s + 1)
y = sympy.inverse_laplace_transform(G*f, s, t)\(\frac{A K_{p} e^{- \frac{t}{\tau_{p}}} \theta\left(t\right)}{\tau_{p}}\)
De forma adimensional se puede escribir:
\[\frac{y (t)}{K_p A} = \mathrm{e}^{- \frac{t}{\tau_p}}\]
Se obtiene la función simétrica a la respuesta a una entrada en escalón, lo que implica que tiene las mismas características.
4.5 Respuesta a una función sinusoidal
Para una entrada sinusoidal de tipo:
\[f (t) = M U (t) \sin \omega t\]
se obtiene la siguiente respuesta:
\[y (s) = \frac{K_p}{\tau_p s + 1} \frac{M \omega}{s^2 + \omega^2}\]
que en tiempo real es:
\[y (t) = \frac{K_p M}{1 + \omega^2 \tau_p^2} \left( \omega \tau_p \mathrm{e}^{- \frac{t}{\tau_p}} \cos \omega t + \sin \omega t \right)\]
Aplicando la propiedad trigonométrica:
\[x \sin \alpha + y \cos \alpha = z \sin (\alpha + \varphi)\]
donde \(z^2 = x^2 + y^2\) y \(\tan \varphi = \frac{y}{x}\), se obtiene:
\[y (t) = \frac{K_p M \omega \tau_p}{1 + \omega^2 \tau_p^2} \mathrm{e}^{- \frac{t}{\tau_p}} + \frac{K_p M}{\sqrt[]{1 + \omega^2 \tau_p^2}} \sin (\omega t + \varphi) \tag{4.5}\]
donde \(\varphi = \mathrm{atan} (- \omega \tau_p)\).
El primer término de la respuesta es un término transitorio, ya que tiende a cero cuando el tiempo tiende a infinito. Este término pierde importancia para tiempos grandes:
\[ y (t) = \overbrace{\underset{\begin{array}{l} \big\downarrow t \to \infty \\ 0 \end{array}}{\frac{K_p M \omega \tau_p}{1 + \omega^2 \tau_p^2} \mathrm{e}^{- \frac{t}{\tau_p}}}}^{\text{Transitorio}} + \overbrace{\frac{K_p M}{\sqrt{1 + \omega^2 \tau_p^2}} \sin (\omega t + \varphi)}^\text{Estacionario} \]
La respuesta obtenida es de tipo sinusoidal con la misma frecuencia de oscilación \(\omega\) que la entrada pero con un desfase \(\varphi\). Además, el desfase está directamente relacionado con la frecuencia angular. Al aumentar el desfase, aumenta el desfase.
En el caso de que la frecuencia angular tienda a infinito, el desfase tiende a \(\frac{\pi}{2}\), que es el desfase máximo.
También se puede obtener un resultado equivalente utilizando SymPy:
using SymPy, Plots
@syms t::real Kp::real=>"K_p" M::real w::real=>"omega" Tp::real=>"tau_p" s
L(f) = sympy.laplace_transform(f, t, s, noconds=true)
iL(f) = sympy.inverse_laplace_transform(f, s, t)
f = L(M*sin(w*t))
G = Kp/(Tp*s+1)
y = iL(G*f)\(K_{p} M \omega \left(\frac{\tau_{p}^{2} e^{- \frac{t \left(\omega^{2} \tau_{p}^{2} + 1\right)}{\omega^{2} \tau_{p}^{3} + \tau_{p}}} \theta\left(t\right)}{\omega^{2} \tau_{p}^{3} + \tau_{p}} + \left(- \frac{\tau_{p} \cos{\left(\omega t \right)}}{\omega^{2} \tau_{p}^{2} + 1} + \frac{\sin{\left(t \left|{\omega}\right| \right)}}{\left(\omega^{2} \tau_{p}^{2} + 1\right) \left|{\omega}\right|}\right) \theta\left(t\right)\right)\)
4.6 Problemas
Problema 4.1 ★
Estudiar la respuesta de un proceso de función de transferencia:
\[G(s) = K_p \frac{\tau s +1}{\tau_p s +1}\]
con las siguientes entradas:
- Función escalón unidad
- Entrada sinusoidal. Estudiar especialmente el comportamiento a tiempos largos, es decir, para \(t \to \infty\).
Solución
- La respuesta de este proceso para una entrada en escalón unidadserá:
\[y (s) = G (s) \frac{1}{s} = K_p \frac{\tau s + 1}{\tau_p s + 1} \frac{1}{s}\]
Para obtener la respuesta en tiempo real hay que realizar la transformada inversa de Laplace:
using SymPy, Plots
@syms t::real Kp::real=>"K_p" T::real=>"tau" Tp::real=>"tau_p" s
G = Kp*(T*s+1)/(Tp*s+1)
f = 1/s
y = G*f
y_t = sympy.inverse_laplace_transform(y, s, t)
y_t = expand(y_t)
y_t = collect(y_t, Kp)
y_t = collect(y_t, exp(-t/Tp))\(K_{p} \left(\left(\frac{\tau \theta\left(t\right)}{\tau_{p}} - \theta\left(t\right)\right) e^{- \frac{t}{\tau_{p}}} + \theta\left(t\right)\right)\)
Es decir, la ecuación anterior se obtiene:
\[y (t) = K_p \left[ 1 - \left( 1 - \frac{\tau}{\tau_p} \right) e^{- \frac{t}{\tau_p}} \right]\]
Se obtiene una ecuación muy similar a la Ecuación 4.4. Evidentemente, el sentido físico e influencia sobre la respuesta del sistema de la ganancia del proceso (\(K_p\)) serán los mismos que para un sistema de primer orden. La influencia de la constante de tiempo del proceso (\(\tau_p\)) también será muy parecida.
La única diferencia es el término \(\tau / \tau_p\). Para ver el significado físico de la constante de tiempo \(\tau\) se puede representar la respuesta para diferentes valores \((K_p = \tau_p = 1)\):
taus = [0, .25, .5, 1, 2, 3]
p = plot()
for tau in taus
plot!(y_t(Kp=>1, Tp=>1, T=>tau), 0, 5, label=tau, lw=2)
end
pSe observa como \(\tau\), para valores inferiores a la unidad, actúa como un adelanto, avanza la respuesta del proceso de primer orden en un valor igual a \(\tau\).
En el caso de una entrada sinusoidal (\(\sin(\omega t)\)), la transformada de Laplace de la respuesta es:
\[y (s) = \frac{K_p (\tau s + 1)}{\tau_p s + 1} \frac{\omega}{s^2 + \omega^2}\]
Recurriendo de nuevo a SymPy para realizar la transformada inversa de Laplace:
w = symbols("omega", real=true)
f_t = sin(w*t)
f = sympy.laplace_transform(f_t, t, s, noconds=true)
y = G*f
y_t = sympy.inverse_laplace_transform(y, s, t)\(K_{p} \omega \left(\left(\frac{\left(\tau - \tau_{p}\right) \cos{\left(\omega t \right)}}{\omega^{2} \tau_{p}^{2} + 1} + \frac{\left(\omega^{2} \tau \tau_{p} + 1\right) \sin{\left(t \left|{\omega}\right| \right)}}{\left(\omega^{2} \tau_{p}^{2} + 1\right) \left|{\omega}\right|}\right) \theta\left(t\right) + \frac{\left(- \tau \tau_{p} + \tau_{p}^{2}\right) e^{- \frac{t \left(\omega^{2} \tau_{p}^{2} + 1\right)}{\omega^{2} \tau_{p}^{3} + \tau_{p}}} \theta\left(t\right)}{\omega^{2} \tau_{p}^{3} + \tau_{p}}\right)\)
y_t = cancel(y_t)
y_t = collect(y_t, exp(t/Tp))
y_t = collect(y_t, sin(w*t))
y_t = collect(y_t, cos(w*t))\(\frac{\left(- K_{p} \omega \tau \left|{\omega}\right| \theta\left(t\right) + K_{p} \omega \tau_{p} \left|{\omega}\right| \theta\left(t\right) + \left(K_{p} \omega^{3} \tau \tau_{p} \sin{\left(t \left|{\omega}\right| \right)} \theta\left(t\right) + K_{p} \omega \sin{\left(t \left|{\omega}\right| \right)} \theta\left(t\right) + \left(K_{p} \omega \tau \left|{\omega}\right| \theta\left(t\right) - K_{p} \omega \tau_{p} \left|{\omega}\right| \theta\left(t\right)\right) \cos{\left(\omega t \right)}\right) e^{\frac{t}{\tau_{p}}}\right) e^{- \frac{t}{\tau_{p}}}}{\omega^{2} \tau_{p}^{2} \left|{\omega}\right| + \left|{\omega}\right|}\)
Es decir:
\[y (t) = \frac{(K_p \omega \tau_p^2 - K_p \omega \tau \tau_p) \mathrm{e}^{- \frac{t}{\tau_p}}}{\tau_p (\omega^2 \tau^2_p + 1)} + \frac{K_p \omega^3 \tau \tau_p + K_p \omega}{\omega (\omega^2 \tau^2_p + 1)} \sin (\omega t) - \frac{K_p \omega \tau_p - K_p \omega \tau}{\omega^2 \tau^2_p + 1} \cos (\omega t)\]
Considerando la propiedad trigonométrica:
\[x \sin \alpha + y \cos \alpha = z \sin (\alpha + \varphi)\]
donde \(z^2 = x^2 + y^2\) y \(\tan \varphi = \frac{y}{x}\), se obtiene:
\[y (t) = \frac{(K_p \omega \tau_p^2 - K_p \omega \tau \tau_p) e^{- \frac{t}\tau_p}}{\tau_p (\omega^2 \tau^2_p + 1)} + \sqrt{\frac{K_p^2 \omega^2 \tau^2 + K^2_p}{\omega^2 \tau^2_p + 1}} \sin (\omega t + \varphi)\]
donde
\(\varphi = \mathrm{atan} \left( \frac{\omega^2 (K_p \omega \tau_p - K_p \omega \tau)^2}{(K_p \omega^3 \tau \tau_p + K_p \omega)^2} \right) = \mathrm{atan} \left( \frac{\omega^2 (\tau_p - \tau)^2}{(\omega^2 \tau \tau_p + 1)^2} \right)\).
Se obtiene una respuesta con una parte transitoria, que rápidamente se anula al aumentar el tiempo:
\[\mathrm{Transitorio} = \lim_{t \to 0} \frac{(K_p \omega \tau_p^2 - K_p \omega \tau \tau_p) e^{- \frac{t}{\tau_p}}}{\tau_p (\omega^2 \tau^2_p + 1)} = 0\]
La parte estacionaria de la respuesta es:
\[y (t) = K_p \sqrt{\frac{\omega^2 \tau^2 + 1}{\omega^2 \tau^2_p + 1}} \sin (\omega t + \varphi)\]
Si se compara la respuesta estacionaria respecto a la entrada se comprueba que tiene la misma frecuencia angular \(\omega\) y que tiene un desfase \(\varphi\), lo que significa que está retrasada. Este desfase depende de la frecuencia angular. La amplitud de la respuesta también depende de la frecuencia angular de la entrada.
Nuevamente la respuesta obtenida es muy similar a la de un proceso de primer orden. Por tanto, la influencia de la ganancia y constante de tiempo del proceso serán similares.
La constante de tiempo \(\tau\) nuevamente actúa como un adelanto, compensando en cierta medida el efecto de la constante de tiempo del proceso. Al aumentar su valor disminuye el desfase y la amplitud disminuye en menor medida que para un proceso de primer orden.
Problema 4.2
Dibujar la respuesta de un sistema de primer orden de constante de tiempo 0.5 y ganancia 1 para las siguientes entradas:
- Un impulso unidad
- Un pulso unidad de duración 5
- Un cambio sinusoidal de amplitud unidad y frecuencia \(\omega = 0.5\).
Problema 4.3 ★
Sea un sistema de primer orden de ganancia unidad y constante de tiempo 0.5. Inicialmente el sistema está en estado estacionario. Se introduce una entrada en rampa cuando el tiempo es igual a 0.
- Desarrollar una expresión que muestre los cambios en el proceso con el tiempo.
- ¿Cuál es la mínima y la máxima diferencia entre la entrada y la salida?
- Dibujar la entrada y la salida en función del tiempo.
Solución
- Hay varias estrategias para desarrollar una expresión que muestre los cambios en el proceso tras una rampa unidad (\(y(t)\)). Se puede resolver directamente la ecuación diferencial o se puede utilizar la función de transferencia y realizar la transformada inversa de Laplace.
Resolución de la ecuación diferencial
Un sistema de primer orden viene descrito por la siguiente ecuación diferencial (Ecuación 4.2):
\[\tau_p \frac{\mathrm{d} y}{\mathrm{d} t} + y = K_p f(t)\]
donde \(f (t)\) es la entrada al sistema, en nuestro caso una rampa unidad (\(f(t) = t\)). Sustituyendo las constantes y la función de entrada, se obtiene la ecuación diferencial a resolver:
\[0.5 \frac{\mathrm{d} y}{\mathrm{d} t} + y = t\]
Se puede resolver la ecuación diferencial utilizando SymPy:
using SymPy
@syms t::real y()
D = Differential(t)
eq = 1//2*D(y(t))+ y(t) ~ t\(y{\left(t \right)} + \frac{\frac{d}{d t} y{\left(t \right)}}{2} = t\)
Se ha supuesto que se están utilizando variables de desviación, lo que supone que \(y (t = 0) = 0\):
ics = Dict(y(0)=>0);Resolviendo la ecuación diferencial, se obtiene:
dsolve(eq, ics=ics)\(y{\left(t \right)} = t - \frac{1}{2} + \frac{e^{- 2 t}}{2}\)
Función de transferencia
La función de transferencia de este proceso es:
\[G = \frac{y (s)}{f (s)} = \frac{1}{0.5 s + 1}\]
y la entrada es una rampa unidad, cuya transformada de Laplace es:
\[f(s) = \frac{1}{s^2}\]
Por tanto la respuesta del proceso es:
\[y(s) = G f(s) = \frac{1}{0.5 s + 1} \frac{1}{s^2}\]
Para obtener la respuesta dependiente del tiempo hay que realizar la transformada inversa de Laplace:
\[y (t) =\mathcal{L}^{- 1} \left( \frac{1}{0.5 s + 1} \frac{1}{s^2} \right)\]
Para realizar la transformada inversa se puede utilizar la técnica de separar en fracciones simples o simplemente utilizar SymPy:
@syms s
y_t = sympy.inverse_laplace_transform(1/(1//2*s+1)*1/s^2, s, t)
expand(y_t)\(t \theta\left(t\right) - \frac{\theta\left(t\right)}{2} + \frac{e^{- 2 t} \theta\left(t\right)}{2}\)
Lógicamente, la respuesta obtenida es igual a la obtenida por el método anterior.
- La diferencia entre la entrada y la salida es:
\[y (t) - f (t) = \frac{e^{- 2 t}}{2} - \frac{1}{2}\]
La función exponencial es continua y decreciente, si se trata de exponentes negativos, para valores de \(t\) mayores que cero, como es el caso. Cuanto \(t = 0\), \(e^{- 2 t} = 1\). Si el tiempo tiende a infinito, \(e^{- 2 t} \underset{t \rightarrow \infty}{\rightarrow} 0\). Por tanto, la diferencia mínima se dará cuando el tiempo es igual a cero, en ese caso la diferencia es de 0. La diferencia máxima se produce cuando el tiempo tiende a infinito, la diferencia es \(y (t) - f (t) \underset{t \rightarrow \infty}{\rightarrow} - \frac{1}{2}\).
- El gráfico de la entrada y la salida en función del tiempo es:
Esta gráfica se ha dibujado a partir de la función respuesta obtenida en el apartado a). En caso de no querer obtener esa función se puede programar el problema con Xcos o Simulink y obtener el resultado mediante métodos numéricos.
A continuación se muestra el programa junto con el resultado obtenido:
using Plots, LaTeXStrings
plot(y_t, 0, 3, lw=2, label=L"y(t)", xlabel=L"t", legend=:topleft)
plot!(t, 0, 3, label=L"t", lw=2)Problema 4.4
En la figura siguiente se muestra un tanque agitado de mezcla, donde \(F_i\) son caudales volumétricos y \(c_i\) son concentraciones:
Asumiendo que todos los caudales de entrada y salida son constantes, pero que las concentraciones de las corrientes de entrada puede cambiar, contestar a las siguientes preguntas:
Demostrar que la composición de la corriente de salida tiene un comportamiento de un sistema de primer orden para cambios en la composición de entrada.
Encontrar las funciones de transferencia entre las composiciones de las corrientes de salida y las de entrada. Dibujar el diagrama de bloques correspondiente.
Definir la ganancia y la constante de tiempo de este sistema.
Problema 4.5
Repetir el problema anterior con el siguiente tanque de mezcla:
