How to map read_variable(res, "FlowActivePowerVariable__Line") to branch data?

[yasirroni]

I use read_variable(res, "FlowActivePowerVariable__Line") to read the results of branch power flow. I use PowerModels to calculate branch flow.

In my use case, the branch name is not in order, something like this:

lines
Dict{Int64, Line} with ... entries:
  5  => Line(5, ...)
  16 => Line(16,...)
  12 => Line(12,...)

Also, some lines are active, some are not. Some lines even has identical f_bus and t_bus.

Using read_variable(res, "FlowActivePowerVariable__Line"), I seem to get unordered but monotonically increasing integer number, not at all similar with the my branch names.

Thus, what is the best idea to "map" or rename the column names of DataFrame returned by read_variable(res, "FlowActivePowerVariable__Line")? Also, since the number of line is reduced, does it merge identical arc or unintentionally deletes one of them?

[jd-lara]

There are a lot of questions about what you are trying to do here. 1. Why are you calculating the Line flows with PowerModels? Which version of PowerSimulations are you using?

Are you solving a Decision model or a simulation? I don't understand how are you getting those results when we call get_variable the return is a Dataframe not a Dict, why does the order matter? . The code you pasted doesn't come from us.

I don't know what kind of system you are running if "Some lines even has identical f_bus and t_bus." that seems to be a data issue.

Please provide more details.

[yasirroni]

Hi @jd-lara, thanks for the reply.

Why are you calculating the Line flows with PowerModels?

Because Sienna didn't solve power flow in PowerSimulations. From the doc, PowerModels power flow is supported.

set_network_model!(
    template_uc,
    NetworkModel(
        NFAPowerModel,  # or DCMPPowerModel,
    ),
)

Which version of PowerSimulations are you using?

version = "0.30.1"

Are you solving a Decision model or a simulation?

A decision model.

problem = DecisionModel(template_uc, sys; optimizer=solver, horizon=hours)
build!(problem; output_dir=mktempdir())
solve!(problem)

I don't understand how are you getting those results when we call get_variable the return is a Dataframe not a Dict, why does the order matter?

Did you mean read_variable(res, "FlowActivePowerVariable__Line")? If yes, I don't think order is matter. The problem is that the column name. You can see that the column name if the DataFrame didn't match the PSY.Line name. Thus, I'm a bit confused. What is the proper way to map the columns of read_variable(res, "FlowActivePowerVariable__Line") to their respective Line name in get_name(line)? currently, I can guess it just by ordering first, but I don't know if this could introduce bug later. In my case, the Line names is not monotonically increasing (for example, not every integer exist: ["1", "2", "5", "10"]).

if "Some lines even has identical f_bus and t_bus." that seems to be a data issue.

This is what other engineer decided when he model transmission expansion planning (and providing data/result). Yes, I can merge multiple lines into one, but require pre-computation. And both PowerModels and MATPOWER support this kind of data representation anyway.

[jd-lara]

Sienna uses the same model as PowerModels inside of the DecisionModel, PowerModels uses the arcs to map the Line names and we keep that reporting mechanism. We need to update the code for the functions add_pm_variable_refs! to reflect the line name in the results instead of the arcs, if you volunteer to update that code we can improve that reporting as suggested.

[jd-lara]

about this: "This is what other engineer decided when he model transmission expansion planning (and providing data/result). Yes, I can merge multiple lines into one, but require pre-computation. And both PowerModels and MATPOWER support this kind of data representation anyway." -> How is this acceptable? Are you asking us to perform a correction to incorrect data? What's the expected outcome for this situation?

[yasirroni]

Hi @jd-lara,

I'm happy to make a PR around that. I think using Arc is not a good idea to pass Line to PowerModels, as I mentioned earlier that PowerModels support multiple parallel lines. In Sienna PowerSystems, it is also supported in data creation without any warning. Thus, maybe,

for line in values(lines)
    println(get_name(get_arc(line)))
end

5 -> 4
6 -> 7
12 -> 11
6 -> 7

Thus, using Line name is a good way to go.

for line in values(lines)
    println(get_name(line))
end

5
16
12
8

As far as I know, PowerModels store the line data in three id format as shown in my comments lanl-ansi/PowerModels.jl#977 (comment), that is (branch_id, f_bus, t_bus).

I might need to read some code to contribute on this add_pm_variable_refs!. Notes that I think it will took me sometime to make this PR.

Can you give me some direction on the way to start? Thank you.

[yasirroni]

How is this acceptable? Are you asking us to perform a correction to incorrect data? What's the expected outcome for this situation?

Just want PowerModels support, that is keeping the Line id separates for data with parallel lines.

See this example for proof of support of parallel lines:

pm_data = PowerModels.parse_file("../../env/lib/python3.13/site-packages/matpower/data/case9.m")
# add new parallel branch with identical br_x and br_r, but lower rate_a
branch_id = "7"
n_branch = length(network["branch"])
branch_id_new = string(n_branch + 1)
network["branch"][branch_id_new] = copy(network["branch"][branch_id])

pm = instantiate_model(pm_data, DCPPowerModel, PowerModels.build_opf)
result = optimize_model!(pm, optimizer=Ipopt.Optimizer)

julia> vars = var(pm)
Dict{Symbol, Any} with 5 entries:
  :p         => Dict{Any, Any}((9, 4, 9)=>-0_p[(9, 9, 4)], (7, 8, 2)=>0_p[(7, 8, 2)…
  :p_dc_cost => Dict{Int64, Any}()
  :pg        => 1-dimensional DenseAxisArray{VariableRef,1,...} with index sets:…
  :p_dc      => 1-dimensional DenseAxisArray{Any,1,...} with index sets:…
  :pg_cost   => Dict{Int64, Any}(2=>850.0000000000001 0_pg[2]² + 120 0_pg[2] + 600,…

julia> vars[:p]
Dict{Any, Any} with 20 entries:
  (9, 4, 9)  => -0_p[(9, 9, 4)]
  (7, 8, 2)  => 0_p[(7, 8, 2)]
  (8, 8, 9)  => 0_p[(8, 8, 9)]
...
  (3, 6, 5)  => -0_p[(3, 5, 6)]
  (10, 8, 2) => 0_p[(10, 8, 2)]
  (8, 9, 8)  => -0_p[(8, 8, 9)]
  (7, 2, 8)  => -0_p[(7, 8, 2)]
  (9, 9, 4)  => 0_p[(9, 9, 4)]
...
  (2, 5, 4)  => -0_p[(2, 4, 5)]

We can see different variable for (x, 8, 2). I've tested they both works on AC, DC, and NFA power flow properly (considering each rate_a, br_r, and br_x).

[jd-lara]

Hi @jd-lara,

I'm happy to make a PR around that. I think using Arc is not a good idea to pass Line to PowerModels, as I mentioned earlier that PowerModels support multiple parallel lines. In Sienna PowerSystems, it is also supported in data creation without any warning. Thus, maybe,

for line in values(lines)
    println(get_name(get_arc(line)))
end

5 -> 4
6 -> 7
12 -> 11
6 -> 7

Thus, using Line name is a good way to go.

for line in values(lines)
    println(get_name(line))
end

5
16
12
8

As far as I know, PowerModels store the line data in three id format as shown in my comments lanl-ansi/PowerModels.jl#977 (comment), that is (branch_id, f_bus, t_bus).

I might need to read some code to contribute on this add_pm_variable_refs!. Notes that I think it will took me sometime to make this PR.

Can you give me some direction on the way to start? Thank you.

its being a while since I looked at the code but the place to start is get_branch_to_pm when we create the PM data model under the hood. If lines are in parallel the id and the arcs are reflected there.

In this function

function get_branch_to_pm(
    ix::Int,
    branch::PSY.ACBranch,
    ::Type{StaticBranchUnbounded},
    ::Type{<:PM.AbstractPowerModel},
)
    PM_branch = Dict{String, Any}(
        "br_r" => PSY.get_r(branch),
        "shift" => 0.0,
        "br_x" => PSY.get_x(branch),
        "g_to" => 0.0,
        "g_fr" => 0.0,
        "b_fr" => PSY.get_b(branch).from,
        "f_bus" => PSY.get_number(PSY.get_arc(branch).from),
        "br_status" => Float64(PSY.get_available(branch)),
        "t_bus" => PSY.get_number(PSY.get_arc(branch).to),
        "b_to" => PSY.get_b(branch).to,
        "index" => ix,
        "angmin" => PSY.get_angle_limits(branch).min,
        "angmax" => PSY.get_angle_limits(branch).max,
        "transformer" => false,
        "tap" => 1.0,
    )
    return PM_branch
end

and in here

function get_branches_to_pm(
    sys::PSY.System,
    network_model::NetworkModel{S},
    ::Type{T},
    branch_template::BranchModelContainer,
    start_idx = 0,
) where {T <: PSY.ACBranch, S <: PM.AbstractPowerModel}
    PM_branches = Dict{String, Any}()
    PMmap_br = Dict{PM_MAP_TUPLE, T}()

    network_reduction = get_network_reduction(network_model)
    removed_branches_names = PNM.get_removed_branches(network_reduction)
    for (d, device_model) in branch_template
        comp_type = get_component_type(device_model)
        if comp_type <: TwoTerminalHVDCTypes
            continue
        end
        !(comp_type <: T) && continue
        start_idx += length(PM_branches)
        for (i, branch) in enumerate(get_available_components(device_model, sys))
            if PSY.get_name(branch) ∈ removed_branches_names
                @debug "Skipping branch $(PSY.get_name(branch)) since it is radial"
                continue
            end
            ix = i + start_idx
            PM_branches["$(ix)"] =
                get_branch_to_pm(ix, branch, get_formulation(device_model), S)
            if PM_branches["$(ix)"]["br_status"] == true
                f = PM_branches["$(ix)"]["f_bus"]
                t = PM_branches["$(ix)"]["t_bus"]
                PMmap_br[(from_to = (ix, f, t), to_from = (ix, t, f))] = branch <- here the index is being written correctly to the map. 
            end
        end
    end
    return PM_branches, PMmap_br
end

you can see we pass the index so it shouldn't get over written.

Possibly there could be a bug in this function that performs the reverse map

function add_pm_variable_refs!(
    container::OptimizationContainer,
    d_class::Type,
    device_types::Vector,
    pm_map::Dict,
    pm_variable_map::Dict,
    pm_variable_types::Base.KeySet,
    time_steps::UnitRange{Int},
)
    for d_type in Set(device_types)
        devices = [d for d in pm_map if typeof(d[2]) == d_type]
        for (pm_v, ps_v) in pm_variable_map[d_class]
            if pm_v in pm_variable_types
                for dir in fieldnames(typeof(ps_v))
                    var_type = getfield(ps_v, dir)
                    var_type === nothing && continue
                    var_container = add_variable_container!(
                        container,
                        var_type,
                        d_type,
                        [PSY.get_name(d[2]) for d in devices],
                        time_steps,
                    )
                    for t in time_steps, (pm_d, d) in devices
                        var = PM.var(container.pm, t, pm_v, getfield(pm_d, dir))
                        var_container[PSY.get_name(d), t] = var
                    end
                end
            end
        end
    end
    return
end