From b4a29bbb97a91c4258085576f163f46de64ca5bf Mon Sep 17 00:00:00 2001 From: Henri Saudubray Date: Thu, 17 Apr 2025 18:19:08 +0200 Subject: [PATCH] feat (doc): some notes --- doc/notes.typ | 97 +++++++++++++++++++++++++++++++++++++++++++++ src/lib/hsim/sim.ml | 52 ++++++++++++++++++++---- 2 files changed, 142 insertions(+), 7 deletions(-) create mode 100644 doc/notes.typ diff --git a/doc/notes.typ b/doc/notes.typ new file mode 100644 index 0000000..2e48a1c --- /dev/null +++ b/doc/notes.typ @@ -0,0 +1,97 @@ +#import "@preview/cetz:0.3.2" +#import "@preview/cetz-plot:0.1.1" + +#set page(paper: "a4") +#set par(justify: true) + += Notes on hybrid system simulation + +== The problem + +The initial problem is _transparent assertions_. We want to simulate assertions +with their own solvers, because simulating assertions with the same solver as +the parent system changes the final results, when assertions are supposed to be +transparent. + +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 +a parameter, and *recursive*, in that the simulation of an assertion is itself +the simulation of a hybrid system with assertions. + +The main idea is as follows: a solver can be seen as a synchronous function (an +inner state, a step function, and a reset function). +```ocaml +type ('p, 'a, 'b) dnode = DNode : + { s: 's; step: 's -> 'a -> 'b * 's; reset : 'p -> 's -> 's } -> + ('p, 'a, 'b) dnode +``` +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 +reach, and returns as output an actual time reached and a dense solution. +```ocaml +type ('y, 'yder) csolver = + (('y, 'yder) ivp, time, time * (time -> 'y)) dnode +``` +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 +and optional zero-crossing. +```ocaml +type ('y, 'zin, 'zout) zsolver = + (('y, 'zout) zc, time * (time -> 'y), time * 'zin option) dnode +``` +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, +dense solution and optional zero-crossing. +```ocaml +type ('y, 'yder, 'zin, 'zout) solver = + (('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 +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 +node, operating on streams of functions defined on successive intervals. + +== Output format + +In "lazy" mode, the output has format ```ocaml 'b signal = 'b value option```. + +In "greedy" mode, the output has format ```ocaml 'b value list```. + +Indeed, 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: + +#align(center, ```ocaml +f = concat { start = 0.0; length = 1.0; u = fun t -> t } + { start = 1.0; length = 0.0; u = fun t -> 2.0 } +```) + +#columns(2, [ + #align(center, cetz.canvas({ + cetz-plot.plot.plot( + size: (1, 1), x-tick-step: 1, y-tick-step: 1, axis-style: "left", { + cetz-plot.plot.add(((0,0),(1,1))) + cetz-plot.plot.add(((1,2),), mark: "o", mark-size: .03)})})) + #colbreak() + #linebreak() + $ f(t) = cases( + t & "if" 0.0 <= t <= 1.0, + 2.0 & "if" t = 1.0 + partial, + bot & "otherwise", + ) $ +]) + +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. diff --git a/src/lib/hsim/sim.ml b/src/lib/hsim/sim.ml index 40bc906..c7222c7 100644 --- a/src/lib/hsim/sim.ml +++ b/src/lib/hsim/sim.ml @@ -60,6 +60,10 @@ module LazySim (S : SimState) = let (h, f, z), sstate = solver.step sstate stop in let mstate = model.cset mstate (f h) in let h' = input.start +. h in + let fout t = + model.fout mstate (input.u (now +. t)) (f (now +. t)) in + let out = + { start = input.start +. now; length = h -. now; u = fout } in let state = match z with | None -> let status = @@ -70,11 +74,6 @@ module LazySim (S : SimState) = let status = S.running ~mode:Discrete ~now:h' status in let mstate = model.zset mstate z in S.update ~status ~mstate ~sstate state in - let fout t = - model.fout mstate (input.u (now +. t)) (f (now +. t)) in - let out = - { start = input.start +. now; length = h -. now; u = fout } - in Some out, state in let reset m s = let mstate = model.reset m (S.mstate s) in @@ -92,20 +91,59 @@ module LazySim (S : SimState) = module GreedySim (S : SimState) = struct + (* TODO: greedy simulation *) + let sim (HNode model : ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode) (DNode solver : ('y, 'yder, 'zin, 'zout) solver) : ('p, 'a, 'b) sim = let state = S.init ~mstate:model.state ~sstate:solver.state in let rec step state input = + let status = S.status state and mstate = S.mstate state + and sstate = S.sstate state in match input, S.is_running state with | Some input, _ -> 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 state = S.update ~status state in - None, state + step state None | None, false -> None, state - | None, true -> assert false + | None, true -> + let input = S.input state and now = S.now state + and stop = S.stop state in + match S.mode state with + | Discrete -> + let o, mstate = model.step mstate (input.u now) in + let state = + let h = model.horizon mstate in + if h <= 0.0 then S.update ~mstate state + else if now >= stop then raise Common.Utils.TODO + else if model.jump mstate then + let y = model.cget mstate in + let fder t = model.fder mstate (offset input now t) in + let fzer t = model.fzer mstate (offset input now t) in + let ivp = { fder; stop = stop -. now; init = y } in + let zc = { yc = y; fzer } in + let sstate = solver.reset (ivp, zc) sstate in + let status = S.running ~mode:Continuous status in + S.update ~status ~mstate ~sstate state + else + let status = S.running ~mode:Continuous status in + S.update ~status state in + let start = input.start +. now in + Some { start; length = 0.0; u = fun _ -> o }, state + | Continuous -> + let (h, f, z), sstate = solver.step sstate stop in + let mstate = model.cset mstate (f h) in + let h' = input.start +. h in + let fout t = + model.fout mstate (input.u (now +. t)) (f (now +. t)) in + let out = + { start = input.start +. now; length = h -. now; u = fout } in + match z with + | None -> + let status = + if h >= stop then S.running ~mode:Discrete ~now:h' status in let reset = assert false in DNode { state; step; reset }