diff --git a/README.md b/README.md
index 543d1f7a..ea8d3f70 100644
--- a/README.md
+++ b/README.md
@@ -184,13 +184,12 @@ The following reformulation methods are currently supported:
 
     All variables must be included in exactly one partition. For manual partitioning, ensure each variable appears in exactly one group. For automatic partitioning, variables are divided as evenly as possible among the specified number of partitions.
 
-6. [Cutting Planes](https://pubsonline.informs.org/doi/10.1287/ijoc.2015.0669): This method iteratively generates cutting planes using a separation problem and a relaxed Big-M formulation, then applies a final reformulation method. The `CuttingPlanes` struct is created with the following arguments:
+6. [Cutting Planes](https://pubsonline.informs.org/doi/10.1287/ijoc.2015.0669): This method iteratively generates cutting planes using a separation problem and a relaxed Big-M formulation. The `CuttingPlanes` struct is created with the following arguments:
 
     - `optimizer`: Optimizer to use when solving the separation and relaxed Big-M subproblems. This is a required value.
     - `max_iter`: Maximum number of cutting plane iterations. Default: `3`.
     - `seperation_tolerance`: Convergence tolerance for the separation problem objective. Default: `1e-6`.
-    - `final_reform_method`: Reformulation method to apply after cutting plane iterations. Default: `BigM()`.
-    - `M_value`: Big-M value to use in the relaxed Big-M reformulation during iterations. Default: `1e9`.
+    - `M_value`: Big-M value to use in the relaxed Big-M reformulation. Default: `1e9`.
 
 ## Infinite-Dimensional GDP
 To model disjunctions, logical variables, and logical constraints with infinite-dimensional optimization problems (e.g., dynamic and stochastic optimization), DisjunctiveProgramming is also compatible with [InfiniteOpt.jl](https://github.com/infiniteopt/InfiniteOpt.jl). For this, the syntax is largely the same, users simply need to import `InfiniteOpt` and use `InfiniteGDPModel`. They also can use `InfiniteLogical` to declare infinite logical variables as shown below:
diff --git a/docs/src/index.md b/docs/src/index.md
index aa1902e3..d342515a 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -61,7 +61,7 @@ A [`GDPModel`](@ref) is a `JuMP Model` with a [`GDPData`](@ref) field in the mod
 - `Logical Constraints`: Selector (cardinality) or proposition (Boolean) constraints describing the relationships between the logical variables.
 - `Disjunct Constraints`: Constraints associated with each disjunct in the model.
 - `Disjunctions`: Disjunction constraints.
-- `Solution Method`: The reformulation technique or solution method. Currently, supported methods include Big-M, Hull, and Indicator Constraints.
+- `Solution Method`: The reformulation technique or solution method. Currently, supported methods include Big-M, Multiple Big-M, Hull, P-Split, Indicator Constraints, and Logic-based Outer Approximation (LOA).
 - `Reformulation Variables`: List of JuMP variables created when reformulating a GDP model into a MIP model.
 - `Reformulation Constraints`: List of constraints created when reformulating a GDP model into a MIP model.
 - `Ready to Optimize`: Flag indicating if the model can be optimized.
@@ -163,9 +163,15 @@ The following reformulation methods are currently supported:
 
 1. [Big-M](https://optimization.cbe.cornell.edu/index.php?title=Disjunctive_inequalities#Big-M_Reformulation[1][2]): The [`BigM`](@ref) struct is used.
 
-2. [Hull](https://optimization.cbe.cornell.edu/index.php?title=Disjunctive_inequalities#Convex-Hull_Reformulation[1][2]): The [`Hull`](@ref) struct is used.
+2. Multiple Big-M: A tighter big-M reformulation that computes a separate big-M value for each disjunct constraint (rather than a single global value) by solving auxiliary mini-models. This is invoked with the [`MBM`](@ref) struct, which requires an optimizer.
 
-3. [Indicator](https://jump.dev/JuMP.jl/stable/manual/constraints/#Indicator-constraints): This method reformulates each disjunct constraint into an indicator constraint with the Boolean reformulation counterpart of the Logical variable used to define the disjunct constraint. This is invoked with [`Indicator`](@ref).
+3. [Hull](https://optimization.cbe.cornell.edu/index.php?title=Disjunctive_inequalities#Convex-Hull_Reformulation[1][2]): The [`Hull`](@ref) struct is used.
+
+4. P-Split: A partitioned reformulation that splits the variables into groups and applies a P-split formulation to each group, giving a relaxation between Big-M and Hull in tightness. This is invoked with the [`PSplit`](@ref) struct.
+
+5. [Indicator](https://jump.dev/JuMP.jl/stable/manual/constraints/#Indicator-constraints): This method reformulates each disjunct constraint into an indicator constraint with the Boolean reformulation counterpart of the Logical variable used to define the disjunct constraint. This is invoked with [`Indicator`](@ref).
+
+6. Logic-based Outer Approximation (LOA): An iterative solution method for nonlinear GDPs, rather than a single-shot reformulation. It alternates between a primal NLP (the model reformulated by an inner method — `BigM`, `MBM`, or `Hull` — with the disjunct binaries fixed at the current selection) and a master MILP that accumulates outer-approximation and no-good cuts until the bound meets the incumbent. This is invoked with the [`LOA`](@ref) struct, e.g. `optimize!(m, gdp_method = LOA(Ipopt.Optimizer))`. Unlike the reformulation methods above, `LOA` runs the full algorithm on every `optimize!` call and cannot be passed to `reformulate_model`.
 
 ## Release Notes
 
diff --git a/ext/InfiniteDisjunctiveProgramming.jl b/ext/InfiniteDisjunctiveProgramming.jl
index 5f77dce7..7c9e6734 100644
--- a/ext/InfiniteDisjunctiveProgramming.jl
+++ b/ext/InfiniteDisjunctiveProgramming.jl
@@ -253,7 +253,7 @@ function _interpolate_at(
     args::NTuple{N, <:Real}
     ) where {N}
     # lower-corner cell index per dimension
-    idx_lo = ntuple(d -> 
+    idx_lo = ntuple(d ->
         clamp(searchsortedlast(grids[d], args[d]),1, length(grids[d]) - 1), N
     )
     # max over the 2^N corners; bit d of k picks lower or upper
@@ -324,7 +324,10 @@ function DP.copy_and_reformulate(
     end
     sub = DP.GDPSubmodel(sub_copy, decision_vars, fwd_map)
     JuMP.set_optimizer(sub.model, method.optimizer)
-    JuMP.set_silent(sub.model)
+    # NOTE: previously called JuMP.set_silent(sub.model) here, but
+    # that hides the Gurobi B&B trace even when the caller passes
+    # OutputFlag=1 / LogToConsole=1 via optimizer_with_attributes.
+    # The caller is now responsible for silencing if desired.
     return sub
 end
 
@@ -374,4 +377,133 @@ function DP.add_cut(
     return
 end
 
+################################################################################
+#                       LOA FOR INFINITEMODEL
+################################################################################
+# LOA override for InfiniteModel: after the inner reformulation, build
+# the transformation backend and hand the LOA loop the backend model —
+# a plain scalar JuMP model whose variables, measures, and constraints
+# are already transcribed per support. The loop runs entirely on base
+# code and solves the backend directly, which never invalidates it, so
+# the committed fixes and warm starts survive the final `optimize!`.
+
+# Per-support transcribed refs of an InfiniteOpt object, as a Vector
+# (finite objects give one element).
+_transcribed_refs(x::AbstractArray) = vec(x)
+_transcribed_refs(x) = [x]
+
+# Index a per-support vector, broadcasting a single-element (finite)
+# vector across all supports.
+_at(values::Vector, k::Int) =
+    length(values) == 1 ? values[1] : values[k]
+
+# Per-support transcribed binary references for an indicator: the
+# transcribed underlying binary, complemented for a complement-form
+# indicator.
+function _transcribed_binary_refs(model::InfiniteOpt.InfiniteModel, indicator)
+    binary_ref = DP._indicator_to_binary(model)[indicator]
+    binaries = _transcribed_refs(InfiniteOpt.transformation_variable(
+        DP._underlying_binary(binary_ref)))
+    binary_ref isa JuMP.GenericAffExpr || return binaries
+    return [1.0 - binary for binary in binaries]
+end
+
+# Transcribe the inner Hull disaggregation map per support. A finite
+# variable under an infinite indicator keys its single disaggregated
+# copy by each per-support binary reference, mirroring the per-support
+# disjunct cuts that look it up.
+function _transcribed_disaggregation_map(model::InfiniteOpt.InfiniteModel)
+    result = Dict{Tuple{JuMP.VariableRef,
+        Union{JuMP.VariableRef, JuMP.AffExpr}}, JuMP.VariableRef}()
+    for ((variable, indicator), disaggregated) in DP._disaggregation_map(model)
+        binary_refs = _transcribed_binary_refs(model, indicator)
+        variables = _transcribed_refs(
+            InfiniteOpt.transformation_variable(variable))
+        disaggregateds = _transcribed_refs(InfiniteOpt.transformation_variable(
+            disaggregated))
+        n = max(length(binary_refs), length(variables), length(disaggregateds))
+        for k in 1:n
+            result[(_at(variables, k), _at(binary_refs, k))] =
+                _at(disaggregateds, k)
+        end
+    end
+    return result
+end
+
+function DP.build_loa_problem(
+    model::InfiniteOpt.InfiniteModel,
+    method::DP.LOA
+    )
+    InfiniteOpt.build_transformation_backend!(model)
+    nlp = InfiniteOpt.transformation_model(model)
+
+    binaries = JuMP.VariableRef[]
+    for (_, binary_ref) in DP._indicator_to_binary(model)
+        append!(binaries, _transcribed_refs(
+            InfiniteOpt.transformation_variable(
+                DP._underlying_binary(binary_ref))))
+    end
+    unique!(binaries)
+
+    disjunct_constraints = Tuple{Union{JuMP.VariableRef, JuMP.AffExpr},
+        JuMP.AbstractJuMPScalar, _MOI.AbstractScalarSet}[]
+    for (_, disjunction) in DP._disjunctions(model)
+        for indicator in disjunction.constraint.indicators
+            haskey(DP._indicator_to_constraints(model), indicator) ||
+                continue
+            binary_refs = _transcribed_binary_refs(model, indicator)
+            for cref in DP._indicator_to_constraints(model)[indicator]
+                cref isa DP.DisjunctConstraintRef || continue
+                constraint = DP._disjunct_constraints(model)[
+                    JuMP.index(cref)].constraint
+                constraint isa JuMP.ScalarConstraint || continue
+                DP._is_linear_F(typeof(constraint.func)) && continue
+                funcs = _transcribed_refs(InfiniteOpt.transformation_expression(
+                    constraint.func))
+                for k in 1:max(length(binary_refs), length(funcs))
+                    push!(disjunct_constraints,
+                        (_at(binary_refs, k), _at(funcs, k), constraint.set))
+                end
+            end
+        end
+    end
+
+    global_constraints = Tuple{JuMP.AbstractJuMPScalar,
+        _MOI.AbstractScalarSet}[]
+    reform_set = Set(DP._reformulation_constraints(model))
+    for (F, S) in JuMP.list_of_constraint_types(model)
+        F === InfiniteOpt.GeneralVariableRef && continue
+        DP._is_linear_F(F) && continue
+        for cref in JuMP.all_constraints(model, F, S)
+            cref in reform_set && continue
+            constraint = JuMP.constraint_object(cref)
+            constraint isa JuMP.ScalarConstraint || continue
+            for func in _transcribed_refs(InfiniteOpt.transformation_expression(
+                constraint.func))
+                push!(global_constraints, (func, constraint.set))
+            end
+        end
+    end
+
+    return DP._LOAProblem(nlp, binaries,
+        disjunct_constraints, global_constraints,
+        _transcribed_disaggregation_map(model))
+end
+
+# Inject the outcome into the transcription backend and mark it ready so
+# InfiniteOpt value/objective queries return it. No binary restore (the
+# relaxed binaries are backend-transcribed and discarded on rebuild) and no
+# solver stash (a later reformulation mutates the InfiniteModel, invalidates
+# the backend, and rebuilds it fresh with the real optimizer).
+function DP._load_solution(
+    model::InfiniteOpt.InfiniteModel,
+    problem::DP._LOAProblem,
+    outcome::NamedTuple,
+    method::DP.LOA
+    )
+    DP._inject_solution(problem.nlp, problem, outcome)
+    InfiniteOpt.set_transformation_backend_ready(model, true)
+    return
+end
+
 end
diff --git a/src/DisjunctiveProgramming.jl b/src/DisjunctiveProgramming.jl
index 85f03455..6f1537e7 100644
--- a/src/DisjunctiveProgramming.jl
+++ b/src/DisjunctiveProgramming.jl
@@ -28,6 +28,7 @@ include("print.jl")
 include("extension_api.jl")
 include("utilities.jl")
 include("psplit.jl")
+include("loa.jl")
 
 # Define additional stuff that should not be exported
 const _EXCLUDE_SYMBOLS = [Symbol(@__MODULE__), :eval, :include]
diff --git a/src/cuttingplanes.jl b/src/cuttingplanes.jl
index dbc818e5..63692258 100644
--- a/src/cuttingplanes.jl
+++ b/src/cuttingplanes.jl
@@ -8,25 +8,17 @@ function collect_cutting_planes_vars(model::JuMP.AbstractModel)
     return collect_all_vars(model)
 end
 
-# Extract solution from a solved model (in-place). Extensions
-# override for models where values live on a backend.
+# Read primal values from a solved model. Returns
+# `Dict{var, Vector{value}}` — per-support shape uniformly: finite
+# models trivially have one "support" (length-1 Vector), the
+# InfiniteOpt extension overrides this dispatch to populate
+# multi-support Vectors. Skips fixed vars.
 function extract_solution(model::JuMP.AbstractModel)
     dvars = collect_cutting_planes_vars(model)
     V = eltype(dvars)
     T = JuMP.value_type(typeof(model))
     return Dict{V, Vector{T}}(
-        v => [JuMP.value(v)] for v in dvars)
-end
-
-# Extract solution from a GDPSubmodel (SEP path).
-function extract_solution(sub::GDPSubmodel)
-    V = eltype(sub.decision_vars)
-    T = JuMP.value_type(typeof(sub.model))
-    sol = Dict{V, Vector{T}}()
-    for var in sub.decision_vars
-        sol[var] = JuMP.value.(sub.fwd_map[var])
-    end
-    return sol
+        v => [JuMP.value(v)] for v in dvars if !JuMP.is_fixed(v))
 end
 
 # Set quadratic separation objective: min Σ (x_k - rBM_k)².
diff --git a/src/datatypes.jl b/src/datatypes.jl
index bdcd4c2e..8bbbce16 100644
--- a/src/datatypes.jl
+++ b/src/datatypes.jl
@@ -430,7 +430,7 @@ mutable struct _Hull{V <: JuMP.AbstractVariableRef, T} <: AbstractReformulationM
     function _Hull(method::Hull{T}, vrefs::Set{V}) where {T, V <: JuMP.AbstractVariableRef}
         new{V, T}(
             method.value,
-            Dict{V, Vector{V}}(vref => V[] for vref in vrefs), 
+            Dict{V, Vector{V}}(vref => V[] for vref in vrefs),
             Dict{Tuple{V, Union{V, JuMP.GenericAffExpr{T, V}}}, V}()
         )
     end
@@ -445,24 +445,20 @@ A type for using the cutting planes approach for disjunctive constraints.
 - `optimizer::O`: Optimizer to use when solving mini-models (required).
 - `max_iter::Int`: Number of iterations (default = `3`).
 - `seperation_tolerance::T`: Tolerance for the separation problem (default = `1e-6`).
-- `final_reform_method::AbstractReformulationMethod`: Final reformulation 
-method to use after cutting planes (default = `BigM()`).
 - `M_value::T`: Big-M value to use in the final reformulation (default = `1e9`).
 """
 struct CuttingPlanes{O, T} <: AbstractReformulationMethod
     optimizer::O;
     max_iter::Int
     seperation_tolerance::T
-    final_reform_method::AbstractReformulationMethod
     M_value::T
     function CuttingPlanes(
         optimizer::O;
         max_iter::Int = 3,
         seperation_tolerance::T = 1e-6,
-        final_reform_method = BigM(),
         M_value::T = 1e9
     ) where {O, T}
-        new{O, T}(optimizer, max_iter, seperation_tolerance, final_reform_method, M_value)
+        new{O, T}(optimizer, max_iter, seperation_tolerance, M_value)
     end
 end
 
@@ -473,25 +469,24 @@ end
 """
     GDPSubmodel{M, V, W}
 
-A unified submodel wrapper used by MBM and cutting plane
-reformulations. It encapsulates a flat JuMP optimization
-submodel built from a single disjunct's feasible region,
-along with mappings back to the original model's variables.
+A JuMP submodel paired with the parent model's decision variables
+and a map from each parent variable to its submodel representation.
 
 ## Fields
-- `model::M`: The JuMP submodel representing a disjunct's
-   feasible region (constraints and variable bounds).
-- `decision_vars::Vector{V}`: Ordered decision variables in
-   the submodel, matching the original model's ordering.
-- `fwd_map::Dict{V, Vector{W}}`: Forward map from original
-   model variables to their submodel counterparts.
+- `model::M`: The JuMP submodel.
+- `decision_vars::Vector{V}`: Parent-model decision variables in
+   the order used by the submodel.
+- `fwd_map::Dict{V, W}`: Map from each parent-model variable to
+   its submodel counterpart. `W` can be a single variable
+   reference (1:1 copy) or a collection of references (e.g. one
+   per transcription support).
 """
 struct GDPSubmodel{M <: JuMP.AbstractModel,
                    V <: JuMP.AbstractVariableRef,
-                   W <: JuMP.AbstractVariableRef}
+                   W}
     model::M
     decision_vars::Vector{V}
-    fwd_map::Dict{V, Vector{W}}
+    fwd_map::Dict{V, W}
 end
 
 """
@@ -567,6 +562,146 @@ A type for using indicator constraint approach for linear disjunctive constraint
 """
 struct Indicator <: AbstractReformulationMethod end
 
+################################################################################
+#                              LOA
+################################################################################
+"""
+    LOA{O, P, R, T} <: AbstractSolutionMethod
+
+Logic-based Outer Approximation solver for GDP models. Iterates a primary
+NLP (original model reformulated by `inner_method`, binaries fixed per
+iteration) and a master MILP accumulating OA and no-good cuts.
+`inner_method` is `BigM` (default), `MBM`, or `Hull`.
+
+LOA is a solution algorithm, not a reformulation: it is invoked with
+`optimize!(model, gdp_method = LOA(...))`, always re-runs on `optimize!`,
+and is not accepted by [`reformulate_model`](@ref). After the run the
+results are loaded into `model`, so the standard JuMP queries answer
+directly: `termination_status` reports the LOA-level outcome
+(`LOCALLY_SOLVED` on convergence or combination exhaustion, the binding
+limit otherwise), `objective_bound` and `relative_gap` report the master
+bound, `raw_status` carries the run summary, and every feasible
+combination visited is exposed as a solution-pool result
+(`value(x; result = k)`, best first). Duals are not available, matching
+other MI(N)LP solvers.
+
+## Fields
+- `nlp_optimizer::O`: solver for the primary NLP.
+- `mip_optimizer::P`: solver for the master MILP (default `nlp_optimizer`).
+- `inner_method::R`: NLP reformulation — `BigM`, `MBM`, or `Hull`.
+- `max_iter::Int`: max iterations after set-covering seeding.
+- `set_cover_max_iter::Int`: max iterations of the set-covering seed
+  loop.
+- `M_value::T`: big-M for the disjunct OA cut gating term.
+- `max_slack::T`: upper bound per slack variable.
+- `oa_penalty::T`: penalty on slacks in the master objective.
+- `use_nlpf::Bool`: solve the NLPF feasibility subproblem when the
+  primary NLP is infeasible (default `true`). When `false`, an
+  infeasible combination contributes only its no-good cut, no OA cuts.
+- `convergence_tol::Float64`: relative gap tolerance for the early stop.
+- `slack_tol::Float64`: max total slack for which the bound still counts
+  as converged (positive slack = nonconvex crossing, keep iterating).
+- `iteration_time_limit::Float64`: budget (s) for the iteration loop.
+- `time_limit::Float64`: overall budget (s) for the algorithm; the
+  final loading solve is not capped (default 3600; `Inf` disables).
+"""
+struct LOA{O, P, R, T} <: AbstractSolutionMethod
+    nlp_optimizer::O
+    mip_optimizer::P
+    inner_method::R
+    max_iter::Int
+    set_cover_max_iter::Int
+    M_value::T
+    max_slack::T
+    oa_penalty::T
+    use_nlpf::Bool
+    convergence_tol::Float64
+    slack_tol::Float64
+    iteration_time_limit::Float64
+    time_limit::Float64
+    function LOA(
+        nlp_optimizer::O;
+        mip_optimizer::P = nlp_optimizer,
+        max_iter::Int = 10,
+        set_cover_max_iter::Int = 8,
+        M_value::T = 1e9,
+        max_slack::T = 1e3,
+        oa_penalty::T = 1e3,
+        inner_method::R = BigM(M_value),
+        use_nlpf::Bool = true,
+        convergence_tol::Float64 = 1e-6,
+        slack_tol::Float64 = 1e-4,
+        iteration_time_limit::Float64 = Inf,
+        time_limit::Float64 = 3600.0
+        ) where {O, P, R <: AbstractReformulationMethod, T}
+        R <: Union{BigM, MBM, Hull} || error(
+            "LOA inner_method must be BigM, MBM, or Hull (got $R). " *
+            "PSplit is not yet supported.")
+        new{O, P, R, T}(nlp_optimizer, mip_optimizer, inner_method,
+            max_iter, set_cover_max_iter, M_value, max_slack, oa_penalty,
+            use_nlpf, convergence_tol, slack_tol,
+            iteration_time_limit, time_limit)
+    end
+end
+
+# The problem the LOA loop operates on: the model solved as the NLP
+# subproblem, the binary variables backing the indicators, the
+# nonlinear disjunct `(binary_ref, function, set)` triples (`binary_ref`
+# is the binary or its `1 - y` complement expression), the nonlinear
+# global `(function, set)` pairs, and the `(variable, binary_ref) ->
+# disaggregated variable` map from an inner Hull reformulation
+# (empty for Big-M / MBM). Built by `build_loa_problem`. The
+# set-covering seed generates its combinations on the fly from the
+# master, so none are stored here.
+struct _LOAProblem{M <: JuMP.AbstractModel, V <: JuMP.AbstractVariableRef, T}
+    nlp::M
+    binaries::Vector{V}
+    disjunct_constraints::Vector{Tuple{Union{V, JuMP.GenericAffExpr{T, V}},
+        JuMP.AbstractJuMPScalar, _MOI.AbstractScalarSet}}
+    global_constraints::Vector{Tuple{JuMP.AbstractJuMPScalar,
+        _MOI.AbstractScalarSet}}
+    disaggregation_map::Dict{Tuple{V, Union{V, JuMP.GenericAffExpr{T, V}}}, V}
+
+    # Inner constructor deriving the value type `T` from the NLP model.
+    # `T` only appears inside `Union` field types, so the default
+    # constructors could be called without binding it (Aqua
+    # unbound_args); this replaces them.
+    function _LOAProblem(
+        nlp::M,
+        binaries::Vector{V},
+        disjunct_constraints::Vector,
+        global_constraints::Vector,
+        disaggregation_map
+        ) where {M <: JuMP.AbstractModel, V <: JuMP.AbstractVariableRef}
+        return new{M, V, JuMP.value_type(M)}(nlp, binaries,
+            disjunct_constraints, global_constraints, disaggregation_map)
+    end
+end
+
+# The LOA master MILP plus the NLP-to-master variable map.
+# `objective` is the NLP's objective, linearized per iteration for
+# the objective cuts; `disaggregator` is the master `_Hull` for Hull
+# cuts, `nothing` for Big-M / MBM.
+struct _LOAMaster{M <: JuMP.AbstractModel, VM, OF, AO, DG}
+    model::M
+    variable_map::VM
+    objective_sense::_MOI.OptimizationSense
+    objective::OF
+    alpha_oa::AO
+    disaggregator::DG
+end
+
+# The optimizer left in a model's backend slot after an LOA solve: a
+# `MockOptimizer` preloaded with the run's results so every standard
+# JuMP query answers with no solve, plus the composed solver name (the
+# mock's own is hardcoded to "Mock"). `model.jl` detects this type to
+# swap the real solver back before a later reformulation solve.
+struct _LOAResultCache{O <: _MOI.Utilities.MockOptimizer} <:
+        _MOI.AbstractOptimizer
+    mock::O
+    name::String
+end
+
 ################################################################################
 #                              GDP Data
 ################################################################################
@@ -593,6 +728,10 @@ mutable struct GDPData{M <: JuMP.AbstractModel, V <: JuMP.AbstractVariableRef, C
     # Helpful metadata for most reformulations (not just one of them)
     variable_bounds::Dict{V, Tuple{T, T}}
 
+    # Hull disaggregations recorded during reformulation, consumed by
+    # solution algorithms (LOA): (variable, indicator) => disaggregated
+    disaggregation_map::Dict{Tuple{V, LogicalVariableRef{M}}, V}
+
     # Reformulation variables and constraints
     reformulation_variables::Vector{V}
     reformulation_constraints::Vector{C}
@@ -600,6 +739,9 @@ mutable struct GDPData{M <: JuMP.AbstractModel, V <: JuMP.AbstractVariableRef, C
     # Solution data
     solution_method::Union{Nothing, AbstractSolutionMethod}
     ready_to_optimize::Bool
+    # The optimizer an LOA solve displaced (it solves the NLP in place),
+    # stashed so a later reformulation re-solve can put it back.
+    displaced_optimizer::Any
 
     # Default constructor
     function GDPData{M, V, C}() where {M <: JuMP.AbstractModel, V <: JuMP.AbstractVariableRef, C}
@@ -613,10 +755,12 @@ mutable struct GDPData{M <: JuMP.AbstractModel, V <: JuMP.AbstractVariableRef, C
             Dict{LogicalVariableRef{M}, Vector{Union{DisjunctConstraintRef{M}, DisjunctionRef{M}}}}(),
             Dict{Union{DisjunctConstraintRef{M}, DisjunctionRef{M}}, LogicalVariableRef{M}}(),
             Dict{V, Tuple{T, T}}(),
+            Dict{Tuple{V, LogicalVariableRef{M}}, V}(),
             Vector{V}(),
             Vector{C}(),
             nothing,
             false,
+            nothing,
         )
     end
 end
diff --git a/src/hull.jl b/src/hull.jl
index 779369e9..8824dc0f 100644
--- a/src/hull.jl
+++ b/src/hull.jl
@@ -42,6 +42,8 @@ function _disaggregate_variable(
     #temp storage
     push!(method.disjunction_variables[vref], dvref)
     method.disjunct_variables[vref, bvref] = dvref
+    #record the disaggregation for solution algorithms (LOA)
+    _disaggregation_map(model)[(vref, lvref)] = dvref
     #create bounding constraints
     dvname = JuMP.name(dvref)
     lbname = isempty(dvname) ? "" : "$(dvname)_lower_bound"
diff --git a/src/loa.jl b/src/loa.jl
new file mode 100644
index 00000000..4465e359
--- /dev/null
+++ b/src/loa.jl
@@ -0,0 +1,1026 @@
+################################################################################
+#                    LOGIC-BASED OUTER APPROXIMATION (LOA)
+################################################################################
+
+################################################################################
+#                            SENSE HANDLERS
+################################################################################
+_penalty_sign(::Val{_MOI.MIN_SENSE}) = 1
+_penalty_sign(::Val{_MOI.MAX_SENSE}) = -1
+_worst_objective(::Val{_MOI.MIN_SENSE}) = Inf
+_worst_objective(::Val{_MOI.MAX_SENSE}) = -Inf
+_is_better(::Val{_MOI.MIN_SENSE}, new, best) = new < best
+_is_better(::Val{_MOI.MAX_SENSE}, new, best) = new > best
+_gap(::Val{_MOI.MIN_SENSE}, best, bound) = best - bound
+_gap(::Val{_MOI.MAX_SENSE}, best, bound) = bound - best
+# A valid overall bound is the weaker of two partial bounds.
+_loosest_bound(::Val{_MOI.MIN_SENSE}, best, bound) = min(best, bound)
+_loosest_bound(::Val{_MOI.MAX_SENSE}, best, bound) = max(best, bound)
+
+# One progress record per NLP solve, for external convergence traces.
+function _log_loa_progress(t_start::Float64, best_objective, master_bound)
+    bound = master_bound === nothing ? NaN : master_bound
+    @info "LOA progress: elapsed=$(time() - t_start) " *
+        "incumbent=$best_objective bound=$bound"
+    return
+end
+
+################################################################################
+#                            MAIN LOOP
+################################################################################
+# LOA optimize hook: build the problem, seed with set-covering
+# combinations, iterate the master/NLP loop, then load the results into
+# the model by injecting every feasible combination found (best first).
+function _optimize_hook(model::JuMP.AbstractModel, method::LOA; kwargs...)
+    _clear_reformulations(model)
+    reformulate_model(model, method.inner_method)
+
+    problem = build_loa_problem(model, method)
+    # The master copies the NLP while its binaries still carry
+    # integrality; the NLP is then relaxed once and each iteration
+    # overwrites the binary fixes in place.
+    master = _build_loa_master(problem, method)
+    _relax_binaries(problem)
+    nlp = problem.nlp
+    displaced_optimizer = _install_nlp_optimizer(nlp, method.nlp_optimizer,
+        gdp_data(model).displaced_optimizer)
+    JuMP.set_silent(nlp)
+    t_start = time()
+    overall_deadline = t_start + method.time_limit
+    loop_deadline =
+        min(t_start + method.iteration_time_limit, overall_deadline)
+    original_time_limit = JuMP.time_limit_sec(nlp)
+    sense_token = Val(master.objective_sense)
+    best_objective = _worst_objective(sense_token)
+    best_result = nothing
+    master_bound = nothing
+    feasible_results = NamedTuple[]
+    master_failure = nothing
+    cuts_added = false
+    converged = false
+
+    # Set-covering seed. This mimics Pyomo GDPopt's set-covering
+    # initialization: borrow the master as a covering MILP (only its
+    # objective changes), let it pick a combination that activates the
+    # nonlinear disjuncts still lacking a linearization, solve the NLP
+    # there, and seed the resulting OA and no-good cuts. Each NLP
+    # warm-starts from the last feasible primal.
+    previous_result = nothing
+    cover_disjuncts = _cover_disjuncts(problem)
+    needs_cover = trues(length(cover_disjuncts))
+    num_covered = 0
+    for iteration in 1:method.set_cover_max_iter
+        (iteration == 1 || any(needs_cover)) || break
+        time() < loop_deadline || break
+        # Swap in the covering objective, solve, read off the combination,
+        # then restore the OA objective (with its accumulated slack
+        # penalties) before emitting cuts against it.
+        oa_objective = JuMP.objective_function(master.model)
+        JuMP.@objective(master.model, Max,
+            _cover_objective(master, cover_disjuncts, needs_cover,
+                num_covered))
+        _cap_remaining_time(master.model, loop_deadline)
+        JuMP.optimize!(master.model)
+        solved = JuMP.is_solved_and_feasible(master.model)
+        combination = solved ? _extract_combination(problem, master) : nothing
+        JuMP.set_objective_sense(master.model, master.objective_sense)
+        JuMP.set_objective_function(master.model, oa_objective)
+        if !solved
+            master_failure = JuMP.termination_status(master.model)
+            break
+        end
+        _set_nlp_warm_start(previous_result)
+        result = _solve_nlp(problem, combination, method;
+            deadline = loop_deadline)
+        avoid_combination(master.model, combination, master.variable_map)
+        _add_oa_cuts(problem, master, result, method)
+        cuts_added = true
+        if result.feasible
+            push!(feasible_results, result)
+            previous_result = result
+            if _is_better(sense_token, result.objective, best_objective)
+                best_objective = result.objective
+                best_result = result
+            end
+        end
+        _log_loa_progress(t_start, best_objective, master_bound)
+        # A disjunct counts as covered only once it is active in a
+        # feasible NLP. An infeasible combination still contributes its
+        # no-good cut above but leaves the coverage targets untouched.
+        if result.feasible
+            for i in eachindex(cover_disjuncts)
+                needs_cover[i] || continue
+                _disjunct_active(result.combination, cover_disjuncts[i]) &&
+                    (needs_cover[i] = false)
+            end
+            num_covered = count(!, needs_cover)
+        end
+    end
+
+    # Master/NLP loop: `alpha_oa` is the bound, the NLP refines the
+    # incumbent. Exit on convergence, infeasible master, max_iter, or time.
+    for _ in 1:method.max_iter
+        time() < loop_deadline || break
+        _cap_remaining_time(master.model, loop_deadline)
+        JuMP.optimize!(master.model)
+        if !JuMP.is_solved_and_feasible(master.model)
+            master_failure = JuMP.termination_status(master.model)
+            break
+        end
+        master_bound = JuMP.value(master.alpha_oa)
+        if best_result !== nothing
+            gap = _gap(sense_token, best_objective, master_bound)
+            total_slack = abs(JuMP.objective_value(master.model) -
+                master_bound) / method.oa_penalty
+            tol = method.convergence_tol * max(abs(best_objective), 1.0)
+            if gap <= tol && total_slack <= method.slack_tol
+                converged = true
+                break
+            end
+        end
+        combination = _extract_combination(problem, master)
+        _set_nlp_warm_start(previous_result)
+        result = _solve_nlp(problem, combination, method;
+            deadline = loop_deadline)
+        avoid_combination(master.model, combination, master.variable_map)
+        _add_oa_cuts(problem, master, result, method)
+        cuts_added = true
+        if result.feasible
+            push!(feasible_results, result)
+            previous_result = result
+            if _is_better(sense_token, result.objective, best_objective)
+                best_objective = result.objective
+                best_result = result
+            end
+        end
+        _log_loa_progress(t_start, best_objective, master_bound)
+    end
+
+    _restore_time_limit(nlp, original_time_limit)
+    status = _loa_termination_status(best_result !== nothing, converged,
+        master_failure, cuts_added, time() >= loop_deadline)
+    message = _loa_status_message(status, best_objective, best_result,
+        master_bound, sense_token)
+    best_result === nothing ? (@warn message) : (@info message)
+    # Must be read before the injection replaces the NLP's optimizer.
+    solver = "LOA(" * string(nameof(typeof(method.inner_method))) *
+        ", nlp = " * JuMP.solver_name(nlp) *
+        ", master = " * JuMP.solver_name(master.model) * ")"
+    sort!(feasible_results; by = result -> result.objective,
+        rev = master.objective_sense == _MOI.MAX_SENSE)
+    # No-good cuts mean the master bound only covers unvisited
+    # combinations, so the reported bound folds in the incumbent
+    # (`raw_status` above keeps the unfolded master bound).
+    bound = if best_result === nothing
+        master_bound
+    elseif _infeasible_master(master_failure)
+        best_objective
+    elseif master_bound === nothing
+        nothing
+    else
+        _loosest_bound(sense_token, best_objective, master_bound)
+    end
+    outcome = (results = feasible_results, status = status,
+        raw_status = message, bound = bound,
+        solve_time = time() - t_start, solver_name = solver,
+        displaced_optimizer = displaced_optimizer)
+    # Injection is not a reformulation; a later method must rebuild.
+    _set_solution_method(model, method)
+    _set_ready_to_optimize(model, false)
+    _load_solution(model, problem, outcome, method)
+    return
+end
+
+# Master statuses that report infeasibility. A local claim
+# (LOCALLY_INFEASIBLE, e.g. from an NLP-based MILP solver like Juniper)
+# counts for combination exhaustion, but only a global certificate
+# justifies the global INFEASIBLE report.
+_infeasible_master(::Nothing) = false
+_infeasible_master(status::_MOI.TerminationStatusCode) = status in
+    (_MOI.INFEASIBLE, _MOI.INFEASIBLE_OR_UNBOUNDED,
+        _MOI.LOCALLY_INFEASIBLE)
+_proven_infeasible_master(::Nothing) = false
+_proven_infeasible_master(status::_MOI.TerminationStatusCode) = status in
+    (_MOI.INFEASIBLE, _MOI.INFEASIBLE_OR_UNBOUNDED)
+
+# LOA-level termination status, reported the way a native solver would:
+# convergence and combination exhaustion (master infeasible with an
+# incumbent in hand) are locally solved (a local NLP solver cannot
+# certify global optimality), a master proven infeasible before any cut
+# proves the linear relaxation infeasible, exhaustion without an
+# incumbent is only a local claim, and otherwise the binding limit is
+# reported (`OTHER_LIMIT` for a master solve that failed without an
+# infeasibility status, e.g. on its own time cap).
+function _loa_termination_status(
+    found_incumbent::Bool,
+    converged::Bool,
+    master_failure::Union{Nothing, _MOI.TerminationStatusCode},
+    cuts_added::Bool,
+    timed_out::Bool
+    )
+    (converged || (_infeasible_master(master_failure) && found_incumbent)) &&
+        return _MOI.LOCALLY_SOLVED
+    _proven_infeasible_master(master_failure) && !cuts_added &&
+        return _MOI.INFEASIBLE
+    _infeasible_master(master_failure) && return _MOI.LOCALLY_INFEASIBLE
+    timed_out && return _MOI.TIME_LIMIT
+    master_failure === nothing || return _MOI.OTHER_LIMIT
+    return _MOI.ITERATION_LIMIT
+end
+
+# Human-readable run summary: logged at the end of the hook and stored
+# as the model's `RawStatusString`. The bound is rigorous only for a
+# convex inner problem; on a nonconvex one it can cross the incumbent
+# (negative gap).
+function _loa_status_message(
+    status::_MOI.TerminationStatusCode,
+    best_objective::Real,
+    best_result,
+    master_bound,
+    sense_token::Val
+    )
+    best_result === nothing &&
+        return "LOA finished [$status]: no feasible incumbent found."
+    master_bound === nothing &&
+        return "LOA finished [$status]: incumbent $best_objective " *
+            "(master produced no bound; best seed kept)."
+    gap = _gap(sense_token, best_objective, master_bound)
+    relative = abs(best_objective) > 1e-10 ?
+        gap / abs(best_objective) : gap
+    label = status == _MOI.LOCALLY_SOLVED ? "converged" : "limit hit"
+    return "LOA finished [$label]: incumbent $best_objective, master " *
+        "bound $master_bound, gap $gap (relative $relative)."
+end
+
+################################################################################
+#                       SET-COVERING INITIALIZATION
+################################################################################
+# The nonlinear disjuncts to cover: one entry per distinct indicator that
+# owns a nonlinear disjunct constraint, carrying that indicator's
+# `binary_ref`. Keyed by `(underlying binary, active value)` so the two
+# disjuncts of a single-binary disjunction (`y` and its `1 - y`
+# complement) stay distinct. Purely linear disjuncts need no cover: the
+# inner reformulation already places them in the master exactly.
+function _cover_disjuncts(problem::_LOAProblem{M, V, T}) where {M, V, T}
+    seen = Set{Tuple{V, Bool}}()
+    disjuncts = Union{V, JuMP.GenericAffExpr{T, V}}[]
+    for (binary_ref, _, _) in problem.disjunct_constraints
+        key = (_underlying_binary(binary_ref),
+            _underlying_value(binary_ref, true))
+        key in seen && continue
+        push!(seen, key)
+        push!(disjuncts, binary_ref)
+    end
+    return disjuncts
+end
+
+# The set-covering objective (master space), mimicking Pyomo GDPopt:
+# maximize active disjuncts weighted `num_covered + 1` if still uncovered
+# else `1`, so one uncovered disjunct outweighs every covered one and each
+# solve must activate a new disjunct when the logic allows. Empty
+# `disjuncts` gives the zero expression: a constant objective that still
+# seeds one feasible combination (needed to bound `alpha_oa`).
+function _cover_objective(
+    master::_LOAMaster,
+    disjuncts,
+    needs_cover,
+    num_covered::Int
+    )
+    T = JuMP.value_type(typeof(master.model))
+    V = JuMP.variable_ref_type(typeof(master.model))
+    expr = JuMP.GenericAffExpr{T, V}(zero(T))
+    for i in eachindex(disjuncts)
+        weight = needs_cover[i] ? num_covered + 1 : 1
+        activation = _remap_indicator_to_binary(disjuncts[i],
+            master.variable_map)
+        JuMP.add_to_expression!(expr, T(weight), activation)
+    end
+    return expr
+end
+
+################################################################################
+#                       FINITE PROBLEM CONSTRUCTION
+################################################################################
+# The concrete variable behind an indicator's binary reference (the
+# underlying variable of a `1 - y` complement expression), and the value
+# it takes when the indicator is active.
+_underlying_binary(binary_ref::JuMP.AbstractVariableRef) = binary_ref
+_underlying_binary(binary_ref::JuMP.GenericAffExpr) =
+    only(keys(binary_ref.terms))
+_underlying_value(::JuMP.AbstractVariableRef, active::Bool) = active
+_underlying_value(::JuMP.GenericAffExpr, active::Bool) = !active
+
+"""
+    build_loa_problem(
+        model::JuMP.AbstractModel,
+        method::LOA
+        )::_LOAProblem
+
+Build the problem the LOA loop operates on from `model` (already
+reformulated by the LOA inner method): the NLP subproblem, the binary
+variables backing the indicators, the nonlinear disjunct
+`(binary_ref, function, set)` triples, the nonlinear global
+`(function, set)` pairs, and the Hull disaggregation map recorded in
+the GDP data during reformulation. The NLP is `model` itself; the
+InfiniteOpt extension overloads this to build the problem from the
+transcribed backend instead.
+
+## Returns
+- `_LOAProblem`: the problem.
+"""
+function build_loa_problem(model::JuMP.AbstractModel, method::LOA)
+    V = JuMP.variable_ref_type(typeof(model))
+    T = JuMP.value_type(typeof(model))
+    binary_map = _indicator_to_binary(model)
+
+    binaries = V[_underlying_binary(binary_ref)
+        for (_, binary_ref) in binary_map]
+    unique!(binaries)
+
+    disjunct_constraints = Tuple{Union{V, JuMP.GenericAffExpr{T, V}},
+        JuMP.AbstractJuMPScalar, _MOI.AbstractScalarSet}[]
+    for (_, disjunction) in _disjunctions(model)
+        for indicator in disjunction.constraint.indicators
+            haskey(_indicator_to_constraints(model), indicator) || continue
+            binary_ref = binary_map[indicator]
+            for cref in _indicator_to_constraints(model)[indicator]
+                cref isa DisjunctConstraintRef || continue
+                constraint = _disjunct_constraints(model)[
+                    JuMP.index(cref)].constraint
+                constraint isa JuMP.ScalarConstraint || continue
+                _is_linear_F(typeof(constraint.func)) && continue
+                push!(disjunct_constraints,
+                    (binary_ref, constraint.func, constraint.set))
+            end
+        end
+    end
+
+    global_constraints = Tuple{JuMP.AbstractJuMPScalar,
+        _MOI.AbstractScalarSet}[]
+    reform_set = Set(_reformulation_constraints(model))
+    for (F, S) in JuMP.list_of_constraint_types(model)
+        F === V && continue
+        _is_linear_F(F) && continue
+        for cref in JuMP.all_constraints(model, F, S)
+            cref in reform_set && continue
+            constraint = JuMP.constraint_object(cref)
+            constraint isa JuMP.ScalarConstraint || continue
+            push!(global_constraints, (constraint.func, constraint.set))
+        end
+    end
+
+    binary_disaggregations =
+        Dict{Tuple{V, Union{V, JuMP.GenericAffExpr{T, V}}}, V}(
+            (variable, binary_map[indicator]) => disaggregated
+            for ((variable, indicator), disaggregated)
+                in _disaggregation_map(model))
+    return _LOAProblem(model, binaries, disjunct_constraints,
+        global_constraints, binary_disaggregations)
+end
+
+# Relax the binaries to continuous in [0, 1]; each NLP solve then
+# overwrites their fixes in place.
+function _relax_binaries(problem::_LOAProblem)
+    for binary in problem.binaries
+        JuMP.unset_binary(binary)
+        JuMP.set_lower_bound(binary, 0.0)
+        JuMP.set_upper_bound(binary, 1.0)
+    end
+    return
+end
+
+# Undo `_relax_binaries` and the loop's last combination fix: drop the
+# fix and the relaxation bounds, then restore integrality so the model
+# is handed back with its binaries exactly as the reformulation created
+# them. (A fixed binary carries no bounds — `fix(; force = true)`
+# deleted them — while an unfixed one still carries the relaxation's.)
+function _restore_binaries(problem::_LOAProblem)
+    for binary in problem.binaries
+        JuMP.is_fixed(binary) && JuMP.unfix(binary)
+        JuMP.has_lower_bound(binary) && JuMP.delete_lower_bound(binary)
+        JuMP.has_upper_bound(binary) && JuMP.delete_upper_bound(binary)
+        JuMP.set_binary(binary)
+    end
+    return
+end
+
+################################################################################
+#                         OPTIMIZER HANDOVER
+################################################################################
+# Install the NLP solver, returning the optimizer it displaces. LOA
+# solves the NLP on the model itself, so that one is the user's own. A
+# re-solve displaces the last run's cache instead, so `previous` (the
+# optimizer already stashed) carries through untouched.
+function _install_nlp_optimizer(
+    nlp::JuMP.AbstractModel,
+    nlp_optimizer,
+    previous
+    )
+    backend = JuMP.backend(nlp)
+    displaced = if backend isa _MOI.Utilities.CachingOptimizer &&
+            backend.optimizer !== nothing
+        JuMP.unsafe_backend(nlp)
+    else
+        nothing
+    end
+    JuMP.set_optimizer(nlp, nlp_optimizer)
+    return displaced isa _LOAResultCache ? previous : displaced
+end
+
+# Model types without a MOI backend of their own never carry a cache.
+function _has_loa_cache(model::JuMP.GenericModel)
+    backend = JuMP.backend(model)
+    return backend isa _MOI.Utilities.CachingOptimizer &&
+        backend.optimizer isa _LOAResultCache
+end
+_has_loa_cache(::JuMP.AbstractModel) = false
+
+# Restore the displaced optimizer, emptied so the re-solve copies the
+# model in. Without one, drop the cache rather than serve its results.
+function _restore_displaced_optimizer(model::JuMP.AbstractModel)
+    displaced = gdp_data(model).displaced_optimizer
+    if displaced === nothing
+        _MOI.Utilities.drop_optimizer(JuMP.backend(model))
+    else
+        _MOI.empty!(displaced)
+        JuMP.set_optimizer(model, () -> displaced)
+    end
+    return
+end
+
+################################################################################
+#                         MASTER CONSTRUCTION
+################################################################################
+# True for linear constraint function types (variable refs / affine).
+# Nonlinear functions enter the master as OA cuts, not at copy time.
+_is_linear_F(::Type{<:JuMP.AbstractVariableRef}) = true
+_is_linear_F(::Type{<:JuMP.GenericAffExpr}) = true
+_is_linear_F(::Type{<:AbstractVector{<:JuMP.AbstractVariableRef}}) = true
+_is_linear_F(::Type{<:AbstractVector{<:JuMP.GenericAffExpr}}) = true
+_is_linear_F(::Type) = false
+
+# Build the LOA master MILP: copy the NLP's variables and linear
+# constraints, install the `alpha_oa` objective auxiliary, and record
+# the NLP-to-master variable map. Nonlinear constraints are not
+# copied; they enter later as OA cuts per NLP solve.
+function _build_loa_master(problem::_LOAProblem, method::LOA)
+    nlp = problem.nlp
+    objective = JuMP.objective_function(nlp)
+    objective_sense = JuMP.objective_sense(nlp)
+    variable_type = JuMP.variable_ref_type(typeof(nlp))
+
+    master = _copy_model(nlp)
+    JuMP.set_optimizer(master, method.mip_optimizer)
+    JuMP.set_silent(master)
+
+    variable_map = copy_variables_onto_model(master, nlp)
+
+    for (F, S) in JuMP.list_of_constraint_types(nlp)
+        F === variable_type && continue
+        _is_linear_F(F) || continue
+        for constraint_ref in JuMP.all_constraints(nlp, F, S)
+            constraint = JuMP.constraint_object(constraint_ref)
+            new_func = replace_variables_in_constraint(
+                constraint.func, variable_map)
+            JuMP.@constraint(master, new_func in constraint.set)
+        end
+    end
+
+    alpha_oa = JuMP.@variable(master, base_name = "alpha_oa")
+    JuMP.@objective(master, objective_sense, alpha_oa)
+
+    disaggregator = _build_disaggregator(problem, variable_map,
+        method.inner_method)
+    return _LOAMaster(master, variable_map, objective_sense, objective,
+        alpha_oa, disaggregator)
+end
+
+# Master-space Hull disaggregator for the convex-hull cut emitter
+# (`nothing` for Big-M / MBM): remap the `(variable, binary_ref) ->
+# disaggregated variable` map into master refs.
+_build_disaggregator(::_LOAProblem, ::Dict, ::Union{BigM, MBM}) = nothing
+function _build_disaggregator(
+    problem::_LOAProblem,
+    variable_map::Dict{V, V},
+    inner::Hull
+    ) where {V <: JuMP.AbstractVariableRef}
+    hull = _Hull(inner, Set{valtype(variable_map)}())
+    for ((variable, binary_ref), disaggregated) in problem.disaggregation_map
+        hull.disjunct_variables[(variable_map[variable],
+            _remap_indicator_to_binary(binary_ref, variable_map))] =
+            variable_map[disaggregated]
+    end
+    return hull
+end
+
+################################################################################
+#                          NLP SUBPROBLEM
+################################################################################
+function _cap_remaining_time(target::JuMP.AbstractModel, deadline::Float64)
+    isfinite(deadline) || return
+    JuMP.set_time_limit_sec(target, max(0.0, deadline - time()))
+    return
+end
+_restore_time_limit(model::JuMP.AbstractModel, ::Nothing) =
+    JuMP.unset_time_limit_sec(model)
+_restore_time_limit(model::JuMP.AbstractModel, seconds::Real) =
+    JuMP.set_time_limit_sec(model, seconds)
+
+# Fix each binary at its combination value, overwriting any previous
+# fix in place.
+function _fix_combination(combination::AbstractDict)
+    for (binary, value) in combination
+        JuMP.fix(binary, value ? 1.0 : 0.0; force = true)
+    end
+    return
+end
+
+# Primal values keyed by the NLP's variables, skipping fixed ones
+# (the combination binaries and any user-fixed variables).
+function _extract_solution(nlp::JuMP.AbstractModel)
+    V = JuMP.variable_ref_type(typeof(nlp))
+    T = JuMP.value_type(typeof(nlp))
+    return Dict{V, T}(v => JuMP.value(v)
+        for v in JuMP.all_variables(nlp) if !JuMP.is_fixed(v))
+end
+
+# Solve the NLP at a fixed combination. If infeasible, fall
+# through to NLPF (a slacked version that always solves) so the master
+# still gets a linearization site, not just a no-good cut. NLPF can be
+# switched off via `use_nlpf = false`.
+function _solve_nlp(
+    problem::_LOAProblem,
+    combination,
+    method::LOA;
+    deadline::Float64 = Inf
+    )
+    nlp = problem.nlp
+    _cap_remaining_time(nlp, deadline)
+    _fix_combination(combination)
+    JuMP.optimize!(nlp, ignore_optimize_hook = true)
+    if JuMP.is_solved_and_feasible(nlp)
+        return (combination = combination,
+            linearization_point = _extract_solution(nlp),
+            objective = JuMP.objective_value(nlp), feasible = true)
+    end
+
+    # Primary NLP infeasible — try NLPF (unless disabled, in which case
+    # the combination contributes only its no-good cut).
+    if method.use_nlpf
+        nlpf = _solve_nlpf(problem, combination, method; deadline = deadline)
+        nlpf === nothing || return nlpf
+    end
+
+    V = JuMP.variable_ref_type(typeof(nlp))
+    T = JuMP.value_type(typeof(nlp))
+    return (combination = combination,
+        linearization_point = Dict{V, T}(),
+        objective = Inf, feasible = false)
+end
+
+################################################################################
+#                       NLPF (FEASIBILITY SUBPROBLEM)
+################################################################################
+# When the primary NLP is infeasible, copy the NLP (fixes included),
+# slack every scalar inequality with one nonnegative `u`, minimize `u`,
+# and return the primal as the linearization site (equalities/bounds
+# exact, inequalities violated by at most `u`).
+function _solve_nlpf(
+    problem::_LOAProblem,
+    combination,
+    method::LOA;
+    deadline::Float64 = Inf
+    )
+    copy, ref_map = JuMP.copy_model(problem.nlp)
+    JuMP.set_optimizer(copy, method.nlp_optimizer)
+    JuMP.set_silent(copy)
+    _cap_remaining_time(copy, deadline)
+
+    var_type = JuMP.variable_ref_type(typeof(copy))
+    u = JuMP.@variable(copy, lower_bound = 0.0, base_name = "_nlpf_u")
+
+    for (F, S) in JuMP.list_of_constraint_types(copy)
+        F === var_type && continue
+        _nlpf_should_slack(S) || continue
+        for cref in JuMP.all_constraints(copy, F, S)
+            con = JuMP.constraint_object(cref)
+            JuMP.delete(copy, cref)
+            new_func = _nlpf_slacked_func(con.func, u, con.set)
+            JuMP.@constraint(copy, new_func in con.set)
+        end
+    end
+
+    JuMP.@objective(copy, Min, u)
+
+    JuMP.optimize!(copy, ignore_optimize_hook = true)
+    # Use the primal only at a genuine feasible point; a solver can report
+    # `has_values` with a nonfeasible/NaN primal that would poison the cut.
+    JuMP.is_solved_and_feasible(copy) || return nothing
+
+    V = JuMP.variable_ref_type(typeof(problem.nlp))
+    T = JuMP.value_type(typeof(problem.nlp))
+    linearization_point = Dict{V, T}(v => JuMP.value(ref_map[v])
+        for v in JuMP.all_variables(problem.nlp) if !JuMP.is_fixed(v))
+    return (combination = combination,
+        linearization_point = linearization_point,
+        objective = Inf, feasible = false)
+end
+
+_nlpf_should_slack(::Type{<:_MOI.LessThan}) = true
+_nlpf_should_slack(::Type{<:_MOI.GreaterThan}) = true
+_nlpf_should_slack(::Type) = false
+
+_nlpf_slacked_func(func, u, ::_MOI.LessThan) = func - u
+_nlpf_slacked_func(func, u, ::_MOI.GreaterThan) = func + u
+
+# Iter-to-iter NLP warm start: seed the next NLP from the last FEASIBLE
+# primal (no-op on the first seed and after an NLPF fall-through, whose
+# primal is slack-distorted).
+function _set_nlp_warm_start(previous)
+    previous === nothing && return
+    for (variable, value) in previous.linearization_point
+        JuMP.set_start_value(variable, value)
+    end
+    return
+end
+
+################################################################################
+#                       LOA RESULT CACHE (MOI LAYER)
+################################################################################
+# `_LOAResultCache` delegates the MOI interface to its inner mock. Only
+# `SolverName` (the mock hardcodes "Mock"), the bound attributes (a
+# raw `KeyError` from the mock's generic storage when the master
+# produced no bound), and `optimize!` are intercepted.
+_MOI.is_empty(cache::_LOAResultCache) = _MOI.is_empty(cache.mock)
+_MOI.empty!(cache::_LOAResultCache) = _MOI.empty!(cache.mock)
+# "Solving" the cache would clear JuMP's dirty flag and serve the stored
+# results again, so an edited model would report a stale point as fresh.
+function _MOI.optimize!(cache::_LOAResultCache)
+    return error("This model holds the results of an LOA solve and " *
+        "cannot be re-solved directly (`ignore_optimize_hook = true`). " *
+        "Call `optimize!(model, gdp_method = ...)` to rebuild and solve.")
+end
+_MOI.copy_to(cache::_LOAResultCache, src::_MOI.ModelLike) =
+    _MOI.copy_to(cache.mock, src)
+_MOI.add_variable(cache::_LOAResultCache) = _MOI.add_variable(cache.mock)
+function _MOI.add_constraint(
+    cache::_LOAResultCache,
+    func::_MOI.AbstractFunction,
+    set::_MOI.AbstractSet
+    )
+    return _MOI.add_constraint(cache.mock, func, set)
+end
+_MOI.delete(cache::_LOAResultCache, index::_MOI.Index) =
+    _MOI.delete(cache.mock, index)
+_MOI.is_valid(cache::_LOAResultCache, index::_MOI.Index) =
+    _MOI.is_valid(cache.mock, index)
+function _MOI.modify(
+    cache::_LOAResultCache,
+    ci::_MOI.ConstraintIndex,
+    change::_MOI.AbstractFunctionModification
+    )
+    return _MOI.modify(cache.mock, ci, change)
+end
+function _MOI.modify(
+    cache::_LOAResultCache,
+    attr::_MOI.ObjectiveFunction,
+    change::_MOI.AbstractFunctionModification
+    )
+    return _MOI.modify(cache.mock, attr, change)
+end
+function _MOI.supports_constraint(
+    cache::_LOAResultCache,
+    F::Type{<:_MOI.AbstractFunction},
+    S::Type{<:_MOI.AbstractSet}
+    )
+    return _MOI.supports_constraint(cache.mock, F, S)
+end
+function _MOI.supports(
+    cache::_LOAResultCache,
+    attr::Union{_MOI.AbstractModelAttribute,
+        _MOI.AbstractOptimizerAttribute}
+    )
+    return _MOI.supports(cache.mock, attr)
+end
+function _MOI.supports(
+    cache::_LOAResultCache,
+    attr::Union{_MOI.AbstractVariableAttribute,
+        _MOI.AbstractConstraintAttribute},
+    IdxType::Type{<:_MOI.Index}
+    )
+    return _MOI.supports(cache.mock, attr, IdxType)
+end
+function _MOI.get(
+    cache::_LOAResultCache,
+    attr::Union{_MOI.AbstractModelAttribute,
+        _MOI.AbstractOptimizerAttribute}
+    )
+    return _MOI.get(cache.mock, attr)
+end
+function _MOI.get(
+    cache::_LOAResultCache,
+    attr::_MOI.AbstractVariableAttribute,
+    index::_MOI.VariableIndex
+    )
+    return _MOI.get(cache.mock, attr, index)
+end
+function _MOI.get(
+    cache::_LOAResultCache,
+    attr::_MOI.AbstractConstraintAttribute,
+    index::_MOI.ConstraintIndex
+    )
+    return _MOI.get(cache.mock, attr, index)
+end
+function _MOI.set(
+    cache::_LOAResultCache,
+    attr::Union{_MOI.AbstractModelAttribute,
+        _MOI.AbstractOptimizerAttribute},
+    value
+    )
+    return _MOI.set(cache.mock, attr, value)
+end
+function _MOI.set(
+    cache::_LOAResultCache,
+    attr::_MOI.AbstractVariableAttribute,
+    index::_MOI.VariableIndex,
+    value
+    )
+    return _MOI.set(cache.mock, attr, index, value)
+end
+function _MOI.set(
+    cache::_LOAResultCache,
+    attr::_MOI.AbstractConstraintAttribute,
+    index::_MOI.ConstraintIndex,
+    value
+    )
+    return _MOI.set(cache.mock, attr, index, value)
+end
+_MOI.get(cache::_LOAResultCache, ::_MOI.SolverName) = cache.name
+# A real solver throws `GetAttributeNotAllowed` when no bound exists.
+function _MOI.get(
+    cache::_LOAResultCache,
+    attr::Union{_MOI.ObjectiveBound, _MOI.RelativeGap}
+    )
+    haskey(cache.mock.model_attributes, attr) ||
+        throw(_MOI.GetAttributeNotAllowed(attr,
+            "The LOA master did not produce a bound."))
+    return _MOI.get(cache.mock, attr)
+end
+# No duals from an MI(N)LP solve; matches `dual_status = NO_SOLUTION`.
+# Both methods are needed: the caching layer swallows the first throw
+# and recomputes `DualObjectiveValue` from per-constraint duals.
+_MOI.get(cache::_LOAResultCache, attr::_MOI.DualObjectiveValue) =
+    throw(_MOI.GetAttributeNotAllowed(attr, "LOA does not provide duals."))
+function _MOI.get(
+    cache::_LOAResultCache,
+    attr::_MOI.ConstraintDual,
+    index::_MOI.ConstraintIndex
+    )
+    return throw(_MOI.GetAttributeNotAllowed(attr,
+        "LOA does not provide duals."))
+end
+
+################################################################################
+#                       LOAD SOLUTION BY INJECTION
+################################################################################
+# Complete primal for one pool result. The binaries were fixed during
+# the loop (absent from the linearization point), so the combination
+# supplies them.
+function _result_primal(problem::_LOAProblem, result)
+    nlp = problem.nlp
+    V = JuMP.variable_ref_type(typeof(nlp))
+    T = JuMP.value_type(typeof(nlp))
+    primal = Dict{V, T}()
+    for v in JuMP.all_variables(nlp)
+        if haskey(result.linearization_point, v)
+            primal[v] = result.linearization_point[v]
+        elseif JuMP.is_fixed(v)
+            primal[v] = JuMP.fix_value(v)
+        else
+            primal[v] = zero(T)
+        end
+    end
+    for (binary, active) in result.combination
+        primal[binary] = active ? one(T) : zero(T)
+    end
+    return primal
+end
+
+# Attach a `_LOAResultCache` preloaded with the outcome so standard
+# JuMP result queries answer with no solve (an empty pool injects only
+# the statuses). `is_model_dirty = false` re-opens the result queries
+# JuMP guards behind OptimizeNotCalled.
+function _inject_solution(
+    nlp::JuMP.AbstractModel,
+    problem::_LOAProblem,
+    outcome::NamedTuple
+    )
+    T = JuMP.value_type(typeof(nlp))
+    primals = [_result_primal(problem, result)
+        for result in outcome.results]
+    # The mock is the one MOI optimizer that accepts written results.
+    # Serve stored objectives rather than re-evaluating at the primal.
+    JuMP.set_optimizer(nlp, () -> _LOAResultCache(
+        _MOI.Utilities.MockOptimizer(
+            _MOI.Utilities.UniversalFallback(_MOI.Utilities.Model{T}()),
+            T; eval_objective_value = false),
+        outcome.solver_name); add_bridges = false)
+    # Attach copies the model in and creates the `optimizer_index` map.
+    _MOI.Utilities.attach_optimizer(JuMP.backend(nlp))
+    # Result attributes are read-only through the caching layer, so they
+    # are written on the inner optimizer directly.
+    cache = JuMP.unsafe_backend(nlp)
+    _MOI.set(cache, _MOI.TerminationStatus(), outcome.status)
+    _MOI.set(cache, _MOI.RawStatusString(), outcome.raw_status)
+    _MOI.set(cache, _MOI.SolveTimeSec(), Float64(outcome.solve_time))
+    _MOI.set(cache, _MOI.ResultCount(), length(primals))
+    if outcome.bound !== nothing
+        _MOI.set(cache, _MOI.ObjectiveBound(), outcome.bound)
+        isempty(primals) || _MOI.set(cache, _MOI.RelativeGap(),
+            abs(outcome.bound - outcome.results[1].objective) /
+                abs(outcome.results[1].objective))
+    end
+    indices = Dict(variable => JuMP.optimizer_index(variable)
+        for variable in JuMP.all_variables(nlp))
+    for k in eachindex(primals)
+        _MOI.set(cache, _MOI.PrimalStatus(k), _MOI.FEASIBLE_POINT)
+        _MOI.set(cache, _MOI.ObjectiveValue(k),
+            outcome.results[k].objective)
+        for (variable, value) in primals[k]
+            _MOI.set(cache, _MOI.VariablePrimal(k),
+                indices[variable], value)
+        end
+    end
+    # No public API for this; must stay the last mutation.
+    nlp.is_model_dirty = false
+    return
+end
+
+# Finalize a finite GDP. Dispatchable so the InfiniteOpt extension
+# overrides it.
+function _load_solution(
+    model::JuMP.AbstractModel,
+    problem::_LOAProblem,
+    outcome::NamedTuple,
+    method::LOA
+    )
+    # Drop the solver first: an NLP-only solver fails JuMP's ZeroOne
+    # supports check during the binary restore, even detached.
+    _MOI.Utilities.drop_optimizer(JuMP.backend(problem.nlp))
+    _restore_binaries(problem)
+    _inject_solution(problem.nlp, problem, outcome)
+    gdp_data(model).displaced_optimizer = outcome.displaced_optimizer
+    return
+end
+
+################################################################################
+#                       COMBO EXTRACTION (master -> NLP)
+################################################################################
+# Read each binary's master value, rounded to `Bool` (a MILP solver
+# returns values within its integer tolerance, never exactly 0/1).
+function _extract_combination(problem::_LOAProblem, master::_LOAMaster)
+    return Dict(binary => round(Bool, JuMP.value(master.variable_map[binary]))
+        for binary in problem.binaries)
+end
+
+################################################################################
+#                          OA CUT EMISSION
+################################################################################
+# Whether the disjunct behind `binary_ref` is active in `combination`
+# (a complement reference inverts its underlying binary's value).
+_disjunct_active(combination::AbstractDict,
+    binary_ref::JuMP.AbstractVariableRef) = combination[binary_ref]
+_disjunct_active(combination::AbstractDict, binary_ref::JuMP.GenericAffExpr) =
+    !combination[_underlying_binary(binary_ref)]
+
+# Emit all OA cuts for one NLP result: the objective cut, a slacked row
+# per nonlinear global, and a gated cut per active nonlinear disjunct
+# constraint (Big-M / MBM or Hull gating per the inner method).
+function _add_oa_cuts(
+    problem::_LOAProblem,
+    master::_LOAMaster,
+    result::NamedTuple,
+    method::LOA
+    )
+    isempty(result.linearization_point) && return
+    sense_token = Val(master.objective_sense)
+    penalty_sign = _penalty_sign(sense_token)
+    linearization = _linearize_at(master.objective,
+        result.linearization_point, master.variable_map)
+    _add_objective_cut(sense_token, master, linearization, method)
+    for (func, set) in problem.global_constraints
+        linearization = _linearize_at(func, result.linearization_point,
+            master.variable_map)
+        _add_global_oa_row(master, linearization, set, method, penalty_sign)
+    end
+    for (binary_ref, func, set) in problem.disjunct_constraints
+        _disjunct_active(result.combination, binary_ref) || continue
+        linearization = _linearize_at(func, result.linearization_point,
+            master.variable_map)
+        _emit_disjunct_oa_cut(method.inner_method, set, master,
+            _remap_indicator_to_binary(binary_ref, master.variable_map),
+            linearization, method, penalty_sign)
+    end
+    return
+end
+
+# Slacked objective cut. MIN: `lin <= alpha_oa + sigma`; MAX symmetric.
+# The slack (like the disjunct/global cuts) keeps a nonconvex objective
+# linearization from making the master infeasible.
+function _add_objective_cut(
+    sense_token::Val,
+    master::_LOAMaster,
+    lin,
+    method::LOA
+    )
+    penalty_sign = _penalty_sign(sense_token)
+    slack = _penalized_slack(master, method, penalty_sign)
+    _add_objective_cut_body(sense_token, master, lin, slack)
+end
+_add_objective_cut_body(::Val{_MOI.MIN_SENSE}, master, lin, slack) =
+    JuMP.@constraint(master.model, lin <= master.alpha_oa + slack)
+_add_objective_cut_body(::Val{_MOI.MAX_SENSE}, master, lin, slack) =
+    JuMP.@constraint(master.model, lin >= master.alpha_oa - slack)
+
+# Fresh penalized slack added to the master objective. Shared by the
+# disjunct, global, and objective cuts so a nonconvex (invalid)
+# linearization can't make the master infeasible.
+function _penalized_slack(master::_LOAMaster, method::LOA, penalty_sign::Int)
+    slack = JuMP.@variable(master.model,
+        lower_bound = 0.0, upper_bound = method.max_slack)
+    JuMP.set_objective_function(master.model,
+        JuMP.objective_function(master.model) +
+            penalty_sign * method.oa_penalty * slack)
+    return slack
+end
+
+# Extract RHS from an MOI set.
+_set_rhs(s::Union{_MOI.LessThan, _MOI.GreaterThan, _MOI.EqualTo}) =
+    _MOI.constant(s)
+
+# The `<= 0` directions of an OA cut for `set`: `lin - rhs` for LessThan,
+# `rhs - lin` for GreaterThan, both for EqualTo / Interval. The caller
+# adds the gating, slack, and (for Hull) disaggregation per direction.
+_oa_cut_terms(set::_MOI.LessThan, lin) = (lin - _set_rhs(set),)
+_oa_cut_terms(set::_MOI.GreaterThan, lin) = (_set_rhs(set) - lin,)
+_oa_cut_terms(set::_MOI.EqualTo, lin) =
+    (lin - _MOI.constant(set), _MOI.constant(set) - lin)
+_oa_cut_terms(set::_MOI.Interval, lin) = (lin - set.upper, set.lower - lin)
+
+# Big-M / MBM disjunct cut: each direction slacked and gated by
+# `M(1 - binary)`.
+function _emit_disjunct_oa_cut(
+    ::Union{BigM, MBM},
+    set,
+    master::_LOAMaster,
+    binary_ref,
+    linearization,
+    method::LOA,
+    penalty_sign::Int
+    )
+    slack = _penalized_slack(master, method, penalty_sign)
+    for term in _oa_cut_terms(set, linearization)
+        JuMP.@constraint(master.model,
+            term - slack <= method.M_value * (1 - binary_ref))
+    end
+    return
+end
+
+# Convex-hull disjunct cut: `disaggregate_expression` rewrites each
+# direction into disaggregated space (variables -> per-disjunct copies,
+# constant scaled by the binary), so the cut switches off at `y = 0` with
+# no big-M. It is linear in its argument, so the two-sided sets need no
+# special casing.
+function _emit_disjunct_oa_cut(
+    ::Hull,
+    set,
+    master::_LOAMaster,
+    binary_ref,
+    linearization,
+    method::LOA,
+    penalty_sign::Int
+    )
+    slack = _penalized_slack(master, method, penalty_sign)
+    for term in _oa_cut_terms(set, linearization)
+        body = disaggregate_expression(
+            master.model, term, binary_ref, master.disaggregator)
+        JuMP.@constraint(master.model, body - slack <= 0)
+    end
+    return
+end
+
+# Slacked global OA row(s): each direction carries a penalized slack so a
+# nonconvex linearization can't make the master infeasible.
+function _add_global_oa_row(
+    master::_LOAMaster,
+    lin,
+    set,
+    method::LOA,
+    penalty_sign::Int
+    )
+    slack = _penalized_slack(master, method, penalty_sign)
+    for term in _oa_cut_terms(set, lin)
+        JuMP.@constraint(master.model, term <= slack)
+    end
+    return
+end
diff --git a/src/mbm.jl b/src/mbm.jl
index 23c5096e..cb726fc1 100644
--- a/src/mbm.jl
+++ b/src/mbm.jl
@@ -508,7 +508,8 @@ function replace_variables_in_constraint(
     W = _var_ref_type(T, var_map)
     new_aff = zero(JuMP.GenericAffExpr{C, W})
     for (var, coef) in fun.terms
-        JuMP.add_to_expression!(new_aff, coef, var_map[var])
+        JuMP.add_to_expression!(new_aff, coef,
+            replace_variables_in_constraint(var, var_map))
     end
     new_aff.constant = new_aff.constant + fun.constant
     return new_aff
@@ -522,7 +523,9 @@ function replace_variables_in_constraint(
     new_quad = zero(JuMP.GenericQuadExpr{C, W})
     for (vars, coef) in fun.terms
         JuMP.add_to_expression!(new_quad,
-            coef * var_map[vars.a] * var_map[vars.b])
+            coef *
+            replace_variables_in_constraint(vars.a, var_map) *
+            replace_variables_in_constraint(vars.b, var_map))
     end
     new_aff = replace_variables_in_constraint(fun.aff, var_map)
     JuMP.add_to_expression!(new_quad, new_aff)
diff --git a/src/model.jl b/src/model.jl
index 0f26a4f0..c1d23950 100644
--- a/src/model.jl
+++ b/src/model.jl
@@ -39,17 +39,42 @@ end
 
 # Define what should happen to solve a GDPModel
 # See https://github.com/jump-dev/JuMP.jl/blob/9ea1df38fd320f864ab4c93c78631d0f15939c0b/src/JuMP.jl#L718-L745
+# The kwarg form is the installed hook; it redispatches on the method so
+# solution algorithms (e.g. LOA) can own their entire solve.
 function _optimize_hook(
-    model::JuMP.AbstractModel; 
+    model::JuMP.AbstractModel;
     gdp_method::AbstractSolutionMethod = BigM(),
     kwargs...
-    ) # can add more kwargs if wanted
-    if !_ready_to_optimize(model) || _solution_method(model) != gdp_method
-        reformulate_model(model, gdp_method)
+    )
+    return _optimize_hook(model, gdp_method; kwargs...)
+end
+
+# Reformulation methods: reformulate if stale, then solve in place.
+function _optimize_hook(
+    model::JuMP.AbstractModel,
+    method::AbstractReformulationMethod;
+    kwargs...
+    )
+    # A prior LOA solve left its result cache in the slot; put the
+    # optimizer it displaced back, unless the user has since set another.
+    _has_loa_cache(model) && _restore_displaced_optimizer(model)
+    gdp_data(model).displaced_optimizer = nothing
+    if !_ready_to_optimize(model) || _solution_method(model) != method
+        reformulate_model(model, method)
     end
     return JuMP.optimize!(model; ignore_optimize_hook = true, kwargs...)
 end
 
+# Fallback for solution algorithms without an optimize hook.
+function _optimize_hook(
+    model::JuMP.AbstractModel,
+    method::AbstractSolutionMethod;
+    kwargs...
+    )
+    error("Solution method `$(nameof(typeof(method)))` is not " *
+        "supported; LOA is currently the only solution algorithm.")
+end
+
 ################################################################################
 #                              GDP DATA
 ################################################################################
@@ -84,6 +109,7 @@ _constraint_to_indicator(model::JuMP.AbstractModel) = gdp_data(model).constraint
 _reformulation_variables(model::JuMP.AbstractModel) = gdp_data(model).reformulation_variables
 _reformulation_constraints(model::JuMP.AbstractModel) = gdp_data(model).reformulation_constraints
 _variable_bounds(model::JuMP.AbstractModel) = gdp_data(model).variable_bounds
+_disaggregation_map(model::JuMP.AbstractModel) = gdp_data(model).disaggregation_map
 _solution_method(model::JuMP.AbstractModel) = gdp_data(model).solution_method # Get the current solution method
 _ready_to_optimize(model::JuMP.AbstractModel) = gdp_data(model).ready_to_optimize # Determine if the model is ready to call `optimize!` without a optimize hook
 
diff --git a/src/reformulate.jl b/src/reformulate.jl
index 9a3b65cb..264fe294 100644
--- a/src/reformulate.jl
+++ b/src/reformulate.jl
@@ -2,12 +2,20 @@
 #                              REFORMULATE
 ################################################################################
 """
-    reformulate_model(model::JuMP.AbstractModel, method::AbstractSolutionMethod = BigM())
+    reformulate_model(
+        model::JuMP.AbstractModel,
+        method::AbstractReformulationMethod = BigM()
+        )
 
 Reformulate a `GDPModel` using the specified `method`. Prior to reformulation,
-all previous reformulation variables and constraints are deleted.
+all previous reformulation variables and constraints are deleted. Solution
+algorithms such as [`LOA`](@ref) are not reformulation methods and are run
+via `optimize!(model, gdp_method = ...)` instead.
 """
-function reformulate_model(model::JuMP.AbstractModel, method::AbstractSolutionMethod = BigM())
+function reformulate_model(
+    model::JuMP.AbstractModel,
+    method::AbstractReformulationMethod = BigM()
+    )
     #clear all previous reformulations
     _clear_reformulations(model)
     #reformulate
@@ -25,6 +33,7 @@ function _clear_reformulations(model::JuMP.AbstractModel)
     empty!(gdp_data(model).reformulation_constraints)
     empty!(gdp_data(model).reformulation_variables)
     empty!(gdp_data(model).variable_bounds)
+    empty!(gdp_data(model).disaggregation_map)
     return
 end
 
diff --git a/src/utilities.jl b/src/utilities.jl
index bca1b748..2d66e377 100644
--- a/src/utilities.jl
+++ b/src/utilities.jl
@@ -8,6 +8,23 @@ function _copy_model(
     return M()
 end
 
+# OVERRIDABLE. Copy `source`'s decision variables into `target` (with
+# bounds and integrality) and return a `source → target` ref map.
+# Constraints and objective are NOT copied — the caller adds whichever
+# subset it wants. InfiniteOpt overrides to additionally copy
+# parameters, derivatives, and parameter functions.
+function copy_variables_onto_model(
+    target::JuMP.AbstractModel,
+    source::JuMP.AbstractModel
+    )
+    V = JuMP.variable_ref_type(typeof(source))
+    ref_map = Dict{V, V}()
+    for variable in JuMP.all_variables(source)
+        ref_map[variable] = variable_copy(target, variable)
+    end
+    return ref_map
+end
+
 """
     copy_and_reformulate(model, decision_vars, reform_method, method)
 
@@ -62,6 +79,121 @@ function reformulate_and_relax(
     return sub, undo_relax
 end
 
+"""
+    extract_solution(sub::GDPSubmodel)
+
+Read the primal solution of `sub.model` after a solve, keyed by the
+parent-model decision variables via `sub.fwd_map`. Shape follows
+`fwd_map` values: `Vector`-valued fwd_maps (CP/MBM) yield per-support
+`Vector`s; scalar fwd_maps (LOA feas) yield scalars.
+"""
+function extract_solution(sub::GDPSubmodel)
+    return Dict(
+        var => JuMP.value.(sub.fwd_map[var]) for var in sub.decision_vars)
+end
+
+################################################################################
+#                          INDICATOR FIXING
+################################################################################
+"""
+    fix_indicator(model, indicator::LogicalVariableRef, value::Bool)
+    fix_indicator(binary_ref, value::Bool)
+
+Fix a logical indicator's binary backing variable to `value` (true →
+1.0, false → 0.0). The 3-arg form is the user-facing API: pass the
+model and the `LogicalVariableRef`. The 2-arg form takes the binary
+backing reference directly (or its complement-form expression
+`1 - other_binary` when the indicator was declared as a complement).
+
+Mirrors [`relax_logical_vars`](@ref) for selectively fixing rather
+than relaxing.
+"""
+fix_indicator(model::JuMP.AbstractModel,
+    indicator::LogicalVariableRef, value::Bool) =
+    fix_indicator(_indicator_to_binary(model)[indicator], value)
+fix_indicator(binary_ref::JuMP.AbstractVariableRef, value::Bool) =
+    JuMP.fix(binary_ref, value ? 1.0 : 0.0; force = true)
+function fix_indicator(
+    binary_ref::JuMP.GenericAffExpr, value::Bool
+    )
+    underlying, _ = only(binary_ref.terms)
+    JuMP.fix(underlying, value ? 0.0 : 1.0; force = true)
+    return
+end
+
+"""
+    unfix_indicator(model, indicator::LogicalVariableRef)
+    unfix_indicator(binary_ref)
+
+Undo [`fix_indicator`](@ref). No-op if not currently fixed.
+"""
+unfix_indicator(model::JuMP.AbstractModel,
+    indicator::LogicalVariableRef) =
+    unfix_indicator(_indicator_to_binary(model)[indicator])
+function unfix_indicator(binary_ref)
+    JuMP.is_fixed(binary_ref) && JuMP.unfix(binary_ref)
+    return
+end
+function unfix_indicator(binary_ref::JuMP.GenericAffExpr)
+    underlying = only(keys(binary_ref.terms))
+    JuMP.is_fixed(underlying) && JuMP.unfix(underlying)
+    return
+end
+
+################################################################################
+#                          NO-GOOD CUTS
+################################################################################
+"""
+    avoid_combination(model, combination)
+    avoid_combination(model, combination, binary_map)
+
+Add a no-good cut to `model` excluding `combination` from any future
+solution. The added constraint
+    Σ_{j active} (1 - y_j) + Σ_{j inactive} y_j ≥ 1
+forces at least one indicator to differ from `combination`.
+
+`combination` maps `LogicalVariableRef` (or a binary variable) →
+`Bool` (whether it was active). The 3-arg form lets you supply an
+explicit `binary_map` (defaulting to `_indicator_to_binary(model)`)
+translating each combination key to the binary used in the cut, e.g.
+when the binaries belong to a copy of the original model (an LOA
+master).
+
+Returns the constraint reference of the added cut.
+"""
+avoid_combination(model::JuMP.AbstractModel, combination) =
+    avoid_combination(
+        model, combination, _indicator_to_binary(model))
+function avoid_combination(model::JuMP.AbstractModel,
+    combination, binary_map
+    )
+    V = JuMP.variable_ref_type(typeof(model))
+    T = JuMP.value_type(typeof(model))
+    cut = JuMP.GenericAffExpr{T, V}(zero(T))
+    for (indicator, active) in combination
+        haskey(binary_map, indicator) || continue
+        add_no_good_terms(cut, binary_map[indicator], active)
+    end
+    return JuMP.@constraint(model, cut >= one(T))
+end
+
+"""
+    add_no_good_terms(cut, binary_ref, active::Bool)
+
+Append one indicator's contribution to a no-good cut expression:
+adds `1 - y_j` if the indicator was active in the excluded
+combination, or `y_j` otherwise. Used by [`avoid_combination`](@ref).
+"""
+function add_no_good_terms(cut, binary_ref, active::Bool)
+    if active
+        JuMP.add_to_expression!(cut, -1.0, binary_ref)
+        JuMP.add_to_expression!(cut, 1.0)
+    else
+        JuMP.add_to_expression!(cut, 1.0, binary_ref)
+    end
+    return
+end
+
 ################################################################################
 #                          LOGICAL VARIABLE RELAXATION
 ################################################################################
@@ -193,6 +325,7 @@ function copy_gdp_data(
     # - indicator_to_constraints
     # - constraint_to_indicator
     # - variable_bounds
+    # - disaggregation_map
     # SINGLE VALUES (copy directly)
     # - solution_method
     # - ready_to_optimize
@@ -301,6 +434,11 @@ function copy_gdp_data(
         new_gdp.variable_bounds[var_map[v]] = bounds
     end
 
+    # Copying the disaggregation map
+    for ((v, lv_ref), dv) in old_gdp.disaggregation_map
+        new_gdp.disaggregation_map[(var_map[v], lv_map[lv_ref])] = var_map[dv]
+    end
+
     # Copying solution method and ready to optimize
     new_gdp.solution_method = old_gdp.solution_method
     new_gdp.ready_to_optimize = old_gdp.ready_to_optimize
@@ -399,5 +537,101 @@ function _remap_constraint_to_indicator(
     ::Dict{DisjunctConstraintRef{M}, DisjunctConstraintRef{M}},
     disj_map::Dict{DisjunctionRef{M}, DisjunctionRef{M}}
 ) where {M <: JuMP.AbstractModel}
-    return disj_map[con_ref]   
+    return disj_map[con_ref]
+end
+
+################################################################################
+#                    LINEARIZATION & EXPRESSION CONVERSION
+################################################################################
+# First-order Taylor for a single var or affine expression mapped into
+# master space.
+function _linearize_at(
+    variable::JuMP.AbstractVariableRef,
+    ::AbstractDict,
+    ref_map::AbstractDict
+    )
+    target = ref_map[variable]
+    return JuMP.GenericAffExpr{Float64, typeof(target)}(0.0, target => 1.0)
+end
+function _linearize_at(
+    func::JuMP.GenericAffExpr,
+    ::AbstractDict,
+    ref_map::AbstractDict
+    )
+    V = valtype(ref_map)
+    T = JuMP.value_type(V) <: Number ? JuMP.value_type(V) : Float64
+    result = JuMP.GenericAffExpr{T, V}(T(func.constant))
+    for (variable, coefficient) in func.terms
+        JuMP.add_to_expression!(result, coefficient, ref_map[variable])
+    end
+    return result
+end
+
+# Convert JuMP expression trees to Julia Expr with
+# MOI.VariableIndex leaves for MOI.Nonlinear evaluation.
+function _to_nlp_expr(expr::JuMP.GenericNonlinearExpr, idx::Dict)
+    args = Any[_to_nlp_expr(a, idx) for a in expr.args]
+    return Expr(:call, expr.head, args...)
+end
+function _to_nlp_expr(expr::JuMP.GenericAffExpr, idx::Dict)
+    parts = Any[expr.constant]
+    for (var, coef) in expr.terms
+        push!(parts, Expr(:call, :*, coef, _MOI.VariableIndex(idx[var])))
+    end
+    length(parts) == 1 && return parts[1]
+    return Expr(:call, :+, parts...)
+end
+function _to_nlp_expr(expr::JuMP.GenericQuadExpr, idx::Dict)
+    parts = Any[_to_nlp_expr(expr.aff, idx)]
+    for (pair, coef) in expr.terms
+        push!(parts, Expr(:call, :*, coef,
+            _MOI.VariableIndex(idx[pair.a]),
+            _MOI.VariableIndex(idx[pair.b])))
+    end
+    length(parts) == 1 && return parts[1]
+    return Expr(:call, :+, parts...)
+end
+function _to_nlp_expr(var::JuMP.AbstractVariableRef, idx::Dict)
+    return _MOI.VariableIndex(idx[var])
+end
+_to_nlp_expr(x::Number, ::Dict) = x
+
+# First-order Taylor linearization of a quadratic or nonlinear
+# expression at point xk via MOI.Nonlinear reverse-mode AD.
+function _linearize_at(
+    func::Union{JuMP.GenericQuadExpr, JuMP.GenericNonlinearExpr},
+    xk::Dict,
+    ref_map
+    )
+    vars = JuMP.AbstractVariableRef[]
+    _interrogate_variables(v -> push!(vars, v), func)
+    unique!(vars)
+    isempty(vars) && return JuMP.AffExpr(JuMP.value(v -> 0.0, func))
+
+    n = length(vars)
+    T = JuMP.value_type(typeof(JuMP.owner_model(vars[1])))
+    idx = Dict(vars[i] => i for i in 1:n)
+    nlp = _MOI.Nonlinear.Model()
+    _MOI.Nonlinear.set_objective(nlp, _to_nlp_expr(func, idx))
+    ord = [_MOI.VariableIndex(i) for i in 1:n]
+    evaluator = _MOI.Nonlinear.Evaluator(
+        nlp, _MOI.Nonlinear.SparseReverseMode(), ord)
+    _MOI.initialize(evaluator, [:Grad])
+
+    xk_vec = [get(xk, v, zero(T)) for v in vars]
+    f_xk = _MOI.eval_objective(evaluator, xk_vec)
+    grad = zeros(T, n)
+    _MOI.eval_objective_gradient(evaluator, grad, xk_vec)
+
+    constant = T(f_xk)
+    for i in 1:n
+        constant -= grad[i] * xk_vec[i]
+    end
+    V = typeof(ref_map[vars[1]])
+    result = JuMP.GenericAffExpr{T, V}(constant)
+    for i in 1:n
+        iszero(grad[i]) && continue
+        JuMP.add_to_expression!(result, grad[i], ref_map[vars[i]])
+    end
+    return result
 end
diff --git a/test/constraints/cuttingplanes.jl b/test/constraints/cuttingplanes.jl
index 9f3a7931..04ecba59 100644
--- a/test/constraints/cuttingplanes.jl
+++ b/test/constraints/cuttingplanes.jl
@@ -5,15 +5,13 @@ function test_CuttingPlanes_datatype()
     @test method.optimizer == HiGHS.Optimizer
     @test method.max_iter == 3
     @test method.seperation_tolerance == 1e-6
-    @test method.final_reform_method isa BigM
     @test method.M_value == 1e9
 
     method = CuttingPlanes(HiGHS.Optimizer;max_iter=10,
-    seperation_tolerance=1e-4, final_reform_method=Indicator(), M_value=1e6
+    seperation_tolerance=1e-4, M_value=1e6
     )
     @test method.max_iter == 10
     @test method.seperation_tolerance == 1e-4
-    @test method.final_reform_method isa Indicator
     @test method.M_value == 1e6
 end
 
diff --git a/test/constraints/loa.jl b/test/constraints/loa.jl
new file mode 100644
index 00000000..4ab77659
--- /dev/null
+++ b/test/constraints/loa.jl
@@ -0,0 +1,1020 @@
+using HiGHS, Ipopt, Juniper
+
+function test_loa_datatype()
+    method = LOA(HiGHS.Optimizer)
+    @test method.nlp_optimizer == HiGHS.Optimizer
+    @test method.mip_optimizer == HiGHS.Optimizer
+    @test method.max_iter == 10
+    @test method.set_cover_max_iter == 8
+    @test method.M_value == 1e9
+    @test method.max_slack == 1000.0
+    @test method.oa_penalty == 1000.0
+    @test method.use_nlpf == true
+    @test method.convergence_tol == 1e-6
+    @test method.slack_tol == 1e-4
+    @test method.inner_method isa BigM
+    @test method.inner_method.value == 1e9
+
+    method = LOA(HiGHS.Optimizer; max_iter = 50, M_value = 1e6,
+        max_slack = 500.0, oa_penalty = 200.0, use_nlpf = false,
+        convergence_tol = 1e-4, slack_tol = 1e-3)
+    @test method.max_iter == 50
+    @test method.M_value == 1e6
+    @test method.max_slack == 500.0
+    @test method.oa_penalty == 200.0
+    @test method.use_nlpf == false
+    @test method.convergence_tol == 1e-4
+    @test method.slack_tol == 1e-3
+    @test method.inner_method isa BigM
+    @test method.inner_method.value == 1e6
+
+    method = LOA(HiGHS.Optimizer; inner_method = MBM(HiGHS.Optimizer))
+    @test method.inner_method isa MBM
+end
+
+function test_cover_disjuncts()
+    # Only disjuncts owning a nonlinear constraint need covering: the two
+    # `x^2` disjuncts do, the two linear `x` disjuncts do not.
+    model = GDPModel()
+    @variable(model, -5 <= x <= 5)
+    @variable(model, Y[1:2], Logical)
+    @constraint(model, x^2 <= 3, Disjunct(Y[1]))
+    @constraint(model, x^2 >= 5, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @variable(model, W[1:2], Logical)
+    @constraint(model, x <= 1, Disjunct(W[1]))
+    @constraint(model, x >= -1, Disjunct(W[2]))
+    @disjunction(model, W)
+
+    DP.reformulate_model(model, BigM(1e9))
+    method = LOA(HiGHS.Optimizer)
+    problem = DP.build_loa_problem(model, method)
+    disjuncts = DP._cover_disjuncts(problem)
+
+    # Two nonlinear disjuncts to cover; the linear W disjuncts omitted.
+    @test length(disjuncts) == 2
+    binary_map = DP._indicator_to_binary(model)
+    underlying = Set(DP._underlying_binary(dref) for dref in disjuncts)
+    @test underlying == Set([DP._underlying_binary(binary_map[Y[1]]),
+        DP._underlying_binary(binary_map[Y[2]])])
+end
+
+function test_oa_cut_terms()
+    # The `<= 0` cut directions per set, at a scalar linearization value
+    # of 5: LessThan/GreaterThan give one signed direction, EqualTo and
+    # Interval give both.
+    @test DP._oa_cut_terms(MOI.LessThan(2.0), 5.0) == (3.0,)
+    @test DP._oa_cut_terms(MOI.GreaterThan(2.0), 5.0) == (-3.0,)
+    @test DP._oa_cut_terms(MOI.EqualTo(2.0), 5.0) == (3.0, -3.0)
+    @test DP._oa_cut_terms(MOI.Interval(1.0, 4.0), 5.0) == (1.0, -4.0)
+end
+
+function test_is_linear_F()
+    # Scalar and vector variable-ref / affine function types are linear
+    # (copied into the master at build time); anything else, e.g. a
+    # quadratic, is nonlinear and enters only as an OA cut.
+    @test DP._is_linear_F(JuMP.VariableRef)
+    @test DP._is_linear_F(JuMP.AffExpr)
+    @test DP._is_linear_F(Vector{JuMP.VariableRef})
+    @test DP._is_linear_F(Vector{JuMP.AffExpr})
+    @test !DP._is_linear_F(JuMP.QuadExpr)
+end
+
+function test_no_good_cut()
+    model = GDPModel()
+    @variable(model, x)
+    @variable(model, Y[1:2], Logical)
+    @constraint(model, x <= 3, Disjunct(Y[1]))
+    @constraint(model, x >= 5, Disjunct(Y[2]))
+    @disjunction(model, Y)
+
+    DP.reformulate_model(model, BigM(1e9))
+    method = LOA(HiGHS.Optimizer)
+    problem = DP.build_loa_problem(model, method)
+    master = DP._build_loa_master(problem, method)
+    master_model = master.model
+
+    binary_map = DP._indicator_to_binary(model)
+    combo = Dict(
+        DP._underlying_binary(binary_map[Y[1]]) =>
+            DP._underlying_value(binary_map[Y[1]], true),
+        DP._underlying_binary(binary_map[Y[2]]) =>
+            DP._underlying_value(binary_map[Y[2]], false))
+
+    num_cons_before = length(JuMP.all_constraints(
+        master_model;
+        include_variable_in_set_constraints = false))
+    cref = DP.avoid_combination(master.model, combo, master.variable_map)
+    num_cons_after = length(JuMP.all_constraints(
+        master_model;
+        include_variable_in_set_constraints = false))
+
+    @test num_cons_after == num_cons_before + 1
+    # The cut is (1 - y1) + y2 >= 1, i.e. normalized -y1 + y2 >= 0:
+    # excludes exactly the (Y1 active, Y2 inactive) combination.
+    y1 = master.variable_map[binary_map[Y[1]]]
+    y2 = master.variable_map[binary_map[Y[2]]]
+    @test JuMP.normalized_coefficient(cref, y1) == -1.0
+    @test JuMP.normalized_coefficient(cref, y2) == 1.0
+    @test JuMP.normalized_rhs(cref) == 0.0
+end
+
+function test_loa_hook_simple()
+    model = GDPModel(HiGHS.Optimizer)
+    set_silent(model)
+    @variable(model, 0 <= x <= 10)
+    @variable(model, Y[1:2], Logical)
+    @constraint(model, x <= 3, Disjunct(Y[1]))
+    @constraint(model, x <= 7, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @objective(model, Max, x)
+
+    method = LOA(HiGHS.Optimizer)
+    optimize!(model, gdp_method = method)
+
+    # LOA marks the model as needing a rebuild before a reformulation solve.
+    @test !DP._ready_to_optimize(model)
+    @test DP._solution_method(model) === method
+    # The incumbent is injected, not solved: queries return x = 7 via Y[2].
+    @test objective_value(model) ≈ 7.0 atol = 1e-6
+    @test value(x) ≈ 7.0 atol = 1e-6
+
+    # Injection leaves a clean model: the incumbent is loaded with no fixed
+    # binaries and no pin constraints.
+    @test !is_fixed(DP._indicator_to_binary(model)[Y[2]])
+end
+
+function test_loa_reoptimize()
+    # After an LOA solve the model must re-solve cleanly with an ordinary
+    # reformulation method (the injected mock optimizer is swapped back).
+    model = GDPModel(HiGHS.Optimizer)
+    set_silent(model)
+    @variable(model, 0 <= x <= 10)
+    @variable(model, Y[1:2], Logical)
+    @constraint(model, x <= 3, Disjunct(Y[1]))
+    @constraint(model, x <= 7, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @objective(model, Max, x)
+
+    optimize!(model, gdp_method = LOA(HiGHS.Optimizer))
+    @test value(x) ≈ 7.0 atol = 1e-6
+
+    optimize!(model, gdp_method = Hull())
+    @test termination_status(model) == MOI.OPTIMAL
+    @test objective_value(model) ≈ 7.0 atol = 1e-4
+end
+
+function test_loa_reoptimize_respects_user_optimizer()
+    # If the user sets their own optimizer after LOA, a re-solve must keep
+    # it (the stashed solver only restores when the injected mock is still
+    # attached). Regression: an unconditional restore clobbered the choice.
+    model = GDPModel(HiGHS.Optimizer)
+    set_silent(model)
+    @variable(model, 0 <= x <= 10)
+    @variable(model, Y[1:2], Logical)
+    @constraint(model, x <= 3, Disjunct(Y[1]))
+    @constraint(model, x <= 7, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @objective(model, Max, x)
+
+    optimize!(model, gdp_method = LOA(HiGHS.Optimizer))
+    set_optimizer(model, HiGHS.Optimizer)
+    set_attribute(model, "time_limit", 123.0)
+    set_silent(model)
+    optimize!(model, gdp_method = BigM())
+    # A clobbering restore would rebuild a fresh solver, losing the setting.
+    @test get_attribute(model, "time_limit") == 123.0
+    @test value(x) ≈ 7.0 atol = 1e-6
+end
+
+function test_loa_solve_simple()
+    # Simple GDP: max x s.t. (x <= 3) OR (x <= 7), 0 <= x <= 10
+    # Optimal: x = 7 (select second disjunct)
+    model = GDPModel(HiGHS.Optimizer)
+    set_silent(model)
+    @variable(model, 0 <= x <= 10)
+    @variable(model, Y[1:2], Logical)
+    @constraint(model, x <= 3, Disjunct(Y[1]))
+    @constraint(model, x <= 7, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @objective(model, Max, x)
+
+    optimize!(model, gdp_method = LOA(HiGHS.Optimizer))
+    @test termination_status(model) == MOI.LOCALLY_SOLVED
+    @test objective_value(model) ≈ 7.0 atol=1e-4
+end
+
+function test_loa_solve_simple_with_mbm()
+    # Same GDP as test_loa_solve_simple, but inner_method = MBM so the
+    # NLP model is built with per-constraint tight Ms instead of BigM.
+    # Optimum unchanged: x = 7 via the second disjunct.
+    model = GDPModel(HiGHS.Optimizer)
+    set_silent(model)
+    @variable(model, 0 <= x <= 10)
+    @variable(model, Y[1:2], Logical)
+    @constraint(model, x <= 3, Disjunct(Y[1]))
+    @constraint(model, x <= 7, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @objective(model, Max, x)
+
+    optimize!(model, gdp_method = LOA(HiGHS.Optimizer;
+        inner_method = MBM(HiGHS.Optimizer)))
+    @test termination_status(model) == MOI.LOCALLY_SOLVED
+    @test objective_value(model) ≈ 7.0 atol=1e-4
+end
+
+function test_loa_solve_two_disjunctions()
+    # Two disjunctions: max x + z
+    # D1: (x <= 3) OR (x <= 7)
+    # D2: (z <= 2) OR (z <= 5)
+    # Optimal: x = 7, z = 5
+    model = GDPModel(HiGHS.Optimizer)
+    set_silent(model)
+    @variable(model, 0 <= x <= 10)
+    @variable(model, 0 <= z <= 10)
+    @variable(model, Y[1:2], Logical)
+    @variable(model, W[1:2], Logical)
+    @constraint(model, x <= 3, Disjunct(Y[1]))
+    @constraint(model, x <= 7, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @constraint(model, z <= 2, Disjunct(W[1]))
+    @constraint(model, z <= 5, Disjunct(W[2]))
+    @disjunction(model, W)
+    @objective(model, Max, x + z)
+
+    optimize!(model, gdp_method = LOA(HiGHS.Optimizer))
+    @test termination_status(model) == MOI.LOCALLY_SOLVED
+    @test objective_value(model) ≈ 12.0 atol=1e-4
+end
+
+function test_loa_not_reformulation()
+    # LOA is a solution algorithm: reformulate_model rejects it.
+    method = LOA(HiGHS.Optimizer)
+    model = GDPModel(HiGHS.Optimizer)
+    @test_throws MethodError DP.reformulate_model(model, method)
+end
+
+function test_loa_nonlinear_global()
+    # max x s.t. x^2 <= 25 (global), (x <= 3) ∨ (x <= 8), 0 <= x <= 10.
+    # Disjunct Y[2] permits x up to 8 but the global x^2 <= 25 bounds
+    # x to 5. Verifies that the global-cut pass emits the
+    # linearization of the global into the master without breaking
+    # the loop. Result must hit the global-binding optimum.
+    ipopt = optimizer_with_attributes(Ipopt.Optimizer,
+        "print_level" => 0, "sb" => "yes")
+    juniper = optimizer_with_attributes(Juniper.Optimizer,
+        "nl_solver" => ipopt, "log_levels" => [])
+    model = GDPModel(juniper)
+    set_silent(model)
+    @variable(model, 0 <= x <= 10)
+    @constraint(model, x^2 <= 25)
+    @variable(model, Y[1:2], Logical)
+    @constraint(model, x <= 3, Disjunct(Y[1]))
+    @constraint(model, x <= 8, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @objective(model, Max, x)
+    optimize!(model,
+        gdp_method = LOA(juniper; mip_optimizer = HiGHS.Optimizer))
+    @test termination_status(model) in (MOI.OPTIMAL, MOI.LOCALLY_SOLVED)
+    @test objective_value(model) ≈ 5.0 atol = 1e-3
+end
+
+function test_loa_nonlinear_equality_global()
+    # max x s.t. x^2 == 25 (global nonlinear equality), (x <= 3) ∨
+    # (x <= 8), 0 <= x <= 10. The Y[1] seed is NLP-infeasible (x <= 3
+    # contradicts x = 5), so NLPF slacks the inequality while keeping
+    # the equality exact. The equality emits BOTH cut directions into
+    # the master. Optimum: x = 5 via Y[2].
+    ipopt = optimizer_with_attributes(Ipopt.Optimizer,
+        "print_level" => 0, "sb" => "yes")
+    model = GDPModel(ipopt)
+    set_silent(model)
+    @variable(model, 0 <= x <= 10)
+    @constraint(model, x^2 == 25)
+    @variable(model, Y[1:2], Logical)
+    @constraint(model, x <= 3, Disjunct(Y[1]))
+    @constraint(model, x <= 8, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @objective(model, Max, x)
+    optimize!(model,
+        gdp_method = LOA(ipopt; mip_optimizer = HiGHS.Optimizer))
+    @test termination_status(model) in (MOI.OPTIMAL, MOI.LOCALLY_SOLVED)
+    @test objective_value(model) ≈ 5.0 atol = 1e-3
+    @test value(x) ≈ 5.0 atol = 1e-3
+    @test value(Y[2]) ≈ 1.0 atol = 1e-6
+end
+
+function test_loa_nonlinear_equality_disjunct()
+    # Nonlinear equality inside a disjunct: Y1: x <= 3, Y2: x^2 == 64.
+    # The disjunct cut emits both gated directions. Optimum: x = 8.
+    ipopt = optimizer_with_attributes(Ipopt.Optimizer,
+        "print_level" => 0, "sb" => "yes")
+    juniper = optimizer_with_attributes(Juniper.Optimizer,
+        "nl_solver" => ipopt, "log_levels" => [])
+    model = GDPModel(juniper)
+    set_silent(model)
+    @variable(model, 0 <= x <= 10)
+    @variable(model, Y[1:2], Logical)
+    @constraint(model, x <= 3, Disjunct(Y[1]))
+    @constraint(model, x^2 == 64, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @objective(model, Max, x)
+    optimize!(model,
+        gdp_method = LOA(juniper; mip_optimizer = HiGHS.Optimizer))
+    @test termination_status(model) in (MOI.OPTIMAL, MOI.LOCALLY_SOLVED)
+    @test objective_value(model) ≈ 8.0 atol = 1e-3
+    @test value(Y[2]) ≈ 1.0 atol = 1e-6
+end
+
+function test_loa_nonlinear_interval_disjunct()
+    # Nonlinear Interval constraint inside a disjunct: Y1: x <= 3,
+    # Y2: 36 <= x^2 <= 64. The disjunct cut emits both gated
+    # directions. Optimum: x = 8 via Y2.
+    ipopt = optimizer_with_attributes(Ipopt.Optimizer,
+        "print_level" => 0, "sb" => "yes")
+    juniper = optimizer_with_attributes(Juniper.Optimizer,
+        "nl_solver" => ipopt, "log_levels" => [])
+    model = GDPModel(juniper)
+    set_silent(model)
+    @variable(model, 0 <= x <= 10)
+    @variable(model, Y[1:2], Logical)
+    @constraint(model, x <= 3, Disjunct(Y[1]))
+    @constraint(model, 36 <= x^2 <= 64, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @objective(model, Max, x)
+    optimize!(model,
+        gdp_method = LOA(juniper; mip_optimizer = HiGHS.Optimizer))
+    @test termination_status(model) in (MOI.OPTIMAL, MOI.LOCALLY_SOLVED)
+    @test objective_value(model) ≈ 8.0 atol = 1e-3
+    @test value(Y[2]) ≈ 1.0 atol = 1e-6
+end
+
+function test_loa_restores_prior_time_limit()
+    # A time limit set by the user before LOA must survive the loop's
+    # per-solve caps when only `iteration_time_limit` is finite (the
+    # `_restore_time_limit(::Real)` path).
+    model = GDPModel(HiGHS.Optimizer)
+    set_silent(model)
+    @variable(model, 0 <= x <= 10)
+    @variable(model, Y[1:2], Logical)
+    @constraint(model, x <= 3, Disjunct(Y[1]))
+    @constraint(model, x <= 7, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @objective(model, Max, x)
+    set_time_limit_sec(model, 90.0)
+    optimize!(model, gdp_method = LOA(HiGHS.Optimizer;
+        iteration_time_limit = 60.0, time_limit = Inf))
+    @test time_limit_sec(model) == 90.0
+    @test termination_status(model) == MOI.LOCALLY_SOLVED
+    @test objective_value(model) ≈ 7.0 atol = 1e-4
+end
+
+function test_loa_restores_displaced_optimizer()
+    # LOA solves the NLP on the model itself, so it displaces the
+    # model's own optimizer. A later reformulation solve must get that
+    # optimizer back, not LOA's NLP solver.
+    model = GDPModel(HiGHS.Optimizer)
+    set_silent(model)
+    set_optimizer_attribute(model, "mip_rel_gap", 0.25)
+    @variable(model, 0 <= x <= 10)
+    @variable(model, Y[1:2], Logical)
+    @constraint(model, x <= 3, Disjunct(Y[1]))
+    @constraint(model, x <= 7, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @objective(model, Max, x)
+    displaced = JuMP.unsafe_backend(model)
+
+    optimize!(model, gdp_method = LOA(HiGHS.Optimizer))
+    @test DP.gdp_data(model).displaced_optimizer === displaced
+
+    # The reformulation solve restores it, attributes intact.
+    optimize!(model, gdp_method = BigM())
+    @test JuMP.unsafe_backend(model) === displaced
+    @test get_optimizer_attribute(model, "mip_rel_gap") == 0.25
+    @test DP.gdp_data(model).displaced_optimizer === nothing
+    @test termination_status(model) == MOI.OPTIMAL
+    @test objective_value(model) ≈ 7.0 atol = 1e-4
+end
+
+function test_loa_repeated_solve_keeps_displaced_optimizer()
+    # A second LOA solve displaces the first one's result cache, not the
+    # model's own optimizer. The stash must still hold the latter.
+    model = GDPModel(HiGHS.Optimizer)
+    set_silent(model)
+    @variable(model, 0 <= x <= 10)
+    @variable(model, Y[1:2], Logical)
+    @constraint(model, x <= 3, Disjunct(Y[1]))
+    @constraint(model, x <= 7, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @objective(model, Max, x)
+    displaced = JuMP.unsafe_backend(model)
+
+    optimize!(model, gdp_method = LOA(HiGHS.Optimizer))
+    optimize!(model, gdp_method = LOA(HiGHS.Optimizer))
+    @test DP.gdp_data(model).displaced_optimizer === displaced
+
+    optimize!(model, gdp_method = BigM())
+    @test JuMP.unsafe_backend(model) === displaced
+    @test termination_status(model) == MOI.OPTIMAL
+    @test objective_value(model) ≈ 7.0 atol = 1e-4
+end
+
+function test_loa_no_optimizer_to_restore()
+    # With no optimizer to put back, the result cache must be dropped
+    # rather than left to serve its stored results as a fresh solve.
+    model = GDPModel()
+    @variable(model, 0 <= x <= 10)
+    @variable(model, Y[1:2], Logical)
+    @constraint(model, x <= 3, Disjunct(Y[1]))
+    @constraint(model, x <= 7, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @objective(model, Max, x)
+
+    optimize!(model, gdp_method = LOA(HiGHS.Optimizer))
+    @test DP.gdp_data(model).displaced_optimizer === nothing
+    @test objective_value(model) ≈ 7.0 atol = 1e-4
+
+    @test_throws JuMP.NoOptimizer optimize!(model, gdp_method = BigM())
+    set_optimizer(model, HiGHS.Optimizer)
+    set_silent(model)
+    optimize!(model, gdp_method = BigM())
+    @test objective_value(model) ≈ 7.0 atol = 1e-4
+end
+
+function test_loa_cache_refuses_direct_solve()
+    # `ignore_optimize_hook = true` would "solve" the result cache and
+    # clear JuMP's dirty flag, reporting a stale point as freshly
+    # solved. It must error instead.
+    model = GDPModel(HiGHS.Optimizer)
+    set_silent(model)
+    @variable(model, 0 <= x <= 10)
+    @variable(model, Y[1:2], Logical)
+    @constraint(model, x <= 3, Disjunct(Y[1]))
+    @constraint(model, x <= 7, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @objective(model, Max, x)
+    optimize!(model, gdp_method = LOA(HiGHS.Optimizer))
+    @test_throws ErrorException optimize!(model;
+        ignore_optimize_hook = true)
+
+    # Still true once the model has been edited underneath the cache.
+    @constraint(model, x <= 1)
+    @test !has_values(model)
+    @test_throws ErrorException optimize!(model;
+        ignore_optimize_hook = true)
+end
+
+function test_loa_limit_hit_report()
+    # Stop the main loop on `max_iter = 1` after the master has produced
+    # a bound: the report must label the run "limit hit" (not converged),
+    # and the single loop iteration still finds the off-diagonal optimum.
+    model = GDPModel(HiGHS.Optimizer)
+    set_silent(model)
+    @variable(model, 0 <= x <= 10)
+    @variable(model, 0 <= z <= 10)
+    @variable(model, Y[1:2], Logical)
+    @variable(model, W[1:2], Logical)
+    @constraint(model, x <= 4, Disjunct(Y[1]))
+    @constraint(model, x >= 6, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @constraint(model, z <= 4, Disjunct(W[1]))
+    @constraint(model, z >= 6, Disjunct(W[2]))
+    @disjunction(model, W)
+    @objective(model, Min, x - z)
+    method = LOA(HiGHS.Optimizer; max_iter = 1)
+    @test_logs (:info, r"limit hit") match_mode = :any begin
+        optimize!(model, gdp_method = method)
+    end
+    @test termination_status(model) == MOI.ITERATION_LIMIT
+    @test primal_status(model) == MOI.FEASIBLE_POINT
+    @test objective_value(model) ≈ -10.0 atol = 1e-4
+end
+
+function test_loa_complement_indicator_nonlinear_disjunct()
+    # Regression: complement-form indicators carry `1 - y_base` (an
+    # AffExpr) as their binary reference. When the complement disjunct has a
+    # nonlinear constraint, the disjunct-cut pass feeds that AffExpr
+    # straight into the OA cut gating term `M(1 - binary)`, which must
+    # accept an AffExpr binary (not just a plain variable ref).
+    #
+    # Setup: 0 <= x <= 10, Y2 ≡ ¬Y1.
+    #   Disjunct(Y1): x <= 3     (linear)
+    #   Disjunct(Y2): x^2 <= 64  (nonlinear, optimum: x = 8 with Y2)
+    ipopt = optimizer_with_attributes(Ipopt.Optimizer,
+        "print_level" => 0, "sb" => "yes")
+    juniper = optimizer_with_attributes(Juniper.Optimizer,
+        "nl_solver" => ipopt, "log_levels" => [])
+    model = GDPModel(juniper)
+    set_silent(model)
+    @variable(model, 0 <= x <= 10)
+    @variable(model, Y1, Logical)
+    @variable(model, Y2, Logical, logical_complement = Y1)
+    @constraint(model, x <= 3, Disjunct(Y1))
+    @constraint(model, x^2 <= 64, Disjunct(Y2))
+    @disjunction(model, [Y1, Y2])
+    @objective(model, Max, x)
+    optimize!(model,
+        gdp_method = LOA(juniper; mip_optimizer = HiGHS.Optimizer))
+    @test termination_status(model) in (MOI.OPTIMAL, MOI.LOCALLY_SOLVED)
+    @test objective_value(model) ≈ 8.0 atol = 1e-3
+end
+
+function test_loa_nlpf_infeasible_disjunct()
+    # Y1 disjunct constraint x^2 >= 200 is NLP-infeasible against the
+    # variable bound x in [0, 10] (max x^2 = 100). The primary NLP at
+    # the Y1=true seed therefore fails and NLPF (V&G 1990 eq. 8) must
+    # kick in: slack `u` on the GreaterThan, minimize u, return the
+    # x = 10 / u = 100 primal as the linearization site for the OA cut
+    # on x^2 >= 200. With the master's per-cut penalized slack
+    # absorbing the resulting (invalid in isolation) linearization and
+    # the no-good cut forbidding Y1, LOA should still converge to the
+    # Y2=true optimum x = 5.
+    ipopt = optimizer_with_attributes(Ipopt.Optimizer,
+        "print_level" => 0, "sb" => "yes")
+    juniper = optimizer_with_attributes(Juniper.Optimizer,
+        "nl_solver" => ipopt, "log_levels" => [])
+    model = GDPModel(juniper)
+    set_silent(model)
+    @variable(model, 0 <= x <= 10)
+    @variable(model, Y[1:2], Logical)
+    @constraint(model, x^2 >= 200, Disjunct(Y[1]))
+    @constraint(model, x <= 5, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @objective(model, Max, x)
+    optimize!(model,
+        gdp_method = LOA(juniper; mip_optimizer = HiGHS.Optimizer))
+    @test termination_status(model) in (MOI.OPTIMAL, MOI.LOCALLY_SOLVED)
+    @test objective_value(model) ≈ 5.0 atol = 1e-3
+    @test value(Y[2]) ≈ 1.0 atol = 1e-6
+    @test value(Y[1]) ≈ 0.0 atol = 1e-6
+end
+
+function test_loa_nlpf_disabled()
+    # `use_nlpf = false` skips the NLPF fall-through entirely: an
+    # NLP-infeasible combination returns the bare infeasible sentinel
+    # (empty linearization point), so it contributes only its no-good
+    # cut. Same setup as test_loa_nlpf_infeasible_disjunct, where NLPF
+    # WOULD succeed (x = 10, u = 100) and return a nonempty
+    # linearization point — the empty point proves NLPF was skipped.
+    ipopt = optimizer_with_attributes(Ipopt.Optimizer,
+        "print_level" => 0, "sb" => "yes")
+    model = GDPModel()
+    @variable(model, 0 <= x <= 10)
+    @variable(model, Y[1:2], Logical)
+    @constraint(model, x^2 >= 200, Disjunct(Y[1]))
+    @constraint(model, x <= 5, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @objective(model, Max, x)
+
+    DP.reformulate_model(model, BigM(1e9))
+    method = LOA(ipopt; use_nlpf = false)
+    problem = DP.build_loa_problem(model, method)
+    DP._relax_binaries(problem)
+    JuMP.set_optimizer(problem.nlp, method.nlp_optimizer)
+    JuMP.set_silent(problem.nlp)
+
+    binary_map = DP._indicator_to_binary(model)
+    combination = Dict(
+        DP._underlying_binary(binary_map[Y[1]]) => true,
+        DP._underlying_binary(binary_map[Y[2]]) => false)
+
+    result = DP._solve_nlp(problem, combination, method)
+    @test result.feasible == false
+    @test result.objective == Inf
+    @test isempty(result.linearization_point)
+
+    # End-to-end: with NLPF off, LOA still reaches the Y2 optimum x = 5
+    # off no-good cuts alone.
+    juniper = optimizer_with_attributes(Juniper.Optimizer,
+        "nl_solver" => ipopt, "log_levels" => [])
+    model = GDPModel(juniper)
+    set_silent(model)
+    @variable(model, 0 <= x <= 10)
+    @variable(model, Y[1:2], Logical)
+    @constraint(model, x^2 >= 200, Disjunct(Y[1]))
+    @constraint(model, x <= 5, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @objective(model, Max, x)
+    optimize!(model, gdp_method = LOA(juniper;
+        mip_optimizer = HiGHS.Optimizer, use_nlpf = false))
+    @test termination_status(model) in (MOI.OPTIMAL, MOI.LOCALLY_SOLVED)
+    @test objective_value(model) ≈ 5.0 atol = 1e-3
+    @test value(Y[2]) ≈ 1.0 atol = 1e-6
+end
+
+function test_loa_solve_nlp_infeasible_fallback()
+    # `_solve_nlp` fallback: both the primary NLP and NLPF are
+    # infeasible at the fixed combination, so it returns the infeasible
+    # sentinel (objective = Inf, feasible = false, empty linearization
+    # point, combination echoed back unchanged).
+    #
+    # The global x^2 == 400 is unsatisfiable under 0 <= x <= 10 (needs
+    # x = 20). It is NONLINEAR, so it never enters the master and the
+    # primary NLP is the one that catches it as infeasible. It is an
+    # EQUALITY, so NLPF does not slack it and NLPF is infeasible too.
+    # This is what distinguishes the path from
+    # test_loa_no_feasible_incumbent, whose LINEAR equality makes the
+    # master infeasible first, short-circuiting the NLP solve entirely.
+    ipopt = optimizer_with_attributes(Ipopt.Optimizer,
+        "print_level" => 0, "sb" => "yes")
+    model = GDPModel()
+    @variable(model, 0 <= x <= 10)
+    @constraint(model, x^2 == 400)
+    @variable(model, Y[1:2], Logical)
+    @constraint(model, x <= 3, Disjunct(Y[1]))
+    @constraint(model, x <= 7, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @objective(model, Max, x)
+
+    DP.reformulate_model(model, BigM(1e9))
+    method = LOA(ipopt)
+    problem = DP.build_loa_problem(model, method)
+    DP._relax_binaries(problem)
+    JuMP.set_optimizer(problem.nlp, method.nlp_optimizer)
+    JuMP.set_silent(problem.nlp)
+
+    binary_map = DP._indicator_to_binary(model)
+    combination = Dict(
+        DP._underlying_binary(binary_map[Y[1]]) => false,
+        DP._underlying_binary(binary_map[Y[2]]) => true)
+
+    result = DP._solve_nlp(problem, combination, method)
+    @test result.feasible == false
+    @test result.objective == Inf
+    @test isempty(result.linearization_point)
+    @test result.combination === combination
+end
+
+function test_loa_sense_primitives()
+    # Sense-dispatched primitives the main loop reads through. The
+    # finite/infinite solve tests cover every branch except the
+    # minimize-sense gap (no minimizing model reaches a master bound in
+    # those tests), so pin all six here for both senses directly.
+    minv = Val(MOI.MIN_SENSE)
+    maxv = Val(MOI.MAX_SENSE)
+    @test DP._penalty_sign(minv) == 1
+    @test DP._penalty_sign(maxv) == -1
+    @test DP._worst_objective(minv) == Inf
+    @test DP._worst_objective(maxv) == -Inf
+    @test DP._is_better(minv, 3.0, 5.0)
+    @test !DP._is_better(minv, 5.0, 3.0)
+    @test DP._is_better(maxv, 5.0, 3.0)
+    @test !DP._is_better(maxv, 3.0, 5.0)
+    @test DP._gap(minv, 5.0, 3.0) == 2.0
+    @test DP._gap(maxv, 3.0, 5.0) == 2.0
+end
+
+function test_loa_iteration_loop()
+    # Force the master/NLP refinement loop to run its body rather than
+    # exhaust the master on the seeds or converge on the first master
+    # solve. Two disjunctions over disjoint x/z half-lines put the
+    # optimum on an OFF-diagonal combination the set-covering seeds (the
+    # diagonal (Y1,W1) and (Y2,W2)) never try. After seeding, the master
+    # stays feasible and its bound beats the seed incumbent, so LOA
+    # enters the loop body, extracts the off-diagonal combination, solves
+    # its NLP, and updates the incumbent before converging. `Min` sense
+    # also drives the minimize-sense gap path.
+    #   min x - z
+    #   D1: (x <= 4) [Y1]  v  (x >= 6) [Y2]
+    #   D2: (z <= 4) [W1]  v  (z >= 6) [W2]
+    # Optimum: Y1 & W2 -> x = 0, z = 10, obj = -10 (off-diagonal).
+    model = GDPModel(HiGHS.Optimizer)
+    set_silent(model)
+    @variable(model, 0 <= x <= 10)
+    @variable(model, 0 <= z <= 10)
+    @variable(model, Y[1:2], Logical)
+    @variable(model, W[1:2], Logical)
+    @constraint(model, x <= 4, Disjunct(Y[1]))
+    @constraint(model, x >= 6, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @constraint(model, z <= 4, Disjunct(W[1]))
+    @constraint(model, z >= 6, Disjunct(W[2]))
+    @disjunction(model, W)
+    @objective(model, Min, x - z)
+    optimize!(model, gdp_method = LOA(HiGHS.Optimizer))
+    @test termination_status(model) == MOI.LOCALLY_SOLVED
+    @test objective_value(model) ≈ -10.0 atol = 1e-4
+    @test value(Y[1]) ≈ 1.0 atol = 1e-6
+    @test value(W[2]) ≈ 1.0 atol = 1e-6
+end
+
+function test_loa_time_limits()
+    # Exercise the wall-clock budget paths. A finite `time_limit`
+    # (overall cap) drives `_cap_remaining_time` during the subproblem
+    # solves; a finite `iteration_time_limit` (loop-only budget) covers
+    # the post-loop time-limit restore from no prior limit (the
+    # `::Nothing` path). Both limits are generous: they govern the
+    # path taken, not the result.
+    function build_model()
+        model = GDPModel(HiGHS.Optimizer)
+        set_silent(model)
+        @variable(model, 0 <= x <= 10)
+        @variable(model, Y[1:2], Logical)
+        @constraint(model, x <= 3, Disjunct(Y[1]))
+        @constraint(model, x <= 7, Disjunct(Y[2]))
+        @disjunction(model, Y)
+        @objective(model, Max, x)
+        return model
+    end
+
+    model = build_model()
+    optimize!(model, gdp_method = LOA(HiGHS.Optimizer; time_limit = 60.0))
+    @test termination_status(model) == MOI.LOCALLY_SOLVED
+    @test objective_value(model) ≈ 7.0 atol = 1e-4
+
+    model = build_model()
+    optimize!(model, gdp_method = LOA(HiGHS.Optimizer;
+        iteration_time_limit = 60.0, time_limit = Inf))
+    @test termination_status(model) == MOI.LOCALLY_SOLVED
+    @test objective_value(model) ≈ 7.0 atol = 1e-4
+end
+
+function test_loa_no_feasible_incumbent()
+    # Every disjunct combination is NLP-infeasible: the global x == 20
+    # contradicts the x <= 10 bound under any selection. NLPF does not
+    # slack equalities, so NLPF also fails (no values) and `_solve_nlp`
+    # returns its infeasible fallback. No seed or loop NLP yields a
+    # feasible incumbent, so `best_result` stays `nothing`, the master
+    # is infeasible (no bound), and LOA exits on the "no feasible
+    # incumbent" warning from `_loa_status_message`.
+    model = GDPModel(HiGHS.Optimizer)
+    set_silent(model)
+    @variable(model, 0 <= x <= 10)
+    @constraint(model, x == 20)
+    @variable(model, Y[1:2], Logical)
+    @constraint(model, x <= 3, Disjunct(Y[1]))
+    @constraint(model, x <= 7, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @objective(model, Max, x)
+    method = LOA(HiGHS.Optimizer)
+    @test_logs (:warn, r"no feasible incumbent") match_mode = :any begin
+        optimize!(model, gdp_method = method)
+    end
+    @test !DP._ready_to_optimize(model)
+    # The failed run still answers the standard result API the way a
+    # native solver does: the linear equality makes the very first
+    # master infeasible (a genuine certificate), the pool is empty, and
+    # value/bound queries throw the standard errors.
+    @test termination_status(model) == MOI.INFEASIBLE
+    @test result_count(model) == 0
+    @test !has_values(model)
+    @test occursin("no feasible incumbent", raw_status(model))
+    @test_throws MOI.ResultIndexBoundsError value(x)
+    @test_throws MOI.GetAttributeNotAllowed objective_bound(model)
+end
+
+function test_loa_hull_linear()
+    # Same GDP as test_loa_solve_simple, but inner_method = Hull so the
+    # NLP and master carry disaggregated variables and the disjunct
+    # cuts are convex-hull cuts rather than big-M. Optimum unchanged:
+    # x = 7 via the second disjunct.
+    model = GDPModel(HiGHS.Optimizer)
+    set_silent(model)
+    @variable(model, 0 <= x <= 10)
+    @variable(model, Y[1:2], Logical)
+    @constraint(model, x <= 3, Disjunct(Y[1]))
+    @constraint(model, x <= 7, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @objective(model, Max, x)
+    optimize!(model, gdp_method = LOA(HiGHS.Optimizer;
+        inner_method = Hull()))
+    @test termination_status(model) == MOI.LOCALLY_SOLVED
+    @test objective_value(model) ≈ 7.0 atol = 1e-4
+end
+
+function test_loa_hull_two_disjunctions()
+    # Two independent disjunctions under Hull; optimum x = 7, z = 5.
+    model = GDPModel(HiGHS.Optimizer)
+    set_silent(model)
+    @variable(model, 0 <= x <= 10)
+    @variable(model, 0 <= z <= 10)
+    @variable(model, Y[1:2], Logical)
+    @variable(model, W[1:2], Logical)
+    @constraint(model, x <= 3, Disjunct(Y[1]))
+    @constraint(model, x <= 7, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @constraint(model, z <= 2, Disjunct(W[1]))
+    @constraint(model, z <= 5, Disjunct(W[2]))
+    @disjunction(model, W)
+    @objective(model, Max, x + z)
+    optimize!(model, gdp_method = LOA(HiGHS.Optimizer;
+        inner_method = Hull()))
+    @test termination_status(model) == MOI.LOCALLY_SOLVED
+    @test objective_value(model) ≈ 12.0 atol = 1e-4
+end
+
+function test_loa_hull_nonlinear_global()
+    # max x s.t. x^2 <= 25 (global), (x <= 3) ∨ (x <= 8) under Hull.
+    # Linear disjuncts enter the master as exact Hull perspectives; the
+    # global nonlinear enters via global OA cuts. Optimum: x = 5.
+    ipopt = optimizer_with_attributes(Ipopt.Optimizer,
+        "print_level" => 0, "sb" => "yes")
+    juniper = optimizer_with_attributes(Juniper.Optimizer,
+        "nl_solver" => ipopt, "log_levels" => [])
+    model = GDPModel(juniper)
+    set_silent(model)
+    @variable(model, 0 <= x <= 10)
+    @constraint(model, x^2 <= 25)
+    @variable(model, Y[1:2], Logical)
+    @constraint(model, x <= 3, Disjunct(Y[1]))
+    @constraint(model, x <= 8, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @objective(model, Max, x)
+    optimize!(model, gdp_method = LOA(juniper;
+        mip_optimizer = HiGHS.Optimizer, inner_method = Hull()))
+    @test termination_status(model) in (MOI.OPTIMAL, MOI.LOCALLY_SOLVED)
+    @test objective_value(model) ≈ 5.0 atol = 1e-3
+end
+
+function test_loa_hull_nonlinear_disjunct()
+    # Nonlinear disjunct constraint under Hull: this is the case big-M
+    # and Hull genuinely differ on. 0 <= x <= 10, Y1: x <= 3,
+    # Y2: x^2 <= 64. The convex-hull OA cut disaggregates the
+    # linearization of x^2. Optimum: x = 8 via Y2.
+    ipopt = optimizer_with_attributes(Ipopt.Optimizer,
+        "print_level" => 0, "sb" => "yes")
+    juniper = optimizer_with_attributes(Juniper.Optimizer,
+        "nl_solver" => ipopt, "log_levels" => [])
+    model = GDPModel(juniper)
+    set_silent(model)
+    @variable(model, 0 <= x <= 10)
+    @variable(model, Y[1:2], Logical)
+    @constraint(model, x <= 3, Disjunct(Y[1]))
+    @constraint(model, x^2 <= 64, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @objective(model, Max, x)
+    optimize!(model, gdp_method = LOA(juniper;
+        mip_optimizer = HiGHS.Optimizer, inner_method = Hull()))
+    @test termination_status(model) in (MOI.OPTIMAL, MOI.LOCALLY_SOLVED)
+    @test objective_value(model) ≈ 8.0 atol = 1e-3
+end
+
+function test_loa_hull_complement_nonlinear()
+    # Hull with a complement indicator Y2 ≡ ¬Y1 and a nonlinear
+    # disjunct under Y2. The disaggregator is keyed by the complement
+    # binary `1 - y1` (an AffExpr); the cut must disaggregate against
+    # it. Optimum: x = 8 via Y2.
+    ipopt = optimizer_with_attributes(Ipopt.Optimizer,
+        "print_level" => 0, "sb" => "yes")
+    juniper = optimizer_with_attributes(Juniper.Optimizer,
+        "nl_solver" => ipopt, "log_levels" => [])
+    model = GDPModel(juniper)
+    set_silent(model)
+    @variable(model, 0 <= x <= 10)
+    @variable(model, Y1, Logical)
+    @variable(model, Y2, Logical, logical_complement = Y1)
+    @constraint(model, x <= 3, Disjunct(Y1))
+    @constraint(model, x^2 <= 64, Disjunct(Y2))
+    @disjunction(model, [Y1, Y2])
+    @objective(model, Max, x)
+    optimize!(model, gdp_method = LOA(juniper;
+        mip_optimizer = HiGHS.Optimizer, inner_method = Hull()))
+    @test termination_status(model) in (MOI.OPTIMAL, MOI.LOCALLY_SOLVED)
+    @test objective_value(model) ≈ 8.0 atol = 1e-3
+end
+
+function test_loa_hull_nested_disaggregation_map()
+    # The inner disjunction (over y) is a disjunct of the outer (gated by
+    # z[1]). Hull's nested redispatch must record the inner
+    # disaggregations in the GDP data map, else LOA-Hull emits the
+    # nested cut on the aggregated variable.
+    model = GDPModel()
+    @variable(model, 0 <= x <= 10)
+    @variable(model, y[1:2], Logical)
+    @variable(model, z[1:2], Logical)
+    @constraint(model, x <= 3, Disjunct(y[1]))
+    @constraint(model, x <= 8, Disjunct(y[2]))
+    @disjunction(model, y, Disjunct(z[1]))
+    @constraint(model, x <= 1, Disjunct(z[2]))
+    @disjunction(model, z)
+    DP.reformulate_model(model, Hull())
+    disaggregations = DP._disaggregation_map(model)
+    @test haskey(disaggregations, (x, y[1]))
+    @test haskey(disaggregations, (x, y[2]))
+end
+
+function test_loa_result_attributes()
+    # Post-LOA the model answers the full standard JuMP result API from
+    # the injected result cache: LOA-level status, composed solver name,
+    # master bound and gap, wall-clock time, and one solution-pool
+    # result per feasible combination (best first). Both disjuncts are
+    # nonlinear so set covering must visit both combinations, making
+    # the pool size deterministic.
+    ipopt = optimizer_with_attributes(Ipopt.Optimizer,
+        "print_level" => 0, "sb" => "yes")
+    model = GDPModel(ipopt)
+    set_silent(model)
+    @variable(model, 0 <= x <= 10)
+    @constraint(model, cap, x <= 9)
+    @variable(model, Y[1:2], Logical)
+    @constraint(model, x^2 <= 9, Disjunct(Y[1]))
+    @constraint(model, x^2 <= 49, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @objective(model, Max, x)
+    optimize!(model,
+        gdp_method = LOA(ipopt; mip_optimizer = HiGHS.Optimizer))
+
+    @test termination_status(model) == MOI.LOCALLY_SOLVED
+    @test primal_status(model) == MOI.FEASIBLE_POINT
+    @test dual_status(model) == MOI.NO_SOLUTION
+    @test_throws MOI.GetAttributeNotAllowed dual_objective_value(model)
+    @test_throws MOI.GetAttributeNotAllowed dual(cap)
+    @test startswith(solver_name(model), "LOA(BigM")
+    @test occursin("Ipopt", solver_name(model))
+    @test occursin("HiGHS", solver_name(model))
+    @test occursin("LOA finished", raw_status(model))
+    @test objective_bound(model) ≈ 7.0 atol = 1e-3
+    @test relative_gap(model) <= 1e-3
+    @test solve_time(model) >= 0.0
+    # Solution pool: both combinations feasible, best first.
+    @test result_count(model) == 2
+    @test objective_value(model) ≈ 7.0 atol = 1e-3
+    @test value(x) ≈ 7.0 atol = 1e-3
+    @test value(Y[2]) ≈ 1.0 atol = 1e-6
+    @test primal_status(model; result = 2) == MOI.FEASIBLE_POINT
+    @test objective_value(model; result = 2) ≈ 3.0 atol = 1e-3
+    @test value(x; result = 2) ≈ 3.0 atol = 1e-3
+    @test primal_status(model; result = 3) == MOI.NO_SOLUTION
+    # ConstraintPrimal is served by evaluating at the injected primal.
+    @test value(cap) ≈ 7.0 atol = 1e-3
+    # solution_summary exercises nearly every query in one call.
+    summary = sprint(show, solution_summary(model))
+    @test occursin("LOA(", summary)
+    @test occursin("LOCALLY_SOLVED", summary)
+    # The binaries come back binary, unfixed, and bound-free.
+    bins = [v for (_, v) in DP._indicator_to_binary(model)
+        if v isa JuMP.VariableRef]
+    @test all(JuMP.is_binary, bins)
+    @test all(!JuMP.is_fixed, bins)
+    @test all(!JuMP.has_lower_bound, bins)
+end
+
+function test_reformulate_after_loa_restores_binaries()
+    # LOA relaxes the logical binaries and fixes them per iteration; at
+    # exit it must hand them back binary and unfixed so a later
+    # reformulation of the same model works unchanged. Regression: a
+    # nested disjunction's tight-M looked up a relaxed binary's bounds
+    # and threw a KeyError.
+    model = GDPModel(HiGHS.Optimizer)
+    set_silent(model)
+    @variable(model, 1 <= x <= 9)
+    @variable(model, Y[1:2], Logical)
+    @variable(model, W[1:2], Logical)
+    @constraint(model, 1 <= x <= 2, Disjunct(W[1]))
+    @constraint(model, 2 <= x <= 3, Disjunct(W[2]))
+    @constraint(model, x >= 8, Disjunct(Y[2]))
+    @disjunction(model, [W[1], W[2]], Disjunct(Y[1]))
+    @disjunction(model, [Y[1], Y[2]])
+    @objective(model, Max, x)
+
+    optimize!(model, gdp_method = LOA(HiGHS.Optimizer))
+    bins = [v for (_, v) in DP._indicator_to_binary(model)
+        if v isa JuMP.VariableRef]
+    @test all(JuMP.is_binary, bins)           # restored at LOA exit
+    @test all(!JuMP.is_fixed, bins)
+
+    @test optimize!(model, gdp_method = BigM()) isa Nothing
+    @test termination_status(model) == MOI.OPTIMAL
+    @test objective_value(model) ≈ 9
+end
+
+@testset "LOA" begin
+    test_loa_datatype()
+    test_cover_disjuncts()
+    test_oa_cut_terms()
+    test_is_linear_F()
+    test_no_good_cut()
+    test_loa_hook_simple()
+    test_loa_reoptimize()
+    test_loa_reoptimize_respects_user_optimizer()
+    test_loa_solve_simple()
+    test_loa_solve_simple_with_mbm()
+    test_loa_solve_two_disjunctions()
+    test_loa_not_reformulation()
+    test_loa_nonlinear_global()
+    test_loa_nonlinear_equality_global()
+    test_loa_nonlinear_equality_disjunct()
+    test_loa_nonlinear_interval_disjunct()
+    test_loa_restores_prior_time_limit()
+    test_loa_restores_displaced_optimizer()
+    test_loa_repeated_solve_keeps_displaced_optimizer()
+    test_loa_no_optimizer_to_restore()
+    test_loa_cache_refuses_direct_solve()
+    test_loa_limit_hit_report()
+    test_loa_complement_indicator_nonlinear_disjunct()
+    test_loa_nlpf_infeasible_disjunct()
+    test_loa_nlpf_disabled()
+    test_loa_solve_nlp_infeasible_fallback()
+    test_loa_sense_primitives()
+    test_loa_iteration_loop()
+    test_loa_time_limits()
+    test_loa_no_feasible_incumbent()
+    test_loa_hull_linear()
+    test_loa_hull_two_disjunctions()
+    test_loa_hull_nonlinear_global()
+    test_loa_hull_nonlinear_disjunct()
+    test_loa_hull_complement_nonlinear()
+    test_loa_hull_nested_disaggregation_map()
+    test_loa_result_attributes()
+    test_reformulate_after_loa_restores_binaries()
+end
diff --git a/test/constraints/mbm.jl b/test/constraints/mbm.jl
index 11f164da..ab9e8287 100644
--- a/test/constraints/mbm.jl
+++ b/test/constraints/mbm.jl
@@ -52,7 +52,7 @@ function test__replace_variables_quad_numeric_map()
     @test result3.aff.terms[y] ≈ 3.0
 end
 
-function test_replace_variables_in_constraint()
+function testreplace_variables_in_constraint()
     model = Model()
     sub_model = Model()
     @variable(model, x[1:3])
@@ -809,7 +809,7 @@ end
     test_mbm()
     test__var_ref_type_numeric_map()
     test__replace_variables_quad_numeric_map()
-    test_replace_variables_in_constraint()
+    testreplace_variables_in_constraint()
     test_prepare_max_M_objective()
     test_raw_M()
     test_maximize_M()
diff --git a/test/extensions/InfiniteDisjunctiveProgramming.jl b/test/extensions/InfiniteDisjunctiveProgramming.jl
index 47bb884b..5fe4186a 100644
--- a/test/extensions/InfiniteDisjunctiveProgramming.jl
+++ b/test/extensions/InfiniteDisjunctiveProgramming.jl
@@ -513,6 +513,76 @@ function test_mbm_finite_and_integer_var()
         [MOI.OPTIMAL, MOI.LOCALLY_SOLVED]
 end
 
+function test_mbm_infinite_simple()
+    model = InfiniteGDPModel(HiGHS.Optimizer)
+    set_silent(model)
+    model = InfiniteGDPModel(HiGHS.Optimizer)
+    set_silent(model)
+    K = 4
+    @infinite_parameter(model, t ∈ [0, 1], num_supports = K)
+    @variable(model, 0 <= x <= 10, Infinite(t))
+    @variable(model, Y[1:2], InfiniteLogical(t))
+    @constraint(model, x >= 5, Disjunct(Y[1]))
+    @constraint(model, x <= 3, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    JuMP.fix(Y[2], true)  # force disj 2 active
+    @objective(model, Min, ∫(x, t))
+    DP.reformulate_model(model, BigM(10.0))
+    set_optimizer(model, HiGHS.Optimizer)
+    set_silent(model)
+    optimize!(model, ignore_optimize_hook = true)
+    sol = DP.extract_solution(model)
+    @test haskey(sol, x)
+    @test length(sol[x]) == K
+    @test all(v -> isapprox(v, 0.0; atol=1e-6), sol[x])
+end
+
+# add_cut adds one pointwise-sum cut to the transformation backend and
+# marks the backend ready so the next optimize! does NOT re-transcribe.
+function test_add_cut_infinite()
+    model = InfiniteGDPModel(HiGHS.Optimizer)
+    set_silent(model)
+    K = 3
+    @infinite_parameter(model, t ∈ [0, 1], num_supports = K)
+    @variable(model, 0 <= x <= 10, Infinite(t))
+    @variable(model, Y[1:2], InfiniteLogical(t))
+    @constraint(model, x >= 5, Disjunct(Y[1]))
+    @constraint(model, x <= 3, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    DP.reformulate_model(model, BigM(10.0))
+    InfiniteOpt.build_transformation_backend!(model)
+    transcribed = InfiniteOpt.transformation_model(model)
+    n_before = JuMP.num_constraints(transcribed;
+        count_variable_in_set_constraints = false)
+    rBM_sol = Dict(x => [1.0, 2.0, 3.0])
+    sep_sol = Dict(x => [0.5, 1.5, 2.5])
+    DP.add_cut(model, [x], rBM_sol, sep_sol)
+    n_after = JuMP.num_constraints(transcribed;
+        count_variable_in_set_constraints = false)
+    @test n_after == n_before + 1
+    # set_transformation_backend_ready(true) — next optimize! should
+    # reuse without re-transcribing (otherwise our cut would be lost)
+    @test InfiniteOpt.transformation_backend_ready(model)
+end
+
+# MBM with finite + integer variables in an InfiniteModel.
+function test_mbm_finite_and_integer_var()
+    model = InfiniteGDPModel(HiGHS.Optimizer)
+    set_silent(model)
+    @infinite_parameter(model, t ∈ [0, 1], num_supports = 10)
+    @variable(model, 0 <= x <= 10, Infinite(t))
+    @variable(model, 0 <= w <= 5, Int)
+    @variable(model, Y[1:2], InfiniteLogical(t))
+    @constraint(model, x + w >= 5, Disjunct(Y[1]))
+    @constraint(model, x + w <= 3, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @objective(model, Min, ∫(x, t) + w)
+    @test optimize!(model,
+        gdp_method = MBM(HiGHS.Optimizer)) isa Nothing
+    @test termination_status(model) in
+        [MOI.OPTIMAL, MOI.LOCALLY_SOLVED]
+end
+
 function test_mbm_infinite_simple()
     model = InfiniteGDPModel(HiGHS.Optimizer)
     set_silent(model)
@@ -782,13 +852,353 @@ function test_CuttingPlanes_multiparameter()
 
     @objective(model, Min, ∫(∫(x, t), s))
 
-    # Should not throw
-    @test optimize!(model,
-        gdp_method = CuttingPlanes(
-            HiGHS.Optimizer; max_iter = 5)
-    ) isa Nothing
+    optimize!(model,
+        gdp_method = CuttingPlanes(HiGHS.Optimizer; max_iter = 5))
     @test termination_status(model) in
         [MOI.OPTIMAL, MOI.LOCALLY_SOLVED]
+    # Y[2] (x <= 3) is the minimizer: x = 0 across (t, s) → ∫∫x = 0.
+    @test objective_value(model) ≈ 0.0 atol = 1e-4
+    @test all(isapprox.(value(x), 0.0; atol = 1e-6))
+end
+
+function test_loa_infinite_nonlinear_global()
+    # max ∫x dt s.t. x(t)^2 <= 25 (global, per-support nonlinear),
+    #   (x <= 3) ∨ (x <= 8), 0 <= x <= 10 over t ∈ [0, 1].
+    # Disjunct Y[2] permits x up to 8 but the global x^2 <= 25 caps
+    # x at 5. The per-support global transcribes to an `AbstractArray`
+    # of scalar constraints, so this exercises the per-support fan-out
+    # of the transcribed global constraints. Without the global cut the
+    # master would allow x = 8 and report 8.0; the binding optimum
+    # is ∫5 dt = 5.
+    ipopt = optimizer_with_attributes(Ipopt.Optimizer,
+        "print_level" => 0, "sb" => "yes")
+    model = InfiniteGDPModel(ipopt)
+    set_silent(model)
+    @infinite_parameter(model, t ∈ [0, 1], num_supports = 10)
+    @variable(model, 0 <= x <= 10, Infinite(t))
+    @constraint(model, x^2 <= 25)
+    @variable(model, Y[1:2], InfiniteLogical(t))
+    @constraint(model, x <= 3, Disjunct(Y[1]))
+    @constraint(model, x <= 8, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @objective(model, Max, ∫(x, t))
+    optimize!(model,
+        gdp_method = LOA(ipopt; mip_optimizer = HiGHS.Optimizer))
+    @test termination_status(model) in
+        (MOI.OPTIMAL, MOI.LOCALLY_SOLVED)
+    @test objective_value(model) ≈ 5.0 atol = 1e-2
+    @test all(value(Y[2]))
+    @test all(isapprox.(value(x), 5.0; atol = 1e-2))
+end
+
+function test_loa_infinite_complement_indicator()
+    # Regression: `_indicator_to_binary(model)[Y2]` returns the AffExpr
+    # `1 - binary(Y1)` for a logical-complement indicator. Fixing and
+    # gating must resolve it onto the underlying binary's transcription
+    # with the value inverted.
+    #
+    # max ∫ x dt with 0 ≤ x ≤ 10 over t ∈ [0,1]:
+    #   Y1: x ≤ 3       Y2 ≡ ¬Y1: x ≤ 8
+    # Y2 is the maximizer (x = 8 across t), giving ∫8 dt = 8.
+    ipopt = optimizer_with_attributes(Ipopt.Optimizer,
+        "print_level" => 0, "sb" => "yes")
+    model = InfiniteGDPModel(ipopt)
+    set_silent(model)
+    @infinite_parameter(model, t ∈ [0, 1], num_supports = 5)
+    @variable(model, 0 <= x <= 10, Infinite(t))
+    @variable(model, Y1, Logical)
+    @variable(model, Y2, Logical, logical_complement = Y1)
+    @constraint(model, x <= 3, Disjunct(Y1))
+    @constraint(model, x <= 8, Disjunct(Y2))
+    @disjunction(model, [Y1, Y2])
+    @objective(model, Max, ∫(x, t))
+    optimize!(model,
+        gdp_method = LOA(ipopt; mip_optimizer = HiGHS.Optimizer))
+    @test termination_status(model) in (MOI.OPTIMAL, MOI.LOCALLY_SOLVED)
+    @test objective_value(model) ≈ 8.0 atol = 1e-2
+end
+
+function test_loa_infinite_complement_nonlinear_disjunct()
+    # Regression: a finite (complement-form AffExpr) indicator gating a
+    # NONLINEAR constraint on an infinite variable. The constraint
+    # transcribes to one scalar row per support while the binary
+    # reference is a single expression, so the fan-out must broadcast
+    # it across the constraint's supports.
+    #
+    # max ∫ x dt with 0 ≤ x ≤ 10 over t ∈ [0,1]:
+    #   Y1: x ≤ 3        Y2 ≡ ¬Y1: x² ≤ 64
+    # Y2 (complement) is the maximizer (x = 8 across t), giving 8.
+    ipopt = optimizer_with_attributes(Ipopt.Optimizer,
+        "print_level" => 0, "sb" => "yes")
+    juniper = optimizer_with_attributes(Juniper.Optimizer,
+        "nl_solver" => ipopt, "log_levels" => [])
+    model = InfiniteGDPModel(juniper)
+    set_silent(model)
+    @infinite_parameter(model, t ∈ [0, 1], num_supports = 5)
+    @variable(model, 0 <= x <= 10, Infinite(t))
+    @variable(model, Y1, Logical)
+    @variable(model, Y2, Logical, logical_complement = Y1)
+    @constraint(model, x <= 3, Disjunct(Y1))
+    @constraint(model, x^2 <= 64, Disjunct(Y2))
+    @disjunction(model, [Y1, Y2])
+    @objective(model, Max, ∫(x, t))
+    optimize!(model,
+        gdp_method = LOA(juniper; mip_optimizer = HiGHS.Optimizer))
+    @test termination_status(model) in (MOI.OPTIMAL, MOI.LOCALLY_SOLVED)
+    @test objective_value(model) ≈ 8.0 atol = 1e-2
+end
+
+function test_loa_infinite_multidim_parameter()
+    # Regression: multi-D dependent parameter group. Transcription
+    # expands a variable over ξ[1:2] across its joint supports; the
+    # transcribed binaries and constraint rows must stay aligned.
+    #
+    # max w with 0 ≤ x ≤ 10 over ξ[1:2] ∈ [0,1]², w ≤ x at every joint
+    # support:
+    #   Y[1]: x ≤ 3       Y[2]: x ≤ 8
+    # Y[2] is the maximizer (x = 8 across ξ allows w = 8).
+    ipopt = optimizer_with_attributes(Ipopt.Optimizer,
+        "print_level" => 0, "sb" => "yes")
+    model = InfiniteGDPModel(ipopt)
+    set_silent(model)
+    @infinite_parameter(model, ξ[1:2] ∈ [0, 1], num_supports = 4)
+    @variable(model, 0 <= x <= 10, Infinite(ξ))
+    @variable(model, 0 <= w <= 10)
+    @constraint(model, w <= x)
+    @variable(model, Y[1:2], InfiniteLogical(ξ))
+    @constraint(model, x <= 3, Disjunct(Y[1]))
+    @constraint(model, x <= 8, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @objective(model, Max, w)
+    optimize!(model,
+        gdp_method = LOA(ipopt; mip_optimizer = HiGHS.Optimizer))
+    @test termination_status(model) in (MOI.OPTIMAL, MOI.LOCALLY_SOLVED)
+    @test objective_value(model) ≈ 8.0 atol = 1e-2
+end
+
+function test_loa_infinite_nlpf_infeasible_disjunct()
+    # InfiniteOpt LOA: Y[1]'s per-support constraint x >= 200 is
+    # NLP-infeasible against the bound x in [0, 10]. With infinite
+    # indicators `Y[1:2], InfiniteLogical(t)`, the infeasible
+    # Y[1]-everywhere seed makes the primary NLP fail, so NLPF runs: it
+    # copies the transcribed NLP (per-support binary fixes included),
+    # slacks the inequality, and returns a linearization point. The
+    # disjunct is deliberately LINEAR:
+    # the model is then convex, so LOA's master bound is rigorous and
+    # it converges to Y[2] active everywhere (x = 5), giving ∫x dt = 5.
+    # (A nonconvex infeasible disjunct such as x^2 >= 200 would make LOA
+    # a heuristic and let the local NLP solver report a constraint-
+    # violating point as solved — out of scope for the NLPF path here.)
+    ipopt = optimizer_with_attributes(Ipopt.Optimizer,
+        "print_level" => 0, "sb" => "yes")
+    juniper = optimizer_with_attributes(Juniper.Optimizer,
+        "nl_solver" => ipopt, "log_levels" => [])
+    model = InfiniteGDPModel(juniper)
+    set_silent(model)
+    @infinite_parameter(model, t ∈ [0, 1], num_supports = 3)
+    @variable(model, 0 <= x <= 10, Infinite(t))
+    @variable(model, Y[1:2], InfiniteLogical(t))
+    @constraint(model, x >= 200, Disjunct(Y[1]))
+    @constraint(model, x <= 5, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @objective(model, Max, ∫(x, t))
+    optimize!(model,
+        gdp_method = LOA(juniper; mip_optimizer = HiGHS.Optimizer))
+    @test termination_status(model) in
+        (MOI.OPTIMAL, MOI.LOCALLY_SOLVED)
+    @test objective_value(model) ≈ 5.0 atol = 1e-2
+end
+
+function test_loa_infinite_aggregate_global()
+    # min ∫x dt s.t. ∫(x^2, t) <= 4 (aggregate global), x >= y,
+    #   (y >= 1) ∨ (y >= 3), 0 <= x, y <= 10 over t ∈ [0, 1].
+    # The aggregate global transcribes to a single scalar (the
+    # measure is expanded), exercising the single-row branch of
+    # the transcribed global constraints. Y[1] (y >= 1) is the cheaper
+    # disjunct: x = 1 satisfies x >= y and ∫x^2 = 1 <= 4, giving
+    # objective ∫1 dt = 1.
+    ipopt = optimizer_with_attributes(Ipopt.Optimizer,
+        "print_level" => 0, "sb" => "yes")
+    model = InfiniteGDPModel(ipopt)
+    set_silent(model)
+    @infinite_parameter(model, t ∈ [0, 1], num_supports = 10)
+    @variable(model, 0 <= x <= 10, Infinite(t))
+    @variable(model, 0 <= y <= 10)
+    @constraint(model, ∫(x^2, t) <= 4)
+    @constraint(model, x >= y)
+    @variable(model, Y[1:2], Logical)
+    @constraint(model, y >= 1, Disjunct(Y[1]))
+    @constraint(model, y >= 3, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @objective(model, Min, ∫(x, t))
+    optimize!(model,
+        gdp_method = LOA(ipopt; mip_optimizer = HiGHS.Optimizer))
+    @test termination_status(model) in
+        (MOI.OPTIMAL, MOI.LOCALLY_SOLVED)
+    @test objective_value(model) ≈ 1.0 atol = 1e-2
+end
+
+function test_loa_infinite_aggregate_disjunct()
+    # An aggregate (measure) constraint INSIDE a disjunct: it
+    # transcribes to a single scalar row that must be gated by the
+    # indicator's binary. max ∫x dt with Y[1]: ∫(x^2, t) <= 4 and
+    # Y[2]: x <= 0. Y[1] is optimal: constant x = 2 gives ∫x^2 = 4
+    # and ∫x dt = 2 (Y[2] caps the objective at 0). Big-M tightening
+    # is disabled: it needs bounds of every constraint variable and a
+    # `MeasureRef` has none.
+    ipopt = optimizer_with_attributes(Ipopt.Optimizer,
+        "print_level" => 0, "sb" => "yes")
+    model = InfiniteGDPModel(ipopt)
+    set_silent(model)
+    @infinite_parameter(model, t ∈ [0, 1], num_supports = 10)
+    @variable(model, 0 <= x <= 10, Infinite(t))
+    @variable(model, Y[1:2], Logical)
+    @constraint(model, ∫(x^2, t) <= 4, Disjunct(Y[1]))
+    @constraint(model, x <= 0, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @objective(model, Max, ∫(x, t))
+    optimize!(model, gdp_method = LOA(ipopt;
+        mip_optimizer = HiGHS.Optimizer, inner_method = BigM(1e4, false)))
+    @test termination_status(model) in (MOI.OPTIMAL, MOI.LOCALLY_SOLVED)
+    @test objective_value(model) ≈ 2.0 atol = 1e-2
+    @test value(Y[1])
+end
+
+function test_loa_infinite_iteration_loop()
+    # Force the InfiniteModel LOA main loop to fix combinations that
+    # vary per support. Two disjunctions over disjoint x/z half-lines
+    # put the optimum on an off-diagonal combination the set-covering
+    # seeds never try, so the master stays feasible and the loop body
+    # runs with master-extracted per-support values.
+    #   min ∫(x - z) dt
+    #   D1: (x <= 4) [Y1]  ∨  (x >= 6) [Y2]
+    #   D2: (z <= 4) [W1]  ∨  (z >= 6) [W2]
+    # Optimum: Y1 & W2 everywhere → x = 0, z = 10, ∫(x - z) = -10.
+    model = InfiniteGDPModel(HiGHS.Optimizer)
+    set_silent(model)
+    @infinite_parameter(model, t ∈ [0, 1], num_supports = 3)
+    @variable(model, 0 <= x <= 10, Infinite(t))
+    @variable(model, 0 <= z <= 10, Infinite(t))
+    @variable(model, Y[1:2], InfiniteLogical(t))
+    @variable(model, W[1:2], InfiniteLogical(t))
+    @constraint(model, x <= 4, Disjunct(Y[1]))
+    @constraint(model, x >= 6, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @constraint(model, z <= 4, Disjunct(W[1]))
+    @constraint(model, z >= 6, Disjunct(W[2]))
+    @disjunction(model, W)
+    @objective(model, Min, ∫(x - z, t))
+    optimize!(model, gdp_method = LOA(HiGHS.Optimizer))
+    @test termination_status(model) in
+        (MOI.OPTIMAL, MOI.LOCALLY_SOLVED)
+    @test objective_value(model) ≈ -10.0 atol = 1e-3
+    @test all(value(Y[1]))
+    @test all(value(W[2]))
+    @test all(isapprox.(value(x), 0.0; atol = 1e-6))
+    @test all(isapprox.(value(z), 10.0; atol = 1e-6))
+end
+
+function test_loa_infinite_hull_linear()
+    # InfiniteOpt LOA with inner_method = Hull. Linear disjuncts enter
+    # the master as exact per-support Hull perspectives over
+    # disaggregated infinite variables. max ∫x dt, Y1: x <= 3,
+    # Y2: x <= 8. Optimum: x = 8 across t → ∫8 dt = 8.
+    model = InfiniteGDPModel(HiGHS.Optimizer)
+    set_silent(model)
+    @infinite_parameter(model, t ∈ [0, 1], num_supports = 5)
+    @variable(model, 0 <= x <= 10, Infinite(t))
+    @variable(model, Y[1:2], InfiniteLogical(t))
+    @constraint(model, x <= 3, Disjunct(Y[1]))
+    @constraint(model, x <= 8, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @objective(model, Max, ∫(x, t))
+    optimize!(model, gdp_method = LOA(HiGHS.Optimizer;
+        inner_method = Hull()))
+    @test termination_status(model) in (MOI.OPTIMAL, MOI.LOCALLY_SOLVED)
+    @test objective_value(model) ≈ 8.0 atol = 1e-2
+end
+
+function test_loa_infinite_hull_nonlinear_disjunct()
+    # The infinite Hull case big-M differs from: a NONLINEAR constraint
+    # on an infinite variable inside a disjunct, under a plain infinite
+    # indicator. The convex-hull OA cut is built per support over the
+    # transcribed disaggregated infinite variable. max ∫x dt,
+    # Y1: x <= 3, Y2: x^2 <= 64. Optimum: x = 8 via Y2.
+    ipopt = optimizer_with_attributes(Ipopt.Optimizer,
+        "print_level" => 0, "sb" => "yes")
+    juniper = optimizer_with_attributes(Juniper.Optimizer,
+        "nl_solver" => ipopt, "log_levels" => [])
+    model = InfiniteGDPModel(juniper)
+    set_silent(model)
+    @infinite_parameter(model, t ∈ [0, 1], num_supports = 5)
+    @variable(model, 0 <= x <= 10, Infinite(t))
+    @variable(model, Y[1:2], InfiniteLogical(t))
+    @constraint(model, x <= 3, Disjunct(Y[1]))
+    @constraint(model, x^2 <= 64, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @objective(model, Max, ∫(x, t))
+    optimize!(model, gdp_method = LOA(juniper;
+        mip_optimizer = HiGHS.Optimizer, inner_method = Hull()))
+    @test termination_status(model) in (MOI.OPTIMAL, MOI.LOCALLY_SOLVED)
+    @test objective_value(model) ≈ 8.0 atol = 1e-2
+    @test all(value(Y[2]))
+    @test all(isapprox.(value(x), 8.0; atol = 1e-2))
+end
+
+function test_loa_infinite_hull_complement_nonlinear()
+    # Hull with a complement indicator Y2 ≡ ¬Y1 gating a nonlinear
+    # constraint on an infinite variable. The per-support disaggregator
+    # is keyed by the transcribed complement binary `1 - y1_k`
+    # (an AffExpr). max ∫x dt, Y1: x <= 3, Y2: x^2 <= 64. Optimum 8.
+    ipopt = optimizer_with_attributes(Ipopt.Optimizer,
+        "print_level" => 0, "sb" => "yes")
+    juniper = optimizer_with_attributes(Juniper.Optimizer,
+        "nl_solver" => ipopt, "log_levels" => [])
+    model = InfiniteGDPModel(juniper)
+    set_silent(model)
+    @infinite_parameter(model, t ∈ [0, 1], num_supports = 5)
+    @variable(model, 0 <= x <= 10, Infinite(t))
+    @variable(model, Y1, Logical)
+    @variable(model, Y2, Logical, logical_complement = Y1)
+    @constraint(model, x <= 3, Disjunct(Y1))
+    @constraint(model, x^2 <= 64, Disjunct(Y2))
+    @disjunction(model, [Y1, Y2])
+    @objective(model, Max, ∫(x, t))
+    optimize!(model, gdp_method = LOA(juniper;
+        mip_optimizer = HiGHS.Optimizer, inner_method = Hull()))
+    @test termination_status(model) in (MOI.OPTIMAL, MOI.LOCALLY_SOLVED)
+    @test objective_value(model) ≈ 8.0 atol = 1e-2
+end
+
+function test_loa_infinite_hull_finite_var()
+    # A FINITE variable disaggregated inside a nonlinear disjunct
+    # constraint gated by an INFINITE indicator. The fan-out slices the
+    # binary per support, so the disaggregator must key the finite
+    # variable's single copy by each per-support binary (driven off the
+    # binary, not the variable). max ∫x dt + w over t ∈ [0,1]:
+    #   Y1: x^2 + w^2 <= 25      Y2: x <= 0, w <= 3
+    # Y1 is active (Y2 caps the objective at 3); x(t) = w = sqrt(12.5),
+    # so the optimum is ∫sqrt(12.5) dt + sqrt(12.5) = 2 sqrt(12.5)
+    # ≈ 7.0711.
+    ipopt = optimizer_with_attributes(Ipopt.Optimizer,
+        "print_level" => 0, "sb" => "yes")
+    juniper = optimizer_with_attributes(Juniper.Optimizer,
+        "nl_solver" => ipopt, "log_levels" => [])
+    model = InfiniteGDPModel(juniper)
+    set_silent(model)
+    @infinite_parameter(model, t ∈ [0, 1], num_supports = 4)
+    @variable(model, 0 <= x <= 10, Infinite(t))
+    @variable(model, 0 <= w <= 10)
+    @variable(model, Y[1:2], InfiniteLogical(t))
+    @constraint(model, x^2 + w^2 <= 25, Disjunct(Y[1]))
+    @constraint(model, x <= 0, Disjunct(Y[2]))
+    @constraint(model, w <= 3, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @objective(model, Max, ∫(x, t) + w)
+    optimize!(model, gdp_method = LOA(juniper;
+        mip_optimizer = HiGHS.Optimizer, inner_method = Hull()))
+    @test termination_status(model) in (MOI.OPTIMAL, MOI.LOCALLY_SOLVED)
+    @test objective_value(model) ≈ 7.0711 atol = 1e-2
 end
 
 @testset "InfiniteDisjunctiveProgramming" begin
@@ -855,4 +1265,22 @@ end
         test_CuttingPlanes_multiparameter()
     end
 
+    # LOA runs on the transcribed backend, a plain JuMP model, so it
+    # needs no `JuMP.copy_model(::InfiniteModel)` from the InfiniteOpt
+    # fork (unlike MBM).
+    @testset "LOA" begin
+        test_loa_infinite_nonlinear_global()
+        test_loa_infinite_aggregate_global()
+        test_loa_infinite_aggregate_disjunct()
+        test_loa_infinite_complement_indicator()
+        test_loa_infinite_complement_nonlinear_disjunct()
+        test_loa_infinite_multidim_parameter()
+        test_loa_infinite_nlpf_infeasible_disjunct()
+        test_loa_infinite_iteration_loop()
+        test_loa_infinite_hull_linear()
+        test_loa_infinite_hull_nonlinear_disjunct()
+        test_loa_infinite_hull_complement_nonlinear()
+        test_loa_infinite_hull_finite_var()
+    end
+
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index 06e8813d..ca34f816 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -18,6 +18,7 @@ include("constraints/mbm.jl")
 include("constraints/bigm.jl")
 include("constraints/psplit.jl")
 include("constraints/cuttingplanes.jl")
+include("constraints/loa.jl")
 include("constraints/hull.jl")
 include("constraints/fallback.jl")    
 include("constraints/disjunction.jl")
diff --git a/test/solve.jl b/test/solve.jl
index c3ca9441..ee89e6e0 100644
--- a/test/solve.jl
+++ b/test/solve.jl
@@ -62,7 +62,15 @@ function test_linear_gdp_example(m, use_complements = false)
     @test value(Y[2])
     @test !value(W[1])
     @test !value(W[2])
-  
+
+    @test optimize!(m, gdp_method = LOA(HiGHS.Optimizer)) isa Nothing
+    @test termination_status(m) == MOI.LOCALLY_SOLVED
+    @test objective_value(m) ≈ 11 atol=1e-3
+    @test value.(x) ≈ [9,2] atol=1e-3
+    @test !value(Y[1])
+    @test value(Y[2])
+    @test !value(W[1])
+    @test !value(W[2])
 
     m_copy, ref_map = JuMP.copy_model(m)
     lv_map = DP.copy_gdp_data(m, m_copy, ref_map)
@@ -130,11 +138,20 @@ function test_quadratic_gdp_example(use_complements = false) #psplit does not wo
 
     @test optimize!(m, gdp_method = PSplit(2,m)) isa Nothing
     @test termination_status(m) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED]
-    @test objective_value(m) ≈ 6.1237 atol=1e-3  
-    @test value.(x) ≈ [4.0825, 2.0412] atol=1e-3 
-    @test !value(Y[1]) 
+    @test objective_value(m) ≈ 6.1237 atol=1e-3
+    @test value.(x) ≈ [4.0825, 2.0412] atol=1e-3
+    @test !value(Y[1])
     @test value(Y[2])
-    @test !value(W[1]) 
+    @test !value(W[1])
+    @test !value(W[2])
+
+    @test optimize!(m, gdp_method = LOA(optimizer)) isa Nothing
+    @test termination_status(m) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED]
+    @test objective_value(m) ≈ 6.1237 atol=1e-3
+    @test value.(x) ≈ [4.0825, 2.0412] atol=1e-3
+    @test !value(Y[1])
+    @test value(Y[2])
+    @test !value(W[1])
     @test !value(W[2])
 end
 
@@ -158,6 +175,12 @@ function test_generic_model(m)
     # TODO add meaningful tests to check the constraints/variables
 end
 
+function test_unsupported_solution_method()
+    m = GDPModel(HiGHS.Optimizer)
+    @test_throws ErrorException optimize!(
+        m, gdp_method = DummySolutionMethod())
+end
+
 @testset "Solve Linear GDP" begin
     test_linear_gdp_example(GDPModel(HiGHS.Optimizer))
     test_linear_gdp_example(GDPModel(HiGHS.Optimizer), true)
@@ -167,4 +190,5 @@ end
         )
     test_quadratic_gdp_example()
     test_generic_model(GDPModel{Float32}(mockoptimizer))
+    test_unsupported_solution_method()
 end
diff --git a/test/utilities.jl b/test/utilities.jl
index 2b87e6e1..aa12a1cd 100644
--- a/test/utilities.jl
+++ b/test/utilities.jl
@@ -59,6 +59,7 @@ end
 # Prepare helpful test types 
 struct BadVarRef <: JuMP.AbstractVariableRef end
 struct DummyReformulation <: AbstractReformulationMethod end
+struct DummySolutionMethod <: AbstractSolutionMethod end
 
 # Define types/methods to test using DP with extension models
 struct MyVar{I} <: JuMP.AbstractVariable
@@ -238,6 +239,125 @@ function test_get_constant_variable()
     @test DP.get_constant(x) == 0.0
 end
 
+function test_linearize_nonlinear_exp()
+    # exp(x) + y at (1, 2):
+    # f = e + 2, ∇f = [e, 1]
+    # linear: e*(x-1) + 1*(y-2) + (e+2) = e*x + y
+    model = GDPModel()
+    @variable(model, x)
+    @variable(model, y)
+    func = @expression(model, exp(x) + y)
+    xk = Dict{JuMP.AbstractVariableRef, Float64}(
+        x => 1.0, y => 2.0)
+    id_map = Dict(x => x, y => y)
+    lin = DP._linearize_at(func, xk, id_map)
+    @test JuMP.constant(lin) ≈ 0.0 atol = 1e-8
+    @test JuMP.coefficient(lin, x) ≈ exp(1.0) atol = 1e-8
+    @test JuMP.coefficient(lin, y) ≈ 1.0 atol = 1e-8
+end
+
+function test_linearize_nonlinear_sin()
+    # sin(x) at x = π/6:
+    # f = 0.5, f' = cos(π/6) = √3/2
+    # linear: 0.5 + (√3/2)(x - π/6)
+    model = GDPModel()
+    @variable(model, x)
+    func = @expression(model, sin(x))
+    xk = Dict{JuMP.AbstractVariableRef, Float64}(
+        x => π / 6)
+    id_map = Dict(x => x)
+    lin = DP._linearize_at(func, xk, id_map)
+    expected_const = 0.5 - (√3 / 2) * (π / 6)
+    @test JuMP.constant(lin) ≈ expected_const atol = 1e-8
+    @test JuMP.coefficient(lin, x) ≈ √3 / 2 atol = 1e-8
+end
+
+function test_linearize_nonlinear_multivar()
+    # exp(x) * sin(y) at (1, π/2):
+    # f = e*1 = e, ∂f/∂x = e*sin(π/2) = e, ∂f/∂y = e*cos(π/2) = 0
+    # linear: e + e*(x-1) + 0*(y-π/2) = e*x
+    model = GDPModel()
+    @variable(model, x)
+    @variable(model, y)
+    func = @expression(model, exp(x) * sin(y))
+    xk = Dict{JuMP.AbstractVariableRef, Float64}(
+        x => 1.0, y => π / 2)
+    id_map = Dict(x => x, y => y)
+    lin = DP._linearize_at(func, xk, id_map)
+    @test JuMP.constant(lin) ≈ 0.0 atol = 1e-8
+    @test JuMP.coefficient(lin, x) ≈ exp(1.0) atol = 1e-8
+    @test JuMP.coefficient(lin, y) ≈ 0.0 atol = 1e-8
+end
+
+function test_to_nlp_expr()
+    model = GDPModel()
+    @variable(model, x)
+    @variable(model, y)
+    idx = Dict(x => 1, y => 2)
+
+    # NonlinearExpr
+    nl = @expression(model, exp(x))
+    e = DP._to_nlp_expr(nl, idx)
+    @test e == Expr(:call, :exp, MOI.VariableIndex(1))
+
+    # AffExpr
+    aff = @expression(model, 2x + 3y + 1)
+    e = DP._to_nlp_expr(aff, idx)
+    @test e isa Expr
+    @test e.head == :call && e.args[1] == :+
+
+    # Number
+    @test DP._to_nlp_expr(42, idx) == 42
+end
+
+function test_fix_unfix_indicator()
+    # `fix_indicator` drives the backing binary: a plain indicator fixes
+    # it at the value, a complement indicator fixes the underlying
+    # binary at the inverse. `unfix_indicator` undoes both forms.
+    model = GDPModel()
+    @variable(model, x)
+    @variable(model, Y1, Logical)
+    @variable(model, Y2, Logical, logical_complement = Y1)
+    @constraint(model, x <= 3, Disjunct(Y1))
+    @constraint(model, x >= 5, Disjunct(Y2))
+    @disjunction(model, [Y1, Y2])
+    DP.reformulate_model(model, BigM())
+    y1 = DP._indicator_to_binary(model)[Y1]
+
+    DP.fix_indicator(model, Y1, true)
+    @test JuMP.is_fixed(y1) && JuMP.fix_value(y1) == 1.0
+    DP.fix_indicator(model, Y1, false)
+    @test JuMP.fix_value(y1) == 0.0
+    DP.unfix_indicator(model, Y1)
+    @test !JuMP.is_fixed(y1)
+
+    DP.fix_indicator(model, Y2, true)     # complement: underlying -> 0
+    @test JuMP.fix_value(y1) == 0.0
+    DP.fix_indicator(model, Y2, false)
+    @test JuMP.fix_value(y1) == 1.0
+    DP.unfix_indicator(model, Y2)
+    @test !JuMP.is_fixed(y1)
+end
+
+function test_avoid_combination_default_map()
+    # The 2-arg form resolves binaries via the model's own
+    # indicator-to-binary map. Excluding (Y1 active, Y2 inactive) adds
+    # (1 - y1) + y2 >= 1, i.e. normalized -y1 + y2 >= 0.
+    model = GDPModel()
+    @variable(model, x)
+    @variable(model, Y[1:2], Logical)
+    @constraint(model, x <= 3, Disjunct(Y[1]))
+    @constraint(model, x >= 5, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    DP.reformulate_model(model, BigM())
+    binary_map = DP._indicator_to_binary(model)
+
+    cref = DP.avoid_combination(model, Dict(Y[1] => true, Y[2] => false))
+    @test JuMP.normalized_coefficient(cref, binary_map[Y[1]]) == -1.0
+    @test JuMP.normalized_coefficient(cref, binary_map[Y[2]]) == 1.0
+    @test JuMP.normalized_rhs(cref) == 0.0
+end
+
 @testset "Utility Functions" begin
     test_all_variables()
     test_collect_all_vars()
@@ -245,4 +365,10 @@ end
     test_get_constant_quadratic()
     test_get_constant_number()
     test_get_constant_variable()
+    test_linearize_nonlinear_exp()
+    test_linearize_nonlinear_sin()
+    test_linearize_nonlinear_multivar()
+    test_to_nlp_expr()
+    test_fix_unfix_indicator()
+    test_avoid_combination_default_map()
 end