Solvers

• 1: The MNA Solver
• 2: Power Flow Solvers
• 3: Alternative Solver Implementation
• 4: State-Space Extraction

1 - The MNA Solver

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

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

Setting up

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

Which assembly function runs depends on the network:

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

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

Stepping

solve does the same four things every step.

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

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

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

Iterative components

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

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

Linear backends

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

Instrumentation

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

Source

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

2 - Power Flow Solvers

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

Solver Implementations

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

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

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

sim.set_pf_solver_use_sparse(True)

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

Generator Reactive Power Limits

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

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

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

sim.set_pf_solver_enforce_q_limits(True)

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

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

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

3 - Alternative Solver Implementation

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

DAESolver

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

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

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

ODESolver and ODEintSolver

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

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

DiakopticsSolver

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

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

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

Linear backends under MNA

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

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

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

Source

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

4 - State-Space Extraction

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

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

Main classes

The implementation is organized around three main parts:

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

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

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

Usage

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

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

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

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

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

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

Modal analysis of the extracted model

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

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

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

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

Examples

The feature is demonstrated in:

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

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