Merge pull request #31 from invenia/sm/ptdf

Add PTDF and LODF calculation functions

Project.toml

@@ -1,14 +1,16 @@
 name = "FullNetworkSystems"
 uuid = "877b7152-b508-43dc-81fb-72341a693988"
 authors = ["Invenia Technical Computing Corporation"]
-version = "1.6.0"
+version = "1.7.0"
 
 [deps]
 AxisKeys = "94b1ba4f-4ee9-5380-92f1-94cde586c3c5"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 Dictionaries = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 InlineStrings = "842dd82b-1e85-43dc-bf29-5d0ee9dffc48"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
 AxisKeys = "0.2"

src/FullNetworkSystems.jl

@@ -5,6 +5,8 @@ using Dates
 using Dictionaries
 using DocStringExtensions
 using InlineStrings
+using LinearAlgebra
+using SparseArrays
 
 export System, SystemDA, SystemRT
 export Zone, Generator, Bus, Branch
@@ -14,7 +16,7 @@ export gens_per_zone, branches_by_breakpoints, get_datetimes
 export get_zones, get_buses, get_generators, get_branches, get_lines, get_transformers
 export get_regulation_requirements, get_operating_reserve_requirements, get_good_utility_requirements
 export get_gens_per_bus, get_loads_per_bus, get_incs_per_bus, get_decs_per_bus, get_psls_per_bus
-export get_ptdf, get_lodfs
+export compute_ptdf, get_ptdf, retrieve_ptdf, compute_lodf, get_lodfs
 export get_initial_commitment, get_initial_downtime, get_initial_uptime
 export get_increments, get_decrements, get_virtuals, get_price_sensitive_loads
 export get_availability, get_must_run
@@ -25,6 +27,8 @@ export get_commitment, get_regulation_commitment
 
 include("system.jl")
 include("accessors.jl")
+include("block_inv.jl")
+include("matrices.jl")
 include("deprecated.jl")
 
 end

src/accessors.jl

@@ -49,6 +49,9 @@ get_loads_per_bus(system::System) = system.loads_per_bus
 
 "Returns the power transfer distribution factor of the system."
 get_ptdf(system::System) = system.ptdf
+"Returns the power transfer distribution factor of the system. Calculates it if missing from system"
+retrieve_ptdf(system::System) = coalesce(get_ptdf(system), compute_ptdf(system))
+
 "Returns a `Dictionary` of the line outage distribution factor matrices for the `System` indexed by contingencies."
 get_lodfs(system::System) = system.lodfs
 

src/block_inv.jl

@@ -0,0 +1,86 @@
+# Used to compute PTDF, but it generic code for inverting a large matrix
+# Could/should be open sourced. See:
+# https://github.com/JuliaLinearAlgebra/GenericLinearAlgebra.jl/pull/46
+
+function _block_inv(
+    A::AbstractMatrix,
+    B::AbstractMatrix,
+    C::AbstractMatrix,
+    D_inv::AbstractMatrix,
+)
+    B_D_inv = B * D_inv
+    # Compute -B_D_inv * C + A and store it in A
+    BLAS.gemm!('N', 'N', -1.0, B_D_inv, C, 1.0, A)
+    A = inv(A)
+    B = A * B_D_inv
+    D_inv_C = D_inv * C
+    # Compute -D_inv_C * A and store it in C
+    mul!(C, -D_inv_C, A)
+    # Compute D_inv_C * B + D_inv and store it in D_inv
+    BLAS.gemm!('N', 'N', 1.0, D_inv_C, B, 1.0, D_inv)
+    return A, -B, C, D_inv
+end
+
+@views function _partition_big_mat(mat::AbstractMatrix; block_size::Int=13_000)
+    A = mat[1:block_size, 1:block_size]
+    B = mat[1:block_size, (block_size + 1):end]
+    C = mat[(block_size + 1):end, 1:block_size]
+    D = mat[(block_size + 1):end, (block_size + 1):end]
+    return A, B, C, D
+end
+
+function _blocks_big_mat(
+    mat::T; block_size::Int=13_000
+) where T<:AbstractMatrix{F} where F
+
+    # SubMat is the type that `_partition_big_mat` returns
+    SubMat = SubArray{F, 2, T, Tuple{UnitRange{Int}, UnitRange{Int}}, false}
+    mat_blocks = Tuple{SubMat, SubMat, SubMat, SubMat}[]
+    D = mat
+
+    while true
+        A, B, C, D = _partition_big_mat(D; block_size=block_size)
+        pushfirst!(mat_blocks, (A, B, C, D))
+        size(D, 1) <= block_size && break
+    end
+
+    return mat_blocks
+end
+
+"""
+    big_mat_inv(mat::AbstractMatrix; block_size::Int=13_000) -> AbstractMatrix
+
+Receives a matrix that is supposed to be inverted. If the size of the matrix is larger than
+the defined `block_size`, it first partitions the matrix into smaller blocks until the
+matrices that are supposed to be inverted have size less than `block_size`.
+The partitioned matrix would look like: `mat = [A B; C D]` where the size of A is guaranteed
+to be smaller than the `block_size`. If matrix D is larger than `block_size`, it
+gets partitioned `D = [A1 B1;C1 D1]` and this process continues until all Ais and Dis are
+smaller than `block_size`.
+
+The default `block_size` is set to be `13_000` as we have empirically observed that, for
+matrices smaller than this size, the built-in `inv` can efficiently handle the inversion.
+This was set when doing the calculation of admittance matrix inverse in MISO and depending
+on the application, this number can be adjusted.
+
+
+Staring from the right bottom corner of the partitioned matrix, we use block inversion
+matrix lemma (https://en.wikipedia.org/wiki/Block_matrix) iteratively until the full matrix
+inversion is computed.
+"""
+function big_mat_inv(mat::AbstractMatrix; block_size::Int=13_000)
+    # If the matrix is smaller than the specified block size, just do regular inversion
+    size(mat, 1) <= block_size && return inv(mat)
+    # partition the matrix into smaller blocks.
+    blocks = _blocks_big_mat(mat, block_size=block_size)
+    # iteratively calculating the matrix inversion of each block
+    A, B, C, D = popfirst!(blocks)
+    A, B, C, D = _block_inv(A, B, C, inv(D))
+    num_blocks = length(blocks)
+    for bl_ in 1:num_blocks
+        inverted_mat = [A B; C D]
+        A, B, C, D = popfirst!(blocks)
+        A, B, C, D = _block_inv(A, B, C, inverted_mat)
+    end
+    return [A B; C D]
+end

src/matrices.jl

@@ -0,0 +1,211 @@
+"""
+    compute_ptdf(system::System; block_size, reference_bus_index) -> KeyedArray
+    compute_ptdf(buses::Buses, branches::Branches; block_size, reference_bus_index) -> KeyedArray
+
+Takes a system, or data for that system, representing a `M` branch, `N` bus grid
+and returns the `M * N` DC-Power Transfer Distribution Factor (DC-PTDF) matrix of the network.
+
+For a ~15,000 bus system with aggregated borders, this is expected to take ~1 minute.
+
+# Keywords
+- `block_size=13_000`: Block size to be used when partitioning a big matrix for inversion.
+- `reference_bus=first(keys(buses))`: The name of the reference bus.
+
+# Output
+- `::KeyedArray`: The PTDF matrix; the axes contain the branch and bus names.
+
+!!! note
+    The input data must have no isolated components or islands.
+"""
+function compute_ptdf(system::System; kwargs...)
+    return compute_ptdf(get_buses(system), get_branches(system); kwargs...)
+end
+
+function compute_ptdf(
+    buses::Buses,
+    branches::Branches;
+    block_size=13_000,
+    reference_bus=nothing,
+)
+    bus_names = collect(keys(buses))
+    reference_bus_index = _reference_bus(reference_bus, bus_names)
+
+    incid_matrix = _incidence(buses, branches)
+    n_branches, n_buses = size(incid_matrix)
+
+    # Remove column related to reference bus from incidence matrix
+    incid_matrix = incid_matrix[:, Not(reference_bus_index)]
+
+    B_fl_tilde = sparse(diagm(_series_susceptance(branches))) * incid_matrix
+    B_bus_tilde_inv = big_mat_inv(
+        Matrix(incid_matrix' * B_fl_tilde),
+        block_size=block_size
+    )
+    ptdf_matrix = B_fl_tilde * B_bus_tilde_inv
+
+    # Add reference bus column back, filled with zeros
+    @views ptdf_matrix = hcat(
+        ptdf_matrix[:, 1:(reference_bus_index - 1)],
+        zeros(n_branches),
+        ptdf_matrix[:, reference_bus_index:end],
+    )
+
+    return KeyedArray(ptdf_matrix, (collect(keys(branches)), bus_names))
+end
+
+function _reference_bus(reference_bus, bus_names)
+    reference_bus === nothing && return 1
+
+    idx = findfirst(==(reference_bus), bus_names)
+    idx === nothing && throw(ArgumentError("Reference bus '$reference_bus' not found."))
+    return idx
+end
+
+"""
+    _series_susceptance(branches) -> Vector{Float64}
+
+Calculates the susceptance of the elements in the branch Dictionary The calculation is
+different depending if the element is a line (no tap) or transformer (tap present).
+"""
+function _series_susceptance(branches)
+    susceptance = map(_branch_susceptance, branches)
+    return collect(susceptance)
+end
+
+function _branch_susceptance(b)::Float64
+    if b.tap === missing
+        return -1 / b.reactance
+    end
+
+    return imag(1 / ((b.resistance + b.reactance * 1im) * (b.tap * exp(b.angle * 1im))))
+end
+
+"""
+    _incidence(buses, branches) -> SparseMatrix
+
+Returns the sparse edge-node incidence matrix related to the buses and branches used as
+inputs. Matrix axes correspond to `(keys(branches), keys(buses))`
+"""
+function _incidence(buses, branches)
+    n_buses = length(buses)
+    n_branches = length(branches)
+
+    # Define the mapping of buses/branches to the incidence/PTDF matrix
+    bus_lookup = _make_ax_ref(buses)
+
+    # Compute incidence matrix
+    A_to = sparse(
+        1:n_branches,
+        [bus_lookup[b.to_bus] for b in branches],
+        fill(-1, n_branches),
+        n_branches,
+        n_buses
+    )
+    A_from = sparse(
+        1:n_branches,
+        [bus_lookup[b.from_bus] for b in branches],
+        fill(1, n_branches),
+        n_branches,
+        n_buses
+    )
+    incid_matrix = A_to + A_from
+
+    return incid_matrix
+end
+
+function _make_ax_ref(ax::Dictionary)
+    return Dictionary(keys(ax), 1:length(ax))
+end
+
+"""
+    compute_lodf(system, branch_names_out) -> KeyedArray
+    compute_lodf(system::System, ptdf_matrix, branch_names_out) -> KeyedArray
+    compute_lodf(buses, branches, ptdf, branch_names_out) -> KeyedArray
+
+Returns the `M*O` DC-Line Outage Distribution Factor (DC-LODF) matrix of the network.
+
+**Important Note:** In the current implementation, we use `lodf` only if the contingency
+scenario does not have any line coming in service. We can also use this function if we want
+to ignore the lines coming in service.
+
+# Inputs
+- `buses::Buses`
+- `branches::Branches`
+- `ptdf_matrix`: The pre-calculated PTDF matrix of the system
+- `branch_names_out`: The names of the branches that are going out in the contingency scenario.
+
+# Output
+- The LODF matrix as a `KeyedArray`. The axes are the branch names and `branch_names_out`.
+
+!!! note
+    The resulting LODF matrix is sensitive to the input PTDF matrix. Using a thresholded
+    PTDF as input might lead to imprecisions in constrast to using the full PTDF.
+"""
+function compute_lodf(system::System, branch_names_out)
+    ptdf_matrix = get_ptdf(system)
+    ismissing(ptdf_matrix) && throw(ArgumentError("System PTDF is missing."))
+
+    return compute_lodf(system, ptdf_matrix, branch_names_out)
+end
+
+function compute_lodf(system::System, ptdf_matrix, branch_names_out)
+    buses = get_buses(system)
+    branches = get_branches(system)
+
+    return compute_lodf(buses, branches, ptdf_matrix, branch_names_out)
+end
+
+function compute_lodf(buses::Buses, branches::Branches, ptdf_matrix, branch_names_out)
+    branch_out_names = collect(filter(in(branch_names_out), keys(branches)))
+    branches_out = getindices(branches, branch_out_names)
+
+    if length(branch_out_names) < length(unique(branch_names_out))
+        @debug("Some of the lines to go out were not found in the line data.")
+    end
+
+    if isempty(branches_out)
+        @debug(
+            "All the lines to go out are already out of service.
+            You can ignore this contingency."
+        )
+        return KeyedArray(Matrix{Float64}(undef, 0, 0), (String[], Int[]))
+    end
+
+    incid_out = _incidence(buses, branches_out)
+
+    branch_names = collect(keys(branches))
+    branch_lookup = _make_ax_ref(branches)
+
+    # Our monitored lines are all the lines
+    ptdf_mo = ptdf_matrix.data * incid_out'
+    # Indices of the branches going out
+    ind_br_out = [branch_lookup[b] for b in branch_out_names]
+    ptdf_oo = ptdf_mo[ind_br_out, :]
+    lodf_matrix = ptdf_mo * inv(I - ptdf_oo)
+    # Discard any name that wasn't matched, and ensure the order is in line with the PSSE
+    lodf_matrix = KeyedArray(lodf_matrix, (branch_names, branch_out_names))
+
+    # If a monitored line is going out, manually correct LODF values so that the
+    # post-contingency flow is zero
+    for br in branch_out_names
+        if br in branch_names
+            _correct_lodf!(lodf_matrix, br)
+        end
+    end
+
+    return lodf_matrix
+end
+
+"""
+    _correct_lodf!(lodf_matrix::KeyedArray, br)
+
+Sets the LODF row corresponding to branch `br` to zero, except for the element `(br, br)`,
+which is set to -1. This is to ensure the post-contingency flow on a line that is going out
+and is also monitored is set to zero.
+"""
+function _correct_lodf!(lodf_matrix::KeyedArray, br)
+    lodf_matrix(br, :) .= zeros(size(lodf_matrix(br, :)))
+    lodf_matrix[Key(br), Key(br)] = -1.0
+
+    return lodf_matrix
+end

test/block_inv.jl

@@ -0,0 +1,25 @@
+# This file contains tests related to the matrix block inversion procedure
+
+@testset "Block matrix inversion" begin
+    big_mat_inv = FullNetworkSystems.big_mat_inv
+    @testset "fallback to `inv`" begin
+        n = 1000
+        M = randn(n, n)
+        # default block size should be >1000, so should fallback to `inv` here
+        @test big_mat_inv(M) == inv(M)
+    end
+
+    @testset "use block algorithm" begin
+        n = 1000
+        for _ in 1:3
+            M = randn(n, n)
+            @test inv(M) ≈ big_mat_inv(M; block_size=500) rtol=1e-3
+        end
+
+        n = 2000
+        for _ in 1:3
+            M = randn(n, n)
+            @test inv(M) ≈ big_mat_inv(M; block_size=1800) rtol=1e-3
+        end
+    end
+end