diff --git a/Project.toml b/Project.toml index 79f3213b8..5670f304e 100644 --- a/Project.toml +++ b/Project.toml @@ -19,6 +19,7 @@ JuMP = "4076af6c-e467-56ae-b986-b466b2749572" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" +PowerFlows = "94fada2c-fd9a-4e89-8d82-81405f5cb4f6" PowerModels = "c36e90e8-916a-50a6-bd94-075b64ef4655" PowerNetworkMatrices = "bed98974-b02a-5e2f-9fe0-a103f5c450dd" PowerSystems = "bcd98974-b02a-5e2f-9ee0-a103f5c450dd" @@ -45,6 +46,7 @@ Logging = "1" MathOptInterface = "1" PowerModels = "^0.21" PowerNetworkMatrices = "^0.11" +PowerFlows = "0.7" PowerSystems = "^4.4" PrettyTables = "2" ProgressMeter = "^1.5" diff --git a/src/PowerSimulations.jl b/src/PowerSimulations.jl index 0a36addcd..fcc20d905 100644 --- a/src/PowerSimulations.jl +++ b/src/PowerSimulations.jl @@ -238,6 +238,10 @@ export PieceWiseLinearCostVariable export TimeDurationOn export TimeDurationOff export PowerOutput +export PowerFlowVoltageAngle +export PowerFlowVoltageMagnitude +export PowerFlowLineReactivePower +export PowerFlowLineActivePower # Constraints export AbsoluteValueConstraint @@ -342,6 +346,7 @@ import LinearAlgebra import JSON3 import PowerSystems import InfrastructureSystems +import PowerFlows import PowerNetworkMatrices import PowerNetworkMatrices: PTDF, VirtualPTDF export PTDF @@ -447,6 +452,7 @@ const MOI = MathOptInterface const MOIU = MathOptInterface.Utilities const MOPFM = MOI.FileFormats.Model const PNM = PowerNetworkMatrices +const PFS = PowerFlows const TS = TimeSeries ################################################################################ @@ -490,6 +496,7 @@ include("core/results_by_time.jl") # Order Required include("operation/problem_template.jl") +include("core/power_flow_data_wrapper.jl") include("core/optimization_container.jl") include("core/store_common.jl") include("initial_conditions/initial_condition_chronologies.jl") @@ -583,6 +590,7 @@ include("network_models/pm_translator.jl") include("network_models/network_slack_variables.jl") include("network_models/area_balance_model.jl") include("network_models/hvdc_networks.jl") +include("network_models/power_flow_evaluation.jl") include("initial_conditions/initialization.jl") diff --git a/src/core/auxiliary_variables.jl b/src/core/auxiliary_variables.jl index c2baca428..50adb7ba9 100644 --- a/src/core/auxiliary_variables.jl +++ b/src/core/auxiliary_variables.jl @@ -13,4 +13,35 @@ Auxiliary Variable for Thermal Generation Models that solve for power above min """ struct PowerOutput <: AuxVariableType end +""" +Auxiliary Variables that are calculated using a `PowerFlowEvaluationModel` +""" +abstract type PowerFlowAuxVariableType <: AuxVariableType end + +""" +Auxiliary Variable for the bus angle results from power flow evaluation +""" +struct PowerFlowVoltageAngle <: PowerFlowAuxVariableType end + +""" +Auxiliary Variable for the bus voltage magnitued results from power flow evaluation +""" +struct PowerFlowVoltageMagnitude <: PowerFlowAuxVariableType end + +""" +Auxiliary Variable for the line reactive flow from power flow evaluation +""" +struct PowerFlowLineReactivePower <: PowerFlowAuxVariableType end + +""" +Auxiliary Variable for the line active flow from power flow evaluation +""" +struct PowerFlowLineActivePower <: PowerFlowAuxVariableType end + convert_result_to_natural_units(::Type{PowerOutput}) = true +convert_result_to_natural_units(::Type{PowerFlowLineReactivePower}) = true +convert_result_to_natural_units(::Type{PowerFlowLineActivePower}) = true + +"Whether the auxiliary variable is calculated using a `PowerFlowEvaluationModel`" +is_from_power_flow(::Type{<:AuxVariableType}) = false +is_from_power_flow(::Type{<:PowerFlowAuxVariableType}) = true diff --git a/src/core/network_model.jl b/src/core/network_model.jl index eaaa0154f..e017a2af1 100644 --- a/src/core/network_model.jl +++ b/src/core/network_model.jl @@ -34,6 +34,7 @@ mutable struct NetworkModel{T <: PM.AbstractPowerModel} duals::Vector{DataType} radial_network_reduction::PNM.RadialNetworkReduction reduce_radial_branches::Bool + power_flow_evaluation::Vector{PFS.PowerFlowEvaluationModel} subsystem::Union{Nothing, String} modeled_branch_types::Vector{DataType} @@ -44,6 +45,7 @@ mutable struct NetworkModel{T <: PM.AbstractPowerModel} reduce_radial_branches = false, subnetworks = Dict{Int, Set{Int}}(), duals = Vector{DataType}(), + power_flow_evaluation = Vector{PFS.PowerFlowEvaluationModel}[], ) where {T <: PM.AbstractPowerModel} _check_pm_formulation(T) new{T}( @@ -54,6 +56,7 @@ mutable struct NetworkModel{T <: PM.AbstractPowerModel} duals, PNM.RadialNetworkReduction(), reduce_radial_branches, + power_flow_evaluation, nothing, Vector{DataType}(), ) @@ -70,6 +73,7 @@ get_reference_buses(m::NetworkModel{T}) where {T <: PM.AbstractPowerModel} = collect(keys(m.subnetworks)) get_subnetworks(m::NetworkModel) = m.subnetworks get_bus_area_map(m::NetworkModel) = m.bus_area_map +get_power_flow_evaluation(m::NetworkModel) = m.power_flow_evaluation has_subnetworks(m::NetworkModel) = !isempty(m.bus_area_map) get_subsystem(m::NetworkModel) = m.subsystem diff --git a/src/core/optimization_container.jl b/src/core/optimization_container.jl index b94438f9b..ea295a561 100644 --- a/src/core/optimization_container.jl +++ b/src/core/optimization_container.jl @@ -78,6 +78,7 @@ mutable struct OptimizationContainer <: IS.Optimization.AbstractOptimizationCont built_for_recurrent_solves::Bool metadata::IS.Optimization.OptimizationContainerMetadata default_time_series_type::Type{<:PSY.TimeSeriesData} + power_flow_evaluation_data::Vector{PowerFlowEvaluationData} end function OptimizationContainer( @@ -120,6 +121,7 @@ function OptimizationContainer( false, IS.Optimization.OptimizationContainerMetadata(), T, + Vector{PowerFlowEvaluationData}[], ) end @@ -157,6 +159,8 @@ get_jump_model(container::OptimizationContainer) = container.JuMPmodel get_metadata(container::OptimizationContainer) = container.metadata get_optimizer_stats(container::OptimizationContainer) = container.optimizer_stats get_parameters(container::OptimizationContainer) = container.parameters +get_power_flow_evaluation_data(container::OptimizationContainer) = + container.power_flow_evaluation_data get_resolution(container::OptimizationContainer) = get_resolution(container.settings) get_settings(container::OptimizationContainer) = container.settings get_time_steps(container::OptimizationContainer) = container.time_steps @@ -170,6 +174,12 @@ is_synchronized(container::OptimizationContainer) = set_time_steps!(container::OptimizationContainer, time_steps::UnitRange{Int64}) = container.time_steps = time_steps +function reset_power_flow_is_solved!(container::OptimizationContainer) + for pf_e_data in get_power_flow_evaluation_data(container) + pf_e_data.is_solved = false + end +end + function has_container_key( container::OptimizationContainer, ::Type{T}, @@ -749,6 +759,7 @@ function build_impl!( @debug "Total operation count $(PSI.get_jump_model(container).operator_counter)" _group = LOG_GROUP_OPTIMIZATION_CONTAINER + add_power_flow_data!(container, get_power_flow_evaluation(transmission_model), sys) check_optimization_container(container) return end @@ -800,6 +811,8 @@ function solve_impl!(container::OptimizationContainer, system::PSY.System) end end + # Order is important because if a dual is needed then it could move the results to the + # temporary primal container _, optimizer_stats.timed_calculate_aux_variables = @timed calculate_aux_variables!(container, system) @@ -809,9 +822,7 @@ function solve_impl!(container::OptimizationContainer, system::PSY.System) _, optimizer_stats.timed_calculate_dual_variables = @timed calculate_dual_variables!(container, system, is_milp(container)) - status = RunStatus.SUCCESSFULLY_FINALIZED - - return status + return RunStatus.SUCCESSFULLY_FINALIZED end function compute_conflict!(container::OptimizationContainer) @@ -1636,8 +1647,24 @@ function deserialize_key(container::OptimizationContainer, name::AbstractString) end function calculate_aux_variables!(container::OptimizationContainer, system::PSY.System) - aux_vars = get_aux_variables(container) - for key in keys(aux_vars) + aux_var_keys = keys(get_aux_variables(container)) + pf_aux_var_keys = filter(is_from_power_flow ∘ get_entry_type, aux_var_keys) + non_pf_aux_var_keys = setdiff(aux_var_keys, pf_aux_var_keys) + # We should only have power flow aux vars if we have power flow evaluators + @assert isempty(pf_aux_var_keys) || !isempty(get_power_flow_evaluation_data(container)) + + reset_power_flow_is_solved!(container) + # Power flow-related aux vars get calculated once per power flow + for (i, pf_e_data) in enumerate(get_power_flow_evaluation_data(container)) + @debug "Processing power flow $i" + solve_powerflow!(pf_e_data, container) + for key in pf_aux_var_keys + calculate_aux_variable_value!(container, key, system) + end + end + + # Other aux vars get calculated once at the end + for key in non_pf_aux_var_keys calculate_aux_variable_value!(container, key, system) end return RunStatus.SUCCESSFULLY_FINALIZED @@ -1914,3 +1941,14 @@ function get_time_series_initial_values!( ) return ts_values end + +lookup_value(container::OptimizationContainer, key::VariableKey) = + get_variable(container, key) +lookup_value(container::OptimizationContainer, key::ParameterKey) = + calculate_parameter_values(get_parameter(container, key)) +lookup_value(container::OptimizationContainer, key::AuxVarKey) = + get_aux_variable(container, key) +lookup_value(container::OptimizationContainer, key::ExpressionKey) = + get_expression(container, key) +lookup_value(container::OptimizationContainer, key::ConstraintKey) = + get_constraint(container, key) diff --git a/src/core/parameters.jl b/src/core/parameters.jl index 2946dd68f..fc850b571 100644 --- a/src/core/parameters.jl +++ b/src/core/parameters.jl @@ -47,6 +47,7 @@ function add_component_name!(attr::TimeSeriesAttributes, name::String, uuid::Str return end +get_component_names(attr::TimeSeriesAttributes) = keys(attr.component_name_to_ts_uuid) _get_ts_uuid(attr::TimeSeriesAttributes, name) = attr.component_name_to_ts_uuid[name] struct VariableValueAttributes{T <: OptimizationContainerKey} <: ParameterAttributes diff --git a/src/core/power_flow_data_wrapper.jl b/src/core/power_flow_data_wrapper.jl new file mode 100644 index 000000000..08da670cd --- /dev/null +++ b/src/core/power_flow_data_wrapper.jl @@ -0,0 +1,22 @@ +mutable struct PowerFlowEvaluationData{T <: PFS.PowerFlowContainer} + power_flow_data::T + """ + Records which PSI keys are read as input to the power flow and how the data are mapped. + For `PowerFlowData`, values are `Dict{String, Int64}` mapping component name to matrix + index of bus; for `SystemPowerFlowContainer`, values are Dict{Union{String, Int64}, + Union{String, Int64}} mapping component name/bus number to component name/bus number. + """ + input_key_map::Dict{<:OptimizationContainerKey, <:Any} + is_solved::Bool +end + +function PowerFlowEvaluationData(power_flow_data::T) where {T <: PFS.PowerFlowContainer} + return PowerFlowEvaluationData{T}( + power_flow_data, + Dict{OptimizationContainerKey, Nothing}(), + false, + ) +end + +get_power_flow_data(ped::PowerFlowEvaluationData) = ped.power_flow_data +get_input_key_map(ped::PowerFlowEvaluationData) = ped.input_key_map diff --git a/src/initial_conditions/initialization.jl b/src/initial_conditions/initialization.jl index cff6e81ed..025ac623d 100644 --- a/src/initial_conditions/initialization.jl +++ b/src/initial_conditions/initialization.jl @@ -12,6 +12,8 @@ function get_initial_conditions_template(model::OperationModel) network_model.radial_network_reduction = get_radial_network_reduction(get_network_model(model.template)) network_model.subnetworks = get_subnetworks(get_network_model(model.template)) + # Initialization does not support PowerFlow evaluation + network_model.power_flow_evaluation = Vector{PFS.PowerFlowEvaluationModel}[] bus_area_map = get_bus_area_map(get_network_model(model.template)) if !isempty(bus_area_map) diff --git a/src/network_models/power_flow_evaluation.jl b/src/network_models/power_flow_evaluation.jl new file mode 100644 index 000000000..5768b1241 --- /dev/null +++ b/src/network_models/power_flow_evaluation.jl @@ -0,0 +1,352 @@ +function _add_aux_variables!( + container::OptimizationContainer, + component_map::Dict{Type{<:AuxVariableType}, <:Set{<:Tuple{DataType, Any}}}, +) + for (var_type, components) in pairs(component_map) + component_types = unique(first.(components)) + for component_type in component_types + component_names = [v for (k, v) in components if k <: component_type] + sort!(component_names) + add_aux_variable_container!( + container, + var_type(), + component_type, + component_names, + get_time_steps(container), + ) + end + end +end + +# Trait that determines what keys serve as input to each type of power flow, if they exist +pf_input_keys(::PFS.PowerFlowData) = + [ActivePowerVariable, PowerOutput, ActivePowerTimeSeriesParameter] +pf_input_keys(::PFS.PSSEExporter) = + [ActivePowerVariable, PowerOutput, ActivePowerTimeSeriesParameter, + PowerFlowVoltageAngle, PowerFlowVoltageMagnitude] + +# Maps the StaticInjection component type by name to the +# index in the PowerFlow data arrays going from Bus number to bus index +function _make_temp_component_map(pf_data::PFS.PowerFlowData, sys::PSY.System) + temp_component_map = Dict{DataType, Dict{String, Int}}() + available_injectors = PSY.get_components(PSY.get_available, PSY.StaticInjection, sys) + bus_lookup = PFS.get_bus_lookup(pf_data) + for comp in available_injectors + comp_type = typeof(comp) + bus_dict = get!(temp_component_map, comp_type, Dict{String, Int}()) + bus_number = PSY.get_number(PSY.get_bus(comp)) + bus_dict[get_name(comp)] = bus_lookup[bus_number] + end + return temp_component_map +end + +_get_temp_component_map_lhs(comp::PSY.Component) = PSY.get_name(comp) +_get_temp_component_map_lhs(comp::PSY.Bus) = PSY.get_number(comp) + +# Creates dicts of components by type +function _make_temp_component_map(::PFS.SystemPowerFlowContainer, sys::PSY.System) + temp_component_map = + Dict{DataType, Dict{Union{String, Int64}, String}}() + # TODO don't hardcode the types here, handle get_available more elegantly, likely `ComponentSelector` use case + relevant_components = vcat( + collect.([ + PSY.get_components(PSY.get_available, PSY.StaticInjection, sys), + PSY.get_components(Union{PSY.Bus, PSY.Branch}, sys)], + )..., + ) + for comp_type in unique(typeof.(relevant_components)) + # NOTE we avoid using bus numbers here because PSY.get_bus(system, number) is O(n) + temp_component_map[comp_type] = + Dict( + _get_temp_component_map_lhs(c) => PSY.get_name(c) for + c in relevant_components if c isa comp_type + ) + end + return temp_component_map +end + +function _make_pf_input_map!( + pf_e_data::PowerFlowEvaluationData, + container::OptimizationContainer, + sys::PSY.System, +) + pf_data = get_power_flow_data(pf_e_data) + temp_component_map = _make_temp_component_map(pf_data, sys) + map_type = valtype(temp_component_map) # Dict{String, Int} for PowerFlowData, Dict{Union{String, Int64}, String} for SystemPowerFlowContainer + input_keys = pf_input_keys(pf_data) + + # Second map that persists to store the bus index that the variable + # has to be added/substracted to in the power flow data dictionary + pf_data_opt_container_map = Dict{OptimizationContainerKey, map_type}() + added_injection_types = DataType[] + for (key, val) in Iterators.flatten([ + get_variables(container), + get_aux_variables(container), + get_parameters(container), + ]) + # Skip irrelevant keys + (get_entry_type(key) in input_keys) || continue + + comp_type = get_component_type(key) + # Skip types that have already been handled (prefer variable over aux variable, aux variable over parameter) + (comp_type in added_injection_types) && continue + push!(added_injection_types, comp_type) + + name_bus_ix_map = map_type() + comp_names = + (key isa ParameterKey) ? get_component_names(get_attributes(val)) : axes(val)[1] + for comp_name in comp_names + name_bus_ix_map[comp_name] = temp_component_map[comp_type][comp_name] + end + pf_data_opt_container_map[key] = name_bus_ix_map + end + + pf_e_data.input_key_map = pf_data_opt_container_map + return +end + +# Trait that determines what branch aux vars we can get from each PowerFlowContainer +branch_aux_vars(::PFS.ACPowerFlowData) = + [PowerFlowLineActivePower, PowerFlowLineReactivePower] +branch_aux_vars(::PFS.ABAPowerFlowData) = [PowerFlowLineActivePower] +branch_aux_vars(::PFS.PTDFPowerFlowData) = [PowerFlowLineActivePower] +branch_aux_vars(::PFS.vPTDFPowerFlowData) = [PowerFlowLineActivePower] +branch_aux_vars(::PFS.PSSEExporter) = DataType[] + +# Same for bus aux vars +bus_aux_vars(::PFS.ACPowerFlowData) = [PowerFlowVoltageAngle, PowerFlowVoltageMagnitude] +bus_aux_vars(::PFS.ABAPowerFlowData) = [PowerFlowVoltageAngle] +bus_aux_vars(::PFS.PTDFPowerFlowData) = DataType[] +bus_aux_vars(::PFS.vPTDFPowerFlowData) = DataType[] +bus_aux_vars(::PFS.PSSEExporter) = DataType[] + +_get_branch_component_tuples(pfd::PFS.PowerFlowData) = + zip(PFS.get_branch_type(pfd), keys(PFS.get_branch_lookup(pfd))) + +_get_branch_component_tuples(pfd::PFS.SystemPowerFlowContainer) = + [(typeof(c), get_name(c)) for c in PSY.get_components(PSY.Branch, PFS.get_system(pfd))] + +_get_bus_component_tuples(pfd::PFS.PowerFlowData) = + tuple.(PSY.ACBus, keys(PFS.get_bus_lookup(pfd))) # get_bus_type returns a ACBusTypes, not the DataType we need here + +_get_bus_component_tuples(pfd::PFS.SystemPowerFlowContainer) = + [ + (typeof(c), PSY.get_number(c)) for + c in PSY.get_components(PSY.Bus, PFS.get_system(pfd)) + ] + +function add_power_flow_data!( + container::OptimizationContainer, + evaluators::Vector{PFS.PowerFlowEvaluationModel}, + sys::PSY.System, +) + container.power_flow_evaluation_data = Vector{PowerFlowEvaluationData}() + sizehint!(container.power_flow_evaluation_data, length(evaluators)) + # For each output key, what components are we working with? + branch_aux_var_components = + Dict{Type{<:AuxVariableType}, Set{Tuple{<:DataType, String}}}() + bus_aux_var_components = Dict{Type{<:AuxVariableType}, Set{Tuple{<:DataType, <:Int}}}() + for evaluator in evaluators + @info "Building PowerFlow evaluator using $(evaluator)" + pf_data = PFS.make_power_flow_container(evaluator, sys; + time_steps = length(get_time_steps(container))) + pf_e_data = PowerFlowEvaluationData(pf_data) + my_branch_aux_vars = branch_aux_vars(pf_data) + my_bus_aux_vars = bus_aux_vars(pf_data) + + my_branch_components = _get_branch_component_tuples(pf_data) + for branch_aux_var in my_branch_aux_vars + to_add_to = get!( + branch_aux_var_components, + branch_aux_var, + Set{Tuple{<:DataType, String}}(), + ) + push!.(Ref(to_add_to), my_branch_components) + end + + my_bus_components = _get_bus_component_tuples(pf_data) + for bus_aux_var in my_bus_aux_vars + to_add_to = + get!(bus_aux_var_components, bus_aux_var, Set{Tuple{<:DataType, <:Int}}()) + push!.(Ref(to_add_to), my_bus_components) + end + push!(container.power_flow_evaluation_data, pf_e_data) + end + + _add_aux_variables!(container, branch_aux_var_components) + _add_aux_variables!(container, bus_aux_var_components) + + # Make the input maps after adding aux vars so output of one power flow can be input of another + for pf_e_data in get_power_flow_evaluation_data(container) + _make_pf_input_map!(pf_e_data, container, sys) + end + return +end + +# How to update the PowerFlowData given a component type. A bit duplicative of code in PowerFlows.jl. +_update_pf_data_component!( + pf_data::PFS.PowerFlowData, + ::Type{<:PSY.StaticInjection}, + index, + t, + value, +) = (pf_data.bus_activepower_injection[index, t] += value) +_update_pf_data_component!( + pf_data::PFS.PowerFlowData, + ::Type{<:PSY.ElectricLoad}, + index, + t, + value, +) = (pf_data.bus_activepower_withdrawals[index, t] += value) + +function _write_value_to_pf_data!( + pf_data::PFS.PowerFlowData, + container::OptimizationContainer, + key::OptimizationContainerKey, + component_map) + result = lookup_value(container, key) + for (device_name, index) in component_map + injection_values = result[device_name, :] + for t in get_time_steps(container) + value = jump_value(injection_values[t]) + _update_pf_data_component!(pf_data, get_component_type(key), index, t, value) + end + end + return +end + +function update_pf_data!( + pf_e_data::PowerFlowEvaluationData{<:PFS.PowerFlowData}, + container::OptimizationContainer, +) + pf_data = get_power_flow_data(pf_e_data) + PFS.clear_injection_data!(pf_data) + input_map = get_input_key_map(pf_e_data) + for (key, component_map) in input_map + _write_value_to_pf_data!(pf_data, container, key, component_map) + end + return +end + +_update_component( + ::Type{<:Union{ActivePowerVariable, PowerOutput, ActivePowerTimeSeriesParameter}}, + comp::PSY.Component, + value, +) = (comp.active_power = value) +# Sign is flipped for loads (TODO can we rely on some existing function that encodes this information?) +_update_component( + ::Type{<:Union{ActivePowerVariable, PowerOutput, ActivePowerTimeSeriesParameter}}, + comp::PSY.ElectricLoad, + value, +) = (comp.active_power = -value) +_update_component(::Type{PowerFlowVoltageAngle}, comp::PSY.Component, value) = + comp.angle = value +_update_component(::Type{PowerFlowVoltageMagnitude}, comp::PSY.Component, value) = + comp.magnitude = value + +function update_pf_system!( + sys::PSY.System, + container::OptimizationContainer, + input_map::Dict{<:OptimizationContainerKey, <:Any}, + time_step::Int, +) + for (key, component_map) in input_map + result = lookup_value(container, key) + for (device_id, device_name) in component_map + injection_values = result[device_id, :] + comp = PSY.get_component(get_component_type(key), sys, device_name) + val = jump_value(injection_values[time_step]) + _update_component(get_entry_type(key), comp, val) + end + end +end + +""" +Update a `PowerFlowEvaluationData` containing a `PowerFlowContainer` that does not +`supports_multi_period` using a single `time_step` of the `OptimizationContainer`. To +properly keep track of outer step number, time steps must be passed in sequentially, +starting with 1. +""" +function update_pf_data!( + pf_e_data::PowerFlowEvaluationData{PFS.PSSEExporter}, + container::OptimizationContainer, + time_step::Int, +) + pf_data = get_power_flow_data(pf_e_data) + input_map = get_input_key_map(pf_e_data) + update_pf_system!(PFS.get_system(pf_data), container, input_map, time_step) + if !isnothing(pf_data.step) + outer_step, _ = pf_data.step + # time_step == 1 means we have rolled over to a new outer step + # (TODO it works but seems a little brittle, consider redesigning) + (time_step == 1) && (outer_step += 1) + pf_data.step = (outer_step, time_step) + end + return +end + +"Fetch the most recently solved `PowerFlowEvaluationData`" +function latest_solved_power_flow_evaluation_data(container::OptimizationContainer) + datas = get_power_flow_evaluation_data(container) + return datas[findlast(x -> x.is_solved, datas)] +end + +function solve_powerflow!( + pf_e_data::PowerFlowEvaluationData, + container::OptimizationContainer) + pf_data = get_power_flow_data(pf_e_data) + if PFS.supports_multi_period(pf_data) + update_pf_data!(pf_e_data, container) + PFS.solve_powerflow!(pf_data) + else + for t in get_time_steps(container) + update_pf_data!(pf_e_data, container, t) + PFS.solve_powerflow!(pf_data) + end + end + pf_e_data.is_solved = true + return +end + +# Currently nothing to write back to the optimization container from a PSSEExporter +calculate_aux_variable_value!(::OptimizationContainer, + ::AuxVarKey{T, <:Any} where {T <: PowerFlowAuxVariableType}, + ::PSY.System, ::PowerFlowEvaluationData{PFS.PSSEExporter}) = nothing + +_get_pf_result(::Type{PowerFlowVoltageAngle}, pf_data::PFS.PowerFlowData) = + PFS.get_bus_angles(pf_data) +_get_pf_result(::Type{PowerFlowVoltageMagnitude}, pf_data::PFS.PowerFlowData) = + PFS.get_bus_magnitude(pf_data) +_get_pf_result(::Type{PowerFlowLineActivePower}, pf_data::PFS.PowerFlowData) = + PFS.get_branch_flow_values(pf_data) +# TODO implement method for PowerFlowLineReactivePower -- I don't think we have a PowerFlowData field for this? +# _fetch_pf_result(pf_data::PFS.PowerFlowData, ::Type{PowerFlowLineActivePower}) = ... + +_get_pf_lookup(::Type{<:PSY.Bus}, pf_data::PFS.PowerFlowData) = PFS.get_bus_lookup(pf_data) +_get_pf_lookup(::Type{<:PSY.Branch}, pf_data::PFS.PowerFlowData) = + PFS.get_branch_lookup(pf_data) + +function calculate_aux_variable_value!(container::OptimizationContainer, + key::AuxVarKey{T, U}, + system::PSY.System, pf_e_data::PowerFlowEvaluationData{<:PFS.PowerFlowData}, +) where {T <: PowerFlowAuxVariableType, U} + @debug "Updating $key from PowerFlowData" + pf_data = get_power_flow_data(pf_e_data) + src = _get_pf_result(T, pf_data) + lookup = _get_pf_lookup(U, pf_data) + dest = get_aux_variable(container, key) + for component_id in axes(dest, 1) # these are bus numbers or branch names + dest[component_id, :] = src[lookup[component_id], :] + end + return +end + +function calculate_aux_variable_value!(container::OptimizationContainer, + key::AuxVarKey{T, <:Any} where {T <: PowerFlowAuxVariableType}, + system::PSY.System) + pf_e_data = latest_solved_power_flow_evaluation_data(container) + pf_data = get_power_flow_data(pf_e_data) + # Skip the aux vars that the current power flow isn't meant to update + (key in branch_aux_vars(pf_data) || key in bus_aux_vars(pf_data)) && return + calculate_aux_variable_value!(container, key, system, pf_e_data) +end diff --git a/src/services_models/transmission_interface.jl b/src/services_models/transmission_interface.jl index b76d994ae..797300be0 100644 --- a/src/services_models/transmission_interface.jl +++ b/src/services_models/transmission_interface.jl @@ -55,13 +55,13 @@ function add_constraints!(container::OptimizationContainer, model::ServiceModel{T, ConstantMaxInterfaceFlow}, ) where {T <: PSY.TransmissionInterface} expr = get_expression(container, InterfaceTotalFlow(), T) - interfaces, timesteps = axes(expr) + interfaces, time_steps = axes(expr) constraint_container_ub = lazy_container_addition!( container, InterfaceFlowLimit(), T, interfaces, - timesteps; + time_steps; meta = "ub", ) constraint_container_lb = lazy_container_addition!( @@ -69,12 +69,12 @@ function add_constraints!(container::OptimizationContainer, InterfaceFlowLimit(), T, interfaces, - timesteps; + time_steps; meta = "lb", ) int_name = PSY.get_name(interface) min_flow, max_flow = PSY.get_active_power_flow_limits(interface) - for t in timesteps + for t in time_steps constraint_container_ub[int_name, t] = JuMP.@constraint(get_jump_model(container), expr[int_name, t] <= max_flow) constraint_container_lb[int_name, t] = diff --git a/src/simulation/simulation_problem_results.jl b/src/simulation/simulation_problem_results.jl index faa327d54..f5c1607c4 100644 --- a/src/simulation/simulation_problem_results.jl +++ b/src/simulation/simulation_problem_results.jl @@ -79,13 +79,16 @@ function set_results_timestamps!( result.results_timestamps = results_timestamps end -list_result_keys(res::SimulationProblemResults, ::AuxVarKey) = list_aux_variable_keys(res) -list_result_keys(res::SimulationProblemResults, ::ConstraintKey) = list_dual_keys(res) +list_result_keys(res::SimulationProblemResults, ::AuxVarKey) = + list_aux_variable_keys(res) +list_result_keys(res::SimulationProblemResults, ::ConstraintKey) = + list_dual_keys(res) list_result_keys(res::SimulationProblemResults, ::ExpressionKey) = list_expression_keys(res) list_result_keys(res::SimulationProblemResults, ::ParameterKey) = list_parameter_keys(res) -list_result_keys(res::SimulationProblemResults, ::VariableKey) = list_variable_keys(res) +list_result_keys(res::SimulationProblemResults, ::VariableKey) = + list_variable_keys(res) get_cached_results(res::SimulationProblemResults, ::Type{<:AuxVarKey}) = get_cached_aux_variables(res) diff --git a/src/utils/jump_utils.jl b/src/utils/jump_utils.jl index 129414c0d..8d7e04ef2 100644 --- a/src/utils/jump_utils.jl +++ b/src/utils/jump_utils.jl @@ -106,9 +106,9 @@ function _to_matrix( array::SparseAxisArray{T, N, K}, columns, ) where {T, N, K <: NTuple{N, Any}} - timesteps = Set{Int}(k[N] for k in keys(array.data)) - data = Matrix{Float64}(undef, length(timesteps), length(columns)) - for (ix, col) in enumerate(columns), t in timesteps + time_steps = Set{Int}(k[N] for k in keys(array.data)) + data = Matrix{Float64}(undef, length(time_steps), length(columns)) + for (ix, col) in enumerate(columns), t in time_steps data[t, ix] = array.data[(col..., t)] end return data diff --git a/test/Project.toml b/test/Project.toml index e15b719fb..f82a28131 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -15,6 +15,7 @@ JuMP = "4076af6c-e467-56ae-b986-b466b2749572" Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" Memento = "f28f55f0-a522-5efc-85c2-fe41dfb9b2d9" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" +PowerFlows = "94fada2c-fd9a-4e89-8d82-81405f5cb4f6" PowerModels = "c36e90e8-916a-50a6-bd94-075b64ef4655" PowerNetworkMatrices = "bed98974-b02a-5e2f-9fe0-a103f5c450dd" PowerSimulations = "e690365d-45e2-57bb-ac84-44ba829e73c4" @@ -33,5 +34,4 @@ UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" [compat] Ipopt = "=1.4.0" -PowerSystemCaseBuilder = "^1.2.0" julia = "^1.6" diff --git a/test/includes.jl b/test/includes.jl index 647e724e9..ba8c7baeb 100644 --- a/test/includes.jl +++ b/test/includes.jl @@ -8,6 +8,7 @@ using HydroPowerSimulations import PowerSystemCaseBuilder: PSITestSystems using PowerNetworkMatrices using StorageSystemsSimulations +using PowerFlows # Test Packages using Test diff --git a/test/runtests.jl b/test/runtests.jl index 9c0568af4..2778d93da 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -11,28 +11,31 @@ Aqua.test_unbound_args(PowerSimulations) const LOG_FILE = "power-simulations-test.log" const DISABLED_TEST_FILES = [ -# "test_basic_model_structs.jl", -# "test_device_branch_constructors.jl", -# "test_device_hydro_generation_constructors.jl", -# "test_device_load_constructors.jl", -# "test_device_hybrid_generation_constructors.jl", -# "test_device_renewable_generation_constructors.jl", -# "test_device_storage_constructors.jl", -# "test_device_thermal_generation_constructors.jl", -# "test_jump_model_utils.jl", -# "test_model_decision.jl", -# "test_problem_template.jl", -# "test_model_emulation.jl", -# "test_network_constructors.jl", -# "test_services_constructor.jl", -# "test_simulation_models.jl", -# "test_simulation_sequence.jl", -# "test_simulation_build.jl", -# "test_initialization_problem.jl", -# "test_simulation_execute.jl", -# "test_simulation_results.jl", -# "test_simulation_results_export.jl", -# "test_simulation_store.jl", +# "test_basic_model_structs.jl" +# "test_device_branch_constructors.jl" +# "test_device_hvdc.jl" +# "test_device_load_constructors.jl" +# "test_device_renewable_generation_constructors.jl" +# "test_device_thermal_generation_constructors.jl" +# "test_formulation_combinations.jl" +# "test_ic_reconciliation.jl" +# "test_initialization_problem.jl" +# "test_model_decision.jl" +# "test_model_emulation.jl" +# "test_network_constructors.jl" +# "test_print.jl" +# "test_problem_template.jl" +# "test_recorder_events.jl" +# "test_services_constructor.jl" +# "test_simulation_build.jl" +# "test_simulation_execute.jl" +# "test_simulation_models.jl" +# "test_simulation_partitions.jl" +# "test_simulation_results.jl" +# "test_simulation_results_export.jl" +# "test_simulation_sequence.jl" +# "test_simulation_store.jl" +# "test_utils.jl" ] LOG_LEVELS = Dict( diff --git a/test/test_basic_model_structs.jl b/test/test_basic_model_structs.jl index a5c37eb5c..02cd696a8 100644 --- a/test/test_basic_model_structs.jl +++ b/test/test_basic_model_structs.jl @@ -6,6 +6,11 @@ end @testset "NetworkModel Tests" begin @test_throws ArgumentError NetworkModel(PM.AbstractPowerModel) + @test NetworkModel( + PTDFPowerModel; + use_slacks = true, + power_flow_evaluation = [DCPowerFlow(), PSSEExportPowerFlow(:v33, "exports", true)], + ) isa NetworkModel end #= diff --git a/test/test_network_constructors.jl b/test/test_network_constructors.jl index 2b90bb648..95c0bf452 100644 --- a/test/test_network_constructors.jl +++ b/test/test_network_constructors.jl @@ -632,7 +632,6 @@ end # These models are easier to test due to their lossless nature @testset "StandardPTDF/DCPPowerModel Radial Branches Test" begin new_sys = PSB.build_system(PSITestSystems, "c_sys5_radial") - for net_model in [DCPPowerModel, PTDFPowerModel] template_uc = template_unit_commitment(; network = NetworkModel(net_model; diff --git a/test/test_simulation_results.jl b/test/test_simulation_results.jl index d3ec74477..5ab275dc8 100644 --- a/test/test_simulation_results.jl +++ b/test/test_simulation_results.jl @@ -155,21 +155,31 @@ function run_simulation( export_path; in_memory = false, system_to_file = true, + uc_network_model = nothing, + ed_network_model = nothing, ) template_uc = get_template_basic_uc_simulation() template_ed = get_template_nomin_ed_simulation() + isnothing(uc_network_model) && ( + uc_network_model = + NetworkModel(CopperPlatePowerModel; duals = [CopperPlateBalanceConstraint]) + ) + isnothing(ed_network_model) && ( + ed_network_model = + NetworkModel( + CopperPlatePowerModel; + duals = [CopperPlateBalanceConstraint], + use_slacks = true, + ) + ) set_device_model!(template_ed, InterruptiblePowerLoad, StaticPowerLoad) set_network_model!( template_uc, - NetworkModel(CopperPlatePowerModel; duals = [CopperPlateBalanceConstraint]), + uc_network_model, ) set_network_model!( template_ed, - NetworkModel( - CopperPlatePowerModel; - duals = [CopperPlateBalanceConstraint], - use_slacks = true, - ), + ed_network_model, ) models = SimulationModels(; decision_models = [ @@ -977,3 +987,59 @@ end test_decision_problem_results(results, sys_ed, sys_uc, in_memory) test_emulation_problem_results(results, in_memory) end + +function load_pf_export(root, export_subdir) + raw_path, md_path = get_psse_export_paths(export_subdir) + sys = System(joinpath(root, raw_path), JSON3.read(joinpath(root, md_path), Dict)) + # TODO I think the necessity of this might speak to a unit misinterpretation somewhere + set_units_base_system!(sys, "NATURAL_UNITS") + return sys +end + +@testset "Test power flow in the loop" begin + file_path = mktempdir(; cleanup = true) + export_path = mktempdir(; cleanup = true) + pf_path = mktempdir(; cleanup = true) + c_sys5_hy_uc = PSB.build_system(PSITestSystems, "c_sys5_hy_uc") + c_sys5_hy_ed = PSB.build_system(PSITestSystems, "c_sys5_hy_ed") + sim = run_simulation( + c_sys5_hy_uc, + c_sys5_hy_ed, + file_path, + export_path; + ed_network_model = NetworkModel( + CopperPlatePowerModel; + duals = [CopperPlateBalanceConstraint], + use_slacks = true, + power_flow_evaluation = + [DCPowerFlow(), PSSEExportPowerFlow(:v33, pf_path, true)], + ), + ) + results = SimulationResults(sim) + results_ed = get_decision_problem_results(results, "ED") + thermal_results = first( + values( + PSI.read_results_with_keys(results_ed, + [PSI.VariableKey(ActivePowerVariable, ThermalStandard)]), + ), + ) + first_result = first(thermal_results) + last_result = last(thermal_results) + + @test length(filter(x -> isdir(joinpath(pf_path, x)), readdir(pf_path))) == 48 * 12 + first_export = load_pf_export(pf_path, "export_1_1") + last_export = load_pf_export(pf_path, "export_48_12") + + # Test that the active powers written to the first and last exports line up with the real simulation results + for gen_name in get_name.(get_components(ThermalStandard, c_sys5_hy_ed)) + this_first_result = first_result[gen_name] + this_first_exported = + get_active_power(get_component(ThermalStandard, first_export, gen_name)) + @test isapprox(this_first_result, this_first_exported) + + this_last_result = last_result[gen_name] + this_last_exported = + get_active_power(get_component(ThermalStandard, last_export, gen_name)) + @test isapprox(this_last_result, this_last_exported) + end +end