Merge 6e4b4b2 into ef19e15

CHANGELOG.md

@@ -2,11 +2,13 @@ PowerModels.jl Change Log
 =========================
 
 ### Staged
+- Added support for matrix-based analysis of basic network data (#728)
 - Added support for `relax_integrality` in `run_model`
 - Added `export_pti` to write a PSSE file (#752)
 - Added `parse_files` to create a PM-multinetwork from multiples files
-- Add support for convering matpower ramp rates into per-unit (#561)
+- Added support for convering matpower ramp rates into per-unit (#561)
 - Fixed bug in dual reporting in `constraint_power_balance_ls` (#741)
+- Fixed sign convetion for power injection in `calc_bus_injection`
 
 ### v0.17.3
 - Added a to file variant of `export_matpower`

docs/make.jl

@@ -16,7 +16,8 @@ makedocs(
             "Storage Model" => "storage.md",
             "Switch Model" => "switch.md",
             "Multi Networks" => "multi-networks.md",
-            "Utilities" => "utilities.md"
+            "Utilities" => "utilities.md",
+            "Basic Data Utilities" => "basic-data-utilities.md"
         ],
         "Library" => [
             "Network Formulations" => "formulations.md",

docs/src/basic-data-utilities.md

@@ -0,0 +1,112 @@
+# Basic Data Utilities
+
+By default PowerModels uses a data model that captures the bulk of the features
+of realistic transmission network datasets such as, inactive devices, breakers
+and HVDC lines. However, these features preclude popular matrix-based
+analysis of power network datasets such as incidence, admittance, and power
+transfer distribution factor (PTDF) matrices. To support these types of
+analysis PowerModels introduces the concept of a _basic networks_, which are
+network datasets that satisfy the properties required to interpret the system
+in a matrix form.
+
+The `make_basic_network` is provided to ensure that a given network dataset
+satisfies the properties required for a matrix interpretation (the specific
+requirements are outlined in the function documentation block). If the given
+dataset does not satisfy the properties, `make_basic_network` transforms the
+dataset to enforce them.
+
+```@docs
+make_basic_network
+```
+
+The standard procedure for loading basic network data is as follows,
+```julia
+data = make_basic_network(parse_file("<path to network data file>"))
+```
+modifications to the original network data file are indicated by logging
+messages in the terminal.
+
+!!! tip
+    If `make_basic_network` results in significant changes to a dataset,
+    `export_file` can be used to inspect and modify the new derivative dataset
+    that conforms to the basic network requirements.
+
+
+## Matrix-Based Data
+
+Using a basic network dataset the following functions can be used to extract
+key power system quantities in vectors and matrix forms. The prefix `_basic_`
+distinguishes these functions from similar tools that operate on any type of
+PowerModels data, including those that are not amenable to a vector/matrix
+format.
+
+```@docs
+calc_basic_bus_voltage
+calc_basic_bus_injection
+calc_basic_branch_series_impedance
+calc_basic_incidence_matrix
+calc_basic_admittance_matrix
+calc_basic_susceptance_matrix
+calc_basic_branch_susceptance_matrix
+calc_basic_ptdf_matrix
+calc_basic_ptdf_row
+```
+
+!!! warning
+    Several variants of the real-valued susceptance matrix are possible.
+    PowerModels uses the version based on inverse of branch series impedance,
+    that is `imag(inv(r + x im))`. One may observe slightly different results
+    when compared to tools that use other variants such as `1/x`.
+
+
+## Matrix-Based Computations
+
+Matrix-based network data can be combined to compute a number of useful
+quantities. For example, by combining the incidence matrix and the series
+impedance one can drive the susceptance and branch susceptance matrices as follows,
+
+```julia
+import LinearAlgebra: Diagonal
+
+bz = calc_basic_branch_series_impedance(data)
+A  = calc_basic_incidence_matrix(data)
+
+Y  = imag(Diagonal(inv.(bz)))
+B  = A'*Y*A    # equivalent to calc_basic_susceptance_matrix
+BB = (A'*Y)'   # equivalent to calc_basic_branch_susceptance_matrix
+```
+
+The bus voltage angles can be combined with the susceptance and branch susceptance
+matrices to observe how power flows through the network as follows,
+
+```julia
+va = angle.(calc_basic_bus_voltage(data))
+B  = calc_basic_susceptance_matrix(data)
+BB = calc_basic_branch_susceptance_matrix(data)
+
+bus_injection =  -B * va
+branch_power  = -BB * va
+```
+
+In the inverse operation, bus injection values can be combined with a PTDF matrix to compute branch flow values as follows,
+
+```julia
+bi   = real(calc_basic_bus_injection(data))
+PTDF = calc_basic_ptdf_matrix(data)
+
+branch_power = PTDF * bi
+```
+
+Finally, the following function provides a tool to solve a DC power flow on
+basic network data using Julia's native linear equation solver,
+
+```@docs
+compute_basic_dc_pf
+```
+
+!!! tip
+    By default PowerModels uses Julia's SparseArrays to ensure the best
+    performance of matrix operations on large power network datasets.
+    The function `Matrix(sparse_array)` can be used to covert a sparse matrix
+    into a full matrix when that is preferred.
+

src/PowerModels.jl

@@ -43,6 +43,7 @@ include("io/pti.jl")
 include("io/psse.jl")
 
 include("core/data.jl")
+include("core/data_basic.jl")
 include("core/ref.jl")
 include("core/base.jl")
 include("core/types.jl")

src/core/admittance_matrix.jl

@@ -295,8 +295,8 @@ function calc_bus_injection(data::Dict{String,<:Any})
     for (i,bus) in data["bus"]
         if bus["bus_type"] != 4
             bvals = bus_values[bus["index"]]
-            p_delta = - bvals["pg"] + bvals["ps"] + bvals["pd"]
-            q_delta = - bvals["qg"] + bvals["qs"] + bvals["qd"]
+            p_delta = bvals["pg"] - bvals["ps"] - bvals["pd"]
+            q_delta = bvals["qg"] - bvals["qs"] - bvals["qd"]
         else
             p_delta = NaN
             q_delta = NaN
@@ -333,7 +333,7 @@ function solve_theta(am::AdmittanceMatrix, bus_injection::Vector{Float64})
     end
     bi[am.ref_idx] = 0.0
 
-    theta = m \ -bi
+    theta = m \ bi
 
     return theta
 end

src/core/data.jl

@@ -2197,7 +2197,7 @@ function _standardize_cost_terms!(components::Dict{String,<:Any}, comp_order::In
             comp["ncost"] = comp_order
             #println("std gen cost: $(comp["cost"])")
 
-            Memento.warn(_LOGGER, "Updated $(cost_comp_name) $(comp["index"]) cost function with order $(length(current_cost)) to a function of order $(comp_order): $(comp["cost"])")
+            Memento.info(_LOGGER, "updated $(cost_comp_name) $(comp["index"]) cost function with order $(length(current_cost)) to a function of order $(comp_order): $(comp["cost"])")
             push!(modified, comp["index"])
         end
     end
@@ -2585,21 +2585,24 @@ end
 ""
 function _select_largest_component!(data::Dict{String,<:Any})
     ccs = calc_connected_components(data)
-    Memento.info(_LOGGER, "found $(length(ccs)) components")
 
-    ccs_order = sort(collect(ccs); by=length)
-    largest_cc = ccs_order[end]
+    if length(ccs) > 1
+        Memento.info(_LOGGER, "found $(length(ccs)) components")
 
-    Memento.info(_LOGGER, "largest component has $(length(largest_cc)) buses")
+        ccs_order = sort(collect(ccs); by=length)
+        largest_cc = ccs_order[end]
 
-    for (i,bus) in data["bus"]
-        if bus["bus_type"] != 4 && !(bus["index"] in largest_cc)
-            bus["bus_type"] = 4
-            Memento.info(_LOGGER, "deactivating bus $(i) due to small connected component")
+        Memento.info(_LOGGER, "largest component has $(length(largest_cc)) buses")
+
+        for (i,bus) in data["bus"]
+            if bus["bus_type"] != 4 && !(bus["index"] in largest_cc)
+                bus["bus_type"] = 4
+                Memento.info(_LOGGER, "deactivating bus $(i) due to small connected component")
+            end
         end
-    end
 
-    correct_reference_buses!(data)
+        correct_reference_buses!(data)
+    end
 end
 
 
@@ -2770,3 +2773,164 @@ function _cc_dfs(i, neighbors, component_lookup, touched)
         end
     end
 end
+
+
+
+"""
+given a network data dict and a mapping of current-bus-ids to new-bus-ids
+modifies the data dict to reflect the proposed new bus ids.
+"""
+function update_bus_ids!(data::Dict{String,<:Any}, bus_id_map::Dict{Int,Int}; injective=true)
+    if _IM.ismultinetwork(data)
+        for (i,nw_data) in data["nw"]
+            _update_bus_ids!(nw_data, bus_id_map; injective=injective)
+        end
+    else
+         _update_bus_ids!(data, bus_id_map; injective=injective)
+    end
+end
+
+
+function _update_bus_ids!(data::Dict{String,<:Any}, bus_id_map::Dict{Int,Int}; injective=true)
+    # verify bus id map is injective
+    if injective
+        new_bus_ids = Set{Int}()
+        for (i,bus) in data["bus"]
+            new_id = get(bus_id_map, bus["index"], bus["index"])
+            if !(new_id in new_bus_ids)
+                push!(new_bus_ids, new_id)
+            else
+                Memento.error(_LOGGER, "bus id mapping given to update_bus_ids has an id clash on new bus id $(new_id)")
+            end
+        end
+    end
+
+
+    # start renumbering process
+    renumbered_bus_dict = Dict{String,Any}()
+
+    for (i,bus) in data["bus"]
+        new_id = get(bus_id_map, bus["index"], bus["index"])
+        bus["index"] = new_id
+        bus["bus_i"] = new_id
+        renumbered_bus_dict["$new_id"] = bus
+    end
+
+    data["bus"] = renumbered_bus_dict
+
+
+    # update bus numbering in dependent components
+    for (i, load) in data["load"]
+        load["load_bus"] = get(bus_id_map, load["load_bus"], load["load_bus"])
+    end
+
+    for (i, shunt) in data["shunt"]
+        shunt["shunt_bus"] = get(bus_id_map, shunt["shunt_bus"], shunt["shunt_bus"])
+    end
+
+    for (i, gen) in data["gen"]
+        gen["gen_bus"] = get(bus_id_map, gen["gen_bus"], gen["gen_bus"])
+    end
+
+    for (i, strg) in data["storage"]
+        strg["storage_bus"] = get(bus_id_map, strg["storage_bus"], strg["storage_bus"])
+    end
+
+
+    for (i, switch) in data["switch"]
+        switch["f_bus"] = get(bus_id_map, switch["f_bus"], switch["f_bus"])
+        switch["t_bus"] = get(bus_id_map, switch["t_bus"], switch["t_bus"])
+    end
+
+    branches = []
+    if haskey(data, "branch")
+        append!(branches, values(data["branch"]))
+    end
+
+    if haskey(data, "ne_branch")
+        append!(branches, values(data["ne_branch"]))
+    end
+
+    for branch in branches
+        branch["f_bus"] = get(bus_id_map, branch["f_bus"], branch["f_bus"])
+        branch["t_bus"] = get(bus_id_map, branch["t_bus"], branch["t_bus"])
+    end
+
+    for (i, dcline) in data["dcline"]
+        dcline["f_bus"] = get(bus_id_map, dcline["f_bus"], dcline["f_bus"])
+        dcline["t_bus"] = get(bus_id_map, dcline["t_bus"], dcline["t_bus"])
+    end
+end
+
+
+
+"""
+given a network data dict merges buses that are connected by closed switches
+converting the dataset into a pure bus-branch model.
+"""
+function resolve_swithces!(data::Dict{String,<:Any})
+    if _IM.ismultinetwork(data)
+        for (i,nw_data) in data["nw"]
+            _resolve_swithces!(nw_data)
+        end
+    else
+         _resolve_swithces!(data)
+    end
+end
+
+""
+function _resolve_swithces!(data::Dict{String,<:Any})
+    if length(data["switch"]) <= 0
+        return
+    end
+
+    bus_sets = Dict{Int,Set{Int}}()
+
+    switch_status_key = pm_component_status["switch"]
+    switch_status_value = pm_component_status_inactive["switch"]
+
+    for (i,switch) in data["switch"]
+        if switch[switch_status_key] != switch_status_value && switch["state"] == 1
+            if !haskey(bus_sets, switch["f_bus"])
+                bus_sets[switch["f_bus"]] = Set{Int}([switch["f_bus"]])
+            end
+            if !haskey(bus_sets, switch["t_bus"])
+                bus_sets[switch["t_bus"]] = Set{Int}([switch["t_bus"]])
+            end
+
+            merged_set = Set{Int}([bus_sets[switch["f_bus"]]..., bus_sets[switch["t_bus"]]...])
+            bus_sets[switch["f_bus"]] = merged_set
+            bus_sets[switch["t_bus"]] = merged_set
+        end
+    end
+
+    bus_id_map = Dict{Int,Int}()
+    for bus_set in Set(values(bus_sets))
+        bus_min = minimum(bus_set)
+        Memento.info(_LOGGER, "merged buses $(join(bus_set, ",")) in to bus $(bus_min) based on switch status")
+        for i in bus_set
+            if i != bus_min
+                bus_id_map[i] = bus_min
+            end
+        end
+    end
+
+    update_bus_ids!(data, bus_id_map, injective=false)
+
+    for (i, branch) in data["branch"]
+        if branch["f_bus"] == branch["t_bus"]
+            Memento.warn(_LOGGER, "switch removal resulted in both sides of branch $(i) connect to bus $(branch["f_bus"]), deactivating branch")
+            branch[pm_component_status["branch"]] = pm_component_status_inactive["branch"]
+        end
+    end
+
+    for (i, dcline) in data["dcline"]
+        if dcline["f_bus"] == dcline["t_bus"]
+            Memento.warn(_LOGGER, "switch removal resulted in both sides of dcline $(i) connect to bus $(branch["f_bus"]), deactivating dcline")
+            branch[pm_component_status["dcline"]] = pm_component_status_inactive["dcline"]
+        end
+    end
+
+    Memento.info(_LOGGER, "removed $(length(data["switch"])) switch components")
+    data["switch"] = Dict{String,Any}()
+end