[DNMY] Support in the loop PowerFlow evaluation

Project.toml

@@ -18,6 +18,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"
@@ -35,15 +36,16 @@ Dates = "1"
 Distributed = "1"
 DocStringExtensions = "~v0.9"
 HDF5 = "~0.17"
-InfrastructureSystems = "2"
+InfrastructureSystems = "^2.4"
 InteractiveUtils = "1"
 JSON3 = "1"
 JuMP = "1"
 LinearAlgebra = "1"
 Logging = "1"
 MathOptInterface = "1"
 PowerModels = "^0.21"
-PowerNetworkMatrices = "^0.11"
+PowerNetworkMatrices = "^0.12"
+PowerFlows = "0.7"
 PowerSystems = "^4.5"
 PrettyTables = "2"
 ProgressMeter = "^1.5"

src/PowerSimulations.jl

@@ -240,6 +240,10 @@ export RateofChangeConstraintSlackDown
 export TimeDurationOn
 export TimeDurationOff
 export PowerOutput
+export PowerFlowVoltageAngle
+export PowerFlowVoltageMagnitude
+export PowerFlowLineReactivePowerFromTo, PowerFlowLineReactivePowerToFrom
+export PowerFlowLineActivePowerFromTo, PowerFlowLineActivePowerToFrom
 
 # Constraints
 export AbsoluteValueConstraint
@@ -345,6 +349,7 @@ import LinearAlgebra
 import JSON3
 import PowerSystems
 import InfrastructureSystems
+import PowerFlows
 import PowerNetworkMatrices
 import PowerNetworkMatrices: PTDF, VirtualPTDF
 export PTDF
@@ -456,6 +461,7 @@ const MOI = MathOptInterface
 const MOIU = MathOptInterface.Utilities
 const MOPFM = MOI.FileFormats.Model
 const PNM = PowerNetworkMatrices
+const PFS = PowerFlows
 const TS = TimeSeries
 
 ################################################################################
@@ -499,6 +505,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")
@@ -593,6 +600,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")
 

src/core/auxiliary_variables.jl

@@ -13,4 +13,51 @@ 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 in the from -> to direction from power flow evaluation
+"""
+struct PowerFlowLineReactivePowerFromTo <: PowerFlowAuxVariableType end
+
+"""
+Auxiliary Variable for the line reactive flow in the to -> from direction from power flow evaluation
+"""
+struct PowerFlowLineReactivePowerToFrom <: PowerFlowAuxVariableType end
+
+"""
+Auxiliary Variable for the line active flow in the from -> to direction from power flow evaluation
+"""
+struct PowerFlowLineActivePowerFromTo <: PowerFlowAuxVariableType end
+
+"""
+Auxiliary Variable for the line active flow in the to -> from direction from power flow evaluation
+"""
+struct PowerFlowLineActivePowerToFrom <: PowerFlowAuxVariableType end
+
 convert_result_to_natural_units(::Type{PowerOutput}) = true
+convert_result_to_natural_units(
+    ::Type{
+        <:Union{
+            PowerFlowLineReactivePowerFromTo, PowerFlowLineReactivePowerToFrom,
+            PowerFlowLineActivePowerFromTo, PowerFlowLineActivePowerToFrom,
+        },
+    },
+) = true
+
+"Whether the auxiliary variable is calculated using a `PowerFlowEvaluationModel`"
+is_from_power_flow(::Type{<:AuxVariableType}) = false
+is_from_power_flow(::Type{<:PowerFlowAuxVariableType}) = true

src/core/network_model.jl

@@ -8,6 +8,10 @@ function _check_pm_formulation(::Type{T}) where {T <: PM.AbstractPowerModel}
     end
 end
 
+_maybe_flatten_pfem(pfem::Vector{PFS.PowerFlowEvaluationModel}) = pfem
+_maybe_flatten_pfem(pfem::PFS.PowerFlowEvaluationModel) =
+    PFS.flatten_power_flow_evaluation_model(pfem)
+
 """
 Establishes the model for the network specified by type.
 
@@ -34,6 +38,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 +49,10 @@ mutable struct NetworkModel{T <: PM.AbstractPowerModel}
         reduce_radial_branches = false,
         subnetworks = Dict{Int, Set{Int}}(),
         duals = Vector{DataType}(),
+        power_flow_evaluation::Union{
+            PFS.PowerFlowEvaluationModel,
+            Vector{PFS.PowerFlowEvaluationModel},
+        } = PFS.PowerFlowEvaluationModel[],
     ) where {T <: PM.AbstractPowerModel}
         _check_pm_formulation(T)
         new{T}(
@@ -54,6 +63,7 @@ mutable struct NetworkModel{T <: PM.AbstractPowerModel}
             duals,
             PNM.RadialNetworkReduction(),
             reduce_radial_branches,
+            _maybe_flatten_pfem(power_flow_evaluation),
             nothing,
             Vector{DataType}(),
         )
@@ -70,6 +80,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
 

src/core/optimization_container.jl

@@ -82,6 +82,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(
@@ -124,6 +125,7 @@ function OptimizationContainer(
         false,
         IS.Optimization.OptimizationContainerMetadata(),
         T,
+        Vector{PowerFlowEvaluationData}[],
     )
 end
 
@@ -161,6 +163,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
@@ -174,6 +178,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},
@@ -753,6 +763,9 @@ function build_impl!(
     @debug "Total operation count $(PSI.get_jump_model(container).operator_counter)" _group =
         LOG_GROUP_OPTIMIZATION_CONTAINER
 
+    TimerOutputs.@timeit BUILD_PROBLEMS_TIMER "Power Flow Initialization" begin
+        add_power_flow_data!(container, get_power_flow_evaluation(transmission_model), sys)
+    end
     check_optimization_container(container)
     return
 end
@@ -804,6 +817,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)
 
@@ -813,9 +828,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)
@@ -1647,8 +1660,26 @@ 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))
+
+    TimerOutputs.@timeit RUN_SIMULATION_TIMER "Power Flow Evaluation" begin
+        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
+    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
@@ -1925,3 +1956,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)

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)
 function _get_ts_uuid(attr::TimeSeriesAttributes, name)
     return attr.component_name_to_ts_uuid[name]
 end

src/core/power_flow_data_wrapper.jl

@@ -0,0 +1,24 @@
+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.
+    The Symbol is a category of data: `:active_power`, `:reactive_power`, etc. The
+    `OptimizationContainerKey` is a source of that data in the `OptimizationContainer`. For
+    `PowerFlowData`, leaf values are `Dict{String, Int64}` mapping component name to matrix
+    index of bus; for `SystemPowerFlowContainer`, leaf values are Dict{Union{String, Int64},
+    Union{String, Int64}} mapping component name/bus number to component name/bus number.
+    """
+    input_key_map::Dict{Symbol, <:Dict{<:OptimizationContainerKey, <:Any}}
+    is_solved::Bool
+end
+
+function PowerFlowEvaluationData(power_flow_data::T) where {T <: PFS.PowerFlowContainer}
+    return PowerFlowEvaluationData{T}(
+        power_flow_data,
+        Dict{Symbol, Dict{OptimizationContainerKey, <:Any}}(),
+        false,
+    )
+end
+
+get_power_flow_data(ped::PowerFlowEvaluationData) = ped.power_flow_data
+get_input_key_map(ped::PowerFlowEvaluationData) = ped.input_key_map

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)