Simulation¶

class pathsim.simulation.Simulation(blocks=None, connections=None, events=None, dt=0.01, dt_min=1e-16, dt_max=None, Solver=<class 'pathsim.solvers.ssprk22.SSPRK22'>, tolerance_fpi=1e-10, iterations_max=200, log=True, diagnostics=False, **solver_kwargs)[source]¶

Bases: object

Class that performs transient analysis of the dynamical system, defined by the blocks and connecions. It manages all the blocks and connections and the timestep update.

The global system equation is evaluated by fixed point iteration, so the information from each timestep gets distributed within the entire system and is available for all blocks at all times.

The minimum number of fixed-point iterations ‘iterations_min’ is set to ‘None’ by default and then the length of the longest internal signal path (with passthrough) is used as the estimate for minimum number of iterations needed for the information to reach all instant time blocks in each timestep. Dont change this unless you know that the actual path is shorter or something similar that prohibits instant time information flow.

Convergence check for the fixed-point iteration loop with ‘tolerance_fpi’ is based on max absolute error (max-norm) to previous iteration and should not be touched.

Multiple numerical integrators are implemented in the ‘pathsim.solvers’ module. The default solver is a fixed timestep 2nd order Strong Stability Preserving Runge Kutta (SSPRK22) method which is quite fast and has ok accuracy, especially if you are forced to take small steps to cover the behaviour of forcing functions. Adaptive timestepping and implicit integrators are also available.

Manages an event handling system based on zero crossing detection. Uses ‘Event’ objects to monitor solver states of stateful blocks and applys transformations on the state in case an event is detected.

Example

This is how to setup a simple system simulation using the ‘Simulation’ class:

import numpy as np

from pathsim import Simulation, Connection
from pathsim.blocks import Source, Integrator, Scope

src = Source(lambda t: np.cos(2*np.pi*t))
itg = Integrator()
sco = Scope(labels=["source", "integrator"])

sim = Simulation(
    blocks=[src, itg, sco],
    connections=[
        Connection(src[0], itg[0], sco[0]),
        Connection(itg[0], sco[1])
        ],
    dt=0.01
    )

sim.run(4)
sim.plot()

Parameters:

blocks (list[Block]) – blocks that define the system
connections (list[Connection]) – connections that connect the blocks
events (list[Event]) – list of event trackers (zero crossing detection, schedule, etc.)
dt (float) – transient simulation timestep in time units, default see ´SIM_TIMESTEP´ in ´_constants.py´
dt_min (float) – lower bound for transient simulation timestep, default see ´SIM_TIMESTEP_MIN´ in ´_constants.py´
dt_max (float) – upper bound for transient simulation timestep, default see ´SIM_TIMESTEP_MAX´ in ´_constants.py´
Solver (Solver) – ODE solver class for numerical integration from ´pathsim.solvers´, default is ´pathsim.solvers.ssprk22.SSPRK22´ (2nd order expl. Runge Kutta)
tolerance_fpi (float) – absolute tolerance for convergence of algebraic loops and internal optimizers of implicit ODE solvers, default see ´SIM_TOLERANCE_FPI´ in ´_constants.py´
iterations_max (int) – maximum allowed number of iterations for implicit ODE solver optimizers and algebraic loop solver, default see ´SIM_ITERATIONS_MAX´ in ´_constants.py´
log (bool | string) – flag to enable logging, default see ´LOG_ENABLE´ in ´_constants.py´ (alternatively a path to a log file can be specified)
solver_kwargs (dict) – additional parameters for numerical solvers such as absolute (´tolerance_lte_abs´) and relative (´tolerance_lte_rel´) tolerance, defaults are defined in ´_constants.py´

time¶

global simulation time, starting at ´0.0´

Type:: float

graph¶

internal graph representation for fast system funcion evluations using DAG with algebraic depths

Type:: Graph

boosters¶

list of boosters (fixed point accelerators) that wrap algebraic loop closing connections assembled from the system graph

Type:: None | list[ConnectionBooster]

engine¶

global integrator (ODE solver) instance serving as a dummy to get attributes and access to intermediate evaluation stages

Type:: Solver

logger¶

global simulation logger

Type:: logging.Logger

_blocks_dyn¶

blocks with internal ´Solver´ instances (stateful)

Type:: list[Block]

_blocks_evt¶

blocks with internal events (discrete time, eventful)

Type:: list[Block]

_active¶

flag for setting the simulation as active, used for interrupts

Type:: bool

property size¶

Get size information of the simulation, such as total number of blocks and dynamic states, with recursive retrieval from subsystems

Returns:: sizes – size of simulation (number of blocks) and number of internal states (from internal engines)
Return type:: tuple[int]

plot(*args, **kwargs)[source]¶

Plot the simulation results by calling all the blocks that have visualization capabilities such as the ‘Scope’ and ‘Spectrum’.

This is a quality of life method. Blocks can be visualized individually due to the object oriented nature, but it might be nice to just call the plot metho globally and look at all the results at once. Also works for models loaded from an external file.

Parameters:

args (tuple) – args for the plot methods
kwargs (dict) – kwargs for the plot method

save_checkpoint(path, recordings=True)[source]¶

Save simulation state to checkpoint files (JSON + NPZ).

Creates two files: {path}.json (structure/metadata) and {path}.npz (numerical data). Blocks and events are keyed by type and insertion order for deterministic cross-instance matching.

Parameters:

path (str) – base path without extension
recordings (bool) – include scope/spectrum recording data (default: True)

load_checkpoint(path)[source]¶

Load simulation state from checkpoint files (JSON + NPZ).

Restores simulation time and all block/event states from a previously saved checkpoint. Matching is based on block/event type and insertion order, so the simulation must be constructed with the same block types in the same order.

Parameters:: path (str) – base path without extension

add_block(block)[source]¶

Adds a new block to the simulation, initializes its local solver instance and collects internal events of the new block.

This works dynamically for running simulations.

Parameters:: block (Block) – block to add to the simulation

remove_block(block)[source]¶

Removes a block from the simulation.

This works dynamically for running simulations. The graph is lazily rebuilt on the next simulation update.

Parameters:: block (Block) – block to remove from the simulation

add_connection(connection)[source]¶

Adds a new connection to the simulation and checks if the new connection overwrites any existing connections.

This works dynamically for running simulations.

Parameters:: connection (Connection) – connection to add to the simulation

remove_connection(connection)[source]¶

Removes a connection from the simulation.

This works dynamically for running simulations. The graph is lazily rebuilt on the next simulation update.

Parameters:: connection (Connection) – connection to remove from the simulation

add_event(event)[source]¶

Checks and adds a new event to the simulation.

This works dynamically for running simulations.

Parameters:: event (Event) – event to add to the simulation

remove_event(event)[source]¶

Removes an event from the simulation.

This works dynamically for running simulations.

Parameters:: event (Event) – event to remove from the simulation

reset(time=0.0)[source]¶

Reset the blocks to their initial state and the global time of the simulation.

For recording blocks such as ‘Scope’, their recorded data is also reset.

Resets linearization automatically, since resetting the blocks resets their internal operators.

Afterwards the system function is evaluated with ‘_update’ to update the block inputs and outputs.

Parameters:: time (float) – simulation time for reset

linearize()[source]¶

Linearize the full system in the current simulation state at the current simulation time.

This is achieved by linearizing algebraic and dynamic operators of the internal blocks. See definition of the ‘Block’ class.

Before linearization, the global system function is evaluated to get the blocks into the current simulation state. This is only really relevant if no solving attempt has been happened before.

delinearize()[source]¶: Revert the linearization of the full system.

steadystate(reset=False)[source]¶

Find steady state solution (DC operating point) of the system by switching all blocks to steady state solver, solving the fixed point equations, then switching back.

The steady state solver forces all the temporal derivatives, i.e. the right hand side equation (including external inputs) of the engines of dynamic blocks to zero.

Note

This is really a sort of pseudo-steady-state solve. It does NOT compute the limit \(t\rightarrow\infty\) but rather forces all time derivatives to zero at a given moment in time.

This means, for a given t it computes the block states x such that:

\[0 = f(x, t)\]

instead of the real steady state:

\[\lim_{t \rightarrow \infty} x(t)\]

Parameters:: reset (bool) – reset the simulation before solving for steady state (default False)

estimate_initial_step(order=None)[source]¶

Estimate a good initial timestep from the system state and dynamics.

Implements the automatic initial step size selection of Hairer, Norsett and Wanner: a first guess is formed from the scaled norms of the state and its derivative, then refined with an estimate of the second derivative obtained from an explicit Euler probe. The scaling uses the absolute and relative error tolerances of the integration engine.

Note

Blocks expose their state derivative through the ‘Block.derivative’ method. Blocks that do not define it (stateless or with a custom dynamic path) are skipped, and if no derivative information is available the configured timestep is returned unchanged.

Parameters:: order (int, None) – order of the integration scheme, defaults to the order of the engine
Returns:: dt – estimated initial timestep
Return type:: float

References

pss(period, transient=0.0, iterations_max=50, tolerance=1e-08, dt=None, reset=False)[source]¶

Find the periodic steady state (limit cycle) of the system by shooting.

For a system with period T the periodic steady state is the initial condition \(x_0\) that is invariant under the period map \(P\) which integrates the system over one period,

\[P(x_0) = x_0, \qquad P(x_0) = x(t_0 + T;\, x_0)\]

The fixed point is found by an Anderson accelerated shooting iteration on \(x_0\), where each shot integrates one period with fixed timesteps. Data sampling is suppressed during the shooting iterations, so recording blocks stay clean; after convergence the engines are left at the limit-cycle initial state and the simulation time at the cycle start, ready for a final recording run over one period.

Note

This is a matrix-free shooting method (no monodromy matrix) and is best suited to single-step integrators. An optional transient run can be used to approach the limit cycle before shooting. Event handling is not considered during the shooting iterations.

Parameters:

period (float) – period of the steady state
transient (float) – duration of an optional transient run to approach the limit cycle
iterations_max (int) – maximum number of shooting iterations
tolerance (float) – convergence tolerance on the shooting residual
dt (float, None) – fixed timestep for the period integration, defaults to self.dt
reset (bool) – reset the simulation before solving (default False)

Returns:

success (bool) – whether the shooting iteration converged
iterations (int) – number of shooting iterations used
residual (float) – shooting residual at the last iteration

References

timestep_fixed_explicit(dt=None)[source]¶

Advances the simulation by one timestep ‘dt’ for explicit fixed step solvers.

Deprecated since version 1.0.0: Use timestep() instead.

Parameters:

dt (float) – timestep

Returns:

success (bool) – indicator if the timestep was successful
max_error (float) – maximum local truncation error from integration
scale (float) – rescale factor for timestep
total_evals (int) – total number of system evaluations
total_solver_its (int) – total number of implicit solver iterations

timestep_fixed_implicit(dt=None)[source]¶

Advances the simulation by one timestep ‘dt’ for implicit fixed step solvers.

Deprecated since version 1.0.0: Use timestep() instead.

Parameters:

dt (float) – timestep

Returns:

success (bool) – indicator if the timestep was successful
max_error (float) – maximum local truncation error from integration
scale (float) – rescale factor for timestep
total_evals (int) – total number of system evaluations
total_solver_its (int) – total number of implicit solver iterations

timestep_adaptive_explicit(dt=None)[source]¶

Advances the simulation by one timestep ‘dt’ for explicit adaptive solvers.

Deprecated since version 1.0.0: Use timestep() instead.

Parameters:

dt (float) – timestep

Returns:

success (bool) – indicator if the timestep was successful
max_error (float) – maximum local truncation error from integration
scale (float) – rescale factor for timestep
total_evals (int) – total number of system evaluations
total_solver_its (int) – total number of implicit solver iterations

timestep_adaptive_implicit(dt=None)[source]¶

Advances the simulation by one timestep ‘dt’ for implicit adaptive solvers.

Deprecated since version 1.0.0: Use timestep() instead.

Parameters:

dt (float) – timestep

Returns:

success (bool) – indicator if the timestep was successful
max_error (float) – maximum local truncation error from integration
scale (float) – rescale factor for timestep
total_evals (int) – total number of system evaluations
total_solver_its (int) – total number of implicit solver iterations

timestep(dt=None, adaptive=True)[source]¶

Advances the transient simulation by one timestep ‘dt’.

Automatic behavior selection based on selected Solver and adaptive flag:

Explicit solvers: Uses _update() for system evaluation
Implicit solvers: Uses _solve() for implicit update equation
Adaptive solvers (with adaptive=True): Reverts timestep if error too large or event not close
Fixed solvers (or adaptive=False): Always completes timestep, resolves events at detected time

If discrete events are detected, they are handled according to stepping mode:

Fixed stepping: Events resolved at interpolated time within step
Adaptive stepping: Events approached via timestep rescaling (secant method)

Parameters:

dt (float) – timestep size for transient simulation
adaptive (bool) – explicitly enable/disable adaptive timestepping; when False, adaptive solvers are forced to take fixed steps without error control (default True)

Returns:

success (bool) – indicator if the timestep was successful
max_error (float) – maximum local truncation error from integration
scale (float) – rescale factor for timestep
total_evals (int) – total number of system evaluations
total_solver_its (int) – total number of implicit solver iterations

step(dt=None, adaptive=True)[source]¶: Wraps ‘Simulation.timestep’ for backward compatibility

collect()[source]¶

Collect all current simulation results from the internal

Deprecated since version 1.0.0: its against pathsims philosophy

recording blocks

Returns:: results
Return type:: dict

stop()[source]¶: Set the flag for active simulation to ‘False’, intended to be called from the outside (for example by events) to interrupt the timestepping loop in ‘run’.

run(duration=10, reset=False, adaptive=True)[source]¶

Perform multiple simulation timesteps for a given ‘duration’.

Tracks the total number of block evaluations (proxy for function calls, although larger, since one function call of the system equation consists of many block evaluations) and the total number of solver iterations for implicit solvers.

Additionally the progress of the simulation is tracked by a custom ‘ProgressTracker’ class that is a dynamic generator and interfaces the logging system.

Parameters:

duration (float) – simulation time (in time units)
reset (bool) – reset the simulation before running (default False)
adaptive (bool) – use adaptive timesteps if solver is adaptive (default True)

Returns:

stats – stats of simulation run tracked by the ‘ProgressTracker’

Return type:

dict

run_streaming(duration=10, reset=False, adaptive=True, tickrate=10, func_callback=None)[source]¶

Perform simulation with streaming output at a fixed wall-clock rate.

This method runs the simulation as fast as possible while yielding intermediate results at a fixed rate defined by ‘tickrate’. Useful for real-time visualization and UI updates.

The progress is tracked and logged using the ‘ProgressTracker’ class.

Parameters:

duration (float) – simulation time (in time units)
reset (bool) – reset the simulation before running (default False)
adaptive (bool) – use adaptive timesteps if solver is adaptive (default True)
tickrate (float) – output rate in Hz, i.e., yields per second of wall-clock time (default 10)
func_callback (callable | None) – callback function that is called at every tick, can be used for data extraction, its return value is yielded by this generator

Yields:

result – The return value of the ‘func_callback’ callable.

run_realtime(duration=10, reset=False, adaptive=True, tickrate=30, speed=1.0, func_callback=None)[source]¶

Perform simulation paced to wall-clock time.

This method runs the simulation synchronized to real time, optionally scaled by ‘speed’. The simulation advances to match elapsed wall-clock time, yielding results at the rate defined by ‘tickrate’.

Useful for interactive simulations, hardware-in-the-loop testing, or when simulation should match real-world timing.

The progress is tracked and logged using the ‘ProgressTracker’ class.

Parameters:

duration (float) – simulation time (in time units)
reset (bool) – reset the simulation before running (default False)
adaptive (bool) – use adaptive timesteps if solver is adaptive (default True)
tickrate (float) – output rate in Hz, i.e., yields per second of wall-clock time (default 30)
speed (float) – time scaling factor where 1.0 is real-time, 2.0 is twice as fast, 0.5 is half speed (default 1.0)
func_callback (callable | None) – callback function that is called at every tick, can be used for data extraction, its return value is yielded by this generator

Yields:

result – The return value of the ‘func_callback’ callable.