1 - State-Space Nodal
The state-space nodal (SSN) method represents a component by its own continuous
state-space model and couples it to the network through the nodal admittance
matrix. A component is described by
$$\frac{d\mathbf{x}}{dt} = \mathbf{A}\mathbf{x} + \mathbf{B}\mathbf{u},
\qquad
\mathbf{y} = \mathbf{C}\mathbf{x} + \mathbf{D}\mathbf{u},$$
where $\mathbf{x}$ is the internal state, and the input $\mathbf{u}$ and output
$\mathbf{y}$ are the terminal quantities exchanged with the network (a voltage and
the corresponding current). Trapezoidal discretisation of $(\mathbf{A}, \mathbf{B})$
yields a discrete model $(\mathbf{A}_d, \mathbf{B}_d)$ and a Norton equivalent: a
constant conductance $\mathbf{W}$ stamped into the system matrix plus a history
current source recomputed each step from the previous state and input. Because the
component is solved simultaneously with the network in the same nodal system, SSN
is numerically robust without the parasitic snubbers that delayed
current-injection schemes require.
This builds directly on Nodal Analysis and is
the companion of
State-Space Extraction, which
recovers a state-space model from an MNA simulation rather than starting from one.
Shift to the Dynamic-Phasor Envelope
In the Dynamic Phasors (shifted-frequency)
domain a physical quantity is written as a slowly varying complex envelope on a
carrier $\omega_s$,
$$x(t) = \operatorname{Re}\{\, X(t)\, e^{j \omega_s t} \,\},$$
so that differentiation maps to
$\frac{dx}{dt} \rightarrow \frac{dX}{dt} + j \omega_s X$. The physical
state-space model therefore becomes, for the envelope,
$$\frac{d\mathbf{X}}{dt} = (\mathbf{A} - j \omega_s \mathbf{I})\,\mathbf{X}
+ \mathbf{B}\,\mathbf{U}.$$
DPsim does not solve the dynamic-phasor system as a complex matrix: every
complex admittance $g = g_r + j g_i$ is stamped as a real $2 \times 2$ block
$\left[\begin{smallmatrix} g_r & -g_i \ g_i & g_r \end{smallmatrix}\right]$, with the real
and imaginary node parts living in separate halves of a real-valued system. The
DP-SSN model matches this by carrying the envelope in real/imaginary-split form
$[\mathbf{X}_r; \mathbf{X}_i]$ and writing the carrier shift explicitly as a real
augmented model,
$$\mathbf{A}_\text{aug} =
\begin{bmatrix} \mathbf{A} & \omega_s \mathbf{I} \\ -\omega_s \mathbf{I} & \mathbf{A} \end{bmatrix},
\qquad
\mathbf{B}_\text{aug} = \begin{bmatrix} \mathbf{B} & \mathbf{0} \\ \mathbf{0} & \mathbf{B} \end{bmatrix}.$$
The same real trapezoidal discretisation used by the EMT SSN models then applies
unchanged, the Norton conductance comes out as a real block that stamps directly
into the augmented network, and the interface reads and writes the real and
imaginary node halves directly with no complex assembly. The $\pm \omega_s$
coupling block is algebraically exact; placing a rotating-frame term inside the
real state matrix follows the established SSN practice for rotating-machine
models.
Components
The single-phase dynamic-phasor SSN models are:
Full_Serial_RLC, a series resistor-inductor-capacitor one-port with a
hand-derived state-space model, used as the reference component.GenericTwoTerminalVTypeSSN and GenericTwoTerminalITypeSSN, which accept a
user-supplied $(\mathbf{A}, \mathbf{B}, \mathbf{C}, \mathbf{D})$ and build the
V-type (voltage input, current output) or I-type (current input, voltage
output) stamping accordingly.
All three reproduce the classical dynamic-phasor stamping of the same circuit,
and the reconstructed time-domain waveform matches the EMT and EMT-SSN results
within discretisation error.
Three-Phase Components
The same real-augmented model extends per phase to DP::Ph3. The
$3 \times 3$ $\mathbf{A}$, $\mathbf{B}$, $\mathbf{C}$, $\mathbf{D}$ matrices
are general, so off-diagonal entries can couple the phases together:
Full_Serial_RLC, the three-phase series RLC one-port.GenericTwoTerminalVTypeSSN and GenericTwoTerminalITypeSSN, the
three-phase generic V-type and I-type components.
As in the single-phase case, all three reproduce the classical three-phase
dynamic-phasor stamping exactly, and the reconstructed time-domain waveform
matches the EMT and EMT-SSN results within discretisation error once
corrected for the RMS-to-peak scaling that EMT::Ph3 sources apply and
DP::Ph3 sources do not, since the DP envelope already carries the complex
amplitude directly. The notebooks below only exercise the symmetrical,
diagonal case; coupling between phases is not covered by existing tests.
Validation and Examples
Two notebooks accompany the single-phase models, both on a single-carrier
series RLC one-port. examples/Notebooks/Circuits/DP_generalizedSSN_RLC.ipynb
validates the DP-SSN models against the classical dynamic-phasor stamping and
the EMT and EMT-SSN waveforms. examples/Notebooks/Circuits/DP_SSN_RLC_accuracy.ipynb
studies the time-step and frequency-dependent accuracy against a small-step
EMT reference.
The three-phase analogues,
examples/Notebooks/Circuits/DP_Ph3_generalizedSSN_RLC.ipynb and
examples/Notebooks/Circuits/DP_Ph3_SSN_RLC_accuracy.ipynb, repeat both
studies on DP::Ph3 circuits, including a current-driven network with generic
V-type and I-type components and a three-phase fault transient.
Further Reading
- C. Dufour, J. Mahseredjian, and J. Bélanger, A Combined State-Space Nodal Method for the Simulation of Power System Transients, IEEE Transactions on Power Delivery, vol. 26, no. 2, pp. 928–935, 2011. https://doi.org/10.1109/TPWRD.2010.2090364
- C. Dufour and D. S. Nasrallah, State-space-nodal rotating machine models with improved numerical stability, IECON 2016 – 42nd Annual Conference of the IEEE Industrial Electronics Society, 2016. https://doi.org/10.1109/IECON.2016.7793690
- A. A. Kida, A. C. S. Lima, F. A. Moreira, J. R. Martí, and J. Tarazona, Inaccuracies due to the frequency warping in simulation of electrical systems using combined state–space nodal analysis, Electric Power Systems Research, vol. 223, art. 109657, 2023. https://doi.org/10.1016/j.epsr.2023.109657
2 - State-Space Extraction
Discrete-time state-space model is extracted from the EMT MNA simulation model. The extracted model is assembled in the form:
$$\mathbf{x}[k+1]
=
\mathbf{A}_{d,\mathrm{local}} \mathbf{x}[k]
+
\mathbf{B}_{d,\mathrm{MNA}} \mathbf{x}_{\mathrm{MNA}}[k+1]$$
$$\mathbf{Y} \mathbf{x}_{\mathrm{MNA}}[k+1]
=
\mathbf{C}_{d,\mathrm{MNA}} \mathbf{x}[k]$$
where:
- $\mathbf{x}$ is the extraction-state vector,
- $\mathbf{x}_{\mathrm{MNA}}$ is the full MNA unknown vector,
- $\mathbf{Y}$ is the active MNA system matrix,
- $\mathbf{A}_{d,\mathrm{local}}$ contains local component state-transition contributions,
- $\mathbf{B}_{d,\mathrm{MNA}}$ maps MNA unknowns to the state update,
- $\mathbf{C}_{d,\mathrm{MNA}}$ maps extraction states to MNA current injections.
Eliminating the MNA unknown vector gives the global discrete-time state matrix
$$\mathbf{A}_{d}
=
\mathbf{A}_{d,\mathrm{local}}
+
\mathbf{B}_{d,\mathrm{MNA}}
\operatorname{solve}
\left(
\mathbf{Y},
\mathbf{C}_{d,\mathrm{MNA}}
\right)$$
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
4 - Nodal Analysis
A circuit with $b$ branches has $2b$ unknowns since there are $b$ voltages and $b$ currents.
Hence, $2b$ linear independent equations are required to solve the circuit.
If the circuit has $n$ nodes and $b$ branches, it has
- Kirchoff’s current law (KCL) equations
- Kirchoff’s voltage law (KVL) equations
- Characteristic equations (Ohm’s Law)
There are only $n-1$ KCLs since the nth equation is a linear combination of the remaining $n-1$.
At the same time, it can be demonstrated that if we can imagine a very high number of closed paths in the network, only $b-n+1$ are able to provide independent KVLs.
Finally there are $b$ characteristic equations, describing the behavior of the branch, making a total of $2b$ linear independent equations.
The nodal analysis method reduces the number of equations that need to be solved simultaneously.
$n-1$ voltage variables are defined and solved, writing $n-1$ KCL based equations.
A circuit can be solved using Nodal Analysis as follows
- Select a reference node (mathematical ground) and number the remaining $n-1$ nodes, that are the independent voltage variables
- Represent every branch current $i$ as a function of node voltage variables $v$ with the general expression $i = g(v)$
- Write $n-1$ KCL based equations in terms of node voltage variable.
The resulting equations can be written in matrix form and have to be solved for $v$.
$$\boldsymbol{Y} \boldsymbol{v} = \boldsymbol{i}$$
5 - Powerflow
The power flow problem is about the calculation of voltage magnitudes and angles for one set of buses.
The solution is obtained from a given set of voltage magnitudes and power levels for a specific model of the network configuration.
The power flow solution exhibits the voltages and angles at all buses and real and reactive flows can be deduced from the same.
Power System Model
Power systems are modeled as a network of buses (nodes) and branches (lines).
To a network bus, components such a generator, load, and transmission substation can be connected.
Each bus in the network is fully described by the following four electrical quantities:
- $\vert V_{k} \vert$: the voltage magnitude
- $\theta_{k}$: the voltage phase angle
- $P_{k}$: the active power
- $Q_{k}$: the reactive power
There are three types of networks buses: VD bus, PV bus and PQ bus.
Depending on the type of the bus, two of the four electrical quantities are specified as shown in the table below.
| 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:
$$-\textbf{J}(\vec{x}) \Delta \vec{x} = \vec{f} (\vec{x})$$
where $\Delta \vec{x}$ is the correction of the solution vector and $\textbf{J}(\vec{x})$ is the Jacobian matrix.
The solution vector $\vec{x}$ represents the voltage $\vec{V}$ by polar or cartesian quantities.
The mismatch function $\vec{f}$ will either represent the power mismatch $\Delta \vec{S}$ in terms of
$$\left [ \begin{array}{c} \Delta \vec{P} \\ \Delta \vec{Q} \end{array} \right ]$$
or the current mismatch $\Delta \vec{I}$ in terms of
$$\left [ \begin{array}{c} \Delta \vec{I_{real}} \\ \Delta \vec{I_{imag}} \end{array} \right ]$$
where the vectors split the complex quantities into real and imaginary parts.
Futhermore, the solution vector $\vec{x}$ will represent $\vec{V}$ either by polar coordinates
$$\left [ \begin{array}{c} \vec{\delta} \\ \vert \vec{V} \vert \end{array} \right ]$$
or rectangular coordinates
$$\left [ \begin{array}{c} \vec{V_{real}} \\ \vec{V_{imag}} \end{array} \right ]$$
This results in four different formulations of the powerflow problem:
- with power mismatch function and polar coordinates
- with power mismatch function and rectangular coordinates
- with current mismatch function and polar coordinates
- with current mismatch function and rectangular coordinates
To solve the problem using NR, we need to formulate $\textbf{J} (\vec{x})$ and $\vec{f} (\vec{x})$ for each powerflow problem formulation.
Of these four, DPsim currently implements the power mismatch function with polar
coordinates, detailed below. It is the formulation used by both the dense
(PFSolverPowerPolar) and sparse (PFSolverPowerPolarSparse) solvers; the other three
are not implemented.
Powerflow Problem with Power Mismatch Function and Polar Coordinates
The injected power at a node $k$ is given by:
$$S_{k} = V_{k} I _{k}^{*}$$
The current injection into any bus $k$ may be expressed as:
$$I_{k} = \sum_{j=1}^{N} Y_{kj} V_{j}$$
Substitution yields:
$$\begin{align}
S_{k} &= V_{k} \left ( \sum_{j=1}^{N} Y_{kj} V_{j} \right )^{*} \nonumber \\
&= V_{k} \sum_{j=1}^{N} Y_{kj}^{*} V_{j} ^{*} \nonumber
\end{align}$$
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:
$$\begin{align}
S_{k} &= V_{k} \sum_{j=1}^{N} Y_{kj}^{*} V_{j}^{*} \nonumber \\
&= \vert V_{k} \vert \angle \theta_{k} \sum_{j=1}^{N} (G_{kj} + jB_{kj})^{*} ( \vert V_{j} \vert \angle \theta_{j})^{*} \nonumber \\
&= \vert V_{k} \vert \angle \theta_{k} \sum_{j=1}^{N} (G_{kj} - jB_{kj}) ( \vert V_{j} \vert \angle - \theta_{j}) \nonumber \\
&= \sum_{j=1} ^{N} \vert V_{k} \vert \angle \theta_{k} ( \vert V_{j} \vert \angle - \theta_{j}) (G_{kj} - jB_{kj}) \nonumber \\
&= \sum_{j=1} ^{N} \left ( \vert V_{k} \vert \vert V_{j} \vert \angle (\theta_{k} - \theta_{j}) \right ) (G_{kj} - jB_{kj}) \nonumber \\
&= \sum_{j=1} ^{N} \vert V_{k} \vert \vert V_{j} \vert \left ( cos(\theta_{k} - \theta_{j}) + jsin(\theta_{k} - \theta_{j}) \right ) (G_{kj} - jB_{kj})
\end{align}$$
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:
$$\begin{align}
{P}_{k} = \sum_{j=1}^{N} \vert V_{k} \vert \vert V_{j} \vert \left ( G_{kj}cos(\theta_{k} - \theta_{j}) + B_{kj} sin(\theta_{k} - \theta_{j}) \right ) \\
{Q}_{k} = \sum_{j=1}^{N} \vert V_{k} \vert \vert V_{j} \vert \left ( G_{kj}sin(\theta_{k} - \theta_{j}) - B_{kj} cos(\theta_{k} - \theta_{j}) \right )
\end{align}$$
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$:
$$\begin{align}
\vec{x} = \left[ \begin{array}{c} \vec{\theta} \\ \vert \vec{V} \vert \\ \end{array} \right ]
= \left[ \begin{array}{c} \theta_{2} \\ \theta_{3} \\ \vdots \\ \theta_{N} \\ \vert V_{N_{PV+1}} \vert \\ \vert V_{N_{PV+2}} \vert \\ \vdots \\ \vert V_{N} \vert \end{array} \right]
\end{align}$$
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:
$$\begin{align}
P_{k} = P_{k} (\vec{x}) \Rightarrow P_{k}(\vec{x}) - P_{k} &= 0 \quad \quad k = 2,...,N \\
Q_{k} = Q_{k} (\vec{x}) \Rightarrow Q_{k} (\vec{x}) - Q_{k} &= 0 \quad \quad k = N_{PV}+1,...,N
\end{align}$$
We now define the mismatch vector $\vec{f} (\vec{x})$ as:
$$\begin{align}
\vec{f} (\vec{x}) = \left [ \begin{array}{c} f_{1}(\vec{x}) \\ \vdots \\ f_{N-1}(\vec{x}) \\ ------ \\ f_{N}(\vec{x}) \\ \vdots \\ f_{2N-N_{PV} -1}(\vec{x}) \end{array} \right ]
= \left [ \begin{array}{c} P_{2}(\vec{x}) - P_{2} \\ \vdots \\ P_{N}(\vec{x}) - P_{N} \\ --------- \\ Q_{N_{PV}+1}(\vec{x}) - Q_{N_{PV}+1} \\ \vdots \\ Q_{N}(\vec{x}) - Q_{N} \end{array} \right]
= \left [ \begin{array}{c} \Delta P_{2} \\ \vdots \\ \Delta P_{N} \\ ------ \\ \Delta Q_{N_{PV}+1} \\ \vdots \\ \Delta Q_{N} \end{array} \right ]
= \vec{0}
\end{align}$$
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.
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:
$$\begin{align}
J_{jk}^{P \theta} &= \frac{\partial P_{j} (\vec{x} ) } {\partial \theta_{k}} = \vert V_{j} \vert \vert V_{k} \vert \left ( G_{jk} sin(\theta_{j} - \theta_{k}) - B_{jk} cos(\theta_{j} - \theta_{k} ) \right ) \\
J_{jj}^{P \theta} &= \frac{\partial P_{j}(\vec{x})}{\partial \theta_{j}} = -Q_{j} (\vec{x} ) - B_{jj} \vert V_{j} \vert ^{2} \\
J_{jk}^{Q \theta} &= \frac{\partial Q_{j}(\vec{x})}{\partial \theta_{k}} = - \vert V_{j} \vert \vert V_{k} \vert \left ( G_{jk} cos(\theta_{j} - \theta_{k}) + B_{jk} sin(\theta_{j} - \theta_{k}) \right ) \\
J_{jj}^{Q \theta} &= \frac{\partial Q_{j}(\vec{x})}{\partial \theta_{k}} = P_{j} (\vec{x} ) - G_{jj} \vert V_{j} \vert ^{2} \\
J_{jk}^{PV} &= \frac{\partial P_{j} (\vec{x} ) } {\partial \vert V_{k} \vert } = \vert V_{j} \vert \left ( G_{jk} cos(\theta_{j} - \theta_{k}) + B_{jk} sin(\theta_{j} - \theta_{k}) \right ) \\
J_{jj}^{PV} &= \frac{\partial P_{j}(\vec{x})}{\partial \vert V_{j} \vert } = \frac{P_{j} (\vec{x} )}{\vert V_{j} \vert} + G_{jj} \vert V_{j} \vert \\
J_{jk}^{QV} &= \frac{\partial Q_{j} (\vec{x} ) } {\partial \vert V_{k} \vert } = \vert V_{j} \vert \left ( G_{jk} sin(\theta_{j} - \theta_{k}) + B_{jk} cos(\theta_{j} - \theta_{k}) \right ) \\
J_{jj}^{QV} &= \frac{\partial Q_{j}(\vec{x})}{\partial \vert V_{j} \vert } = \frac{Q_{j} (\vec{x} )}{\vert V_{j} \vert} - B_{jj} \vert V_{j} \vert \\
\end{align}$$
The formulas above use the voltage magnitude $\vert V_k \vert$ as the unknown. The DPsim
implementation instead uses the relative voltage increment $\Delta \vert V_k \vert / \vert V_k \vert$,
which scales every voltage-magnitude column ($J^{PV}$, $J^{QV}$) by $\vert V_k \vert$. For
example $J_{jj}^{PV}$ becomes $P_j(\vec{x}) + G_{jj} \vert V_j \vert^2$. This is paired with
the multiplicative voltage update $\vert V_k \vert \leftarrow \vert V_k \vert (1 + \Delta \vert V_k \vert / \vert V_k \vert)$,
so the solution is identical; only the scaling of the voltage columns differs.
The linear system of equations that is solved in every Newton iteration can be written in matrix form as follows:
$$\begin{align}
-\left [ \begin{array}{cccccc}
\frac{\partial \Delta P_{2} }{\partial \theta_{2}} & \cdots & \frac{\partial \Delta P_{2} }{\partial \theta_{N}} &
\frac{\partial \Delta P_{2} }{\partial \vert V_{N_{G+1}} \vert} & \cdots & \frac{\partial \Delta P_{2} }{\partial \vert V_{N} \vert} \\
\vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\
\frac{\partial \Delta P_{N} }{\partial \theta_{2}} & \cdots & \frac{\partial \Delta P_{N}}{\partial \theta_{N}} &
\frac{\partial \Delta P_{N}}{\partial \vert V_{N_{G+1}} \vert } & \cdots & \frac{\partial \Delta P_{N}}{\partial \vert V_{N} \vert} \\
\frac{\partial \Delta Q_{N_{G+1}} }{\partial \theta_{2}} & \cdots & \frac{\partial \Delta Q_{N_{G+1}} }{\partial \theta_{N}} &
\frac{\partial \Delta Q_{N_{G+1}} }{\partial \vert V_{N_{G+1}} \vert } & \cdots & \frac{\partial \Delta Q_{N_{G+1}} }{\partial \vert V_{N} \vert} \\
\vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\
\frac{\partial \Delta Q_{N}}{\partial \theta_{2}} & \cdots & \frac{\partial \Delta Q_{N}}{\partial \theta_{N}} &
\frac{\partial \Delta Q_{N}}{\partial \vert V_{N_{G+1}} \vert } & \cdots & \frac{\partial \Delta Q_{N}}{\partial \vert V_{N} \vert}
\end{array} \right ]
\left [ \begin{array}{c} \Delta \theta_{2} \\ \vdots \\ \Delta \theta_{N} \\ \Delta \vert V_{N_{G+1}} \vert \\ \vdots \\ \Delta \vert V_{N} \vert \end{array} \right ]
= \left [ \begin{array}{c} \Delta P_{2} \\ \vdots \\ \Delta P_{N} \\ \Delta Q_{N_{G+1}} \\ \vdots \\ \Delta Q_{N} \end{array} \right ]
\end{align}$$
Solution of the Problem
The solution update formula is given by:
$$\begin{align}
\vec{x}^{(i+1)} = \vec{x}^{(i)} + \Delta \vec{x}^{(i)} = \vec{x}^{(i)} - \textbf{J}^{-1} \vec{f} (\vec{x}^{(i)})
\end{align}$$
To sum up, the NR algorithm, for application to the power flow problem is:
- Set the iteration counter to $i=1$. Use the initial solution $V_{i} = 1 \angle 0^{\circ}$
- Compute the mismatch vector $\vec{f}({\vec{x}})$ using the power flow equations
- Perform the following stopping criterion tests:
- If $\vert \Delta P_{i} \vert < \epsilon_{P}$ for all type PQ and PV buses and
- If $\vert \Delta Q_{i} \vert < \epsilon_{Q}$ for all type PQ
- Then go to step 6
- Otherwise, go to step 4.
- Evaluate the Jacobian matrix $\textbf{J}^{(i)}$ and compute $\Delta \vec{x}^{(i)}$.
- Compute the update solution vector $\vec{x}^{(i+1)}$. Return to step 3.
- Stop.
Convergence and Step Control
The iteration is governed by two parameters:
- Tolerance (default $10^{-8}$): the run is converged once every entry of the
mismatch vector satisfies $\vert f_{i}(\vec{x}) \vert <$ tolerance (an infinity-norm
test over all $\Delta P$ and $\Delta Q$ components).
- Maximum iterations (default $20$): an upper bound on the number of Newton steps
per power flow solve.
To improve robustness far from the solution, the full Newton step is scaled by a
single factor $\alpha \in (0, 1]$ rather than damped component-wise. The factor is the
largest value that keeps the per-step changes within fixed bounds:
$$\alpha = \min \left( 1,\ \frac{\Delta\theta_{max}}{\max_k \vert \Delta\theta_k \vert},\ \frac{\Delta V_{max}}{\max_k \vert \Delta V_k / V_k \vert} \right)$$
with $\Delta\theta_{max} = 0.2\ \text{rad}$ and $\Delta V_{max} = 0.1\ \text{pu}$. Because
the whole step is scaled by one factor, the Newton search direction is preserved, so
$\alpha = 1$ near the solution and quadratic convergence is retained; $\alpha < 1$ only
bounds large early steps. Voltage magnitudes are updated multiplicatively
($V_k \leftarrow V_k (1 + \alpha, \Delta V_k / V_k)$), consistent with the relative
voltage increment used in the Jacobian.
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.