#########################################################################################

##

## MAIN SIMULATION ENGINE

## (simulation.py)

##

## This module contains the simulation class that manages

## the blocks, connections, events and specific simulation methods.

##

#########################################################################################

# IMPORTS ===============================================================================

import json

import warnings

import numpy as np

import time

import datetime

import logging

from pathsim import __version__

from ._constants import (

SIM_TIMESTEP,

SIM_TIMESTEP_MIN,

SIM_TIMESTEP_MAX,

SIM_TOLERANCE_FPI,

SIM_ITERATIONS_MAX,

LOG_ENABLE

)

from .optim.booster import ConnectionBooster

from .utils.graph import Graph

from .utils.analysis import Timer

from .utils.deprecation import deprecated

from .utils.portreference import PortReference

from .utils.progresstracker import ProgressTracker

from .utils.diagnostics import Diagnostics, ConvergenceTracker, StepTracker

from .utils.logger import LoggerManager

from .solvers import SSPRK22, SteadyState

from .blocks._block import Block

from .events._event import Event

from .connection import Connection

# TRANSIENT SIMULATION CLASS ============================================================

class Simulation:

"""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:

.. code-block:: python

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´

Attributes

----------

time : float

global simulation time, starting at ´0.0´

graph : Graph

internal graph representation for fast system funcion evluations

using DAG with algebraic depths

boosters : None | list[ConnectionBooster]

list of boosters (fixed point accelerators) that wrap algebraic

loop closing connections assembled from the system graph

engine : Solver

global integrator (ODE solver) instance serving as a dummy to

get attributes and access to intermediate evaluation stages

logger : logging.Logger

global simulation logger

_blocks_dyn : list[Block]

blocks with internal ´Solver´ instances (stateful)

_blocks_evt : list[Block]

blocks with internal events (discrete time, eventful)

_active : bool

flag for setting the simulation as active, used for interrupts

"""

def __init__(

self,

blocks=None,

connections=None,

events=None,

dt=SIM_TIMESTEP,

dt_min=SIM_TIMESTEP_MIN,

dt_max=SIM_TIMESTEP_MAX,

Solver=SSPRK22,

tolerance_fpi=SIM_TOLERANCE_FPI,

iterations_max=SIM_ITERATIONS_MAX,

log=LOG_ENABLE,

diagnostics=False,

**solver_kwargs

):

#system definition

self.blocks = []

self.connections = []

self.events = []

#simulation timestep and bounds

self.dt = dt

self.dt_min = dt_min

self.dt_max = dt_max

#numerical integrator to be used (class definition)

self.Solver = Solver

#numerical integrator instance

self.engine = Solver()

#internal system graph -> initialized later

self.graph = None

self._graph_dirty = False

#internal algebraic loop solvers -> initialized later

self.boosters = None

#error tolerance for fixed point loop and implicit solver

self.tolerance_fpi = tolerance_fpi

#additional solver parameters

self.solver_kwargs = solver_kwargs

#iterations for fixed-point loop

self.iterations_max = iterations_max

#enable logging flag

self.log = log

#initial simulation time

self.time = 0.0

#collection of blocks with internal ODE solvers

self._blocks_dyn = []

#collection of blocks with internal events

self._blocks_evt = []

#flag for setting the simulation active

self._active = True

#convergence trackers for the three solver loops

self._loop_tracker = ConvergenceTracker()

self._solve_tracker = ConvergenceTracker()

self._step_tracker = StepTracker()

#diagnostics snapshot (None when disabled)

self.diagnostics = Diagnostics() if diagnostics else None

#diagnostics history (list of snapshots per timestep)

self._diagnostics_history = [] if diagnostics == "history" else None

#initialize logging

logger_mgr = LoggerManager(

enabled=bool(self.log),

output=self.log if isinstance(self.log, str) else None,

level=logging.INFO,

date_format='%H:%M:%S'

)

self.logger = logger_mgr.get_logger("simulation")

self.logger.info(f"LOGGING (log: {self.log})")

#prepare and add blocks (including internal events)

if blocks is not None:

for block in blocks:

self.add_block(block)

#check and add connections

if connections is not None:

for connection in connections:

self.add_connection(connection)

#check and add events

if events is not None:

for event in events:

self.add_event(event)

#check if blocks from connections are in simulation

self._check_blocks_are_managed()

#assemble the system graph for simulation

self._assemble_graph()

def __contains__(self, other):

"""Check if blocks, connections or events are

already part of the simulation

Paramters

---------

other : obj

object to check if its part of simulation

Returns

-------

bool

"""

return (

other in self.blocks or

other in self.connections or

other in self.events

)

def __bool__(self):

"""Boolean evaluation of Simulation instances

Returns

-------

active : bool

is the simulation active

"""

return self._active

# methods for access to metadata ----------------------------------------------

@property

def size(self):

"""Get size information of the simulation, such as total number

of blocks and dynamic states, with recursive retrieval from subsystems

Returns

-------

sizes : tuple[int]

size of simulation (number of blocks) and number

of internal states (from internal engines)

"""

total_n, total_nx = 0, 0

for block in self.blocks:

n, nx = block.size

total_n += n

total_nx += nx

return total_n, total_nx

# visualization ---------------------------------------------------------------

def plot(self, *args, **kwargs):

"""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

"""

for block in self.blocks:

if block: block.plot(*args, **kwargs)

# checkpoint methods ----------------------------------------------------------

@staticmethod

def _checkpoint_key(type_name, type_counts):

"""Generate a deterministic checkpoint key from block/event type

and occurrence index (e.g. 'Integrator_0', 'Scope_1').

Parameters

----------

type_name : str

class name of the block or event

type_counts : dict

running counter per type name, mutated in place

Returns

-------

key : str

deterministic checkpoint key

"""

idx = type_counts.get(type_name, 0)

type_counts[type_name] = idx + 1

return f"{type_name}_{idx}"

def save_checkpoint(self, path, recordings=True):

"""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)

"""

#strip extension if provided

if path.endswith('.json') or path.endswith('.npz'):

path = path.rsplit('.', 1)[0]

#simulation metadata

checkpoint = {

"version": "1.0.0",

"pathsim_version": __version__,

"created": datetime.datetime.now(datetime.timezone.utc).isoformat(),

"simulation": {

"time": self.time,

"dt": self.dt,

"dt_min": self.dt_min,

"dt_max": self.dt_max,

"solver": self.Solver.__name__,

"tolerance_fpi": self.tolerance_fpi,

"iterations_max": self.iterations_max,

},

"blocks": [],

"events": [],

}

npz_data = {}

#checkpoint all blocks (keyed by type + insertion index)

type_counts = {}

for block in self.blocks:

key = self._checkpoint_key(block.__class__.__name__, type_counts)

b_json, b_npz = block.to_checkpoint(key, recordings=recordings)

b_json["_key"] = key

checkpoint["blocks"].append(b_json)

npz_data.update(b_npz)

#checkpoint external events (keyed by type + insertion index)

type_counts = {}

for event in self.events:

key = self._checkpoint_key(event.__class__.__name__, type_counts)

e_json, e_npz = event.to_checkpoint(key)

e_json["_key"] = key

checkpoint["events"].append(e_json)

npz_data.update(e_npz)

#write files

with open(f"{path}.json", "w", encoding="utf-8") as f:

json.dump(checkpoint, f, indent=2, ensure_ascii=False)

np.savez(f"{path}.npz", **npz_data)

def load_checkpoint(self, path):

"""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

"""

#strip extension if provided

if path.endswith('.json') or path.endswith('.npz'):

path = path.rsplit('.', 1)[0]

#read files

with open(f"{path}.json", "r", encoding="utf-8") as f:

checkpoint = json.load(f)

npz = np.load(f"{path}.npz", allow_pickle=False)

try:

#version check

cp_version = checkpoint.get("pathsim_version", "unknown")

if cp_version != __version__:

warnings.warn(

f"Checkpoint was saved with PathSim {cp_version}, "

f"current version is {__version__}"

)

#restore simulation state

sim_data = checkpoint["simulation"]

self.time = sim_data["time"]

self.dt = sim_data["dt"]

self.dt_min = sim_data["dt_min"]

self.dt_max = sim_data["dt_max"]

#solver type check

if sim_data["solver"] != self.Solver.__name__:

warnings.warn(

f"Checkpoint solver '{sim_data['solver']}' differs from "

f"current solver '{self.Solver.__name__}'"

)

#index checkpoint blocks by key

block_data = {b["_key"]: b for b in checkpoint.get("blocks", [])}

#restore blocks by type + insertion order

type_counts = {}

for block in self.blocks:

key = self._checkpoint_key(block.__class__.__name__, type_counts)

if key in block_data:

block.load_checkpoint(key, block_data[key], npz)

else:

warnings.warn(

f"Block '{key}' not found in checkpoint"

)

#index checkpoint events by key

event_data = {e["_key"]: e for e in checkpoint.get("events", [])}

#restore external events by type + insertion order

type_counts = {}

for event in self.events:

key = self._checkpoint_key(event.__class__.__name__, type_counts)

if key in event_data:

event.load_checkpoint(key, event_data[key], npz)

else:

warnings.warn(

f"Event '{key}' not found in checkpoint"

)

finally:

npz.close()

# adding system components ----------------------------------------------------

def add_block(self, block):

"""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

"""

#check if block already in block list

if block in self.blocks:

_msg = f"block {block} already part of simulation"

self.logger.error(_msg)

raise ValueError(_msg)

#initialize numerical integrator of block with parent

block.set_solver(self.Solver, self.engine, **self.solver_kwargs)

#add to dynamic list if solver was initialized

if block.engine:

self._blocks_dyn.append(block)

#add to eventful list if internal events

if block.events:

self._blocks_evt.append(block)

#add block to global blocklist

self.blocks.append(block)

#mark graph for rebuild

if self.graph:

self._graph_dirty = True

def remove_block(self, block):

"""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

"""

#check if block is in block list

if block not in self.blocks:

_msg = f"block {block} not part of simulation"

self.logger.error(_msg)

raise ValueError(_msg)

#remove from global blocklist

self.blocks.remove(block)

#remove from dynamic list

if block in self._blocks_dyn:

self._blocks_dyn.remove(block)

#remove from eventful list

if block in self._blocks_evt:

self._blocks_evt.remove(block)

#mark graph for rebuild

if self.graph:

self._graph_dirty = True

def add_connection(self, connection):

"""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

"""

#check if connection already in connection list

if connection in self.connections:

_msg = f"{connection} already part of simulation"

self.logger.error(_msg)

raise ValueError(_msg)

#add connection to global connection list

self.connections.append(connection)

#mark graph for rebuild

if self.graph:

self._graph_dirty = True

def remove_connection(self, connection):

"""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

"""

#check if connection is in connection list

if connection not in self.connections:

_msg = f"{connection} not part of simulation"

self.logger.error(_msg)

raise ValueError(_msg)

#remove from global connection list

self.connections.remove(connection)

#mark graph for rebuild

if self.graph:

self._graph_dirty = True

def add_event(self, event):

"""Checks and adds a new event to the simulation.

This works dynamically for running simulations.

Parameters

----------

event : Event

event to add to the simulation

"""

#check if event already in event list

if event in self.events:

_msg = f"{event} already part of simulation"

self.logger.error(_msg)

raise ValueError(_msg)

#add event to global event list

self.events.append(event)

def remove_event(self, event):

"""Removes an event from the simulation.

This works dynamically for running simulations.

Parameters

----------

event : Event

event to remove from the simulation

"""

#check if event is in event list

if event not in self.events:

_msg = f"{event} not part of simulation"

self.logger.error(_msg)

raise ValueError(_msg)

#remove from global event list

self.events.remove(event)

# system assembly -------------------------------------------------------------

def _assemble_graph(self):

"""Build the internal graph representation for fast system function

evaluation and algebraic loop resolution.

"""

#reset all block inputs to clear stale values from removed connections

for block in self.blocks:

block.inputs.reset()

#time the graph construction

with Timer(verbose=False) as T:

self.graph = Graph(self.blocks, self.connections)

self._graph_dirty = False

#create boosters for loop closing connections

if self.graph.has_loops:

self.boosters = [

ConnectionBooster(conn) for conn in self.graph.loop_closing_connections()

]

#log block summary

num_dynamic = len(self._blocks_dyn)

num_static = len(self.blocks) - num_dynamic

num_eventful = len(self._blocks_evt)

self.logger.info(

f"BLOCKS (total: {len(self.blocks)}, dynamic: {num_dynamic}, "

f"static: {num_static}, eventful: {num_eventful})"

)

#log graph info

self.logger.info(

"GRAPH (nodes: {}, edges: {}, alg. depth: {}, loop depth: {}, runtime: {})".format(

*self.graph.size, *self.graph.depth, T

)

)

# topological checks ----------------------------------------------------------

def _check_blocks_are_managed(self):

"""Check whether the blocks that are part of the connections are

in the simulation block set ('self.blocks') and therefore managed

by the simulation.

If not, there will be a warning in the logging.

"""

# Collect connection blocks

conn_blocks = set()

for conn in self.connections:

conn_blocks.update(conn.get_blocks())

# Check subset actively managed

for blk in conn_blocks:

if blk not in self.blocks:

self.logger.warning(

f"{blk} in 'connections' but not in 'blocks'!"

)

# solver management -----------------------------------------------------------

def _set_solver(self, Solver=None, **solver_kwargs):

"""Initialize all blocks with solver for numerical integration

and tolerance for local truncation error ´tolerance_lte´.

If blocks already have solvers, change the numerical integrator

to the ´Solver´ class.

Parameters

----------

Solver : Solver

numerical solver definition from ´pathsim.solvers´

solver_kwargs : dict

additional parameters for numerical solvers

"""

#update global solver class

if Solver is not None:

self.Solver = Solver

#update solver parmeters

self.solver_kwargs.update(solver_kwargs)

#initialize dummy engine to get solver attributes

self.engine = self.Solver()

#iterate all blocks and set integration engines with tolerances

self._blocks_dyn = []

for block in self.blocks:

block.set_solver(self.Solver, self.engine, **self.solver_kwargs)

#add dynamic blocks to list

if block.engine:

self._blocks_dyn.append(block)

#logging message

self.logger.info(

"SOLVER (dyn. blocks: {}) -> {} (adaptive: {}, explicit: {})".format(

len(self._blocks_dyn),

self.engine,

self.engine.is_adaptive,

self.engine.is_explicit

)

)

# resetting -------------------------------------------------------------------

def reset(self, time=0.0):

"""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

"""

self.logger.info(f"RESET (time: {time})")

#set active again

self._active = True

#reset simulation time

self.time = time

#reset integration engine

self.engine.reset()

#reset all blocks to initial state

for block in self.blocks:

block.reset()

#reset all event managers

for event in self.events:

event.reset()

#reset convergence trackers and diagnostics

self._loop_tracker.reset()

self._solve_tracker.reset()

self._step_tracker.reset()

if self.diagnostics is not None:

self.diagnostics = Diagnostics()

if self._diagnostics_history is not None:

self._diagnostics_history.clear()

#evaluate system function

self._update(self.time)

# linearization ---------------------------------------------------------------

def linearize(self):

"""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.

"""

#evaluate system function at current time

self._update(self.time)

#linearize all internal blocks and time it

with Timer(verbose=False) as T:

for block in self.blocks:

block.linearize(self.time)

self.logger.info(f"LINEARIZED (runtime: {T})")

def delinearize(self):

"""Revert the linearization of the full system."""

for block in self.blocks:

block.delinearize()

self.logger.info("DELINEARIZED")

# event system helpers --------------------------------------------------------

def _get_active_events(self):

"""Generator that yields all active events from simulation

and internal block events.

"""

for event in self.events:

if event:

yield event

for block in self._blocks_evt:

for event in block.events:

if event:

yield event

def _estimate_events(self, t):

"""Estimate the time until the next.

Parameters

----------

t : float

evaluation time for event estimation

Returns

-------

float | None

esimated time until next event (delta)

"""

dt_evt_min = None

#check external events

for event in self._get_active_events():

#get the estimate

dt_evt = event.estimate(self.time)

#no estimate available

if dt_evt is None: continue

#smaller than min

if dt_evt_min is None or dt_evt < dt_evt_min:

dt_evt_min = dt_evt

#return time until next event or None

return dt_evt_min

def _detected_events(self, t):

"""Check for possible (active) events and return them chronologically,

sorted by their timestep ratios (closest to the initial point in time).

Parameters

----------

t : float

evaluation time for event function

Returns

-------

detected : list[Event]

list of detected events within timestep

"""

#iterate all event managers

detected_events = []

for event in self._get_active_events():

#check if an event is detected

detected, close, ratio = event.detect(t)

#event was detected during the timestep

if detected:

detected_events.append([event, close, ratio])

#return detected events sorted by ratio

return sorted(detected_events, key=lambda e: e[-1])

# solving system equations ----------------------------------------------------

def _update(self, t):

"""Distribute information within the system by evaluating the directed acyclic graph

(DAG) formed by the algebraic passthroughs of the blocks and resolving algebraic loops

through accelerated fixed-point iterations.

Effectively evaluates the right hand side function of the global

system ODE/DAE

.. math::

\\begin{equnarray}

\\dot{x} &= f(x, t) \\\\

0 &= g(x, t)

\\end{equnarray}

by converging the whole system (´f´ and ´g´) to a fixed-point at a given point

in time ´t´.

If no algebraic loops are present in the system, convergence is

guaranteed after the first stage (evaluation of the DAG in '_dag').

Otherwise, accelerated fixed-point iterations ('_loops') are performed as a second

stage on the DAGs (broken cycles) of blocks that are part of or tainted by upstream

algebraic loops.

Parameters

----------

t : float

evaluation time for system function

"""

#lazy graph rebuild if dirty

if self._graph_dirty:

self._assemble_graph()

self._graph_dirty = False

#evaluate DAG

self._dag(t)

#algebraic loops -> solve them

if self.graph.has_loops:

self._loops(t)