restore ff code

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"

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")
 

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
-=#

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

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