Documentation

• 1: User Guide
• 1.1: Install
• 1.2: Examples
• 1.3: Logging Results
• 1.4: Real-Time
• 1.5: Co-simulation and Interfaces
• 2: Tutorials
• 2.1: Python Tutorials
• 2.1.1: Your First Simulation
• 2.1.2: Adding Dynamics
• 2.1.3: A Network, and Where It Starts From
• 2.1.4: Applying a Fault
• 2.1.5: The Same Circuit in Two Domains
• 2.1.6: Adding a Machine
• 2.1.7: Exchanging Data With Another Tool
• 3: Developer Guide
• 3.1: Architecture and Conventions
• 3.1.1: Architecture
• 3.1.2: Build
• 3.1.3: Coding Conventions
• 3.2: Attributes and Scheduling
• 3.2.1: Attributes
• 3.2.2: Attribute Usage Guidelines
• 3.2.3: Interfaces
• 3.2.4: Task Scheduling
• 3.2.5: Adding Tasks to a Component
• 3.3: Writing a Model
• 3.3.1: Component and Solver Initialization
• 3.3.2: Real-Time Execution
• 3.3.3: Interfacing with the MNA Solver
• 3.3.4: Subcomponent Handling
• 3.3.5: Add New Model
• 3.3.6: Create New Simulation
• 3.3.7: Debugging
• 3.3.8: Logger Implementation
• 3.4: Solvers
• 3.4.1: The MNA Solver
• 3.4.2: Power Flow Solvers
• 3.4.3: Alternative Solver Implementation
• 3.4.4: State-Space Extraction
• 3.5: Model Implementations
• 3.5.1: Reduced Order Generator Implementation
• 3.5.2: Switch and Load Implementation
• 3.5.3: Source Implementation
• 3.5.4: Injection and Compensation Implementation
• 3.5.5: SSN Component Implementation
• 3.5.6: Signal Component Implementation
• 3.5.7: DP Ph1 Averaged VSI Implementation
• 3.5.8: DP Ph3 Averaged VSI Implementation
• 3.5.9: EMT Ph3 Averaged VSI Implementation
• 3.5.10: EMT Ph3 Grid-Forming VSI Implementation
• 4: Concepts
• 4.1: Nodal Analysis
• 4.2: Dynamic Phasors
• 4.3: Powerflow
• 4.4: State-Space Nodal
• 4.5: State-Space Nodal Across Domains
• 4.6: State-Space Extraction
• 4.7: Modal Analysis
• 4.8: Alternative Solution Methods
• 4.9: Models
• 4.9.1: RLC-Elements
• 4.9.2: Branches
• 4.9.3: Transformer
• 4.9.4: Ideal Transformer Model
• 4.9.5: Switches
• 4.9.6: Loads
• 4.9.7: Sources
• 4.9.8: Network Injection and Compensation
• 4.9.9: Synchronous Generator
• 4.9.9.1: Reduced Order Machine Models
• 4.9.10: Synchronous Generator Regulators
• 4.9.10.1: Exciters and Power System Stabiliser
• 4.9.10.2: Turbines and Governors
• 4.9.11: Power Electronics
• 4.9.11.1: EMT Ph3 Averaged Voltage Source Inverter
• 4.9.11.2: DP Ph1 Averaged Voltage Source Inverter
• 4.9.11.3: DP Ph3 Averaged Voltage Source Inverter
• 4.9.11.4: EMT Ph3 Grid-Forming Inverter
• 4.9.12: Converter Control
• 4.9.13: State-Space Nodal Components
• 4.9.14: Signal Processing Blocks
• 5: Reference
• 5.1: Model Availability
• 5.2: Signal Models
• 5.3: State-Space Extraction Support
• 6: Contributing
• 6.1: LLM Pull Request Review
• 7: How to Cite

1 - User Guide

For readers who want to run simulations: installation, building from source, and the tasks that come up while setting up and running a study.

To change DPsim itself rather than use it, see the developer guide. For the mathematics behind what the solvers compute, see concepts.

1.1 - Install

DPsim is a Python module and C++ library for dynamic power system simulation.

The quickest route to a result is the Python module: install it, then work through the tutorials, which build up one idea at a time from a source and a resistor. If you would rather read a finished study than build one, the example notebooks run complete scenarios and plot them.

Building from source is only needed for a platform without a published wheel, or to work on DPsim itself; see build.

Try it without installing

The example notebooks run in the browser with no local installation:

Python package

DPsim is published on PyPI and installs like any other Python package:

python3 -m venv venv
source venv/bin/activate
pip install dpsim

If you prefer conda, the equivalent is:

conda create -n dpsim python=3.13
conda activate dpsim
pip install dpsim

Docker

You need Docker installed first. The prepared image on Docker Hub bundles the module together with a JupyterLab session:

docker run -p 8888:8888 sogno/dpsim

Then open http://localhost:8888/lab?token=3adaa57df44cea75e60c0169e1b2a98ae8f7de130481b5bc.

Note that the image pins that access token in its startup command, so it is the same for everyone who runs the image. Publish the port on localhost only, as above, and do not expose it to an untrusted network.

To build the image yourself rather than pulling it, see build.

1.2 - Examples

If you are starting out, start with the tutorials instead. They work through one idea at a time in order, each as a complete script. This page is the inventory of what else the repository carries, which is the right thing once you know what you are looking for.

The examples come in both languages. The Python notebooks run a scenario and plot the result in one place, so they are the better way to see a complete study. The C++ examples are the better reference for using DPsim as a library, and are what the real time and co-simulation scenarios are written in.

Where each kind fits

Notebooks

Run them in the browser with no local installation:

Locally, they live under examples/Notebooks and need the Python package on the path, as described in the build section.

Understanding_DP.ipynb is worth reading early if dynamic phasors are new to you, since it builds the intuition the concepts section then formalises.

Continuing from a tutorial

Each tutorial ends at the point where the notebooks take over, so a category is usually the natural next step from the rung that introduced the idea.

C++ examples

Under dpsim/examples/cxx, built as part of a normal build and produced as executables in the build directory.

The target name is the source file name without its extension, regardless of which directory the source sits in, and the executables are written flat into the build tree. So to build and run a single example from your build directory:

cmake --build . --target DP_VS_RL1
./dpsim/examples/cxx/DP_VS_RL1

Co-simulation

The dpsim-villas examples exchange data with other simulators or with hardware through VILLASnode. These need WITH_VILLAS enabled and are therefore not available on Windows. See interfaces for how the coupling works and real time for running them against the wall clock.

1.3 - Logging Results

A simulation records nothing unless asked. Every quantity you want afterwards has to be named before the run, because the solver keeps only what the current step needs and discards the rest.

dpsimpy.Logger is the one to reach for. It is the CSV data logger and it is what nearly every example and notebook uses. There is also a real-time data logger, which is the right choice under a real-time timer for the reason given below, and an interface for sending results somewhere other than a file.

Only attributes can be logged

So the question is never “how do I log this value” but “is this value an attribute”. If it is, one line records it. If it is not, no logger option will reach it, and the answer is to expose it as an attribute in the model, which is a code change described under attributes.

print_attribute_list() on any object prints what it publishes, which is the reliable way to find out. Anything absent from that list cannot be logged.

The same fact explains a convenience: because attributes are the unit, a derived quantity that is itself an attribute is logged the same way as a terminal voltage, with no special handling.

Registering an attribute

logger = dpsimpy.Logger("my_simulation")
logger.log_attribute("n1.v", "v", n1)
logger.log_attribute("load.i_intf", "i_intf", load)

sim.add_logger(logger)

The three arguments are the column name you want in the output, the name of the attribute on the object, and the object itself. The first is yours to choose and is what you will key on when reading the file back; the second must match an attribute the object actually publishes.

Common attribute names are v on a node, and i_intf and v_intf on a component for the current through it and the voltage across it. Components publish their own states as well: a machine offers w_r, delta and Te, a converter its control states.

Naming an attribute that does not exist fails when the logger is set up, not at the end of the run, so a typo costs a second rather than a simulation.

Where the file goes

By default the file is logs/<logger name>.csv under the working directory, where the name is the one given to the Logger constructor.

Logger.set_log_dir changes that directory and Logger.get_log_dir reports it. Calling it is what produces the nested logs/<something>/<name>.csv layout the example notebooks use, so a script that does not call it gets the flat form. Setting it is worth doing when one script runs several simulations, since two loggers with the same name otherwise write to the same file.

The same setting governs the diagnostic text log, which is why the two land side by side.

The file is plain CSV with a time column first and one column per logged attribute, in the order they were registered. A three-phase quantity becomes three columns suffixed _0, _1, _2, and a complex quantity in an envelope domain is written as a complex value that the reading side parses back.

Reading it back

import villas.dataprocessing.readtools as rt

results = rt.read_timeseries_dpsim("logs/my_simulation.csv")
voltage = results["n1.v"]

The keys are the column names you chose. Each value carries time and values arrays. In an envelope domain the values are complex, so abs() gives the magnitude, and frequency_shift_list recovers the waveform as shown in comparing domains.

The first row is written before the first solve, so it holds the state the simulation started from rather than a result. A plot that appears to begin at zero usually begins at that row.

What it costs

Logging is opt-in and costs nothing for what you do not ask for, but what you do ask for is written every step. A run of 100 000 steps logging 50 attributes writes five million values, and on a large network the file, not the solve, becomes the slow part.

Three things help. Log the attributes you will actually look at rather than everything available. Prefer a specific attribute over a whole matrix; the rows_max and cols_max arguments to log_attribute cap how much of a matrix quantity is written. And note that a large time step reduces the file in exact proportion, so a study that only needs the envelope does not need the step of one that needs the waveform.

Data logging against diagnostic logging

The word covers two unrelated things and they are easy to confuse.

What this page describes is the data logger: numerical results, CSV, opt-in per attribute. The other is the diagnostic log, the text file recording what the solver did, controlled by dpsimpy.LogLevel and passed to component constructors. Raising a component’s log level makes it describe its own initialization and stamping in prose; it has no effect on the CSV.

They land in the same logs/ directory side by side, one as .csv and the other as .log, which is why they get mistaken for each other. A component constructed with LogLevel.debug produces a great deal of text and no additional results.

The real-time data logger mentioned at the top is a data logger like the first, not a third kind of log. It records the same results and differs only in buffering them in memory and writing at the end, because a disk write inside a real-time step has no bound on how long it takes. See real-time.

1.4 - Real-Time

Normally a simulation runs as fast as it can: a one-second study finishes in whatever time the solver needs. In a real-time simulation, one second of simulated time takes one second of wall clock time, no faster and no slower.

When you need it

Only when something outside the simulation has its own clock. A controller running on real hardware, another simulator you are coupled to, or a person turning a dial all move at their own pace, and the simulation has to move with them. If nothing outside is waiting, running in real time only makes the study slower.

This is why co-simulation and real-time usually appear together.

What it changes

Two things. The simulation waits at the end of each step until the wall clock catches up, so a run takes as long as the time it simulates. And a step that takes longer than its own duration is an overrun: the simulation has missed its deadline and can no longer claim to be in step with the outside world.

Overruns are the whole difficulty. Everything else about real-time execution is arranging for them not to happen: keeping the step’s work bounded, and keeping the operating system from interrupting it.

What it takes

A time step that comfortably exceeds the work done in it, models that do nothing slow inside the step, and a host that will not interrupt at the wrong moment. Millisecond steps are undemanding; microsecond steps need a tuned kernel and careful models, and are where most of the effort goes.

The requirements on the host and on the models are in real-time execution. With that tuning, steps as low as 5 us synchronised to an FPGA through VILLASnode have been achieved.

Running a Real-Time Simulation

Before running a simulation, you can run the following commands as root:

echo "evacuating cores"
tuna isolate -c 9,11,13,15

echo "disabling RT throttling"
echo -1 > /proc/sys/kernel/sched_rt_runtime_us

echo "stopping systemd services"
systemctl stop polkit
systemctl stop containerd
systemctl stop crond
systemctl stop chronyd

As a reference, real-time simulation examples are provided in the dpsim/examples/cxx and dpsim-villas/examples/cxx folder of the DPsim repository.

To benefit from the PREEMPT_RT feature and the isolated cores, the simulation has to be started using the chrt command to set the scheduling policy and priority, and the taskset command to pin the process to the isolated cores.

• chrt man-page
• taskset man-page

In the following example, we set the FIFO scheduling policy with the highest priority (99) and pin the execution of the simulation to CPU cores 9,11,13,15 which have been reserved previously (see above).

# the simple RT_DP_CS_R_1 simulation
taskset -c 9,11,13,15 chrt -f 99 build/dpsim/examples/cxx/RT_DP_CS_R_1

# Cosimulation using VILLASnode, FPGA synchronized time step, and exchanging data via Aurora interface.
# Here we need sudo, to interact with the FPGA. We disable logging (log=false) and set the time step to 50 us (-t 0.00005).
sudo taskset -c 9,11,13,15 chrt -f 99 build/dpsim-villas/examples/cxx/FpgaCosim3PhInfiniteBus -o log=false -t 0.00005 -d 10

1.5 - Co-simulation and Interfaces

Interfaces can be used to exchange simulation signals between a DPsim simulation and other soft- or hardware, for example an MQTT-broker or an FPGA. Simulation signals in the form of Attributes can be imported or exported once per simulation time step. Interfaces are subclasses of Interface and implement the methods addExport and addImport, which add dependencies to the passed attribute that forward the attribute value from or to the interface. This way, attributes that are imported are read from the interface before they are used in any DPsim component. Attributes that are exported are written to the interface after they are set by a DPsim component.

Where the boundary is, and where to read further

This page documents only DPsim’s side of the interface: which attributes are exchanged, when they are read and written relative to the time step, and how the simulation is synchronized. Everything on the other side of the JSON configuration belongs to VILLASnode and is documented there rather than here, so the two should be read together.

The parts most often needed are these:

• Node types is the reference for the type key and its per-type options. Which protocols are available, and what each one requires, is decided here rather than in DPsim.
• Nodes covers the surrounding configuration structure that the node object sits in.
• Hooks describe the processing that can be applied to samples in transit, such as scaling, limiting or statistics. Anything that can be done with a hook does not need to be done in the simulation, which is usually the better place for it.

2 - Tutorials

Worked simulations in order of difficulty, starting from something trivial. Each one adds exactly one new idea to the one before it, and each is a complete script you can run rather than a fragment to assemble.

Every script on these pages was run before it was written up, and the numbers quoted are the numbers it produced. Where a result is surprising, the page says so rather than leaving you to wonder whether you typed something wrong.

Two tracks

Python is the track to start with, and the difference is not a matter of taste. With the Python package installed you edit a script and run it; there is no build step between a change and a result, and the loop is a few seconds long. That is the right way to learn what the simulator does.

A C++ track covers the same ground for anyone embedding the solver in an application or writing a new model, where the Python API is not the interface being used. It costs a compile and link on every change, so it is the wrong place to learn the concepts and the right place to work once you know them. It also requires a working build of DPsim itself rather than only the installed package; see build. The C++ examples are the material it is being built from.

The two tracks describe the same simulator and share the concept pages behind them. Only the calling code differs, so nothing learned in one track has to be relearned in the other.

If a tutorial is not what you want

The User Guide covers installation and individual features, Concepts has the mathematics behind the models, the Developer Guide covers changing DPsim rather than using it, and Reference has the generated API and the model availability tables.

2.1 - Python Tutorials

Each tutorial starts from the one before and adds exactly one new thing. Work through them in order; each is a complete runnable script rather than a fragment.

1. Your first simulation. A source and a resistor. The shape of a script, and how to read a result back.
2. Adding dynamics. An inductor, the transient it produces, and how to choose a time step.
3. A network, and where it starts from. A line between two buses, and initializing the dynamic run from a powerflow.
4. Applying a fault. Switching during a run, and why clearing a fault needs more care than applying one.
5. The same circuit in two domains. Waveforms against envelopes, and what the envelope buys.
6. Adding a machine. A synchronous generator, initializing it correctly, and what the model order changes.
7. Exchanging data with another tool. Handing a value out of a running simulation, and where the boundary lies.

A rung on converters and their control belongs between the last two and is not written yet.

What you need

DPsim importable from Python. If it is not, see install or build.

Reading results back and plotting them uses the data processing package the example notebooks also use. It is separate from DPsim and is imported as villas.dataprocessing.

2.1.1 - Your First Simulation

This page goes from nothing to a plotted result. It assumes DPsim is installed and importable; if it is not, start with install or build.

The circuit is deliberately trivial, a voltage source feeding a resistor, so that nothing in it distracts from the shape of the script. Every real simulation has the same five parts in the same order.

The whole script

import dpsimpy
import villas.dataprocessing.readtools as rt
from villas.dataprocessing.timeseries import TimeSeries as ts

name = "first_simulation"

# 1. Nodes
gnd = dpsimpy.dp.SimNode.gnd
n1 = dpsimpy.dp.SimNode("n1")

# 2. Components
src = dpsimpy.dp.ph1.VoltageSource("src")
src.V_ref = complex(100, 0)

load = dpsimpy.dp.ph1.Resistor("load")
load.R = 10.0

# 3. Connections
src.connect([gnd, n1])
load.connect([n1, gnd])

# 4. Topology
system = dpsimpy.SystemTopology(50, [gnd, n1], [src, load])

# 5. Logging, then run
logger = dpsimpy.Logger(name)
logger.log_attribute("n1.v", "v", n1)
logger.log_attribute("load.i_intf", "i_intf", load)

sim = dpsimpy.Simulation(name)
sim.set_domain(dpsimpy.Domain.DP)
sim.set_system(system)
sim.set_time_step(1e-3)
sim.set_final_time(0.1)
sim.add_logger(logger)
sim.run()

Running it prints solver progress and writes logs/first_simulation.csv.

What each part is doing

Nodes come first because components connect to them, not to each other. SimNode.gnd is the reference node and is shared; every network needs it. Nodes are chosen from a domain namespace, dpsimpy.dp here, and a node from one domain cannot be connected to a component from another.

Components are created, then configured through their attributes. src.V_ref = complex(100, 0) sets the source reference as a complex phasor, because this is the dynamic phasor domain and a voltage is an envelope rather than an instantaneous value. In EMT the same field would carry a different meaning; see dynamic phasors.

Logging is opt-in, and takes three steps in order. Nothing is recorded unless a logger asks for it.

logger = dpsimpy.Logger(name)                          # 1. create it; the name becomes the file name
logger.log_attribute("n1.v", "v", n1)                  # 2. register each attribute you want
logger.log_attribute("load.i_intf", "i_intf", load)
sim.add_logger(logger)                                 # 3. attach it to the simulation, before run()

log_attribute takes the column name you want, the name of the attribute on the object, and the object itself. "v" on a node is its voltage; "i_intf" on a component is the current through it. The first argument is yours to choose and is the key you will use when reading the file back; the second must be an attribute the object actually publishes, and print_attribute_list() on the object shows what that is.

The order is what makes it work. A logger registers attributes before it is attached, and it must be attached before run(), because the column header is written from whatever is registered when the first row is written. A logger created but never passed to add_logger produces no file at all, which is the usual reason for a run that appears to have logged nothing.

Reading the results

results = rt.read_timeseries_dpsim("logs/" + name + ".csv")

print(sorted(results.keys()))
# ['load.i_intf', 'n1.v']

v = results["n1.v"]
print(v.time[1], abs(v.values[1]))
# 0.001 100.0

The keys are the column names given to log_attribute. Each value is a time series with time and values arrays; in a dynamic phasor simulation the values are complex, and abs() gives the envelope magnitude.

Note that the first sample is zero. The log is written before the first solve, so row zero is the state the simulation started from rather than a result. From t = 0.001 onwards this circuit sits at exactly 100 V and 10 A, which is what a 100 V source across 10 Ω should give.

Plotting

import villas.dataprocessing.plottools as pt

pt.plot_timeseries(1, results["n1.v"].abs())
pt.plot_timeseries(2, results["load.i_intf"].abs())

The first argument is a figure number, so repeated calls with the same number overlay curves on one axis. .abs() is needed for complex results; plotting a complex series directly is not meaningful.

Recovering the waveform from a dynamic phasor result

A dynamic phasor result is an envelope, not a waveform. To compare it against an instantaneous result, shift it back onto the carrier:

emt = ts.frequency_shift_list(results, 50)
print(sorted(emt.keys()))
# ['load.i_intf_shift', 'n1.v_shift']

The script

The complete script for this page is 01_first_simulation.py under examples/Python/Tutorials. The numbers quoted above are the numbers it prints, so if the two ever disagree the page is the one that is wrong.

Next

This circuit has no dynamics at all: it is a source and a resistor, so it reaches its final value in one step and stays there. The next step is to add an element that stores energy, which is where the choice of time step starts to matter and where a result becomes worth plotting.

The examples work through larger networks, and the models used here are derived under sources and RLC elements.

2.1.2 - Adding Dynamics

The circuit in your first simulation reaches its final value in a single step. A resistor is a purely algebraic element: it stores no energy, so the circuit has no state variable and its response to a change is instantaneous.

An inductor stores energy in its magnetic field, and its current cannot change instantaneously. That current becomes a state variable, the circuit becomes first order, and it acquires a transient worth looking at and a reason to care about the time step.

The circuit

A source, a resistor and an inductor in series. Only the inductor is new; everything else is the same shape as before.

import dpsimpy
import villas.dataprocessing.readtools as rt

name = "rl_circuit"

gnd = dpsimpy.dp.SimNode.gnd
n1 = dpsimpy.dp.SimNode("n1")
n2 = dpsimpy.dp.SimNode("n2")

src = dpsimpy.dp.ph1.VoltageSource("src")
src.V_ref = complex(100, 0)

r = dpsimpy.dp.ph1.Resistor("r")
r.R = 10.0

l = dpsimpy.dp.ph1.Inductor("l")
l.L = 0.05

src.connect([gnd, n1])
r.connect([n1, n2])
l.connect([n2, gnd])

system = dpsimpy.SystemTopology(50, [gnd, n1, n2], [src, r, l])

logger = dpsimpy.Logger(name)
logger.log_attribute("i_l", "i_intf", l)

sim = dpsimpy.Simulation(name)
sim.set_domain(dpsimpy.Domain.DP)
sim.set_system(system)
sim.set_time_step(1e-4)
sim.set_final_time(0.05)
sim.add_logger(logger)
sim.run()

current = rt.read_timeseries_dpsim("logs/" + name + ".csv")["i_l"]

A second node appears because the resistor and the inductor meet somewhere, and that junction is a node like any other. Components connect to nodes, never directly to each other, so a series chain of two elements always needs the node between them.

What to expect before running it

Two numbers are worth working out first, because they are what the result should be checked against.

The steady-state current follows from the impedance at the system frequency,

and the transient decays with the time constant $\tau = L/R = 5$ ms, so the circuit settles after roughly five of those, about 25 ms. The simulation runs for 50 ms, comfortably past that.

Getting 5.37 A at the end is the check that the circuit was built as intended. A wrong connection order or a missing component usually shows up here rather than as an error.

Why the time step matters now

Run the same circuit twice, once at 0.1 ms and once at 5 ms, and compare the inductor current:

At 5 ms the two disagree by nearly a factor of two. By 20 ms they agree to better than a tenth of a percent, and both end at the right steady-state value.

The rule that follows: choose the step against the fastest thing you need to see, not against the duration of the run or the value you expect at the end. Here the fastest thing is $\tau = 5$ ms, and 0.1 ms resolves it with room to spare.

What the values mean in this domain

The current is complex, and abs() gives the magnitude of the envelope rather than an instantaneous current. In this domain the 50 Hz oscillation is not in the numbers at all: it has been moved into the carrier and handled analytically, which is why a 0.1 ms step is generous here and would be merely adequate for the same circuit solved as a waveform.

That difference is the subject of a later step. For now it is enough to know that a flat line in a dynamic phasor result means a steady sinusoid, not a constant.

The script

The complete script for this page is 02_adding_dynamics.py under examples/Python/Tutorials. The numbers quoted above are the numbers it prints, so if the two ever disagree the page is the one that is wrong.

Next

The circuit still has one source and one branch. Next is a network with a line between two buses, where the state the simulation starts from stops being obvious.

The elements used here are derived under RLC elements, and the trapezoidal companion models behind them under nodal analysis.

2.1.3 - A Network, and Where It Starts From

The circuits so far started from nothing and settled. That is fine for a resistor and an inductor, and useless for a network: a real system is already running when you start looking at it, and the transient you care about is the one caused by an event, not by switching the whole grid on.

This tutorial builds a two-bus network, solves its steady state with a powerflow, and starts the dynamic simulation from that solution.

Part one: the powerflow

A powerflow is a different kind of simulation. It has no time step in any meaningful sense; it solves the algebraic steady state, iterating until the bus voltages are consistent with the specified powers.

import dpsimpy
import villas.dataprocessing.readtools as rt

Vnom = 20e3

n1pf = dpsimpy.sp.SimNode("n1", dpsimpy.PhaseType.Single)
n2pf = dpsimpy.sp.SimNode("n2", dpsimpy.PhaseType.Single)

slack = dpsimpy.sp.ph1.NetworkInjection("slack")
slack.set_parameters(voltage_set_point=Vnom)
slack.set_base_voltage(Vnom)
slack.modify_power_flow_bus_type(dpsimpy.PowerflowBusType.VD)

line = dpsimpy.sp.ph1.PiLine("line")
line.set_parameters(R=0.5, L=0.5 / 314, C=50e-6)
line.set_base_voltage(Vnom)

load = dpsimpy.sp.ph1.Load("load")
load.set_parameters(active_power=100e3, reactive_power=50e3, nominal_voltage=Vnom)
load.modify_power_flow_bus_type(dpsimpy.PowerflowBusType.PQ)

slack.connect([n1pf])
line.connect([n1pf, n2pf])
load.connect([n2pf])

system_pf = dpsimpy.SystemTopology(50, [n1pf, n2pf], [slack, line, load])

logger_pf = dpsimpy.Logger("pf")
logger_pf.log_attribute("v1", "v", n1pf)
logger_pf.log_attribute("v2", "v", n2pf)

sim_pf = dpsimpy.Simulation("pf")
sim_pf.set_system(system_pf)
sim_pf.set_domain(dpsimpy.Domain.SP)
sim_pf.set_solver(dpsimpy.Solver.NRP)
sim_pf.set_solver_component_behaviour(dpsimpy.SolverBehaviour.Initialization)
sim_pf.do_init_from_nodes_and_terminals(False)
sim_pf.set_time_step(0.1)
sim_pf.set_final_time(0.1)
sim_pf.add_logger(logger_pf)
sim_pf.run()

Four things here are new and none of them are optional.

The powerflow is built in the static phasor domain, dpsimpy.sp, whatever domain the dynamic run will use. The solver is set to Solver.NRP, the Newton-Raphson powerflow solver, rather than the default nodal solver.

modify_power_flow_bus_type is what makes the problem solvable. Every bus must declare which two of its four quantities are known: VD fixes voltage magnitude and angle, and there must be exactly one such bus, the slack, which absorbs whatever mismatch remains. PQ fixes active and reactive power, which is what a load specifies. Without these the powerflow has no boundary conditions.

set_base_voltage is required on components because the solver works in per unit, and do_init_from_nodes_and_terminals(False) tells the components not to try to initialize themselves from node voltages that do not exist yet, since establishing those voltages is the job this simulation is doing.

Running it gives 20 000 V at the slack and 20 075 V at the load bus. The load bus sitting above nominal is not an error: the line’s shunt capacitance supplies more reactive power at this load than the series impedance drops.

Part two: the dynamic run

The dynamic network is built separately, in the domain the simulation will actually use, and then takes its initial state from the powerflow solution.

n1 = dpsimpy.dp.SimNode("n1", dpsimpy.PhaseType.Single)
n2 = dpsimpy.dp.SimNode("n2", dpsimpy.PhaseType.Single)

slack_d = dpsimpy.dp.ph1.NetworkInjection("slack")
slack_d.set_parameters(V_ref=complex(Vnom, 0))

line_d = dpsimpy.dp.ph1.PiLine("line")
line_d.set_parameters(
    series_resistance=0.5,
    series_inductance=0.5 / 314,
    parallel_capacitance=50e-6,
)

load_d = dpsimpy.dp.ph1.RXLoad("load")
load_d.set_parameters(active_power=100e3, reactive_power=50e3, volt=Vnom)

slack_d.connect([n1])
line_d.connect([n1, n2])
load_d.connect([n2])

system_dp = dpsimpy.SystemTopology(50, [n1, n2], [slack_d, line_d, load_d])
system_dp.init_with_powerflow(systemPF=system_pf, domain=dpsimpy.Domain.DP)

logger = dpsimpy.Logger("dyn")
logger.log_attribute("v2", "v", n2)

sim = dpsimpy.Simulation("dyn")
sim.set_system(system_dp)
sim.set_domain(dpsimpy.Domain.DP)
sim.set_time_step(1e-3)
sim.set_final_time(0.05)
sim.add_logger(logger)
sim.run()

init_with_powerflow matches the two networks by node name and copies the solved voltages across, which is why the node names must agree between the two topologies. It is the only line connecting the two halves.

The result is the point of the whole exercise. The dynamic run starts at 20 074.8 V and ends at 20 074.9 V: it begins in steady state rather than settling into one. Without the powerflow it would start from zero and spend the first several cycles charging the line, and any event applied during that period would be mixed in with a startup transient that has nothing to do with the system.

Parameter names differ between domains

The same component takes different keyword names in different domains. The static phasor line takes R, L and C; the dynamic phasor line takes series_resistance, series_inductance and parallel_capacitance. The load is Load with nominal_voltage in the powerflow and RXLoad with volt in the dynamic run.

The script

The complete script for this page is 03_two_bus_network.py under examples/Python/Tutorials. The numbers quoted above are the numbers it prints, so if the two ever disagree the page is the one that is wrong.

Next

The network now starts where it should, so an event applied to it produces a clean response. Next is applying one: a fault, using a switch.

The powerflow method is described under powerflow, and the loads and lines used here under loads and branches.

2.1.4 - Applying a Fault

The network from the previous tutorial starts in steady state, so anything that happens to it now is a response to the event rather than to startup. This tutorial applies a fault at the load bus, clears it, and looks at what the clearing does.

Scheduling an event

A switch is an ordinary component. What makes it a fault is that its state is changed partway through the run by an event.

fault = dpsimpy.dp.ph1.Switch("fault")
fault.set_parameters(open_resistance=1e9, closed_resistance=10.0)
fault.open()
fault.connect([gnd, n2])

# ... build the topology including `fault` ...

sim.add_event(dpsimpy.event.SwitchEvent(0.1, fault, True))   # apply
sim.add_event(dpsimpy.event.SwitchEvent(0.2, fault, False))  # clear

The switch is created open and connected between the load bus and ground, so closing it puts a 10 Ω path to ground at that bus. SwitchEvent takes the time, the switch, and the state to move to: True closes, False opens.

The switch must be in the topology’s component list like anything else. A switch that is created, connected and given events but left out of the list produces a run with no fault and no error.

Set a time step small enough to resolve the event. At 0.1 ms the fault instant is captured within one step; a millisecond step would smear it.

What happens

The load bus sits at 20 080 V, drops to 19 090 V while the fault is on, and recovers afterwards. The drop is modest because a 10 Ω fault on a 20 kV bus is not a solid short and the source is stiff.

The interesting part is the instant of clearing:

The bus voltage rings between 34 kV and 9 kV within a millisecond, a 70% overshoot on a network that was in steady state a moment earlier.

Why, and what to do about it

This is not the physical response of the circuit. Opening the switch asks the simulation to interrupt the current flowing through the line inductance within one time step, and an inductor current cannot change instantaneously. With the trapezoidal companion model the result is a numerical oscillation that decays slowly, as explained under switches.

A real breaker does not do this, because an arc forms across the opening contacts and dissipates the stored energy over a short but finite interval. The variable-resistance switch reproduces that: it raises its resistance over several steps rather than in one.

Two changes are needed, and the second is easy to forget:

fault = dpsimpy.dp.ph1.varResSwitch("fault")
fault.set_parameters(open_resistance=1e9, closed_resistance=10.0)
fault.open()
fault.set_init_parameters(1e-4)   # must match the simulation time step

With the same fault at the same instant:

The oscillation is gone. The voltage rises smoothly to a 20 943 V peak, a 4% overshoot rather than 70%, and settles back to its pre-fault value.

Use the plain switch for switching that does not interrupt inductive current, and the variable-resistance switch for faults, particularly at a machine terminal or a transformer winding. The cost is that the system matrix changes on every step of the transition rather than once, so each of those steps needs a refactorisation.

The script

The complete script for this page is 04_applying_a_fault.py under examples/Python/Tutorials. The numbers quoted above are the numbers it prints, so if the two ever disagree the page is the one that is wrong.

Next

The results so far have been envelopes. Next is running the same circuit as instantaneous waveforms and comparing the two, which is where the domains stop being an abstraction.

2.1.5 - The Same Circuit in Two Domains

Every result so far has been an envelope, and the pages have said that an envelope is not a waveform without showing what the difference costs. This tutorial runs the same circuit both ways and puts the two on one axis.

The circuit is the RL branch from adding dynamics, unchanged.

Building the same circuit twice

Only the namespaces differ. dpsimpy.emt.ph1 instead of dpsimpy.dp.ph1, and Domain.EMT instead of Domain.DP:

def build(domain, ns, ph1, name, dt):
    gnd = ns.SimNode.gnd
    n1 = ns.SimNode("n1")
    n2 = ns.SimNode("n2")

    src = ph1.VoltageSource("src")
    src.set_parameters(
        V_ref=complex(100, 0),
        f_src=(50.0 if domain == dpsimpy.Domain.EMT else 0.0),
    )

    r = ph1.Resistor("r"); r.set_parameters(R=10.0)
    l = ph1.Inductor("l"); l.set_parameters(L=0.05)

    src.connect([gnd, n1]); r.connect([n1, n2]); l.connect([n2, gnd])

    system = dpsimpy.SystemTopology(50, [gnd, n1, n2], [src, r, l])
    logger = dpsimpy.Logger(name)
    logger.log_attribute("i_l", "i_intf", l)

    sim = dpsimpy.Simulation(name)
    sim.set_domain(domain); sim.set_system(system)
    sim.set_time_step(dt); sim.set_final_time(0.06)
    sim.add_logger(logger)
    sim.run()
    return rt.read_timeseries_dpsim("logs/%s.csv" % name)

emt = build(dpsimpy.Domain.EMT, dpsimpy.emt, dpsimpy.emt.ph1, "cmp_emt", 5e-5)["i_l"]
dp = build(dpsimpy.Domain.DP, dpsimpy.dp, dpsimpy.dp.ph1, "cmp_dp", 1e-3)
dp_shift = ts.frequency_shift_list(dp, 50)["i_l_shift"]

The one parameter that means different things

f_src is not the same quantity in the two domains, and this is the single easiest way to get a wrong answer here.

The comparison

The dashed line is the dynamic phasor result shifted back onto the 50 Hz carrier. It lies on the EMT waveform. The third curve is the envelope magnitude itself, which is what the DP simulation actually computed: the smooth rise to 5.37 A that the oscillation is riding on.

Twenty times fewer steps, and the same answer to four significant figures.

Why this works, and when it does not

The saving is not because the envelope model is coarser. For a single carrier the transform is exact. The 50 Hz oscillation has been moved out of the integrated quantity and into a coefficient handled analytically, so the step size is set by how fast the envelope changes rather than by the carrier. Here the envelope settles with a 5 ms time constant, and 1 ms resolves it comfortably.

What an envelope domain cannot represent is content outside the band it retains around the carrier. A harmonic, a fast switching transient, or a wideband disturbance is simply absent. That is the trade, and it is the reason both domains exist rather than one being better.

The transform itself is derived under dynamic phasors.

Reading it back

frequency_shift_list appends _shift to every key, so the shifted series is i_l_shift and not i_l. Asking for the original name after shifting raises a KeyError that reads like a missing signal.

The script

The complete script for this page is 05_comparing_domains.py under examples/Python/Tutorials. The numbers quoted above are the numbers it prints, so if the two ever disagree the page is the one that is wrong.

Next

The circuits so far have been passive. Next is a synchronous machine, where the model order becomes a choice and initialization from a powerflow stops being optional.

2.1.6 - Adding a Machine

Everything so far has been passive. A synchronous machine brings two things that no previous tutorial needed: it has mechanical state, so it can swing, and it has to be told the operating point it starts from rather than deducing it.

The network is a machine feeding a strong grid through a line, with a fault applied at the machine terminal for 100 ms.

Machine parameters

A machine is specified by operational parameters rather than winding data:

gen = dpsimpy.dp.ph1.SynchronGenerator4OrderVBR("gen")
gen.set_operational_parameters_per_unit(
    nom_power=555e6, nom_voltage=24e3, nom_frequency=60.0, H=3.7,
    Ld=1.81, Lq=1.76, L0=0.15,
    Ld_t=0.3, Lq_t=0.65, Td0_t=8.0, Tq0_t=1.0,
)

The inductances are in per unit on the machine’s own base and the time constants in seconds. H is the inertia constant, and it sets how fast the machine can accelerate: a low H swings further for the same disturbance.

Each model order takes a different parameter set, and the difference is exactly the states it keeps. The third order model omits Lq_t and Tq0_t because it has no q-axis rotor state at all. The sixth order model adds the subtransient set, Ld_s, Lq_s, Td0_s, Tq0_s and Taa. Passing the wrong set raises a TypeError listing the accepted signature. The equations are derived under reduced order machine models.

The machine base and the network around it

The machine parameters are per unit on the machine’s own base, while the line is given in ohms, so the two only make sense together. The base impedance follows from the machine rating,

A line reactance of a few tenths of an ohm is therefore a few tenths per unit, which is an ordinary transmission connection. The same line specified as 20 mH would be 7.5 Ω, above 7 per unit, and no machine delivers rated power through that.

Initializing the machine

The network is initialized from a powerflow exactly as in the two-bus tutorial. The machine additionally needs its own operating point:

system_dp.init_with_powerflow(systemPF=system_pf, domain=dpsimpy.Domain.DP)

vterm = n1.initial_single_voltage()          # from the powerflow, magnitude and angle
gen.set_initial_values(
    init_complex_electrical_power=complex(300e6, 0),
    init_mechanical_power=300e6,
    init_complex_terminal_voltage=vterm,
)

The difference is measurable. With the angle assumed zero, the rotor speed oscillates by ±0.33% for the first half second, and a fault applied at 0.5 s lands on top of a transient that has nothing to do with it. Taking the voltage from the node, the pre-fault speed is flat to 5 × 10⁻⁶ pu and the only thing in the result is the fault.

That the mechanical power equals the scheduled electrical power is the other half of the same condition: if they disagree, the machine accelerates or decelerates from the first step.

What the model order changes

The same fault, applied to the same machine at the same instant, with three model orders:

All three peak at the clearing instant and all three recover, but the third order model swings noticeably further. That is not numerical: omitting the q-axis rotor removes damping that is physically present, so the third order machine is optimistic about how far it swings and pessimistic about how well it settles.

The practical reading is that model order is a statement about which phenomena you intend to capture. For a first-swing stability question the fourth order model is the usual choice. The subtransient orders matter when the first cycles after the fault are the subject rather than the envelope of the swing.

The script

The complete script for this page is 06_a_machine.py under examples/Python/Tutorials. The numbers quoted above are the numbers it prints, so if the two ever disagree the page is the one that is wrong.

Next

The machine is a source of energy with its own dynamics. Next is a converter, where the dynamics are in the control rather than in a rotor.

2.1.7 - Exchanging Data With Another Tool

Everything so far has ended with a CSV read after the run. This tutorial hands a value out while the simulation runs, which is what co-simulation, hardware in the loop and any live coupling are built on.

The far side here is a file, so nothing external has to be running. A file is a poor co-simulation partner and an excellent first one: the mechanism is identical to an MQTT broker or an FPGA, and only the configuration changes.

What is different about this run

sim = dpsimpy.RealTimeSimulation(name)
...
sim.run(1)

It is a RealTimeSimulation, not a Simulation, and run takes a start delay in seconds. An exchange is paced by the wall clock rather than by how fast the solver can go, because the other side is a real thing running in real time. A one-second simulation takes one second.

That also means the results carry wall-clock timestamps rather than simulation time, which is visible in the output below.

Configuring the far side

interface_config = {
    "type": "file",
    "format": "csv",
    "uri": "logs/exchanged.csv",
    "out": {"flush": True},
}

interface = dpsimpyvillas.InterfaceVillas(
    name="dpsim-file", config=json.dumps(interface_config)
)

This dictionary is not DPsim configuration. It is a VILLASnode node description, passed through as JSON, and its keys are documented by VILLASnode rather than here. Changing type from file to mqtt and giving a broker address is the entire difference between writing to disk and publishing to a broker; nothing in the simulation changes.

Choosing what crosses the boundary

interface.export_attribute(boundary.attr("i_intf").derive_coeff(0, 0), 0)
sim.add_interface(interface)

Two things are worth reading slowly.

The exported quantity is an attribute, the same unit the logger works in. attr("i_intf") takes the whole interface current, which is a matrix, and derive_coeff(0, 0) selects one element of it. Without that you would be handing across a matrix where the far side expects a number.

What comes out

The file is written in VILLASnode’s sample format, not as a DPsim result CSV:

# secs,nsecs,offset,sequence,signal0
1785540392,258011384,nan,0,5.00000000000000000+0.00000000000000000i

A wall-clock timestamp in seconds and nanoseconds, an offset, a sequence number, then one column per exported signal. Over one second at a 10 ms step this run wrote 101 rows.

The offset column is nan here, and that is correct rather than a misconfiguration. It reports the delay between when a sample was created and when it was received, so it only has a value on an incoming path, where those two instants genuinely differ and the difference is the transport latency. These samples originate in the simulation and go straight out, so there is no receive event to measure against and the column is empty by construction.

The sequence number is what the far side uses to detect a missed sample. The queueless interface relies on it directly, which is why it reserves the first imported signal for a sequence counter.

Reading in the other direction

import_attribute is the mirror of export_attribute and makes an incoming value drive something in the simulation, typically the reference of a controlled source. Two options change the timing: blockOnRead halts the simulation at the start of every step until a new value arrives, and syncOnSimulationStart holds the whole run until the far side has produced its first value. Both are described under co-simulation.

Without either, the simulation reads whatever arrived most recently and carries on, which is the right behaviour when the far side is slower and the wrong one when the exchange must be lock-step.

The script

The complete script for this page is 08_exchanging_data.py under examples/Python/Tutorials. It needs a build with VILLASnode available.

Where to read further

DPsim documents its own side of the boundary: which attributes cross, and when they are read and written relative to the step. Everything on the other side of that JSON belongs to VILLASnode and is documented there. The co-simulation page collects the links.

The theory of what a delay across a coupling costs is under the ideal transformer model, which is the same argument whether the two sides are two solvers or two machines.

3 - Developer Guide

These pages describe how DPsim is built rather than what it computes. They are the background for adding a component, changing a solver, or working out why a simulation behaves as it does. For the physics and the numerical methods, see concepts.

Two ideas run through the codebase and are worth reading first. Attributes expose component parameters and state to the logger, the Python bindings and the scheduler. Tasks carry declared attribute dependencies, and those declarations are what the scheduler uses to order and parallelise a timestep.

3.1 - Architecture and Conventions

Start here. These pages describe the shape of the library, the modules it divides into and the class hierarchy underneath them, together with the conventions any change is expected to follow.

3.1.1 - Architecture

Modules and Dependencies

The figure below shows the main components of the DPsim library and their dependencies on other software projects. All functionality is implemented in the C++ core, which can be used standalone or together with the Python interface. The Python interface is a thin wrapper of the C++ core. Jupyter notebooks can either use the DPsim Python interface to run simulations or call executables implemented in C++. The data analysis and plotting is always done in Python using common libraries like Matplotlib. To collect the simulation results from within Python, one can use the villas-dataprocessing Python package.

Another approach to get data in or out of DPsim is the VILLASnode interface, which does not depend on Python at all. The main purpose of the VILLASnode interface is to exchange data during the simulation runtime, for example, in real-time simulation experiments. The data could be send to other simulators, hardware or other software components like databases. Storing the data in databases can be another way of managing (also offline) simulation results if the Python CSV method is not desireable.

The CIM reader is based on the CIM++ library and provides a comfortable alternative to defining the grid manually in C++ or Python. In principle, it calls the same functions to create elements, which are also used in the C++ defined example scenarios, but automatically. DPsim also provides a way to visualize the defined networks before simulation.

The main solver of DPsim is currently the MNA solver because it enables a rather deterministic computation time per simulation time step, which is necessary for real-time simulation. Apart from that, it is also well established in offline circuit simulation. The only dependency of the MNA solver is the linear algebra library Eigen. For some component models, it is possible to use the Sundials ODE solver in combination with the MNA solver. In that case, the component is solved by the ODE solver whereas the network is still handled by the MNA solver. A DAE solver is currently under development. Its main purpose will be offline simulation, for example, to provide reference results where simulation runtime and real-time execution are not relevant.

The component models depend mostly on the Eigen library. Even if components are used in combination with Sundials ODE / DAE solvers, we try to keep the specific functions required by these solvers independent of the Sundials package.

Class Hierarchy

The Simulation class holds references to instances of Interface, Solver, Logger and SystemTopology. For a simulation scenario, the minimum description would include a SystemTopology and a solver type. The Solver instance is then created by the Simulation.

An important function of the Simulation is to collect all tasks, which have to be executed during the simulation. These tasks include computation steps of the individual power system component models as well as read and write tasks of the interfaces and logging variables etc. Before the scheduling is done, Simulation calls getTasks() to retrieve the tasks from three of these classes: Solver, Interface and Logger. The power system component and signal-model tasks are collected by the Solver instances and relayed to the Simulation, while interfaces and loggers contribute their own tasks directly.

All power system element classes inherit from the IdentifiedObject class. This class corresponds with the IdentifiedObject of the IEC61970 CIM and has a uid and name attribute as well.

The next layer of specialization includes information on the topological connection between network elements. An electrical bus and network nodes in general are represented by the TopologiclaNode class. The connection of electrical components, TopologicalPowerComp, is managed via terminals of type TopologicalTerminal. These three types describe the electrical connections of the network, which are bidirectional and include voltages and currents. The signal type elements, TopologicalSignalComp, can only have unidirectional components, which are not expressed using node and terminals. Instead, the attribute system is used to define signal type connections.

3.1.2 - Build

All builds start from a checkout of the repository. To build and read the code, cloning over HTTPS needs no account:

git clone https://github.com/sogno-platform/dpsim.git
cd dpsim

If you intend to contribute, clone your own fork over SSH instead, since contributions are accepted from forks only and pushing needs an authenticated remote:

git clone git@github.com:<your-user>/dpsim.git
cd dpsim
git remote add upstream https://github.com/sogno-platform/dpsim.git

The container route below is the most reproducible, because the image already carries every dependency at the version CI uses. The native routes need those dependencies installed by hand.

Container based

The commands below use docker, but the images are ordinary OCI images, so podman works as a drop-in replacement throughout. On Fedora and Rocky, podman is usually the one already installed. Substitute podman for docker in every command if you prefer it.

The repository ships a development image with all required dependencies:

docker build -t sogno/dpsim:dev -f packaging/Docker/Dockerfile.dev .

Alternatively, pull the prebuilt image instead of building it:

docker pull sogno/dpsim:dev

Then start an interactive session with the working copy mounted into the container:

docker run -it -p 8888:8888 -v $(pwd):/dpsim --privileged sogno/dpsim:dev bash

The -p option maps port 8888 so a JupyterLab instance inside the container is reachable from the host. The --privileged option is required for debug builds. On Windows, the current directory is spelled differently:

docker run -it -p 8888:8888 -v ${pwd}:/dpsim --privileged sogno/dpsim:dev bash

Inside the container, the C++ and Python libraries build as follows:

cd /dpsim
mkdir build && cd build
cmake ..
cmake --build . --target dpsimpy

Targets that are not built by default have to be named explicitly, for example:

cmake --build . --target dpsimpy dpsimpyvillas

To build everything:

cmake --build .

Optional features are enabled through the CMake options defined in the CMakeLists.txt files, for example:

cmake .. -DWITH_GSL=ON

To use the freshly built Python package without installing it, put both the compiled extension and the pure Python package on the path:

cd /dpsim/build
export PYTHONPATH=$(pwd):$(pwd)/../python/src

This is the setup most contributors work with, since it picks up a rebuild immediately without any reinstall step.

Do not use pip install -e . for this. An editable install only links the pure Python sources; dpsimpy is a compiled extension, so edits to the C++ are not picked up and you keep running whatever binary was built at install time. The failure is silent, since the import still succeeds and simply gives you stale behaviour. Either rebuild and rely on PYTHONPATH as above, or reinstall the package after every C++ change.

To summarise the three ways to get DPsim, in increasing order of involvement: pip install dpsim for a released Linux wheel, a native build plus PYTHONPATH for development, and make install to place a build system wide.

If you develop inside a conda environment, the equivalent is to register the same two directories from within the active environment. This needs conda-build installed:

cd /dpsim/build
conda develop $(pwd) && conda develop $(pwd)/../python/src

Note that this writes into the environment, so it becomes specific to your setup.

To run JupyterLab against it:

cd /dpsim
jupyter lab --ip="0.0.0.0" --allow-root --no-browser

To install DPsim system wide instead:

cd /dpsim/build
sudo make install

CMake for Linux

The authoritative dependency list is whatever the Dockerfiles install, since that is what CI builds against. See packaging/Docker/Dockerfile.dev for the Fedora set, and install-fedora-deps.sh or install-ubuntu-deps.sh for scripts that install them.

Both libcimpp and villas-node are optional. Neither needs to be built from source, though the images do not yet take the same route for both.

libcimpp publishes prebuilt .deb and .rpm packages per CIM version as release assets. The Fedora and Debian images install those directly, while the Rocky image still builds it from source:

# Pick the package matching your distribution and the CIM version you need.
wget https://github.com/sogno-platform/libcimpp/releases/download/release%2Fv2.2.0/libcimpp_CGMES_2.4.15_16FEB2016-2.2.0-Linux.deb
sudo apt-get install -y ./libcimpp_CGMES_2.4.15_16FEB2016-2.2.0-Linux.deb
sudo ldconfig

VILLASnode is served from the package repositories at https://packages.fein-aachen.org, which carry both debian/ and redhat/. Note that the images currently still build it from source, pinned to a specific commit, so the packaged version is the more convenient route for a local build but is not what CI exercises.

Building either from source remains supported, and the deps scripts above do that, which is what you want when you need a specific commit rather than a release.

Sundials is only needed for the DAE solver. If your distribution does not package it, the version CI uses is:

git clone --branch v3.2.1 --recurse-submodules --depth 1 https://github.com/LLNL/sundials.git
mkdir -p sundials/build && cd sundials/build
cmake .. -DCMAKE_BUILD_TYPE=Release
make -j$(nproc) install

Cloning, building and installing then work exactly as in the container section above.

CMake for Windows

Windows is built in CI on windows-latest, so the recipe below mirrors what .github/workflows/build_test_windows.yaml runs. You need Visual Studio with the C++ desktop development workload, CMake and Git for Windows. For Python support, install Python 3 and add it to your PATH. Let CMake pick the default generator rather than naming a Visual Studio version, so the build follows whichever Visual Studio you have.

For the C++ libraries only:

mkdir build
cd build
cmake -DWITH_PYBIND=OFF ..
cmake --build . --target dpsim --target dpsim-models --parallel

For the Python bindings, install pybind11 first:

pip install pybind11[global]
mkdir build
cd build
cmake -DWITH_PYBIND=ON ..
cmake --build . --target dpsimpy --parallel

If CMake rejects the spdlog dependency because of its minimum policy version, add -DCMAKE_POLICY_VERSION_MINIMUM=3.5, which is what CI currently does as a workaround.

The dpsim-villas library is not available on Windows, since it requires VILLASnode, which does not build there. WITH_VILLAS therefore stays off and the dpsimpyvillas target does not exist, so co-simulation examples cannot be built on Windows. The CIM reader is likewise not part of the CI Windows build, as libcimpp is not installed there.

CMake for macOS

macOS is not covered by CI, so treat this as a starting point rather than a supported path. Install the dependencies with Homebrew:

/bin/bash -c "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/HEAD/install.sh)"
brew install gcc git cmake graphviz python3 gsl eigen spdlog

Then build as in the container section. Building on Apple Silicon is known to fail while building libcimpp, see issue #609. Configure with -DWITH_CIM=OFF if you do not need the CIM reader.

Python package

Wheels are produced by cibuildwheel in the publish_to_pypi workflow, currently for manylinux x86_64 and CPython 3.9 through 3.13. To build a source distribution locally:

python3 -m build --sdist

Nix

DPsim can be built using Nix, a declarative package manager for reproducible builds. The following steps require a working single-user or multi-user installation of Nix, but not necessarily NixOS.

DPsim uses the Flakes feature, which has to be enabled:

echo "experimental-features=nix-command flakes" > ~/.config/nix/nix.conf

Building DPsim, including all its dependencies:

nix build github:sogno-platform/dpsim

The result is placed in the result folder of the current directory. For development, a local environment can be set up with:

nix develop github:sogno-platform/dpsim

The Flake reference above can be replaced by a local path such as . when the repository is already checked out.

Documentation

The Python and C++ references are generated by separate CMake targets. Both are also built and published by the documentation workflow on every push to master.

Python

Install Sphinx or use the Docker image, then:

mkdir -p build && cd build
cmake ..
make docs

The result is generated in build/docs/sphinx/html/. Note that this target requires the Python bindings, so it is only available when configured with -DWITH_PYBIND=ON.

C++

Install Doxygen or use the Docker image, then:

mkdir -p build && cd build
cmake ..
make docs_cxx

The result is generated in build/docs/doxygen/html/.

Website

The surrounding website is a Hugo site under docs/hugo. It needs the Hugo version pinned in the documentation workflow, since the theme does not build with arbitrary versions:

cd docs/hugo
npm ci
hugo --minify

3.1.3 - Coding Conventions

Conventions that apply across the codebase. For the process of getting a change merged, see contributing.

This is a summary of general guidelines for the development of DPsim.

Scaling of Voltages and Currents

Voltage quantities are expressed either as phase-to-phase RMS values (denominated as RMS3PH) or as phase-to-ground peak values (denominated as PEAK1PH):

• Initialisation quantities (e.g. initialSingleVoltage of SimPowerComp) as RMS3PH values
• Simulation quantities in both SP and DP domain (e.g. mIntfVoltage of DP::Ph1::PiLine) as RMS3PH values
• Simulation quantities in the EMT domain (e.g. mIntfVoltage of EMT::Ph3::Transformer) as PEAK1PH values

Current quantities are expressed either as RMS or as PEAK values:

• Simulation quantities in both SP and DP domain (e.g. mIntfCurrent of DP::Ph1::PiLine) as RMS values
• Simulation quantities in the EMT domain (e.g. mIntfCurrent of EMT::Ph3::Transformer) as PEAK values

Logging

Debug or trace should be the default log level for information that might be nice to have but not necessary for every simulation case.

Calls to the logger that might occur during simulation must use spdlog macros, like SPDLOG_LOGGER_INFO.

3.2 - Attributes and Scheduling

Attributes are the unit of state in DPsim, and they are not only a way to expose a value: the scheduler builds the execution order from the dependencies that components declare over them. The two subjects are one subject, which is why they sit together.

This is also what logging and the co-simulation interfaces operate on, so a quantity that is not an attribute cannot be recorded or exchanged.

3.2.1 - Attributes

In DPsim, an attribute is a special kind of variable which usually stores a scalar or matrix value used in the simulation. Examples for attributes are the voltage of a node, the reference current of a current source, or the left and right vectors of the MNA matrix system. In general, attributes are instances of the Attribute<T> class, but they are usually stored and accessed through a custom smart pointer of type const AttributeBase::Ptr (which expands to const AttributePointer<AttributeBase>).

Through the template parameter T of the Attribute<T> class, attributes can have different value types, most commonly Real, Complex, Matrix, or MatrixComp. Additionally, attributes can fall into one of two categories: Static attributes have a fixed value which can only be changed explicitly through the attribute’s set-method or through a mutable reference obtained through get. Dynamic attributes on the other hand can dynamically re-compute their value from other attributes every time they are read. This can for example be used to create a scalar attribute of type Real whose value always contains the magnitude of another, different attribute of type Complex.

Any simulation component or class which inherits from IdentifiedObject contains an instance of an AttributeList. This list can be used to store all the attributes present in this component and later access them via a String instead of having to use the member variable directly. For reasons of code clarity and runtime safety, the member variables should still be used whenever possible.

Creating and Storing Attributes

Normally, a new attribute is created by using the create or createDynamic method of an AttributeList object. These two methods will create a new attribute of the given type and insert it into the AttributeList under the given name. After the name, create can take an additional parameter of type T which will be used as the initial value for this attribute. Afterwards, a pointer to the attribute is returned which can then be stored in a component’s member variable. Usually this is done in the component’s constructor in an initialization list:

/// Component class Base::Ph1::PiLine

public:
  // Definition of attributes
  const Attribute<Real>::Ptr mSeriesRes;
  const Attribute<Real>::Ptr mSeriesInd;
  const Attribute<Real>::Ptr mParallelCap;
  const Attribute<Real>::Ptr mParallelCond;

// Component constructor: Initializes the attributes in the initialization list
Base::Ph1::PiLine(CPS::AttributeList::Ptr attributeList) :
  mSeriesRes(attributeList->create<Real>("R_series")),
  mSeriesInd(attributeList->create<Real>("L_series")),
  mParallelCap(attributeList->create<Real>("C_parallel")),
  mParallelCond(attributeList->create<Real>("G_parallel")) { };

When a class has no access to an AttributeList object (for example the Simulation class), attributes can instead be created through the make methods on AttributeStatic<T> and AttributeDynamic<T>:

// Simulation class
Simulation::Simulation(String name, Logger::Level logLevel) :
  mName(AttributeStatic<String>::make(name)),
  mFinalTime(AttributeStatic<Real>::make(0.001)),
  mTimeStep(AttributeStatic<Real>::make(0.001)),
  mSplitSubnets(AttributeStatic<Bool>::make(true)),
  mSteadyStateInit(AttributeStatic<Bool>::make(false)),
  //...
{
  // ...
}

Working with Static Attributes

As stated above, the value of a static attribute can only be changed through the attribute’s set-method or by writing its value through a mutable reference obtained by calling get. This means that the value will not change between consecutive reads. Because of the performance benefits static attributes provide over dynamic attributes, attributes should be static whenever possible.

The value of a static attribute can be read by using the attribute’s get-function (i.e. attr->get) or by applying the * operator on the already dereferenced pointer (i.e. **attr), which is overloaded to also call the get function. Both methods return a mutable reference to the attribute’s value of type T&:

AttributeBase::Ptr attr = AttributeStatic<Real>::make(0.001);
Real read1 = attr->get(); //read1 = 0.001
Real read2 = **attr; //read2 = 0.001
Real& read3 = **attr; //read3 = 0.001

The value of an attribute can be changed by either writing to the mutable reference obtained from get, or by calling the set-method:

AttributeBase::Ptr attr = AttributeStatic<Real>::make(0.001);
Real read1 = **attr; //read1 = 0.001
**attr = 0.002;
Real read2 = **attr; //read2 = 0.002
attr->set(0.003);
Real read3 = **attr; //read3 = 0.003

Working with Dynamic Attributes

In general, dynamic attributes can be accessed via the same get and set-methods described above for static attributes. However, dynamic attributes can additionally have dependencies on other attributes which affect the behavior of these methods. Usually, this is used to dynamically compute the attribute’s value from the value of another attribute. In the simplest case, a dynamic attribute can be set to reference another (static or dynamic) attribute using the setReference-method. After this method has been called, the dynamic attribute’s value will always reflect the value of the attribute it references:

AttributeBase::Ptr attr1 = AttributeStatic<Real>::make(0.001);
AttributeBase::Ptr attr2 = AttributeDynamic<Real>::make();

attr2->setReference(attr1);

Real read1 = **attr2; //read1 = 0.001
**attr1 = 0.002;
Real read2 = **attr2; //read2 = 0.002

When working with references between multiple dynamic attributes, the direction in which the references are defined can be important: References should always be set in such a way that the reference relationships form a one-way chain. Only the last attribute in such a reference chain (which itself does not reference anything) should be modified by external code (i.e. through mutable references or the set-method). This ensures that changes are always reflected in all attributes in the chain. For example, the following setup might lead to errors because it overwrites an existing reference:

// Overwriting an existing reference relationship
AttributeBase::Ptr A = AttributeDynamic<Real>::make();
AttributeBase::Ptr B = AttributeDynamic<Real>::make();
AttributeBase::Ptr C = AttributeDynamic<Real>::make();

B->setReference(A); // Current chain: B -> A
B->setReference(C); // Current chain: B -> C, reference on A is overwritten

**C = 0.1; // Change will not be reflected in A

Correct implementation:

AttributeBase::Ptr A = AttributeDynamic<Real>::make();
AttributeBase::Ptr B = AttributeDynamic<Real>::make();
AttributeBase::Ptr C = AttributeDynamic<Real>::make();

B->setReference(A); // Current chain: B -> A
C->setReference(B); // Current chain: C -> B -> A

**A = 0.1; // Updating the last attribute in the chain will update A, B, and C

Aside from setting references, it is also possible to completely recompute a dynamic attribute’s value every time it is read. This can for example be used to create attributes which reference a single matrix coefficient of another attribute, or which represent the magnitude or phase of a complex attribute. Dynamic attributes which depend on one other attribute in this way are also called derived attributes, and they can be created by calling one of the various derive... methods on the original attribute:

AttributeBase::Ptr attr1 = AttributeStatic<Complex>::make(Complex(3, 4));
AttributeBase::Ptr attr2 = attr1->deriveMag();

Real read1 = **attr2; // read1 = 5
**attr1 = Complex(1, 0);
Real read2 = **attr2; // read2 = 1

There is also a general derive-method which can take a custom getter and setter lambda function for computing the derived attribute from its dependency. For more complex cases involving dependencies on multiple attributes, the AttributeDynamic class has a method called addTask which can be used to add arbitrary computation tasks which are executed when the attribute is read or written to. For more information, check the method comments in Attribute.h.

Using Attributes for Logging and Interfacing

When setting up a simulation, there are some methods which require an instance of AttributeBase::Ptr as a parameter. Examples for this are the logger methods (e.g. DataLogger::logAttribute) and interface methods (e.g. InterfaceVillas::exportAttribute). To obtain the required attribute pointer, one can either directly access the public member variables of the component the attribute belongs to, or use the component’s attribute(String name) method which will look up the attribute in the component’s AttributeList:

auto r1 = DP::Ph1::Resistor::make("r_1");
r1->setParameters(5);

auto logger = DataLogger::make("simName");
// Access the attribute through the member variable
logger->logAttribute("i12", r1->mIntfCurrent);

auto intf = std::make_shared<InterfaceVillas>(config);
// Access the attribute through the AttributeList
intf->exportAttribute(r1->attribute('i_intf'), 0, true, true);

// Access the attribute through the member variable and use deriveCoeff to convert it to a scalar value
intf->exportAttribute(r1->mIntfVoltage->deriveCoeff<Complex>(0, 0), 0, true);

When creating a simulation in Python, the component’s member variables are usually not accessible, so the attr-method has to be used for all accesses:

# dpsim-mqtt.py
intf = dpsimpyvillas.InterfaceVillas(name='dpsim-mqtt', config=mqtt_config)
intf.import_attribute(evs.attr('V_ref'), 0, True)
intf.export_attribute(r12.attr('i_intf').derive_coeff(0, 0), 0)

Using Attributes to Schedule Tasks

Attributes are also used to determine dependencies of tasks on data, which is information required by the scheduler. For the usual MNAPreStep and MNAPostStep tasks, these dependencies are configured in the mnaAddPreStepDependencies and mnaAddPostStepDependencies methods:

void DP::Ph1::Inductor::mnaAddPostStepDependencies(
    AttributeBase::List &prevStepDependencies, AttributeBase::List &attributeDependencies,
    AttributeBase::List &modifiedAttributes, Attribute<Matrix>::Ptr &leftVector
  ) {
    attributeDependencies.push_back(leftVector);
  modifiedAttributes.push_back(mIntfVoltage);
  modifiedAttributes.push_back(mIntfCurrent);
}

Here, the MNA post step depends on the solution vector of the system, leftVector, and modifies mIntfVoltage and mIntfCurrent. Therefore, this task needs to be scheduled after the system solution that computes leftVector and before tasks that require the voltage and current interface vectors of the inductance, e.g. the task logging these values.

3.2.2 - Attribute Usage Guidelines

This page gives practical rules for deciding when a model variable should be a DPsim attribute. For details on the attribute mechanism itself, see attributes.

Rule of Thumb

Use an attribute if the value must be visible to DPsim infrastructure, for example logging, interfaces, Python access, string-based lookup, or scheduling.

Otherwise, prefer a normal C++ member variable or a local variable.

Quick Decision Checklist

Before adding a new attribute, ask:

1. Does it need to be logged?
2. Does it need to be imported or exported?
3. Does it need to be accessed from Python or by name?
4. Is it used as a scheduler dependency?
5. Is it an externally relevant model input, output, state, or setpoint?
6. Is a normal C++ variable insufficient?

If the answer to all questions is no, do not make it an attribute.

Use an Attribute For

Use an attribute if the value:

• should be logged
• should be imported or exported through an interface
• should be accessed from Python or generic code by name
• is read or modified by scheduled tasks
• is an externally relevant model input, output, state, or setpoint
• is a derived view of another attribute, for example one matrix coefficient

Typical examples are interface voltages and currents, source references, controller setpoints, and values exchanged through VILLASnode.

Do Not Use an Attribute For

Prefer a normal C++ variable if the value:

• is only used inside one method
• is a temporary intermediate result
• is a cached coefficient or solver helper
• is a fixed implementation detail
• duplicates another existing attribute
• never needs logging, interface access, Python access, or scheduling

Do not create attributes for every variable in the model equations.

Choose the Simplest Attribute Type

If decided that a value should be an attribute, choose the simplest suitable attribute type.

Prefer a static attribute when the value is stored directly by the component:

const Attribute<Real>::Ptr mPower;

MyComponent::MyComponent(const String& name)
    : IdentifiedObject(name),
      mPower(mAttributes->create<Real>("P", 0.0)) {}

Use dynamic, referenced, or derived attributes only when the value must depend on another attribute.

For example, use a derived attribute when exporting one coefficient of a matrix or vector attribute:

intf->exportAttribute(component->mIntfCurrent->deriveCoeff<Complex>(0, 0), 0, true);

Avoid long chains of dynamic, referenced, or derived attributes unless they are really needed.

Access Attributes in C++

In model code, prefer typed attribute members:

mPower->set(power);
const Real power = **mPower;

Use string-based lookup mainly in generic code, logging, interfaces, tests, or Python-style access:

logger->logAttribute("P", component->mPower);
intf->exportAttribute(component->attribute("P"), 0, true);

When assigning a new value, prefer set() if update tasks should be triggered or if the assignment should be explicit. Direct mutable access through **attribute can be used for simple static attributes, but it does not express this intent as clearly.

Examples

Private member variable: no need to use an attribute

This value is stored as a member because it is used by several functions of the class. It is still internal to the implementation: it does not need to be logged, imported or exported, used by the scheduler, or accessed by name.

class MyComponent : public IdentifiedObject {
private:
  Real mConductance = 0.0; // used internally by several methods
};

Externally visible output: use an attribute

This value is a model output. It may be useful for logging, plotting, interfaces, tests, or Python access, so it should be registered as an attribute.

class MyComponent : public IdentifiedObject {
public:
  const Attribute<Real>::Ptr mPower;

  explicit MyComponent(const String& name)
      : IdentifiedObject(name),
        mPower(mAttributes->create<Real>("P", 0.0)) {}

  void updatePower(Real power) {
    mPower->set(power);
  }
};

Runtime setpoint: use an attribute

This value is a model input or setpoint. It may be changed from outside the component, for example through Python, an interface, or a test setup.

class MySource : public IdentifiedObject {
public:
  const Attribute<Real>::Ptr mVoltageRef;

  explicit MySource(const String& name)
      : IdentifiedObject(name),
        mVoltageRef(mAttributes->create<Real>("V_ref", 0.0)) {}

  void setParameters(Real voltageRef) {
    mVoltageRef->set(voltageRef);
  }
};

Derived scalar value: use a derived attribute

This value is not stored separately. It is a scalar view of an existing vector or matrix attribute, which avoids duplicating data manually.

intf->exportAttribute(component->mIntfCurrent->deriveCoeff<Complex>(0, 0), 0, true);

3.2.3 - Interfaces

Why you would exchange data at all, and what it costs, is under co-simulation. This page is how it is configured and what it does to the task graph.

Choosing and configuring an interface

The VILLASnode interface is designed to make use of the various node types and protocols supported by the VILLASframework. By utilizing the nodes provided by VILLASnode, it can be configured to import and export attributes using a wide range of protocols. There are two interface implementations for VILLASnode: InterfaceVillas, which is queued, and InterfaceVillasQueueless.

To create and configure one of the VILLASnode interface instance, create a new shared pointer of type InterfaceVillas or InterfaceVillasQueueless and supply it with a configuration string in the first constructor argument. This configuration must be a valid JSON object containing the settings for the VILLASnode node type that should be used for data import and export. This means that the JSON contains a type key describing what node type to use, as well as any additional configuration options required for this node type. The valid configuration keys can be found in the VILLASnode documentation.

After the object is created, the exportAttribute and importAttribute methods can be used to set up the data exchange between the DPsim simulation and the configured node. The attributes given as the first parameter to these methods are attributes belonging to components in the simulation which should be read or updated by the interface. As an example, for exporting and importing attributes via the MQTT protocol, the VILLASnode interfaces can be configured as follows:

Using C++:

// JSON configuration adhering to the VILLASnode documentation
std::string mqttConfig = R"STRING({
    "type": "mqtt",
    "format": "json",
    "host": "mqtt",
    "in": {
        "subscribe": "/mqtt-dpsim"
    },
    "out": {
        "publish": "/dpsim-mqtt"
    }
})STRING";

// Creating a new InterfaceVillas object
std::shared_ptr<InterfaceVillas> intf = std::make_shared<InterfaceVillas>(mqttConfig);

// Configuring the InterfaceVillas to import and export attributes
intf->importAttribute(evs->mVoltageRef, 0, true, true);
intf->exportAttribute(r12->mIntfCurrent->deriveCoeff<Complex>(0, 0), 1, true, "v_load");

Using Python:

# JSON configuration adhering to the VILLASnode documentation
mqtt_config = '''{
        "type": "mqtt",
        "format": "json",
        "host": "mqtt",
        "in": {
            "subscribe": "/mqtt-dpsim"
        },
        "out": {
            "publish": "/dpsim-mqtt"
        }
}'''

# Creating a new InterfaceVillas object
intf = dpsimpyvillas.InterfaceVillas(name='dpsim-mqtt', config=mqtt_config)

# Configuring the InterfaceVillas to import and export attributes
intf.import_attribute(evs.attr('V_ref'), 0, True)
intf.export_attribute(r12.attr('i_intf').derive_coeff(0, 0), 0)

Adding an Interface to the Simulation

After a new interface has been created and configured, it can be added to a simulation using the Simulation::addInterface method:

// Create and configure simulation
RealTimeSimulation sim(simName);
sim.setSystem(sys);
sim.setTimeStep(timeStep);
sim.setFinalTime(10.0);

// Create and configure interface
auto intf = //...

// Add interface to simulation
sim.addInterface(intf);

Adding an interface also adds two tasks to the simulation, one before the step and one after, so an imported value is in place before anything reads it and an exported one is sent after everything that could change it. The transfer itself happens on separate threads, so a slow far side does not hold up the solver. How that is arranged, and why it matters, is under how an interface is scheduled.

Synchronizing the Simulation with the Environment

To allow for synchronizing the DPsim simulation with external services, the Interface class provides some additional configuration options in the importAttribute and exportAttribute methods. For imports, setting the blockOnRead parameter will completely halt the simulation at the start of every time step until a new value for this attribute was read from the environment. Additionally, the syncOnSimulationStart parameter can be set for every import to indicate that this attribute is used to synchronize the start of the simulation. When a simulation contains any interfaces importing attributes which have syncOnSimulationStart set, the Simulation::sync will be called before the first time step. This method will:

• write out all attributes configured for export to the environment
• block until all attributes with syncOnSimulationStart set have been read from the environment at least once
• write out all exported attributes again

Note that this setting operates independently of the blockOnRead flag. This means that with both flags set, the simulation will block again after the synchronization at the start of the first time step until another value is received for the attribute in question.

The two tasks

Adding an interface adds a PreStep and a PostStep task.

PreStep is declared to modify every attribute imported from the environment, so the scheduler places it before any task that depends on those attributes. An imported value is therefore in place before anything reads it.

PostStep is declared to depend on every attribute exported to the environment, so it runs after anything that might modify them.

Why PostStep declares a modified attribute

PostStep modifies nothing in the simulation: it only sends values outward. The scheduler prunes tasks whose outputs nothing needs, so a task that modifies nothing is dropped, and the export would silently never happen.

To prevent that, PostStep is declared to modify Scheduler::external. That attribute exists to make a task reachable when its real effect is outside the simulation.

Task execution is not the moment of transfer

When these tasks execute is not when the data actually crosses the boundary. The interface spawns a reader thread and a writer thread and communicates with them over a lock-free queue.

The consequence is the useful part: a slow import or export does not block the solver. The simulation hands a value to the queue and continues. That is what makes an interface to a slow or unreliable far side usable at all, and it is also why a value read this step may have been produced some time ago.

Blocking is opt-in through blockOnRead and syncOnSimulationStart on the import, described on the co-simulation page. Those are the only ways the exchange paces the simulation.

3.2.4 - Task Scheduling

Within each simulation timestep, DPsim executes a set of tasks: discrete units of computation contributed by components, the solver, interfaces, and loggers. Before the first timestep the scheduler collects all tasks, resolves their data dependencies into a directed acyclic graph, and produces an ordered schedule. That schedule is then replayed on every timestep with no further graph analysis.

Tasks

The Task base class

Every task is an instance of a class that inherits from CPS::Task (dpsim-models/include/dpsim-models/Task.h). Each subclass implements one member function:

virtual void execute(Real time, Int timeStepCount) = 0;

To participate in scheduling, a task declares its data dependencies through three attribute lists that are populated in the task’s constructor:

All three lists hold AttributeBase::Ptr objects, the same pointers used throughout the component model. See Attributes for details on the attribute system.

The component text logger (CPS::Logger, backed by spdlog) is a separate mechanism used for human-readable debug and diagnostic output. It is not part of the scheduling system and can print any value regardless of whether it is an attribute.

For practical rules on when a variable should be an attribute versus a plain member variable, see Attribute Usage Guidelines.

Common component task conventions

The names below are component and solver conventions, not scheduler-level concepts. The scheduler only sees the attribute dependencies a task declares; it has no notion of a “PreStep” or “PostStep” and never orders tasks by these names.

MNA components typically define two task classes per component:

This is a common pattern rather than a fixed rule; the exact work each task does is component-specific. Signal-domain components (regulators, governors, control blocks) define their own task list via getTasks(); many separate previous-step state handling from output updates, for example a PreStep that copies state from the previous step and a Step that updates the block outputs.

The solver itself contributes a task that solves the MNA system; individual components do not depend on it by name, they depend on leftVector instead (see below).

Building the schedule

Task collection

Simulation::prepSchedule() collects all tasks before the first timestep from three top-level sources:

• Solvers: each solver contributes its task list via Solver::getTasks(). For MNA solvers this list bundles:
  • the matrix-solve task,
  • MNA component pre-/post-step tasks from MNASimPowerComp::mnaTasks() (built during solver initialization via mnaAddPreStepDependencies() / mnaAddPostStepDependencies()),
  • signal-domain component tasks returned by SimSignalComp::getTasks(),
  • optional solver-side tasks, such as state-space extraction, when enabled.
• Interfaces: each interface contributes its own tasks via Interface::getTasks(). These typically depend on the attributes exchanged with external systems.
• Loggers: each logger contributes a logging task via Logger::getTask(), depending on the logged attributes so values are written after the producing tasks have run.

All tasks are placed in a flat Task::List and handed to the scheduler.

Dependency resolution

Scheduler::resolveDeps() (dpsim/src/Scheduler.cpp) translates the attribute-level declarations into directed edges between tasks. For every attribute in mModifiedAttributes, it finds all tasks that list that attribute in their mAttributeDependencies and adds an edge:

Task A modifies attr_X and task B depends on it, so the edge runs A to B and the scheduler must place A first.

A special Root sentinel task is inserted as a sink for all mPrevStepDependencies entries. Its role is explained in the pruning step below.

Topological sort and pruning

Scheduler::topologicalSort() first runs a backward breadth-first search (BFS) from Root, marking every task that transitively contributes to a simulation output. Tasks not reachable in this pass are dropped from the schedule because they produce data no downstream consumer reads in the current timestep.

Kahn’s algorithm then processes the remaining tasks in dependency order and appends them to the schedule. The result is a flat, ordered list in which every task appears after all of its current-step predecessors.

The Root sentinel matters here: it holds a reference to an external attribute updated by an interface or by the solver, so the backward BFS reaches it and keeps every task that writes previous-timestep state, even when that output is only consumed in the next timestep.

Level scheduling

For parallel execution the ordered list is converted into levels by Scheduler::levelSchedule(). Each task is assigned to the level one greater than the highest-level task it depends on:

Tasks within the same level have no data dependencies between them and can execute in parallel. The scheduler guarantees that all tasks in level k finish before any task in level k+1 starts.

Scheduler variants

The scheduler is chosen at Simulation construction time; SequentialScheduler is the default.

Per-timestep execution

Scheduler::step(time, timeStepCount) is called once per timestep. For the sequential scheduler:

for (auto& task : mSchedule)
    task->execute(time, timeStepCount);

Parallel schedulers distribute tasks across threads within each level and synchronize with a barrier before advancing to the next level.

This page describes how the scheduler works. For how to give a component tasks of its own, see adding tasks to a component.

3.2.5 - Adding Tasks to a Component

How to attach tasks to a component. For how the scheduler consumes them, see scheduling.

Signal components

Signal components inherit from SimSignalComp and return their tasks from getTasks(). The usual pattern is to define inner Task classes whose constructors populate the dependency lists, then instantiate them in getTasks():

class MyComponent : public SimSignalComp {
public:
    const Attribute<Real>::Ptr mInput;       // written by upstream component
    const Attribute<Real>::Ptr mOutput;      // read by downstream component
    const Attribute<Real>::Ptr mOutputPrev;  // state carried across timesteps

    class PreStep : public Task {
    public:
        explicit PreStep(MyComponent& comp)
            : Task(**comp.mName + ".PreStep"), mComp(comp) {
            mPrevStepDependencies.push_back(mComp.mOutput);
            mModifiedAttributes.push_back(mComp.mOutputPrev);
        }
        void execute(Real time, Int timeStepCount) override {
            **mComp.mOutputPrev = **mComp.mOutput;
        }
    private:
        MyComponent& mComp;
    };

    class Step : public Task {
    public:
        explicit Step(MyComponent& comp)
            : Task(**comp.mName + ".Step"), mComp(comp) {
            mAttributeDependencies.push_back(mComp.mInput);
            mModifiedAttributes.push_back(mComp.mOutput);
        }
        void execute(Real time, Int timeStepCount) override {
            mComp.signalStep(time, timeStepCount);
        }
    private:
        MyComponent& mComp;
    };

    Task::List getTasks() override {
        return { std::make_shared<PreStep>(*this),
                 std::make_shared<Step>(*this) };
    }
};

PreStep uses mPrevStepDependencies for mOutput because it reads the value produced last timestep, not the value that Step will produce this timestep. Using mAttributeDependencies here would create a same-step dependency on Step and force PreStep after Step, which is backwards.

MNA power components

MNA components inherit from MNASimPowerComp<VarType>. Instead of getTasks(), they implement two hook functions that MNASimPowerComp calls when it builds the MnaPreStep and MnaPostStep tasks during solver initialization.

void DP::Ph1::MyComponent::mnaAddPreStepDependencies(
    AttributeBase::List& prevStepDependencies,
    AttributeBase::List& attributeDependencies,
    AttributeBase::List& modifiedAttributes) {

    prevStepDependencies.push_back(mIntfCurrent);   // read from previous step
    modifiedAttributes.push_back(mRightVector);     // stamp right-hand side
}

void DP::Ph1::MyComponent::mnaAddPostStepDependencies(
    AttributeBase::List& prevStepDependencies,
    AttributeBase::List& attributeDependencies,
    AttributeBase::List& modifiedAttributes,
    Attribute<Matrix>::Ptr& leftVector) {

    attributeDependencies.push_back(leftVector);  // wait for matrix solve
    modifiedAttributes.push_back(mIntfVoltage);
    modifiedAttributes.push_back(mIntfCurrent);
}

PostStep must always list leftVector in attributeDependencies. This creates the edge from the solver’s matrix-solve task to every component’s PostStep, ensuring the solution vector is available before voltages and currents are extracted.

Dependency declaration checklist

• Every attribute read inside execute() must appear in mAttributeDependencies or mPrevStepDependencies.
• Every attribute written inside execute() must appear in mModifiedAttributes.
• State carried from the previous timestep goes in mPrevStepDependencies, not mAttributeDependencies.
• MnaPostStep must list leftVector in attributeDependencies.
• No attribute should appear in both mAttributeDependencies and mPrevStepDependencies for the same task.

3.3 - Writing a Model

The path from an empty file to a working component: what to declare, which hooks the solver calls and in what order, how a component built from other components is assembled, and how to debug one that runs but produces the wrong answer.

For the equations a model should implement, see Concepts. For worked examples of finished models, see model implementations.

3.3.1 - Component and Solver Initialization

Initialization is the phase between constructing the system topology and running the first timestep. Its job is to size the system matrices, derive initial state from power-flow results, register MNA tasks, and stamp static conductances. Two constraints drive its structure:

• The system matrix size depends on the total number of simulation nodes, including virtual nodes declared by composite components and their sub-components. All virtual nodes must therefore be known before the matrices are allocated.
• Component parameter values (impedances, initial phasors) depend on terminal voltages and powers, which are only available after a power-flow solve.

These two constraints impose an ordering that is captured in the solver’s initialization sequence.

MNA Solver Initialization Sequence

MnaSolver::initialize() executes the following steps in order.

Step 1 — Identify topology objects

identifyTopologyObjects() iterates over SystemTopology::mComponents and sorts each component into one of four lists:

Ground nodes are excluded here.

Step 2 — Create sub-components (pre-pass)

Before the matrix can be sized, every composite component’s sub-component tree must be fully instantiated so that all virtual nodes are visible. The solver calls createSubComponents() recursively on every MNA component:

• Only sub-components newly registered by this call are recursed into, because eagerly-constructed sub-components (created in the constructor before connect() has run) are not yet safe to recurse into.
• This step is a pre-pass only — it must not set parameter values derived from terminal data or frequency.

For details on the three-stage composite lifecycle (createSubComponents, initializeParentFromNodesAndTerminals, mnaCompInitialize), see Subcomponent Handling.

Step 3 — Collect virtual nodes and assign indices

collectVirtualNodes() visits every component and calls virtualNodes() to collect all virtual SimNode objects, then appends them to the solver’s node list. assignMatrixNodeIndices() then assigns a contiguous integer index to every simulation node (real and virtual), which determines the row/column layout of the system matrices.

After this step the matrix size is fixed.

Step 4 — Allocate empty matrices

createEmptyVectors() and createEmptySystemMatrix() allocate the left-side vector, right-side vector, system matrix (dense or sparse depending on the solver variant), and switch-variant copies. For sparse solvers, mBaseSystemMatrix and mLuFactorizations are also allocated here, with one variant per switch combination.

Step 5 — Initialize components (initializeComponents)

This step has three sub-passes over the component lists.

5a — Power components: initializeFromNodesAndTerminals

For every SimPowerComp<VarType> in mMNAComponents and mMNAIntfVariableComps:

1. checkForUnconnectedTerminals() validates connectivity.
2. If mInitFromNodesAndTerminals is set (the default), initializeFromNodesAndTerminals(mSystem.mSystemFrequency) is called.

This is where components read their terminal voltages and powers and derive physical parameters (impedances, initial phasor values, per-unit quantities). For composite components initializeFromNodesAndTerminals() is final in CompositePowerComp and sequences the three lifecycle stages automatically; non-composite power components override it directly.

5b — Signal components: initialize(omega, timeStep)

Each SimSignalComp in mSimSignalComps receives initialize(mSystem.mSystemOmega, mTimeStep). This is the hook for signal-domain components (regulators, governors, PSS blocks) to allocate their state buffers, set initial values, and wire up attribute connections.

5c — MNA components: mnaInitialize

Each MNA component (including switches) receives mnaInitialize(omega, timeStep, leftVector). In MNASimPowerComp this method:

1. Clears and re-registers MNAPreStep / MNAPostStep tasks according to the hasPreStep / hasPostStep flags.
2. Initializes mRightVector to zero with the correct size.
3. Calls mnaCompInitialize(omega, timeStep, leftVector) on the component.

In mnaCompInitialize, component classes call updateMatrixNodeIndices() and perform any one-time MNA setup that requires the final node layout (e.g. allocating per-component history vectors sized to the system).

Nodes are initialized last via SimNode::initialize(), which zeros the node voltage.

Step 6 — Optional steady-state initialization

If mSteadyStateInit is set, steadyStateInitialization() iterates the MNA solve until the phasor solution converges. The flag mIsInInitialization is set to true for this sub-phase so that components can distinguish initialization solves from simulation solves via mBehaviour (see below).

Step 7 — Set simulation behaviour

After initialization solves are complete, the solver calls setBehaviour(TopologicalPowerComp::Behaviour::MNASimulation) on every TopologicalPowerComp and setBehaviour(SimSignalComp::Behaviour::Simulation) on every SimSignalComp.

The Behaviour enum (defined in TopologicalPowerComp) has three values:

Components that need different behaviour between initialization and simulation check mBehaviour in their pre/post-step methods or in mnaCompPreStep.

Step 8 — Initialize system matrices (initializeSystem)

initializeSystem() selects one of three paths:

• Parallel frequencies (initializeSystemWithParallelFrequencies): stamps each frequency into a separate thread.
• Variable matrix (initializeSystemWithVariableMatrix): used by MnaSolverSysRecomp; saves static switch matrices as base matrices and adds variable elements on top.
• Precomputed matrices (initializeSystemWithPrecomputedMatrices): the common path. Calls switchedMatrixStamp() for each switch combination, which iterates over all static MNA components and calls mnaApplySystemMatrixStamp() and mnaApplyRightSideVectorStamp(). LU factorizations are computed for each variant.

After this step the solver is ready to execute timesteps.

Component Class Hierarchy and Init Hooks

The following diagram shows which initialization methods live in which class, and the override points for component authors.

Methods marked * are the virtual override points for component authors. Methods marked final must not be overridden; the base class sequences them correctly.

Component Method Contracts

The table below summarizes which initialization method has which responsibilities. A tick means the operation belongs in that method; a cross means it must not appear there.

Common pitfalls

• Accessing terminals in the constructor or createSubComponents: terminal data (initial voltage, connected power) is not yet populated. The topology is set up but power-flow has not run.
• Accessing mFrequencies(0,0) in createSubComponents: the system frequency matrix is set on SimPowerComp via initialize(Matrix) which only runs later. Use the frequency argument passed to initializeParentFromNodesAndTerminals or the omega argument in mnaCompInitialize.
• Zero-valued shunt branches: a capacitor or reactor with zero admittance injects a zero row/column into the system matrix, which makes the LU factorization singular. Guard with a strict > 0 check and omit the branch rather than inserting a zero stamp.

Composite Component Initialization Sequence

The following diagram shows how the solver and a composite component interact during initialization. For further details see Subcomponent Handling.

PFSolver Initialization

PFSolver::initialize() follows a simpler sequence because it operates only on single-phase SP components with no sub-component tree and does not need a createSubComponents pre-pass.

PFSolver::setSolverAndComponentBehaviour() is the equivalent of Step 7 for the MNA solver: it calls setBehaviour(Behaviour::PFSimulation) or setBehaviour(Behaviour::Initialization) on all components to allow them to switch stamping modes.

Known Design Issues (issue #59)

The following areas were identified in GitHub issue #59 as needing improvement.

SimPowerComp::initialize(Matrix frequencies) naming clash

SimPowerComp<T>::initialize(Matrix frequencies) is called by the solver to propagate frequency information down the component tree.

• For power components, use initializeFromNodesAndTerminals() or initializeParentFromNodesAndTerminals().
• For signal components, use the initialize(Real omega, Real timeStep) hook provided by SimSignalComp.
• For anything else (e.g. setting up state-space matrices), add a named helper called from one of the above.

The base implementation of SimPowerComp::initialize(Matrix) should be renamed to something that cannot be accidentally overridden (e.g. propagateFrequencies()), and an override guard should be added to catch accidental overrides.

Sub-component construction in constructors

Some components create and register sub-components eagerly in their constructor before connect() has been called on those sub-components. This works today because the solver’s createSubComponents pre-pass skips already-registered sub-components, but it couples topology creation to object construction and makes components harder to reason about. The long-term goal is to migrate all sub-component construction to createSubComponents(), giving a clear rule: the constructor only allocates and the topology stage wires.

Signal component initialize not sequenced with power flow

Signal components receive initialize(omega, timeStep) after initializeFromNodesAndTerminals on power components but before the MNA tasks are registered. If a signal component’s initial state depends on the power-flow solution (e.g. an exciter initializing to match the generator terminal voltage), it must read the relevant attribute values directly — there is no formal mechanism today to express this dependency in the initialization sequence. A future improvement would be to give signal components access to the settled power-flow solution before their initialize is called.

3.3.2 - Real-Time Execution

Why you would run in real time, and how to start such a run, is under real-time simulation. This page is what has to be true of the host and of the models for a deadline to be met.

DPsim runs in real time on any system, but without tuning the smallest reliable step is nowhere near microseconds, because operating system noise and other processes interfere. With the tuning below, steps as low as 5 us synchronised to an FPGA through VILLASnode have been achieved.

Operating System and Kernel

A kernel built with PREEMPT_RT improves latency when issuing system calls and enables the FIFO scheduler that avoids preemption during the run.

This used to mean tracking down an out-of-tree patch set. It no longer does: PREEMPT_RT was merged into the mainline Linux kernel in 6.12, so a recent kernel can be built with it directly and a growing number of distributions ship or package one. Check what you already have before installing anything:

uname -v | grep -q PREEMPT_RT && echo "already real-time" || echo "not a PREEMPT_RT kernel"

If you need one, most distributions still offer a binary package. On Rocky Linux:

sudo dnf --enablerepo=rt install kernel-rt kernel-rt-devel

More aggressive tuning can involve isolating a set of cores for exclusive use by the real-time simulation. This way, the kernel will not schedule any processes on these cores. Add the kernel parameters isolcpus and nohz_full using, for example, grubby:

sudo grubby --update-kernel=ALL --args="isolcpus=9,11,13,15 nohz_full=9,11,13,15"

Something similar, but less invasive and non-permanent can be achieved using tuna:

sudo tuna isolate -c 9,11,13,15

To avoid real-time throttling to cause overruns disable this feature:

sudo bash -c "echo -1 > /proc/sys/kernel/sched_rt_runtime_us"

Note that this is not persistent when rebooting.

Simulation Model Tuning

Real time capable models cannot issue any system calls during simulation as the context switch to the kernel introduces unacceptable latencies. This means models cannot allocate memory, use mutexes or other interrupt-driven synchronization primitives, read or write data from files. You should turn off logging, when time steps in the low milliseconds are desired. There is a RealTimeDataLogger that can be used to output simulation results in these cases. Note however, that this logger pre-allocated the memory required for all of the logging required during simulations. Your machine may run out of memory, when the simulation is long or you log too many signals.

You can increase the performance of your simulation by adding the -flto and -march=native compiler flags:

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 8801cbe8d..4a2843269 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -79,7 +79,7 @@ include(CheckSymbolExists)
 check_symbol_exists(timerfd_create sys/timerfd.h HAVE_TIMERFD)
 check_symbol_exists(getopt_long getopt.h HAVE_GETOPT)
 if(CMAKE_BUILD_TYPE STREQUAL "Release" OR CMAKE_BUILD_TYPE STREQUAL "RelWithDebInfo")
-       add_compile_options(-Ofast)
+       add_compile_options(-Ofast -flto -march=native)
 endif()

 # Get version info and buildid from Git

Where the time step comes from

By default the simulation paces itself against the host clock, which is enough for almost everything. Synchronising the step to an external source instead is only necessary when the accuracy of the step itself matters at the nanosecond level, which in practice means hardware in the loop against equipment with its own clock.

That distinction is worth making before reaching for it: locking to an external source constrains the whole run and is not a general improvement, only the answer to a specific requirement.

Writing a model that can hold a deadline

A real-time capable model must issue no system calls during simulation: the context switch into the kernel costs more than the deadline allows.

Turn logging off when the step is in the low milliseconds. RealTimeDataLogger exists for the cases that still need results: it buffers in memory and writes at the end rather than touching the disk inside the step.

3.3.3 - Interfacing with the MNA Solver

The various solver classes based on MNASolver are used to perform Nodal Analysis during a DPsim simulation. For components to be able to influence the input variables of the MNA, they have to implement certain methods defined in the MNAInterface interface class. While it is possible to individually implement MNAInterface for every component, the behavior of many components can be unified in a common base class. This base class is called MNASimPowerComp<T>. Currently, it is the only class which directly implements MNAInterface and in turn all MNA components inherit from this class. Much like the CompositePowerComp class for Composite Components, the MNASimPowerComp class provides some common behavior for all MNA components, e.g. the creation and registration of the MNAPreStep and MNAPostStep tasks. Additionally, MNASimPowerComp provides a set of virtual methods prefixed mnaComp... which can be implemented by the child component classes to provide their own MNA behavior. These methods are:

virtual void mnaCompInitialize(Real omega, Real timeStep, Attribute<Matrix>::Ptr leftVector);
virtual void mnaCompApplySystemMatrixStamp(SparseMatrixRow& systemMatrix);
virtual void mnaCompApplyRightSideVectorStamp(Matrix& rightVector);
virtual void mnaCompUpdateVoltage(const Matrix& leftVector);
virtual void mnaCompUpdateCurrent(const Matrix& leftVector);
virtual void mnaCompPreStep(Real time, Int timeStepCount);
virtual void mnaCompPostStep(Real time, Int timeStepCount, Attribute<Matrix>::Ptr &leftVector);
virtual void mnaCompAddPreStepDependencies(AttributeBase::List &prevStepDependencies, AttributeBase::List &attributeDependencies, AttributeBase::List &modifiedAttributes);
virtual void mnaCompAddPostStepDependencies(AttributeBase::List &prevStepDependencies, AttributeBase::List &attributeDependencies, AttributeBase::List &modifiedAttributes, Attribute<Matrix>::Ptr &leftVector);
virtual void mnaCompInitializeHarm(Real omega, Real timeStep, std::vector<Attribute<Matrix>::Ptr> leftVector);
virtual void mnaCompApplySystemMatrixStampHarm(SparseMatrixRow& systemMatrix, Int freqIdx);
virtual void mnaCompApplyRightSideVectorStampHarm(Matrix& sourceVector);
virtual void mnaCompApplyRightSideVectorStampHarm(Matrix& sourceVector, Int freqIdx);

MNASimPowerComp provides empty default implementations for all of these methods, so component classes are not forced to implement any of them.

Controlling Common Base Class Behavior

Child component classes can control the behavior of the base class through the constructor arguments of MNASimPowerComp. The two boolean variables hasPreStep and hasPostStep can be used to control whether the MNAPreStep and MNAPostStep tasks will be created and registered. If these tasks are created, the mnaCompPreStep / mnaCompPostStep and mnaCompAddPreStepDependencies / mnaCompAddPostStepDependencies methods will be called during the component’s lifecycle. If the tasks are not created, these methods are superfluous and should not be implemented in the child class.

Currently, the MNASimPowerComp base class only exhibits additional behavior over the mnaComp... methods in the mnaInitialize method. In this method, the list of MNA tasks is cleared, and the new tasks are added according to the hasPreStep and hasPostStep parameters. Additionally, the right vector attribute mRightVector required by MNAInterface is set to a zero-vector with its length equal to that of the system leftVector. If this behavior is not desired, e.g. for resistors which have no influence on the system right vector, the right vector can be re-set to have zero size in the mnaCompInitialize method:

void DP::Ph1::Resistor::mnaCompInitialize(Real omega, Real timeStep, Attribute<Matrix>::Ptr leftVector) {
  updateMatrixNodeIndices();

  **mRightVector = Matrix::Zero(0, 0);
  //...
}

For all other MNA methods, the MNASimPowerComp base class will just call the associated mnaComp... method. For more details, take a look at the implementations in MNASimPowerComp.cpp.

3.3.4 - Subcomponent Handling

In DPsim, there are many components which can be broken down into individual subcomponents. Examples are the PiLine, consisting of an inductor, three resistors, and two capacitors, or the NetworkInjection which contains a voltage source. On the C++ class level, these subcomponents are represented by member variables within the larger component class. In this guide, all components which have subcomponents are called composite components.

Creating Composite Components

While normal components are usually subclasses of SimPowerComp<T> or MNASimPowerComp<T>, there exists a special base class for composite components called CompositePowerComp<T>. This class provides multiple methods and parameters for configuring how the subcomponents should be handled with respect to the MNAPreStep and MNAPostStep tasks. The main idea here is that the subcomponents do not register their own MNA tasks, but instead their MNA methods like mnaPreStep and mnaPostStep are called explicitly in the tasks of the composite component. In the constructor of CompositePowerComp<T>, the parameters hasPreStep and hasPostStep can be set to automatically create and register a MNAPreStep or MNAPostStep task that will call the mnaCompPreStep or mnaCompPostStep method on execution. Additionally, all subcomponents should be registered as soon as they are created using the addMNASubComponent-method. This method takes multiple parameters defining how and in what order the subcomponent’s pre- and post- steps should be called, as well as if the subcomponent should be stamped into the system rightVector.

Initialization lifecycle

Composite components are initialized in three stages, each with a defined role. (These are distinct from the electrical phases A/B/C of a three-phase component.)

1. Topology stage (createSubComponents()). Decides which sub-components exist and how they are wired: make_shared, connect() to network/virtual nodes, and addMNASubComponent(). This runs in a pre-pass before the MNA system matrix is sized, so any virtual nodes owned by sub-components are visible to collectVirtualNodes(). Because it runs before power-flow results or the simulation frequency are guaranteed to be available, createSubComponents() must not read terminal data (initialSingleVoltage(), singleActivePower(), …), system frequency (mFrequencies(0,0)), or compute any power-/impedance-derived value. It must be idempotent — guard the body with mSubCompCreated (a protected field inherited from CompositePowerComp).
2. Parameterization stage (initializeParentFromNodesAndTerminals(Real frequency)). Sets the values the sub-components created in stage 1 will use. This is where terminal reads, frequency-dependent impedance/admittance calculations, and setParameters() calls on sub-components belong. The simulation frequency is passed in as a direct argument, so there is no need to access mFrequencies(0,0). This is the hook concrete composites must implement — do not override initializeFromNodesAndTerminals() directly; the base class owns that method and calls this hook at the right time.
3. MNA-init stage (mnaCompInitialize()). Unchanged; already recurses into sub-components.

CompositePowerComp<VarType>::initializeFromNodesAndTerminals() is final and sequences these stages:

void initializeFromNodesAndTerminals(Real frequency) final {
  createSubComponents();                            // idempotent safety net for paths
                                                     //   that reach this composite without
                                                     //   the solver's pre-pass having run
  initializeParentFromNodesAndTerminals(frequency); // parent derives values,
                                                     //   setParameters() on subs
  for (auto subComp : mSubComponents) {
    subComp->initialize(mFrequencies);              // propagate frequencies down
    subComp->initializeFromNodesAndTerminals(frequency);
  }
}

The loop re-enters this same final wrapper for any sub-component that is itself a composite, so the whole tree initializes correctly without each level manually calling initialize()/initializeFromNodesAndTerminals() on its children.

A sub-component whose very existence (not just its value) depends on a parameterization-stage value — e.g. picking an inductor vs. a capacitor based on the sign of computed reactive power — cannot be registered in createSubComponents(). Create and register it directly inside initializeParentFromNodesAndTerminals() instead. This is safe because the MNA-registered sub-component list is not consumed until MnaSolver::initialize() finishes the parameterization stage for all components. The one constraint: the late-registered sub-component must not introduce new virtual nodes — those must be declared in the constructor or setParameters(), before collectVirtualNodes() runs.

// DP_Ph1_PiLine.cpp
DP::Ph1::PiLine::PiLine(String uid, String name, Logger::Level logLevel)
  : Base::Ph1::PiLine(mAttributes),
  // Call the constructor of CompositePowerComp and enable automatic pre- and post-step creation
  CompositePowerComp<Complex>(uid, name, true, true, logLevel)
{
  //...
}

void DP::Ph1::PiLine::createSubComponents() {
  if (mSubCompCreated)
    return;
  mSubCompCreated = true;

  // Create series sub components
  mSubSeriesResistor = std::make_shared<DP::Ph1::Resistor>(**mName + "_res", mLogLevel);

  // Setup mSubSeriesResistor... (only from values already known from this
  // component's own setParameters()/constructor/Attributes - no terminal or
  // frequency reads here)

  // Register the resistor as a subcomponent. The resistor's pre- and post-step will be called before the pre- and post-step of the parent,
  // and the resistor does not contribute to the `rightVector`.
  addMNASubComponent(mSubSeriesResistor, MNA_SUBCOMP_TASK_ORDER::TASK_BEFORE_PARENT, MNA_SUBCOMP_TASK_ORDER::TASK_BEFORE_PARENT, false);

  mSubSeriesInductor = std::make_shared<DP::Ph1::Inductor>(**mName + "_ind", mLogLevel);

  // Setup mSubSeriesInductor...

  // Register the inductor as a subcomponent. The inductor's pre- and post-step will be called before the pre- and post-step of the parent,
  // and the inductor does contribute to the `rightVector`.
  addMNASubComponent(mSubSeriesInductor, MNA_SUBCOMP_TASK_ORDER::TASK_BEFORE_PARENT, MNA_SUBCOMP_TASK_ORDER::TASK_BEFORE_PARENT, true);
  //...
}

void DP::Ph1::PiLine::initializeParentFromNodesAndTerminals(Real frequency) {
  //...
  // Frequency-dependent values go here, not in createSubComponents().
  Real omega = 2. * PI * frequency;
  Complex impedance = {**mSeriesRes, omega * **mSeriesInd};
  //...
}

Orchestrating MNA Method Calls

By choosing which methods to override in the composite component class, subcomponent handling can either be offloaded to the CompositePowerComp base class or manually implemented in the new component class. By default, CompositePowerComp provides all methods demanded by MNAInterface in such a way that the subcomponents’ MNA-methods are properly called. To also allow for the composite component class to perform further actions in these MNA-methods, there exist multiple methods prefixed with mnaParent, e.g. mnaParentPreStep or mnaParentAddPostStepDependencies. These parent methods will usually be called after the respective method has been called on the subcomponents. For the mnaPreStep and mnaPostStep methods, this behavior can be set explicitly in the addMNASubComponent method.

If a composite component requires a completely custom implementation of some MNA-method, e.g. for skipping certain subcomponents or for calling the subcomponent’s methods in a different order, the composite component class can still override the original MNA-method with the mnaComp prefix instead of the mnaParent prefix. This will prevent the CompositePowerComp base class from doing any subcomponent handling in this specific MNA-method, so the subcomponent method calls have to be performed explicitly if desired. Given this, the following two implementations of the mnaAddPreStepDependencies method are equivalent:

void DP::Ph1::PiLine::mnaParentAddPreStepDependencies(AttributeBase::List &prevStepDependencies, AttributeBase::List &attributeDependencies, AttributeBase::List &modifiedAttributes) {
  // Only add the dependencies of the composite component, the subcomponent's dependencies are handled by the base class
  prevStepDependencies.push_back(mIntfCurrent);
  prevStepDependencies.push_back(mIntfVoltage);
  modifiedAttributes.push_back(mRightVector);
}

void DP::Ph1::PiLine::mnaCompAddPreStepDependencies(AttributeBase::List &prevStepDependencies, AttributeBase::List &attributeDependencies, AttributeBase::List &modifiedAttributes) {
  // Manually add pre-step dependencies of subcomponents
  for (auto subComp : mSubcomponentsMNA) {
    subComp->mnaAddPreStepDependencies(prevStepDependencies, attributeDependencies, modifiedAttributes);
  }

  // Add pre-step dependencies of component itself
  prevStepDependencies.push_back(mIntfCurrent);
  prevStepDependencies.push_back(mIntfVoltage);
  modifiedAttributes.push_back(mRightVector);
}

3.3.5 - Add New Model

This page walks through adding a component model, using a three phase dynamic phasor inductor as the example.

Where the code lives

Component models live in the dpsim-models subproject, which builds the CPS library. Headers and sources are separate trees, both organised by domain:

dpsim-models
 |- include
 |   \ dpsim-models
 |       |- Base            shared base classes, one per component family
 |       |- DP              dynamic phasor
 |       |- EMT             electromagnetic transient
 |       |- SP              static phasor
 |       \ Signal           domain independent control and signal models
 \- src
     |- Base
     |- DP
     |- EMT
     |- SP
     \ Signal

Namespaces follow the same shape, with the phase count nested inside the domain:

CPS::{DP,EMT,SP}::{Ph1,Ph3}::{Name}
CPS::Signal::{Name}

File names encode the same information, so the example model needs two files:

• dpsim-models/include/dpsim-models/DP/DP_Ph3_Inductor.h
• dpsim-models/src/DP/DP_Ph3_Inductor.cpp

declaring the class CPS::DP::Ph3::Inductor.

Choosing a base class

DPsim supports several solvers, and each requires certain member functions on the component. Which ones you implement is determined by the interfaces you inherit rather than by the solver itself.

For an MNA component, derive from MNASimPowerComp<VarType>, with Complex as the variable type in the DP and SP domains and Real in EMT. The MNA hooks are declared on MNAInterface, which MNASimPowerComp implements, so that is where to look for the full set.

If the model is naturally expressed as several existing components wired together rather than as a single stamp, derive from CompositePowerComp and add subcomponents instead. The pi-line is a worked example. See subcomponents.

If the component is better described by its own state-space model coupled to the network, see state-space nodal for that alternative.

Attributes

Every component exposes its parameters and state through attributes, declared in the class and registered in the constructor. Attributes are what make a value visible to the logger, to the Python bindings and to the task scheduler.

How to declare, read and derive attributes is described under attributes and attribute usage.

Tasks and step functions

Pre-step and post-step functions are registered as tasks, and the scheduler derives the order in which they may run from the attributes each task reads and writes. Declaring those dependencies correctly matters: a task that modifies an attribute without declaring it can be scheduled in the wrong order, or in parallel with a reader.

How tasks are built and how the dependency graph is derived is described under scheduling.

Registering the new component

A new component is not picked up automatically. Three places have to be updated:

• dpsim-models/src/CMakeLists.txt, adding the new source file to the list, for example DP/DP_Ph3_Inductor.cpp
• dpsim-models/include/dpsim-models/Components.h, adding the header so that including that one file gives access to every component
• dpsim/src/pybind/DPComponents.cpp, or the matching EMTComponents.cpp, SPComponents.cpp or SignalComponents.cpp, to expose the class to Python

The Python binding follows the existing pattern in those files:

py::class_<CPS::DP::Ph1::Resistor, std::shared_ptr<CPS::DP::Ph1::Resistor>,
           CPS::SimPowerComp<CPS::Complex>>(mDPPh1, "Resistor", py::multiple_inheritance())
    .def("set_parameters", &CPS::DP::Ph1::Resistor::setParameters, "R"_a);

Name the arguments using the "R"_a form shown above. Without it the Python signature and the generated reference fall back to positional placeholders such as arg0, and callers cannot use keyword arguments.

Initialization

Components are initialized either from power flow results or from explicitly set initial values, and the solver calls into the component in a defined order. Do not add an initialize(Real) overload of your own for user-facing initialization; use the documented hooks instead.

See initialization for the full sequence, and note the scaling conventions in guidelines, since initialization quantities are RMS3PH while EMT simulation quantities are PEAK1PH.

3.3.6 - Create New Simulation

Here, we will show the implementation of a new simulation scenario defined in C++, which is using DPsim as a library.

Directory Structure

In the end, your directory structure should look like as follows:

my-project
  |- CMakeLists.txt
  |- source
      |- my-scenario.cpp
  |- dpsim (as submodule)

CMake File

Your CMakeLists.txt could look like this:

cmake_minimum_required(VERSION 3.14)
project(my-project CXX)

add_subdirectory(dpsim)

add_executable(my-scenario source/my-scenario.cpp)
target_link_libraries(my-scenario dpsim)

Build the Project