Add tolerance-driven bilinear approximation config API

Project.toml

@@ -27,8 +27,7 @@ PowerFlows = "94fada2c-fd9a-4e89-8d82-81405f5cb4f6"
 [sources]
 InfrastructureSystems = {rev = "IS4", url = "https://github.com/Sienna-Platform/InfrastructureSystems.jl"}
 PowerSystems = {rev = "psy6", url = "https://github.com/Sienna-Platform/PowerSystems.jl"}
-# TODO: move IOM back to rev = "main" once ac/canonical-key-component-type merges.
-InfrastructureOptimizationModels = {rev = "ac/canonical-key-component-type", url = "https://github.com/Sienna-Platform/InfrastructureOptimizationModels.jl"}
+InfrastructureOptimizationModels = {rev = "main", url = "https://github.com/Sienna-Platform/InfrastructureOptimizationModels.jl"}
 PowerNetworkMatrices = {rev = "psy6", url = "https://github.com/Sienna-Platform/PowerNetworkMatrices.jl"}
 
 [extensions]
@@ -40,7 +39,7 @@ DocStringExtensions = "~0.8, ~0.9"
 InfrastructureSystems = "3"
 InteractiveUtils = "1.11.0"
 JuMP = "^1.28"
-PowerNetworkMatrices = "^0.19, ^0.20, ^0.22"
+PowerNetworkMatrices = "^0.22"
 PowerSystems = "5.3"
 PrettyTables = "3"
 ProgressMeter = "1.11.0"

src/PowerOperationsModels.jl

@@ -231,6 +231,7 @@ include("core/constraints.jl")
 include("core/auxiliary_variables.jl")
 include("core/parameters.jl")
 include("core/formulations.jl")
+include("core/bilinear_configs.jl")
 include("core/network_formulations.jl")
 include("core/problem_template.jl")
 include("core/feedforward_interface.jl")
@@ -567,7 +568,7 @@ export HydroWaterFactorModel
 export HydroWaterModelReservoir
 export HydroTurbineBilinearDispatch
 export HydroTurbineWaterLinearDispatch
-export HydroTurbineBin2BilinearDispatch
+export HydroTurbineMILPBilinearDispatch
 export HydroTurbineWaterLinearCommitment
 export HydroEnergyModelReservoir
 export HydroTurbineEnergyDispatch

src/ac_transmission_models/branch_constructor.jl

@@ -1675,15 +1675,6 @@ end
 ####################### Two-Terminal VSC HVDC Construct ####################
 ############################################################################
 
-# Quadratic / bilinear approximation traits — same scheme used by the MT
-# converter formulations.
-_quad_config(::Type{HVDCTwoTerminalVSCNLP}) = IOM.NoQuadApproxConfig()
-_quad_config(::Type{HVDCTwoTerminalVSCLP}) =
-    IOM.SolverSOS2QuadConfig(DEFAULT_INTERPOLATION_LENGTH)
-_bilinear_config(::Type{HVDCTwoTerminalVSCNLP}) = IOM.NoBilinearApproxConfig()
-_bilinear_config(::Type{HVDCTwoTerminalVSCLP}) =
-    IOM.Bin2Config(IOM.SolverSOS2QuadConfig(DEFAULT_INTERPOLATION_LENGTH))
-
 function construct_device!(
     container::OptimizationContainer,
     sys::PSY.System,
@@ -1740,31 +1731,26 @@ function construct_device!(
         (min = -_vsc_cable_i_max(d), max = _vsc_cable_i_max(d)) for d in devices
     ]
 
-    quad_cfg, bilin_cfg = _quad_config(F), _bilinear_config(F)
+    quad_cfg, bilin_cfg =
+        _build_converter_configs(F, device_model, vcat(v_f_bounds, v_t_bounds), i_bounds)
 
-    v_f_sq_expr = IOM._add_quadratic_approx!(
-        quad_cfg, container, PSY.TwoTerminalVSCLine,
-        line_names, time_steps, v_f_var, v_f_bounds, "v_f_sq",
-    )
-    v_t_sq_expr = IOM._add_quadratic_approx!(
-        quad_cfg, container, PSY.TwoTerminalVSCLine,
-        line_names, time_steps, v_t_var, v_t_bounds, "v_t_sq",
-    )
+    # The converter loss terms read `i_sq`; build it once and reuse it for both
+    # terminal bilinears.
     i_sq_expr = IOM._add_quadratic_approx!(
         quad_cfg, container, PSY.TwoTerminalVSCLine,
         line_names, time_steps, i_var, i_bounds, "i_sq",
     )
 
-    IOM._add_bilinear_approx!(
-        bilin_cfg, container, PSY.TwoTerminalVSCLine,
+    _add_converter_bilinear!(
+        bilin_cfg, quad_cfg, container, PSY.TwoTerminalVSCLine,
         line_names, time_steps,
-        v_f_sq_expr, i_sq_expr, v_f_var, i_var,
+        v_f_var, i_var, i_sq_expr,
         v_f_bounds, i_bounds, "vi_ft",
     )
-    IOM._add_bilinear_approx!(
-        bilin_cfg, container, PSY.TwoTerminalVSCLine,
+    _add_converter_bilinear!(
+        bilin_cfg, quad_cfg, container, PSY.TwoTerminalVSCLine,
         line_names, time_steps,
-        v_t_sq_expr, i_sq_expr, v_t_var, i_var,
+        v_t_var, i_var, i_sq_expr,
         v_t_bounds, i_bounds, "vi_tf",
     )
 

src/common_models/quadratic_converter_loss.jl

@@ -91,3 +91,145 @@ function _add_abs_value_constraints!(
     end
     return
 end
+
+#########################################
+####### Tolerance-driven configs ########
+#########################################
+#
+# The MILP/LP converter-loss formulations size their `v·I` (bilinear) and `I²`
+# (quadratic) surrogates from the same `DeviceModel` attributes used by
+# `HydroTurbineMILPBilinearDispatch` — `"bilinear_approximation"`,
+# `"bilinear_quadratic_method"`, and the `"bilinear_relative_tolerance"` /
+# `"bilinear_absolute_tolerance"` pair — bridged to IOM configs through
+# `_build_bilinear_config` (src/core/bilinear_configs.jl). A relative tolerance
+# is scaled to absolute by the term magnitude (`_resolve_tolerance`): the
+# product `max|v|·max|i|` for the bilinear, `max|i|²` for the standalone `I²`.
+# The NLP formulations keep both terms exact.
+#
+# The converters reuse the standalone `I²` we build for the loss term instead of
+# letting the bilinear recompute it. That works because the squares-based schemes
+# (`bin2`, `hybs`, `none`) accept a precomputed `(xsq, ysq)`, so we pass the
+# loss's `i_sq` straight through; the discretization-based schemes (`nmdt`,
+# `dnmdt`) never build `I²` at all, so there is nothing to duplicate and we use
+# the raw `(x_var, y_var)` form. `_add_converter_bilinear!` centralizes that
+# branch.
+
+# Worst-case domain width across devices, used to size the tolerance-driven
+# discretizations. Errors if the width is non-finite (missing/infinite limits).
+function _max_delta(bounds)
+    delta = maximum(b.max - b.min for b in bounds)
+    isfinite(delta) || error(
+        "Converter bilinear approximation requires finite variable bounds to " *
+        "size the discretization, but got a non-finite domain width ($(delta)). " *
+        "Check the device voltage/current limits.",
+    )
+    return delta
+end
+
+# Build (quad_cfg, bilin_cfg) for a converter-loss formulation from the device
+# `v_bounds`/`i_bounds`. MILP/LP read the attributes, derive the worst-case
+# domain widths and term magnitudes, and size the discretization; NLP keeps the
+# loss terms exact (no approximation).
+function _build_converter_configs(
+    ::Type{F},
+    model::DeviceModel,
+    v_bounds,
+    i_bounds,
+) where {F <: Union{QuadraticLossConverterMILP, HVDCTwoTerminalVSCLP}}
+    method = get_attribute(model, "bilinear_approximation")
+    quad_method = get_attribute(model, "bilinear_quadratic_method")
+    abs_tol = get_attribute(model, "bilinear_absolute_tolerance")
+    rel_tol = get_attribute(model, "bilinear_relative_tolerance")
+    method == "none" &&
+        return (IOM.NoQuadApproxConfig(), IOM.NoBilinearApproxConfig())
+    v_delta, i_delta = _max_delta(v_bounds), _max_delta(i_bounds)
+    # Bilinear v·i sized against the product magnitude max|v|·max|i|.
+    tol_prod = _resolve_tolerance(abs_tol, rel_tol, _max_abs(v_bounds) * _max_abs(i_bounds))
+    bilin_cfg = _build_bilinear_config(method, quad_method, tol_prod, v_delta, i_delta)
+    quad_cfg = _converter_quad_config(bilin_cfg, quad_method, abs_tol, rel_tol, i_bounds)
+    return (quad_cfg, bilin_cfg)
+end
+
+function _build_converter_configs(
+    ::Type{F},
+    ::DeviceModel,
+    _v_bounds,
+    _i_bounds,
+) where {F <: Union{QuadraticLossConverterNLP, HVDCTwoTerminalVSCNLP}}
+    return (IOM.NoQuadApproxConfig(), IOM.NoBilinearApproxConfig())
+end
+
+# Quad config for the standalone loss `I²`. For bin2/hybs the bilinear's inner
+# quad is reused — the bin2/hybs tolerance bound assumes the squares share that
+# inner quad (see bilinear_approximations/bin2.jl). For nmdt/dnmdt the bilinear
+# uses a discretization and never builds `I²`, so the loss `I²` is sized on its
+# own from the quad method and tolerance over the `I` domain (scaled by the `I²`
+# magnitude max|i|² when the tolerance is relative).
+_converter_quad_config(
+    bilin_cfg::Union{IOM.Bin2Config, IOM.HybSConfig},
+    ::String,
+    _abs_tol,
+    _rel_tol,
+    _i_bounds,
+) = bilin_cfg.quad_config
+
+function _converter_quad_config(
+    ::Union{IOM.NMDTBilinearConfig, IOM.DNMDTBilinearConfig},
+    quad_method::String,
+    abs_tol,
+    rel_tol,
+    i_bounds,
+)
+    Q = _quad_config_type(quad_method)
+    tol_sq = _resolve_tolerance(abs_tol, rel_tol, _max_abs(i_bounds)^2)
+    depth = IOM.tolerance_depth(Q; tolerance = tol_sq, max_delta = _max_delta(i_bounds))
+    return Q(; depth = depth)
+end
+
+# Add the bilinear `x·y` approximation, reusing the precomputed `ysq` (= the
+# loss `i_sq`) for the squares-based schemes and building `xsq` internally; the
+# discretization-based schemes ignore `ysq`/`quad_cfg` and take the raw form
+# (so no `xsq` is created — no model bloat).
+function _add_converter_bilinear!(
+    bilin_cfg::Union{IOM.Bin2Config, IOM.HybSConfig, IOM.NoBilinearApproxConfig},
+    quad_cfg,
+    container::OptimizationContainer,
+    ::Type{C},
+    names,
+    time_steps,
+    x_var,
+    y_var,
+    ysq,
+    x_bounds,
+    y_bounds,
+    meta::String,
+) where {C <: IS.InfrastructureSystemsComponent}
+    xsq = IOM._add_quadratic_approx!(
+        quad_cfg, container, C, names, time_steps,
+        x_var, x_bounds, meta * "_xsq",
+    )
+    return IOM._add_bilinear_approx!(
+        bilin_cfg, container, C, names, time_steps,
+        xsq, ysq, x_var, y_var, x_bounds, y_bounds, meta,
+    )
+end
+
+function _add_converter_bilinear!(
+    bilin_cfg::Union{IOM.NMDTBilinearConfig, IOM.DNMDTBilinearConfig},
+    quad_cfg,
+    container::OptimizationContainer,
+    ::Type{C},
+    names,
+    time_steps,
+    x_var,
+    y_var,
+    ysq,
+    x_bounds,
+    y_bounds,
+    meta::String,
+) where {C <: IS.InfrastructureSystemsComponent}
+    return IOM._add_bilinear_approx!(
+        bilin_cfg, container, C, names, time_steps,
+        x_var, y_var, x_bounds, y_bounds, meta,
+    )
+end

src/core/bilinear_configs.jl

@@ -0,0 +1,150 @@
+############################ Validation helpers ############################################
+
+function _validate_tolerance(tolerance::Float64)
+    (isfinite(tolerance) && tolerance > 0) || throw(
+        ArgumentError(
+            "bilinear approximation `tolerance` must be finite and > 0, got $tolerance",
+        ),
+    )
+    return tolerance
+end
+
+"""
+Characteristic magnitude of a variable over its per-device `bounds`
+(`max|x|` across all devices). Used to turn a relative tolerance into an
+absolute one — see [`_resolve_tolerance`](@ref).
+"""
+_max_abs(bounds) = maximum(max(abs(b.min), abs(b.max)) for b in bounds)
+
+"""
+Resolve the absolute bilinear/quadratic approximation tolerance from the
+`absolute` and `relative` attribute values. A relative tolerance is scaled to
+absolute by the characteristic product/term magnitude `scale`
+(`τ_abs = relative · scale`). Each argument is a positive number or `nothing`;
+the discretization must satisfy every tolerance that is set, so the effective
+absolute tolerance is the smallest of those provided. At least one must be set.
+"""
+function _resolve_tolerance(absolute, relative, scale::Float64)
+    tols = Float64[]
+    isnothing(absolute) || push!(tols, Float64(absolute))
+    isnothing(relative) || push!(tols, Float64(relative) * scale)
+    isempty(tols) && throw(
+        ArgumentError(
+            "at least one of `bilinear_absolute_tolerance` or " *
+            "`bilinear_relative_tolerance` must be set (both are unset)",
+        ),
+    )
+    return _validate_tolerance(minimum(tols))
+end
+
+function _quad_config_type(method::String)
+    if method == "solver_sos2"
+        return IOM.SolverSOS2QuadConfig
+    elseif method == "manual_sos2"
+        return IOM.ManualSOS2QuadConfig
+    elseif method == "sawtooth"
+        return IOM.SawtoothQuadConfig
+    elseif method == "nmdt"
+        return IOM.NMDTQuadConfig
+    elseif method == "dnmdt"
+        return IOM.DNMDTQuadConfig
+    else
+        error(
+            "Unsupported bilinear quadratic method \"$(method)\". " *
+            "Supported: \"solver_sos2\", \"manual_sos2\", \"sawtooth\", " *
+            "\"nmdt\", \"dnmdt\".",
+        )
+    end
+end
+
+const _BIN2_QUAD_METHODS = ("solver_sos2", "manual_sos2", "sawtooth", "nmdt", "dnmdt")
+const _HYBS_QUAD_METHODS = ("solver_sos2", "manual_sos2", "sawtooth")
+
+function _validate_quad_method(method::String, scheme::String, supported)
+    method in supported || error(
+        "Unsupported bilinear quadratic method \"$(method)\" for bilinear " *
+        "approximation \"$(scheme)\". Supported: " *
+        join(("\"$(m)\"" for m in supported), ", ") * ".",
+    )
+    return method
+end
+
+############################ Translation to IOM configs ####################################
+
+"""
+Build the IOM bilinear config consumed by `IOM._add_bilinear_approx!` from
+string attribute values, sizing the discretization from `tolerance` and the
+per-device domain widths (`delta_x`, `delta_y`).
+
+`method` selects the bilinear approximation scheme: `"bin2"`, `"hybs"`,
+`"nmdt"`, `"dnmdt"`, or `"none"`. `quad_method` selects the inner quadratic
+PWL method used by the `"bin2"` and `"hybs"` schemes (`"solver_sos2"`,
+`"manual_sos2"`, `"sawtooth"`, and — for `"bin2"` only — `"nmdt"`, `"dnmdt"`);
+it is ignored by the other schemes.
+
+Each IOM `tolerance_depth` / `tolerance_epigraph_depth` helper inverts its
+method's worst-case-gap bound and allocates the error budget across the inner
+quadratic, so POM never sizes the inner quad by hand — it just builds the inner
+quad at the returned `depth` (with the IOM-default `epigraph_depth`).
+
+Errors when `method` or `quad_method` is unrecognized, when `quad_method` is
+invalid for the selected scheme, or when `tolerance` is non-finite or ≤ 0.
+"""
+function _build_bilinear_config(
+    method::String,
+    quad_method::String,
+    tolerance::Float64,
+    delta_x::Float64,
+    delta_y::Float64,
+)
+    method == "none" && return IOM.NoBilinearApproxConfig()
+    _validate_tolerance(tolerance)
+    if method == "bin2"
+        _validate_quad_method(quad_method, method, _BIN2_QUAD_METHODS)
+        Q = _quad_config_type(quad_method)
+        depth = IOM.tolerance_depth(
+            IOM.Bin2Config{Q};
+            tolerance = tolerance,
+            max_delta_x = delta_x,
+            max_delta_y = delta_y,
+        )
+        return IOM.Bin2Config(Q(; depth))
+    elseif method == "hybs"
+        _validate_quad_method(quad_method, method, _HYBS_QUAD_METHODS)
+        Q = _quad_config_type(quad_method)
+        depth = IOM.tolerance_depth(
+            IOM.HybSConfig{Q};
+            tolerance = tolerance,
+            max_delta_x = delta_x,
+            max_delta_y = delta_y,
+        )
+        epigraph_depth = IOM.tolerance_epigraph_depth(
+            IOM.HybSConfig{Q};
+            tolerance = tolerance,
+            max_delta_x = delta_x,
+            max_delta_y = delta_y,
+        )
+        return IOM.HybSConfig(Q(; depth); epigraph_depth)
+    elseif method == "nmdt"
+        depth = IOM.tolerance_depth(
+            IOM.NMDTBilinearConfig;
+            tolerance = tolerance,
+            max_delta_x = delta_x,
+            max_delta_y = delta_y,
+        )
+        return IOM.NMDTBilinearConfig(; depth)
+    elseif method == "dnmdt"
+        depth = IOM.tolerance_depth(
+            IOM.DNMDTBilinearConfig;
+            tolerance = tolerance,
+            max_delta_x = delta_x,
+            max_delta_y = delta_y,
+        )
+        return IOM.DNMDTBilinearConfig(; depth)
+    else
+        error(
+            "Unsupported bilinear approximation \"$(method)\". " *
+            "Supported: \"bin2\", \"hybs\", \"nmdt\", \"dnmdt\", \"none\".",
+        )
+    end
+end

src/core/definitions.jl

@@ -73,7 +73,6 @@ const OBJECTIVE_FUNCTION_NEGATIVE = -1.0
 const INITIALIZATION_PROBLEM_HORIZON_COUNT = 3
 # The DEFAULT_RESERVE_COST value is used to avoid degeneracy of the solutions, reserve cost isn't provided.
 const DEFAULT_RESERVE_COST = 1.0
-const DEFAULT_INTERPOLATION_LENGTH = 4
 const KiB = 1024
 const MiB = KiB * KiB
 const GiB = MiB * KiB