dnmt

.claude/claude.md

@@ -45,7 +45,6 @@ src/
     add_jump_expressions.jl   # JuMP expression builders
     add_param_container.jl    # Parameter container builders
     add_constraint_dual.jl    # Constraint & dual helpers
-    add_pwl_methods.jl        # Piecewise-linear variable/constraint methods
     constraint_helpers.jl     # Generic constraint utilities
     range_constraint.jl       # Min/max range constraints
     duration_constraints.jl   # Min up/down time constraints
@@ -82,6 +81,17 @@ src/
     initial_conditions_update_in_memory_store.jl
     model_numerical_analysis_utils.jl
     optimization_debugging.jl # Debug/diagnostic tools
+  quadratic_approximations/  # PWL approximation of x² (univariate)
+    pwl_utils.jl              # Shared breakpoint generator, _square helper
+    incremental.jl            # Incremental PWL formulation (δ/z variables)
+    solver_sos2.jl            # Solver-native SOS2 quadratic approximation
+    manual_sos2.jl            # Manual SOS2 (binary adjacency) quadratic approximation
+    sawtooth.jl               # Sawtooth relaxation approximation
+    epigraph.jl               # Epigraph-based approximation
+  bilinear_approximations/  # Approximation of bilinear products x·y
+    mccormick.jl              # McCormick envelopes for bilinear terms
+    bilinear.jl               # Bin2 separable bilinear approximation
+    hybs.jl                   # HybS hybrid separable approximation
   utils/                  # General-purpose utilities
     jump_utils.jl             # JuMP helper functions
     dataframes_utils.jl       # DataFrame manipulation
@@ -142,6 +152,12 @@ scripts/formatter/        # Code formatting (JuliaFormatter)
   concrete formulations call into.
 - **`src/objective_function/`** translates cost curves into JuMP objective terms. Each cost
   curve type has its own file.
+- **`src/quadratic_approximations/`** implements PWL approximation methods for x²:
+  SOS2 (solver and manual), sawtooth, epigraph, plus the incremental formulation
+  and shared breakpoint utilities.
+- **`src/bilinear_approximations/`** implements approximation methods for bilinear
+  products x·y: Bin2 separable decomposition, HybS hybrid relaxation, and McCormick
+  envelopes.
 - **`src/operation/`** implements `DecisionModel` and `EmulationModel` — the two main model
   types — plus serialization, output stores, and the problem template.
 - **`src/utils/`** is for pure utility functions with no domain coupling.

docs/make.jl

@@ -20,7 +20,8 @@ pages = OrderedDict(
         "Developers" => ["Developer Guidelines" => "reference/developer_guidelines.md",
         "Internals" => "reference/internal.md"],
         "Public API" => "reference/public.md",
-        "Quadratic Approximations" => "reference/quadratic_approximations.md",
+        "Piecewise Linear Cost Functions" => "explanation/piecewise_linear_cost_functions.md",
+        "Quadratic Approximations" => "explanation/quadratic_approximations.md",
     ],
 )
 

docs/src/explanation/piecewise_linear_cost_functions.md

@@ -0,0 +1,188 @@
+# Piecewise Linear Cost Functions
+
+```@meta
+CurrentModule = InfrastructureOptimizationModels
+```
+
+This package provides two formulations for representing piecewise linear (PWL) cost
+functions in optimization models. Both formulations express an operating point and its
+cost as a point on a piecewise linear curve connecting breakpoints
+``(P_0, C_0), (P_1, C_1), \ldots, (P_n, C_n)``, but they approach the problem differently.
+
+## Lambda Formulation (Convex Combination)
+
+Used by `CostCurve{PiecewisePointCurve}` and `FuelCurve{PiecewisePointCurve}`.
+
+The lambda formulation assigns a weight ``\lambda_i`` to each **breakpoint**. The
+operating point and cost are expressed as a weighted average of the breakpoint values.
+
+### Formulation
+
+Given ``n + 1`` breakpoints with power levels ``P_0, \ldots, P_n`` and costs
+``C(P_0), \ldots, C(P_n)``:
+
+**Variables:**
+
+```math
+\lambda_i \in [0, 1], \quad i = 0, 1, \ldots, n
+```
+
+**Constraints:**
+
+```math
+\begin{aligned}
+p &= \sum_{i=0}^{n} \lambda_i \, P_i && \text{(linking)} \\
+C(p) &= \sum_{i=0}^{n} \lambda_i \, C(P_i) && \text{(cost)} \\
+\sum_{i=0}^{n} \lambda_i &= u && \text{(normalization, } u \text{ = on-status)} \\
+\text{at most two adjacent } &\lambda_i \text{ may be nonzero} && \text{(adjacency)}
+\end{aligned}
+```
+
+### Adjacency Enforcement
+
+The adjacency condition keeps the operating point on a single line segment. For
+**convex** cost curves (where each segment is more expensive than the last), the
+optimizer naturally satisfies this condition — no extra constraints are needed.
+
+For **non-convex** cost curves, an
+[SOS2 constraint](https://en.wikipedia.org/wiki/Special_ordered_set) is added to
+enforce adjacency explicitly. This requires solver support for SOS2, and effectively
+introduces additional branching in the solver.
+
+### Compact Form
+
+When the power variable represents output above minimum generation
+(PowerAboveMinimumVariable), a compact variant adjusts the linking constraint to
+include a ``P_{\min}`` offset:
+
+```math
+\sum_{i=0}^{n} \lambda_i \, P_i = p + P_{\min} \cdot u
+```
+
+### Variables and Constraints
+
+| Container Type                                       | Description                                  |
+|:---------------------------------------------------- |:-------------------------------------------- |
+| [`PiecewiseLinearCostVariable`](@ref)                | Lambda weights ``\lambda_i \in [0, 1]``      |
+| [`PiecewiseLinearCostConstraint`](@ref)              | Links power variable to weighted breakpoints |
+| [`PiecewiseLinearCostNormalizationConstraint`](@ref) | Ensures ``\sum \lambda_i = u``               |
+
+### Implementation
+
+| Function                        | Role                                                           |
+|:------------------------------- |:-------------------------------------------------------------- |
+| `_add_pwl_variables!`           | Creates ``n + 1`` lambda variables per time step               |
+| `_add_pwl_constraint_standard!` | Adds linking and normalization constraints                     |
+| `_add_pwl_constraint_compact!`  | Compact form with ``P_{\min}`` offset                          |
+| `_add_pwl_term!`                | Orchestrates variables, constraints, SOS2, and cost expression |
+| `add_pwl_sos2_constraint!`      | Adds SOS2 adjacency for non-convex curves                      |
+
+## Delta Formulation (Incremental / Block Offers)
+
+Used by `MarketBidCost` and `ImportExportCost` offer curves.
+
+The delta formulation assigns a variable ``\delta_k`` to each **segment**, representing
+how much of that segment has been used. The operating point is the sum of segment
+contributions from the first breakpoint.
+
+### Formulation
+
+Given ``n`` segments with breakpoints ``P_0, \ldots, P_n`` and marginal costs (slopes)
+``m_k = \frac{C(P_k) - C(P_{k-1})}{P_k - P_{k-1}}``:
+
+**Variables:**
+
+```math
+\delta_k \geq 0, \quad k = 1, \ldots, n
+```
+
+**Constraints:**
+
+```math
+\begin{aligned}
+p &= \sum_{k=1}^{n} \delta_k + P_{\min,\text{offset}} && \text{(linking)} \\
+\delta_k &\leq P_k - P_{k-1} && \text{(segment capacity)} \\
+C(p) &= \sum_{k=1}^{n} m_k \, \delta_k && \text{(cost)}
+\end{aligned}
+```
+
+### Convexity Advantage
+
+For convex offer curves (``m_1 \leq m_2 \leq \cdots \leq m_n``), the fill-order
+condition is automatically satisfied. A cost-minimizing optimizer fills cheap segments
+before expensive ones, so no SOS2 constraints or binary variables are needed. This is
+the common case in power systems, where generator marginal costs are typically
+increasing.
+
+### Variables and Constraints
+
+| Container Type                                           | Description                            |
+|:-------------------------------------------------------- |:-------------------------------------- |
+| [`PiecewiseLinearBlockIncrementalOffer`](@ref)           | Delta variables for incremental offers |
+| [`PiecewiseLinearBlockDecrementalOffer`](@ref)           | Delta variables for decremental offers |
+| [`PiecewiseLinearBlockIncrementalOfferConstraint`](@ref) | Linking + capacity for incremental     |
+| [`PiecewiseLinearBlockDecrementalOfferConstraint`](@ref) | Linking + capacity for decremental     |
+
+### Implementation
+
+| Function                           | Role                                                     |
+|:---------------------------------- |:-------------------------------------------------------- |
+| `add_pwl_variables!`               | Creates ``n`` delta variables per time step              |
+| `_add_pwl_constraint!`             | Adds linking and segment capacity constraints            |
+| `add_pwl_term!`                    | Orchestrates variables, constraints, and cost expression |
+| `get_pwl_cost_expression`          | Builds ``\sum m_k \, \delta_k`` cost expression          |
+| `add_pwl_block_offer_constraints!` | Low-level block-offer constraint builder                 |
+
+## Comparison
+
+|                      | Lambda (``\lambda``)               | Delta (``\delta``)                  |
+|:-------------------- |:---------------------------------- |:----------------------------------- |
+| **Thinks about**     | Breakpoints (the dots)             | Segments (the lines)                |
+| **Variables**        | ``n + 1`` (one per breakpoint)     | ``n`` (one per segment)             |
+| **Output equation**  | ``p = \sum \lambda_i \, P_i``      | ``p = \sum \delta_k + P_{\min}``    |
+| **Cost equation**    | ``C = \sum \lambda_i \, C(P_i)``   | ``C = \sum m_k \, \delta_k``        |
+| **Adjacency rule**   | Must be enforced explicitly (SOS2) | Often automatic (convex case)       |
+| **Binary variables** | Usually needed (non-convex)        | Often not needed (convex case)      |
+| **Used by**          | `CostCurve`, `FuelCurve`           | `MarketBidCost`, `ImportExportCost` |
+
+## When Does the Choice Matter?
+
+For **convex** cost functions — where each segment costs more than the last — the delta
+formulation is typically simpler and faster. The optimizer fills cheap segments first on
+its own, so no extra adjacency constraints are needed. This is the common case in power
+systems.
+
+For **non-convex** cost functions — where a segment might be cheaper than the one before
+it — both formulations need additional constraints (SOS2 or binary variables) to force
+the correct ordering. In that case the lambda formulation is the traditional choice.
+
+## Incremental Interpolation (General PWL Approximation)
+
+In addition to cost functions, the package provides a general-purpose incremental
+(delta) method for approximating arbitrary nonlinear functions as piecewise linear. This
+is used for PWL approximation of constraints (e.g., loss functions, quadratic terms).
+
+### Formulation
+
+Given breakpoints ``(x_1, y_1), \ldots, (x_{n+1}, y_{n+1})`` where ``y_i = f(x_i)``,
+and interpolation variables ``\delta_i \in [0, 1]`` with binary ordering variables
+``z_i \in \{0, 1\}``:
+
+```math
+\begin{aligned}
+x &= x_1 + \sum_{i=1}^{n} \delta_i \, (x_{i+1} - x_i) \\
+y &= y_1 + \sum_{i=1}^{n} \delta_i \, (y_{i+1} - y_i) \\
+z_i &\geq \delta_{i+1}, \quad z_i \leq \delta_i \quad \text{(ordering)}
+\end{aligned}
+```
+
+The ordering constraints ensure segments are filled sequentially: ``\delta_1`` must be
+filled before ``\delta_2`` can begin.
+
+### Implementation
+
+| Function                                             | Role                                        |
+|:---------------------------------------------------- |:------------------------------------------- |
+| `_get_breakpoints_for_pwl_function`                  | Generates equally-spaced breakpoints        |
+| `add_sparse_pwl_interpolation_variables!`            | Creates ``\delta`` and ``z`` variables      |
+| `_add_generic_incremental_interpolation_constraint!` | Adds interpolation and ordering constraints |

docs/src/reference/quadratic_approximations.md → docs/src/explanation/quadratic_approximations.md

@@ -99,7 +99,7 @@ linear constraints per component per time step. Same approximation quality as So
 ## Sawtooth
 
 Approximates ``x^2`` using the recursive sawtooth MIP formulation from
-Beach, Burlacu, Hager, and Hildebrand (2024). This method requires only ``O(\log(1/\varepsilon))``
+Beach, Burlacu, Hager, and Hildebrand (2023). This method requires only ``O(\log(1/\varepsilon))``
 binary variables to achieve error ``\varepsilon``, compared to ``O(1/\sqrt{\varepsilon})``
 for the SOS2 methods.
 
@@ -162,3 +162,97 @@ constraints per component per time step. The approximation interpolates ``x^2``
 To match the number of breakpoints between methods, set ``S = 2^L``. At equal breakpoint
 count the approximation quality is identical, but the sawtooth uses ``L = \log_2 S``
 binary variables instead of ``S``.
+
+## Error Scaling
+
+Both SOS2 and sawtooth produce the same PWL interpolation of ``x^2`` when they use the
+same number of uniform breakpoints. With ``n`` uniform segments on ``[0, 1]``, the
+interpolant ``F_n(x)`` satisfies the classical error bound (Barmann et al., 2023):
+
+```math
+0 \leq F_n(x) - x^2 \leq \frac{1}{4n^2} \quad \text{for all } x \in [0, 1]
+```
+
+The maximum error is attained at the midpoint of each segment. Since both methods
+interpolate at the same breakpoints, the pointwise error is identical — the difference
+lies entirely in how efficiently the segments are encoded.
+
+### Same number of binary variables
+
+If we budget ``L`` binary/integer variables for each method, SOS2 gets ``L`` segments
+while sawtooth gets ``2^L`` segments. The error gap grows exponentially:
+
+| ``L`` | SOS2 segments | SOS2 error            | Sawtooth segments | Sawtooth error        | Ratio |
+|:----- |:------------- |:--------------------- |:----------------- |:--------------------- |:----- |
+| 1     | 1             | ``2.50\times10^{-1}`` | 2                 | ``6.25\times10^{-2}`` | 4x    |
+| 2     | 2             | ``6.25\times10^{-2}`` | 4                 | ``1.56\times10^{-2}`` | 4x    |
+| 3     | 3             | ``2.78\times10^{-2}`` | 8                 | ``3.91\times10^{-3}`` | 7x    |
+| 4     | 4             | ``1.56\times10^{-2}`` | 16                | ``9.77\times10^{-4}`` | 16x   |
+| 5     | 5             | ``1.00\times10^{-2}`` | 32                | ``2.44\times10^{-4}`` | 41x   |
+| 6     | 6             | ``6.94\times10^{-3}`` | 64                | ``6.10\times10^{-5}`` | 114x  |
+| 7     | 7             | ``5.10\times10^{-3}`` | 128               | ``1.53\times10^{-5}`` | 334x  |
+| 8     | 8             | ``3.91\times10^{-3}`` | 256               | ``3.81\times10^{-6}`` | 1024x |
+
+SOS2 error decays **polynomially** as ``O(1/L^2)``; sawtooth error decays
+**exponentially** as ``O(4^{-L})``.
+
+### Same number of segments
+
+When both methods use the same number of uniform segments ``n``, they produce identical
+PWL interpolations, so the approximation error is the same. The difference is in
+formulation size: SOS2 needs ``O(n)`` binary variables, sawtooth needs only
+``\log_2(n)``.
+
+| Segments ``n`` | Error                 | SOS2 vars | Sawtooth vars | Var ratio |
+|:-------------- |:--------------------- |:--------- |:------------- |:--------- |
+| 2              | ``6.25\times10^{-2}`` | 1         | 1             | 1.0x      |
+| 4              | ``1.56\times10^{-2}`` | 3         | 2             | 1.5x      |
+| 8              | ``3.91\times10^{-3}`` | 7         | 3             | 2.3x      |
+| 16             | ``9.77\times10^{-4}`` | 15        | 4             | 3.8x      |
+| 32             | ``2.44\times10^{-4}`` | 31        | 5             | 6.2x      |
+| 64             | ``6.10\times10^{-5}`` | 63        | 6             | 10.5x     |
+| 128            | ``1.53\times10^{-5}`` | 127       | 7             | 18.1x     |
+| 256            | ``3.81\times10^{-6}`` | 255       | 8             | 31.9x     |
+
+## Extension to Bilinear Terms
+
+Bilinear terms ``z = xy`` arise throughout optimization (energy systems, pooling problems,
+gas networks). The standard univariate reformulation uses the identity:
+
+```math
+xy = \frac{1}{4}\left[(x + y)^2 - (x - y)^2\right]
+```
+
+Each squared term is approximated independently with one of the methods above. With both
+terms at the same refinement level, the bilinear error satisfies:
+
+```math
+\varepsilon_{xy} \leq \frac{1}{2} \, \varepsilon_{\text{quad}}
+```
+
+All scaling relationships from the univariate case carry over with a constant factor of
+``1/2``. The exponential gap between SOS2 and sawtooth persists.
+
+## Practical Considerations
+
+**SOS2** has the advantage of being natively supported by commercial solvers (Gurobi,
+CPLEX) with specialized branching rules. Both the solver SOS2 and manual SOS2
+formulations produce sharp LP relaxations for convex functions.
+
+**Sawtooth** introduces auxiliary continuous variables and big-M-type constraints, which
+may interact less favorably with presolve and cutting planes. However, Beach et al. (2023)
+show that the sawtooth relaxation is **sharp** (its LP relaxation equals the convex hull)
+and **hereditarily sharp** (sharpness is preserved at every node in the branch-and-bound
+tree). This strong theoretical property, combined with exponentially fewer binary
+variables, can lead to significant solver performance gains for problems requiring high
+approximation accuracy.
+
+Whether the tighter approximation at a given variable budget outweighs the structural
+advantages of SOS2 depends on the specific problem and solver.
+
+## References
+
+  - Beach, B., Burlacu, R., Barmann, A., Hager, L., Kleinert, T. (2023). *Enhancements of discretization approaches for non-convex MIQCQPs*. Journal of Global Optimization.
+  - Barmann, A., Burlacu, R., Hager, L., Kleinert, T. (2023). *On piecewise linear approximations of bilinear terms: structural comparison of univariate and bivariate MIP formulations*. Journal of Global Optimization, 85, 789-819.
+  - Yarotsky, D. (2017). *Error bounds for approximations with deep ReLU networks*. Neural Networks, 94, 103-114.
+  - Huchette, J.A. (2018). *Advanced mixed-integer programming formulations: methodology, computation, and application*. PhD thesis, MIT.

docs/src/reference/developer_guidelines.md

@@ -1,6 +1,6 @@
 # Developer Guidelines
 
-In order to contribute to `PowerSystems.jl` repository please read the following sections of
+In order to contribute to `InfrastructureOptimizationModels.jl` repository please read the following sections of
 [`InfrastructureSystems.jl`](https://github.com/NREL-Sienna/InfrastructureSystems.jl)
 documentation in detail: