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,6 @@ export HydroWaterFactorModel
 export HydroWaterModelReservoir
 export HydroTurbineBilinearDispatch
 export HydroTurbineWaterLinearDispatch
-export HydroTurbineBin2BilinearDispatch
 export HydroTurbineWaterLinearCommitment
 export HydroEnergyModelReservoir
 export HydroTurbineEnergyDispatch
@@ -810,15 +810,13 @@ export HVDCTwoTerminalLossless
 export HVDCTwoTerminalDispatch
 export HVDCTwoTerminalPiecewiseLoss
 export HVDCTwoTerminalLCC
-export HVDCTwoTerminalVSCNLP
-export HVDCTwoTerminalVSCLP
+export HVDCTwoTerminalVSC
 
 # Converter Formulations
 export LosslessConverter
 export LinearLossConverter
 export AbstractQuadraticLossConverter
-export QuadraticLossConverterMILP
-export QuadraticLossConverterNLP
+export QuadraticLossConverter
 
 # DC Line Formulations
 export DCLosslessLine

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,36 +1731,31 @@ 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",
     )
 
-    _register_pq_sq_expressions!(
-        container, devices, line_names, time_steps, device_model,
+    _register_vsc_apparent_power_squares!(
+        bilin_cfg, container, devices, line_names, time_steps, device_model,
         network_model,
     )
 
@@ -1784,9 +1770,12 @@ function construct_device!(
     add_constraints!(
         container, HVDCVSCConverterPowerConstraint, devices, device_model, network_model,
     )
-    _maybe_add_reactive_power_constraints!(
-        container, devices, device_model, network_model,
-        HVDCVSCApparentPowerLimitConstraint,
+    _add_vsc_apparent_power_limit!(
+        bilin_cfg,
+        container,
+        devices,
+        device_model,
+        network_model,
     )
 
     add_constraint_dual!(container, sys, device_model)

src/common_models/quadratic_converter_loss.jl

@@ -1,8 +1,9 @@
 # Shared helpers for quadratic / two-term converter losses
 #   loss(I) = a * I^2 + b * |I| + c
-# Used by multi-terminal InterconnectingConverter formulations
-# (QuadraticLossConverterMILP, QuadraticLossConverterNLP) and two-terminal
-# TwoTerminalVSCLine formulations (HVDCTwoTerminalVSCLP, HVDCTwoTerminalVSCNLP).
+# Used by the multi-terminal InterconnectingConverter formulation
+# (QuadraticLossConverter) and the two-terminal TwoTerminalVSCLine formulation
+# (HVDCTwoTerminalVSC). Both default to the exact (NLP) case and opt into MILP
+# approximations via the `"bilinear_approximation"` attribute.
 #
 # `|I|` is represented by an LP surrogate: a single non-negative variable
 # `CurrentAbsoluteValueVariable` bounded below by both `i` and `-i`. The
@@ -91,3 +92,111 @@ function _add_abs_value_constraints!(
     end
     return
 end
+
+# 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.
+function _build_converter_configs(
+    ::Type{F},
+    model::DeviceModel,
+    v_bounds,
+    i_bounds,
+) where {F <: Union{QuadraticLossConverter, HVDCTwoTerminalVSC}}
+    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
+
+# 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