diff --git a/src/PowerOperationsModels.jl b/src/PowerOperationsModels.jl index c9a0957..5a7967e 100644 --- a/src/PowerOperationsModels.jl +++ b/src/PowerOperationsModels.jl @@ -14,6 +14,7 @@ import ProgressMeter import PowerSystems import PowerSystems: get_component import Serialization +import SparseArrays import TimerOutputs import InteractiveUtils: methodswith diff --git a/src/ac_transmission_models/AC_branches.jl b/src/ac_transmission_models/AC_branches.jl index bf58fe5..ea24559 100644 --- a/src/ac_transmission_models/AC_branches.jl +++ b/src/ac_transmission_models/AC_branches.jl @@ -755,26 +755,48 @@ function add_constraints!( end function _make_flow_expressions!( - jump_model::JuMP.Model, name::String, time_steps::UnitRange{Int}, - ptdf_col::AbstractVector{Float64}, + ptdf_col::Vector{Float64}, nodal_balance_expressions::Matrix{JuMP.AffExpr}, ) @debug "Making Flow Expression on thread $(Threads.threadid()) for branch $name" + nz_idx = [i for i in eachindex(ptdf_col) if abs(ptdf_col[i]) > PTDF_ZERO_TOL] + hint = length(nz_idx) expressions = Vector{JuMP.AffExpr}(undef, length(time_steps)) for t in time_steps - expressions[t] = JuMP.@expression( - jump_model, - sum( - ptdf_col[i] * nodal_balance_expressions[i, t] for - i in 1:length(ptdf_col) + acc = IOM.get_hinted_aff_expr(hint) + @inbounds for i in nz_idx + JuMP.add_to_expression!(acc, ptdf_col[i], nodal_balance_expressions[i, t]) + end + expressions[t] = acc + end + return name, expressions +end + +function _make_flow_expressions!( + name::String, + time_steps::UnitRange{Int}, + ptdf_col::SparseArrays.SparseVector{Float64, Int}, + nodal_balance_expressions::Matrix{JuMP.AffExpr}, +) + @debug "Making Flow Expression on thread $(Threads.threadid()) for branch $name" + nz_idx = SparseArrays.nonzeroinds(ptdf_col) + nz_val = SparseArrays.nonzeros(ptdf_col) + hint = length(nz_idx) + expressions = Vector{JuMP.AffExpr}(undef, length(time_steps)) + for t in time_steps + acc = IOM.get_hinted_aff_expr(hint) + @inbounds for k in eachindex(nz_idx) + JuMP.add_to_expression!( + acc, + nz_val[k], + nodal_balance_expressions[nz_idx[k], t], ) - ) + end + expressions[t] = acc end return name, expressions - # change when using the not concurrent version - # return expressions end function _add_expression_to_container!( @@ -789,7 +811,6 @@ function _add_expression_to_container!( name = PSY.get_name(reduction_entry) if name in branches branch_flow_expr[name, :] .= _make_flow_expressions!( - jump_model, name, time_steps, ptdf_col, @@ -812,7 +833,6 @@ function _add_expression_to_container!( for name in names if name in branches branch_flow_expr[name, :] .= _make_flow_expressions!( - jump_model, name, time_steps, ptdf_col, @@ -837,7 +857,6 @@ function _add_expression_to_container!( for name in names if name in branches branch_flow_expr[name, :] .= _make_flow_expressions!( - jump_model, name, time_steps, ptdf_col, @@ -873,13 +892,10 @@ function add_expressions!( time_steps, ) - jump_model = get_jump_model(container) - - tasks = map(collect(name_to_arc_map)) do pair + tasks = map(name_to_arc_map) do pair (name, (arc, _)) = pair ptdf_col = ptdf[arc, :] Threads.@spawn _make_flow_expressions!( - jump_model, name, time_steps, ptdf_col, diff --git a/src/ac_transmission_models/branch_constructor.jl b/src/ac_transmission_models/branch_constructor.jl index 6d24f9a..568b072 100644 --- a/src/ac_transmission_models/branch_constructor.jl +++ b/src/ac_transmission_models/branch_constructor.jl @@ -147,13 +147,15 @@ function construct_device!( ) end - add_branch_parameters!( - container, - DynamicBranchRatingTimeSeriesParameter, - devices, - device_model, - network_model, - ) + if haskey(get_time_series_names(device_model), DynamicBranchRatingTimeSeriesParameter) + add_branch_parameters!( + container, + DynamicBranchRatingTimeSeriesParameter, + devices, + device_model, + network_model, + ) + end if haskey( get_time_series_names(device_model), @@ -260,7 +262,7 @@ function construct_device!( add_constraints!( container, - PostContingencyEmergencyRateLimitConstrain, + PostContingencyEmergencyRateLimitConstraint, branches, branches_outages, device_model, diff --git a/src/common_models/add_parameters.jl b/src/common_models/add_parameters.jl index 5ccb0db..e3387bc 100644 --- a/src/common_models/add_parameters.jl +++ b/src/common_models/add_parameters.jl @@ -27,7 +27,7 @@ function add_parameters!( if get_rebuild_model(get_settings(container)) && has_container_key(container, T, D) return end - _add_parameters!(container, T(), devices, model) + _add_parameters!(container, T, devices, model) return end @@ -41,7 +41,7 @@ function add_parameters!( has_container_key(container, T, U, PSY.get_name(service)) return end - _add_parameters!(container, T(), service, model) + _add_parameters!(container, T, service, model) return end @@ -59,7 +59,7 @@ function add_branch_parameters!( if get_rebuild_model(get_settings(container)) && has_container_key(container, T, D) return end - _add_time_series_parameters!(container, T(), network_model, devices, model) + _add_time_series_parameters!(container, T, network_model, devices, model) return end @@ -69,7 +69,7 @@ end function _add_parameters!( container::OptimizationContainer, - param::T, + param::Type{T}, devices::U, model::DeviceModel{D, W}, ) where { @@ -87,7 +87,7 @@ end function _check_dynamic_branch_rating_ts( ts::AbstractArray, - ::T, + ::Type{T}, device::PSY.Device, model::DeviceModel{D, W}, ) where {D <: PSY.Component, T <: TimeSeriesParameter, W <: AbstractDeviceFormulation} @@ -121,7 +121,7 @@ _size_wrapper(::Tuple) = () function _add_time_series_parameters!( container::OptimizationContainer, - param::T, + param::Type{T}, devices, model::DeviceModel{D, W}, ) where {D <: PSY.Component, T <: TimeSeriesParameter, W <: AbstractDeviceFormulation} @@ -131,25 +131,48 @@ function _add_time_series_parameters!( end time_steps = get_time_steps(container) - ts_name = _get_time_series_name(T(), first(devices), model) + ts_name = _get_time_series_name(T, first(devices), model) device_names = String[] devices_with_time_series = D[] initial_values = Dict{String, AbstractArray}() + # device name -> ts_uuid cache so the second loop below doesn't re-query IS. + device_ts_uuids = Dict{String, String}() + model_interval = get_interval(get_settings(container)) + is_ts_interval = _to_is_interval(model_interval) + model_resolution = get_resolution(get_settings(container)) + is_ts_resolution = _to_is_resolution(model_resolution) @debug "adding" T D ts_name ts_type _group = IOM.LOG_GROUP_OPTIMIZATION_CONTAINER for device::D in devices if !PSY.has_time_series(device, ts_type, ts_name) - @info "Time series $(ts_type):$(ts_name) for $D, $(PSY.get_name(device)) not found skipping parameter addition." + @debug "Time series $(ts_type):$(ts_name) for $D, $(PSY.get_name(device)) not found. Skipping parameter addition for this device." continue end - push!(device_names, PSY.get_name(device)) + device_name = PSY.get_name(device) + push!(device_names, device_name) push!(devices_with_time_series, device) - ts_uuid = string(IS.get_time_series_uuid(ts_type, device, ts_name)) + ts_uuid = string( + IS.get_time_series_uuid( + ts_type, + device, + ts_name; + resolution = is_ts_resolution, + interval = is_ts_interval, + ), + ) + device_ts_uuids[device_name] = ts_uuid if !(ts_uuid in keys(initial_values)) initial_values[ts_uuid] = - IOM.get_time_series_initial_values!(container, ts_type, device, ts_name) + IOM.get_time_series_initial_values!( + container, + ts_type, + device, + ts_name; + interval = model_interval, + resolution = model_resolution, + ) _check_dynamic_branch_rating_ts(initial_values[ts_uuid], param, device, model) end end @@ -195,7 +218,7 @@ function _add_time_series_parameters!( IOM.add_component_name!( IOM.get_attributes(param_container), device_name, - string(IS.get_time_series_uuid(ts_type, device, ts_name)), + device_ts_uuids[device_name], ) end return @@ -207,7 +230,7 @@ end function _add_time_series_parameters!( container::OptimizationContainer, - param::T, + ::Type{T}, network_model::NetworkModel{<:AbstractPTDFModel}, devices, model::DeviceModel{D, W}, @@ -223,13 +246,24 @@ function _add_time_series_parameters!( all_branch_maps_by_type = PNM.get_all_branch_maps_by_type(net_reduction_data) # TODO: Temporary workaround to get the name where we assume all the names are the same across devices. - ts_name = _get_time_series_name(T(), first(devices), model) + ts_name = _get_time_series_name(T, first(devices), model) + model_interval = get_interval(get_settings(container)) + ts_interval = model_interval device_name_axis, ts_uuid_axis = - get_branch_argument_parameter_axes(net_reduction_data, devices, ts_type, ts_name) + get_branch_argument_parameter_axes( + net_reduction_data, + devices, + ts_type, + ts_name; + interval = ts_interval, + ) if isempty(device_name_axis) @info "No devices with time series $ts_name found for $D devices. Skipping parameter addition." return end + # name -> ts_uuid cache built from the axis pair so the per-branch loop below + # doesn't re-query IS.get_time_series_uuid for each branch. + branch_ts_uuids = Dict{String, String}(zip(device_name_axis, ts_uuid_axis)) additional_axes = () param_container = add_param_container!( container, @@ -251,8 +285,7 @@ function _add_time_series_parameters!( if device_with_time_series === nothing continue end - ts_uuid = - string(IS.get_time_series_uuid(ts_type, device_with_time_series, ts_name)) + ts_uuid = branch_ts_uuids[name] has_entry, tracker_container = search_for_reduced_branch_parameter!( reduced_branch_tracker, @@ -267,7 +300,8 @@ function _add_time_series_parameters!( container, ts_type, device_with_time_series, - ts_name, + ts_name; + interval = ts_interval, ) ts_vals = _unwrap_for_param.(Ref(param_instance), raw_ts_vals, Ref(additional_axes)) @@ -298,30 +332,34 @@ function _add_time_series_parameters!( IOM.add_component_name!( IOM.get_attributes(param_container), name, - string(IS.get_time_series_uuid(ts_type, device_with_time_series, ts_name)), + ts_uuid, ) end return end -_get_time_series_name(::T, ::PSY.Component, model::DeviceModel) where {T <: ParameterType} = +_get_time_series_name( + ::Type{T}, + ::PSY.Component, + model::DeviceModel, +) where {T <: ParameterType} = get_time_series_names(model)[T] -_get_time_series_name(::StartupCostParameter, device::PSY.Component, ::DeviceModel) = +_get_time_series_name(::Type{StartupCostParameter}, device::PSY.Component, ::DeviceModel) = IS.get_name(PSY.get_start_up(PSY.get_operation_cost(device))) -_get_time_series_name(::ShutdownCostParameter, device::PSY.Component, ::DeviceModel) = +_get_time_series_name(::Type{ShutdownCostParameter}, device::PSY.Component, ::DeviceModel) = IS.get_name(PSY.get_shut_down(PSY.get_operation_cost(device))) _get_time_series_name( - ::IncrementalCostAtMinParameter, + ::Type{IncrementalCostAtMinParameter}, device::PSY.Device, ::DeviceModel, ) = IS.get_name(PSY.get_incremental_initial_input(PSY.get_operation_cost(device))) _get_time_series_name( - ::DecrementalCostAtMinParameter, + ::Type{DecrementalCostAtMinParameter}, device::PSY.Device, ::DeviceModel, ) = @@ -331,33 +369,37 @@ _get_time_series_name( # _get_expected_time_series_eltype — for ObjectiveFunctionParameter ################################################################################# -_get_expected_time_series_eltype(::T) where {T <: ParameterType} = Float64 -_get_expected_time_series_eltype(::StartupCostParameter) = NTuple{3, Float64} +_get_expected_time_series_eltype(::Type{T}) where {T <: ParameterType} = Float64 +_get_expected_time_series_eltype(::Type{StartupCostParameter}) = NTuple{3, Float64} ################################################################################# # _param_to_vars — lookup: ObjectiveFunctionParameter → variable types ################################################################################# -_param_to_vars(::FuelCostParameter, ::AbstractDeviceFormulation) = (ActivePowerVariable,) -_param_to_vars(::StartupCostParameter, ::AbstractThermalFormulation) = (StartVariable,) -_param_to_vars(::StartupCostParameter, ::ThermalMultiStartUnitCommitment) = +_param_to_vars(::Type{FuelCostParameter}, ::Type{<:AbstractDeviceFormulation}) = + (ActivePowerVariable,) +_param_to_vars(::Type{StartupCostParameter}, ::Type{<:AbstractThermalFormulation}) = + (StartVariable,) +_param_to_vars(::Type{StartupCostParameter}, ::Type{ThermalMultiStartUnitCommitment}) = MULTI_START_VARIABLES -_param_to_vars(::ShutdownCostParameter, ::AbstractThermalFormulation) = (StopVariable,) -_param_to_vars(::AbstractCostAtMinParameter, ::AbstractDeviceFormulation) = (OnVariable,) +_param_to_vars(::Type{ShutdownCostParameter}, ::Type{<:AbstractThermalFormulation}) = + (StopVariable,) +_param_to_vars(::Type{<:AbstractCostAtMinParameter}, ::Type{<:AbstractDeviceFormulation}) = + (OnVariable,) _param_to_vars( ::Union{ - IncrementalPiecewiseLinearSlopeParameter, - IncrementalPiecewiseLinearBreakpointParameter, + Type{IncrementalPiecewiseLinearSlopeParameter}, + Type{IncrementalPiecewiseLinearBreakpointParameter}, }, - ::AbstractDeviceFormulation, + ::Type{<:AbstractDeviceFormulation}, ) = (PiecewiseLinearBlockIncrementalOffer,) _param_to_vars( ::Union{ - DecrementalPiecewiseLinearSlopeParameter, - DecrementalPiecewiseLinearBreakpointParameter, + Type{DecrementalPiecewiseLinearSlopeParameter}, + Type{DecrementalPiecewiseLinearBreakpointParameter}, }, - ::AbstractDeviceFormulation, + ::Type{<:AbstractDeviceFormulation}, ) = (PiecewiseLinearBlockDecrementalOffer,) @@ -367,7 +409,7 @@ _param_to_vars( calc_additional_axes( ::OptimizationContainer, - ::T, + ::Type{T}, ::U, ::DeviceModel{D, W}, ) where { @@ -378,7 +420,7 @@ calc_additional_axes( calc_additional_axes( ::OptimizationContainer, - ::T, + ::Type{T}, ::U, ::ServiceModel{D, W}, ) where { @@ -405,7 +447,7 @@ _unwrap_for_param(::ParameterType, ts_elem, expected_axs) = ts_elem function _add_parameters!( container::OptimizationContainer, - param::T, + param::Type{T}, devices::U, model::DeviceModel{D, W}, ) where { @@ -425,7 +467,7 @@ function _add_parameters!( device_names = String[] active_devices = D[] for device in devices - ts_name = _get_time_series_name(T(), device, model) + ts_name = _get_time_series_name(T, device, model) if PSY.has_time_series(device, ts_type, ts_name) push!(ts_names, ts_name) push!(device_names, PSY.get_name(device)) @@ -444,19 +486,27 @@ function _add_parameters!( container, T, D, - _param_to_vars(T(), W()), + _param_to_vars(T, W), SOSStatusVariable.NO_VARIABLE, false, - _get_expected_time_series_eltype(T()), + _get_expected_time_series_eltype(T), device_names, additional_axes..., time_steps, ) param_instance = T() + model_interval = get_interval(get_settings(container)) + ts_interval = model_interval for (ts_name, device_name, device) in zip(ts_names, device_names, active_devices) raw_ts_vals = - IOM.get_time_series_initial_values!(container, ts_type, device, ts_name) + IOM.get_time_series_initial_values!( + container, + ts_type, + device, + ts_name; + interval = ts_interval, + ) ts_vals = _unwrap_for_param.(Ref(param_instance), raw_ts_vals, Ref(additional_axes)) @assert all(_size_wrapper.(ts_vals) .== Ref(length.(additional_axes))) for step in time_steps @@ -484,7 +534,7 @@ end function _add_parameters!( container::OptimizationContainer, - ::T, + ::Type{T}, service::U, model::ServiceModel{U, V}, ) where {T <: TimeSeriesParameter, U <: PSY.Service, V <: AbstractServiceFormulation} @@ -495,9 +545,18 @@ function _add_parameters!( ts_name = get_time_series_names(model)[T] time_steps = get_time_steps(container) name = PSY.get_name(service) - ts_uuid = string(IS.get_time_series_uuid(ts_type, service, ts_name)) + model_interval = get_interval(get_settings(container)) + ts_interval = model_interval + ts_uuid = string( + IS.get_time_series_uuid( + ts_type, + service, + ts_name; + interval = _to_is_interval(ts_interval), + ), + ) @debug "adding" T U _group = IOM.LOG_GROUP_OPTIMIZATION_CONTAINER - additional_axes = calc_additional_axes(container, T(), [service], model) + additional_axes = calc_additional_axes(container, T, [service], model) parameter_container = add_param_container!(container, T, U, ts_type, @@ -511,7 +570,7 @@ function _add_parameters!( IOM.set_subsystem!(IOM.get_attributes(parameter_container), IOM.get_subsystem(model)) jump_model = get_jump_model(container) - ts_vector = IOM.get_time_series(container, service, T, name) + ts_vector = IOM.get_time_series(container, service, T, name; interval = ts_interval) multiplier = get_multiplier_value(T, service, V) for t in time_steps IOM.set_multiplier!(parameter_container, multiplier, name, t) @@ -527,7 +586,7 @@ end function _add_parameters!( container::OptimizationContainer, - ::T, + ::Type{T}, key::VariableKey{U, D}, model::DeviceModel{D, W}, devices::V, @@ -574,7 +633,7 @@ end function _add_parameters!( container::OptimizationContainer, - ::T, + ::Type{T}, key::VariableKey{U, D}, model::DeviceModel{D, W}, devices::V, @@ -624,7 +683,7 @@ end function _add_parameters!( container::OptimizationContainer, - ::T, + ::Type{T}, key::VariableKey{U, D}, model::DeviceModel{D, W}, devices::V, @@ -672,7 +731,7 @@ end function _add_parameters!( container::OptimizationContainer, - ::T, + ::Type{T}, key::AuxVarKey{U, D}, model::DeviceModel{D, W}, devices::V, @@ -720,7 +779,7 @@ end function _add_parameters!( container::OptimizationContainer, - ::T, + ::Type{T}, devices::V, model::DeviceModel{D, W}, ) where { @@ -770,7 +829,7 @@ end function _add_parameters!( container::OptimizationContainer, - ::T, + ::Type{T}, key::VariableKey{U, S}, model::ServiceModel{S, W}, devices::V, diff --git a/src/core/constraints.jl b/src/core/constraints.jl index c9007d7..e74ade2 100644 --- a/src/core/constraints.jl +++ b/src/core/constraints.jl @@ -128,9 +128,6 @@ The specified constraint is formulated as: """ struct FeedforwardLowerBoundConstraint <: ConstraintType end struct FeedforwardEnergyTargetConstraint <: ConstraintType end -struct FlowActivePowerConstraint <: ConstraintType end #not being used -struct FlowActivePowerFromToConstraint <: ConstraintType end #not being used -struct FlowActivePowerToFromConstraint <: ConstraintType end #not being used """ Struct to create the constraint that set the flow limits through a PhaseShiftingTransformer. @@ -146,9 +143,6 @@ struct FlowLimitConstraint <: ConstraintType end struct FlowLimitFromToConstraint <: ConstraintType end struct FlowLimitToFromConstraint <: ConstraintType end -struct FlowReactivePowerConstraint <: ConstraintType end #not being used -struct FlowReactivePowerFromToConstraint <: ConstraintType end #not being used -struct FlowReactivePowerToFromConstraint <: ConstraintType end #not being used """ Struct to create the constraints that set the power balance across a lossy HVDC two-terminal line. @@ -251,7 +245,7 @@ The specified constraint is formulated as: ``` """ struct FlowRateConstraint <: ConstraintType end -struct PostContingencyEmergencyRateLimitConstrain <: PostContingencyConstraintType end +struct PostContingencyEmergencyRateLimitConstraint <: PostContingencyConstraintType end """ Struct to create the constraint for branch flow rate limits from the 'from' bus to the 'to' bus. diff --git a/src/core/definitions.jl b/src/core/definitions.jl index aa6823c..ef53c43 100644 --- a/src/core/definitions.jl +++ b/src/core/definitions.jl @@ -54,6 +54,7 @@ const JumpSupportedLiterals = # Settings constants const UNSET_HORIZON = Dates.Millisecond(0) const UNSET_RESOLUTION = Dates.Millisecond(0) +const UNSET_INTERVAL = Dates.Millisecond(0) const UNSET_INI_TIME = Dates.DateTime(0) # Tolerance of comparisons @@ -63,6 +64,7 @@ const BALANCE_SLACK_COST = 1e6 const CONSTRAINT_VIOLATION_SLACK_COST = 2e5 const SERVICES_SLACK_COST = 1e5 const COST_EPSILON = 1e-3 +const PTDF_ZERO_TOL = 1e-9 const MISSING_INITIAL_CONDITIONS_TIME_COUNT = 999.0 const MAX_START_STAGES = 3 diff --git a/src/network_models/instantiate_network_model.jl b/src/network_models/instantiate_network_model.jl index f30cb18..ca0342e 100644 --- a/src/network_models/instantiate_network_model.jl +++ b/src/network_models/instantiate_network_model.jl @@ -263,6 +263,7 @@ function IOM.instantiate_network_model!( @info "Applying both radial and degree two reductions" ptdf = PNM.VirtualPTDF( sys; + tol = PTDF_ZERO_TOL, network_reductions = PNM.NetworkReduction[ PNM.RadialReduction(; irreducible_buses = irreducible_buses), PNM.DegreeTwoReduction(; @@ -277,6 +278,7 @@ function IOM.instantiate_network_model!( end ptdf = PNM.VirtualPTDF( sys; + tol = PTDF_ZERO_TOL, network_reductions = PNM.NetworkReduction[PNM.RadialReduction(; irreducible_buses = irreducible_buses, )], @@ -285,12 +287,13 @@ function IOM.instantiate_network_model!( @info "Applying degree two reduction" ptdf = PNM.VirtualPTDF( sys; + tol = PTDF_ZERO_TOL, network_reductions = PNM.NetworkReduction[PNM.DegreeTwoReduction(; irreducible_buses = irreducible_buses, )], ) else - ptdf = PNM.VirtualPTDF(sys) + ptdf = PNM.VirtualPTDF(sys; tol = PTDF_ZERO_TOL) end model.PTDF_matrix = ptdf model.network_reduction = deepcopy(ptdf.network_reduction_data) @@ -371,6 +374,7 @@ function IOM.instantiate_network_model!( @info "Applying both radial and degree two reductions" ptdf = PNM.VirtualPTDF( sys; + tol = PTDF_ZERO_TOL, network_reductions = PNM.NetworkReduction[ PNM.RadialReduction(; irreducible_buses = irreducible_buses), PNM.DegreeTwoReduction(; irreducible_buses = irreducible_buses), @@ -383,6 +387,7 @@ function IOM.instantiate_network_model!( end ptdf = PNM.VirtualPTDF( sys; + tol = PTDF_ZERO_TOL, network_reductions = PNM.NetworkReduction[PNM.RadialReduction(; irreducible_buses = irreducible_buses, )], @@ -391,12 +396,13 @@ function IOM.instantiate_network_model!( @info "Applying degree two reduction" ptdf = PNM.VirtualPTDF( sys; + tol = PTDF_ZERO_TOL, network_reductions = PNM.NetworkReduction[PNM.DegreeTwoReduction(; irreducible_buses = irreducible_buses, )], ) else - ptdf = PNM.VirtualPTDF(sys) + ptdf = PNM.VirtualPTDF(sys; tol = PTDF_ZERO_TOL) end model.PTDF_matrix = ptdf model.network_reduction = deepcopy(ptdf.network_reduction_data) diff --git a/src/network_models/security_constrained_models.jl b/src/network_models/security_constrained_models.jl index 1bdf4d6..b7f7d4c 100644 --- a/src/network_models/security_constrained_models.jl +++ b/src/network_models/security_constrained_models.jl @@ -3,7 +3,7 @@ Min and max limits for post-contingency branch flows for Abstract Branch Formula """ function get_min_max_limits( branch::PSY.ACTransmission, - ::Type{<:PostContingencyEmergencyRateLimitConstrain}, + ::Type{<:PostContingencyEmergencyRateLimitConstraint}, ::Type{<:AbstractBranchFormulation}, ::NetworkModel{<:AbstractPTDFModel}, ) @@ -20,7 +20,7 @@ Add branch post-contingency rate limit constraints for ACBranch considering LODF """ function add_constraints!( container::OptimizationContainer, - cons_type::Type{PostContingencyEmergencyRateLimitConstrain}, + cons_type::Type{PostContingencyEmergencyRateLimitConstraint}, branches::IS.FlattenIteratorWrapper{PSY.ACTransmission}, branches_outages::Vector{T}, device_model::DeviceModel{T, U}, @@ -86,7 +86,7 @@ function add_constraints!( limits = get_min_max_limits( branch, - PostContingencyEmergencyRateLimitConstrain, + PostContingencyEmergencyRateLimitConstraint, U, network_model, ) diff --git a/src/services_models/reserves.jl b/src/services_models/reserves.jl index 671005b..c8eccb0 100644 --- a/src/services_models/reserves.jl +++ b/src/services_models/reserves.jl @@ -131,6 +131,17 @@ function add_reserve_variables!( return end +function _sum_reserve_variables( + vars::AbstractArray{<:JuMP.AbstractVariableRef}, + extra::Int, +) + acc = IOM.get_hinted_aff_expr(length(vars) + extra) + for v in vars + JuMP.add_to_expression!(acc, v) + end + return acc +end + ################################## Reserve Requirement Constraint ########################## function add_constraints!( container::OptimizationContainer, @@ -158,35 +169,35 @@ function add_constraints!( get_variable(container, ActivePowerReserveVariable, SR, service_name) use_slacks = get_use_slacks(model) - ts_vector = IOM.get_time_series(container, service, "requirement") + ts_vector = IOM.get_time_series( + container, + service, + "requirement"; + interval = get_interval(get_settings(container)), + ) use_slacks && (slack_vars = reserve_slacks!(container, service)) requirement = PSY.get_requirement(service) jump_model = get_jump_model(container) + extra = use_slacks ? 1 : 0 if built_for_recurrent_solves(container) param_container = get_parameter(container, RequirementTimeSeriesParameter, SR, service_name) param = get_parameter_column_refs(param_container, service_name) for t in time_steps - if use_slacks - resource_expression = JuMP.@expression( - jump_model, sum(@view reserve_variable[:, t]) + slack_vars[t]) - else - resource_expression = JuMP.@expression( - jump_model, sum(@view reserve_variable[:, t])) - end + resource_expression = + _sum_reserve_variables(@view(reserve_variable[:, t]), extra) + use_slacks && + JuMP.add_to_expression!(resource_expression, slack_vars[t]) constraint[service_name, t] = JuMP.@constraint(jump_model, resource_expression >= param[t] * requirement) end else for t in time_steps - if use_slacks - resource_expression = JuMP.@expression( - jump_model, sum(@view reserve_variable[:, t]) + slack_vars[t]) - else - resource_expression = JuMP.@expression( - jump_model, sum(@view reserve_variable[:, t])) - end + resource_expression = + _sum_reserve_variables(@view(reserve_variable[:, t]), extra) + use_slacks && + JuMP.add_to_expression!(resource_expression, slack_vars[t]) constraint[service_name, t] = JuMP.@constraint( jump_model, resource_expression >= ts_vector[t] * requirement @@ -224,7 +235,12 @@ function add_constraints!( var_r = get_variable(container, ActivePowerReserveVariable, SR, service_name) jump_model = get_jump_model(container) requirement = PSY.get_requirement(service) - ts_vector = IOM.get_time_series(container, service, "requirement") + ts_vector = IOM.get_time_series( + container, + service, + "requirement"; + interval = get_interval(get_settings(container)), + ) param_container = get_parameter(container, RequirementTimeSeriesParameter, SR, service_name) param = get_parameter_column_refs(param_container, service_name) @@ -274,14 +290,11 @@ function add_constraints!( requirement = PSY.get_requirement(service) jump_model = get_jump_model(container) + extra = use_slacks ? 1 : 0 for t in time_steps - resource_expression = JuMP.GenericAffExpr{Float64, JuMP.VariableRef}() - JuMP.add_to_expression!(resource_expression, - JuMP.@expression(jump_model, sum(@view reserve_variable[:, t]))) - # consider a for loop - if use_slacks - JuMP.add_to_expression!(resource_expression, slack_vars[t]) - end + resource_expression = + _sum_reserve_variables(@view(reserve_variable[:, t]), extra) + use_slacks && JuMP.add_to_expression!(resource_expression, slack_vars[t]) constraint[service_name, t] = JuMP.@constraint(jump_model, resource_expression >= requirement) end diff --git a/src/static_injector_models/thermal_generation.jl b/src/static_injector_models/thermal_generation.jl index 85ca712..dfd3b2c 100644 --- a/src/static_injector_models/thermal_generation.jl +++ b/src/static_injector_models/thermal_generation.jl @@ -94,6 +94,7 @@ initial_condition_variable(::InitialTimeDurationOff, d::PSY.ThermalGen, ::Abstra function proportional_cost(container::OptimizationContainer, cost::PSY.ThermalGenerationCost, S::Type{OnVariable}, T::PSY.ThermalGen, U::Type{<:AbstractThermalFormulation}, t::Int) return onvar_cost(container, cost, S, T, U, t) + PSY.get_constant_term(PSY.get_vom_cost(PSY.get_variable(cost))) + PSY.get_fixed(cost) end +is_time_variant_term(::OptimizationContainer, ::PSY.ThermalGenerationCost, ::OnVariable, ::PSY.ThermalGen, ::AbstractThermalFormulation, t::Int) = false is_time_variant_term(::OptimizationContainer, ::PSY.ThermalGenerationCost, ::Type{OnVariable}, ::Type{<:PSY.ThermalGen}, ::Type{<:AbstractThermalFormulation}, t::Int) = false diff --git a/src/twoterminal_hvdc_models/branch_constructor.jl b/src/twoterminal_hvdc_models/branch_constructor.jl index c1bd087..60aab5f 100644 --- a/src/twoterminal_hvdc_models/branch_constructor.jl +++ b/src/twoterminal_hvdc_models/branch_constructor.jl @@ -244,7 +244,7 @@ function construct_device!( add_constraints!( container, - PostContingencyEmergencyRateLimitConstrain, + PostContingencyEmergencyRateLimitConstraint, branches, branches_outages, device_model, diff --git a/src/utils/psy_utils.jl b/src/utils/psy_utils.jl index 43ad022..b42e7af 100644 --- a/src/utils/psy_utils.jl +++ b/src/utils/psy_utils.jl @@ -1,3 +1,9 @@ +_to_is_interval(interval::Dates.Millisecond) = + interval == UNSET_INTERVAL ? nothing : interval + +_to_is_resolution(resolution::Dates.Millisecond) = + resolution == UNSET_RESOLUTION ? nothing : resolution + function get_available_reservoirs(sys::PSY.System) return PSY.get_components( x -> (PSY.get_available(x)), diff --git a/test/performance/performance_test.jl b/test/performance/performance_test.jl index 7b15857..e061741 100644 --- a/test/performance/performance_test.jl +++ b/test/performance/performance_test.jl @@ -6,6 +6,7 @@ using InfrastructureOptimizationModels const IOM = InfrastructureOptimizationModels using PowerSystems const PSY = PowerSystems +import InfrastructureSystems as IS using Logging using PowerSystemCaseBuilder using PowerNetworkMatrices @@ -53,12 +54,31 @@ function set_device_models!(template::OperationsProblemTemplate, uc::Bool = true end try + # Build both systems, then merge the 5-minute SingleTimeSeries from the + # realization system onto the DA system so a single System carries raw + # 1-hour and 5-minute data that can be transformed separately per model. sys_rts_da = build_system(PSISystems, "modified_RTS_GMLC_DA_sys") sys_rts_rt = build_system(PSISystems, "modified_RTS_GMLC_RT_sys") - for sys in [sys_rts_da, sys_rts_rt] - g = get_component(ThermalStandard, sys, "121_NUCLEAR_1") - set_must_run!(g, true) + # Drop the transform that PSB pre-baked so we can attach new per-resolution + # transforms and leave both static series intact. + PSY.transform_single_time_series!( + sys_rts_da, + Hour(48), + Hour(24); + resolution = Hour(1), + delete_existing = true, + ) + PSY.transform_single_time_series!( + sys_rts_da, + Hour(1), + Minute(15); + resolution = Minute(5), + delete_existing = false, + ) + + for g in get_components(ThermalStandard, sys_rts_da) + get_name(g) == "121_NUCLEAR_1" && set_must_run!(g, true) end for i in 1:2 @@ -66,7 +86,6 @@ try NetworkModel( PTDFPowerModel; use_slacks = true, - PTDF_matrix = PTDF(sys_rts_da), duals = [CopperPlateBalanceConstraint], ), ) @@ -76,7 +95,6 @@ try NetworkModel( PTDFPowerModel; use_slacks = true, - PTDF_matrix = PTDF(sys_rts_rt), duals = [CopperPlateBalanceConstraint], ), ) @@ -93,17 +111,23 @@ try optimizer_solve_log_print = false, direct_mode_optimizer = true, check_numerical_bounds = false, + horizon = Hour(48), + interval = Hour(24), + resolution = Hour(1), ) ed = DecisionModel( template_ed, - sys_rts_rt; + sys_rts_da; name = "ED", optimizer = optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.01, "log_to_console" => false), initialize_model = true, check_numerical_bounds = false, + horizon = Hour(48), + interval = Hour(24), + resolution = Hour(1), ) # Build diff --git a/test/test_model_decision.jl b/test/test_model_decision.jl index 87edcc7..0e4ccd2 100644 --- a/test/test_model_decision.jl +++ b/test/test_model_decision.jl @@ -214,6 +214,46 @@ end end end +@testset "Test Locational Marginal Prices between DC lossless with PowerModels vs PTDFPowerModel" begin + networks = [DCPPowerModel, PTDFPowerModel] + sys = PSB.build_system(PSITestSystems, "c_sys5") + ptdf = VirtualPTDF(sys) + # These are the duals of interest for the test + dual_constraint = [[NodalBalanceActiveConstraint], [CopperPlateBalanceConstraint]] + LMPs = [] + for (ix, network) in enumerate(networks) + template = get_template_dispatch_with_network( + NetworkModel(network; PTDF_matrix = ptdf, duals = dual_constraint[ix]), + ) + if network == PTDFPowerModel + set_device_model!( + template, + DeviceModel(PSY.Line, StaticBranch; duals = [FlowRateConstraint]), + ) + end + model = DecisionModel(template, sys; optimizer = HiGHS_optimizer) + @test build!(model; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED + res = OptimizationProblemOutputs(model) + + # These tests require results to be working + if network == PTDFPowerModel + push!(LMPs, abs.(psi_ptdf_lmps(res, ptdf))) + else + duals = read_dual( + res, + NodalBalanceActiveConstraint, + ACBus; + table_format = TableFormat.WIDE, + ) + duals = abs.(duals[:, propertynames(duals) .!== :DateTime]) + push!(LMPs, duals[!, sort(propertynames(duals))]) + end + end + @test isapprox(LMPs[1], LMPs[2], atol = 100.0) +end + @testset "Test OptimizationProblemOutputs interfaces" begin sys = PSB.build_system(PSITestSystems, "c_sys5_re") template = get_template_dispatch_with_network( diff --git a/test/test_multi_interval.jl b/test/test_multi_interval.jl new file mode 100644 index 0000000..8d15796 --- /dev/null +++ b/test/test_multi_interval.jl @@ -0,0 +1,136 @@ +@testset "Multi-interval validation errors" begin + c_sys5 = PSB.build_system(PSITestSystems, "c_sys5"; add_single_time_series = true) + transform_single_time_series!(c_sys5, Hour(24), Hour(12); delete_existing = false) + transform_single_time_series!(c_sys5, Hour(24), Hour(24); delete_existing = false) + + template = get_thermal_dispatch_template_network() + + # Without interval kwarg on a multi-interval system, should error + @test_throws IS.ConflictingInputsError DecisionModel( + template, + c_sys5; + optimizer = HiGHS_optimizer, + horizon = Hour(24), + ) + + # With a non-existent interval on a system that has no SingleTimeSeries to auto-transform from, should error. + c_sys5_forecasts = PSB.build_system(PSITestSystems, "c_sys5"; add_forecasts = true) + @test_throws IS.ConflictingInputsError DecisionModel( + template, + c_sys5_forecasts; + optimizer = HiGHS_optimizer, + horizon = Hour(24), + interval = Hour(6), + ) +end + +@testset "DecisionModel with explicit interval" begin + c_sys5 = PSB.build_system(PSITestSystems, "c_sys5"; add_single_time_series = true) + transform_single_time_series!(c_sys5, Hour(24), Hour(12); delete_existing = false) + transform_single_time_series!(c_sys5, Hour(24), Hour(24); delete_existing = false) + + template = get_thermal_dispatch_template_network() + + # Build with interval = 24h + model_24h = DecisionModel( + template, + c_sys5; + optimizer = HiGHS_optimizer, + horizon = Hour(24), + interval = Hour(24), + ) + @test build!(model_24h; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + @test solve!(model_24h) == IOM.RunStatus.SUCCESSFULLY_FINALIZED + + # Build with interval = 12h on the same system + model_12h = DecisionModel( + template, + c_sys5; + optimizer = HiGHS_optimizer, + horizon = Hour(24), + interval = Hour(12), + ) + @test build!(model_12h; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + @test solve!(model_12h) == IOM.RunStatus.SUCCESSFULLY_FINALIZED + + # Verify models have the correct interval in store params + @test IOM.get_interval(model_24h) == Dates.Millisecond(Hour(24)) + @test IOM.get_interval(model_12h) == Dates.Millisecond(Hour(12)) +end + +@testset "Auto-transform SingleTimeSeries with interval" begin + c_sys5 = PSB.build_system(PSITestSystems, "c_sys5"; add_single_time_series = true) + template = get_thermal_dispatch_template_network() + + model = DecisionModel( + template, + c_sys5; + optimizer = HiGHS_optimizer, + horizon = Hour(24), + interval = Hour(24), + ) + @test build!(model; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED + + # Second model with a different interval reuses the same system and auto-transforms. + model2 = DecisionModel( + template, + c_sys5; + optimizer = HiGHS_optimizer, + horizon = Hour(24), + interval = Hour(12), + ) + @test build!(model2; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + @test solve!(model2) == IOM.RunStatus.SUCCESSFULLY_FINALIZED +end + +@testset "RTS system shared across two intervals - build only" begin + sys_rts = PSB.build_system(PSISystems, "modified_RTS_GMLC_DA_sys") + # Clear any pre-existing transform so we can attach two fresh intervals. + PSY.transform_single_time_series!(sys_rts, Hour(24), Hour(24); delete_existing = true) + PSY.transform_single_time_series!(sys_rts, Hour(24), Hour(12); delete_existing = false) + + template = get_template_standard_uc_simulation() + set_network_model!(template, NetworkModel(CopperPlatePowerModel)) + + model_24h = DecisionModel( + template, + sys_rts; + name = "UC_24h", + optimizer = HiGHS_optimizer, + horizon = Hour(24), + interval = Hour(24), + ) + @test build!(model_24h; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + @test IOM.get_interval(model_24h) == Dates.Millisecond(Hour(24)) + + model_12h = DecisionModel( + template, + sys_rts; + name = "UC_12h", + optimizer = HiGHS_optimizer, + horizon = Hour(24), + interval = Hour(12), + ) + @test build!(model_12h; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + @test IOM.get_interval(model_12h) == Dates.Millisecond(Hour(12)) + + # Same underlying system, different intervals selected per model. + @test get_system(model_24h) === get_system(model_12h) +end + +@testset "Single interval system works without interval kwarg" begin + # Backward compatibility: existing single-interval systems work without changes + c_sys5 = PSB.build_system(PSITestSystems, "c_sys5") + template = get_thermal_dispatch_template_network() + model = DecisionModel(template, c_sys5; optimizer = HiGHS_optimizer) + @test build!(model; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED +end diff --git a/test/test_network_constructors.jl b/test/test_network_constructors.jl new file mode 100644 index 0000000..018de5c --- /dev/null +++ b/test/test_network_constructors.jl @@ -0,0 +1,74 @@ +function add_dummy_time_series_data!(sys) + # Attach dummy data so the problem builds: + dummy_data = Dict( + DateTime("2020-01-01T08:00:00") => [5.0, 6, 7, 7, 7], + DateTime("2020-01-01T08:30:00") => [9.0, 9, 9, 9, 8], + DateTime("2020-01-01T09:00:00") => [6.0, 6, 5, 5, 4], + ) + resolution = Dates.Minute(5) + dummy_forecast = Deterministic("max_active_power", dummy_data, resolution) + load = collect(get_components(StandardLoad, sys))[1] + add_time_series!(sys, load, dummy_forecast) + return sys +end + +# Regression test for https://github.com/NREL-Sienna/PowerSimulations.jl/issues/1594 +# Combines a NetworkModel with radial + degree-two reductions, a Line DeviceModel +# with a filter_function, and a request for FlowRateConstraint duals. Before the +# fix in src/devices_models/devices/common/add_constraint_dual.jl, the dual +# container was sized along PSY.get_name.(devices) — every device passing the +# filter — while the FlowRateConstraint container was sized along the +# post-reduction axis from get_branch_argument_constraint_axis. The resulting +# axis mismatch raised DimensionMismatch in process_duals during dual +# extraction. Building a model is enough to detect the regression: after the +# fix, axes(dual)[1] must equal axes(constraint)[1] for every meta. +@testset "FlowRateConstraint duals with branch filter and network reductions" begin + sys = build_system(PSITestSystems, "case11_network_reductions") + add_dummy_time_series_data!(sys) + nr = NetworkReduction[RadialReduction(), DegreeTwoReduction()] + ptdf = PTDF(sys; network_reductions = nr) + + template = OperationsProblemTemplate( + NetworkModel(PTDFPowerModel; + PTDF_matrix = ptdf, + duals = [CopperPlateBalanceConstraint], + reduce_radial_branches = PNM.has_radial_reduction(ptdf.network_reduction_data), + reduce_degree_two_branches = PNM.has_degree_two_reduction( + ptdf.network_reduction_data, + ), + use_slacks = false), + ) + # Mirror the filter shape from issue #1594: a voltage threshold that selects + # all lines in this all-230 kV system. The filter is registered (so the + # filter_function code path runs) but does not exclude any branch from a + # series chain, so reductions still drop lines from the constraint axis. + set_device_model!( + template, + DeviceModel( + Line, + StaticBranch; + duals = [FlowRateConstraint], + attributes = Dict( + "filter_function" => + x -> PSY.get_base_voltage(PSY.get_from(PSY.get_arc(x))) >= 230.0, + ), + ), + ) + set_device_model!(template, Transformer2W, StaticBranch) + ps_model = DecisionModel(template, sys; optimizer = HiGHS_optimizer) + @test build!(ps_model; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + + container = get_optimization_container(ps_model) + # The unfiltered Line set has 12 entries; full reduction leaves 6 entries + # in the constraint axis. The dual container must use the same 6 entries. + for meta in ("lb", "ub") + cons_key = ConstraintKey(FlowRateConstraint, Line, meta) + cons = get_constraint(container, cons_key) + dual = get_duals(container)[cons_key] + @test axes(dual)[1] == axes(cons)[1] + @test length(axes(cons)[1]) < + length(collect(get_components(Line, sys))) + @test "4-5-i_1" in axes(cons)[1] + end +end diff --git a/test/test_utils/model_checks.jl b/test/test_utils/model_checks.jl index 5b275ca..cad573e 100644 --- a/test/test_utils/model_checks.jl +++ b/test/test_utils/model_checks.jl @@ -102,27 +102,41 @@ function psi_checksolve_test(model::DecisionModel, status, expected_output, tol @test isapprox(obj_value, expected_output, atol = tol) end -# currently unused. we've removed get_system on Outputs, so would need to pass -# the system separately if we want to use this. - -#= -function psi_ptdf_lmps(outputs::OptimizationProblemOutputs, ptdf) - cp_duals = - read_dual(outputs, IOM.ConstraintKey(CopperPlateBalanceConstraint, PSY.System)) +function psi_ptdf_lmps(res::OptimizationProblemOutputs, ptdf) + cp_duals = read_dual( + res, + IOM.ConstraintKey(CopperPlateBalanceConstraint, PSY.System); + table_format = TableFormat.WIDE, + ) λ = Matrix{Float64}(cp_duals[:, propertynames(cp_duals) .!= :DateTime]) - flow_duals = read_dual(outputs, IOM.ConstraintKey(NetworkFlowConstraint, PSY.Line)) - μ = Matrix{Float64}(flow_duals[:, PNM.get_branch_ax(ptdf)]) - - buses = get_components(Bus, get_system(outputs)) + flow_ub = read_dual( + res, + IOM.ConstraintKey(FlowRateConstraint, PSY.Line, "ub"); + table_format = TableFormat.WIDE, + ) + flow_lb = read_dual( + res, + IOM.ConstraintKey(FlowRateConstraint, PSY.Line, "lb"); + table_format = TableFormat.WIDE, + ) + arcs = PNM.get_arc_axis(ptdf) + nr = PNM.get_network_reduction_data(ptdf) + branch_names = [PSY.get_name(nr.direct_branch_map[arc]) for arc in arcs] + μ = + Matrix{Float64}(flow_ub[:, branch_names]) .+ + Matrix{Float64}(flow_lb[:, branch_names]) + + buses = get_components(Bus, IOM.get_source_data(res)) lmps = OrderedDict() for bus in buses - lmps[get_name(bus)] = μ * ptdf[:, get_number(bus)] + bus_number = get_number(bus) + ptdf_col = [ptdf[arc, bus_number] for arc in arcs] + lmps[get_name(bus)] = μ * ptdf_col end lmp = λ .+ DataFrames.DataFrame(lmps) return lmp[!, sort(propertynames(lmp))] end -=# function check_variable_unbounded( model::DecisionModel,