Chapter 3 Compilation and Simulation

This chapter explains the basic principles behind the compilation and execution of Zélus programs. Understanding these details is useful for linking Zélus and OCaml code, either by importing OCaml values into Zélus or by exporting compiled Zélus code into OCaml.

3.1 Combinatorial Functions

Combinatorial functions are defined by the form:

[ let ] [ fun ] ident pattern = result-expr

They are compiled directly into OCaml functions without any particular transformations.

3.2 Sequential Nodes

Synchronous, possibly stateful, functions are defined by the form:

[ let ] [ node ] ident pattern = result-expr

They are translated as follows:

A record type ident is introduced to represent the current value of the node’s internal state. Such records contains fields for each unit delay and -> within the node, and also for every other node instantiated within it. Values of this type are thus trees of nested records that reflect the structure of instantiations.
A step function ident_step is introduced. Calling ident_step self pattern, where self is the current state of the node (of type ident) and pattern passes the current value of the input, returns the current output of the node and directly updates the state record self.
A reset function ident_reset is also introduced. Calling ident_reset self resets the state self to its initial value.
Finally, there is also an allocation function ident_alloc () which returns an initial state (of type ident).

3.3 Hybrid Nodes and Interaction with a Solver

Figure 3.1: Basic structure of the hybrid simulation algorithm

Hybrid nodes are also compiled into several functions. Their execution requires run-time support to coordinate phases of integration by a numeric solver with the execution of discrete reactions. Before outlining the form of compilation, we describe the principles behind the simulation loop provided by the run-time.

The simulation loop of a hybrid system can be formalized as a synchronous function that defines four streams lx(n), y(n), z(n) and t(n), with n ∈ ℕ. Here, t(n) is the increasing sequence of simulation times at which the solver stops;¹ lx(n) is the value at time t(n) ∈ ℝ of the continuous state variables, that is, the variables defined by their derivatives in the original model; y(n) is the value at time t(n) of the model’s discrete state; and z(n) signals any zero-crossings occurring at instant t(n).

The simulation loop function has two modes: discrete (D) and continuous (C). The two modes and the transitions between them are sketched in figure 3.1. The function begins in D, where all computations that may change the discrete state or that have side effects are executed. Several discrete steps may be executed without advancing the simulation time t in what is termed an event cascade. Simulation time is advanced in C where a numeric solver is invoked to approximate the values of ODEs within the model. The function loops back to the D mode whenever a zero-crossing occurs (section 2.2.2) or a timer horizon (section 2.2.3) is reached. We now formalize the details of the two modes.

The Continuous Mode (C)

In the continuous mode C, a numeric solver computes an approximation of the solution of the set of (active) ODEs within a model. We represent the solver abstractly as a function solve(f)(g) parameterized by two other functions f and g:

x′(τ) = f(y(n), τ, x(τ)) yields the derivatives of continuous state variables x at τ ∈ ℝ, given the discrete states y(n) and an approximation to the continuous states x(τ);
upz(τ) = g(y(n), τ, x(τ)) computes the value of the zero-crossing expressions upz given the same arguments as f.

For a given f and g, the continuous mode C computes s, lx, z, and t such that:

(lx, z, t, s)(n+1) = solve(f)(g)(s, y, lx, t, step)(n)

where

s(n)

is the internal state of the solver at instant t(n) ∈ ℝ. Calling solve(f)(g) updates the state to s(n+1).

lx(n)

is the value of x at instant t(n), that is, lx(n) = x(t(n)); lx is a discrete-time signal while x is a continuous-time signal.

t(n+1)

is bounded by the horizon given to the solver t(n) + step(n), that is,

t(n) ≤ t(n+1) ≤ t(n) + step(n)

x

is a solution of an ODE, parameterized by the current discrete state y(n), that is:

∀ T ∈ [t(n), t(n+1) + margin] x(T) = lx(n) +

∫

t(n)

f (y(n), τ, x(τ))) dτ

Integration is performed on an interval strictly larger than [t(n),t(n+1)] to account for possible zero-crossings at instant t(n+1). The constant margin is considered as a parameter of the solver. In practice, solve(f)(g) computes a finite sample of approximated values on the interval [t(n), t(n+1) + margin]. The value of lx(n+1) is

lx(n+1) = x(t(n+1))

z(n+1)

signals any zero-crossings detected at time t(n+1). It comprises k ∈ ℕ boolean elements z(n+1)(i) such that:

z(n+1)(i) =

(∀ T ∈

⎡
⎣

t(n), t(n+1)

⎞
⎠

upz(T)(i) < 0)

∧

∃ m ≤ margin. (∀ T ∈ [t(n+1), t(n+1)+m] upz(T)(i) ≥ 0)

This definition assumes that the solver also stops whenever a zero-crossing expression passes through zero from positive to negative. An event occurs with a transition to the discrete phase (D) when a zero-crossing occurs or the horizon t(n)+step(n) is reached:

event = z(n+1)(0) ∨ … ∨ z(n+1)(k) ∨ t(n+1) = t(n) + step(n)

If the solver raises an error (for example, a division by zero occurs, or a suitable approximation cannot be found), we consider that the simulation fails.

The Discrete Mode (D)

All discrete changes occur in this mode. It is entered when an event occurs during integration. During a discrete phase, the function next defines y, lx, step, encore, z, and t:

(y, lx, step, encore)(n+1)	= next (y, lx, z, t)(n)
z(n+1)	= false
t(n+1)	= t(n)

where:

y(n+1): is the new discrete state (outside of mode D, y(n+1) = y(n));
lx(n+1): is the new continuous state, which may be directly changed in the discrete mode;
step(n+1): is the new step size;
encore(n+1): is true if an additional discrete step must be performed. Function next can decide to trigger another discrete event instantaneously causing an event cascade.
t(n): (the simulation time) is unchanged during a discrete phase.

The initial values for y(0), lx(0) and s(0) are given by an initialization function init. Finally, solve(f)(g) may decide to reset its internal state when the continuous state changes. If init_solve(lx(n), s(n)) initializes the solver state, then we have:

reinit	= (lx(n+1) ≠ lx(n))
s(n+1)	= if reinit then init_solve(lx(n+1), s(n)) else s(n)

These definitions give a synchronous interpretation of the simulation loop, seen as a stream function that iteratively computes the sequences lx, y and t at instant n+1 relative to their values at instant n and an internal state. By writing solve(f)(g), we abstract the actual choices of the integration method and zero-crossing detection algorithms. A more detailed description of solve(f)(g) would be possible (for example, as an automaton with two states: one that integrates, one that detects zero-crossings) but these details have no influence on the code generation problem which is independent of such simulation details.

Hybrid functions are defined by the form:

[ let ] [ hybrid ] ident pattern = result-expr

They are compiled into six functions:

A record type ident is introduced to represent the current value of the function’s internal state. In addition to fields for each unit delay, -> operator, and node instantiation within the function body, it also contains fields for continuous state variables (those defined with der) and zero-crossing expressions (those registered with up()).
A step function ident_step is introduced. Calling ident_step self pattern, where self is the current state of the node (of type ident) and pattern passes the current value of the input, returns the current output of the node and directly updates the state record self.
The abstract function f, in the above formalization, is realized by the function ident_derivatives self pattern that computes the current values of derivatives from the current state and stores them in the internal state self.
The abstract function g, in the above formalization, is realized by the function ident_crossings self pattern that computes the current values of zero-crossing expressions and stores them in the internal state self.
A reset function ident_reset is also introduced. Calling ident_reset self resets the state self to its initial value.
Finally, there is also an allocation function ident_alloc () which returns an initial state (of type ident).

1: They are called major time steps in Simulink.