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