Meu fluxo de trabalho em Julia: Legibilidade VS. Desempenho

Olá!

Após o meu último post, recebi muitas perguntas como: “O que é Julia? Por que você está usando Julia? Vale a pena?”

Bem, vou tentar responder as duas últimas perguntas mostrando o meu fluxo de trabalho: como eu codifico e que tipo de ganhos eu obtenho (para a primeira pergunta eu vou indicar o website oficial 🙂 ).

Requisitos necessários

É necessário um pouco de conhecimento de C e Julia para que você possa entender totalmente este post.

Meu caso de uso

Muitas das minhas atividades no INPE exigem análises que só podem ser realizadas por algoritmos de otimização ou simulações de Monte Carlo. Por exemplo, durante os estudos da Pré-Fase A no Centro de Projeto Integrado de Missões Espaciais (CPRIME), às vezes utilizamos algoritmos evolutivos para obter as melhores soluções para o projeto de um satélite. Também é comum utilizar simulações de Monte Carlo para encontrar os piores cenários possíveis, por exemplo, o estado inicial do satélite (atitude e velocidade angular) que produz o maior tempo para adquirir o apontamento inicial para o Sol. Assim, alguns algoritmos que codificamos devem ser executados milhares, até milhões, de vezes. Então, o desempenho é algo que realmente importa para mim.

Por outro lado, somos principalmente engenheiros que estão acostumados com o MATLAB. Os algoritmos são escritos para resolver equações e simular sistemas dinâmicos. Portanto, quanto mais próximo o algoritmo estiver das equações matemáticas, melhor.

A seguir, descreverei meu fluxo de trabalho: como eu escrevo o código em Julia desde a validação inicial até a versão que deve ser utilizada dentro de um otimizador. Eu tentei muitas abordagens ao longo dos anos, mas esta foi a que forneceu o menor tempo para obter os resultados.

Exemplo numérico

Para ilustrar, vamos supor que queremos um algoritmo para executar a atualização de covariância do filtro de Kalman de um sistema dinâmico. Não dê muita importância a essas equações, apenas selecionei um cenário aleatório para ilustrar melhor o que estamos fazendo.

Inicialização

\(\begin{aligned}\mathbf{F}_{k-1} &= \exp\left[\begin{array}{c c c} \mathbf{w}^{\times} & -\mathbf{I}_3 & \mathbf{0}_{3 \times 12} \\ \mathbf{0}_{3 \times 3} & \lambda \cdot \mathbf{I}_3 & \mathbf{0}_{3 \times 12} \\ \mathbf{0}_{12 \times 3} & \mathbf{0}_{12 \times 3} & \mathbf{0}_{12 \times 12 } \end{array} \right] \\ \mathbf{H}_k &= \left[ \begin{array}{c c c c c c} \mathbf{s}^{\times} & \mathbf{0}_{3 \times 3} & \mathbf{s}^{\times} & \mathbf{I}_3 & \mathbf{0}_{3 \times 3} & \mathbf{0}_{3 \times 3} \\ \mathbf{B}^{\times} & \mathbf{0}_{3 \times 3} & \mathbf{0}_{3 \times 3} & \mathbf{0}_{3 \times 3} & \mathbf{B}^{\times} & \mathbf{I}_3 \end{array} \right] \\ \mathbf{P}_p &= \mathbf{I}_{18} \\ \mathbf{P}_u &= \mathbf{I}_{18} \end{aligned}\)

Loop de simulação

Para \(k = 2 \) até \(60000\), faça:

\(\begin{aligned} \mathbf{P}_p &= \mathbf{F}_{k-1} \mathbf{P}_p \mathbf{F}_{k-1}^{T} + \mathbf{Q} \\ \mathbf{K}_k &= \mathbf{P}_p\mathbf{H}_k^T(\mathbf{R} + \mathbf{H}_{k} \mathbf{P}_p \mathbf{H}_{k}^{T}) \\ \mathbf{P}_u &= (\mathbf{I}_{18} – \mathbf{K}_k\mathbf{H}_k)\mathbf{P}_p(\mathbf{I}_{18}-\mathbf{K}_k\mathbf{H}_k)^T + \mathbf{K}_k\mathbf{R}\mathbf{K}_k^T \end{aligned} \)

Toda a vez que \(\mathbf{P}_u \) for modificada, guarde \(\mbox{trace}(\mathbf{P}_u) \).

Nota

\(\mathbf{w}^{\times} \), \(\mathbf{s}^{\times} \) e \(\mathbf{B}^{\times} \) são matrizes constantes, \(\mathbf{0}_{n \times m} \) é uma matriz \(n \times m\) preenchida com 0s, e \(\mathbf{I}_n\) é uma matrix identidade \(n \times n\).

Primeiro: escreva o código rápido, o mais legível possível

A primeira coisa é codificar uma versão inicial do algoritmo que você quer, sem se preocupar muito com o desempenho. Nesta fase, a legibilidade é fundamental, porque facilitará a depuração.

Nossa primeira versão do código será:

using LinearAlgebra

# Constants
# ==============================================================================

const O3x3   = zeros(3,3)
const O3x12  = zeros(3,12)
const O12x18 = zeros(12,18)

function kalman_filter()

    # Constants
    # ==========================================================================

    λ = 1/100
    Q = 1e-20I
    R = 1e-2I

    # Dynamic model
    # ==========================================================================

    wx = Float64[ 0 -1  0;
                  1  0  0;
                  0  0  0;]

    Ak_1 = [  wx    -I   O3x12;
             O3x3  λ*I   O3x12;
                  O12x18      ;]

    Fk_1 = exp(Ak_1)

    # Measurement model
    # ==========================================================================

    sx = Float64[ 0  0  0;
                  0  0 -1;
                  0  1  0;]

    Bx = Float64[ 0  0 -1;
                  0  0  0;
                  1  0  0;]

    Hk = [sx O3x3   sx    I O3x3 O3x3;
          Bx O3x3 O3x3 O3x3   Bx    I;]

    # Kalman filter initialization
    # ==========================================================================

    Pu = Matrix{Float64}(I,18,18)

    # Simulation
    # ==========================================================================

    result = Vector{Float64}(undef, 60000)

    result[1] = tr(Pu)

    for k = 2:60000
        Pp = Fk_1*Pu*Fk_1' + Q
        K  = Pp*Hk'*pinv(R + Hk*Pp*Hk')
        Pu = (I - K*Hk)*Pp*(I - K*Hk)' + K*R*K'

        result[k] = tr(Pu)
    end

    return result
end

É possível ver quão perto do equacionamento matemático esse código está? Mesmo os símbolos como λ podem ser utilizados em Julia, que suporta qualquer caractere UTF-8 no nome das variáveis. Agora, imagine algo com centenas de equações. É mais fácil validar o código e encontrar problemas quando você escrever usando uma linguagem computacional parecida com a linguagem matemática.

Informação

O I em Julia não é uma matriz identidade 3 x 3 como nas equações acima. Ele é chamado de operador UniformScalling. Em nosso código, ele se transforma em uma matriz identidade com a dimensão necessária para a operação fazer sentido, isto é 3 x 3.

Nesta primeira versão do código, temos duas tarefas: 1) ver se funciona corretamente; e 2) ver se o desempenho é bom o suficiente para as nossas necessidades. Sabemos que após 60.000 iterações, o traço da matriz de covariância \(P_u\) é aproximadamente 7,00089. Em Julia, podemos executar o código no REPL e ver se isso é verdade:

julia> result = kalman_filter();

julia> result[end]
7.000894875771719

Legal 🎊🎉! Pelo menos o código faz o que é necessário. O próximo passo é medir o desempenho e verificar se é o suficiente. Isso é muito importante: Julia é uma linguagem muito rápida para mim não porque é mais rápida que C, mas porque o tempo entre a escrita do código inicial e o resultado final é muito menor. Às vezes você simplesmente não precisa de um código com grande desempenho e quer codificar o mais rápido possível, Julia pode fazer isso conforme foi mostrado. Por outro lado, algumas vezes você realmente precisa de velocidade, Julia pode fazer isso também como será mostrado a seguir.

A maneira correta de avaliar o desempenho de um código em Julia é colocá-lo dentro de uma função (já fizemos isso) e usar o excelente pacote BenchmarkTools.jl:

julia> using BenchmarkTools

julia> @btime kalman_filter()
  1.296 s (2640271 allocations: 2.83 GiB)

Agora sabemos que o tempo médio de execução da função (após a etapa de compilação) é de 1,296s e que foi alocado 2,83 GiB de memória. Neste momento, eu me pergunto: “Essa velocidade é suficiente para o que eu preciso fazer?” Por exemplo, se eu precisar executar esse código apenas uma vez, posso conviver com seu estado atual. No entanto, se eu precisar que esse código seja utilizado dentro de um otimizador sendo executado 1.000.000 de vezes, eu realmente desejaria otimizá-lo.

Segundo: Modifique seu código para melhorar o desempenho, se necessário

Há algumas coisas que você precisa tomar cuidado para aumentar o desempenho de um código em Julia. A documentação é muito boa e eu realmente sugiro que você a leia. Para começar, você deve se preocupar com três coisas:

  • Estabilidade do tipo: Julia deve saber o tipo das variáveis no código a partir apenas do tipo das entradas. Isso pode ser verificado com a macro @code_warntype.
  • Alocações: o processo para alocar memória é lento, você deve realizar o mínimo de alocações possíveis para melhorar o desempenho.
  • Evite cálculos duplicados: os compiladores estão ficando cada vez mais inteligentes e às vezes eles conseguem identificar cálculos duplicados, mas não faz mal ajudá-los!

Podemos ver facilmente que nosso código é tipo-estável (type stable) ao chamar @code_warntype kalman_filter(). Você verá que nenhuma variável pode ter tipos distintos. Eu não vou me aprofundar aqui porque este é um assunto para um post inteiro. Para maiores informações, consulte a documentação.

Existem alguns métodos para verificar quais partes do código estão alocando. O método “oficial” é iniciar Julia com a opção --track-allocation=user, executar a função e, em seguida, analisar o arquivo .mem que será criado. Isso é descrito aqui. No entanto, existe um método que eu, pessoalmente, acho melhor. Ele é baseado no pacote TimerOutputs.jl. Esse pacote fornece uma maneira de analisar apenas linhas específicas. No nosso caso, devemos nos preocupar com o código executado dentro do loop for. O arquivo anterior deve ser modificado da seguinte maneira:

using LinearAlgebra
using TimerOutputs

...

    reset_timer!()
    for k = 2:60000
        @timeit "Pp" Pp = Fk_1*Pu*Fk_1' + Q
        @timeit "K"  K  = Pp*Hk'*pinv(R + Hk*Pp*Hk')
        @timeit "Pu" Pu = (I - K*Hk)*Pp*(I - K*Hk)' + K*R*K'

        result[k] = tr(Pu)
    end

    print_timer()

    return result
end

Agora, quando você executar kalman_filter();, você verá na tela:

 ──────────────────────────────────────────────────────────────────
                           Time                   Allocations
                   ──────────────────────   ───────────────────────
 Tot / % measured:      2.01s / 69.0%           2.86GiB / 99.2%

 Section   ncalls     time   %tot     avg     alloc   %tot      avg
 ──────────────────────────────────────────────────────────────────
 K          60.0k    673ms  48.5%  11.2μs   1.08GiB  38.3%  19.0Ki
 Pu         60.0k    511ms  36.9%  8.52μs   1.29GiB  45.4%  22.5KiB
 Pp         60.0k    203ms  14.7%  3.39μs    472MiB  16.3%  8.06KiB
 ─────────────────────────────────────────────────────────────────

OK, agora temos algo para analisar. Vamos realizar quatro modificações:

  • Essas três variáveis Pp, K, e Pu estão sendo alocadas a cada iteração. Uma coisa que devemos fazer é alocá-las antes do loop e dizer para Julia utilizar o mesmo espaço para armazenar os novos valores. A última ação pode ser feita utilizando o operador .= que executa um broadcasted assignment. Em outras palavras, cada elemento da variável do lado esquerdo será igual ao elemento correspondente da expressão do lado direito, o que não realiza alocações. Também é bom substituir o operador + pelo operador .+ pela mesma razão, uma vez que a soma de matrizes é uma operação executada elemento por elemento.
  • A variável UniformScaling será trocada por matrizes pré-alocadas. Isso também ajudou a reduzir as alocações.
  • Note que (I - K*Hk) está sendo calculado duas vezes dentro do loop. Nós usaremos uma variável auxiliar para evitar isso. É importante que essa variável auxiliar também seja pré-alocada antes do loop pela mesma razão que discutimos no primeiro ponto.
  • Como podemos garantir que todos os elementos das matrizes dentro do for são acessados por índices válidos, então podemos utilizar a macro @inbounds como descrito na documentação. Isso evita operações de verificação de dimensões, melhorando o desempenho.

Dado isto, nosso código agora deve ser:

using LinearAlgebra

# Constants
# ==============================================================================

const O3x3   = zeros(3,3)
const O3x12  = zeros(3,12)
const O12x18 = zeros(12,18)

function kalman_filter()

    # Constants
    # ==========================================================================

    I6  = Matrix{Float64}(I,6,6)
    I18 = Matrix{Float64}(I,18,18)
    λ   = 1/100
    Q   = 1e-20I18
    R   = 1e-2I6

    # Dynamic model
    # ==========================================================================

    wx = Float64[ 0 -1  0;
                  1  0  0;
                  0  0  0;]

    Ak_1 = [  wx     -I   O3x12;
             O3x3   λ*I   O3x12;
                   O12x18     ;]

    Fk_1 = exp(Ak_1)

    # Measurement model
    # ==========================================================================

    sx = Float64[ 0  0  0;
                  0  0 -1;
                  0  1  0;]

    Bx = Float64[ 0  0 -1;
                  0  0  0;
                  1  0  0;]

    Hk = [sx O3x3   sx    I O3x3 O3x3;
          Bx O3x3 O3x3 O3x3   Bx    I;]

    # Kalman filter initialization
    # ==========================================================================

    Pu = Matrix{Float64}(I,18,18)
    Pp = similar(Pu)

    # Simulation
    # ==========================================================================

    result = Vector{Float64}(undef, 60000)

    result[1] = tr(Pu)

    # Auxiliary variables.
    K    = zeros(18,6)
    aux1 = zeros(18,18)

    @inbounds for k = 2:60000
        Pp   .= Fk_1*Pu*Fk_1' .+ Q
        K    .= Pp*Hk'*pinv(R .+ Hk*Pp*Hk')
        aux1 .= I18 .- K*Hk
        Pu   .= aux1*Pp*aux1' .+ K*R*K'

        result[k] = tr(Pu)
    end

    return result
end

O benchmark dessa nova versão é:

julia> @btime kalman_filter();
  1.109 s (2340283 allocations: 2.06 GiB)

Acabamos de reduzir o tempo computacional em 14% com apenas pequenas modificações! Se essa função for utilizada em um otimizador que a executa 1.000.000 de vezes, economizaremos mais de 2 dias.

Terceiro: Indo ao extremo!

Se é realmente necessário reduzir o tempo computacional, então ainda podemos melhorar algumas coisas aqui. No entanto, alerta de spoiler, isso prejudicará a legibilidade. Por isso, aconselho seguir este caminho somente depois que seu código estiver testado e validado e se você realmente precisar de mais velocidade.

Olhando para o último benchmark, ainda vemos que o código está alocando 2GiB. Isso vem das matrizes temporárias que Julia deve alocar para armazenar multiplicações intermediárias. Por exemplo, na multiplicação Fk_1*Pu*Fk_1' .+ Q, é necessário alocar pelo menos um array temporário para armazenar Fk_1*Pu ou Pu*Fk_1'. O compilador não pode (ainda?) adivinhar que isso pode ser alocado antes do loop ao invés de alocá-lo a cada iteração. Portanto, se realmente precisamos reduzir o tempo de computação, devemos reduzir drasticamente essas alocações.

A abordagem é bastante semelhante para a segunda fase: criar matrizes temporárias fora do loop e usá-las para armazenar as operações intermediárias. No entanto, aqui não podemos usar o operador de broadcast de multiplicação .*, pois ele executará uma multiplicação de elementos ao invés de um produto de matrizes. Dessa forma, a solução é usar a seguinte função:

mul!(Y, A, B)

que calcula o produto matricial A*B e armazena na matriz Y. Essa função também retorna a matriz Y para facilitar a concatenação de operações. Finalmente, o código modificado é:

using LinearAlgebra

# Constants
# ==============================================================================

const O3x3   = zeros(3,3)
const O3x12  = zeros(3,12)
const O12x18 = zeros(12,18)

function kalman_filter()

    # Constants
    # ==========================================================================

    I6  = Matrix{Float64}(I,6,6)
    I18 = Matrix{Float64}(I,18,18)
    λ   = 1/100
    Q   = 1e-20I18
    R   = 1e-2I6

    # Dynamic model
    # ==========================================================================

    wx = Float64[ 0 -1  0;
                  1  0  0;
                  0  0  0;]

    Ak_1 = [  wx     -I   O3x12;
             O3x3   λ*I   O3x12;
                   O12x18     ;]

    Fk_1 = exp(Ak_1)

    # Measurement model
    # ==========================================================================

    sx = Float64[ 0  0  0;
                  0  0 -1;
                  0  1  0;]

    Bx = Float64[ 0  0 -1;
                  0  0  0;
                  1  0  0;]

    Hk = [sx O3x3   sx    I O3x3 O3x3;
          Bx O3x3 O3x3 O3x3   Bx    I;]

    # Kalman filter initialization
    # ==========================================================================

    Pu = Matrix{Float64}(I,18,18)
    Pp = similar(Pu)

    # Simulation
    # ==========================================================================

    result = Vector{Float64}(undef, 60000)

    result[1] = tr(Pu)

    # Auxiliary variables.
    K    = zeros(18,6)
    aux1 = zeros(18,18)
    aux2 = zeros(18,18)
    aux3 = zeros(18,6)
    aux4 = zeros(6,6)
    aux5 = zeros(18,18)
    aux6 = zeros(18,18)

    for k = 2:60000
        # Pp = Fk_1*Pu*Fk_1' + Q
        mul!(aux2, mul!(aux1, Fk_1, Pu), Fk_1')
        @. Pp = aux2 + Q

        # K = Pp*Hk'*pinv(R + Hk*Pp*Hk')
        mul!(aux4, Hk, mul!(aux3, Pp, Hk'))
        mul!(K, aux3, pinv(R + aux4))

        # Pu = (I - K*Hk)*Pp*(I - K*Hk)' + K*R*K'
        mul!(aux1, K, Hk)
        @. aux2 = I18 - aux1
        mul!(aux6, mul!(aux5, aux2, Pp), aux2')
        mul!(aux5, mul!(aux3, K, R), K')
        @. Pu = aux6 + aux5

        result[k] = tr(Pu)
    end

    return result
end

Bem … isso é … feio 😅 No entanto, vamos fazer um benchmark:

julia> @btime kalman_filter()
  848.334 ms (1680299 allocations: 915.17 MiB)

Wow! Esse código é 23% mais rápido que a versão anterior e 34% mais rápido que a primeira versão. Comparado com esta última, a nova versão nos pouparia mais de 5 dias se executássemos o algoritmo 1.000.000 de vezes. O trade-off aqui é ter um código menos legível, matematicamente falando.

Informação

Julia tem excelentes capacidades de meta-programação. Pode-se facilmente construir uma macro que irá expandir as multiplicações para usar as variáveis auxiliares, mantendo a legibilidade. No entanto, este é um tópico um pouco avançado para cobrir por agora.

OK Ronan, e quanto ao C ++?

Essa é a pergunta mais frequente que me fazem quando falo de Julia. Naturalmente, o C ++ pode ser mais rápido no tempo de execução. No entanto, como Julia é uma linguagem dinâmica (parece interpretada), é mais fácil codificar e testar os resultados. Essa é uma das razões pelas quais linguagens como MATLAB e Python foram criadas. Assim, se você considerar o tempo total (desenvolvimento, teste e execução), na minha experiência, você pode vencer o C++ usando Julia (a menos que você seja um mágico em C++, o que a maioria dos engenheiros que eu conheço, inclusive eu, não são).

De qualquer forma, vamos apenas escrever o mesmo algoritmo em C++ para ver o quão longe estamos em termos de desempenho. Observe que estas são minhas habilidades em C++. Tenho quase certeza de que esse código pode ser otimizado (alguém pode fazer isso pra mim, por favor?). No entanto, eu não sou um programador experiente em Julia também. Assim, tenho a sensação de que o código em Julia também pode ser otimizado (alguém, também, pode fazer isso pra mim, por favor?).

A versão em C++ utiliza a biblioteca Eigen para as operações matriciais. Caso contrário, será um código muito, muito longo. O código pode ser visto a seguir:

#include <iostream>
#include "Eigen/Core"
#include "Eigen/QR"
#include "unsupported/Eigen/MatrixFunctions"

using namespace Eigen;

int main(int argc, char *argv[])
{
    double l = 1.0/100.0;
    Matrix<double, 18, 18> Q;
    Matrix<double,  6,  6> R;
    Matrix<double, 18, 18> Ak_1;
    Matrix<double, 18, 18> Fk_1;
    Matrix<double,  6, 18> Hk;
    Matrix<double, 18, 18> Pp;
    Matrix<double, 18, 18> Pu;
    Matrix<double, 18,  6> K;

    VectorXd result(60000);

    // Auxiliary variables.
    MatrixXd I6(6,6);
    MatrixXd I18(18,18);

    MatrixXd aux1(6,6);
    MatrixXd aux2(18,18);
    MatrixXd pinv(6,6);

    // Constants
    // =========================================================================

    I6  = MatrixXd::Identity(6,6);
    I18 = MatrixXd::Identity(18,18);
    Q   = I18*1e-20;
    R   = I6*1e-2;

    // Dynamic model
    // =========================================================================

    Ak_1 << 0, -1,  0, -1,  0,  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
            1,  0,  0,  0, -1,  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
            0,  0,  0,  0,  0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
            0,  0,  0,  l,  0,  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
            0,  0,  0,  0,  l,  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
            0,  0,  0,  0,  0,  l, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
            0,  0,  0,  0,  0,  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
            0,  0,  0,  0,  0,  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
            0,  0,  0,  0,  0,  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
            0,  0,  0,  0,  0,  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
            0,  0,  0,  0,  0,  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
            0,  0,  0,  0,  0,  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
            0,  0,  0,  0,  0,  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
            0,  0,  0,  0,  0,  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
            0,  0,  0,  0,  0,  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
            0,  0,  0,  0,  0,  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
            0,  0,  0,  0,  0,  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
            0,  0,  0,  0,  0,  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;

    // WARNING: Matrix exponential is unsupported in Eigen 3.3.7.
    Fk_1 = Ak_1.exp();

    // Measurement model
    // =========================================================================

    Hk << 0,  0,  0,   0, 0, 0,  0,  0,  0,   1, 0, 0,   0,  0,  0,  0, 0, 0,
          0,  0, -1,   0, 0, 0,  0,  0, -1,   0, 1, 0,   0,  0,  0,  0, 0, 0,
          0,  1,  0,   0, 0, 0,  0,  1,  0,   0, 0, 1,   0,  0,  0,  0, 0, 0,
          0,  0, -1,   0, 0, 0,  0,  0,  0,   0, 0, 0,   0,  0, -1,  1, 0, 0,
          0,  0,  0,   0, 0, 0,  0,  0,  0,   0, 0, 0,   0,  0,  0,  0, 1, 0,
          1,  0,  0,   0, 0, 0,  0,  0,  0,   0, 0, 0,   1,  0,  0,  0, 0, 1;

    // Kalman filter initialization
    // =========================================================================

    Pp.setIdentity();
    Pu.setIdentity();

    // Simulation
    // =========================================================================

    result[0] = Pu.trace();

    for(int k = 2; k <= 60000; k++) {
        Pp   = Fk_1*Pu*Fk_1.transpose() + Q;
        aux1 = R + Hk*Pp*Hk.transpose();
        pinv = aux1.completeOrthogonalDecomposition().pseudoInverse();
        K    = Pp*Hk.transpose()*pinv;
        aux2 = I18 - K*Hk;
        Pu   = aux2*Pp*aux2.transpose() + K*R*K.transpose();

        result[k-1] = Pu.trace();
    }

    return 0;
}

Na minha máquina, o tempo médio de execução utilizando -O3 é:

$ time ./a.out
./a.out  0,83s user 0,00s system 99% cpu 0,836 total

que é muito, muito próximo da versão mais rápida do código em Julia.

Informação

Alguns podem argumentar que seria mais rápido se eu usasse -march=native. De fato, é! No entanto, para uma comparação justa, eu teria também que recompilar a imagem de sistema de Julia com a mesma flag, mas, infelizmente, não tenho tempo para fazê-lo agora. Este é um procedimento fácil se você usar o pacote PackageCompiler.jl.

Conclusão

Espero que você tenha entendido como Julia está me ajudando a atingir meus objetivos mais rapidamente. Agora, temos uma linguagem fácil para codificar e testar algoritmos e que, com pequenas adaptações, pode atingir o desempenho de C. Eu não tenho mais necessidade de traduzir todos os algoritmos que precisam de alto desempenho para uma linguagem compilada, o que acaba sendo muito demorado para um engenheiro que não possui conhecimento avançadíssimo de C / C++.

O pacote SatelliteToolbox.jl para Linguagem Julia

Olá!

Neste post, gostaria de apresentar o SatelliteToolbox.jl, que é um pacote para a linguagem Julia com muitas opções para analisar missões espaciais. Ele é usado diariamente no Instituto Nacional de Pesquisas Espaciais (INPE). Primeiramente, apresento um breve histórico sobre o pacote e, em seguida, mostro algumas análises interessantes que podem ser feitas com ele.

Histórico

Em 2013, ingressei no INPE como Engenheiro de Sistemas Espaciais Júnior. Fui designado para a divisão de sistemas espaciais onde tive que trabalhar com o subsistema de controle de atitude e órbita (AOCS). Como eu tinha apenas um conhecimento intermediário sobre órbitas, decidi me aprofundar nesse assunto codificando algoritmos e comparando os resultados com a herança do INPE e com a literatura no meu tempo livre.

O primeiro passo foi selecionar a linguagem! Em meu Doutorado, utilizei o MATLAB para simular um sistema de navegação inercial, mas as simulações de Monte Carlo eram tão lentas que precisei reescrever muitas partes em C usando CMEX. Por outro lado, na minha pesquisa de pós-doutorado, onde estudei estimação em sistemas distribuídos não-lineares, decidi usar FORTRAN (com os padrões de 2008 para ter pelo menos um código legível …) para que a velocidade de execução não fosse um problema. Sim, o desempenho foi muito bom, mas demorei muito tempo para codificar. Então ouvi falar de uma nova linguagem que prometia o melhor dos dois mundos: algo que se assemelhava a uma linguagem interpretada com a velocidade de uma linguagem compilada! E foi assim que conheci Julia.

Naquela época (usando a v0.2, eu acho), Julia era uma linguagem realmente nova. Mas decidi aceitar as dificuldades e tentar codificar meus algoritmos com ela. De qualquer forma, era apenas um pequeno projeto pessoal para aprender mais sobre órbitas. Eu enfrentei muitos bugs e até tive que usar o branch master (pre-v0.3) devido a alguns problemas e funcionalidades que não existiam, mas foi divertido 🙂

Depois de alguns anos (e várias re-escritas devido a mudanças na linguagem), a v0.4 da linguagem Julia foi lançada. Neste momento, dada a quantidade de código que eu tinha e o estado da linguagem, comecei a ver que esse monte de algoritmos poderia, de fato, ser utilizado para auxiliar minhas atividades no INPE. Por isso, decidi criar um pacote privado, chamado SatToolbox.jl, para organizar tudo o que fiz.

Depois de algum tempo, esse pequeno projeto pessoal acabou se tornando o núcleo de um simulador do conceito operacional de missões espaciais chamado Forplan, que está sendo desenvolvido no Centro de Projeto Integrado de Missões Espaciais (CPRIME) do INPE. Dado o bom feedback que recebi, decidi renomear o pacote para SatelliteToolbox.jl e lançá-lo como um pacote oficial da linguagem Julia em março de 2018.

Neste post, gostaria de descrever brevemente o SatelliteToolbox.jl e como ele pode ser utilizado para algumas análises relativas a missões espaciais. Todos os recursos disponíveis podem ser vistos na documentação (disponível apenas em inglês). Uma lista breve dos algoritmos implementados até o momento deste post (na v0.5.0) é:

  • Modelos atmosféricos terrestre:
  • Modelos do campo geomagnético:
  • Índices espaciais:
    • Capacidade de se obter automaticamente diversos índices espaciais, como F10.7, Ap, Kp, etc.
  • Funções para realizar análises gerais relacionadas com órbitas, como conversão de anomalias, cálculo de perturbações, etc.
  • Propagadores orbitais:
    • Two body;
    • J2;
    • J4; e
    • SGP4/SDP4.
  • Funções para converter entre sistemas de referência ECI e ECEF:
    • Toda a teoria IAU-76/FK5 é suportada. Então, a conversão entre quaisquer sistemas de referência a seguir está disponível:
      • ITRF: International Terrestrial Reference Frame;
      • PEF: Pseudo-Earth Fixed reference frame;
      • MOD: Mean-Of-Date reference frame;
      • TOD: True-Of-Date reference frame;
      • GCRF: Geocentric Celestial Reference Frame;
      • J2000: J2000 reference frame;
      • TEME: True Equator, Mean Equinox reference frame.
    • Toda a teoria IAU-2006/2010 é suportada. Então, a conversão entre quaisquer sistemas de referência a seguir está disponível:
      • ITRF: International Terrestrial Reference Frame;
      • TIRS: Terrestrial Intermediate Reference Frame;
      • CIRS: Celestial Intermediate Reference Frame;
      • GCRF: Geocentric Celestial Reference Frame.
  • Funções para converter entre referências Geocêntricas e Geodésicas (WGS-84).

A seguir, eu forneço alguns exemplos de como o SatelliteToolbox.jl pode ser utilizado para analisar missões espaciais.

Instalação

A primeira coisa a se fazer (assumindo que você já instalou Julia, que pode ser obtida aqui) é instalar o pacote. Isso pode ser feito digitando:

julia> using Pkg
julia> Pkg.add("SatelliteToolbox") 

Informação

Se o pacote já estiver instalado, então tenha certeza que você possui no mínimo a versão v0.6.0. Você pode atualizar todos os seus pacotes através do comando Pkg.update().

Para carregar o pacote, que deve ser feito toda a vez que Julia é iniciada, digite:

julia> using SatelliteToolbox 

Informação

Aqui e a seguir, o comando que você deve digitar é o que está apresentado depois de julia>, como você verá no Julia REPL. Todo o resto é o que você deve ver impresso na tela.

Exemplos

Agora, mostrarei algumas análises que podem ser feitas com as funções que estão disponíveis.

Requisitos necessários

Para manter este post pequeno, assumirei que você possui conhecimento sobre linguagem Julia e um pouco de experiência em assuntos relacionados a satélites e órbitas.

Ano Novo na ISS

Vamos ver como podemos calcular onde os astronautas a bordo da ISS estavam durante o Ano Novo em Greenwich! A primeira coisa que devemos fazer é obter a informação sobre a órbita da ISS. Nesse caso, devemos obter o TLE (Two-Line Element), que é um formato de dados consistindo em duas linhas com 70 carácteres cada e que contém toda a informação relacionada à órbita. O TLE a seguir foi obtido do site Celestrak em 4 de Janeiro de 2019, às 12:25 (horário de Brasília).

ISS (ZARYA)             
1 25544U 98067A 19004.25252738 .00000914 00000-0 21302-4 0 9994
2 25544 51.6417 96.7089 0002460 235.6509 215.6919 15.5373082014978

Esse TLE deve ser carregado em uma variável em Julia. Existem inúmeros métodos para se fazer isso utilizando o SatelliteToolbox.jl. Aqui, utilizaremos uma string especial:

julia> iss_tle = tle"""
       ISS (ZARYA)
       1 25544U 98067A   19004.25252738  .00000914  00000-0  21302-4 0  9994
       2 25544  51.6417  96.7089 0002460 235.6509 215.6919 15.53730820149783
       """[1]
                             TLE
    ==========================================================
                            Name: ISS (ZARYA)
                Satellite number: 25544
        International designator: 98067A
                    Epoch (Year): 19
                     Epoch (Day):   4.25252738
              Epoch (Julian Day): 2458487.75253
              Element set number: 999
                     Inclination:  51.64170000 deg
                            RAAN:  96.70890000 deg
             Argument of perigee: 235.65090000 deg
                    Mean anomaly: 215.69190000 deg
                 Mean motion (n):  15.53730820 revs/day
               Revolution number: 14978

                              B*: 0.000021 1/[er]

                        1   d
                       ---.--- n: 0.000009 rev/day²
                        2  dt

                        1   d²
                       ---.--- n: 0.000000 rev/day³
                        6  dt²
    ==========================================================

Esse código carrega o primeiro TLE especificado dentro da string entre tle"""...""" para a variável iss_tle.

Agora, devemos inicializar um propagador orbital utilizando o TLE carregado. Nesse caso, vamos usar o SGP4:

julia> orbp = init_orbit_propagator(Val{:sgp4}, iss_tle)
OrbitPropagatorSGP4{Float64}(Orbit{Float64,Float64,Float64,Float64,Float64,Float64,Float64}(2.45848775252738e6, 6784.486114511487, 0.000246, 0.9013176963271557, 1.687888720981944, 4.112884090287905, 3.764246850766715), SGP4_GravCte{Float64}
  R0: Float64 6378.137
  XKE: Float64 0.07436685316871385
  J2: Float64 0.00108262998905
  J3: Float64 -2.53215306e-6
  J4: Float64 -1.61098761e-6
, SGP4_Structure{Float64}
  epoch: Float64 2.45848775252738e6
  n_0: Float64 0.06779429624677841
  e_0: Float64 0.000246
  i_0: Float64 0.9013176963271557
  Ω_0: Float64 1.687888720981944
  ω_0: Float64 4.112884090287905
  M_0: Float64 3.764533824882357
  bstar: Float64 2.1302e-5
  Δt: Float64 0.0
  a_k: Float64 1.0637096874073868
  e_k: Float64 0.000246
  i_k: Float64 0.9013176963271557
  Ω_k: Float64 1.687888720981944
  ω_k: Float64 4.112884090287905
  M_k: Float64 3.764533824882357
  n_k: Float64 0.06778673761247853
  all_0: Float64 1.0637096874073868
  nll_0: Float64 0.06778673761247853
  AE: Float64 1.0
  QOMS2T: Float64 1.880276800610929e-9
  β_0: Float64 0.9999999697419996
  ξ: Float64 19.424864323113187
  η: Float64 0.005082954423839129
  sin_i_0: Float64 0.7841453225081564
  θ: Float64 0.6205772419196336
  θ²: Float64 0.3851161131885794
  A_30: Float64 2.53215306e-6
  k_2: Float64 0.000541314994525
  k_4: Float64 6.0412035375e-7
  C1: Float64 4.150340425449004e-10
  C3: Float64 0.0052560300138783985
  C4: Float64 7.530189312128724e-7
  C5: Float64 0.0005696111334271365
  D2: Float64 1.4236674016273006e-17
  D3: Float64 7.305590907411524e-25
  D4: Float64 4.371134237708994e-32
  dotM: Float64 0.06779430410299993
  dotω: Float64 4.494429738092806e-5
  dotΩ1: Float64 -6.037619137612582e-5
  dotΩ: Float64 -6.0409004140717795e-5
  algorithm: Symbol sgp4
  sgp4_gc: SGP4_GravCte{Float64}
  sgp4_ds: SatelliteToolbox.SGP4_DeepSpace{Float64}
)

A variável orbp agora possui a estrutura do propagador orbital do tipo SGP4 com a órbita especificada pelo TLE iss_tle. Esse TLE foi gerado no dia Juliano 2458487.75253 (2019-01-04 06:03:38.592 +0000). Portanto, devemos retro-propagar a órbita para o instante desejado 2019-01-01 00:00:00.000 +0000 (Ano Novo em Greenwich). Isso pode ser feito através da função propagate_to_epoch! como a seguir:

julia> o,r_teme,v_teme = propagate_to_epoch!(orbp, DatetoJD(2019,1,1,0,0,0))
(Orbit{Float64,Float64}(2.4584845e6, 6.784512486615914e6, 0.00024606332568975184, 0.901078684552016, 1.971145471942967, 3.902165730735552, 0.4001307036976068), [4.61152e6, -9.76729e5, -4.88282e6], [-998.41, 7209.56, -2387.48])

A função propagate_to_epoch! retorna três valores. O primeiro o são os elementos orbitais osculantes no instante da propagação, o segundo r_teme é o vetor posição e o último v_teme é o vetor velocidade. Esses vetores são representados no mesmo sistema de referência utilizado para descrever os elementos orbitais quando o propagador foi inicializado. Como estamos utilizando o TLE, então esses vetores estão representados no sistema de referência TEME (True Equator, Mean Equinox).

TEME é um sistema de referência inercial (ECI, Earth-Centered Inertial). Então, devemos converter o vetor posição para um sistema de referência fixo à Terra (ECEF, Earth-Centered, Earth-Fixed) para que possamos calcular a posição da ISS (latitude, longitude e altitude) no instante desejado. O SatelliteToolbox.jl possui toda a teoria IAU-76/FK5 relacionada com a conversão entre sistemas de referência. Neste exemplo, converteremos o TEME para o International Terrestrial Reference Frame (ITRF) para um cálculo mais preciso. Para esse tipo de conversão, deve-se obter os dados da orientação da Terra (EOP, Earth Orientation Data) que são fornecidos pelo IERS. O SatelliteToolbox.jl pode facilmente carregar e usar esses dados através do comando:

julia> eop = get_iers_eop()
[ Info: Downloading file 'EOP_IAU1980.TXT' from 'https://datacenter.iers.org/data/latestVersion/223_EOP_C04_14.62-NOW.IAU1980223.txt'.

Aviso

Os dados EOP são medidos e publicados pelo IERS. Isso significa que o sistema de referência ITRF não pode ser utilizado aqui para prever qual será a posição da ISS, uma vez que os dados EOP não estarão disponíveis. Nesse caso, se um cálculo preciso não é necessário, então pode-se utilizar o sistema de referência Pseudo-Earth Fixed (PEF). A conversão entre o TEME e o PEF utilizando a teoria IAU-76/FK5 não requer dados externos.

A DCM (matriz de cosenos diretores) que gira o TEME para se alinhar ao ITRF é calculada por:

julia> D_ITRF_TEME = rECItoECEF(TEME(), ITRF(), DatetoJD(2019,1,1,0,0,0), eop)
3×3 StaticArrays.SArray{Tuple{3,3},Float64,2,9}:
 -0.17984      0.983696     6.88009e-7
 -0.983696    -0.17984     -1.32013e-6
 -1.17487e-6  -9.14204e-7   1.0

Então, o vetor posição representado no ITRF é:

julia> r_itrf = D_ITRF_TEME*r_teme
3-element StaticArrays.SArray{Tuple{3},Float64,1,3}:
 -1.7901372282879825e6
 -4.360672084341918e6
 -4.882826426845337e6

Finalmente, considerando o elipsóide de referência WGS-84, a latitude, longitude e altitude da ISS durante o Ano Novo em Greenwich pode ser obtida pela função ECEFtoGeodetic conforme mostrado a seguir:

julia> lat,lon,h = ECEFtoGeodetic(r_itrf)
(-0.8061562370049091, -1.9603374908330662, 419859.0733353887)

julia> rad2deg(lat)
-46.18935000852941

julia> rad2deg(lon)
-112.3190646460004

julia> h/1000
419.8590733353887

ou seja, latitude 46.189° S, longitude 112.319° W e altitude 419.859 km. Isso está de acordo com o histórico disponível no website I.S.S. Tracker:

Informação

A pequena diferença pode ser explicada pelo TLE que utilizamos para calcular a posição da ISS e o que é utilizado pelo histórico do I.S.S. Tracker. Uma vez que o TLE possui erros e a propagação orbital os aumenta, é melhor obter um TLE gerado o mais próximo possível do instante desejado. Isso não foi buscado aqui.

Perfil de densidade atmosférica

Neste segundo exemplo, usaremos as funções disponíveis no SatelliteToolbox.jl para calcular o perfil de densidade atmosférica. Existem muitos modelos disponíveis na literatura. O SatelliteToolbox.jl implementa quatro deles: o modelo atmosférico exponencial, o Jacchia-Roberts 1971, o Jacchia-Bowman 2008 e o NRLMSISE-00. Todos eles, exceto o primeiro, requerem como entrada alguns índices espaciais, como o F10.7, que mede a atividade do Sol, e o Ap, que mede a atividade geomagnética. O SatelliteToolbox.jl está preparado para obter todos os arquivos necessários da Internet a fim de que esses índices possam ser facilmente calculados. Isso pode ser feito por:

julia> init_space_indices(wdcfiles_newest_year = 2018)
[ Info: Downloading file 'DTCFILE.TXT' from 'http://sol.spacenvironment.net/jb2008/indices/DTCFILE.TXT'.
[ Info: Downloading file 'fluxtable.txt' from 'ftp://ftp.geolab.nrcan.gc.ca/data/solar_flux/daily_flux_values/fluxtable.txt'.
[ Info: Downloading file 'SOLFSMY.TXT' from 'http://sol.spacenvironment.net/jb2008/indices/SOLFSMY.TXT'.
[ Info: Downloading file 'kp2017.wdc' from 'ftp://ftp.gfz-potsdam.de/pub/home/obs/kp-ap/wdc/kp2017.wdc'.
[ Info: Downloading file 'kp2016.wdc' from 'ftp://ftp.gfz-potsdam.de/pub/home/obs/kp-ap/wdc/kp2016.wdc'.
[ Info: Downloading file 'kp2018.wdc' from 'ftp://ftp.gfz-potsdam.de/pub/home/obs/kp-ap/wdc/kp2018.wdc'.

Informação

A keyword wdcfiles_newest_year não é necessária, mas foi utilizada aqui para evitar um erro dado que o arquivo kp2019.wdc não estava disponível quando esse tutorial foi escrito. Para maiores informações, veja a documentação.

Vamos calcular o perfil de densidade atmosférica de 100 km a 1000 km (passo de 1 km), utilizando todos os quatro modelos, em 2018-11-1 00: 00: 00 +0000 sobre a cidade de São José dos Campos, SP, Brasil (Latitude 23.2237 ° S, longitude 45,9009 ° W).

O modelo atmosférico exponencial, o mais simples, não depende nem dos índices espaciais nem da localização, apenas da altitude. Assim, o perfil atmosférico é calculado por:

julia> at_exp = expatmosphere.(100e3:1e3:1000e3)
901-element Array{Float64,1}:
 5.297e-7
 4.4682006197154693e-7
 3.7690800034030025e-7
 3.1793478585921236e-7
 2.6818886298003355e-7
 2.2622647607479199e-7
 1.9082976790512217e-7
 1.6097143424840973e-7
 1.357849088663833e-7
 1.1453921350666081e-7
 9.661e-8
 8.418342953216043e-8
 7.335524073901479e-8
 6.391983997068123e-8
 5.5698078292918074e-8
 4.8533850005678735e-8
 4.229112868106303e-8
 3.685138444423763e-8
 3.211133345951797e-8
 ⋮
 3.316262334792707e-15
 3.2979959953424034e-15
 3.279830268908571e-15
 3.2617646013035953e-15
 3.2437984413923887e-15
 3.225931241075577e-15
 3.2081624552727763e-15
 3.190491541905968e-15
 3.172917961882957e-15
 3.155441179080928e-15
 3.1380606603300916e-15
 3.1207758753974133e-15
 3.1035862969704427e-15
 3.086491400641224e-15
 3.069490664890299e-15
 3.0525835710707956e-15
 3.0357696033926073e-15
 3.019e-15

Cada elemento é a densidade atmosférica [kg / m³] relacionada a uma altitude.

Informação

Aqui, utilizamos o operador de broadcast . para calcular a densidade em todo o intervalo de altitudes utilizando apenas uma linha de código. Para maiores informações, veja a documentação.

Para o modelo Jacchia-Robert 2008, devemos especificar a latitude geodésica [rad], a longitude [rad] e a altitude [m]. Observe que, como já inicializamos os índices espaciais, todas as informações necessárias serão obtidas automaticamente:

julia> at_jb2008 = jb2008.(DatetoJD(2018,11,1,0,0,0), deg2rad(-23.2237), deg2rad(-45.9009), 100e3:1e3:1000e3)
901-element Array{JB2008_Output{Float64},1}:
 JB2008_Output{Float64}
  nN2: Float64 8.646783160590725e18
  nO2: Float64 2.002012321529147e18
  nO: Float64 6.354156624129256e17
  nAr: Float64 1.0339387366523766e17
  nHe: Float64 1.4269240166433797e14
  nH: Float64 0.9667293309147541
  rho: Float64 5.323059860159937e-7
  T_exo: Float64 677.6952552196058
  Tz: Float64 193.02497135461914

 JB2008_Output{Float64}
  nN2: Float64 7.196278316404481e18
  nO2: Float64 1.6373098215597222e18
  nO: Float64 5.865508118864621e17
  nAr: Float64 8.604946802614022e16
  nHe: Float64 1.1875563628018653e14
  nH: Float64 0.9651359622807867
  rho: Float64 4.430112278422747e-7
  T_exo: Float64 677.6952552196058
  Tz: Float64 195.66319833921864

...

 JB2008_Output{Float64}
  nN2: Float64 7.63695901782906
  nO2: Float64 0.005854548029499879
  nO: Float64 1.9770627164533935e7
  nAr: Float64 4.496252494092111e-9
  nHe: Float64 9.905504449767151e10
  nH: Float64 2.2106945558924915e11
  rho: Float64 1.0288388821434309e-15
  T_exo: Float64 677.6952552196058
  Tz: Float64 677.6523510137463

Cada elemento é uma instância da estrutura JB2008_Output que contém a densidade dos constituintes da atmosfera em [kg / m³] relacionadas a uma altitude.

O modelo NRLMSISE-00 requer as mesmas informações, mas em uma ordem diferente. Mais uma vez, como já inicializamos os índices espaciais, todas as informações necessárias são obtidas automaticamente:

julia> at_nrlmsise00 = nrlmsise00.(DatetoJD(2018,11,1,0,0,0), 100e3:1e3:1000e3, deg2rad(-23.2237), deg2rad(-45.9009))
901-element Array{NRLMSISE00_Output{Float64},1}:
 NRLMSISE00_Output{Float64}
  den_N: Float64 3.225647667164233e11
  den_N2: Float64 1.1558415665482785e19
  den_O: Float64 4.649965500403523e17
  den_aO: Float64 4.631659520454273e-37
  den_O2: Float64 2.6263326718789934e18
  den_H: Float64 2.533162671436194e13
  den_He: Float64 1.2320073447340945e14
  den_Ar: Float64 1.160809744681819e17
  den_Total: Float64 6.968049043353933e-7
  T_exo: Float64 1027.3184649
  T_alt: Float64 215.25904311781903
  flags: NRLMSISE00_Flags

 NRLMSISE00_Output{Float64}
  den_N: Float64 3.5691748620576013e11
  den_N2: Float64 1.0012594360559639e19
  den_O: Float64 4.6065649467733504e17
  den_aO: Float64 1.8947221849916303e-36
  den_O2: Float64 2.2358683099960123e18
  den_H: Float64 2.35052621777078e13
  den_He: Float64 1.121318459050076e14
  den_Ar: Float64 9.845632305957742e16
  den_Total: Float64 6.029280387245405e-7
  T_exo: Float64 1027.3184649
  T_alt: Float64 213.7198609515809
  flags: NRLMSISE00_Flags

...

 NRLMSISE00_Output{Float64}
  den_N: Float64 5.097414213511122e6
  den_N2: Float64 106.44260908421272
  den_O: Float64 7.592117961201309e7
  den_aO: Float64 2.0421422183370042e9
  den_O2: Float64 0.07891050457188623
  den_H: Float64 1.3014187084557657e11
  den_He: Float64 8.255445331597499e10
  den_Ar: Float64 1.442732462296107e-7
  den_Total: Float64 8.205713083292234e-16
  T_exo: Float64 724.4998315669409
  T_alt: Float64 724.4998315398782
  flags: NRLMSISE00_Flags

Cada elemento é uma instância da estrutura NRLMSISE00_Output que contém a densidade dos constituintes da atmosfera em [kg / m³] relacionadas a uma altitude.

NRLMSISE-00 online

Os valores calculados aqui podem ser comparados à versão online do modelo NRLMSISE-00. Observe, no entanto, que a versão online no momento em que este tutorial foi escrito permite apenas o cálculo entre 1960/02/14 e 2018/03/17.

O modelo Jacchia-Roberts 1971 não possui suporte à obtenção automática dos índices espaciais ainda. Portanto, precisaremos fazer isso manualmente. São necessários três índices: o F10.7 diário, a média do F10.7 (janela de 81 dias, centrada no tempo de entrada) e o índice geomagnético Kp (com um atraso de 3 horas). Essa informação pode ser obtida por:

julia> F107 = get_space_index(F10(), DatetoJD(2018,11,1,0,0,0))
65.8

julia> F107m = get_space_index(F10M(), DatetoJD(2018,11,1,0,0,0); window = 81)
68.29135802469136

julia> kp = get_space_index(Kp(), DatetoJD(2018,11,1,0,0,0)-3/24)
0.875

Assim, o perfil atmosférico computado por JR1971 é obtido por:

julia> at_jr1971 = jr1971.(DatetoJD(2018,11,1,0,0,0), deg2rad(-23.2237), deg2rad(-45.9009), 100e3:1e3:1000e3, F107, F107m, kp)
901-element Array{JR1971_Output{Float64},1}:
 JR1971_Output{Float64}
  nN2: Float64 5.7192521805880885e19
  nO2: Float64 1.0370130013611293e19
  nO: Float64 1.2248930040184422e19
  nAr: Float64 4.797323831280962e17
  nHe: Float64 3.150115795435398e15
  nH: Float64 0.0
  rho: Float64 3.4060767884871413e-6
  T_exo: Float64 657.1677377132266
  Tz: Float64 191.2125970557249

 JR1971_Output{Float64}
  nN2: Float64 8.080629730659746e18
  nO2: Float64 1.6344357654348047e18
  nO: Float64 1.0616020508369499e18
  nAr: Float64 9.003667794720594e16
  nHe: Float64 7.366460213846506e13
  nH: Float64 0.0
  rho: Float64 4.96920849944813e-7
  T_exo: Float64 657.1677377132266
  Tz: Float64 193.37386958205954

 ...

 JR1971_Output{Float64}
  nN2: Float64 4.812412406159134
  nO2: Float64 0.00284582346906836
  nO: Float64 2.7544317744479913e7
  nAr: Float64 1.3825250583995572e-9
  nHe: Float64 1.0969262471305798e11
  nH: Float64 3.6326221208226086e11
  rho: Float64 1.3378408900651963e-15
  T_exo: Float64 669.4417264661663
  Tz: Float64 669.441725137747

Cada elemento é uma instância da estrutura JR1971_Output que contém a densidade dos constituintes da atmosfera em [kg / m³] relacionadas a uma altitude.

Finalmente, utilizando o pacote PyPlot.jl, os perfis atmosféricos (altitude vs. densidade) em escala semi-log podem ser plotados usando:

julia> using PyPlot
julia> figure()
julia> h = 100:1:1000
julia> semilogx(at_exp, h, map(x->x.rho, at_jb2008), h, map(x->x.den_Total,at_nrlmsise00), h, map(x->x.rho,at_jr1971), h)
julia> legend(["Exp.", "JB2008", "NRLMSISE-00", "JR1971"])
julia> xlabel("Density [kg/m^3]")
julia> ylabel("Altitude [km]")
julia> title("Atmospheric Density, 2018-11-01 00:00:00+0000")
julia> grid()

que leva a:

Atmospheric density: 2018-11-01 00:00:00+0000

Para mais informações sobre as muitas opções para calcular a densidade atmosférica, consulte a documentação.

Atenção

Se a data desejada não estiver disponível nos arquivos dos índices espaciais, um erro será gerado ao calcular a densidade atmosférica.

Conclusão

Espero que este tutorial tenha ajudado você a entender um pouco de como o pacote SatelliteToolbox.jl pode ser usado para realizar análises relacionadas a satélites e órbitas. Se você tiver alguma dúvida, por favor, fique à vontade para deixar um comentário abaixo!

Atualização 2019-03-24: Atualização do tutorial para considerar a versão v0.6.0 do SatelliteToolbox.jl.

Atualização 2018-01-08: Correção do algoritmo que plota os dados, pois a variável h não estava corretamente definida (obrigado Bernard_GODARD).

Mudanças na v3 do OrdinaryDiffEq.jl

Olá!

Eu sei, um longo, longo tempo passou desde a minha última postagem. Eu gostaria de me desculpar! Estive muito ocupado devido a uma quantidade enorme de trabalho a ser feito no INPE.

De qualquer forma, esse será uma postagem pequena. Apenas gostaria de atualizar o código-fonte da última lição de tal forma que rode corretamente utilizando Julia v0.6 e OrdinaryDiffEq.jl v3.

Como eu disse, o pacote OrdinaryDiffEq.jl, que é parte do projeto JuliaDiffEq, é maravilhoso. Na minha humilde opinião, é um dos melhores solvers de equações diferenciais de código aberto disponíveis. O desenvolvimento dessa suíte está ocorrendo a uma taxa muito elevada. Isso é bom, pois o pacote tem ficado melhor a cada dia. O lado ruim é que a API também se altera muito rápido e às vezes seu código precisa ser modificado para continuar a funcionar nas novas versões.

Experiência necessária

Isso é apenas uma atualização do código-fonte apresentado na minha última postagem. Portanto, recomendo fortemente a lê-la primeiro:

Mudanças na v3 do OrdinaryDiffEq.jl

Desde que eu postei meu último tutorial, a API do OrdinaryDiffEq.jl mudou (na verdade, a API é definida em DiffEqBase.jl, mas ele será atualizado quando você atualizar o OrdinaryDiffEq.jl). Existem três modificações que devem ser feitas para que o código-fonte anterior possa funcionar corretamente com Julia v0.6 e OrdinaryDiffEq.jl v3.0.3.

Interface de Array de Dados Customizada

Anteriormente, uma interface de array de dados customizada podia ser criada através de uma nova estrutura que herdava o tipo abstrato DEDataArray{T}. Isso se alterou e agora o novo tipo abstrato correto é DEDataVector{T}. Então, nossa definição do SimWorkspace deve ser:

type SimWorkspace{T} <: DEDataVector{T}
    x::Array{T,1}
    a::T
end

Assinatura da equação dinâmica

Na versão utilizada no último tutorial, a assinatura da função dinâmica era:

function f(t,u,du)

Devido a algumas mudanças internas, a assinatura agora é:

function f(du,u,t,p)

p é um vetor de parâmetros que não é utilizado nesse exemplo, mas simplifica muitas coisas. Se você quiser saber mais informações, por favor veja http://docs.juliadiffeq.org/stable/tutorials/ode_example.html#Defining-Parameterized-Functions-1.

Assinatura da função de condição do callback

A assinatura da função de condição do callback se alterou de:

function condition_control_loop(t,u,integrator)

para:

function condition_control_loop(u,t,integrator)

Código-fonte atualizado

Finalmente, o código-fonte atualizado do último tutorial pode ser visto a seguir. Ele irá fornecer exatamente o mesmo resultado.

using DiffEqBase
using OrdinaryDiffEq

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

# Global variable to store the acceleration.
# This solution was selected just for the sake of simplification. Don't use
# global variables!
a = 0.0

# Parameters.
r = 10.0   # Reference position.
k = 0.3    # Proportional gain.
d = 0.8    # Derivative gain.

# Configuration of the simulation.
tf     = 30.0             # Final simulation time.
tstops = collect(0:1:tf)  # Instants that the control loop will be computed.
u0     = [0.0; 0.0]       # Initial state.

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

# Dynamic equation.
function dyn_eq(du,u,p,t)
    du .= [0 1; 0 0]*u + [0;1]*a
end

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

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

# Affect function.
function control_loop!(integrator)
    global a

    p = integrator.u[1]
    v = integrator.u[2]

    a = k*(r-p) + d*(0.0-v)
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)
nothing

Conclusões

Esse tutorial apenas atualiza o código-fonte do último de tal forma que se possa executá-lo utilizando Julia v0.6 e OrdinaryDiffEq.jl v3.0.3.

Estou planejando para o meu próxima tutorial um exemplo simples de como codificar filtros de Kalman em linguagem Julia!

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:

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

Old way to create a custom type in OrdinaryDiffEq
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 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`:

type MyDataArray{T} <: DEDataArray{T}
    x::Array{T,1}
    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:

type SimWorkspace{T} <: DEDataArray{T}
    x::Array{T,1}
    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(t,u,du)
    du[1] = u[2]
    du[2] = u.a
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 user_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 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`:

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 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.

Source code

using DiffEqBase
using OrdinaryDiffEq

###############################################################################
#                                    Types
###############################################################################

type SimWorkspace{T} <: DEDataArray{T}
    x::Array{T,1}
    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(t,u,du)
    du[1] = u[2]
    du[2] = u.a
end

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

# Condition function.
function condition_control_loop(t,u,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 user_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)

Utilizando julia para simular sistemas compostos de partes contínuas e discretas

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. Note que a teoria permanece a mesma!

Olá!

Nesse primeiro post, eu gostaria de explicar como eu contornei um problema que eu sempre tive relacionado com a simulação de sistemas dinâmicos que possuem partes discretas e partes contínuas. Note que com certeza esse é o território do MathWorks® Simulink. Não existe dúvida que o Simulink é a ferramenta mais rápida de prototipagem de tais tipos de aplicação em que você possui um modelo dinâmico contínuo modificado por partes discretas que são amostradas em intervalos de tempo específicos.

Entretanto, existem algumas situações em que você não consegue utilizar o Simulink (não ter uma licença é a mais óbvia…). No meu caso de uso específico, eu estou codificando o simulador do Subsistema Controle de Atitude e Órbita (AOCS, em inglês) para o satélite Amazonia-1. O Simulink me deu alguns problemas relacionados com a simulação de algumas descontinuidades, como, por exemplo, sombras de algumas partes do satélite nos sensores solares. Verifiquei que eu necessitava de mais controle no solver do que o Simulink poderia me oferecer. Portanto, decidi transcrever o simulador para uma linguagem de programação “padrão”. Após um extensiva pesquisa, a linguagem julia e o pacote OrdinaryDiffEq me forneceram a melhor solução para simular o meu sistema. Portanto, neste post, gostaria de compartilhar como eu fui capaz de resolver o problema de simular sistemas com partes discretas e contínuas utilizando essas ferramentas.

Experiência necessária

Tentarei manter esse tutorial o mais simples possível. No entanto, ajudaria se você possuir um conhecimento mínimo sobre:

  1. Linguagem julia;
  2. Sistemas de controle; e
  3. Equações diferenciais ordinárias.

Definição do problema

Na maioria dos sistemas de controle, você possui quatro partes:

  1. A planta, que é o que você quer controlar;
  2. Os sensores, que medem quantidades relacionadas com os processos que você quer controlar;
  3. Os atuadores, que agirão e modificarão o estado da sua planta; e
  4. O controlador, que irá ler a saída dos sensores e comandará os atuadores para modificar a planta dado um conjunto de especificações.

No entanto, hoje em dia o controlador é normalmente embutido em um microcomputador, que possui intervalos pré-definidos para executar suas funções. Em outras palavras, apenas em intervalos de tempo específicos (chamados de instantes de amostragem) o controlador irá ler a saídas dos sensores, atualizar a lei de controle, e enviar novos comandos para os atuadores.

Para simular tais sistemas, você deverá considerar dois tipos de modelos: contínuos e discretos. A planta, por exemplo, é um modelo contínuo, uma vez que é normalmente modelada por um conjunto de equações diferenciais e geralmente varia para todo `t`. O controlador, conforme mencionado anteriormente, apenas atualiza seu estado interno em um conjunto de instantes determinados.

Academicamente falando, a simulação de um sistema desse tipo pode ser resolvida utilizando a transformada Z para converter o modelo contínuo da planta em um modelo discreto. Dessa forma, todos os modelos serão discretos e poderão ser facilmente simulados. Esse método possui muitos problemas. O maior deles, na minha opinião, é que ele funciona apenas para sistemas lineares e pode ser muito difícil de ser aplicado em plantas complexas com muitas variáveis de estado.

Portanto, o problema é: como eu posso utilizar uma linguagem de programação “padrão” para simular um sistema em que a parte contínua (a planta) é controlada por um processo discreto (o controlador)? Eu irei mostrar a seguir uma solução possível utilizando julia e o pacote OrdinaryDiffEq.

Sistema dinâmico proposto

Eu decidi propor um sistema dinâmico simples para ser controlado. Isso ajudará a tornar o tutorial mais interessante e fácil de entender. O problema que selecionei é um simples controle de uma partícula em uma dimensão definido como se segue:

Controlar a posição de uma partícula que se move em uma dimensão. Você poderá alterar a aceleração da partícula e você possuí acesso às medidas de posição e velocidade.

Essa partícula pode ser modelada por um conjunto de duas equações diferenciais descrevendo a derivada temporal de sua posição e velocidade:

\(\begin{array}{ccc} \dot{p} &=& v \\ \dot{v} &=& a \end{array} ,\)

onde \(p \) é a posição da partícula, \(v \) é a velocidade da partícula e \(a \) é a aceleração da partícula, que podemos controlar.

Então, definimos o vetor de estado como:

\(\mathbf{x} = \left[\begin{array}{c} p \\ v \end{array}\right] . \)

Finalmente, nosso sistema dinâmico pode ser escrito como:

\(\dot{\mathbf{x}} = \left[\begin{array}{cc} 0 & 1 \\ 0 & 0 \end{array}\right]\mathbf{x} + \left[\begin{array}{c} 0 \\ 1 \end{array}\right] \cdot a . \)

Nota: Esse é precisamente um daqueles sistemas que podem ser facilmente simulado utilizando a abordagem da transformada Z quando você possui um controlador discreto. Mas a ideia é manter esse tutorial simples de tal forma que você possa entender o conteito para estendê-lo a cenários mais complexos.

Controlador

Nós queremos controlar a aceleração de uma partícula de tal forma que ela vá para um ponto específico no espaço, chamado \(r \). Nesse caso, eu irei propor um controlador PD (muito) simples. A lei de controle será:

\(a = 0.3 \cdot (r-p) + 0.8 \cdot (0.0-v) \)

O controle deverá ser executado em um processo discreto que é atualizado a cada 1s.

Nota: Para aqueles que não estão familiarizados com sistemas de controle, apenas assuma que calculando a aceleração conforme mostrado aqui, nós podemos levar a partícula para a posição desejada \(r \).

Nas próximas seções, eu descreverei como você poderá usar julia e o pacote OrdinaryDiffEq para simular esse sistema em uma solução muito elegante. Você verá que será capaz de facilmente estender essa ideia para lidar com sistemas muito mais complexos.

Por que eu escolhi julia?!

No meu instituto, o INPE, eu tenho que responder quase todo dia o porquê eu escolhi julia como a linguagem de programação para realizar meu trabalho. É um pouco difícil introduzir novas tecnologias especialmente quando você está lidando com áreas que pessoas ainda utilizam FORTRAN devido a todos aqueles códigos legados.

Eu comecei a utilizar julia desde a versão 0.2 (estamos atualmente na versão 0.5). As razões do porquê eu realmente gostei dessa linguagem são:

  1. Ela parece uma linguagem interpretada no sentido que você não precisa se preocupar com os tipos das variáveis, o que deixa um programador MATLAB bastante confortável;
  2. Você pode muito facilmente chamar funções escritas em C / FORTRAN, isso é importante porque modelos como IGRF e MSIS são disponibilizados em FORTRAN;
  3. O desempenho das simulações quando a vetorização é impossível (ou muito difícil) é muito melhor que o de outras linguagens interpretadas (como MATLAB, Octave, Python, etc.);
  4. Você não precisa se preocupar em vetorizar o seu código;
  5. Programação paralela é muito fácil de se obter; e
  6. A sintaxe é muito elegante (ok, eu sei, isso é uma opinião pessoal 🙂 ).

Claro que existem muitos outros pontos positivos sobre julia, por favor visite o website http://www.julialang.org para mais informações.

Além disso, para o tipo de aplicação que eu estou focando nesse tutorial, o pacote OrdinaryDiffEq desenvolvido e mantido por Chris Rackauckas é único. Chris desenvolvou uma interface maravilhosa para se interagir com os solvers que deixa você fazer praticamente tudo. Conclusão: você possui todos aqueles solvers que você precisa (Runge-Kutta, Dormand-Prince, etc.) podendo interagir muito profundamente com o algoritmo para adaptá-lo às suas necessidades específicas. Devo admitir que eu não conheço nenhum outro pacote que forneça uma solução tão elegante para nenhuma outra linguagem (por favor, me avise nos comentários se eu estou me esquecendo de algo!).

Configuração

Primeiramente, faça o download do pacote julia apropriado para o seu sistema operacional em http://www.julialang.org (note que eu testei tudo aqui utilizando o openSUSE Tumbleweed e julia 0.5).

Após a instalação, abra julia e instale o OrdinaryDiffEq:

Pkg.init()
Pkg.add("OrdinaryDiffEq")

Para mais informações sobre esses comandos para instalar o pacote, por favor leia a documentação aqui: http://docs.julialang.org/en/stable/manual/packages/

Simulando o sistema

A primeira coisa é carregar os módulos que vamos precisar:

Pkg.init()
Pkg.add("OrdinaryDiffEq")

Resolver uma EDO utilizando o OrdinaryDiffEq é muito similar aos outros solver disponíveis no MATLAB, por exemplo `ode45`. Você deverá definir uma função dinâmica, selecionar o estado inicial, ajustar algumas opções, e chamar o solver. Então, vamos definir a função dinâmica.

No OrdinaryDiffEq, a função dinâmica deverá possuir a seguinte assinatura:

function f(t,u,du) 

onde `t` é o tempo, `u` é o vetor de estado e `du` é a derivada temporal do vetor de estado. Note que `u` e `du` possuem a mesma dimensão do estado inicial. Essa função deverá computar `du`, mas o usuário deve ser bastante cuidadoso para não alterar a referência ao vetor! Caso contrário, o solver não irá utilizar o valor computado.

Para o nosso problema proposto, a função dinâmica será:

function dyn_eq(t,u,du)
    du .= [0 1;
           0 0]*u + [0;1]*a
end

Nota 1: Nesse exemplo, a variável de aceleração será global. Isso tornará as coisas simples, mas é um padrão de codificação muito (MUITO) ruim. Entretanto, a maneira correta de se evitar essa variável global é um pouco mais complicada e será tratada em um tutorial futuro.

Nota 2: O sinal `.=` significa que os valores no lado esquerdo serão copiados para o vetor `du` sem a alteração da referência.

Ok, nosso modelo do sistema dinâmico (a parte contínua do nosso problema) está pronto! A seguir, nós iremos codificar a lei de controle discreta.

Controlador discreto utilizando callbacks

Callbacks são um tipo de interface com o solver que faz o pacote OrdinaryDiffEq ser bastante poderoso. Você define duas funções: uma para verificar se um evento ocorreu (função condição) e outra para ser executada caso esse evento tenha ocorrido (função efeito). Dessa forma, a cada instante de amostragem, o solver chama a função condição e, se ele disser que o evento ocorreu, então a função efeito é executada.

O OrdinaryDiffEq provê dois tipos de Callbacks: contínuos e discretos. Nesse caso, estamos interessados em callbacks discretos em que você consegue definir precisamente antes da execução quais serão os instantes de amostragem em que a função efeito será chamada. Para maiores informações, por favor veja a documentação em http://docs.juliadiffeq.org/latest/features/callback_functions.html

Vamos definir primeiro a função condição para verificar se nós estamos em um intervalo de amostragem em que a lei de controle deve ser atualizada:

tf = 30.0
tstops = collect(0:1:tf)
function condition_control_loop(t,u,integrator)
    (t in tstops)
end

A variável `tf` define o tempo de simulação e o array `tstops` define os intervalos de amostragem do loop de controle. A função `condition_control_loop`, nossa função condição, simplesmente verifica se `t` é um dos nossos intervalos de amostragem. Se essa função retornar `true`, então a função efeito será executada. Note que também é possível definir uma espécie de callback baseado no vetor de estado `u`, mas isto não é o nosso caso.

Não se preocupe com a variável `integrator`. Ela basicamente define uma interface com o integrador que deixa você verificar e alterar muitas e muitas coisas. Eu irei fornecer mais informações em tutoriais futuros. Se você desejar saber mais sobre isso, verifique a documentação em http://docs.juliadiffeq.org/latest/basics/integrator.html

Agora, vamos definir nossa função efeito que, nesse caso, é nosso controlador:

r = 10.0
k = 0.3
d = 0.8
function control_loop!(integrator)
    global a

    p = integrator.u[1]
    v = integrator.u[2]

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

onde `r` é a posição de referência e `k` e `d` são, respectivamente, os termos proporcionais e derivativos do controlador PD. Nesse caso, o loop de controle apenas necessita calcular a nova aceleração, que está escrita na variável global 😕 `a`, utilizando a posição e velocidade atual da partícula.

Agora, nós definiremos o callback:

cb = DiscreteCallback(condition_control_loop, control_loop!)

Executando o solver

Finalmente, precisamos apenas definir nosso estado inicial e chamar o solver do pacote OrdinaryDiffEq para solucionar o problema para nós:

u0   = [0.0; 0.0]
prob = ODEProblem(dyn_eq, u0, (0.0, tf))
sol  = solve(prob, Tsit5(), callback = cb, tstops=tstops)

`Tsit5()`significa que estamos utilizando o método Tsitouras 5/4 Runge-Kutta. Outras soluções possíveis são `DP5()` para o método Dormand-Prince’s 5/4 Runge-Kutta, ou `BS3()` para o método Bogacki-Shampine 3/2. Todos esses métodos possuem passo de amostragem adaptativo, que é uma boa escolha para simular sistemas dinâmicos, mas você também poderá escolher métodos com passo fixo. Para maiores informações, veja a documentação em http://docs.juliadiffeq.org/latest/solvers/ode_solve.html

Após isto, a solução do problema será armazenada na variável `sol`. O vetor de tempo contendo todos os passos de amostragem aceitos pelo solver é `sol.t` e os vetores de estado em cada um desses instantes são armazenados em `sol.u`.

Se você estiver utilizando o PyPlot, você poderá plotar os resultados utilizando:

plot(sol.t, sol.u)

O resultado esperado dessa simulação é:

Note que pode-se ver facilmente os pontos em que o loop de controle é calculado, onde imediatamente se altera a velocidade da partícula.

Conclusão

Nesse primeiro tutorial, eu mostrei como eu utilizei julia para simular sistemas que possuem partes contínuas e discretas. O conceito foi introduzido por um problema muito simples em que a posição de uma partícula deve ser controlada utilizando um loop de controle discreto.

Eu espero que você possa ter entendido o conceito e o quão fácil será estendê-lo para simulações muito mais complexas. 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.

Source code

using DiffEqBase
using OrdinaryDiffEq

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

# Global variable to store the acceleration.
# This solution was selected just for the sake of simplification. Don't use
# global variables!
a = 0.0

# Parameters.
r = 10.0   # Reference position.
k = 0.3    # Proportional gain.
d = 0.8    # Derivative gain.

# Configuration of the simulation.
tf     = 30.0             # Final simulation time.
tstops = collect(0:1:tf)  # Instants that the control loop will be computed.
u0     = [0.0; 0.0]       # Initial state.

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

# Dynamic equation.
function dyn_eq(t,u,du)
    du .= [0 1; 0 0]*u + [0;1]*a
end

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

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

# Affect function.
function control_loop!(integrator)
    global a

    p = integrator.u[1]
    v = integrator.u[2]

    a = k*(r-p) + d*(0.0-v)
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)
nothing

Introduzindo… meu novo website!

Olá!

Eu finalmente decidi mudar o meu website do servidor do INPE (http://www.inpe.br/ete/dse/ronan) para o meu próprio servidor. Isso ocorreu porque eu precisava de mais liberdade para instalar um CMS, que não é possível atualmente no meu antigo servidor.

Hoje, eu finalmente terminei a configuração e agora vocês podem esperar alguns posts nas próximas semanas sobre o meu trabalho no INPE (o Instituto Nacional de Pesquisas Espaciais). Além disso, vou tentar disponibilizar tutoriais sobre modelagem e simulação de sistemas espaciais, que atualmente é uma das mais minhas maiores responsabilidades, bem como informações relacionadas a processamento de sinais, que foi a área do meu Doutorado.

A linguagem padrão de cada post será o Inglês, mas prometo me esforçar para traduzir cada um deles também para o Português.

Espero que gostem! Fiquem ligados 😉