DPsim is a solver library for dynamic power system simulation.
It supports both the electromagnetic transient (EMT) and dynamic phasor (DP) domain for dynamic simulation.
A powerflow solver is included for standalone usage or to initialize dynamic simulations.
It provides a Python module which can be embedded in any Python 3 application / scripts.
The simulation core is implemented in highly-efficient C++ code.
It supports real-time execution with time-steps down to 50 uS.
It can load models in the IEC61970 Common Information Model (CIM) / Common Grid Model Exchange Standard (CGMES) XML format.
It can be interfaced to a variety of protocols and interfaces via VILLASnode.
Connect
Using or want to use DPsim? Find out more here:
GitHub Discussions - Ask questions, share ideas, and get community support. This is our main point of contact.
LF Energy Slack - Chat with other users and developers and get help in the #sogno or #sogno-dpsim channel.
For email inquiries, please contact the Institute for Automation of Complex Power Systems (ACS), which coordinates DPsim development: post_acs@eonerc.rwth-aachen.de.
How to contribute
If you want to get more involved with DPsim, we welcome contributions of all kinds, including code, documentation, examples, models, bug reports, feature requests, and reviews.
Please open a Pull Request or issue on GitHub, or start a discussion there to propose ideas and get feedback from the community.
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.
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
constAttribute<Real>::PtrmSeriesRes;constAttribute<Real>::PtrmSeriesInd;constAttribute<Real>::PtrmParallelCap;constAttribute<Real>::PtrmParallelCond;// Component constructor: Initializes the attributes in the initialization list
Base::Ph1::PiLine(CPS::AttributeList::PtrattributeList):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(Stringname,Logger::LevellogLevel):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&:
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:
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::PtrA=AttributeDynamic<Real>::make();AttributeBase::PtrB=AttributeDynamic<Real>::make();AttributeBase::PtrC=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::PtrA=AttributeDynamic<Real>::make();AttributeBase::PtrB=AttributeDynamic<Real>::make();AttributeBase::PtrC=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:
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:
autor1=DP::Ph1::Resistor::make("r_1");r1->setParameters(5);autologger=DataLogger::make("simName");// Access the attribute through the member variable
logger->logAttribute("i12",r1->mIntfCurrent);autointf=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:
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:
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 - State-Space Extraction
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 EMT 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.
Supported scope
State-space extraction is available for real-valued EMT simulations with the direct MNA solver.
Supported components with extraction states are:
EMT 3-phase inductors,
EMT 3-phase capacitors,
EMT 3-phase two-terminal V-type SSN components.
Supported algebraic components without extraction states are:
EMT 3-phase resistors,
EMT 3-phase voltage sources.
Other components are rejected explicitly when state-space extraction is enabled.
Usage
In C++, state-space extraction can be enabled on a simulation as follows:
The examples compare equivalent EMT 3-phase RLC systems built from native MNA components and generic two-terminal V-type SSN components.
1.5 - Interfaces
Interfaces can be used to exchange simulation signals between a DPsim simulation and other soft- or hardware, for example an MQTT-broker or an FPGA.
Simulation signals in the form of Attributes can be imported or exported once per simulation time step.
Interfaces are subclasses of Interface and implement the methods addExport and addImport, which add dependencies to the passed attribute that forward the attribute value from or to the interface.
This way, attributes that are imported are read from the interface before they are used in any DPsim component.
Attributes that are exported are written to the interface after they are set by a DPsim component.
Interfacing with VILLASnode
This feature requires the compilation of DPsim with the WITH_VILLAS feature flag. For use of the VILLASnode interface in python, the dpsimpyvillas target has to built in addition to the normal dpsimpy package.
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: InterfaceVillasQueued and InterfaceVillasQueueless.
InterfaceVillasQueued uses a ring buffer to store signal data between DPsim and VILLASnode to allow the protocol used in VILLASnode to operate at a different rate and non-synchronized to the DPsim time step.
InterfaceVillasQueueless uses direct communication with a VILLASnode node type implementing a specific protocol without using a buffer, thus enabling significantly lower latency communication.
With InterfaceVillasQueueless, the protocol operates at the time step of DPsim, i.e., an attribute update directly triggers a write() call to the connected VILLASnode node type.
InterfaceVillasQueued should be used when using non- or soft real-time protocols or communication mediums, such as MQTT or connections via the internet.
InterfaceVillasQueueless should be used when communicating using reliable, low latency, real-time protocols, e.g., with FPGAs, via dedicated fibre networks, or with local real-time applications.
To create and configure one of the VILLASnode interface instance, create a new shared pointer of type InterfaceVillasQueued 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.
Note for InterfaceVillasQueueless:
The queueless interface expects the first input signal in the VILLASnode configuration to be a sequence number that is incremented every time step.
If the value of the sequence number is not incremented by one between two consecutive time steps, an overrun is detected.
Because logging outputs can cause large delays and overruns should not occur spuriously, the interface only reports a warning when a large number of overruns occur.
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::stringmqttConfig=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<InterfaceVillasQueued>intf=std::make_shared<InterfaceVillasQueued>(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 documentationmqtt_config='''{
"type": "mqtt",
"format": "json",
"host": "mqtt",
"in": {
"subscribe": "/mqtt-dpsim"
},
"out": {
"publish": "/dpsim-mqtt"
}
}'''# Creating a new InterfaceVillas objectintf=dpsimpyvillas.InterfaceVillas(name='dpsim-mqtt',config=mqtt_config)# Configuring the InterfaceVillas to import and export attributesintf.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
RealTimeSimulationsim(simName);sim.setSystem(sys);sim.setTimeStep(timeStep);sim.setFinalTime(10.0);// Create and configure interface
autointf=//...
// 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.
1.6 - 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:
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:
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.7 - 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.)
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).
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.
MNA-init stage (mnaCompInitialize()). Unchanged; already recurses into sub-components.
CompositePowerComp<VarType>::initializeFromNodesAndTerminals() is final and sequences these stages:
voidinitializeFromNodesAndTerminals(Realfrequency)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(autosubComp: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.
Value derivation in the parameterization stage must be numerically safe at degenerate inputs. Sub-component values come from terminal power, voltage, and frequency, so a zero rated power, zero reactance, or zero capacitance can produce a division by zero and inject NaN/inf into the system matrix or source vector, which then persists for the rest of the simulation. Guard such expressions: use admittance form (Y = jwC, branch current V * Y) rather than impedance form (V / Z), and create optional shunt branches only when their value is strictly positive — a zero-valued shunt is just an open circuit. A degenerate value that slips through will produce a wrong but non-crashing topology, so an explicit check with a log message is better than assuming the value is nonzero.
// DP_Ph1_PiLine.cpp
DP::Ph1::PiLine::PiLine(Stringuid,Stringname,Logger::LevellogLevel):Base::Ph1::PiLine(mAttributes),// Call the constructor of CompositePowerComp and enable automatic pre- and post-step creation
CompositePowerComp<Complex>(uid,name,true,true,logLevel){//...
}voidDP::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);//...
}voidDP::Ph1::PiLine::initializeParentFromNodesAndTerminals(Realfrequency){//...
// Frequency-dependent values go here, not in createSubComponents().
Realomega=2.*PI*frequency;Compleximpedance={**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:
voidDP::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);}
voidDP::Ph1::PiLine::mnaCompAddPreStepDependencies(AttributeBase::List&prevStepDependencies,AttributeBase::List&attributeDependencies,AttributeBase::List&modifiedAttributes){// Manually add pre-step dependencies of subcomponents
for(autosubComp: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);}
2 - Getting Started
How to install, build and run the DPsim project.
2.1 - Real-Time
This page describes recommended techniques to optimize the host operating system for real-time execution of DPsim.
DPsim supports real-time execution, where the simulation time equals the real or wall clock time, on any system.
However, without modifying system parameters the minimum time step reliably achievable will not be in the range of microseconds.
This is due to operating system noise and other processes interfering with the simulation.
With proper tuning, we have achieved real-time time steps as low as 5 us synchronized to a FPGA using VILLASnode.
Synchronizing the time step to an external source is only necessary, when very high time step accuracy, with maximum deviations in the nanoseconds, is required.
Operating System and Kernel
Using a Linux kernel with the PREEMPT_RT feature improves latency when issuing system calls and enables the FIFO scheduler that lets us avoid preemption during the real-time simulation.
Most distributions offer a binary package for a PREEMPT_RT enabled kernel. For example on Rocky Linux:
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:
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
Running a Real-Time Simulation
Before running a simulation, you can run the following commands as root:
As a reference, real-time simulation examples are provided in the dpsim/examples/cxx and dpsim-villas/examples/cxx folder of the DPsim repository.
To benefit from the PREEMPT_RT feature and the isolated cores, the simulation has to be started using the chrt command to set the scheduling policy and priority, and the taskset command to pin the process to the isolated cores.
In the following example, we set the FIFO scheduling policy with the highest priority (99) and pin the execution of the simulation to CPU cores 9,11,13,15 which have been reserved previously (see above).
# the simple RT_DP_CS_R_1 simulationtaskset -c 9,11,13,15 chrt -f 99 build/dpsim/examples/cxx/RT_DP_CS_R_1
# Cosimulation using VILLASnode, FPGA synchronized time step, and exchanging data via Aurora interface.# Here we need sudo, to interact with the FPGA. We disable logging (log=false) and set the time step to 50 us (-t 0.00005).sudo taskset -c 9,11,13,15 chrt -f 99 build/dpsim-villas/examples/cxx/FpgaCosim3PhInfiniteBus -o log=false -t 0.00005 -d 10
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
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
exportPYTHONPATH=$(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
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.
Building DPsim, including all its dependencies can be done by running:
nix build github:sogno-platform/dpsim
The build result will be available within the result folder of your current directory.
For development purposes, a local development environment can be setup by you running:
nix develop github:sogno-platform/dpsim
Please note, that the Flake reference above (github:sogno-platform/dpsim) can be substituted by a local path (.), in case you have locally checked out the DPsim repo.
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.3 - 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
DPsim can be installed in Linux as a Python module:
pip install dpsim
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.13
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 in Linux:
pip install dpsim
Some examples also need additional packages to process the data or show graphics. Please have a look at the import section of them in case you run into any errors.
From Source
To build and install DPsim from the source files, please refer to the build section.
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 - State-Space Extraction
Discrete-time state-space model is extracted from the EMT MNA simulation model. The extracted model is assembled in the form:
The resulting matrix $\mathbf{A}_{d}$ describes the homogeneous discrete-time dynamics of the EMT MNA simulation model at the current operating point and system-matrix configuration.
References
J. A. Hollman and J. R. Marti, Step-by-step eigenvalue analysis with EMTP discrete-time solutions, IEEE Transactions on Power Systems, 2010. https://doi.org/10.1109/TPWRS.2009.2039810
Y. Han, H. Sun, B. Huang, S. Qin, M. Mu, and Y. Yu, “Discrete-Time State-Space Construction Method for SSO Analysis of Renewable Power Generation Integrated AC/DC Hybrid System,” IEEE Transactions on Power Systems, 2022. https://doi.org/10.1109/TPWRS.2021.3115248
3.2 - 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.3 - 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$.
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.
Bus Type
Known
Unknown
$VD$
$\vert V_{k} \vert, \theta_{k}$
$P_{k}, Q_{k}$
$PV$
$P_{k}, \vert V_{k} \vert$
$Q_{k}, \theta_{k}$
$PQ$
$P_{k}, Q_{k}$
$\vert V_{k} \vert, \theta_{k}$
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:
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
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
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:
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:
Please follow the build instructions to checkout your code and install the basic dependencies and tools.
4.1 - Attribute Usage Guidelines
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:
Does it need to be logged?
Does it need to be imported or exported?
Does it need to be accessed from Python or by name?
Is it used as a scheduler dependency?
Is it an externally relevant model input, output, state, or setpoint?
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:
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.
classMyComponent:publicIdentifiedObject{private:RealmConductance=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.
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.3 - 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.
Pull Requests
There are no strict formal requirements besides the following:
Developer Certificate of Origin (DCO)
We require a Developer Certificate of Origin. See more here.
Code Formatting with pre-commit
We enforce code formatting automatically using pre-commit. Please run pre-commit install the first time you clone the repository to run pre-commit before each commit automatically. If you forgot to do this, you will need to use the command pre-commit run --all-files one time to format your changes.
Development in Forks Only
We accept contributions made in forks only. The main repository is not intended for contributor-specific branches.
SPDX Headers
Use SPDX headers to indicate copyright and licensing information, especially when introducing new files to the codebase. For example:
/* Author: John Smith <John.Smith@example.com>
* SPDX-FileCopyrightText: 2025 Example.com
* SPDX-License-Identifier: MPL-2.0
*/
Creating New Releases (info for maintainers)
DPsim currently uses to Semantic 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:
Update setup.cfg
Update CMakeLists.txt
Update sonar-project.properties
Update CHANGELOG.md and include all the unreleaed changes in the list
Create a new tag with an increased version number (can be done during the Release in GitHub)
Python Packages
Due to the creation of a new tag, a new PyPi package will be deployed automatically.
Only Linux packages are currently available, other platforms will be supported in the future.
Container Images
To release an updated Docker image, the container workflow needs to be triggered manually.
If a Pull Request changes a container image, this is not updated automatically in the container image register.
5 - Models
Mathematical description of the models implemented in DPsim.
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 - Power Electronics
DPsim provides several averaged power-electronic inverter models for simulation using EMT, DP, and SP network modeling domains.
Three-Phase Averaged Voltage Source Inverter with State-Space Nodal Interface
The EMT::Ph3::AvVoltSourceInverterStateSpace model represents a grid-following averaged voltage source inverter in the EMT domain.
It is implemented as a variable state-space nodal component and can therefore be directly stamped into the MNA system.
The model includes a PLL, filtered active/reactive power measurement, outer power control, inner current control, and an LC filter with coupling resistance to the grid node.
In DPSim, synchronous generator control systems are solved separately from the electric network. The outputs of the electric network (active and reactive power, node voltages, branch currents and rotor speed of synchronous generators) at time $k- \Delta t$ are used as the input of the controllers to calculate their states at time $k$. Because of the relatively slow response of the controllers, the error in the network solution due to the time delay $\Delta t$ introduced by this approach is negligible.
Exciter
DC1 type model is the standard IEEE type DC1 exciter, whereas the other model is a simplified version of the IEEE DC1 type model. The inputs of the exciters are the magnitude of the terminal voltage of the generator connected to the exciter $v_h$ and the voltage reference $v_{ref}$, which is defined as a variable since other devices such as over-excitation limiters or power system stabilizers (PSS) modify such reference with additional signals. At the moment, no over-excitation limiters have been implemented in DPSim so that the reference voltage is given by:
$$
v_{ref}(t) = v_{ref,0} + v_{pss}(t)
$$
where $v_{ref,0}$ is initialized after the power flow computations and $v_{pss}(t)$ is the output of the (optional) PSS connected to the exciter. The output of the exciter systems is the induced emf by the field current at $t=k + \Delta t$: $v_{ef}(k + \Delta t)$ (sometimes the alternative notation $e_{fd}(k + \Delta t)$ is used).
IEEE Type DC1 exciter model
Fig. 1: Control diagram of the IEEE Type DC1 exciterAdapted from: Milano, Frequency Variations in Power Systems
This model is used to represent field controlled dc commutator exciters with continuously acting voltage regulators (especially the direct-acting rheostatic, rotating
amplifier, and magnetic amplifier types). The control diagram of this exciter is depicted in Fig. 1 and it is described by the following set of differential equations:
$$
T_{b} \frac{d}{dt} v_{b}(t) = v_{ref} - v_{R}(t) - v_{f}(t) - v_{b}(t),
$$
$$
T_{a} \frac{d}{dt} v_{a}(t) = K_{a} v_{in}(t) - v_{a}(t),
$$
$$
T_{f} \frac{d}{dt} v_{f}(t) - K_{f} \frac{d}{dt} v_{ef}(t) = -v_{f}(t),
$$
$$
T_{ef} \frac{d}{dt} v_{ef}(t) = v_{a}(t) - (K_{ef} + sat(t)) v_{ef}(t),
$$
where $v_h$ is the module of the machine’s terminal voltage, and $v_{in}$ is the amplifier input signal, which for the IEEE Type DC1 is given by:
$$
v_{in}(t) = T_{c} \frac{d}{dt} v_b(t) + v_b(t).
$$
The ceiling function approximates the saturation of the excitation winding:
$$
sat(t) = A_{ef} e^{(B_{ef} | v_{ef}(t) | )}
$$
The set of differential equations are discretized using forward euler in order to solve it numerically, which leads to the following set of algebraic equations:
$$
v_R(k + \Delta t) = v_R(k) + \frac{\Delta t}{T_R} ( v_h(k) - v_R(k) ),
$$
$$
v_b(k + \Delta t) = v_b(k)(1 - \frac{\Delta t}{T_b}) + \frac{\Delta t}{T_b} ( v_{ref}(k) - v_R(k) - v_f(k)),
$$
$$
v_{in}(k + \Delta t) = \Delta t \cdot \frac{T_c}{T_b} (v_{ref}(k) - v_R(k) - v_{f}(k) - v_b(k)) + v_b(k+1),
$$
$$
v_a(k + \Delta t) = v_a(k) + \frac{\Delta t}{T_a} ( v_{in}(k) K_a - v_a(k) ),
$$
$$
v_f(k + \Delta t) = (1 - \frac{\Delta t}{T_f}) v_f(k) + \frac{\Delta t K_f}{T_f T_{ef}} ( v_{a}(k) - (K_{ef} + sat(k)) v_{ef}(k) ),
$$
$$
v_{ef}(k + \Delta t) = v_{ef}(k) + \frac{\Delta t}{T_{ef}} ( v_{a}(k) - (sat(k) + K_{ef}) v_{ef}(k)),
$$
$$
sat(k) = A_{ef} e^{(B_{ef} | v_{ef}(k) | )}
$$
Since the values of all variables for $t=k$ are known, $v_{ef}(k+1)$ can be easily calculated using the discretised equations, which is carried out in the preStep function of the generator connected to each exciter.
The initial values of all variables, which are used in the first simulation step, are calculated assuming that the simulation starts in the steady. This is equivalent to assume that all derivative are equal to zero, which leads to:
$$
v_R(k=0) = v_h(k=0),
$$
$$
v_f(k=0) = 0
$$
$$
v_a(k=0) = K_{ef} v_{ef}(k=0) + A_{ef} e^{B_{ef} |v_{ef} (k=0)|} v_{ef}(k=0),
$$
$$
v_{in}(k=0) = \frac{v_a(k=0)}{K_a},
$$
$$
v_b(k=0) = v_{in}(k=0),
$$
$$
v_{ref}(t=0) = v_{in}(t=0) + v_b(t=0),
$$
where $v_h(k=0)$, $v_{ef}(k=0)$ are calculated after the power flow analysis and after the initialization of synchronous machines (see section initialization of SG).
Simplified IEEE Type DC1 exciter model (DC1Simp)
Fig. 2: Control diagram of the IEEE Type DC1 exciterAdapted from: Milano, Power System Modelling and Scripting
Because the time constants $T_b$ and $T_c$ of the IEEE Type DC1 exciter model are frequently small enough to be neglected, in DPSim a simplified model of this exciter which neglect these time constants is also implemented. The control diagram of this exciter is depicted in Fig. 2 and it is described by the following set of differential equations:
$$
T_R \frac{d}{dt} v_R(t) = v_h(t) - v_R(t)
$$
$$
T_a \frac{d}{dt} v_a(t) = - v_a(t) + K_a v_{in}(t)
$$
$$
T_f \frac{d}{dt} v_f(t) - K_f \frac{d}{dt} v_{ef}(t) = -v_f(t),
$$
$$
T_e \frac{d}{dt} v_{ef}(t) = v_a(t) - v_{ef}(t) (sat(t) + K_{ef})
$$
where $v_h$ is the module of the machine’s terminal voltage, and $v_{in}$ is the amplifier input signal, which is given by:
$$
v_{in}(t) = v_{ref} (t) - v_R(t) - v_f(t)
$$
The set of differential equations are discretized using forward euler in order to solve it numerically, which leads to the following set of algebraic equations:
$$
v_R(k + \Delta t) = v_R(k) + \frac{\Delta t}{T_R} ( v_h(k) - v_R(k) ),
$$
$$
v_{in}(k) = v_{ref}(k) - v_R(k) - v_f(k),
$$
$$
v_a(k + \Delta t) = v_a(k) + \frac{\Delta t}{T_a} ( v_{in}(k) K_a - v_a(k) ),
$$
$$
v_f(k + \Delta t) = (1 - \frac{\Delta t}{T_f}) v_f(k) + \frac{\Delta t K_f}{T_f T_{ef}} ( v_{a}(k) - (K_{ef} + sat(k)) v_{ef}(k) ),
$$
$$
v_{ef}(k + \Delta t) = v_{ef}(k) + \frac{\Delta t}{T_{ef}} ( v_{a}(k) - (sat(k) + K_{ef}) v_{ef}(k)),
$$
$$
sat(k) = A_{ef} e^{(B_{ef} | v_{ef}(k) | )}
$$
Since the values of all variables for $t=k$ are known, $v_{ef}(k+1)$ can be easily calculated using the discretised equations, which is carried out in the preStep function of the generator connected to each exciter.
The initial values of all variables, which are used in the first simulation step, are calculated assuming that the simulation starts in the steady. This is equivalent to assume that all derivative are equal to zero, which leads to:
$$
v_R(k=0) = v_h(k=0),
$$
$$
v_f(k=0) = 0,
$$
$$
v_a(k=0) = K_{ef} v_{ef}(k=0) + A_{ef} e^{B_{ef} |v_{ef} (k=0)|} v_{ef}(k=0),
$$
$$
v_{in}(k=0) = \frac{v_a(k=0)}{K_a},
$$
$$
v_{ref}(t=0) = v_R(t=0) + v_{in}(t=0),
$$
where $v_h(k=0)$, $v_{ef}(k=0)$ are calculated using the power flow analysis and after the initialization of synchronous machines (see section initialization of SG).
Static Exciter
Fig. 3: Control diagram of the Static ExciterAdapted from [6]
The control diagram of this is depicted in Fig. 3. It can be observed as a simplified version of the DC1 type exciter which is composed only by the regulator, the amplifier and an optional transducer. To discretize the lead-lag compensator using forward euler it is better to split this block into two parallel blocks as depicted in Fig. 4.
Fig. 4: Control diagram of the Static Exciter
where:
$$
C_{a} = \frac{T_{a}}{T_{b}}, \quad C_{b} = \frac{T_{b}-T_{a}}{T_{b}}.
$$
and it is described by the following set of differential equations:
$$
T_{R} \frac{d}{dt} v_{r}(t) = v_{h}(t) - v_{r}(t)
$$
$$
T_{b} \frac{d}{dt} x_{b}(t) = v_{in}(t) - x_{b}(t)
$$
$$
T_{e} \frac{d}{dt} e_{fd}(t) = K_{a} v_{e}(t) - e_{fd}(t),
$$
Then, the set of differential equations are discretized using forward euler in order to solve it numerically, which leads to the following set of algebraic equations:
To consider the saturation of $e_{fd}$ there are two different implementations, which is automatically selected depending of value of the parameter $K_{bc}$:
where $e^{*}_{fd}$ represents the output of the exciter.
Anti-windup ($K_{bc}>0$): for controllers with an integral component, i.e. also for PID controllers, the so-called “windup effect” can occur when using the standard saturation function. A strategy for limiting the anti-windup effect is shown in Fig. 5.
Fig. 5: Control diagram of the Static Exciter with anti windup strategy
which means that the input of the differential equation describing $e_{fd}$, $v_{e}$, takes now the following form:
The initial values of all variables, which are used in the first simulation step, are calculated assuming that the simulation starts in the steady. This is equivalent to assume that all derivative are equal to zero, which leads to:
PSS is a controller of synchronous generators used to enhance damping of electromechanical oscillations. The PSS1A implemented in DPSim accepts three optional input signals: rotor speed $\omega$, active power $P$, and terminal voltage magnitude $V_h$. The combined input signal is:
$$
s(t) = K_w \omega(t) + K_p P(t) + K_v V_h(t)
$$
Setting $K_p = K_v = 0$ recovers the speed-only special case. The PSS output $v_{pss}$ at time $t=k$ is a signal used as the input of the AVR to calculate the field voltage at $t=k+\Delta t$, $v_{fd}(k+\Delta t)$. At present, only one PSS is implemented in DPSim which is a simplified version of the IEEE PSS1A type model.
IEEE PSS1A type PSS
Fig. 6: Control diagram of the PSS Type 1 (speed input only;the implementation also accepts active power $K_p P$ and terminal voltage $K_v V_h$).Adapted from: Milano, Power System Modelling and Scripting
The control diagram of this PSS is depicted in Fig. 6. It includes a washout filter and two lead-lag blocks and is described by the following set of differential equations:
$$
T_w \frac{d}{dt} v_1(t) = -(s(t) + v_1(t)),
$$
$$
T_2 \frac{d}{dt} v_2(t) = (1 - \frac{T_1}{T_2})(s(t) + v_1(t)) - v_2(t),
$$
$$
T_4 \frac{d}{dt} v_3(t) = (1 - \frac{T_3}{T_4})\left(v_2(t) + \frac{T_1}{T_2}(s(t) + v_1(t))\right) - v_3(t),
$$
$$
v_{pss}(t) = v_3(t) + \frac{T_3}{T_4}\left(v_2(t) + \frac{T_1}{T_2}(s(t) + v_1(t))\right),
$$
where $s(t) = K_w \omega(t) + K_p P(t) + K_v V_h(t)$ is the combined input signal and $v_{pss}(t)$ is the output signal used to modify the reference voltage of the AVR.
The set of differential equations are discretized using forward euler in order to solve it numerically, which leads to the following set of algebraic equations:
$$
v_1(k + \Delta t) = v_1(k) - \frac{\Delta t}{T_w} (s(k) + v_1(k)),
$$
$$
v_2(k + \Delta t) = v_2(k) + \frac{\Delta t}{T_2} \left((1-\frac{T_1}{T_2})(s(k) + v_1(k)) - v_2(k)\right),
$$
$$
v_3(k + \Delta t) = v_3(k) + \frac{\Delta t}{T_4} \left((1-\frac{T_3}{T_4})\left(v_2(k) + \frac{T_1}{T_2}(s(k) + v_1(k))\right) - v_3(k)\right),
$$
$$
v_{pss}(k) = v_3(k) + \frac{T_3}{T_4} \left(v_2(k) + \frac{T_1}{T_2} (s(k) + v_1(k))\right)
$$
Since the values of all variables for $t=k$ are known, $v_{pss}(k)$ can be easily calculated using the discretised equations, which is carried out in the preStep function of the generator connected to each exciter. Then, $v_{pss}(k)$ is used as input of the AVR to calculate the field voltage at time $k+1$. The values $v_1(k+1)$, $v_2(k+1)$, $v_3(k+1)$ are stored and used to calculate the PSS output of the next time step.
The initial values of all variables, which are used in the first simulation step, are calculated assuming that the simulation starts in steady state. This is equivalent to assuming that all derivatives are equal to zero, which leads to:
$$
v_1(k=0) = -s(k=0),
$$
$$
v_2(k=0) = (1 - \frac{T_1}{T_2})(s(k=0) + v_1(k=0)),
$$
$$
v_3(k=0) = (1 - \frac{T_3}{T_4})\left(v_2(k=0) + \frac{T_1}{T_2}(s(k=0) + v_1(k=0))\right),
$$
$$
v_{pss}(k=0) = v_3(k=0) + \frac{T_3}{T_4}\left(v_2(k=0) + \frac{T_1}{T_2}(s(k=0) + v_1(k=0))\right),
$$
where $s(k=0) = K_w \omega(k=0) + K_p P(k=0) + K_v V_h(k=0)$ is evaluated after the power flow analysis and initialization of synchronous machines (see section initialization of SG). In steady state $\omega(k=0) = 1.0$ (pu), and if $K_p = K_v = 0$ then $v_2 = v_3 = v_{pss} = 0$.
Turbine Governor Models
In DPsim there are two types of Turbine Governor implementations. The Turbine Governor Type 1 implements both the turbine and the governor in one component. In contrast, Steam Turbine and Steam Turbine Governor are implemented as two separate classes and their objects are created independently. Steam/Hydro Turbine and Steam/Hydro Turbine Governor are two blocks that must be connected in series.
The input of the turbine governor models is the mechanical omega at time $t=k-\Delta t$ and the output is the mechanical power at time $t=k$. This variable is then used by the SG to predict the mechanical omega at time $t=k+\Delta t$.
Turbine Governor Type 1
Fig. 7: Control diagram of the turbine governor type 1Source: Milano, Power System Modelling and Scripting
This model includes a governor, a servo and a reheat block. The control diagram of this governor is depicted in Fig. 7 and it is described by the following set of differential equations:
$$
p_{in}(t) = p_{ref} + \frac{1}{R} (\omega_{ref} - \omega(t)),
$$
$$
T_s \frac{d}{dt} x_{g1}(t) = p_{in}(t) - x_{g1}(t),
$$
$$
T_c \frac{d}{dt} x_{g2}(t) = \left(1 - \frac{T_3}{T_c}\right) x_{g1}(t) - x_{g2}(t),
$$
$$
T_5 \frac{d}{dt} x_{g3}(t) = \left(1 - \frac{T_4}{T_5}\right) \left(x_{g2}(t) + \frac{T_3}{T_c} x_{g1}(t)\right) - x_{g3}(t),
$$
$$
\tau_m(t) = x_{g3}(t) + \frac{T_4}{T_5} \left(x_{g2}(t) + \frac{T_3}{T_c} x_{g1}(t)\right),
$$
where $\omega(t)$ is the input signal and $\tau_m(t)$ is the output signal of the governor.
The differential equations are discretized using the forward Euler method, which leads to the following set of algebraic equations:
$$
p_{in}(k-\Delta t) = p_{ref} + \frac{1}{R} (\omega_{ref} - \omega(k-\Delta t)),
$$
$$
x_{g1}(k) = x_{g1}(k-\Delta t) + \frac{\Delta t}{T_s} \left(p_{in}(k-\Delta t) - x_{g1}(k-\Delta t)\right),
$$
$$
x_{g2}(k) = x_{g2}(k-\Delta t) + \frac{\Delta t}{T_c} \left(\left(1 - \frac{T_3}{T_c}\right) x_{g1}(k-\Delta t) - x_{g2}(k-\Delta t)\right),
$$
$$
x_{g3}(k) = x_{g3}(k-\Delta t) + \frac{\Delta t}{T_5} \left(\left(1 - \frac{T_4}{T_5}\right) \left(x_{g2}(k-\Delta t) + \frac{T_3}{T_c} x_{g1}(k-\Delta t)\right) - x_{g3}(k-\Delta t)\right),
$$
$$
\tau_m(k) = x_{g3}(k) + \frac{T_4}{T_5} \left(x_{g2}(k) + \frac{T_3}{T_c} x_{g1}(k)\right).
$$
Since all variables at $t=k-\Delta t$ are known, $\tau_m(k)$ is computed in the preStep of the generator and used to approximate the mechanical equations at time $k+\Delta t$.
Steam Governor
Fig. 8: Control diagram of the steam turbine governorAdapted from [6]
The control diagram of this model is depicted in Fig. 8. This model receives as input the frequency deviation $\Delta\omega = \omega_{ref} - \omega$ from the nominal frequency (normally $50,\text{Hz}$ or $60,\text{Hz}$) and produces the valve opening signal $p_{gv}$ for the turbine. $p_{ref}$ is the mechanical power produced at nominal frequency. The governor implements a lead-lag controller $\frac{K(1+sT_2)}{(1+sT_1)}$ where $K=1/R$ and $R$ is the droop coefficient, followed by a PT1 integrator with embedded rate limiters and an anti-windup loop. To avoid unnecessary dead-beat behaviour, complex transfer functions with more than one pole and zero are decomposed via partial fraction expansion into parallel PT1 elements, as shown in Fig. 9.
Fig. 9: Control diagram of the steam turbine governor after partial-fraction decomposition
Analogous to the static exciter model, the integrator uses an anti-windup strategy as shown in Fig. 10.
Fig. 10: Control diagram of the steam turbine governor with anti-windup strategy
If $T_1 = 0$ the $p_1(k)$ equation is skipped and $p(k)$ is instead:
$$
p(k-\Delta t) = \frac{1}{R} \left(\Delta \omega(k-\Delta t) + \frac{T_{2}}{\Delta t} \left(\Delta \omega(k-\Delta t) - \Delta \omega(k-2\Delta t)\right)\right).
$$
Assuming the simulation starts in steady state (all derivatives zero, $\Delta\omega(0)=0$), the initial values are:
$$
p_{1}(t=0) = 0, \quad p(t=0) = 0, \quad p_{ref} = p_{gv}^{*}(t=0) = p_{gv}(t=0).
$$
Steam Turbine
Fig. 11: Control diagram of the steam turbineAdapted from [6]
The steam turbine receives the valve opening signal $p_{gv}$ from the Steam Governor and outputs mechanical power $p_m$ to the synchronous generator. It is divided into high-pressure (HP), intermediate-pressure (IP), and low-pressure (LP) stages, each modelled as a first-order lag with time constants $T_{CH}$, $T_{RH}$, $T_{CO}$ respectively. Setting a time constant to zero disables that lag element. The total mechanical power is a weighted sum of each stage: $F_{HP} + F_{IP} + F_{LP} = 1$ must hold. The forward-Euler discretised equations are:
Fig. 12: Control diagram of a hydro turbine governorAdapted from [6]
The Hydro Turbine Governor receives the frequency deviation $\Delta\omega = \omega_{ref} - \omega$ as input and produces the valve/gate opening signal $p_{gv}$ for the turbine. $p_{ref}$ is the mechanical power produced at nominal frequency. The controller transfer function is $K\frac{1+sT_2}{(1+sT_1)(1+sT_3)}$, where $K=\frac{1}{R}$ and $R$ is the droop coefficient. The transfer function is decomposed into two parallel PT1 blocks as shown in Fig. 13.
Fig. 13: Control diagram of a hydro turbine governor after partial-fraction decomposition
Assuming the simulation starts in steady state (all derivatives zero, $\Delta\omega(t=0)=0$), the initial values are:
$$
x_{1}(t=0) = 0, \quad x_{2}(t=0) = 0, \quad p_{ref} = p_{gv}(t=0).
$$
Hydro Turbine
Fig. 14: Control diagram of a hydro turbineAdapted from [6]
The Hydro Turbine receives the gate opening signal $p_{gv}$ from the Hydro Turbine Governor and outputs mechanical power $p_m$ to the synchronous generator. The transfer function is specified by the water starting time $T_W$ and can be represented as the sum of two parallel blocks as shown in Fig. 15.
Fig. 15: Control diagram of a hydro turbine after decompositionAdapted from [6]
Assuming the simulation starts in steady state (all derivatives zero), the initial values are:
$$
x_{1}(t=0) = p_{gv}(t=0), \quad p_{m}(t=0) = p_{gv}(t=0).
$$
References
[1] “IEEE Recommended Practice for Excitation System Models for Power System Stability Studies,” in IEEE Std 421.5-2016 (Revision of IEEE Std 421.5-2005) , vol., no., pp.1-207, 26 Aug. 2016, doi: 10.1109/IEEESTD.2016.7553421.
[2] F. Milano, “Power system modelling and scripting,” in Power System Modelling and Scripting. London: Springer-Verlag, 2010, ISBN: 978-3-642-13669-6. doi: 10.1007/978-3-642-13669-6.
[3] F. Milano, A. Manjavacas, “Frequency Variations in Power Systems: Modeling, State Estimation, and Control”. ISBN: 978-1-119-55184-3.
[4] F. Milano, “Power System Analysis Toolbox: Documentation for PSAT”, ISBN: 979-8573500560.
[6] A. Roehder, B. Fuchs, J. Massman, M. Quester, A. Schnettler, “Transmission system stability assessment within an integrated grid development process”.
5.3 - 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
$$\begin{array}{c|c c c}
~ & j & k & l \cr
\hline
j & & & -1 \cr
k & & & T \cr
l & 1 & -T & 0
\end{array}
\begin{pmatrix}
v_j \cr
v_k \cr
i_{l} \cr
\end{pmatrix}
=
\begin{pmatrix}
\cr
\cr
0\cr
\end{pmatrix}$$
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.4 - Branches
RX-Line
PI-Line
Transformer
5.5 - Induction Machine
5.6 - 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$.
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:
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
$$\begin{align}
\langle v_L \rangle(t) &= \Big \langle L \frac{\mathrm{d} i_L}{\mathrm{d}t} \Big \rangle(t) \nonumber \\
&= L \frac{\mathrm{d}}{dt} \langle i_L \rangle(t) + j \omega L \ \langle i_L \rangle(t)
\end{align}$$
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
$\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:
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]:
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.
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:
Description of typical simulation and development 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
Extending the simulator with new component or control models.
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
Using DPsim for a new simulation scenario.
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:
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
How to contribute to DPsim.
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.
