Developer Guide

• 1: Architecture and Conventions
• 1.1: Architecture
• 1.2: Build
• 1.3: Coding Conventions
• 2: Attributes and Scheduling
• 2.1: Attributes
• 2.2: Attribute Usage Guidelines
• 2.3: Interfaces
• 2.4: Task Scheduling
• 2.5: Adding Tasks to a Component
• 3: Writing a Model
• 3.1: Component and Solver Initialization
• 3.2: Real-Time Execution
• 3.3: Interfacing with the MNA Solver
• 3.4: Subcomponent Handling
• 3.5: Add New Model
• 3.6: Create New Simulation
• 3.7: Debugging
• 3.8: Logger Implementation
• 4: Solvers
• 4.1: The MNA Solver
• 4.2: Power Flow Solvers
• 4.3: Alternative Solver Implementation
• 4.4: State-Space Extraction
• 5: Model Implementations
• 5.1: Reduced Order Generator Implementation
• 5.2: Switch and Load Implementation
• 5.3: Source Implementation
• 5.4: Injection and Compensation Implementation
• 5.5: SSN Component Implementation
• 5.6: Signal Component Implementation
• 5.7: DP Ph1 Averaged VSI Implementation
• 5.8: DP Ph3 Averaged VSI Implementation
• 5.9: EMT Ph3 Averaged VSI Implementation
• 5.10: EMT Ph3 Grid-Forming VSI Implementation

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.

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.

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

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.

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.

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.

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

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.

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.

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

The build process is similar to the one of DPsim:

cd my-project
mkdir build && cd build
cmake ..
make my-scenario

3.7 - Debugging

Mixed Python C++ Debugging

Prerequisites

Your vscode launch.json should have two configurations, one to launch the python process and one to attach gdb:

{
    "version": "0.2.0",
    "configurations": [
        {
            "name": "Python: Current File",
            "type": "python",
            "request": "launch",
            "program": "${file}",
            "console": "integratedTerminal",
            "stopOnEntry": true,
            "env": {"PYTHONPATH": "${workspaceFolder}/build${pathSeparator}${env:PYTHONPATH}"}
        },
        {
            "name": "(gdb) Attach",
            "type": "cppdbg",
            "request": "attach",
            "program": "/usr/bin/python",
            "processId": "${command:pickProcess}",
            "MIMode": "gdb",
            "setupCommands": [
                {
                    "description": "Enable pretty-printing for gdb",
                    "text": "-enable-pretty-printing",
                    "ignoreFailures": true
                }
            ]
        }
    ]
}

The python debugger will stop on entry (“stopOnEntry”: true). Make sure to adapt your PYTHONPATH variable if necessary.

The C++ code has to be build in debug mode

cmake .. -DCMAKE_BUILD_TYPE=Debug

Attaching C++ Debugger

• open the python example to be debugged
• go to the debug menu and select / run the “Python: Current File” configuration
• the python debugger should stop at entry
• set C++ breakpoints
• go to the debug menu and run the “(gdb) Attach” configuration
• select a process… choose the python process with the “—adapter-access-token” part
• you can view the whole description when you hover over the process with the mouse
• press play to continue Python debugging… the c++ debugger will stop at the next breakpoint

You can automate this by using the vscode extension “Python C++ Debugger” and by adding this configuration to the launch.json above:

{
    "name": "Python C++ Debugger",
    "type": "pythoncpp",
    "request": "launch",
    "pythonConfig": "custom",
    "pythonLaunchName": "Python: Current File",
    "cppConfig": "default (gdb) Attach"
}

This will automatically run both debuggers and select the current process.

It can take a while before the debugger hits the C++ breakpoints.

C++ Debugging

Use the following launch.json for vscode and set the program path:

{
    "version": "0.2.0",
    "configurations": [
        {
            "name": "(gdb) Launch",
            "type": "cppdbg",
            "request": "launch",
            "program": "${workspaceFolder}/dpsim/build/Examples/Cxx/example",
            "args": [],
            "stopAtEntry": true,
            "cwd": "${workspaceFolder}",
            "environment": [],
            "externalConsole": false,
            "MIMode": "gdb",
            "setupCommands": [
                {
                    "description": "Enable pretty-printing for gdb",
                    "text": "-enable-pretty-printing",
                    "ignoreFailures": true
                }
            ]
        }
    ]
}

3.8 - Logger Implementation

Using the loggers is covered under logging results. This page covers the classes.

Distinct things share the name

Only the data loggers have anything to do with results. CPS::Logger is a different subsystem that happens to share the word, and conflating the two is the most common confusion here. A component constructed with Logger::Level::debug writes prose about its own initialization and contributes nothing to any CSV.

DataLogger

Holds a map from column name to attribute and appends a row per step. log(Real time, Int timeStepCount) returns early when the logger is disabled or when timeStepCount % mDownsampling != 0, so down-sampling is a modulo on the step counter rather than a time comparison, and it is exact regardless of step size.

The header is written lazily on the first row, by testing mLogFile.tellp() == 0. That means the column set is fixed by whatever was registered before the first log call; registering an attribute afterwards would produce rows that no longer match the header.

Values are written with std::scientific in fixed-width columns, which is what makes the output readable as a table and also what makes it larger than a minimal CSV would be.

The constructor takes (name, enabled, downsampling). The Python binding exposes only the name, so enabled and downsampling are unreachable from Python. A binding that took all three would make down-sampling available to notebook users, who currently have only the time step.

RealTimeDataLogger

Exists because writing to disk inside a real-time step is not acceptable: the file system offers no bound on how long a write takes, and one slow write overruns the step. It preallocates mAttributeData from either the final time and step size or an explicit row count, fills it during the run, and writes at the end.

The preallocation is the point, and it is also the constraint: the row count must be known before the run, so a real-time simulation of indefinite length needs a different arrangement.

DataLoggerInterface

The abstract seam. Implement it to send results somewhere other than a file, which is what the co-simulation interfaces do rather than logging and re-reading. Simulation::addLogger accepts anything implementing it.

Scheduling

A logger contributes a task like any other component, so the scheduler places it by its declared attribute dependencies. A logged attribute therefore keeps alive the task that produces it, which has a consequence worth knowing: logging an attribute can change which tasks the scheduler considers reachable. A model whose results change when a logger is added is exhibiting a missing dependency declaration elsewhere, not a logging bug. See adding tasks to a component.

Source

dpsim/src/DataLogger.cpp, dpsim/src/RealTimeDataLogger.cpp, dpsim/include/dpsim/DataLoggerInterface.h, and dpsim-models/include/dpsim-models/Logger.h for the unrelated diagnostic logger.

4 - Solvers

The MNA solver is the default and the one almost every simulation uses; a component’s side of it is under interfacing with the MNA solver. The remaining pages cover the powerflow solver, the alternatives to nodal analysis, the linear algebra backends, and the state-space model that can be recovered from a running simulation.

The methods themselves are derived under alternative solution methods.

4.1 - The MNA Solver

MnaSolver<VarType> is the solver almost every simulation uses. The method it implements is derived under nodal analysis; what a component must provide to take part is under interfacing with the MNA solver. This page is the solver itself.

Note the spelling: the class is MnaSolver even though the file is MNASolver.h.

Setting up

initialize runs the sequence described under component and solver initialization: identify the topology objects, create sub-components, collect virtual nodes, assignMatrixNodeIndices, size the matrices, initialize the components, then assemble.

Which assembly function runs depends on the network:

• initializeSystemWithPrecomputedMatrices when the switch combinations are few enough to enumerate. Every combination gets its own factorised matrix up front, so a switching event becomes a lookup rather than a refactorisation.
• initializeSystemWithVariableMatrix when a component changes its own stamp continuously and enumeration is impossible.
• initializeSystemWithParallelFrequencies for a harmonic study, where several frequencies are solved side by side.

resolveSystemMatrixRecomputationMode chooses between them when the mode is Auto; SystemMatrixRecomputationMode::Enabled and Disabled force it either way.

Stepping

solve does the same four things every step.

It zeroes the right-hand side and sums the stamps the components’ pre-step tasks produced, which is why a component that fails to declare its dependencies can find its contribution missing rather than wrong. It calls updateSwitchStatus, which produces an index into the precomputed matrices. It solves through the linear solver for that index. Then it hands the solution to the components' post-step tasks.

The switch index is the point of the precomputed strategy: with the factorisations already built, a switching event costs a different lookup rather than new numerical work. That is what makes a network with frequent switching affordable, and it is why the number of switches is bounded in practice, since the enumeration grows as two to the power of that number.

solveWithSystemMatrixRecomputation is the other path. It asks hasVariableComponentChanged each step and rebuilds and refactorises only when something reports a change, which is the expensive but general case used by variable components such as the SSN models.

Iterative components

After the solve, the solver checks whether any synchronous generator reports requiresIteration. If so it repeats the solve step until none does, which is how the predictor-corrector and two-stage machine models reach the implicit solution rather than its explicit approximation. Models that do not request iteration cost nothing here.

This loop is the reason a machine model can be iterative without the whole solver being iterative.

Linear backends

The solver does not implement its own factorisation; it selects an adapter through MnaSolverFactory. The choices and their tuning are described under alternative solvers, which also covers the ordering and partial-refactorisation options that matter most when the matrix changes every step.

Instrumentation

Solver::mLogSolveTimes records the wall-clock duration of each solve into mSolveTimes, which is the measurement to use when comparing backends or step sizes rather than timing the whole run.

Source

dpsim/src/MNASolver.cpp, dpsim/src/MNASolverDirect.cpp, and dpsim/include/dpsim/MNASolverFactory.h.

4.2 - Power Flow Solvers

What DPsim implements. For the underlying formulation, the mismatch function and the Jacobian, see power flow.

Solver Implementations

DPsim ships two implementations of the Newton-Raphson power flow solver with power mismatch and polar coordinates. Both produce identical results (to round-off); they differ only in how the Jacobian is stored and factorized:

• PFSolverPowerPolar (dense): assembles a dense Jacobian and computes a fresh factorization every Newton iteration. This is the default.
• PFSolverPowerPolarSparse (sparse): assembles the Jacobian into a sparse matrix whose sparsity pattern is fixed (derived once from the network admittance matrix). The symbolic factorization (ordering) is analyzed once and reused; only the numeric values are recomputed each Newton iteration. The first iteration of every power flow solve does a full factorization with pivoting, and subsequent iterations refactorize while reusing that ordering (via KLU when available). This scales better on large, sparse grids.

The dense solver is used by default. To opt in to the sparse solver:

sim.set_pf_solver_use_sparse(True)

The flag is ignored and the dense solver is used if DPsim was built without a sparse linear solver. The benchmark notebook examples/Notebooks/Grids/PF_Sparse_vs_Dense.ipynb runs a range of network sizes both ways, verifies the converged voltages match, and compares run time.

Generator Reactive Power Limits

A PV bus assumes its generator can produce whatever reactive power the Newton-Raphson solution asks for, holding $\vert V_k \vert$ at its setpoint. Real generators cannot: Q is bounded by $Q_{min}$ and $Q_{max}$. DPsim can enforce these bounds with a bidirectional PV↔PQ outer loop:

1. Run the inner Newton-Raphson solve to convergence (as described above).
2. For every PV bus, compute the generator’s actual reactive output. If it exceeds $Q_{max}$ or falls below $Q_{min}$, pin the injection at the violated limit and convert the bus to PQ.
3. For every bus pinned this way in an earlier pass, check whether the constraint has relaxed: if $\vert V_k \vert$ has moved past its original setpoint in the releasing direction, restore voltage control and convert the bus back to PV.
4. Repeat from step 1 until no bus switches, an outer-iteration cap is hit, or a per-bus switch counter trips (an anti-oscillation guard, since a bus can otherwise toggle PV↔PQ indefinitely near the boundary).

Enforcement is opt-in and defaults off, so a system with no limits configured behaves exactly as before:

sim.set_pf_solver_enforce_q_limits(True)

$Q_{min}$/$Q_{max}$ are set per generator via SynchronGenerator.set_parameters(..., q_limit_max=..., q_limit_min=...); the defaults are $\pm\infty$ (unlimited). Generators sharing a bus have their limits summed. The two limits are enforced independently, with no assumption about sign or relative magnitude: asymmetric bounds (e.g. $Q_{max}=150$ MVAr, $Q_{min}=-30$ MVAr) and same-sign bounds (e.g. a generator restricted to $Q \in [20, 150]$ MVAr, always producing, or $Q \in [-150, -20]$ MVAr, always absorbing) are both enforced correctly.

Limitation: no P-dependent capability curve. $Q_{min}$ and $Q_{max}$ are constants set once per generator, not a function of active power output $P$. A real synchronous generator’s reactive capability is a “D-curve” bounded by three physically distinct mechanisms: the stator (armature) current limit $\sqrt{P^2+Q^2} \le S_{rated}$, the rotor (field) current / heating limit on the over-excited (Q-providing) side, and the under-excitation limiter (UEL) / steady-state stability limit on the under-excited (Q-absorbing) side. All three tighten as $P$ approaches rated output, and the over- and under-excited bounds come from unrelated physical limits, so the true feasible region is neither symmetric in $Q$ nor independent of $P$. DPsim does not model this curve; a generator’s $Q$ headroom is the same regardless of how much $P$ it is producing at the time. Flat per-generator limits are a common baseline in power-flow tools generally, so this is not a regression, but a P-dependent capability curve is not currently implemented.

The notebook examples/Notebooks/Grids/PF_Generator_Qlimits.ipynb validates the switching behavior on a small hand-wired case (binding and non-binding limits, dense vs. sparse agreement).

4.3 - Alternative Solver Implementation

The methods are derived under alternative solution methods. This page covers the code and the configuration.

DAESolver

Wraps the IDA integrator from Sundials. initialize builds the state and derivative vectors, registers each component’s residual function, then creates the solver with IDACreate, passes the solver instance as user data so the residual callbacks can reach it, and sets scalar relative and absolute tolerances with IDASStolerances.

Components take part by implementing DAEInterface and contributing their residual. The offset vector recorded at the top of the file defines how each component’s block is laid out within the global residual.

Two practical notes. The solver chooses its own steps, so a run’s cost is not predictable and it cannot be used under a real-time timer. And several std::cout calls remain in the initialization path, so it prints to standard output independently of the logger.

ODESolver and ODEintSolver

ODESolver wraps CVODE from Sundials for a single component, sizing the problem from mOdePreState and attaching a dense linear solver. ODEintSolver does the same job with boost’s odeint, calling comp->odeint(y, ydot, t).

Both integrate one component across a network step while the network itself stays on its fixed step, so the coupling is staggered and first-order accurate regardless of the inner integrator’s order.

DiakopticsSolver

Constructed with the system and an explicit list of components to tear, which must implement MNATearInterface. system.splitSubnets performs the partition, initSubnets builds the per-subnetwork node and component lists, and mNodeSubnetMap records which subnetwork owns each node.

createTearMatrices is specialised per value type, and the sizes differ in a way worth noting: the Real specialisation allocates tearComponents * phaseMultiplier, while the Complex one allocates twice that, because a complex quantity is carried as a real-augmented pair. The phase multiplier is 3 when the subnetwork phase type is ABC and 1 otherwise, taken from the first node of the system.

The removed-branch system is dense and small. A comment in the source notes that the reduction could still be sped up by exploiting the block diagonal structure of the inverse, so the present implementation is correct rather than optimal.

Linear backends under MNA

DirectLinearSolverConfiguration tunes the chosen backend, and not every option applies to every one:

• SCALING_METHOD: none, sum or max
• FILL_IN_REDUCTION_METHOD: AMD, AMD_NV, AMD_RA or COLAMD. The NV and RA variants take the set of time-varying entries into account when ordering, which is what makes partial refactorization effective for a network with switching elements.
• PARTIAL_REFACTORIZATION_METHOD: none, factorization path, or refactorization restart. This is the lever that matters when a switch or a variable component changes the matrix every step.
• USE_BTF: block triangular form on or off

The defaults are chosen for a general network. The combination of an ordering that knows about varying entries with partial refactorization is what makes repeated switching affordable, and it is inert if the matrix never changes.

Source

Under dpsim/src/: DAESolver.cpp, ODESolver.cpp, ODEintSolver.cpp, DiakopticsSolver.cpp, and the six *Adapter.cpp files.

4.4 - State-Space Extraction

The method itself, what the extracted model means and where it is valid, is derived under state-space extraction. This page covers enabling it and reading the result.

State-space extraction is optional and can be enabled through the Simulation API. During simulation setup, the MNA solver creates an MNAStateSpaceExtractor. During the solver task flow, a state-space extraction task uses the active direct linear solver to update the extracted discrete-time state matrix.

Main classes

The implementation is organized around three main parts:

• MNAStateSpaceExtractor assembles and stores the extracted discrete-time state matrix.
• MNAStateSpaceContributor represents the state-space contribution of one supported component.
• MNAStateSpaceContributorFactory creates contributors for supported MNA components.

The extractor is owned by the MNA solver. Component contributors are created during solver initialization and are used to stamp the local matrices needed for the MNA-coupled state-space formulation.

For the components that support extraction in each domain, see state-space extraction support.

Usage

In C++, state-space extraction can be enabled as follows. The example below uses EMT Ph3; replace Domain::EMT with Domain::DP for DP Ph1:

Simulation sim("Example");
sim.setDomain(Domain::EMT);
sim.setSolverType(Solver::Type::MNA);
sim.doStateSpaceExtraction(true);
sim.run();

const auto &extractor = sim.getStateSpaceExtractor();
const Matrix &Ad = extractor.getDiscreteStateMatrix();

In Python, the corresponding API is shown below. Replace dpsimpy.Domain.EMT with dpsimpy.Domain.DP for DP Ph1:

sim = dpsimpy.Simulation("Example")
sim.set_domain(dpsimpy.Domain.EMT)
sim.set_solver(dpsimpy.Solver.MNA)
sim.do_state_space_extraction(True)
sim.run()

extractor = sim.get_state_space_extractor()
Ad = extractor.get_discrete_state_matrix()

Modal analysis of the extracted model

StateSpaceModalAnalysis is constructed from an MNAStateSpaceExtractor and computes the modes of whatever the extractor last produced. The method is described under modal analysis.

update() runs Eigen::EigenSolver on the discrete state matrix and throws if it does not converge. It then maps each discrete eigenvalue to the continuous plane with 2 / dt * (z - 1) / (z + 1) and keeps both sets, retrievable through getDiscreteEigenvalues and getContinuousEigenvalues.

Participation factors are the elementwise product of the right eigenvectors with the transpose of the left ones. They require inverting the right eigenvector matrix, so update() throws with an explicit message when that matrix is singular. That happens for a defective state matrix, which is a property of the system rather than a numerical problem; the eigenvalues are still valid in that case, only the participation factors are unavailable.

setAnalysisFrame selects between StateSpaceAnalysisFrame::Native, which analyses the states as the components hold them, and GlobalDQ0, which transforms into one common frame first. The second needs setGlobalDq0Frame(omega, theta0). getStateNames returns names matching the frame in use, so a participation factor can be attributed to a named state rather than to an index.

Examples

The feature is demonstrated in:

• EMT Ph3 RLC extraction
• EMT Ph3 composite extraction
• DP Ph1 RLC extraction
• DP Ph1 composite extraction

Equivalent Python notebooks are available in examples/Notebooks/StateSpace.

5 - Model Implementations

One page per model family, covering the class hierarchy, how the component interfaces with the solver, its attributes and state layout, and the traps in configuring it. The equations behind each are under Concepts, which names no class; these pages name nothing else.

Which domains implement which model is in model availability, generated from the headers.

5.1 - Reduced Order Generator Implementation

The equations are derived under reduced order machine models. This page covers only their arrangement in code.

Class hierarchy

Base::ReducedOrderSynchronGenerator<VarType> holds everything independent of domain and of order: the per unit base values, the operational parameters, the mechanical states, the controller attachments and the discretisation coefficients. It is templated on Real for EMT and Complex for DP and SP, which is why the axis frame quantities appear twice, as mVdq0/mIdq0 in the real specialisation and mVdq/mIdq in the complex one.

Each domain then provides a ReducedOrderSynchronGeneratorVBR layer holding the frame transform, and each order a concrete class. The order is recorded in mSGOrder, which selects which coefficients are computed.

Network interface

setModelAsNortonSource chooses between the two interface forms. The default is the Norton equivalent, in which the machine contributes only to the right hand side vector and requests no virtual nodes. The Thevenin form requests two virtual nodes instead. Both represent the same model; the Norton form is cheaper because it leaves the system matrix untouched between steps [Wang2010].

Coefficients

calculateAuxiliarConstants computes the discretisation coefficients once, since they depend only on the parameters and the step size. The member names map to the symbols on the theory page as follows.

The naming looks wrong at first and is not. Zd_t is built from $L_q - L_q’$ and Zq_t from $L_d - L_d’$, because each is named for the axis whose coefficient it feeds rather than for the parameters it is assembled from. That follows the physics: the d-axis internal voltage arises from q-axis rotor flux and decays with $T_{q0}’$, so mAd_t correctly combines $L_q - L_q’$ with $T_{q0}’$ and multiplies the q-axis current.

Read a coefficient’s use rather than its assignment line before concluding an axis is swapped.

Step sequence

mnaCompPreStep runs before the network solve and does three things in order. It advances the controllers, saving mEf_prev and mMechTorque_prev first because the trapezoidal history terms need the previous values. It calls stepInPerUnit, which updates the frame transforms from mThetaMech, recomputes the axis frame state from the terminal quantities, and evaluates the history voltage into mEh_vbr. It then stamps the result into the right hand side vector.

Each concrete order implements only specificInitialization and stepInPerUnit. Everything else is inherited.

Initialization

Initialization runs from the powerflow solution, not from user supplied states. The base class computes the load angle as the phase of $V + j L_q I$, projects the terminal voltage and current onto the axis frame, and derives the field voltage from the no-load relation. Only then does specificInitialization set the order specific states, which is why a concrete class can assume mVdq and mIdq are already populated.

Attached controllers are initialized afterwards from the machine’s own initial values, so an exciter or governor never needs its own operating point.

Controllers

Excitation, governor, turbine and power system stabilizer attach through the base class and are optional, guarded by mHasExciter, mHasGovernorAndTurbine, mHasTurbineGovernor and mHasPSS. The stabilizer output feeds the exciter within the same step, and the governor output feeds the turbine, so the order of the calls in mnaCompPreStep is load bearing.

Source code

• Base class header, implementation
• Concrete orders under dpsim-models/src/{SP,DP,EMT}/ named SynchronGenerator<N>OrderVBR
• Availability per domain is listed under model availability

References

• [Wang2010] IEEE Xplore document 5411963. Cited in the machine model pages as the basis for interfacing a machine to a nodal solver through a current source that leaves the system matrix unchanged.

5.2 - Switch and Load Implementation

The models are derived under switches and loads. This page covers only their arrangement in code. Availability per domain is in model availability.

Switches

Switch implements Base::Ph1::Switch and stamps one admittance chosen by mIsClosed, using MNAStampUtils::stampAdmittance so the grounded-terminal cases are handled centrally. SeriesSwitch folds a series resistance into the same branch.

varResSwitch additionally implements MNAVariableCompInterface, which is what allows it to change the system matrix during a run. Its hasParameterChanged is called each step and drives the transition:

• Opening multiplies the resistance by mDeltaResOpen each step until it passes the target open value, then clamps to it and reports the transition finished.
• Closing uses mDeltaResClosed, which is 0, so the first step takes the resistance to zero, the clamp catches it and sets the closed value. Closing is therefore immediate by construction, not by a separate code path.

Its initializeFromNodesAndTerminals carries a comment saying it is not used.

Loads

RXLoad is a CompositePowerComp. In initializeFromNodesAndTerminals it converts the powers to element values and builds sub-components:

• a resistor, only if the active power is non-zero
• an inductor if the reactance is positive, a capacitor if negative, and nothing if the reactive power is zero

PQLoadCS wraps a current source and sets its reference in updateSetPoint from conj(S / mNomVoltage). The nominal voltage, not the terminal voltage, is deliberate; the line using the terminal voltage is present but commented out. Changing it would make the component nonlinear and require an iterative solve.

Shunt takes a conductance and a susceptance directly and additionally carries per-unit attributes, since it is the form the powerflow solver consumes.

Source

• Switches: {SP,DP,EMT}_Ph{1,3}_Switch, DP_Ph3_SeriesSwitch, EMT_Ph3_SeriesSwitch, {DP,SP}_Ph1_varResSwitch under dpsim-models/src/
• Loads: DP_Ph1_RXLoad, EMT_Ph3_RXLoad, DP_Ph1_RXLoadSwitch, DP_Ph1_PQLoadCS, SP_Ph1_Load, {SP,DP}_Ph1_Shunt, EMT_Ph3_Shunt

5.3 - Source Implementation

The models are derived under sources. This page covers only the code.

Ideal sources

VoltageSource requests one virtual node, which carries the source current as the extra unknown, and stamps the constraint rows that fix the terminal voltage difference. CurrentSource requests none and contributes only to the right hand side.

ControlledVoltageSource and ControlledCurrentSource are the same components with their reference supplied as an attribute rather than a parameter, so another component or an interface can drive them. The reference is read during the pre-step, which is why it is the previous step’s value.

VoltageSourceNorton

Stamps directly rather than through a virtual node. mnaCompApplySystemMatrixStamp adds mConductance to both diagonal entries and subtracts it from the two off-diagonal entries, guarded by terminalNotGrounded, and mnaCompApplyRightSideVectorStamp sets the equivalent current mIntfVoltage / mResistance with opposite signs at the two terminals.

mConductance is computed in setParameters as 1 / resistance, so calling setParameters is mandatory before the run and a zero resistance is a division by zero rather than an ideal source.

VoltageSourceRamp

A composite wrapping a VoltageSource whose reference it rewrites each step in updateState(time). Three regimes: before mSwitchTime the reference is unchanged; during mRampTime the added voltage is interpolated linearly while the added frequency is blended by a raised sine 0.5 + 0.5 * sin(pi * t / T - pi/2); afterwards both are fully applied.

The two are blended differently on purpose. A linear frequency interpolation applied as a phase offset would step the phase at both ends of the ramp; the raised sine has zero derivative at both ends, so the frequency contribution enters and leaves smoothly. The consequence is that the instantaneous frequency during the ramp is not the linear interpolation between the two values, and reading mAddSrcFreq as “the frequency at the midpoint” is wrong.

Note also that the added frequency term is applied as mAddSrcFreq * time, using absolute simulation time rather than time since the switch, so the phase contribution depends on when in the run the ramp occurs.

ProfileVoltageSource

Holds a std::filesystem::path, a sample vector and an index, and reads the file in readFromFile at construction. It implements DAEInterface in addition to the MNA hooks.

The samples are stepped by index rather than interpolated against simulation time, so the profile’s sample rate and the simulation step must match for the waveform to have the intended duration. It is bound in Python and constructing it with a file that is not a readable sample list raises rather than crashing, which is covered by a test.

Source

Under dpsim-models/src/{SP,DP,EMT}/. Availability per domain is in model availability.

5.4 - Injection and Compensation Implementation

The models are derived under network injection and compensation. This page covers only the code.

NetworkInjection

A CompositePowerComp wrapping a single VoltageSource sub-component. It owns no equations of its own; it exists so that the external network is a named component rather than a bare source, and so that the driving waveform can be swapped without changing the network description.

setParameters is overloaded by the kind of generator wanted behind it: a constant phasor for a fixed source, a start frequency with a rate of change for a ramp, and an initial phasor with a modulation frequency for a modulated one. Which overload is called determines which SignalGenerator the sub-source is given; see signal component implementation.

Because the source is ideal, adding an impedance to represent a finite short circuit level is the caller’s job. Nothing in the component does it.

SVC

Not composite. It computes a susceptance each step and realises it by reconfiguring an internal reactive element, so it implements the variable-component interface and forces a refactorisation whenever the value changes.

updateSusceptance performs both lags with the trapezoidal rule, using precomputed constants Fac1 = dt / (2 Tr), Fac2 = dt Kr / (2 Tr) and Fac3 = dt / (2 Tm). The measurement lag is applied first, then the error is formed in per unit against mNomVolt, then the susceptance follows from the previous value and the present and previous error.

The result is clamped to mBMax and mBMin before use, and the internal element is only rebuilt when the value actually changed. The sign of the clamped susceptance selects which element is formed: positive gives an inductance 1 / (omega * B * mBN), negative a capacitance B * mBN / (-omega). mBN is the base susceptance, so B is per unit.

SolidStateTransformer

A CompositePowerComp that represents each side as a current source rather than as a coupled winding pair. setParameters(nomV1, nomV2, Pref, Q1ref, Q2ref) takes the two nominal voltages and three power set points; the active power is common to both sides, while the reactive powers are set per side.

Values are held in per unit internally, so the nominal voltages are the base rather than a turns ratio. There is no magnetising branch, no leakage impedance and no angle dependence, which is the representation the concept page describes and not an omission.

Source

• DP_Ph1_NetworkInjection, and the SP and EMT variants alongside it
• DP_Ph1_SVC
• SP_Ph1_SolidStateTransformer

Availability per domain is in model availability.

5.5 - SSN Component Implementation

The method is derived under state-space nodal components and state-space nodal. This page covers only the code.

Base hierarchy

SSNComp holds the continuous matrices, the discrete pair, the equivalent admittance mW, the history vector mYHist and the state attribute x. Two branches specialise it by which quantity is the input:

• VTypeSSNComp takes voltage in and gives current out, so it stamps as an admittance
• ITypeSSNComp is the dual

Terminal-count layers sit on top: TwoTerminalVTypeSSNComp, TwoTerminalITypeSSNComp and FourTerminalVTypeSSNComp handle the mapping from terminal quantities to the model input, and the Variable layers add re-forming of the model between steps. Every EMT::Ph3 base sets PhaseType::ABC in its constructor, which the concrete components rely on.

What a component provides

A fixed-model component only calls SSNComp::setParameters(A, B, C, D) with its chosen state, input and output. EMT::Ph3::SSN::Inductor is the whole pattern:

Matrix aMatrix = Matrix::Zero(3, 3);   // x = i_abc
Matrix bMatrix = inductance.inverse(); // u = v_abc
Matrix cMatrix = Matrix::Identity(3, 3);
Matrix dMatrix = Matrix::Zero(3, 3);
SSNComp::setParameters(aMatrix, bMatrix, cMatrix, dMatrix);

The base does the rest: recomputeDiscreteModel calls Math::calculateStateSpaceTrapezoidalMatrices and sets mW = mC * mdB + mD, calculateHistoryVector returns mC * (mdA * x + mdB * u), and the post step updates the state from the old and new input.

A varying component additionally overrides updateStateSpaceModel (a no-op for linear components) and, for the variable bases, updateComponentParameters to report whether the model changed. Only when it reports a change is the system matrix refactorised.

Domain differences

The formulation differs by domain, and so does the code path. The theory is under SSN across domains.

A component supplies the same real (A, B, C, D) in either domain. EMT::SSNComp discretises them directly. DP::SSNComp does not: buildAugmentedA(omega) assembles the real-augmented 2n x 2n matrix with A on both diagonal blocks and +wI / -wI off-diagonal, buildAugmentedB places B on both diagonal blocks, and the result goes through the same Math::calculateStateSpaceTrapezoidalMatrices helper as EMT. The discrete blocks are then folded back into complex form as topLeft + j * bottomLeft, which is the inverse of the [[P, -Q], [Q, P]] representation. mW and the history vector are complex as a result.

recomputeDiscreteModel therefore takes omega in DP and takes no argument in EMT. A component that hardcodes a frequency here rather than using the value handed to mnaCompInitialize is wrong at any other system frequency.

Frame metadata

getLocalAbcStateBlocks returns nothing by default and should be overridden only for states that genuinely form physical abc triples. It is consumed by tooling that reasons about the state vector in the phase frame, and declaring a block that is not one produces wrong groupings rather than an error.

Initialization

calculateSteadyStateStateFromInput evaluates (jωI − A)⁻¹ B u, which requires the continuous model to be set first. Components with real control states cannot use the default initializeFromNodesAndTerminals on the mixed base; see DP Ph1 averaged VSI implementation for that case and for the requirement that the state matrix be handed over already carrier shifted.

The components

Fixed models: SSN_Full_Serial_RLC, SSN_Capacitor, SSN_Inductor, SSNTypeV2T, SSNTypeI2T. Varying models: SSN_Variable_Serial_RLC, PiecewiseLinearInductor, and the inverter models under power electronics. The Generic two- and four-terminal classes take the matrices from the caller instead of forming them, so they are the route to an SSN component without writing C++. Availability per domain is in model availability.

5.6 - Signal Component Implementation

The models are derived under converter control and signal processing blocks. This page covers only the code.

Base and scheduling

Signal blocks derive from SimSignalComp and take no part in the nodal solve. They contribute tasks through getTasks() rather than through the MNA hooks, and the scheduler orders them from the attribute dependencies those tasks declare. A block that reads an attribute without declaring it may still produce the right answer, by luck of ordering, and then change behaviour when an unrelated component is added; see adding tasks to a component.

Most blocks follow a two-task shape: a PreStep that copies the current values into the previous ones, and a Step that computes the new state and output. The split exists so that a value consumed by another block within the same step is unambiguous about which timestep it belongs to.

The mInputPrev / mInputCurr pattern

Blocks that integrate with the trapezoidal rule need both the present and the previous input, so they carry mInputPrev, mInputCurr, mStatePrev, mStateCurr and the matching output pair. The PreStep task performs the shift. Integrator::signalStep is the whole pattern:

**mStateCurr = **mStatePrev + mTimeStep / 2.0 * **mInputCurr
                            + mTimeStep / 2.0 * **mInputPrev;
**mOutputCurr = **mStateCurr;

VCO::signalStep deliberately does not, using state + dt * input, because it accumulates an angle rather than integrating a control signal.

Every one of these blocks needs setSimulationParameters(timestep) before the run, since the step size appears directly in the update. Blocks that expose setInitialValues must also have it called, or they start from zero rather than from the operating point.

State-space blocks

PLL is written as an explicit state-space block rather than as arithmetic, setting

mA << 0, mKi, 0, 0;
mB << 1, mKp, 0, 1;
mC << 1, 0, 0, 1;
mD << 0, 0, 0, 0;

The first input is the nominal frequency and is held constant, which is how the feed-forward term enters. Writing it this way means the block can be discretised by the same helpers as anything else rather than by hand.

FIRFilter

FIRFilter keeps a circular buffer and a write index, and step sums mFilter[i] * mSignal[...] over the filter length before advancing the index. It contributes a single Step task. Filter coefficients are supplied by the caller; nothing validates their length against the buffer or checks that they sum to a sensible gain.

Generators

SignalGenerator is the abstract base; SineWaveGenerator, DCGenerator, CosineFMGenerator and FrequencyRampGenerator are the concrete ones, and all expose their value through a sigOut attribute that a source component references.

Source

Under dpsim-models/src/Signal/. Availability is in model availability; these blocks are domain independent and appear there as a list rather than a matrix.

5.7 - DP Ph1 Averaged VSI Implementation

The equations are derived under DP Ph1 averaged voltage source inverter. This page covers only their arrangement in code.

Class and base

DP::Ph1::AvVoltSourceInverterStateSpace is final and derives from DP::Ph1::MixedVTypeVariableSSNComp. The mixed base is what makes the model possible in this domain: eight of the twelve states are real baseband control states and only the last four are the real and imaginary parts of the two carrier-band envelopes, so the component cannot use the plain complex SSN base.

State layout

The state order is fixed by a private StateIndex enum, which the linearization indexes directly.

The base does not impose this ordering. It is told only how many real and how many complex states there are, and sizes the packed real vector as realStateCount + 2 * complexStateCount. The three-phase model orders its states the other way round, envelopes first and controls after, and is equally valid. What the base does require is that the derived class hand it a state matrix that is already carrier shifted: the steady-state solve assumes it, and a model that supplies an unshifted matrix initializes to the wrong operating point rather than failing.

The default initializeFromNodesAndTerminals throws unless realStateCount is zero, so any model with real control states, which includes this one, must override it.

Tracking Psi rather than the raw PLL angle keeps the tracked quantity bounded. The raw angle grows without limit, which costs relinearization accuracy as a run gets longer.

Parameters

initializeFromNodesAndTerminals derives the initial state from the connected node voltage, so the operating point comes from the powerflow rather than from user supplied states.

Source and examples

• Source: header, implementation
• C++ example
• Python notebook
• The state-space nodal interfacing this model relies on is described under state-space nodal

5.8 - DP Ph3 Averaged VSI Implementation

The equations are derived under DP Ph3 averaged voltage source inverter. This page covers only their arrangement in code.

Class and base

DP::Ph3::AvVoltSourceInverterStateSpace is final and derives from DP::Ph1::MixedVTypeVariableSSNComp, the same mixed base as the single-phase model. Per-phase complex quantities are carried as std::array<Complex, 3>.

State layout

Twenty states by default, ordered envelopes first and controls afterwards, or twenty-two with the optional negative-sequence loop enabled.

This is the reverse of the single-phase ordering, which places controls first. The base does not care: it is given only the counts of real and complex states and sizes the packed real vector as realStateCount + 2 * complexStateCount. What it does require is a state matrix that is already carrier shifted, since the steady-state solve assumes it.

The last two states are the difference from the single-phase model beyond the per-phase filter. Three independent phase envelopes admit a negative-sequence component that a single positive-sequence envelope cannot represent, so the controller carries its own negative-sequence integrator pair.

Enabling the negative-sequence loop

The constructor takes an enableNegSeqControl flag, false by default. The two references $i_{nd,\mathrm{ref}}$ and $i_{nq,\mathrm{ref}}$ are the last two arguments of setParameters and default to zero, which makes the loop a suppressor rather than an injector. The measured $i_{rc,nd}$ and $i_{rc,nq}$ are exposed as the irc_n_d and irc_n_q attributes, and stay at zero while the loop is disabled.

The two integrators are appended after the control block rather than inserted next to the other control states, so enabling the flag leaves every envelope and positive-sequence control index unchanged. Code indexing into the state vector therefore does not need to know about the flag.

The theory behind the loop is derived under DP Ph3 averaged VSI.

Source and examples

• Source: header, implementation
• C++ example
• Python notebook
• Single-phase counterpart: DP Ph1 averaged VSI implementation

5.9 - EMT Ph3 Averaged VSI Implementation

The equations are derived under EMT Ph3 averaged voltage source inverter. This page covers only their arrangement in code.

Class and base

EMT::Ph3::AvVoltSourceInverterStateSpace is final and derives from EMT::Ph3::TwoTerminalVTypeVariableSSNComp. Unlike the dynamic phasor ports of this model, every state here is real, so it uses the plain variable state-space nodal base rather than the mixed one. The base sets PhaseType::ABC in its constructor, which this component relies on.

State layout

Fourteen real states, controls first and filter states afterwards.

The first state is the raw PLL angle. The dynamic phasor ports track the deviation from the nominal carrier phase instead, because there the angle is compared against a carrier and an unboundedly growing value costs relinearization accuracy. In EMT there is no carrier to drift against, so the raw angle is used directly.

Six real filter states here correspond to two complex envelopes in the single-phase dynamic phasor model and six in the three-phase one. That correspondence is the practical statement of what the envelope transform buys.

Source and examples

• Source: header, implementation
• C++ example
• Python notebook
• Dynamic phasor counterparts: DP Ph1, DP Ph3

5.10 - EMT Ph3 Grid-Forming VSI Implementation

The equations and the linearization mathematics are derived under EMT Ph3 grid-forming voltage source inverter. This page covers only their arrangement in code.

Class and base

EMT::Ph3::SSN_GFM is final and derives from EMT::Ph3::TwoTerminalVTypeVariableSSNComp. All seventeen states are real.

State layout

Omega and Theta being states rather than inputs is what makes this grid forming: the converter carries its own frequency and angle instead of tracking a measured one through a PLL.

Numerical linearization

The Jacobians are not written out by hand. calculateNumericalJacobians forms all four by central differences of the nonlinear state and output functions, so a change to the control equations needs no matching change to any matrix code.

The perturbation for column $j$ is absoluteStep + relativeStep * max(1, |x_j|), defaulting to 1e-8 and 1e-6 and adjustable at runtime. The max(1, ...) floor means the step is effectively absolute for small states and relative for large ones, which keeps the difference well conditioned across states whose magnitudes differ by orders.

Because the model is time varying, the state-space form and its stamp are recomputed every step rather than cached. That is the cost of this approach and the reason it is used only where the control is genuinely nonlinear.

Source and examples

• Source: header, implementation
• C++ example, IEEE 9-bus mixed machine and inverter
• Python notebook
• Grid-following counterpart: EMT Ph3 averaged VSI implementation