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

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…

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

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

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

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.

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
``````

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!

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