Estendendo as Soluções de EDOs em Julia Através da Criação de Data Arrays Customizados para as Simulações

Tutorial /Julia
7 min.
Atualizado em:

Olá!

Nesse post, eu mostrarei como podemo utilizar um Data Array para estender as possibilidade do pacote OrdinaryDiffEq.jl da linguagem julia para simular sistemas dinâmicos.

Experiência necessária

Esse tutorial provê informações adicionais para resolver o problema introduzido no meu último tutorial , então é muito importante que você o leia antes.

Definição do problema

No meu último post , eu mostrei como o pacote OrdinaryDiffEq.jl da linguagem julia pode ser utilizado para simular sistema de controle que possuem partes discretas e contínuas. No entanto, existiam dois grande problemas:

  1. Eu precisei utilizar uma variável global para passar a informação da parte discreta para a função dinâmica (naquele caso, a informação era a aceleração);
  2. Eu teria que utilizar alguns hacks muito ruins para plotar a aceleração comandada (criar um vetor global para guardar os valores computados).

Como mencionei anteriormente, a solução ótima para contornar esses problemas é era de alguma forma difícil de se obter. Felizmente, Chris Rackauckas , que é o mantenedor do pacote OrdinaryDiffEq.jl, introduziu o conceito de Data Arrays, que simplifica muito a solução para aqueles problemas. Note que você precisará da v1.4.0 ou superior do pacote OrdinaryDiffEq para ter acesso à nova interface Data Array.

Nesse tutorial, o problema que queremos resolver é: como criar parâmetros customizados em que a parte discreta do sistema modifica e que influenciará na função dinâmica?

Uma pequena história

Quando eu comecei a usar o OrdinaryDiffEq.jl para simular o subsistema de controle e atitude do satélite Amazonia-1, eu verifiquei que eu precisaria de dois tipos de variáveis:

  1. Aquelas que compõe o vetor de estados e que serão integradas pelo solver; e
  2. Alguns parâmetros que o loop de controle pode modificar e que irá influenciar na solução dinâmica, como, por exemplo, o torque comandado para as rodas de reação e magneto-torquers.

Entretanto, naquela época, não existia um meio fácil de se conseguir isso. O único método é criar uma estrutura (type) customizada em julia para guardar os parâmetros que são externos ao vetor de estado. O problema é que você precisava implementar muitos métodos para que essa nova estrutura fosse compatível com a interface do OrdinaryDiffEq.jl. A seguir, você poderá ver o que deveria ser feito para que se definisse um parâmetro chamado f1. Esse código foi utilizado na issue aberta por mim no GitHub: https://github.com/JuliaDiffEq/DifferentialEquations.jl/issues/117

Método antigo para se criar um tipo customizado no OrdinaryDiffEq.jl

importall Base
import RecursiveArrayTools.recursivecopy!

using DifferentialEquations

type SimType{T} <: AbstractArray{T,1}
    x::Array{T,1}

    f1::T
end

function SimType{T}(::Type{T}, length::Int64)
    return SimType{T}(zeros(length), T(0))
end

function SimType{T}(S::SimType{T})
    return SimType(T, length(S.x))
end

function SimType{T}(S::SimType{T}, dims::Dims)
    return SimType(T, prod(dims))
end

similar{T}(A::SimType{T}) = begin
    R = SimType(A)
    R.f1 = A.f1
    R
end

similar{T}(A::SimType{T}, ::Type{T}, dims::Dims) = begin
    R = SimType(A, dims)
    R.f1 = A.f1
    R
end

done(A::SimType, i::Int64) = done(A.x,i)
eachindex(A::SimType)      = eachindex(A.x)
next(A::SimType, i::Int64) = next(A.x,i)
start(A::SimType)          = start(A.x)

length(A::SimType) = length(A.x)
ndims(A::SimType)  = ndims(A.x)
size(A::SimType)   = size(A.x)

recursivecopy!(B::SimType, A::SimType) = begin
    recursivecopy!(B.x, A.x)
    B.f1 = A.f1
end

getindex( A::SimType,    i::Int...) = (A.x[i...])
setindex!(A::SimType, x, i::Int...) = (A.x[i...] = x)

Então, Chris Rackauckas decidiu criar uma nova interface dentro do pacote OrdinaryDiffEq.jl para que você possa facilmente construir suas próprias estruturas customizadas para serem utilizadas nos solvers.

Nota: Na verdade, a definição da nova interface foi colocada no pacote DiffEqBase , que é uma dependência do OrdinaryDiffEq.jl.

Interface Data Type

Essa interface permite que você defina suas próprias estruturas a serem utilizadas pelo solver. Para fazer isso, você simplesmente precisa definir um novo tipo em julia (type) que herda o tipo DEDataVector:

mutable struct MyDataArray{T} <: DEDataVector{T}
    x::Vector{T}
    a::T
    b::Symbol
end

onde deve existir um array x que define o vetor de estados e qualquer outro parâmetro que você irá necessitar na sua simulação. Dado isso, você poderá acessar quaisquer parâmetros dentro do solver e das funções de callback utilizando u.a, u.b, etc. Você poderá também acessar o i-ésimo elemento do vetor de estados utilizando u[i]. Em outras palavras, você não precisará da notação pesada u.x[i].

Isso praticamente substitui todo aquele código que eu mostrei anteriormente 👍.

Resolvendo o nosso problema

Então, agora nós podemos remover a variável global a (🎉🎊🎉🎊🎉🎊) definindo o seguinte tipo:

mutable struct SimWorkspace{T} <: DEDataVector{T}
    x::Vector{T}
    a::T
end

Dentro de nossa função dinâmica, nós podemos agora obter a aceleração comandada utilizando u.a:

# Dynamic equation.
function dyn_eq!(du, u, t, p)
    du[1] = u[2]
    du[2] = u.a
    return nothing
end

Agora, nós precisamos apenas modificar a função efeito do nosso callback para escrever nessa variável:

# Affect function.
function control_loop!(integrator)
    p = integrator.u[1]
    v = integrator.u[2]

    a = k*(r-p) + d*(0.0-v)

    for c in full_cache(integrator)
        c.a = a
    end
end

Você deverá ter notado um for bem estranho aqui. Esse problema é criado pela estrutura interna do solucionadores de equações diferenciais ordinárias. Por exemplo, o Algoritmo Runge-Kutta de 4ª ordem necessita quatro caches chamados k_1, k_2, k_3 e k_4, que são basicamente a função dinâmica computada em diferentes pontos do espaço de estados, sendo utilizados para obter a solução no próximo instante de amostragem. Nesse caso, o parâmetro a deve ser definido em todos esses caches para fornecer a solução correta. Isso é precisamente o que aquele for está fazendo. Ele está passando por todos os caches no solver e definindocorretamente o parâmetro a com a aceleração computada pela função efeito desse callback. Note que o vetor de estado atual integrator.u está dentro da lista full_cache(integrator) e, consequentemente, é também definido dentro daquele loop. O lado bom de se utilizar esse for é que você não precisa se preocupar com quantos caches o solver selecionado possui.

Finalmente, nós precisamos apenas definir nosso estado inicial utilizando o tipo SimWorkspace:

u0 = SimWorkspace(zeros(2), 0.0)

onde o primeiro parâmetro está inicializando o vetor de estados com zeros (posição e velocidade) e o segundo parâmetro está inicializando a aceleração. Após a execução do programa, nós podemos acessar a aceleração comandada em cada instante de amostragem que foi salvo utilizando sol[i].a.

Se você estiver utilizando o PyPlot , então o vetor de estado e a aceleração podem ser plotados através de:

plot(sol.t, sol.u, sol.t, map(x->x.a, sol.u))

que deverá fornecer o seguinte resultado:

Se você quiser alterar o estado inicial, então você deverá apenas modificar a variável u0 antes de chamar o solver:

u0 = SimWorkspace([5.0; -1.0], 0.0)

ou

u0 = SimWorkspace(zeros(2), 0.0)
u0[1] = 5.0
u0[2] = -1.0

que deverá fornecer o seguinte resultado:

Conclusões

Nesse tutorial, eu mostrei como você pode utilizar a interface Data Array do pacote OrdinaryDiffEq.jl para criar espaços de trabalho customizados para suas simulaçãoes. Portanto, agora você sabe como criar parâmetros para serem utilizados no solver que não pertencem ao vetor de estado.

Eu forneci um exemplo através da extensão da solução apresentada em um tutorial anterior , onde nós precisamos de uma variável global para lidar com esse tipo de parâmetro. Note que, embora o exemplo apresentado aqui seja muito simples, esse conceito pode ser facilmente estendido para lhe fornecer uma solução simples e elegante para simular sistemas muito mais complexos.

Se não ficou claro, por favor, me avise nos comentários. Eu irei utilizar o seu conselho para melhorar esse e meus futuros tutoriais!

A seguir, você encontrará o código fonte desse tutorial.

using DiffEqBase
using OrdinaryDiffEq

###############################################################################
#                                 Structures
###############################################################################

mutable struct SimWorkspace{T} <: DEDataVector{T}
    x::Vector{T}
    a::T
end

###############################################################################
#                                  Variables
###############################################################################

# Parameters.
r  = 10.0   # Reference position.
k  = 0.3    # Proportional gain.
d  = 0.8    # Derivative gain.
vl = 2.0    # Velocity limit.

# Configuration of the simulation.
tf     = 30.0             # Final simulation time.
tstops = collect(0:1:tf)  # Instants that the control loop will be computed.

# Initial state of the simulation.
#                 |-- State Vector --| |-- Parameter a -- |
u0 = SimWorkspace(      zeros(2),               0.0        )

###############################################################################
#                                  Functions
###############################################################################

# Dynamic equation.
function dyn_eq(du, u, t, p)
    du[1] = u[2]
    du[2] = u.a
    return nothing
end

# CALLBACK: Control loop.
# =======================

# Condition function.
function condition_control_loop(u, t, integrator)
    (t in tstops)
end

# Affect function.
function control_loop!(integrator)
    p = integrator.u[1]
    v = integrator.u[2]

    a = k*(r-p) + d*(0.0-v)

    for c in full_cache(integrator)
        c.a = a
    end
end

cb = DiscreteCallback(condition_control_loop, control_loop!)

###############################################################################
#                                Solver
###############################################################################

prob = ODEProblem(dyn_eq, u0, (0.0, tf))
sol = solve(prob, Tsit5(), callback = cb, tstops = tstops)
#Julia
#OrdinaryDiffEq.jl
#Simulação
comments powered by Disqus