AC Power Flow

JuliaGrid uses standard network components and the Unified Branch Model for power flow analysis, allowing load profiles, generator capacities, voltage specifications, contingency analysis, and planning to be defined efficiently.

To begin, let us generate the PowerSystem type, as illustrated by the following example:

@labels(Integer)

@power(MW, MVAr, MVA)
@voltage(pu, deg, V)

system = powerSystem()

addBus!(system; label = 1, type = 3)
addBus!(system; label = 2, type = 1, active = 21.7, reactive = 12.7)
addBus!(system; label = 3, type = 1, active = 11.2, reactive = -3.0)
addBus!(system; label = 4, type = 2, conductance = 2.1, susceptance = 1.2)

addBranch!(system; from = 1, to = 2, resistance = 0.02, reactance = 0.06)
addBranch!(system; from = 1, to = 3, resistance = 0.05, reactance = 0.21)
addBranch!(system; from = 2, to = 3, resistance = 0.13, reactance = 0.26)
addBranch!(system; from = 3, to = 4, reactance = 0.17, susceptance = 0.2, conductance = 1e-4)

addGenerator!(system; bus = 1)
addGenerator!(system; bus = 3, active = 40.0, reactive = 42.4)

To review, we can conceptualize the bus/branch model as the graph denoted by $\mathcal{G} = (\mathcal{N}, \mathcal{E})$, where we have the set of buses $\mathcal{N} = \{1, \dots, n\}$, and the set of branches $\mathcal{E} \subseteq \mathcal{N} \times \mathcal{N}$ within the power system:

julia> 𝒩 = collect(keys(system.bus.label))4-element Vector{Int64}:
 1
 2
 3
 4
julia> ℰ = hcat([𝒩[system.branch.layout.from] 𝒩[system.branch.layout.to]])4×2 Matrix{Int64}:
 1  2
 1  3
 2  3
 3  4

Nodal Network Equations

As previously demonstrated in the section on the Nodal Network Equations, we observe the system of equations:

\[ \mathbf{\bar {I}} = \mathbf{Y} \mathbf{\bar {V}}.\]

The complex current injection at the bus $i \in \mathcal{N}$ is defined as:

\[ \bar{I}_{i} = \cfrac{S_{i}^*}{\bar{V}_{i}^*}.\]

Thus, for any given bus $i \in \mathcal{N}$, we can express it as:

\[ \cfrac{S_{i}^*}{\bar{V}_{i}^*} = \sum_{j = 1}^n Y_{ij} \bar {V}_j.\]

The complex power injection denoted by $S_i$ consists of both the active power $P_i$ and reactive power $Q_i$. This relationship can be represented as follows:

\[ \cfrac{P_i - \text{j}Q_i}{\bar{V}_{i}} = \sum_{j = 1}^n Y_{ij} \bar {V}_j.\]

Recognizing that $Y_{ij} = G_{ij} + \text{j}B_{ij}$, $\bar{V}_{i} = V_i \text{e}^{\text{j}\theta_{i}}$, $\bar{V}_{j} = V_j \text{e}^{\text{j}\theta_{j}}$, and by defining $\theta_{ij} = \theta_{i} - \theta_{j}$, we can break down the above equation into its real and imaginary parts, resulting in two equations that describe bus $i \in \mathcal{N}$ as follows:

\[ \begin{aligned} {P}_{i} &={V}_{i}\sum\limits_{j=1}^n (G_{ij}\cos\theta_{ij}+B_{ij}\sin\theta_{ij})V_j\\ {Q}_{i} &={V}_{i}\sum\limits_{j=1}^n (G_{ij}\sin\theta_{ij}-B_{ij}\cos\theta_{ij})V_j. \end{aligned}\]

As demonstrated by the above equations, the bus $i \in \mathcal{N}$ contains four unknown variables: the active power injection ${P}_{i}$, reactive power injection ${Q}_{i}$, bus voltage magnitude ${V}_{i}$, and bus voltage angle ${\theta}_{i}$.

To solve these equations, it is necessary to specify two known variables. Although any two variables can be selected mathematically, the choice is determined by the devices that are connected to a particular bus. The standard options are listed in the table below, and these options are used to define the bus types [9].

Consequently, JuliaGrid operates with sets $\mathcal{N}_{\text{pq}}$ and $\mathcal{N}_{\text{pv}}$ that contain demand and generator buses, respectively, and exactly one slack bus in the set $\mathcal{N}_{\text{sb}}$. The bus types are stored in the variable:

julia> system.bus.layout.type4-element Vector{Int8}:
 3
 1
 1
 2

It should be noted that JuliaGrid cannot handle systems with multiple slack buses. Additionally, when using functions such as newtonRaphson, fastNewtonRaphsonBX, fastNewtonRaphsonXB, and gaussSeidel, the bus type can be modified as discussed in the section on Bus Type Modification.

Furthermore, active power injections ${P}_{i}$ and reactive power injections ${Q}_{i}$ can be expressed as:

\[ \begin{aligned} P_{i} &= P_{\text{p}i} - P_{\text{d}i} \\ Q_{i} &= Q_{\text{p}i} - Q_{\text{d}i}, \end{aligned}\]

where ${P}_{\text{d}i}$ and ${Q}_{\text{d}i}$ denote the active and reactive power demanded at the bus $i \in \mathcal{N}$, while ${P}_{\text{p}i}$ and ${Q}_{\text{p}i}$ correspond to the active and reactive power produced by the generators at the bus $i \in \mathcal{N}$.

To provide a more comprehensive understanding, it is important to note that each bus $i \in \mathcal{N}$ has the capacity to host multiple generators. This scenario can be conceptualized by introducing the set $\mathcal{S}_i$, which encompasses all generators connected to bus $i \in \mathcal{N}$. With this perspective, we can calculate the values of ${P}_{\text{p}i}$ and ${Q}_{\text{p}i}$ as follows:

\[ \begin{aligned} P_{\text{p}i} &= \sum_{k \in \mathcal{S}_i} P_{\text{g}k}\\ Q_{\text{p}i} &= \sum_{k \in \mathcal{S}_i} Q_{\text{g}k}, \end{aligned}\]

where $P_{\text{g}k}$ and $Q_{\text{g}k}$ represent the active and reactive power outputs of the $k$-th generator within the set $\mathcal{S}_i$.

As a way to summarize, the power injection vectors, represented as $\mathbf{P} = [P_i]$ and $\mathbf{Q} = [Q_i]$ can be computed based on the following variables and expressions:

julia> 𝐏 = system.bus.supply.active - system.bus.demand.active4-element Vector{Float64}:
  0.0
 -0.217
  0.28800000000000003
  0.0
julia> 𝐐 = system.bus.supply.reactive - system.bus.demand.reactive4-element Vector{Float64}:
  0.0
 -0.127
  0.454
  0.0

Positive active or reactive power values, $P_i > 0$ or $Q_i > 0$, signify that power is being supplied into the power system from the specific bus. This indicates that the generators connected to this bus are producing more power than what the connected load is consuming. Conversely, negative values, $P_i < 0$ or $Q_i < 0$, indicate that the bus is drawing in active or reactive power from the power system. This suggests that the load's demand is exceeding the output from the generators.

Newton-Raphson Method

The Newton-Raphson method is commonly used in AC power flow calculations due to its quadratic rate of convergence. It provides an accurate approximation of the roots of the system of nonlinear equations:

\[ \mathbf{f}(\mathbf{x}) = \mathbf{0}.\]

This, in turn, allows for the determination of unknown voltage magnitudes and angles of buses, represented by the state vector $\mathbf x = [\mathbf x_\text{a}, \mathbf x_\text{m}]^T$. The state vector comprises two components:

• $\mathbf x_\text{a} \in \mathbb{R}^{n-1}$, which holds the bus voltage angles of demand and generator buses, represented by $\mathbf x_\text{a} = [\theta_i]$, where $i \in \mathcal{N}_{\text{pq}} \cup \mathcal{N}_{\text{pv}}$,
• $\mathbf x_\text{m} \in \mathbb{R}^{n_{\text{pq}}}$, which holds the bus voltage magnitudes of demand buses, represented by $\mathbf x_\text{m} = [V_i]$, where $i \in \mathcal{N}_{\text{pq}}$, and $n_{\text{pq}} = |\mathcal{N}_{\text{pq}}|$.

Knowing the voltage magnitudes and angles for certain types of buses is a consequence of the structure of the state vector $\mathbf x$. Specifically, the voltage magnitude and angle at the slack bus are known, as well as the voltage magnitude at generator buses.

As detailed in the Nodal Network Equations section of this manual, the expressions for active and reactive power injection are as follows:

\[ \begin{aligned} {P}_{i} &={V}_{i}\sum\limits_{j=1}^n (G_{ij}\cos\theta_{ij}+B_{ij}\sin\theta_{ij})V_j\\ {Q}_{i} &={V}_{i}\sum\limits_{j=1}^n (G_{ij}\sin\theta_{ij}-B_{ij}\cos\theta_{ij})V_j. \end{aligned}\]

Using the above equations, we can define the active power injection function for demand and generator buses:

\[ f_{P_i}(\mathbf x) = {V}_{i}\sum\limits_{j=1}^n (G_{ij}\cos\theta_{ij}+B_{ij}\sin\theta_{ij})V_j - {P}_{i} = 0, \;\;\; \forall i \in \mathcal{N}_{\text{pq}} \cup \mathcal{N}_{\text{pv}},\]

and the reactive power injection function for demand buses:

\[ f_{Q_i}(\mathbf x) = {V}_{i}\sum\limits_{j=1}^n (G_{ij}\sin\theta_{ij}-B_{ij}\cos\theta_{ij})V_j - {Q}_{i} = 0, \;\;\; \forall i \in \mathcal{N}_{\text{pq}}.\]

The active and reactive mismatches, often denoted as $\Delta P_i(\mathbf x)$ and $\Delta Q_i(\mathbf x)$, respectively, are defined as the functions $f_{P_i}(\mathbf x)$ and $f_{Q_i}(\mathbf x)$. The first terms on the right-hand side represent power injections at a bus, while the second term is constant and is obtained based on the active and reactive powers of the generators that supply a bus and active and reactive powers demanded by consumers at the same bus. Therefore, the Newton-Raphson method solves the system of nonlinear equations:

\[ \mathbf{f(x)} = \begin{bmatrix} \mathbf{f}_{\text{P}}(\mathbf x) \\ \mathbf{f}_{\text{Q}}(\mathbf x) \end{bmatrix} = \mathbf 0,\]

where the first $n - 1$ equations correspond to demand and generator buses, and the last $n_{\text{pq}}$ equations correspond to demand buses.

Initialization

To compute the voltage magnitudes and angles of buses using the Newton-Raphson method in JuliaGrid, you must first execute the acModel! function to set up the system, followed by initializing the Newton-Raphson method using the newtonRaphson function. The following code snippet demonstrates this process:

acModel!(system)
analysis = newtonRaphson(system)

This results in the creation of the starting vectors of bus voltage magnitudes $\mathbf{V}^{(0)}$ and angles $\bm{\Theta}^{(0)}$, as shown below:

julia> 𝐕⁽⁰⁾ = analysis.voltage.magnitude4-element Vector{Float64}:
 1.0
 1.0
 1.0
 1.0
julia> 𝚯⁽⁰⁾ = analysis.voltage.angle4-element Vector{Float64}:
 0.0
 0.0
 0.0
 0.0

Here, we utilize a "flat start" approach in our method. It is important to keep in mind that when dealing with initial conditions in this manner, the Newton-Raphson method may encounter difficulties.

Iterative Process

To implement the Newton-Raphson method, the iterative approach based on the Taylor series expansion, JuliaGrid provides the mismatch! and solve! functions. These functions are utilized to carry out the Newton-Raphson method iteratively until a stopping criterion is reached, as demonstrated in the following code snippet:

for iteration = 1:100
    stopping = mismatch!(system, analysis)
    if all(stopping .< 1e-8)
        break
    end
    solve!(system, analysis)
end

The mismatch! function calculates the mismatch in active power injection for demand and generator buses and the mismatch in reactive power injection for demand buses at each iteration $\nu = \{1, 2, \dots\}$. The equations used for these computations are:

\[ f_{P_i}(\mathbf x^{(\nu-1)}) = {V}_{i}^{(\nu-1)}\sum\limits_{j=1}^n (G_{ij}\cos\theta_{ij}^{(\nu-1)}+B_{ij}\sin\theta_{ij}^{(\nu-1)}){V}_{j}^{(\nu-1)} - {P}_{i}, \;\;\; \forall i \in \mathcal{N}_{\text{pq}} \cup \mathcal{N}_{\text{pv}},\]

as well as the reactive power injection mismatch for demand buses:

\[ f_{Q_i}(\mathbf x^{(\nu-1)}) = {V}_{i}^{(\nu)}\sum\limits_{j=1}^n (G_{ij}\sin\theta_{ij}^{(\nu-1)}-B_{ij}\cos\theta_{ij}^{(\nu-1)}){V}_{j}^{(\nu-1)} - {Q}_{i}, \;\;\; \forall i \in \mathcal{N}_{\text{pq}}.\]

The resulting vector from these calculations is stored in the mismatch variable of the ACPowerFlow abstract type and can be accessed through the following line of code:

julia> 𝐟 = analysis.method.mismatch6-element Vector{Float64}:
  1.5543122344752192e-15
  4.579669976578771e-16
 -1.810613397539215e-16
  3.497202527569243e-15
  1.1345091532888318e-15
  0.0

In addition to computing the mismatches in active and reactive power injection, the mismatch! function also returns the maximum absolute values of these mismatches. These maximum values are used as termination criteria for the iteration loop if both are less than a predefined stopping criterion $\epsilon$:

\[ \max \{|f_{P_i}(\mathbf x^{(\nu-1)})|,\; \forall i \in \mathcal{N}_{\text{pq}} \cup \mathcal{N}_{\text{pv}} \} < \epsilon \\ \max \{|f_{Q_i}(\mathbf x^{(\nu-1)})|,\; \forall i \in \mathcal{N}_{\text{pq}} \} < \epsilon.\]

Next, the function solve! computes the increments of bus voltage angle and magnitude at each iteration using:

\[ \mathbf{\Delta} \mathbf{x}^{(\nu-1)} = -\mathbf{J}(\mathbf{x}^{(\nu-1)})^{-1} \mathbf{f}(\mathbf{x}^{(\nu-1)}),\]

where $\mathbf{\Delta} \mathbf{x} = [\mathbf \Delta \mathbf x_\text{a}, \mathbf \Delta \mathbf x_\text{m}]^T$ consists of the vector of bus voltage angle increments $\mathbf \Delta \mathbf x_\text{a} \in \mathbb{R}^{n-1}$ and bus voltage magnitude increments $\mathbf \Delta \mathbf x_\text{m} \in \mathbb{R}^{n_{\text{pq}}}$, and $\mathbf{J}(\mathbf{x}) \in \mathbb{R}^{n_{\text{u}} \times n_{\text{u}}}$ is the Jacobian matrix, $n_{\text{u}} = n + n_{\text{pq}} - 1$.

These values are stored in the ACPowerFlow abstract type and can be accessed after each iteration:

julia> 𝚫𝐱 = analysis.method.increment6-element Vector{Float64}:
  7.92699879268019e-11
 -3.5508852576600275e-10
 -2.6233074858150384e-9
  9.233062009151665e-10
  5.1042499891742e-9
  7.485907667424543e-9
julia> 𝐉 = analysis.method.jacobian6×6 SparseArrays.SparseMatrixCSC{Float64, Int64} with 28 stored entries:
 18.4159   -3.36093     ⋅         6.36094  -1.56596      ⋅
 -3.38112  15.4035    -7.11439   -1.65565   3.10891     0.0233733
   ⋅       -7.11439    7.11439     ⋅       -0.0238268   0.0233733
 -6.83212   1.7057      ⋅        18.0563   -3.08559      ⋅
  1.66532  -2.81034   -0.025953  -3.36147  14.9752     -6.40723
   ⋅        0.025953  -0.025953    ⋅       -6.53154     6.40723
julia> 𝐋 = analysis.method.factorization.L6×6 SparseArrays.SparseMatrixCSC{Float64, Int64} with 17 stored entries:
  1.0           ⋅            ⋅           ⋅           ⋅         ⋅
 -0.00400889   1.0           ⋅           ⋅           ⋅         ⋅
   ⋅            ⋅           1.0          ⋅           ⋅         ⋅
 -0.465213     0.00308853  -0.177715    1.0          ⋅         ⋅
   ⋅            ⋅          -0.371291    0.0618339   1.0        ⋅
 -0.00178072  -0.444182     0.0918455  -0.346344   -0.203397  1.0
julia> 𝐔 = analysis.method.factorization.U6×6 SparseArrays.SparseMatrixCSC{Float64, Int64} with 17 stored entries:
 0.498347  0.00163724   ⋅        -0.498347       ⋅         -0.00166901
  ⋅        0.493224     ⋅         4.33681e-19    ⋅         -0.502794
  ⋅         ⋅          0.619986  -0.113149      0.214146   -0.0527194
  ⋅         ⋅           ⋅         0.250011     -0.0158957   0.092718
  ⋅         ⋅           ⋅          ⋅            0.688866   -0.12927
  ⋅         ⋅           ⋅          ⋅             ⋅          0.299378

The JuliaGrid implementation of the AC power flow follows a specific order to store the increment $\mathbf{\Delta} \mathbf{x}$ and mismatch $\mathbf{f(x)}$ vectors. The first $n-1$ elements of both vectors correspond to the demand and generator buses in the same order as they appear in the input data. The first $n-1$ elements of the increment vector $\mathbf{\Delta} \mathbf{x}$ correspond to the voltage angle increments $\mathbf \Delta \mathbf x_\text{a}$, while the first $n-1$ elements of the mismatch vector $\mathbf{f(x)}$ correspond to the mismatch in active power injections $\mathbf{f}_{\text{P}}(\mathbf x)$.

The last $n_{\text{pq}}$ elements of the increment $\mathbf{\Delta} \mathbf{x}$ and mismatch $\mathbf{f(x)}$ vectors correspond to the demand buses in the order they appear in the input data. For the increment vector $\mathbf{\Delta} \mathbf{x}$, it matches the bus voltage magnitude increments $\mathbf \Delta \mathbf x_\text{m}$, while for the mismatch vector $\mathbf{f(x)}$, it matches the mismatch in reactive power injections $\mathbf{f}_{\text{Q}}(\mathbf x)$.

These specified orders dictate the row and column order of the Jacobian matrix $\mathbf{J}(\mathbf{x})$.

Finally, the function solve! adds the computed increment term to the previous solution to obtain a new solution:

\[ \mathbf {x}^{(\nu)} = \mathbf {x}^{(\nu-1)} + \mathbf \Delta \mathbf {x}^{(\nu-1)}.\]

The bus voltage magnitudes $\mathbf{V} = [V_i]$ and angles $\bm{\Theta} = [\theta_i]$ are then updated based on the obtained solution $\mathbf {x}$. It is important to note that only the voltage magnitudes related to demand buses and angles related to demand and generator buses are updated; not all values are updated. Therefore, the final solution obtained by JuliaGrid is stored in the following vectors:

julia> 𝐕 = analysis.voltage.magnitude4-element Vector{Float64}:
 1.0
 1.0058448714519173
 1.0892355535521518
 1.1103697460384185
julia> 𝚯 = analysis.voltage.angle4-element Vector{Float64}:
  0.0
 -0.00644900951642222
 -0.00046072472071607
 -0.0041086656422506945

Jacobian Matrix

To complete the tutorial on the Newton-Raphson method, we will now describe the Jacobian matrix and provide the equations involved in its evolution. Without loss of generality, we assume that the slack bus is the first bus, followed by the set of demand buses and the set of generator buses:

\[ \begin{aligned} \mathcal{N}_{\text{sb}} &= \{ 1 \} \\ \mathcal{N}_{\text{pq}} &= \{2, \dots, m\} \\ \mathcal{N}_{\text{pv}} &= \{m + 1,\dots, n\}, \end{aligned}\]

where $\mathcal{N} = \mathcal{N}_{\text{sb}} \cup \mathcal{N}_{\text{pq}} \cup \mathcal{N}_{\text{pv}}$. Therefore, we can express:

\[ \begin{aligned} \mathbf x_\text{a} &= [\theta_2,\dots,\theta_n]^T, \;\;\;\;\;\; \mathbf \Delta \mathbf x_\text{a} = [\Delta \theta_2,\dots,\Delta \theta_n]^T \\ \mathbf x_\text{m} &= [V_2,\dots,V_{m}]^T, \;\;\; \mathbf \Delta \mathbf x_\text{m} = [\Delta V_2,\dots,\Delta V_{m}]^T. \end{aligned}\]

The Jacobian matrix $\mathbf{J(x^{(\nu)})} \in \mathbb{R}^{n_{\text{u}} \times n_{\text{u}}}$ is:

\[ \mathbf{J(x^{(\nu)})}= \left[ \begin{array}{ccc|ccc} \cfrac{\mathrm \partial{{f_{P_2}}(\mathbf x^{(\nu)})}} {\mathrm \partial \theta_2} & \cdots & \cfrac{\mathrm \partial{{f_{P_2}}(\mathbf x^{(\nu)})}}{\mathrm \partial \theta_{n}} & \cfrac{\mathrm \partial{{f_{P_2}}(\mathbf x^{(\nu)})}}{\mathrm \partial V_{2}} &\cdots & \cfrac{\mathrm \partial{{f_{P_2}}(\mathbf x^{(\nu)})}}{\mathrm \partial V_{m}}\\ \;\vdots & \\ \cfrac{\mathrm \partial{{f_{P_n}}(\mathbf x^{(\nu)})}} {\mathrm \partial \theta_2} & \cdots & \cfrac{\mathrm \partial{{f_{P_n}}(\mathbf x^{(\nu)})}}{\mathrm \partial \theta_{n}} & \cfrac{\mathrm \partial{{f_{P_n}}(\mathbf x^{(\nu)})}}{\mathrm \partial V_{2}} &\cdots & \cfrac{\mathrm \partial{{f_{P_n}}(\mathbf x^{(\nu)})}}{\mathrm \partial V_{m}} \\[10pt] \hline \\ \cfrac{\mathrm \partial{{f_{Q_2}}(\mathbf x^{(\nu)})}} {\mathrm \partial \theta_2} & \cdots & \cfrac{\mathrm \partial{{f_{Q_2}}(\mathbf x^{(\nu)})}}{\mathrm \partial \theta_{n}} & \cfrac{\mathrm \partial{{f_{Q_2}}(\mathbf x^{(\nu)})}}{\mathrm \partial V_{2}} &\cdots & \cfrac{\mathrm \partial{{f_{Q_2}}(\mathbf x^{(\nu)})}}{\mathrm \partial V_{m}}\\ \;\vdots & \\ \cfrac{\mathrm \partial{{f_{Q_{m}}}(\mathbf x^{(\nu)})}} {\mathrm \partial \theta_2} & \cdots & \cfrac{\mathrm \partial{{f_{Q_{m}}}(\mathbf x^{(\nu)})}}{\mathrm \partial \theta_{n}} & \cfrac{\mathrm \partial{{f_{Q_{m}}}(\mathbf x^{(\nu)})}}{\mathrm \partial V_{2}} &\cdots & \cfrac{\mathrm \partial{{f_{Q_{m}}}(\mathbf x^{(\nu)})}}{\mathrm \partial V_{m}} \end{array} \right].\]

The Jacobian matrix can be expressed using four block matrices:

\[ \mathbf{J(x^{(\nu)})} = \begin{bmatrix} \mathbf{J_{11}(x^{(\nu)})} &\mathbf{J_{12}(x^{(\nu)})} \\ \mathbf{J_{21}(x^{(\nu)})} & \mathbf{J_{22}(x^{(\nu)})} \end{bmatrix},\]

where diagonal elements of the Jacobian sub-matrices are defined as follows:

\[ \begin{aligned} \cfrac{\mathrm \partial{{f_{P_i}}(\mathbf x^{(\nu)})}} {\mathrm \partial \theta_{i}} &= {V}_{i}^{(\nu)}\sum\limits_{j=1}^n (-G_{ij}\sin\theta_{ij}^{(\nu)}+B_{ij}\cos\theta_{ij}^{(\nu)}){V}_{j}^{(\nu)} - B_{ii}({V}_{i}^{(\nu)})^2\\ \cfrac{\mathrm \partial{{f_{P_i}}(\mathbf x^{(\nu)})}}{\mathrm \partial V_{i}^{(\nu)}} &= \sum\limits_{j=1}^n (G_{ij}\cos\theta_{ij}^{(\nu)}+B_{ij}\sin\theta_{ij}^{(\nu)}){V}_{j}^{(\nu)} + G_{ii}{V}_{i}^{(\nu)}\\ \cfrac{\mathrm \partial{{f_{Q_i}}(\mathbf x^{(\nu)})}} {\mathrm \partial \theta_{i}} &= {V}_{i}^{(\nu)} \sum\limits_{j=1}^n (G_{ij}\cos\theta_{ij}^{(\nu)}+B_{ij}\sin\theta_{ij}^{(\nu)}){V}_{j}^{(\nu)} - G_{ii}({V}_{i}^{(\nu)})^2\\ \cfrac{\mathrm \partial{{f_{Q_i}}(\mathbf x^{(\nu)})}}{\mathrm \partial V_{i}} &= \sum\limits_{j=1}^n (G_{ij}\sin\theta_{ij}^{(\nu)}-B_{ij}\cos\theta_{ij}^{(\nu)}){V}_{j}^{(\nu)} - B_{ii}{V}_{i}^{(\nu)}, \end{aligned}\]

while non-diagonal elements of the Jacobian sub-matrices are:

\[ \begin{aligned} \cfrac{\mathrm \partial{{f_{P_i}}(\mathbf x^{(\nu)})}}{\mathrm \partial \theta_{j}} &= (G_{ij}\sin\theta_{ij}^{(\nu)}-B_{ij}\cos\theta_{ij}^{(\nu)}){V}_{i}^{(\nu)}{V}_{j}^{(\nu)}\\ \cfrac{\mathrm \partial{{f_{P_i}}(\mathbf x^{(\nu)})}}{\mathrm \partial V_{j}^{(\nu)}} &= (G_{ij}\cos\theta_{ij}^{(\nu)}+B_{ij}\sin\theta_{ij}^{(\nu)}){V}_{i}^{(\nu)}\\ \cfrac{\mathrm \partial{{f_{Q_i}}(\mathbf x^{(\nu)})}}{\mathrm \partial \theta_{j}} &= -(G_{ij}\cos\theta_{ij}^{(\nu)} + B_{ij}\sin\theta_{ij}^{(\nu)}){V}_{i}^{(\nu)}{V}_{j}^{(\nu)}\\ \cfrac{\mathrm \partial{{f_{Q_i}}(\mathbf x^{(\nu)})}}{\mathrm\partial V_{j}} &= (G_{ij}\sin\theta_{ij}^{(\nu)}-B_{ij}\cos\theta_{ij}^{(\nu)}){V}_{i}^{(\nu)}. \end{aligned}\]

Fast Newton-Raphson Method

Although the fast Newton-Raphson method may converge more slowly than the traditional Newton-Raphson method, the shorter solution time for the updates often compensates for this slower convergence, resulting in a shorter overall solution time. This is particularly true for systems that are not heavily loaded, where a shorter overall solution time is almost always achieved. It is important to note that if the algorithm converges, it will converge to a correct solution [6].

The fast Newton-Raphson method involves decoupling the power flow equations. Namely, the Newton-Raphson method is based on the equations:

\[ \begin{bmatrix} \mathbf{J_{11}(x)} &\mathbf{J_{12}(x)} \\ \mathbf{J_{21}(x)} & \mathbf{J_{22}(x)} \end{bmatrix} \begin{bmatrix} \mathbf \Delta \mathbf x_\text{a} \\ \mathbf \Delta \mathbf x_\text{m} \end{bmatrix} + \begin{bmatrix} \mathbf{f}_{\text{P}}(\mathbf x) \\ \mathbf{f}_{\text{Q}}(\mathbf x) \end{bmatrix} = \mathbf 0,\]

where the iteration index has been omitted for simplicity. However, in transmission grids, there exists a strong coupling between active powers and voltage angles, as well as between reactive powers and voltage magnitudes. To achieve decoupling, two conditions must be satisfied: first, the resistance values $r_{ij}$ of the branches must be small compared to their reactance values $x_{ij}$, and second, the angle differences must be small, i.e., $\theta_{ij} \approx 0$ [10]. Therefore, starting from the above equation, we have:

\[ \begin{bmatrix} \mathbf{J_{11}(x)} & \mathbf{0} \\ \mathbf{0} & \mathbf{J_{22}(x)} \end{bmatrix} \begin{bmatrix} \mathbf \Delta \mathbf x_\text{a} \\ \mathbf \Delta \mathbf x_\text{m} \end{bmatrix} + \begin{bmatrix} \mathbf{f}_{\text{P}}(\mathbf x) \\ \mathbf{f}_{\text{Q}}(\mathbf x) \end{bmatrix} = \mathbf 0,\]

which gives the decoupled system as follows:

\[ \begin{aligned} \mathbf{f}_{\text{P}}(\mathbf x) &= -\mathbf{J_{11}(x)} \mathbf \Delta \mathbf x_\text{a} \\ \mathbf{f}_{\text{Q}}(\mathbf x) &= -\mathbf{J_{22}(x)} \mathbf \Delta \mathbf x_\text{m}. \end{aligned}\]

To examine the problem, it is helpful to express it as:

\[ \begin{aligned} {f}_{P_2}(\mathbf x) &= -\Delta \theta_2\cfrac{\mathrm \partial{{f_{P_2}}(\mathbf x)}} {\mathrm \partial \theta_2} - \cdots - \Delta \theta_n \cfrac{\mathrm \partial{{f_{P_2}}(\mathbf x)}}{\mathrm \partial \theta_{n}} \\ & \vdots \\ {f}_{P_n}(\mathbf x) &= -\Delta \theta_2\cfrac{\mathrm \partial{{f_{P_n}}(\mathbf x)}} {\mathrm \partial \theta_2} - \cdots - \Delta \theta_n \cfrac{\mathrm \partial{{f_{P_i}}(\mathbf x)}}{\mathrm \partial \theta_{n}}\\ {f}_{Q_2}(\mathbf x) &= - \Delta V_2 \cfrac{\mathrm \partial{{f_{Q_2}}(\mathbf x)}}{\mathrm \partial V_{2}} - \cdots - \Delta V_{n_{\text{pq}}} \cfrac{\mathrm \partial{{f_{Q_2}}(\mathbf x)}}{\mathrm \partial V_{m}}\\ & \vdots \\ {f}_{Q_{m}}(\mathbf x) &= - \Delta V_2 \cfrac{\mathrm \partial{{f_{Q_{m}}}(\mathbf x)}}{\mathrm \partial V_{2}} - \cdots - \Delta V_{m} \cfrac{\mathrm \partial{{f_{Q_{m}}}(\mathbf x)}}{\mathrm \partial V_{m}}. \end{aligned}\]

Firstly, the second part of the expressions is expanded as follows:

\[ \begin{aligned} {f}_{Q_2}(\mathbf x) &= -\cfrac{\Delta V_2}{V_2}V_2 \cfrac{\mathrm \partial{{f_{Q_2}}(\mathbf x)}}{\mathrm \partial V_{2}} - \cdots - \cfrac{\Delta V_{m}}{V_{m}} V_{m} \cfrac{\mathrm \partial{{f_{Q_2}}(\mathbf x)}}{\mathrm \partial V_{m}}\\ & \vdots \\ {f}_{Q_{m}}(\mathbf x) &= - \cfrac{\Delta V_2}{V_2}V_2 \cfrac{\mathrm \partial{{f_{Q_{m}}}(\mathbf x)}}{\mathrm \partial V_{2}} - \cdots - \cfrac{\Delta V_{m}}{V_{m}} V_{m} \cfrac{\mathrm \partial{{f_{Q_{m}}}(\mathbf x)}}{\mathrm \partial V_{m}}. \end{aligned}\]

Next, the Jacobian elements are derived. To achieve this, we can use the expressions defined for the Newton-Raphson method. For demand buses, the above expansions are applied as:

\[ \begin{aligned} \cfrac{\mathrm \partial{{f_{P_i}}(\mathbf x)}} {\mathrm \partial \theta_{i}} &= {V}_{i}\sum\limits_{j=1}^n (-G_{ij}\sin\theta_{ij}+B_{ij}\cos\theta_{ij}){V}_{j} - B_{ii}{V}_{i}^2\\ \cfrac{\mathrm \partial{{f_{P_i}}(\mathbf x)}}{\mathrm \partial \theta_{j}} &= (G_{ij}\sin\theta_{ij}-B_{ij}\cos\theta_{ij}){V}_{i}{V}_{j}\\ V_i \cfrac{\mathrm \partial{{f_{Q_i}}(\mathbf x)}}{\mathrm \partial V_{i}} &= V_i\sum\limits_{j=1}^n (G_{ij}\sin\theta_{ij}-B_{ij}\cos\theta_{ij}){V}_{j} - B_{ii}{V}_{i}^2\\ V_j \cfrac{\mathrm \partial{{f_{Q_i}}(\mathbf x)}}{\mathrm\partial V_{j}} &= (G_{ij}\sin\theta_{ij}-B_{ij}\cos\theta_{ij}) V_i V_j. \end{aligned}\]

As the definition of reactive power is given by the equation:

\[ {Q}_{i} ={V}_{i}\sum\limits_{j=1}^n (G_{ij}\sin\theta_{ij}-B_{ij}\cos\theta_{ij})V_j,\]

the Jacobian elements can be expressed in the following manner:

\[ \begin{aligned} \cfrac{\mathrm \partial{{f_{P_i}}(\mathbf x)}} {\mathrm \partial \theta_{i}} &= -Q_i - B_{ii}{V}_{i}^2\\ \cfrac{\mathrm \partial{{f_{P_i}}(\mathbf x)}}{\mathrm \partial \theta_{j}} &= (G_{ij}\sin\theta_{ij}-B_{ij}\cos\theta_{ij}) V_i V_j\\ V_i \cfrac{\mathrm \partial{{f_{Q_i}}(\mathbf x)}}{\mathrm \partial V_{i}} &= Q_i - B_{ii} V_i^2\\ V_j \cfrac{\mathrm \partial{{f_{Q_i}}(\mathbf x)}}{\mathrm\partial V_{j}} &= (G_{ij}\sin\theta_{ij} - B_{ij}\cos\theta_{ij}) V_i V_j. \end{aligned}\]

The decoupled model is established through the following approximations:

\[ \begin{aligned} \sin(\theta_{ij}) \approx 0 \\ \cos(\theta_{ij}) \approx 1 \\ Q_i << B_{ii}V_i^2. \end{aligned}\]

Thus, when the approximations are made, the Jacobian elements are simplified, resulting in the decoupled model where the Jacobian elements are:

\[ \begin{aligned} \cfrac{\mathrm \partial{{f_{P_i}}(\mathbf x)}} {\mathrm \partial \theta_{i}} &= -B_{ii}{V}_{i}^2\\ \cfrac{\mathrm \partial{{f_{P_i}}(\mathbf x)}} {\mathrm \partial \theta_{j}} &= -B_{ij}{V}_{i}{V}_{j}\\ V_i \cfrac{\mathrm \partial{{f_{Q_i}}(\mathbf x)}} {\mathrm \partial V_{i}} &= -B_{ii}{V}_{i}^2\\ V_j \cfrac{\mathrm \partial{{f_{Q_i}}(\mathbf x)}}{\mathrm\partial V_{j}} &= -B_{ij}{V}_{i}{V}_{j}. \end{aligned}\]

Thus, the initial system of equations becomes:

\[ \begin{aligned} {f}_{P_2}(\mathbf x) &= B_{22} \Delta \theta_2 {V}_{2}^2 + \cdots + B_{2n} \Delta \theta_n {V}_{2}{V}_{n} \\ & \vdots \\ {f}_{P_n}(\mathbf x) &= B_{n2} \Delta \theta_2 {V}_{2}{V}_{n} + \cdots + B_{nn} \Delta \theta_n {V}_{n}^2 \\ {f}_{Q_2}(\mathbf x) &= B_{22} \cfrac{\Delta V_2}{V_2} {V}_{2}^2 + \cdots + B_{2m} \cfrac{\Delta V_{m}}{V_{m}} {V}_{2}V_{m} \\ & \vdots \\ {f}_{Q_{m}}(\mathbf x) &= B_{m2} \cfrac{\Delta V_2}{V_2} {V}_{2}V_{m} + \cdots + B_{mm} \cfrac{\Delta V_{m}}{V_{m}} V_{m}^2. \end{aligned}\]

Using $V_j \approx 1$, wherein $V_i^2 = V_iV_j, j=i$, the first part of the equations can be simplified to:

\[ \begin{aligned} {f}_{P_2}(\mathbf x) &= B_{22} \Delta \theta_2 {V}_{2} + \cdots + B_{2n} \Delta \theta_n {V}_{2}\\ & \vdots \\ {f}_{P_n}(\mathbf x) &= B_{n2} \Delta \theta_2 {V}_{n} + \cdots + B_{nn} \Delta \theta_n {V}_{n}. \end{aligned}\]

Similarly, the second part of the equations can be simplified to:

\[ \begin{aligned} {f}_{Q_2}(\mathbf x) &= B_{22} {V}_{2} \Delta V_2 + \cdots + B_{2m} V_2 \Delta V_{m} \\ & \vdots \\ {f}_{Q_{m}}(\mathbf x) &= B_{m2} V_{m} \Delta V_2 + \cdots + B_{mm} V_{m} \Delta V_{m}. \end{aligned}\]

The fast Newton-Raphson method is ultimately based on the system of equations presented below:

\[ \begin{aligned} \cfrac{{f}_{P_2}(\mathbf x)}{{V}_{2}} &= B_{22} \Delta \theta_2 + \cdots + B_{2n} \Delta \theta_n \\ & \vdots \\ \cfrac{{f}_{P_n}(\mathbf x)}{{V}_{n}} &= B_{n2} \Delta \theta_2 + \cdots + B_{nn} \Delta \theta_n \\ \cfrac{{f}_{Q_2}(\mathbf x)}{{V}_{2}} &= B_{22} \Delta V_2 + \cdots + B_{2m} \Delta V_{m} \\ & \vdots \\ \cfrac{{f}_{Q_{m}}(\mathbf x)}{V_{m}} &= B_{m2} \Delta V_2 + \cdots + B_{mm} \Delta V_{m}. \end{aligned}\]

This system can be written as:

\[ \begin{aligned} \mathbf{h}_{\text{P}}(\mathbf x) &= \mathbf{B}_1 \mathbf \Delta \mathbf x_\text{a} \\ \mathbf{h}_{\text{Q}}(\mathbf x) &= \mathbf{B}_2 \mathbf \Delta \mathbf x_\text{m}. \end{aligned}\]

One of the main advantages of this approach is that the Jacobian matrices $\mathbf{B}_1$ and $\mathbf{B}_2$ are constant and need only be formed once. Furthermore, this method can be used to define both the XB and BX versions of the fast Newton-Raphson method.

XB Version

The matrix $\mathbf{B}_1$ is formed by neglecting the resistance $r_{ij}$, shunt susceptance $\Im \{ y_{\text{sh}i} \}$, charging susceptance $\Im \{ y_{\text{s}ij} \}$, and transformer tap ratio magnitude $\tau_{ij}$. The matrix $\mathbf{B}_2$ is constructed by disregarding the transformer phase shift angle $\phi_{ij}$. This approach corresponds to the standard fast Newton-Raphson method and is known to exhibit exceptional convergence properties in typical scenarios [10].

To initialize the XB version of the fast Newton-Raphson method, one can utilize the following code snippet:

acModel!(system)
analysis = fastNewtonRaphsonXB(system)

BX Version

The matrix $\mathbf{B}_1$ ignores the shunt susceptance$\Im \{ y_{\text{sh}i} \}$, charging susceptance $\Im \{ y_{\text{s}ij} \}$, and transformer tap ratio magnitude $\tau_{ij}$. The matrix $\mathbf{B}_2$ ignores the resistance $r_{ij}$ and transformer phase shift angle $\phi_{ij}$. In usual cases, the iteration count for the BX version is comparable to the XB scheme. However, for systems with high $r_{ij}/x_{ij}$ ratios, the BX scheme requires considerably fewer iterations than the XB scheme to solve the power flow [10].

To initialize the BX version of the fast Newton-Raphson method, you can use the following code:

acModel!(system)
analysis = fastNewtonRaphsonBX(system)

Initialization

When a user creates the fast Newton-Raphson method in JuliaGrid, the Jacobian matrices $\mathbf{B}_1$ and $\mathbf{B}_2$ are formed to correspond to the active and reactive power equations, respectively:

julia> 𝐁₁ = analysis.method.active.jacobian3×3 SparseArrays.SparseMatrixCSC{Float64, Int64} with 7 stored entries:
 -18.0769     3.07692    ⋅
   3.07692  -13.4657    5.88235
    ⋅         5.88235  -5.88235
julia> 𝐁₂ = analysis.method.reactive.jacobian3×3 SparseArrays.SparseMatrixCSC{Float64, Int64} with 7 stored entries:
 -20.5128     3.84615    ⋅
   3.84615  -14.3904    5.88235
    ⋅         5.88235  -5.77035

Additionally, during this stage, JuliaGrid generates the starting vectors for bus voltage magnitudes $\mathbf{V}^{(0)}$ and angles $\bm{\Theta}^{(0)}$ as demonstrated below:

julia> 𝐕⁽⁰⁾ = analysis.voltage.magnitude4-element Vector{Float64}:
 1.0
 1.0
 1.0
 1.0
julia> 𝚯⁽⁰⁾ = analysis.voltage.angle4-element Vector{Float64}:
 0.0
 0.0
 0.0
 0.0

Iterative Process

JuliaGrid offers the mismatch! and solve! functions to implement the fast Newton-Raphson method iterations. These functions are used iteratively until a stopping criterion is met, as shown in the code snippet below:

for iteration = 1:100
    stopping = mismatch!(system, analysis)
    if all(stopping .< 1e-8)
        break
    end
    solve!(system, analysis)
end

The functions $\mathbf{f}_{\text{P}}(\mathbf x)$ and $\mathbf{f}_{\text{Q}}(\mathbf x)$ remain free of approximations, with only the calculation of the state variable increments affected [6]. As a result, we still use the following equations to compute the mismatches:

\[ \begin{aligned} f_{P_i}(\mathbf x) &= {V}_{i}\sum\limits_{j=1}^n (G_{ij}\cos\theta_{ij}+B_{ij}\sin\theta_{ij})V_j - {P}_{i} = 0, \;\;\; \forall i \in \mathcal{N}_{\text{pq}} \cup \mathcal{N}_{\text{pv}}\\ f_{Q_i}(\mathbf x) &= {V}_{i}\sum\limits_{j=1}^n (G_{ij}\sin\theta_{ij}-B_{ij}\cos\theta_{ij})V_j - {Q}_{i} = 0, \;\;\; \forall i \in \mathcal{N}_{\text{pq}}. \end{aligned}\]

Therefore, the mismatch! function calculates the mismatch in active power injection for demand and generator buses and the mismatch in reactive power injection for demand buses at each iteration $\nu = \{1, 2, \dots\}$:

\[ \begin{aligned} h_{P_i}(\mathbf {x}^{(\nu-1)}) &= \sum\limits_{j=1}^n (G_{ij}\cos\theta_{ij}^{(\nu-1)}+B_{ij}\sin\theta_{ij}^{(\nu-1)}){V}_{j}^{(\nu-1)} - \cfrac{{P}_{i}}{{V}_{i}^{(\nu-1)}}, \;\;\; \forall i \in \mathcal{N}_{\text{pq}} \cup \mathcal{N}_{\text{pv}} \\ h_{Q_i}(\mathbf {x}^{(\nu-1)}) &= \sum\limits_{j=1}^n (G_{ij}\sin\theta_{ij}^{(\nu-1)}-B_{ij}\cos\theta_{ij}^{(\nu-1)}){V}_{j}^{(\nu-1)} - \cfrac{{Q}_{i}}{{V}_{i}^{(\nu-1)}}, \;\;\; \forall i \in \mathcal{N}_{\text{pq}}. \end{aligned}\]

The resulting vectors from these calculations are stored in the ACPowerFlow abstract type and can be accessed through the following:

julia> 𝐡ₚ = analysis.method.active.increment3-element Vector{Float64}:
  1.0521034860626955e-9
 -1.8809173701904967e-10
 -2.280102039884849e-10
julia> 𝐡ₒ = analysis.method.reactive.increment3-element Vector{Float64}:
 9.293374091829092e-11
 4.4771575718049794e-12
 4.400285291193953e-12

In addition to computing the mismatches in active and reactive power injection, the mismatch! function also returns the maximum absolute values of these mismatches. These maximum values are used as termination criteria for the iteration loop if both are less than a predefined stopping criterion $\epsilon$:

\[ \max \{|h_{P_i}(\mathbf x^{(\nu)})|,\; \forall i \in \mathcal{N}_{\text{pq}} \cup \mathcal{N}_{\text{pv}} \} < \epsilon \\ \max \{|h_{Q_i}(\mathbf x^{(\nu)})|,\; \forall i \in \mathcal{N}_{\text{pq}} \} < \epsilon.\]

Next, the function solve! computes the bus voltage angle increments:

\[ \mathbf \Delta \mathbf x_\text{a}^{(\nu-1)} = \mathbf{B}_1^{-1} \mathbf{h}_{\text{P}}(\mathbf x^{(\nu-1)}).\]

To obtain the voltage angle increments, JuliaGrid initially performs LU factorization on the Jacobian matrix $\mathbf{B}_1 = \mathbf{L}_1\mathbf{U}_1$. This factorization is executed only once and is utilized in each iteration of the algorithm:

julia> 𝐋₁ = analysis.method.active.factorization.L3×3 SparseArrays.SparseMatrixCSC{Float64, Int64} with 5 stored entries:
  1.0         ⋅         ⋅
   ⋅         1.0        ⋅
 -0.524625  -0.160564  1.0
julia> 𝐔₁ = analysis.method.active.factorization.U3×3 SparseArrays.SparseMatrixCSC{Float64, Int64} with 5 stored entries:
 -0.5    ⋅         0.5
   ⋅   -0.854545   0.145455
   ⋅     ⋅        -0.314811

The vector of increments that corresponds to the active power equations can be accessed using:

julia> 𝚫𝐱ₐ = analysis.method.active.increment3-element Vector{Float64}:
  1.0521034860626955e-9
 -1.8809173701904967e-10
 -2.280102039884849e-10

The solution is then updated as follows:

\[ \mathbf x_\text{a}^{(\nu)} = \mathbf x_\text{a}^{(\nu-1)} + \mathbf \Delta \mathbf x_\text{a}^{(\nu-1)}.\]

It is important to note that only the voltage angles related to demand and generator buses are updated, while the vector of bus voltage angles of all buses is stored:

julia> 𝚯 = analysis.voltage.angle4-element Vector{Float64}:
  0.0
 -0.006449009470429831
 -0.00046072472210989724
 -0.004108665646915064

After calculating the update for voltage angles, to calculate the magnitude updates the fast Newton-Raphson method then solves the equation:

\[ \mathbf \Delta \mathbf x_\text{m}^{(\nu-1)} = \mathbf{B}_2^{-1} \mathbf{h}_{\text{Q}}(\mathbf x^{(\nu)}).\]

Similarly to the previous instance, JuliaGrid initially executes LU factorization on the Jacobian matrix $\mathbf{B}_2 = \mathbf{L}_2\mathbf{U}_2$. However, it provides the flexibility for users to opt for QR factorization instead. This factorization occurs only once and is utilized in each iteration of the fast Newton-Raphson algorithm:

julia> 𝐋₂ = analysis.method.reactive.factorization.L3×3 SparseArrays.SparseMatrixCSC{Float64, Int64} with 5 stored entries:
  1.0         ⋅         ⋅
   ⋅         1.0        ⋅
 -0.492513  -0.189366  1.0
julia> 𝐔₂ = analysis.method.reactive.factorization.U3×3 SparseArrays.SparseMatrixCSC{Float64, Int64} with 5 stored entries:
 -0.495194    ⋅         0.504806
   ⋅        -0.842105   0.157895
   ⋅          ⋅        -0.318121

The vector of increments that corresponds to the reactive power equations can be accessed using:

julia> 𝚫𝐱ₘ = analysis.method.active.increment3-element Vector{Float64}:
  1.0521034860626955e-9
 -1.8809173701904967e-10
 -2.280102039884849e-10

Finally, the solution is updated as follows:

\[ \mathbf x_\text{m}^{(\nu)} = \mathbf x_\text{m}^{(\nu-1)} + \mathbf \Delta \mathbf x_\text{m}^{(\nu-1)}.\]

Again, it is important to note that only the voltage magnitudes of demand buses are updated, while the vector of bus voltage magnitude for all buses is stored:

julia> 𝐕 = analysis.voltage.magnitude4-element Vector{Float64}:
 1.0
 1.005844871456561
 1.0892355535531821
 1.1103697460394555

Gauss-Seidel Method

As elaborated in the Nodal Network Equations section of this manual, each bus is associated with the balance equation expressed as:

\[ \sum_{j = 1}^n Y_{ij} \bar {V}_j = \cfrac{P_i - \text{j}Q_i}{\bar{V}_{i}}, \;\;\; \forall i \in \mathcal{N}.\]

In its expanded form, this can be written as:

\[ \begin{aligned} Y_{11} & \bar{V}_{1} + \cdots+ Y_{1n}\bar{V}_{n} = \frac{{P}_{1} - j{Q}_{1}}{\bar{V}_{1}^*} \\ \; \vdots & \\ Y_{n1} & \bar{V}_{1} + \cdots+ Y_{nn}\bar{V}_{n} = \frac{{P}_{n} - j{Q}_{n}}{\bar{V}_{n}^*}. \end{aligned}\]

While the Gauss-Seidel method directly solves the system of equations, it suffers from very slow convergence, which increases almost linearly with the system size, necessitating numerous iterations to obtain the desired solution [11]. Moreover, the convergence time of the Gauss-Seidel method increases significantly for large-scale systems and can face convergence issues for systems with high active power transfers. Nevertheless, power flow programs utilize both the Gauss-Seidel and Newton-Raphson methods in a complementary manner. Specifically, the Gauss-Seidel method is employed to obtain a quick approximate solution from a "flat start", while the Newton-Raphson method is utilized to obtain the final accurate solution [8].

The Gauss-Seidel method is usually applied to a system of $n$ complex equations, where one represents the slack bus. Consequently, one equation can be eliminated, resulting in a power flow problem with $n-1$ equations.

Initialization

JuliaGrid provides a way to utilize the Gauss-Seidel method for solving the AC power flow problem and determining the magnitudes and angles of bus voltages. To use this method, we need to execute the acModel! function first to set up the system and then initialize the Gauss-Seidel method using the gaussSeidel function. The code snippet below demonstrates this process:

acModel!(system)
analysis = gaussSeidel(system)

This results in the creation of the starting vectors of bus voltage magnitudes $\mathbf{V}^{(0)}$ and angles $\bm{\Theta}^{(0)}$, as shown below:

julia> 𝐕⁽⁰⁾ = analysis.voltage.magnitude4-element Vector{Float64}:
 1.0
 1.0
 1.0
 1.0
julia> 𝚯⁽⁰⁾ = analysis.voltage.angle4-element Vector{Float64}:
 0.0
 0.0
 0.0
 0.0

Iterative Process

JuliaGrid offers the mismatch! and solve! functions to implement the Gauss-Seidel method iterations. These functions are used iteratively until a stopping criterion is met, as shown in the code snippet below:

for iteration = 1:300
    stopping = mismatch!(system, analysis)
    if all(stopping .< 1e-8)
        break
    end
    solve!(system, analysis)
end

In contrast to the Newton-Raphson and fast Newton-Raphson methods, the Gauss-Seidel method does not require the calculation of the mismatch in active and reactive power injection at each iteration. Instead, the mismatch! function is used solely to verify the convergence criteria. At each iteration $\nu = \{1, 2, \dots\}$, we calculate the active power injection mismatch for demand and generator buses, as shown below:

\[ {f}_{P_i}(\mathbf x^{(\nu-1)}) = \Re\{\bar{V}_i^{(\nu - 1)} \bar{I}_i^{*(\nu - 1)}\} - P_i, \;\;\; \forall i \in \mathcal{N}_{\text{pq}} \cup \mathcal{N}_{\text{pv}}.\]

We also compute the reactive power injection mismatch for demand buses, given by:

\[ {f}_{Q_i}(\mathbf x^{(\nu-1)}) = \Im\{\bar{V}_i^{(\nu - 1)} \bar{I}_i^{*(\nu - 1)}\} - Q_i, \;\;\; \forall i \in \mathcal{N}_{\text{pq}}.\]

However, these mismatches are not stored, as they are only used to obtain the maximum absolute values of these mismatches. The maximum values of these mismatches are used as termination criteria for the iteration loop if both are less than a predefined stopping criterion $\epsilon$, as shown below:

\[ \max \{|{f}_{P_i}(\mathbf x^{(\nu-1)})|,\; \forall i \in \mathcal{N}_{\text{pq}} \cup \mathcal{N}_{\text{pv}} \} < \epsilon \\ \max \{|{f}_{Q_i}(\mathbf x^{(\nu-1)})|,\; \forall i \in \mathcal{N}_{\text{pq}} \} < \epsilon.\]

After initializing complex bus voltages $\bar{V}_i^{(0)}$ for all buses in the power system, the function solve! proceeds to compute the voltages for demand buses using the Gauss-Seidel method:

\[ \bar{V}_{i}^{(\nu)} = \cfrac{1}{{Y}_{ii}} \Bigg(\cfrac{{P}_{i} - j{Q}_{i}}{\bar{V}_{i}^{*(\nu-1)}} - \sum\limits_{\substack{j = 1}}^{i - 1} {Y}_{ij}\bar{V}_{j}^{(\nu)} - \sum\limits_{\substack{j = i + 1}}^{n} {Y}_{ij}\bar{V}_{j}^{(\nu-1)}\Bigg), \;\;\; \forall i \in \mathcal{N}_{\text{pq}}.\]

The next step is to determine the solution for generator buses in two stages: first, the reactive power injection is calculated, and then the bus complex voltage is updated using the following equations:

\[ \begin{aligned} Q_i^{(\nu)} &= -\Im \left\{ \bar{V}_{i}^{*(\nu)} \sum\limits_{j=1}^n {Y}_{ij}\bar{V}_{j}^{(\nu)}\right\}, \;\;\; \forall i \in \mathcal{N}_{\text{pv}} \\ \bar{V}_{i}^{(\nu )} &:= \cfrac{1}{{Y}_{ii}} \Bigg(\cfrac{{P}_{i} - j{Q}_{i}^{(\nu)}}{\bar{V}_{i}^{*(\nu )}}- \sum\limits_{\substack{j = 1,\;j \neq i}}^{n} {Y}_{ij}\bar{V}_{j}^{(\nu)} \Bigg), \;\;\; \forall i \in \mathcal{N}_{\text{pv}}. \end{aligned}\]

The obtained voltage magnitude may not be equal to the magnitude specified for the generator bus, so a voltage correction step is necessary:

\[ \bar{V}_{i}^{(\nu)} := {V}_{i}^{(0)} \cfrac{\bar{V}_{i}^{(\nu)}}{{V}_{i}^{(\nu)}}, \;\;\; \forall i \in \mathcal{N}_{\text{pv}}.\]

JuliaGrid stores the final results in vectors that contain all bus voltage magnitudes and angles:

julia> 𝐕 = analysis.voltage.magnitude4-element Vector{Float64}:
 1.0
 1.005844871851792
 1.0892355545361385
 1.1103697470414973
julia> 𝚯 = analysis.voltage.angle4-element Vector{Float64}:
  0.0
 -0.006449009444866342
 -0.0004607248321877602
 -0.004108665753722359

Power Analysis

Once the computation of voltage magnitudes and angles at each bus is completed, various electrical quantities can be determined. JuliaGrid offers the power! function, which enables the calculation of powers associated with buses, branches, and generators. Here is an example code snippet demonstrating its usage:

power!(system, analysis)

The function stores the computed powers in the rectangular coordinate system. It calculates the following powers related to buses, branches, and generators:

Power Injections

Active and reactive power injections are stored as the vectors $\mathbf{P} = [P_i]$ and $\mathbf{Q} = [Q_i]$, respectively, and can be retrieved using the following commands:

julia> 𝐏 = analysis.power.injection.active4-element Vector{Float64}:
 -0.025304583936174074
 -0.21699999730489078
  0.28800000125457986
 -2.4077102614260874e-18
julia> 𝐐 = analysis.power.injection.reactive4-element Vector{Float64}:
 -0.5224650282662511
 -0.12699999649491905
  0.4540000073996279
  9.892532361792298e-21

Generator Power Injections

The power! function in JuliaGrid also computes the active and reactive power injections from the generators at each bus. The active power supplied by the generators to the buses can be calculated by summing the given generator active powers in the input data, except for the slack bus, which can be determined as:

\[ P_{\text{p}i} = P_i + P_{\text{d}i},\;\;\; i \in \mathcal{N}_{\text{sb}},\]

where $P_{\text{d}i}$ represents the active power demanded by consumers at the slack bus. The active power injections from the generators at each bus are stored as the vector, denoted by $\mathbf{P}_{\text{p}} = [P_{\text{p}i}]$, can be obtained using:

julia> 𝐏ₚ = analysis.power.supply.active4-element Vector{Float64}:
 -0.025304583936174074
  0.0
  0.4
  0.0

The calculation of reactive power injection from the generators at generator or slack buses can be achieved using the subsequent equation:

\[ Q_{\text{p}i} = Q_i + Q_{\text{d}i},\;\;\; \forall i \in \mathcal{N}_{\text{pv}} \cup \mathcal{N}_{\text{sb}},\]

where $Q_{\text{d}i}$ represents the reactive power demanded by consumers at the corresponding bus. Further, the reactive power injected by the generators at buses from $\mathcal{N}_{\text{pq}}$ can be calculated by summing the given generator reactive powers in the input data. The vector of these reactive power injections by the generators to the buses, denoted by $\mathbf{Q}_{\text{p}} = [Q_{\text{p}i}]$, can be retrieved using the following command:

julia> 𝐐ₚ = analysis.power.supply.reactive4-element Vector{Float64}:
 -0.5224650282662511
  0.0
  0.424
  0.0

Power at Bus Shunt Elements

Active and reactive powers associated with the shunt elements at each bus are represented by the vectors $\mathbf{P}_{\text{sh}} = [{P}_{\text{sh}i}]$ and $\mathbf{Q}_{\text{sh}} = [{Q}_{\text{sh}i}]$. To retrieve these powers in JuliaGrid, use the following commands:

julia> 𝐏ₛₕ = analysis.power.shunt.active4-element Vector{Float64}:
 0.0
 0.0
 0.0
 0.02589134047804497
julia> 𝐐ₛₕ = analysis.power.shunt.reactive4-element Vector{Float64}:
 -0.0
 -0.0
 -0.0
 -0.014795051701739982

Power Flows

The resulting active and reactive power flows at each from-bus end are stored as the vectors $\mathbf{P}_{\text{i}} = [P_{ij}]$ and $\mathbf{Q}_{\text{i}} = [Q_{ij}],$ respectively, and can be retrieved using the following commands:

julia> 𝐏ᵢ = analysis.power.from.active4-element Vector{Float64}:
  0.0681800941327122
 -0.09348467806888605
 -0.14919987912917734
  0.026012308231465508
julia> 𝐐ᵢ = analysis.power.from.reactive4-element Vector{Float64}:
 -0.11979262337047558
 -0.4026724048957746
 -0.24793254773638929
 -0.25400850645944056

The vectors of active and reactive power flows at the to-bus end are stored as $\mathbf{P}_{\text{j}} = [P_{ji}]$ and $\mathbf{Q}_{\text{j}} = [Q_{ji}]$, respectively, and can be retrieved using the following code:

julia> 𝐏ⱼ = analysis.power.to.active4-element Vector{Float64}:
 -0.06780011817571369
  0.10202890060379553
  0.15995879241931948
 -0.02589134047804497
julia> 𝐐ⱼ = analysis.power.to.reactive4-element Vector{Float64}:
 0.12093255124147115
 0.4385581395423945
 0.26945037431667357
 0.014795051701739979

Power at Branch Shunt Elements

Active and reactive powers associated with the branch shunt elements at each branch are represented by the vectors $\mathbf{P}_{\text{s}} = [P_{\text{s}ij}]$ and $\mathbf{Q}_{\text{s}} = [Q_{\text{s}ij}]$. We can retrieve these values using the following code:

julia> 𝐏ₛ = analysis.power.charging.active4-element Vector{Float64}:
 0.0
 0.0
 0.0
 0.00012096775342053239
julia> 𝐐ₛ = analysis.power.charging.reactive4-element Vector{Float64}:
 -0.0
 -0.0
 -0.0
 -0.2419355068410648

Power at Branch Series Elements

Active and reactive powers associated with the branch series element at each branch are represented by the vectors $\mathbf{P}_{\text{l}} = [P_{\text{l}ij}]$ and $\mathbf{Q}_{\text{l}} = [Q_{\text{l}ij}]$. We can retrieve these values using the following code:

julia> 𝐏ₗ = analysis.power.series.active4-element Vector{Float64}:
  0.0003799759569985204
  0.008544222534909472
  0.010758913290142138
 -1.0842021724855044e-19
julia> 𝐐ₗ = analysis.power.series.reactive4-element Vector{Float64}:
 0.001139927870995561
 0.03588573464661977
 0.021517826580284276
 0.0027220520833633236

Generator Power Outputs

To obtain the output active powers of each generator connected to bus $i \in \mathcal{N}_{\text{pq}} \cup \mathcal{N}_{\text{pv}}$, the given active power in the input data is utilized. For the generator connected to the slack bus, the output active power is determined using the equation:

\[ P_{\text{g}i} = P_i + P_{\text{d}i},\;\;\; i \in \mathcal{N}_{\text{sb}}.\]

In the case of multiple generators connected to the slack bus, the first generator in the input data is assigned the obtained value of $P_{\text{g}i}$. Then, this amount of power is reduced by the output active power of the other generators.

To retrieve the vector of active power outputs of generators, denoted as $\mathbf{P}_{\text{g}} = [P_{\text{g}i}]$, $i \in \mathcal{S}$, where the set $\mathcal{S}$ represents the set of generators, users can utilize the following command:

julia> 𝐏ₒ = analysis.power.generator.active2-element Vector{Float64}:
 -0.025304583936174074
  0.4

The output reactive powers of each generator located at the bus is obtained as:

\[ Q_{\text{g}i} = Q_i + Q_{\text{d}i},\;\;\; i \in \mathcal{N}.\]

If there are multiple generators at the same bus, the reactive power is allocated proportionally among the generators based on their reactive power capabilities.

To retrieve the vector of reactive power outputs of generators, denoted as $\mathbf{Q}_{\text{g}} = [Q_{\text{g}i}]$, $i \in \mathcal{S}$, users can utilize:

julia> 𝐐ₒ = analysis.power.generator.reactive2-element Vector{Float64}:
 -0.5224650282662511
  0.4240000073996279

Current Analysis

JuliaGrid offers the current! function, which enables the calculation of currents associated with buses and branches. Here is an example code snippet demonstrating its usage:

current!(system, analysis)

The function stores the computed currents in the polar coordinate system. It calculates the following currents related to buses and branches:

Current Injections

In JuliaGrid, complex current injections are stored in the vector of magnitudes denoted as $\mathbf{I} = [I_i]$ and the vector of angles represented as $\bm{\psi} = [\psi_i]$. You can retrieve them using the following commands:

julia> 𝐈 = analysis.current.injection.magnitude4-element Vector{Float64}:
 0.5230774586325025
 0.24997084766089062
 0.4935966411616148
 2.168404344971009e-18
julia> 𝛙 = analysis.current.injection.angle4-element Vector{Float64}:
  1.619191576605966
  2.605637769016233
 -1.0059543720614923
  3.141592653589793

Current Flows

To obtain the vectors of magnitudes $\mathbf{I}_{\text{i}} = [I_{ij}]$ and angles $\bm{\psi}_{\text{i}} = [\psi_{ij}]$ for the resulting complex current flows, you can use the following commands:

julia> 𝐈ᵢ = analysis.current.from.magnitude4-element Vector{Float64}:
 0.1378361267952858
 0.41338172516233695
 0.28768189283066065
 0.23441849266339332
julia> 𝛙ᵢ = analysis.current.from.angle4-element Vector{Float64}:
 1.0533688181418361
 1.7989158163819308
 2.1060717428167246
 1.4682841231815587

Similarly, we can obtain the vectors of magnitudes $\mathbf{I}_{\text{j}} = [I_{ji}]$ and angles $\bm{\psi}_{\text{j}} = [\psi_{ji}]$ of the resulting complex current flows using the following code:

julia> 𝐈ⱼ = analysis.current.to.magnitude4-element Vector{Float64}:
 0.1378361267952858
 0.41338172516233695
 0.28768189283066065
 0.026856261289684834
julia> 𝛙ⱼ = analysis.current.to.angle4-element Vector{Float64}:
 -2.088223835447957
 -1.3426768372078626
 -1.0355209107730685
 -2.626555205096993

Current at Branch Series Elements

To obtain the vectors of magnitudes $\mathbf{I}_{\text{l}} = [I_{\text{l}ij}]$ and angles $\bm{\psi}_{\text{l}} = [\psi_{\text{l}ij}]$ of the resulting complex current flows, one can use the following code:

julia> 𝐈ₗ = analysis.current.series.magnitude4-element Vector{Float64}:
 0.13783612679528548
 0.41338172516233634
 0.2876818928306608
 0.12653881253831878
julia> 𝛙ₗ = analysis.current.series.angle4-element Vector{Float64}:
 1.0533688181418344
 1.798915816381931
 2.1060717428167246
 1.3809084790963575