Calculations — power-grid-model documentation (2024)

Power flow

Power flow is a “what-if” based grid calculation that will calculate the node voltages and the power flow through the branches, based on assumed load/generation profiles.Some typical use-cases are network planning and contingency analysis.

Input:

• Network data: topology + component attributes

• Assumed load/generation profile

Output:

• Node voltage magnitude and angle

• Power flow through branches

See Power flow algorithms for detailed documentation on the calculation methods.

State estimation

State estimation is a statistical calculation method that determines the most probable state of the grid, based onnetwork data and measurements. Here, measurements can be power flow or voltage values with certain kind of uncertainty, which wereeither measured, estimated or forecasted.

Input:

• Network data: topology + component attributes

• Power flow / voltage measurements with uncertainty

Output:

• Node voltage magnitude and angle

• Power flow through branches

• Deviation between measurement values and estimated state

In order to perform a state estimation, the system should be observable. If the system is not observable, the calculation willraise a sparse matrix error; the matrix reprensentation of the system is too sparse to solve and particular so when it issingular. In short, meeting the requirement of observability indicates that the system is either an overdetermined system (when the number of measurements is larger than the number ofunknowns) or a balanced system (when the number of measurements is equal to the number of unknowns). For each node, there are two unknowns, u and u_angle, so the following conditions should be met:

\[ \begin{eqnarray} n_{measurements} & >= & n_{unknowns} \end{eqnarray}\]

Where

\[ \begin{eqnarray} n_{unknowns} & = & 2 & \cdot & n_{nodes} \end{eqnarray}\]

The number of measurements can be found by taking the sum of the following:

• number of nodes with a voltage sensor with magnitude only

• two times the number of nodes with a voltage sensor with magnitude and angle

• two times the number of nodes without appliances connected

• two times the number of nodes where all connected appliances are measured by a power sensor

• two times the number of branches with a power sensor

Short circuit calculations

Short circuit calculation is carried out to analyze the worst case scenario when a fault has occurred.The currents flowing through branches and node voltages are calculated.Some typical use-cases are selection or design of components like conductors or breakers and power system protection, e.g. relay co-ordination.

Input:

• Network data: topology + component attributes

• Fault type and impedance.

• In the API call: choose between minimum and maximum voltage scaling to calculate the minimum or maximum short circuit currents (according to IEC 60909).

Output:

• Node voltage magnitude and angle

• Current flowing through branches and fault.

Newton-Raphson power flow

Algorithm call: CalculationMethod.newton_raphson

This is the traditional method for power flow calculations. This method uses a Taylor series, ignoring the higher orderterms, to solve the nonlinear set of equations iteratively:

\[ \begin{eqnarray} f(x) & = y \end{eqnarray}\]

Where:

\[\begin{split} \begin{eqnarray} x = \begin{bmatrix} \delta \\ U \end{bmatrix} = \begin{bmatrix} \delta_2 \\ \vdots \\ \delta_N \\ U_2 \\ \vdots \\ U_N \end{bmatrix} \quad\text{and}\quad y = \begin{bmatrix} P \\ Q \end{bmatrix} = \begin{bmatrix} P_2 \\ \vdots \\ P_N \\ Q_2 \\ \vdots \\ Q_N \end{bmatrix} \quad\text{and}\quad f(x) = \begin{bmatrix} P(x) \\ Q(x) \end{bmatrix} = \begin{bmatrix} P_{2}(x) \\ \vdots \\ P_{N}(x) \\ Q_{2}(x) \\ \vdots \\ Q_{N}(x) \end{bmatrix} \end{eqnarray}\end{split}\]

As can be seen in the equations above \(\delta_1\) and \(V_1\) are omitted, because they are known for the slack bus.In each iteration \(i\) the following equation is solved:

\[ \begin{eqnarray} J(i) \Delta x(i) & = \Delta y(i) \end{eqnarray}\]

Where

\[ \begin{eqnarray} \Delta x(i) & = x(i+1) - x(i) \quad\text{and}\quad \Delta y(i) & = y - f(x(i)) \end{eqnarray}\]

\(J\) is the Jacobian, a matrix with all partialderivatives of \(\dfrac{\partial P}{\partial \delta}\), \(\dfrac{\partial P}{\partial U}\), \(\dfrac{\partial Q}{\partial \delta}\)and \(\dfrac{\partial Q}{\partial U}\).

For each iteration the following steps are executed:

• Compute \(\Delta y(i)\)

• Compute the Jacobian \(J(i)\)

• Using LU decomposition, solve \(J(i) \Delta x(i) = \Delta y(i)\) for \(\Delta x(i)\)

• Compute \(x(i+1)\) from \(\Delta x(i) = x(i+1) - x(i)\)

Iterative current power flow

Algorithm call: CalculationMethod.iterative_current

This algorithm is a Jacobi-like method for powerflow analysis.It has a linear convergence rate as opposed to the quadratic convergence rate in the Newton-Raphson method. This means that more iterations is needed to achieve similar convergence.Additionally, Newton-Raphson will be more robust converging in case of greater meshed configurations. Nevertheless, the iterative current algorithm will be faster most of the time.

The algorithm is as follows:

1. Build \(Y_{bus}\) matrix

2. Initialization of \(U_N^0\) to \(1\) plus the intrinsic phase shift of transformers

3. Calculate injected currents: \(I_N^i\) for \(i^{th}\) iteration. The injected currents are calculated as per ZIP model of loads and generation using \(U_N\).\( \begin{eqnarray} I_N = \overline{S_{Z}} \cdot U_{N} + \overline{(\frac{S_{I}}{U_{N}})} \cdot |U_{N}| + \overline{(\frac{S_{P}}{U_N})} \end{eqnarray}\)

4. Solve linear equation: \(YU_N^i = I_N^i\)

5. Check convergence: If maximum voltage deviation from the previous iteration is greater than the tolerance setting (ie. \(u^{(i-1)}_\sigma > u_\epsilon\)), then go back to step 3.

The iterative current algorithm only needs to calculate injected currents before solving linear equations.This is more straightforward than calculating the Jacobian, which was done in the Newton-Raphson algorithm.

Factorizing the matrix of linear equation is the most computationally heavy task.The \(Y_{bus}\) matrix here does not change across iterations which means it only needs to be factorized once to solve the linear equations in all iterations.The \(Y_{bus}\) matrix also remains unchanged in certain batch calculations like timeseries calculations.

Linear power flow

Algorithm call: CalculationMethod.linear

This is an approximation method where we assume that all loads and generations are of constant impedance type regardless of their actual LoadGenType.By doing so, we obtain huge performance benefits as the computation required is equivalent to a single iteration of the iterative methods.It will be more accurate when most of the load/generation types are of constant impedance or the actual node voltages are close to 1 p.u.When all the load/generation types are of constant impedance, the Linear method will be the fastest without loss of accuracy.Therefore power-grid-model will use this method regardless of the input provided by the user in this case.

The algorithm is as follows:

1. Assume injected currents by loads \(I_N=0\) since we model loads/generation as impedance.Admittance of each load/generation is calculated from rated power as \(y=-\overline{S}\)

2. Build \(Y{bus}\). Add the admittances of loads/generation to the diagonal elements.

3. Solve linear equation: \(YU_N = I_N\) for \(U_N\)

Linear current power flow

Algorithm call: CalculationMethod.linear_current

This algorithm is essentially a single iteration of iterative current power flow.

This approximation method will give better results when most of the load/generation types resemble constant current.Similar to iterative current, batch calculations like timeseries will also be faster.Same reason applies here: the \(Y_{bus}\) matrix does not change across batches and as a result the same factorization could be used.

In practical grids most loads and generations correspond to the constant power type. Linear current would give a better approximation than Linear in such case. This is because we approximate the load as current instead of impedance.There is a correlation in voltage error of approximation with respect to the actual voltage for all approximations. They are most accurate when the actual voltages are close to 1 p.u. and the error increases as we deviate from this level.When we approximate the load as impedance at 1 p.u., the voltage error has quadratic relation to the actual voltage. When it is approximated as a current at 1 p.u., the voltage error is only linearly dependent in comparison.

Iterative linear state estimation

Algorithm call: CalculationMethod.iterative_linear

Linear WLS requires all measurements to be linear. This is only possible when all measurements are phasor unit measurements,which is not realistic in distribution grids. Therefore, traditional measurements are linearized prior to running the algorithm:

• Bus voltage: Linear WLS requires a voltage phasor including a phase angle. Given that the phase shift in the distribution grid is very small,it is assumed that the angle shift is zero plus the intrinsic phase shift of transformers. For a certain bus i, the voltagemagnitude measured at that bus is translated into a voltage phasor, where \(\theta_i\) is the intrinsic transformer phase shift:

\[ \begin{eqnarray} \underline{U}_i = U_i \cdot e^{j \theta_i} \end{eqnarray}\]

• Branch/shunt power flow: Linear WLS requires a complex current phasor. To make this translation, the voltage at the terminal shouldalso be measured, otherwise the nominal voltage with zero angle is used as an estimation. With the measured (linearized) voltagephasor, the current phasor is calculated as follows:

\[ \begin{eqnarray} \underline{I} = (\underline{S}/\underline{U})^* \end{eqnarray}\]

• Bus power injection: Linear WLS requires a complex current phasor. Similar as above, if the bus voltage is not measured,the nominal voltage with zero angle will be used as an estimation. The current phasor is calculated as follows:

\[ \begin{eqnarray} \underline{I} = (\underline{S}/\underline{U})^* \end{eqnarray}\]

The aggregated apparent power flow is considered as a single measurement, with variance \(\sigma_S^2 = \sigma_P^2 + \sigma_Q^2\).

The assumption made in the linearization of measurements introduces a system error to the algorithm, because the phase shifts ofbus voltages are ignored in the input measurement data. This error is corrected by applying an iterative approach to the linear WLSalgorithm. In each iteration, the resulted voltage phase angle will be applied as the phase shift of the measured voltage phasorfor the next iteration:

• Initialization: let \(\underline{U}^{(k)}\) be the column vector of the estimated voltage phasor in the k-th iteration. Let Bus \(s\)be the slack bus, which is connected to the external network (source). \(\underline{U}^{(0)}\) is initialized as follows:

  • For bus \(i\), if there is no voltage measurement, assign \(\underline{U}^{(0)} = e^{j \theta_i}\), where \(\theta_i\) is the intrinsictransformer phase shift between Bus \(i\) and Bus \(s\).

  • For bus \(i\), if there is a voltage measurement, assign \(\underline{U}^{(0)} = U_{meas,i}e^{j \theta_i}\), where \(U_{meas,i}\) isthe measured voltage magnitude.

• In iteration \(k\), follow the steps below:

  • Linearize the voltage measurements by using the phase angle of the calculated voltages of iteration \(k-1\).

  • Linearize the complex power flow measurements by using either the linearized voltage measurement if the bus is measured, orthe voltage phasor result from iteration \(k-1\).

  • Compute the temporary new voltage phasor \(\underline{\tilde{U}}^{(k)}\) using the pre-factorized matrix. See also Matrix-prefactorization

  • Normalize the voltage phasor angle by setting the angle of the slack bus to zero:

  • If the maximum deviation between \(\underline{U}^{(k)}\) and \(\underline{U}^{(k-1)}\) is smaller than the error tolerance \(\epsilon\),stop the iteration. Otherwise, continue until the maximum number of iterations is reached.

In the iteration process, the phase angle of voltages at each bus is updated using the last iteration;the system error of the phase shift converges to zero. Because the matrix is pre-built andpre-factorized, the computation cost of each iteration is much smaller than for the Newton-Raphsonmethod, where the Jacobian matrix needs to be constructed and factorized each time.

Newton-Raphson state estimation

Algorithm call: CalculationMethod.newton_raphson

The Newton-Raphson state estimation considers the problem as a system of real, non-linear equations.It iteratively solves the first order Taylor expansion of the system. It does not make any assumptions on linearizationIt could be more accurate and converge in less iterations than the iterative linear approach.However, the Jacobian matrix needs to be calculated every iteration and cannot be prefactorized, which makes it slower on average.

The Newton-Raphson method considers all measurements to be independent real measurements.I.e., \(\sigma_P\), \(\sigma_Q\) and \(\sigma_U\) are all independent values.The rationale behind to calculation is similar to that of the Newton-Raphson for power flow.Consequently, the iteration process differs slightly from that of iterative linear state estimation, as shown below.

• Initialization: let \(\boldsymbol{U}^{(k)}\) be the column vector of the estimated voltage magnitude and \(\boldsymbol{\theta}^{(k)}\) the column vector of theestimated voltage angle in k-th iteration. Let Bus \(s\) be the slack bus which is connected to external network (source).We initialize \(\boldsymbol{U}^{(0)}\) and \(\boldsymbol{\theta}^{(k)}\) as follows:

  • For Bus \(i\), if there is no voltage measurement, assign \(U_i^{(0)} = 1\) and \(\theta_i^{(0)} = 0\),where \(\theta_i\) is the intrinsic transformer phase shift between Bus \(i\) and Bus \(s\).

  • For Bus \(i\), if there is a voltage measurement, \(U_i^{(0)} = U_{\text{mea},i}\) and \(\theta_i^{(0)} = \theta_i\),where \(U_{\text{mea},i}\) is the measured voltage magnitude.

• In iteration \(k\), follow the steps below:

  • Construct and (LU-)factorize the Jacobian matrix for \(\boldsymbol{x}^{(k)}\), using the network parameters.

  • Compute the temporary new voltage magnitude and angle result, \(\widetilde{\boldsymbol{U}}^{(k)}\) and \(\widetilde{\boldsymbol{\theta}}^{(k)}\),using the prefactorized matrix. See also Matrix-prefactorization

  • Normalize the result voltage phasor angle by setting angle of slack bus to zero:\(\underline{U}_i^{(k)} = \widetilde{\underline{U}}_i^{(k)} \times |\widetilde{\underline{U}}_s^{(k)}| / \widetilde{\underline{U}}_s^{(k)}\).

  • If the maximum deviation between \(\underline{\boldsymbol{U}}^{(k)}\) and \(\underline{\boldsymbol{U}}^{(k-1)}\) is smaller than the tolerance \(\epsilon\),stop the iteration, otherwise continue until the maximum number of iterations is reached.Note: we’re using the phasor here: \(\underline{\boldsymbol{U}} = \boldsymbol{U}e^{j\boldsymbol{\theta}}\).

As for the iterative linear approach, during iterations, phase angles of voltage at each bus are updated using ones from the previous iteration.The system error of the phase shift converges to zero.

IEC 60909 short circuit calculation

Algorithm call: CalculationMethod.iec60909

The assumptions used for calculations in power-grid-model are aligned to the ones mentioned in IEC 60909.

• The state of the grid with respect to loads and generations are ignored for the short circuit calculation. (Note: Shunt admittances are included in calculation.)

• The pre-fault voltage is considered in the calculation and is calculated based on the grid parameters and topology. (Excl. loads and generation)

• The calculations are assumed to be time-independent. (Voltages are sine throughout with the fault occurring at a zero crossing of the voltage, the complexity of rotating machines and harmonics are neglected, etc.)

• To account for the different operational conditions, a voltage scaling factor of c is applied to the voltage source while running short circuit calculation function.The factor c is determined by the nominal voltage of the node that the source is connected to and the API option to calculate the minimum or maximum short circuit currents.The table to derive c according to IEC 60909 is shown below.

There are 4 types of fault situations that can occur in the grid, along with the following possible combinations of the associated phases:

• Three-phase to ground: abc

• Single phase to ground: a, b, c

• Two phase: ab, bc, ac

• Two phase to ground: ab, bc, ac

Power flow with automatic tap changing

Some of the most important regulators in the grid affect the tap position of transformers.These Transformer tap regulators try to regulate a control voltage \(U_{\text{control}}\) such that it is within a specified voltage band.The \(U_{\text{control}}\) may be compensated for the voltage drop during transport.Power flow calculations that take the behavior of these regulators into account may be toggled by providing one of the following strategies to the tap_changing_strategy option.

Example: dependent batch update

# 3 scenarios, 3 objects (lines)# for each scenario, only one line is specifiedline_update = initialize_array('update', 'line', (3, 1))# set the mutations for each scenario: disable one of the three linesfor component_update, component_id in zip(line_update, (3, 5, 8)): component_update['id'] = component_id component_update['from_status'] = 0 component_update['to_status'] = 0non_independent_update_data = {'line': line_update}

Example: full batch update data

# 3 scenarios, 3 objects (lines)# for each scenario, all lines are specifiedline_update = initialize_array('update', 'line', (3, 3))# use broadcasting to specify the default stateline_update['id'] = [[3, 5, 8]]line_update['from_status'] = 1line_update['to_status'] = 1# set the mutations for each scenario: disable one of the three linesfor component_idx, scenario in enumerate(line_update): component = scenario[component_idx] component['from_status'] = 0 component['to_status'] = 0independent_update_data = {'line': line_update}