3  Como abordar la dinámica de un sistema

Introducción

Este es un capítulo clave del temario, ya que nos va a proporcionar las herramientas con las que podremos afrontar con éxito los siguientes temas. En este capítulo confluyen varias de las asignaturas que habéis estado viendo en los últimos años:

• Matemáticas: Utilizaremos operadores matemáticos, números complejos, ecuaciones diferenciales, límites, métodos numéricos y propiedades trigonométricas. (Pista: los números complejos y las oscilaciones suelen ir de la mano, como ya habéis visto, por ejemplo, en electricidad con la corriente alterna.)

• Operaciones básicas: Utilizaremos balances macroscópicos con el matiz interesante de que los realizaremos en estado no estacionario. Es una diferencia importante respecto a cuando se utilizan para diseño, donde se suelen plantear las operaciones en estado estacionario.

El objetivo de este capítulo es poder plantear modelos matemáticos de operaciones básicas que nos permitan poder simular su comportamiento, su dinámica, en el tiempo. Seguiremos esta secuencia:

1. Detallar el funcionamiento de la operación: Cuanto más conozcamos la operación, mejor será nuestro modelo matemático y las simulaciones serán más fiables

2. Plantear el modelo matemático de la operación utilizando balances no estacionarios, es decir, se tratará de ecuaciones diferenciales

3. Resolución de las ecuaciones diferenciales:

   • Solución analítica (si la hay): Resolver ecuaciones diferenciales, aunque sean ordinarias y no muy complejas, (como es nuestro caso) no es sencillo. Para simplificar la tarea recurriremos a una herramienta matemática que facilita mucho su resolución: la transformada de Laplace.
   • Solución mediante métodos numéricos: Cuando no sea posible la obtención de soluciones analíticas, utilizaremos una herramienta de modernización y simulación informática que nos permitirá resolver numéricamente las ecuaciones diferenciales

Advertencia

En esta unidad se introducen muchos conceptos que, aunque no son difíciles una vez se dominan, requieren de práctica y un cierto pensamiento abstracto.

Hay que dedicar tiempo suficiente para practicar y asimilar las herramientas que se presentan, ya que su uso es clave en los temas siguientes.

La unidad se compone de los siguientes apartados:

Ejemplo de dinámica: Se explican los pasos a seguir para desarrollar el modelo matemático de una operación básica. Consideraremos que el modelo está terminado al plantear la ecuación o ecuaciones diferenciales.

Transformada de Laplace: Introducción de la transformada de Laplace como una herramienta útil para resolver ecuaciones diferenciales. Se introducirá la definición del operador y sus principales propiedades.

Función de transferencia: La función de transferencia nos permitirá trabajar con los modelos matemáticos sin tener que estar realizan la transformada de Laplace de manera continua.

Funciones singulares]: De cara a plantear las simulaciones es bueno conocer unas cuantas funciones matemáticas que, aunque son idealizaciones, son representativas de comportamientos del mundo real.

Inversión de la transformada de Laplace: Realizar la transformada de Laplace es sencillo, en cambio, realizar la operación inversa es mucho más complejo. Es una situación parecida a derivar e integrar. Las dos operaciones son importantes, pero resulta más complejo integrar que derivar. El problema es que es un paso imprescindible para obtener la solución de la ecuación diferencial planteada.

Expansión en fracciones simples: Técnicas para realizar la transformada inversa de Laplace.

Tablas de transformadas de Laplace: Se presentan las transformadas de Laplace y transformadas inversas de las funciones más habituales.

Tip

Aunque resolver las ecuaciones diferenciales, realizar las transformadas de Laplace y las transformadas inversas no es complicado, es un tarea tediosa.

Utilizaremos de manera habitual SymPy, un sistema algebraico computacional (CAS), que nos eliminará gran parte del trabajo. Es recomendable saberlo utilizar, pero no es imprescindible para superar la asignatura.

3.1 Un ejemplo de dinámica de un sistema. ¿Qué se desea conocer?

Como ejemplo se considera un sistema dinámico de nivel. En este caso se pretende encontrar el comportamiento con el tiempo de un depósito en función del caudal de entrada. Para ello habrá que buscar un modelo matemático que describa dicho comportamiento. Se recomienda seguir el siguiente procedimiento:

1. Se realiza el diagrama de flujo del sistema: En ese diagrama se marcan todas las variables implicadas. En este caso se trata de un depósito que se descarga por gravedad. Las variables implicadas son:

   • Sección del depósito: \(A(t)\)

   • Resistencia de la tubería al paso del fluido: \(R\)

   • Caudales volumétricos de entrada y salida: \(q_1 (t)\) y \(q_2 (t)\)

   • Nivel del depósito: \(h(t)\)

   • Densidad del fluido: \(\rho (t)\)

Figura 3.1: El nivel del depósito depende de los caudales de entrada y salida del mismo.

2. Planteamiento del modelo matemático: Un modelo matemático es un conjunto de ecuaciones que relacionan entre sí las variables del sistema. Se basan en ecuaciones de estado, leyes de equilibrio, ecuaciones cinéticas y balances de materia, energía y cantidad de movimiento.

   Para el ejemplo, al aplicar el Balance Macroscópico de Materia:

   \[ \frac{d A (t) h (t) \rho (t)}{d t} = \rho (t) q_1 (t) - \rho(t) q_2 (t) \tag{3.1}\]

   En el momento de plantear el modelo matemático hay que explicitar todas las hipótesis o simplificaciones realizadas. En este caso se han realizado las siguientes suposiciones:

   • La densidad (\(\rho\)) y el area del depósito (\(A\)) son constantes e independientes del tiempo.

   • El caudal de salida del depósito depende del nivel y de la resistencia de la tubería al paso del fluido de este modo: \(q_2 (t) = \frac{h (t)}{R}\).

   Teniendo en cuenta las hipótesis anteriores la Ecuación 3.1 queda como:

   \[A \frac{d h}{d t} = q_1 (t) - \frac{h (t)}{R} \tag{3.2}\]

3. Definición de las variables de desviación: Debido a que lo realmente importante en control es conocer cuánto se ha desplazado el sistema respecto al estado estacionario se definen unas variables de desviación.

   Para el ejemplo, el balance macroscópico de materia en estado estacionario es:

   \[0 = q_{1 e} - \frac{h_e}{R} \tag{3.3}\]

   donde \(_e\) indica estado estacionario.

   Restando ambos balances macroscópicos de materia (Ecuación 3.1 y Ecuación 3.3):

   \[A \frac{d (h (t) - h_e)}{d t} = q_1 (t) - q_{1 e} - \frac{h(t) - h_e}{R} \tag{3.4}\]

   Se definen las siguientes variables de desviación:

   \[\begin{aligned} H (t) &\equiv h (t) - h_e & \\ Q_1 (t) &\equiv q_1 (t) - q_{1 e} \end{aligned}\]

   Las mayúsculas indican que se tratan de variables de desviación.

   Sustituyendo en la Ecuación 3.4 se obtiene:

   \[A \frac{d H (t)}{d t} = Q_1 (t) - \frac{H (t)}{R} \tag{3.5}\]

   Esta ecuación es el modelo matemático que representa la respuesta dinámica del tanque a cambios en el caudal de entrada. Resolviendo esta ecuación diferencial se puede saber cómo varía el nivel con el tiempo según cambia el caudal de entrada.

3.2 La transformada de Laplace como herramienta útil

La transformada de Laplace \(\bar{f} (s)\) de una función \(f (t)\) se define como:

\[\mathcal{L} [f (t)] \equiv \bar{f} (s) = \int_0^{\infty} f (t) e^{- s t} \mathrm{d}t\]

donde \(t \in \mathbb{R}\) y \(s \in \mathbb{C}\).

El uso de transformadas de Laplace ofrece un método simple y elegante de resolver ecuaciones diferenciales como las que se obtienen en los modelos matemáticos de los procesos alimentarios.

Entre las diferentes propiedades de las transformadas de Laplace cabe destacar:

1. Es un operador lineal:

   \[\mathcal{L} [a_1 f_1 (t) + a_2 f_2 (t)] = a_1 \mathcal{L} [f_1 (t)] + a_2 \mathcal{L} [f_2 (t)]\]

2. La transformada de una derivada es:

   \[\mathcal{L} \left[ \frac{d f (t)}{d t} \right] = s \bar{f} (s) - f (0)\]

   Es importante resaltar que una ecuación diferencial ordinaria de primer orden pasa a ser una ecuación lineal de primer grado.

   La transformada de la segunda derivada es:

   \[\mathcal{L} \left[ \frac{d^2 f (t)}{d t^2} \right] = s^2 \bar{f} (s) - s f (0) - f' (0)\]

   Generalizando:

   \[\mathcal{L} \left[ \frac{d^{(n)} f (t)}{d t^n} \right] = s^n \bar{f} (s) - s^{n - 1} f (0) - s^{n - 2} f' (0) - \ldots - s f^{(n - 2)} (0) - f^{(n - 1)} (0)\]

   Esta propiedad permite resolver de manera sencilla ecuaciones diferenciales ya que las convierte en ecuaciones algebraicas. El Problema 2.3 muestra cómo se puede utilizar la transformada de Laplace para resolver una ecuación diferencial.

3. La transformada de Laplace de una integral es:

   \[\mathcal{L} \left[ \int_0^t f (t) d t \right] = \frac{1}{s} \bar{f} (s)\]

4. Translación de la transformada:

   \[\mathcal{L} [e^{- \alpha t} f (t)] = \bar{f} (\alpha + s)\]

   Translación de la función:

   \[\mathcal{L} [f (t - t_0) U(t-t_0)] = e^{- t_0 s} \bar{f}(s) \tag{3.6}\]

5. Teorema del valor final:

   \[ \mathcal{L} \left[ \underset{t \rightarrow \infty}{\lim} f (t) \right] = \underset{s \rightarrow 0}{\lim} [s \bar{f}(s)] \tag{3.7}\]

Ejemplo

La aplicación de la transformada de Laplace a una función es sencilla disponiendo de las tablas de transformadas de Laplace y del conocimiento de las propiedades anteriores. En el Prob. 3.1 se muestra como se aplica.

3.3 La función de transferencia. Álgebra de funciones de transferencia

Para un proceso sencillo, como el del nivel del depósito, se puede plantear un esquema sencillo que describa en cierta medida al sistema:

Figura 3.2: Representación simplificada de un proceso.

Para el ejemplo del nivel del depósito, \(f (t) = Q_1 (t)\) y \(y (t) = H (t)\). Es decir, en este caso el depósito es el proceso, la salida es la variación del nivel del depósito y la entrada al proceso es el caudal de entrada.

En el caso de trabajar utilizando las transformadas de Laplace de las funciones de entrada y salida, se puede representar la dinámica del proceso mediante el uso de la función de transferencia. La función de transferencia \(G (s)\) liga la entrada y la salida del sistema:

\[G (s) = \frac{\mathcal{L} [y (t)]}{\mathcal{L} [f (t)]} = \frac{\bar{y} (s)}{\bar{f} (t)}\]

donde \(\bar{y} (s)\) es la transformada de Laplace de la respuesta del proceso definida utilizando variables de desviación y \(\bar{f} (s)\) es la transformada de la función de desviación de la entrada.

Figura 3.3: La función de transferencia del proceso es \(G(s)\), \(\bar{f}(s)\) y \(\bar{y}(s)\) son la transformada de Laplace de las funciones de entrada y salida del proceso.

Para obtener la función de transferencia del ejemplo se realiza la transformada de Laplace de la ecuación diferencial Ecuación 3.5:

\[\mathcal{L} \left[ A \frac{d H (t)}{d t} \right] =\mathcal{L} \left[ Q_1 (t) - \frac{H (t)}{R} \right]\]

Se realizan las transformadas considerando que el operador es lineal y conociendo la transformada de una derivada:

\[A s \bar{H} (s) = \overline{Q_1} (s) - \frac{\bar{H} (s)}{R}\]

Se supone que para \(t = 0\) el sistema está todavía en estado estacionario, las variables de desviación son nulas.

Operando la ecuación anterior para encontrar la función de transferencia:

\[G (s) = \frac{\bar{H} (s)}{\overline{Q_1} (s)} = \frac{R}{R A s + 1}\]

Este es la función de transferencia típica de un sistema de primer orden.

Para los procesos más complicados, con diagramas de bloques más complejos, se recurre al álgebra de funciones de transferencia.

En la situación en la que se presente un conjunto de procesos en paralelo, que se representan como una serie de bloques en paralelo, como los de la figura siguiente:

Figura 3.4: Reducción de tres bloques en paralelo.

Se puede encontrar un bloque equivalente que reúna las tres funciones de transferencia \(G_1\), \(G_2\) y \(G_3\) en una sola función de transferencia \(G\):

Figura 3.5: Reducción de un grupo de \(n\) bloques en serie.

Buscar un bloque equivalente tiene la ventaja de que solo será necesario trabajar con una función de transferencia en lugar de con tres, si no se necesita conocer las funciones \(X_1\), \(X_2\) y \(X_3\).

En este caso obtener la función de transferencia equivalente es una tarea sencilla:

\[G (s) = \frac{Y (s)}{X (s)} = \frac{Y_1 (s) + Y_2 (s) + Y_3 (s)}{X (s)} = G_1 (s) + G_2 (s) + G_3 (s)\]

Generalizando:

\[G (s) = \sum_i G_i (s)\]

Para un diagrama de bloques en serie:

\[G (s) = \frac{X_n}{X_0} = \frac{X_1}{X_0} \frac{X_2}{X_1} \cdots \frac{X_n}{X_{n - 1}} = G_1 G_2 \ldots G_n\]

Generalizando:

\[G (s) = \prod_i G_i (s) \tag{3.8}\]

Ejemplo

El Problema 3.6 muestra un ejemplo de sistema en serie, se trata de dos depósitos de almacenamiento en serie. La curiosidad es que se plantea el sistema con interacción y el sistema de tanques sin interacción, para el que no es aplicable de manera directa la Ecuación 3.8.

3.4 Transformadas de algunas funciones singulares

A continuación se muestran las transformadas de algunas funciones con las que se trabajará frecuentemente más adelante ya que pueden ser asimiladas como las perturbaciones más frecuentes.

Si no se dice lo contrario todas estas funciones se definen para que su valor sea nulo a tiempo menor que cero.

3.4.1 Función escalón

Es una función cuyo valor para tiempos menores que cero es nulo y que alcanza el valor \(M\) para tiempo mayor que 0:

Figura 3.6: Función escalón de altura \(M\).

Esta función se define como:

\[f (t) = \left\{\begin{array}{ll} 0 & t < 0\\ M & t > 0 \end{array}\right.\]

La transformada de Laplace de esta función es:

\[\mathcal{L} [f (t)] = \frac{M}{s}\]

Si \(M\) es igual a 1 se habla de la función escalón unidad, \(U(t)\), o de la función de Heaviside, \(\theta(t)\).

Encontrar la transformada inversa de Laplace es muy sencillo utilizando SymPy:

# Definición de las variables que vamos a utilizar
@syms M t s

# Definimos un escalón unidad de altura M
f = M*Heaviside(t)

# Calculamos la transformada inversa de Laplace
sympy.laplace_transform(f, t, s, noconds=true)

\(\frac{M}{s}\)

En el cálculo anterior, se ha utilizado la función de Heaviside como función escalón unidad. También se puede realizar el cálculo definiendo la función escalón como una constante \(M\) y asumiendo de manera implícita que ese valor es solo para \(t > 0\):

# Definimos un escalón unidad de altura M
f = M

# Calculamos la transformada inversa de Laplace
sympy.laplace_transform(f, t, s, noconds=true)

\(\frac{M}{s}\)

Aunque la mayor parte de las veces las funciones de entrada sucede a tiempo 0, no es raro que aparezcan desplazadas en el tiempo, que sufran un retraso \(t_0\).

Figura 3.7: Función pulso unidad retrasada en el tiempo un valor \(t_0\).

En el caso de que la función tenga un retraso \(t_0\):

\[f (t) = \left\{\begin{array}{ll} 0 & t < t_0\\ M & t > t_0 \end{array}\right.\]

O lo que es lo mismo:

\[f (t - t_0) = \left\{\begin{array}{ll} 0 & t < 0\\ M & t > 0 \end{array}\right.\]

Por tanto, aplicando la propiedad de la transformada de Laplace, (ecuación {eq}translacion), la transformada de Laplace será:

\[\mathcal{L} [f (t - t_0)] = \frac{M}{s} \mathrm{e}^{- s t_0}\]

3.4.2 Función pulso

Se trata de una función pulso con área \(A = M t_0\). A continuación, en la figura de la derecha se muestra un pulso de altura unidad, \(M=1\). El pulso (línea roja) es la suma de un pulso unidad (línea verde) al que se le suma un pulso de altura -1 retrasado un tiempo igual a la anchura del pulso, es decir, \(t_0\):

Figura 3.8: Un pulso, figura de la derecha, es la resta de dos funciones escalón de igual altura. El escalón negativo (azul) anula al escalón verde y forma el pulso (naranja).

La función pulso se define como:

\[f (t) = \left\{\begin{array}{cc} 0 & t < 0\\ M & 0 < t < t_0\\ 0 & t > t_0 \end{array}\right.\]

Utilizando la definición del escalón unidad también se puede escribir como:

\[f (t) = M [U (t) - U (t - t_0)]\]

Por tanto, la transformada de Laplace será:

\[\mathcal{L} [f (t)] = \bar{f} (s) = M \left( \frac{1}{s} - \frac{\mathrm{e}^{- s t_0}}{s} \right) = \frac{M}{s} (1 - \mathrm{e}^{- s t_0})\]

3.4.3 Función impulso

Se trata de un pulso de área \(A\) tal que \(M \rightarrow \infty\) y \(t_0 \rightarrow 0\):

Figura 3.9: Función impulso de área \(A\).

La transformada de Laplace de esta función es:

\[\mathcal{L} [f (t)] = \bar{f} (s) = A\]

En el caso particular de que el área sea 1 se habla de la función delta de Dirac \(\delta (t)\).

A continuación se calcula la transformada de Laplace de un impulso de altura A, \(f(t) = A \delta(t)\):

@syms A

sympy.laplace_transform(A*sympy.DiracDelta(t), t, s, noconds=true)

\(A\)

Se puede comprobar fácilmente que el impulso es la derivada de la función escalón.

Ejemplo

En el Problema 3.4 se puede comprobar las diferencias y similitudes en la respuesta de un proceso a una entrada en escalón y en impulso.

3.4.4 Función rampa

Se trata de una función lineal de pendiente \(M\):

Figura 3.10: Función rampa de pendiente \(M\).

Esta función se define como:

\[f (t) = \left\{\begin{array}{ll} 0 & t < 0\\ M t & t > 0 \end{array}\right. = M t U (t)\]

La transformada de Laplace es:

\[\mathcal{L} [M t U (t)] = \frac{M}{s^2}\]

sympy.laplace_transform(M*t, t, s, noconds=true)

\(\frac{M}{s^{2}}\)

3.4.5 Funciones trigonométricas

La función seno es:

Figura 3.11: Función seno de frecuencia \(\omega\).

Se define la función como:

\[f (t) = \left\{\begin{array}{cl} 0 & t < 0\\ M \sin (\omega t) & t > 0 \end{array}\right. = M U (t) \sin (\omega t)\]

donde \(M\) es la amplitud y \(\omega\) es la frecuencia angular, expresada normalmente como rad/s.

La transformada de Laplace de la función seno es:

\[\mathcal{L} [M \sin (\omega t)] = \frac{M \omega}{s^2 + \omega^2}\]

y la de la función coseno:

\[\mathcal{L} [M \cos (\omega t)] = \frac{M s}{s^2 + \omega^2}\]

@syms w=>"omega"
sympy.laplace_transform(M*sin(w*t), t, s, noconds=true)

\(\frac{M \omega}{\omega^{2} + s^{2}}\)

3.5 Inversión de transformadas. De vuelta al tiempo real

Continuando con el ejemplo se estudiará la salida del sistema para una entrada de tipo escalón unidad:

\[\bar{H} (s) = \left( \frac{R}{R A s + 1} \right) \overline{Q_1} (s) = \left( \frac{R}{R A s + 1} \right) \frac{1}{s}\]

Mediante el operador transformada inversa de Laplace (\(\mathcal{L}^{- 1}\)) se obtiene la salida en tiempo real. Para ello hay que descomponer la función a invertir en partes asimilables a las que se encuentran en las tablas de transformadas de Laplace (Sección 3.7):

\[\bar{H} (s) = \left( \frac{R}{R A s + 1} \right) \frac{1}{s} = \frac{a}{s} + \frac{b}{s + \frac{1}{A R}}\]

Donde \(a\) y \(b\) son dos variables a determinar. Obviamente, \(a = R\) y \(b = - R\). Por tanto,

\[H (t) = \left( R - R \mathrm{e}^{- \frac{t}{R A}} \right) U (t) = R U (t) \left( 1 - \mathrm{e}^{- \frac{t}{R A}} \right) = R U (t) \left( 1 - \mathrm{e}^{- \frac{t}{\tau}} \right)\]

donde \(\tau = R A\) es la constante de tiempo y tiene dimensiones de tiempo.

También se puede realizar la transformada inversa de Laplace de una manera muy simple utilizando SymPy:

using SymPy, Plots, LaTeXStrings

# Definición de los símbolos necesarios
@syms R::real A::real t::real
@syms s

# Definición de la función de transferencia, G(s)
G = R/(R*A*s+1)

# Definición de la función de entrada, f(s), 
# en este caso un escalón unidad
f = 1/s

# Cálculo de la respuesta, H(t)
sympy.inverse_laplace_transform(G*f, s, t)

\(R \left(\theta\left(t\right) - e^{- \frac{t}{A R}} \theta\left(t\right)\right)\)

Cuanto mayor es \(\tau\) más lenta es la respuesta, más tarda el sistema en alcanzar el estado estacionario. Se comprueba que cuanto menor es la sección del tanque más rápida es la respuesta. Si \(\tau\) es grande se dice que el sistema presenta una gran inercia.

Taus = [1//2, 1, 2, 4, 8 ]
f = 1/s
Yt = []
fig_inv = plot(lw=2,
    xlabel=L"t", ylabel=L"y(t)", legend=:bottomright)

for τ in Taus
    G = 1/(τ*s+1)
    y = G*f
    yt = sympy.inverse_laplace_transform(y, s, t)
    push!(Yt, yt)
end

for i in 1:5
    tau = float(Taus[i])
    plot!(Yt[i], -0.5, 10, lw=2, label=L"\tau_%$i = %$tau")
end

display(fig_inv)

Figura 3.12: Influencia en la respuesta de la constante de tiempo.

Ejemplo

La técnica propuesta en este capítulo para obtener modelos matemáticos se puede utilizar para modelos de mayor complejidad, como el que se obtiene en la resolución del Problema 3.7.

3.6 Expansión en fracciones parciales

Advertencia

En Man (2009) hay un método más sencillo que no requiere realizar derivadas.

Un método de inversión de transformadas de Laplace es la expansión en fracciones parciales. Para invertir la función \(\bar{f} (s)\) se reordena como una suma de funciones simples:

\[\bar{f} (s) = \overline{f_1} (s) + \overline{f_2} (s) + \ldots + \overline{f_s} (s)\]

Una vez realizada la descomposición se realiza la inversión de transformadas:

\[f (t) =\mathcal{L}^{- 1} [\overline{f_1} (s)] +\mathcal{L}^{- 1} [\overline{f_2} (s)] + \ldots +\mathcal{L}^{- 1} [\overline{f_s} (s)]\]

Normalmente las funciones a invertir en control de procesos aparecen como fracciones de polinomios en s del tipo:

\[\bar{f} (s) = \frac{\bar{z} (s)}{\bar{p} (s)}\]

donde z(s) es un polinomio de orden m y p(s) es un polinomio de orden n.

Para realizar la inversión de la transformada de Laplace se factoriza el denominador:

\[\bar{f} (s) = \frac{z (s)}{(s - p_1) (s - p_2) \ldots (s - p_n)}\]

donde \(p_i\) son las raíces (ceros) del polinomio p(s).

Si todos los \(p_i\) son diferentes, se puede expresar f(s) como una suma de n términos:

\[\bar{f} (s) = \frac{A}{s - p_1} + \frac{B}{s - p_2} + \frac{C}{s - p_3} + \cdots + \frac{W}{s - p_n}\]

Los numeradores de la ecuación anterior se evalúan de la siguiente manera:

\[\begin{aligned} & A = \lim_{s \to p_1} [(s - p_1) \bar{f} (s)] & \\ & B = \lim_{s \to p_2} [(s - p_2) \bar{f} (s)] & \\ & \vdots & \\ & W = \lim_{s \to p_n} [(s - p_n) \bar{f} (s)] & \end{aligned}\]

Si existen raíces del denominador de f(s) repetidas, de nuevo se expresa f(s) como un producto de fracciones simples. Por ejemplo, en el caso de que una raíz se repita dos veces:

\[\bar{f} (s) = \frac{z (s)}{(s - p_1)^2 (s - p_3) (s - p_4) \ldots (s - p_n)}\]

la descomposición en fracciones simples es:

\[\bar{f} (s) = \frac{A}{(s - p_1)^2} + \frac{B}{s - p_1} + \frac{C}{s - p_3} + \cdots + \frac{W}{s - p_n} \tag{3.9}\]

Si la raíz se repite 3 veces (orden 3) la expansión sería así:

\[\begin{aligned} & \bar{f} (s) = \frac{z (s)}{(s - p_1)^3 (s - p_4) (s - p_5) \ldots (s - p_n)} & \nonumber\\ & \bar{f} (s) = \frac{A}{(s - p_1)^3} + \frac{B}{(s - p_1)^2} + \frac{C}{s - p_1} + \frac{D}{s - p_4} + \cdots + \frac{W}{s - p_n} & \end{aligned} \tag{3.10}\]

Los numeradores de la Ecuación 3.9 se calculan de la siguiente manera:

\[\begin{aligned} & A = \lim_{s \to p_1} [(s - p_1)^2 \bar{f} (s)] & \\ & B = \lim_{s \to p_1} \left\{ \frac{\mathrm{d}}{\mathrm{d}s} [(s - p_1)^2 \bar{f} (s)] \right\} & \\ & C = \lim_{s \to p_3} [(s - p_3) \bar{f} (s)] & \\ & \vdots & \end{aligned}\]

Para encontrar el numerador C de la Ecuación 3.10 se debe tomar la segunda derivada. Generalizando al término \(A_j\) de una raíz de orden N en \(p_1\):

\[A_j = \lim_{s \to p_1} \left\{ \frac{\mathrm{d}^{m - 1}}{\mathrm{d}s^{m - 1}} [(s - p_1)^N \bar{f} (s)] \right\} \frac{1}{(m - 1) !}\]

Aunque su utilización es un tanto tediosa existen casos, como el apartado d) del problema 3.2, que muy difícilmente se puede resolver sin utilizar la técnica de expansión en fracciones parciales. Esta técnica tiene la ventaja de ser muy sistemática.

La expansión en fracciones simples también se puede realizar con SymPy utilizando la instrucción apart:

using SymPy
@syms s t
apart(3*s+3/(2*s^2+s+1))

\(3 s + \frac{3}{2 s^{2} + s + 1}\)

3.7 Transformadas de Laplace

\(\mathbf{f (t), t > 0}\)	\(\mathbf{\bar{f}(s)}\)
Impulso unidad, \(\delta(t_0)\)	1
Pulso unidad, \(\delta_A(t)\)	\(\frac{1}{A} \frac{1 - e^{- sA}}{s}\)
Escalón unidad	\(\frac{1}{s}\)
Rampa, \(f(t) = t\)	\(\frac{1}{s^2}\)
\(t^2\)	\(\frac{2!}{s^3}\)
\(t^n\)	\(\frac{n!}{s^{n + 1}}\)
\(e^{- a t}\)	\(\frac{1}{s + a}\)
\(t^n e^{- a t}\)	\(\frac{n!}{(s + a)^{n + 1}}\)
\(\sin (\omega t)\)	\(\frac{\omega}{s^2 + \omega^2}\)
\(\cos (\omega t)\)	\(\frac{s}{s^2 + \omega^2}\)
\(\sinh (\omega t)\)	\(\frac{\omega}{s^2 - \omega^2}\)
\(\cosh (\omega t)\)	\(\frac{s}{s^2 - \omega^2}\)
\(e^{- a t} \sin (\omega t)\)	\(\frac{\omega}{(s + a)^2 + \omega^2}\)
\(e^{- a t} \cos (\omega t)\)	\(\frac{s + a}{(s + a)^2 + \omega^2}\)

Tabla 3.1: Transformadas de Laplace de funciones seleccionadas.

\(\mathbf{\bar{f} (s)}\)	\(\mathbf{f (t)}\)
\(\frac{1}{(s + a) (s + b)}\)	\(\frac{e^{- a t} - e^{- b t}}{b - a}\)
\(\frac{1}{(s + a) (s + b) (s + c)}\)	\(\frac{e^{- a t}}{(b - a) (c - a)} + \frac{e^{- b t}}{(c - b) (a - b)} + \frac{e^{- a t}}{(a - c) (b - c)}\)
\(\frac{s + a}{(s + b) (s + c)}\)	\(\frac{1}{c - b} [(a - b) e^{- b t} -(a - c) e^{- c t}]\)
\(\frac{a}{(s + b)^2}\)	\(a t e^{- b t}\)
\(\frac{a}{(s + b)^3}\)	\(\frac{a}{2} t^2 e^{- b t}\)
\(\frac{a}{(s + b)^{n + 1}}\)	\(\frac{a}{n!} t^n b^{- b t}\)
\(\frac{1}{s (a s + 1)}\)	\(1 - e^{- t / a}\)
\(\frac{1}{s (a s + 1)^2}\)	\(1 - \frac{a + t}{a} e^{- t / a}\)
\(\frac{\omega^2}{s (s^2 + 2 \zeta \omega s + \omega^2)}\)	\(1 + \frac{e^{-\zeta \omega t}}{ \sqrt[]{1 - \zeta^2}} \sin \left( \omega \sqrt[]{1- \zeta^2} t - \varphi \right)\) donde \(\cos \varphi = - \zeta\)
\(\frac{s}{(1 + a s) (s^2 + \omega^2)}\)	\(- \frac{1}{1 + a^2 \omega^2}e^{- t / a} + \frac{1}{\sqrt[]{1 + a^2 \omega^2}} \cos (\omega t +\varphi)\) donde \(\varphi = \arctan (a \omega)\)
\(\frac{s}{(s^2 + \omega^2)^2}\)	\(\frac{1}{2 \omega} \sin (\omega t)\)
\(\frac{1}{(s + a) [(s + b)^2 + \omega^2]}\)	\(\frac{e^{- a t}}{(a - b)^2+ \omega^2} + \frac{e^{- b t} \sin (\omega t + \varphi)}{\omega\sqrt[]{(a - b)^2 + \omega^2}}\) donde \(\varphi = \arctan \left(\frac{\omega}{a - b} \right)\)

Tabla 3.2: Transformada inversa de Laplace de funciones seleccionadas.

3.8 Problemas

Problema 3.1 ★

Hallar la transformada de Laplace de las siguientes funciones:

1. \(f (t) = \mathrm{e}^{- 2 t} \cos (t)\)

2. \(f (t) = \mathrm{e}^{- 4 t} + \theta(t-2) \sin (t - 2) + t^2 \mathrm{e}^{- 2 t}\)

3. \(f (t) = 2 \mathrm{e}^{- t} \cos (10 t) - t^4 + 6 \mathrm{e}^{- (t - 10)} \theta(t-10)\)

Solución

1. Aplicando la propiedad 4. Traslación de la transformada:

\[\mathcal{L} (\mathrm{e}^{- 2 t} \cos t) = \frac{s + 2}{s^2 + 4 s + 5}\]

También se puede realizar la transformada de Laplace utilizando SymPy. En primer lugar cargamos la librería:

using SymPy

A continuación definimos las variables \(s\) y \(t\). Las definimos por separado para indicar que \(t\) es una variable real, mientras que \(s\) es una variable compleja:

@syms s
@syms t::real

(t,)

Finalmente calculamos la transformada de Laplace:

sympy.laplace_transform(exp(-2*t)*cos(t), t, s)

((s + 2)/((s + 2)^2 + 1), -2, True)

SymPy ofrece más información de la que necesitamos, ya que solo necesitamos la función transformada. Podemos obtenerla añadiendo el parámetro noconds:

sympy.laplace_transform(exp(-2*t)*cos(t), t, s, noconds=true)

\(\frac{s + 2}{\left(s + 2\right)^{2} + 1}\)

Puede ser conveniente definir un función para simplificar la realización de la transformada de Laplace, ya que siempre vamos a utilizar las variables \(s\) y \(t\):

L(f) = sympy.laplace_transform(f, t, s, noconds=true)

L(exp(-2*t)*cos(t))

\(\frac{s + 2}{\left(s + 2\right)^{2} + 1}\)

2. El operador transformada de Laplace es lineal:

\[\mathcal{L} (\mathrm{e}^{- 4 t} + \sin (t - 2) + t^2 \mathrm{e}^{- 2 t}) =\mathcal{L} (\mathrm{e}^{- 4 t}) +\mathcal{L} (\sin (t - 2)) +\mathcal{L} (t^2 \mathrm{e}^{- 2 t})\]

La transformada de Laplace del primer sumando es directa a partir de las tablas:

\[\mathcal{L} [\mathrm{e}^{- 4 t}] = \frac{1}{s + 4}\]

Para calcular el segundo de los sumandos es necesario aplicar, de nuevo, la propiedad 4.:

\[\mathcal{L} [\sin (t - 2)] = \mathrm{e}^{- 2 s} \mathcal{L} [\sin (t)] = \mathrm{e}^{- 2 s} \frac{1}{s^2 + 1}\]

Por último, el tercer sumando se obtiene de la consulta de las tablas:

\[\mathcal{L} [t^2 \mathrm{e}^{- 2 t}] = \frac{2!}{(s + 2)^{2 + 1}} = \frac{2}{(s + 2)^3}\]

Por tanto,

\[\bar{f} (s) = \frac{1}{s + 4} + \mathrm{e}^{- 2 s} \frac{1}{s^2 + 1} + \frac{2}{(s + 2)^3}\]

Con SymPy obtenemos el mismo resultado:

L(exp(-4*t) + sin(t-2)*Heaviside(t-2) + t^2*exp(-2*t))

\(\frac{e^{- 2 s}}{s^{2} + 1} + \frac{1}{s + 4} + \frac{2}{\left(s + 2\right)^{3}}\)

3. De nuevo, aplicando la linealidad del operador transformada:

\[\begin{aligned} & \mathcal{L} [2 \mathrm{e}^{- t} \cos (10 t)] = 2 \frac{s + 1}{(s + 1)^2 + 10^2} & \\ & \mathcal{L} [t^4] = \frac{4!}{s^5} = \frac{24}{s^5} & \\ & \mathcal{L} [6 \mathrm{e}^{- (t - 10)}] = 6 \mathrm{e}^{- 10 s} \frac{1}{s + 1} & \end{aligned}\]

Por tanto,

\[\bar{f} (s) = 2 \frac{s + 1}{(s + 1)^2 + 10^2} - \frac{24}{s^5} + 6 \mathrm{e}^{- 10 s} \frac{1}{s + 1}\]

Utilizando SymPy:

L(2*exp(-t)*cos(10*t)) - L(t^4) + L(6*exp(t-10)*Heaviside(t-10))

\(\frac{2 \left(s + 1\right)}{\left(s + 1\right)^{2} + 100} + \frac{6 e^{- 10 s}}{s - 1} - \frac{24}{s^{5}}\)

Problema 3.2 ★

Hallar la transformada inversa de:

1. \(f (s) = \frac{1}{\frac{s}{3} + 1}\)

2. \(f (s) = \frac{s}{(s + 1) (s^2 + 1)}\)

3. \(f (s) = \frac{3 s + 4}{s^3 (s + 2)}\)

4. \(f (s) = \frac{s^2}{(s^2 + 1)^2}\)

Solución

1. Para realizar la transformada inversa de Laplace se puede consultar las tablas, donde se encuentra que:

\[\mathcal{L} (\mathrm{e}^{-at}) = \frac{1}{s + a}\]

Por tanto, para resolver el ejercicio se operará de la siguiente manera:

\[f (s) = \frac{1}{\frac{s}{3} + 1} = \frac{3}{s + 3}\]

Por tanto, realizando la transformada inversa de Laplace y considerando que el operador es asociativo:

\[f (t) =\mathcal{L}^{- 1} \left( \frac{3}{s + 3} \right) = 3 \mathrm{e}^{- 3 t}\]

También se puede resolver el problema utilizando SymPy, como siempre empezaremos cargando la librería y defiendo las variables \(t\) y \(s\):

using SymPy
@syms s
@syms t::real;

La función para realizar la transformada inversa de Laplace es inverse_laplace_transform y la sintaxis es sencilla y similar a laplace_transform:

sympy.inverse_laplace_transform(1/(s/3 + 1), s, t)

\(3 e^{- 3 t} \theta\left(t\right)\)

Al igual que en el problema 3.01 definimos una función L() para que hacer que fuera más sencillo el calcular la transformada de Laplace, podemos hacer lo mismo para la transformada inversa de Laplace. Definimos una nueva función llamada invL, su funcionamiento es trivial:

invL(f) = sympy.inverse_laplace_transform(f, s, t)

invL(1/(s/3 + 1))

\(3 e^{- 3 t} \theta\left(t\right)\)

2. Para poder realizar la transformada inversa de Laplace es necesario realizar la descomposición de la función en fracciones simples:

\[f (s) = \frac{s}{(s + 1) (s^2 + 1)} = \frac{A}{s + 1} + \frac{B (s + 1)}{s^2 + 1}\]

Realizando la suma y simplificando el denominador se obtiene:

\[s = A (s^2 + 1) + B (s^2 + 2 s + 1)\]

Para encontrar \(A\) y \(B\) hay que resolver el siguiente sistema:

\[\left\{\begin{array}{l} 0 = A + B\\ 1 = 2 B \end{array}\right.\]

Por tanto, \(A = - \frac{1}{2}\) y \(B = \frac{1}{2}\). Como consecuencia, aplicando la propiedad asociativa:

\[\begin{split} f (t) =&\mathcal{L}^{- 1} \left( \frac{s}{(s + 1) (s^2 + 1)} \right) =\mathcal{L}^{- 1} \left( \frac{- \frac{1}{2}}{s + 1} + \frac{\frac{1}{2}(s + 1)}{s^2 + 1} \right) =\\ &\mathcal{L}^{- 1} \left( \frac{- \frac{1}{2}}{s + 1} \right) +\mathcal{L}^{- 1} \left( \frac{\frac{1}{2} (s + 1)}{s^2 + 1} \right) \end{split}\]

La transformada inversa de Laplace del primer sumando es directa:

\[\mathcal{L}^{- 1} \left( \frac{- \frac{1}{2}}{s + 1} \right) = - \frac{1}{2} \mathrm{e}^{- t}\]

La transformada inversa de Laplace del segundo sumando requiere descomponer la fracción para encontrar funciones como las que aparecen en las tablas:

\[\mathcal{L}^{- 1} \left( \frac{\frac{1}{2} (s + 1)}{s^2 + 1} \right) = \frac{1}{2} \left[ \mathcal{L}^{- 1} \left( \frac{s}{s^2 + 1} \right) + \mathcal{L}^{- 1} \left( \frac{1}{s^2 + 1} \right) \right] = \frac{1}{2} (\cos t + \sin t)\]

El resultado del problema es:

\[f (t) = - \frac{1}{2} \mathrm{e}^{- t} + \frac{1}{2} (\cos t + \sin t)\]

Utilizando SymPy podemos realizar la descomposición en fracciones simples:

f = s/((s+1)*(s^2+1))
apart(f)

\(\frac{s + 1}{2 \left(s^{2} + 1\right)} - \frac{1}{2 \left(s + 1\right)}\)

Tambien podemos realizar la transformada inversa de Laplace:

invL(f)

\(\left(\frac{\sin{\left(t \right)}}{2} + \frac{\cos{\left(t \right)}}{2}\right) \theta\left(t\right) - \frac{e^{- t} \theta\left(t\right)}{2}\)

3. Para realizar la transformada inversa de Laplace es necesario descomponer f(s) en fracciones simples:

\[f (s) = \frac{3 s + 4}{s^3 (s + 2)} = \frac{A s^2 + B s + C}{s^3} + \frac{D}{s + 2}\]

Tras realizar la suma y simplificar el denominador de la ecuación se obtiene:

\[3 s + 4 = A s^3 + s^2 (2 A + B) + s (2 B + C) + 2 C + D s^3\]

Para encontrar \(A\), \(B\), \(C\) y \(D\) hay que resolver el siguiente sistema de ecuaciones:

\[\left\{\begin{array}{l} 0 = A + D\\ 0 = 2 A + B\\ 3 = 2 B + C\\ 4 = 2 C \end{array}\right.\]

La solución es \(A = - \frac{1}{4}\), \(B = \frac{1}{2}\), \(C = 2\) y \(D = \frac{1}{4}\). Por tanto,

\[f (s) = - \frac{1}{4} + \frac{\frac{1}{2}}{s^2} + \frac{2}{s^3} + \frac{\frac{1}{4}}{s + 2}\]

La realización de la transformada inversa de Laplace es trivial:

\[f (t) = - \frac{1}{4} + \frac{1}{2} t + t^2 + \frac{1}{4} \mathrm{e}^{- 2 t}\]

Podemos comprobar los resultados con SymPy. En primer lugar, realizaremos la descomposición de fracciones simples:

apart((3*s+4)/(s^3*(s+2)))

\(\frac{1}{4 \left(s + 2\right)} - \frac{1}{4 s} + \frac{1}{2 s^{2}} + \frac{2}{s^{3}}\)

La transformada inversa de Laplace:

invL((3*s+4)/(s^3*(s+2)))

\(\frac{t \theta\left(t\right)}{2} + \frac{\left(\left(4 t^{2} - 2 t + 1\right) e^{2 t} - 1\right) e^{- 2 t} \theta\left(t\right)}{4} - \frac{\theta\left(t\right)}{4} + \frac{e^{- 2 t} \theta\left(t\right)}{4}\)

4. De nuevo, antes de realizar la transformada inversa de Laplace hay que realizar la descomposición en fracciones simples. Para realizarla se seguirá en este caso el método que aparece en el Tema 2. Se buscan las raíces del denominador de f(s), estas raíces son dobles y son i y -i. Por tanto la descomposición en raíces simples, teniendo en cuenta que las raíces son dobles, es:

\[f (s) = \frac{s^2}{(s^2 + 1)^2} = \frac{s^2}{(s - \mathrm{i})^2 (s + \mathrm{i})^2} = \frac{A}{(s - \mathrm{i})^2} + \frac{B}{s - \mathrm{i}} + \frac{C}{(s + \mathrm{i})^2} + \frac{D}{s + \mathrm{i}}\]

Según el método, para encontrar \(A\):

\[A = \lim_{s \rightarrow \mathrm{i}} \left[ (s - \mathrm{i})^2 \frac{s^2}{(s - \mathrm{i})^2 (s + \mathrm{i})^2} \right] = \frac{1}{4}\]

Para \(B\), \(C\) y \(D\):

\[\begin{multline} B = \lim_{s \rightarrow \mathrm{i}} \left\{ \frac{\mathrm{d}}{\mathrm{d}s} \left[ (s - \mathrm{i})^2 \frac{s^2}{(s - \mathrm{i})^2 (s + \mathrm{i})^2} \right] \right\} =\lim_{s \rightarrow \mathrm{i}} \left[ \frac{2 s (s + \mathrm{i})^2 - s^2 2 (s + 1)}{(s + \mathrm{i})^4} \right] = \\ \frac{2 \mathrm{i}(- 4) - (- 1) 2 (2 \mathrm{i})}{16} = \frac{- 8 \mathrm{i}+ 4 \mathrm{i}}{16} = - \frac{\mathrm{i}}{4} \end{multline}\]

\[C = \lim_{s \rightarrow \mathrm{i}} \left[ (s + i)^2 \frac{s^2}{(s - \mathrm{i})^2 (s + \mathrm{i})^2} \right] = \frac{1}{4}\]

\[\begin{multline} D = \lim_{s \rightarrow \mathrm{i}} \left\{ \frac{\mathrm{d}}{\mathrm{d}s} \left[ (s + \mathrm{i})^2 \frac{s^2}{(s - \mathrm{i})^2 (s + \mathrm{i})^2} \right] \right\} = \lim_{s \rightarrow \mathrm{i}} \left[ \frac{2 s (s - \mathrm{i})^2 - s^2 2 (s - 1)}{(s - \mathrm{i})^4} \right] =\\ \frac{2 (- \mathrm{i}) (- 4) - (- 1) 2 (-2\mathrm{i})}{16} = \frac{\mathrm{i}}{4} \end{multline}\]

Por tanto, la descomposición en fracciones simples queda de la manera siguiente:

\[\frac{s^2}{(s - \mathrm{i})^2 (s + \mathrm{i})^2} = \frac{\frac{1}{4}}{(s -\mathrm{i})^2} - \frac{\frac{\mathrm{i}}{4}}{s - \mathrm{i}} + \frac{\frac{1}{4}}{(s + \mathrm{i})^2} + \frac{\frac{\mathrm{i}}{4}}{s + \mathrm{i}}\]

La realización de la transformada inversa de Laplace es directa:

\[F (t) = \frac{1}{4} t \mathrm{e}^{\mathrm{i}t} - \frac{\mathrm{i}}{4} \mathrm{e}^{\mathrm{i} t} + \frac{1}{4} t \mathrm{e}^{- \mathrm{i}t} + \frac{\mathrm{i}}{4} \mathrm{e}^{- \mathrm{i} t}\]

Para obtener una respuesta en tiempo real util es necesario eliminar los números complejos de la función. Para ello, hay que considerar la escritura exponencial de un número complejo:

\[\alpha = \beta + \mathrm{i}\gamma \Longrightarrow \mathrm{e}^{\alpha t} = \mathrm{e}^{\beta t} (\cos \gamma t + \mathrm{i}\sin \gamma t)\]

De manera que:

\[\mathrm{e}^{\mathrm{i}t} = \cos t + \mathrm{i}\sin t\] \[\mathrm{e}^{- \mathrm{i}t} = \cos (- t) + \mathrm{i}\sin (- t) = \cos t - \mathrm{i} \sin t\]

Sustituyendo en F(t):

\[F (t) = \frac{1}{4} t (\cos t + \mathrm{i}\sin t) - \frac{\mathrm{i}}{4}(\cos t + \mathrm{i}\sin t) + \frac{1}{4} t (\cos t - \mathrm{i}\sin(t) + \frac{\mathrm{i}}{4} (\cos t - \mathrm{i}\sin t)\]

Operando:

\[F(t) = \frac{1}{2} t \cos t + \frac{1}{2} \sin t\]

También se puede resolver el problema utilizando SymPy.

Podemos realizar la descomposición en factores simples:

f = s^2/(s^2+1)^2

apart(f)

\(\frac{1}{s^{2} + 1} - \frac{1}{\left(s^{2} + 1\right)^{2}}\)

La transformada inversa de Laplace:

invL(f)

\(\frac{\left(t \cos{\left(t \right)} + \sin{\left(t \right)}\right) \theta\left(t\right)}{2} + \sin{\left(t \right)} \theta\left(t\right)\)

Problema 3.3 ★

Usando la técnica de la transformada de Laplace, encontrar las respuestas transitoria y estacionaria del sistema descrito por la ecuación diferencial siguiente:

\[\frac{\mathrm{d}^2 y}{\mathrm{d}t^2} + 3 \frac{\mathrm{d}y}{\mathrm{d}t} + 2 y = 1\]

con las condiciones iniciales \(y (0) = y'(0) = 1\).

Solución

La transformada de una derivada de una función y de la derivada segunda es:

\[\mathcal{L} \left( \frac{\mathrm{d}f(t)}{\mathrm{d}t} \right) = s f(s) - f(0)\]

\[\mathcal{L} \left( \frac{\mathrm{d}^2 f (t)}{\mathrm{d}t^2} \right) = - \left. \frac{\mathrm{d}}{\mathrm{d}t} f (t) \right|_{t = 0} + s^2 f (s) - f (0) s\]

Realizando la transformada de Laplace a la ecuación diferencial a resolver y considerando las condiciones iniciales se obtiene:

\[s^{2_{}} \bar{y} - s - 1 + 3 (s \bar{y} - 1) + 2 \bar{y} = \frac{1}{s}\]

Despejando \(\bar{y}\):

\[\bar{y} = \frac{s^2 + 4 s + 1}{s (s + 1) (s + 2)}\]

Para poder realizar la transformada inversa de Laplace y poder obtener y(t) hay que realizar primero la descomposición en fracciones simples de \(\bar{y} (s)\):

\[\frac{s^2 + 4 s + 1}{s (s + 1) (s + 2)} = \frac{A}{s} + \frac{B}{s + 1} + \frac{C}{s + 2}\]

Sumando las fracciones simples y simplificando el denominador se obtiene la siguiente ecuación:

\[s^2 + 4 s + 1 = A (s + 1) (s + 2) + Bs (s + 2) + Cs (s + 1)\]

Operando se obtiene el sistema de ecuaciones siguiente:

\[\left\{\begin{array}{l} (1) s^2 = (A + B + C) s^2\\ (4) s = (3 A + 3 B + C) s\\ 1 = 2 A \end{array}\right.\]

La solución de este sistema es \(A = \frac{1}{2}\), \(B = 2\) y \(C = -\frac{3}{2}\). Por tanto,

\[\bar{y} = \frac{\frac{1}{2}}{s} + \frac{2}{s + 1} - \frac{\frac{3}{2}}{s + 2}\]

Realizando la transformada inversa de Laplace, consultando las tablas, se obtiene:

\[y (t) = \frac{1}{2} U (t) + 2 \mathrm{e}^{- t} - \frac{3}{2} \mathrm{e}^{- 2 t}\]

donde \(U(t)\) es el escalón unidad.

La parte estacionaria de \(y (t)\) es \(\frac{1}{2} U (t)\) ya que:

\[\lim_{t \rightarrow \infty} y (t) = \frac{1}{2}\]

El resto de la solución es la respuesta transitoria, es decir, la parte de la solución dependiente del tiempo.

Resolución con SymPy

En primer lugar, cargaremos la librería y definiremos las variables \(s\) y \(t\). Como es habitual, especificaremos que el tiempo es una variable real. Además, definiremos \(y\) como una función, ya que \(y(t)\) será nuestra función incógnita:

using SymPy

@syms t::real s y()::real

D = Differential(t);

El siguiente paso es definir las condiciones iniciales, \(y(0) = 1\) y \(\frac{\mathrm{d} y(0)}{\mathrm{d} t} = 1\). Observad la sintaxis:

ic = Dict(y(0)=>1, D(y)(0)=>1);

Para reducir errores es conveniente definir la ecuación diferencial y la llamaremos deq. Para definir una ecuación en SymPy se utiliza la función Eq(lhs, rhs) donde lhs indica la parte izquierda de la igualdad y rhs la parte derecha:

deq = Eq(diff(y(t), t, 2) + 3diff(y(t),t) + 2y(t), 1)

\(2 y{\left(t \right)} + 3 \frac{d}{d t} y{\left(t \right)} + \frac{d^{2}}{d t^{2}} y{\left(t \right)} = 1\)

Ya estamos en condiciones de resolver la ecuación diferencial mediante la instrucción dsolve(), que utiliza internamente la transformada inversa de Laplace. La solución de la ecuación la guardamos en la variable sol:

sol = dsolve(deq, y(t), ics=ic)

\(y{\left(t \right)} = \frac{1}{2} + 2 e^{- t} - \frac{3 e^{- 2 t}}{2}\)

El resultado es una igualdad, en el caso de que nos interese la parte derecha:

sol.rhs

\(\frac{1}{2} + 2 e^{- t} - \frac{3 e^{- 2 t}}{2}\)

Si es la parte izquierda:

sol.lhs

\(y{\left(t \right)}\)

Para encontrar la respuesta estacionaria hay que calcular el límite cuando el tiempo tiende a infinito, ya que entonces habrá desaparecido la influencia de la respuesta transitoria:

sympy.limit(sol.rhs, t, oo)

\(\frac{1}{2}\)

Problema 3.4 ★

La respuesta a un escalón unidad de un sistema viene dada por:

\[y (t) = 1 - \frac{7}{3} \mathrm{e}^{- t} + \frac{3}{2} \mathrm{e}^{- 2 t} - \frac{1}{6} \mathrm{e}^{- 4 t}\]

1. ¿Cuál es la función de transferencia de dicho sistema?

2. ¿Cuál sería la respuesta en tiempo real si la entrada fuera un impulso unidad?

Solución

1. Debemos encontrar aquella función de transferencia G(s) que para una entrada de tipo escalón unidad (\(x (s) = \frac{1}{s}\)) produzca una respuesta como la del enunciado del problema. La transformada de Laplace de y(t) es y(s):

Por tanto la función de transferencia buscada será:

\[G (s) = \frac{y (s)}{x (s)}\]

Para poder encontrar la función de transferencia hay que hacer la transformada de Laplace a la función y(t):

\[y (s) =\mathcal{L} [y (t)] = \frac{1}{s} - \frac{\frac{7}{3}}{s + 1} + \frac{\frac{3}{2}}{s + 2} - \frac{\frac{1}{6}}{s + 4} = \frac{s + 8}{s (s + 1) (s + 2) (s + 4)}\]

Por tanto,

\[G (s) = \frac{s + 8}{(s + 1) (s + 2) (s + 4)}\]

2. La función impulso es la derivada de la función escalón. En consecuencia, la respuesta en tiempo real a una entrada en impulso no será más que derivar y(t) respecto al tiempo:

\[[y (t)]_{\mathrm{impulso}} = \frac{\mathrm{d}}{\mathrm{d}t} [y (t)]_{\mathrm{escalón}} = \frac{7}{3} \mathrm{e}^{- t} - 3 \mathrm{e}^{- 2 t} + \frac{2}{3} \mathrm{e}^{- 4 t}\]

Problema 3.5 ★

Resolver las siguientes ecuaciones diferenciales:

1. \(\frac{\mathrm{d}^2 y}{\mathrm{d}t^2} + 4 y = 3\) con \(y (0) = y' (0) = 1\)

2. \(\frac{\mathrm{d}y}{\mathrm{d}t} + 2 y = 5 \sin (3 t)\) con \(y (0) = 1\)

Solución

1. Para resolver la ecuación diferencial, en primer lugar se realizará la transformada de Laplace:

\[\begin{aligned} & \mathcal{L} \left[ \frac{\mathrm{d}^2 y}{\mathrm{d}t^2} + 4 y \right] =\mathcal{L} [3] & \\ & s^2 \bar{y} - s - 1 + 4 \bar{y} = \frac{3}{s} & \end{aligned}\]

A continuación se despeja \(\bar{y}\):

\[\bar{y} = \frac{s}{s^2 + 4} + \frac{1}{s^2 + 4} + \frac{3}{s (s^2 + 4)}\]

Para obtener la respuesta en tiempo real hay que realizar la transformada inversa. La transformada inversa de los dos primeros sumandos es directa a partir de la tabla de transformadas de Laplace. La transformada inversa del tercer de los sumandos no es directa, por lo que es necesario descomponerlo en fracciones simples:

\[\frac{3}{s (s^2 + 4)} = \frac{As}{s^2 + 4} + \frac{B}{s} = \frac{3}{4} \left( \frac{1}{s} - \frac{s}{s^2 + 4} \right)\]

donde:

\[\begin{aligned} & 3 = As^2 + B (s^2 + 4) & \\ & \left\{\begin{array}{l} 0 = A + B\\ 3 = 4 B \end{array}\right. & \end{aligned}\]

La realización de las transformadas de Laplace es ahora directa. La solución a la ecuación diferencial es:

\[y (t) = \cos (2 t) + \frac{1}{2} \sin (2 t) + \frac{3}{4} (1 - \cos(2 t))\]

También se puede resolver la ecuación utilizando SymPy. A continuación se muestra la serie de instrucciones necesarias para introducir las condiciones iniciales y la ecuación:

using SymPy

@syms s t::real y()

D = Differential(t)

# Condiciones iniciales
ic = Dict(D(y)(0)=>1, y(0)=>1)

# Ecuación diferencial
deq = D(D(y))(t) + 4*y(t) ~ 3

# Resolución de la ecuación diferencial
dsolve(deq, y(t), ics=ic)

\(y{\left(t \right)} = \frac{\sin{\left(2 t \right)}}{2} + \frac{\cos{\left(2 t \right)}}{4} + \frac{3}{4}\)

2. La resolución de la ecuación se puede realizar aplicando la transformada de Laplace o directamente utilizando SymPy:

# Condiciones iniciales
ic2 = Dict(y(0)=>1);

# Ecuación diferencial
deq2 = D(y(t)) + 2y(t) ~ 5sin(3t)

\(2 y{\left(t \right)} + \frac{d}{d t} y{\left(t \right)} = 5 \sin{\left(3 t \right)}\)

dsolve(deq2, y(t), ics=ic2)

\(y{\left(t \right)} = \frac{10 \sin{\left(3 t \right)}}{13} - \frac{15 \cos{\left(3 t \right)}}{13} + \frac{28 e^{- 2 t}}{13}\)

Problema 3.6 ★

Determinar las funciones de transferencia \(\frac{h_2 (s)}{F_i (s)}\) de los sistemas siguientes:

El fluido tiene una densidad constante y las resistencias son lineales (\(F = \frac{h}{R}\)).

Solución

La primera figura muestra un sistema de tanques sin interacción ya que el nivel del segundo tanque no influye en el nivel del primer tanque o en el caudal \(F_1\). En cambio, el segundo sistema se trata de un sistema de tanques con interacción ya que el caudal de salida del primer tanque dependerá tanto del nivel del tanque 1 como del nivel del tanque 2.

Tanques sin interacción

En ambos casos el objetivo del problema es encontrar la siguiente función de transferencia (las primas indican que se tratan de variables de desviación):

\[G (s) = \frac{h_2' (s)}{F_i' (s)}\]

representada en el siguiente diagrama:

Si el sistema no presenta interacciones el bloque anterior es equivalente a un sistema en el que hay dos procesos en serie como el que se muestra a continuación:

en el que cada uno de los bloques representa a un depósito. Las funciones de transferencia de estos bloques serán:

\[\begin{aligned} G_1 (s) = \frac{h_1' (s)}{F_i' (s)} & & \\ G_2 (s) = \frac{h_2' (s)}{h_1' (s)} & & \end{aligned}\]

Al tratarse de dos procesos en serie:

\[G (s) = G_1 (s) G_2 (s)\]

Para calcular \(G_1\) hay que obtener el modelo matemático que represente la variación con el tiempo del nivel de un depósito en función del caudal de entrada. Se trata del mismo modelo que el del ejemplo del capítulo 2, por tanto:

\[G_1 = \frac{K_1}{\tau_1 s + 1} = \frac{h'_1}{F'_i}\]

donde \(K_1 = R_1\) y \(\tau_1 = A_1 R_1\).

De igual manera para el segundo de los depósitos:

\[G_2 = \frac{K_2}{\tau_2 s + 1} = \frac{h'_2}{F'_1}\]

donde \(K_2 = R_2\), \(\tau_2 = A_2 R_2\) y \(F_1' = \frac{h'_1}{R_1}\).

Sustituyendo \(F'_1\) y operando se encuentra la función de transferencia buscada:

\[G = \frac{K_2}{(\tau_1 s + 1) (\tau_2 s + 1)}\]

Se puede apreciar que se trata de un sistema de segundo orden. La respuesta de este sistema será sobreamortiguada (\(\zeta > 1\)), ya que las raíces del denominador son números reales.

En la figura siguiente se muestra la variación con el tiempo del nivel del depósito para un cambio en el caudal de entrada en forma de escalón unidad. Se observa que la dinámica es más lenta a medida que aumenta el número de tanques:

Tanques con interacción

Para obtener el modelo matemático se realizarán los Balances Macroscópicos de Materia a los tanques (se ha supuesto que la densidad es constante y que no depende del tiempo):

\[\begin{aligned} \text{Tanque 1: }A_1 \frac{\mathrm{d}h_1}{\mathrm{d}t} &= F_i - F_1\\ \text{Tanque 2: } A_2 \frac{\mathrm{d}h_2}{\mathrm{d}t} &= F_1 - F_2 \end{aligned}\]

Los caudales de salida de los tanques dependen de las resistencias lineales:

\[\begin{aligned} F_1 &= \frac{h_1 - h_2}{R_1}\\ F_2 &= \frac{h_2}{R_2} \end{aligned}\]

El caudal de salida del primer tanque depende, lógicamente, de la diferencia de alturas entre los dos tanques ya que la presión ejercida por la columna de fluido depende de esa diferencia de nivel. En el caso de que \(h_2\) fuera mayor que \(h_1\) el fluido circularía en sentido contrario. Esta ecuación marca la interacción entre los tanques e impide que las ecuaciones de los balances de materia se puedan resolver de manera independiente. Las dos ecuaciones diferenciales se deben resolver de manera simultanea.

Sustituyendo los caudales en los balances de materia se encuentra:

\[\begin{aligned} A_1 R_1 \frac{\mathrm{d}h_1}{\mathrm{d}t} + h_1 - h_2 = R_1 F_i & \\ A_2 R_2 \frac{\mathrm{d}h_2}{\mathrm{d}t} + \left( 1 + \frac{R_2}{R_1} \right) h_2 - \frac{R_2}{R_1} h_1 = 0 \end{aligned}\]

Tal como se ha comentado las dos ecuaciones diferenciales están ligadas (las incógnitas \(h_1 (t)\) y \(h_2 (t)\) aparecen en ambas ecuaciones), se deben resolver simultáneamente. Esta es la diferencia fundamental con los sistemas sin interacción.

A continuación se realizan los balances de materia en estado estacionario:

\[\begin{aligned} h_{1, e} - h_{2, e} = R_1 F_i\\ \left( 1 + \frac{R_2}{R_1} \right) h_{2, e} - \frac{R_2}{R_1} h_{1, e} = 0 \end{aligned}\]

donde los subíndices \(e\) indican los valores en estado estacionario.

Restando los balances de materia y los balances en estado estacionario:

\[\begin{aligned} A_1 R_1 \frac{\mathrm{d}h'_1}{\mathrm{d}t} + h'_1 - h'_2 = R_1 F'_i\\ A_2 R_2 \frac{\mathrm{d}h'_2}{\mathrm{d}t} + \left( 1 + \frac{R_2}{R_1} \right) h'_2 - \frac{R_2}{R_1} h'_1 = 0 \end{aligned}\]

tomando como variables de desviación:

\[\begin{aligned} h'_1 &= h_1 - h_{1, e}\\ h'_2 &= h_2 - h_{2, e}\\ F'_i &= F_i - F_{i, e} \end{aligned}\]

Haciendo la transformada deLaplace:

\[\begin{aligned} (A_1 R_1 s + 1) h'_1 (s) - h'_2 (s) = R_1 F'_i (s)\\ \left[ A_2 R_2 s + \left( 1 + \frac{R_2}{R_1} \right) \right] h'_2 (s) - \frac{R_2}{R_1} h'_1 (s) = 0 \end{aligned}\]

Resolviendo el sistema de ecuaciones:

\[\begin{aligned} h'_1 &= \frac{(\tau_2 R_1) s + R_1 + R_2}{\tau_1 \tau_2 s^2 + (\tau_1 + \tau_2 + A_1 R_2) s + 1} F_i'\\ h'_2 &= \frac{R_2}{\tau_1 \tau_2 s^2 + (\tau_1 + \tau_2 + A_1 R_2) s + 1} F_i' \end{aligned}\]

donde \(\tau_1 = A_1 R_1\) y \(\tau_2 = A_2 R_2\).

La respuesta es igual a la de los tanques sin interacción excepto por el término \(A_1 R_2\) del denominador. Es el factor de interacción, ya que cuanto mayor sea mayor es la interacción entre los tanques.

La respuesta, de nuevo, es sobreamortiguada pero más lenta que la del sistema sin interacción. Para el sistema de tanques sin interacción, la constante de tiempo es:

\[\tau_{\mathrm{sin}} = \sqrt{\tau_1 + \tau_2}\]

Para el sistema de tanques con interacción:

\[\tau_{\mathrm{con}} = \sqrt{\tau_1 + \tau_2 + A_1 R_2}\]

Al ser el factor de interacción positivo es evidente que \(\tau_{\mathrm{con}} > \tau_{\mathrm{sin}}\).

Problema 3.7 ★

En un reactor tanque agitado se lleva a acabo una reacción irreversible de primer orden, \(A \rightarrow B\), con un coeficiente cinético k = 1 h^-1. El caudal de alimentación es de 10 l/min, con \(c_{A_0}\) = 0.5 mol/l. El volumen del tanque es de 1 m³. Determinar el efecto causado por una subida instantánea de la concentración a 0.7 mol/l. Indicar cual hubiese sido el transitorio, si la perturbación se hubiera producido en forma de impulso unidad. Cuáles serán la ganancia y la constante de tiempo del sistema?

Solución

A partir del enunciado del problema se puede realizar el siguiente diagrama de flujo:

donde:

\[\left\{\begin{array}{l} q = 10 \text{l/min} = 0.6\ \mathrm{ m^3/h}\\ V = 1 \ \mathrm{m^3} \end{array}\right.\]

El subíndice A indica que es la concentración del componente A y el subíndice i indica que se trata de la concentración de entrada.

La reacción química que se produce en el reactor es la siguiente:

\[A \rightarrow B\]

cuya constante cinética es \(k = 1 \text{ h}^{-1}\).

La concentración de entrada en el reactor antes de que se produzca el cambio de la concentración en escalón es \(c_{A i_0} = c_{A i} (t = 0) = 0.5\ \mathrm{mol/l}\).

A continuación hay que plantear el balance macroscópico de materia en estado no estacionario (BMM):

\[V \frac{\mathrm{d}c_A}{\mathrm{d}t} = q c_{A i} - q c_A - k c_A V\]

donde se ha supuesto que los caudales volumétricos (q), las densidades (\(\rho\)) y el volumen del reactor (V) son constantes e independientes del tiempo.

El balance macroscópico de materia  en estado estacionario (BMMe) del sistema es:

\[0 = q c_{A i_0} - q c_{A e} - k c_{A_{} e} V\]

donde el subíndice e indica que se trata de las concentraciones en estado estacionario. Este balance representa la situación antes de que se produzca ningún cambio.

Operando la ecuación anterior:

\[q c_{A i_0} = c_{A e} (q + k V)\]

El tiempo de residencia del reactor es:

\[\theta = \frac{V}{q}\]

Sustituyendo y operando:

\[c_{A e} = \frac{c_{A i_0}}{1 + k \theta}\]

Sustituyendo los calores de las variables se obtiene que la concentración de salida del reactor del componente A en estadio estacionario, es decir, para tiempos menores o iguales a cero es:

\[c_{A e} = 0.1875 \text{ mol/l}\]

El siguiente paso es definir las variables de desviación para la entrada:

\[C_{A i} = c_{A i} - c_{A i_0}\]

y para la salida:

\[C_A = c_A - c_{A e}\]

Restando el BMMe al BMM se obtiene:

\[V \frac{\mathrm{d}c_A}{\mathrm{d}t} = q (c_{A i} - c_{A i_0}) - q (c_A - c_{A e}) - k (c_A - c_{A e}) V\]

Sustituyendo las variables de desviación y operando:

\[V \frac{\mathrm{d}C_A}{\mathrm{d}t} = q C_{A i} - C_A (q + k V)\]

Esta última ecuación es el modelo matemático del proceso planteado. Con objeto de resolver la ecuación diferencial se aplica la transformada de Laplace a la ecuación anterior:

\[V s \overline{C_A} = q \overline{C_{A i}} - \overline{C_A} (q + k V)\]

Para poder responder a las preguntas propuestas por el problema hay que plantear la siguiente función de transferencia:

\[G = \frac{\overline{C_A}}{\overline{C_{A i}}}\]

Combinando las dos ecuaciones anteriores se obtiene:

\[G = \frac{q}{V s + q + k V} = \frac{\frac{q}{q + k V}}{\frac{V}{q + k V} s + 1} = \frac{\frac{1}{1 + k \theta}}{\frac{1}{\theta^{- 1} + k} s + 1}\]

Lógicamente la función de transferencia es la propia de un sistema de primer orden. La ganancia del proceso será:

\[K = \frac{1}{1 + k \theta} = 0.375\]

Y la constante de tiempo es:

\[\tau = \frac{1}{\theta^{- 1} + k} = 0.625 h = 37.5 \min\]

El efecto de un cambio súbito de concentración de 0.5 a 0.7 mol/l, es decir, para una entrada en escalón de altura 0.2 mol/l es:

\[C_A = 0.375 \cdot 0.2 \text{ mol/l} \left( 1 - \mathrm{e}^{- \frac{t}{37.5 \min}} \right) = 0.075 \text{ mol/l} \left( 1 - \mathrm{e}^{- \frac{t}{37.5 \min}} \right)\]

Deshaciendo las variables de desviación se obtiene la ecuación que describirá la variación de la concentración del componente A en función del tiempo:

\[c_A = 0.075 \text{ mol/l} \left( 1 - \mathrm{e}^{- \frac{t}{37.5 \min}} \right) + 0.1875 \text{ mol/l}\]

Para una entrada que sea un impulso unidad la respuesta es:

\[C_A = \frac{K A}{\tau} e^{- \frac{t}{\tau}} = 0.6 \text{ mol/l } \mathrm{e}^{- \frac{t}{37.5 \min}}\]

Problema 3.8 ★

Considerar el tanque encamisado de la figura utilizado como precalentador. Suponiendo que la capacitancia de la pared del tanque es despreciable y que la temperatura en el interior del mismo es uniforme.

Deducir:

1. La función de transferencia para variaciones de \(T_1\) y/o \(T_s\).

2. El diagrama de bloques del sistema.

3. La respuesta del sistema para una entrada en escalón unidad en \(T_1\) (\(T_s\) cte).

Solución

1. En este problema se plantea un proceso con dos entradas (\(T_1\) y \(T_s\)) y una respuesta (\(T_2\)). En primer lugar se plantea el Balance Macroscópico de Energía:

\[m C_p \frac{\mathrm{d}T_2}{\mathrm{d}t} = w C_p (T_1 - T_2) + U A (T_s - T_2) \tag{3.11}\]

donde \(m\) es la cantidad de líquido a calentar en el interior del tanque, \(C_p\) es el calor específico, \(U\) es el coeficiente global de transferencia de calor y \(A\) es el área de intercambio. Como el enunciado del problema no indica que exista ningún sistema de control de nivel se asume que los caudales de las corrientes de entrada y de salida son iguales (\(w_1 = w_2 = w\)).

Operando y agrupando constantes en la Ecuación 3.11 se obtiene:

\[\tau \frac{\mathrm{d}T_2}{\mathrm{d}t} + T_2 = k_1 T_1 + k_2 T_s \tag{3.12}\]

donde:

\[\begin{aligned} & \tau = \frac{m C_p}{w C_p + U A} & \\ & k_1 = \frac{w C_p}{w C_p + U A} & \\ & k_2 = \frac{U A}{w C_p + U A} & \end{aligned}\]

A continuación se realiza el Balance Macroscópico de Energía en estado estacionario:

\[0 = k_1 T_{1, e} + k_2 T_{s, e} - T_{2, e} \tag{3.13}\]

Restando la Ecuación 3.12 y Ecuación 3.13, definiendo variables de desviación y realizando la transformada de Laplace a la ecuación resultante se obtiene:

\[\begin{aligned} & \tau s \bar{T}'_2 + \bar{T}_2' = k_1 \bar{T}_1' + k_2 \bar{T}_s' & \\ & \bar{T}_2' = \frac{k_1}{\tau s + 1} \bar{T}_1' + \frac{k_2}{\tau s + 1} \bar{T}_s' & \end{aligned}\]

donde las prima indica que se trata de variables de desviación.

2. El diagrama de bloques de este proceso es:

donde:

\[\begin{aligned} & G_1 = \frac{k_1}{\tau s + 1} & \\ & G_2 = \frac{k_2}{\tau s + 1} & \end{aligned}\]

3. La variación de la temperatura de la corriente de salida del tanque para un cambio en escalón unidad de \(T_1\) será:

\[\bar{T}_2' = \frac{k_1}{\tau s + 1} \frac{1}{s} + 0\]

Realizando la transformada inversa de Laplace y deshaciendo las variables de desviación se obtiene:

\[T_2 = T_{2, e} + k_1 \left( 1 - \mathrm{e}^{- \frac{t}{\tau}} \right)\]

Man, Yiu-Kwong. 2009. «An Improved Heaviside Approach to Partial Fraction Expansion and Its Applications». International Journal of Mathematical Education in Science and Technology 40 (6): 808-14. https://doi.org/10.1080/00207390902825310.