Where is the "rate_a" in PowerModels.jl from

The data provided by PowerModels.jl is inconsistent with my local raw Matpower file.
The latter is all zero (seems to mean that there is no limit), but the former has a vector of numbers that I can’t tell where they come from. (Moreover, what is its unit (in MVA? or in p.u. (which is unlikely)))

The code is

import PowerModels
R = let
    PowerModels.silence()
    t = PowerModels.make_basic_network(PowerModels.parse_file("data/case118.m"))
    A = PowerModels.calc_basic_incidence_matrix(t)
    B = size(A, 1)
    [t["branch"]["$b"]["rate_a"] for b in 1:B]
end

make_basic_network calls calc_thermal_limits!:

If you just want to read the raw data, you probably want to use parse_file instead of make_basic_network:

t = PowerModels.parse_file("data/case118.m")

So, it is not an honest rate_a in physical sense, but rather an inference conditional on other related data (e.g. voltage, angles).
(will this calculated data be useful in practice? I wonder.)

See Basic Data Utilities · PowerModels, in particular “all branches have explicit thermal limits” and “users requiring any of the features listed above for their analysis should use the non-basic PowerModels routines.”

PowerModels is designed around some very specific formulations. They’re not for everyone.

You are correct, in that the derived rate_a will never be a binding constraint in this model, but it is a valid upper bound on the amount of apparent power that can flow across the line. That still provides some useful information I think.

@oscar gives good advise. make_basic_network was intended primarily as a tool for teaching. If you don’t want the properties that it enforces then work with the data output by parse_file.

To answer your original question, “Where is the “rate_a” in PowerModels.jl from” if you do not call,
PowerModels.silence() then you should see some logging messages about how rate_a values were missing so they were replaced with different ones.

Side note - Are you familiar with PGLib-OPF and PGLib.jl? These datasets have been cleaned up, so you should not see as many warnings during data reading process.

2 Likes

I’ll try to learn them, thanks!

Update: I’ve worked on this network pglib_opf_case118_ieee.m

PowerFlowWithJuMP
import SparseArrays, PGLib; RawDict = PGLib.pglib("pglib_opf_case118_ieee.m");
import JuMP, Gurobi; GRB_ENV = Gurobi.Env();
import LinearAlgebra.dot as dot

function rectify_dir!(bft)
    for r in 1:size(bft, 1)
        if bft[r, 1] > bft[r, 2]
            bft[r, 1], bft[r, 2] = bft[r, 2], bft[r, 1]
        elseif bft[r, 1] == bft[r, 2]
            error("self-branch")
        end
    end
end; function assert_is_connected(bft)
    B = size(bft, 1)
    pb = Vector(1:B) # primal bs
    sn = [1] # subnet
    for ite in 1:B # only du a full iteration as a conservative choice
        progress = false
        for (i, b) in enumerate(pb) # take out branch b
            if bft[b, 1] in sn
                on = bft[b, 2] # the other node
                on ∉ sn && push!(sn, on)
            elseif bft[b, 2] in sn
                on = bft[b, 1] # the other node
                on ∉ sn && push!(sn, on)
            else
                continue
            end
            popat!(pb, i); progress = true; break # fathom branch b
        end
        if progress == false
            error("ite = $ite. All the rest branches are not connected to the subnet being investigated. Check the current subnet")
        end
        1:maximum(bft) ⊆ sn && return  # The graph is proved to be connected
    end
    error("here shouldn't be reached")
end; function bft_2_A(bft)
    B, N = size(bft, 1), maximum(bft)
    return SparseArrays.sparse([Vector(1:B); Vector(1:B)], vec(bft), [-ones(Int, B); ones(Int, B)], B, N)
end; function adjnode(n)
    nv = Int[]
    for ci in findall(x -> x == n, bft)
        b, c = ci.I # row(= branch), col
        on = bft[b, 3 - c]
        on ∉ nv && push!(nv, on)
    end
    sort(nv)
end; function get_b(n, m) # 🟠 use this after `rectify_dir!(bft)` had been done
    n > m && ((n, m) = (m, n))
    for r in 1:size(bft, 1)
        (bft[r, 1] == n && bft[r, 2] == m) && return r
    end
    error("no such branch exist")
end;
function bpf(b, Dvec = Dref, pvec = pg) # dispatch function of power flow on branch b
    b == 9 && return -pvec[findfirst(x -> x == 10, Gnode)::Int] # b8[9] = -3
    b == 7 && return bpf(9) # b8[7] = -1
    b == 94 && return bpf(93) # b8[94] = -1
    b == 127 && return -bpf(126) # b8[127] = -1
    b == 134 && return -pvec[findfirst(x -> x == 87, Gnode)::Int] # b8[134] = -3
    b == 176 && return -pvec[findfirst(x -> x == 111, Gnode)::Int] # b8[176] = -3
    b == 184 && return Dvec[findfirst(x -> x == 117, Dnode)::Int] # b8[184] = -2
    dup_bvec = [67, 76, 99, 86, 124, 139, 142]
    b ∈ dup_bvec && return false # b8[dup_bvec] .= 0; b8[dup_bvec .- 1] .= 2
    f = findfirst(x -> x == b, F2B)
    return pf[f] # the JuMP variable
end
B, N = 186, 118; # fixed for this network
G = 54; # ✅ in case118, there is no such case that 2 generators connecting to the same bus
D = 99; # demands are data
bft = let # load a rectified bft
    bft = [RawDict["branch"]["$r"][c == 1 ? "f_bus" : "t_bus"] for r in 1:length(RawDict["branch"]), c in 1:2]; # raw bft
    assert_is_connected(bft);
    rectify_dir!(bft);
    bft
end; A = bft_2_A(bft);
Dnode = [RawDict["load"]["$d"]["load_bus"] for d in 1:D]; # function n2d(n) return findfirst(x -> x == n, Dnode)::Int end;
Dref  = [RawDict["load"]["$d"]["pd"]       for d in 1:D]; # demand vector (std reference)
Gnode = [ RawDict["gen"]["$g"]["gen_bus"]  for g in 1:G]; # function n2g(n) return findfirst(x -> x == n, Gnode)::Int end;
Cgen = let
    Cgen = zeros(G)
    for g in 1:G
        v = RawDict["gen"]["$g"]["cost"]
        isempty(v) || (Cgen[g] = v[1] / 1000)
    end
    Cgen
end;
n8 = let # elements being
    # [1]: a bus with load only, and whose degree is larger than 2
    # [2]: a bus with gen only, and whose degree is larger than 2
    # [3]: a bus with both power injection and load
    # [-1]: a bus with load only, and whose degree is 1, we won't add KCL constr here
    # [-3]: a bus with gen only, and whose degree is 1, we won't add KCL constr here
    # [-2]: a bare bus whose degree is 2, we won't add KCL constr here
    # [0]: a bare joint bus
    n8 = zeros(Int8, N)
    n8[Dnode] .+= 1
    n8[Gnode] .+= 2
    n8[117] = -1
    n8[[10, 87, 111]] .= -3
    n8[[9, 63, 81]] .= -2
    n8
end; b8 = let # Note this procedure is specific to case118, therefore executing once is sufficient
    # elements being
    # [2]: this line represents a 2-parallel branch, works as a normal line
    # [1]: a normal line
    # [0]: this line was banned (pf ≡ 0) due to parallel redundancy, we won't allocate pf variable here
    # [-1]: power flow of this line hinges on one of its adjacent line, we won't allocate pf variable here
    # [-2]: power flow of this line hinges on the corresponding pure load, we won't allocate pf variable here
    # [-3]: power flow of this line hinges on the corresponding pure gen, we won't allocate pf variable here
    b8 = ones(Int8, size(bft, 1))
    b8[9] = -3
    b8[7] = -1
    b8[94] = -1
    b8[127] = -1
    b8[134] = -3
    b8[176] = -3
    b8[184] = -2
    dup_bvec = [67, 76, 99, 86, 124, 139, 142]
    b8[dup_bvec] .= 0; b8[dup_bvec .- 1] .= 2
    b8
end; F2B = findall(s -> s ∈ 1:2, b8); F = length(F2B);
model = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)); JuMP.@variable(model, pg[g = 1:G] >= 0);
JuMP.@variable(model, pf[f = 1:F]);
let # KCL constrs
    for n in findall(s -> s == 0, n8)
        bv, sgv = SparseArrays.findnz(view(A, :, n))
        JuMP.@constraint(model, dot(sgv, bpf.(bv)) == 0)
    end
    for n in findall(s -> s == 1, n8)
        bv, sgv = SparseArrays.findnz(view(A, :, n))
        JuMP.@constraint(model, dot(sgv, bpf.(bv)) == Dref[findfirst(x -> x == n, Dnode)::Int])
    end
    for n in findall(s -> s == 3, n8)
        bv, sgv = SparseArrays.findnz(view(A, :, n))
        JuMP.@constraint(model, dot(sgv, bpf.(bv)) + pg[findfirst(x -> x == n, Gnode)::Int] == Dref[findfirst(x -> x == n, Dnode)::Int])
    end
    for n in findall(s -> s == 2, n8)
        bv, sgv = SparseArrays.findnz(view(A, :, n))
        JuMP.@constraint(model, dot(sgv, bpf.(bv)) + pg[findfirst(x -> x == n, Gnode)::Int] == 0)
    end
end
JuMP.@objective(model, Min, dot(Cgen, pg))
JuMP.optimize!(model); JuMP.assert_is_solved_and_feasible(model; allow_local = false);
pg_val = JuMP.value.(pg) # output power of all generators
bpf_val = @. JuMP.value(bpf(1:B)) # power flow on all branches

  • /api/*.m - heavily loaded test cases (i.e. binding thermal limit constraints)
  • /sad/*.m - small phase angle difference cases (i.e. binding phase angle difference constraints)

I can understand the first thermal limit, but not the other one.
Why is small phase angle difference considered “binding”? (I thought large is binding?)
What is meant by “phase angle difference constraint”?

What does the “removing” in these logging mean?
I can still access them normally.

julia> import PGLib

julia> d = PGLib.pglib("pglib_opf_case118_ieee.m");
[info | PowerModels]: removing 3 cost terms from generator 32: Float64[]
# omit many lines
[info | PowerModels]: removing 1 cost terms from generator 30: [2575.8442, 0.0]
[info | PowerModels]: removing 3 cost terms from generator 3: Float64[]

julia> d["gen"]["30"]["cost"]
2-element Vector{Float64}:
 2575.8442
    0.0