From e0f0873374f13dd47a6be2ae87c0d7f9c17eec7e Mon Sep 17 00:00:00 2001 From: Carleton Coffrin Date: Tue, 2 Feb 2021 10:24:40 -0700 Subject: [PATCH] Basic-Data Idea (#756) * basic data support functions * add tests for bus number updates and switch resolution * add basic data tests * skip select_largest_component logic when only a single component is found * add basic data docs * update basic data docs * update calc_bus_injection to match standard conventions, injection is positive * comment preliminary basic jacobian implementation * add note to changelog --- CHANGELOG.md | 4 +- docs/make.jl | 3 +- docs/src/basic-data-utilities.md | 112 +++++++ src/PowerModels.jl | 1 + src/core/admittance_matrix.jl | 6 +- src/core/data.jl | 186 +++++++++++- src/core/data_basic.jl | 498 +++++++++++++++++++++++++++++++ src/prob/pf.jl | 10 +- test/data-basic.jl | 203 +++++++++++++ test/data.jl | 42 ++- test/runtests.jl | 2 + 11 files changed, 1045 insertions(+), 22 deletions(-) create mode 100644 docs/src/basic-data-utilities.md create mode 100644 src/core/data_basic.jl create mode 100644 test/data-basic.jl diff --git a/CHANGELOG.md b/CHANGELOG.md index f9b48708f..680bf1210 100644 --- a/CHANGELOG.md +++ b/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` diff --git a/docs/make.jl b/docs/make.jl index 90dfd9704..628ee1b58 100644 --- a/docs/make.jl +++ b/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", diff --git a/docs/src/basic-data-utilities.md b/docs/src/basic-data-utilities.md new file mode 100644 index 000000000..592d595cd --- /dev/null +++ b/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("")) +``` +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. + diff --git a/src/PowerModels.jl b/src/PowerModels.jl index ea5dd12b0..ba3d12192 100644 --- a/src/PowerModels.jl +++ b/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") diff --git a/src/core/admittance_matrix.jl b/src/core/admittance_matrix.jl index 681847c3f..196321d3b 100644 --- a/src/core/admittance_matrix.jl +++ b/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 diff --git a/src/core/data.jl b/src/core/data.jl index d54d8d23e..8791c120e 100644 --- a/src/core/data.jl +++ b/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 diff --git a/src/core/data_basic.jl b/src/core/data_basic.jl new file mode 100644 index 000000000..e94a88397 --- /dev/null +++ b/src/core/data_basic.jl @@ -0,0 +1,498 @@ +# tools for working with the "basic" versions of PowerModels data dict + +""" +given a powermodels data dict produces a new data dict that conforms to the +following basic network model requirements. +- no dclines +- no switches +- no inactive components +- all components are numbered from 1-to-n +- the network forms a single connected component +- there exactly one phase angle reference bus +- generation cost functions are quadratic +- all branches have explicit thermal limits +- phase shift on all transformers is set to 0.0 +- bus shunts have 0.0 conductance values +users requiring any of the features listed above for their analysis should use +the non-basic PowerModels routines. +""" +function make_basic_network(data::Dict{String,<:Any}) + if _IM.ismultinetwork(data) + Memento.error(_LOGGER, "make_basic_network does not support multinetwork data") + end + + # make a copy of data so that modifications do not change the input data + data = deepcopy(data) + + # TODO transform PWL costs into linear costs + for (i,gen) in data["gen"] + if get(gen, "cost_model", 2) != 2 + Memento.error(_LOGGER, "make_basic_network only supports network data with polynomial cost functions, generator $(i) has a piecewise linear cost function") + end + end + standardize_cost_terms!(data, order=2) + + # set conductance to zero on all shunts + for (i,shunt) in data["shunt"] + if !isapprox(shunt["gs"], 0.0) + Memento.warn(_LOGGER, "setting conductance on shunt $(i) from $(shunt["gs"]) to 0.0") + shunt["gs"] = 0.0 + end + end + + # ensure that branch components always have a rate_a value + calc_thermal_limits!(data) + + # set phase shift to zero on all branches + for (i,branch) in data["branch"] + if !isapprox(branch["shift"], 0.0) + Memento.warn(_LOGGER, "setting phase shift on branch $(i) from $(branch["shift"]) to 0.0") + branch["shift"] = 0.0 + end + end + + # ensure single connected component + select_largest_component!(data) + + # ensure that components connected in inactive buses are also inactive + propagate_topology_status!(data) + + # ensure there is exactly one reference bus + ref_buses = Set{Int}() + for (i,bus) in data["bus"] + if bus["bus_type"] == 3 + push!(ref_buses, bus["index"]) + end + end + if length(ref_buses) > 1 + Memento.warn(_LOGGER, "network data specifies $(length(ref_buses)) reference buses") + for ref_bus_id in ref_buses + data["bus"]["$(ref_bus_id)"]["bus_type"] = 2 + end + ref_buses = Set{Int}() + end + if length(ref_buses) == 0 + gen = _biggest_generator(data["gen"]) + @assert length(gen) > 0 + gen_bus = gen["gen_bus"] + ref_bus = data["bus"]["$(gen_bus)"] + ref_bus["bus_type"] = 3 + Memento.warn(_LOGGER, "setting bus $(gen_bus) as reference based on generator $(gen["index"])") + end + + # remove switches by merging buses + resolve_swithces!(data) + + # switch resolution can result in new parallel branches + correct_branch_directions!(data) + + # set remaining unsupported components as inactive + dcline_status_key = pm_component_status["dcline"] + dcline_inactive_status = pm_component_status_inactive["dcline"] + for (i,dcline) in data["dcline"] + dcline[dcline_status_key] = dcline_inactive_status + end + + # remove inactive components + for (comp_key, status_key) in pm_component_status + comp_count = length(data[comp_key]) + status_inactive = pm_component_status_inactive[comp_key] + data[comp_key] = _filter_inactive_components(data[comp_key], status_key=status_key, status_inactive_value=status_inactive) + if length(data[comp_key]) < comp_count + Memento.info(_LOGGER, "removed $(comp_count - length(data[comp_key])) inactive $(comp_key) components") + end + end + + # re-number non-bus component ids + for comp_key in keys(pm_component_status) + if comp_key != "bus" + data[comp_key] = _renumber_components(data[comp_key]) + end + end + + # renumber bus ids + bus_ordered = sort([bus for (i,bus) in data["bus"]], by=(x) -> x["index"]) + + bus_id_map = Dict{Int,Int}() + for (i,bus) in enumerate(bus_ordered) + bus_id_map[bus["index"]] = i + end + update_bus_ids!(data, bus_id_map) + + data["basic_network"] = true + + return data +end + +""" +given a component dict returns a new dict where inactive components have been +removed. +""" +function _filter_inactive_components(comp_dict::Dict{String,<:Any}; status_key="status", status_inactive_value=0) + filtered_dict = Dict{String,Any}() + + for (i,comp) in comp_dict + if comp[status_key] != status_inactive_value + filtered_dict[i] = comp + end + end + + return filtered_dict +end + +""" +given a component dict returns a new dict where components have been renumbered +from 1-to-n ordered by the increasing values of the orginal component id. +""" +function _renumber_components(comp_dict::Dict{String,<:Any}) + renumbered_dict = Dict{String,Any}() + + comp_ordered = sort([comp for (i,comp) in comp_dict], by=(x) -> x["index"]) + + for (i,comp) in enumerate(comp_ordered) + comp = deepcopy(comp) + comp["index"] = i + renumbered_dict["$i"] = comp + end + + return renumbered_dict +end + + + +""" +given a basic network data dict, returns a complex valued vector of bus voltage +values in rectangular coordinates as they appear in the network data. +""" +function calc_basic_bus_voltage(data::Dict{String,<:Any}) + if !get(data, "basic_network", false) + Memento.warn(_LOGGER, "calc_basic_bus_voltage requires basic network data and given data may be incompatible. make_basic_network can be used to transform data into the appropriate form.") + end + + b = [bus for (i,bus) in data["bus"] if bus["bus_type"] != 4] + bus_ordered = sort(b, by=(x) -> x["index"]) + + return [bus["vm"]*cos(bus["va"]) + bus["vm"]*sin(bus["va"])im for bus in bus_ordered] +end + +""" +given a basic network data dict, returns a complex valued vector of bus power +injections as they appear in the network data. +""" +function calc_basic_bus_injection(data::Dict{String,<:Any}) + if !get(data, "basic_network", false) + Memento.warn(_LOGGER, "calc_basic_bus_injection requires basic network data and given data may be incompatible. make_basic_network can be used to transform data into the appropriate form.") + end + + bi_dict = calc_bus_injection(data) + bi_vect = [bi_dict[1][i] + bi_dict[2][i]im for i in 1:length(data["bus"])] + + return bi_vect +end + +""" +given a basic network data dict, returns a complex valued vector of branch +series impedances. +""" +function calc_basic_branch_series_impedance(data::Dict{String,<:Any}) + if !get(data, "basic_network", false) + Memento.warn(_LOGGER, "calc_basic_branch_series_impedance requires basic network data and given data may be incompatible. make_basic_network can be used to transform data into the appropriate form.") + end + + b = [branch for (i,branch) in data["branch"] if branch["br_status"] != 0] + branch_ordered = sort(b, by=(x) -> x["index"]) + + return [branch["br_r"] + branch["br_x"]im for branch in branch_ordered] +end + + +""" +given a basic network data dict, returns a sparse integer valued incidence +matrix with one row for each branch and one column for each bus in the network. +In each branch row a +1 is used to indicate the _from_ bus and -1 is used to +indicate _to_ bus. +""" +function calc_basic_incidence_matrix(data::Dict{String,<:Any}) + if !get(data, "basic_network", false) + Memento.warn(_LOGGER, "calc_basic_incidence_matrix requires basic network data and given data may be incompatible. make_basic_network can be used to transform data into the appropriate form.") + end + + I = Int64[] + J = Int64[] + V = Int64[] + + b = [branch for (i,branch) in data["branch"] if branch["br_status"] != 0] + branch_ordered = sort(b, by=(x) -> x["index"]) + for (i,branch) in enumerate(branch_ordered) + push!(I, i); push!(J, branch["f_bus"]); push!(V, 1) + push!(I, i); push!(J, branch["t_bus"]); push!(V, -1) + end + + return sparse(I,J,V) +end + +""" +given a basic network data dict, returns a sparse complex valued admittance +matrix with one row and column for each bus in the network. +""" +function calc_basic_admittance_matrix(data::Dict{String,<:Any}) + if !get(data, "basic_network", false) + Memento.warn(_LOGGER, "calc_basic_admittance_matrix requires basic network data and given data may be incompatible. make_basic_network can be used to transform data into the appropriate form.") + end + + am = calc_admittance_matrix(data) + + # conj can be removed once #734 is resolved + return conj.(am.matrix) +end + + +""" +given a basic network data dict, returns a sparse real valued susceptance +matrix with one row and column for each bus in the network. +This susceptance matrix reflects the imaginary part of an admittance +matrix that only considers the branch series impedance. +""" +function calc_basic_susceptance_matrix(data::Dict{String,<:Any}) + if !get(data, "basic_network", false) + Memento.warn(_LOGGER, "calc_basic_susceptance_matrix requires basic network data and given data may be incompatible. make_basic_network can be used to transform data into the appropriate form.") + end + + am = calc_susceptance_matrix(data) + + # -1.0 can be removed once #734 is resolved + return -1.0*am.matrix +end + + +""" +given a basic network data dict, returns a sparse real valued branch susceptance +matrix with one row for each branch and one column for each bus in the network. +Multiplying the branch susceptance matrix by bus phase angels yields a vector +active power flow values for each branch. +""" +function calc_basic_branch_susceptance_matrix(data::Dict{String,<:Any}) + if !get(data, "basic_network", false) + Memento.warn(_LOGGER, "calc_basic_branch_susceptance_matrix requires basic network data and given data may be incompatible. make_basic_network can be used to transform data into the appropriate form.") + end + + I = Int64[] + J = Int64[] + V = Float64[] + + b = [branch for (i,branch) in data["branch"] if branch["br_status"] != 0] + branch_ordered = sort(b, by=(x) -> x["index"]) + for (i,branch) in enumerate(branch_ordered) + g,b = calc_branch_y(branch) + push!(I, i); push!(J, branch["f_bus"]); push!(V, b) + push!(I, i); push!(J, branch["t_bus"]); push!(V, -b) + end + + return sparse(I,J,V) +end + +""" +given a basic network data dict, computes real valued vector of bus voltage +phase angles by solving a dc power flow. +""" +function compute_basic_dc_pf(data::Dict{String,<:Any}) + if !get(data, "basic_network", false) + Memento.warn(_LOGGER, "compute_basic_dc_pf requires basic network data and given data may be incompatible. make_basic_network can be used to transform data into the appropriate form.") + end + + num_bus = length(data["bus"]) + ref_bus_id = reference_bus(data)["index"] + + bi = real(calc_basic_bus_injection(data)) + + sm = calc_basic_susceptance_matrix(data) + + for i in 1:num_bus + if i == ref_bus_id + # TODO improve scaling of this value + sm[i,i] = 1.0 + else + if !iszero(sm[ref_bus_id,i]) + sm[ref_bus_id,i] = 0.0 + end + end + end + bi[ref_bus_id] = 0.0 + + theta = -sm \ bi + + return theta +end + + + +""" +given a basic network data dict, returns a real valued ptdf matrix with one +row for each branch and one column for each bus in the network. +Multiplying the ptdf matrix by bus injection values yields a vector +active power flow values on each branch. +""" +function calc_basic_ptdf_matrix(data::Dict{String,<:Any}) + if !get(data, "basic_network", false) + Memento.warn(_LOGGER, "calc_basic_ptdf_matrix requires basic network data and given data may be incompatible. make_basic_network can be used to transform data into the appropriate form.") + end + + num_bus = length(data["bus"]) + num_branch = length(data["branch"]) + + # -1.0 can be removed once #734 is resolved + b_inv = calc_susceptance_matrix_inv(data).matrix + + ptdf = zeros(num_branch, num_bus) + for (i,branch) in data["branch"] + branch_idx = branch["index"] + bus_fr = branch["f_bus"] + bus_to = branch["t_bus"] + g,b = calc_branch_y(branch) + for n in 1:num_bus + ptdf[branch_idx, n] = -b*(b_inv[bus_fr, n] - b_inv[bus_to, n]) + end + end + + return ptdf +end + +""" +given a basic network data dict and a branch index returns a row of the ptdf +matrix reflecting that branch. +""" +function calc_basic_ptdf_row(data::Dict{String,<:Any}, branch_index::Int) + if !get(data, "basic_network", false) + Memento.warn(_LOGGER, "calc_basic_ptdf_row requires basic network data and given data may be incompatible. make_basic_network can be used to transform data into the appropriate form.") + end + + if branch_index < 1 || branch_index > length(data["branch"]) + Memento.error(_LOGGER, "branch index of $(branch_index) is out of bounds, valid values are $(1)-$(length(data["branch"]))") + end + branch = data["branch"]["$(branch_index)"] + g,b = calc_branch_y(branch) + + num_bus = length(data["bus"]) + num_branch = length(data["branch"]) + + am = calc_susceptance_matrix(data) + + if_fr = injection_factors_va(am, branch["f_bus"]) + if_to = injection_factors_va(am, branch["t_bus"]) + + ptdf_column = zeros(num_bus) + for n in 1:num_bus + # -1.0 can be removed once #734 is resolved + ptdf_column[n] = -b*(get(if_fr, n, 0.0) - get(if_to, n, 0.0)) + end + + return ptdf_column +end + + +#= +""" +given a basic network data dict, returns a sparse real valued Jacobian matrix +of the ac power flow problem. The power variables are ordered by p and then q +while voltage values are ordered by voltage angle and then voltage magnitude. +""" +function calc_basic_jacobian_matrix(data::Dict{String,<:Any}) + if !get(data, "basic_network", false) + Memento.warn(_LOGGER, "calc_basic_jacobian_matrix requires basic network data and given data may be incompatible. make_basic_network can be used to transform data into the appropriate form.") + end + + num_bus = length(data["bus"]) + + vm = [bus["vm"] for (i,bus) in data["bus"]] + va = [bus["va"] for (i,bus) in data["bus"]] + am = calc_admittance_matrix(data) + + + neighbors = [Set{Int}([i]) for i in 1:num_bus] + I, J, V = findnz(am.matrix) + for nz in eachindex(V) + push!(neighbors[I[nz]], J[nz]) + push!(neighbors[J[nz]], I[nz]) + end + + + J0_I = Int64[] + J0_J = Int64[] + J0_V = Float64[] + + for i in 1:num_bus + f_i_r = i + f_i_i = i + num_bus + + for j in neighbors[i] + x_j_fst = j + num_bus + x_j_snd = j + + push!(J0_I, f_i_r); push!(J0_J, x_j_fst); push!(J0_V, 0.0) + push!(J0_I, f_i_r); push!(J0_J, x_j_snd); push!(J0_V, 0.0) + push!(J0_I, f_i_i); push!(J0_J, x_j_fst); push!(J0_V, 0.0) + push!(J0_I, f_i_i); push!(J0_J, x_j_snd); push!(J0_V, 0.0) + end + end + + J = sparse(J0_I, J0_J, J0_V) + + + for i in 1:num_bus + f_i_r = i + f_i_i = i + num_bus + + for j in neighbors[i] + x_j_fst = j + num_bus + x_j_snd = j + + bus_type = data["bus"]["$(j)"]["bus_type"] + if bus_type == 1 + if i == j + y_ii = am.matrix[i,i] + J[f_i_r, x_j_fst] = 2*real(y_ii)*vm[i] + sum( real(am.matrix[i,k]) * vm[k] * cos(va[i] - va[k]) - imag(am.matrix[i,k]) * vm[k] * sin(va[i] - va[k]) for k in neighbors[i] if k != i) + J[f_i_r, x_j_snd] = vm[i] * sum( real(am.matrix[i,k]) * vm[k] * -sin(va[i] - va[k]) - imag(am.matrix[i,k]) * vm[k] * cos(va[i] - va[k]) for k in neighbors[i] if k != i) + + J[f_i_i, x_j_fst] = 2*imag(y_ii)*vm[i] + sum( imag(am.matrix[i,k]) * vm[k] * cos(va[i] - va[k]) + real(am.matrix[i,k]) * vm[k] * sin(va[i] - va[k]) for k in neighbors[i] if k != i) + J[f_i_i, x_j_snd] = vm[i] * sum( imag(am.matrix[i,k]) * vm[k] * -sin(va[i] - va[k]) + real(am.matrix[i,k]) * vm[k] * cos(va[i] - va[k]) for k in neighbors[i] if k != i) + else + y_ij = am.matrix[i,j] + J[f_i_r, x_j_fst] = vm[i] * (real(y_ij) * cos(va[i] - va[j]) - imag(y_ij) * sin(va[i] - va[j])) + J[f_i_r, x_j_snd] = vm[i] * vm[j] * (real(y_ij) * sin(va[i] - va[j]) - imag(y_ij) * -cos(va[i] - va[j])) + + J[f_i_i, x_j_fst] = vm[i] * (imag(y_ij) * cos(va[i] - va[j]) + real(y_ij) * sin(va[i] - va[j])) + J[f_i_i, x_j_snd] = vm[i] * vm[j] * (imag(y_ij) * sin(va[i] - va[j]) + real(y_ij) * -cos(va[i] - va[j])) + end + elseif bus_type == 2 + if i == j + J[f_i_r, x_j_fst] = 0.0 + J[f_i_i, x_j_fst] = 1.0 + + y_ii = am.matrix[i,i] + J[f_i_r, x_j_snd] = vm[i] * sum( real(am.matrix[i,k]) * vm[k] * -sin(va[i] - va[k]) - imag(am.matrix[i,k]) * vm[k] * cos(va[i] - va[k]) for k in neighbors[i] if k != i) + J[f_i_i, x_j_snd] = vm[i] * sum( imag(am.matrix[i,k]) * vm[k] * -sin(va[i] - va[k]) + real(am.matrix[i,k]) * vm[k] * cos(va[i] - va[k]) for k in neighbors[i] if k != i) + else + J[f_i_r, x_j_fst] = 0.0 + J[f_i_i, x_j_fst] = 0.0 + + y_ij = am.matrix[i,j] + J[f_i_r, x_j_snd] = vm[i] * vm[j] * (real(y_ij) * sin(va[i] - va[j]) - imag(y_ij) * -cos(va[i] - va[j])) + J[f_i_i, x_j_snd] = vm[i] * vm[j] * (imag(y_ij) * sin(va[i] - va[j]) + real(y_ij) * -cos(va[i] - va[j])) + end + elseif bus_type == 3 + if i == j + J[f_i_r, x_j_fst] = 1.0 + J[f_i_r, x_j_snd] = 0.0 + J[f_i_i, x_j_fst] = 0.0 + J[f_i_i, x_j_snd] = 1.0 + end + else + @assert false + end + end + end + + return J +end +=# + diff --git a/src/prob/pf.jl b/src/prob/pf.jl index 2dc6289b6..d04365a43 100644 --- a/src/prob/pf.jl +++ b/src/prob/pf.jl @@ -158,10 +158,10 @@ function instantiate_pf_data(data::Dict{String,<:Any}) gen_bus = data["bus"]["$(gen["gen_bus"])"] if gen["gen_status"] != 0 if gen_bus["bus_type"] == 3 - p_delta[gen_bus["index"]] += gen["pg"] - q_delta[gen_bus["index"]] += gen["qg"] + p_delta[gen_bus["index"]] -= gen["pg"] + q_delta[gen_bus["index"]] -= gen["qg"] elseif gen_bus["bus_type"] == 2 - q_delta[gen_bus["index"]] += gen["qg"] + q_delta[gen_bus["index"]] -= gen["qg"] else @assert false end @@ -192,8 +192,8 @@ function instantiate_pf_data(data::Dict{String,<:Any}) bus_type_idx = Int[data["bus"]["$(bus_id)"]["bus_type"] for bus_id in am.idx_to_bus] - p_delta_base_idx = Float64[p_delta[bus_id] for bus_id in am.idx_to_bus] - q_delta_base_idx = Float64[q_delta[bus_id] for bus_id in am.idx_to_bus] + p_delta_base_idx = Float64[-p_delta[bus_id] for bus_id in am.idx_to_bus] + q_delta_base_idx = Float64[-q_delta[bus_id] for bus_id in am.idx_to_bus] p_inject_idx = [0.0 for bus_id in am.idx_to_bus] q_inject_idx = [0.0 for bus_id in am.idx_to_bus] diff --git a/test/data-basic.jl b/test/data-basic.jl new file mode 100644 index 000000000..1ef6102ba --- /dev/null +++ b/test/data-basic.jl @@ -0,0 +1,203 @@ +# Tests of basic data transformation and utilities + + +@testset "test basic network transformation" begin + + @testset "7-bus with inactive components" begin + data = make_basic_network(PowerModels.parse_file("../test/data/matpower/case7_tplgy.m")) + + @test length(data["bus"]) == 4 + @test length(data["load"]) == 2 + @test length(data["shunt"]) == 2 + @test length(data["gen"]) == 1 + @test length(data["branch"]) == 2 + @test length(data["dcline"]) == 0 + @test length(data["switch"]) == 0 + @test length(data["storage"]) == 0 + + result = run_opf(data, ACPPowerModel, ipopt_solver) + @test isapprox(result["objective"], 1036.52; atol=1e0) + end + + @testset "5-bus with switches components" begin + data = make_basic_network(PowerModels.parse_file("../test/data/matpower/case5_sw.m")) + + @test length(data["bus"]) == 4 + @test length(data["load"]) == 3 + @test length(data["shunt"]) == 1 + @test length(data["gen"]) == 5 + @test length(data["branch"]) == 5 + @test length(data["dcline"]) == 0 + @test length(data["switch"]) == 0 + @test length(data["storage"]) == 0 + + result = run_opf(data, ACPPowerModel, ipopt_solver) + @test isapprox(result["objective"], 16551.7; atol=1e0) + end + +end + + +@testset "test basic network functions" begin + + @testset "basic bus injection" begin + data = make_basic_network(PowerModels.parse_file("../test/data/matpower/case14.m")) + + result = run_opf(data, DCPPowerModel, ipopt_solver) + update_data!(data, result["solution"]) + + bi = calc_basic_bus_injection(data) + + @test isapprox(real(sum(bi)), 0.0; atol=1e-6) + + + result = run_opf(data, ACPPowerModel, ipopt_solver) + update_data!(data, result["solution"]) + + bi = calc_basic_bus_injection(data) + + @test isapprox(sum(bi), 0.092872 - 0.0586887im; atol=1e-6) + end + + @testset "basic incidence matrix" begin + data = make_basic_network(PowerModels.parse_file("../test/data/matpower/case5_sw.m")) + + I = calc_basic_incidence_matrix(data) + + @test size(I, 1) == length(data["branch"]) + @test size(I, 2) == length(data["bus"]) + @test sum(I) == 0 + @test sum(abs.(I)) == 2*length(data["branch"]) + end + + @testset "basic admittance matrix" begin + data = make_basic_network(PowerModels.parse_file("../test/data/matpower/case5_sw.m")) + + AM = calc_basic_admittance_matrix(data) + + @test size(AM, 1) == length(data["bus"]) + @test size(AM, 2) == length(data["bus"]) + @test isapprox(sum(AM), 0.0151187 + 0.00624668im; atol=1e-6) + end + + @testset "basic branch series impedance and susceptance matrix" begin + data = make_basic_network(PowerModels.parse_file("../test/data/matpower/case5_sw.m")) + + bz = calc_basic_branch_series_impedance(data) + A = calc_basic_incidence_matrix(data) + + # docs example + Y = imag(LinearAlgebra.Diagonal(inv.(bz))) + SM_1 = A'*Y*A + + SM_2 = calc_basic_susceptance_matrix(data) + + @test isapprox(SM_1, SM_2; atol=1e-6) + + + result = run_opf(data, DCPPowerModel, ipopt_solver) + update_data!(data, result["solution"]) + + va = angle.(calc_basic_bus_voltage(data)) + B = calc_basic_susceptance_matrix(data) + + # docs example + bus_injection = -B * va + @test isapprox(bus_injection, real(calc_basic_bus_injection(data)); atol=1e-6) + end + + @testset "basic bus voltage and branch power matrix" begin + data = make_basic_network(PowerModels.parse_file("../test/data/matpower/case5_sw.m")) + + va = angle.(calc_basic_bus_voltage(data)) + + bz = calc_basic_branch_series_impedance(data) + A = calc_basic_incidence_matrix(data) + + Y = imag(LinearAlgebra.Diagonal(inv.(bz))) + BB_1 = (A'*Y)' + + BB_2 = calc_basic_branch_susceptance_matrix(data) + + @test isapprox(BB_1, BB_2; atol=1e-6) + + # docs example + branch_power = -BB_1*va + + @test length(branch_power) == length(data["branch"]) + + for (i,branch) in data["branch"] + g,b = calc_branch_y(branch) + pf = -b*(va[branch["f_bus"]] - va[branch["t_bus"]]) + @test isapprox(branch_power[branch["index"]], pf; atol=1e-6) + end + println() + end + + @testset "basic dc power flow" begin + data = make_basic_network(PowerModels.parse_file("../test/data/matpower/case5_sw.m")) + solution = compute_dc_pf(data) + + va = compute_basic_dc_pf(data) + + for (i,val) in enumerate(va) + @test isapprox(solution["bus"]["$(i)"]["va"], val; atol=1e-6) + end + end + + @testset "basic ptdf matrix" begin + data = make_basic_network(PowerModels.parse_file("../test/data/matpower/case5_sw.m")) + + P = calc_basic_ptdf_matrix(data) + + @test size(P, 1) == length(data["branch"]) + @test size(P, 2) == length(data["bus"]) + @test isapprox(sum(P), 0.9894736; atol=1e-6) + + + result = run_opf(data, DCPPowerModel, ipopt_solver) + update_data!(data, result["solution"]) + + bi = real(calc_basic_bus_injection(data)) + # accounts for vm = 1.0 assumption + for (i,shunt) in data["shunt"] + if !isapprox(shunt["gs"], 0.0) + bi[shunt["shunt_bus"]] += shunt["gs"] + end + end + + va = angle.(calc_basic_bus_voltage(data)) + + # docs example + branch_power = P*bi + + for (i,branch) in data["branch"] + g,b = calc_branch_y(branch) + pf = -b*(va[branch["f_bus"]] - va[branch["t_bus"]]) + @test isapprox(branch_power[branch["index"]], pf; atol=1e-6) + end + end + + @testset "basic ptdf columns" begin + data = make_basic_network(PowerModels.parse_file("../test/data/matpower/case5_sw.m")) + + P = calc_basic_ptdf_matrix(data) + + for i in 1:length(data["branch"]) + row = calc_basic_ptdf_row(data, i) + @test isapprox(P[i,:], row; atol=1e-6) + end + end + + # @testset "basic jacobian matrix" begin + # data = make_basic_network(PowerModels.parse_file("../test/data/matpower/case5_sw.m")) + + # J = calc_basic_jacobian_matrix(data) + + # @test size(J, 1) == 2*length(data["bus"]) + # @test size(J, 2) == 2*length(data["bus"]) + # @test isapprox(sum(J), 10.0933; atol=1e-4) + # end + +end + diff --git a/test/data.jl b/test/data.jl index 797eb37b9..fa6783595 100644 --- a/test/data.jl +++ b/test/data.jl @@ -803,4 +803,44 @@ end end end -end \ No newline at end of file +end + + + +@testset "test renumber bus ids" begin + @testset "5-bus with dcline" begin + data = PowerModels.parse_file("../test/data/matpower/case5_dc.m") + result_1 = run_opf(data, ACPPowerModel, ipopt_solver) + + update_bus_ids!(data, Dict(1 => 10, 2=>20, 3=>30, 4=>40, 5=>50)) + + result_2 = run_opf(data, ACPPowerModel, ipopt_solver) + + @test isapprox(result_1["objective"], result_2["objective"]; atol=1e-6) + end + + @testset "5-bus with switches" begin + data = PowerModels.parse_file("../test/data/matpower/case5_sw.m") + result_1 = PowerModels._run_opf_sw(data, DCPPowerModel, ipopt_solver) + + update_bus_ids!(data, Dict(1 => 10, 2=>20, 3=>30, 4=>40, 10=>100)) + + result_2 = PowerModels._run_opf_sw(data, DCPPowerModel, ipopt_solver) + + @test isapprox(result_1["objective"], result_2["objective"]; atol=1e-6) + end +end + +@testset "test resolve switches" begin + @testset "5-bus with switches" begin + data = PowerModels.parse_file("../test/data/matpower/case5_sw.m") + resolve_swithces!(data) + + @test length(data["switch"]) == 0 + @test length(data["bus"]) == 4 + + result = run_opf(data, ACPPowerModel, ipopt_solver) + + @test isapprox(result["objective"], 16641.20; atol=1e0) + end +end diff --git a/test/runtests.jl b/test/runtests.jl index ef6efb9ea..a8c984cd9 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -44,6 +44,8 @@ include("common.jl") include("data.jl") + include("data-basic.jl") + include("model.jl") include("am.jl")