diff --git a/Project.toml b/Project.toml index 85dbbb0c35..fd1cf39fb9 100644 --- a/Project.toml +++ b/Project.toml @@ -45,13 +45,13 @@ JuMP = "1" LinearAlgebra = "1" Logging = "1" MathOptInterface = "1" -PowerModels = "~0.20" -PowerNetworkMatrices = "^0.9" +PowerModels = "^0.20" +PowerNetworkMatrices = "^0.10" PowerSystems = "^3" PrettyTables = "2" ProgressMeter = "^1.5" SHA = "0.7" Serialization = "1" -TimeSeries = "~0.23" +TimeSeries = "~0.23, ~0.24" TimerOutputs = "~0.5" julia = "^1.6" diff --git a/src/core/definitions.jl b/src/core/definitions.jl index b51e41e7db..142967f063 100644 --- a/src/core/definitions.jl +++ b/src/core/definitions.jl @@ -53,9 +53,6 @@ const KiB = 1024 const MiB = KiB * KiB const GiB = MiB * KiB -const UNSUPPORTED_POWERMODELS = - [PM.SOCBFPowerModel, PM.SOCBFConicPowerModel, PM.IVRPowerModel] - const PSI_NAME_DELIMITER = "__" const M_VALUE = 1e6 diff --git a/src/core/network_model.jl b/src/core/network_model.jl index 058d512f77..caa2d8873f 100644 --- a/src/core/network_model.jl +++ b/src/core/network_model.jl @@ -20,7 +20,7 @@ Establishes the model for a particular device specified by type. - `use_slacks::Bool`: Adds slacks to the network modelings - `PTDF::PTDF`: PTDF Array calculated using PowerNetworkMatrices - `duals::Vector{DataType}`: Constraint types to calculate the duals - + - `reduce_radial_branches::Bool`: Skips modeling radial branches in the system to reduce problem size # Example ptdf_array = PTDF(system) @@ -32,21 +32,34 @@ mutable struct NetworkModel{T <: PM.AbstractPowerModel} subnetworks::Dict{Int, Set{Int}} bus_area_map::Dict{PSY.ACBus, Int} duals::Vector{DataType} + radial_network_reduction::PNM.RadialNetworkReduction + reduce_radial_branches::Bool function NetworkModel( ::Type{T}; use_slacks = false, PTDF_matrix = nothing, + reduce_radial_branches = false, subnetworks = Dict{Int, Set{Int}}(), duals = Vector{DataType}(), ) where {T <: PM.AbstractPowerModel} _check_pm_formulation(T) - new{T}(use_slacks, PTDF_matrix, subnetworks, Dict{PSY.ACBus, Int}(), duals) + new{T}( + use_slacks, + PTDF_matrix, + subnetworks, + Dict{PSY.ACBus, Int}(), + duals, + PNM.RadialNetworkReduction(), + reduce_radial_branches, + ) end end get_use_slacks(m::NetworkModel) = m.use_slacks get_PTDF_matrix(m::NetworkModel) = m.PTDF_matrix +get_reduce_radial_branches(m::NetworkModel) = m.reduce_radial_branches +get_radial_network_reduction(m::NetworkModel) = m.radial_network_reduction get_duals(m::NetworkModel) = m.duals get_network_formulation(::NetworkModel{T}) where {T} = T get_reference_buses(m::NetworkModel{T}) where {T <: PM.AbstractPowerModel} = @@ -62,6 +75,15 @@ function add_dual!(model::NetworkModel, dual) return end +function check_radial_branch_reduction_compatibility( + ::Type{T}, +) where {T <: PM.AbstractPowerModel} + if T ∈ INCOMPATIBLE_WITH_RADIAL_BRANCHES_POWERMODELS + error("Network Model $T is not compatible with radial branch reduction") + end + return +end + function instantiate_network_model( model::NetworkModel{T}, sys::PSY.System, @@ -69,6 +91,10 @@ function instantiate_network_model( if isempty(model.subnetworks) model.subnetworks = PNM.find_subnetworks(sys) end + if model.reduce_radial_branches + check_radial_branch_reduction_compatibility(T) + model.radial_network_reduction = PNM.RadialNetworkReduction(sys) + end return end @@ -90,7 +116,12 @@ end function instantiate_network_model(model::NetworkModel{StandardPTDFModel}, sys::PSY.System) if get_PTDF_matrix(model) === nothing @info "PTDF Matrix not provided. Calculating using PowerNetworkMatrices.PTDF" - model.PTDF_matrix = PNM.PTDF(sys) + model.PTDF_matrix = + PNM.PTDF(sys; reduce_radial_branches = model.reduce_radial_branches) + end + if model.reduce_radial_branches + @assert !isempty(model.PTDF_matrix.radial_network_reduction) + model.radial_network_reduction = model.PTDF_matrix.radial_network_reduction end get_PTDF_matrix(model).subnetworks model.subnetworks = deepcopy(get_PTDF_matrix(model).subnetworks) @@ -107,19 +138,27 @@ function _assign_subnetworks_to_buses( ) where {T <: Union{CopperPlatePowerModel, StandardPTDFModel}} subnetworks = model.subnetworks temp_bus_map = Dict{Int, Int}() - for bus in get_available_components(PSY.ACBus, sys) + radial_network_reduction = PSI.get_radial_network_reduction(model) + for bus in PSI.get_available_components(PSY.ACBus, sys) bus_no = PSY.get_number(bus) + mapped_bus_no = PNM.get_mapped_bus_number(radial_network_reduction, bus) if haskey(temp_bus_map, bus_no) model.bus_area_map[bus] = temp_bus_map[bus_no] + continue else + bus_mapped = false for (subnet, bus_set) in subnetworks - if bus_no ∈ bus_set + if mapped_bus_no ∈ bus_set temp_bus_map[bus_no] = subnet model.bus_area_map[bus] = subnet + bus_mapped = true break end end end + if !bus_mapped + error("Bus $(PSY.summary(bus)) not mapped to any reference bus") + end end return end diff --git a/src/core/optimization_container.jl b/src/core/optimization_container.jl index be85ae1a9c..9686a1d06d 100644 --- a/src/core/optimization_container.jl +++ b/src/core/optimization_container.jl @@ -431,9 +431,14 @@ function _make_system_expressions!( subnetworks::Dict{Int, Set{Int}}, dc_bus_numbers::Vector{Int}, ::Type{<:PM.AbstractPowerModel}, + bus_reduction_map::Dict{Int64, Set{Int64}}, ) time_steps = get_time_steps(container) - ac_bus_numbers = collect(Iterators.flatten(values(subnetworks))) + if isempty(bus_reduction_map) + ac_bus_numbers = collect(Iterators.flatten(values(subnetworks))) + else + ac_bus_numbers = collect(keys(bus_reduction_map)) + end container.expressions = Dict( ExpressionKey(ActivePowerBalance, PSY.ACBus) => _make_container_array(ac_bus_numbers, time_steps), @@ -450,9 +455,14 @@ function _make_system_expressions!( subnetworks::Dict{Int, Set{Int}}, dc_bus_numbers::Vector{Int}, ::Type{<:PM.AbstractActivePowerModel}, + bus_reduction_map::Dict{Int64, Set{Int64}}, ) time_steps = get_time_steps(container) - ac_bus_numbers = collect(Iterators.flatten(values(subnetworks))) + if isempty(bus_reduction_map) + ac_bus_numbers = collect(Iterators.flatten(values(subnetworks))) + else + ac_bus_numbers = collect(keys(bus_reduction_map)) + end container.expressions = Dict( ExpressionKey(ActivePowerBalance, PSY.ACBus) => _make_container_array(ac_bus_numbers, time_steps), @@ -467,6 +477,7 @@ function _make_system_expressions!( subnetworks::Dict{Int, Set{Int}}, ::Vector{Int}, ::Type{CopperPlatePowerModel}, + bus_reduction_map::Dict{Int64, Set{Int64}}, ) time_steps = get_time_steps(container) subnetworks_ref_buses = collect(keys(subnetworks)) @@ -482,9 +493,14 @@ function _make_system_expressions!( subnetworks::Dict{Int, Set{Int}}, dc_bus_numbers::Vector{Int}, ::Type{T}, + bus_reduction_map::Dict{Int64, Set{Int64}}, ) where {T <: Union{PTDFPowerModel, StandardPTDFModel}} time_steps = get_time_steps(container) - ac_bus_numbers = sort!(collect(Iterators.flatten(values(subnetworks)))) + if isempty(bus_reduction_map) + ac_bus_numbers = collect(Iterators.flatten(values(subnetworks))) + else + ac_bus_numbers = collect(keys(bus_reduction_map)) + end subnetworks = collect(keys(subnetworks)) container.expressions = Dict( ExpressionKey(ActivePowerBalance, PSY.System) => @@ -492,7 +508,9 @@ function _make_system_expressions!( ExpressionKey(ActivePowerBalance, PSY.DCBus) => _make_container_array(dc_bus_numbers, time_steps), ExpressionKey(ActivePowerBalance, PSY.ACBus) => - _make_container_array(ac_bus_numbers, time_steps), + # Bus numbers are sorted to guarantee consistency in the order between the + # containers + _make_container_array(sort!(ac_bus_numbers), time_steps), ) return end @@ -502,9 +520,10 @@ function initialize_system_expressions!( ::Type{T}, subnetworks::Dict{Int, Set{Int}}, system::PSY.System, + bus_reduction_map::Dict{Int64, Set{Int64}}, ) where {T <: PM.AbstractPowerModel} dc_bus_numbers = [PSY.get_number(b) for b in PSY.get_components(PSY.DCBus, system)] - _make_system_expressions!(container, subnetworks, dc_bus_numbers, T) + _make_system_expressions!(container, subnetworks, dc_bus_numbers, T, bus_reduction_map) return end @@ -520,7 +539,7 @@ function build_impl!( transmission, transmission_model.subnetworks, sys, - ) + transmission_model.radial_network_reduction.bus_reduction_map) # Order is required for device_model in values(template.devices) diff --git a/src/devices_models/device_constructors/branch_constructor.jl b/src/devices_models/device_constructors/branch_constructor.jl index ec1b7d6d9d..616500cafb 100644 --- a/src/devices_models/device_constructors/branch_constructor.jl +++ b/src/devices_models/device_constructors/branch_constructor.jl @@ -106,7 +106,7 @@ function construct_device!( container, devices, device_model, - get_network_formulation(network_model), + network_model, ) end return @@ -273,7 +273,7 @@ function construct_device!( container, devices, device_model, - get_network_formulation(network_model), + network_model, ) add_constraint_dual!(container, sys, device_model) return @@ -331,7 +331,7 @@ function construct_device!( ) where {T <: PSY.ACBranch} devices = get_available_components(T, sys, get_attribute(model, "filter_function")) - branch_rate_bounds!(container, devices, model, get_network_formulation(network_model)) + branch_rate_bounds!(container, devices, model, network_model) add_constraints!(container, RateLimitConstraintFromTo, devices, model, network_model) add_constraints!(container, RateLimitConstraintToFrom, devices, model, network_model) @@ -356,7 +356,7 @@ function construct_device!( ) where {T <: PSY.ACBranch} devices = get_available_components(T, sys, get_attribute(model, "filter_function")) - branch_rate_bounds!(container, devices, model, get_network_formulation(network_model)) + branch_rate_bounds!(container, devices, model, network_model) add_constraint_dual!(container, sys, model) return end diff --git a/src/devices_models/devices/AC_branches.jl b/src/devices_models/devices/AC_branches.jl index 1176f83700..e14215cec1 100644 --- a/src/devices_models/devices/AC_branches.jl +++ b/src/devices_models/devices/AC_branches.jl @@ -104,11 +104,18 @@ function branch_rate_bounds!( container::OptimizationContainer, devices::IS.FlattenIteratorWrapper{B}, ::DeviceModel{B, <:AbstractBranchFormulation}, - ::Type{<:PM.AbstractDCPModel}, + network_model::NetworkModel{<:PM.AbstractDCPModel}, ) where {B <: PSY.ACBranch} var = get_variable(container, FlowActivePowerVariable(), B) + + radial_network_reduction = get_radial_network_reduction(network_model) + radial_branches_names = PNM.get_radial_branches(radial_network_reduction) + for d in devices name = PSY.get_name(d) + if name ∈ radial_branches_names + continue + end for t in get_time_steps(container) JuMP.set_upper_bound(var[name, t], PSY.get_rate(d)) JuMP.set_lower_bound(var[name, t], -1.0 * PSY.get_rate(d)) @@ -121,15 +128,23 @@ function branch_rate_bounds!( container::OptimizationContainer, devices::IS.FlattenIteratorWrapper{B}, ::DeviceModel{B, <:AbstractBranchFormulation}, - ::Type{<:PM.AbstractPowerModel}, + network_model::NetworkModel{<:PM.AbstractPowerModel}, ) where {B <: PSY.ACBranch} vars = [ get_variable(container, FlowActivePowerFromToVariable(), B), get_variable(container, FlowActivePowerToFromVariable(), B), ] + + time_steps = get_time_steps(container) + radial_network_reduction = get_radial_network_reduction(network_model) + radial_branches_names = PNM.get_radial_branches(radial_network_reduction) + for d in devices name = PSY.get_name(d) - for t in get_time_steps(container), var in vars + if name ∈ radial_branches_names + continue + end + for t in time_steps, var in vars JuMP.set_upper_bound(var[name, t], PSY.get_rate(d)) JuMP.set_lower_bound(var[name, t], -1.0 * PSY.get_rate(d)) end @@ -168,14 +183,55 @@ function add_constraints!( container::OptimizationContainer, cons_type::Type{RateLimitConstraint}, devices::IS.FlattenIteratorWrapper{T}, - model::DeviceModel{T, U}, - ::NetworkModel{X}, + ::DeviceModel{T, U}, + network_model::NetworkModel{V}, ) where { T <: PSY.ACBranch, U <: AbstractBranchFormulation, - X <: PM.AbstractActivePowerModel, + V <: PM.AbstractActivePowerModel, } - add_range_constraints!(container, cons_type, FlowActivePowerVariable, devices, model, X) + time_steps = get_time_steps(container) + radial_network_reduction = get_radial_network_reduction(network_model) + if isempty(radial_network_reduction) + device_names = [PSY.get_name(d) for d in devices] + else + device_names = PNM.get_meshed_branches(radial_network_reduction) + end + + con_lb = + add_constraints_container!( + container, + cons_type(), + T, + device_names, + time_steps; + meta = "lb", + ) + con_ub = + add_constraints_container!( + container, + cons_type(), + T, + device_names, + time_steps; + meta = "ub", + ) + + array = get_variable(container, FlowActivePowerVariable(), T) + + for device in devices + ci_name = PSY.get_name(device) + if ci_name ∈ PNM.get_radial_branches(radial_network_reduction) + continue + end + limits = get_min_max_limits(device, RateLimitConstraint, U) # depends on constraint type and formulation type + for t in time_steps + con_ub[ci_name, t] = + JuMP.@constraint(container.JuMPmodel, array[ci_name, t] <= limits.max) + con_lb[ci_name, t] = + JuMP.@constraint(container.JuMPmodel, array[ci_name, t] >= limits.min) + end + end return end @@ -215,7 +271,7 @@ function add_constraints!( cons_type::Type{RateLimitConstraintFromTo}, devices::IS.FlattenIteratorWrapper{B}, ::DeviceModel{B, <:AbstractBranchFormulation}, - ::NetworkModel{T}, + network_model::NetworkModel{T}, ) where {B <: PSY.ACBranch, T <: PM.AbstractPowerModel} rating_data = [(PSY.get_name(h), PSY.get_rate(h)) for h in devices] @@ -231,7 +287,13 @@ function add_constraints!( ) constraint = get_constraint(container, cons_type(), B) + radial_network_reduction = get_radial_network_reduction(network_model) + radial_branches_names = PNM.get_radial_branches(radial_network_reduction) + for r in rating_data + if r[1] ∈ radial_branches_names + continue + end for t in time_steps constraint[r[1], t] = JuMP.@constraint( get_jump_model(container), @@ -250,7 +312,7 @@ function add_constraints!( cons_type::Type{RateLimitConstraintToFrom}, devices::IS.FlattenIteratorWrapper{B}, ::DeviceModel{B, <:AbstractBranchFormulation}, - ::NetworkModel{T}, + network_model::NetworkModel{T}, ) where {B <: PSY.ACBranch, T <: PM.AbstractPowerModel} rating_data = [(PSY.get_name(h), PSY.get_rate(h)) for h in devices] @@ -266,7 +328,13 @@ function add_constraints!( ) constraint = get_constraint(container, cons_type(), B) + radial_network_reduction = get_radial_network_reduction(network_model) + radial_branches_names = PNM.get_radial_branches(radial_network_reduction) + for r in rating_data + if r[1] ∈ radial_branches_names + continue + end for t in time_steps constraint[r[1], t] = JuMP.@constraint( get_jump_model(container), diff --git a/src/devices_models/devices/common/add_to_expression.jl b/src/devices_models/devices/common/add_to_expression.jl index 54aa14a590..c9c7e1cb68 100644 --- a/src/devices_models/devices/common/add_to_expression.jl +++ b/src/devices_models/devices/common/add_to_expression.jl @@ -79,7 +79,7 @@ function add_to_expression!( ::Type{U}, devices::IS.FlattenIteratorWrapper{V}, model::DeviceModel{V, W}, - ::NetworkModel{X}, + network_model::NetworkModel{X}, ) where { T <: SystemBalanceExpressions, U <: TimeSeriesParameter, @@ -89,11 +89,12 @@ function add_to_expression!( } param_container = get_parameter(container, U(), V) multiplier = get_multiplier_array(param_container) + radial_network_reduction = get_radial_network_reduction(network_model) for d in devices, t in get_time_steps(container) - bus_number = PSY.get_number(PSY.get_bus(d)) + bus_no = PNM.get_mapped_bus_number(radial_network_reduction, PSY.get_bus(d)) name = PSY.get_name(d) _add_to_jump_expression!( - get_expression(container, T(), PSY.ACBus)[bus_number, t], + get_expression(container, T(), PSY.ACBus)[bus_no, t], get_parameter_column_refs(param_container, name)[t], multiplier[name, t], ) @@ -107,7 +108,7 @@ function add_to_expression!( ::Type{U}, devices::IS.FlattenIteratorWrapper{V}, ::DeviceModel{V, W}, - ::NetworkModel{X}, + network_model::NetworkModel{X}, ) where { T <: ActivePowerBalance, U <: OnStatusParameter, @@ -116,13 +117,13 @@ function add_to_expression!( X <: PM.AbstractPowerModel, } parameter = get_parameter_array(container, U(), V) - + radial_network_reduction = get_radial_network_reduction(network_model) for d in devices, t in get_time_steps(container) - bus_number = PSY.get_number(PSY.get_bus(d)) + bus_no = PNM.get_mapped_bus_number(radial_network_reduction, PSY.get_bus(d)) name = PSY.get_name(d) mult = get_expression_multiplier(U(), T(), d, W()) _add_to_jump_expression!( - get_expression(container, T(), PSY.ACBus)[bus_number, t], + get_expression(container, T(), PSY.ACBus)[bus_no, t], parameter[name, t], mult, ) @@ -139,7 +140,7 @@ function add_to_expression!( ::Type{U}, devices::IS.FlattenIteratorWrapper{V}, ::DeviceModel{V, W}, - ::NetworkModel{X}, + network_model::NetworkModel{X}, ) where { T <: SystemBalanceExpressions, U <: VariableType, @@ -149,11 +150,12 @@ function add_to_expression!( } variable = get_variable(container, U(), V) expression = get_expression(container, T(), PSY.ACBus) + radial_network_reduction = get_radial_network_reduction(network_model) for d in devices, t in get_time_steps(container) name = PSY.get_name(d) - bus_number = PSY.get_number(PSY.get_bus(d)) + bus_no = PNM.get_mapped_bus_number(radial_network_reduction, PSY.get_bus(d)) _add_to_jump_expression!( - expression[bus_number, t], + expression[bus_no, t], variable[name, t], get_variable_multiplier(U(), V, W()), ) @@ -215,8 +217,9 @@ function add_to_expression!( var = get_variable(container, U(), V) nodal_expr = get_expression(container, T(), PSY.ACBus) sys_expr = get_expression(container, T(), PSY.System) + radial_network_reduction = get_radial_network_reduction(network_model) for d in devices - bus_no_to = PSY.get_number(PSY.get_arc(d).to) + bus_no_to = PNM.get_mapped_bus_number(radial_network_reduction, PSY.get_arc(d).to) ref_bus_from = get_reference_bus(network_model, PSY.get_arc(d).from) ref_bus_to = get_reference_bus(network_model, PSY.get_arc(d).to) for t in get_time_steps(container) @@ -250,8 +253,10 @@ function add_to_expression!( var = get_variable(container, U(), V) nodal_expr = get_expression(container, T(), PSY.ACBus) sys_expr = get_expression(container, T(), PSY.System) + radial_network_reduction = get_radial_network_reduction(network_model) for d in devices - bus_no_from = PSY.get_number(PSY.get_arc(d).from) + bus_no_from = + PNM.get_mapped_bus_number(radial_network_reduction, PSY.get_arc(d).from) ref_bus_to = get_reference_bus(network_model, PSY.get_arc(d).to) ref_bus_from = get_reference_bus(network_model, PSY.get_arc(d).from) for t in get_time_steps(container) @@ -352,12 +357,14 @@ function add_to_expression!( } variable = get_variable(container, U(), V) expression = get_expression(container, T(), PSY.ACBus) + radial_network_reduction = get_radial_network_reduction(network_model) for d in devices name = PSY.get_name(d) - bus_number = PSY.get_number(PSY.get_arc(d).from) + bus_no_ = PSY.get_number(PSY.get_arc(d).from) + bus_no = PNM.get_mapped_bus_number(radial_network_reduction, bus_no_) for t in get_time_steps(container) _add_to_jump_expression!( - expression[bus_number, t], + expression[bus_no, t], variable[name, t], get_variable_multiplier(U(), V, W()), ) @@ -375,7 +382,7 @@ function add_to_expression!( ::Type{U}, devices::IS.FlattenIteratorWrapper{V}, ::DeviceModel{V, W}, - ::NetworkModel{X}, + network_model::NetworkModel{X}, ) where { T <: ActivePowerBalance, U <: FlowActivePowerToFromVariable, @@ -385,12 +392,14 @@ function add_to_expression!( } variable = get_variable(container, U(), V) expression = get_expression(container, T(), PSY.ACBus) + radial_network_reduction = get_radial_network_reduction(network_model) for d in devices name = PSY.get_name(d) - bus_number = PSY.get_number(PSY.get_arc(d).to) + bus_no_ = PSY.get_number(PSY.get_arc(d).to) + bus_no = PNM.get_mapped_bus_number(radial_network_reduction, bus_no_) for t in get_time_steps(container) _add_to_jump_expression!( - expression[bus_number, t], + expression[bus_no, t], variable[name, t], get_variable_multiplier(U(), V, W()), ) @@ -405,7 +414,7 @@ function add_to_expression!( ::Type{U}, devices::IS.FlattenIteratorWrapper{V}, ::DeviceModel{V, W}, - ::NetworkModel{X}, + network_model::NetworkModel{X}, ) where { T <: SystemBalanceExpressions, U <: OnVariable, @@ -415,11 +424,13 @@ function add_to_expression!( } variable = get_variable(container, U(), V) expression = get_expression(container, T(), PSY.ACBus) + radial_network_reduction = get_radial_network_reduction(network_model) for d in devices, t in get_time_steps(container) name = PSY.get_name(d) - bus_number = PSY.get_number(PSY.get_bus(d)) + bus_no_ = PSY.get_number(PSY.get_bus(d)) + bus_no = PNM.get_mapped_bus_number(radial_network_reduction, bus_no_) _add_to_jump_expression!( - expression[bus_number, t], + expression[bus_no, t], variable[name, t], get_variable_multiplier(U(), d, W()), ) @@ -574,10 +585,12 @@ function add_to_expression!( multiplier = get_multiplier_array(param_container) sys_expr = get_expression(container, T(), PSY.System) nodal_expr = get_expression(container, T(), PSY.ACBus) + radial_network_reduction = get_radial_network_reduction(network_model) for d in devices name = PSY.get_name(d) device_bus = PSY.get_bus(d) - bus_no = PSY.get_number(device_bus) + bus_no_ = PSY.get_number(device_bus) + bus_no = PNM.get_mapped_bus_number(radial_network_reduction, bus_no_) ref_bus = get_reference_bus(network_model, device_bus) param = get_parameter_column_refs(param_container, name) for t in get_time_steps(container) @@ -605,9 +618,11 @@ function add_to_expression!( parameter = get_parameter_array(container, U(), V) sys_expr = get_expression(container, T(), PSY.System) nodal_expr = get_expression(container, T(), PSY.ACBus) + radial_network_reduction = get_radial_network_reduction(network_model) for d in devices, t in get_time_steps(container) name = PSY.get_name(d) - bus_no = PSY.get_number(PSY.get_bus(d)) + bus_no_ = PSY.get_number(PSY.get_bus(d)) + bus_no = PNM.get_mapped_bus_number(radial_network_reduction, bus_no_) mult = get_expression_multiplier(U(), T(), d, W()) device_bus = PSY.get_bus(d) ref_bus = get_reference_bus(network_model, device_bus) @@ -637,10 +652,11 @@ function add_to_expression!( variable = get_variable(container, U(), V) sys_expr = get_expression(container, T(), PSY.System) nodal_expr = get_expression(container, T(), PSY.ACBus) + radial_network_reduction = get_radial_network_reduction(network_model) for d in devices name = PSY.get_name(d) device_bus = PSY.get_bus(d) - bus_no = PSY.get_number(device_bus) + bus_no = PNM.get_mapped_bus_number(radial_network_reduction, device_bus) ref_bus = get_reference_bus(network_model, device_bus) for t in get_time_steps(container) _add_to_jump_expression!( @@ -675,11 +691,11 @@ function add_to_expression!( variable = get_variable(container, U(), V) sys_expr = get_expression(container, T(), PSY.System) nodal_expr = get_expression(container, T(), PSY.ACBus) + radial_network_reduction = get_radial_network_reduction(network_model) for d in devices name = PSY.get_name(d) - device_bus = PSY.get_bus(d) - bus_no = PSY.get_number(device_bus) - ref_bus = get_reference_bus(network_model, device_bus) + bus_no = PNM.get_mapped_bus_number(radial_network_reduction, PSY.get_bus(d)) + ref_bus = get_reference_bus(network_model, PSY.get_bus(d)) for t in get_time_steps(container) _add_to_jump_expression!( sys_expr[ref_bus, t], @@ -715,16 +731,20 @@ function add_to_expression!( } var = get_variable(container, U(), V) expression = get_expression(container, T(), PSY.ACBus) + radial_network_reduction = get_radial_network_reduction(network_model) for d in devices + bus_no_from = + PNM.get_mapped_bus_number(radial_network_reduction, PSY.get_arc(d).from) + bus_no_to = PNM.get_mapped_bus_number(radial_network_reduction, PSY.get_arc(d).to) for t in get_time_steps(container) flow_variable = var[PSY.get_name(d), t] _add_to_jump_expression!( - expression[PSY.get_number(PSY.get_arc(d).from), t], + expression[bus_no_from, t], flow_variable, -1.0, ) _add_to_jump_expression!( - expression[PSY.get_number(PSY.get_arc(d).to), t], + expression[bus_no_to, t], flow_variable, 1.0, ) @@ -752,9 +772,11 @@ function add_to_expression!( var = get_variable(container, U(), V) nodal_expr = get_expression(container, T(), PSY.ACBus) sys_expr = get_expression(container, T(), PSY.System) + radial_network_reduction = get_radial_network_reduction(network_model) for d in devices - bus_no_from = PSY.get_number(PSY.get_arc(d).from) - bus_no_to = PSY.get_number(PSY.get_arc(d).to) + bus_no_from = + PNM.get_mapped_bus_number(radial_network_reduction, PSY.get_arc(d).from) + bus_no_to = PNM.get_mapped_bus_number(radial_network_reduction, PSY.get_arc(d).to) ref_bus_from = get_reference_bus(network_model, PSY.get_arc(d).from) ref_bus_to = get_reference_bus(network_model, PSY.get_arc(d).to) for t in get_time_steps(container) @@ -826,16 +848,20 @@ function add_to_expression!( ) where {T <: ActivePowerBalance, U <: PhaseShifterAngle, V <: PhaseAngleControl} var = get_variable(container, U(), PSY.PhaseShiftingTransformer) expression = get_expression(container, T(), PSY.ACBus) + radial_network_reduction = get_radial_network_reduction(network_model) for d in devices + bus_no_from = + PNM.get_mapped_bus_number(radial_network_reduction, PSY.get_arc(d).from) + bus_no_to = PNM.get_mapped_bus_number(radial_network_reduction, PSY.get_arc(d).to) for t in get_time_steps(container) flow_variable = var[PSY.get_name(d), t] _add_to_jump_expression!( - expression[PSY.get_number(PSY.get_arc(d).from), t], + expression[bus_no_from, t], flow_variable, -get_variable_multiplier(U(), d, V()), ) _add_to_jump_expression!( - expression[PSY.get_number(PSY.get_arc(d).to), t], + expression[bus_no_to, t], flow_variable, get_variable_multiplier(U(), d, V()), ) @@ -1079,7 +1105,8 @@ function add_to_expression!( variable = get_variable(container, U(), PSY.ACBus) expression = get_expression(container, T(), PSY.ACBus) @assert_op length(axes(variable, 1)) == length(axes(expression, 1)) - for t in get_time_steps(container), n in axes(variable, 1) + # We uses axis here to avoid double addition of the slacks to the aggregated buses + for t in get_time_steps(container), n in axes(expression, 1) _add_to_jump_expression!( expression[n, t], variable[n, t], @@ -1102,8 +1129,8 @@ function add_to_expression!( } variable = get_variable(container, U(), PSY.ACBus, "P") expression = get_expression(container, T(), PSY.ACBus) - bus_numbers = PSY.get_number.(get_available_components(PSY.ACBus, sys)) - for t in get_time_steps(container), n in bus_numbers + # We uses axis here to avoid double addition of the slacks to the aggregated buses + for t in get_time_steps(container), n in axes(expression, 1) _add_to_jump_expression!( expression[n, t], variable[n, t], @@ -1126,8 +1153,8 @@ function add_to_expression!( } variable = get_variable(container, U(), PSY.ACBus, "Q") expression = get_expression(container, T(), PSY.ACBus) - bus_numbers = PSY.get_number.(get_available_components(PSY.ACBus, sys)) - for t in get_time_steps(container), n in bus_numbers + # We uses axis here to avoid double addition of the slacks to the aggregated buses + for t in get_time_steps(container), n in axes(expression, 1) _add_to_jump_expression!( expression[n, t], variable[n, t], diff --git a/src/initial_conditions/initialization.jl b/src/initial_conditions/initialization.jl index 80ff0f1fd7..e5222f6bfc 100644 --- a/src/initial_conditions/initialization.jl +++ b/src/initial_conditions/initialization.jl @@ -5,7 +5,12 @@ function get_initial_conditions_template(model::OperationModel) get_network_formulation(model.template); use_slacks = get_use_slacks(get_network_model(model.template)), PTDF_matrix = get_PTDF_matrix(get_network_model(model.template)), + reduce_radial_branches = get_reduce_radial_branches( + get_network_model(model.template), + ), ) + network_model.radial_network_reduction = + get_radial_network_reduction(get_network_model(model.template)) network_model.subnetworks = get_subnetworks(get_network_model(model.template)) bus_area_map = get_bus_area_map(get_network_model(model.template)) if !isempty(bus_area_map) diff --git a/src/network_models/network_constructor.jl b/src/network_models/network_constructor.jl index e2d4a5d850..7f53454439 100644 --- a/src/network_models/network_constructor.jl +++ b/src/network_models/network_constructor.jl @@ -202,15 +202,14 @@ function construct_network!( end if get_use_slacks(model) - add_variables!(container, SystemBalanceSlackUp, sys, T) - add_variables!(container, SystemBalanceSlackDown, sys, T) + add_variables!(container, SystemBalanceSlackUp, sys, model) + add_variables!(container, SystemBalanceSlackDown, sys, model) add_to_expression!( container, ActivePowerBalance, SystemBalanceSlackUp, sys, model, - T, ) add_to_expression!( container, @@ -218,7 +217,6 @@ function construct_network!( SystemBalanceSlackDown, sys, model, - T, ) add_to_expression!( container, @@ -226,7 +224,6 @@ function construct_network!( SystemBalanceSlackUp, sys, model, - T, ) add_to_expression!( container, @@ -234,7 +231,6 @@ function construct_network!( SystemBalanceSlackDown, sys, model, - T, ) objective_function!(container, PSY.ACBus, model) end diff --git a/src/network_models/network_slack_variables.jl b/src/network_models/network_slack_variables.jl index 619aeab047..82e1638e49 100644 --- a/src/network_models/network_slack_variables.jl +++ b/src/network_models/network_slack_variables.jl @@ -37,7 +37,13 @@ function add_variables!( U <: PM.AbstractActivePowerModel, } time_steps = get_time_steps(container) - bus_numbers = PSY.get_number.(get_available_components(PSY.ACBus, sys)) + radial_network_reduction = get_radial_network_reduction(network_model) + if isempty(radial_network_reduction) + bus_numbers = PSY.get_number.(get_available_components(PSY.ACBus, sys)) + else + bus_numbers = collect(keys(PNM.get_bus_reduction_map(radial_network_reduction))) + end + variable = add_variable_container!(container, T(), PSY.ACBus, bus_numbers, time_steps) for t in time_steps, n in bus_numbers variable[n, t] = JuMP.@variable( @@ -59,7 +65,12 @@ function add_variables!( U <: PM.AbstractPowerModel, } time_steps = get_time_steps(container) - bus_numbers = PSY.get_number.(get_available_components(PSY.ACBus, sys)) + radial_network_reduction = get_radial_network_reduction(network_model) + if isempty(radial_network_reduction) + bus_numbers = PSY.get_number.(get_available_components(PSY.ACBus, sys)) + else + bus_numbers = collect(keys(PNM.get_bus_reduction_map(radial_network_reduction))) + end variable_active = add_variable_container!(container, T(), PSY.ACBus, "P", bus_numbers, time_steps) variable_reactive = diff --git a/src/network_models/pm_translator.jl b/src/network_models/pm_translator.jl index 627a5b7cf0..ee86e0544b 100644 --- a/src/network_models/pm_translator.jl +++ b/src/network_models/pm_translator.jl @@ -381,7 +381,7 @@ end function get_branches_to_pm( sys::PSY.System, - ::Type{S}, + network_model::NetworkModel{S}, ::Type{T}, branch_template::BranchModelContainer, start_idx = 0, @@ -389,6 +389,8 @@ function get_branches_to_pm( PM_branches = Dict{String, Any}() PMmap_br = Dict{PM_MAP_TUPLE, T}() + radial_network_reduction = get_radial_network_reduction(network_model) + radial_branches_names = PNM.get_radial_branches(radial_network_reduction) for (d, device_model) in branch_template comp_type = get_component_type(device_model) if comp_type <: TwoTerminalHVDCTypes @@ -398,6 +400,10 @@ function get_branches_to_pm( start_idx += length(PM_branches) filter_func = get_attribute(device_model, "filter_function") for (i, branch) in enumerate(get_available_components(comp_type, sys, filter_func)) + if PSY.get_name(branch) ∈ radial_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) @@ -413,7 +419,7 @@ end function get_branches_to_pm( sys::PSY.System, - ::Type{S}, + network_model::NetworkModel{S}, ::Type{T}, branch_template::BranchModelContainer, start_idx = 0, @@ -442,7 +448,7 @@ end function get_branches_to_pm( ::PSY.System, - ::Type{PTDFPowerModel}, + network_model::NetworkModel{PTDFPowerModel}, ::Type{T}, branch_template::BranchModelContainer, start_idx = 0, @@ -485,13 +491,13 @@ end function pass_to_pm(sys::PSY.System, template::ProblemTemplate, time_periods::Int) ac_lines, PMmap_ac = get_branches_to_pm( sys, - get_network_formulation(template), + get_network_model(template), PSY.ACBranch, template.branches, ) two_terminal_dc_lines, PMmap_dc = get_branches_to_pm( sys, - get_network_formulation(template), + get_network_model(template), TwoTerminalHVDCTypes, template.branches, length(ac_lines), diff --git a/src/network_models/powermodels_interface.jl b/src/network_models/powermodels_interface.jl index cd5f059944..ea0b6bb991 100644 --- a/src/network_models/powermodels_interface.jl +++ b/src/network_models/powermodels_interface.jl @@ -6,6 +6,14 @@ ################################################################################# # Model Definitions +const UNSUPPORTED_POWERMODELS = + [PM.SOCBFPowerModel, PM.SOCBFConicPowerModel, PM.IVRPowerModel] + +const INCOMPATIBLE_WITH_RADIAL_BRANCHES_POWERMODELS = [ + PM.SparseSDPWRMPowerModel, + PTDFPowerModel, +] + function instantiate_nip_expr_model(data::Dict{String, Any}, model_constructor; kwargs...) return PM.instantiate_model(data, model_constructor, instantiate_nip_expr; kwargs...) end @@ -288,17 +296,25 @@ function powermodels_network!( ) where {S <: PM.AbstractPowerModel} time_steps = get_time_steps(container) pm_data, PM_map = pass_to_pm(sys, template, time_steps[end]) - buses = get_available_components(PSY.ACBus, sys) - for t in time_steps, bus in buses - pm_data["nw"]["$(t)"]["bus"]["$(bus.number)"]["inj_p"] = + network_model = get_network_model(template) + radial_network_reduction = get_radial_network_reduction(network_model) + if isempty(radial_network_reduction) + ac_bus_numbers = PSY.get_number.(get_available_components(PSY.ACBus, sys)) + else + bus_reduction_map = PNM.get_bus_reduction_map(radial_network_reduction) + ac_bus_numbers = collect(keys(bus_reduction_map)) + end + + for t in time_steps, bus_no in ac_bus_numbers + pm_data["nw"]["$(t)"]["bus"]["$bus_no"]["inj_p"] = container.expressions[ExpressionKey(ActivePowerBalance, PSY.ACBus)][ - bus.number, + bus_no, t, ] - pm_data["nw"]["$(t)"]["bus"]["$(bus.number)"]["inj_q"] = + pm_data["nw"]["$(t)"]["bus"]["$bus_no"]["inj_q"] = container.expressions[ExpressionKey(ReactivePowerBalance, PSY.ACBus)][ - bus.number, + bus_no, t, ] end @@ -321,10 +337,19 @@ function powermodels_network!( pm_data, PM_map = pass_to_pm(sys, template, time_steps[end]) buses = get_available_components(PSY.ACBus, sys) - for t in time_steps, bus in buses - pm_data["nw"]["$(t)"]["bus"]["$(PSY.get_number(bus))"]["inj_p"] = + network_model = get_network_model(template) + radial_network_reduction = get_radial_network_reduction(network_model) + if isempty(radial_network_reduction) + ac_bus_numbers = PSY.get_number.(get_available_components(PSY.ACBus, sys)) + else + bus_reduction_map = PNM.get_bus_reduction_map(radial_network_reduction) + ac_bus_numbers = collect(keys(bus_reduction_map)) + end + + for t in time_steps, bus_no in ac_bus_numbers + pm_data["nw"]["$(t)"]["bus"]["$bus_no"]["inj_p"] = container.expressions[ExpressionKey(ActivePowerBalance, PSY.ACBus)][ - PSY.get_number(bus), + bus_no, t, ] # pm_data["nw"]["$(t)"]["bus"]["$(bus.number)"]["inj_q"] = 0.0 diff --git a/src/operation/operation_model_interface.jl b/src/operation/operation_model_interface.jl index 9a1e01c0c0..faa323ae6b 100644 --- a/src/operation/operation_model_interface.jl +++ b/src/operation/operation_model_interface.jl @@ -415,7 +415,7 @@ end function serialize_optimization_model(model::OperationModel, save_path::String) serialize_jump_optimization_model( - get_optimization_container(model), + get_jump_model(get_optimization_container(model)), save_path, ) return diff --git a/test/test_network_constructors.jl b/test/test_network_constructors.jl index 6e4860a9d3..519b19d28f 100644 --- a/test/test_network_constructors.jl +++ b/test/test_network_constructors.jl @@ -1,29 +1,30 @@ # Note to devs. Use GLPK or Cbc for models with linear constraints and linear cost functions # Use OSQP for models with quadratic cost function and linear constraints and ipopt otherwise +const NETWORKS_FOR_TESTING = [ + (PM.ACPPowerModel, fast_ipopt_optimizer), + (PM.ACRPowerModel, fast_ipopt_optimizer), + (PM.ACTPowerModel, fast_ipopt_optimizer), + #(PM.IVRPowerModel, fast_ipopt_optimizer), #instantiate_ivp_expr_model not implemented + (PM.DCPPowerModel, fast_ipopt_optimizer), + (PM.DCMPPowerModel, fast_ipopt_optimizer), + (PM.NFAPowerModel, fast_ipopt_optimizer), + (PM.DCPLLPowerModel, fast_ipopt_optimizer), + (PM.LPACCPowerModel, fast_ipopt_optimizer), + (PM.SOCWRPowerModel, fast_ipopt_optimizer), + (PM.SOCWRConicPowerModel, scs_solver), + (PM.QCRMPowerModel, fast_ipopt_optimizer), + (PM.QCLSPowerModel, fast_ipopt_optimizer), + #(PM.SOCBFPowerModel, fast_ipopt_optimizer), # not implemented + (PM.BFAPowerModel, fast_ipopt_optimizer), + #(PM.SOCBFConicPowerModel, fast_ipopt_optimizer), # not implemented + (PM.SDPWRMPowerModel, scs_solver), + (PM.SparseSDPWRMPowerModel, scs_solver), + (PTDFPowerModel, fast_ipopt_optimizer), +] + @testset "All PowerModels models construction" begin - networks = [ - (PM.ACPPowerModel, fast_ipopt_optimizer), - (PM.ACRPowerModel, fast_ipopt_optimizer), - (PM.ACTPowerModel, fast_ipopt_optimizer), - #(PM.IVRPowerModel, fast_ipopt_optimizer), #instantiate_ivp_expr_model not implemented - (PM.DCPPowerModel, fast_ipopt_optimizer), - (PM.DCMPPowerModel, fast_ipopt_optimizer), - (PM.NFAPowerModel, fast_ipopt_optimizer), - (PM.DCPLLPowerModel, fast_ipopt_optimizer), - (PM.LPACCPowerModel, fast_ipopt_optimizer), - (PM.SOCWRPowerModel, fast_ipopt_optimizer), - (PM.SOCWRConicPowerModel, scs_solver), - (PM.QCRMPowerModel, fast_ipopt_optimizer), - (PM.QCLSPowerModel, fast_ipopt_optimizer), - #(PM.SOCBFPowerModel, fast_ipopt_optimizer), # not implemented - (PM.BFAPowerModel, fast_ipopt_optimizer), - #(PM.SOCBFConicPowerModel, fast_ipopt_optimizer), # not implemented - (PM.SDPWRMPowerModel, scs_solver), - (PM.SparseSDPWRMPowerModel, scs_solver), - (PTDFPowerModel, fast_ipopt_optimizer), - ] c_sys5 = PSB.build_system(PSITestSystems, "c_sys5") - for (network, solver) in networks + for (network, solver) in NETWORKS_FOR_TESTING template = get_thermal_dispatch_template_network( NetworkModel(network; PTDF_matrix = PTDF(c_sys5)), ) @@ -696,3 +697,86 @@ end ) @test all(isapprox.(sum(zone_2_gen .+ zone_2_load; dims = 2), 0.0; atol = 1e-3)) end + +# These models are easier to test due to their lossless nature +@testset "StandardPTDF/DCPPowerModel Radial Branches Test" begin + new_sys = PSB.build_system(PSITestSystems, "c_sys5_radial") + + for net_model in [DCPPowerModel, StandardPTDFModel] + template_uc = template_unit_commitment(; + network = NetworkModel(net_model; + reduce_radial_branches = true, + use_slacks = false, + ), + ) + thermal_model = ThermalStandardUnitCommitment + set_device_model!(template_uc, ThermalStandard, thermal_model) + + ##### Solve Reduced Model #### + solver = GLPK_optimizer + uc_model_red = DecisionModel( + template_uc, + new_sys; + optimizer = solver, + name = "UC_RED", + store_variable_names = true, + ) + + @test build!(uc_model_red; output_dir = mktempdir(; cleanup = true)) == + PSI.BuildStatus.BUILT + solve!(uc_model_red) + + res_red = ProblemResults(uc_model_red) + + flow_lines = read_variable(res_red, "FlowActivePowerVariable__Line") + line_names = DataFrames.names(flow_lines)[2:end] + + ##### Solve Original Model #### + template_uc_orig = template_unit_commitment(; + network = NetworkModel(net_model; + reduce_radial_branches = false, + use_slacks = false, + ), + ) + set_device_model!(template_uc_orig, ThermalStandard, thermal_model) + + uc_model_orig = DecisionModel( + template_uc_orig, + new_sys; + optimizer = solver, + name = "UC_ORIG", + store_variable_names = true, + ) + + @test build!(uc_model_orig; output_dir = mktempdir(; cleanup = true)) == + PSI.BuildStatus.BUILT + solve!(uc_model_orig) + + res_orig = ProblemResults(uc_model_orig) + + flow_lines_orig = read_variable(res_orig, "FlowActivePowerVariable__Line") + + for line in line_names + @test isapprox(flow_lines[!, line], flow_lines_orig[!, line]) + end + end +end + +@testset "All PowerModels models construction with reduced radial branches" begin + new_sys = PSB.build_system(PSITestSystems, "c_sys5_radial") + for (network, solver) in NETWORKS_FOR_TESTING + if network ∈ PSI.INCOMPATIBLE_WITH_RADIAL_BRANCHES_POWERMODELS + continue + end + template = get_thermal_dispatch_template_network( + NetworkModel(network; + PTDF_matrix = PTDF(new_sys), + reduce_radial_branches = true, + use_slacks = true), + ) + ps_model = DecisionModel(template, new_sys; optimizer = solver) + @test build!(ps_model; output_dir = mktempdir(; cleanup = true)) == + PSI.BuildStatus.BUILT + @test ps_model.internal.container.pm !== nothing + end +end diff --git a/test/test_utils/mock_operation_models.jl b/test/test_utils/mock_operation_models.jl index a41128efb5..759684109c 100644 --- a/test/test_utils/mock_operation_models.jl +++ b/test/test_utils/mock_operation_models.jl @@ -117,6 +117,7 @@ function mock_construct_device!( PSI.get_network_formulation(template), PSI.get_network_model(template).subnetworks, PSI.get_system(problem), + Dict{Int64, Set{Int64}}(), ) if PSI.validate_available_devices(model, PSI.get_system(problem)) PSI.construct_device!(