[WIP] Add tolerance constructor to approximation configs

src/bilinear_approximations/bin2.jl

@@ -22,14 +22,83 @@ abstract type BilinearApproxConfig end
 Config for Bin2 bilinear approximation using z = ½((x+y)² − x² − y²).
 
 # Fields
-- `quad_config::QuadraticApproxConfig`: quadratic method used for x², y², and (x+y)²
+- `quad_config::Q`: quadratic method used for x², y², and (x+y)²
 - `add_mccormick::Bool`: whether to add reformulated McCormick cuts through separable variables (default true)
+
+The Q type parameter lets tolerance helpers dispatch on the inner quad method;
+see `tolerance_depth(::Type{Bin2Config{Q}}; …)`.
 """
-struct Bin2Config <: BilinearApproxConfig
-    quad_config::QuadraticApproxConfig
+struct Bin2Config{Q <: QuadraticApproxConfig} <: BilinearApproxConfig
+    quad_config::Q
     add_mccormick::Bool
+
+    Bin2Config(
+        quad_config::Q;
+        add_mccormick::Bool = true,
+    ) where {Q <: QuadraticApproxConfig} =
+        new{Q}(quad_config, add_mccormick)
+end
+
+# --- Tolerance helpers ---
+#
+# Notation: Δx, Δy are domain lengths (Δx = x_max − x_min); ε denotes errors.
+#
+# Bilinear identity:  xy = ½((x+y)² − x² − y²).
+# Approximation:      z  = ½(z_p − z_x − z_y), where z_• is the inner-quad
+# approximation of •² and (•)² for • ∈ {x, y, x+y}.
+#
+# Let ε_x = x² − z_x, ε_y = y² − z_y, ε_p = (x+y)² − z_p be the per-term
+# inner-quad errors. For one-sided-over inner quads (Sawtooth, SolverSOS2,
+# ManualSOS2), each ε_• ∈ [0, ε_•^max] where the bound scales as Δ²·c at
+# depth L (c is the inner quad's per-unit error coefficient):
+#   ε_x^max = Δx²·c,  ε_y^max = Δy²·c,  ε_p^max = (Δx+Δy)²·c.
+#
+# Substituting into z − xy yields  z − xy = ½(ε_x + ε_y − ε_p).
+# With each ε_• ∈ [0, ε_•^max], the range of z − xy is
+#   ½(0 + 0 − ε_p^max)  ≤  z − xy  ≤  ½(ε_x^max + ε_y^max − 0),
+# so |z − xy| ≤ max(½ε_p^max, ½(ε_x^max + ε_y^max)).
+#
+# Now (Δx+Δy)² = Δx² + 2ΔxΔy + Δy² ≥ Δx² + Δy² (since Δx, Δy ≥ 0), so
+# ε_p^max ≥ ε_x^max + ε_y^max, and the max collapses to ½ε_p^max.
+# To hit user-target τ, ask the inner quad for ε_p^max ≤ 2τ on Δx+Δy.
+
+"""
+    tolerance_depth(::Type{Bin2Config{Q}}; tolerance, max_delta_x, max_delta_y)::Int
+
+Inner-quad depth such that Bin2's worst-case overestimation gap is ≤ `tolerance`.
+Derivation: see the comment block above. Forwards to
+`tolerance_depth(Q; tolerance = 2·τ, max_delta = Δx + Δy)`.
+
+Defined for one-sided-over inner quads: `SawtoothQuadConfig`, `SolverSOS2QuadConfig`,
+`ManualSOS2QuadConfig`, `NMDTQuadConfig`, `DNMDTQuadConfig`. `EpigraphQuadConfig`
+is excluded — it is one-sided-under, so the sign of `ε_p` flips and the bound
+above no longer applies; an Epigraph inner quad can drive `z` arbitrarily far
+from `xy` under MIN/MAX objectives.
+
+**Caveat for NMDT/DNMDT inner Q**: these are only one-sided-over when their
+`epigraph_depth = 0`. With `epigraph_depth > 0`, the inner result becomes free
+in `[epigraph(x), nmdt(x)]`, which crosses `x²` and breaks the derivation. Pass
+NMDT/DNMDT inner Qs with `epigraph_depth = 0` only.
+"""
+function tolerance_depth(
+    ::Type{Bin2Config{Q}};
+    tolerance::Float64,
+    max_delta_x::Float64,
+    max_delta_y::Float64,
+) where {
+    Q <: Union{
+        SawtoothQuadConfig,
+        SolverSOS2QuadConfig,
+        ManualSOS2QuadConfig,
+        NMDTQuadConfig,
+        DNMDTQuadConfig,
+    },
+}
+    return tolerance_depth(Q;
+        tolerance = 2 * tolerance,
+        max_delta = max_delta_x + max_delta_y,
+    )
 end
-Bin2Config(quad_config::QuadraticApproxConfig) = Bin2Config(quad_config, true)
 
 # --- Unified bilinear approximation dispatch ---
 

src/bilinear_approximations/hybs.jl

@@ -9,21 +9,113 @@ struct HybSBoundConstraint <: ConstraintType end
 """
 Config for HybS (Hybrid Separable) bilinear approximation.
 
-Combines Bin2 lower bound and Bin3 upper bound with shared quadratic for x², y²
-and LP-only epigraph for (x+y)², (x−y)².
+Adds two inequalities that sandwich `z ≈ x·y`:
+- lower: `z ≥ ½(z_p1 − z_x − z_y)` where `z_p1 ≤ (x+y)²` is an epigraph (LP-only)
+  lower bound and `z_x ≥ x²`, `z_y ≥ y²` come from the inner quadratic `Q`.
+- upper: `z ≤ ½(z_x + z_y − z_p2)` where `z_p2 ≤ (x−y)²` is an epigraph (LP-only)
+  lower bound on the cross-difference.
 
 # Fields
-- `quad_config::QuadraticApproxConfig`: quadratic method used for the shared x² and y² terms
-- `epigraph_depth::Int`: depth for the epigraph Q^{L1} LP-only approximation of cross-terms (x±y)²
+- `quad_config::Q`: quadratic method for the shared x² and y² terms
+- `epigraph_depth::Int`: depth for the epigraph approximation of cross-terms (x±y)²
 - `add_mccormick::Bool`: whether to add standard McCormick envelope cuts on the product variable (default false)
+
+`Q` must be a one-sided-over quadratic approximation (so `z_x ≥ x²` and
+`z_y ≥ y²` hold). The reason is asymmetry: `z_x` and `z_y` appear with sign `−`
+in the lower bound and sign `+` in the upper bound, so if `z_x` could
+*under*-estimate `x²` then the `−z_x` term in the lower can drive `lower > xy`
+and the sandwich is invalid. This rules out `EpigraphQuadConfig` (one-sided
+under) and the two-sided `NMDTQuadConfig` / `DNMDTQuadConfig` for the tolerance
+helpers below; only `SawtoothQuadConfig`, `SolverSOS2QuadConfig`, and
+`ManualSOS2QuadConfig` are supported. See `tolerance_depth(::Type{HybSConfig{Q}};
+…)` and `tolerance_epigraph_depth(::Type{HybSConfig{Q}}; …)`.
 """
-struct HybSConfig <: BilinearApproxConfig
-    quad_config::QuadraticApproxConfig
+struct HybSConfig{Q <: QuadraticApproxConfig} <: BilinearApproxConfig
+    quad_config::Q
     epigraph_depth::Int
     add_mccormick::Bool
+
+    HybSConfig(
+        quad_config::Q;
+        epigraph_depth::Int,
+        add_mccormick::Bool = false,
+    ) where {Q <: QuadraticApproxConfig} =
+        new{Q}(quad_config, epigraph_depth, add_mccormick)
+end
+
+# --- Tolerance helpers ---
+#
+# Notation: Δx, Δy are domain lengths; ε denotes errors. Let
+#   ε_q = max(x² − z_x, y² − z_y, 0)     inner-quad one-sided-over error
+#   ε_e = max((x+y)² − z_p1, (x−y)² − z_p2, 0)   epigraph one-sided-under error
+# both ≥ 0 by construction.
+#
+# Where do the lower/upper expressions come from? Plug z_p1 = (x+y)² − ε_p1,
+# z_x = x² + ε_x, z_y = y² + ε_y into the lower-bound inequality:
+#   z ≥ ½(z_p1 − z_x − z_y)
+#     = ½((x+y)² − ε_p1 − x² − ε_x − y² − ε_y)
+#     = ½(2xy − ε_p1 − ε_x − ε_y)
+#     = xy − ½(ε_p1 + ε_x + ε_y)
+# So lower − xy ∈ [−ε_q − ½ε_e, 0]. Similarly substituting into
+# z ≤ ½(z_x + z_y − z_p2) gives upper − xy ∈ [0, ε_q + ½ε_e]. Combining,
+#   |z − xy| ≤ ε_q + ½ε_e.
+#
+# The total error is ε_q + ½ε_e. To meet τ, the helpers below allocate
+#   ε_q ≤ τ/2  →  inner Q at tolerance τ/2 over max_delta = max(Δx, Δy)
+#   ε_e ≤ τ    →  epigraph at tolerance τ   over max_delta = Δx + Δy
+# so the sum is ≤ τ/2 + ½·τ = τ. The 50/50 split between the two error
+# sources is an arbitrary design choice — any other allocation (e.g. 30/70)
+# that keeps the sum ≤ τ would also be valid; 50/50 is chosen for simplicity.
+#
+# Restricted to one-sided-over Q (see struct docstring for the asymmetry
+# argument). Other Q raise MethodError.
+
+"""
+    tolerance_depth(::Type{HybSConfig{Q}}; tolerance, max_delta_x, max_delta_y)::Int
+
+Inner-quad depth for HybS at target tolerance `τ`. Returns the smallest depth
+whose inner-quad error on `[ax, ax+Δx]` (and `[ay, ay+Δy]`) is `≤ τ/2`, which
+satisfies the `ε_q ≤ τ/2` half of the HybS budget.
+
+Only defined for `Q ∈ {SawtoothQuadConfig, SolverSOS2QuadConfig,
+ManualSOS2QuadConfig}` — the one-sided-over inner quads for which the HybS
+sandwich is valid. Other Q raise `MethodError`.
+
+See also `tolerance_epigraph_depth(::Type{HybSConfig{Q}}; …)` for the second knob.
+"""
+function tolerance_depth(
+    ::Type{HybSConfig{Q}};
+    tolerance::Float64,
+    max_delta_x::Float64,
+    max_delta_y::Float64,
+) where {Q <: Union{SawtoothQuadConfig, SolverSOS2QuadConfig, ManualSOS2QuadConfig}}
+    return tolerance_depth(Q;
+        tolerance = tolerance / 2,
+        max_delta = max(max_delta_x, max_delta_y),
+    )
+end
+
+"""
+    tolerance_epigraph_depth(::Type{HybSConfig{Q}}; tolerance, max_delta_x, max_delta_y)::Int
+
+Epigraph depth for HybS at target tolerance `τ`. Returns the smallest depth
+whose epigraph error on the cross-term range `Δx + Δy` is `≤ τ`, which
+satisfies the `½ε_e ≤ τ/2` half of the HybS budget.
+
+Only defined for `Q ∈ {SawtoothQuadConfig, SolverSOS2QuadConfig,
+ManualSOS2QuadConfig}`. Other Q raise `MethodError`.
+"""
+function tolerance_epigraph_depth(
+    ::Type{HybSConfig{Q}};
+    tolerance::Float64,
+    max_delta_x::Float64,
+    max_delta_y::Float64,
+) where {Q <: Union{SawtoothQuadConfig, SolverSOS2QuadConfig, ManualSOS2QuadConfig}}
+    return tolerance_depth(EpigraphQuadConfig;
+        tolerance = tolerance,
+        max_delta = max_delta_x + max_delta_y,
+    )
 end
-HybSConfig(quad_config::QuadraticApproxConfig, epigraph_depth::Int) =
-    HybSConfig(quad_config, epigraph_depth, false)
 
 # --- Unified HybS dispatch methods ---
 
@@ -145,7 +237,7 @@ function _add_bilinear_approx!(
     end
 
     # --- Epigraph Q^{L1} lower bound for (x+y)² and (x−y)² (no binaries) ---
-    epi_cfg = EpigraphQuadConfig(config.epigraph_depth)
+    epi_cfg = EpigraphQuadConfig(; depth = config.epigraph_depth)
     zp1_expr = _add_quadratic_approx!(
         epi_cfg,
         container, C, names, time_steps,

src/bilinear_approximations/nmdt.jl

@@ -10,19 +10,71 @@ Config for double-NMDT bilinear approximation (discretizes both x and y).
 
 # Fields
 - `depth::Int`: number of binary discretization levels L for both x and y
+
+The worst-case relaxation gap is `Δx·Δy·2^{-2L-2}`. See
+`tolerance_depth(::Type{DNMDTBilinearConfig}; …)` to derive `depth` from a target
+tolerance.
 """
 struct DNMDTBilinearConfig <: BilinearApproxConfig
     depth::Int
+
+    DNMDTBilinearConfig(; depth::Int) = new(depth)
+end
+
+"""
+    tolerance_depth(::Type{DNMDTBilinearConfig}; tolerance, max_delta_x, max_delta_y)::Int
+
+Smallest DNMDT bilinear depth `L` whose worst-case relaxation gap on
+`[ax, ax+Δx] × [ay, ay+Δy]` falls within `tolerance`. Inverts
+`Δx·Δy·2^{-2L-2} ≤ τ`:
+```
+L = ⌈(log₂(Δx·Δy/τ) − 2) / 2⌉
+```
+clamped to `L ≥ 1`.
+"""
+function tolerance_depth(
+    ::Type{DNMDTBilinearConfig};
+    tolerance::Float64,
+    max_delta_x::Float64,
+    max_delta_y::Float64,
+)
+    return _ceil_positive((log2(max_delta_x * max_delta_y / tolerance) - 2) / 2)
 end
 
 """
 Config for single-NMDT bilinear approximation (discretizes x only).
 
 # Fields
 - `depth::Int`: number of binary discretization levels L for x
+
+The worst-case relaxation gap is `Δx·Δy·2^{-L-2}`. See
+`tolerance_depth(::Type{NMDTBilinearConfig}; …)` to derive `depth` from a target
+tolerance.
 """
 struct NMDTBilinearConfig <: BilinearApproxConfig
     depth::Int
+
+    NMDTBilinearConfig(; depth::Int) = new(depth)
+end
+
+"""
+    tolerance_depth(::Type{NMDTBilinearConfig}; tolerance, max_delta_x, max_delta_y)::Int
+
+Smallest NMDT bilinear depth `L` whose worst-case relaxation gap on
+`[ax, ax+Δx] × [ay, ay+Δy]` falls within `tolerance`. Inverts
+`Δx·Δy·2^{-L-2} ≤ τ`:
+```
+L = ⌈log₂(Δx·Δy/τ) − 2⌉
+```
+clamped to `L ≥ 1`.
+"""
+function tolerance_depth(
+    ::Type{NMDTBilinearConfig};
+    tolerance::Float64,
+    max_delta_x::Float64,
+    max_delta_y::Float64,
+)
+    return _ceil_positive(log2(max_delta_x * max_delta_y / tolerance) - 2)
 end
 
 # --- DNMDT bilinear approximation ---

src/quadratic_approximations/common.jl

@@ -9,6 +9,14 @@ struct QuadraticExpression <: ExpressionType end
 "Abstract supertype for quadratic approximation method configurations."
 abstract type QuadraticApproxConfig end
 
+"""
+    _ceil_positive(x::Float64)::Int
+
+Smallest integer ≥ x, clamped to ≥ 1. Used by every `tolerance_depth` helper
+to convert a real-valued depth bound (e.g. `log₂(Δ²/τ)/2`) into a usable depth.
+"""
+_ceil_positive(x::Float64)::Int = max(1, ceil(Int, x))
+
 """
     _normed_variable!(container, C, names, time_steps, x_var, bounds, meta)
 
@@ -48,7 +56,11 @@ function _normed_variable!(
         IS.@assert_op b.max > b.min
         lx = b.max - b.min
         result = result_expr[name, t] = JuMP.AffExpr(0.0)
-        add_linear_to_jump_expression!(result, x_var[name, t], 1.0 / lx, -b.min / lx)
+        # add_proportional + add_constant accept any AbstractJuMPScalar for x_var,
+        # so this works for both VariableRef and AffExpr inputs (the latter is
+        # used by Bin2 to normalize the (x+y) expression).
+        add_proportional_to_jump_expression!(result, x_var[name, t], 1.0 / lx)
+        add_constant_to_jump_expression!(result, -b.min / lx)
     end
     return result_expr
 end

src/quadratic_approximations/epigraph.jl

@@ -18,9 +18,35 @@ Config for epigraph (Q^{L1}) LP-only lower-bound quadratic approximation.
 
 # Fields
 - `depth::Int`: number of tangent-line breakpoints (2^depth + 1 tangent lines); pure LP, zero binary variables
+
+The worst-case underestimation gap is `Δ²·2^{-2L-2}`. Per-segment derivation:
+adjacent tangents at `a` and `a+h` (h = Δ/2^L) meet at the midpoint `a+h/2`
+with value `a² + ah`; true `x² = a² + ah + h²/4`, so the gap is `h²/4 =
+Δ²·2^{-2L-2}`. See `tolerance_depth(::Type{EpigraphQuadConfig}; …)` to derive
+`depth` from a target tolerance.
 """
 struct EpigraphQuadConfig <: QuadraticApproxConfig
     depth::Int
+
+    EpigraphQuadConfig(; depth::Int) = new(depth)
+end
+
+"""
+    tolerance_depth(::Type{EpigraphQuadConfig}; tolerance, max_delta)::Int
+
+Smallest epigraph depth `L` whose worst-case underestimation gap on `[a, a+Δ]`
+falls within `tolerance`. Inverts the closed-form bound `Δ²·2^{-2L-2} ≤ τ`:
+```
+L = ⌈(log₂(Δ²/τ) − 2) / 2⌉
+```
+clamped to `L ≥ 1`.
+"""
+function tolerance_depth(
+    ::Type{EpigraphQuadConfig};
+    tolerance::Float64,
+    max_delta::Float64,
+)
+    return _ceil_positive((log2(max_delta^2 / tolerance) - 2) / 2)
 end
 
 """

src/quadratic_approximations/manual_sos2.jl

@@ -16,12 +16,40 @@ Config for manual binary-variable SOS2 quadratic approximation.
 # Fields
 - `depth::Int`: number of PWL segments (breakpoints = depth + 1)
 - `pwmcc_segments::Int`: number of piecewise McCormick cut partitions; 0 to disable (default 4)
+
+The worst-case PWL overestimation gap is `Δ²/(4·depth²)`. `pwmcc_segments` is
+an LP-relaxation tightener (adds piecewise-McCormick cuts on the concave
+relaxation surface); it does **not** change the MIP-optimal worst-case error,
+only the LP relaxation given to branch-and-bound. See
+`tolerance_depth(::Type{ManualSOS2QuadConfig}; …)` to derive `depth` from a target
+tolerance.
 """
 struct ManualSOS2QuadConfig <: QuadraticApproxConfig
     depth::Int
     pwmcc_segments::Int
+
+    ManualSOS2QuadConfig(; depth::Int, pwmcc_segments::Int = 4) =
+        new(depth, pwmcc_segments)
+end
+
+"""
+    tolerance_depth(::Type{ManualSOS2QuadConfig}; tolerance, max_delta)::Int
+
+Smallest SOS2 segment count `d` whose worst-case PWL gap on `[a, a+Δ]` falls
+within `tolerance`. Inverts `Δ²/(4·d²) ≤ τ`:
+```
+d = ⌈Δ / (2·√τ)⌉
+```
+clamped to `d ≥ 1`. `pwmcc_segments` does not enter the error bound, so it is
+left to the constructor.
+"""
+function tolerance_depth(
+    ::Type{ManualSOS2QuadConfig};
+    tolerance::Float64,
+    max_delta::Float64,
+)
+    return _ceil_positive(max_delta / (2 * sqrt(tolerance)))
 end
-ManualSOS2QuadConfig(depth::Int) = ManualSOS2QuadConfig(depth, 4)
 
 """
     _add_quadratic_approx!(config::ManualSOS2QuadConfig, container, C, names, time_steps, x_var, bounds, meta)