Add hydro formulation with bilinear approximation

Project.toml

@@ -28,7 +28,6 @@ PowerFlowsExt = "PowerFlows"
 
 [sources]
 InfrastructureSystems = {url = "https://github.com/NREL-Sienna/InfrastructureSystems.jl", rev = "IS4"}
-InfrastructureOptimizationModels = {url = "https://github.com/NREL-Sienna/InfrastructureOptimizationModels.jl", rev = "main"}
 
 [compat]
 Dates = "1"

src/PowerOperationsModels.jl

@@ -141,8 +141,8 @@ import InfrastructureOptimizationModels:
     _include_constant_min_gen_power_in_constraint,
     add_variable_cost_to_objective!,
     _vom_offer_direction,
-    _add_pwl_constraint!,
-    add_pwl_term!,
+    add_pwl_constraint_delta!,
+    add_pwl_term_delta!,
     get_output_offer_curves,
     # Internal utilities used by market bid overrides and proportional_cost
     is_time_variant,
@@ -175,7 +175,7 @@ import InfrastructureOptimizationModels:
     set_service_model!,
     finalize_template!,
     make_empty_jump_model_with_settings,
-    _set_model!
+    set_model!
 
 using InfrastructureOptimizationModels # TODO: use explicit imports.
 
@@ -527,6 +527,7 @@ export HydroWaterFactorModel
 export HydroWaterModelReservoir
 export HydroTurbineBilinearDispatch
 export HydroTurbineWaterLinearDispatch
+export HydroTurbineBin2BilinearDispatch
 export HydroEnergyModelReservoir
 export HydroTurbineEnergyDispatch
 export HydroTurbineEnergyCommitment

src/common_models/market_bid_overrides.jl

@@ -175,7 +175,7 @@ function add_variable_cost_to_objective!(
     ::ImportExportSourceModel,
 )
     isnothing(get_output_offer_curves(cost_function)) && return
-    add_pwl_term!(
+    add_pwl_term_delta!(
         IncrementalOffer(),
         container,
         component,
@@ -194,7 +194,7 @@ function add_variable_cost_to_objective!(
     ::ImportExportSourceModel,
 )
     isnothing(get_input_offer_curves(cost_function)) && return
-    add_pwl_term!(
+    add_pwl_term_delta!(
         DecrementalOffer(),
         container,
         component,
@@ -221,7 +221,7 @@ function add_variable_cost_to_objective!(
     if !(isnothing(get_output_offer_curves(cost_function)))
         error("Component $(component_name) is not allowed to participate as a supply.")
     end
-    add_pwl_term!(
+    add_pwl_term_delta!(
         DecrementalOffer(),
         container,
         component,
@@ -241,7 +241,7 @@ _vom_offer_direction(::AbstractControllablePowerLoadFormulation) = DecrementalOf
 """
 PWL block offer constraints for ORDC (ReserveDemandCurve).
 """
-function _add_pwl_constraint!(
+function add_pwl_constraint_delta!(
     container::OptimizationContainer,
     component::T,
     ::U,
@@ -273,7 +273,7 @@ end
 """
 PWL cost terms for StepwiseCostReserve (AbstractServiceFormulation).
 """
-function add_pwl_term!(
+function add_pwl_term_delta!(
     container::OptimizationContainer,
     component::T,
     cost_data::PSY.CostCurve{PSY.PiecewiseIncrementalCurve},
@@ -300,7 +300,7 @@ function add_pwl_term!(
     slopes = IS.get_y_coords(data)
     break_points = PSY.get_x_coords(data)
     for t in time_steps
-        pwl_vars = add_pwl_variables!(
+        pwl_vars = add_pwl_variables_delta!(
             container,
             PiecewiseLinearBlockIncrementalOffer,
             T,
@@ -309,9 +309,9 @@ function add_pwl_term!(
             length(slopes);
             upper_bound = Inf,
         )
-        add_pwl_constraint!(container, component, U(), break_points, pwl_vars, t)
+        add_pwl_constraint_delta!(container, component, U(), break_points, pwl_vars, t)
         pwl_cost_expressions[t] =
-            get_pwl_cost_expression(pwl_vars, slopes, multiplier * dt)
+            get_pwl_cost_expression_delta(pwl_vars, slopes, multiplier * dt)
     end
     return pwl_cost_expressions
 end

src/core/formulations.jl

@@ -321,6 +321,11 @@ Formulation type to add injection variables for a HydroTurbine connected to rese
 """
 struct HydroTurbineBilinearDispatch <: AbstractHydroDispatchFormulation end
 
+"""
+Formulation type to add injection variables for a HydroTurbine connected to reservoirs using a bilinear model (with water flow variables) [`PowerSystems.HydroGen`](@extref). Uses a linearized approximation.
+"""
+struct HydroTurbineBin2BilinearDispatch <: AbstractHydroDispatchFormulation end
+
 """
 Formulation type to add injection variables for a HydroTurbine connected to reservoirs using a linear model [`PowerSystems.HydroGen`](@extref).
 The model assumes a shallow reservoir. The head for the conversion between water flow and power can be approximated as a linear function of the water flow on which the head elevation is always the intake elevation.

src/core/problem_template.jl

@@ -140,7 +140,7 @@ function set_device_model!(
     template::OperationsProblemTemplate,
     model::DeviceModel{D},
 ) where {D <: IS.InfrastructureSystemsComponent}
-    _set_model!(template.devices, model)
+    set_model!(template.devices, model)
     return
 end
 
@@ -152,7 +152,7 @@ function set_device_model!(
     template::OperationsProblemTemplate,
     model::DeviceModel{D},
 ) where {D <: PSY.Branch}
-    _set_model!(template.branches, model)
+    set_model!(template.branches, model)
     return
 end
 
@@ -191,15 +191,15 @@ function set_service_model!(
     service_name::String,
     model::ServiceModel{T, <:AbstractServiceFormulation},
 ) where {T <: PSY.Service}
-    _set_model!(template.services, (service_name, Symbol(T)), model)
+    set_model!(template.services, (service_name, Symbol(T)), model)
     return
 end
 
 function set_service_model!(
     template::OperationsProblemTemplate,
     model::ServiceModel{<:PSY.Service, <:AbstractServiceFormulation},
 )
-    _set_model!(template.services, model)
+    set_model!(template.services, model)
     return
 end
 

src/energy_storage_models/storage_models.jl

@@ -1916,7 +1916,7 @@ function _add_variable_cost_to_objective!(
     @debug "Market Bid" _group = LOG_GROUP_COST_FUNCTIONS component_name
     incremental_cost_curves = PSY.get_incremental_offer_curves(cost_function)
     if !isnothing(incremental_cost_curves)
-        add_pwl_term!(
+        add_pwl_term_delta!(
             IncrementalOffer(),
             container,
             component,
@@ -1942,7 +1942,7 @@ function _add_variable_cost_to_objective!(
     @debug "Market Bid" _group = LOG_GROUP_COST_FUNCTIONS component_name
     decremental_cost_curves = PSY.get_decremental_offer_curves(cost_function)
     if !isnothing(decremental_cost_curves)
-        add_pwl_term!(
+        add_pwl_term_delta!(
             DecrementalOffer(),
             container,
             component,

src/network_models/instantiate_network_model.jl

@@ -106,6 +106,29 @@ function _get_irreducible_buses_due_to_dlrs(
     return collect(irreducible_buses)
 end
 
+function _get_unmodeled_branch_types(
+    branch_models::BranchModelContainer,
+    sys::PSY.System,
+)
+    unmodeled = DataType[]
+    for d in PSY.get_existing_device_types(sys)
+        if d <: PSY.ACTransmission && !haskey(branch_models, Symbol(d))
+            push!(unmodeled, d)
+        end
+    end
+    return unmodeled
+end
+
+function _validate_network_and_branches(
+    model::NetworkModel,
+    branch_models::BranchModelContainer,
+    sys::PSY.System,
+)
+    unmodeled = _get_unmodeled_branch_types(branch_models, sys)
+    IOM._check_branch_network_compatibility(model, unmodeled)
+    return
+end
+
 #################################################################################
 # Generic fallback for AbstractPowerModel (Ybus-based models: ACP, ACR, etc.)
 #################################################################################
@@ -116,6 +139,7 @@ function IOM.instantiate_network_model!(
     number_of_steps::Int,
     sys::PSY.System,
 ) where {T <: AbstractPowerModel}
+    _validate_network_and_branches(model, branch_models, sys)
     if isempty(model.subnetworks)
         model.subnetworks = PNM.find_subnetworks(sys)
     end
@@ -181,6 +205,7 @@ function IOM.instantiate_network_model!(
     number_of_steps::Int,
     sys::PSY.System,
 )
+    _validate_network_and_branches(model, branch_models, sys)
     PNM.populate_branch_maps_by_type!(model.network_reduction)
     empty!(model.reduced_branch_tracker)
     IOM.set_number_of_steps!(model.reduced_branch_tracker, number_of_steps)
@@ -197,6 +222,7 @@ function IOM.instantiate_network_model!(
     number_of_steps::Int,
     sys::PSY.System,
 )
+    _validate_network_and_branches(model, branch_models, sys)
     if isempty(model.subnetworks)
         model.subnetworks = PNM.find_subnetworks(sys)
     end
@@ -225,7 +251,7 @@ function IOM.instantiate_network_model!(
         model,
         branch_models,
     )
-    IOM._check_branch_network_compatibility(model, branch_models, sys)
+    _validate_network_and_branches(model, branch_models, sys)
     if IOM.get_PTDF_matrix(model) === nothing || !isempty(irreducible_buses)
         if IOM.get_PTDF_matrix(model) !== nothing
             @warn "Provided PTDF Matrix is being ignored since irreducible buses were identified because of DLRs. Recalculating PTDF Matrix with PowerNetworkMatrices.PTDF and the identified irreducible buses."
@@ -334,7 +360,7 @@ function IOM.instantiate_network_model!(
         model,
         branch_models,
     )
-    IOM._check_branch_network_compatibility(model, branch_models, sys)
+    _validate_network_and_branches(model, branch_models, sys)
     if IOM.get_PTDF_matrix(model) === nothing || !isempty(irreducible_buses)
         if IOM.get_PTDF_matrix(model) !== nothing
             @warn "Provided PTDF Matrix is being ignored since irreducible buses were identified because of DLRs. Recalculating PTDF Matrix with PowerNetworkMatrices.PTDF and the identified irreducible buses."

src/services_models/reserves.jl

@@ -543,7 +543,7 @@ function _add_reserves_variable_cost_to_objective!(
     end
 
     pwl_cost_expressions =
-        add_pwl_term!(container, component, variable_cost, T(), U())
+        add_pwl_term_delta!(container, component, variable_cost, T(), U())
     for t in time_steps
         add_to_expression!(
             container,

src/static_injector_models/hydro_generation.jl

@@ -488,7 +488,7 @@ function add_variables!(
     T <: HydroTurbineFlowRateVariable,
     U <: Union{Vector{D}, IS.FlattenIteratorWrapper{D}},
     W <: Union{Vector{E}, IS.FlattenIteratorWrapper{E}},
-    X <: Union{HydroTurbineBilinearDispatch, HydroTurbineWaterLinearDispatch},
+    X <: Union{HydroTurbineBilinearDispatch, HydroTurbineWaterLinearDispatch, HydroTurbineBin2BilinearDispatch},
 } where {
     D <: PSY.HydroTurbine,
     E <: PSY.HydroReservoir,
@@ -1835,6 +1835,80 @@ function add_constraints!(
     return
 end
 
+"""
+This function define the relationship between turbined flow and power produced with a linear approximation for the bilinear product.
+"""
+function add_constraints!(
+    container::OptimizationContainer,
+    sys::PSY.System,
+    ::Type{TurbinePowerOutputConstraint},
+    devices::IS.FlattenIteratorWrapper{V},
+    model::DeviceModel{V, W},
+    ::NetworkModel{X},
+) where {
+    V <: PSY.HydroTurbine,
+    W <: HydroTurbineBin2BilinearDispatch,
+    X <: PM.AbstractPowerModel,
+}
+    time_steps = get_time_steps(container)
+    base_power = get_model_base_power(container)
+    names = PSY.get_name.(devices)
+    constraint =
+        add_constraints_container!(
+            container,
+            TurbinePowerOutputConstraint(),
+            V,
+            names,
+            time_steps,
+        )
+    power = get_variable(container, ActivePowerVariable(), V)
+    flow = get_variable(container, HydroTurbineFlowRateVariable(), V)
+    head = get_variable(container, HydroReservoirHeadVariable(), PSY.HydroReservoir)
+    for d in devices
+        name = PSY.get_name(d)
+        conversion_factor = PSY.get_conversion_factor(d)
+        reservoirs = filter(PSY.get_available, PSY.get_connected_head_reservoirs(sys, d))
+        powerhouse_elevation = PSY.get_powerhouse_elevation(d)
+
+        fh_prod = IOM._add_bilinear_approx!(
+            IOM.Bin2Config(IOM.SolverSOS2QuadConfig(4)),
+            container,
+            V,
+            PSY.get_name.(reservoirs),
+            time_steps,
+            flow[name, :, :],
+            head,
+            [
+                (
+                    min = get_variable_lower_bound(HydroTurbineFlowRateVariable(), d, W()),
+                    max = get_variable_upper_bound(HydroTurbineFlowRateVariable(), d, W())
+                ) for _=1:length(reservoirs)
+            ],
+            [
+                (
+                    min = get_variable_lower_bound(HydroReservoirHeadVariable(), res, W()),
+                    max = get_variable_upper_bound(HydroReservoirHeadVariable(), res, W())
+                ) for res in reservoirs
+            ],
+            "$(get_name(d))_FlowHeadProduct"
+        )
+
+        for t in time_steps
+            constraint[name, t] = JuMP.@constraint(
+                container.JuMPmodel,
+                power[name, t] ==
+                GRAVITATIONAL_CONSTANT * WATER_DENSITY * conversion_factor *
+                sum(
+                    fh_prod[PSY.get_name(res), t]
+                    - powerhouse_elevation * flow[name, PSY.get_name(res), t]
+                    for res in reservoirs
+                ) / (1e6 * base_power)
+            )
+        end
+    end
+    return
+end
+
 ############################################################################
 ############################### Expressions ################################
 ############################################################################

src/static_injector_models/hydrogeneration_constructor.jl

@@ -1809,7 +1809,7 @@ function construct_device!(
     network_model::NetworkModel{S},
 ) where {
     H <: PSY.HydroTurbine,
-    D <: Union{HydroTurbineBilinearDispatch, HydroTurbineWaterLinearDispatch},
+    D <: Union{HydroTurbineBilinearDispatch, HydroTurbineWaterLinearDispatch, HydroTurbineBin2BilinearDispatch},
     S <: AbstractActivePowerModel,
 }
     devices = get_available_components(model, sys)
@@ -1873,7 +1873,7 @@ function construct_device!(
     network_model::NetworkModel{S},
 ) where {
     H <: PSY.HydroTurbine,
-    D <: Union{HydroTurbineBilinearDispatch, HydroTurbineWaterLinearDispatch},
+    D <: Union{HydroTurbineBilinearDispatch, HydroTurbineWaterLinearDispatch, HydroTurbineBin2BilinearDispatch},
     S <: AbstractActivePowerModel,
 }
     devices = get_available_components(model, sys)

src/static_injector_models/thermal_generation.jl

@@ -1559,7 +1559,7 @@ Add PWL cost terms for ThermalDispatchNoMin formulation.
 Rejects non-convex or negative-slope PWL data since ThermalDispatchNoMin
 cannot use SOS-2 formulations.
 """
-function IOM.add_pwl_term!(
+function IOM.add_pwl_term_lambda!(
     container::IOM.OptimizationContainer,
     component::T,
     cost_function::Union{
@@ -1621,7 +1621,7 @@ function IOM.add_pwl_term!(
     temp_cost_function =
         IOM.create_temporary_cost_function_in_system_per_unit(cost_function, data)
     for t in time_steps
-        IOM.add_pwl_variables!(container, T, name, t, data)
+        IOM.add_pwl_variables_lambda!(container, T, name, t, data)
         power_var = IOM.get_variable(container, U(), T)[name, t]
         IOM._add_pwl_constraint_standard!(
             container,
@@ -1632,7 +1632,7 @@ function IOM.add_pwl_term!(
             power_var,
         )
         pwl_cost =
-            IOM.get_pwl_cost_expression(
+            IOM.get_pwl_cost_expression_lambda(
                 container,
                 component,
                 t,

test/Project.toml

@@ -32,7 +32,6 @@ TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
 
 [sources]
-InfrastructureOptimizationModels = {rev = "main", url = "https://github.com/NREL-Sienna/InfrastructureOptimizationModels.jl"}
 InfrastructureSystems = {rev = "IS4", url = "https://github.com/NREL-Sienna/InfrastructureSystems.jl"}
 
 [compat]