πŸ“ Linear analysis

Linear analysis in VeraGrid provides a fast approximation of the power flow and power transfer impacts using sensitivity factors instead of solving the full non-linear power flow every time. It is especially useful for screening studies, contingency analysis, time-series approximations, and transfer sensitivity calculations.

The most important outputs are:

  • PTDF: Power Transfer Distribution Factors for branches HvdcLine and VSC.

  • LODF: Line Outage Distribution Factors for branches HvdcLine and VSC.

  • Pf: Lossless Power flows through the branches.

  • Pbus: Modified injections to match the PTDF formation. (i.e. with distributed slack, the SLACK and PV nodes get modified)

These factors describe how injections or outages propagate through branch flows.

It is not intended to replace the non-linear power flow in all situations. Its main value is speed and interpretability.

Interpretation Notes

  • Large PTDF absolute values indicate strong sensitivity of a branch to injections at a bus.

  • Large LODF absolute values indicate strong flow redistribution after the outage of another branch.

  • Linear branch loadings are very useful for screening, but not all overloads found linearly will match the exact non-linear AC solution.

  • Time-series linear analysis is often good for ranking and trend analysis, even when exact MW values differ somewhat from Newton-Raphson.

Options

The linear analysis options are:

  • distribute_slack: whether the active-power mismatch is distributed across PV and slack buses instead of concentrated at one reference

  • correct_values: whether the analysis should apply internal corrections for out-of-range values

  • ptdf_threshold: PTDF sparsification threshold

  • lodf_threshold: LODF sparsification threshold

For most workflows, distribute_slack and correct_values are the most relevant settings.

Theory

Linear analysis uses a lossless active-power approximation to produce two main sensitivity matrices:

  • PTDF, the Power Transfer Distribution Factor matrix, which maps bus active-power injections to branch active-power flows.

  • LODF, the Line Outage Distribution Factor matrix, which maps the pre-contingency flow of an outaged branch to the flow changes on the remaining branches.

The formulas below are applied independently to each connected island. A slack correction must never transfer active-power mismatch between electrically disconnected islands.

Network Matrices

For one island, let:

  • n be the number of buses in the island.

  • m be the number of passive branches with impedance in the island.

  • P \in \mathbb{R}^{n} be the bus active-power injection vector.

  • P_f \in \mathbb{R}^{m} be the branch active-power flow vector, oriented from the branch from bus to the branch to bus.

  • \theta \in \mathbb{R}^{n} be the bus voltage-angle vector.

Define the branch-bus incidence matrix A \in \mathbb{R}^{m \times n} as:

A_{k i} =
\begin{cases}
1, & i = F_k \\
-1, & i = T_k \\
0, & \text{otherwise}
\end{cases}

where F_k and T_k are the from and to buses of branch k.

With this orientation, the linear branch flow model is:

P_f = B_f \theta

and the nodal active-power balance is:

P = A^{\mathsf{T}} P_f = B_{\mathrm{bus}} \theta

where:

B_f = \operatorname{diag}(b) A

and:

B_{\mathrm{bus}} = A^{\mathsf{T}} \operatorname{diag}(b) A

Here b_k is the active-power sensitivity of branch k with respect to its angle difference. In the simplest DC model, b_k = 1 / x_k. In the actual numerical implementation, transformer tap effects and other branch parameters are included in the produced B_f and B_{\mathrm{bus}} matrices.

Single-Reference PTDF

Because B_{\mathrm{bus}} is singular for a connected lossless island, one bus angle is used as the reference. Let s be the selected reference bus and let \mathcal{N} be the set of non-reference buses.

The reduced nodal equation is:

P_{\mathcal{N}} = B_{\mathrm{bus}, \mathcal{N}\mathcal{N}} \theta_{\mathcal{N}}

so:

\theta_{\mathcal{N}} = B_{\mathrm{bus}, \mathcal{N}\mathcal{N}}^{-1} P_{\mathcal{N}}

The single-reference PTDF matrix H^{(0)} \in \mathbb{R}^{m \times n} is then built as:

H^{(0)}_{:, \mathcal{N}} = B_{f, :, \mathcal{N}} B_{\mathrm{bus}, \mathcal{N}\mathcal{N}}^{-1}

and:

H^{(0)}_{:, s} = 0

Thus, for a raw injection vector P:

P_f = H^{(0)} P

In single-slack mode, this flow is balanced by assigning the island active-power mismatch to the reference slack bus. For nodal-balance checks, define:

\Delta P = \mathbf{1}^{\mathsf{T}} P

and:

w_i^{\mathrm{slack}} =
\begin{cases}
1, & i = s \\
0, & i \ne s
\end{cases}

The effective injection vector is:

P^{\mathrm{eff}} = P - w^{\mathrm{slack}} \Delta P

or, component-wise:

P_i^{\mathrm{eff}} = P_i \quad \text{for } i \ne s

P_s^{\mathrm{eff}} = P_s - \sum_{i=1}^{n} P_i = - \sum_{i \ne s} P_i

Therefore, when distribute_slack = False, the PTDF production matrix is not modified:

D_{\mathrm{PTDF}} = I

H = H^{(0)}

but the Kirchhoff balance check must compare branch-flow divergence against P^{\mathrm{eff}}, not necessarily against the raw injection vector P.

Distributed-Slack PTDF

When distribute_slack = True, the PTDF is modified so that the balancing injection is distributed only over buses that can regulate active power in the linear model.

The numerical bus type convention is:

Code

Bus type

1

PQ

2

PV

3

Slack

The latest distributed-slack convention uses PV and Slack buses only. Define the participating set:

\mathcal{G} = \left\{ i \;\middle|\; \operatorname{bus\_type}_i = 2 \;\lor\; \operatorname{bus\_type}_i = 3 \right\}

The participation vector w \in \mathbb{R}^{n} is:

w_i =
\begin{cases}
\dfrac{1}{|\mathcal{G}|}, & i \in \mathcal{G} \\
0, & i \notin \mathcal{G}
\end{cases}

so that:

\mathbf{1}^{\mathsf{T}} w = 1

The distributed-slack transaction matrix is:

D = I - w \mathbf{1}^{\mathsf{T}}

or, entry-wise:

D_{ij} = \delta_{ij} - w_i

where \delta_{ij} is the Kronecker delta. Only the rows of PV and Slack buses are modified, because only those buses have nonzero participation weights.

The effective injection vector is:

P^{\mathrm{eff}} = D P

therefore:

P^{\mathrm{eff}} = P - w \left( \mathbf{1}^{\mathsf{T}} P \right)

and:

\sum_{i=1}^{n} P_i^{\mathrm{eff}} = 0

The distributed-slack PTDF is produced by right-multiplying the single-reference PTDF by D:

H = H^{(0)} D

The branch flows are then computed from the raw injection vector as:

P_f = H P = H^{(0)} D P = H^{(0)} P^{\mathrm{eff}}

Each PTDF column is therefore a balanced transaction. For a unit injection at bus j:

D e_j = e_j - w

so column j of the distributed-slack PTDF represents an injection at bus j balanced by withdrawals over the PV and Slack buses according to w.

This replaces the older all-bus convention:

D_{ij}^{\mathrm{old}} =
\begin{cases}
1, & i = j \\
-\dfrac{1}{n - 1}, & i \ne j
\end{cases}

The older convention assigned the balancing power to all other buses, including PQ buses, and excluded the injection bus itself from the balancing set. The new convention uses a fixed participation vector based on bus_types, so PQ buses do not absorb active-power mismatch.

PTDF Production Summary

For each island, the PTDF production can be summarized as:

H^{(0)}_{:, \mathcal{N}} = B_{f, :, \mathcal{N}} B_{\mathrm{bus}, \mathcal{N}\mathcal{N}}^{-1}

H^{(0)}_{:, s} = 0

If distribute_slack = False:

H = H^{(0)}

If distribute_slack = True:

D = I - w \mathbf{1}^{\mathsf{T}}

H = H^{(0)} D

The resulting PTDF maps bus injections to branch flows:

P_f = H P

For distributed slack, this is equivalent to:

P_f = H^{(0)} P^{\mathrm{eff}}

with:

P^{\mathrm{eff}} = D P

Nodal-Balance Check

With the branch-flow orientation defined by A, branch-flow divergence is:

\operatorname{div}(P_f) = A^{\mathsf{T}} P_f

The linear Kirchhoff balance check is:

A^{\mathsf{T}} P_f - P^{\mathrm{eff}} \approx 0

For distribute_slack = True:

P^{\mathrm{eff}} = D P

For distribute_slack = False with one reference slack bus s:

P^{\mathrm{eff}} = P - w^{\mathrm{slack}} \left( \mathbf{1}^{\mathsf{T}} P \right)

This distinction is important: the vector used to compute flows and the vector used to check nodal balance are not always the same raw vector. If the original P comes from an AC operating point, it may contain active losses or mismatch that must be absorbed by the slack convention in the lossless linear model.

LODF Production

The LODF matrix is produced from branch-to-branch PTDFs.

For branch k, define the endpoint transaction vector:

a_k = e_{F_k} - e_{T_k}

where e_{F_k} and e_{T_k} are unit vectors at the from and to buses of branch k. Collecting all endpoint transaction vectors gives:

A^{\mathsf{T}} = \begin{bmatrix} a_1 & a_2 & \cdots & a_m \end{bmatrix}

The branch-to-branch PTDF matrix is:

K = H A^{\mathsf{T}}

so:

K_{\ell k} = H_{\ell, :} \left( e_{F_k} - e_{T_k} \right)

K_{\ell k} is the sensitivity of branch \ell to a unit transaction from the from bus of branch k to the to bus of branch k.

For the outage of branch k, the LODF denominator is:

d_k = 1 - K_{kk}

For non-islanding outages, the LODF entries are:

\operatorname{LODF}_{\ell k} = \frac{K_{\ell k}}{1 - K_{kk}} \quad \text{for } \ell \ne k

and the diagonal is set to:

\operatorname{LODF}_{kk} = -1

The post-contingency flow approximation for outage k is then:

P_{f, \ell}^{\mathrm{post}(k)} = P_{f, \ell}^{\mathrm{pre}} + \operatorname{LODF}_{\ell k} P_{f, k}^{\mathrm{pre}}

for \ell \ne k, and:

P_{f, k}^{\mathrm{post}(k)} = 0

If:

1 - K_{kk} \approx 0

then the outaged branch is island-forming or numerically close to island-forming. In that case, the standard LODF value is not finite, because removing the branch changes network connectivity. Production code should mark or correct those values according to the selected correction policy.

The branch endpoint transaction a_k is balanced because:

\mathbf{1}^{\mathsf{T}} a_k = 0

Therefore, for the corrected distributed-slack matrix:

D a_k = a_k

This means the LODF calculation is independent of the distributed-slack participation vector, provided the PTDF uses the consistent formulation:

H = H^{(0)} D

The slack convention matters for mapping arbitrary bus injections P to flows, but endpoint-to-endpoint branch outage transactions are already balanced.

Avoiding Double Correction

There are two valid usage patterns.

If the PTDF already includes the slack-distribution matrix:

H = H^{(0)} D

then compute flows from the raw injection vector:

P_f = H P

and check nodal balance with:

P^{\mathrm{eff}} = D P

Alternatively, if using the single-reference PTDF externally:

H = H^{(0)}

then first correct the injection vector:

P^{\mathrm{eff}} = D P

and then compute:

P_f = H^{(0)} P^{\mathrm{eff}}

Do not apply the slack matrix to both the PTDF and the input vector, because that applies the same slack correction twice.

SRAP: Automatic Power Reduction System

The Automatic Power Reduction System, or SRAP, is a linear sensitivity-based mechanism used in contingency analysis to determine whether an overload can be dismissed by corrective generation re-dispatch.

The idea is simple:

  • a contingency produces an overload on a monitored branch

  • a set of generators has known sensitivities with respect to that overload

  • available upward or downward redispatch is combined with those sensitivities

  • if the overload can be removed within the allowed corrective margin, the contingency can be dismissed

Example:

Imagine that a line overloads by 5 MW after a contingency. Three plants are identified as significant for relieving that overload:

  • plant 1: generating 80 MW, with PTDF sensitivity 0.11

  • plant 2: generating 30 MW, with PTDF sensitivity 0.09

  • plant 3: generating 50 MW, with PTDF sensitivity 0.07

Assume the SRAP limit is 90 MW of redispatch.

We build arrays ordered by sensitivity:

  • sensitivity = [0.11, 0.09, 0.07]

  • p_available = [80, 30, 30]

Then we compute the contribution f = sensitivity * p_available.

If the resulting corrective capability is larger than the overload, the contingency can be considered solvable by SRAP and may be dismissed from the critical set.

A contingency study with SRAP activated can be run as follows:

con_options = ContingencyAnalysisOptions()
con_options.use_srap = True
con_options.engine = ContingencyEngine.Linear

con_drv = ContingencyAnalysisDriver(
    grid=grid,
    options=con_options,
    engine=EngineType.VeraGrid
)

con_drv.run()

API

Snapshot Linear Analysis

This is the standard single-state linear study.

Using the simplified API:

import os
import VeraGridEngine as gce

folder = os.path.join('..', 'Grids_and_profiles', 'grids')
fname = os.path.join(folder, 'IEEE 5 Bus.xlsx')

main_circuit = gce.open_file(fname)

results = gce.linear_power_flow(grid=main_circuit)

print("Bus results:\n", results.get_bus_df())
print("Branch results:\n", results.get_branch_df())
print("PTDF:\n", results.mdl(gce.ResultTypes.PTDF).to_df())
print("LODF:\n", results.mdl(gce.ResultTypes.LODF).to_df())

Using the driver directly:

import os
import VeraGridEngine as gce

folder = os.path.join('..', 'Grids_and_profiles', 'grids')
fname = os.path.join(folder, 'IEEE 5 Bus.xlsx')

main_circuit = gce.open_file(fname)

options_ = gce.LinearAnalysisOptions(distribute_slack=False, correct_values=True)

drv = gce.LinearAnalysisDriver(grid=main_circuit, options=options_)
drv.run()

print("Bus results:\n", drv.results.get_bus_df())
print("Branch results:\n", drv.results.get_branch_df())
print("PTDF:\n", drv.results.mdl(gce.ResultTypes.PTDF).to_df())
print("LODF:\n", drv.results.mdl(gce.ResultTypes.LODF).to_df())

Output:

Bus results:
        Vm   Va       P    Q
Bus 0  1.0  0.0  2.1000  0.0
Bus 1  1.0  0.0 -3.0000  0.0
Bus 2  1.0  0.0  0.2349  0.0
Bus 3  1.0  0.0 -0.9999  0.0
Bus 4  1.0  0.0  4.6651  0.0

Branch results:
                  Pf   loading
Branch 0-1  2.497192  0.624298
Branch 0-3  1.867892  0.832394
Branch 0-4 -2.265084 -0.828791
Branch 1-2 -0.502808 -0.391900
Branch 2-3 -0.267908 -0.774300
Branch 3-4 -2.400016 -1.000006

PTDF:
               Bus 0     Bus 1     Bus 2  Bus 3     Bus 4
Branch 0-1  0.193917 -0.475895 -0.348989    0.0  0.159538
Branch 0-3  0.437588  0.258343  0.189451    0.0  0.360010
Branch 0-4  0.368495  0.217552  0.159538    0.0 -0.519548
Branch 1-2  0.193917  0.524105 -0.348989    0.0  0.159538
Branch 2-3  0.193917  0.524105  0.651011    0.0  0.159538
Branch 3-4 -0.368495 -0.217552 -0.159538    0.0 -0.480452

LODF:
            Branch 0-1  Branch 0-3  Branch 0-4  Branch 1-2  Branch 2-3  Branch 3-4
Branch 0-1   -1.000000    0.344795    0.307071   -1.000000   -1.000000   -0.307071
Branch 0-3    0.542857   -1.000000    0.692929    0.542857    0.542857   -0.692929
Branch 0-4    0.457143    0.655205   -1.000000    0.457143    0.457143    1.000000
Branch 1-2   -1.000000    0.344795    0.307071   -1.000000   -1.000000   -0.307071
Branch 2-3   -1.000000    0.344795    0.307071   -1.000000   -1.000000   -0.307071
Branch 3-4   -0.457143   -0.655205    1.000000   -0.457143   -0.457143   -1.000000

Time-Series Linear Analysis

The time-series driver uses the same linear sensitivity idea to approximate branch flows over many time steps.

This is useful when:

  • the grid topology does not change significantly

  • you want a fast approximation of flow evolution over time

  • you want to compare linear and non-linear flow tracking

The time-series driver is LinearAnalysisTimeSeriesDriver.

import VeraGridEngine as gce

grid = gce.open_file("IEEE39_1W.veragrid")

drv = gce.LinearAnalysisTimeSeriesDriver(grid=grid)
drv.run()

print(drv.results.Sf)
print(drv.results.loading)

Benchmark

Linear Algebra Frameworks Benchmark

IEEE 39 1-year time series

The experiment measures the time taken by the time-series simulation using different linear algebra solvers.

The power flow tolerance is set to 1e-4.

The time in seconds taken using each solver is:

KLU

LAPACK

ILU

SuperLU

Pardiso

Test 1

82.03

82.10

81.79

82.88

93.23

Test 2

80.22

80.84

81.71

81.37

95.29

Test 3

79.53

82.32

82.75

80.98

92.62

Test 4

80.06

82.66

82.14

80.17

97.60

Test 5

80.07

80.51

81.94

80.03

93.39

Average

80.38

81.68

82.07

81.09

94.42

2869 Pegase 1-week time series

The experiment measures the time taken by the time-series simulation using different linear algebra solvers.

The power flow tolerance is set to 1e-4.

The time in seconds taken using each solver is:

KLU

LAPACK

ILU

SuperLU

Pardiso

Test 1

2.46

2.50

2.52

2.48

2.54

Test 2

2.35

2.31

2.36

2.32

2.59

Test 3

2.40

2.42

2.46

2.46

2.46

Test 4

2.33

2.31

2.34

2.33

2.42

Test 5

2.31

2.32

2.45

2.33

2.51

Average

2.37

2.37

2.43

2.39

2.51

From these tests, the solvers are roughly equivalent for this type of simulation, except Pardiso, which performs worse than the others in these specific benchmarks.

Linear vs Non-Linear Comparison

One of the most common uses of the linear model is to compare it against a Newton-Raphson time-series power flow.

import os
from matplotlib import pyplot as plt
import VeraGridEngine as gce

plt.style.use('fivethirtyeight')

folder = os.path.join('..', 'Grids_and_profiles', 'grids')
fname = os.path.join(folder, 'IEEE39_1W.veragrid')
main_circuit = gce.open_file(fname)

ptdf_driver = gce.LinearAnalysisTimeSeriesDriver(grid=main_circuit)
ptdf_driver.run()

pf_options_ = gce.PowerFlowOptions(solver_type=gce.SolverType.NR)
ts_driver = gce.PowerFlowTimeSeriesDriver(grid=main_circuit, options=pf_options_)
ts_driver.run()

fig = plt.figure(figsize=(30, 6))
ax1 = fig.add_subplot(131)
ax1.set_title('Newton-Raphson based flow')
ax1.plot(ts_driver.results.Sf.real)
ax1.set_ylabel('MW')
ax1.set_xlabel('Time')

ax2 = fig.add_subplot(132)
ax2.set_title('PTDF based flow')
ax2.plot(ptdf_driver.results.Sf.real)
ax2.set_ylabel('MW')
ax2.set_xlabel('Time')

ax3 = fig.add_subplot(133)
ax3.set_title('Difference')
diff = ts_driver.results.Sf.real - ptdf_driver.results.Sf.real
ax3.plot(diff)
ax3.set_ylabel('MW')
ax3.set_xlabel('Time')

fig.set_tight_layout(tight=True)

plt.show()

PTDF flows comparison.png

This comparison is important because it shows where the linear model is accurate enough for screening and where the full non-linear power flow is still necessary.

Registered Result Properties

LinearAnalysisResults registered properties

The snapshot linear analysis result stores sensitivity matrices and the base linearized operating point.

Property

Type

Description

branch_names

StrVec

Names aligned with branch-indexed result arrays.

bus_names

StrVec

Names aligned with bus-indexed result arrays.

bus_types

IntVec

Bus type code used by the solved numerical model.

PTDF

Mat

Power transfer distribution factor matrix.

LODF

Mat

Line outage distribution factor matrix.

HvdcDF

Mat

HVDC power-transfer sensitivity matrix.

HvdcODF

Mat

HVDC outage distribution factor matrix.

VscDF

Mat

VSC power-transfer sensitivity matrix.

VscODF

Mat

VSC outage distribution factor matrix.

Sf

Vec

Complex branch power flow at the from side.

Sbus

Vec

Complex bus power injection.

voltage

CxVec

Complex bus voltage solution.

loading

Vec

Branch loading result.

The result model also exposes tabular views through mdl() for:

  • ResultTypes.PTDF

  • ResultTypes.LODF

  • ResultTypes.HvdcPTDF

  • ResultTypes.HvdcODF

  • ResultTypes.VscPTDF

  • ResultTypes.VscODF

  • ResultTypes.BranchActivePowerFrom

  • ResultTypes.BranchLoading

    • ResultTypes.BranchActivePower

LinearAnalysisTimeSeriesResults registered properties

The time-series linear analysis result stores linearized voltages, injections, flows, losses, and loading over time.

Property

Type

Description

branch_names

StrVec

Names aligned with branch-indexed result arrays.

bus_names

StrVec

Names aligned with bus-indexed result arrays.

bus_types

IntVec

Bus type code used by the solved numerical model.

voltage

CxMat

Complex bus voltage solution.

Sf

CxMat

Complex branch power flow at the from side.

S

CxMat

Complex bus power result matrix.

losses

CxMat

Complex branch losses.

loading

CxMat

Branch loading result.