[DNMY] Support in the loop PowerFlow evaluation

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"
@@ -47,6 +48,7 @@ Logging = "1"
 MathOptInterface = "1"
 PowerModels = "^0.20"
 PowerNetworkMatrices = "^0.10"
+PowerFlows = "0.6"
 PowerSystems = "^3"
 PrettyTables = "2"
 ProgressMeter = "^1.5"

src/PowerSimulations.jl

@@ -238,6 +238,10 @@ export LowerBoundFeedForwardSlack
 export TimeDurationOn
 export TimeDurationOff
 export PowerOutput
+export PowerFlowAngle
+export PowerFlowVoltage
+export PowerFlowLineReactivePower
+export PowerFlowLineActivePower
 
 # Constraints
 export AbsoluteValueConstraint
@@ -339,6 +343,7 @@ import LinearAlgebra
 import JSON3
 import PowerSystems
 import InfrastructureSystems
+import PowerFlows
 import PowerNetworkMatrices
 import PowerNetworkMatrices: PTDF, VirtualPTDF
 export PTDF
@@ -399,6 +404,7 @@ const MOI = MathOptInterface
 const MOIU = MathOptInterface.Utilities
 const MOPFM = MOI.FileFormats.Model
 const PNM = PowerNetworkMatrices
+const PFS = PowerFlows
 const TS = TimeSeries
 
 ################################################################################
@@ -530,6 +536,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/powerflow_evaluation.jl")
 
 include("initial_conditions/initialization.jl")
 

src/core/auxiliary_variables.jl

@@ -31,7 +31,29 @@ Auxiliary Variable for Thermal Generation Models that solve for power above min
 """
 struct PowerOutput <: AuxVariableType end
 
+"""
+Auxiliary Variable for the bus angle results from power flow evaluation
+"""
+struct PowerFlowAngle <: AuxVariableType end
+
+"""
+Auxiliary Variable for the bus voltage magnitued results from power flow evaluation
+"""
+struct PowerFlowVoltage <: AuxVariableType end
+
+"""
+Auxiliary Variable for the line reactive flow from power flow evaluation
+"""
+struct PowerFlowLineReactivePower <: AuxVariableType end
+
+"""
+Auxiliary Variable for the line active flow from power flow evaluation
+"""
+struct PowerFlowLineActivePower <: AuxVariableType end
+
 should_write_resulting_value(::Type{<:AuxVariableType}) = true
 
 convert_result_to_natural_units(::Type{<:AuxVariableType}) = false
 convert_result_to_natural_units(::Type{PowerOutput}) = true
+convert_result_to_natural_units(::Type{PowerFlowLineReactivePower}) = true
+convert_result_to_natural_units(::Type{PowerFlowLineActivePower}) = true

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
+    powerflow_evaluation::Union{Nothing, PFS.PowerFlowEvaluationModel}
 
     function NetworkModel(
         ::Type{T};
@@ -42,6 +43,7 @@ mutable struct NetworkModel{T <: PM.AbstractPowerModel}
         reduce_radial_branches = false,
         subnetworks = Dict{Int, Set{Int}}(),
         duals = Vector{DataType}(),
+        powerflow_evaluation = nothing,
     ) where {T <: PM.AbstractPowerModel}
         _check_pm_formulation(T)
         new{T}(
@@ -52,6 +54,7 @@ mutable struct NetworkModel{T <: PM.AbstractPowerModel}
             duals,
             PNM.RadialNetworkReduction(),
             reduce_radial_branches,
+            powerflow_evaluation,
         )
     end
 end
@@ -66,6 +69,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_powerflow_evaluation(m::NetworkModel) = m.powerflow_evaluation
 has_subnetworks(m::NetworkModel) = !isempty(m.bus_area_map)
 
 function add_dual!(model::NetworkModel, dual)

src/core/optimization_container.jl

@@ -110,6 +110,7 @@ mutable struct OptimizationContainer <: AbstractModelContainer
     built_for_recurrent_solves::Bool
     metadata::OptimizationContainerMetadata
     default_time_series_type::Type{<:PSY.TimeSeriesData}
+    power_flow_data::Union{PFS.PowerFlowData, Nothing}
 end
 
 function OptimizationContainer(
@@ -154,6 +155,7 @@ function OptimizationContainer(
         false,
         OptimizationContainerMetadata(),
         T,
+        nothing,
     )
 end
 
@@ -650,8 +652,8 @@ 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_powerflow_evaluation(transmission_model), sys)
     check_optimization_container(container)
-
     return
 end
 

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.powerflow_evaluation = nothing
     bus_area_map = get_bus_area_map(get_network_model(model.template))
     if !isempty(bus_area_map)
         network_model.bus_area_map = get_bus_area_map(get_network_model(model.template))

src/network_models/powerflow_evaluation.jl

@@ -0,0 +1,72 @@
+function add_power_flow_data!(::OptimizationContainer, ::Nothing, ::PSY.System)
+    # NO OP function
+end
+
+function _add_branches_aux_variables!(
+    container::OptimizationContainer,
+    vars::Vector{DataType},
+    branch_types::Vector{DataType},
+    branch_names::Vector{String},
+)
+    time_steps = get_time_steps(container)
+    for var_type in vars
+        for D in Set(branch_types)
+            add_aux_variable_container!(
+                container,
+                var_type,
+                D,
+                branch_names[branch_types == D],
+                time_steps,
+            )
+        end
+    end
+    return
+end
+
+function _add_buses_aux_variables!(
+    container::OptimizationContainer,
+    vars::Vector{DataType},
+    bus_names::Vector{String},
+)
+    time_steps = get_time_steps(container)
+    for var_type in vars
+        add_aux_variable_container!(
+            container,
+            var_type,
+            D,
+            bus_names,
+            time_steps,
+        )
+    end
+    return
+end
+
+function add_power_flow_data!(
+    container::OptimizationContainer,
+    evaluator::T,
+    sys::PSY.System,
+) where {T <: Union{PFS.PTDFDCPowerFlow, PFS.vPTDFDCPowerFlow}}
+    @info "Building PowerFlow evaluator using $(evaluator)"
+    container.power_flow_data =
+        PFS.PowerFlowData(evaluator, sys; time_steps = length(get_time_steps(container)))
+end
+
+function add_power_flow_data!(
+    container::OptimizationContainer,
+    evaluator::PFS.DCPowerFlow,
+    sys::PSY.System,
+)
+    @info "Building PowerFlow evaluator using $(evaluator)"
+    container.power_flow_data =
+        PFS.PowerFlowData(evaluator, sys; time_steps = length(get_time_steps(container)))
+end
+
+function add_power_flow_data!(
+    container::OptimizationContainer,
+    evaluator::PFS.ACPowerFlow,
+    sys::PSY.System,
+)
+    @info "Building PowerFlow evaluator using $(evaluator)"
+    container.power_flow_data =
+        PFS.PowerFlowData(evaluator, sys; time_steps = length(get_time_steps(container)))
+end

src/services_models/transmission_interface.jl

@@ -33,26 +33,26 @@ 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!(
         container,
         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] =

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

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"

test/test_network_constructors.jl

@@ -21,6 +21,7 @@ const NETWORKS_FOR_TESTING = [
     (PM.SparseSDPWRMPowerModel, scs_solver),
 ]
 
+
 @testset "All PowerModels models construction" begin
     c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
     for (network, solver) in NETWORKS_FOR_TESTING
@@ -652,7 +653,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;