All the explanations on this section are implemented here.

Canonical Newton-Raphson

The Newton-Raphson method is the standard power flow method tough at schools. GridCal implements slight but important modifications of this method that turns it into a more robust, industry-standard algorithm. The Newton-Raphson method is the first order Taylor approximation of the power flow equation.

The expression to update the voltage solution in the Newton-Raphson algorithm is the following:

\textbf{V}_{t+1} = \textbf{V}_t + \textbf{J}^{-1}(\textbf{S}_0 - \textbf{S}_{calc})


  • \textbf{V}_t: Voltage vector at the iteration t (current voltage).
  • \textbf{V}_{t+1}: Voltage vector at the iteration t+1 (new voltage).
  • \textbf{J}: Jacobian matrix.
  • \textbf{S}_0: Specified power vector.
  • \textbf{S}_{calc}: Calculated power vector using \textbf{V}_t.

In matrix form we get:

\textbf{J}_{11} & \textbf{J}_{12} \\
\textbf{J}_{21} & \textbf{J}_{22} \\
\Delta \textbf{P}\\
\Delta \textbf{Q}\\

Jacobian in power equations

The Jacobian matrix is the derivative of the power flow equation for a given voltage set of values.

\textbf{J} =
\textbf{J}_{11} & \textbf{J}_{12} \\
\textbf{J}_{21} & \textbf{J}_{22} \\


  • J11 = Re\left(\frac{\partial \textbf{S}}{\partial \theta}[pvpq, pvpq]\right)
  • J12 = Re\left(\frac{\partial \textbf{S}}{\partial |V|}[pvpq, pq]\right)
  • J21 = Im\left(\frac{\partial \textbf{S}}{\partial \theta}[pq, pvpq]\right)
  • J22 = Im\left(\frac{\partial \textbf{S}}{\partial |V|}[pq pq]\right)


  • \textbf{S} = \textbf{V} \cdot \left(\textbf{I} + \textbf{Y}_{bus} \times \textbf{V} \right)^*

Here we introduced two complex-valued derivatives:

  • \frac{\partial S}{\partial |V|} = V_{diag} \cdot \left(Y_{bus} \times V_{diag,norm} \right)^* + I_{diag}^* \cdot V_{diag,norm}
  • \frac{\partial S}{\partial \theta} =  1j \cdot V_{diag} \cdot \left(I_{diag} - Y_{bus} \times V_{diag} \right)^*


  • V_{diag}: Diagonal matrix formed by a voltage solution.
  • V_{diag,norm}: Diagonal matrix formed by a voltage solution, where every voltage is divided by its module.
  • I_{diag}: Diagonal matrix formed by the current injections.
  • Y_{bus}: And of course, this is the circuit full admittance matrix.

This Jacobian form can be used for other methods.


In 1982 S. Iwamoto and Y. Tamura present a method [1] where the Jacobian matrix J is only computed at the beginning, and the iteration control parameter µ is computed on every iteration. In GridCal, J and µ are computed on every iteration getting a more robust method on the expense of a greater computational effort.

The Iwamoto and Tamura’s modification to the Newton-Raphson algorithm consist in finding an optimal acceleration parameter µ that determines the length of the iteration step such that, the very iteration step does not affect negatively the solution process, which is one of the main drawbacks of the Newton-Raphson method:

\textbf{V}_{t+1} = \textbf{V}_t + \mu \cdot \textbf{J}^{-1}\times (\textbf{S}_0 - \textbf{S}_{calc})

Here µ is the Iwamoto optimal step size parameter. To compute the parameter µ we must do the following:

\textbf{J'} = Jacobian(\textbf{Y}, \textbf{dV})

The matrix \textbf{J'} is the Jacobian matrix computed using the voltage derivative numerically computed as the voltage increment \textbf{dV}= \textbf{V}_{t} - \textbf{V}_{t-1} (voltage difference between the current and the previous iteration).

\textbf{dx} = \textbf{J}^{-1} \times  (\textbf{S}_0 - \textbf{S}_{calc})

\textbf{a} = \textbf{S}_0 - \textbf{S}_{calc}

\textbf{b} = \textbf{J} \times \textbf{dx}

\textbf{c} = \frac{1}{2} \textbf{dx} \cdot (\textbf{J'} \times \textbf{dx})

g_0 = -\textbf{a} \cdot \textbf{b}

g_1 = \textbf{b} \cdot \textbf{b} + 2  \textbf{a} \cdot \textbf{c}

g_2 = -3  \textbf{b} \cdot \textbf{c}

g_3 = 2  \textbf{c} \cdot \textbf{c}

G(x) = g_0 + g_1 \cdot x + g_2 \cdot x^2 + g_3 \cdot x^3

µ = solve(G(x), x_0=1)

There will be three solutions to the polynomial G(x). Only the last solution will be real, and therefore it is the only valid value for µ. The polynomial can be solved numerically using 1 as the seed.

[1]Iwamoto, S., and Y. Tamura. “A load flow calculation method for ill-conditioned power systems.”IEEE transactions on power apparatus and systems 4 (1981): 1736-1743.

Newton-Raphson in Current Equations

Newton-Raphson in current equations is similar to the regular Newton-Raphson algorithm presented before, but the mismatch is computed with the current instead of power.

The Jacobian is then:

Re\left\{\left[\frac{\partial I}{\partial \theta}\right]\right\}_{(pqpv, pqpv)} &
Re\left\{\left[\frac{\partial I}{\partial Vm}\right]\right\}_{(pqpv, pq)} \\
Im\left\{\left[\frac{\partial I}{\partial \theta}\right]\right\}_{(pq, pqpv)} &
Im\left\{\left[\frac{\partial I}{\partial Vm}\right]\right\}_{(pq,pq)}


\left[\frac{\partial I}{\partial Vm}\right] = [Y] \times [E_{diag}]

\left[\frac{\partial I}{\partial \theta}\right] = 1j \cdot [Y] \times [V_{diag}]

The mismatch is computed as increments of current:

F = \left[
 Re\{\Delta I\} \\
 Im\{\Delta I\}


[\Delta I] = \left( \frac{S_{specified}}{V} \right)^*  - ([Y] \times [V] - [I^{specified}])

The steps of the algorithm are equal to the the algorithm presented in Newton-Raphson.