diff --git a/Project.toml b/Project.toml
index a5320aa..f3fbb43 100644
--- a/Project.toml
+++ b/Project.toml
@@ -20,7 +20,7 @@ Reexport = "1"
 julia = "1.10"
 Juniper = "0.9.3"
 Ipopt = "1.9.0"
-InfiniteOpt = "0.6"
+InfiniteOpt = "0.6.3"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
diff --git a/ext/InfiniteDisjunctiveProgramming.jl b/ext/InfiniteDisjunctiveProgramming.jl
index 426f7ac..5f77dce 100644
--- a/ext/InfiniteDisjunctiveProgramming.jl
+++ b/ext/InfiniteDisjunctiveProgramming.jl
@@ -9,8 +9,8 @@ import DisjunctiveProgramming as DP
 ################################################################################
 function DP.InfiniteGDPModel(args...; kwargs...)
     return DP.GDPModel{
-        InfiniteOpt.InfiniteModel, 
-        InfiniteOpt.GeneralVariableRef, 
+        InfiniteOpt.InfiniteModel,
+        InfiniteOpt.GeneralVariableRef,
         InfiniteOpt.InfOptConstraintRef
         }(args...; kwargs...)
 end
@@ -26,7 +26,7 @@ end
 ################################################################################
 DP.InfiniteLogical(prefs...) = DP.Logical(InfiniteOpt.Infinite(prefs...))
 
-_is_parameter(vref::InfiniteOpt.GeneralVariableRef) = 
+_is_parameter(vref::InfiniteOpt.GeneralVariableRef) =
     _is_parameter(InfiniteOpt.dispatch_variable_ref(vref))
 _is_parameter(::InfiniteOpt.DependentParameterRef) = true
 _is_parameter(::InfiniteOpt.IndependentParameterRef) = true
@@ -171,17 +171,207 @@ function DP.disaggregate_expression(
 end
 
 ################################################################################
-#                              ERROR MESSAGES
+#                          MBM FOR INFINITEMODEL
+################################################################################
+
+# Copy the InfiniteModel, strip everything but VariableInfo bounds,
+# add back the selected disjunct constraints, transcribe, and return
+# only the other disjunct's constraints plus variable bounds.
+function DP.copy_model_with_constraints(
+    model::InfiniteOpt.InfiniteModel,
+    constraints::Vector{<:DP.DisjunctConstraintRef},
+    method::DP._MBM
+    )
+    # Filter out every source constraint at copy time instead of
+    # copying then deleting. Equivalent end state, fewer allocations.
+    mini, ref_map = JuMP.copy_model(
+        model; filter_constraints = cref -> false
+        )
+
+    for cref in constraints
+        con = JuMP.constraint_object(cref)
+        T = one(JuMP.value_type(typeof(mini)))
+        JuMP.@constraint(mini, ref_map[con.func] * T in con.set)
+    end
+
+    InfiniteOpt.build_transformation_backend!(mini)
+    transcribed = InfiniteOpt.transformation_model(mini)
+    JuMP.set_optimizer(transcribed, method.optimizer)
+    JuMP.set_silent(transcribed)
+
+    # fwd_map needs every ref reachable from disjunct constraints —
+    # decision vars + parameters + parameter functions so the
+    # objective substitution in `prepare_max_M_objective` can look up
+    # any term it sees.
+    decision_vars = DP.collect_all_vars(model)
+    fwd_map = Dict{InfiniteOpt.GeneralVariableRef,
+        Vector{InfiniteOpt.GeneralVariableRef}}()
+    for v in decision_vars
+        fwd_map[v] = [ref_map[v]]
+    end
+    for p in InfiniteOpt.all_parameters(model)
+        fwd_map[p] = [ref_map[p]]
+    end
+    for pf in InfiniteOpt.all_parameter_functions(model)
+        fwd_map[pf] = [ref_map[pf]]
+    end
+    return DP.GDPSubmodel(mini, decision_vars, fwd_map)
+end
+
+function DP.prepare_max_M_objective(
+    ::InfiniteOpt.InfiniteModel,
+    obj::JuMP.ScalarConstraint{T, S},
+    sub::DP.GDPSubmodel
+    ) where {T, S <: _MOI.LessThan}
+    flat_map = Dict(v => ws[1] for (v, ws) in sub.fwd_map)
+    obj_func = DP.replace_variables_in_constraint(obj.func, flat_map)
+    return obj_func - obj.set.upper
+end
+
+function DP.prepare_max_M_objective(
+    ::InfiniteOpt.InfiniteModel,
+    obj::JuMP.ScalarConstraint{T, S},
+    sub::DP.GDPSubmodel
+    ) where {T, S <: _MOI.GreaterThan}
+    flat_map = Dict(v => ws[1] for (v, ws) in sub.fwd_map)
+    obj_func = DP.replace_variables_in_constraint(obj.func, flat_map)
+    return obj.set.lower - obj_func
+end
+
+# Constant interpolation
+function _interpolate(
+    grids::NTuple{N, AbstractVector{<:Real}},
+    values::AbstractArray{<:Real, N}
+    ) where {N}
+    # mimic the call form of Interpolations.jl's interpolation
+    return (args...) -> _interpolate_at(grids, values, args)
+end
+
+function _interpolate_at(
+    grids::NTuple{N, AbstractVector{<:Real}},
+    values::AbstractArray{<:Real, N},
+    args::NTuple{N, <:Real}
+    ) where {N}
+    # lower-corner cell index per dimension
+    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
+    return maximum(
+        values[ntuple(d -> idx_lo[d] +((k >> (d - 1)) & 1), N)...]
+        for k in 0:(2^N - 1)
+        )
+end
+
+# Transcribe mini_expr, solve per support on the transcribed JuMP
+# model, and aggregate to a scalar if uniform, else to a parameter
+# function on main.
+function DP.raw_M(
+    sub::DP.GDPSubmodel{<:InfiniteOpt.InfiniteModel},
+    mini_expr::JuMP.AbstractJuMPScalar,
+    method::DP._MBM
+    )
+    objectives = InfiniteOpt.transformation_expression(mini_expr)
+    transcribed = InfiniteOpt.transformation_model(sub.model)
+    inner_sub = DP.GDPSubmodel(transcribed,JuMP.VariableRef[],
+        Dict{JuMP.VariableRef, Vector{JuMP.VariableRef}}()
+        )
+    M_vals = Array{typeof(method.default_M)}(undef, size(objectives))
+    for I in eachindex(objectives)
+        m = DP.raw_M(inner_sub, objectives[I], method)
+        m === nothing && return nothing
+        M_vals[I] = m
+    end
+    all(==(first(M_vals)), M_vals) && return first(M_vals)
+    mini_prefs = InfiniteOpt.parameter_refs(mini_expr)
+    reverse_map = Dict(ws[1] => v for (v, ws) in sub.fwd_map)
+    prefs = Tuple(reverse_map[p] for p in mini_prefs)
+    main = JuMP.owner_model(first(prefs))
+    grids = Tuple(InfiniteOpt.supports(p) for p in prefs)
+    param_func = InfiniteOpt.build_parameter_function(
+        error, _interpolate(grids, M_vals), prefs)
+    return InfiniteOpt.add_parameter_function(main, param_func)
+end
+
+################################################################################
+#                    CUTTING PLANES FOR INFINITEMODEL
 ################################################################################
-function DP.reformulate_model(::InfiniteOpt.InfiniteModel, ::DP.MBM)
-    error("The `MBM` method is not supported for `InfiniteModel`." *
-          "Please use `BigM`, `Hull`, `Indicator`, or `PSplit` instead.")
+
+# Build CP subproblem: reformulate the InfiniteModel in-place, transcribe,
+# copy, and wrap in GDPSubmodel with forward variable map.
+function DP.copy_and_reformulate(
+    model::InfiniteOpt.InfiniteModel,
+    decision_vars::Vector{InfiniteOpt.GeneralVariableRef},
+    reform_method::DP.AbstractReformulationMethod,
+    method::DP.CuttingPlanes
+    )
+    DP.reformulate_model(model, reform_method)
+    InfiniteOpt.build_transformation_backend!(model)
+    transcribed = InfiniteOpt.transformation_model(model)
+    transcription_fwd = Dict{InfiniteOpt.GeneralVariableRef,
+        Vector{JuMP.VariableRef}}()
+    for v in DP.collect_all_vars(model)
+        transcription_var = InfiniteOpt.transformation_variable(v)
+        var_prefs = InfiniteOpt.parameter_refs(v)
+        transcription_fwd[v] = isempty(var_prefs) ?
+            [transcription_var] : vec(transcription_var)
+    end
+    sub_copy, copy_map = JuMP.copy_model(transcribed)
+    fwd_map = Dict{InfiniteOpt.GeneralVariableRef, Vector{JuMP.VariableRef}}()
+    for v in decision_vars
+        haskey(transcription_fwd, v) || continue
+        fwd_map[v] = [copy_map[transcribed_var] for transcribed_var in transcription_fwd[v]]
+    end
+    sub = DP.GDPSubmodel(sub_copy, decision_vars, fwd_map)
+    JuMP.set_optimizer(sub.model, method.optimizer)
+    JuMP.set_silent(sub.model)
+    return sub
 end
 
-function DP.reformulate_model(::InfiniteOpt.InfiniteModel, ::DP.CuttingPlanes)
-    error("The `CuttingPlanes` method is not supported for `InfiniteModel`." *
-          "Please use `BigM`, `Hull`, `Indicator`, or `PSplit` instead.")
+# Read per-support values from the transformation backend.
+function DP.extract_solution(model::InfiniteOpt.InfiniteModel)
+    dvars = DP.collect_cutting_planes_vars(model)
+    V = eltype(dvars)
+    T = JuMP.value_type(typeof(model))
+    sol = Dict{V, Vector{T}}()
+    for v in dvars
+        transcription_var = InfiniteOpt.transformation_variable(v)
+        var_prefs = InfiniteOpt.parameter_refs(v)
+        sol[v] = isempty(var_prefs) ? [JuMP.value(transcription_var)] :
+            JuMP.value.(vec(transcription_var))
+    end
+    return sol
 end
 
+# Add a pointwise-sum cut directly to the transformation backend and mark
+# it ready so the next optimize! doesn't re-transcribe and wipe the cut.
+function DP.add_cut(
+    model::InfiniteOpt.InfiniteModel,
+    decision_vars::Vector{InfiniteOpt.GeneralVariableRef},
+    rBM_sol::Dict{<:JuMP.AbstractVariableRef, <:Vector{<:Number}},
+    sep_sol::Dict{<:JuMP.AbstractVariableRef, <:Vector{<:Number}}
+    )
+    transcribed = InfiniteOpt.transformation_model(model)
+    cut_expr = zero(JuMP.GenericAffExpr{
+        JuMP.value_type(typeof(transcribed)),
+        JuMP.variable_ref_type(transcribed)})
+    for var in decision_vars
+        haskey(rBM_sol, var) || continue
+        haskey(sep_sol, var) || continue
+        rbm_vals = rBM_sol[var]
+        sep_vals = sep_sol[var]
+        transcription_var = InfiniteOpt.transformation_variable(var)
+        transcribed_vars = transcription_var isa AbstractArray ?
+            vec(transcription_var) : [transcription_var]
+        for k in eachindex(transcribed_vars)
+            xi = 2 * (sep_vals[k] - rbm_vals[k])
+            JuMP.add_to_expression!(cut_expr, xi, transcribed_vars[k])
+            JuMP.add_to_expression!(cut_expr, -xi * sep_vals[k])
+        end
+    end
+    JuMP.@constraint(transcribed, cut_expr >= 0)
+    InfiniteOpt.set_transformation_backend_ready(model, true)
+    return
 end
 
+end
diff --git a/src/mbm.jl b/src/mbm.jl
index 973c051..23c5096 100644
--- a/src/mbm.jl
+++ b/src/mbm.jl
@@ -242,7 +242,7 @@ end
     prepare_max_M_objective(model, obj::ScalarConstraint, sub::GDPSubmodel)
 
 Convert a constraint into an objective expression for M-value
-maximization. Returns a single JuMP expression to pass to `_raw_M`.
+maximization. Returns a single JuMP expression to pass to `raw_M`.
 """
 function prepare_max_M_objective(
     ::JuMP.AbstractModel,
@@ -250,7 +250,7 @@ function prepare_max_M_objective(
     sub::GDPSubmodel
     ) where {T, S <: _MOI.LessThan}
     flat_map = Dict(v => ws[1] for (v, ws) in sub.fwd_map)
-    expr = -obj.set.upper + _replace_variables_in_constraint(obj.func, flat_map)
+    expr = -obj.set.upper + replace_variables_in_constraint(obj.func, flat_map)
     return expr
 end
 
@@ -260,14 +260,20 @@ function prepare_max_M_objective(
     sub::GDPSubmodel
     ) where {T, S <: _MOI.GreaterThan}
     flat_map = Dict(v => ws[1] for (v, ws) in sub.fwd_map)
-    expr = obj.set.lower - _replace_variables_in_constraint(obj.func, flat_map)
+    expr = obj.set.lower - replace_variables_in_constraint(obj.func, flat_map)
     return expr
 end
 
-# Solve the submodel for a single objective expression.
-# Returns a scalar M value, or nothing if infeasible.
-function _raw_M(
-    sub::GDPSubmodel,
+"""
+    raw_M(sub::GDPSubmodel, objective, method::_MBM)
+
+Maximize `objective` over `sub` to obtain one raw M value for MBM.
+Returns `max(obj_value, 0)` on optimal, `nothing` on infeasible
+(signals the constraint is redundant in the combined region), or
+`method.default_M` otherwise (unbounded, numerical failure, etc).
+"""
+function raw_M(
+    sub::GDPSubmodel{<:JuMP.AbstractModel},
     objective::JuMP.AbstractJuMPScalar,
     method::_MBM
     )
@@ -291,7 +297,7 @@ function _maximize_M(
     method::_MBM
     ) where {T, S <: Union{_MOI.LessThan, _MOI.GreaterThan}}
     sub = _get_submodel(model, constraints, method)
-    return _raw_M(sub,
+    return raw_M(sub,
         prepare_max_M_objective(model, objective, sub), method)
 end
 
@@ -321,8 +327,8 @@ function _maximize_M(
     set_value = objective.set.value
     ge_obj = JuMP.ScalarConstraint(objective.func, MOI.GreaterThan(set_value))
     le_obj = JuMP.ScalarConstraint(objective.func, MOI.LessThan(set_value))
-    raw_lower = _raw_M(sub, prepare_max_M_objective(model, ge_obj, sub), method)
-    raw_upper = _raw_M(sub, prepare_max_M_objective(model, le_obj, sub), method)
+    raw_lower = raw_M(sub, prepare_max_M_objective(model, ge_obj, sub), method)
+    raw_upper = raw_M(sub, prepare_max_M_objective(model, le_obj, sub), method)
     (raw_lower === nothing || raw_upper === nothing) &&
         return nothing
     return [raw_lower, raw_upper]
@@ -341,8 +347,8 @@ function _maximize_M(
         MOI.GreaterThan(set_values[1]))
     le_obj = JuMP.ScalarConstraint(objective.func,
         MOI.LessThan(set_values[2]))
-    raw_lower = _raw_M(sub, prepare_max_M_objective(model, ge_obj, sub), method)
-    raw_upper = _raw_M(sub, prepare_max_M_objective(model, le_obj, sub), method)
+    raw_lower = raw_M(sub, prepare_max_M_objective(model, ge_obj, sub), method)
+    raw_upper = raw_M(sub, prepare_max_M_objective(model, le_obj, sub), method)
     (raw_lower === nothing || raw_upper === nothing) &&
         return nothing
     return [raw_lower, raw_upper]
@@ -361,7 +367,7 @@ function _maximize_M(
     for i in 1:objective.set.dimension
         le_obj = JuMP.ScalarConstraint(
             objective.func[i], MOI.LessThan(zero(val_type)))
-        raw = _raw_M(sub,
+        raw = raw_M(sub,
             prepare_max_M_objective(model, le_obj, sub),
             method)
         raw === nothing && return nothing
@@ -384,7 +390,7 @@ function _maximize_M(
         ge_obj = JuMP.ScalarConstraint(
             objective.func[i],
             MOI.GreaterThan(zero(val_type)))
-        raw = _raw_M(sub,
+        raw = raw_M(sub,
             prepare_max_M_objective(model, ge_obj, sub),
             method)
         raw === nothing && return nothing
@@ -410,10 +416,10 @@ function _maximize_M(
         le_obj = JuMP.ScalarConstraint(
             objective.func[i],
             MOI.LessThan(zero(val_type)))
-        raw_ge = _raw_M(sub,
+        raw_ge = raw_M(sub,
             prepare_max_M_objective(model, ge_obj, sub),
             method)
-        raw_le = _raw_M(sub,
+        raw_le = raw_M(sub,
             prepare_max_M_objective(model, le_obj, sub),
             method)
         (raw_ge === nothing || raw_le === nothing) &&
@@ -456,7 +462,7 @@ function copy_model_with_constraints(
     for cref in constraints
         con = JuMP.constraint_object(cref)
         flat_map = Dict(v => only(ws) for (v, ws) in fwd_map)
-        expr = _replace_variables_in_constraint(con.func, flat_map)
+        expr = replace_variables_in_constraint(con.func, flat_map)
         T = one(JuMP.value_type(typeof(sub_model)))
         JuMP.@constraint(sub_model, expr * T in con.set)
     end
@@ -474,7 +480,7 @@ end
 # Replace variable refs in an expression using a map. Uses AbstractDict
 # because the InfiniteModel MBM path maps decision vars to VariableRefs
 # and parameter functions to Numbers in the same dict (via _build_flat_map).
-function _replace_variables_in_constraint(
+function replace_variables_in_constraint(
     fun::JuMP.AbstractVariableRef,
     var_map::AbstractDict
     )
@@ -495,7 +501,7 @@ function _var_ref_type(
     return V
 end
 
-function _replace_variables_in_constraint(
+function replace_variables_in_constraint(
     fun::T, var_map::AbstractDict
     ) where {T <: JuMP.GenericAffExpr}
     C = JuMP.value_type(T)
@@ -508,7 +514,7 @@ function _replace_variables_in_constraint(
     return new_aff
 end
 
-function _replace_variables_in_constraint(
+function replace_variables_in_constraint(
     fun::T, var_map::AbstractDict
     ) where {T <: JuMP.GenericQuadExpr}
     C = JuMP.value_type(T)
@@ -518,23 +524,23 @@ function _replace_variables_in_constraint(
         JuMP.add_to_expression!(new_quad,
             coef * var_map[vars.a] * var_map[vars.b])
     end
-    new_aff = _replace_variables_in_constraint(fun.aff, var_map)
+    new_aff = replace_variables_in_constraint(fun.aff, var_map)
     JuMP.add_to_expression!(new_quad, new_aff)
     return new_quad
 end
 
-function _replace_variables_in_constraint(fun::Number, var_map::AbstractDict)
+function replace_variables_in_constraint(fun::Number, var_map::AbstractDict)
     return fun
 end
 
-function _replace_variables_in_constraint(fun::T,
+function replace_variables_in_constraint(fun::T,
     var_map::AbstractDict) where {T <: JuMP.GenericNonlinearExpr}
-    new_args = Any[_replace_variables_in_constraint(
+    new_args = Any[replace_variables_in_constraint(
         arg, var_map) for arg in fun.args]
     return T(fun.head, new_args)
 end
 
-function _replace_variables_in_constraint(fun::Vector, var_map::AbstractDict)
-    return [_replace_variables_in_constraint(expr,
+function replace_variables_in_constraint(fun::Vector, var_map::AbstractDict)
+    return [replace_variables_in_constraint(expr,
         var_map) for expr in fun]
 end
diff --git a/src/utilities.jl b/src/utilities.jl
index b0203d7..bca1b74 100644
--- a/src/utilities.jl
+++ b/src/utilities.jl
@@ -30,7 +30,7 @@ function copy_and_reformulate(
     orig_to_copy = Dict{V, V}(
         v => ref_map[v] for v in decision_vars)
     JuMP.@objective(copy, sense,
-        _replace_variables_in_constraint(obj, orig_to_copy)
+        replace_variables_in_constraint(obj, orig_to_copy)
         )
     fwd_map = Dict{V, Vector{V}}(v => [ref_map[v]] for v in decision_vars)
     sub = GDPSubmodel(copy, decision_vars, fwd_map)
@@ -221,7 +221,7 @@ function copy_gdp_data(
         old_con_ref = LogicalConstraintRef(model, idx)
         new_con_ref = LogicalConstraintRef(new_model, idx)
         c = lc_data.constraint
-        expr = _replace_variables_in_constraint(c.func, lv_map)
+        expr = replace_variables_in_constraint(c.func, lv_map)
         new_con = JuMP.build_constraint(error, expr, c.set)
         JuMP.add_constraint(new_model, new_con, lc_data.name)
         lc_map[old_con_ref] = new_con_ref
@@ -233,7 +233,7 @@ function copy_gdp_data(
         old_dc_ref = DisjunctConstraintRef(model, idx)
         old_indicator = old_gdp.constraint_to_indicator[old_dc_ref]
         new_indicator = lv_map[old_indicator]
-        new_expr = _replace_variables_in_constraint(old_constraint.func, 
+        new_expr = replace_variables_in_constraint(old_constraint.func, 
         var_map
         )
         # Update to new_gdp.disjunct_constraints
@@ -247,7 +247,7 @@ function copy_gdp_data(
     # Copying disjunctions
     for (idx, disj_data) in old_gdp.disjunctions
         old_disj = disj_data.constraint
-        new_indicators = [_replace_variables_in_constraint(indicator, lv_map) 
+        new_indicators = [replace_variables_in_constraint(indicator, lv_map) 
         for indicator in old_disj.indicators
             ]
         new_disj = Disjunction(new_indicators, old_disj.nested)
@@ -383,7 +383,7 @@ function _remap_indicator_to_binary(
     bref::JuMP.GenericAffExpr,
     var_map::Dict{V, V}
 ) where {V <: JuMP.AbstractVariableRef}
-    return _replace_variables_in_constraint(bref, var_map)
+    return replace_variables_in_constraint(bref, var_map)
 end
 
 function _remap_constraint_to_indicator(
diff --git a/src/variables.jl b/src/variables.jl
index f96cfe0..e3369db 100644
--- a/src/variables.jl
+++ b/src/variables.jl
@@ -547,7 +547,9 @@ function get_variable_info(vref::JuMP.AbstractVariableRef;
     has_lb::Bool = JuMP.has_lower_bound(vref), 
     has_ub::Bool = JuMP.has_upper_bound(vref), 
     has_fix::Bool = JuMP.is_fixed(vref), 
-    has_start::Bool = JuMP.has_start_value(vref), 
+    has_start::Bool = JuMP.has_start_value(vref) &&
+        !(JuMP.start_value(vref) isa Number &&
+          isnan(JuMP.start_value(vref))),
     has_binary::Bool = JuMP.is_binary(vref), 
     has_integer::Bool = JuMP.is_integer(vref), 
     lower_bound::Union{Number, Function} = has_lb ? JuMP.lower_bound(vref) : 0, 
diff --git a/test/constraints/mbm.jl b/test/constraints/mbm.jl
index d181c08..11f164d 100644
--- a/test/constraints/mbm.jl
+++ b/test/constraints/mbm.jl
@@ -27,7 +27,7 @@ function test__var_ref_type_numeric_map()
     @test DP._var_ref_type(typeof(aff), var_map) == VariableRef
 end
 
-# _replace_variables_in_constraint with QuadExpr where var_map
+# replace_variables_in_constraint with QuadExpr where var_map
 # maps some vars to Numbers. Covers lines 569, 571, 574.
 function test__replace_variables_quad_numeric_map()
     model = Model()
@@ -38,17 +38,17 @@ function test__replace_variables_quad_numeric_map()
 
     # both Number (line 569)
     map1 = Dict{VariableRef, Any}(x[1] => 2.0, x[2] => 3.0)
-    result1 = DP._replace_variables_in_constraint(quad1, map1)
+    result1 = DP.replace_variables_in_constraint(quad1, map1)
     @test result1.aff.constant ≈ 6.0
 
     # ra Number, rb VariableRef (line 571)
     map2 = Dict{VariableRef, Any}(x[1] => 2.0, x[2] => y)
-    result2 = DP._replace_variables_in_constraint(quad1, map2)
+    result2 = DP.replace_variables_in_constraint(quad1, map2)
     @test result2.aff.terms[y] ≈ 2.0
 
     # rb Number, ra VariableRef (line 574)
     map3 = Dict{VariableRef, Any}(x[1] => y, x[2] => 3.0)
-    result3 = DP._replace_variables_in_constraint(quad1, map3)
+    result3 = DP.replace_variables_in_constraint(quad1, map3)
     @test result3.aff.terms[y] ≈ 3.0
 end
 
@@ -64,14 +64,14 @@ function test_replace_variables_in_constraint()
     #Test GenericVariableRef
     new_vars = Dict{AbstractVariableRef, AbstractVariableRef}()
     [new_vars[x[i]] = @variable(sub_model) for i in 1:3]
-    varref = DP._replace_variables_in_constraint(x[1], new_vars)
-    expr1 = DP._replace_variables_in_constraint(
+    varref = DP.replace_variables_in_constraint(x[1], new_vars)
+    expr1 = DP.replace_variables_in_constraint(
         constraint_object(con1).func, new_vars)
-    expr2 = DP._replace_variables_in_constraint(
+    expr2 = DP.replace_variables_in_constraint(
         constraint_object(con2).func, new_vars)
-    expr3 = DP._replace_variables_in_constraint(
+    expr3 = DP.replace_variables_in_constraint(
         constraint_object(con3).func, new_vars)
-    expr4 = DP._replace_variables_in_constraint(
+    expr4 = DP.replace_variables_in_constraint(
         constraint_object(con4).func, new_vars)
     @test expr1 == JuMP.@expression(sub_model, new_vars[x[1]] + 1 - 1)
     @test expr2 == JuMP.@expression(sub_model, new_vars[x[2]]*new_vars[x[1]])
@@ -80,7 +80,7 @@ function test_replace_variables_in_constraint()
     expected = JuMP.@expression(sub_model, sin(new_vars[x[3]]) - 0.0)
     @test JuMP.isequal_canonical(expr3, expected)
     @test expr4 == [new_vars[x[i]] for i in 1:3]
-    @test_throws MethodError DP._replace_variables_in_constraint(
+    @test_throws MethodError DP.replace_variables_in_constraint(
         "String", new_vars)
 end
 
@@ -134,19 +134,19 @@ function test_raw_M()
         DisjunctConstraintRef[con2], mbm)
     obj = DP.prepare_max_M_objective(model,
         constraint_object(con), sub)
-    @test DP._raw_M(sub, obj, mbm) == 0.0
+    @test DP.raw_M(sub, obj, mbm) == 0.0
     set_upper_bound(x, 1)
     sub2 = DP.copy_model_with_constraints(model,
         DisjunctConstraintRef[con], mbm)
     obj2 = DP.prepare_max_M_objective(model,
         constraint_object(con2), sub2)
-    @test DP._raw_M(sub2, obj2, mbm) == 15
+    @test DP.raw_M(sub2, obj2, mbm) == 15
     set_integer(y)
     @constraint(model, con3, y*x == 15,
         Disjunct(Y[1]))
     obj3 = DP.prepare_max_M_objective(model,
         constraint_object(con2), sub2)
-    @test DP._raw_M(sub2, obj3, mbm) == 15
+    @test DP.raw_M(sub2, obj3, mbm) == 15
     # Fresh _MBM after changing bounds
     JuMP.fix(y, 5; force=true)
     mbm2 = DP._MBM(
@@ -155,7 +155,7 @@ function test_raw_M()
         DisjunctConstraintRef[con], mbm2)
     obj4 = DP.prepare_max_M_objective(model,
         constraint_object(con2), sub3)
-    @test DP._raw_M(sub3, obj4, mbm2) == 10
+    @test DP.raw_M(sub3, obj4, mbm2) == 10
     # Infeasible region → nothing
     delete_lower_bound(x)
     mbm3 = DP._MBM(
@@ -164,7 +164,7 @@ function test_raw_M()
         DisjunctConstraintRef[con2], mbm3)
     obj5 = DP.prepare_max_M_objective(model,
         constraint_object(con2), sub4)
-    @test DP._raw_M(sub4, obj5, mbm3) == nothing
+    @test DP.raw_M(sub4, obj5, mbm3) == nothing
 
     # infeasible (x >= 100 but x <= 1)
     set_upper_bound(x, 1)
@@ -175,7 +175,7 @@ function test_raw_M()
         mbm4)
     obj6 = DP.prepare_max_M_objective(model,
         constraint_object(con), sub5)
-    @test DP._raw_M(sub5, obj6, mbm4) == nothing
+    @test DP.raw_M(sub5, obj6, mbm4) == nothing
 
     # Unbounded subproblem → default_M fallback.
     # No lower bound on x means max(5 - x) s.t. x <= 3
@@ -191,7 +191,7 @@ function test_raw_M()
         DisjunctConstraintRef[ub_con1], mbm_ub)
     obj_ub = DP.prepare_max_M_objective(model_ub,
         constraint_object(ub_con2), sub_ub)
-    @test DP._raw_M(sub_ub, obj_ub, mbm_ub) == mbm_ub.default_M
+    @test DP.raw_M(sub_ub, obj_ub, mbm_ub) == mbm_ub.default_M
 end
 
 function test_maximize_M()
@@ -783,6 +783,20 @@ function test_get_variable_info()
     @test info_custom.has_ub == false
 end
 
+# A NaN-valued start is treated as no start so the NaN doesn't
+# propagate to copies or to solver inputs.
+function test_get_variable_info_nan_start()
+    model = GDPModel()
+    @variable(model, x)
+    JuMP.set_start_value(x, NaN)
+    @test JuMP.has_start_value(x) == true
+    @test isnan(JuMP.start_value(x))
+
+    info = DP.get_variable_info(x)
+    @test info.has_start == false
+    @test info.start == 0
+end
+
 @testset "MBM" begin
     test__copy_model()
     test_variable_properties()
@@ -791,6 +805,7 @@ end
     test_variable_copy()
     test__copy_model_with_constraints()
     test_get_variable_info()
+    test_get_variable_info_nan_start()
     test_mbm()
     test__var_ref_type_numeric_map()
     test__replace_variables_quad_numeric_map()
diff --git a/test/extensions/InfiniteDisjunctiveProgramming.jl b/test/extensions/InfiniteDisjunctiveProgramming.jl
index 8ba1e5c..47bb884 100644
--- a/test/extensions/InfiniteDisjunctiveProgramming.jl
+++ b/test/extensions/InfiniteDisjunctiveProgramming.jl
@@ -55,7 +55,7 @@ function test_infinite_logical()
     @test binary_variable(y) isa InfiniteOpt.GeneralVariableRef
 end
 
-function test__is_parameter()
+function test_is_parameter()
     model = InfiniteGDPModel()
     @infinite_parameter(model, t ∈ [0, 1])
     @infinite_parameter(model, s[1:2] ∈ [0, 1], independent = true)
@@ -77,6 +77,34 @@ function test__is_parameter()
     @test IDP._is_parameter(y) == false
 end
 
+# _is_parameter on unwrapped concrete dispatch types
+# (DependentParameterRef, IndependentParameterRef, FiniteParameterRef,
+# ParameterFunctionRef, Any fallback).
+function test_is_parameter_concrete_dispatches()
+    model = InfiniteGDPModel()
+    # Scalar + `independent = true` array both give IndependentParameterRef;
+    # a default array parameter gives DependentParameterRef.
+    @infinite_parameter(model, t ∈ [0, 1])
+    @infinite_parameter(model, s[1:2] ∈ [0, 1], independent = true)
+    @infinite_parameter(model, q[1:2] ∈ [0, 1])
+    @finite_parameter(model, p == 1.0)
+    @variable(model, x, Infinite(t))
+    @parameter_function(model, pf == t -> 2*t)
+    dvr = InfiniteOpt.dispatch_variable_ref
+    # Verify each ref hits the intended dispatch.
+    @test dvr(t) isa InfiniteOpt.IndependentParameterRef
+    @test dvr(s[1]) isa InfiniteOpt.IndependentParameterRef
+    @test dvr(q[1]) isa InfiniteOpt.DependentParameterRef
+    @test dvr(p) isa InfiniteOpt.FiniteParameterRef
+    @test dvr(pf) isa InfiniteOpt.ParameterFunctionRef
+    @test IDP._is_parameter(dvr(t)) == true
+    @test IDP._is_parameter(dvr(s[1])) == true
+    @test IDP._is_parameter(dvr(q[1])) == true
+    @test IDP._is_parameter(dvr(p)) == true
+    @test IDP._is_parameter(dvr(pf)) == true
+    @test IDP._is_parameter(dvr(x)) == false    # Any fallback
+end
+
 function test_requires_disaggregation()
     model = InfiniteGDPModel()
     @infinite_parameter(model, t ∈ [0, 1])
@@ -336,10 +364,232 @@ function test_logical_value()
     @test eltype(val) == Bool
 end
 
-function test_unsupported_methods_error()
+# raw_M against an InfiniteModel where M is constant across supports.
+# Setup: x(t) ∈ [0, 10], disj1: x ≥ 5, disj2: x ≤ 3.
+# For disj1 slack r(x) = 5 - x maximized over disj2's region x ∈ [0, 3]:
+# max(5 - x) = 5 at x = 0. Same at every support ⇒ scalar M = 5.
+function test_raw_M_infinite_scalar()
+    model = InfiniteGDPModel()
+    @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, con, x >= 5, Disjunct(Y[1]))
+    @constraint(model, con2, x <= 3, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    mbm = DP._MBM(MBM(HiGHS.Optimizer), model)
+    sub = DP.copy_model_with_constraints(
+        model, DP.DisjunctConstraintRef[con2], mbm)
+    obj = DP.prepare_max_M_objective(
+        model, JuMP.constraint_object(con), sub)
+    @test length(InfiniteOpt.parameter_refs(obj)) == 1
+    @test DP.raw_M(sub, obj, mbm) == 5.0
+end
+
+# raw_M with a support-varying M. Setup: x(t) ∈ [0, 10], disj1: x ≤ 2t,
+# disj2: x ≥ 0.5. Slack r(x) = x - 2t maximized over x ∈ [0.5, 10]:
+# max(x - 2t) = 10 - 2t. Varies with t ⇒ raw_M returns a pfunc whose
+# raw values at supports are max-of-cell upper bounds for 10 - 2t.
+function test_raw_M_infinite_param_function()
+    model = InfiniteGDPModel()
+    supports = [0.0, 0.25, 0.5, 0.75, 1.0]
+    @infinite_parameter(model, t ∈ [0, 1], supports = supports)
+    @variable(model, 0 <= x <= 10, Infinite(t))
+    @variable(model, Y[1:2], InfiniteLogical(t))
+    @parameter_function(model, f == t -> 2*t)
+    @constraint(model, con, x <= f, Disjunct(Y[1]))
+    @constraint(model, con2, x >= 0.5, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    mbm = DP._MBM(MBM(HiGHS.Optimizer), model)
+    sub = DP.copy_model_with_constraints(
+        model, DP.DisjunctConstraintRef[con2], mbm)
+    obj = DP.prepare_max_M_objective(
+        model, JuMP.constraint_object(con), sub)
+    M = DP.raw_M(sub, obj, mbm)
+    @test M isa InfiniteOpt.GeneralVariableRef
+    raw_fn = InfiniteOpt.raw_function(M)
+    # max-of-corners is conservative: raw_fn(t) ≥ 10 - 2t at supports.
+    for t_val in supports
+        @test raw_fn(t_val) >= 10.0 - 2*t_val - 1e-6
+    end
+end
+
+# Piecewise-constant max-of-corners: returns the maximum value over
+# the 2^n corners of the cell containing the query.
+function test_interpolate()
+    grid1 = [0.0, 1.0, 2.0, 3.0]
+    vals1 = [10.0, 20.0, 40.0, 50.0]
+    f = IDP._interpolate((grid1,), vals1)
+    # At grid points: max over the cell to the right (or last cell).
+    @test f(0.0) == 20.0   # max(vals[1], vals[2])
+    @test f(1.0) == 40.0   # max(vals[2], vals[3])
+    @test f(3.0) == 50.0   # last cell: max(vals[3], vals[4])
+    # Between grid points: max of the surrounding two values.
+    @test f(0.5) == 20.0
+    @test f(1.5) == 40.0
+    @test f(2.25) == 50.0
+    # Out-of-range clamps to the boundary cell.
+    @test f(-1.0) == 20.0
+    @test f(4.0) == 50.0
+
+    # 2D: max over the 4 surrounding corners.
+    gx = [0.0, 1.0, 2.0]
+    gy = [0.0, 10.0]
+    vals2 = [x * y for x in gx, y in gy]   # 3x2 matrix
+    g = IDP._interpolate((gx, gy), vals2)
+    @test g(0.0, 0.0) == 10.0   # corners (0,0)=0, (1,0)=0, (0,10)=0, (1,10)=10
+    @test g(2.0, 10.0) == 20.0  # last cell, max corner is (2,10)=20
+    @test g(0.5, 5.0) == 10.0   # corners 0,0,0,10 -> 10
+    @test g(1.5, 5.0) == 20.0   # corners 0,0,10,20 -> 20
+end
+
+# extract_solution returns per-support values from the transformation
+# backend. Setup: force disj 2 active (x ≤ 3), BigM-reformulate, solve
+# min ∫x ⇒ x = 0 at every support.
+function test_extract_solution_infinite()
+    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)
-    @test_throws ErrorException DP.reformulate_model(model, MBM(HiGHS.Optimizer))
-    @test_throws ErrorException DP.reformulate_model(model, CuttingPlanes(HiGHS.Optimizer))
+    set_silent(model)
+
+    @infinite_parameter(model, t ∈ [0, 1], num_supports = 10)
+    @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)
+
+    @objective(model, Min, ∫(x, t))
+
+    @test optimize!(model, gdp_method = MBM(HiGHS.Optimizer)) isa Nothing
+    @test termination_status(model) in
+        [MOI.OPTIMAL, MOI.LOCALLY_SOLVED]
+    # x=0 with disjunct 2 active (x <= 3) gives min
+    @test objective_value(model) ≈ 0.0 atol = 0.1
+end
+
+function test_mbm_infinite_param_dependent()
+    model = InfiniteGDPModel(HiGHS.Optimizer)
+    set_silent(model)
+
+    @infinite_parameter(model, t ∈ [0, 1], num_supports = 20)
+    @variable(model, -10 <= x <= 10, Infinite(t))
+    @variable(model, Y[1:2], InfiniteLogical(t))
+
+    # Parameter-dependent constraints:
+    # Disjunct 1: x(t) <= 2*t
+    # Disjunct 2: x(t) >= 1 - t
+    @parameter_function(model, f1 == t -> 2*t)
+    @parameter_function(model, f2 == t -> 1 - t)
+    @constraint(model, x <= f1, Disjunct(Y[1]))
+    @constraint(model, x >= f2, Disjunct(Y[2]))
+    @disjunction(model, Y)
+
+    @objective(model, Min, ∫(x, t))
+
+    @test optimize!(model, gdp_method = MBM(HiGHS.Optimizer)) isa Nothing
+    @test termination_status(model) in
+        [MOI.OPTIMAL, MOI.LOCALLY_SOLVED]
+end
+
+function test_mbm_vs_bigm_infinite()
+    # Compare MBM and BigM: should give same
+    # feasible set and optimal value.
+    for method_pair in [
+        (BigM(100), MBM(HiGHS.Optimizer))
+    ]
+        model1 = InfiniteGDPModel(HiGHS.Optimizer)
+        set_silent(model1)
+        @infinite_parameter(model1, t ∈ [0, 1], num_supports = 10)
+        @variable(model1, 0 <= x1 <= 10, Infinite(t))
+        @variable(model1, Y1[1:2], InfiniteLogical(t))
+        @constraint(model1, x1 >= 5, Disjunct(Y1[1]))
+        @constraint(model1, x1 <= 3, Disjunct(Y1[2]))
+        @disjunction(model1, Y1)
+        @objective(model1, Min, ∫(x1, t))
+        optimize!(model1, gdp_method = method_pair[1])
+        obj1 = objective_value(model1)
+
+        model2 = InfiniteGDPModel(HiGHS.Optimizer)
+        set_silent(model2)
+        @infinite_parameter(model2, t2 ∈ [0, 1], num_supports = 10)
+        @variable(model2, 0 <= x2 <= 10, Infinite(t2))
+        @variable(model2, Y2[1:2], InfiniteLogical(t2))
+        @constraint(model2, x2 >= 5, Disjunct(Y2[1]))
+        @constraint(model2, x2 <= 3, Disjunct(Y2[2]))
+        @disjunction(model2, Y2)
+        @objective(model2, Min, ∫(x2, t2))
+        optimize!(model2, gdp_method = method_pair[2])
+        obj2 = objective_value(model2)
+
+        @test obj1 ≈ obj2 atol = 0.5
+    end
 end
 
 function test_methods()
@@ -393,6 +643,154 @@ function test_methods()
     @test value(z) ≈ expected_z atol=tol
 end
 
+function test_mbm_with_derivatives()
+    model = InfiniteGDPModel(HiGHS.Optimizer)
+    set_silent(model)
+
+    @infinite_parameter(model, t ∈ [0, 1], num_supports = 10)
+    @variable(model, -5 <= x <= 5, Infinite(t))
+    @variable(model, Y[1:2], InfiniteLogical(t))
+
+    @constraint(model, ∂(x, t) >= 1, Disjunct(Y[1]))
+    @constraint(model, ∂(x, t) <= -1, Disjunct(Y[2]))
+    @disjunction(model, Y)
+
+    set_upper_bound(∂(x, t), 10)
+    set_lower_bound(∂(x, t), -10)
+
+    @objective(model, Min, ∫(x^2, t))
+
+    juniper = JuMP.optimizer_with_attributes(
+        Juniper.Optimizer,
+        "nl_solver" => JuMP.optimizer_with_attributes(
+            Ipopt.Optimizer, "print_level" => 0),
+        "log_levels" => []
+    )
+    set_optimizer(model, juniper)
+    @test optimize!(model, gdp_method = MBM(juniper)) isa Nothing
+    @test termination_status(model) in
+        [MOI.OPTIMAL, MOI.LOCALLY_SOLVED,
+         MOI.ALMOST_LOCALLY_SOLVED]
+end
+
+function test_CuttingPlanes_infinite_simple()
+    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, Y[1:2], InfiniteLogical(t))
+
+    @constraint(model, x >= 5, Disjunct(Y[1]))
+    @constraint(model, x <= 3, Disjunct(Y[2]))
+    @disjunction(model, Y)
+
+    @objective(model, Min, ∫(x, t))
+
+    # Should not throw
+    @test optimize!(model,
+        gdp_method = CuttingPlanes(
+            HiGHS.Optimizer; max_iter = 5)
+    ) isa Nothing
+    @test termination_status(model) in
+        [MOI.OPTIMAL, MOI.LOCALLY_SOLVED]
+end
+
+function test_CuttingPlanes_infinite_two_disj()
+    model = InfiniteGDPModel(HiGHS.Optimizer)
+    set_silent(model)
+
+    @infinite_parameter(model, t ∈ [0, 1], num_supports = 10)
+    @variable(model, 0 <= x[1:2] <= 10, Infinite(t))
+    @variable(model, W1[1:2], InfiniteLogical(t))
+    @variable(model, W2[1:2], InfiniteLogical(t))
+
+    @constraint(model, x[1] >= 2, Disjunct(W1[1]))
+    @constraint(model, x[1] <= 1, Disjunct(W1[2]))
+    @disjunction(model, W1)
+
+    @constraint(model, x[2] >= 3, Disjunct(W2[1]))
+    @constraint(model, x[2] <= 2, Disjunct(W2[2]))
+    @disjunction(model, W2)
+
+    @objective(model, Min, ∫(x[1] + x[2], t))
+
+    # Compare cutting planes vs BigM
+    optimize!(model,
+        gdp_method = CuttingPlanes(
+            HiGHS.Optimizer; max_iter = 10)
+    )
+    cp_obj = objective_value(model)
+
+    model2 = InfiniteGDPModel(HiGHS.Optimizer)
+    set_silent(model2)
+    @infinite_parameter(model2, t2 ∈ [0, 1], num_supports = 10)
+    @variable(model2, 0 <= x2[1:2] <= 10, Infinite(t2))
+    @variable(model2, V1[1:2], InfiniteLogical(t2))
+    @variable(model2, V2[1:2], InfiniteLogical(t2))
+    @constraint(model2, x2[1] >= 2, Disjunct(V1[1]))
+    @constraint(model2, x2[1] <= 1, Disjunct(V1[2]))
+    @disjunction(model2, V1)
+    @constraint(model2, x2[2] >= 3, Disjunct(V2[1]))
+    @constraint(model2, x2[2] <= 2, Disjunct(V2[2]))
+    @disjunction(model2, V2)
+    @objective(model2, Min, ∫(x2[1] + x2[2], t2))
+    optimize!(model2, gdp_method = BigM())
+    bigm_obj = objective_value(model2)
+
+    @test cp_obj ≈ bigm_obj atol = 1.0
+end
+
+
+
+function test_CuttingPlanes_with_cuts()
+    # Maximization with single-constraint disjuncts where Hull
+    # is strictly tighter than BigM. BigM allows x+y up to
+    # variable bounds (20), Hull limits to max(5,8)=8 — this
+    # forces cuts to tighten the relaxation. Finite var w exercises
+    # the isempty(var_prefs) branch in add_cut.
+    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 <= y <= 10, Infinite(t))
+    @variable(model, 0 <= w <= 10)
+    @variable(model, Y[1:2], InfiniteLogical(t))
+    @constraint(model, x + y <= 5, Disjunct(Y[1]))
+    @constraint(model, x + y <= 8, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @objective(model, Max, ∫(x + y, t) + w)
+    cutting_planes = CuttingPlanes(HiGHS.Optimizer;
+        max_iter = 30, seperation_tolerance = 1e-6)
+    @test optimize!(model, gdp_method = cutting_planes) isa Nothing
+    @test termination_status(model) in
+        [MOI.OPTIMAL, MOI.LOCALLY_SOLVED]
+end
+
+function test_CuttingPlanes_multiparameter()
+    model = InfiniteGDPModel(HiGHS.Optimizer)
+    set_silent(model)
+
+    @infinite_parameter(model, t ∈ [0, 1], num_supports = 5)
+    @infinite_parameter(model, s ∈ [0, 2], num_supports = 4)
+    @variable(model, 0 <= x <= 10, Infinite(t, s))
+    @variable(model, Y[1:2], InfiniteLogical(t, s))
+
+    @constraint(model, x >= 5, Disjunct(Y[1]))
+    @constraint(model, x <= 3, Disjunct(Y[2]))
+    @disjunction(model, Y)
+
+    @objective(model, Min, ∫(∫(x, t), s))
+
+    # Should not throw
+    @test optimize!(model,
+        gdp_method = CuttingPlanes(
+            HiGHS.Optimizer; max_iter = 5)
+    ) isa Nothing
+    @test termination_status(model) in
+        [MOI.OPTIMAL, MOI.LOCALLY_SOLVED]
+end
+
 @testset "InfiniteDisjunctiveProgramming" begin
 
     @testset "Model" begin
@@ -405,7 +803,8 @@ end
 
     @testset "Variables" begin
         test_infinite_logical()
-        test__is_parameter()
+        test_is_parameter()
+        test_is_parameter_concrete_dispatches()
         test_requires_disaggregation()
         test_variable_properties_infiniteopt()
         test_variable_properties_from_expr()
@@ -429,7 +828,17 @@ end
     @testset "Methods" begin
         test_get_constant()
         test_disaggregate_expression_infiniteopt()
-        test_unsupported_methods_error()
+    end
+
+    @testset "MBM" begin
+        test_interpolate()
+        test_raw_M_infinite_scalar()
+        test_raw_M_infinite_param_function()
+        test_mbm_finite_and_integer_var()
+        test_mbm_infinite_simple()
+        test_mbm_infinite_param_dependent()
+        test_mbm_vs_bigm_infinite()
+        test_mbm_with_derivatives()
     end
 
     @testset "Integration" begin
@@ -437,4 +846,13 @@ end
         test_methods()
     end
 
-end
\ No newline at end of file
+    @testset "Cutting Planes" begin
+        test_extract_solution_infinite()
+        test_add_cut_infinite()
+        test_CuttingPlanes_infinite_simple()
+        test_CuttingPlanes_infinite_two_disj()
+        test_CuttingPlanes_with_cuts()
+        test_CuttingPlanes_multiparameter()
+    end
+
+end