Cutting-planes with quadrature weights for InfGDP

ext/InfiniteDisjunctiveProgramming.jl

@@ -303,38 +303,130 @@ function DP.copy_and_reformulate(
     model::InfiniteOpt.InfiniteModel,
     decision_vars::Vector{InfiniteOpt.GeneralVariableRef},
     reform_method::DP.AbstractReformulationMethod,
-    method::DP.CuttingPlanes
+    method::DP._CuttingPlanes
     )
+    # Quadrature weights come from the clean model, before reformulation.
+    for (v, w) in _compute_quadrature_weights(model, decision_vars)
+        method.weights[v] = w
+    end
     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]]
+        transcription_var = InfiniteOpt.transformation_variable(v)
+        tvars = isempty(InfiniteOpt.parameter_refs(v)) ?
+            [transcription_var] : vec(transcription_var)
+        fwd_map[v] = [copy_map[tv] for tv in tvars]
     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
 
+# Collect the objective's measure data (nested measures included),
+# mapped from each measured parameter to its measure data.
+_collect_measure_data(data, expr::Number) = nothing
+function _collect_measure_data(data::Dict, expr::InfiniteOpt.GeneralVariableRef)
+    dispatch = InfiniteOpt.dispatch_variable_ref(expr)
+    dispatch isa InfiniteOpt.MeasureRef || return nothing
+    md = InfiniteOpt.measure_data(expr)
+    prefs = InfiniteOpt.parameter_refs(md)
+    for p in (prefs isa AbstractArray ? prefs : (prefs,))
+        data[p] = md
+    end
+    return _collect_measure_data(data, InfiniteOpt.measure_function(expr))
+end
+function _collect_measure_data(data::Dict, expr::JuMP.GenericAffExpr)
+    for (v, _) in expr.terms
+        _collect_measure_data(data, v)
+    end
+    return nothing
+end
+function _collect_measure_data(data::Dict, expr::JuMP.GenericQuadExpr)
+    _collect_measure_data(data, expr.aff)
+    for (pair, _) in expr.terms
+        _collect_measure_data(data, pair.a)
+        _collect_measure_data(data, pair.b)
+    end
+    return nothing
+end
+function _collect_measure_data(data::Dict, expr::JuMP.GenericNonlinearExpr)
+    for arg in expr.args
+        _collect_measure_data(data, arg)
+    end
+    return nothing
+end
+
+# Whether a (scalar) parameter ranges over a plain interval;
+# dependent groups and distribution parameters return false.
+function _is_interval_parameter(pref)
+    pref isa AbstractArray && return false
+    dispatch = InfiniteOpt.dispatch_variable_ref(pref)
+    return InfiniteOpt.infinite_domain(dispatch) isa
+        InfiniteOpt.IntervalDomain
+end
+
+# Compute quadrature weights on a constraint-free copy of the model:
+# give it the objective sum_v m(v), with m the objective's own
+# measure data per parameter (if unmeasured, a default integral for
+# scalar interval parameters, else a support average - the measures
+# that keep the copy's supports, and so the weights, aligned with
+# the transcription), transcribe, and take each support's weight
+# from the transcribed coefficients. Finite variables get a unit
+# weight.
+function _compute_quadrature_weights(
+    model::InfiniteOpt.InfiniteModel,
+    decision_vars::Vector{InfiniteOpt.GeneralVariableRef}
+    )
+    mini, ref_map = JuMP.copy_model(model; filter_constraints = cref -> false)
+    measure_data = Dict{InfiniteOpt.GeneralVariableRef,
+        InfiniteOpt.AbstractMeasureData}()
+    _collect_measure_data(measure_data, JuMP.objective_function(mini))
+    terms = Any[]
+    for v in decision_vars
+        prefs = InfiniteOpt.parameter_refs(v)
+        isempty(prefs) && continue
+        expr = ref_map[v]
+        for g in prefs
+            mini_g = g isa AbstractArray ? [ref_map[p] for p in g] : ref_map[g]
+            key = mini_g isa AbstractArray ? first(mini_g) : mini_g
+            md = get(measure_data, key, nothing)
+            if md !== nothing
+                expr = InfiniteOpt.measure(expr, md)
+            elseif _is_interval_parameter(mini_g)
+                expr = InfiniteOpt.integral(expr, mini_g)
+            else
+                expr = InfiniteOpt.support_sum(expr, mini_g) /
+                    InfiniteOpt.num_supports(key)
+            end
+        end
+        push!(terms, expr)
+    end
+    JuMP.@objective(mini, Min, sum(terms))
+    InfiniteOpt.build_transformation_backend!(mini)
+    obj = JuMP.objective_function(InfiniteOpt.transformation_model(mini))
+    weights = Dict{InfiniteOpt.GeneralVariableRef, Vector{Float64}}()
+    for v in decision_vars
+        if isempty(InfiniteOpt.parameter_refs(v))
+            weights[v] = [1.0]
+        else
+            tvars = vec(InfiniteOpt.transformation_variable(ref_map[v]))
+            weights[v] = [JuMP.coefficient(obj, tv) for tv in tvars]
+        end
+    end
+    return weights
+end
+
 # 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
+function DP.extract_solution(
+    model::InfiniteOpt.InfiniteModel,
+    reform_state::DP._CuttingPlanes
+    )
+    sol = Dict{InfiniteOpt.GeneralVariableRef, Vector{Float64}}()
+    for v in reform_state.decision_vars
         transcription_var = InfiniteOpt.transformation_variable(v)
         var_prefs = InfiniteOpt.parameter_refs(v)
         sol[v] = isempty(var_prefs) ? [JuMP.value(transcription_var)] :
@@ -343,28 +435,26 @@ function DP.extract_solution(model::InfiniteOpt.InfiniteModel)
     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.
+# Add a quadrature-weighted 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},
+    reform_state::DP._CuttingPlanes,
     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
+    cut_expr = zero(JuMP.AffExpr)
+    for var in reform_state.decision_vars
         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]
+        w = reform_state.weights[var]
         for k in eachindex(transcribed_vars)
-            xi = 2 * (sep_vals[k] - rbm_vals[k])
+            xi = 2 * w[k] * (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

src/cuttingplanes.jl

@@ -8,14 +8,14 @@ function collect_cutting_planes_vars(model::JuMP.AbstractModel)
     return collect_all_vars(model)
 end
 
-# Extract solution from a solved model (in-place). Extensions
+# Extract solution from the solved rBM model (in-place). Extensions
 # override for models where values live on a backend.
-function extract_solution(model::JuMP.AbstractModel)
-    dvars = collect_cutting_planes_vars(model)
-    V = eltype(dvars)
-    T = JuMP.value_type(typeof(model))
+function extract_solution(
+    model::JuMP.AbstractModel,
+    reform_state::_CuttingPlanes{O, T, V}
+    ) where {O, T, V}
     return Dict{V, Vector{T}}(
-        v => [JuMP.value(v)] for v in dvars)
+        v => [JuMP.value(v)] for v in reform_state.decision_vars)
 end
 
 # Extract solution from a GDPSubmodel (SEP path).
@@ -29,9 +29,11 @@ function extract_solution(sub::GDPSubmodel)
     return sol
 end
 
-# Set quadratic separation objective: min Σ (x_k - rBM_k)².
+# Set the separation objective: min Σ_var Σ_k w_k (x_k - rBM_k)²,
+# with w the support weights carried by the CP working state.
 function _set_separation_objective(
     sub::GDPSubmodel,
+    weights::Dict,
     rBM_sol::Dict{<:JuMP.AbstractVariableRef, <:Vector{<:Number}}
     )
     obj_expr = zero(JuMP.GenericQuadExpr{
@@ -41,11 +43,10 @@ function _set_separation_objective(
     for var in sub.decision_vars
         sub_vars = sub.fwd_map[var]
         vals = rBM_sol[var]
-        for k in 1:length(sub_vars)
+        w = weights[var]
+        for k in eachindex(sub_vars)
             JuMP.add_to_expression!(obj_expr,
-                (sub_vars[k] - vals[k]) *
-                (sub_vars[k] - vals[k])
-                )
+                w[k] * (sub_vars[k] - vals[k]) * (sub_vars[k] - vals[k]))
         end
     end
     JuMP.@objective(sub.model, Min, obj_expr)
@@ -54,34 +55,33 @@ end
 
 # Solve the separation problem. Returns (separation_obj, separation_sol).
 function _solve_separation(
-    separation::GDPSubmodel,
+    reform_state::_CuttingPlanes,
     rBM_sol::Dict{<:JuMP.AbstractVariableRef, <:Vector{<:Number}}
     )
-    _set_separation_objective(separation, rBM_sol)
+    separation = reform_state.separation
+    _set_separation_objective(separation, reform_state.weights, rBM_sol)
     JuMP.optimize!(separation.model, ignore_optimize_hook = true)
     separation_obj = JuMP.objective_value(separation.model)
     separation_sol = extract_solution(separation)
     return separation_obj, separation_sol
 end
 
-# Add cut: Σ_var Σ_k 2*(sep_k - rBM_k)*(x_k - sep_k) ≥ 0
+# Add cut: Σ_var Σ_k 2*w_k*(sep_k - rBM_k)*(x_k - sep_k) ≥ 0
 function add_cut(
     model::JuMP.AbstractModel,
-    decision_vars::Vector{<:JuMP.AbstractVariableRef},
+    reform_state::_CuttingPlanes{O, T, V},
     rBM_sol::Dict{<:JuMP.AbstractVariableRef,<:Vector{<:Number}},
     separation_sol::Dict{<:JuMP.AbstractVariableRef,<:Vector{<:Number}}
-    )
-    cut_expr = zero(JuMP.GenericAffExpr{
-        JuMP.value_type(typeof(model)),
-        JuMP.variable_ref_type(model)})
-    for var in decision_vars
+    ) where {O, T, V}
+    cut_expr = zero(JuMP.GenericAffExpr{T, V})
+    for var in reform_state.decision_vars
         rbm_vals = rBM_sol[var]
         sep_vals = separation_sol[var]
-        for k in 1:length(rbm_vals)
-            xi = 2 * (sep_vals[k] - rbm_vals[k])
+        w = reform_state.weights[var]
+        for k in eachindex(rbm_vals)
+            xi = 2 * w[k] * (sep_vals[k] - rbm_vals[k])
             JuMP.add_to_expression!(cut_expr, xi, var)
-            JuMP.add_to_expression!(
-                cut_expr, -xi * sep_vals[k])
+            JuMP.add_to_expression!(cut_expr, -xi * sep_vals[k])
         end
     end
     cref = JuMP.@constraint(model, cut_expr >= 0)
@@ -98,27 +98,29 @@ function reformulate_model(
     method::CuttingPlanes
     )
     _clear_reformulations(model)
-    decision_vars = collect_cutting_planes_vars(model)
+    reform_state = _CuttingPlanes(method, model)
 
     # Build separation subproblem from the clean (unreformulated) model
-    separation = copy_and_reformulate(model, decision_vars, Hull(), method)
-    JuMP.relax_integrality(separation.model)
+    reform_state.separation = copy_and_reformulate(
+        model, reform_state.decision_vars, Hull(), reform_state)
+    JuMP.relax_integrality(reform_state.separation.model)
 
     # rBM: BigM in-place, relax logical vars
-    reformulate_model(model, BigM(method.M_value))
-    JuMP.set_optimizer(model, method.optimizer)
+    reformulate_model(model, BigM(reform_state.M_value))
+    JuMP.set_optimizer(model, reform_state.optimizer)
     JuMP.set_silent(model)
     relaxed_vars = relax_logical_vars(model)
 
     # Cutting plane loop: rBM <-> SEP until convergence
-    for iter in 1:method.max_iter
+    for iter in 1:reform_state.max_iter
         JuMP.optimize!(model, ignore_optimize_hook = true)
-        rBM_sol = extract_solution(model)
-        separation_obj, separation_sol = _solve_separation(separation, rBM_sol)
-        if separation_obj <= method.seperation_tolerance
+        rBM_sol = extract_solution(model, reform_state)
+        separation_obj, separation_sol =
+            _solve_separation(reform_state, rBM_sol)
+        if separation_obj <= reform_state.separation_tolerance
             break
         end
-        add_cut(model, decision_vars, rBM_sol, separation_sol)
+        add_cut(model, reform_state, rBM_sol, separation_sol)
     end
 
     unrelax_logical_vars(relaxed_vars)

src/datatypes.jl

@@ -445,9 +445,10 @@ 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 
+- `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`).
+- `M_value::T`: Big-M value of the relaxed BigM model the loop cuts
+against (default = `1e9`).
 """
 struct CuttingPlanes{O, T} <: AbstractReformulationMethod
     optimizer::O;
@@ -462,7 +463,8 @@ struct CuttingPlanes{O, T} <: AbstractReformulationMethod
         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,
+            final_reform_method, M_value)
     end
 end
 
@@ -494,6 +496,35 @@ struct GDPSubmodel{M <: JuMP.AbstractModel,
     fwd_map::Dict{V, Vector{W}}
 end
 
+# Per-run working state for one cutting planes solve: config copied
+# from `CuttingPlanes` plus the decision variables, the separation
+# submodel, and the support weights the separation objective and the
+# cuts share (unit weights for finite models; extensions store
+# quadrature weights at submodel build time).
+mutable struct _CuttingPlanes{O, T, V <: JuMP.AbstractVariableRef} <:
+        AbstractReformulationMethod
+    optimizer::O
+    max_iter::Int
+    separation_tolerance::T
+    final_reform_method::AbstractReformulationMethod
+    M_value::T
+    decision_vars::Vector{V}
+    weights::Dict{V, Vector{T}}
+    separation::Union{Nothing, GDPSubmodel{<:JuMP.AbstractModel, V}}
+
+    function _CuttingPlanes(
+        method::CuttingPlanes{O, T},
+        model::JuMP.AbstractModel
+        ) where {O, T}
+        dvars = collect_cutting_planes_vars(model)
+        V = eltype(dvars)
+        new{O, T, V}(method.optimizer, method.max_iter,
+            method.seperation_tolerance, method.final_reform_method,
+            method.M_value, dvars,
+            Dict{V, Vector{T}}(v => ones(T, 1) for v in dvars), nothing)
+    end
+end
+
 """
     PSplit <: AbstractReformulationMethod
 
@@ -681,4 +712,4 @@ A `VariableProperties` object with blank info.
 function VariableProperties(expr)
     info = _free_variable_info()
     return VariableProperties(info, "", nothing, nothing)
-end
+end