Add attribute to select approximation method

Project.toml

@@ -29,7 +29,7 @@ PowerFlowsExt = "PowerFlows"
 [sources]
 InfrastructureSystems = {url = "https://github.com/NREL-Sienna/InfrastructureSystems.jl", rev = "IS4"}
 PowerSystems = {url = "https://github.com/NREL-Sienna/PowerSystems.jl", rev = "psy6"}
-InfrastructureOptimizationModels = {url = "https://github.com/NREL-Sienna/InfrastructureOptimizationModels.jl", rev = "lk/pom-test-fixes"}
+InfrastructureOptimizationModels = {url = "https://github.com/Sienna-Platform/InfrastructureOptimizationModels.jl", rev = "ac/tolerance-option"}
 
 [compat]
 Dates = "1"

src/PowerOperationsModels.jl

@@ -537,6 +537,7 @@ export HydroWaterFactorModel
 export HydroWaterModelReservoir
 export HydroTurbineBilinearDispatch
 export HydroTurbineWaterLinearDispatch
+export HydroTurbineMILPBilinearDispatch
 export HydroTurbineWaterLinearCommitment
 export HydroEnergyModelReservoir
 export HydroTurbineEnergyDispatch

src/core/formulations.jl

@@ -326,6 +326,44 @@ Formulation type to add injection variables for a HydroTurbine connected to rese
 """
 struct HydroTurbineBilinearDispatch <: AbstractHydroDispatchFormulation end
 
+"""
+MILP formulation for the turbined-flow × head bilinear product in the hydro
+turbine power-output constraint. Adds injection variables for a HydroTurbine
+connected to reservoirs using a linearized approximation of the bilinear model.
+
+Selects the bilinear approximation scheme via `DeviceModel` attributes (so users
+do not need to depend on `InfrastructureOptimizationModels`). Users pick the
+desired approximation gap via `bilinear_tolerance`; the constraint constructor
+combines that with the per-device flow and head bounds to size each method's
+discretization automatically (no manual segment count).
+
+# Attributes
+- `bilinear_approximation::String` (default `"bin2"`): top-level scheme.
+  Supported: `"bin2"`, `"hybs"`, `"nmdt"`, `"dnmdt"`, `"none"`.
+- `bilinear_quadratic_method::String` (default `"solver_sos2"`): inner quadratic
+  PWL method used when `bilinear_approximation ∈ {"bin2","hybs"}`. Supported:
+  `"solver_sos2"`, `"manual_sos2"`, `"sawtooth"`, `"epigraph"`, `"nmdt"`,
+  `"dnmdt"`, `"none"`.
+- `bilinear_tolerance::Float64` (default `1e-2`): maximum approximation gap
+  for the chosen method. Each IOM tolerance constructor uses it together with
+  the per-device `max_delta` (computed at the call site from variable bounds)
+  to pick a depth via `ceil`.
+- `bilinear_add_mccormick::Union{Bool,Nothing}` (default `nothing`): when
+  `nothing`, defers to the IOM struct default (`Bin2Config` → `true`,
+  `HybSConfig` → `false`). Ignored by `nmdt`/`dnmdt`/`none`.
+- `bilinear_epigraph_depth::Union{Int,Nothing}` (default `nothing`): when
+  `nothing`, defers to the IOM struct default (`NMDTQuadConfig` /
+  `DNMDTQuadConfig` → `3*depth`). `"hybs"` has no IOM default for this field
+  and *must* override.
+
+# Sentinel convention
+`nothing` means "use the IOM constructor's default value." POM does not
+duplicate IOM struct defaults; it just passes through the user's overrides.
+
+See: [`PowerSystems.HydroGen`](@extref).
+"""
+struct HydroTurbineMILPBilinearDispatch <: 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.
@@ -364,6 +402,7 @@ These types share constructors.
 """
 const HydroTurbineWaterFormulation = Union{
     HydroTurbineBilinearDispatch,
+    HydroTurbineMILPBilinearDispatch,
     HydroTurbineWaterLinearDispatch,
     HydroTurbineWaterLinearCommitment,
 }

src/static_injector_models/hydro_generation.jl

@@ -421,6 +421,42 @@ function get_default_attributes(
     return Dict{String, Any}("head_fraction_usage" => 0.0)
 end
 
+"""
+Default `DeviceModel` attributes for `HydroTurbineMILPBilinearDispatch`. The
+returned dictionary picks the bilinear approximation scheme used inside the
+turbine-power constraint; see [`HydroTurbineMILPBilinearDispatch`](@ref) for the
+full attribute reference and the `nothing` sentinel convention.
+
+The default tolerance is `1e-2` paired with `Bin2` + `SolverSOS2`. The actual
+segment count is computed per device at constraint-build time from the
+per-device flow and head ranges, so two systems with very different bounds will
+get appropriately different discretizations from the same tolerance setting.
+"""
+function get_default_attributes(
+    ::Type{T},
+    ::Type{D},
+) where {T <: PSY.HydroTurbine, D <: HydroTurbineMILPBilinearDispatch}
+    return Dict{String, Any}(
+        # Top-level bilinear approximation scheme.
+        # Supported: "bin2", "hybs", "nmdt", "dnmdt", "none".
+        "bilinear_approximation" => "bin2",
+        # Inner quadratic PWL method (used when bilinear_approximation ∈ {"bin2","hybs"}).
+        # Supported: "solver_sos2", "manual_sos2", "sawtooth", "epigraph",
+        # "nmdt", "dnmdt", "none".
+        "bilinear_quadratic_method" => "solver_sos2",
+        # Maximum approximation gap. Combined with the per-device max_delta
+        # to size the discretization automatically.
+        "bilinear_tolerance" => 1e-2,
+        # `nothing` ⇒ defer to IOM struct default
+        # (Bin2Config: true; HybSConfig: false; ignored otherwise).
+        "bilinear_add_mccormick" => nothing,
+        # `nothing` ⇒ defer to IOM struct default
+        # (HybSConfig has no default and must be overridden; NMDT/DNMDT
+        #  quadratic: 3*depth; ignored otherwise).
+        "bilinear_epigraph_depth" => nothing,
+    )
+end
+
 function get_default_attributes(
     ::Type{T},
     ::Type{D},
@@ -1817,6 +1853,214 @@ function add_constraints!(
     return
 end
 
+"""
+Build an `IOM.QuadraticApproxConfig` from attribute values.
+
+`tolerance` is the requested approximation gap (`bilinear_tolerance`).
+`max_delta` is the domain extent the quadratic config will be applied to
+(computed at the call site from per-device variable bounds).
+`epi` is either an `Int` (`bilinear_epigraph_depth` override) or `nothing`,
+in which case IOM's struct default applies — relevant only for `NMDTQuadConfig`
+and `DNMDTQuadConfig`, which fall back to `3*depth` when no override is given.
+
+Errors with a list of supported method strings when `method` is unrecognized.
+"""
+function _build_quadratic_config(
+    method::String,
+    tolerance::Float64,
+    max_delta::Float64,
+    epi,
+)
+    if method == "solver_sos2"
+        return IOM.SolverSOS2QuadConfig(; tolerance, max_delta)
+    elseif method == "manual_sos2"
+        return IOM.ManualSOS2QuadConfig(; tolerance, max_delta)
+    elseif method == "sawtooth"
+        return IOM.SawtoothQuadConfig(; tolerance, max_delta)
+    elseif method == "epigraph"
+        return IOM.EpigraphQuadConfig(; tolerance, max_delta)
+    elseif method == "nmdt"
+        return if epi === nothing
+            IOM.NMDTQuadConfig(; tolerance, max_delta)
+        else
+            IOM.NMDTQuadConfig(; tolerance, max_delta, epigraph_depth = epi)
+        end
+    elseif method == "dnmdt"
+        return if epi === nothing
+            IOM.DNMDTQuadConfig(; tolerance, max_delta)
+        else
+            IOM.DNMDTQuadConfig(; tolerance, max_delta, epigraph_depth = epi)
+        end
+    elseif method == "none"
+        return IOM.NoQuadApproxConfig()
+    else
+        error(
+            "Unsupported bilinear_quadratic_method \"$(method)\". " *
+            "Supported: \"solver_sos2\", \"manual_sos2\", \"sawtooth\", " *
+            "\"epigraph\", \"nmdt\", \"dnmdt\", \"none\".",
+        )
+    end
+end
+
+"""
+Build an `IOM.BilinearApproxConfig` from the `DeviceModel`'s attributes and the
+per-device flow / head domain widths.
+
+Reads `bilinear_approximation`, `bilinear_quadratic_method`,
+`bilinear_tolerance`, `bilinear_add_mccormick`, and `bilinear_epigraph_depth`
+from `model`. Sentinel (`nothing`) attribute values mean "let the IOM
+constructor's default apply" — POM never duplicates an IOM struct default.
+
+For `bin2` / `hybs`, the inner quadratic config is applied to x², y², and
+(x±y)². The widest of those is (x±y), with delta = `max_delta_x + max_delta_y`,
+so that value is passed as the inner quadratic's `max_delta`. This guarantees
+the requested tolerance for all three (the bound on x² and y² will be
+correspondingly tighter, which is fine).
+
+Errors when:
+- `bilinear_approximation` is not one of `"bin2"`, `"hybs"`, `"nmdt"`,
+  `"dnmdt"`, `"none"`.
+- `bilinear_approximation == "hybs"` but `bilinear_epigraph_depth === nothing`
+  (`HybSConfig` has no IOM-side default for `epigraph_depth`).
+
+See [`HydroTurbineMILPBilinearDispatch`](@ref) for the attribute reference.
+"""
+function _build_bilinear_config(
+    model::DeviceModel{<:PSY.HydroTurbine, HydroTurbineMILPBilinearDispatch},
+    max_delta_x::Float64,
+    max_delta_y::Float64,
+)
+    method = get_attribute(model, "bilinear_approximation")
+    tolerance = get_attribute(model, "bilinear_tolerance")
+    add_mc = get_attribute(model, "bilinear_add_mccormick")
+    epi = get_attribute(model, "bilinear_epigraph_depth")
+
+    if method == "bin2"
+        inner_delta = max_delta_x + max_delta_y
+        quad = _build_quadratic_config(
+            get_attribute(model, "bilinear_quadratic_method"),
+            tolerance,
+            inner_delta,
+            epi,
+        )
+        return if add_mc === nothing
+            IOM.Bin2Config(quad)
+        else
+            IOM.Bin2Config(quad; add_mccormick = add_mc)
+        end
+    elseif method == "hybs"
+        epi === nothing && error(
+            "bilinear_approximation = \"hybs\" requires a non-nothing " *
+            "bilinear_epigraph_depth attribute (IOM.HybSConfig has no default).",
+        )
+        inner_delta = max_delta_x + max_delta_y
+        quad = _build_quadratic_config(
+            get_attribute(model, "bilinear_quadratic_method"),
+            tolerance,
+            inner_delta,
+            epi,
+        )
+        return if add_mc === nothing
+            IOM.HybSConfig(quad; epigraph_depth = epi)
+        else
+            IOM.HybSConfig(quad; epigraph_depth = epi, add_mccormick = add_mc)
+        end
+    elseif method == "nmdt"
+        return IOM.NMDTBilinearConfig(; tolerance, max_delta_x, max_delta_y)
+    elseif method == "dnmdt"
+        return IOM.DNMDTBilinearConfig(; tolerance, max_delta_x, max_delta_y)
+    elseif method == "none"
+        return IOM.NoBilinearApproxConfig()
+    else
+        error(
+            "Unsupported bilinear_approximation \"$(method)\" for " *
+            "HydroTurbineMILPBilinearDispatch. " *
+            "Supported: \"bin2\", \"hybs\", \"nmdt\", \"dnmdt\", \"none\".",
+        )
+    end
+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 <: HydroTurbineMILPBilinearDispatch,
+    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)
+
+        flow_lb = get_variable_lower_bound(HydroTurbineFlowRateVariable, d, W)
+        flow_ub = get_variable_upper_bound(HydroTurbineFlowRateVariable, d, W)
+        flow_delta = flow_ub - flow_lb
+
+        head_bounds = [
+            (
+                min = get_variable_lower_bound(HydroReservoirHeadVariable, res, W),
+                max = get_variable_upper_bound(HydroReservoirHeadVariable, res, W),
+            ) for res in reservoirs
+        ]
+        # Worst-case head range across the turbine's reservoirs — gives a
+        # single config that meets the requested tolerance for every pair.
+        head_delta = maximum(b.max - b.min for b in head_bounds)
+
+        flow_bounds = repeat([(min = flow_lb, max = flow_ub)], length(reservoirs))
+
+        fh_prod = IOM._add_bilinear_approx!(
+            _build_bilinear_config(model, flow_delta, head_delta),
+            container,
+            V,
+            PSY.get_name.(reservoirs),
+            time_steps,
+            flow[name, :, :],
+            head,
+            flow_bounds,
+            head_bounds,
+            "$(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,11 @@ _maybe_add_on_variables!(
 _maybe_add_on_variables!(
     ::OptimizationContainer,
     devices,
-    ::Union{Type{HydroTurbineBilinearDispatch}, Type{HydroTurbineWaterLinearDispatch}},
+    ::Union{
+        Type{HydroTurbineBilinearDispatch},
+        Type{HydroTurbineMILPBilinearDispatch},
+        Type{HydroTurbineWaterLinearDispatch},
+    },
 ) = nothing
 
 function construct_device!(

test/runtests.jl

@@ -22,11 +22,15 @@ const DISABLED_TEST_FILES = [  # Can generate with ls -1 test | grep "test_.*.jl
 # "test_device_synchronous_condenser_constructors.jl",
 # "test_device_thermal_generation_constructors.jl",
 # "test_formulation_combinations.jl",
+# "test_import_export_cost.jl",
 # "test_initialization_problem.jl",
+# "test_is_time_variant_proportional.jl",
+# "test_market_bid_cost.jl",
+# "test_mbc_parameter_population.jl",
 # "test_model_decision.jl",
 # "test_multi_interval.jl",
-# "test_network_constructors_with_dlr.jl",
 # "test_network_constructors.jl",
+# "test_network_constructors_with_dlr.jl",
 # "test_problem_template.jl",
 # "test_storage_device_models.jl",
 # "test_transfer_initial_conditions.jl",