odehybrid | An Uncommon Lab

Overview

odehybrid is a set of tools for building simulations of hybrid continuous- and discrete-time dynamical systems in MATLAB. Many systems fall into this category, such as control systems, where the physical process being controlled is continuous but the controller's logic and math are discretely executed and updated at some period.

In general, continuous-time systems are ordinary differential equations (ODEs) and can be written as:

$$\dot{x}(t) = f(t, x(t), \ldots)$$ where $x$ is the continuous state (e.g., a column matrix with position and velocity).

Discrete-time systems are updated at some period, $\Delta t$, as follows:

$$y(t + \Delta t) = g(t, y(t), \ldots)$$ where $y$ is the discrete state (e.g., an acceleration command or mode logic).

These systems interact, such as:

$$\dot{x}(t) = f(t, x(t), y(t), \ldots)$$ $$y(t + \Delta t) = g(t, x(t), y(t), \ldots)$$

By expanding on these simple ideas, one can build large simulations with many systems of systems. But first, let's examine a very quick simulation.

Quick Example

The primary function is called odehybrid, and this works in an analogous manner to MATLAB's built-in ode45. For example, here is a tiny simulation of an unstable continuous-time process with a discrete controller running at 10Hz.

$$ \dot{x}(t) = \begin{bmatrix} 0 & 1 \\ 2 & 0 \end{bmatrix} x(t) + \begin{bmatrix} 0 \\ 1 \end{bmatrix} u(t) $$ $$ u(t + \Delta t) = -\begin{bmatrix} 8 & 4 \end{bmatrix} x(t) $$ where $u$ is held constant from $t$ to $t + \Delta t$ (sample and hold), and $\Delta t$ is always 0.1s (10Hz). The simulation for this is quick:

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); hold on; stairs(tu, u, 'r'); xlabel('Time');    % Plot 'em.
legend('x_1', 'x_2', 'u', 'Location', 'se');                % Label 'em.

There are many, many more examples on the Examples page.

More States

In the example above, the states were simple: a column vector and a scalar. However, odehybrid can handle quite a bit more. For instance, one can have multiple continuous and discrete states with different types, such as matrices, cell arrays, or even structs. Suppressing dependence on time for clarity:

$$ \begin{bmatrix} \dot{x}_1 & \dot{x}_2 & \ldots \end{bmatrix} = f(t, x_1, x_2, \ldots, y_1, y_2, \ldots) $$ $$ \begin{bmatrix} x_1 & x_2 & \ldots & y_1 & y_2 & \ldots \end{bmatrix} \leftarrow g(t, x_1, x_2, \ldots, y_1, y_2, \ldots) $$ That is, $f$ is a function that returns the derivative of all of the continuous states, and $g$ is a function that can update all of the states (both continuous and discrete) at some period, $\Delta t$. In this case, the interface becomes:

% ODE
[xd1, xd2, ...] = f(t, x1, x2, ..., y1, y2, ...) ...;

% Discrete-time update
[x1, x2, ... y1, y2, ...] = g(t, x1, x2, ..., y1, y2, ...) ...;

% Call odehybrid to simulate and return histories of x1, x2, ... y1, y2.
[t, x1, x2, ..., td, y1, y2, ...] = odehybrid(solver, @f, @g, dt, ts, ...
                                              {x1, x2, ...}, .,..
                                              {y2, y2, ...});

See: odehybrid

Continuous Time

Propagating a continuous-time dynamical system over time requires the use of ordinary differential equation solvers. These update the state by taking derivatives (calling $f$) at various points in time so that an accurate picture of the change of $x$ can be determined. A fantastic example of a continuous solver is MATLAB's built-in ode45. This function is a variable-step solver, so one does not specify a time step. Instead, it uses the output of $f$ at various points to determine how large a time step can be justified.

MATLAB has many variable-step solvers built-in, and one can use any of these for the continuous solver, such as the use of ode45 for the example above. However, these are meant to propagate a system for a “long” amount of time, and so they will always take at least a certain number of steps. When a continuous-time process is frequently interrupted by a discrete-time process, the number of steps can be far too large. For this reason, odehybrid includes an adaptive-step ODE solver, rkadapt, which has an indentical interface to ode45. By default, it uses the Runge-Kutta method known as Dormand-Prince 5(4) (fifth order method with fourth order error-checking). The function can further take any explicit Runge-Kutta adaptive-step Butcher tableau as inputs, allowing it to be used for any explicit adaptive-step solver.

For some systems, adaptive-step solvers can be overkill, such as those with dynamics that are slower than the discrete-time update rate or whose dynamics don't change much. In these cases, one can use a faster fixed-step solver. Fixed-step solvers still use numerous outputs from $f$ at various points to inform the state update, but do so with with a specified time step, which eliminates some of the extra calls to $f$ that are necessary with adaptive-step methods. The most common strategy in this class is known as “the” Runge-Kutta fourth order method (though there are many well developed Runge-Kutta fourth-order methods). This is implemented in rk4. Further, the function rkfixed allows one to input any Butcher tableau for an explicit fixed-step Runge-Kutta method.

These can be used as the solver quite simply as:

[...] = odehyrbid(@rkadapt, ...);

For fixed-step solvers, one needs to specify the time step using the odeset options, which are the eighth input to odehybrid (right after the initial states).

options = odeset('MaxStep', 0.1); % Set 0.1s time step for ODE.
[...] = odehyrbid(@rk4, ..., options);
[...] = odehyrbid(@rkfixed, ..., options);

Discrete Time

Discrete updates are much easier to propagate than continuous-time systems. More or less, the discrete function is called at the specified rate, updating the discrete states and potentially the continuous-states. For example, if the discrete period were 1s and the simulation started at 0s., then the simulation would start with a call to the discrete function, updating the states from the initial state. It would call the continuous solver from 0s to 1s, accepting the updated continuous states as a result. It would then call the discrete function with the updated continuous states and the last discrete states. It would then propagate on to 2s and repeat this process until the end time is reached.

odehybrid allows multiple discrete-time update functions at various rates. Each function operates on all states. For example:

% ODE
[xd1, xd2, ...] = f(t, x1, x2, ..., y1, y2, ...) ...;

% Discrete-time updates
[x1, x2, ... y1, y2, ...] = g1(t, x1, x2, ..., y1, y2, ...) ...;
[x1, x2, ... y1, y2, ...] = g2(t, x1, x2, ..., y1, y2, ...) ...;

% Call odehybrid to simulate and return histories of x1, x2, ... y1, y2.
[t, x1, x2, ..., td, y1, y2, ...] = odehybrid(solver, ...
                                              @f, ...
                                              {@g1, @g2}, [dt1, dt2], ...
                                              ts, ...
                                              {x1, x2, ...}, .,..
                                              {y2, y2, ...});

State Histories

The outputs of odehybrid are histories of each state. For instance, if a state were a vector of length $n$, the output would be a matrix of size $s$-by-$n$, where $s$ is the number of samples. If the state is a struct, the output will be an $s$-by-1 struct array. Otherwise, the state history will be an $s$-by-1 cell array, which element of which contains the state at the specified time. This therefore mimics the ode45 output type but extends it for many more types.

A consequence of allowing the discrete updates to update the continuous-states as well is that the states will have two values at every discrete update — the value before the update and the value after. If one goes to interpolate with this data using, e.g., interp1, it won't work, because it won't know how to interpolate across those time steps which have two values. For instance, the histories may look like this:

t = [... 1.5; 1.7;   2; 2; 2.4; 2.6; ...];
x = [...   0; 0.1; 0.2; 1; 1.1; 1.2; ...];

Since x has two values at 2s, normal interpolation won't work. odehybrid includes an interpolation utility, interpd, to allow one to use either the “left” or “right” value for interpolation right at 2s, so it can handle this case just fine.

xi = interpd(t, x, ti, '-'); % "keep left" interpolation

See: interpd

Logging

When the systems become large and complex, simply outputting a state history isn't enough to understand what's going on inside the system. In these cases, one can use the included TimeSeriesLogger. It has a simple interface and can add new signals on the fly.

log = TimeSeriesLogger();
log.add('alpha', t, a); % Add a new point to the signal known as 'alpha'.
...
log.add('alpha', t, a); % Add another point at a later time.

When a logger is the final argument to odehybrid, it's provided to the ODE and the discrete updates. They can then log whatever is necessary:

% Create a logger.
log = TimeSeriesLogger();

% ODE
function [xd1, xd2, ...] = f(t, x1, x2, ..., y1, y2, ... log)
    ... 
    if nargin == n                  % See if a log was passed in.
      log.add('name', t, y2 - x1);  % Log something.
    end
end

% Discrete-time update
[x1, x2, ... y1, y2, ...] = g(t, x1, x2, ..., y1, y2, ... log) ...;

% Call odehybrid to simulate and return histories of x1, x2, ... y1, y2.
[t, x1, x2, ..., td, y1, y2, ...] = odehybrid(solver, @f, @g, dt, ts, ...
                                              {x1, x2, ...}, ...
                                              {y2, y2, ...}, ...
                                              options, ...
                                              log);

Afterwards, one can automatically plot the logged signals.

log.plot();

There are numerous options for logging, accessing a logged signal directly and controlling what gets plotted and in what figure.

One critical thing to note: when calling the ODE function, ($f$), the log will not always be passed in. The reason for this is that an adaptive-step solver may “try out” a time step, see that the time step was too big, and abandon it for a smaller time step. If somewhere were to be logged during this call, it wouldn't make sense (it would be from a discarded time step). For this reason, one must see if the log is provided in one's ODE, such as:

function xd = an_ode_to_joy(t, x, y1, y2, log)
    xd = ...;
    if nargin == 5
      log.add('var', t, some_local_var);
    end
end

Discrete updates will always receive the log as an input, because discrete steps never need to be discarded.

See: TimeSeriesLogger

Honoring the Past

There's quite a bit of good work in MATLAB's ODE functions, with many features exposed via the options input. odehybrid can take these options as well and passes them on to the continuous solver. For instance, one can use events (zero-crossing detection) and an output function (live plot updates) via these options. Because odehybrid works with numerous states, it will map these states to event and output functions, just like they're mapped to the ODE and discrete update functions. That is, an event function will have the following interface:

[z, terminal, direction] = event_fcn(t, x1, x2, ..., y1, y2, ...) ...;

and an output function will have a similar interface:

status = output_fcn(t, x1, x2, ..., y1, y2, ... flag) ...;

See odeset in the MATLAB documentation for more. (Note that the included continuous solvers, rk4, rkfixed, and rkadapt will use the output function, but do not use events.

Docs