GitHub

PTDF matrix

In this tutorial the methods for computing the Power Transfer Distribution Factors (PTDF) are presented. Before diving into this tutorial we encourage the user to load PowerNetworkMatrices, hit the ? key in the REPL terminal and look for the documention of the different PTDF methods available.

Evaluation of the PTDF matrix

The PTDF matrix can be evaluated according to three different approaches:

  • Dense: considers functions for dense matrix multiplication and inversion
  • KLU: considers functions for sparse matrix multiplication and inversion (default)
  • MKLPardiso: uses Intel's MKL Pardiso solver for sparse matrix operations (only available on Intel-based systems running Linux or Windows)

The evaluation of the PTDF matrix can be easily performed starting from importing the system's data and then by simply calling the PTDF method.

julia> using PowerNetworkMatrices
julia> using PowerSystemCaseBuilder
julia> const PNM = PowerNetworkMatricesPowerNetworkMatrices
julia> const PSB = PowerSystemCaseBuilderPowerSystemCaseBuilder
julia> sys = PSB.build_system(PSB.PSITestSystems, "c_sys5");[ Info: Loaded time series from storage file existing=/home/runner/.julia/packages/PowerSystemCaseBuilder/QOo1f/data/serialized_system/8584b9e729c8aa68ee5405660c6258cde1f67ed3b68f114823707912e9a0d16c/c_sys5_time_series_storage.h5 new=/tmp/jl_Gnt86c compression=InfrastructureSystems.CompressionSettings(false, InfrastructureSystems.CompressionTypesModule.CompressionTypes.DEFLATE = 1, 3, true)
julia> ptdf_1 = PTDF(sys);[ Info: Finding subnetworks via iterative union find
julia> get_ptdf_data(ptdf_1)6×5 transpose(::Matrix{Float64}) with eltype Float64:
 -0.368495  -0.217552  -0.159538  0.0  -0.480452
  0.193917  -0.475895  -0.348989  0.0   0.159538
  0.193917   0.524105   0.651011  0.0   0.159538
  0.437588   0.258343   0.189451  0.0   0.36001
  0.193917   0.524105  -0.348989  0.0   0.159538
  0.368495   0.217552   0.159538  0.0  -0.519548

Advanced users might be interested in computing the PTDF matrix starting from either the data contained in the IncidenceMatrix and BA_matrix structures, or by the information related to the branches and buses of the system.

julia> # evaluate the BA_matrix and Incidence_Matrix
       ba_matrix = BA_Matrix(sys);[ Info: Finding subnetworks via iterative union find
julia> a_matrix = IncidenceMatrix(sys);[ Info: Finding subnetworks via iterative union find
julia> # get the PTDF matrix starting from the values of the
       # previously computed matrices
       ptdf_2 = PTDF(a_matrix, ba_matrix);
julia> get_ptdf_data(ptdf_2)6×5 transpose(::Matrix{Float64}) with eltype Float64:
 -0.368495  -0.217552  -0.159538  0.0  -0.480452
  0.193917  -0.475895  -0.348989  0.0   0.159538
  0.193917   0.524105   0.651011  0.0   0.159538
  0.437588   0.258343   0.189451  0.0   0.36001
  0.193917   0.524105  -0.348989  0.0   0.159538
  0.368495   0.217552   0.159538  0.0  -0.519548

Available methods for the computation of the PTDF matrix

As previously mentioned, the PTDF matrix can be evaluated considering different approaches. The method can be selected by specifying the field linear_solver in the PTDF function.

julia> ptdf_dense = PTDF(sys, linear_solver="Dense");[ Info: Finding subnetworks via iterative union find
julia> get_ptdf_data(ptdf_dense)6×5 transpose(::Matrix{Float64}) with eltype Float64:
 -0.368495  -0.217552  -0.159538  0.0  -0.480452
  0.193917  -0.475895  -0.348989  0.0   0.159538
  0.193917   0.524105   0.651011  0.0   0.159538
  0.437588   0.258343   0.189451  0.0   0.36001
  0.193917   0.524105  -0.348989  0.0   0.159538
  0.368495   0.217552   0.159538  0.0  -0.519548
julia> ptdf_klu = PTDF(sys, linear_solver="KLU");[ Info: Finding subnetworks via iterative union find
julia> get_ptdf_data(ptdf_klu)6×5 transpose(::Matrix{Float64}) with eltype Float64:
 -0.368495  -0.217552  -0.159538  0.0  -0.480452
  0.193917  -0.475895  -0.348989  0.0   0.159538
  0.193917   0.524105   0.651011  0.0   0.159538
  0.437588   0.258343   0.189451  0.0   0.36001
  0.193917   0.524105  -0.348989  0.0   0.159538
  0.368495   0.217552   0.159538  0.0  -0.519548
julia> ptdf_mkl = PTDF(sys, linear_solver="MKLPardiso");[ Info: Finding subnetworks via iterative union find
julia> get_ptdf_data(ptdf_mkl)6×5 transpose(::Matrix{Float64}) with eltype Float64:
 -0.368495  -0.217552  -0.159538  0.0  -0.480452
  0.193917  -0.475895  -0.348989  0.0   0.159538
  0.193917   0.524105   0.651011  0.0   0.159538
  0.437588   0.258343   0.189451  0.0   0.36001
  0.193917   0.524105  -0.348989  0.0   0.159538
  0.368495   0.217552   0.159538  0.0  -0.519548

By default the "KLU" method is selected, which appeared to require significant less time and memory with respect to "Dense". Please note that regardless of which method (KLU, Dense, or MKLPardiso) is used, the resulting PTDF matrix is stored as a dense one.

Note on MKLPardiso: The MKLPardiso solver option is only available on Intel-based systems running Linux or Windows. On other platforms (e.g., ARM-based systems or macOS), use KLU or Dense instead.

Evaluating the PTDF matrix considering distributed slack bus

Whenever needed, the PTDF matrix can be computed with a distributed slack bus. To do so, a vector of type Vector{Float64} needs to be defined, specifying the weights for each bus of the system. These weights identify how the load on the slack bus is redistributed accross the system.

julia> # consider equal distribution accross each bus for this example
       buscount = length(PSY.get_available_components(PSY.ACBus, sys));ERROR: UndefVarError: `PSY` not defined in `Main`
Suggestion: check for spelling errors or missing imports.
julia> dist_slack = 1 / buscount * ones(buscount);ERROR: UndefVarError: `buscount` not defined in `Main`
Suggestion: add an appropriate import or assignment. This global was declared but not assigned.
julia> dist_slack_array = dist_slack / sum(dist_slack);ERROR: UndefVarError: `dist_slack` not defined in `Main`
Suggestion: add an appropriate import or assignment. This global was declared but not assigned.

Once the vector of the weights is defined, the PTDF matrix can be computed by defining the input argument dist_slack (empty array Float64[] by default):

julia> ptdf_distr = PTDF(sys, dist_slack=dist_slack_array);ERROR: UndefVarError: `dist_slack_array` not defined in `Main`
Suggestion: add an appropriate import or assignment. This global was declared but not assigned.

The difference between a the matrix computed with and without the dist_slack field defined can be seen as follows:

julia> # with no distributed slack bus
       get_ptdf_data(ptdf_klu)6×5 transpose(::Matrix{Float64}) with eltype Float64:
 -0.368495  -0.217552  -0.159538  0.0  -0.480452
  0.193917  -0.475895  -0.348989  0.0   0.159538
  0.193917   0.524105   0.651011  0.0   0.159538
  0.437588   0.258343   0.189451  0.0   0.36001
  0.193917   0.524105  -0.348989  0.0   0.159538
  0.368495   0.217552   0.159538  0.0  -0.519548
julia> # with distributed slack bus
       get_ptdf_data(ptdf_distr)ERROR: UndefVarError: `ptdf_distr` not defined in `Main`
Suggestion: check for spelling errors or missing imports.

"Sparse" PTDF matrix

The PTDF matrix can be computed in a "sparse" fashion by defining the input argument tol. If this argument is defined, then elements of the PTDF matrix whose absolute values are below the set threshold are dropped. In addition, the matrix will be stored as a sparse one of type SparseArrays.SparseMatrixCSC{Float64, Int} instead of Matrix{Float64}.

By considering an "extreme" value of 0.2 as tol, the PTDF matrix can be computed as follows:

julia> ptdf_sparse = PTDF(sys, tol=0.2);[ Info: Finding subnetworks via iterative union find
julia> get_ptdf_data(ptdf_sparse)6×5 LinearAlgebra.Transpose{Float64, SparseArrays.SparseMatrixCSC{Float64, Int64}} with 15 stored entries:
 -0.368495  -0.217552    ⋅         ⋅   -0.480452
   ⋅        -0.475895  -0.348989   ⋅     ⋅
   ⋅         0.524105   0.651011   ⋅     ⋅
  0.437588   0.258343    ⋅         ⋅    0.36001
   ⋅         0.524105  -0.348989   ⋅     ⋅
  0.368495   0.217552    ⋅         ⋅   -0.519548

NOTE: 0.2 was used for the purpose of this tutorial. In practice much smaller values are used (e.g., 1e-5).

Network Reductions

The PTDF matrix can be computed with network reductions applied to simplify the system topology. Network reductions eliminate certain buses and branches while preserving the electrical characteristics of the network. This can significantly reduce computation time and memory usage for large systems.

Two types of network reductions are supported:

  • RadialReduction: Eliminates radial (leaf) buses that have only one connection
  • DegreeTwoReduction: Eliminates degree-two buses (buses with exactly two connections) by combining their incident branches

For detailed information about these reductions, see the RadialReduction tutorial and DegreeTwoReduction tutorial.

Using Network Reductions with PTDF

To apply network reductions, pass a vector of NetworkReduction objects to the network_reductions keyword argument:

julia> # Apply radial reduction
       ptdf_radial = PTDF(sys, network_reductions=[RadialReduction()]);ERROR: TypeError: in keyword argument network_reductions, expected Vector{NetworkReduction}, got a value of type Vector{RadialReduction}
julia> # Apply degree-two reduction
       ptdf_degree_two = PTDF(sys, network_reductions=[DegreeTwoReduction()]);ERROR: TypeError: in keyword argument network_reductions, expected Vector{NetworkReduction}, got a value of type Vector{DegreeTwoReduction}
julia> # Combine multiple reductions (order matters - RadialReduction first is recommended)
       ptdf_combined = PTDF(sys, network_reductions=[RadialReduction(), DegreeTwoReduction()]);[ Info: Finding subnetworks via iterative union find

Protecting Specific Buses from Reduction

Both reduction types allow you to specify buses that should not be eliminated using the irreducible_buses parameter:

julia> # Protect specific buses from radial reduction
       reduction = RadialReduction(irreducible_buses=[1, 2])RadialReduction([1, 2])
julia> ptdf_protected = PTDF(sys, network_reductions=[reduction]);ERROR: TypeError: in keyword argument network_reductions, expected Vector{NetworkReduction}, got a value of type Vector{RadialReduction}

DegreeTwoReduction Options

The DegreeTwoReduction has an additional option to control whether buses with reactive power injections are reduced:

julia> # Preserve buses with reactive power injections
       reduction = DegreeTwoReduction(reduce_reactive_power_injectors=false)DegreeTwoReduction(Int64[], false)
julia> ptdf_preserve_reactive = PTDF(sys, network_reductions=[reduction]);ERROR: TypeError: in keyword argument network_reductions, expected Vector{NetworkReduction}, got a value of type Vector{DegreeTwoReduction}

Accessing Reduction Information

After computing the PTDF matrix with reductions, you can access information about what was reduced:

julia> ptdf_reduced = PTDF(sys, network_reductions=[RadialReduction(), DegreeTwoReduction()]);[ Info: Finding subnetworks via iterative union find
julia> # Get the reduction data
       reduction_data = get_network_reduction_data(ptdf_reduced)

Combining Reductions with Other Options

Network reductions can be combined with other PTDF options like distributed slack and sparsification:

julia> ptdf_full_options = PTDF(sys,
           linear_solver="KLU",
           dist_slack=dist_slack_array,
           tol=1e-5,
           network_reductions=[RadialReduction(), DegreeTwoReduction()]
       );ERROR: UndefVarError: `dist_slack_array` not defined in `Main`
Suggestion: add an appropriate import or assignment. This global was declared but not assigned.

NOTE: The reference (slack) bus is automatically protected from elimination during reductions.