# Category Archives: Simulation

Posts related to simulation of systems.

# My Julia workflow: Readability vs. Performance

Hi!

After my last post, I got many questions like: “What’s Julia? Why you are using Julia? Is it worth it?”

Well, I will try to answer the last two questions by showing my workflow: how I code and what kind of gains I get (for the first question I will point you to the official website 🙂 )

Required background

It is required a little knowledge of C and Julia so that you can fully understand this post.

## My use case

Many of my activities at INPE require analysis that can only be accomplished by optimization algorithms or Monte Carlo simulations. For example, during the Pre-Phase A studies at the Space Mission Integrated Design Center (CPRIME), we sometimes use evolutionary algorithms to obtain optimal solutions for the spacecraft design. It is also common to use Monte Carlo simulations to find the worst case scenarios, e.g. the initial state of the spacecraft (attitude and angular velocity) that yields the longest time to acquire the Sun pointing. Thus, the algorithms that I use are often required to run a thousand, even a million, times. So performance is something that really matters for me.

On the other hand, we are mostly engineers that are used to MATLAB. The algorithms are written to solve equations and simulate dynamic systems. Hence, the closest the algorithm is to the mathematical equations the better.

In the following, I will describe my workflow: how do I write Julia code from the initial validation up to the version that is meant to be used inside an optimizer. I tried many approaches over the years, but this one was what gave me the shortest time to obtain the results.

## Numerical example

To illustrate, let’s assume that we want an algorithm to perform the Kalman filter covariance update of a dynamic system. Do not give too much importance to those equations, I just selected a random scenario to improve the clarity of what we are doing.

Initialization

\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}

Simulation loop

For $$k = 2$$ to $$60000$$, do:

\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}

Every time $$\mathbf{P}_u$$ is modified, store $$\mbox{trace}(\mathbf{P}_u)$$.

Note

$$\mathbf{w}^{\times}$$, $$\mathbf{s}^{\times}$$, and $$\mathbf{B}^{\times}$$ are constant matrices, $$\mathbf{0}_{n \times m}$$ is a $$n \times m$$ matrix filled with 0s, and $$\mathbf{I}_n$$ is an $$n \times n$$ identify matrix.

## First: Code fast, as readable as possible

The very first thing is to code an initial version of the algorithm you want without worrying too much about the performance. In this phase, readability is paramount, because it will make debugging easier.

Our first version of the code looks like:

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

Is it possible to see how close to the mathematical equation it is? Even the symbols like λ can be used in Julia, that supports any UTF-8 character in the variables name. Now, imagine something with hundreds of equations. It will be way easier to validate the code and find problems when you write using a computational language that is close to the mathematical one.

Information

The I in Julia is not a 3 x 3 identity matrix as in the equations above. It is called the UniformScalling operator. In our code, it turns into an identity matrix with the size it must have to make sense, i.e 3 x 3.

In this very first version of the code, we have two tasks: 1) see if it works properly; and 2) see if the performance is good enough for our needs. We know that after 60,000 iterations the trace of the covariance matrix $$P_u$$ is approximately 7.00089. In Julia, we can run the code in REPL and see if this is true:

julia> result = kalman_filter();

julia> result[end]
7.000894875771719

Good 🎊🎉! At least the code does what it is supposed to. The next thing is to measure the performance and verify if it is enough. This is very important: Julia is a very fast language for me not because it is faster than C, but because the time I spend between the initial code and the final result is way shorter. Sometimes you just do not need a performant code and want to code as fast as possible, Julia can do this as I just show. Sometimes you need speed, Julia can do this as well as I will show.

The correct way to benchmark a Julia code is to put inside a function (we have already done that) and use the excellent package BenchmarkTools.jl.

julia> using BenchmarkTools

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

Now we know that the mean execution time of the function (after the compilation step) is 1.296s and it has allocated 2.83 GiB of memory. At this moment I ask myself: “Is this speed enough for what I need to do?” For example, if I need to run this code just once, then I can live with its current state. However, if I need this code to be called inside an optimizer that will run 1,000,000 times, then I would really want to optimize it.

## Second: Modify your code to improve the performance if necessary

There are some things you need to take care to increase the performance of Julia code. The documentation is very good and I truly suggest you read. For starters, you must look for three things:

• Code type-stability: Julia must know the type of the variables given only the type of the inputs. This can be checked with the macro @code_warntype.
• Allocations: the process to allocate is slow, you want to have the least amount possible of allocations to improve the performance.
• Avoid duplicate computations: the compilers are getting smarter every day and sometimes they got it, but it does not harm to help them!

We can easily see that our code is type stable by calling @code_warntype kalman_filter(). You will see that no variable can have multiple types. I will not go deeper here, because this is a subject for an entire post. For more information, see the documentation.

There are some methods to check which parts of the code are allocating. The “official” one is to start Julia with the option --track-allocation=user, execute the function, and then analyze the .mem file that will be created. This is described here. However, there is a method that I personally find better. It is based on the package TimerOutputs.jl. It provides a way to analyze just specific lines. In our case, we should be concerned about the code executed inside the for-loop. The previous file should be modified as follows:

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

Now, when you run kalman_filter();, you will see on the screen:

 ──────────────────────────────────────────────────────────────────
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, now we have something to work with. We will perform four modifications:

• Those three variables Pp, K, and Pu are being allocated at every iteration. One thing we should do is pre-allocate them before the loop and tell Julia to use the same space to store the new values. The latter can be accomplished by using the operator .= that performs a broadcasted assignment. In other words, each element of the left-side variable will be equal to the corresponding element of the right-side expression, leading to no allocation in this operation. It is also good to replace the + for .+ for the same reason, since summing matrices is an operation performed element-by-element.
• The UniformScaling variable will be exchanged by a pre-allocated matrices. This also helped to reduce the allocations.
• Notice that (I - K*Hk) is being computed twice inside the loop. We will use an auxiliary variable to avoid this. It is important that this auxiliary variable is also pre-allocated before the loop for the same reason we discussed on the first point.
• Since we can ensure that all arrays are accessed by valid indices inside the for, then we can use the @inbounds macro as described in the documentation. This avoid bound checking operations, improving performance.

Given that, our code now should look like:

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

The benchmark using this new version provides:

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

We just decreased the computational time 14% with only minor modifications! If this function is used in an optimizer that calls it 1,000,000 times, then we just save more than 2 days.

## Third: Going extreme!

If it is really necessary to reduce the computational time, then we can still do things here. However, spoiler alert, this will hurt the readability. Hence, I advise going this path only after your code has been tested and validated and if you really need more speed.

Looking at the last benchmark, we still see that the code is allocating 2GiB. This is coming from the temporary arrays Julia must allocate to store in-between multiplications. For example, in the multiplication Fk_1*Pu*Fk_1' .+ Q, it needs to allocate at least one temporary array to store Fk_1*Pu or Pu*Fk_1'. The compiler cannot (yet?) guess that this can be allocated before the loop instead of allocating it at every iteration. Hence, if we really need to reduce the computation time, then we need to drastically reduce those allocations.

The approach is quite similar for the second phase: create temporary matrices outside the loop and use them to store the in-between operations. However, here we cannot use the broadcast multiplication operator .*, since it will perform an element-wise multiplication instead of a matrix product. Thus, the solution is to use the following function:

mul!(Y, A, B)

that computes the matrix product A*B and stores in the matrix Y. This function also returns the matrix Y to make operation concatenation easier. Thus, the modified code is:

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

Well… this is… ugly 😅 However, let’s benchmark it:

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

Wow! It’s 23% faster than the previous version and 34% faster than the first version. Compared to the latter, this code would save us more than 5 days if we perform 1,000,000 evaluations. The trade-off here was to have a less readable code, mathematically speaking.

Information

Julia has superb meta-programming capabilities. One can easily construct a macro that will expand the multiplications to use the auxiliary variables while maintaining the readability. However, this is a somewhat advanced topic to cover for now.

## OK Ronan, what about C++?

This is the most asked question when I talk about Julia. Of course, C++ can be faster in execution time. However, since Julia is a dynamic language (seems interpreted), it is way easier to code and test the results. That’s one of the reasons why languages like MATLAB and Python were created. Thus, if you account for the total time (developing, testing, and executing), in my experience, you can beat C++ using Julia (unless, of course, you are a C++ magician, which most engineers I know, including me, are not).

Anyway, let’s just write the very same algorithm in C++ to see how far we are. Notice that these are my C++ skills. I am almost sure that this code can be optimized (can anyone do this for me, please?). However, I am not a very experienced Julia programmer either. Thus, I have a feeling that the Julia code can also be optimized (can anyone also do this for me, please?).

The C++ version uses Eigen to handle matrices operations. Otherwise, it will be a very, very long code. The code can be seen as follows:

#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;
}

In my machine, the mean execution time using -O3 is:

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

which is very, very close to the fastest Julia version of this code.

Information

Some may argue that it would be faster if I used -march=native. Indeed, it is! However, for a fair comparison, I would have also to recompile Julia system image with the same flag, which, unfortunately, I do not have time now. This is an easy procedure if you use the package PackageCompiler.jl.

## Conclusion

I hope that you understood how Julia is helping me to achieve my objectives faster. We now have a language that is easy to code and test algorithms, and that with minor adaptations can reach C performance. I have no more need to translate all the algorithms I need to be executed really fast to a compiled language, which turns out to be very time consuming for an engineer that does not code very well using C/C++.

# The Satellite Toolbox for Julia

Hi!

In this post, I would like to introduce the SatelliteToolbox.jl, which is a package for Julia language with many options to analyze space missions. It is used on a daily basis at the Brazilian National Institute for Space Research (INPE). First, it is presented a brief history about the package, and then I show some interesting analysis that can be done with it.

## History

In 2013, I joined INPE as a Junior Space Systems Engineer. I was assigned to the space systems division in which I had to work with the attitude and orbit control subsystem (AOCS). Since I had only an intermediate knowledge about orbits, I decided to go into this subject by coding algorithms and comparing the results to the INPE heritage and to the literature in my spare time.

The very first thing was to select the language! In my Ph.D., I used MATLAB to simulate an inertial navigation system, but the Monte Carlo simulations were so slow that I had to rewrite many parts in C using CMEX. On the other hand, in my post-doctoral research where I studied estimation in distributed, non-linear systems, I decided to go with FORTRAN (using the 2008 standards to have at least a readable code…) so that the execution speed is not an issue. Yeah, the performance was very good, but it took me too much time to code. Then I heard about a new language that was promising the best of both worlds: something that resembles an interpreted language with speed of a compiled one! And that’s how I met Julia.

By that time (using v0.2 I think), Julia was a really new language. But I decided to accept the rough edges and try to code my algorithms using it. Anyway, it was just a personal side project to learn more about orbits. I did face many bugs, I had to use master (pre-v0.3) due to some bugs and missing features, but it was fun 🙂

After some years (and huge rewrites due to breaking changes), Julia released its v0.4. In this time, given the amount of code I had and the state of the language, I started to see that this bunch of algorithms can indeed be used for something at INPE to help in my activities. Hence, I decided to create a private package, which was called SatToolbox.jl, to organize everything I have done.

After some time, this little side project turned out to be a very important core for a simulator, called Forplan, of the mission operational concept we started to develop at INPE’s Space Mission Integrated Design Center (CPRIME). Given the good feedback I received, I decided to rename the toolbox to SatelliteToolbox.jl and it was released as an official Julia package in March 2018.

In this post, I would like to briefly describe SatelliteToolbox.jl and how it can be used to analyze space missions. The entire set of the features can be seen in the documentation. A brief, non-exhaustive list of the algorithms implemented by the time of this post (in v0.5.0) is:

• Earth atmospheric models:
• Earth geomagnetic field models:
• Space indices:
• Capability to automatically fetch many space indices, e.g. F10.7, Ap, Kp, etc.
• Functions to perform general analysis related to orbits, e.g. converting anomalies, computing perturbations, etc.
• Orbit propagators:
• Two body;
• J2;
• J4; and
• SGP4/SDP4.
• Functions to convert between ECI and ECEF references frames:
• All the IAU-76/FK5 theory is supported. Hence, the conversion between any of the following frames is available:
• 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.
• All the IAU-2006/2010 theory is supported. Hence, the conversion between any of the following frames is available:
• ITRF: International Terrestrial Reference Frame;
• TIRS: Terrestrial Intermediate Reference Frame;
• CIRS: Celestial Intermediate Reference Frame;
• GCRF: Geocentric Celestial Reference Frame.
• Functions to convert between Geocentric and Geodetic (WGS-84) references.

In the following, I provide a few examples of how SatelliteToolbox.jl can be used to analyze space missions.

## Installation

The very first thing (provided that you have Julia installed, which can be obtained here) is to install the package. This can be done by typing:

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

Information

If the package is already installed, then make sure you have at least the version v0.6.0. You can update all of your packages using the command Pkg.update().

To load the package, which must be done every time Julia is restarted, type:

julia> using SatelliteToolbox 

Information

Here and in the following, the command you have to type is what is presented after julia>, as you will see in the Julia REPL. Everything else is what you should see printed on the screen.﻿

## Examples

Now, I will show some analysis that can be done with the functions that are available.

Necessary background

To keep this post short, I will assume that you have knowledge about Julia and a little background in satellite and orbits.

### New Year in ISS

Let’s see how can we calculate where the astronauts onboard the ISS were during the New Year in Greenwich! The first thing we need to do is to get information about the ISS orbit. In this case, we must obtain the TLE (Two-Line Element) file, which is a data format consisting of two lines with 70 characters each that contain all information related to the orbit. The following TLE was obtained from Celestrak at January 4, 2019, 12:25 +0200.

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

This TLE must be loaded into a variable inside Julia. There are a number of methods to do this using SatelliteToolbox.jl. Here, we will use a special string type:

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²
==========================================================

This code loads the first TLE specified inside the string enclosed by tle"""...""" into the variable iss_tle.

Now, we have to initialize an orbit propagator using the loaded TLE. In this case, we will use the 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}
)

The variable orbp now holds the orbit propagator structure of type SGP4 with the orbit specified by the TLE iss_tle. This TLE was generated at the Julian day 2458487.75253 (2019-01-04 06:03:38.592 +0000). Thus, we have to backpropagate the orbit to the desired instant 2019-01-01 00:00:00.000 +0000 (New Year in Greenwich). This can be accomplished by the function propagate_to_epoch! as follows:

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])

The function propagate_to_epoch! returns three values. The first one o is the osculating orbit elements at the propagation instant, the second one r_teme is the position vector, and the last one v_teme is the velocity vector. Those vectors are represented on the same reference frame that was used to describe the orbit elements when the propagator was initialized. Since we are using the TLE, those vectors are represented in the TEME (True Equator, Mean Equinox) reference frame.

TEME is an Earth-Centered Inertial (ECI) reference frame. Hence, we must convert the position vector to an Earth-Centered, Earth-Fixed (ECEF) reference frame so that we can compute what was the ISS position (latitude, longitude, and altitude) at the desired instant. SatelliteToolbox.jl has the entire IAU-76/FK5 theory related to the conversion between reference frames. For this example, we will convert TEME into the International Terrestrial Reference Frame (ITRF) for more accurate computation. This kind of conversion requires the Earth Orientation Data (EOP) that is provided by IERS. SatelliteToolbox.jl can easily load and use this data as follows:

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

Warning

The EOP data is measured and published by IERS. This means that the ITRF cannot be used here to predict what will be the ISS position in the future since the EOP data will not be available. In this scenario, if high accuracy is not needed, then the Pseudo-Earth Fixed (PEF) reference frame can be used. The conversion between TEME and PEF using IAU-76/FK5 does not require external data.

The DCM (Direction Cosine Matrix) that rotates TEME into alignment with ITRF is computed by:

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

Thus, the position vector represented in ITRF is:

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

Finally, considering the WGS-84 reference ellipsoid, the latitude, longitude, and altitude of the ISS during the New Year in Greenwich can be obtained by the function ECEFtoGeodetic as follows:

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

-46.18935000852941

-112.3190646460004

julia> h/1000
419.8590733353887

i.e., latitude 46.189° S, longitude 112.319° W, altitude 419.859 km. This is in agreement with the historical information on I.S.S. Tracker website:

Information

The minor difference can be explained by the TLE we used to computed the ISS position and the one used by the historical data of I.S.S. Tracker. Since a TLE has errors and the propagation increases them, it is desirable to obtain a TLE as closest as possible to the desired instant, which was not pursued here.

## Atmospheric density profile

In this second example, we will use the built-in functions of SatelliteToolbox.jl to compute the atmospheric density profile. There are many models available in the literature. SatelliteToolbox.jl implements four of them: Exponential atmospheric model, Jacchia-Roberts 1971, Jacchia-Bowman 2008, and NRLMSISE-00. All of them but the former requires as input some space indices, like the F10.7 that measures the Sun activity and the Ap that measures the geomagnetic activity. SatelliteToolbox.jl is prepared to download all the required files from the Internet so that those indices can be easily obtained. This can be accomplished by:

julia> init_space_indices(wdcfiles_newest_year = 2018)
[ Info: Downloading file 'kp2018.wdc' from 'ftp://ftp.gfz-potsdam.de/pub/home/obs/kp-ap/wdc/kp2018.wdc'.

Information

The keyword wdcfiles_newest_year is not necessary, but it is added to avoid an error since the kp2019.wdc file was not available by the time this tutorial was written. For more information, see the documentation.

We will compute the atmospheric density profile from 100 km to 1000 km (steps of 1 km) using all the four models at 2018-11-1 00:00:00+0000 over São José dos Campos, SP, Brazil (Latitude 23.2237° S, Longitude 45.9009° W).

The exponential atmospheric model, the simpler one, depends neither on space indices nor on the location, only on the altitude. Thus, the atmospheric profile is computed by:

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

Each element is the atmospheric density [kg/m³] related to one altitude.

Information

Here we use the Julia broadcast operator . to compute for the entire altitude interval using only one line. For more information, see the documentation.

For the Jacchia-Robert 2008, we must specify the geodetic latitude [rad], longitude [rad], and altitude. Notice that, since we have already initialized the space indices, all the required information will be gathered automatically:

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

Each element is an instance of the structure JB2008_Output that contains the density of the atmospheric species in [kg/m³] related to one altitude.

The NRLMSISE-00 requires the same information, but in a different order. One more time, since we have already initialized the space indices, all the required information is fetched automatically:

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

Each element is an instance of the structure NRLMSISE00_Output that contains the density of the atmospheric species in [kg/m³] related to one altitude.

Online NRLMSISE-00

The values computed here can be compared to the online NRLMSISE-00 model. Notice, however, that the online version by the time this tutorial was written only allows computation between 1960/02/14 and 2018/03/17.

The Jacchia-Roberts 1971 model does not support fetching the space indices automatically yet. Hence, we will need to do this manually. It requires three indices: the daily F10.7, the averaged F10.7 (81-day window, centered on input time), and the Kp geomagnetic index (with a delay of 3 hours). That information can be fetched by:

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

Thus, the atmospheric profile computed by JR1971 is obtained by:

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

Each element is an instance of the structure JR1971_Output that contains the density of the atmospheric species in [kg/m³] related to one altitude.

Finally, using the PyPlot.jl package, the atmospheric profiles (altitude vs. density) in semi-log scale can be plotted using:

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()

Warning

If the desired date is not available on the space indices files, then an error will be thrown when computing the atmospheric density.

## Conclusion

I hope that this tutorial has helped you to understand a little bit how SatelliteToolbox.jl can be used to perform analysis related to satellites and orbits. If you have any question, please, feel free to leave a comment below!

Update 2019-03-24: Update the tutorial to match the version v0.6.0 of SatelliteToolbox.jl.

Update 2018-01-08: Fix algorithm that plots the data, because the variable h was not defined (thanks Bernard_GODARD).

# Changes in OrdinaryDiffEq.jl v3

Hi!

I know, it has been a long, long time since my last post. I would like to apologize! I have been very busy due to a massive amount of work to be done at INPE.

Anyway, this will be a very small post. I just want to update the source code of the last lesson so that it can work under Julia v0.6 and OrdinaryDiffEq.jl v3.

As I said, the package OrdinaryDiffEq.jl, which is part of the project JuliaDiffEq, is amazing. In my humble opinion, one of the best open-source differential equation solvers out there. The development of this suite is occurring at a very fast rate. This is good, because the package is becoming better each day. The downside is that API changes also very fast and sometimes your code needs refactoring to continue to work in the newer versions.

# Necessary background

This is just an update of the source code presented in the last post. Hence, I highly suggest to read it first:

# Changes in OrdinaryDiffEq.jl v3

Since I posted the last tutorial, OrdinaryDiffEq.jl API has changed (actually, the API is defined in DiffEqBase.jl, but it will be updated when you update OrdinaryDiffEq.jl). There are three modifications that we must perform so that the previous source code works correctly with Julia v0.6 and OrdinaryDiffEq.jl v3.0.3.

## Custom Data Array Interface

Previously, a custom data array interface could be created by means of a new structure that inherits the abstract type DEDataArray{T}. This have changed and now the correct abstract type is DEDataVector{T}. Hence, our SimWorkspace definition must be:

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

## Dynamic equation footprint

In the version used for the last tutorial, the footprint of the dynamic function was:

function f(t,u,du)


Due to some internal changes, the footprint now is:

function f(du,u,t,p)

p is a vector of parameters that is unused in this example, but simplifies a lot of things. If you want more information, please see http://docs.juliadiffeq.org/stable/tutorials/ode_example.html#Defining-Parameterized-Functions-1.

## Callback condition function footprint

The callback condition function footprint was also changed from:

function condition_control_loop(t,u,integrator) 

to:

function condition_control_loop(u,t,integrator)

# Updated source code

Finally, the updated source-code of the last tutorial can be seen below. It will provide exactly the same result.

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

# Conclusions

This tutorial just updated the source code of the last one so that it can work using Julia v0.6 and OrdinaryDiffEq.jl v3.0.3.

I am planning for my next tutorial a simple example of how to code Kalman filters in Julia language!

# Extending the ODE Solutions in Julia by Creating Custom Data Arrays for the Simulations

NOTE: The source code of this tutorial will not work under the latest versions of Julia and OrdinaryDiffEq.jl. You can see an updated version of the source code here.

Hi!

In this new post, I will show how we can create a Data Array to extend the possibilities of the julia package OrdinaryDiffEq to simulate dynamic systems.

# Necessary background

This tutorial provides additional information to solve a problem introduced by my last tutorial, then it is very important to read it before:

# Problem definition

In my last post, I showed how I use the julia package OrdinaryDiffEq to simulate control systems that have continuous and discrete parts. However, there were two big problems:

1. I needed to use a global variable to pass the information from the discrete part of the system to the dynamic function (in that case, the information was the acceleration);
2. I would have to do some very bad hacks to plot the computed acceleration (creating a global vector to store the computed values).

As I mentioned, the optimal solution to circumvent those problems is was somewhat difficult to achieve. Fortunately, Chris Rackauckas, which is the maintainer of OrdinaryDiffEq, introduced the concept of Data Arrays, which greatly simplified the solution for those problems. Notice that you will need v1.4.0 or higher of OrdinaryDiffEq in order to have access to the new Data Array interface.

In this tutorial, the problem we want to solve is: how to create custom parameters that the discrete part of the system can modify and that will influence the dynamic function?

# A little history

When I started to use OrdinaryDiffEq to simulate the attitude and orbit control subsystem of Amazonia-1 satellite, I saw that I needed two kinds of variables:

1. The ones that compose the state-vector and that will be integrated by the solver; and
2. Some parameters that the control loop can modify and will influence the dynamic solution, like the commanded torque for the reaction wheels and magneto-torquers.

However, by that time, there was no easy way to achieve this. The only method is to create a custom type in julia to store the parameters that are external to the state-vector. The problem is that you needed to implement yourself many methods in order to make this new type compatible with the OrdinaryDiffEq interface. In the following, you can see what you had to do so that you can define a parameter called f1. This code was used in this issue opened by me at 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)

Hence, Chris Rackauckas decided to create a new interface inside the OrdinaryDiffEq package so that you can easily build your own custom types to be used by the solver.

Note: Actually the definition of the new interface is placed at the package DiffEqBase, which is a dependency of OrdinaryDiffEq.

# Data Array interface

This interface allows you to define your own types to be used by the solver. In order to do this, you just need to define a type that inherits DEDataArray:

type MyDataArray{T} <: DEDataArray{T}
x::Array{T,1}
a::T
b::Symbol
end

where there must be an x array that defines the state-vector and any other parameter that you will need in your simulation. Given that, you can access any of those parameters inside the solver and callback functions using u.a, u.b, etc. You can also access the i-th element of the state vector using u[i]. In other words, you do not need to use the cumbersome notation u.x[i].

This pretty much replaces all the code I showed before 👍.

# Solving our problem

Hence, we now can remove the global variable a (🎉🎊🎉🎊🎉🎊) by defining the following type:

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

Inside our dynamic function, we can now get the commanded acceleration using u.a:

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

Now, we just need to modify the affect function of our callback to write on this variable:

# 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

You might have noticed that there is a strange for here. This is a problem created by the internal structure of the ordinary differential equation solvers. For example, 4-th order Runge-Kutta algorithm needs to use four caches called $$k_1$$, $$k_2$$, $$k_3$$, and $$k_4$$, which are essentially the dynamic function computed at different space-state points, that are used to obtain the solution at the next sampling step. In this case, the parameter a must be set in all those caches to provide the correct solution.  This is precisely what that for is doing. It is looping for all caches in the solver, and correctly setting the parameter a with the acceleration computed by this callback affect function. Notice that the current state-vector integrator.u is inside the list user_cache(integrator) and, consequently,  is also being set inside that loop. The good thing is that by using this for, you do not have to worry about how many caches the selected solver has.

Finally, we just need to create our initial state using the SimWorkspace type:

u0 = SimWorkspace(zeros(2), 0.0)

where the first parameter is initializing the state-vector with zeros (position and velocity), and the second parameter is initializing the acceleration. After the execution, we can access the commanded acceleration at each saved time step using sol[i].a.

If you are using PyPlot, then the state-vector and the acceleration can be plot using:

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

which should provide you the following result:

If you would like to change the initial state, then you just need to modify u0 before calling the solver:

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

or

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

which will provide you the following result:

# Conclusions

In this tutorial, I showed how you can use the Data Array interface of package OrdinaryDiffEq to create your own custom simulation workspaces. Hence, now you know how to create parameters to be used inside the solver that does not belong to the state-space vector.

I provided an example by extending the solution presented in a previous tutorial, where we needed a global variable to handle this kind of parameters. Notice that, although the presented example was very simple, this concept can be easily extended to provide you an easy, elegant way to simulate much more complex systems.

If I was not very clear, please, let me know in the comments! I will use all your advice to improve this and my future tutorials!

In the following, you can find the source code of this 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)

Edit (2017-02-28): Grammar corrections and minor text improvements.

# Using julia to simulate systems composed of continuous and discrete parts

NOTE: The source code of this tutorial will not work under the latest versions of Julia and OrdinaryDiffEq.jl. You can see an updated version of the source code here. Notice that the theory remains the same!

Hi!

In this first post, I would like to explain how I circumvent a problem I always had related to simulation of dynamic systems that has discrete and continuous parts. Notice that this is for sure MathWorks® Simulink territory. There is no doubt that Simulink is the faster prototyping tool for such kind of applications in which you have a continuous dynamic model modified by discrete parts that are sampled at specific time intervals.

However, sometimes there are situations in which you cannot use Simulink (not having a license being the most obvious one…). In my specific use case, I am coding the Attitude and Orbit Control Subsystem (AOCS) simulator for the Amazonia-1 Satellite. Simulink gave me some problems related to simulation of some discontinuities, like shadows of some satellite parts on the coarse solar sensors. It turns out that I needed more control in the solver than Simulink could offer me. Hence, I decided to port the simulator to a “standard” language. After a very long research, julia and the package OrdinaryDiffEq provided me the best solution to simulate my system. Hence, in this post, I would like to share how I was able to solve the problem of simulating systems with discrete and continuous parts using these tools.

# Necessary background

I will try to keep this tutorial as simple as possible. However, it would be helpful if you have a minimum knowledge about:

1. Julia language;
2. Control systems; and
3. Ordinary Differential Equations.

# Problem definition

In most control systems, you have four parts:

1. The plant, which is what you want to control;
2. The sensors, which will measure quantities related to the processes you want to control;
3. The actuators, which will act and modify the state of your plant; and
4. The controller, which will read the sensor outputs and command the actuators to modify the plant given a set of specifications.

However, nowadays the controller is usually embedded on a microcomputer, which has predefined intervals to execute its functions. In other words, only on specific time instants (called sampling steps) the controller will read the sensor output, update the control law, and send new commands to the actuators.

In order to simulate such system, you have to consider two kinds of models: continuous and discrete. The plant, for example, is a continuous model, since it is usually modeled by a set of differential equations and generally varies for every t. The controller, as mentioned before, only updates its internal state on a determined set of time instants.

Academically speaking, the simulation of such system can be solved by using the Z-transform to convert the plant continuous model to a discrete one. Hence, every model will be discrete and can be simulated easily. This method has many problems. The biggest one, in my opinion, is that it only works for linear systems and it can be very difficult to be used in complex plants with many state variables.

Hence, the problem is: how can I use a “standard” language to simulate a system in which you have a continuous part (the plant) that is controlled by a discrete process (the controller)? I will provide in this post one possible solution using julia and the package OrdinaryDiffEq.

### Proposed dynamic system

I decided to propose a very simple dynamic system to be controlled. This will help to make the tutorial more interesting and easy to understand. The problem I selected is a simple particle control in one dimension defined as follows:

Control the position of a particle that moves in one dimension. You can change the particle acceleration and you have access to the position and velocity measurements.

This particle can be modeled by a set of two differential equations describing the time-derivative of its position and velocity:

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

where $$p$$ is the particle position, $$v$$ is the particle velocity, and $$a$$ is the acceleration, which we can control.

Hence, we can define the state vector as:

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

Finally, our dynamic system can be written as:

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

Note: This is precisely one of those systems that can be easily simulated using the Z-transform approach when you have a discrete controller. But the idea is to keep this tutorial simple so that you can understand the idea to extend to more complex scenarios.

### Controller

We want to control the acceleration of a particle so that it goes to a specific point in space, called $$r$$. In this case, I will propose a (very) simple PD controller. The control law will be:

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

The controller must be a discrete process that is updated every 1s.

Note: For those unfamiliar with control systems, just assume that by computing the acceleration as shown here we can drive the particle to the desired position $$r$$.

In the next sections, I will describe how you can use julia and the package OrdinaryDiffEq to simulate this system in a very elegant solution. You will see that you can easily extend this idea to handle much more complex systems.

# Why did I choose julia?!

At my institute, INPE, I have to answer almost on a daily basis why I am using julia as the language to do my work. It is kind difficult to introduce new technologies especially when you are dealing with areas that people still use FORTRAN due to all those legacy codes.

I started to use julia since version 0.2 (we are currently on 0.5). The reasons why I really like this language are:

1. It seems an interpreted language in the sense that you do not need to mind about the types of the variables, which makes a MATLAB programmer very comfortable;
2. You can very easily call C / FORTRAN functions inside your julia code, this is important since models like IGRF and MSIS are made available in FORTRAN;
3. The performance of the simulations when vectorization is impossible (or very difficult) is much better than other interpreted languages (like MATLAB, Octave, Python, etc.);
4. You don’t need to mind to vectorize your code;
5. Parallel programming is now very easy to achieve; and
6. The syntax if very elegant (ok, I know, this is a very personal opinion 🙂 ).

Furthermore, for the kind of application I am focusing on this tutorial, the package OrdinaryDiffEq developed and maintained by Chris Rackauckas is a game changer. Chris developed an amazing interface to interact with the solvers that let you do almost everything. Conclusion: you have all the solvers you need (Runge-Kutta, Dormand-Prince, etc.) with the power to interact deeply with the algorithm to adapt it to your specific needs. I have to admit that I am not aware of any other package that provides such elegant solution to any other language (please, let me know in the comments if I missed something!).

# Setup

First, download the appropriate julia package for you operational system at http://www.julialang.org (notice that I tested everything here using openSUSE Tumbleweed and julia 0.5).

After installation, open julia and install the OrdinaryDiffEq:

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

# Simulating the system

The first thing is to load the modules we need:

using DiffEqBase
using OrdinaryDiffEq

Solving an ODE using OrdinaryDiffEq is very like to other solvers available in MATLAB, e.g. ode45. You need to define your dynamic function, set the initial state, adjust some options, and call the solver. Hence, let’s define our dynamic function.

In OrdinaryDiffEq, the dynamic function must have the following footprint:

function f(t,u,du) 

where t is the time, u is the state-vector, and du is the time-derivative of the state vector. Notice that u and du have the same dimension of the initial state. This function must compute du, but the user must be careful to not change the vector reference! Otherwise, the solver will not use the computed value.

For our proposed problem, the dynamic function will be:

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

Note 1: In this example, the acceleration variable will be global. This will make things easier, but it is a very (VERY) bad programming standard. However, the right way to avoid this global variable is a little bit more complicated and will be handled in a future tutorial.

Note 2: The .= sign means that the values on the left-side will be copied to the vector du without changing the reference.

Ok, the dynamic system model (the continuous part of our problem) is done! In the next, we will code our discrete control law.

### Discrete controller using Callbacks

Callbacks are one type of interface with the solver that makes OrdinaryDiffEq package powerful. You define two functions: one to check if an event occurred (condition function), and other to be executed in case this event has happened (affect function). Hence, at every sampling step, the solver calls the condition function and, if it tells that the event has occurred, then the affect function is executed.

OrdinaryDiffEq provides two kinds of Callbacks: Continuous and Discrete. In this case, we are interested in discrete callbacks in which we can precisely define prior to the execution what are the sampling steps that the affect function will be called. For more information, please see the documentation at http://docs.juliadiffeq.org/latest/features/callback_functions.html

Let’s define first the condition function that will check if we are on a sampling interval in which the control law must be updated:

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

The variable tf defines the simulation time, and the array tstops defines the sampling intervals of the control loop. The function condition_control_loop, our condition function, simply checks if t is one of our sampling instants. If this function return true, then the affect function will be called. Notice that it is also possible to define some kind of callback based on the state-vector u, but it is not our case.

Don’t mind with the integrator variable. It is basically an interface with the integrator that lets you check and change many, many things. I will provide more information in future tutorials. If you want to know more about it, check the documentation at http://docs.juliadiffeq.org/latest/basics/integrator.html

Now, let’s define our affect function that, in this case, is our controller:

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

where r is the reference position, and k and d are respectively the proportional and derivative terms of the PD controller. In this case, the controller loop only needs to compute the new acceleration, which is written in the global variable 😕 a, using the current position and velocity of the particle.

Now, we can define the callback:

cb = DiscreteCallback(condition_control_loop, control_loop!)

# Calling the solver

Finally, we just need to define our initial state and call the solver of OrdinaryDiffEq package to solve the problem for us:

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

Tsit5() means that we are using the Tsitouras 5/4 Runge-Kutta method. Other possible options are DP5() for Dormand-Prince’s 5/4 Runge-Kutta method, or BS3() for Bogacki-Shampine 3/2 method. All these methods have adaptive sampling steps, which is a good choice to simulate dynamic systems, but you can also use methods with fixed steps. For more information, see the documentation at http://docs.juliadiffeq.org/latest/solvers/ode_solve.html

After that, the solution of the problem will be stored in the variable sol. The time vector containing all sampling steps accepted by the solver is sol.t and the state vectors on each of these instants are stored in sol.u.

If you are using PyPlot, you can plot the solution this way:

plot(sol.t, sol.u)

The expected result of this simulation is:

Notice that you can easily see the points in which the control loop is computed, which immediately change the velocity of the particle.

# Conclusion

In this first tutorial, I showed how I use julia to simulate systems that have continuous and discrete parts. The concept was introduced by a very simple problem in which the position of a particle must be controlled using a discrete control loop.

I hope you could understand well the concept and how easy it will be to extend it to much more complex simulations. If I was not very clear, please, let me know in the comments! I will use all your advice to improve this and my future tutorials!

In the following, you can find the source code of this 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

• Edit (2017-02-26): Grammar corrections.