hms/hsim
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

feat (sim): greedy simulation

Browse source
This commit is contained in:
Henri Saudubray 2025-04-22 17:57:10 +02:00
parent b4a29bbb97
commit 3d317f65a0
Signed by: hms
GPG key ID: 7065F57ED8856128
5 changed files with 421 additions and 262 deletions
Download patch file Download diff file Expand all files Collapse all files

145
doc/hsim.tex
View file

@ -1,145 +0,0 @@
\documentclass[a4paper]{article}
\usepackage{fullpage}
\usepackage{listings}
\usepackage{xcolor}
\title{A formalization of a simulation engine for hybrid systems}
\author{}
\date{}
\setlength{\parindent}{0pt}
\lstdefinelanguage{ml}{
basicstyle=\ttfamily,
morekeywords=[1]{
type, float, option, let, fun, with, if, else, then, in, as, match,
},
morekeywords=[2]{
HNode, DNode, Idle, Running, Discrete, Continuous
},
keywordstyle=[1]{\color{blue}},
keywordstyle=[2]{\color{red}},
commentstyle=\itshape,
columns=[l]fullflexible,
sensitive=true,
morecomment=[s]{(*}{*)},
keepspaces=true,
literate=
{'a}{$\alpha$}{1}
{'b}{$\beta$}{1}
{'p}{$\rho$}{1}
{'s}{$\sigma$}{1}
{'y}{$y$}{1}
{'yder}{$\dot{y}$}{1}
{'zin}{$z_{in}$}{1}
{'zout}{$z_{out}$}{1}
{fun\ }{$\lambda$}{1}
{->}{$\to$}{1}
{+.}{$+$}{1}
{-.}{$-$}{1}
{=}{$=$}{1}
{>=}{$\geq$}{1}
{<=}{$\leq$}{1}
}
\lstnewenvironment{ml}{\lstset{language=ml}}{}
\newcommand{\mlf}[1]{\lstinline[language=ml]{#1}}
\begin{document}
\maketitle
A discrete synchronous function, or node, can be seen as a pair of a step and
reset function, which operate on an inner state:
\begin{ml}
type ('p, 'a, 'b) dnode =
DNode : { s : 's;
step : 's -> 'a -> 'b * 's;
reset : 'p -> 's -> 's } -> ('p, 'a, 'b) dnode
\end{ml}
A hybrid node is quite similar: it has an inner state, a step and a reset
function; but the step function is decomposed into multiple distinct elements
for the purpose of the simulation:
\begin{ml}
type ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode =
HNode : { s : 's;
step : 's -> 'a -> 'b * 's;
fder : 's -> 'a -> 'y -> 'yder;
fzer : 's -> 'a -> 'y -> 'zout;
fout : 's -> 'a -> 'y -> 'b;
reset : 'p -> 's -> 's;
horizon : 's -> time;
jump : 's -> bool;
cget : 's -> 'y;
cset : 's -> 'y -> 's;
zset : 's -> 'zin -> 's;
} -> ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode
\end{ml}
\mlf{step} and \mlf{reset} are on discrete steps, and behave in the same way as
for a discrete node. \mlf{fder}, \mlf{fzer} and \mlf{fout} are used during
integration phases with a solver. \mlf{fder} represents the derivative function
in an initial value problem. \mlf{fzer} is the zero-crossing function, which
computes the values of all zero-crossing function we wish to monitor. \mlf{fout}
computes the output of the system at a particular instant.
\begin{figure}
\begin{ml}
let sim (HNode model) (DNode solver) =
let s = { status = Idle; mstate = model.state; sstate = solver.state } in
let step state input =
match input, state.status with
| Some i, _ ->
let status = Running
{ mode = Continuous; input = i; now = 0.0; stop = i.length } in
None, { state with status }
| None, Idle -> None, state
| None, Running ({ mode = Discrete; _ } as r) ->
let o, mstate = model.step state.mstate (r.input.u r.now) in
let state =
let h = model.horizon mstate in
if h <= 0.0 then { state with mstate }
else if r.now >= r.stop then s
else if model.jump mstate then
let y = model.cget mstate in
let fder t = model.fder mstate (offset r.input r.now t) in
let fzer t = model.fzer mstate (offset r.input r.now t) in
let ivp = { fder; stop = r.stop -. r.now; init = y } in
let zc = { yc = y; fzer } in
let sstate = solver.reset (ivp, zc) state.sstate in
let status = Running { r with mode = Continuous } in
{ status; mstate; sstate }
else { state with status = Running { r with mode = Continuous } }
in Some { start = r.now; length = 0.0; u = fun _ -> o }, state
| None, Running ({ mode = Continuous; _ } as r) ->
let (h, f, z), sstate = solver.step state.sstate r.stop in
let mstate = model.cset state.mstate (f h) in
let now = r.input.start +. h in
let state = match z with
| None ->
let status =
if h >= r.stop then Running { r with mode = Discrete; now }
else Running { r with now } in
{ status; mstate; sstate }
| Some z ->
let status = Running { r with mode = Discrete; now } in
{ status; mstate = model.zset mstate z; sstate } in
let fout t =
model.fout mstate (r.input.u (r.now +. t)) (f (r.now +. t)) in
let out =
{ start = r.input.start +. r.now; length = h -. r.now; u = fout } in
Some out, state in
let reset (m, s) { mstate; sstate; _ } =
let mstate = model.reset m mstate in
let sstate = solver.reset s sstate in
{ status = Idle; mstate; sstate } in
DNode { state = s; step; reset }
\end{ml}
\label{fig:ml:sim}
\caption{Hybrid System Simulation}
\end{figure}
\end{document}

363
doc/notes.typ
View file

@ -1,17 +1,36 @@
#import "@preview/cetz:0.3.2" #import "@preview/cetz:0.3.2"
#import "@preview/cetz-plot:0.1.1" #import "@preview/cetz-plot:0.1.1"
#import "@preview/finite:0.4.1"
#set page(paper: "a4") #set page(paper: "a4")
#set par(justify: true) #set par(justify: true)
#let breakl(b) = {
raw(b.text.replace(regex(" *\\\\\n *"), " "), lang: b.lang, block: b.block)
}
#show raw: set text(font: "CaskaydiaCove NF")
#show link: underline
= Notes on hybrid system simulation = Notes on hybrid system simulation
== TODO
Multiple different axes of research:
- abstract over the notion of internal state and of our ability to copy some or
part of it
- make free variables contained in solver closures explicit, and require the
ability to copy them for greedy simulation (without the solver itself knowing
about them?)
- think about the notion of input values with overlapping domains
- think about successive input values carrying information about possible
discontinuities
== The problem == The problem
The initial problem is _transparent assertions_. We want to simulate assertions The overarching problem is _transparent assertions_. We want to simulate
with their own solvers, because simulating assertions with the same solver as assertions with their own solvers, because simulating assertions with the same
the parent system changes the final results, when assertions are supposed to be solver as their parent system changes the final results, when assertions are
transparent. supposed to be transparent.
Assertions are themselves hybrid systems, and can contain their own assertions. Assertions are themselves hybrid systems, and can contain their own assertions.
We thus need a semantics that is both *independent* of the solver, which becomes We thus need a semantics that is both *independent* of the solver, which becomes
@ -25,73 +44,337 @@ type ('p, 'a, 'b) dnode = DNode :
{ s: 's; step: 's -> 'a -> 'b * 's; reset : 'p -> 's -> 's } -> { s: 's; step: 's -> 'a -> 'b * 's; reset : 'p -> 's -> 's } ->
('p, 'a, 'b) dnode ('p, 'a, 'b) dnode
``` ```
An ODE solver is, in fact, a special case of a discrete node, which is An ODE solver is, in fact, a special case of a discrete node, which is
initialized with an initial value problem, takes as input a time horizon to initialized with an initial value problem, takes as input a time horizon to
reach, and returns as output an actual time reached and a dense solution. reach, and returns as output an actual time reached and a dense solution.
```ocaml ```ocaml
type ('y, 'yder) csolver = type ('y, 'yder) ivp = { init : 'y; fder : time -> 'y -> 'yder; stop : time }
(('y, 'yder) ivp, time, time * (time -> 'y)) dnode type ('y, 'yder) csolver = (('y, 'yder) ivp, time, time * (time -> 'y)) dnode
``` ```
Similarly, a zero-crossing solver is initialized with a zero-crossing problem, Similarly, a zero-crossing solver is initialized with a zero-crossing problem,
takes as input a horizon and dense function, and returns an actual time reached takes as input a horizon and dense function, and returns an actual time reached
and optional zero-crossing. and optional zero-crossing.
```ocaml ```ocaml
type ('y, 'zout) zc = { fzer : time -> 'y -> 'zout; yc : 'y }
type ('y, 'zin, 'zout) zsolver = type ('y, 'zin, 'zout) zsolver =
(('y, 'zout) zc, time * (time -> 'y), time * 'zin option) dnode (('y, 'zout) zc, time * (time -> 'y), time * 'zin option) dnode
``` ```
A complete solver is then a composition of these two: it is initialized with A complete solver is then a composition of these two: it is initialized with
both an `ivp` and `zc`, takes in a horizon, and returns an actual time reached, both an `ivp` and `zc`, takes in a horizon, and returns an actual time reached,
dense solution and optional zero-crossing. dense solution and optional zero-crossing.
```ocaml
#breakl(```ocaml
type ('y, 'yder, 'zin, 'zout) solver = type ('y, 'yder, 'zin, 'zout) solver =
(('y, 'yder) ivp * ('y, 'zout) zc, time, time * (time -> 'y) * 'zin option) dnode (('y, 'yder) ivp * ('y, 'zout) zc, time, time * (time -> 'y) * 'zin option) \
``` dnode
The construction of a full solver from an ODE and a zero-crossing solver is ```)
The construction of a full solver from an ODE and a zero-crossing solver is
available in the file `src/lib/hsim/solver.ml`. available in the file `src/lib/hsim/solver.ml`.
The simulation of a hybrid system with a solver can itself be seen as a discrete The simulation of a hybrid system with a solver can itself be seen as a discrete
node, operating on streams of functions defined on successive intervals. node, operating on streams of functions defined on successive intervals.
== Output format ```ocaml
type ('p, 'a, 'b) sim = ('p, 'a signal, 'b signal) dnode
```
In "lazy" mode, the output has format ```ocaml 'b signal = 'b value option```. == The simulation
In "greedy" mode, the output has format ```ocaml 'b value list```. === Modes and output format
Indeed, we cannot concatenate a discrete step (i.e. a function defined on a There are two simulation modes, depending on whether the solver provides a way
singleton interval $[t,t]$) with anything else, as it may be hidden by the to copy its internal state as part of its interface, or more precisely, whether
(inclusive) border of the other function: the dense solution returned by the solver's step function remains valid after we
modify the solver's state (through another integration step or a full reset).
In *_lazy_* mode, the output has format ```ocaml 'b signal = 'b value option```.
The input stream first provides a value for the simulation to run, and must then
provide as many `None` values as needed for the solver to finish integrating as
much as possible (that is, until it reaches a zero-crossing or the horizon
provided by the current input value). This mode is the one used when the solver
modifies its solution in place (for example, #smallcaps([Sundials])). Each step
of the simulation runs _one_ step of the solver, and then runs all subsystems as
much as they need to in order for them to finish using the output of the solver.
The solver output is then not needed anymore, and we can call the solver again
to progress with the simulation.
In *_greedy_* mode, the output has format ```ocaml 'b value list```. The input
stream provides a value for the simulation to run, and the simulation then calls
the solver as much as it needs to until it reaches the end of the possible
integration period (zero-crossing or horizon). We concatenate the list of
results (all functions returned from the solver steps are continuous at their
edges and defined on more than a single instant, and we can thus create a large,
piecewise-defined function from all the results up to this point). The solver
then performs as many discrete steps as needed (i.e. we cascade). These are kept
in the output list _as is_. We cannot concatenate a discrete step (i.e. a
function defined on a singleton interval $[t,t]$) with anything else, as it may
be hidden by the (inclusive) border of the other function. For instance,
#align(center, ```ocaml #align(center, ```ocaml
f = concat { start = 0.0; length = 1.0; u = fun t -> t } f = concat { start = 0.0; length = 1.0; u = fun t -> t }
{ start = 1.0; length = 0.0; u = fun t -> 2.0 } { start = 1.0; length = 0.0; u = fun t -> 2.0 }
```) ```)
#columns(2, [ #columns(2, [#{
#align(center, cetz.canvas({ align(center, cetz.canvas({
cetz-plot.plot.plot( import cetz-plot.plot: *
size: (1, 1), x-tick-step: 1, y-tick-step: 1, axis-style: "left", { plot(size: (1, 0.5), axis-style: "left", x-tick-step: 1, x-label: "time",
cetz-plot.plot.add(((0,0),(1,1))) y-tick-step: 1, y-label: none, {
cetz-plot.plot.add(((1,2),), mark: "o", mark-size: .03)})})) add(((0,0),(1,1)))
#colbreak() add(((1,2),), mark: "o", mark-size: .03)
#linebreak() })
}))
colbreak()
$ f(t) = cases( $ f(t) = cases(
t & "if" 0.0 <= t <= 1.0, t & "if" 0.0 <= t <= 1.0,
2.0 & "if" t = 1.0 + partial, 2.0 & "if" t = 1.0 + partial,
bot & "otherwise", bot & "otherwise",
) $ ) $
}])
Here, what is `f.u 1.0`? We cannot differentiate between the multiple
infinitesimal steps at a given real time `t` using only floating-point numbers
as input. Therefore, we keep the list of inputs separated, apart from the
adjacent continuous functions, which we can concatenate together. Until we reach
the horizon given by the input function, we keep alternating continuous and
discrete steps. The output is a list of functions, which we can interpret as a
stream, and the subsystems (the assertions, for instance) can then be called
repeatedly until they reach the same horizon as the parent system. This has the
advantage of not needing/having optional values in the input and output streams,
and the calling system does not need to guess how many steps the simulation will
need to reach the requested horizon.
=== Steps
A simulation step is either a discrete step or a continuous step.
Discrete steps call the `model.step` function, build a constant output function
defined on a singleton interval `[now, now]`, and end in one of four possible
cases:
- We need to cascade (perform another discrete step)
- We reached the horizon (reset the solver? idle mode?)
- We can still integrate, but there has been a state jump (reset solver)
- We can still integrate, and there has been no state jump (perform a continuous
step)
#columns(4, [#align(center, {
import cetz: *
import cetz-plot.plot: *
let p = palette.new(colors: (blue,blue, blue)).with(stroke: true, fill: true)
canvas({
plot(
size: (2,2), axis-style: "left", x-tick-step: none, x-label: none,
x-ticks: ((1,none),), x-grid: "both", y-tick-step: none, y-label: none,
y-min: 1, y-ticks: ((4, none),), y-grid: "both",
style: p, {
add(((0,2),(.3,2.1),(.7,2.9),(1,3)), line: "spline")
add(((1,4),), mark: "o", mark-size: .03)
add(((1,5),(1.3,4.7),(1.7,3.2),(2,3)), line: "spline")
})
})
colbreak(); canvas({
plot(
size: (2,2), axis-style: "left", x-tick-step: none, x-label: none,
x-ticks: ((3,none),), x-grid: "both", y-tick-step: none, y-label: none,
y-min: 1, y-ticks: ((3, none),), y-grid: "both", {
add(((0,2),(.5,3.2),(1,5.5),(2,4),(3,3)), line: "spline")
})
})
colbreak(); canvas({
plot(
size: (2,2), axis-style: "left", x-tick-step: none, x-label: none,
x-ticks: ((1, none),), x-grid: "both", y-tick-step: none, y-label: none,
y-min: 1, y-ticks: ((5, none),), y-grid: "both", {
add(((0,2),(0.5,2.3),(1,3)), line: "spline")
add(((1,5),(1.5,4.7),(2,4)), line: "spline")
})
})
colbreak(); canvas({
plot(
size: (2,2), axis-style: "left", x-tick-step: none, x-label: none,
x-ticks: ((2,none),), x-grid: "both", y-tick-step: none, y-label: none,
y-min: 1, y-ticks: ((4, none),), y-grid: "both", {
add(((0,2),(0.5,2.3),(1,3),(1.3,5),(2,4),(3,3)), line: "spline")
})
})
})])
Continuous steps call the `solver.step` function, parameterized by `model.fder`
and `model.fzer`, build an output function using `model.fout` defined on the
interval `[now, now + h]`, and end in one of three possible cases:
- We found a zero-crossing (perform a discrete step)
- We found no zero-crossing, but have reached the horizon (perform a discrete
step)
- We found no zero-crossing and have not reached the horizon (perform another
continuous step)
#columns(3, [#align(center, {
import cetz: *
import cetz-plot.plot: *
let p = palette.new(colors: (blue,blue, blue)).with(stroke: true, fill: true)
canvas({
plot(
size: (2,2), axis-style: "left", x-tick-step: none, x-label: none,
x-ticks: ((.7,none),), x-grid: "both", y-tick-step: none, y-label: none,
y-min: -2, y-ticks: ((0,$0$),), y-grid: "both", {
add(((0,-1),(.5,-.6),(.7,0),(1,2.5),(2,1),(3,.5)), line: "spline")
})
})
colbreak(); canvas({
plot(
size: (2,2), axis-style: "left", x-tick-step: none, x-label: none,
x-ticks: ((3,none),), x-grid: "both", y-tick-step: none, y-label: none,
y-min: 1, y-ticks: ((3, none),), y-grid: "both", {
add(((0,2),(.5,3.2),(1,5.5),(2,4),(3,3)), line: "spline")
})
})
colbreak(); canvas({
plot(
size: (2,2), axis-style: "left", x-tick-step: none, x-label: none,
x-ticks: ((2,none),), x-grid: "both", y-tick-step: none, y-label: none,
y-min: 1, y-ticks: ((4, none),), y-grid: "both", {
add(((0,2),(0.5,2.3),(1,3),(1.3,5),(2,4),(3,3)), line: "spline")
})
})
})])
=== Resets
==== Solver reset
When do we need to reset the solver itself?
If we receive a new input that is not continuous w.r.t. the previous input, we
should reset the solver before starting the integration. This implies that an
input value is tagged with information relating it to the previous input value.
If this is not the case, we can reset whenever we receive a new input value,
although this is less satisfactory.
If a discrete step results in a state jump, we should reset the solver before
continuing the integration.
==== Simulation reset
Two possible options for the simulation reset:
- ask for reset parameters for both the model and solver:
```ocaml
'p * (('y, 'yder) ivp * ('y, 'zin, 'zout) zc
```
- ask for only what is needed for the reset of the model, then build a new
initial value problem and zero-crossing problem using the new model state and
reset the solver using them. We can build `fder` and `fzer` quite easily. We
still need to build `init`, `yc` and `time`.
`init` and `yc` need to use `model.cget` in order to obtain a continuous state
element. However, the usual simulation process first starts by executing a
discrete step to initialize the model state. Does `model.cget` make sense
before this first step has been executed?
#figure(finite.automaton(
(D: (D: "cascade", C: "no cascade"),
C: (C: "no zero-crossing", D: "zero-crossing")),
initial: "D", final: (), layout: finite.layout.linear.with(spacing: 3)
)) <automaton>
Tagged values that indicate whether they are discontinuous w.r.t. the previous
input value do not help us here. Indeed, we could consider that the first
value in the stream is always considered discontinuous, which would force a
reset as we receive an input value. However, the initialization of the solver
still depends on the initialization of the model, which (might) need a
discrete step, and we are back to square one.
Additionally, what stop time do we use for the initial state?
If possible, this would have the advantage of making the reset function only
dependent on the model, since the solver is updated to match the model's
simulation requirements.
```ocaml
val sim : ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode ->
('y, 'yder, 'zin, 'zout) solver ->
('p, 'a, 'b) sim
```
== Mathematical model
#link("https://zelus.di.ens.fr/cc2015/fullpaper.pdf")[[CC'15]] defines the
simulation loop as follows:
#block(fill: rgb(230, 230, 230), stroke: black, inset: 5pt, [
The simulation loop of a hybrid system is a _synchronous function_ defining
four streams $t(n)$, $l x(n)$, $y(n)$ and $z(n)$ with $n in bb(N)$.
$t(n) in bb(R)$ is the increasing sequence of instants at which the solver
stops. $l x(n)$ is the value at time $t(n)$ of the _continuous state
variables_, that is, of all variables defined by their derivatives in the
original model. $y(n)$ is the value at time $t(n)$ of the _discrete state_.
$z(n)$ indicates any _zero-crossings_ at instant $t(n)$ on signals monitored
by the solver, that is, any signals that become equal to or pass through zero.
The synchronous function has two modes: the discrete mode ($D$) contains all
computations that may change the discrete state or that have side effects. The
continuous mode ($C$) is where ODEs are solved. The two modes alternate
according to the execution scheme summarized in @automaton.
In the _continuous mode_, the solver computes an approximation of the solution
of the ODEs and monitors a set of expressions for zero-crossings. The solver
is abstracted by a function $s o l v e(f)(g)$ parameterized by $f$ and $g$
where:
- $x'(tau)=f(y(n),tau,x(tau))$ defines the derivatives of continuous state
variables $x$ at instant $tau in bb(R)$;
- $u p z(tau)=g(y(n),tau,x(tau))$ defines the current values of a set of
zero-crossing signals $u p z$, indexed by $i in {1, ..., k}$.
The continuous mode $C$ computes four sequences $s$, $l x$, $z$ and $t$ such
that $ (l x, z, t, s)(n+1) = s o l v e(f)(g)(s,y,l x,t,s t e p)(n) $ where
- $s(n)$ is the internal state of the solver at instant $t(n) in bb(R)$.
Calling $s o l v e(f)(g)$ updates the state to $s(n+1)$.
- $x$ is an approximation of a solution of the ODE
$ x'(tau)=f(y(n),tau,x(tau)) $
It is parameterized by the current discrete state $y(n)$ and initialized at
instant $t(n)$ with the value of $l x(n)$, that is, $x(t(n)) = l x(n)$.
- $l x(n+1)$ is bounded by the horizon $t(n) + s t e p(n)$ that the solver has
been asked to reach, that is
$ t(n) <= t(n+1) <= t(n) + s t e p(n) $
- $z(n+1)$ signals any zero-crossings detected at time $t(n+1)$. An event
occurs with a transition to the discrete mode $D$ when horizon
$t(n) + s t e p(n)$ is reached, or when at least one of the zero-crossing
signals $u p z(i)$, for $i in {1, ..., k}$ crosses zero, which is indicated
by a true value for the corresponding boolean output $z(n+1)(i)$. $
e v e n t = z(n+1)(0) or ... or z(n+1)(k) or (t(n+1) = t(n) + s t e p(n)) $
If the solver raises an error (for example, a division by zero or an inability
to find a solution), we consider that the simulation fails.
All discrete changes occur in the _discrete mode_ $D$. It is entered when an
event is raised during integration. During a discrete phase, the function
$n e x t$ defines $y$, $l x$, $s t e p$, $e n c o r e$, $z$ and $t$:
$ (y, l x, s t e p, e n c o r e)(n+1) & = n e x t(y, l x, z, t)(n) \
z(n+1) & = f a l s e \
t(n+1) & = t(n) $
where
- $y(n+1)$ is the new discrete state; outside of mode $D$, $y(n+1) = y(n)$.
- $l x(n+1)$ is the new continuous state, which may be changed directly in the
discrete mode.
- $s t e p(n+1)$ is the new step size.
- $e n c o r e(n+1)$ is true if an additional discrete step must be performed.
Function $n e x t$ can decide to trigger instantaneously another discrete
event, causing an _event cascade_.
- $t(n)$ (the simulation time) is unchanged during discrete phases.
The initial values for $y(0)$, $l x(0)$ and $s(0)$ are given by an
initialization function $i n i t$. Finally, $s o l v e(f)(g)$ may decide to
reset its internal state if the continuous state changes. If
$i n i t \_ s o l v e$ initializes the solver state, we have
$ r e i n i t & = (l x(n+1) != l x(n)) \
s(n+1) & = "if" r e i n i t "then" i n i t \_ s o l v e(l x(n+1),s(n))
"else" s(n) $
Taken together, the definitions from both modes give a synchronous
interpretation of the simulation loop as a stream function that computes the
sequences $l x$, $y$ and $t$ at instant $n+1$ according to their values at
instant $n$ and an internal state. Writing $s o l v e(f)(g)$ abstracts from
the actual choice of integration method and zero-crossing detection algorithm.
A more detailed description of $s o l v e(f)(g)$ would be possible (for
example, an automaton with two states: one that integrates, and one that
detects zero-crossings) but with no influence on the code generation problem
which must be independent of such simulation details.
Given a program written in a high-level language, we must produce the
functions $i n i t$, $f$, $g$, and $n e x t$. In practice, they are
implemented in an imperative language like C. Code generation for hybrid
models has much in common with code generation for synchronous languages.
]) ])
How do we represent infinitesimal steps in floating-point numbers ?
A possible representation of the above could be
#align(center, ```ocaml
f = [{ start = 0.0; length = 1.0; u = fun t -> t };
{ start = 1.0; length = 0.0; u = fun _ -> 2.0 }]
```)
but returning a list is problematic, however, as it cannot be passed easily as
input to the underlying models. Unless the input can itself be considered as a
list, in which case the simulation simply operates over the whole list,
resetting the solver as needed. This has the advantage of not having to wait
several steps, providing nothing, until the solver is done integrating the
input, before we can provide it with another input.

144
src/lib/hsim/sim.ml
View file

@ -6,19 +6,28 @@ open State
let offset (input : 'a value) (now : time) : time -> 'a = let offset (input : 'a value) (now : time) : time -> 'a =
fun t -> input.u ((now -. input.start) +. t) fun t -> input.u ((now -. input.start) +. t)
let rec compose = function
| [] -> assert false
| [f] -> f
| { start=sl; u=ul; _ } :: l ->
let { start=sr; length=lr; u=ur } = compose l in
let length = sr +. lr -. sl in
{ start=sl; length; u=fun t -> if t <= sr then ur t else ul t }
module LazySim (S : SimState) = module LazySim (S : SimState) =
struct struct
(* TODO: figure out where we initialize the solvers; the initialization (* TODO: figure out where we initialize the solvers; the initialization and
function already supposes an initialized solver state, but could we reset functions already suppose an initialized solver state, but
parameterize [LazySim] with a solver state module that provides its could we parameterize simulation with a solver state module that
own initialization function ? *) provides its own initialization function ? *)
let sim let sim
(HNode model : ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode) (HNode model : ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode)
(DNode solver : ('y, 'yder, 'zin, 'zout) solver) (DNode solver : ('y, 'yder, 'zin, 'zout) solver)
: ('p, 'a, 'b) sim : ('p, 'a, 'b) lazy_sim
= let s = S.init ~mstate:model.state ~sstate:solver.state in = let state = S.init ~mstate:model.state ~sstate:solver.state in
let step state input = let step state input =
let mstate = S.mstate state and sstate = S.sstate state let mstate = S.mstate state and sstate = S.sstate state
and status = S.status state in and status = S.status state in
@ -39,8 +48,9 @@ module LazySim (S : SimState) =
let h = model.horizon mstate in let h = model.horizon mstate in
if h <= 0.0 then S.update ~mstate state if h <= 0.0 then S.update ~mstate state
else if now >= stop then else if now >= stop then
(* TODO: equivalent of [s] (initial state of model and (* TODO: Build an initial state with the initial states for
solvers) ? *) both the solver and model - an equivalent of [s] in
the original version. *)
raise Common.Utils.TODO raise Common.Utils.TODO
else if model.jump mstate then else if model.jump mstate then
let y = model.cget mstate in let y = model.cget mstate in
@ -75,77 +85,97 @@ module LazySim (S : SimState) =
let mstate = model.zset mstate z in let mstate = model.zset mstate z in
S.update ~status ~mstate ~sstate state in S.update ~status ~mstate ~sstate state in
Some out, state in Some out, state in
let reset m s =
let mstate = model.reset m (S.mstate s) in let reset p s =
(* TODO: does [model.cget mstate] make sense before the first discrete
step - can we use [mstate] to reinitialize [sstate] ? *)
let mstate = model.reset p (S.mstate s) in
let y = model.cget mstate in let y = model.cget mstate in
(* TODO: what initial stop time do we use ? *)
let stop = raise Common.Utils.TODO in let stop = raise Common.Utils.TODO in
let ivp = { fder = model.fder mstate; stop; init = y } in let ivp = { fder = model.fder mstate; stop; init = y } in
let zc = { fzer = model.fzer mstate; yc = y } in let zc = { fzer = model.fzer mstate; yc = y } in
let sstate = solver.reset (ivp, zc) (S.sstate s) in let sstate = solver.reset (ivp, zc) (S.sstate s) in
let status = S.idle (S.status s) in let status = S.idle (S.status s) in
S.update ~status ~mstate ~sstate s in S.update ~status ~mstate ~sstate s in
DNode { state = s; step; reset }
DNode { state; step; reset }
end end
module GreedySim (S : SimState) = module GreedySim (S : SimState) =
struct struct
(* TODO: greedy simulation *) (* TODO: greedy simulation: call the solvers and the subsystems as often as
needed until we reach the horizon. *)
let sim let sim
(HNode model : ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode) (HNode model : ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode)
(DNode solver : ('y, 'yder, 'zin, 'zout) solver) (DNode solver : ('y, 'yder, 'zin, 'zout) solver)
: ('p, 'a, 'b) sim : ('p, 'a, 'b) greedy_sim
= let state = S.init ~mstate:model.state ~sstate:solver.state in = let state = S.init ~mstate:model.state ~sstate:solver.state in
let rec step state input = let rec step state input =
let status = S.status state and mstate = S.mstate state let status = S.status state and mstate = S.mstate state
and sstate = S.sstate state in and sstate = S.sstate state in
match input, S.is_running state with if not (S.is_running state) then
| Some input, _ -> let mode = Discrete and now = 0.0 and stop = input.length in
let mode = Discrete and now = 0.0 and stop = input.length in let status = S.running ~mode ~input ~now ~stop (S.status state) in
let status = S.running ~mode ~input ~now ~stop (S.status state) in let state = S.update ~status state in
let state = S.update ~status state in step state input
step state None else let now = S.now state and stop = S.stop state in
| None, false -> None, state match S.mode state with
| None, true -> | Discrete ->
let input = S.input state and now = S.now state let o, mstate = model.step mstate (input.u now) in
and stop = S.stop state in let h = model.horizon mstate in
match S.mode state with let rest, state =
| Discrete -> if h <= 0.0 then step (S.update ~mstate state) input
let o, mstate = model.step mstate (input.u now) in else if now >= stop then [], state
let state = else if model.jump mstate then
let h = model.horizon mstate in let y = model.cget mstate in
if h <= 0.0 then S.update ~mstate state let fder t = model.fder mstate (offset input now t) in
else if now >= stop then raise Common.Utils.TODO let fzer t = model.fzer mstate (offset input now t) in
else if model.jump mstate then let ivp = { fder; stop = stop -. now; init = y } in
let y = model.cget mstate in let zc = { yc = y; fzer } in
let fder t = model.fder mstate (offset input now t) in let sstate = solver.reset (ivp, zc) sstate in
let fzer t = model.fzer mstate (offset input now t) in let status = S.running ~mode:Continuous status in
let ivp = { fder; stop = stop -. now; init = y } in step (S.update ~status ~mstate ~sstate state) input
let zc = { yc = y; fzer } in else
let sstate = solver.reset (ivp, zc) sstate in let status = S.running ~mode:Continuous status in
let status = S.running ~mode:Continuous status in step (S.update ~status state) input in
S.update ~status ~mstate ~sstate state let start = input.start +. now in
else { start; length = 0.0; u = fun _ -> o }::rest, state
let status = S.running ~mode:Continuous status in | Continuous ->
S.update ~status state in let (h, f, z), _sstate = solver.step sstate stop in
let start = input.start +. now in let mstate = model.cset mstate (f h) in
Some { start; length = 0.0; u = fun _ -> o }, state let h' = input.start +. h in
| Continuous -> let fout t =
let (h, f, z), sstate = solver.step sstate stop in model.fout mstate (input.u (now +. t)) (f (now +. t)) in
let mstate = model.cset mstate (f h) in let out =
let h' = input.start +. h in { start = input.start +. now; length = h -. now; u = fout } in
let fout t = match z with
model.fout mstate (input.u (now +. t)) (f (now +. t)) in | None ->
let out = if h >= stop then
{ start = input.start +. now; length = h -. now; u = fout } in let status = S.running ~mode:Discrete ~now:h' status in
match z with let rest, state =
| None -> step (S.update ~status ~mstate ~sstate state) input in
let status = out::rest, state
if h >= stop then S.running ~mode:Discrete ~now:h' status else
in let status = S.running ~now:h' status in
let rest, state =
step (S.update ~status ~mstate ~sstate state) input in
(match rest with
| [] -> [out], state
| f::rest -> compose [out;f] :: rest, state)
| Some z ->
let status = S.running ~mode:Discrete ~now:h' status in
let mstate = model.zset mstate z in
let rest, state =
step (S.update ~status ~mstate ~sstate state) input in
out::rest, state in
let reset = assert false in let reset = assert false in
DNode { state; step; reset } DNode { state; step; reset }
end end

24
src/lib/hsim/state.ml
View file

@ -1,20 +1,6 @@
open Types open Types
(* Model *)
module type ModelState =
sig
type state
end
(* Solvers *)
module type SolverState =
sig
type state
end
(* Simulation *) (* Simulation *)
(** Step mode: discrete or continuous *) (** Step mode: discrete or continuous *)
@ -52,7 +38,7 @@ module type SimState =
val idle : 'a status -> 'a status val idle : 'a status -> 'a status
(** Update the status to running. *) (** Update the status to running. *)
val running : val running :
?mode:mode -> ?input:'a value -> ?mode:mode -> ?input:'a value ->
?now:time -> ?stop:time -> 'a status -> 'a status ?now:time -> ?stop:time -> 'a status -> 'a status
@ -81,7 +67,7 @@ module FunctionalSimState : SimState =
(** Simulation status: (** Simulation status:
- [Idle]: Waiting for input, no activity; - [Idle]: Waiting for input, no activity;
- [Running]: Currently integrating: step [mode], current [input], [now] - [Running]: Currently integrating: step [mode], current [input], [now]
timestamp, and [stop] time. *) timestamp, and [stop] time. *)
type 'a status = type 'a status =
| Idle : 'a status | Idle : 'a status
@ -92,7 +78,7 @@ module FunctionalSimState : SimState =
stop : time; (** How long until we stop. *) stop : time; (** How long until we stop. *)
} -> 'a status } -> 'a status
(** Internal state of the simulation node: model state, solver state and (** Internal state of the simulation node: model state, solver state and
current simulation status. *) current simulation status. *)
type ('a, 'ms, 'ss) state = type ('a, 'ms, 'ss) state =
{ status : 'a status; (** Current simulation status. *) { status : 'a status; (** Current simulation status. *)
@ -137,7 +123,7 @@ module FunctionalSimState : SimState =
match s.status with Running r -> r.input | Idle -> raise Not_running match s.status with Running r -> r.input | Idle -> raise Not_running
let now s = let now s =
match s.status with Running r -> r.now | Idle -> raise Not_running match s.status with Running r -> r.now | Idle -> raise Not_running
let stop s = let stop s =
match s.status with Running r -> r.stop | Idle -> raise Not_running match s.status with Running r -> r.stop | Idle -> raise Not_running
let init ~mstate ~sstate = { status = Idle; mstate; sstate } let init ~mstate ~sstate = { status = Idle; mstate; sstate }
@ -198,7 +184,7 @@ module InPlaceSimState : SimState =
match s.status with Running r -> r.input | Idle -> raise Not_running match s.status with Running r -> r.input | Idle -> raise Not_running
let now s = let now s =
match s.status with Running r -> r.now | Idle -> raise Not_running match s.status with Running r -> r.now | Idle -> raise Not_running
let stop s = let stop s =
match s.status with Running r -> r.stop | Idle -> raise Not_running match s.status with Running r -> r.stop | Idle -> raise Not_running
let init ~mstate ~sstate = { status = Idle; mstate; sstate } let init ~mstate ~sstate = { status = Idle; mstate; sstate }

7
src/lib/hsim/types.mli
View file

@ -81,5 +81,10 @@ type ('y, 'yder, 'zin, 'zout) solver =
(** The simulation of a hybrid system is a synchronous function on streams of (** The simulation of a hybrid system is a synchronous function on streams of
functions. *) functions. *)
type ('p, 'a, 'b) sim = type ('p, 'a, 'b) lazy_sim =
('p, 'a signal, 'b signal) dnode ('p, 'a signal, 'b signal) dnode
(** Greedy simulation takes in an input and computes as many solver and
subsystem steps as needed to reach the input's horizon. *)
type ('p, 'a, 'b) greedy_sim =
('p, 'a value, 'b value list) dnode