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

NOTA: O código-fonte desse tutorial não irá funcionar nas últimas versões de Julia e OrdinaryDiffEq.jl. Você poderá ver uma versão atualizada do código-fonte aqui.

Olá!

Nesse post, eu mostrarei como podemo utilizar um Data Array para estender as possibilidade do pacote OrdinaryDiffEq 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:

https://www.ronanarraes.com/2017/02/utilizando-julia-para-simular-sistemas-compostos-de-partes-continuas-e-discretas/

Definição do problema

No meu último post, eu mostrei como o pacote OrdinaryDiffEq 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, 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 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. 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

Antiga forma de se criar uma estrutura customizada no OrdinariDiffEq

Então, Chris Rackauckas decidiu criar uma nova interface dentro do pacote OrdinaryDiffEq 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.

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 DEDataArray:

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:

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

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

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 definindo corretamente 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 user_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:

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:

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:

ou

que deverá fornecer o seguinte resultado:

Conclusões

Nesse tutorial, eu mostrei como você pode utilizar a interface Data Array do pacote OrdinaryDiffEq 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ário. Eu irei utilizar o seu conselho para melhorar esse e meus futuros tutoriais!

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

Código fonte