Transformada inversa de Laplace numérica

control de procesos
julia
Autor/a

runjaj

Fecha de publicación

21 de mayo de 2023

Es habitual cuando se trabaja el control de procesos utilizar la transformada de Laplace para convertir las ecuaciones diferenciales que describen el comportamiento dinámico de los diferentes componentes del lazo de control en ecuaciones algebráicas.

Aunque no es la mejor opción se puede realizar la transformada inversa de Laplace numéricamente. Para realizar las simulaciones utilizaremos el paquete de Julia InverseLaplace.jl que nos permite utilizar el método de Talbot de manera extremadamente simple.

Descripción del lazo de control

Vamos a plantear una lazo de control sencillo (Figura 1) compuesto por un controlador de tipo proporcional, un elemento final de control, un proceso a controlar y un medidor. Las variables de este lazo de control por retroalimentación son:

  • \(y\), variable controlada
  • \(y_{sp}\), consigna, es decir, el valor deseado de la variable controlada
  • \(y_m\), valor medido de la variable controlada
  • \(\varepsilon\), error
  • \(c\), acción de control
  • \(f\), variable manipulable
Código 
using LaTeXStrings, TikzPictures

TikzPicture(L"""
\draw [draw=black, fill=yellow!10, very thick, dashed] (1.5,-1) rectangle ++(5,2.5) node [below, xshift=-2.5cm]{\textbf{Sistema de control}};;
\node[rectangle, align=right] (entrada) {$y_{sp}(t)$};
\node[circle, , very thick, draw, right=1.5cm of entrada, minimum size=0.5cm] (comp) {};
\node[rectangle, very thick, draw, right=1.5cm of comp, minimum size=1.5cm] (control) {Controlador};
\node[rectangle, very thick, draw, right=1.5cm of control, minimum size=1.5cm, align=center] (efc) {Elem. final\\de control};
\node[rectangle, very thick, draw, right=1.5cm of efc, minimum size=1.5cm] (proc) {Proceso};
\node[rectangle, align=left, right=1.5cm of proc] (salida) {$y(t)$};
\node[rectangle, very thick, draw, below=0.75cm of efc, minimum size=1.5cm] (medidor) {Medidor};
\node [circle, right=.75cm of proc] (output) {};

\draw [->, very thick] (entrada) -- node[above, pos=0.95] {{$+$}} (comp);
\draw [->, very thick] (comp) to node[above] {$\varepsilon(t)$} (control);
\draw [->, very thick] (control) to node[above]{$c(t)$} (efc);
\draw [->, very thick] (efc) to node[above]{$f(t)$} (proc);
\draw [->, very thick] (proc) -- (salida);
\draw [->, very thick] (output.west) |- (medidor);
\draw [->, very thick] (medidor)  -| node[pos=0.3, above right] {$y_m(t)$} node[pos=0.99, left] {{$-$}}   (comp) ;

""", options=raw"font=\sffamily",  preamble=raw"\usetikzlibrary{positioning}", width="400px")

Figura 1: Diagrama de bloques del lazo de control por retroalimentación.

Los diferentes componentes del lazo de control vienen descritos por sus funciones de transferencia (Figura 2):

  • Controlador: \(G_c = K_c\), donde \(K_c\) es la ganancia del controlador
  • Elemento final de control: \(G_v = \frac{2}{s+1}\)
  • Proceso: \(G_p = \frac{4}{5s+1}\)
  • Medidor: \(G_m = \frac{1}{\frac{s}{2}+1}\)
Código 
TikzPicture(L"""
\draw [draw=black, fill=yellow!10, very thick, dashed] (1.5,-1) rectangle ++(4.5,2.5) node [below, xshift=-2.25cm]{\textbf{Sistema de control}};;
\node[rectangle, align=right] (entrada) {$y_{sp}(s)$};
\node[circle, , very thick, draw, right=1.5cm of entrada, minimum size=0.5cm] (comp) {};
\node[rectangle, very thick, draw, right=1.5cm of comp, minimum size=1.5cm] (control) {$G_c(s)$};
\node[rectangle, very thick, draw, right=1.5cm of control, minimum size=1.5cm, align=center] (efc) {$G_f(s)$};
\node[rectangle, very thick, draw, right=1.5cm of efc, minimum size=1.5cm] (proc) {$G_p(s)$};
\node[rectangle, align=left, right=1.5cm of proc] (salida) {$y(s)$};
\node[rectangle, very thick, draw, below=0.75cm of efc, minimum size=1.5cm] (medidor) {$G_m(s)$};
\node [circle, right=.75cm of proc] (output) {};

\draw [->, very thick] (entrada) -- node[above, pos=0.95] {{$+$}} (comp);
\draw [->, very thick] (comp) to node[above] {$\varepsilon(s)$} (control);
\draw [->, very thick] (control) to node[above]{$c(s)$} (efc);
\draw [->, very thick] (efc) to node[above]{$f(s)$} (proc);
\draw [->, very thick] (proc) -- (salida);
\draw [->, very thick] (output.west) |- (medidor);
\draw [->, very thick] (medidor)  -| node[pos=0.3, above right] {$y_m(t)$} node[pos=0.99, left] {{$-$}}   (comp) ;


""", options=raw"font=\sffamily",  preamble=raw"\usetikzlibrary{positioning}", width="400px")

Figura 2: Funciones de transferencia del lazo de control.

Aunque la dinámica de cada uno de los elementos es simple, el conjunto tiene una dinámica compleja. Vamos a estudiar el efecto de la ganancia del control sobre la respuesta del lazo de control.

Para poder simular la dinámica del lazo de control, definiremos la función de transferencia del bloque equivalente (Figura 3):

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

Código 
TikzPicture(L"""
\node[rectangle, align=right] (entrada) {$y_{sp}(s)$};
\node[rectangle, very thick, draw, right=1.5cm of entrada, minimum size=1.5cm] (bloque) {$G(s)$};
\node[rectangle, align=left, right=1.5cm of bloque] (salida) {$y(s)$};

\draw [->, very thick] (entrada) -- (bloque);
\draw [->, very thick] (bloque) to (salida);
""", options=raw"font=\sffamily",  preamble=raw"\usetikzlibrary{positioning}", width="200px")

Figura 3: Bloque equivalente al lazo de contro.

Cálculos con SymPy

Aunque el cálculo no es excesivamente complejo, resulta mucho más rápido, sencillo y libre de errores realizándolo con ayuda de SymPy.jl y Julia.

Empezaremos cargando los paquetes necesarios (SymPy para cálculo algebraico, Plots para representación de gráficos e InverseLaplace para el cálculo numérico de la transformada inversa de Laplace):

using SymPy, Plots, InverseLaplace

Definimos las variables de SymPy que utilzaremos para los cálculos algebráicos:

@syms s t::real Kc::positive=>"K_c"
(s, t, K_c)

Hemos indicado que \(t\) es una variable real, a diferencia de la variable compleja \(s\). La ganancia del controlador es un número real positivo. La parte =>K_ces para que en los resultados de los cálculos se obtenga \(K_c\) y no \(Kc\).

A continuación definimos las funciones de transferencia de los bloques del lazo de control:

Gc = Kc

\(K_{c}\)

Gv = 2/(s+1)

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

Gp = 4/(5s+1)

\(\frac{4}{5 s + 1}\)

Gm = 1/(s/2+1)

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

El cálculo de la función de tranferencia \(G = \frac{y(s)}{y_{sp}(s)}\) es:

G = Gc*Gv*Gp/(1+Gc*Gv*Gp*Gm)

\(\frac{8 K_{c}}{\left(s + 1\right) \left(5 s + 1\right) \left(\frac{8 K_{c}}{\left(\frac{s}{2} + 1\right) \left(s + 1\right) \left(5 s + 1\right)} + 1\right)}\)

Podemos simplificar el resultado:

simplify(G)

\(\frac{8 K_{c} \left(s + 2\right)}{16 K_{c} + \left(s + 1\right) \left(s + 2\right) \left(5 s + 1\right)}\)

La ventaja es que podemos realizar cálculos de manera simbólica. Por ejemplo, podemos calcular el offset, el error residual, para un cambio en la consigna en forma de escalón unidad:

ysp = 1/s

\(\frac{1}{s}\)

El offset es:

Offset = limit(s*ysp - s*G*1/s, s, 0)

\(- \frac{8 K_{c}}{8 K_{c} + 1} + 1\)

Simulación con InverseLaplace

Vamos a simular la respuesta del lazo de control para un cambio en la consigna en forma de escalón unidad para los siguientes valores de consigna: 1/100, 1/10 y 1.

En primer lugar, definiremos la función que simulará la respuesta para \(K_c = 1\):

y_1 = ILT(lambdify(G(Kc=>1)*ysp))
Talbot{SymPy.var"#118#119"{SymPy.var"###292"}}(Nterms=32)

En esta instrucción hay varias partes a destacar:

  • La respuesta del lazo de control se calcula como \(y(s) = G(s) y_{sp}(s)\)
  • G(Kc=>1)sustituye el valor de Kc en G
  • lambdify transforma la función de SymPy en una función de Julia que puede manejar InverseLaplace
  • ILT es la función de InverseLaplace para realizar la transformada inversa de Laplace por el método de Talbot

Definimos las funciones para los otros valores de la ganancia del controlador:

y_10 = ILT(lambdify(G(Kc=>1/10)*1/s))
y_100 = ILT(lambdify(G(Kc=>1/100)*1/s))
Talbot{SymPy.var"#118#119"{SymPy.var"###294"}}(Nterms=32)

Ya solo queda representar las funciones anteriores para ver como cambia la variable controlada para un cambio en la consigna en forma de escalón unidad:

plot(0:.1:20, [t->1, t->y_1.(t), t->y_10.(t),t->y_100.(t)],
    label=[L"$y_{sp}$" L"$K_c=1$" L"$K_c=1/10$" L"$K_c=1/100$"],
    size=(600, 400))

Vemos como obtenemos respuestas sobreamortiguadas para los dos valores menores de la ganancia. Al aumentar la ganancia proporcional logramos una respuesta más rápida y subamortiguada. Tal como era de esperar, el offset disminuye al aumentar \(K_c\).

Con un valor de \(K_c\) de 2.5 llegamos al límite de estabilidad, tenemos una respuesta con oscilaciones sostenidas:

y_25 = ILT(lambdify(G(Kc=>2.5)*1/s))

plot(0:.1:10, [t->1, t->y_25.(t)], 
    label=[L"$y_{sp}$" L"$K_c=2.5$"], size=(600, 400))

Valores más elevados, obtienen respuestas no acotadas. El lazo de control es inestable:

y_35 = ILT(lambdify(G(Kc=>3.5)*1/s))

plot(0:.1:10, [t->1, t->y_35.(t)],
    label=[L"$y_{sp}$" L"$K_c=3.5$"], size=(600, 400))

El problema del cálculo numérico de la tranformada inversa de Laplace es que es bastante inestable y falla para tiempos un poco elevados:

plot(0:.1:15, [t->1, t->y_35.(t)],
    label=[L"$y_{sp}$" L"$K_c=3.5$"], size=(600, 400))

Conclusión

El cálculo numérico de la transferencia de la Laplace es una buena opción para realizar simulaciones sencillas pero con la suficiente complejidad para no poder encontrar una solución analítica. La principal ventaja es su sencillez con Julia.

En el caso de sistemas más complejos, es preferible resolver las ecuaciones diferenciales numéricamente ya que los métodos numéricos son mucho más potentes.