Documentation

• 1: Overview
• 1.1: Architecture
• 1.2: Attributes
• 1.3: Scheduling
• 1.4: Interfacing with the MNA Solver
• 1.5: Subcomponent Handling
• 1.6: Interfaces
• 2: Getting Started
• 2.1: Build
• 2.2: Install
• 2.3: Real-Time
• 3: Concepts
• 3.1: Dynamic Phasors
• 3.2: Nodal Analysis
• 3.3: Powerflow
• 4: Development
• 4.1: Debugging
• 4.2: Guidelines
• 5: Models
• 5.1: Transformer
• 5.2: Branches
• 5.3: Induction Machine
• 5.4: RLC-Elements
• 5.5: Synchronous Generator
• 5.6: VS-Inverter
• 6: Core Tasks
• 6.1: Add New Model
• 6.2: Create New Simulation
• 7: Examples
• 8: Roadmap
• 9: Reference
• 10: Contribution Guidelines

1 - Overview

DPsim is a real-time capable power system simulator that supports dynamic phasor and electromagnetic transient simulation as well as continuous powerflow. It primarily targets large-scale scenarios on commercial off-the-sheld hardware that require deterministic time steps in the range of micro- to milliseconds.

DPsim supports the CIM format as native input for the description of electrical network topologies, component parameters and load flow data, which is used for initialization. For this purpose, CIM++ is integrated in DPsim. Users interact with the C++ simulation kernel via Python bindings, which can be used to script the execution, schedule events, change parameters and retrieve results. Supported by the availability of existing Python frameworks like Numpy, Pandas and Matplotlib, Python scripts have been proven as an easy and flexible way to codify the complete workflow of a simulation from modelling to analysis and plotting, for example in Jupyter notebooks.

The DPsim simulation kernel is implemented in C++ and uses the Eigen linear algebra library. By using a system programming language like C++ and a highly optimized math library, optimal performance and real-time execution can be guaranteed. The integration into the VILLASframework allows DPsim to be used in large-scale co-simulations.

Licensing

The project is released under the terms of the MPL 2.0.

Where should I go next?

• Getting Started: Get started with DPsim
• Examples: Check out some example code!

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 a function, e.g. getTasks(), in order to retrieve the tasks from instances of the four classes mentioned previously. The power system model element tasks are collected by the Solver instances and relayed to the Simulation.

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

1.3 - Scheduling

DPsim implements level scheduling. A task T4 that depends on data modified by task T1 is scheduled to the level following the level of task T1. In the simplest case, all tasks of a level have to be finished before tasks of the next level are started.

The dependencies of tasks on data are determined by referencing the attributes that are read or modified by the task. The scheduler computes the schedule prior to the simulation from the task dependency graph resulting from the tasks’ data dependencies.

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

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

// 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::initializeFromNodesAndTerminals(Real frequency) {
	//...
	// Create series sub components
	mSubSeriesResistor = std::make_shared<DP::Ph1::Resistor>(**mName + "_res", mLogLevel);
	
    // Setup mSubSeriesResistor...

    // 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);
    //...
}

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

1.6 - Interfaces

Interfaces can be used to share data between a DPsim simulation and other, external services, for example an MQTT-broker. For the purpose of this guide, all services receiving and transmitting data besides the running DPsim instance are grouped in the term environment. Therefore, interfaces provide a way for a DPsim simulation to exchange data with the environment. This data is stored in the form of Attributes and can be imported or exported in every simulation time step. Exporting an attribute means that on every time step, the current value of that attribute is read and written out to the environment. Importing an attribute means that on every time step, a new value is read from the environment and the attribute in the simulation is updated to match this value.

Configuring an Interface

On the configuration level, an interface is an instance of the Interface class. Because the base Interface class requires an instance of an InterfaceWorker to construct, it is recommended to not use this base class directly, but instead construct a subclass derived from Interface which internally handles the construction of the InterfaceWorker. Currently, there exists only one such subclass in DPsim which is the InterfaceVillas.

Configuring the InterfaceVillas

This feature requires the compilation of DPsim with the WITH_VILLAS feature flag. For use of the InterfaceVillas in python, the dpsimpyvillas target has to built in addition to the normal dpsimpy package.

The InterfaceVillas is an interface designed to make use of the various node types and protocols supported by the VILLASframework. By utilizing the nodes provided by VILLASnode, the InterfaceVillas can be configured to import and export attributes from and to a wide range of external services. To create and configure an InterfaceVillas instance, create a new shared pointer of type InterfaceVillas 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 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 InterfaceVillas 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. For an explanation of the various parameters, see the code documentation in InterfaceVillas.h. 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 InterfaceVillas 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);

This method will add two new Tasks to the simulation. The interface’s PreStep task is set to modify all attributes that are imported from the environment and is therefore scheduled to execute before any other simulation tasks that depend on these attributes. The interface’s PostStep task is set to depend on all attributes that are exported to the environment and is therefore scheduled to execute after any other simulation tasks that might modify these attributes. To prevent the scheduler from just dropping the PostStep task since it does not modify any attributes and is therefore not seen as relevant to the simulation, the task is set to modify the Scheduler::external attribute. Note that the execution of these tasks might not necessarily coincide with the point in time at which the values are actually written out to or read from the environment. This is because the interface internally spawns two new threads for exchanging data with the environment and then uses a lock-free queue for communication between these reader and writer threads, and the simulation. Because of this, time-intensive import or export operations will not block the main simulation thread unless this is explicitly configured in the interface’s importAttribute and exportAttribute methods.

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.

2 - Getting Started

2.1 - Build

Docker based

Clone the repository

$ git clone git@github.com:sogno-platform/dpsim.git

or using https if you do not have an account

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

In the repository, there is a Docker file with all required dependencies

$ cd dpsim
$ docker build -t sogno/dpsim:dev -f Packaging/Docker/Dockerfile.dev .

Alternatively, the image can be pulled from DockerHub like so

$ docker pull sogno/dpsim:dev

For OS specific instructions on how to install requirements, see the sections below.

Next, run a Docker container

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

The option -p maps the port 8888 of the container to the docker host. This is required to access the jupyter lab instance inside the container. The option --privileged is required for debug builds.

For Windows, you might need to specify the current directory with curly brackets

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

Now, you should be in an interactive session inside the docker container.

The DPsim C++ and Python library without C++ examples or documentation can be built as follows

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

If you need other libraries that are not built by default, you need to target them specifically, for example if you need `dpsimpy´ and ´dpsimpyvillas´:

$ cmake --build . --target dpsimpy dpsimpyvillas

To build everything run

$ cmake --build .

To use other libraries that are installed, use the relevant option defined in the CMakeList.txt files, for example for GSL below, and then build as usual:

$ cmake .. -DWITH_GSL=ON

If you would like to use the Python package, it has to be added to the path. The following command adds the dpsimpy C++/Python package as well as the dpsim pure Python package.

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

If you are using conda or other ways to develop with environments, please keep in mind that this will become specific for your setup. For this case, from within the environment already active:

$ cd /dpsim/build
$ conda develop $(pwd) && conda develop $(pwd)/Source/Python && conda develop $(pwd)/../Source/Python

To run jupyter lab

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

To install dpsim run

$ cd /dpsim/build
$ sudo make install

CMake for Linux

The most recent list of requirements can be found in the Dockerfiles.

Make sure that the required dependencies are installed. The fedora installation script in the DPsim repository is a good place to start from.

Note: There are currently no Debian packages for villas-node and libcimpp16v29a. If you want to use these optional feature, you have to build them manually.

Install Sundials

$ git clone --branch v3.1.1 https://github.com/LLNL/sundials.git
$ mkdir sundials/build
$ pushd sundials/build
$ cmake .. \
    -DBUILD_SHARED_LIBS=ON \
    -DBUILD_STATIC_LIBS=OFF \
    -DEXAMPLES_ENABLE_C=OFF
$ make -j$(nproc) install
$ popd

The following steps to clone, build and install are the same as for the Docker setup.

CMake for Windows

Make sure that the required dependecies are installed:

• Visual Studio 2017 with C++ Desktop development package
• CMake for Windows
• Git for Windows
• For Python support, install Python3, for example, Anaconda, and add Python to your PATH.

Clone the project as explained for Docker.

Open a windows command prompt and navigate into the new DPsim folder. Generate a Visual Studio project with CMake and use it to build the project

$ mkdir build
$ cd build
$ cmake -G "Visual Studio 15 2017 Win64" ..

Open Visual Studio and load the Visual Studio project from the build directory within the DPsim folder.

You can either build the project from within Visual Studio or from the command line by running the following command in the windows command prompt

$ cmake --build .

To install the Python package use Visual Studio and the Release configuration to build the DPsim Python module and then build the INSTALL project.

CMake for macOS

Make sure that the required dependecies are installed

$ /usr/bin/ruby -e "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/master/install)"
$ brew install gcc9 git cmake graphviz python3 gsl eigen spdlog
$ sudo pip3 install numpy

Clone the source as explained for the Docker setup.

Compile

$ mkdir build
$ cmake ..
$ make -j$(sysctl -n hw.ncpu)

To install the generated Python module to your system

$ sudo make install

Python Package for pypi

Follow the previous steps to set up the Docker container.

To build the Python package run

$ python3 setup.py bdist_wheel

Documentation

Python

Install Sphinx or use the Docker image.

Generate the Python documentation by running Sphinx via CMake:

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

The resulting documentation will be generated in Documentation/html/.

C++

Install Doxygen or use the Docker image.

Generate the C++ documentation by running Doxygen via CMake:

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

The resulting documentation will be generated in Documentation/html/Cxx.

2.2 - Install

DPsim is a Python module / C++ library for complex power system simulation. As a Python module, the easiest way to get to know DPsim is via Jupyter Notebooks.

Docker

First, you need to install Docker. Then, you could either build a docker image by yourself as described in the build instructions or download a prepared image from Docker Hub as described in the following.

To start a Jupyter session, run a DPsim Docker container

	$ docker run -p 8888:8888 sogno/dpsim

And access the session by opening the following link: http://localhost:8888/lab?token=3adaa57df44cea75e60c0169e1b2a98ae8f7de130481b5bc

Python

Currently, the pypi packages are not maintained. Until we have updated the packages, please use the docker installation.

Prerequisites

First, you need to make sure that Python is installed and your version is compatible. An easy way to install Python and all required packages is the Anaconda distribution. To get started, install the latest installer for Python 3.x from the downloads section. Then, run the Anaconda Prompt and create a new conda environment:

$ conda create -n dpsim python=3.6

After creating the environment you need to make sure that it is activated. The current environment is displayed at the beginning of the command line in brackets. It should read "(dpsim)…". In case it is not activated, run:

$ activate dpsim

Pip Package Installation

Then, DPsim can be installed as a Python module:

$ pip install dpsim

From Source

To build and install DPsim from the source files, please refer to the build section.

2.3 - Real-Time

This page describes recommended techniques to optimize the host operating system for real-time execution of DPsim.

In principle, real-time execution is supported on all platforms. However, we recommend to use an optimized Linux installation.

Operating System and Kernel

For minimum latency several kernel and driver settings can be optimized.

To get started, we recommend the Redhat Real-time Tuning Guide.

A PREEMPT_RT patched Linux kernel is recommended. Precompiled kernels for Fedora can be found here: http://ccrma.stanford.edu/planetccrma/software/

Use the tuned tool for improving general real-time performance. Please adjust the setting isolated_cpucores according to your hardware and enable the realtime profile as follows:

  $ dnf install tuned-profiles-realtime
  $ echo "realtime" > /etc/tuned/active_profile
  $ echo "isolated_cpucores=6-7" >> /etc/tuned/realtime-variables.conf
  $ systemctl enable tuned && systemctl start tuned
  $ systemctl reboot

Running a real-time simulation

As a reference, real-time simulation examples are provided in the Examples/Cxx folder of the DPsim repository.

In order to run a real-time simulation, the simulation process must be started in a special way in order to change the execution priority, scheduler and CPU affinity. For this purpose the chrt and taskset commands are used. In the following example, we pin the execution of the simulation to CPU cores 6-7 which have been reserved previously by the tuned real-time profile (see above).

  $ taskset --all-tasks --cpu-list 6-7 \
  $ chrt --fifo 99 \
  $ Examples/Cxx/RT_DP_CS_R_1

More details:

• chrt man-page
• taskset man-page

Recommended Hardware

Some proposals for the selection of appropriate server hardware:

• Server-grade CPU, e.g. Intel Xeon. A multi-core system enables true parallel execution of several decoupled systems
• Server-grade network cards, e.g. Intel PRO/1000. These allow offloading of UDP checksumming to the hardware

3 - Concepts

The book introduces the reader to the general concepts implemented in DPsim, a dynamic phasor (DP) real-time simulator, as well as the physical models of the power system components that are used in simulations. The first chapters give an overview of dynamic phasors and nodal analysis which are the two pillars of the main solver implemented in DPsim. The second part describes in detail what are the physical equations for each model and how they are transformed and implemented for dynamic phasor simulations and other domains that are also supported by DPsim.

In order to be able to run a dynamic simulation, DPsim also includes a loadflow solver to compute the initial state of the network if it is not included in the network data. Besides DP simulations, DPsim also comes with EMT models for some components which are used as reference for testing the DP models.

3.1 - Dynamic Phasors

In the power systems community, dynamic phasors were initially introduced for power electronics analysis Sanders1991 as a more general approach than state-space averaging. They were used to construct efficient models for the dynamics of switching gate phenomena with a high level of detail as shown in Mattavelli1999. A few years later, dynamic phasors were also employed for power system simulation as described in Demiray2008. In Strunz2006 the authors combine the dynamic phasor approach with the Electromagnetic Transients Program (EMTP) simulator concept which includes Modified Nodal Analysis (MNA). Further research topics include fault and stability analysis under unbalanced conditions as presented in Stankovic2000 and also rotating machine models have been developed in dynamic phasors Zhang 2007.

Bandpass Signals and Baseband Representation

Although here, dynamic phasors are presented as a power system modelling tool, it should be noted that the concept is also known in other domains, for example, microwave and communications engineering [Maas2003, Suarez2009, Haykin2009, Proakis2001]. In these domains, the approach is often denoted as base band representation or complex envelope. Another common term coming from power electrical engineering is shifted frequency analysis (SFA). In the following, the general approach of dynamic phasors for power system simulation is explained starting from the idea of bandpass signals. This is because the 50 Hz or 60 Hz fundamental and small deviations from it can be seen as such a bandpass signal. Futhermore, higher frequencies, for example, generated by power electronics can be modelled in a similar way.

3.2 - Nodal Analysis

A circuit with $b$ branches has $2b$ unknowns since there are $b$ voltages and $b$ currents. Hence, $2b$ linear independent equations are required to solve the circuit. If the circuit has $n$ nodes and $b$ branches, it has

• Kirchoff’s current law (KCL) equations
• Kirchoff’s voltage law (KVL) equations
• Characteristic equations (Ohm’s Law)

There are only $n-1$ KCLs since the nth equation is a linear combination of the remaining $n-1$. At the same time, it can be demonstrated that if we can imagine a very high number of closed paths in the network, only $b-n+1$ are able to provide independent KVLs. Finally there are $b$ characteristic equations, describing the behavior of the branch, making a total of $2b$ linear independent equations.

The nodal analysis method reduces the number of equations that need to be solved simultaneously. $n-1$ voltage variables are defined and solved, writing $n-1$ KCL based equations. A circuit can be solved using Nodal Analysis as follows

• Select a reference node (mathematical ground) and number the remaining $n-1$ nodes, that are the independent voltage variables
• Represent every branch current $i$ as a function of node voltage variables $v$ with the general expression $i = g(v)$
• Write $n-1$ KCL based equations in terms of node voltage variable.

The resulting equations can be written in matrix form and have to be solved for $v$. $$\boldsymbol{Y} \boldsymbol{v} = \boldsymbol{i}$$

3.3 - Powerflow

The power flow problem is about the calculation of voltage magnitudes and angles for one set of buses. The solution is obtained from a given set of voltage magnitudes and power levels for a specific model of the network configuration. The power flow solution exhibits the voltages and angles at all buses and real and reactive flows can be deduced from the same.

Power System Model

Power systems are modeled as a network of buses (nodes) and branches (lines). To a network bus, components such a generator, load, and transmission substation can be connected. Each bus in the network is fully described by the following four electrical quantities:

• $\vert V_{k} \vert$: the voltage magnitude
• $\theta_{k}$: the voltage phase angle
• $P_{k}$: the active power
• $Q_{k}$: the reactive power

There are three types of networks buses: VD bus, PV bus and PQ bus. Depending on the type of the bus, two of the four electrical quantities are specified as shown in the table below.

Single Phase Power Flow Problem

The power flow problem can be expressed by the goal to bring a mismatch function $\vec{f}$ to zero. The value of the mismatch function depends on a solution vector $\vec{x}$: $$ \vec{f}(\vec{x}) = 0 $$

As $\vec{f}(\vec{x})$ will be nonlinear, the equation system will be solved with Newton-Raphson:

$$-\textbf{J}(\vec{x}) \Delta \vec{x} = \vec{f} (\vec{x})$$

where $\Delta \vec{x}$ is the correction of the solution vector and $\textbf{J}(\vec{x})$ is the Jacobian matrix. The solution vector $\vec{x}$ represents the voltage $\vec{V}$ by polar or cartesian quantities. The mismatch function $\vec{f}$ will either represent the power mismatch $\Delta \vec{S}$ in terms of

$$\left [ \begin{array}{c} \Delta \vec{P} \ \Delta \vec{Q} \end{array} \right ]$$

or the current mismatch $\Delta \vec{I}$ in terms of

$$\left [ \begin{array}{c} \Delta \vec{I_{real}} \ \Delta \vec{I_{imag}} \end{array} \right ]$$

where the vectors split the complex quantities into real and imaginary parts. Futhermore, the solution vector $\vec{x}$ will represent $\vec{V}$ either by polar coordinates

$$\left [ \begin{array}{c} \vec{\delta} \ \vert \vec{V} \vert \end{array} \right ]$$

or rectangular coordinates

$$\left [ \begin{array}{c} \vec{V_{real}} \ \vec{V_{imag}} \end{array} \right ]$$

This results in four different formulations of the powerflow problem:

• with power mismatch function and polar coordinates
• with power mismatch function and rectangular coordinates
• with current mismatch function and polar coordinates
• with current mismatch function and rectangular coordinates

To solve the problem using NR, we need to formulate $\textbf{J} (\vec{x})$ and $\vec{f} (\vec{x})$ for each powerflow problem formulation.

Powerflow Problem with Power Mismatch Function and Polar Coordinates

Formulation of Mismatch Function

The injected power at a node $k$ is given by: $$S_{k} = V_{k} I _{k}^{*}$$

The current injection into any bus $k$ may be expressed as: $$ I_{k} = \sum_{j=1}^{N} Y_{kj} V_{j} $$

Substitution yields:

We may define $G_{kj}$ and $B_{kj}$ as the real and imaginary parts of the admittance matrix element $Y_{kj}$ respectively, so that $Y_{kj} = G_{kj} + jB_{kj}$. Then we may rewrite the last equation:

If we now perform the algebraic multiplication of the two terms inside the parentheses, and collect real and imaginary parts, and recall that $S_{k} = P_{k} + jQ_{k}$, we can express (1) as two equations: one for the real part, $P_{k}$, and one for the imaginary part, $Q_{k}$, according to:

These equations are called the power flow equations, and they form the fundamental building block from which we solve the power flow problem.

We consider a power system network having $N$ buses. We assume one VD bus, $N_{PV}-1$ PV buses and $N-N_{PV}$ PQ buses. We assume that the VD bus is numbered bus $1$, the PV buses are numbered $2,…,N_{PV}$, and the PQ buses are numbered $N_{PV}+1,…,N$. We define the vector of unknown as the composite vector of unknown angles $\vec{\theta}$ and voltage magnitudes $\vert \vec{V} \vert$:

The right-hand sides of equations (2) and (3) depend on the elements of the unknown vector $\vec{x}$. Expressing this dependency more explicitly, we rewrite these equations as:

We now define the mismatch vector $\vec{f} (\vec{x})$ as:

That is a system of nonlinear equations. This nonlinearity comes from the fact that $P_{k}$ and $Q_{k}$ have terms containing products of some of the unknowns and also terms containing trigonometric functions of some the unknowns.

Formulation of Jacobian

As discussed in the previous section, the power flow problem will be solved using the Newton-Raphson method. Here, the Jacobian matrix is obtained by taking all first-order partial derivates of the power mismatch functions with respect to the voltage angles $\theta_{k}$ and magnitudes $\vert V_{k} \vert$ as:

The linear system of equations that is solved in every Newton iteration can be written in matrix form as follows:

Solution of the Problem

The solution update formula is given by:

To sum up, the NR algorithm, for application to the power flow problem is:

1. Set the iteration counter to $i=1$. Use the initial solution $V_{i} = 1 \angle 0^{\circ}$
2. Compute the mismatch vector $\vec{f}({\vec{x}})$ using the power flow equations
3. Perform the following stopping criterion tests:
   • If $\vert \Delta P_{i} \vert < \epsilon_{P}$ for all type PQ and PV buses and
   • If $\vert \Delta Q_{i} \vert < \epsilon_{Q}$ for all type PQ
   • Then go to step 6
   • Otherwise, go to step 4.
4. Evaluate the Jacobian matrix $\textbf{J}^{(i)}$ and compute $\Delta \vec{x}^{(i)}$.
5. Compute the update solution vector $\vec{x}^{(i+1)}$. Return to step 3.
6. Stop.

4 - Development

Environment

We recommend the following development tools:

• Editor: Visual Studio Code
  • Extensions:
    • C/C++
    • Python
    • CMake Tools
    • Docker
    • EditorConfig for VS Code
    • Remote Development
• Docker
• CMake
• Git

Please follow the build instructions to checkout your code and install the basic dependencies and tools.

4.1 - 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
                }
            ]
        }
    ]
}

4.2 - Guidelines

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.

Creating New Releases

Although DPsim currently does not have any conventions on versioning, the periodic creation of new versions can help to mark significant changes and to analyze new portions of code using tools like SonarCloud.

A new version of DPsim has to be indicated as follows:

• Create a new tag with an increased version number
• Update setup.cfg
• Update CMakeLists.txt
• Update sonar-project.properties

Due to the creation of a new tag, a new PyPi package will be deployed automatically. To release an updated Docker image, the container workflow needs to be triggered manually.

5 - Models

The following models are currently available:

• Dynamic phasors
  • inductor, capacitor, resistor
  • current and voltage source
  • load (PQ and Z type)
  • pi-line
  • transmission line (Bergeron)
  • synchronous generator dq-frame full order (Kundur, Krause)
  • inverter averaged
  • inverter with harmonics (comparable to switched model)
  • switch
• EMT
  • inductor, capacitor, resistor
  • current and voltage source
  • load (Z type)
  • pi-line
  • transmission line (Bergeron)
  • synchronous generator dq-frame full order (Kundur, Krause)
  • inverter averaged
  • switch

5.1 - Transformer

2-Winding Transformer

The transformer model is composed of an RL-segment and an ideal transformer. The single line diagram is depicted in the figure below.

If node reduction is not applied, two virtual nodes are created to stamp this model into the system matrix.

Furthermore, the ideal transformer has an additional equation, which requires an extension of the system matrix. The complete matrix stamp for the ideal transformer is

The variable $j$ denotes the high voltage node while $k$ is the low voltage node. $l$ indicates the inserted row and column to accommodate the relation between the two voltages at the ends of the transformer. The transformer ratio is defined as $T = V_{j} / V_{k}$. A phase shift can be introduced if $T$ is considered as a complex number.

5.2 - Branches

RX-Line

PI-Line

Transformer

5.3 - Induction Machine

5.4 - RLC-Elements

EMT Equations and Modified Nodal Analysis

Inductance

An inductance is described by $$ v_j(t) - v_k(t) = v_L(t) = L \frac{\mathrm{d} i_L(t)}{\mathrm{d}t} $$ Integration results in an equation to compute the current at time $t$ from a previous state at $t - \Delta t$. $$ i_L(t) = i_L(t - \Delta t) + \frac{1}{L} \ \int_{t - \Delta t}^{t} v_L(\tau) \ \mathrm{d} \tau $$ There are various methods to discretize this equation in order to solve it numerically. The trapezoidal rule, an implicit second-order method, is commonly applied for circuit simulation: $$ \int_{t - \Delta t}^{t} f(\tau) \ \mathrm{d} \tau \approx \frac{\Delta t}{2}(f(t) + f(t - \Delta t)) $$ Applying the trapezoidal rule to leads to $$ i_L(t) = i_L(t - \Delta t) + \frac{\Delta t}{2L}(v_L(t) + v_L(t - \Delta t)) $$ This can be rewritten in terms of an equivalent conductance and current source and the number of time steps $k$ with size $\Delta t$. $$ i_L(k) = g_L v_L(k) + i_{L,equiv}(k-1) $$ $$ i_{L,equiv}(k-1) = i_L(k-1) + \frac{\Delta t}{2L} v_L(k-1) $$ $$ g_L = \frac{\Delta t}{2L} $$

Hence, components described by differential equations are transformed into a DC equivalent circuit as depicted in the figure below.

Capacitance

The same procedure can be applied to a capacitance. Integration on both side yields $$ i_C(t) = C \frac{\mathrm{d}}{\mathrm{d}t} \ v_C(t) $$ $$ v_C(t) = v_C(t - \Delta t) + \frac{1}{C} \int_{t - \Delta t}^t i_C(\tau) \mathrm{d} \tau $$ Finally, the equivalent circuit is described by a current source and a conductance. $$ i_{C}(k) = g_{C} v_C(k) + i_{C,equiv}(k-1) $$ $$ i_{C,equiv}(k-1) = -i_{C}(k-1) - g_C v_C(k-1) $$ $$ g_{C} = \frac{2C}{\Delta t} $$ This equation set is visualized in the figure below.

Hence, the vector of unknowns $\bm{x}$ and the source vector $\bm{b}$ become time dependent and this leads to the system description: $$ \bm{A} \bm{x}(t) = \bm{b}(t) $$ To simulate the transient behavior of circuits, this linear equation has to be solved repeatedly. As long as the system topology and the time step is fixed, the system matrix is constant.

Extension with Dynamic Phasors

The dynamic phasor concept can be integrated with nodal analysis. The overall procedure does not change but the system equations are rewritten using complex numbers and all variables need to be expressed in terms of dynamic phasors. Therefore, the resistive companion representations of inductances and capacitances have to be adapted as well.

Inductance

In dynamic phasors the integration of the inductance equation yields

$$ \langle i_L \rangle(t) = \langle i_L \rangle(t - \Delta t) + \int_{t - \Delta t}^t \frac{1}{L} \langle v_L \rangle(\tau) - j \omega \ \langle i_L \rangle(\tau) \mathrm{d} \tau $$

Applying the trapezoidal method leads to the finite difference equation: $$ \begin{split} \langle i_L \rangle(k) = \langle i_L \rangle(k-1) + \frac{\Delta t}{2} \bigg[ \frac{1}{L} (\langle v_L \rangle(k) + \langle v_L \rangle(k-1)) - j \omega (\langle i_L \rangle(t) + \langle i_L \rangle(k-1) \bigg] \end{split} $$

Solving this for $\langle i_L \rangle(k)$ results in the \ac{DP} equivalent circuit model: $$ \langle i_L \rangle(k) = \frac{a - jab}{1 + b^2} \langle v_L \rangle(k) + \langle i_{L,equiv} \rangle(k-1) $$ with $$ a = \frac{\Delta t}{2L}, \qquad b = \frac{\Delta t \omega}{2} $$ $$ \langle i_{L,equiv} \rangle(k-1) = \frac{1 - b^2 - j2b}{1 + b^2} \langle i_L \rangle(k-1) + \frac{a - jab}{1 + b^2} \langle v_L \rangle(k-1) $$

Capacitance

Similarly, a capacitance is described by as follows $$ \langle i_C \rangle(k) = C \ \frac{\mathrm{d} \langle v_C \rangle}{\mathrm{d} t} + j \omega C \ \langle v_C \rangle(t) $$ $$ v_C(t) = v_C(t- \Delta t) + \int_{t- \Delta t}^{t} \frac{1}{C} \ i_C(\tau) -j \omega \ v_C(\tau) \ \mathrm{d} \tau $$

Applying the trapezoidal rule for the capacitance equation leads to the finite difference equation:

The DP model for the capacitance is defined by $$ \langle i_C \rangle(k) = \frac{1+jb}{a} \ \langle v_C \rangle(k) + \langle i_{C,equiv} \rangle(k-1) $$ with $$ a = \frac{\Delta t}{2C}, \qquad b = \frac{\Delta t \omega}{2} $$ $$ \langle i_{C,equiv} \rangle(k-1) = - \frac{1-jb}{a} \ \langle v_C \rangle(k-1) - \langle i_C \rangle(k-1) $$

RL-series element

In dynamic phasors the integration of the inductance equation yields $$ \langle v \rangle(t) = L \frac{\mathrm{d}}{dt} \langle i \rangle(t) + j \omega L \ \langle i \rangle(t) + R \ \langle i \rangle(t) $$ $$ \langle i \rangle(t) = \langle i \rangle(t - \Delta t) + \int_{t - \Delta t}^t \frac{1}{L} \langle v \rangle(\tau) - j \omega \ \langle i \rangle(\tau) - \frac{R}{L} \ \langle i \rangle(\tau) \mathrm{d} \tau $$

Applying the trapezoidal method leads to the finite difference equation: $$ \begin{split} \langle i \rangle(k) = \langle i \rangle(k-1) + \frac{\Delta t}{2} \bigg[ \frac{1}{L} (\langle v \rangle(k) + \langle v \rangle(k-1)) - \left( j \omega + \frac{R}{L} \right) (\langle i \rangle(k) + \langle i \rangle(k-1)) \bigg] \end{split} $$

Solving this for $\langle i \rangle(k)$ results in the \ac{DP} equivalent circuit model: $$ \langle i \rangle(k) = \frac{a + Ra^2 - jab}{(1+Ra)^2 + b^2} \langle v \rangle(k) + \langle i_{equiv} \rangle(k-1) $$ with $$ a = \frac{\Delta t}{2L}, \qquad b = \frac{\Delta t \omega}{2} $$ $$ \langle i_{equiv} \rangle(k-1) = \frac{1 - b^2 - j2b + 2Ra + (Ra)^2 - j2Rab}{(1+Ra^2) + b^2} \langle i \rangle(k-1) + \frac{a + Ra^2 - jab}{(1+Ra)^2 + b^2} \langle v \rangle(k-1) $$

5.5 - Synchronous Generator

Two different synchronous machine models are currently available:

• the full order dq0 reference frame model (EMT, DP) [Kundur, Power system stability and control, 1994]
• and the much simpler transient stability model (DP) [Eremia, Handbook of Electrical Power System Dynamics, 2003]

The machine model is interfaced to the nodal analysis network solver through a current source, which only affects the source vector and not the system matrix Wang2010.

Basic Equations

The equations of the stator and rotor voltages are

where

The flux linkage equations are defined as

The inductance matrices are varying with the rotor position $\theta_r$ which varies with time.

The mechanical equations are:

$\theta_r$ is the rotor position, $\omega_r$ is the angular electrical speed, $P$ is the number of poles, $J$ is the moment of inertia, $T_m$ and $T_e$ are the mechanical and electrical torque, respectively. Motor convention is used for all models.

dq0 Reference Frame 9th Order Model

For stator referred variables, the base quantities for per unit are chosen as follows:

• $v_{s base}$ peak value of rated line-to-neutral voltage in V
• $i_{s base}$ peak value of rated line current in A
• $f_{base}$ rated frequency in Hz

The synchronous generator equations in terms of per unit values in the rotor reference frame become:

where

The flux linkages are:

where

with

The mechanical equations in per unit become:

For the simulation, fluxes are chosen as state variables. To avoid the calculation of currents from fluxes using the inverse of the inductance matrix, the equation set needs to be solved for the fluxes analytically. To simplify the calculations, dq axis magnetizing flux linkages are defined [Krause, Analysis of electric machinery and drive systems, 2002]:

Using the flux linkages results in a simpler equation set for the fluxes:

Dynamic Phasor Model

The fundamental dynamic phasors are similar to the dq0 quantities for symmetrical conditions since both yield DC quantities in a rotating reference frame. The network abc dynamic phasor quantities can be converted to dq0 dynamic phasors by applying the symmetrical components transformation and a rotation.

The angle $\delta$ is the orientation of the dq0 reference frame relative to the abc frame.

The winding currents for positive and zero sequence components can be expressed as

In the dynamic phasor case, the equation for $\frac{d}{dt} \langle \lambda_{0s} \rangle_1$ has a frequency shift. To complete the state model, the magnetizing flux linkages are expressed as:

where

The mechanical equations in dynamic phasors are:

Transient Stability Model

5.6 - VS-Inverter

6 - Core Tasks

Each task should give the user

• The prerequisites for this task, if any (this can be specified at the top of a multi-task page if they’re the same for all the page’s tasks. “All these tasks assume that you understand….and that you have already….”).
• What this task accomplishes.
• Instructions for the task. If it involves editing a file, running a command, or writing code, provide code-formatted example snippets to show the user what to do! If there are multiple steps, provide them as a numbered list.
• If appropriate, links to related concept, tutorial, or example pages.

6.1 - Add New Model

Add a Component Model

In this section we will show the implementation of a new component model by means a three-phase dynamic phasor inductor model.

C++ OOP

DPsim implements component models in a sub project called CPowerSystems (CPS) that is located in the models folder. This folder is added to the DPsim CMake project. Every component in DPsim is represented by a C++ class.

DPsim supports different types of solvers (MNA, DAE, NRP). Each solver requires certain member functions in the component class to be implemented. These functions are specified by the solver interface classes: MNAInterface.h, DAEInterface.h, etc.

Directory / Namespace Structure

For the implementation of the new component, we add two new files

• models/Source/DP/DP_Ph3_Inductor.cpp
• models/Include/DP/DP_Ph3_Inductor.h

In these files, we will implement a new C++ class with the name CPS::DP::Ph3::Inductor. The general structure looks as follows.

Directories:

DPsim
 |
 |- Source
 |- Include
  \ models
      |- Source
          |- DP
          |- EMT
          |- Static
            \ Signal
      |- Include
          |- DP
          |- EMT
          |- Static
            \ Signal

Namespaces:

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

Attributes

Each components has a list of attributes, which has to be specified when creating the components class.

TODO: explain attribute system

Tasks for Pre/Post-step Functions

TODO: add example task dependency graph

Adding the new Component to DPsim

After finishing the implementation of the new component, it needs to be added to the following files:

• models/Include/cps/Components.h
• models/Source/CMakeLists.txt
• Sources/Python/Module.cpp

6.2 - 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 could look like this:

cmake_minimum_required(VERSION 3.5)
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

7 - Examples

The DPsim repository includes examples that can be run locally.

8 - Roadmap

Short-term planning for new features is done on the GitHub Project board.

You can also check the Issues List or the Pull Requests on GitHub.

Under Development

• Solver
  • CUDA sparse implementation
  • improve online system matrix computation and refactorization to support nonlinear elements in network solution (NICSLU integration)
  • merge DAE solver branch
• Interfaces
  • reimplement python interface using pybind and expose more models / functionalities
  • add python based examples using the VILLASnode interface
  • support matpower / pypower format for static simulation
• Tests, Examples, CI
  • convert most of the examples to Python and test them against reference results in CI
  • convert more gitlab CI jobs to github actions
  • add IEEE39 system to examples
• Models
  • VBR generator model
  • SVC
  • add tap-change to transfomer

Ideas

• Solver
  • improve integration of diakoptics solver
• Interfaces
  • implement CIM reader in Python using new pybind interface and cimpy library

9 - Reference

The Sphinx documentation describes the Python API.

The Doxygen documentation is automatically generated from the C++ code using Doxygen. It is helpful to understand the general structure of the C++ DPsim core components.

10 - Contribution Guidelines

We’d love to accept your patches and contributions to this project. Please send us a pull request or get in touch with us via mail or slack if you would like to contribute.