# A Simulation Engine

I write a lot of simulations. For me, it's critical that I can pump out a solid simulation quickly and know that it's accurate, flexible, and fast. The best simulation engine out there is probably Simulink from The MathWorks, though others are doing interesting things too, such as SystemModeler from Wolfram. However, these are often geared towards building large, complex simulations in a graphical way. That's not always what I want. Let's consider:

Suppose I were to create a simulation of a satellite's orbit, attitude, and control. I'd probably pop open Simulink and start drawing blocks, especially if I needed the particular strengths of that platform. But now suppose I needed to simulate n interacting smallsats in a formation, but n was an input. Woah, suddenly that's hard to do in a platform where you draw each system as a little box. I would have to manually draw boxes for the 2-sat config, the 3-sat config, 4-sat, .... You can't just spontaneously reconfigure everything like you can with raw code. Plus, I often need the systems to function both in the simulator and outside of it in ways that just aren't possible with these platforms. Finally, if I need to be able to run the simulation tons of times, e.g. to train a controller or do post-flight analysis, these big simulation platforms kind of get in the way. They have serious strengths, and I love using them to their full potential, but when I need something else, I don't want to try to hack it out in a giant platform. I want something lean and easy and fast.

So let's take a step back. If I just needed a simulation of a simple physical state vector (like position and velocity), and I had a dynamics function that produces the time-derivative of that state, such as:

then of course I could simulate this easily with MATLAB's ode45 function. It's light and quick, like what I'm looking for, but it's limited to very small problems and lacks many of the expected features of a "simulator". For instance, what if I have both continuous and discrete processes (like a physical system with a digital controller)? And what if I want want to have multiple, separate states, like the below?

[x_dot, y_dot, z_dot] = f(t, x, y, z);

And what if x were a matrix, y were a struct, and z were a cell array? And what if I want easy-peasy logging? Now we're asking way, way more than our good friend ode45 was designed to handle.

Enter odehybrid.

For my various projects, I started creating my own MATLAB-based simulation platform -- like ode45, but expanded for multiple discrete-time systems, multiple inputs with multiple types, logging, and some custom solvers. This grew along with projects I was working on, always tied to real things that were actually needed in a simulation somewhere, and finally I believe it has filled the gap between the simple continuous-only solvers like ode45 and giant platforms like Simulink or SystemModeler. It's finally ready to be released as an open-source project, and as of the last commit, was called odehybrid, so that's the name that has stuck. Let's see what it does.

First, the ordinary differential equation should have the form:

[x1_dot, x2_dot, ... xn_dot] = ...
f(t, x1, x2, ... xn, y1, y2, ... yn);

where the xk variables are continuous-time states and the yk variables are the discrete-time states. So the ode, f, returns the time-derivative of each of the continuous-time states.

It also works with discrete update equations with the form:

[x1, x2, ... xn, y1, y2, ... yn] = ...
g(x1, x2, ... xn, y1, y2, ... yn);

That is, the discrete update equation can update the discrete-time state and the continuous-time states.

Hang on. What? The continuous stuff can be discretely updated? It can. Consider a satellite changing orbits with a really quick burst of rocket thrust. If we're simulating for years, we don't want to simulate the rocket burst at the 0.01s level necessary to capture the extreme stiffness of the rocket burst. It would be way easier (and not particularly less accurate) to just add the right right when the satellite's controller says to do so. I'll grant this sounds odd at first, but it's super convenient. And if you don't want to change a continuous-time state in a discrete update, just don't change it. It will pass from input to output, and everything will be just like it was.

Ok, cool. For the next step, we want the states to not just be state vectors. That's fine! Even the continuous states needn't be state vectors! Any state can be a matrix with any dimensions, a struct, struct array, or a cell array, each element of which can be any of these types, etc. So, we can keep the position and velocity in a vector, the actuator states in a matrix, the guidance module state in a struct, and the various individual controller states in a cell array.

How does one return the derivative of a struct? Simple, just let each field now mean the time-rate-of-change of that field. E.g., if a state is:

x.p = 1;
x.v = 2;

we simply set up a derivative struct like so:

x_dot.p = x.v;
x_dot.v = -9.81;


Now we just pass these functions over to odehybrid along with the simulation time range, initial states, which solver we'd like to use for the continuous-stuff, and the time step for the discrete update. Here's the interface:

[t, ...                    % Output continuous time-steps
x1, x2, ...,              % Output continuous states
td, ...                   % Output discrete time-steps
y1, y2, ...] = ...        % Output discrete states
odehyrbid( ...
@ode45, ...        % Solver
@f, ...            % ODE
@g, ...            % Discrete update
0.1, ...           % Discrete update rate
[0 100], ...       % Time span
{x1, x2, ...}, ... % Initial continuous states
{y1, y2, ...});    % Initial discrete states


This little interface can do quite a lot. For instance, here's the flight simulator for my autonomous aircraft:

ode = @(varargin) aircraft_dynamics(varargin{:}, ...
options, ...
constants);
de  = @(varargin) control_wrapper(varargin{:}, ...
options, ...
constants);
dt  = get_control_rate(); % Returns a scalar.
ts  = [0 t_end];
xc0 = x;
xd0 = {actuators, environment, nav, control};
options = odeset(...);

% Simulate!
[t, x, td, act, env, nav, ctl] = ...
odehybrid(@rkadapt, ode, de, dt, ts, xc0, xd0);


where I'm using the included rkadapt function for the ODE solver (which is faster for this type of propagation than ode45).

As a bonus, it can run multiple discrete updates and multiple rates. Instead of passing in @g, we'd just pass in {@g1, @g2, ... @gn} and instead of dt, we'd pass in [dt1, dt2, ... dtn].

We can also add a flexible logger to the mix, and odehybrid will take care of passing it to the ode and de functions correctly.

log = TimeSeriesLogger();         % Make a logger.
ode = @(t, x1, ... yn, log) ... ; % Make an argument for it
de  = @(t, x1, ... yn, log) ... ; %   in the ODE and DE.
[...] = odehybrid(..., log);      % Give it to odehybrid.


Here's the last thing: it's open-source. I've benefited from good open-source tools, and this is a good place for me to contribute in my particular knowledge areas.

There's a lot of info over at odehybrid, but for a here-it-is-running scenario, consider this simple simulation of an unstable continuous-time process (something trying to run off to infinity) being stabilized by a discrete-time controller at 10Hz (control folks will recognize this as a PD controller). This is 100% of the code necessary for the simulation.

ode = @(t, x, u) [0 1; 2 0] * x + [0; 1] * u;  % Differential equation
de  = @(t, x, u) deal(x, -[8 4] * x);          % Discrete update equation
dt  = 0.1;                                     % Discrete eq. time step
ts  = [0 5];                                   % From 0 to 5s
x0  = [1; 0];                                  % Initial continuous state
u0  = 0;                                       % Initial discrete state
[t, x, tu, u] = odehybrid(@ode45, ode, de, dt, ts, x0, u0); % Simulate!
plot(t, x, tu, u, '.'); xlabel('Time');                     % Plot 'em.
legend('x_1', 'x_2', 'u', 'Location', 'se');                % Label 'em.


It's available on File Exchange.

If down the road somewhere, you're using odehybrid, and you think, "Hey, this would be better with...," then get in touch! I'll send notes about the organization of the code so you can make the change and contribute to the project, and I'll also update the project when I find a solid addition. The goal is for the project to grow along with actual need in practice, while maintaining the simple and very general-purpose syntax.