diff --git a/Project.toml b/Project.toml index 6d6c548..a3d2d9c 100644 --- a/Project.toml +++ b/Project.toml @@ -20,6 +20,9 @@ Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" +[sources] +InfrastructureSystems = {url = "https://github.com/NREL-Sienna/InfrastructureSystems.jl", rev = "IS4"} + [weakdeps] PowerFlows = "94fada2c-0ca5-4b90-a1fb-4bc5b59ccfc7" diff --git a/src/PowerOperationsModels.jl b/src/PowerOperationsModels.jl index 5a7967e..c096a66 100644 --- a/src/PowerOperationsModels.jl +++ b/src/PowerOperationsModels.jl @@ -111,6 +111,8 @@ import InfrastructureOptimizationModels: add_proportional_cost!, add_proportional_cost_maybe_time_variant!, skip_proportional_cost, + # variable cost + add_variable_cost!, # Network model instantiation (POM extends for concrete network formulations) instantiate_network_model!, # Parameter addition (POM provides concrete implementations) @@ -520,6 +522,10 @@ export HVDCLosses export ConverterDCPower export ConverterCurrentDirection +# Load Variables +export ShiftUpActivePowerVariable +export ShiftDownActivePowerVariable + ######## Hydro Formulations ######## export HydroDispatchRunOfRiver export HydroDispatchRunOfRiverBudget @@ -666,6 +672,11 @@ export DurationConstraint export CommitmentConstraint export StartTypeConstraint export StartupTimeLimitTemperatureConstraint +export ShiftedActivePowerBalanceConstraint +export ShiftUpActivePowerVariableLimitsConstraint +export ShiftDownActivePowerVariableLimitsConstraint +export RealizedShiftedLoadMinimumBoundConstraint +export NonAnticipativityConstraint ################################################################################# # Exports - Expression Types (defined in core/expressions.jl) @@ -732,6 +743,7 @@ export ThermalSecurityConstrainedStandardUnitCommitment export StaticPowerLoad export PowerLoadInterruption export PowerLoadDispatch +export PowerLoadShift # Renewable Formulations export RenewableFullDispatch diff --git a/src/core/constraints.jl b/src/core/constraints.jl index e74ade2..2dc7a17 100644 --- a/src/core/constraints.jl +++ b/src/core/constraints.jl @@ -1047,3 +1047,56 @@ The specified constraints are formulated as: ``` """ struct StorageRegularizationConstraintDischarge <: ConstraintType end + +""" +Struct to create the constraint to balance shifted power over the user-defined time horizons. +For more information check the [`PowerLoadShift`](@ref) formulation. +The specified constraints are formulated as: +```math +\\sum_{t \\in \\text{time horizon}_k } p_t^\\text{shift,up} - p_t^\\text{shift,dn} = 0 , \\quad \\forall k \\text{ time horizons} +``` +""" +struct ShiftedActivePowerBalanceConstraint <: ConstraintType end + +""" +Struct to create the constraint to balance shifted power over the user-defined time horizons. +For more information check the [`PowerLoadShift`](@ref) formulation. +The specified constraints are formulated as: +```math +p_t^\\text{realized} \\ge 0.0 , \\quad \\forall k \\text{ time horizons} +``` +""" +struct RealizedShiftedLoadMinimumBoundConstraint <: ConstraintType end + +""" +Struct to create the non-anticipativity constraint for the [`PowerLoadShift`](@ref) formulation. +This enforces that shift up can only occur after an equal or greater amount of shift down has +already been committed, preventing the optimizer from shifting load up before it has been +shifted down. The constraint is formulated as: + +```math +\\sum_{\\tau=1}^{t} \\left( p_\\tau^\\text{shift,dn} - p_\\tau^\\text{shift,up} \\right) \\ge 0, +\\quad \\forall t \\in \\{1,\\dots,T\\} +``` +""" +struct NonAnticipativityConstraint <: ConstraintType end + +""" +Struct to create the constraint to limit shifted power active power between upper and lower bounds. +For more information check the [`PowerLoadShift`](@ref) formulation. +The specified constraints are formulated as: +```math +0 \\le p_t^\\text{shift, up} \\le P_t^\\text{upper}, \\quad \\forall t \\in \\{1,\\dots,T\\} +``` +""" +struct ShiftUpActivePowerVariableLimitsConstraint <: PowerVariableLimitsConstraint end + +""" +Struct to create the constraint to limit shifted power active power between upper and lower bounds. +For more information check the [`PowerLoadShift`](@ref) formulation. +The specified constraints are formulated as: +```math +0 \\le p_t^\\text{shift, dn} \\le P_t^\\text{lower}, \\quad \\forall t \\in \\{1,\\dots,T\\} +``` +""" +struct ShiftDownActivePowerVariableLimitsConstraint <: PowerVariableLimitsConstraint end diff --git a/src/core/expressions.jl b/src/core/expressions.jl index 6d678fd..7b3284a 100644 --- a/src/core/expressions.jl +++ b/src/core/expressions.jl @@ -11,6 +11,7 @@ struct ComponentReserveDownBalanceExpression <: ExpressionType end struct InterfaceTotalFlow <: ExpressionType end struct PTDFBranchFlow <: ExpressionType end struct PostContingencyNodalActivePowerDeployment <: PostContingencyExpressions end +struct RealizedShiftedLoad <: ExpressionType end ################################################################################# # Hydro Expressions @@ -100,6 +101,7 @@ struct ReserveDeploymentBalanceDownCharge <: StorageReserveChargeExpression end # Method extensions for output writing should_write_resulting_value(::Type{InterfaceTotalFlow}) = true should_write_resulting_value(::Type{PTDFBranchFlow}) = true +should_write_resulting_value(::Type{RealizedShiftedLoad}) = true should_write_resulting_value(::Type{HydroServedReserveUpExpression}) = true should_write_resulting_value(::Type{HydroServedReserveDownExpression}) = true @@ -114,3 +116,4 @@ convert_output_to_natural_units(::Type{InterfaceTotalFlow}) = true convert_output_to_natural_units(::Type{PostContingencyBranchFlow}) = true convert_output_to_natural_units(::Type{PostContingencyActivePowerGeneration}) = true convert_output_to_natural_units(::Type{PTDFBranchFlow}) = true +convert_output_to_natural_units(::Type{RealizedShiftedLoad}) = true diff --git a/src/core/formulations.jl b/src/core/formulations.jl index f212c1a..1159327 100644 --- a/src/core/formulations.jl +++ b/src/core/formulations.jl @@ -81,6 +81,11 @@ Formulation type to enable (continuous) load interruption dispatch """ struct PowerLoadDispatch <: AbstractControllablePowerLoadFormulation end +""" +Formulation type to enable load shifting +""" +struct PowerLoadShift <: AbstractControllablePowerLoadFormulation end + ############################ Regulation Device Formulations ################################ abstract type AbstractRegulationFormulation <: AbstractDeviceFormulation end struct ReserveLimitedRegulation <: AbstractRegulationFormulation end diff --git a/src/core/parameters.jl b/src/core/parameters.jl index a37ba09..83f6df9 100644 --- a/src/core/parameters.jl +++ b/src/core/parameters.jl @@ -30,6 +30,16 @@ Parameter to define active power in time series """ struct ActivePowerInTimeSeriesParameter <: TimeSeriesParameter end +""" +Parameter to define active power in time series +""" +struct ShiftUpActivePowerTimeSeriesParameter <: TimeSeriesParameter end + +""" +Parameter to define active power in time series +""" +struct ShiftDownActivePowerTimeSeriesParameter <: TimeSeriesParameter end + """ Parameter to define requirement time series """ @@ -238,6 +248,8 @@ convert_output_to_natural_units(::Type{ReactivePowerTimeSeriesParameter}) = true convert_output_to_natural_units(::Type{RequirementTimeSeriesParameter}) = true convert_output_to_natural_units(::Type{UpperBoundValueParameter}) = true convert_output_to_natural_units(::Type{LowerBoundValueParameter}) = true +convert_output_to_natural_units(::Type{ShiftUpActivePowerTimeSeriesParameter}) = true +convert_output_to_natural_units(::Type{ShiftDownActivePowerTimeSeriesParameter}) = true convert_output_to_natural_units(::Type{ReservoirLimitParameter}) = true convert_output_to_natural_units(::Type{ReservoirTargetParameter}) = true convert_output_to_natural_units(::Type{EnergyTargetTimeSeriesParameter}) = true diff --git a/src/core/variables.jl b/src/core/variables.jl index ef271e0..660ef5e 100644 --- a/src/core/variables.jl +++ b/src/core/variables.jl @@ -124,6 +124,22 @@ Docs abbreviation: ``\\theta`` """ struct VoltageAngle <: VariableType end +######################################### +###### Power Load Shift Variables ####### +######################################### + +""" +Struct to dispatch the creation of Shifted Up Active Power Variables +Docs abbreviation: ``p^\\text{shift,up}`` +""" +struct ShiftUpActivePowerVariable <: VariableType end + +""" +Struct to dispatch the creation of Shifted Down Active Power Variables +Docs abbreviation: ``p^\\text{shift,dn}`` +""" +struct ShiftDownActivePowerVariable <: VariableType end + ######################################### ####### DC Converter Variables ########## ######################################### diff --git a/src/static_injector_models/electric_loads.jl b/src/static_injector_models/electric_loads.jl index c8c779b..079a7dc 100644 --- a/src/static_injector_models/electric_loads.jl +++ b/src/static_injector_models/electric_loads.jl @@ -24,6 +24,21 @@ get_multiplier_value(::Type{<:TimeSeriesParameter}, d::PSY.ElectricLoad, ::Type{ get_multiplier_value(::Type{ReactivePowerTimeSeriesParameter}, d::PSY.ElectricLoad, ::Type{StaticPowerLoad}) = -1*PSY.get_max_reactive_power(d) get_multiplier_value(::Type{<:TimeSeriesParameter}, d::PSY.ElectricLoad, ::Type{<:AbstractControllablePowerLoadFormulation}) = PSY.get_max_active_power(d) +########################### ShiftablePowerLoad ##################################### + +get_variable_binary(::Type{ShiftUpActivePowerVariable}, ::Type{<:PSY.ElectricLoad}, ::Type{PowerLoadShift}) = false +get_variable_lower_bound(::Type{ShiftUpActivePowerVariable}, d::PSY.ElectricLoad, ::Type{PowerLoadShift}) = 0.0 +get_variable_upper_bound(::Type{ShiftUpActivePowerVariable}, d::PSY.ElectricLoad, ::Type{PowerLoadShift}) = nothing # Unbounded above by default, but can be limited by time series parameters + +get_variable_binary(::Type{ShiftDownActivePowerVariable}, ::Type{<:PSY.ElectricLoad}, ::Type{PowerLoadShift}) = false +get_variable_lower_bound(::Type{ShiftDownActivePowerVariable}, d::PSY.ElectricLoad, ::Type{PowerLoadShift}) = 0.0 +get_variable_upper_bound(::Type{ShiftDownActivePowerVariable}, d::PSY.ElectricLoad, ::Type{PowerLoadShift}) = nothing # Unbounded above by default, but can be limited by time series parameters + +variable_cost(cost::PSY.OperationalCost, ::Type{ShiftUpActivePowerVariable}, ::PSY.ElectricLoad, ::Type{<:AbstractControllablePowerLoadFormulation})=PSY.get_variable(cost) +variable_cost(cost::PSY.OperationalCost, ::Type{ShiftDownActivePowerVariable}, ::PSY.ElectricLoad, ::Type{<:AbstractControllablePowerLoadFormulation})=PSY.get_variable(cost) + +###################################################### + # To avoid ambiguity with default_interface_methods.jl: get_multiplier_value(::Type{<:AbstractPiecewiseLinearBreakpointParameter}, ::PSY.ElectricLoad, ::Type{StaticPowerLoad}) = 1.0 get_multiplier_value(::Type{<:AbstractPiecewiseLinearBreakpointParameter}, ::PSY.ElectricLoad, ::Type{<:AbstractControllablePowerLoadFormulation}) = 1.0 @@ -34,6 +49,8 @@ proportional_cost(cost::Nothing, ::Type{OnVariable}, ::PSY.ElectricLoad, ::Type{ proportional_cost(cost::PSY.OperationalCost, ::Type{OnVariable}, ::PSY.ElectricLoad, ::Type{<:AbstractControllablePowerLoadFormulation})=PSY.get_fixed(cost) objective_function_multiplier(::Type{<:VariableType}, ::Type{<:AbstractControllablePowerLoadFormulation})=OBJECTIVE_FUNCTION_NEGATIVE +objective_function_multiplier(::Type{ShiftUpActivePowerVariable}, ::Type{<:AbstractControllablePowerLoadFormulation})=OBJECTIVE_FUNCTION_NEGATIVE +objective_function_multiplier(::Type{ShiftDownActivePowerVariable}, ::Type{PowerLoadShift})=OBJECTIVE_FUNCTION_POSITIVE #! format: on @@ -81,6 +98,119 @@ function get_default_time_series_names( return Dict{Type{<:TimeSeriesParameter}, String}() end +function get_default_time_series_names( + ::Type{<:PSY.ShiftablePowerLoad}, + ::Type{PowerLoadShift}, +) + return Dict{Type{<:TimeSeriesParameter}, String}( + ActivePowerTimeSeriesParameter => "max_active_power", + ReactivePowerTimeSeriesParameter => "max_active_power", + ShiftUpActivePowerTimeSeriesParameter => "shift_up_max_active_power", + ShiftDownActivePowerTimeSeriesParameter => "shift_down_max_active_power", + ) +end + +####################### Expressions ######################### + +function add_expressions!( + container::OptimizationContainer, + ::Type{T}, + devices::U, + model::DeviceModel{D, W}, +) where { + T <: RealizedShiftedLoad, + U <: Union{Vector{D}, IS.FlattenIteratorWrapper{D}}, + W <: PowerLoadShift, +} where {D <: PSY.ShiftablePowerLoad} + time_steps = get_time_steps(container) + names = PSY.get_name.(devices) + expression = add_expression_container!(container, T, D, names, time_steps) + shift_up = get_variable(container, ShiftUpActivePowerVariable, D) + shift_down = get_variable(container, ShiftDownActivePowerVariable, D) + param_container = get_parameter(container, ActivePowerTimeSeriesParameter, D) + multiplier = get_multiplier_array(param_container) + for t in time_steps, d in devices + name = PSY.get_name(d) + expression[name, t] = + get_parameter_column_refs(param_container, name)[t] * multiplier[name, t] + + shift_up[name, t] - shift_down[name, t] + end + return +end + +function add_to_expression!( + container::OptimizationContainer, + ::Type{T}, + ::Type{U}, + devices::IS.FlattenIteratorWrapper{V}, + device_model::DeviceModel{V, W}, + network_model::NetworkModel{CopperPlatePowerModel}, +) where { + T <: ActivePowerBalance, + U <: RealizedShiftedLoad, + V <: PSY.StaticInjection, + W <: AbstractDeviceFormulation, +} + realized_load = get_expression(container, U, V) + expression = get_expression(container, T, PSY.System) + for d in devices + device_bus = PSY.get_bus(d) + ref_bus = get_reference_bus(network_model, device_bus) + name = PSY.get_name(d) + for t in get_time_steps(container) + JuMP.add_to_expression!( + expression[ref_bus, t], + -1.0, # Realized load enter negative to the balance + realized_load[name, t], + ) + end + end + return +end + +""" +Electric Load implementation to add parameters to PTDF ActivePowerBalance expressions +""" +function add_to_expression!( + container::OptimizationContainer, + ::Type{T}, + ::Type{U}, + devices::IS.FlattenIteratorWrapper{V}, + device_model::DeviceModel{V, W}, + network_model::NetworkModel{X}, +) where { + T <: ActivePowerBalance, + U <: RealizedShiftedLoad, + V <: PSY.ShiftablePowerLoad, + W <: PowerLoadShift, + X <: AbstractPTDFModel, +} + realized_load = get_expression(container, U, V) + sys_expr = get_expression(container, T, _system_expression_type(X)) + nodal_expr = get_expression(container, T, PSY.ACBus) + network_reduction = get_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(network_reduction, bus_no_) + ref_index = _ref_index(network_model, device_bus) + for t in get_time_steps(container) + JuMP.add_to_expression!( + sys_expr[ref_index, t], + -1.0, # Realized load enter negative to the balance + realized_load[name, t], + ) + JuMP.add_to_expression!( + nodal_expr[bus_no, t], + -1.0, # Realized load enter negative to the balance + realized_load[name, t], + ) + end + end + return +end + ####################################### Reactive Power Constraints ######################### """ Reactive Power Constraints on Controllable Loads Assume Constant power_factor @@ -180,6 +310,198 @@ function add_constraints!( return end +function add_constraints!( + container::OptimizationContainer, + T::Type{ShiftedActivePowerBalanceConstraint}, + devices::IS.FlattenIteratorWrapper{V}, + model::DeviceModel{V, W}, + ::NetworkModel{X}, +) where {V <: PSY.ShiftablePowerLoad, W <: PowerLoadShift, X <: PM.AbstractPowerModel} + time_steps = get_time_steps(container) + time_steps_end = time_steps[end] + # Keep this container 2D (name, terminal-time marker) to match standard indexing patterns. + constraint = add_constraints_container!( + container, + T, + V, + PSY.get_name.(devices), + [time_steps_end], + ) + up_variable = get_variable(container, ShiftUpActivePowerVariable, V) + down_variable = get_variable(container, ShiftDownActivePowerVariable, V) + jump_model = get_jump_model(container) + for d in devices + name = PSY.get_name(d) + constraint[name, time_steps_end] = + JuMP.@constraint( + jump_model, + sum(up_variable[name, t] - down_variable[name, t] for t in time_steps) == + 0.0 + ) + end + additional_balance_interval = get_attribute(model, "additional_balance_interval") + if !isnothing(additional_balance_interval) + if !(additional_balance_interval isa Dates.Period) + throw( + IS.InvalidValue( + "The additional_balance_interval attribute must be a Dates.Period, got $(typeof(additional_balance_interval)).", + ), + ) + end + interval_ms = Dates.Millisecond(additional_balance_interval).value + if interval_ms <= 0 + throw( + IS.InvalidValue( + "The additional_balance_interval attribute must be greater than zero.", + ), + ) + end + + resolution = get_resolution(container) + resolution_ms = Dates.Millisecond(resolution).value + + if interval_ms % resolution_ms != 0 + throw( + IS.InvalidValue( + "The additional_balance_interval attribute must be an integer multiple of model resolution (interval_ms = $(interval_ms), resolution_ms = $(resolution_ms)).", + ), + ) + end + + interval_length = interval_ms ÷ resolution_ms + if interval_length > length(time_steps) + throw( + IS.InvalidValue( + "The additional_balance_interval attribute must be less than or equal to the optimization horizon.", + ), + ) + end + + interval_ranges = [ + start_idx:min(start_idx + interval_length - 1, length(time_steps)) for + start_idx in 1:interval_length:length(time_steps) + ] + interval_end_steps = [time_steps[last(interval_range)] for interval_range in interval_ranges] + constraint_aux = add_constraints_container!( + container, + T, + V, + PSY.get_name.(devices), + interval_end_steps; + meta = "additional", + ) + for d in devices + name = PSY.get_name(d) + for interval_range in interval_ranges + end_time_step = time_steps[last(interval_range)] + constraint_aux[name, end_time_step] = JuMP.@constraint( + container.JuMPmodel, + sum( + up_variable[name, t] - down_variable[name, t] for + t in time_steps[interval_range] + ) == 0.0 + ) + end + end + end + return +end + +function add_constraints!( + container::OptimizationContainer, + T::Type{RealizedShiftedLoadMinimumBoundConstraint}, + U::Type{<:ExpressionType}, + devices::IS.FlattenIteratorWrapper{V}, + ::DeviceModel{V, W}, + ::NetworkModel{X}, +) where {V <: PSY.ShiftablePowerLoad, W <: PowerLoadShift, X <: PM.AbstractPowerModel} + time_steps = get_time_steps(container) + constraint = add_constraints_container!( + container, + T, + V, + PSY.get_name.(devices), + time_steps, + ) + realized_load = get_expression(container, U, V) + jump_model = get_jump_model(container) + for d in devices, t in time_steps + name = PSY.get_name(d) + constraint[name, t] = JuMP.@constraint(jump_model, realized_load[name, t] >= 0.0) + end + return +end + +function add_constraints!( + container::OptimizationContainer, + T::Type{NonAnticipativityConstraint}, + devices::IS.FlattenIteratorWrapper{V}, + ::DeviceModel{V, W}, + ::NetworkModel{X}, +) where {V <: PSY.ShiftablePowerLoad, W <: PowerLoadShift, X <: PM.AbstractPowerModel} + time_steps = get_time_steps(container) + constraint = add_constraints_container!( + container, + T, + V, + PSY.get_name.(devices), + time_steps, + ) + up_variable = get_variable(container, ShiftUpActivePowerVariable, V) + down_variable = get_variable(container, ShiftDownActivePowerVariable, V) + jump_model = get_jump_model(container) + for d in devices + name = PSY.get_name(d) + for t in time_steps + constraint[name, t] = JuMP.@constraint( + jump_model, + sum(down_variable[name, τ] - up_variable[name, τ] for τ in 1:t) >= 0.0 + ) + end + end + return +end + +function add_constraints!( + container::OptimizationContainer, + ::Type{ShiftUpActivePowerVariableLimitsConstraint}, + U::Type{<:VariableType}, + devices::IS.FlattenIteratorWrapper{V}, + model::DeviceModel{V, W}, + ::NetworkModel{X}, +) where {V <: PSY.ShiftablePowerLoad, W <: PowerLoadShift, X <: PM.AbstractPowerModel} + add_parameterized_upper_bound_range_constraints( + container, + ShiftUpActivePowerVariableLimitsConstraint, + U, + ShiftUpActivePowerTimeSeriesParameter, + devices, + model, + X, + ) + return +end + +function add_constraints!( + container::OptimizationContainer, + ::Type{ShiftDownActivePowerVariableLimitsConstraint}, + U::Type{<:VariableType}, + devices::IS.FlattenIteratorWrapper{V}, + model::DeviceModel{V, W}, + ::NetworkModel{X}, +) where {V <: PSY.ShiftablePowerLoad, W <: PowerLoadShift, X <: PM.AbstractPowerModel} + add_parameterized_upper_bound_range_constraints( + container, + ShiftDownActivePowerVariableLimitsConstraint, + U, + ShiftDownActivePowerTimeSeriesParameter, + devices, + model, + X, + ) + return +end + ############################## FormulationControllable Load Cost ########################### function add_to_objective_function!( container::OptimizationContainer, @@ -267,3 +589,30 @@ function proportional_cost( ) end end + +########## PowerLoadShift Formulation Costs ############### + +function objective_function!( + container::OptimizationContainer, + devices::IS.FlattenIteratorWrapper{T}, + ::DeviceModel{T, U}, + ::Type{<:PM.AbstractPowerModel}, +) where {T <: PSY.ShiftablePowerLoad, U <: PowerLoadShift} + add_variable_cost!(container, ShiftUpActivePowerVariable, devices, U) + add_variable_cost!(container, ShiftDownActivePowerVariable, devices, U) + return +end + +### Special Method to skip VOM cost on ShiftUpActivePowerVariable ### +function add_variable_cost!( + container::OptimizationContainer, + ::Type{U}, + devices::IS.FlattenIteratorWrapper{T}, + ::Type{V}, +) where {T <: PSY.ShiftablePowerLoad, U <: ShiftUpActivePowerVariable, V <: PowerLoadShift} + for d in devices + op_cost_data = PSY.get_operation_cost(d) + add_variable_cost_to_objective!(container, U, d, op_cost_data, V) + end + return +end diff --git a/src/static_injector_models/load_constructor.jl b/src/static_injector_models/load_constructor.jl index cacaa74..56b5dfc 100644 --- a/src/static_injector_models/load_constructor.jl +++ b/src/static_injector_models/load_constructor.jl @@ -515,3 +515,223 @@ function construct_device!( construct_device!(container, sys, ccs, new_model, network_model) return end + +function construct_device!( + container::OptimizationContainer, + sys::PSY.System, + ::ArgumentConstructStage, + model::DeviceModel{PSY.ShiftablePowerLoad, PowerLoadShift}, + network_model::NetworkModel{<:PM.AbstractPowerModel}, +) + devices = get_available_components(model, sys) + + add_variables!(container, ShiftUpActivePowerVariable, devices, PowerLoadShift) + add_variables!(container, ShiftDownActivePowerVariable, devices, PowerLoadShift) + add_variables!(container, ReactivePowerVariable, devices, PowerLoadShift) + + process_market_bid_parameters!(container, devices, model, false, true) + + if haskey(get_time_series_names(model), ActivePowerTimeSeriesParameter) + add_parameters!(container, ActivePowerTimeSeriesParameter, devices, model) + end + if haskey(get_time_series_names(model), ShiftUpActivePowerTimeSeriesParameter) + add_parameters!(container, ShiftUpActivePowerTimeSeriesParameter, devices, model) + end + if haskey(get_time_series_names(model), ShiftDownActivePowerTimeSeriesParameter) + add_parameters!(container, ShiftDownActivePowerTimeSeriesParameter, devices, model) + end + + # Add realized load expression + add_expressions!(container, RealizedShiftedLoad, devices, model) + + # Add Parameters to expressions + add_to_expression!( + container, + ActivePowerBalance, + RealizedShiftedLoad, + devices, + model, + network_model, + ) + add_to_expression!( + container, + ReactivePowerBalance, + ReactivePowerVariable, + devices, + model, + network_model, + ) + + add_expressions!(container, ProductionCostExpression, devices, model) + add_event_arguments!(container, devices, model, network_model) + return +end + +function construct_device!( + container::OptimizationContainer, + sys::PSY.System, + ::ModelConstructStage, + model::DeviceModel{PSY.ShiftablePowerLoad, PowerLoadShift}, + network_model::NetworkModel{<:PM.AbstractPowerModel}, +) + devices = + get_available_components(model, + sys, + ) + + add_constraints!( + container, + ShiftedActivePowerBalanceConstraint, + devices, + model, + network_model, + ) + add_constraints!( + container, + RealizedShiftedLoadMinimumBoundConstraint, + RealizedShiftedLoad, + devices, + model, + network_model, + ) + add_constraints!( + container, + ShiftUpActivePowerVariableLimitsConstraint, + ShiftUpActivePowerVariable, + devices, + model, + network_model, + ) + add_constraints!( + container, + ShiftDownActivePowerVariableLimitsConstraint, + ShiftDownActivePowerVariable, + devices, + model, + network_model, + ) + add_constraints!( + container, + ReactivePowerVariableLimitsConstraint, + ReactivePowerVariable, + devices, + model, + network_model, + ) + add_constraints!( + container, + NonAnticipativityConstraint, + devices, + model, + network_model, + ) + + add_feedforward_constraints!(container, model, devices) + + objective_function!(container, devices, model, get_network_formulation(network_model)) + add_event_constraints!(container, devices, model, network_model) + add_constraint_dual!(container, sys, model) + return +end + +# AbstractActivePowerModel + PowerLoadShift device model +function construct_device!( + container::OptimizationContainer, + sys::PSY.System, + ::ArgumentConstructStage, + model::DeviceModel{PSY.ShiftablePowerLoad, PowerLoadShift}, + network_model::NetworkModel{<:PM.AbstractActivePowerModel}, +) + devices = get_available_components(model, sys) + + add_variables!(container, ShiftUpActivePowerVariable, devices, PowerLoadShift) + add_variables!(container, ShiftDownActivePowerVariable, devices, PowerLoadShift) + + process_market_bid_parameters!(container, devices, model, false, true) + + if haskey(get_time_series_names(model), ActivePowerTimeSeriesParameter) + add_parameters!(container, ActivePowerTimeSeriesParameter, devices, model) + end + if haskey(get_time_series_names(model), ShiftUpActivePowerTimeSeriesParameter) + add_parameters!(container, ShiftUpActivePowerTimeSeriesParameter, devices, model) + end + if haskey(get_time_series_names(model), ShiftDownActivePowerTimeSeriesParameter) + add_parameters!(container, ShiftDownActivePowerTimeSeriesParameter, devices, model) + end + + # Add realized load expression + add_expressions!(container, RealizedShiftedLoad, devices, model) + + # Add Parameters to expressions + add_to_expression!( + container, + ActivePowerBalance, + RealizedShiftedLoad, + devices, + model, + network_model, + ) + + add_expressions!(container, ProductionCostExpression, devices, model) + add_event_arguments!(container, devices, model, network_model) + return +end + +function construct_device!( + container::OptimizationContainer, + sys::PSY.System, + ::ModelConstructStage, + model::DeviceModel{PSY.ShiftablePowerLoad, PowerLoadShift}, + network_model::NetworkModel{<:PM.AbstractActivePowerModel}, +) + devices = + get_available_components(model, + sys, + ) + + add_constraints!( + container, + ShiftedActivePowerBalanceConstraint, + devices, + model, + network_model, + ) + add_constraints!( + container, + RealizedShiftedLoadMinimumBoundConstraint, + RealizedShiftedLoad, + devices, + model, + network_model, + ) + add_constraints!( + container, + ShiftUpActivePowerVariableLimitsConstraint, + ShiftUpActivePowerVariable, + devices, + model, + network_model, + ) + add_constraints!( + container, + ShiftDownActivePowerVariableLimitsConstraint, + ShiftDownActivePowerVariable, + devices, + model, + network_model, + ) + add_constraints!( + container, + NonAnticipativityConstraint, + devices, + model, + network_model, + ) + + add_feedforward_constraints!(container, model, devices) + + objective_function!(container, devices, model, get_network_formulation(network_model)) + add_event_constraints!(container, devices, model, network_model) + add_constraint_dual!(container, sys, model) + return +end diff --git a/test/test_device_load_constructors.jl b/test/test_device_load_constructors.jl index 38309fa..b96c831 100644 --- a/test/test_device_load_constructors.jl +++ b/test/test_device_load_constructors.jl @@ -239,3 +239,109 @@ end @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED end end + +@testset "PowerLoadShift with NonAnticipativityConstraint" begin + c_sys5_il = + PSB.build_system(PSITestSystems, "c_sys5_il"; add_single_time_series = true) + il_load = first(PSY.get_components(InterruptiblePowerLoad, c_sys5_il)) + + shiftable_load = ShiftablePowerLoad(; + name = "shiftable_load", + available = true, + bus = PSY.get_bus(il_load), + active_power = PSY.get_active_power(il_load), + active_power_limits = (min = 0.0, max = PSY.get_active_power(il_load)), + reactive_power = PSY.get_reactive_power(il_load), + max_active_power = PSY.get_max_active_power(il_load), + max_reactive_power = PSY.get_max_reactive_power(il_load), + base_power = PSY.get_base_power(il_load), + load_balance_time_horizon = 1, + operation_cost = LoadCost(; + variable = CostCurve( + LinearCurve(0.0), + UnitSystem.NATURAL_UNITS, + LinearCurve(1.0), + ), + fixed = 0.0, + ), + ) + PSY.add_component!(c_sys5_il, shiftable_load) + PSY.set_available!(il_load, false) + PSY.copy_time_series!(shiftable_load, il_load) + + tstamps = TimeSeries.timestamp( + PSY.get_time_series_array(SingleTimeSeries, shiftable_load, "max_active_power"), + ) + n = length(tstamps) + up_vals = ones(n) + down_vals = ones(n) + + PSY.add_time_series!( + c_sys5_il, + shiftable_load, + SingleTimeSeries( + "shift_up_max_active_power", + TimeArray(tstamps, up_vals); + scaling_factor_multiplier = PSY.get_max_active_power, + ), + ) + PSY.add_time_series!( + c_sys5_il, + shiftable_load, + SingleTimeSeries( + "shift_down_max_active_power", + TimeArray(tstamps, down_vals); + scaling_factor_multiplier = PSY.get_max_active_power, + ), + ) + + PSY.transform_single_time_series!(c_sys5_il, Hour(24), Hour(24)) + + template = OperationsProblemTemplate( + NetworkModel( + CopperPlatePowerModel; + duals = [CopperPlateBalanceConstraint], + ), + ) + set_device_model!(template, ThermalStandard, ThermalBasicUnitCommitment) + set_device_model!(template, PowerLoad, StaticPowerLoad) + set_device_model!( + template, + DeviceModel( + ShiftablePowerLoad, + PowerLoadShift; + attributes = Dict{String, Any}("additional_balance_interval" => Hour(12)), + ), + ) + + model = DecisionModel( + template, + c_sys5_il; + name = "UC_shiftable", + store_variable_names = true, + optimizer = HiGHS_optimizer, + ) + + @test build!(model; output_dir = mktempdir(; cleanup = true)) == + ModelBuildStatus.BUILT + + @test solve!(model) == RunStatus.SUCCESSFULLY_FINALIZED + + results = OptimizationProblemOutputs(model) + up = read_variable( + results, + "ShiftUpActivePowerVariable__ShiftablePowerLoad"; + table_format = TableFormat.WIDE, + ) + dn = read_variable( + results, + "ShiftDownActivePowerVariable__ShiftablePowerLoad"; + table_format = TableFormat.WIDE, + ) + + # Verify the non-anticipativity constraint holds in the solution: + # the running sum of (shift_down - shift_up) must be >= 0 at every time step. + @test all( + cumsum(dn[!, "shiftable_load"] .- up[!, "shiftable_load"]) .>= -1e-6, + ) +end