Development roadmap for PathSim.
If simulation doesnt converge, we want to know where the problem comes from. Which block for the ode solvers. Which connection for the algebraic loop solvers. We want block wise error and convergence metrics.
Luckiely everything is there already, we just need to expose it through a nice interface.
The adaptive timestep solvers currently use a P controler to adjust the timestep given an estimate for the current truncation error. In the literature and in practice, different kinds of controlers have been explored for that such as I, PI and PID controlers.
It would make sense to implement a new class TruncationErrorControler to modularize this and enable different controler types.
Some runge kutta solvers such as RKBS32 and RKDP54 have the first-same-as-last (FSAL) property, where the last slope of the previous stage is the first slope of the next one. This saves one function evaluation per successful step and therefore increases solver efficiency.
Currently this is not implemented in pathsim, but it would give some additional 5% to 20% speedup for some solvers.
Currently, PathSim is not explicitly typed. Adding type hints makes sense.
Sometimes in big complex systems, individual components have wildly different timescales for their physics. In some cases it makes sense to approximate components with very fast dynamics as being in a steady state at each timestep, such that the component becomes purely algebraic.
To achieve this, the time derivative of the block ode
MATHDISPLAY0ENDMATH
will be forced to zero (transients have settled locally). Then the behavior of the block becomes algebraic and is defined by this algebraic implicit equation
MATHDISPLAY1ENDMATH
that has to be solved for x for given block inputs u and time t.
Implementation ideas:
- define new
pseudosteadystate(t)method for blocks and implement it for blocks that have dynamic operators - add a new
_pssflag to blocks and enablepssas an arg for the blocks that support it - in the timestepping methods of the block check if
_pssis set and handle the logic accordingly
This will work but maybe there is a more elegant solution
- alternatively build a
PeriodicSteadyStateblock that wraps other blocks and solves the algebraic equation through the same interface as theSimulationusing the already existingSteadyStatesolver
I think thats more elegant but might be less efficient
PathSim has a decentralized architecture for the blocks which lends itself to parallelism and asynchronizity. Expensive blocks should compute asynchronously and not make the other blocks wait. With free-threading from Python 3.13, parallelization of the block updates is possible and has been verified with multiprocessing (slow but validation of the concept) for an earlier build.
Near linear scaling of performance with number of cpus for large numbers of blocks can be expected.
Implement a copy method for the blocks, the Subsystem class, and the Simulation.
This should enable convenient copying of standard blocks for defining a system.
Implementing implicit-explicit ode solvers.
Some blocks in large systems may exhibit local stiffness while the coupling to external blocks is non-stiff. In these cases it would be nice to use more stable implicit ode solvers for these blocks while using explicit solvers for the other, non-stiff blocks.
The global solver would remain explicit, while locally, blocks can be flagged as stiff and then use the implicit counterpart. Implicit iterations will only be done locally for that block. Methods that would suit this are:
- ARS443 pair of explicit and diagonally implicit runge kutta methods
- ARK4(3)7L[2]SA pair of explicit and diagonally implicit runge kutta methods with embedded pairs for error control
Using exponential integrators is a way to eliminate stiffness from linear dynamical systems. Many pathsim blocks are pure linear odes such as the StateSpace blocks and its derivates, as well as the Differentiator and the PID.
They are more or less of the following form:
MATHDISPLAY0ENDMATH
Stiffness occurs when the eigenvalues of A are on vastly different scales, i.e. there are fast and slow time constants in the system.
Typpically the ODE is solved from one timestep to the next like this:
MATHDISPLAY1ENDMATH
But since the exact solution can be expressed for linea systems, the problem can be reduced to:
MATHDISPLAY2ENDMATH
It removes the state dependency from the integration. Intuitively it can be explained by superposition due to the linearity. The solution is an exponentially dacaying state added to the input convolved with the impulse response.
The resulting convolution integral can then be solved by some ode solver.
Implementing this as an option for linear ode blocks will help reduce stiffness induced by these kinds of blocks.