diff --git a/Project.toml b/Project.toml index 1aaa99b..c7ad013 100644 --- a/Project.toml +++ b/Project.toml @@ -4,14 +4,14 @@ authors = ["Jose Daniel Lara, Rodrigo Henriquez-Auba, Sourabh Dalvi"] version = "0.3.0" [deps] -Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" +Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" InfrastructureSystems = "2cd47ed4-ca9b-11e9-27f2-ab636a7671f1" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" +MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" PowerSimulations = "e690365d-45e2-57bb-ac84-44ba829e73c4" PowerSystems = "bcd98974-b02a-5e2f-9ee0-a103f5c450dd" Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b" @@ -20,8 +20,8 @@ Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b" DataStructures = "~0.18" DocStringExtensions = "~0.8, ~0.9" InfrastructureSystems = "1" -MathOptInterface = "1" JuMP = "1" +MathOptInterface = "1" PowerSimulations = "^0.23" PowerSystems = "3" julia = "^1.6" diff --git a/src/StorageSystemsSimulations.jl b/src/StorageSystemsSimulations.jl index b7cb766..88c34b9 100644 --- a/src/StorageSystemsSimulations.jl +++ b/src/StorageSystemsSimulations.jl @@ -21,6 +21,10 @@ export ReserveCoverageConstraint export ReserveCoverageConstraintEndOfPeriod export StorageTotalReserveConstraint +# FF +export EnergyTargetFeedforward +export EnergyLimitFeedforward + ################################################################################# # Imports import Logging @@ -67,6 +71,7 @@ include("core/formulations.jl") include("core/variables.jl") include("core/constraints.jl") include("core/expressions.jl") +include("core/parameters.jl") include("core/initial_conditions.jl") include("core/feedforward.jl") diff --git a/src/core/feedforward.jl b/src/core/feedforward.jl index e8db0c1..6a87110 100644 --- a/src/core/feedforward.jl +++ b/src/core/feedforward.jl @@ -1,8 +1,48 @@ +""" +Adds a constraint to enforce a minimum energy level target with a slack variable associated witha penalty term. +""" +struct EnergyTargetFeedforward <: PSI.AbstractAffectFeedforward + optimization_container_key::PSI.OptimizationContainerKey + affected_values::Vector{<:PSI.OptimizationContainerKey} + target_period::Int + penalty_cost::Float64 + function EnergyTargetFeedforward(; + component_type::Type{<:PSY.Component}, + source::Type{T}, + affected_values::Vector{DataType}, + target_period::Int, + penalty_cost::Float64, + meta=PSI.CONTAINER_KEY_EMPTY_META, + ) where {T} + values_vector = Vector{PSI.VariableKey}(undef, length(affected_values)) + for (ix, v) in enumerate(affected_values) + if v <: PSI.VariableType + values_vector[ix] = + PSI.get_optimization_container_key(v(), component_type, meta) + else + error( + "EnergyTargetFeedforward is only compatible with VariableType or ParamterType affected values", + ) + end + end + new( + PSI.get_optimization_container_key(T(), component_type, meta), + values_vector, + target_period, + penalty_cost, + ) + end +end + +PSI.get_default_parameter_type(::EnergyTargetFeedforward, _) = EnergyTargetParameter +PSI.get_optimization_container_key(ff::EnergyTargetFeedforward) = + ff.optimization_container_key + function PSI._add_feedforward_arguments!( container::PSI.OptimizationContainer, model::PSI.DeviceModel, devices::IS.FlattenIteratorWrapper{T}, - ff::PSI.EnergyTargetFeedforward, + ff::EnergyTargetFeedforward, ) where {T <: PSY.Storage} parameter_type = PSI.get_default_parameter_type(ff, T) PSI.add_parameters!(container, parameter_type, ff, model, devices) @@ -43,7 +83,7 @@ function PSI.add_feedforward_constraints!( container::PSI.OptimizationContainer, ::PSI.DeviceModel{T, U}, devices::IS.FlattenIteratorWrapper{T}, - ff::PSI.EnergyTargetFeedforward, + ff::EnergyTargetFeedforward, ) where {T <: PSY.Storage, U <: AbstractStorageFormulation} time_steps = PSI.get_time_steps(container) parameter_type = PSI.get_default_parameter_type(ff, T) @@ -83,58 +123,196 @@ function PSI.add_feedforward_constraints!( return end -#= -# TODO: Check if this is should be restored -function PSI._add_variable_cost_to_objective!( +""" +Adds a constraint to limit the sum of a variable over the number of periods to the source value +""" +struct EnergyLimitFeedforward <: PSI.AbstractAffectFeedforward + optimization_container_key::PSI.OptimizationContainerKey + affected_values::Vector{<:PSI.OptimizationContainerKey} + number_of_periods::Int + function EnergyLimitFeedforward(; + component_type::Type{<:PSY.Component}, + source::Type{T}, + affected_values::Vector{DataType}, + number_of_periods::Int, + meta=PSI.CONTAINER_KEY_EMPTY_META, + ) where {T} + values_vector = Vector{PSI.VariableKey}(undef, length(affected_values)) + for (ix, v) in enumerate(affected_values) + if v <: PSI.VariableType + values_vector[ix] = + PSI.get_optimization_container_key(v(), component_type, meta) + else + error( + "EnergyLimitFeedforward is only compatible with VariableType or ParamterType affected values", + ) + end + end + new( + PSI.get_optimization_container_key(T(), component_type, meta), + values_vector, + number_of_periods, + ) + end +end + +PSI.get_default_parameter_type(::EnergyLimitFeedforward, _) = EnergyLimitParameter +PSI.get_optimization_container_key(ff) = ff.optimization_container_key +get_number_of_periods(ff) = ff.number_of_periods + +@doc raw""" + add_feedforward_constraints(container::OptimizationContainer, + cons_name::Symbol, + param_reference, + var_key::VariableKey) + +Constructs a parameterized integral limit constraint to implement feedforward from other models. +The Parameters are initialized using the upper boundary values of the provided variables. + + +``` sum(variable[var_name, t] for t in 1:affected_periods)/affected_periods <= param_reference[var_name] ``` + +# LaTeX + +`` \sum_{t} x \leq param^{max}`` + +# Arguments +* container::OptimizationContainer : the optimization_container model built in PowerSimulations +* model::DeviceModel : the device model +* devices::IS.FlattenIteratorWrapper{T} : list of devices +* ff::FixValueFeedforward : a instance of the FixValue Feedforward +""" +function PSI.add_feedforward_constraints!( container::PSI.OptimizationContainer, - ::T, - component::U, - op_cost::PSY.MarketBidCost, - ::V, -) where {T <: PSI.EnergySurplusVariable, U <: PSY.Storage, V <: EnergyValueCurve} - component_name = PSY.get_name(component) - @debug "Market Bid" _group = PSI.LOG_GROUP_COST_FUNCTIONS component_name + ::PSI.DeviceModel, + devices::IS.FlattenIteratorWrapper{T}, + ff::EnergyLimitFeedforward, +) where {T <: PSY.Component} time_steps = PSI.get_time_steps(container) - initial_time = PSI.get_initial_time(container) - variable_cost_forecast = PSY.get_variable_cost( - component, - op_cost; - start_time=initial_time, - len=length(time_steps), - ) - variable_cost_forecast_values = TimeSeries.values(variable_cost_forecast) - parameter_container = PSI._get_cost_function_parameter_container( - container, - PSI.CostFunctionParameter(), - component, - T(), - V(), - eltype(variable_cost_forecast_values), - ) - pwl_cost_expressions = - PSI._add_pwl_term!(container, component, variable_cost_forecast_values, T(), V()) - jump_model = PSI.get_jump_model(container) - for t in time_steps - PSI.set_parameter!( - parameter_container, - jump_model, - PSY.get_cost(variable_cost_forecast_values[t]), - # Using 1.0 here since we want to reuse the existing code that adds the mulitpler - # of base power times the time delta. - 1.0, - component_name, - t, - ) - PSI.add_to_expression!( + parameter_type = PSI.get_default_parameter_type(ff, T) + param = PSI.get_parameter_array(container, parameter_type(), T) + multiplier = PSI.get_parameter_multiplier_array(container, parameter_type(), T) + affected_periods = get_number_of_periods(ff) + for var in PSI.get_affected_values(ff) + variable = PSI.get_variable(container, var) + set_name, set_time = JuMP.axes(variable) + IS.@assert_op set_name == [PSY.get_name(d) for d in devices] + IS.@assert_op set_time == time_steps + + if affected_periods > set_time[end] + error( + "The number of affected periods $affected_periods is larger than the periods available $(set_time[end])", + ) + end + no_trenches = set_time[end] รท affected_periods + var_type = PSI.get_entry_type(var) + con_ub = PSI.add_constraints_container!( container, - PSI.ProductionCostExpression, - pwl_cost_expressions[t], - component, - t, + PSI.FeedforwardIntegralLimitConstraint(), + T, + set_name, + 1:no_trenches; + meta="$(var_type)integral", ) - PSI.add_to_objective_variant_expression!(container, pwl_cost_expressions[t]) + + for name in set_name, i in 1:no_trenches + con_ub[name, i] = JuMP.@constraint( + container.JuMPmodel, + sum( + variable[name, t] for + t in (1 + (i - 1) * affected_periods):(i * affected_periods) + ) <= sum( + param[name, t] * multiplier[name, t] for + t in (1 + (i - 1) * affected_periods):(i * affected_periods) + ) + ) + end end + return +end +# TODO: It also needs the add parameters code + +function PSI.update_parameter_values!( + model::PSI.OperationModel, + key::PSI.ParameterKey{T, U}, + input::PSI.DatasetContainer{PSI.InMemoryDataset}, +) where {T <: EnergyLimitParameter, U <: PSY.Generator} + # Enable again for detailed debugging + # TimerOutputs.@timeit RUN_SIMULATION_TIMER "$T $U Parameter Update" begin + optimization_container = PSI.get_optimization_container(model) + # Note: Do not instantite a new key here because it might not match the param keys in the container + # if the keys have strings in the meta fields + parameter_array = PSI.get_parameter_array(optimization_container, key) + parameter_attributes = PSI.get_parameter_attributes(optimization_container, key) + internal = PSI.get_internal(model) + execution_count = internal.execution_count + current_time = PSI.get_current_time(model) + state_values = + PSI.get_dataset_values(input, PSI.get_attribute_key(parameter_attributes)) + component_names, time = axes(parameter_array) + resolution = PSI.get_resolution(model) + interval_time_steps = + Int(PSI.get_interval(model.internal.store_parameters) / resolution) + state_data = PSI.get_dataset(input, PSI.get_attribute_key(parameter_attributes)) + state_timestamps = state_data.timestamps + max_state_index = PSI.get_num_rows(state_data) + + state_data_index = PSI.find_timestamp_index(state_timestamps, current_time) + sim_timestamps = range(current_time; step=resolution, length=time[end]) + old_parameter_values = jump_value.(parameter_array) + # The current method uses older parameter values because when passing the energy output from one stage + # to the next, the aux variable values gets over-written by the lower level model after its solve. + # This approach is a temporary hack and will be replaced in future versions. + for t in time + timestamp_ix = min(max_state_index, state_data_index + 1) + @debug "parameter horizon is over the step" max_state_index > state_data_index + 1 + if state_timestamps[timestamp_ix] <= sim_timestamps[t] + state_data_index = timestamp_ix + end + for name in component_names + # the if statement checks if its the first solve of the model and uses the values stored in the state + # and for subsequent solves uses the state data to update the parameter values for the last set of time periods + # that are equal to the length of the interval i.e. the time periods that dont overlap between each solves. + if execution_count == 0 || t > time[end] - interval_time_steps + # Pass indices in this way since JuMP DenseAxisArray don't support view() + state_value = state_values[name, state_data_index] + if !isfinite(state_value) + error( + "The value for the system state used in $(encode_key_as_string(key)) is not a finite value $(state_value) \ + This is commonly caused by referencing a state value at a time when such decision hasn't been made. \ + Consider reviewing your models' horizon and interval definitions", + ) + end + PSI._set_param_value!(parameter_array, state_value, name, t) + else + # Currently the update method relies on using older parameter values of the EnergyLimitParameter + # to update the parameter for overlapping periods between solves i.e. we ingoring the parameter values + # in the model interval time periods. + state_value = state_values[name, state_data_index] + if !isfinite(state_value) + error( + "The value for the system state used in $(encode_key_as_string(key)) is not a finite value $(state_value) \ + This is commonly caused by referencing a state value at a time when such decision hasn't been made. \ + Consider reviewing your models' horizon and interval definitions", + ) + end + PSI._set_param_value!( + parameter_array, + old_parameter_values[name, t + interval_time_steps], + name, + t, + ) + end + end + end + + IS.@record :execution PSI.ParameterUpdateEvent( + T, + U, + parameter_attributes, + PSI.get_current_timestamp(model), + PSI.get_name(model), + ) return end -=# diff --git a/src/core/parameters.jl b/src/core/parameters.jl new file mode 100644 index 0000000..f21f820 --- /dev/null +++ b/src/core/parameters.jl @@ -0,0 +1,10 @@ +""" +Parameter to define energy limit +""" +struct EnergyLimitParameter <: PSI.VariableValueParameter end +# TODO: Check if EnergyTargetParameter and EnergyLimitParameter should be removed +# This affects feedforwards that can break if not defined +struct EnergyTargetParameter <: PSI.VariableValueParameter end + +convert_result_to_natural_units(::Type{EnergyLimitParameter}) = true +convert_result_to_natural_units(::Type{EnergyTargetParameter}) = true diff --git a/src/storage_models.jl b/src/storage_models.jl index 702b6f2..f191a52 100644 --- a/src/storage_models.jl +++ b/src/storage_models.jl @@ -76,8 +76,8 @@ PSI.variable_cost(cost::PSY.StorageManagementCost, ::PSI.ActivePowerInVariable, ######################## Parameters ################################################## -PSI.get_parameter_multiplier(::PSI.EnergyTargetParameter, ::PSY.Storage, ::AbstractStorageFormulation) = 1.0 - +PSI.get_parameter_multiplier(::EnergyTargetParameter, ::PSY.Storage, ::AbstractStorageFormulation) = 1.0 +PSI.get_parameter_multiplier(::EnergyLimitParameter, ::PSY.Storage, ::AbstractStorageFormulation) = 1.0 #! format: on diff --git a/test/test_storage_simulation.jl b/test/test_storage_simulation.jl index b6d8cd3..c0c0337 100644 --- a/test/test_storage_simulation.jl +++ b/test/test_storage_simulation.jl @@ -110,24 +110,33 @@ end @test run!(model) == RunStatus.SUCCESSFUL end -# TODO: Move this to Hydro? -#= -function test_2_stages_with_storage_ems(in_memory) - template_uc = - get_template_hydro_st_uc(NetworkModel(CopperPlatePowerModel; use_slacks=true)) - template_ed = - get_template_hydro_st_ed(NetworkModel(CopperPlatePowerModel; use_slacks=true)) - set_device_model!(template_ed, InterruptiblePowerLoad, StaticPowerLoad) - c_sys5_hy_uc = PSB.build_system(PSITestSystems, "c_sys5_hy_ems_uc") - c_sys5_hy_ed = PSB.build_system(PSITestSystems, "c_sys5_hy_ems_ed") +@testset "Simulation with 2-Stages EnergyLimitFeedforward with GenericBattery" begin + sys_uc = build_system(PSITestSystems, "c_sys5_bat") + sys_ed = build_system(PSITestSystems, "c_sys5_bat") + + template_uc = get_template_basic_uc_storage_simulation() + template_ed = get_template_dispatch_storage_simulation() + models = SimulationModels(; decision_models=[ - DecisionModel(template_uc, c_sys5_hy_uc; name="UC", optimizer=GLPK_optimizer), - DecisionModel(template_ed, c_sys5_hy_ed; name="ED", optimizer=GLPK_optimizer), + DecisionModel( + template_uc, + sys_uc; + name="UC", + optimizer=GLPK_optimizer, + store_variable_names=true, + ), + DecisionModel( + template_ed, + sys_ed; + name="ED", + optimizer=GLPK_optimizer, + store_variable_names=true, + ), ], ) - sequence_cache = SimulationSequence(; + sequence = SimulationSequence(; models=models, feedforwards=Dict( "ED" => [ @@ -137,31 +146,37 @@ function test_2_stages_with_storage_ems(in_memory) affected_values=[ActivePowerVariable], ), EnergyLimitFeedforward(; - component_type=HydroEnergyReservoir, - source=ActivePowerVariable, - affected_values=[ActivePowerVariable], + component_type=GenericBattery, + source=ActivePowerOutVariable, + affected_values=[ActivePowerOutVariable], number_of_periods=12, ), ], ), ini_cond_chronology=InterProblemChronology(), ) + sim_cache = Simulation(; - name="cache", + name="sim", steps=2, models=models, - sequence=sequence_cache, + sequence=sequence, simulation_folder=mktempdir(; cleanup=true), ) + build_out = build!(sim_cache) @test build_out == PSI.BuildStatus.BUILT - execute_out = execute!(sim_cache; in_memory=in_memory) + + execute_out = execute!(sim_cache) @test execute_out == PSI.RunStatus.SUCCESSFUL -end -@testset "Simulation with 2-Stages with Storage EMS" begin - for in_memory in (true, false) - test_2_stages_with_storage_ems(in_memory) - end + # Test UC Vars are equal to ED params + res = SimulationResults(sim_cache) + res_ed = res.decision_problem_results["ED"] + param_ed = read_realized_parameter(res_ed, "EnergyLimitParameter__GenericBattery") + + res_uc = res.decision_problem_results["UC"] + p_out_bat = read_realized_variable(res_uc, "ActivePowerOutVariable__GenericBattery") + + @test isapprox(param_ed[!, 2], p_out_bat[!, 2] / 100.0; atol=1e-4) end -=# diff --git a/test/test_utils/operations_problems_templates.jl b/test/test_utils/operations_problems_templates.jl index c9f4264..4bb198a 100644 --- a/test/test_utils/operations_problems_templates.jl +++ b/test/test_utils/operations_problems_templates.jl @@ -73,3 +73,39 @@ function get_template_dispatch_with_network(network=StandardPTDFModel) set_device_model!(template, TwoTerminalHVDCLine, HVDCTwoTerminalLossless) return template end + +function get_template_basic_uc_storage_simulation() + template = ProblemTemplate(CopperPlatePowerModel) + set_device_model!(template, ThermalStandard, ThermalBasicUnitCommitment) + set_device_model!(template, RenewableDispatch, RenewableFullDispatch) + set_device_model!(template, PowerLoad, StaticPowerLoad) + device_model = DeviceModel( + GenericBattery, + StorageDispatchWithReserves; + attributes=Dict{String, Any}( + "reservation" => true, + "cycling_limits" => false, + "energy_target" => false, + ), + ) + set_device_model!(template, device_model) + return template +end + +function get_template_dispatch_storage_simulation() + template = ProblemTemplate(CopperPlatePowerModel) + set_device_model!(template, ThermalStandard, ThermalBasicDispatch) + set_device_model!(template, RenewableDispatch, RenewableFullDispatch) + set_device_model!(template, PowerLoad, StaticPowerLoad) + device_model = DeviceModel( + GenericBattery, + StorageDispatchWithReserves; + attributes=Dict{String, Any}( + "reservation" => true, + "cycling_limits" => false, + "energy_target" => false, + ), + ) + set_device_model!(template, device_model) + return template +end