Calculating power transfer distribution factor (PTDF) matrix using PowerModels.jl

Specific Domains Optimization (Mathematical)
1

I want to obtain a PTDF matrix to use in my power system optimisation model using data in PowerModels.jl format (I also raised an issue for this to be integrated into PowerModels.jl itself: https://github.com/lanl-ansi/PowerModels.jl/issues/728).

I have been trying to do this following the method described here in Appendix B2, but I’m getting confused with what the output of PowerModels.calc_susceptance_matrix is - is it the branch susceptance matrix? The nodal susceptance matrix? How do I go from one to the other?

I’m tempted to say that all I need to do is the following:

    network_data = PowerModels.parse_file("network.m")
    B = PowerModels.calc_susceptance_matrix(network_data).matrix
    Binv = PowerModels.calc_susceptance_matrix_inv(network_data).matrix
    PTDF = B*Binv

But I assume this is wrong…

2

I did not have time to check this in detail, so assume there are bugs in there, but I think what you need is something like this,

network_data = PowerModels.parse_file("case30Q.m")
B_inv = PowerModels.calc_susceptance_matrix_inv(network_data).matrix

PTDF = zeros(length(network_data["bus"]), length(network_data["branch"]))
for (b,branch) in network_data["branch"]
    bus_fr = branch["f_bus"]
    bus_to = branch["t_bus"]
    b_val = imag(inv(branch["br_r"] + branch["br_x"]im))
    for n in 1:length(network_data["bus"]) # will break if buses are not numbered 1-to-n
        PTDF[n,branch["index"]] = B_inv[bus_fr, n] - B_inv[bus_to, n]
    end
end

It is important to note that PowerModels supports bus and branch ids that are not 1-to-n and components that can be omitted via status values. PTDF and similar matrices require, all components to be active and the ids to be numbered 1-to-n.

1 Like
3

@jdlara-berkeley, maybe you can provide further details here?

4

Thanks for this, I usually use AxisArrays which might resolve the problem of ids

5

Oh nice. Can you post an example of how you set these up? I tried using them some time ago for a similar goal but was not able to get them working.

6

Sure! Here’s an adapted version of my code for the PTDF calculation:

using PowerModels
using AxisArrays

networkData = PowerModels.parse_file("/home/u0128861/.julia/packages/PowerModels/fEPoB/test/data/matpower/case30.m")

# Create vector of bus and branch names
# If no name / string given, use the index instead
BusNames = sort(
    [get(v, "string", v["index"]) for (k,v) in networkData["bus"]]
)
BranchNames = sort(
    [get(v, "name", v["index"]) for (k,v) in networkData["branch"]]
)
B_inv = PowerModels.calc_susceptance_matrix_inv(networkData)
PTDF = AxisArray(
    zeros(
        length(networkData["bus"]), length(networkData["branch"])
    ),
    Axis{:Bus}(BusNames),
    Axis{:Branch}(BranchNames),
)

# Dictionaries mapping bus / branch indices to their names
busIndices2Names = Dict(
    v["index"] => get(v, "name", v["index"]) for (k,v) in networkData["bus"]
)
branchIndices2Names = Dict(
    v["index"] => get(v, "name", v["index"]) for (k,v) in networkData["branch"]
)

# Compute PTDF matrix
for (b, branch) in networkData["branch"]
    bus_fr = B_inv.bus_to_idx[branch["f_bus"]]
    bus_to = B_inv.bus_to_idx[branch["t_bus"]]
    b_val = imag(inv(branch["br_r"] + branch["br_x"]im))
    branch_name = branchIndices2Names[branch["index"]]
    for (k, bus) in networkData["bus"]
        n = B_inv_bus_to_idx[bus["index"]]
        bus_name = busIndices2Names[bus["index"]]
        PTDF[bus_name, branch_name] = B_inv.matrix[bus_fr, n] - B_inv.matrix[bus_to, n]
    end
end

I noticed you had defined a variable b_val, but this is not used in the PTDF calculation you wrote - where does it go? (As you may have noticed by now, I am not an electrical engineer and I’m quite clueless about power flow…)

For that particular test case, the buses / branches didn’t have names, hence the get(v, "name", v["index"]) calls. I should probably also include checks that a) all buses / branches have names and if not b) that their indices go from 1:n.

Regarding point b), that because if you construct an AxisArray with an axis of integers, it will behave like a normal Array unless you pass the commant atvalue when getting the index:

julia> A = AxisArray([1,2,3], Axis{:String}(["1","2","4"]))
1-dimensional AxisArray{Int64,1,...} with axes:
    :String, ["1", "2", "4"]
And data, a 3-element Array{Int64,1}:
 1
 2
 3

julia> A["4"]
3

julia> A = AxisArray([1,2,3], Axis{:Int}([1,2,4]))
1-dimensional AxisArray{Int64,1,...} with axes:
    :Int, [1, 2, 4]
And data, a 3-element Array{Int64,1}:
 1
 2
 3

julia> A[4]
ERROR: BoundsError: attempt to access 3-element Array{Int64,1} at index [4]
Stacktrace:
 [1] getindex at ./array.jl:744 [inlined]
 [2] getindex(::AxisArray{Int64,1,Array{Int64,1},Tuple{Axis{:Int,Array{Int64,1}}}}, ::Int64) at /home/u0128861/.julia/packages/AxisArrays/IFpjG/src/indexing.jl:39
 [3] top-level scope at none:0

julia> A[atvalue(4)]
3
7

You are right it should probably be something like,

PTDF[bus_name, branch_name] = b_val*(B_inv.matrix[bus_fr, n] - B_inv.matrix[bus_to, n])

But I did not check this in detail. My goal was to show you just roughly how it was done. My recommendation would be to take the precise definition out of a textbook to make sure.

Thanks for the tip about AxisArrays, I had missed the atvalue function!

8

A new basic data feature has been added to PowerModels in v0.17.4 to help support this kind of analysis. For details see,

https://lanl-ansi.github.io/PowerModels.jl/stable/basic-data-utilities/

1 Like