Skip to content

simulator

Hybrid system simulation utilities.

HybridSimulationError

Bases: RuntimeError

Base class for hybrid simulation runtime failures.

InvalidEventSurfaceValueError

Bases: HybridSimulationError

Raised when an event surface returns NaN.

SurfaceEntryError

Bases: HybridSimulationError

Raised when an ERROR surface is zero on location entry.

AmbiguousTransitionError

Bases: HybridSimulationError

Raised when multiple TRIGGER surfaces are zero on location entry.

SimulationProgressError

Bases: HybridSimulationError

Raised when a solver event does not advance physical time.

ensure_state(state)

Validate and coerce a state vector into a 1D array.

Source code in src/flowcean/hybrid/simulator.py
46
47
48
49
50
51
52
def ensure_state(state: Iterable[float]) -> State:
    """Validate and coerce a state vector into a 1D array."""
    array = np.asarray(state, dtype=float)
    if array.ndim != 1:
        message = "State must be a 1D array."
        raise ValueError(message)
    return array

simulate(system, t_span, x0=None, location0=None, *, input_stream=None, capture_inputs=None, capture_derivatives=False, max_jumps=256, rtol=1e-07, atol=1e-09, max_step=None, dense_output=False, sample_times=None, sample_dt=None)

Simulate a hybrid system and return a trace.

Parameters:

Name Type Description Default
system HybridSystem

Hybrid system to simulate.

required
t_span tuple[float, float]

Start and end time for integration.

required
x0 Iterable[float] | None

Optional initial state override.

None
location0 Location | None

Optional initial location override.

None
input_stream InputStream | None

Optional input stream accessor for callbacks.

None
capture_inputs bool | None

Input capture mode. If None, capture iff an input stream is provided.

None
capture_derivatives bool

Whether to re-evaluate Location.dynamics.flow on the returned trace grid and store the sampled derivatives in Trace.dx. This assumes pure flow callbacks under repeated evaluation. Scalar derivative returns are accepted only for single-state systems.

False
max_jumps int

Maximum number of transitions allowed.

256
rtol float

Relative tolerance for the solver.

1e-07
atol float

Absolute tolerance for the solver.

1e-09
max_step float | None

Optional maximum step size.

None
dense_output bool

Whether to build a continuous solution per segment.

False
sample_times Iterable[float] | None

Monotone time grid to sample from the dense solution.

None
sample_dt float | None

Fixed sampling interval to generate a time grid.

None

Returns:

Name Type Description
Trace Trace

The simulation trace with location labels and events.

Source code in src/flowcean/hybrid/simulator.py
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
def simulate(
    system: HybridSystem,
    t_span: tuple[float, float],
    x0: Iterable[float] | None = None,
    location0: Location | None = None,
    *,
    input_stream: InputStream | None = None,
    capture_inputs: bool | None = None,
    capture_derivatives: bool = False,
    max_jumps: int = 256,
    rtol: float = 1e-7,
    atol: float = 1e-9,
    max_step: float | None = None,
    dense_output: bool = False,
    sample_times: Iterable[float] | None = None,
    sample_dt: float | None = None,
) -> Trace:
    """Simulate a hybrid system and return a trace.

    Args:
        system: Hybrid system to simulate.
        t_span: Start and end time for integration.
        x0: Optional initial state override.
        location0: Optional initial location override.
        input_stream: Optional input stream accessor for callbacks.
        capture_inputs: Input capture mode. If ``None``, capture iff an
            input stream is provided.
        capture_derivatives: Whether to re-evaluate ``Location.dynamics.flow``
            on the returned trace grid and store the sampled derivatives in
            ``Trace.dx``. This assumes pure flow callbacks under repeated
            evaluation. Scalar derivative returns are accepted only for
            single-state systems.
        max_jumps: Maximum number of transitions allowed.
        rtol: Relative tolerance for the solver.
        atol: Absolute tolerance for the solver.
        max_step: Optional maximum step size.
        dense_output: Whether to build a continuous solution per segment.
        sample_times: Monotone time grid to sample from the dense solution.
        sample_dt: Fixed sampling interval to generate a time grid.

    Returns:
        Trace: The simulation trace with location labels and events.
    """
    start_location = (
        system.initial_location if location0 is None else location0
    )
    if not isinstance(start_location, Location):
        message = "location0 must be a Location."
        raise TypeError(message)
    location_ids = {id(location) for location in system.locations}
    if id(start_location) not in location_ids:
        message = "location0 must be included in system.locations."
        raise ValueError(message)

    state = ensure_state(x0 if x0 is not None else system.initial_state)
    location = start_location
    effective_input_stream = input_stream or _missing_input_stream
    should_capture = _resolve_capture_inputs(
        capture_inputs=capture_inputs,
        input_stream=input_stream,
    )

    t_segments: list[np.ndarray] = []
    x_segments: list[np.ndarray] = []
    location_segments: list[np.ndarray] = []
    sol_segments: list[Callable[[np.ndarray], np.ndarray] | None] = []
    events: list[Event] = []
    boundaries: list[_Boundary] = []

    t_current = float(t_span[0])
    t_final = float(t_span[1])
    jumps = 0

    sample_grid = _prepare_sample_times(t_span, sample_times, sample_dt)
    needs_dense = dense_output or sample_grid is not None

    initial_entry = _settle_location_entries(
        system,
        location,
        state,
        t_current,
        effective_input_stream,
        first_microstep=0,
        jumps=jumps,
        max_jumps=max_jumps,
    )
    state = initial_entry.state
    location = initial_entry.location
    jumps = initial_entry.jumps
    events.extend(initial_entry.events)
    if initial_entry.events:
        boundaries.append(_Boundary(t_current, state.copy(), location))

    while t_current < t_final:
        transitions = system.transitions_from(location)
        event_fns = _build_event_functions(
            transitions,
            system.parameters,
            location.parameters,
            effective_input_stream,
        )
        segment_start = t_current

        solve_kwargs = {
            "fun": _wrap_flow(
                location,
                system.parameters,
                effective_input_stream,
            ),
            "t_span": (segment_start, t_final),
            "y0": state,
            "events": event_fns or None,
            "rtol": rtol,
            "atol": atol,
            "dense_output": needs_dense,
        }
        if max_step is not None:
            solve_kwargs["max_step"] = max_step

        result = solve_ivp(**solve_kwargs)
        if not result.success:
            message = f"ODE integration failed: {result.message}"
            raise HybridSimulationError(message)

        t_segments.append(result.t)
        x_segments.append(result.y.T)
        location_segments.append(
            np.full(result.t.shape, location, dtype=object),
        )
        sol_segments.append(result.sol)

        if not result.t_events or all(
            len(event_list) == 0 for event_list in result.t_events
        ):
            break

        triggered_index, event_time, event_state = _first_event(
            result.t_events,
            result.y_events,
        )
        if event_time <= segment_start:
            progress_context = (
                f"segment start={segment_start!r}, event time={event_time!r}"
            )
            message = (
                f"An event did not advance physical time ({progress_context})."
            )
            error = SimulationProgressError(message)
            error.add_note(
                "This can result from stateful callbacks, discontinuous event "
                "surfaces, or insufficient floating-point time resolution. "
                "Use deterministic callbacks and continuous event surfaces.",
            )
            raise error

        transition = transitions[triggered_index]
        jumps = _increment_jumps(jumps, max_jumps)
        state, event = _apply_transition(
            transition,
            event_time,
            event_state,
            system.parameters,
            effective_input_stream,
            microstep=0,
        )
        events.append(event)
        location = transition.target

        target_entry = _settle_location_entries(
            system,
            location,
            state,
            event_time,
            effective_input_stream,
            first_microstep=1,
            jumps=jumps,
            max_jumps=max_jumps,
        )
        state = target_entry.state
        location = target_entry.location
        jumps = target_entry.jumps
        events.extend(target_entry.events)
        boundaries.append(_Boundary(event_time, state.copy(), location))
        t_current = event_time

    if sample_grid is None:
        t_all = _concat_segments(t_segments)
        x_all = _concat_segments(x_segments)
        location_objects = _concat_segments(location_segments)
        _apply_boundaries(
            t_all,
            x_all,
            location_objects,
            boundaries,
        )
        unique_times = _unique_time_mask(t_all)
        t_all = t_all[unique_times]
        x_all = x_all[unique_times]
        location_objects = location_objects[unique_times]
        location_all = _location_labels(location_objects)
        u_all = None
        dx_all = None
        if should_capture:
            if input_stream is None:
                message = "Internal error: expected input_stream for capture."
                raise RuntimeError(message)
            u_all = _capture_inputs(t_all, input_stream)
        if capture_derivatives:
            dx_all = _capture_derivatives(
                system=system,
                times=t_all,
                states=x_all,
                locations=location_objects,
                input_stream=effective_input_stream,
            )
        return Trace(
            t=t_all,
            x=x_all,
            location=location_all,
            events=tuple(events),
            u=u_all,
            dx=dx_all,
        )

    rolled = _rollout_segments(
        sample_grid,
        t_segments,
        x_segments,
        location_segments,
        sol_segments,
        boundaries,
    )
    u_all = None
    dx_all = None
    if should_capture:
        if input_stream is None:
            message = "Internal error: expected input_stream for capture."
            raise RuntimeError(message)
        u_all = _capture_inputs(rolled.t, input_stream)
    if capture_derivatives:
        dx_all = _capture_derivatives(
            system=system,
            times=rolled.eval_t,
            states=rolled.x,
            locations=rolled.location,
            input_stream=effective_input_stream,
        )
    location_all = _location_labels(rolled.location)
    return Trace(
        t=rolled.t,
        x=rolled.x,
        location=location_all,
        events=tuple(events),
        u=u_all,
        dx=dx_all,
    )

generate_traces(system, t_span, initial_states, *, input_stream=None, capture_inputs=None, capture_derivatives=False, max_jumps=256, rtol=1e-07, atol=1e-09, max_step=None, dense_output=False, sample_times=None, sample_dt=None)

Simulate a batch of traces for a set of initial states.

The input stream and capture semantics match :func:simulate, including the requirement that capture_derivatives=True assumes pure flow callbacks under repeated evaluation on the returned trace grid. Scalar derivative returns are accepted only for single-state systems.

Source code in src/flowcean/hybrid/simulator.py
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
def generate_traces(
    system: HybridSystem,
    t_span: tuple[float, float],
    initial_states: Iterable[Iterable[float]],
    *,
    input_stream: InputStream | None = None,
    capture_inputs: bool | None = None,
    capture_derivatives: bool = False,
    max_jumps: int = 256,
    rtol: float = 1e-7,
    atol: float = 1e-9,
    max_step: float | None = None,
    dense_output: bool = False,
    sample_times: Iterable[float] | None = None,
    sample_dt: float | None = None,
) -> list[Trace]:
    """Simulate a batch of traces for a set of initial states.

    The input stream and capture semantics match :func:`simulate`, including
    the requirement that ``capture_derivatives=True`` assumes pure flow
    callbacks under repeated evaluation on the returned trace grid. Scalar
    derivative returns are accepted only for single-state systems.
    """
    return [
        simulate(
            system,
            t_span,
            x0=state,
            input_stream=input_stream,
            capture_inputs=capture_inputs,
            capture_derivatives=capture_derivatives,
            max_jumps=max_jumps,
            rtol=rtol,
            atol=atol,
            max_step=max_step,
            dense_output=dense_output,
            sample_times=sample_times,
            sample_dt=sample_dt,
        )
        for state in initial_states
    ]