diff --git a/Project.toml b/Project.toml index 01523f917..53a901c19 100644 --- a/Project.toml +++ b/Project.toml @@ -45,7 +45,7 @@ LinearAlgebra = "1" Logging = "1" MathOptInterface = "1" PowerModels = "^0.21" -PowerNetworkMatrices = "^0.11" +PowerNetworkMatrices = "^0.11.1" PowerSystems = "4" PrettyTables = "2" ProgressMeter = "^1.5" diff --git a/src/PowerSimulations.jl b/src/PowerSimulations.jl index e3a3eb9fb..f4d882d98 100644 --- a/src/PowerSimulations.jl +++ b/src/PowerSimulations.jl @@ -513,7 +513,6 @@ include("feedforward/feedforward_arguments.jl") include("feedforward/feedforward_constraints.jl") include("parameters/add_parameters.jl") -include("parameters/update_parameters.jl") include("simulation/optimization_output_cache.jl") include("simulation/optimization_output_caches.jl") @@ -534,6 +533,10 @@ include("simulation/simulation_internal.jl") include("simulation/simulation.jl") include("simulation/simulation_results_export.jl") include("simulation/simulation_results.jl") +include("operation/operation_model_simulation_interface.jl") +include("parameters/update_container_parameter_values.jl") +include("parameters/update_cost_parameters.jl") +include("parameters/update_parameters.jl") include("devices_models/devices/common/objective_function/common.jl") include("devices_models/devices/common/objective_function/linear_curve.jl") diff --git a/src/core/settings.jl b/src/core/settings.jl index d9dd095cb..063ad1cce 100644 --- a/src/core/settings.jl +++ b/src/core/settings.jl @@ -16,6 +16,7 @@ struct Settings export_pwl_vars::Bool allow_fails::Bool rebuild_model::Bool + export_optimization_model::Bool store_variable_names::Bool check_numerical_bounds::Bool ext::Dict{String, Any} @@ -41,6 +42,7 @@ function Settings( allow_fails::Bool = false, check_numerical_bounds = true, rebuild_model = false, + export_optimization_model = false, store_variable_names = false, ext = Dict{String, Any}(), ) @@ -77,6 +79,7 @@ function Settings( export_pwl_vars, allow_fails, rebuild_model, + export_optimization_model, store_variable_names, check_numerical_bounds, ext, @@ -151,6 +154,7 @@ get_detailed_optimizer_stats(settings::Settings) = settings.detailed_optimizer_s get_direct_mode_optimizer(settings::Settings) = settings.direct_mode_optimizer get_store_variable_names(settings::Settings) = settings.store_variable_names get_rebuild_model(settings::Settings) = settings.rebuild_model +get_export_optimization_model(settings::Settings) = settings.export_optimization_model use_time_series_cache(settings::Settings) = settings.time_series_cache_size > 0 function set_horizon!(settings::Settings, horizon::Dates.TimePeriod) diff --git a/src/operation/decision_model.jl b/src/operation/decision_model.jl index b0718de64..785742ac7 100644 --- a/src/operation/decision_model.jl +++ b/src/operation/decision_model.jl @@ -114,6 +114,7 @@ function DecisionModel{M}( direct_mode_optimizer = false, store_variable_names = false, rebuild_model = false, + export_optimization_model = false, check_numerical_bounds = true, initial_time = UNSET_INI_TIME, time_series_cache_size::Int = IS.TIME_SERIES_CACHE_SIZE_BYTES, @@ -139,6 +140,7 @@ function DecisionModel{M}( check_numerical_bounds = check_numerical_bounds, store_variable_names = store_variable_names, rebuild_model = rebuild_model, + export_optimization_model = export_optimization_model, ) return DecisionModel{M}(template, sys, settings, jump_model; name = name) end @@ -444,7 +446,7 @@ keyword arguments to that function. - `console_level = Logging.Error`: - `file_level = Logging.Info`: - `disable_timer_outputs = false` : Enable/Disable timing outputs - - `serialize::Bool = true`: If true, serialize the model to a file to allow re-execution later. + - `export_optimization_problem::Bool = true`: If true, serialize the model to a file to allow re-execution later. # Examples @@ -459,7 +461,7 @@ function solve!( console_level = Logging.Error, file_level = Logging.Info, disable_timer_outputs = false, - serialize = true, + export_optimization_problem = true, kwargs..., ) build_if_not_already_built!( @@ -500,7 +502,7 @@ function solve!( current_time, ) end - if serialize + if export_optimization_problem TimerOutputs.@timeit RUN_OPERATION_MODEL_TIMER "Serialize" begin serialize_problem(model; optimizer = optimizer) serialize_optimization_model(model) @@ -558,19 +560,3 @@ function solve!( end return get_run_status(model) end - -function update_parameters!( - model::DecisionModel, - decision_states::DatasetContainer{InMemoryDataset}, -) - cost_function_unsynch(get_optimization_container(model)) - for key in keys(get_parameters(model)) - update_parameter_values!(model, key, decision_states) - end - if !is_synchronized(model) - update_objective_function!(get_optimization_container(model)) - obj_func = get_objective_expression(get_optimization_container(model)) - set_synchronized_status!(obj_func, true) - end - return -end diff --git a/src/operation/emulation_model.jl b/src/operation/emulation_model.jl index c635cf35b..cd42d60ec 100644 --- a/src/operation/emulation_model.jl +++ b/src/operation/emulation_model.jl @@ -454,6 +454,30 @@ function update_model!( return end +""" +Update parameter function an OperationModel +""" +function update_parameter_values!( + model::EmulationModel, + key::ParameterKey{T, U}, + input::DatasetContainer{InMemoryDataset}, +) where {T <: ParameterType, U <: PSY.Component} + # Enable again for detailed debugging + # TimerOutputs.@timeit RUN_SIMULATION_TIMER "$T $U Parameter Update" begin + optimization_container = get_optimization_container(model) + update_container_parameter_values!(optimization_container, model, key, input) + parameter_attributes = get_parameter_attributes(optimization_container, key) + IS.@record :execution ParameterUpdateEvent( + T, + U, + parameter_attributes, + get_current_timestamp(model), + get_name(model), + ) + #end + return +end + function update_model!(model::EmulationModel) update_model!(model, get_store(model), InterProblemChronology()) return @@ -507,7 +531,7 @@ keyword arguments to that function. - `export_problem_results::Bool`: If true, export OptimizationProblemResults DataFrames to CSV files. - `output_dir::String`: Required if the model is not already built, otherwise ignored - `enable_progress_bar::Bool`: Enables/Disable progress bar printing - - `serialize::Bool`: If true, serialize the model to a file to allow re-execution later. + - `export_optimization_model::Bool`: If true, serialize the model to a file to allow re-execution later. # Examples @@ -522,7 +546,7 @@ function run!( console_level = Logging.Error, file_level = Logging.Info, disable_timer_outputs = false, - serialize = true, + export_optimization_model = true, kwargs..., ) build_if_not_already_built!( @@ -555,7 +579,7 @@ function run!( run_impl!(model; kwargs...) set_run_status!(model, RunStatus.SUCCESSFULLY_FINALIZED) end - if serialize + if export_optimization_model TimerOutputs.@timeit RUN_OPERATION_MODEL_TIMER "Serialize" begin optimizer = get(kwargs, :optimizer, nothing) serialize_problem(model; optimizer = optimizer) diff --git a/src/operation/operation_model_interface.jl b/src/operation/operation_model_interface.jl index 1eb0e8ecc..d2157ee96 100644 --- a/src/operation/operation_model_interface.jl +++ b/src/operation/operation_model_interface.jl @@ -104,13 +104,22 @@ end function solve_impl!(model::OperationModel) container = get_optimization_container(model) + model_name = get_name(model) + ts = get_current_timestamp(model) + output_dir = get_output_dir(model) + + if get_export_optimization_model(get_settings(model)) + model_output_dir = joinpath(output_dir, "optimization_model_exports") + mkpath(model_output_dir) + tss = replace("$(ts)", ":" => "_") + model_export_path = joinpath(model_output_dir, "exported_$(model_name)_$(tss).json") + serialize_optimization_model(container, model_export_path) + end + status = solve_impl!(container, get_system(model)) set_run_status!(model, status) if status != RunStatus.SUCCESSFULLY_FINALIZED settings = get_settings(model) - model_name = get_name(model) - ts = get_current_timestamp(model) - output_dir = get_output_dir(model) infeasible_opt_path = joinpath(output_dir, "infeasible_$(model_name).json") @error("Serializing Infeasible Problem at $(infeasible_opt_path)") serialize_optimization_model(container, infeasible_opt_path) @@ -432,16 +441,6 @@ function serialize_optimization_model(model::OperationModel) return end -function update_model!(model::OperationModel, source, ini_cond_chronology) - TimerOutputs.@timeit RUN_SIMULATION_TIMER "Parameter Updates" begin - update_parameters!(model, get_decision_states(source)) - end - TimerOutputs.@timeit RUN_SIMULATION_TIMER "Ini Cond Updates" begin - update_initial_conditions!(model, source, ini_cond_chronology) - end - return -end - function instantiate_network_model(model::OperationModel) template = get_template(model) network_model = get_network_model(template) diff --git a/src/operation/operation_model_simulation_interface.jl b/src/operation/operation_model_simulation_interface.jl new file mode 100644 index 000000000..3ec50b2c9 --- /dev/null +++ b/src/operation/operation_model_simulation_interface.jl @@ -0,0 +1,31 @@ +function update_model!(model::OperationModel, source::SimulationState, ini_cond_chronology) + TimerOutputs.@timeit RUN_SIMULATION_TIMER "Parameter Updates" begin + update_parameters!(model, source) + end + TimerOutputs.@timeit RUN_SIMULATION_TIMER "Ini Cond Updates" begin + update_initial_conditions!(model, source, ini_cond_chronology) + end + return +end + +function update_parameters!(model::EmulationModel, state::SimulationState) + data = get_decision_states(state) + update_parameters!(model, data) + return +end + +function update_parameters!( + model::DecisionModel, + simulation_state::SimulationState, +) + cost_function_unsynch(get_optimization_container(model)) + for key in keys(get_parameters(model)) + update_parameter_values!(model, key, simulation_state) + end + if !is_synchronized(model) + update_objective_function!(get_optimization_container(model)) + obj_func = get_objective_expression(get_optimization_container(model)) + set_synchronized_status!(obj_func, true) + end + return +end diff --git a/src/parameters/update_container_parameter_values.jl b/src/parameters/update_container_parameter_values.jl new file mode 100644 index 000000000..5784c0964 --- /dev/null +++ b/src/parameters/update_container_parameter_values.jl @@ -0,0 +1,455 @@ +function _update_parameter_values!( + ::AbstractArray{T}, + ::NoAttributes, + args..., +) where {T <: Union{Float64, JuMP.VariableRef}} end + +######################## Methods to update Parameters from Time Series ##################### +function _set_param_value!(param::JuMPVariableMatrix, value::Float64, name::String, t::Int) + fix_parameter_value(param[name, t], value) + return +end + +function _set_param_value!( + param::DenseAxisArray{Vector{NTuple{2, Float64}}}, + value::Vector{NTuple{2, Float64}}, + name::String, + t::Int, +) + param[name, t] = value + return +end + +function _set_param_value!(param::JuMPFloatArray, value::Float64, name::String, t::Int) + param[name, t] = value + return +end + +function _update_parameter_values!( + parameter_array::AbstractArray{T}, + attributes::TimeSeriesAttributes{U}, + ::Type{V}, + model::DecisionModel, + ::DatasetContainer{InMemoryDataset}, +) where { + T <: Union{JuMP.VariableRef, Float64}, + U <: PSY.AbstractDeterministic, + V <: PSY.Component, +} + initial_forecast_time = get_current_time(model) # Function not well defined for DecisionModels + horizon = get_time_steps(get_optimization_container(model))[end] + ts_name = get_time_series_name(attributes) + multiplier_id = get_time_series_multiplier_id(attributes) + subsystem = get_subsystem(attributes) + template = get_template(model) + if isempty(subsystem) + device_model = get_model(template, V) + else + device_model = get_model(template, V, subsystem) + end + components = get_available_components(device_model, get_system(model)) + ts_uuids = Set{String}() + for component in components + ts_uuid = string(IS.get_time_series_uuid(U, component, ts_name)) + if !(ts_uuid in ts_uuids) + ts_vector = get_time_series_values!( + U, + model, + component, + ts_name, + multiplier_id, + initial_forecast_time, + horizon, + ) + for (t, value) in enumerate(ts_vector) + if !isfinite(value) + error("The value for the time series $(ts_name) is not finite. \ + Check that the data in the time series is valid.") + end + _set_param_value!(parameter_array, value, ts_uuid, t) + end + push!(ts_uuids, ts_uuid) + end + end +end + +function _update_parameter_values!( + parameter_array::AbstractArray{T}, + attributes::TimeSeriesAttributes{U}, + service::V, + model::DecisionModel, + ::DatasetContainer{InMemoryDataset}, +) where { + T <: Union{JuMP.VariableRef, Float64}, + U <: PSY.AbstractDeterministic, + V <: PSY.Service, +} + initial_forecast_time = get_current_time(model) # Function not well defined for DecisionModels + horizon = get_time_steps(get_optimization_container(model))[end] + ts_name = get_time_series_name(attributes) + ts_uuid = string(IS.get_time_series_uuid(U, service, ts_name)) + ts_vector = get_time_series_values!( + U, + model, + service, + get_time_series_name(attributes), + get_time_series_multiplier_id(attributes), + initial_forecast_time, + horizon, + ) + for (t, value) in enumerate(ts_vector) + if !isfinite(value) + error("The value for the time series $(ts_name) is not finite. \ + Check that the data in the time series is valid.") + end + _set_param_value!(parameter_array, value, ts_uuid, t) + end +end + +function _update_parameter_values!( + parameter_array::AbstractArray{T}, + attributes::TimeSeriesAttributes{U}, + ::Type{V}, + model::EmulationModel, + ::DatasetContainer{InMemoryDataset}, +) where {T <: Union{JuMP.VariableRef, Float64}, U <: PSY.SingleTimeSeries, V <: PSY.Device} + initial_forecast_time = get_current_time(model) + template = get_template(model) + device_model = get_model(template, V) + components = get_available_components(device_model, get_system(model)) + ts_name = get_time_series_name(attributes) + ts_uuids = Set{String}() + for component in components + ts_uuid = string(IS.get_time_series_uuid(U, component, ts_name)) + if !(ts_uuid in ts_uuids) + # Note: This interface reads one single value per component at a time. + value = get_time_series_values!( + U, + model, + component, + get_time_series_name(attributes), + get_time_series_multiplier_id(attributes), + initial_forecast_time, + )[1] + if !isfinite(value) + error("The value for the time series $(ts_name) is not finite. \ + Check that the data in the time series is valid.") + end + _set_param_value!(parameter_array, value, ts_uuid, 1) + push!(ts_uuids, ts_uuid) + end + end + return +end + +function _update_parameter_values!( + parameter_array::AbstractArray{T}, + attributes::VariableValueAttributes, + ::Type{<:PSY.Device}, + model::DecisionModel, + state::DatasetContainer{InMemoryDataset}, +) where {T <: Union{JuMP.VariableRef, Float64}} + current_time = get_current_time(model) + state_values = get_dataset_values(state, get_attribute_key(attributes)) + component_names, time = axes(parameter_array) + model_resolution = get_resolution(model) + state_data = get_dataset(state, get_attribute_key(attributes)) + state_timestamps = state_data.timestamps + max_state_index = get_num_rows(state_data) + if model_resolution < state_data.resolution + t_step = 1 + else + t_step = model_resolution ÷ state_data.resolution + end + state_data_index = find_timestamp_index(state_timestamps, current_time) + sim_timestamps = range(current_time; step = model_resolution, length = time[end]) + for t in time + timestamp_ix = min(max_state_index, state_data_index + t_step) + @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 + # 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(get_attribute_key(attributes))) 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 + _set_param_value!(parameter_array, state_value, name, t) + end + end + return +end + +function _update_parameter_values!( + parameter_array::AbstractArray{T}, + attributes::VariableValueAttributes, + ::PSY.Reserve, + model::DecisionModel, + state::DatasetContainer{InMemoryDataset}, +) where {T <: Union{JuMP.VariableRef, Float64}} + current_time = get_current_time(model) + state_values = get_dataset_values(state, get_attribute_key(attributes)) + component_names, time = axes(parameter_array) + model_resolution = get_resolution(model) + state_data = get_dataset(state, get_attribute_key(attributes)) + state_timestamps = state_data.timestamps + max_state_index = get_num_rows(state_data) + if model_resolution < state_data.resolution + t_step = 1 + else + t_step = model_resolution ÷ state_data.resolution + end + state_data_index = find_timestamp_index(state_timestamps, current_time) + sim_timestamps = range(current_time; step = model_resolution, length = time[end]) + for t in time + timestamp_ix = min(max_state_index, state_data_index + t_step) + @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 + # 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(get_attribute_key(attributes))) 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 + _set_param_value!(parameter_array, state_value, name, t) + end + end + return +end + +function _update_parameter_values!( + parameter_array::AbstractArray{T}, + attributes::VariableValueAttributes{VariableKey{OnVariable, U}}, + ::Type{U}, + model::DecisionModel, + state::DatasetContainer{InMemoryDataset}, +) where {T <: Union{JuMP.VariableRef, Float64}, U <: PSY.Device} + current_time = get_current_time(model) + state_values = get_dataset_values(state, get_attribute_key(attributes)) + component_names, time = axes(parameter_array) + model_resolution = get_resolution(model) + state_data = get_dataset(state, get_attribute_key(attributes)) + state_timestamps = state_data.timestamps + max_state_index = get_num_rows(state_data) + if model_resolution < state_data.resolution + t_step = 1 + else + t_step = model_resolution ÷ state_data.resolution + end + state_data_index = find_timestamp_index(state_timestamps, current_time) + + sim_timestamps = range(current_time; step = model_resolution, length = time[end]) + for t in time + timestamp_ix = min(max_state_index, state_data_index + t_step) + @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 + # Pass indices in this way since JuMP DenseAxisArray don't support view() + value = round(state_values[name, state_data_index]) + if !isfinite(value) + error( + "The value for the system state used in $(encode_key_as_string(get_attribute_key(attributes))) is not a finite value $(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 + if 0.0 > value || value > 1.0 + error( + "The value for the system state used in $(encode_key_as_string(get_attribute_key(attributes))): $(value) is out of the [0, 1] range", + ) + end + _set_param_value!(parameter_array, value, name, t) + end + end + return +end + +function _update_parameter_values!( + parameter_array::AbstractArray{T}, + attributes::VariableValueAttributes, + ::Type{<:PSY.Component}, + model::EmulationModel, + state::DatasetContainer{InMemoryDataset}, +) where {T <: Union{JuMP.VariableRef, Float64}} + current_time = get_current_time(model) + state_values = get_dataset_values(state, get_attribute_key(attributes)) + component_names, _ = axes(parameter_array) + state_data = get_dataset(state, get_attribute_key(attributes)) + state_timestamps = state_data.timestamps + state_data_index = find_timestamp_index(state_timestamps, current_time) + for name in component_names + # Pass indices in this way since JuMP DenseAxisArray don't support view() + _set_param_value!(parameter_array, state_values[name, state_data_index], name, 1) + end + return +end + +function _update_parameter_values!( + parameter_array::AbstractArray{T}, + attributes::VariableValueAttributes{VariableKey{OnVariable, U}}, + ::Type{<:PSY.Component}, + model::EmulationModel, + state::DatasetContainer{InMemoryDataset}, +) where {T <: Union{JuMP.VariableRef, Float64}, U <: PSY.Component} + current_time = get_current_time(model) + state_values = get_dataset_values(state, get_attribute_key(attributes)) + component_names, _ = axes(parameter_array) + state_data = get_dataset(state, get_attribute_key(attributes)) + state_timestamps = state_data.timestamps + state_data_index = find_timestamp_index(state_timestamps, current_time) + for name in component_names + # Pass indices in this way since JuMP DenseAxisArray don't support view() + value = round(state_values[name, state_data_index]) + if !isfinite(value) + error( + "The value for the system state used in $(encode_key_as_string(get_attribute_key(attributes))) is not a finite value $(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 + if 0.0 > value || value > 1.0 + error( + "The value for the system state used in $(encode_key_as_string(get_attribute_key(attributes))): $(value) is out of the [0, 1] range", + ) + end + _set_param_value!(parameter_array, value, name, 1) + end + return +end + +function _update_parameter_values!( + ::AbstractArray{T}, + ::VariableValueAttributes, + ::Type{<:PSY.Component}, + ::EmulationModel, + ::EmulationModelStore, +) where {T <: Union{JuMP.VariableRef, Float64}} + error("The emulation model has parameters that can't be updated from its results") + return +end + +""" +Update parameter function an OperationModel +""" +function update_container_parameter_values!( + optimization_container::OptimizationContainer, + model::OperationModel, + key::ParameterKey{T, U}, + input::DatasetContainer{InMemoryDataset}, +) where {T <: ParameterType, U <: PSY.Component} + # Enable again for detailed debugging + # TimerOutputs.@timeit RUN_SIMULATION_TIMER "$T $U Parameter Update" begin + # 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 = get_parameter_array(optimization_container, key) + parameter_attributes = get_parameter_attributes(optimization_container, key) + _update_parameter_values!(parameter_array, parameter_attributes, U, model, input) + return +end + +function update_container_parameter_values!( + optimization_container::OptimizationContainer, + model::OperationModel, + key::ParameterKey{T, U}, + input::DatasetContainer{InMemoryDataset}, +) where {T <: ObjectiveFunctionParameter, U <: PSY.Component} + # 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 = get_parameter_array(optimization_container, key) + # Multiplier is only needed for the objective function since `_update_parameter_values!` also updates the objective function + parameter_multiplier = get_parameter_multiplier_array(optimization_container, key) + parameter_attributes = get_parameter_attributes(optimization_container, key) + _update_parameter_values!( + parameter_array, + parameter_multiplier, + parameter_attributes, + U, + model, + input, + ) + return +end + +function update_container_parameter_values!( + optimization_container::OptimizationContainer, + model::OperationModel, + key::ParameterKey{T, U}, + input::DatasetContainer{InMemoryDataset}, +) where {T <: ObjectiveFunctionParameter, U <: PSY.Service} + # 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 = get_parameter_array(optimization_container, key) + # Multiplier is only needed for the objective function since `_update_parameter_values!` also updates the objective function + parameter_multiplier = get_parameter_multiplier_array(optimization_container, key) + parameter_attributes = get_parameter_attributes(optimization_container, key) + _update_parameter_values!( + parameter_array, + parameter_multiplier, + parameter_attributes, + U, + model, + input, + ) + return +end + +function update_container_parameter_values!( + optimization_container::OptimizationContainer, + model::OperationModel, + key::ParameterKey{FixValueParameter, U}, + input::DatasetContainer{InMemoryDataset}, +) where {U <: PSY.Component} + # 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 = get_parameter_array(optimization_container, key) + parameter_attributes = get_parameter_attributes(optimization_container, key) + _update_parameter_values!(parameter_array, parameter_attributes, U, model, input) + _fix_parameter_value!(optimization_container, parameter_array, parameter_attributes) + return +end + +function update_container_parameter_values!( + optimization_container::OptimizationContainer, + model::OperationModel, + key::ParameterKey{FixValueParameter, U}, + input::DatasetContainer{InMemoryDataset}, +) where {U <: PSY.Service} + # 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 = get_parameter_array(optimization_container, key) + parameter_attributes = get_parameter_attributes(optimization_container, key) + service = PSY.get_component(U, get_system(model), key.meta) + @assert service !== nothing + _update_parameter_values!(parameter_array, parameter_attributes, U, model, input) + _fix_parameter_value!(optimization_container, parameter_array, parameter_attributes) + return +end + +function update_container_parameter_values!( + optimization_container::OptimizationContainer, + model::OperationModel, + key::ParameterKey{T, U}, + input::DatasetContainer{InMemoryDataset}, +) where {T <: ParameterType, U <: PSY.Service} + # 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 = get_parameter_array(optimization_container, key) + parameter_attributes = get_parameter_attributes(optimization_container, key) + service = PSY.get_component(U, get_system(model), key.meta) + @assert service !== nothing + _update_parameter_values!(parameter_array, parameter_attributes, service, model, input) + return +end diff --git a/src/parameters/update_cost_parameters.jl b/src/parameters/update_cost_parameters.jl new file mode 100644 index 000000000..600eca675 --- /dev/null +++ b/src/parameters/update_cost_parameters.jl @@ -0,0 +1,125 @@ +function _update_parameter_values!( + parameter_array::DenseAxisArray, + parameter_multiplier::JuMPFloatArray, + attributes::CostFunctionAttributes, + ::Type{V}, + model::DecisionModel, + ::DatasetContainer{InMemoryDataset}, +) where {V <: PSY.Component} + initial_forecast_time = get_current_time(model) # Function not well defined for DecisionModels + time_steps = get_time_steps(get_optimization_container(model)) + horizon = time_steps[end] + container = get_optimization_container(model) + @assert !is_synchronized(container) + template = get_template(model) + device_model = get_model(template, V) + components = get_available_components(device_model, get_system(model)) + + for component in components + if _has_variable_cost_parameter(component) + name = PSY.get_name(component) + ts_vector = PSY.get_variable_cost( + component, + PSY.get_operation_cost(component); + start_time = initial_forecast_time, + len = horizon, + ) + variable_cost_forecast_values = TimeSeries.values(ts_vector) + for (t, value) in enumerate(variable_cost_forecast_values) + if attributes.uses_compact_power + # TODO implement this + value, _ = _convert_variable_cost(value) + end + # TODO removed an apparently unused block of code here? + _set_param_value!(parameter_array, value, name, t) + update_variable_cost!( + container, + parameter_array, + parameter_multiplier, + attributes, + component, + t, + ) + end + end + end + return +end + +_has_variable_cost_parameter(component::PSY.Component) = + _has_variable_cost_parameter(PSY.get_operation_cost(component)) +_has_variable_cost_parameter(::PSY.MarketBidCost) = true +_has_variable_cost_parameter(::T) where {T <: PSY.OperationalCost} = false + +function _update_pwl_cost_expression( + container::OptimizationContainer, + ::Type{T}, + component_name::String, + time_period::Int, + cost_data::PSY.PiecewiseLinearData, +) where {T <: PSY.Component} + pwl_var_container = get_variable(container, PieceWiseLinearCostVariable(), T) + resolution = get_resolution(container) + dt = Dates.value(resolution) / MILLISECONDS_IN_HOUR + gen_cost = JuMP.AffExpr(0.0) + slopes = PSY.get_slopes(cost_data) + upb = get_breakpoint_upper_bounds(cost_data) + for i in 1:length(cost_data) + JuMP.add_to_expression!( + gen_cost, + slopes[i] * upb[i] * dt * pwl_var_container[(component_name, i, time_period)], + ) + end + return gen_cost +end + +function update_variable_cost!( + container::OptimizationContainer, + parameter_array::JuMPFloatArray, + parameter_multiplier::JuMPFloatArray, + attributes::CostFunctionAttributes{Float64}, + component::T, + time_period::Int, +) where {T <: PSY.Component} + resolution = get_resolution(container) + dt = Dates.value(resolution) / MILLISECONDS_IN_HOUR + base_power = get_base_power(container) + component_name = PSY.get_name(component) + cost_data = parameter_array[component_name, time_period] # TODO is this a new-style cost? + if iszero(cost_data) + return + end + mult_ = parameter_multiplier[component_name, time_period] + variable = get_variable(container, get_variable_type(attributes)(), T) + gen_cost = variable[component_name, time_period] * mult_ * cost_data * base_power * dt + add_to_objective_variant_expression!(container, gen_cost) + set_expression!(container, ProductionCostExpression, gen_cost, component, time_period) + return +end + +function update_variable_cost!( + container::OptimizationContainer, + parameter_array::DenseAxisArray{Vector{NTuple{2, Float64}}}, + parameter_multiplier::JuMPFloatArray, + ::CostFunctionAttributes{Vector{NTuple{2, Float64}}}, + component::T, + time_period::Int, +) where {T <: PSY.Component} + component_name = PSY.get_name(component) + cost_data = parameter_array[component_name, time_period] + if all(iszero.(last.(cost_data))) + return + end + mult_ = parameter_multiplier[component_name, time_period] + gen_cost = + _update_pwl_cost_expression( + container, + T, + component_name, + time_period, + PSY.PiecewiseLinearData(cost_data), + ) + add_to_objective_variant_expression!(container, mult_ * gen_cost) + set_expression!(container, ProductionCostExpression, gen_cost, component, time_period) + return +end diff --git a/src/parameters/update_parameters.jl b/src/parameters/update_parameters.jl index 958cc8e27..53e3f4d35 100644 --- a/src/parameters/update_parameters.jl +++ b/src/parameters/update_parameters.jl @@ -1,470 +1,15 @@ -function _update_parameter_values!( - ::AbstractArray{T}, - ::NoAttributes, - args..., -) where {T <: Union{Float64, JuMP.VariableRef}} end - -######################## Methods to update Parameters from Time Series ##################### -function _set_param_value!(param::JuMPVariableMatrix, value::Float64, name::String, t::Int) - fix_parameter_value(param[name, t], value) - return -end - -function _set_param_value!( - param::DenseAxisArray{Vector{NTuple{2, Float64}}}, - value::Vector{NTuple{2, Float64}}, - name::String, - t::Int, -) - param[name, t] = value - return -end - -function _set_param_value!(param::JuMPFloatArray, value::Float64, name::String, t::Int) - param[name, t] = value - return -end - -function _update_parameter_values!( - parameter_array::AbstractArray{T}, - attributes::TimeSeriesAttributes{U}, - ::Type{V}, - model::DecisionModel, - ::DatasetContainer{InMemoryDataset}, -) where { - T <: Union{JuMP.VariableRef, Float64}, - U <: PSY.AbstractDeterministic, - V <: PSY.Component, -} - initial_forecast_time = get_current_time(model) # Function not well defined for DecisionModels - horizon = get_time_steps(get_optimization_container(model))[end] - ts_name = get_time_series_name(attributes) - multiplier_id = get_time_series_multiplier_id(attributes) - subsystem = get_subsystem(attributes) - template = get_template(model) - if isempty(subsystem) - device_model = get_model(template, V) - else - device_model = get_model(template, V, subsystem) - end - components = get_available_components(device_model, get_system(model)) - ts_uuids = Set{String}() - for component in components - ts_uuid = string(IS.get_time_series_uuid(U, component, ts_name)) - if !(ts_uuid in ts_uuids) - ts_vector = get_time_series_values!( - U, - model, - component, - ts_name, - multiplier_id, - initial_forecast_time, - horizon, - ) - for (t, value) in enumerate(ts_vector) - if !isfinite(value) - error("The value for the time series $(ts_name) is not finite. \ - Check that the data in the time series is valid.") - end - _set_param_value!(parameter_array, value, ts_uuid, t) - end - push!(ts_uuids, ts_uuid) - end - end -end - -function _update_parameter_values!( - parameter_array::AbstractArray{T}, - attributes::TimeSeriesAttributes{U}, - service::V, - model::DecisionModel, - ::DatasetContainer{InMemoryDataset}, -) where { - T <: Union{JuMP.VariableRef, Float64}, - U <: PSY.AbstractDeterministic, - V <: PSY.Service, -} - initial_forecast_time = get_current_time(model) # Function not well defined for DecisionModels - horizon = get_time_steps(get_optimization_container(model))[end] - ts_name = get_time_series_name(attributes) - ts_uuid = string(IS.get_time_series_uuid(U, service, ts_name)) - ts_vector = get_time_series_values!( - U, - model, - service, - get_time_series_name(attributes), - get_time_series_multiplier_id(attributes), - initial_forecast_time, - horizon, - ) - for (t, value) in enumerate(ts_vector) - if !isfinite(value) - error("The value for the time series $(ts_name) is not finite. \ - Check that the data in the time series is valid.") - end - _set_param_value!(parameter_array, value, ts_uuid, t) - end -end - -function _update_parameter_values!( - parameter_array::AbstractArray{T}, - attributes::TimeSeriesAttributes{U}, - ::Type{V}, - model::EmulationModel, - ::DatasetContainer{InMemoryDataset}, -) where {T <: Union{JuMP.VariableRef, Float64}, U <: PSY.SingleTimeSeries, V <: PSY.Device} - initial_forecast_time = get_current_time(model) - template = get_template(model) - device_model = get_model(template, V) - components = get_available_components(device_model, get_system(model)) - ts_name = get_time_series_name(attributes) - ts_uuids = Set{String}() - for component in components - ts_uuid = string(IS.get_time_series_uuid(U, component, ts_name)) - if !(ts_uuid in ts_uuids) - # Note: This interface reads one single value per component at a time. - value = get_time_series_values!( - U, - model, - component, - get_time_series_name(attributes), - get_time_series_multiplier_id(attributes), - initial_forecast_time, - )[1] - if !isfinite(value) - error("The value for the time series $(ts_name) is not finite. \ - Check that the data in the time series is valid.") - end - _set_param_value!(parameter_array, value, ts_uuid, 1) - push!(ts_uuids, ts_uuid) - end - end - return -end - -function _update_parameter_values!( - parameter_array::AbstractArray{T}, - attributes::VariableValueAttributes, - ::Type{<:PSY.Device}, - model::DecisionModel, - state::DatasetContainer{InMemoryDataset}, -) where {T <: Union{JuMP.VariableRef, Float64}} - current_time = get_current_time(model) - state_values = get_dataset_values(state, get_attribute_key(attributes)) - component_names, time = axes(parameter_array) - model_resolution = get_resolution(model) - state_data = get_dataset(state, get_attribute_key(attributes)) - state_timestamps = state_data.timestamps - max_state_index = get_num_rows(state_data) - if model_resolution < state_data.resolution - t_step = 1 - else - t_step = model_resolution ÷ state_data.resolution - end - state_data_index = find_timestamp_index(state_timestamps, current_time) - sim_timestamps = range(current_time; step = model_resolution, length = time[end]) - for t in time - timestamp_ix = min(max_state_index, state_data_index + t_step) - @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 - # 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(get_attribute_key(attributes))) 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 - _set_param_value!(parameter_array, state_value, name, t) - end - end - return -end - -function _update_parameter_values!( - parameter_array::AbstractArray{T}, - attributes::VariableValueAttributes, - ::PSY.Reserve, - model::DecisionModel, - state::DatasetContainer{InMemoryDataset}, -) where {T <: Union{JuMP.VariableRef, Float64}} - current_time = get_current_time(model) - state_values = get_dataset_values(state, get_attribute_key(attributes)) - component_names, time = axes(parameter_array) - model_resolution = get_resolution(model) - state_data = get_dataset(state, get_attribute_key(attributes)) - state_timestamps = state_data.timestamps - max_state_index = get_num_rows(state_data) - if model_resolution < state_data.resolution - t_step = 1 - else - t_step = model_resolution ÷ state_data.resolution - end - state_data_index = find_timestamp_index(state_timestamps, current_time) - sim_timestamps = range(current_time; step = model_resolution, length = time[end]) - for t in time - timestamp_ix = min(max_state_index, state_data_index + t_step) - @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 - # 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(get_attribute_key(attributes))) 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 - _set_param_value!(parameter_array, state_value, name, t) - end - end - return -end - -function _update_parameter_values!( - parameter_array::AbstractArray{T}, - attributes::VariableValueAttributes{VariableKey{OnVariable, U}}, - ::Type{U}, - model::DecisionModel, - state::DatasetContainer{InMemoryDataset}, -) where {T <: Union{JuMP.VariableRef, Float64}, U <: PSY.Device} - current_time = get_current_time(model) - state_values = get_dataset_values(state, get_attribute_key(attributes)) - component_names, time = axes(parameter_array) - model_resolution = get_resolution(model) - state_data = get_dataset(state, get_attribute_key(attributes)) - state_timestamps = state_data.timestamps - max_state_index = get_num_rows(state_data) - if model_resolution < state_data.resolution - t_step = 1 - else - t_step = model_resolution ÷ state_data.resolution - end - state_data_index = find_timestamp_index(state_timestamps, current_time) - - sim_timestamps = range(current_time; step = model_resolution, length = time[end]) - for t in time - timestamp_ix = min(max_state_index, state_data_index + t_step) - @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 - # Pass indices in this way since JuMP DenseAxisArray don't support view() - value = round(state_values[name, state_data_index]) - if !isfinite(value) - error( - "The value for the system state used in $(encode_key_as_string(get_attribute_key(attributes))) is not a finite value $(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 - if 0.0 > value || value > 1.0 - error( - "The value for the system state used in $(encode_key_as_string(get_attribute_key(attributes))): $(value) is out of the [0, 1] range", - ) - end - _set_param_value!(parameter_array, value, name, t) - end - end - return -end - -function _update_parameter_values!( - parameter_array::AbstractArray{T}, - attributes::VariableValueAttributes, - ::Type{<:PSY.Component}, - model::EmulationModel, - state::DatasetContainer{InMemoryDataset}, -) where {T <: Union{JuMP.VariableRef, Float64}} - current_time = get_current_time(model) - state_values = get_dataset_values(state, get_attribute_key(attributes)) - component_names, _ = axes(parameter_array) - state_data = get_dataset(state, get_attribute_key(attributes)) - state_timestamps = state_data.timestamps - state_data_index = find_timestamp_index(state_timestamps, current_time) - for name in component_names - # Pass indices in this way since JuMP DenseAxisArray don't support view() - _set_param_value!(parameter_array, state_values[name, state_data_index], name, 1) - end - return -end - -function _update_parameter_values!( - parameter_array::AbstractArray{T}, - attributes::VariableValueAttributes{VariableKey{OnVariable, U}}, - ::Type{<:PSY.Component}, - model::EmulationModel, - state::DatasetContainer{InMemoryDataset}, -) where {T <: Union{JuMP.VariableRef, Float64}, U <: PSY.Component} - current_time = get_current_time(model) - state_values = get_dataset_values(state, get_attribute_key(attributes)) - component_names, _ = axes(parameter_array) - state_data = get_dataset(state, get_attribute_key(attributes)) - state_timestamps = state_data.timestamps - state_data_index = find_timestamp_index(state_timestamps, current_time) - for name in component_names - # Pass indices in this way since JuMP DenseAxisArray don't support view() - value = round(state_values[name, state_data_index]) - @assert 0.0 <= value <= 1.0 - if !isfinite(value) - error( - "The value for the system state used in $(encode_key_as_string(get_attribute_key(attributes))) is not a finite value $(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 - _set_param_value!(parameter_array, value, name, 1) - end - return -end - -function _update_parameter_values!( - ::AbstractArray{T}, - ::VariableValueAttributes, - ::Type{<:PSY.Component}, - ::EmulationModel, - ::EmulationModelStore, -) where {T <: Union{JuMP.VariableRef, Float64}} - error("The emulation model has parameters that can't be updated from its results") - return -end - -""" -Update parameter function an OperationModel -""" -function update_container_parameter_values!( - optimization_container::OptimizationContainer, - model::OperationModel, - key::ParameterKey{T, U}, - input::DatasetContainer{InMemoryDataset}, -) where {T <: ParameterType, U <: PSY.Component} - # Enable again for detailed debugging - # TimerOutputs.@timeit RUN_SIMULATION_TIMER "$T $U Parameter Update" begin - # 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 = get_parameter_array(optimization_container, key) - parameter_attributes = get_parameter_attributes(optimization_container, key) - _update_parameter_values!(parameter_array, parameter_attributes, U, model, input) - return -end - -function update_container_parameter_values!( - optimization_container::OptimizationContainer, - model::OperationModel, - key::ParameterKey{T, U}, - input::DatasetContainer{InMemoryDataset}, -) where {T <: ObjectiveFunctionParameter, U <: PSY.Component} - # 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 = get_parameter_array(optimization_container, key) - # Multiplier is only needed for the objective function since `_update_parameter_values!` also updates the objective function - parameter_multiplier = get_parameter_multiplier_array(optimization_container, key) - parameter_attributes = get_parameter_attributes(optimization_container, key) - _update_parameter_values!( - parameter_array, - parameter_multiplier, - parameter_attributes, - U, - model, - input, - ) - return -end - -function update_container_parameter_values!( - optimization_container::OptimizationContainer, - model::OperationModel, - key::ParameterKey{T, U}, - input::DatasetContainer{InMemoryDataset}, -) where {T <: ObjectiveFunctionParameter, U <: PSY.Service} - # 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 = get_parameter_array(optimization_container, key) - # Multiplier is only needed for the objective function since `_update_parameter_values!` also updates the objective function - parameter_multiplier = get_parameter_multiplier_array(optimization_container, key) - parameter_attributes = get_parameter_attributes(optimization_container, key) - _update_parameter_values!( - parameter_array, - parameter_multiplier, - parameter_attributes, - U, - model, - input, - ) - return -end - -function update_container_parameter_values!( - optimization_container::OptimizationContainer, - model::OperationModel, - key::ParameterKey{FixValueParameter, U}, - input::DatasetContainer{InMemoryDataset}, -) where {U <: PSY.Component} - # 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 = get_parameter_array(optimization_container, key) - parameter_attributes = get_parameter_attributes(optimization_container, key) - _update_parameter_values!(parameter_array, parameter_attributes, U, model, input) - _fix_parameter_value!(optimization_container, parameter_array, parameter_attributes) - return -end - -function update_container_parameter_values!( - optimization_container::OptimizationContainer, - model::OperationModel, - key::ParameterKey{FixValueParameter, U}, - input::DatasetContainer{InMemoryDataset}, -) where {U <: PSY.Service} - # 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 = get_parameter_array(optimization_container, key) - parameter_attributes = get_parameter_attributes(optimization_container, key) - _update_parameter_values!( - parameter_array, - parameter_attributes, - FixValueParameter, - model, - input, - ) - _fix_parameter_value!(optimization_container, parameter_array, parameter_attributes) - return -end - -function update_container_parameter_values!( - optimization_container::OptimizationContainer, - model::OperationModel, - key::ParameterKey{T, U}, - input::DatasetContainer{InMemoryDataset}, -) where {T <: ParameterType, U <: PSY.Service} - # 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 = get_parameter_array(optimization_container, key) - parameter_attributes = get_parameter_attributes(optimization_container, key) - service = PSY.get_component(U, get_system(model), key.meta) - @assert service !== nothing - _update_parameter_values!(parameter_array, parameter_attributes, service, model, input) - return -end - """ Update parameter function an OperationModel """ function update_parameter_values!( model::OperationModel, key::ParameterKey{T, U}, - input::DatasetContainer{InMemoryDataset}, + simulation_state::SimulationState, ) where {T <: ParameterType, U <: PSY.Component} # Enable again for detailed debugging # TimerOutputs.@timeit RUN_SIMULATION_TIMER "$T $U Parameter Update" begin optimization_container = get_optimization_container(model) + input = get_decision_states(simulation_state) update_container_parameter_values!(optimization_container, model, key, input) parameter_attributes = get_parameter_attributes(optimization_container, key) IS.@record :execution ParameterUpdateEvent( @@ -498,7 +43,7 @@ end function update_parameter_values!( model::OperationModel, key::ParameterKey{FixValueParameter, T}, - input::DatasetContainer{InMemoryDataset}, + simulation_state::SimulationState, ) where {T <: PSY.Service} # Enable again for detailed debugging # TimerOutputs.@timeit RUN_SIMULATION_TIMER "$T $U Parameter Update" begin @@ -509,7 +54,14 @@ function update_parameter_values!( parameter_attributes = get_parameter_attributes(optimization_container, key) service = PSY.get_component(T, get_system(model), key.meta) @assert service !== nothing - _update_parameter_values!(parameter_array, parameter_attributes, service, model, input) + input = get_decision_states(simulation_state) + _update_parameter_values!( + parameter_array, + parameter_attributes, + service, + model, + input, + ) _fix_parameter_value!(optimization_container, parameter_array, parameter_attributes) IS.@record :execution ParameterUpdateEvent( FixValueParameter, @@ -521,129 +73,3 @@ function update_parameter_values!( #end return end - -function _update_parameter_values!( - parameter_array::DenseAxisArray, - parameter_multiplier::JuMPFloatArray, - attributes::CostFunctionAttributes, - ::Type{V}, - model::DecisionModel, - ::DatasetContainer{InMemoryDataset}, -) where {V <: PSY.Component} - initial_forecast_time = get_current_time(model) # Function not well defined for DecisionModels - time_steps = get_time_steps(get_optimization_container(model)) - horizon = time_steps[end] - container = get_optimization_container(model) - @assert !is_synchronized(container) - template = get_template(model) - device_model = get_model(template, V) - components = get_available_components(device_model, get_system(model)) - - for component in components - if _has_variable_cost_parameter(component) - name = PSY.get_name(component) - ts_vector = PSY.get_variable_cost( - component, - PSY.get_operation_cost(component); - start_time = initial_forecast_time, - len = horizon, - ) - variable_cost_forecast_values = TimeSeries.values(ts_vector) - for (t, value) in enumerate(variable_cost_forecast_values) - if attributes.uses_compact_power - # TODO implement this - value, _ = _convert_variable_cost(value) - end - # TODO removed an apparently unused block of code here? - _set_param_value!(parameter_array, value, name, t) - update_variable_cost!( - container, - parameter_array, - parameter_multiplier, - attributes, - component, - t, - ) - end - end - end - return -end - -_has_variable_cost_parameter(component::PSY.Component) = - _has_variable_cost_parameter(PSY.get_operation_cost(component)) -_has_variable_cost_parameter(::PSY.MarketBidCost) = true -_has_variable_cost_parameter(::T) where {T <: PSY.OperationalCost} = false - -function _update_pwl_cost_expression( - container::OptimizationContainer, - ::Type{T}, - component_name::String, - time_period::Int, - cost_data::PSY.PiecewiseLinearData, -) where {T <: PSY.Component} - pwl_var_container = get_variable(container, PieceWiseLinearCostVariable(), T) - resolution = get_resolution(container) - dt = Dates.value(resolution) / MILLISECONDS_IN_HOUR - gen_cost = JuMP.AffExpr(0.0) - slopes = PSY.get_slopes(cost_data) - upb = get_breakpoint_upper_bounds(cost_data) - for i in 1:length(cost_data) - JuMP.add_to_expression!( - gen_cost, - slopes[i] * upb[i] * dt * pwl_var_container[(component_name, i, time_period)], - ) - end - return gen_cost -end - -function update_variable_cost!( - container::OptimizationContainer, - parameter_array::JuMPFloatArray, - parameter_multiplier::JuMPFloatArray, - attributes::CostFunctionAttributes{Float64}, - component::T, - time_period::Int, -) where {T <: PSY.Component} - resolution = get_resolution(container) - dt = Dates.value(resolution) / MILLISECONDS_IN_HOUR - base_power = get_base_power(container) - component_name = PSY.get_name(component) - cost_data = parameter_array[component_name, time_period] # TODO is this a new-style cost? - if iszero(cost_data) - return - end - mult_ = parameter_multiplier[component_name, time_period] - variable = get_variable(container, get_variable_type(attributes)(), T) - gen_cost = variable[component_name, time_period] * mult_ * cost_data * base_power * dt - add_to_objective_variant_expression!(container, gen_cost) - set_expression!(container, ProductionCostExpression, gen_cost, component, time_period) - return -end - -function update_variable_cost!( - container::OptimizationContainer, - parameter_array::DenseAxisArray{Vector{NTuple{2, Float64}}}, - parameter_multiplier::JuMPFloatArray, - ::CostFunctionAttributes{Vector{NTuple{2, Float64}}}, - component::T, - time_period::Int, -) where {T <: PSY.Component} - component_name = PSY.get_name(component) - cost_data = parameter_array[component_name, time_period] - if all(iszero.(last.(cost_data))) - return - end - mult_ = parameter_multiplier[component_name, time_period] - gen_cost = - _update_pwl_cost_expression( - container, - T, - component_name, - time_period, - PSY.PiecewiseLinearData(cost_data), - ) - add_to_objective_variant_expression!(container, mult_ * gen_cost) - set_expression!(container, ProductionCostExpression, gen_cost, component, time_period) - return -end diff --git a/test/test_model_emulation.jl b/test/test_model_emulation.jl index f824c7ec4..968807331 100644 --- a/test/test_model_emulation.jl +++ b/test/test_model_emulation.jl @@ -223,7 +223,7 @@ end model; executions = 10, output_dir = mktempdir(; cleanup = true), - serialize = serialize, + export_optimization_model = serialize, ) == PSI.RunStatus.SUCCESSFULLY_FINALIZED end end diff --git a/test/test_simulation_build.jl b/test/test_simulation_build.jl index d6d4293ca..db4b9e6be 100644 --- a/test/test_simulation_build.jl +++ b/test/test_simulation_build.jl @@ -320,3 +320,80 @@ end ) @test !isempty(c) end + +@testset "Test FixValue Feedforwards" begin + template_uc = get_template_basic_uc_simulation() + set_network_model!(template_uc, NetworkModel(PTDFPowerModel; use_slacks = true)) + set_device_model!(template_uc, DeviceModel(Line, StaticBranchUnbounded)) + set_service_model!(template_uc, ServiceModel(VariableReserve{ReserveUp}, RangeReserve)) + template_ed = + get_template_nomin_ed_simulation(NetworkModel(PTDFPowerModel; use_slacks = true)) + set_device_model!(template_ed, DeviceModel(Line, StaticBranchUnbounded)) + set_service_model!(template_ed, ServiceModel(VariableReserve{ReserveUp}, RangeReserve)) + c_sys5_hy_uc = PSB.build_system(PSITestSystems, "c_sys5_uc"; add_reserves = true) + c_sys5_hy_ed = PSB.build_system(PSITestSystems, "c_sys5_ed"; add_reserves = true) + models = SimulationModels(; + decision_models = [ + DecisionModel( + template_uc, + c_sys5_hy_uc; + name = "UC", + optimizer = HiGHS_optimizer, + initialize_model = false, + ), + DecisionModel( + template_ed, + c_sys5_hy_ed; + name = "ED", + optimizer = ipopt_optimizer, + initialize_model = false, + ), + ], + ) + + sequence = SimulationSequence(; + models = models, + feedforwards = Dict( + "ED" => [ + SemiContinuousFeedforward(; + component_type = ThermalStandard, + source = OnVariable, + affected_values = [ActivePowerVariable], + ), + FixValueFeedforward(; + component_type = VariableReserve{ReserveUp}, + source = ActivePowerReserveVariable, + affected_values = [ActivePowerReserveVariable], + ), + ], + ), + ini_cond_chronology = InterProblemChronology(), + ) + + sim = Simulation(; + name = "reserve_feedforward", + steps = 2, + models = models, + sequence = sequence, + simulation_folder = mktempdir(; cleanup = true), + ) + build_out = build!(sim) + @test build_out == PSI.SimulationBuildStatus.BUILT + ed_power_model = PSI.get_simulation_model(PSI.get_models(sim), :ED) + c = PSI.get_parameter( + PSI.get_optimization_container(ed_power_model), + FixValueParameter(), + VariableReserve{ReserveUp}, + "Reserve1", + ) + @test !isempty(c.multiplier_array) + @test !isempty(c.parameter_array) + c = PSI.get_parameter( + PSI.get_optimization_container(ed_power_model), + FixValueParameter(), + VariableReserve{ReserveUp}, + "Reserve11", + ) + @test !isempty(c.multiplier_array) + @test !isempty(c.parameter_array) +end diff --git a/test/test_simulation_execute.jl b/test/test_simulation_execute.jl index 6bc458176..22b8cbb25 100644 --- a/test/test_simulation_execute.jl +++ b/test/test_simulation_execute.jl @@ -1,4 +1,5 @@ -function test_single_stage_sequential(in_memory, rebuild) +function test_single_stage_sequential(in_memory, rebuild, export_model) + tmp_dir = mktempdir(; cleanup = true) template_ed = get_template_nomin_ed_simulation() c_sys = PSB.build_system(PSITestSystems, "c_sys5_uc") models = SimulationModels([ @@ -8,6 +9,7 @@ function test_single_stage_sequential(in_memory, rebuild) name = "ED", optimizer = ipopt_optimizer, rebuild_model = rebuild, + export_optimization_model = export_model, ), ]) test_sequence = @@ -20,20 +22,29 @@ function test_single_stage_sequential(in_memory, rebuild) steps = 2, models = models, sequence = test_sequence, - simulation_folder = mktempdir(; cleanup = true), + simulation_folder = tmp_dir, ) build_out = build!(sim_single) @test build_out == PSI.SimulationBuildStatus.BUILT execute_out = execute!(sim_single; in_memory = in_memory) @test execute_out == PSI.RunStatus.SUCCESSFULLY_FINALIZED + return tmp_dir end @testset "Single stage sequential tests" begin for in_memory in (true, false), rebuild in (true, false) - test_single_stage_sequential(in_memory, rebuild) + test_single_stage_sequential(in_memory, rebuild, false) end end +@testset "Test model export at each solve" begin + folder = test_single_stage_sequential(true, false, true) + test_path = + joinpath(folder, "consecutive", "problems", "ED", "optimization_model_exports") + @test ispath(test_path) + @test length(readdir(test_path)) == 2 +end + function test_2_stage_decision_models_with_feedforwards(in_memory) template_uc = get_template_basic_uc_simulation() template_ed = get_template_nomin_ed_simulation()