diff --git a/mpisppy/cgbase.py b/mpisppy/cgbase.py index 613edc123..768968fdd 100644 --- a/mpisppy/cgbase.py +++ b/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("") diff --git a/mpisppy/opt/dualcg.py b/mpisppy/opt/dualcg.py index de403b786..a29370f3a 100644 --- a/mpisppy/opt/dualcg.py +++ b/mpisppy/opt/dualcg.py @@ -130,6 +130,13 @@ 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 @@ -137,8 +144,6 @@ def update_dual_center(self): 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("") diff --git a/mpisppy/phbase.py b/mpisppy/phbase.py index 59a942120..e67696b28 100644 --- a/mpisppy/phbase.py +++ b/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 diff --git a/mpisppy/spopt.py b/mpisppy/spopt.py index f81eb475a..60a0efbdf 100644 --- a/mpisppy/spopt.py +++ b/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}") diff --git a/mpisppy/tests/test_cg_main.py b/mpisppy/tests/test_cg_main.py index cdb499fd6..12dfeb845 100644 --- a/mpisppy/tests/test_cg_main.py +++ b/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() diff --git a/mpisppy/tests/test_outer_bound_only.py b/mpisppy/tests/test_outer_bound_only.py index 41cec5d2e..59bf1cd96 100644 --- a/mpisppy/tests/test_outer_bound_only.py +++ b/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()