diff --git a/pyomo/contrib/pyros/CHANGELOG.txt b/pyomo/contrib/pyros/CHANGELOG.txt index 58b051a5e2c..dd9a46df066 100644 --- a/pyomo/contrib/pyros/CHANGELOG.txt +++ b/pyomo/contrib/pyros/CHANGELOG.txt @@ -8,6 +8,7 @@ PyROS 1.3.14 12 May 2026 ------------------------------------------------------------------------------- - Add `CartesianProductSet` class to facilitate representations of Cartesian products of continuous uncertainty sets +- Add caching for uncertainty set exact coordinate value bounds ------------------------------------------------------------------------------- diff --git a/pyomo/contrib/pyros/pyros.py b/pyomo/contrib/pyros/pyros.py index 46e2fb35583..d62a97ad834 100644 --- a/pyomo/contrib/pyros/pyros.py +++ b/pyomo/contrib/pyros/pyros.py @@ -386,10 +386,13 @@ def solve( ) model_data = ModelData(original_model=model, timing=TimingData(), config=None) - with time_code( - timing_data_obj=model_data.timing, - code_block_name="main", - is_main_timer=True, + with ( + uncertainty_set._cache_manager(), + time_code( + timing_data_obj=model_data.timing, + code_block_name="main", + is_main_timer=True, + ), ): kwds.update( dict( diff --git a/pyomo/contrib/pyros/tests/test_grcs.py b/pyomo/contrib/pyros/tests/test_grcs.py index 4f13da7ed4f..410593c82b0 100644 --- a/pyomo/contrib/pyros/tests/test_grcs.py +++ b/pyomo/contrib/pyros/tests/test_grcs.py @@ -3985,7 +3985,7 @@ def solve(self, model, *args, **kwargs): class TestPyROSUnavailableSubsolvers(unittest.TestCase): """ - Check that appropriate exceptionsa are raised if + Check that appropriate exceptions are raised if PyROS is invoked with unavailable subsolvers. """ @@ -5232,5 +5232,242 @@ def test_discrete_set_subsolver_error_recovery(self, name, sec_con_UB): ) +# @SolverFactory.register("slow_solver") +class SlowSolver: + """ + Solver which sleeps for a specified time before solving. + """ + + def __init__(self, sleep_time, sub_solver): + self.sleep_time = sleep_time + self.sub_solver = sub_solver + + self.options = Bunch() + + def available(self, exception_flag=True): + return True + + def license_is_valid(self): + return True + + def __enter__(self): + return self + + def __exit__(self, et, ev, tb): + pass + + def solve(self, model, **kwargs): + """ + Sleep, then solve a model. + + Parameters + ---------- + model : ConcreteModel + Model of interest. + + Returns + ------- + results : SolverResults + Solver results. + """ + + # ensure only one active objective + active_objs = [ + obj for obj in model.component_data_objects(Objective, active=True) + ] + assert len(active_objs) == 1 + + # sleep for specified time + time.sleep(self.sleep_time) + + print("I slept. Now, I will solve.") + # invoke subsolver + results = self.sub_solver.solve(model, **kwargs) + + return results + + +class CustomExactBoundsUncertaintySet(BoxSet): + """ + A box set, modified such that the `parameter_bounds` + attribute calculation involves globally solving the coordinate + value bounding problems. + """ + + def __init__(self, bounds, sleep_time, cache): + super().__init__(bounds) + self.sleep_time = sleep_time + self.cache = cache + + @property + def parameter_bounds(self): + """ + Solve bounding problems to calculate exact parameter bounds. + """ + solver = SlowSolver( + sub_solver=SolverFactory("ipopt"), sleep_time=self.sleep_time + ) + bounds = self._compute_exact_parameter_bounds(solver=solver) + + if not self.cache: + self._cache.clear() + + return bounds + + +@unittest.skipUnless(ipopt_available, "IPOPT is not available.") +class TestPyROSCacheUncertaintySetBounds(unittest.TestCase): + """ + Test behavior of PyROS solver with respect to caching of + an uncertainty set's exact coordinate value bounds. + """ + + def test_pyros_cache_creation(self): + """ + Check that PyROS creates a cache for storing computed exact parameter bounds. + """ + m = build_leyffer_two_cons() + + # Define the uncertainty set + interval = CustomExactBoundsUncertaintySet( + bounds=[(0.25, 2)], sleep_time=0, cache=True + ) + + # Instantiate the PyROS solver + pyros_solver = SolverFactory("pyros") + + # Define subsolvers utilized in the algorithm + local_subsolver = SolverFactory("ipopt") + global_subsolver = SolverFactory("ipopt") + + # check cache exists + self.assertTrue(hasattr(interval, "_cache")) + + # Call the PyROS solver + pyros_solver.solve( + model=m, + first_stage_variables=[m.x1], + second_stage_variables=[m.x2], + uncertain_params=[m.u], + uncertainty_set=interval, + local_solver=local_subsolver, + global_solver=global_subsolver, + options={ + "objective_focus": ObjectiveType.worst_case, + "solve_master_globally": True, + }, + ) + + # check cache has been cleared after solve + self.assertTrue(hasattr(interval, "_cache")) + self.assertEqual( + interval._cache, {}, msg="Did not clear uncertainty set cache after solve." + ) + + def test_pyros_cache_time(self): + """ + Check that the PyROS solve time is improved when + the uncertainty set's exact coordinate value bounds are cached. + """ + m = build_leyffer_two_cons() + + # Define the uncertainty set + interval_cache = CustomExactBoundsUncertaintySet( + bounds=[(0.25, 2)], sleep_time=0.1, cache=True + ) + interval_no_cache = CustomExactBoundsUncertaintySet( + bounds=[(0.25, 2)], sleep_time=0.1, cache=False + ) + + # Instantiate the PyROS solver + pyros_solver = SolverFactory("pyros") + + # Define subsolvers utilized in the algorithm + local_subsolver = SolverFactory("ipopt") + global_subsolver = SolverFactory("ipopt") + + # Call the PyROS solver + results_cache = pyros_solver.solve( + model=m, + first_stage_variables=[m.x1], + second_stage_variables=[m.x2], + uncertain_params=[m.u], + uncertainty_set=interval_cache, + local_solver=local_subsolver, + global_solver=global_subsolver, + options={ + "objective_focus": ObjectiveType.worst_case, + "solve_master_globally": True, + }, + ) + + results_no_cache = pyros_solver.solve( + model=m, + first_stage_variables=[m.x1], + second_stage_variables=[m.x2], + uncertain_params=[m.u], + uncertainty_set=interval_no_cache, + local_solver=local_subsolver, + global_solver=global_subsolver, + options={ + "objective_focus": ObjectiveType.worst_case, + "solve_master_globally": True, + }, + ) + + # caching should always result in less time, + # as not caching reruns the slow solver multiple times + self.assertGreater(results_no_cache.time, results_cache.time) + + def test_pyros_solve_assert_bounds_cache_empty(self): + """ + Check that calling PyROS on an uncertainty set + with a nonempty bounds cache results in an exception. + """ + m = build_leyffer_two_cons() + + # Define the uncertainty set + interval = CustomExactBoundsUncertaintySet( + bounds=[(25, 200)], sleep_time=0, cache=True + ) + self.assertAlmostEqual(interval.parameter_bounds[0][0], 25, places=2) + self.assertAlmostEqual(interval.parameter_bounds[0][1], 200, places=2) + + # change set attributes, leading to outdated parameter bounds + interval.bounds = [(0.25, 2)] + self.assertAlmostEqual(interval.parameter_bounds[0][0], 25, places=2) + self.assertAlmostEqual(interval.parameter_bounds[0][1], 200, places=2) + + # Instantiate the PyROS solver + pyros_solver = SolverFactory("pyros") + + # Define subsolvers utilized in the algorithm + local_subsolver = SolverFactory("ipopt") + global_subsolver = SolverFactory("ipopt") + + # Solve with PyROS + exc_str = ( + r"Nonempty cache for CustomExactBoundsUncertaintySet " + r"exact parameter bounds." + ) + with self.assertRaisesRegex(AssertionError, exc_str): + pyros_solver.solve( + model=m, + first_stage_variables=[m.x1], + second_stage_variables=[m.x2], + uncertain_params=[m.u], + uncertainty_set=interval, + local_solver=local_subsolver, + global_solver=global_subsolver, + options={ + "objective_focus": ObjectiveType.worst_case, + "solve_master_globally": True, + }, + ) + + self.assertAlmostEqual(interval._cache[0, minimize], 25, places=2) + self.assertAlmostEqual(interval._cache[0, maximize], 200, places=2) + + if __name__ == "__main__": unittest.main() diff --git a/pyomo/contrib/pyros/tests/test_uncertainty_sets.py b/pyomo/contrib/pyros/tests/test_uncertainty_sets.py index 5923a5c3ba5..9542dc875fe 100644 --- a/pyomo/contrib/pyros/tests/test_uncertainty_sets.py +++ b/pyomo/contrib/pyros/tests/test_uncertainty_sets.py @@ -24,7 +24,7 @@ ) from pyomo.environ import SolverFactory -from pyomo.core.base import ConcreteModel, Param, Var +from pyomo.core.base import ConcreteModel, Param, Var, maximize, minimize, UnitInterval from pyomo.core.expr import RangedExpression from pyomo.core.expr.compare import assertExpressionsEqual @@ -3736,6 +3736,37 @@ def parameter_bounds(self, val): self._parameter_bounds = val +class CustomDomainUncertaintySet(CustomUncertaintySet): + """ + Test simple custom uncertainty set with specified uncertain parameter domains. + """ + + def __init__(self, dim): + self._dim = dim + self._parameter_bounds = [(0, 1)] * self.dim + + def set_as_constraint(self, uncertain_params=None, block=None): + blk, param_var_list, conlist, aux_vars = ( + _setup_standard_uncertainty_set_constraint_block( + block=block, + uncertain_param_vars=uncertain_params, + dim=self.dim, + num_auxiliary_vars=None, + ) + ) + conlist.add(sum(param_var_list) <= 1) + for var in param_var_list: + conlist.add(0 <= var) + var.domain = UnitInterval + + return UncertaintyQuantification( + block=blk, + uncertainty_cons=list(conlist.values()), + uncertain_param_vars=param_var_list, + auxiliary_vars=aux_vars, + ) + + class TestCustomUncertaintySet(unittest.TestCase): """ Test for a custom uncertainty set subclass. @@ -3768,6 +3799,33 @@ def test_compute_exact_parameter_bounds(self): custom_set._compute_exact_parameter_bounds(baron), [(-1, 1)] * 2 ) + @unittest.skipUnless(baron_available, "BARON is not available") + def test_solve_exact_bounds_optimization(self): + """ + Test parameter bounds computations are cached. + """ + baron = SolverFactory("baron") + custom_set = CustomUncertaintySet(dim=2) + + # check cache exists + self.assertTrue(hasattr(custom_set, "_cache")) + + # check bounds calculation + self.assertEqual( + custom_set._solve_exact_bounds_optimization( + custom_set._create_bounding_model(), 0, minimize, baron + ), + -1.0, + ) + + # check cache access + self.assertIs( + custom_set._solve_exact_bounds_optimization( + custom_set._create_bounding_model(), 0, minimize, baron + ), + custom_set._cache[0, minimize], + ) + @unittest.skipUnless(baron_available, "BARON is not available") def test_solve_feasibility(self): """ @@ -3837,6 +3895,15 @@ def test_is_nonempty(self): with self.assertRaisesRegex(ValueError, exc_str): custom_set.is_nonempty(config=CONFIG) + def test_fbbt_values(self): + """ + Test that `_fbbt_parameter_bounds` returns correct values. + """ + CONFIG = pyros_config() + custom_set = CustomDomainUncertaintySet(dim=2) + + self.assertEqual(custom_set._fbbt_parameter_bounds(config=CONFIG), [(0, 1)] * 2) + @unittest.skipUnless(baron_available, "BARON is not available") def test_is_coordinate_fixed(self): """ diff --git a/pyomo/contrib/pyros/uncertainty_sets.py b/pyomo/contrib/pyros/uncertainty_sets.py index bf230d1ffea..0b22847bcb6 100644 --- a/pyomo/contrib/pyros/uncertainty_sets.py +++ b/pyomo/contrib/pyros/uncertainty_sets.py @@ -17,6 +17,7 @@ """ import abc +import contextlib import math import functools import itertools @@ -524,6 +525,31 @@ def parameter_bounds(self): """ raise NotImplementedError + @property + def _cache(self): + """ + dict : Cache for the bounds defining the minimum bounding box + of `self`. Each key is a 2-tuple containing the positional + index of a coordinate and an + :class:`~pyomo.common.enums.ObjectiveSense` object + indicating the type of bound (`minimize` for a lower bound, + `maximize for an upper bound) being specified for the + coordinate. + """ + try: + return self.__cache + except AttributeError: + self.__cache = {} + return self.__cache + + @contextlib.contextmanager + def _cache_manager(self): + assert ( + not self._cache + ), f"Nonempty cache for {self.__class__.__name__} exact parameter bounds." + yield self + self._cache.clear() + def _create_bounding_model(self): """ Make uncertain parameter value bounding problems (optimize @@ -667,6 +693,7 @@ def validate(self, config): ValueError If nonemptiness check or boundedness check fails. """ + # perform validation checks if not self.is_nonempty(config=config): raise ValueError(f"Nonemptiness check failed for uncertainty set {self}.") @@ -790,44 +817,106 @@ def _compute_exact_parameter_bounds(self, solver, index=None): if index is None: index = [(True, True)] * self.dim - # create bounding model and get all objectives + # create bounding model bounding_model = self._create_bounding_model() - objs_to_optimize = bounding_model.param_var_objectives.items() param_bounds = [] - for idx, obj in objs_to_optimize: - # activate objective for corresponding dimension - obj.activate() + for idx in range(self.dim): bounds = [] - - # solve for lower bound, then upper bound - # solve should be successful for i, sense in enumerate((minimize, maximize)): - # check if the LB or UB should be solved if not index[idx][i]: bounds.append(None) continue - obj.sense = sense - res = solver.solve(bounding_model, load_solutions=False) - if check_optimal_termination(res): - bounding_model.solutions.load_from(res) - else: - raise ValueError( - "Could not compute " - f"{'lower' if sense == minimize else 'upper'} " - f"bound in dimension {idx + 1} of {self.dim}. " - f"Solver status summary:\n {res.solver}." + # NOTE: variables are not initialized for the first + # exact bounds optimization, and are initialized to + # the solution of the most recent iteration at each + # subsequent evaluation. + bounds.append( + self._solve_exact_bounds_optimization( + bounding_model, idx, sense, solver ) - bounds.append(value(obj)) + ) # add parameter bounds for current dimension param_bounds.append(tuple(bounds)) + return param_bounds + + def _solve_exact_bounds_optimization(self, bounding_model, index, sense, solver): + """ + Compute an exact lower or upper bound for a specified + coordinate of the points contained in `self`, by solving + a bounding model. + + For efficiency, the result is cached if the bounding model + is solved successfully. Further, if the cache already contains + an entry corresponding to the coordinate bound of interest, + then that entry is returned and solution of the bounding model + is skipped. + + Parameters + ---------- + bounding_model : ConcreteModel + Bounding model, with an indexed minimization sense + Objective with name 'param_var_objectives' consisting + of `N` entries, all of which have been deactivated. + index : int + The positional index for the coordinate of interest. + sense : ~pyomo.common.ObjectiveSense + Optimization sense for the bounding model objective. + This also indicates the type of bound (lower or upper) + to be computed. + Select `minimize` to compute the lower bound or + `maximize` to compute the upper bound. + solver : ~pyomo.opt.base.solvers.OptSolver + Optimizer to invoke on the bounding model. + + Returns + ------- + bound : float + A value of the lower or upper bound for + the corresponding dimension at the specified index. + + Raises + ------ + ValueError + If there was an unsuccessful attempt to solve + the bounding model. + """ + # we use saved optimization results for a given index + # for either `maximize` UB or `minimize` LB if they exist + if (index, sense) in self._cache: + return self._cache[index, sense] + + # select objective corresponding to specified index + obj = bounding_model.param_var_objectives[index] + obj.activate() + + # optimize with specified sense + obj.sense = sense + try: + res = solver.solve(bounding_model, load_solutions=False) + finally: # ensure sense is minimize when done, deactivate obj.sense = minimize obj.deactivate() - return param_bounds + if check_optimal_termination(res): + bounding_model.solutions.load_from(res) + else: + raise ValueError( + "Could not compute " + f"{'lower' if sense == minimize else 'upper'} " + f"bound in dimension {index + 1} of {self.dim}. " + f"Solver status summary:\n {res.solver}." + ) + + bound = value(obj) + + # store in cache + self._cache[index, sense] = bound + + return bound def _fbbt_parameter_bounds(self, config): """ @@ -859,9 +948,7 @@ def _fbbt_parameter_bounds(self, config): f"{fbbt_infeasible_con_exception!r}" ) - param_bounds = [ - (var.lower, var.upper) for var in bounding_model.param_vars.values() - ] + param_bounds = [(var.lb, var.ub) for var in bounding_model.param_vars.values()] return param_bounds @@ -3991,6 +4078,19 @@ def parameter_bounds(self): return [] return parameter_bounds + @contextlib.contextmanager + def _cache_manager(self): + with contextlib.ExitStack() as stack: + # Verify this (CartesianProductSet's) cache is empty + stack.enter_context(super()._cache_manager()) + for uset in self._all_sets: + # Verify all component caches are empty + stack.enter_context(uset._cache_manager()) + yield self + # This will re-enter when this context manager is exited, + # which will exit the stack context, triggering all the + # context managers entered above to exit. + def _iterate_over_all_sets(self): """ Iterate over the sets contained in `self`.