Como dibujar un diagrama de Bode con Julia

control de procesos
julia
Autor/a

runjaj

Fecha de publicación

5 de julio de 2023

Es frecuente recurrir a librerías especializadas en control de procesos, como puede ser Control System Toolbox de Matlab o ControlSystems.jl para representar los diagramas de Bode.

Vamos a ver que podemos representar estos diagramas utilizando Julia sin ninguna librería especializada. Esta entrada es la explicación larga de este toot:

¿Se puede dibujar un diagrama de #Bode en un toot con #julialang ?

using Plots
using DSP:unwrap!
G(s) = 3/(2s^2+s+1)*exp(-.5s)
x(ω) = real(G(im*ω))
y(ω) = imag(G(im*ω))
AR(ω) = √(x(ω)^2 + y(ω)^2)
φ(ω) = atan(y(ω)/x(ω))*180/pi
ω = 10.0.^(-1:0.001:1);
plot(ω, AR, xscale=:log10, yscale=:log10, xlabel="ω", ylabel="AR", legend=false)
φᵤ = unwrap!([φ(i) for i in ω]);
plot(ω, φᵤ, xscale=:log10, xlabel="ω", ylabel="φ", legend=false)

¡Se puede!

#control #procesos

Un ejemplo simple

Como ejemplo representaremos el diagrama de Bode de esta función de transferencia:

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

Empezamos definiendo la función de transferencia del ejemplo:

G(s) = (3s+1)/(s^2+4s+1)

A continuación hacemos el cambio de variable \(s=\mathrm{i}\omega\) y encontramos la parte real y la parte imaginaria:

\[\begin{align} x(\omega) &= Re(G(\mathrm{i}\omega))\\ y(\omega) &= Im(G(\mathrm{i} \omega)) \end{align}\]

donde \(\mathrm{i} = \sqrt{-1}\) y \(\omega\) es la frecuencia angular.

En Julia, simplemente tenemos que reproducir estas funciones:

x(ω) = real(G(im*ω))
y(ω) = imag(G(im*ω))

La razón de amplitudes es:

\[RA(\omega) = \sqrt{x^2(\omega)+y^2(\omega)}\]

y el desfase es:

\[\phi(\omega) = \arctan\frac{y(\omega)}{x(\omega)}\].

RA(ω) = sqrt(x(ω)^2 + y(ω)^2)
φ(ω) = atan(y(ω)/x(ω))*180/pi

Representaremos el diagrama de Bode para valores de \(\omega\) entre 10-2 y 102, calcularemos 100 puntos:

ω = 10 .^LinRange(-2, 2, 100)
100-element Vector{Float64}:
   0.010000000000000002
   0.010974987654930556
   0.012045035402587823
   0.013219411484660288
   0.014508287784959401
   0.015922827933410922
   0.017475284000076828
   0.019179102616724886
   0.0210490414451202
   0.023101297000831605
   0.025353644939701114
   0.02782559402207126
   0.030538555088334154
   ⋮
  35.93813663804626
  39.442060594376564
  43.28761281083057
  47.50810162102798
  52.14008287999685
  57.2236765935022
  62.80291441834253
  68.92612104349695
  75.64633275546291
  83.02175681319744
  91.11627561154896
 100.0

Ya solo queda representar la razón de amplitudes frente a la frecuencia angular. Normalmente utilizo la librería Plots.jl para representar gráficos, pero quiero probar Makie.jl:

using CairoMakie

fig_RA = Figure()
ax_RA = Axis(fig_RA[1,1],
    xlabel = "ω", xscale = log10,
    ylabel = "RA", yscale = log10,
    xminorgridvisible=true, xminorticksvisible=true,
    xminorticks = IntervalsBetween(10),
    yminorgridvisible=true, yminorticksvisible=true,
    yminorticks = IntervalsBetween(10))
lines!(ax_RA, ω, RA.(ω))
fig_RA

A continuación representaremos el desfase:

fig_φ = Figure()
ax_φ = Axis(fig_φ[1,1],
    xlabel = "ω", xscale = log10,
    ylabel = "φ",
    xminorgridvisible=true, xminorticksvisible=true,
    xminorticks = IntervalsBetween(10),
    yminorgridvisible=true, yminorticksvisible=true,
    yminorticks = IntervalsBetween(4))
lines!(ax_φ, ω, φ.(ω))
fig_φ

Normalmente se dibujan las dos figuras juntas:

fig = Figure()
ax_RA = Axis(fig[1,1],
    xlabel = "ω", xscale = log10,
    ylabel = "RA", yscale = log10,
    xminorgridvisible=true, xminorticksvisible=true,
    xminorticks = IntervalsBetween(10),
    yminorgridvisible=true, yminorticksvisible=true,
    yminorticks = IntervalsBetween(10))
ax_φ = Axis(fig[2,1],
    xlabel = "ω", xscale = log10,
    ylabel = "φ",
    xminorgridvisible=true, xminorticksvisible=true,
    xminorticks = IntervalsBetween(10),
    yminorgridvisible=true, yminorticksvisible=true,
    yminorticks = IntervalsBetween(4))
lines!(ax_RA, ω, RA.(ω))
lines!(ax_φ, ω, φ.(ω))
fig

Solución a los dolores de cabeza del arco tangente

Veamos un caso más complejo:

G(s) = 20*(s+1)/(s*(s+5)*(s^2+2s+10))

x(ω) = real(G(im*ω))
y(ω) = imag(G(im*ω))

RA(ω) = sqrt(x(ω)^2 + y(ω)^2)
φ(ω) = atan(y(ω)/x(ω))*180/pi

ω = 10 .^LinRange(-2, 2, 1000)

fig = Figure()
ax_RA = Axis(fig[1,1],
    xlabel = "ω", xscale = log10,
    ylabel = "RA", yscale = log10,
    xminorgridvisible=true, xminorticksvisible=true,
    xminorticks = IntervalsBetween(10),
    yminorgridvisible=true, yminorticksvisible=true,
    yminorticks = IntervalsBetween(10))
ax_φ = Axis(fig[2,1],
    xlabel = "ω", xscale = log10,
    ylabel = "φ",
    xminorgridvisible=true, xminorticksvisible=true,
    xminorticks = IntervalsBetween(10),
    yminorgridvisible=true, yminorticksvisible=true,
    yminorticks = IntervalsBetween(4))
lines!(ax_RA, ω, RA.(ω))
lines!(ax_φ, ω, φ.(ω))
fig

Tenemos un problema y es que para un valor de \(\omega\) de un poco más de 2 el desfase ha dado una vuelta, lo que provoca esa discontinuidad. Lo más sencillo para solucionar este problema es utilizar la función unwrap!() de la librería DSP.jl.

using DSP:unwrap!

φᵤ = unwrap!([φ(i) for i in ω])

fig_φ = Figure()
ax_φ = Axis(fig_φ[1,1],
    xlabel = "ω", xscale = log10,
    ylabel = "φ",
    xminorgridvisible=true, xminorticksvisible=true,
    xminorticks = IntervalsBetween(10),
    yminorgridvisible=true, yminorticksvisible=true,
    yminorticks = IntervalsBetween(4))
lines!(ax_φ, ω, φᵤ)
fig_φ

Es importante destacar que hemos aumentado en número de puntos para los que hemos calculado el desfase de 100 a 1000 para que el algoritmo detecte correctamente las vueltas.

Una función para no escribir tanto

Podemos crear una función que reuna todo el código anterior:

function bode(G; w_min=1e-2, w_max=1e2, puntos=100)
    x(ω) = real(G(im*ω))
    y(ω) = imag(G(im*ω))
    
    RA(ω) = sqrt(x(ω)^2 + y(ω)^2)
    φ(ω) = atan(y(ω)/x(ω))*180/pi
    
    ω = 10 .^LinRange(log10(w_min), log10(w_max), puntos)
    
    φᵤ = unwrap!([φ(i) for i in ω])
    
    fig = Figure()
    ax_RA = Axis(fig[1,1],
        xlabel = "ω", xscale = log10,
        ylabel = "RA", yscale = log10,
        xminorgridvisible=true, xminorticksvisible=true,
        xminorticks = IntervalsBetween(10),
        yminorgridvisible=true, yminorticksvisible=true,
        yminorticks = IntervalsBetween(10))
    ax_φ = Axis(fig[2,1],
        xlabel = "ω", xscale = log10,
        ylabel = "φ",
        xminorgridvisible=true, xminorticksvisible=true,
        xminorticks = IntervalsBetween(10),
        yminorgridvisible=true, yminorticksvisible=true,
        yminorticks = IntervalsBetween(4))
    lines!(ax_RA, ω, RA.(ω))
    lines!(ax_φ, ω, φᵤ)
    fig
end
bode (generic function with 1 method)

Comprobemos que funciona:

G(s) = 3/(2s^2+s+1)*exp(-.5s)

bode(G; w_min=0.1, w_max=10, puntos=500)

Conclusión

Como vemos crear una función que dibuje los diagramas de Bode (o de Nyquist, si nos apetece) no es complicado. Julia nos permite aplicar la teoría de manera prácticamente directa, ¡lo que es fantástico!

A partir de hay podemos crear nuevas funciones para crear nuestra propia caja de herramientas. ¿Ventajas? Conocemos exactamente lo que hace el código, además de aprender. ¿Inconvenientes? Es fácil caer en soluciones naïf que no tengan en cuenta algunos casos o que el código sea un tanto “chapucero”.