Obtaining incidence + PTDF matrix from network_data

[junglegobs]

This is more of a question rather than an actual issue (though the lack of documentation could be considered an issue).

I'm using PowerModels to read data from matpower files, and I want to use it to extract the incidence and PTDF matrices. I assumed that there would be inbuilt functions for these, and I just wanted to make sure that this is not the case before I go on to write my own functions to do this. If people think this would be useful I could even write this into PowerModels.jl itself?

[ccoffrin]

There is not a pre-built function for the incidence matrix, but I will consider adding it. It would be along this line (given a PM data dict),

    I = Int64[]
    J = Int64[]
    V = Float64[]

    for (i,branch) in data["branch"]
        f_bus = branch["f_bus"]
        t_bus = branch["t_bus"]
        if branch["br_status"] != 0 && <check that buses are active (i.e. not status 4)>
            f_bus = bus_to_idx[f_bus]
            t_bus = bus_to_idx[t_bus]
            push!(I, f_bus); push!(J, i); push!(V,  1)
            push!(I, t_bus); push!(J, i); push!(V,  1)
        end
    end

    M = sparse(I,J,V)

There is no function to compute the PTDF matrix explicitly, but it can be derived from the Admittance matrix and its inverse, which functions do exist for. More details here, https://lanl-ansi.github.io/PowerModels.jl/stable/utilities/#Lazy-Line-Flow-Limits-1

Two additional remarks,

• I did recently learn that PM Admittance matircies are the complex conjugate of the standard form. This will be corrected in a future release.
• Computing the PTDF for large networks (say 10,000+ buses) can be very slow becouse the inverse is a dense matrix.

[junglegobs]

Thanks! Regarding the large PTDF calculation, that's good to know but not a problem in my case at least.

[junglegobs]

If I could possibly bother you with a question: I'm new to the Matpower format and PowerModels.jl, and I'm confused as to the difference between keys of data["branch"] and the "index" of a branch. As far as I can tell these are and must be the same, except that the keys are Strings and not Ints? I ask since the above pseudo code does not work since i is a string (for the example case3.m).

[ccoffrin]

Sure.

As far as I can tell these are and must be the same, except that the keys are Strings and not Ints?

That is correct. The key of every component type in PowerModels are Ints; these are found in comp["index"]. The reason that the keys are strings in the data Dict is so that the data can be serialized to a JSON file, which does not support integer keys.

Here is a more precise example of building an incidence matrix,

using SparseArrays

I = Int64[]
J = Int64[]
V = Float64[]

for (i,branch) in data["branch"]
    idx = branch["index"]
    f_bus = branch["f_bus"]
    t_bus = branch["t_bus"]
    if branch["br_status"] != 0 # && <check that buses are active (i.e. not status 4)>
        push!(I, f_bus); push!(J, idx); push!(V,  1)
        push!(I, t_bus); push!(J, idx); push!(V,  1)
    end
end

A = sparse(I,J,V)

I'll also note that the Julia forums here, https://discourse.julialang.org/, provide the best form for asking general questions about PowerModels.

[YutingMou]

I coded a small script for calculating PTDF and get the results very close to that of MatPower.

using PowerModels
using DataFrames

network_data = PowerModels.parse_file("matpower/case14.m")  #  load the case
bus_reference = 1 ## could adapt by checking bus["bus_type"] == 3
bus_data = network_data["bus"]
bus_data = sort(bus_data, by=key_str->parse(Int, key_str))
branch_data = network_data["branch"]
branch_data = sort(branch_data, by=key_str->parse(Int, key_str))


B =  zeros(length(bus_data), length(bus_data))
for (key ,branch) in branch_data
    bus_fr = branch["f_bus"]
    bus_to = branch["t_bus"]
    B[bus_fr, bus_to] = B[bus_fr, bus_to] - 1/branch["br_x"]
    B[bus_to, bus_fr] = B[bus_to, bus_fr] - 1/branch["br_x"]
    B[bus_fr, bus_fr] = B[bus_fr, bus_fr] + 1/branch["br_x"]
    B[bus_to, bus_to] = B[bus_to, bus_to] + 1/branch["br_x"]
end
B_aux = B
# zero row and col for refbus, then put 1 in diag so we can invert it
B_aux[:,bus_reference] = zeros(length(bus_data),1)
B_aux[bus_reference,:] = zeros(1,length(bus_data))
B_aux[bus_reference,bus_reference] = 1.0
B_inv = inv(B_aux)
B_inv[bus_reference, bus_reference]=0;

PTDF = zeros(length(branch_data), length(bus_data))
for (key ,branch) in branch_data
    bus_fr = branch["f_bus"]
    bus_to = branch["t_bus"]
    for bus_id in 1:length(bus_data) # will break if buses are not numbered 1-to-n
        PTDF[branch["index"],bus_id] = 1/branch["br_x"] *(
            (B_inv[bus_fr, bus_id] - B_inv[bus_fr, bus_reference])
            - (B_inv[bus_to, bus_id] - B_inv[bus_to, bus_reference])
            )
    end
end

Reference: http://www.ee.umn.edu/class/ee8725/HW%204%20ANSWER%20KEY.pdf

[junglegobs]

Thanks @YutingMou 👍

[ccoffrin]

@junglegobs and @YutingMou, the basic data utilities have been added to address these points.

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

I hope you find the new features helpful!