Data.jl

#* Functions for handling the Network Data Dictionary (NDD) of PowerModels.jl
#*------------------------------------------------------------------------------

#=
Builds a network data dictionary conform with PowerModels.jl from several files containing certain input data:
model - whether to calculate AC- or DC-PF (:ac or :dc)
pgfile - .RAW or .m file that contains the detailed grid data (e.g. generator and load information etc.)
busfile - .csv file that contains geographic locations of all buses
branchfile - optional .csv file that contains transmission line lengths
nddfile - .jld2 file in which the final network data dictionary will be saved
=#

"""
    build_network_data(file::String; kwargs...)

Reads a .RAW or .m `file` containing the data of a power grid and calculates the power flow according to the chosen model (see below). Returns a network data dictionary conform with the [PowerModels Network Data Format](https://lanl-ansi.github.io/PowerModels.jl/stable/network-data/) that also contains the calculated power flow solution.
"""
function build_network_data(
        file::String;
        locfile = "", lenfile = "", outfile = "",
        pf_model = :ac,
        kwargs...
    )

    ### Read grid data and calculate operation point
    network_data = calc_pf(file, pf_model=pf_model; kwargs...)
    
    ### Add additional data to network_data, if given
    if isempty(locfile) == false
        add_locs!(network_data, locfile) # add geographic bus locations
    end
    if isempty(lenfile) == false
        add_tl_lengths!(network_data, lenfile)
    end

    add_tl_voltages!(network_data) # add transmission line voltages

    if isempty(outfile) == false
        save(outfile, "network_data", network_data) # save network data
    end

    return network_data
end

#*------------------------------------------------------------------------------

#=
Reads the geographic bus locations (longitude and latitude coordinates in degrees) from a .csv file and adds them to the network data dictionary (NDD). Returns the resulting NDD, in which each bus has a "bus_lat" and "bus_lon" entry. The columns containg the longitude and latitude coordinates in the .csv file should be called "SubLongitude" and "SubLatitude", respectively.
=#
function add_locs!(network_data::Dict{String,<:Any}, csvfile::String)
    pos = Dict{String,Any}(
        "baseMVA"  => 100.0,
        "bus" => Dict{String,Any}(), 
        "per_unit" => true
    )

    BusData::DataFrame = CSV.File(csvfile) |> DataFrame!
    N = size(BusData, 1) # number of buses
    sizehint!(pos["bus"], N)

    for bus in eachrow(BusData)
        pos["bus"][string(bus[:Number])] = Dict(
            "bus_lon" => bus[:SubLongitude], 
            "bus_lat" => bus[:SubLatitude]
        )
    end
    update_data!(network_data, pos) # add positions to network_data
    
    return network_data
end

#*------------------------------------------------------------------------------

#=
Reads transmission line lengths from a .csv file and adds them to the network data dictionary.
=#
function add_tl_lengths!(network_data::Dict{String,<:Any}, csvfile::String)
    len = Dict{String,Any}(
        "baseMVA"  => 100.0,
        "branch" => Dict{String,Any}(), 
        "per_unit" => true
    )

    BranchData::DataFrame = CSV.File(csvfile) |> DataFrame!
    N = size(BranchData, 1) # number of branches
    sizehint!(len["branch"], N)

    for branch in eachrow(BranchData)
        from, to = branch[:BusNumFrom], branch[:BusNumTo]
        indices = findall(
            b -> ((b["f_bus"], b["t_bus"]) == (from, to) 
            || (b["f_bus"], b["t_bus"]) == (to, from)), network_data["branch"]
        )
        for index in indices
            len["branch"][index] = Dict("length" => branch[:GICLineDistanceKm])
        end
    end
    update_data!(network_data, len) # add positions to network_data
    
    return network_data
end

#*------------------------------------------------------------------------------

#=
Determines voltage levels of transmission lines and adds them to the network data dictionary.
=#
function add_tl_voltages!(network_data::Dict{String,<:Any})
    voltages = Dict{String,Any}(
        "baseMVA"  => 100.0,
        "branch" => Dict{String,Any}(), 
        "per_unit" => true
    )
    L = length(network_data["branch"])
    sizehint!(voltages["branch"], L)

    for (i, branch) in network_data["branch"]
        if branch["transformer"] == true
            ### Transformers have a transmission line voltage of 0 kV
            voltages["branch"][i] = Dict(
                "tl_voltage" => 0. 
            )
        else
            ### Transmission lines get base voltage of from bus
            from = branch["f_bus"]
            voltages["branch"][i] = Dict(
                "tl_voltage" => network_data["bus"][string(from)]["base_kv"]
            )
        end
    end
    update_data!(network_data, voltages)

    return network_data
end

#*------------------------------------------------------------------------------

#=
Returns the bus types of all buses in the NDD in form an ordered dictionary. The different bus types are "Generator", "Load and generator", "Load" and "Empty bus". The ordering is fixed so that the dictionary can be used to plot the buses in a specific order.
=#
#? How to call "empty" buses? Load and generator buses in a substation are connected to these empty buses? Is there only one empty bus per substation?
function get_bustypes(network_data::Dict{String,<:Any})
    ### Active loads
    load = unique(
        [l["load_bus"] for l in collect(values(network_data["load"]))
        if l["status"] == 1]
    )

    ### Active generators
    gen = unique(
        [g["gen_bus"] for g in collect(values(network_data["gen"]))
        if g["gen_status"] == 1]
    )
    
    ### Buses with both load and generation
    load_and_gen = [i for i in load if i in gen]
    
    ### Slack bus
    slack = [
        b["index"] for b in collect(values(network_data["bus"])) 
        if b["bus_type"] == 3
    ]

    ### Remove buses with load and generation from pure loads and generators
    filter!(i -> i ∉ load_and_gen, load)
    filter!(i -> i ∉ load_and_gen, gen)
    ### Remove slack bus from other arrays
    filter!(i -> i ∉ slack, load)
    filter!(i -> i ∉ slack, gen)
    filter!(i -> i ∉ slack, load_and_gen)
    
    ### Empty buses with neither load nor generation
    empty = unique(
        [b["index"] for b in collect(values(network_data["bus"]))
        if b["index"] ∉ vcat(load, gen, load_and_gen, slack)]
    )

    ### Check whether the right number of bus types was obtained
    sum = (
        length(load) + length(gen) + length(load_and_gen) + length(empty) + length(slack)
    )
    N::Int64 = length(network_data["bus"]) # number of buses
    @assert sum == N "$sum bus types were obtained instead of $N"

    return OrderedDict(
        "Generator" => gen,
        "Load and generator" => load_and_gen,
        "Load" => load,
        "Slack" => slack,
        "Empty bus" => empty       
    )
end

#*------------------------------------------------------------------------------

#=
Returns the (string) indices of transmission lines in the network data dictionary that would be considered as underground lines depending on the used maximum length of underground lines and minimum MW-load served. If mode = :sum, the sum of MW-loads of both ends of a line has to exceed min_MW_load, and if mode = :single, one end alone has to exceed min_MW_load.
=#
function get_underground_tl(
        network_data::Dict{String,<:Any}, 
        max_length = 12.875, # in km
        min_MW_load = 2.; # in per-unit (assumed 100 MW base)
        mode = :sum # :sum or :single
    )

    MW_loads = get_MW_loads(network_data) # dictionary of MW loads

    underground_tl = []
    if mode == :single # one end alone has to exceed min_MW_load
        filter!(load -> last(load) >= min_MW_load, MW_loads)
        large_MW_loads = collect(keys(MW_loads))
        ### Go through branches and identify underground transmission lines
        for (i, branch) in network_data["branch"]
            if branch["transformer"] == false && branch["length"] < max_length
                from, to = branch["f_bus"], branch["t_bus"]
                if from in large_MW_loads || to in large_MW_loads
                    push!(underground_tl, i)
                end
            end
        end
    elseif mode == :sum # both ends together have to exceed min_MW_load
        ### Go through branches and identify underground transmission lines
        for (i, branch) in network_data["branch"]
            if branch["transformer"] == false && branch["length"] < max_length
                from, to = branch["f_bus"], branch["t_bus"]
                sum_MW_load = 0.
                if haskey(MW_loads, from)
                    sum_MW_load += MW_loads[from]
                end
                if haskey(MW_loads, to)
                    sum_MW_load += MW_loads[to]
                end
                if sum_MW_load >= min_MW_load
                    push!(underground_tl, i)
                end
            end
        end
    else
        throw(ArgumentError("Unknown mode $mode."))
    end

    return underground_tl
end

#=
Returns a dictionary with load bus indices (Int64) as keys and their respective MW load as values.
=#
function get_MW_loads(network_data::Dict{String,<:Any})
    load_dict = Dict{Int64,Float64}()
    for (i, load) in network_data["load"]
        sub_dict = Dict(load["load_bus"] => load["pd"])
        merge!(+, load_dict, sub_dict)
    end
    return load_dict
end

function get_total_load(network_data::Dict{String,<:Any})
    total_MW_load = 0.
    total_Mvar_load = 0.
    total_MVA_load = 0.

    for load in [l for l in values(network_data["load"]) if l["status"] == 1]
        total_MW_load += load["pd"]
        total_Mvar_load += load["qd"]
        total_MVA_load += sqrt(load["pd"]^2 + load["qd"]^2)
    end

    return total_MW_load, total_Mvar_load, total_MVA_load
end

#*------------------------------------------------------------------------------

# Returns an array of tuples that each contain information about specific branches of interest. The tuples contain the respective branch index, the loading of the branch and a boolean value set to `1`, if the branch is a transformer, or `0`, if the branch is a transmission line. 
"""
    get_branches(network_data::Dict{String,<:Any}; kwargs...)

Returns the indices of branches in `network_data` that are loaded (flow/capacity) by at least `min_loading` (see **keyword arguments**).

# Keyword Arguments

`pf_type` can be set to `"MVA-loading"` (default), `"MW-loading"` or `"Mvar-loading"` and specifies whether the loading due to apparent (MVA), active (MW) or reactive power (Mvar) should be considered. If `network_data["pf_model"] => :dc` meaning that it was last updated using a DC power flow solution, `pf_type` is automatically set to `"MW-loading"`.

`min_loading` is a `Float64` (default `0.0`) that specifies the minimum loading of branches to be returned (according to `pf_type`). Only branches with a loading equal or above `min_loading` will be returned. Hence, `min_loading = 1.0` can be used to obtain overloaded branches.

`br_status` can be set to `1` (default) or `0` and specifies whether active (`1`) or inactive branches (`0`) should be returned.
"""
function get_branches(
        network_data::Dict{String,<:Any}; 
        pf_type = "MVA-loading", # or "MW-loading", "Mvar-loading"
        min_loading = 0.0, # minimum loading of branches
        br_status = 1 # active or inactive branches
    )

    ### Check whether network_data contains a DC-PF solution
    if network_data["pf_model"] == :dc
        pf_type = "MW-loading" # use branch loadings due to active power
    end

    # branches = [
    #     (i, br[pf_type], br["transformer"]) for (i, br) in network_data["branch"] if br["br_status"] == 1 && br[pf_type] >= min_loading
    # ]
    if br_status == 1
        branches = [
            br["index"] for br in values(network_data["branch"]) if 
            br["br_status"] == br_status && br[pf_type] >= min_loading
        ]
    elseif br_status == 0
        branches = [
            br["index"] for br in values(network_data["branch"]) if 
            br["br_status"] == br_status
        ]
    else
        throw(ArgumentError("Unknown branch status $(br_status)!"))
    end

    return branches
end

"""
    get_cc_data(network_data::Dict{String,<:Any}, cc::Set{Int64})
"""
function get_cc_data(
        network_data::Dict{String,<:Any}, cc::Set{Int64}; name = "cc"
    )

    cc_network_data = Dict{String,Any}(
        "baseMVA" => network_data["baseMVA"],
        "per_unit" => network_data["per_unit"],
        "source_type" => network_data["source_type"],
        "name" => name,
        "source_version" => network_data["source_version"],
        "bus" => Dict{String,Any}(),
        "gen" => Dict{String,Any}(),
        "load" => Dict{String,Any}(),
        "shunt" => Dict{String,Any}(),
        "branch" => Dict{String,Any}(),
        "dcline" => Dict{String,Any}(),
        "storage" => Dict{String,Any}(),
        "switch" => Dict{String,Any}()
    )
    
    for b in cc
        cc_network_data["bus"]["$b"] = network_data["bus"]["$b"]
        
    end
end

#*------------------------------------------------------------------------------

#=
Deactivates branches that have been overloaded due to a power flow redistribution (flow/capacity > 1.). The branches are identified by their string indices in branch_ids. The status of these branches is set to 0 in network_data. Parallel branches are not automatically deactivated. Afterwards PowerModels function simplify_network! is run in order to remove dangling buses etc. that might occur due to the failures.
=#
function disable_branches!(
        network_data::Dict{String,<:Any}, 
        branch_ids::Array{Int64,1}
    )
    ### Go through branches to deactivate
    for id in branch_ids
        network_data["branch"][id]["br_status"] = 0 # deactivate branch
    end

    simplify_network!(network_data) # remove dangling buses etc.

    return nothing
end

"""
    deactivate_overloads!(network_data::Dict{String,<:Any})
"""
function deactivate_overloads!(network_data::Dict{String,<:Any})
    ### Get overloaded active branches and deactivate them
    overloaded_br = get_branches(network_data, min_loading=1.)
    if isempty(overloaded_br) == false
        info(LOGGER, 
            "Deactivating the following $(length(overloaded_br)) overloaded branches: $overloaded_br"
        )
    else
        info(LOGGER, "No branches are overloaded.")
    end
    for br in overloaded_br
        network_data["branch"]["$br"]["br_status"] = 0
        network_data["branch"]["$br"]["failure_reason"] = "overload"
    end

    simplify_network!(network_data)

    return overloaded_br
end

#=
Function that deactivates (overhead) transmission lines identified by their indices (string) in tl_ids that have been destroyed by a hurricane. For any given line index all parallel lines also become deactivated. The status of these lines is set to 0 in network_data. Afterwards PowerModels function simplify_network! is run in order to remove dangling buses etc. that might occur due to the failures.
=#
function destroy_tl!(
        network_data::Dict{String,<:Any}, 
        tl_ids::Array # string-ID's of destroyed transmission lines
    )

    ### Go through transmission lines to destroy
    for id in tl_ids
        ## Deactivate id and all parallel lines
        from = network_data["branch"][id]["f_bus"]
        to = network_data["branch"][id]["t_bus"]
        ## Find all parallel branches
        ids = findall(
            b -> (b["f_bus"], b["t_bus"]) == (from, to), 
            network_data["branch"]
        )
        for idx in ids
            network_data["branch"][idx]["br_status"] = 0
            network_data["branch"][idx]["failure_reason"] = "wind"
        end
    end

    simplify_network!(network_data) # remove dangling buses etc.

    return nothing
end

#=
Deactivates all components (generator, loads, etc.) in a connected component.
=#
function _deactivate_cc!(
        network_data::Dict{String,<:Any}, 
        active_cc_gens::Array{Int64,1}, 
        active_cc_loads::Array{Int64,1}
    )

    ### Deactivate all active generators in connected component
    for g in active_cc_gens
        network_data["gen"]["$g"]["gen_status"] = 0
    end

    ### Deactivate all active loads in connected component
    for l in active_cc_loads
        network_data["load"]["$l"]["status"] = 0
    end

    ### Deactivate rest of connected component (branches, dangling buses etc.)
    simplify_network!(network_data)
    
    return nothing
end

#*------------------------------------------------------------------------------

"""
    check_voltage_bounds(network_data::Dict{String,Any}; kwargs...)

Checks whether active buses in `network_data` have voltage magnitudes outside of their tolerance interval (usually ranging from 0.9 p.u. to 1.1 p.u.) and returns a dictionary with the bus indices as keys and the respective active loads connected to the bus as values (PQ-generators are excluded). ITCPG's logger prints a `warn` message everytime a voltage violation is detected.

# Keyword Arguments

`limit` can be set to `"lower"`, `"upper"` or `"both"` (default) and specifies which voltage bound should be considered. If e.g. `limit = "lower"` is chosen, only buses with under-voltages are returned.

`filter_active_loads` can be set to `true` or `false` (default) and specifies whether only buses with active loads connected to them should be returned (useful for under-voltage load shedding). 

See also: [`run_outer_ac_pf!`](@ref), [`run_uvls!`](@ref)
"""
function check_voltage_bounds(
        network_data::Dict{String,<:Any}; 
        limit = "both", # considered voltage bounds
        filter_active_loads = false # whether to filter buses with active loads
    )

    vv_buses = Dict{String,Array{String,1}}()
    ### Go through buses and check their voltage magnitudes
    for (key, bus) in network_data["bus"]
        vm = bus["vm"]
        if limit == "both" || limit == "lower"
            if vm < bus["vmin"] && bus["bus_type"] != 4
                warn(LOGGER, 
                    "Bus $key (type $(network_data["bus"][key]["bus_type"])) has a voltage magnitude of |V|=$vm that is below its lower bound of $(bus["vmin"])" 
                )
                ### Get string indices of active loads connected to the bus
                connected_loads = [
                    l_id for (l_id, l) in network_data["load"] if 
                    l["load_bus"] == parse(Int64, key) && l["status"] == 1 && 
                    haskey(l, "is_gen") == false # exclude PQ-generators
                ]
                vv_buses[key] = connected_loads
            end
        end
        if limit == "both" || limit == "upper"
            if vm > bus["vmax"] && bus["bus_type"] != 4
                warn(LOGGER, 
                    "Bus $key (type $(network_data["bus"][key]["bus_type"])) has a voltage magnitude of |V|=$vm that is above its upper bound of $(bus["vmax"])" 
                )
                ### Get string indices of active loads connected to the bus
                connected_loads = [
                    l_id for (l_id, l) in network_data["load"] if 
                    l["load_bus"] == parse(Int64, key) && l["status"] == 1 && 
                    haskey(l, "is_gen") == false # exclude PQ-generators
                ]
                vv_buses[key] = connected_loads
            end
        end
    end

    ### Check whether to exclude buses without active loads
    if filter_active_loads == true
        filter!(x -> isempty(last(x)) == false, vv_buses)
    end

    return vv_buses
end

function check_gen_bounds(network_data::Dict{String,<:Any})
    for (i, gen) in network_data["gen"]
        if gen["gen_status"] == 1
            if gen["pg"] > gen["pmax"] || gen["pg"] < gen["pmin"]
                warn(LOGGER, 
                    "Active power dispatch of generator $i lies outside of interval [$(gen["pmin"]), $(gen["pmax"])]: $(gen["pg"])"
                )
            end
            if gen["qg"] > gen["qmax"] || gen["qg"] < gen["qmin"]
                warn(LOGGER,
                    "Reactive power dispatch of generator $i lies outside of interval [$(gen["qmin"]), $(gen["qmax"])]: $(gen["qg"])"
                )
            end
        end
    end

    return nothing
end

function check_gen_bounds(
        network_data::Dict{String,<:Any},
        cc::Set{Int64}, # set of bus indices in connected component
        cc_i::Int64; # index/number of connected components
        power = "active" # "active": active power, "both": reactive power too
    )

    cc_slack = [
        b["index"] for b in values(network_data["bus"]) if 
        b["bus_type"] == 3 && b["index"] ∈ cc
    ]
    
    if length(cc_slack) != 1
        warn(LOGGER, "(Sub-)Grid #$cc_i has $(length(cc_slack)) slack buses!")
    end

    for (i, gen) in network_data["gen"]
        if gen["gen_bus"] ∈ cc && gen["gen_status"] == 1
            if gen["gen_bus"] ∈ cc_slack
                info(LOGGER, 
                    "Active power dispatch of slack bus $i in #$cc_i with interval [$(gen["pmin"]), $(gen["pmax"])]: $(gen["pg"])"
                )
            else
                if gen["pg"] > gen["pmax"] || gen["pg"] < gen["pmin"]
                    warn(LOGGER, "Active power dispatch of generator $i in #$cc_i lies outside of interval [$(gen["pmin"]), $(gen["pmax"])]: $(gen["pg"])")
                end
                if power == "both"
                    if gen["qg"] > gen["qmax"] || gen["qg"] < gen["qmin"]
                        warn(LOGGER,
                            "Reactive power dispatch of generator $i lies outside of interval [$(gen["qmin"]), $(gen["qmax"])]: $(gen["qg"])"
                        )
                    end
                end
            end
        end
    end

    return nothing
end