My first tutorial about julia: julia and mixed continuous/discrete control systems

General Usage
, ,
1

Hi guys!

Yesterday I wrote my first tutorial about julia. It is a very simple one, but describe a very big problem that I was having for quite sometime. Every time I needed to simulate a control system in which the control loop is discrete, I had to use Simulink. However, sometimes I needed more control of the solver than Simulink can give me. Normally, I ended up coding my own Runge-Kutta and doing the stuff I need.

However, few months ago, @ChrisRackauckas created an amazing interface in the package OrdinaryDiffEq that let you do almost everything inside the solver. This was a game changer for me. Now I have access to all those ODE solvers (with adaptive time step) with the possibility to interact so deeply with them that is possible to solve all the problems I had.

Julia and OrdinaryDiffEq are being currently used by me at INPE (the National Institute for Space Research at Brazil) to code the Attitude and Orbit Control Subsystem simulator of the Amazonia-1 Satellite, and to build a simulator to validate the operation concept of space missions.

I am really trying to show to everyone how good julia can be to scientific computing (I am even giving 0.5 point to my students if the homework is coded in julia :smiley: ). Hence, I decided to make a blog in which I will post little tutorials about my use cases and how julia solved my problems.

If you want to read:

http://www.ronanarraes.com/2017/02/using-julia-to-simulate-systems-composed-of-continuous-and-discrete-parts/

Any suggestions are welcome :slight_smile:

Thanks,
Ronan

6 Likes
2

Hey,
As always, I have a few suggestions :slight_smile:. For the dynamical equation, you should be using the in-place form, i.e. instead of returning du, it should update du:

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

If you want to return the output, you’d need to use the other signature

function dyn_eq(t,u)
    du = [0 1;
          0 0].*u + [0;1]*a # Note this will return du anyways because it's the last line
end

but the first is usually better for performance. Let me know if you think this needs to be documented more because this is a pretty crucial point. I'd at least update the post to have one of these two because I'm not sure the mixed signature works (in all cases, maybe there are times where it works accidentally!).

Also, for a problem like this where you aren't worried about saving the value of `a` at each timepoint, you should consider using a `ParameterizedFunction`. Here's how it is on master:

https://gist.github.com/ChrisRackauckas/49ab10bfbf872c4cfedc06f25aef06aa

I just realized that this will need a tag to work (I thought I tagged this part already!), so it's waiting on:

https://github.com/JuliaLang/METADATA.jl/pull/7847

This will do most of what your array stuff from before did, but it won't save what the value of `a` was. If the saving is important, you either need to do it manually or add the value to `u` like you had before. However, I plan on making using custom arrays much easier than it was before. I'll open an issue to get some ideas and and try to get this together ASAP since I think it's a nice feature.

Also, at the bottom of the latest docs

http://docs.juliadiffeq.org/latest/

You'll see a new section on `Modeling Tools`. Right now I have some packages (soon to be registered) which make it easier to solve some domain-specific problems, such as those which are multi-scale, based on financial markets, and systems biological models. If there's anything you think would help building control systems models, we could add a control systems domain and build some helper functions.
1 Like
3

Hi @ChrisRackauckas,

Thanks for your suggestions! I modified the post to use the cached version of dynamic equation. However, the parameterized functions I will let to another tutorial so that this one is kept simple :slight_smile:

This kind of parameterized functions is very good! In my case, it can avoid creating types for the state-vector for more simple systems.