Clear outer_bound when a solve produces none, instead of keeping a stale one

mpisppy/cgbase.py

@@ -555,21 +555,37 @@ def add_columns_to_mp_from_results(self, all_results):
                 reduced cost, scenario cost, and variable values.
 
         Returns:
-            float: The sum of reduced costs
+            float or None: The sum of reduced costs, or None if any subproblem
+                reported no reduced cost.
+
+        build_columns_from_subproblem_solutions leaves red_cost at None when a
+        subproblem's solve produced no outer bound. Skipping such a term and
+        returning the sum of the rest would overstate the bound -- the caller
+        adds it to rmp_obj_val and treats the result as an outer bound -- so the
+        whole sum is reported unavailable instead. The columns themselves are
+        still added: they are useful regardless of whether a bound came with
+        them.
         """
         sum_redcosts=0
+        bound_available=True
         for rank_results in all_results:
             if rank_results is not None:
                 if isinstance(rank_results, list):
                     for result in rank_results:
                         sname, red_cost, scen_cost, xvec = result
                         self.add_column_for_scenario(sname, scen_cost, xvec)
-                        sum_redcosts+=red_cost
+                        if red_cost is None:
+                            bound_available=False
+                        else:
+                            sum_redcosts+=red_cost
                 else:
                     sname, red_cost, scen_cost, xvec = rank_results
                     self.add_column_for_scenario(sname, scen_cost, xvec)
-                    sum_redcosts+=red_cost
-        return sum_redcosts                
+                    if red_cost is None:
+                        bound_available=False
+                    else:
+                        sum_redcosts+=red_cost
+        return sum_redcosts if bound_available else None
 
     def _build_columns_from_xhat_list(self, xhat_list):
         """
@@ -674,22 +690,32 @@ def iterk_loop(self):
                 sum_redcosts=self.add_columns_to_mp_from_results(all_results)
                 self.add_columns_to_mp_from_results(all_results_xhat_recent) 
                 self.add_columns_to_mp_from_results(all_results_xfeas) 
-                self.outer_bound_candidate = self.rmp_obj_val + sum_redcosts
+                # A None sum means some subproblem produced no reduced cost, so
+                # rmp_obj_val + sum(reduced costs) is not an outer bound this
+                # iteration. Leave the incumbent best bound and the convergence
+                # metric where they are rather than computing from a partial
+                # sum; the next iteration may well have every term.
+                if sum_redcosts is None:
+                    self.outer_bound_candidate = None
+                else:
+                    self.outer_bound_candidate = self.rmp_obj_val + sum_redcosts
 
                 # Update bounds and convergence metric
                 if self.problem_is_lp:
                     self.best_solution_obj_val= self.rmp_obj_val
 
-                if self.best_bound_obj_val is None:
-                    self.best_bound_obj_val = self.outer_bound_candidate
-                elif self.is_minimizing:
-                    self.best_bound_obj_val = max(self.outer_bound_candidate, self.best_bound_obj_val)
-                else:
-                    self.best_bound_obj_val = min(self.outer_bound_candidate, self.best_bound_obj_val)
-                if self.is_minimizing:
-                    self.conv=(self.rmp_obj_val-self.best_bound_obj_val)/abs(self.best_bound_obj_val)
-                else:
-                    self.conv=(self.best_bound_obj_val-self.rmp_obj_val)/abs(self.best_bound_obj_val)
+                if self.outer_bound_candidate is not None:
+                    if self.best_bound_obj_val is None:
+                        self.best_bound_obj_val = self.outer_bound_candidate
+                    elif self.is_minimizing:
+                        self.best_bound_obj_val = max(self.outer_bound_candidate, self.best_bound_obj_val)
+                    else:
+                        self.best_bound_obj_val = min(self.outer_bound_candidate, self.best_bound_obj_val)
+                if self.best_bound_obj_val is not None:
+                    if self.is_minimizing:
+                        self.conv=(self.rmp_obj_val-self.best_bound_obj_val)/abs(self.best_bound_obj_val)
+                    else:
+                        self.conv=(self.best_bound_obj_val-self.rmp_obj_val)/abs(self.best_bound_obj_val)
 
                 if dprogress:
                     print("")

mpisppy/opt/dualcg.py

@@ -130,15 +130,20 @@ def update_dual_center(self):
         Update the regularization center if the outer bound candidate improved.
         """
 
+        # No bound this iteration means there is nothing to compare against and
+        # nothing to recenter on. Checked first: the comparison below would
+        # otherwise raise on None once a center has been established, which the
+        # old ordering only avoided while the center was still unset.
+        if self.outer_bound_candidate is None:
+            return
+
         if self.dual_center_outer_bound is not None:
             if self.is_minimizing:
                 improved = self.outer_bound_candidate >= self.dual_center_outer_bound
             else:
                 improved = self.outer_bound_candidate <= self.dual_center_outer_bound
             if not improved:
                 return
-        elif self.outer_bound_candidate is None:
-            return  
 
         m = self.mp  
 
@@ -241,23 +246,31 @@ def iterk_loop(self):
                 sum_redcosts=self.add_columns_to_mp_from_results(all_results)
                 self.add_columns_to_mp_from_results(all_results_xhat_recent) 
                 self.add_columns_to_mp_from_results(all_results_xfeas) 
-                self.outer_bound_candidate = self.rmp_obj_val + sum_redcosts
+                # See cgbase.add_columns_to_mp_from_results: None means some
+                # subproblem reported no reduced cost, so there is no outer
+                # bound to form from a partial sum this iteration.
+                if sum_redcosts is None:
+                    self.outer_bound_candidate = None
+                else:
+                    self.outer_bound_candidate = self.rmp_obj_val + sum_redcosts
                 if not hasattr(self, "dual_center_outer_bound"):
                     self.dual_center_outer_bound = None
 
                 self.update_dual_center()
 
                 # Update outer bound and convergence metric.
-                if self.best_bound_obj_val is None:
-                    self.best_bound_obj_val = self.outer_bound_candidate
-                elif self.is_minimizing:
-                    self.best_bound_obj_val = max(self.outer_bound_candidate, self.best_bound_obj_val)
-                else:
-                    self.best_bound_obj_val = min(self.outer_bound_candidate, self.best_bound_obj_val)
-                if self.is_minimizing:
-                    self.conv=(self.rmp_obj_val-self.best_bound_obj_val)/abs(self.best_bound_obj_val)
-                else:
-                    self.conv=(self.best_bound_obj_val-self.rmp_obj_val)/abs(self.best_bound_obj_val)
+                if self.outer_bound_candidate is not None:
+                    if self.best_bound_obj_val is None:
+                        self.best_bound_obj_val = self.outer_bound_candidate
+                    elif self.is_minimizing:
+                        self.best_bound_obj_val = max(self.outer_bound_candidate, self.best_bound_obj_val)
+                    else:
+                        self.best_bound_obj_val = min(self.outer_bound_candidate, self.best_bound_obj_val)
+                if self.best_bound_obj_val is not None:
+                    if self.is_minimizing:
+                        self.conv=(self.rmp_obj_val-self.best_bound_obj_val)/abs(self.best_bound_obj_val)
+                    else:
+                        self.conv=(self.best_bound_obj_val-self.rmp_obj_val)/abs(self.best_bound_obj_val)
 
                 if dprogress:
                     print("")

mpisppy/phbase.py

@@ -547,8 +547,10 @@ def post_solve_bound(self, solver_options=None, verbose=False):
                 If True, displays verbose output. Default False.
 
         Returns:
-            float:
-                An outer bound on the optimal objective function value.
+            float or None:
+                An outer bound on the optimal objective function value, or
+                None if any subproblem solve produced no bound (the
+                expectation cannot be formed if a scenario is missing one).
 
         Note:
             This function overwrites current variable values. This is only

mpisppy/spopt.py

@@ -381,11 +381,19 @@ def _vb(msg):
                         raise
 
                 if outer_bound is None:
-                    # Leave outer_bound at its previous value rather than
-                    # publish whatever the solver left in the results object.
-                    # The stale bound is still a valid outer bound (any
-                    # Lagrangian dual value bounds the original problem),
-                    # which is what the pre-outer_bound_only path did here.
+                    # This solve produced no bound, so say so. Keeping the
+                    # previous value would be wrong, even though each stale
+                    # value was a valid bound when it was computed: Ebound
+                    # sums p_s * outer_bound_s across scenarios, and a
+                    # Lagrangian bound is only valid as a *sum* when every
+                    # scenario uses weights from one generation satisfying
+                    # sum_s p_s W_s = 0. Mixing a stale bound with fresh ones
+                    # gives sum_s p_s L_s(W'_s) <= OPT + (sum_s p_s W'_s)^T xbar*,
+                    # whose trailing term does not vanish -- so the result is
+                    # not an outer bound at all, and nothing downstream can
+                    # tell. None instead lets Ebound's missing-bound check fire
+                    # and the spoke decline to send.
+                    s._mpisppy_data.outer_bound = None
                     if gripe:
                         print (f"[{self._get_cylinder_name()}] No outer bound for scenario {s.name}")
                         if results is not None:
@@ -399,6 +407,10 @@ def _vb(msg):
 
             elif sputils.not_good_enough_results(results):
                 s._mpisppy_data.solution_available = False
+                # Same reasoning as the outer_bound_only branch above: a failed
+                # solve computed no bound, and leaving the previous one in place
+                # lets Ebound form a mixed-generation sum that is not a bound.
+                s._mpisppy_data.outer_bound = None
 
                 if gripe:
                     print (f"[{self._get_cylinder_name()}] Solve failed for scenario {s.name}")

mpisppy/tests/test_cg_main.py

@@ -19,6 +19,7 @@
 from mpisppy.tests.examples.sizes.sizes import scenario_creator as sizes_creator, \
                                                scenario_denouement as sizes_denouement
 from mpisppy.tests.utils import get_solver, limit_solver_threads
+from mpisppy.cgbase import CGBase
 import mpisppy.MPI as mpi
 
 solver_available, solver_name, persistent_available, persistent_solver_name = get_solver(persistent_OK=False)
@@ -244,5 +245,65 @@ def test_sizes_iter0_creates_initial_columns(self):
             self.assertGreater(cg.next_col[sname], 0)
 
 
+
+class TestRedCostSumming(unittest.TestCase):
+    """add_columns_to_mp_from_results with a subproblem that reported no bound.
+
+    build_columns_from_subproblem_solutions already leaves red_cost at None
+    when a subproblem's solve produced no outer bound -- it tests the value for
+    finiteness precisely because that can happen. The consumer used to add it
+    unconditionally, which raised
+
+        TypeError: unsupported operand type(s) for +=: 'int' and 'NoneType'
+
+    the moment the case actually occurred (seen in CI, where the size-limited
+    cplex build fails solves that succeed on a full license).
+
+    Silently skipping the missing term would be worse than the crash: the caller
+    adds the sum to rmp_obj_val and treats the result as an outer bound, so
+    dropping a term overstates it. The sum reports itself unavailable instead.
+
+    No solver needed -- the method only calls self.add_column_for_scenario, so a
+    recorder stands in for the CG object.
+    """
+
+    class _Recorder:
+        def __init__(self):
+            self.added = []
+
+        def add_column_for_scenario(self, sname, scen_cost, xvec):
+            self.added.append(sname)
+
+    def _sum(self, all_results):
+        rec = self._Recorder()
+        total = CGBase.add_columns_to_mp_from_results(rec, all_results)
+        return total, rec.added
+
+    def test_all_bounds_present_sums_them(self):
+        results = [[("s0", 1.5, 10.0, {}), ("s1", 2.5, 20.0, {})]]
+        total, added = self._sum(results)
+        self.assertEqual(total, 4.0)
+        self.assertEqual(added, ["s0", "s1"])
+
+    def test_a_missing_bound_makes_the_sum_unavailable(self):
+        results = [[("s0", 1.5, 10.0, {}), ("s1", None, 20.0, {})]]
+        total, added = self._sum(results)
+        self.assertIsNone(total)
+        # The columns are still added: they are useful whether or not a bound
+        # came with them.
+        self.assertEqual(added, ["s0", "s1"])
+
+    def test_single_tuple_per_rank_is_handled_too(self):
+        # Ranks may report one tuple rather than a list of them.
+        self.assertEqual(self._sum([("s0", 3.0, 10.0, {})])[0], 3.0)
+        self.assertIsNone(self._sum([("s0", None, 10.0, {})])[0])
+
+    def test_none_rank_results_are_skipped(self):
+        total, added = self._sum([None, [("s0", 1.0, 10.0, {})]])
+        self.assertEqual(total, 1.0)
+        self.assertEqual(added, ["s0"])
+
+
+
 if __name__ == '__main__':
     unittest.main()

mpisppy/tests/test_outer_bound_only.py

@@ -19,6 +19,7 @@
 import unittest
 
 import numpy as np
+from pyomo.opt import SolverResults, SolverStatus, TerminationCondition
 
 import mpisppy.opt.ph
 import mpisppy.phbase
@@ -190,5 +191,101 @@ def test_normal_solve_still_loads_solution(self):
             self.assertTrue(np.isfinite(s._mpisppy_data.inner_bound))
 
 
+
+class _NoBoundPlugin:
+    """A solver plugin that comes back with no usable bound.
+
+    Only what solve_one touches: an options mapping and a solve() returning a
+    results object whose termination condition is one of the outcomes
+    no_outer_bound_results screens out. Not a persistent solver, so solve_one
+    takes the ordinary path.
+    """
+
+    def __init__(self, termination_condition):
+        self.options = {}
+        self._tc = termination_condition
+
+    def solve(self, s, **kwargs):
+        results = SolverResults()
+        results.solver.status = SolverStatus.warning
+        results.solver.termination_condition = self._tc
+        return results
+
+
+class TestStaleBoundNotRetained(unittest.TestCase):
+    """A solve that produces no bound must clear the subproblem's bound.
+
+    Keeping the previous value looks harmless -- each stale value really was a
+    valid outer bound for the weights it was computed with -- but Ebound sums
+    p_s * outer_bound_s across scenarios, and a Lagrangian bound is valid as a
+    *sum* only when every scenario uses weights from a single generation
+    satisfying sum_s p_s W_s = 0:
+
+        sum_s p_s L_s(W'_s) <= OPT + (sum_s p_s W'_s)^T xbar*
+
+    Mixing one scenario's stale bound with the others' fresh ones leaves that
+    trailing term in place, so the result is not an outer bound at all. Nothing
+    downstream can detect it: Ebound's missing-bound check looks for None, and a
+    retained number is not None. The hub then latches the best outer bound it
+    has seen, so a single too-good value is never corrected by later honest
+    ones.
+
+    No solver needed: the failing solve is stubbed.
+    """
+
+    def _ph_with_prior_bounds(self, value=-100.0):
+        """A PH object whose subproblems already carry bounds from an earlier
+        iteration's weights."""
+        ph = _make_ph()
+        for s in ph.local_scenarios.values():
+            s._mpisppy_data.outer_bound = value
+        return ph
+
+    def test_bound_only_infeasible_clears_the_bound(self):
+        ph = self._ph_with_prior_bounds()
+        name = list(ph.local_scenarios)[0]
+        s = ph.local_scenarios[name]
+        s._solver_plugin = _NoBoundPlugin(TerminationCondition.infeasible)
+        ph.solve_one(None, name, s, gripe=False, need_solution=False,
+                     outer_bound_only=True)
+        self.assertIsNone(s._mpisppy_data.outer_bound)
+
+    def test_bound_only_unbounded_clears_the_bound(self):
+        # An unbounded subproblem has Lagrangian value -inf, so a retained
+        # finite number over-claims regardless of any weight mixing.
+        ph = self._ph_with_prior_bounds()
+        name = list(ph.local_scenarios)[0]
+        s = ph.local_scenarios[name]
+        s._solver_plugin = _NoBoundPlugin(TerminationCondition.unbounded)
+        ph.solve_one(None, name, s, gripe=False, need_solution=False,
+                     outer_bound_only=True)
+        self.assertIsNone(s._mpisppy_data.outer_bound)
+
+    def test_ebound_declines_rather_than_mixing_generations(self):
+        # The end-to-end shape of the bug: one scenario keeps a bound from the
+        # previous weights while the other two get fresh ones. Ebound must
+        # refuse, not return a finite number that is not a bound.
+        ph = self._ph_with_prior_bounds()
+        names = list(ph.local_scenarios)
+        stale = ph.local_scenarios[names[0]]
+        stale._solver_plugin = _NoBoundPlugin(TerminationCondition.infeasible)
+        ph.solve_one(None, names[0], stale, gripe=False, need_solution=False,
+                     outer_bound_only=True)
+        for n in names[1:]:
+            ph.local_scenarios[n]._mpisppy_data.outer_bound = -50.0
+        self.assertIsNone(ph.Ebound())
+
+    def test_failed_ordinary_solve_also_clears_the_bound(self):
+        # The same exposure exists off the outer_bound_only path: the
+        # not_good_enough_results branch used to leave outer_bound untouched.
+        ph = self._ph_with_prior_bounds()
+        name = list(ph.local_scenarios)[0]
+        s = ph.local_scenarios[name]
+        s._solver_plugin = _NoBoundPlugin(TerminationCondition.infeasible)
+        ph.solve_one(None, name, s, gripe=False, need_solution=False)
+        self.assertIsNone(s._mpisppy_data.outer_bound)
+        self.assertFalse(s._mpisppy_data.solution_available)
+
+
 if __name__ == "__main__":
     unittest.main()