Simulación de un lazo de control con ModelingToolkit

control de procesos
julia
Autor/a

runjaj

Fecha de publicación

22 de mayo de 2023

En la entrada anterior dedicada al cálculo numérico de la transformada de Laplace decía que es preferible simular las ecuaciones diferenciales directamente a utilizar el método de Talbot para estimar numéricamente la transformada de Laplace. La última opción solo la usaría para problemas sencillos en los que estamos trabajando con funciones de transferencia y nos puede dar pereza reescribir el problema.

Me quedé con ganas de plantear la simulación de esa entrada utilizando las ecuaciones diferenciales y métodos numéricos más potentes.

Planteamiento del modelo matemático

El sistema a controlar es el lazo de control por retroalimentación siguiente:

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.

El comparador viene descrito por la siguiente ecuación:

\[\varepsilon(t) = y_{sp}(t)-y_m(t)\]

El controlador es un controlador proporcional:

\[c(t) = K_c \varepsilon(t)\]

Nota

Estamos utilizando variables de desviación sobre el estado inicial para que todas las condiciones iniciales sean 0.

La válvula, proceso y medidor vienen definidos por procesos de primer orden. Los sistemas de primer orden vienen descritos por una ecuación diferencial de primer orden del tipo:

\[\frac{\mathrm{d}y(t)}{\mathrm{d}t} = \frac{K f(t) - y(t)}{\tau}\]

donde \(f(t)\) es la entrada y \(y(t)\) es la respuesta de ese proceso de primer orden.

Las ecuaciones quedan así:

\[\begin{align} \frac{\mathrm{d}f(t)}{\mathrm{d}t} &= \frac{K_v c(t) - f(t)}{\tau_v}\\ \frac{\mathrm{d}y(t)}{\mathrm{d}t} &= \frac{K_p f(t) - y(t)}{\tau_p}\\ \frac{\mathrm{d}y_m(t)}{\mathrm{d}t} &= \frac{K_m y(t) - y_m(t)}{\tau_m} \end{align}\]

Los parámetros de los componentes el lazo de contro son:

Elemento Ganancia (\(K\)) Cte. de tiempo (\(\tau\))
Válvula 2 1
Proceso 4 5
Medidor 1 0.5

Escritura de las ecuaciones en ModelingTookit

Para resolver el sistema de ecuaciones anterior con Julia utilizaremos ModelingToolkit.jl. Es una auténtica maravilla, potente y sencillo de utilizar.

Empezaremos cargando las librerías a utilizar:

using ModelingToolkit, OrdinaryDiffEq, Plots

OrdinaryDiffEq es la librería de resolución numérica de ecuaciones diferenciales que utilizará ModelingToolkit para resolver nuestro sistema de ecuaciones.

A continuación, vamos a definir con el macro @variables la variable independiente, el tiempo, y las funciones que componen nuestro sistema:

@variables t err(t) c(t) f(t) y(t) ym(t) ysp(t)

\[ \begin{equation} \left[ \begin{array}{c} t \\ \mathrm{err}\left( t \right) \\ c\left( t \right) \\ f\left( t \right) \\ y\left( t \right) \\ \mathrm{ym}\left( t \right) \\ \mathrm{ysp}\left( t \right) \\ \end{array} \right] \end{equation} \]

También definiremos el único parámetro que va a tener nuestro modelo matemático, la ganancia del controlador \(K_c\):

@parameters Kc

\[ \begin{equation} \left[ \begin{array}{c} Kc \\ \end{array} \right] \end{equation} \]

También vamos a definir el vector con las condiciones iniciales:

u0 = [
    ysp => 0,
    err => 0,
    c => 0,
    f => 0,
    y => 0,
    ym => 0
]
6-element Vector{Pair{Num, Int64}}:
 ysp(t) => 0
 err(t) => 0
   c(t) => 0
   f(t) => 0
   y(t) => 0
  ym(t) => 0

Para poder escribir las ecuaciones diferenciales, tenemos que definir el operador diferencial:

D = Differential(t)
(::Differential) (generic function with 2 methods)

Ya estamos en condiciones de escribir las ecuaciones que describe los diferentes elementos del lazo de control:

eqs = [
    ysp ~ 1,
    err ~ ysp - ym,
    c ~ Kc*err,
    D(f) ~ (2*c-f)/1,
    D(y) ~ (4*f-y)/5,
    D(ym) ~ (y-ym)/.5
]

\[ \begin{align} \mathrm{ysp}\left( t \right) =& 1 \\ \mathrm{err}\left( t \right) =& - \mathrm{ym}\left( t \right) + \mathrm{ysp}\left( t \right) \\ c\left( t \right) =& Kc \mathrm{err}\left( t \right) \\ \frac{\mathrm{d} f\left( t \right)}{\mathrm{d}t} =& 2 c\left( t \right) - f\left( t \right) \\ \frac{\mathrm{d} y\left( t \right)}{\mathrm{d}t} =& \frac{1}{5} \left( - y\left( t \right) + 4 f\left( t \right) \right) \\ \frac{\mathrm{d} \mathrm{ym}\left( t \right)}{\mathrm{d}t} =& 2 \left( - \mathrm{ym}\left( t \right) + y\left( t \right) \right) \end{align} \]

El siguiente paso es definir el modelo matemático a resolver:

@named model = ODESystem(eqs, t)

\[ \begin{align} \frac{\mathrm{d} f\left( t \right)}{\mathrm{d}t} =& 2 c\left( t \right) - f\left( t \right) \\ \frac{\mathrm{d} y\left( t \right)}{\mathrm{d}t} =& \frac{1}{5} \left( - y\left( t \right) + 4 f\left( t \right) \right) \\ \frac{\mathrm{d} \mathrm{ym}\left( t \right)}{\mathrm{d}t} =& 2 \left( - \mathrm{ym}\left( t \right) + y\left( t \right) \right) \\ \mathrm{ysp}\left( t \right) =& 1 \\ \mathrm{err}\left( t \right) =& - \mathrm{ym}\left( t \right) + \mathrm{ysp}\left( t \right) \\ c\left( t \right) =& Kc \mathrm{err}\left( t \right) \end{align} \]

Queda un sistema de 6 ecuaciones, pero que ModelingToolkit tiene la capacidad de simplificar algebraicamente nuestras ecuaciones:

sys = structural_simplify(model)

\[ \begin{align} \frac{\mathrm{d} f\left( t \right)}{\mathrm{d}t} =& 2 c\left( t \right) - f\left( t \right) \\ \frac{\mathrm{d} y\left( t \right)}{\mathrm{d}t} =& \frac{1}{5} \left( - y\left( t \right) + 4 f\left( t \right) \right) \\ \frac{\mathrm{d} \mathrm{ym}\left( t \right)}{\mathrm{d}t} =& 2 \left( - \mathrm{ym}\left( t \right) + y\left( t \right) \right) \end{align} \]

Reduce nuestro modelo matemático a la mitad de ecuaciones. Ya solo nos queda dar valor a la ganancia del controlador para poder realizar la simulación:

p = [Kc => 0.5]
1-element Vector{Pair{Num, Float64}}:
 Kc => 0.5

Resolución numérica del sistema de ecuaciones

Definimos los tiempos para los que deseamos realizar la simulación:

tspan = (0.0, 20.0)
(0.0, 20.0)

Con todos estos elementos ya podemos definir el problema a resolver numéricamente:

prob = ODEProblem(sys, u0, tspan, p, jac=true)
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 20.0)
u0: 3-element Vector{Float64}:
 0.0
 0.0
 0.0

Con jac=true estamos indicando a ModelingToolikt que calcule el jabobiano del sistema algebraicamente.

Finalmente resolvemos numéricamente las ecuaciones:

sol = solve(prob, Tsit5())
retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 31-element Vector{Float64}:
  0.0
  9.999999999999999e-5
  0.0010999999999999998
  0.011099999999999997
  0.04624583951849052
  0.0977638858285024
  0.1599124918811575
  0.24131464158918597
  0.343878035614621
  0.4751152468831308
  0.6402600175761708
  0.8477799735646785
  1.106623290785601
  ⋮
  5.374376168512642
  6.198252551930584
  7.0770366453045765
  8.068608853826673
  9.241895015867271
 10.450692712559249
 11.702739922882195
 13.131150188939367
 14.782946016540626
 16.338933052029752
 18.3644329136901
 20.0
u: 31-element Vector{Vector{Float64}}:
 [0.0, 0.0, 0.0]
 [9.999500016665579e-5, 3.999840004133117e-9, 2.66645334351957e-13]
 [0.0010993952216748256, 4.837871004855705e-7, 3.546211559846735e-10]
 [0.011038621304711685, 4.906580332721622e-5, 3.614801384679568e-7]
 [0.04519250266842801, 0.0008398314074215098, 2.542003397358568e-5]
 [0.09313138616054098, 0.003677220121103394, 0.00023056226937609458]
 [0.14774348325545827, 0.009599705594843893, 0.0009609855197739362]
 [0.2142205737543654, 0.021171009359756812, 0.0030999158764550123]
 [0.290283022039332, 0.04129471045726331, 0.008291135620539708]
 [0.3758878635889408, 0.07486996672686692, 0.01979937085817074]
 [0.46620314969512494, 0.12739063064101813, 0.04284422935087423]
 [0.554310443742814, 0.2056056163251664, 0.08544397963801113]
 [0.628234716061676, 0.3152556694820271, 0.1577397149405979]
 ⋮
 [0.08023115355704097, 0.8717630071442579, 0.920086781037451]
 [0.10605597062873548, 0.7938617681132257, 0.8393409257643297]
 [0.1607992008452688, 0.7520980503849545, 0.776576969666896]
 [0.21117672392097075, 0.7528877521710114, 0.7532819902459517]
 [0.23002168592875502, 0.7833169554716718, 0.7710454140832245]
 [0.21700826568014261, 0.8081467794600239, 0.7990592527650897]
 [0.19900621299314872, 0.8125445572336623, 0.8114956783572311]
 [0.19252188878726254, 0.8037093935564069, 0.8068202542155899]
 [0.19780248833026234, 0.7969326628055511, 0.7983512554754443]
 [0.20165643143442172, 0.7980095875303407, 0.7974623317898745]
 [0.20076312091406062, 0.8006572068110831, 0.8002140714383854]
 [0.19963739469437597, 0.8005487082082924, 0.8006547756483468]

y representar los resultados:

plot(sol, size=(500,400))

No hay ningún problema si queremos representar otras variables:

plot(sol, idxs=[ysp, y, c], size=(500,400))

Para acabar vamos a comparar las diferentes respuestas para los siguientes valores de ganancia del controlador: 0.01, 0.1, 1 y 2.5.

Como el problema a resolver numéricamente prob ya está creado, no es necesario que repitamos el cálculo. Podemos modificar el problema con la instrucción remake. Como vamos a repetir la simulación para varios valores de \(K_c\), definiremos una función:

function simul(Kc_val)
    p = [Kc => Kc_val]
    newprob = remake(prob; p=p)
    solve(newprob, Tsit5())
end
simul (generic function with 1 method)

Con un bucle podemos representar todas las soluciones en un mismo gráfico para poder compararlas:

plot(size=(500,400)
)
for i in [0.01, 0.1, 1, 2.5]
    plot!(simul(i), idxs=[y], label=i)
end
plot!()

Finalmente comprobamos que no tenemos los problemas de estabilidad del método numérico cuando el tiempo es un poco largo:

plot(simul(3.5), size=(500,400))

Conclusiones

Hemos visto que ModelingToolkit nos permite simular un lazo de control de manera muy simple gracias al uso que hace del cálculo simbólico y de los mejores algoritmos disponibles en la resolución de ecuaciones diferenciales.

Hemos planteado el sistema de ecuaciones diferenciales, lo que no ha resultado difícil. Una manera alternativa, quizás más sencilla, es utilizar la ModelingToolkitStandardLibrary que nos permitirá de manera muy simple definir elementos de primer orden, segundo orden, controladores PID, etc., pero esto queda para otro día.