Fixed Bug with Upper Bound in PlasmoBenders

.github/workflows/test.yaml

@@ -29,6 +29,8 @@ jobs:
       - name: Run tests (PlasmoBenders)
         env:
           PKG: PlasmoBenders
+          # Set 1 thread for macOS, 2 for others
+          JULIA_NUM_THREADS: ${{ matrix.os == 'macos-latest' && 1 || 2 }}
         run: julia --color=yes .ci/ci.jl
 
   schwarz:
@@ -48,4 +50,5 @@ jobs:
       - name: Run tests (PlasmoSchwarz)
         env:
           PKG: PlasmoSchwarz
+          JULIA_NUM_THREADS: 1
         run: julia --color=yes .ci/ci.jl

lib/PlasmoBenders/Project.toml

@@ -1,7 +1,7 @@
 name = "PlasmoBenders"
 uuid = "491f1417-53b2-48aa-b9da-44cdd6c031b7"
 authors = ["David Cole"]
-version = "0.3.1"
+version = "0.3.2"
 
 [deps]
 Plasmo = "d3f7391f-f14a-50cc-bbe4-76a32d1bad3c"

lib/PlasmoBenders/src/Benders.jl

@@ -164,6 +164,7 @@ mutable struct BendersAlgorithm{T, V} <: AbstractPBAlgorithm{T, V}
     primal_iters::Dict{T, Matrix{Float64}}
     phis::Dict{T, Vector{Float64}}
     phis_LR::Dict{T, Vector{Float64}}
+    subgraph_objectives::Dict{T, Float64}
 
     time_forward_pass::Float64
     time_backward_pass::Float64
@@ -228,7 +229,8 @@ mutable struct BendersAlgorithm{T, V} <: AbstractPBAlgorithm{T, V}
         optimizer.primal_iters = Dict{T, Matrix{Float64}}()
         optimizer.phis = Dict{T, Vector{Float64}}()
         optimizer.phis_LR = Dict{T, Vector{Float64}}()
-
+        optimizer.subgraph_objectives = Dict{T, Float64}()
+
         optimizer.time_forward_pass = 0.
         optimizer.time_backward_pass = 0.
         optimizer.time_init = 0.
@@ -470,6 +472,7 @@ function BendersAlgorithm(
 
         # Add start object
         push!(optimizer.solve_order, root_object)
+        optimizer.subgraph_objectives[root_object] = 0.0
 
         _init_ext!(optimizer)
 
@@ -479,6 +482,7 @@ function BendersAlgorithm(
         while length(optimizer.ext["search_next"]) > 0
             search_next = optimizer.ext["search_next"][1]
             parent_object = optimizer.parent_objects[search_next]
+            optimizer.subgraph_objectives[search_next] = 0.0
             ############### Add complicating variables ##############
             # Get the linking constraints between last and current node
             _add_complicating_variables!(optimizer, parent_object, search_next, get_add_slacks(optimizer), get_slack_penalty(optimizer))
@@ -720,9 +724,9 @@ function _forward_pass!(optimizer::BendersAlgorithm)
     if get_regularize(optimizer)
         get_regularize_lbs(optimizer)[root_object] = lb[1]
 
-        _regularize_pass!(optimizer, root_object, ub)
+        _regularize_pass!(optimizer, root_object)
     else
-        _add_to_upper_bound!(optimizer, root_object, ub)
+        _save_subproblem_objective!(optimizer, root_object)
         # Save primal information to upcoming objects
         next_objects = optimizer.solve_order_dict[root_object]
         for object in next_objects
@@ -737,6 +741,8 @@ function _forward_pass!(optimizer::BendersAlgorithm)
         optimizer.last_solutions[root_object] = _get_object_last_solutions(root_object)
     end
 
+    ub[1] = optimizer.subgraph_objectives[root_object]
+
     ############# Solve each successive object #################
     if get_parallelize_benders(optimizer)
         _optimize_in_forward_pass_multithread!(optimizer, ub)

lib/PlasmoBenders/src/regularize.jl

@@ -40,8 +40,7 @@ end
 
 function _regularize_pass!(
     optimizer::BendersAlgorithm{T},
-    object,
-    ub
+    object::T
 ) where {T <: Union{Plasmo.OptiGraph, Plasmo.RemoteOptiGraph}}
     next_objects = optimizer.solve_order_dict[object]
     if length(next_objects) > 0
@@ -65,7 +64,7 @@ function _regularize_pass!(
 
         if termination_status(object) != MOI.INFEASIBLE
             obj_val_minus_theta = value(object, original_objective) - _theta_value(optimizer, object)
-            ub[1] += obj_val_minus_theta
+            optimizer.subgraph_objectives[object] = obj_val_minus_theta
             get_regularize_ubs(optimizer)[object] = obj_val_minus_theta
             for next_object in next_objects
                 comp_vars = optimizer.comp_vars[next_object]
@@ -102,7 +101,7 @@ function _regularize_pass!(
             delete!(object_dictionary(object), :_reg_con)
         end
     else
-        _add_to_upper_bound!(optimizer, object, ub)
+        _save_subproblem_objective!(optimizer, object)
         get_regularize_ubs(optimizer)[object] = JuMP.objective_value(object)
 
         if !optimizer.is_MIP #TODO: This can probably be moved to _save_forward_pass_solutions function; doesn't need to be here

lib/PlasmoBenders/src/solution.jl

@@ -96,8 +96,10 @@ end
 
 function _optimize_in_forward_pass_multithread!(optimizer::BendersAlgorithm{OptiGraph}, ub) #TODO: Type these functions
     Threads.@threads for i in 2:(length(optimizer.solve_order))
-        _forward_pass_iteration!(optimizer, i, ub)
+        _forward_pass_iteration!(optimizer, i)
     end
+    # sum the objective values from each thread and add to upper bound
+    ub[1] += sum(optimizer.subgraph_objectives[optimizer.solve_order[i]] for i in 2:(length(optimizer.solve_order)))
 end
 
 function _optimize_in_forward_pass_multithread!(optimizer::BendersAlgorithm{RemoteOptiGraph}, ub)
@@ -111,9 +113,9 @@ function _optimize_in_forward_pass_multithread!(optimizer::BendersAlgorithm{Remo
 
             optimizer.time_subproblem_solves += t_solve
             if !is_feasible
-                ub[1] = Inf
+                optimizer.subgraph_objectives[next_object] = Inf
             else
-                ub[1] += obj_val
+                optimizer.subgraph_objectives[next_object] = obj_val # TODO: if we parallelize beyond two stages, we need to subtract the value of theta here
             end
 
             optimizer.feasibility_map[next_object] = is_feasible
@@ -129,6 +131,8 @@ function _optimize_in_forward_pass_multithread!(optimizer::BendersAlgorithm{Remo
             end
         end
     end
+    # sum the objective values from each thread and add to upper bound
+    ub[1] += sum(optimizer.subgraph_objectives[optimizer.solve_order[i]] for i in 2:(length(optimizer.solve_order)))
 end
 
 function _optimize_remote_graph_forward(optimizer::BendersAlgorithm, next_object::RemoteOptiGraph, i::Int)
@@ -184,11 +188,13 @@ end
 
 function _optimize_in_forward_pass!(optimizer, ub) 
     for i in 2:(length(optimizer.solve_order))
-        _forward_pass_iteration!(optimizer, i, ub)
+        _forward_pass_iteration!(optimizer, i)
     end
+    # sum the objective values from each thread and add to upper bound
+    ub[1] += sum(optimizer.subgraph_objectives[optimizer.solve_order[i]] for i in 2:(length(optimizer.solve_order)))
 end
 
-function _forward_pass_iteration!(optimizer, i, ub)
+function _forward_pass_iteration!(optimizer, i)
     next_object = optimizer.solve_order[i]
 
     comp_vars = optimizer.comp_vars[next_object]
@@ -208,20 +214,21 @@ function _forward_pass_iteration!(optimizer, i, ub)
     object_termination_status = _check_termination_status(next_object, i; add_slacks_bool=get_add_slacks(optimizer), feasibility_cuts_bool=get_feasibility_cuts(optimizer))
 
     if object_termination_status
-        _save_forward_pass_solutions(optimizer, next_object, ub)
+        _save_forward_pass_solutions(optimizer, next_object)
         optimizer.feasibility_map[next_object] = true
     else
         # Need to do feasibility cuts; the check termination status function already tested
         # that feasibility_cuts was true
         println("Subgraph $i in forward pass was infeasible; using feasibility_cuts")
-        _save_feasibility_cut_data(optimizer, next_object, ub)
+        _save_feasibility_cut_data(optimizer, next_object)
+        optimizer.subgraph_objectives[next_object] = Inf
         optimizer.feasibility_map[next_object] = false
     end
 
     PlasmoBenders._unfix_variables(next_object, var_copies)
 end
 
-function _save_forward_pass_solutions(optimizer, next_object, ub)
+function _save_forward_pass_solutions(optimizer, next_object)
     # Add to the upper bound; if it's not the last object, subtract the cost-to-go from upper bound
     obj_val = JuMP.objective_value(next_object)
 
@@ -235,10 +242,10 @@ function _save_forward_pass_solutions(optimizer, next_object, ub)
     if get_regularize(optimizer)
         get_regularize_lbs(optimizer)[next_object] = obj_val
 
-        _regularize_pass!(optimizer, next_object, ub)
+        _regularize_pass!(optimizer, next_object)
     else
         # Save primal information to upcoming objects
-        _add_to_upper_bound!(optimizer, next_object, ub)
+        _save_subproblem_objective!(optimizer, next_object)
         next_objects = optimizer.solve_order_dict[next_object]
         for object in next_objects
             next_comp_vars = optimizer.comp_vars[object]
@@ -262,12 +269,10 @@ function _save_forward_pass_solutions(optimizer, next_object, ub)
     end
 end
 
-function _save_feasibility_cut_data(optimizer, next_object, ub)
+function _save_feasibility_cut_data(optimizer, next_object)
     # See JuMP documentation for implementation of this method:
     # https://jump.dev/JuMP.jl/stable/tutorials/algorithms/benders_decomposition/#Feasibility-cuts
 
-    ub[1] = Inf
-
     # Save primal information to upcoming objects
     if !optimizer.is_MIP
         obj_val = JuMP.dual_objective_value(next_object)
@@ -694,15 +699,15 @@ function _add_initial_relaxed_cuts!(
     optimizer.ext["link_var_mapping"] = link_vars_mapping
 end
 
-function _add_to_upper_bound!(
+function _save_subproblem_objective!(
     optimizer::BendersAlgorithm{T}, 
-    object::T, 
-    ub
+    object::T
 ) where {T <: Plasmo.AbstractOptiGraph}
     obj_val = JuMP.objective_value(object)
-    ub[1] += obj_val
     if length(optimizer.solve_order_dict[object]) > 0
         theta_val = _theta_value(optimizer, object)
-        ub[1] -= theta_val
+        optimizer.subgraph_objectives[object] = obj_val - theta_val
+    else
+        optimizer.subgraph_objectives[object] = obj_val
     end
 end