diff --git a/doc/OnlineDocs/explanation/solvers/pyros/uncertainty_sets.rst b/doc/OnlineDocs/explanation/solvers/pyros/uncertainty_sets.rst index 2188ac14fec..cb948b271de 100644 --- a/doc/OnlineDocs/explanation/solvers/pyros/uncertainty_sets.rst +++ b/doc/OnlineDocs/explanation/solvers/pyros/uncertainty_sets.rst @@ -44,6 +44,7 @@ subclasses are enumerated below: ~pyomo.contrib.pyros.uncertainty_sets.BoxSet ~pyomo.contrib.pyros.uncertainty_sets.BudgetSet ~pyomo.contrib.pyros.uncertainty_sets.CardinalitySet + ~pyomo.contrib.pyros.uncertainty_sets.CartesianProductSet ~pyomo.contrib.pyros.uncertainty_sets.DiscreteScenarioSet ~pyomo.contrib.pyros.uncertainty_sets.EllipsoidalSet ~pyomo.contrib.pyros.uncertainty_sets.FactorModelSet @@ -76,6 +77,9 @@ subclasses are provided below. * - :class:`~pyomo.contrib.pyros.uncertainty_sets.CardinalitySet` - :math:`\begin{array}{l} q^{0} \in \mathbb{R}^{n}, \\ \hat{q} \in \mathbb{R}_{+}^{n}, \\ \Gamma \in [0, n] \end{array}` - :math:`\left\{ q \in \mathbb{R}^{n} \middle| \begin{array}{l} \exists\,\xi \in [0, 1]^n\,:\\ \quad \,q = q^{0} + \hat{q} \circ \xi \\ \quad \displaystyle \sum_{i=1}^{n} \xi_{i} \leq \Gamma \end{array} \right\}` + * - :class:`~pyomo.contrib.pyros.uncertainty_sets.CartesianProductSet` + - :math:`\begin{array}{l} \mathcal{Q}_{1} \subset \mathbb{R}^{n_1}, \\ \mathcal{Q}_{2} \subset \mathbb{R}^{n_2}, \\ \vdots \\ \mathcal{Q}_{m} \subset \mathbb{R}^{n_m} \end{array}` + - :math:`\displaystyle \mathcal{Q}_{1} \times \mathcal{Q}_{2} \times \cdots \times \mathcal{Q}_{m}` * - :class:`~pyomo.contrib.pyros.uncertainty_sets.DiscreteScenarioSet` - :math:`q^{1}, q^{2},\dots , q^{S} \in \mathbb{R}^{n}` - :math:`\{q^{1}, q^{2}, \dots , q^{S}\}` diff --git a/pyomo/contrib/pyros/__init__.py b/pyomo/contrib/pyros/__init__.py index c067fb7b10d..0061c94c60d 100644 --- a/pyomo/contrib/pyros/__init__.py +++ b/pyomo/contrib/pyros/__init__.py @@ -20,4 +20,5 @@ BoxSet, IntersectionSet, AxisAlignedEllipsoidalSet, + CartesianProductSet, ) diff --git a/pyomo/contrib/pyros/tests/test_grcs.py b/pyomo/contrib/pyros/tests/test_grcs.py index 609c84cf2cc..4f13da7ed4f 100644 --- a/pyomo/contrib/pyros/tests/test_grcs.py +++ b/pyomo/contrib/pyros/tests/test_grcs.py @@ -69,6 +69,8 @@ _setup_standard_uncertainty_set_constraint_block, AxisAlignedEllipsoidalSet, BoxSet, + CardinalitySet, + CartesianProductSet, DiscreteScenarioSet, FactorModelSet, Geometry, @@ -460,6 +462,90 @@ def test_two_stg_model_discrete_set(self): ) +@unittest.skipUnless(ipopt_available, "IPOPT is not available.") +class TestPyROSSolveCartesianProductSet(unittest.TestCase): + """ + Test PyROS successfully solves model with cartesian product uncertainty. + """ + + def build_model(self): + m = ConcreteModel() + m.q = Param(range(4), initialize=0, mutable=True) + m.x = Var(initialize=0, bounds=(0, 100)) + m.obj = Objective(expr=m.x, sense=minimize) + m.con = Constraint(expr=m.x >= sum(m.q.values())) + return m + + def test_solve_with_cartesian_product_set(self): + """ + Test solve with cartesian product uncertainty. + """ + m = self.build_model() + cpset = CartesianProductSet( + [ + BoxSet([[0, 1]]), + FactorModelSet( + origin=[0, 0], number_of_factors=1, beta=1, psi_mat=[[1], [3]] + ), + CardinalitySet(origin=[0], positive_deviation=[0.5], gamma=1), + ] + ) + results = SolverFactory("pyros").solve( + model=m, + first_stage_variables=[m.x], + second_stage_variables=[], + uncertain_params=m.q, + uncertainty_set=cpset, + local_solver=SolverFactory("ipopt"), + global_solver=SolverFactory("ipopt"), + options={ + "objective_focus": ObjectiveType.worst_case, + "solve_master_globally": True, + }, + ) + + # check successful termination + self.assertEqual( + results.pyros_termination_condition, + pyrosTerminationCondition.robust_optimal, + msg="Did not identify robust optimal solution to problem instance.", + ) + # need only 2 iterations: + # - first gives nominally optimal x, which is robust infeasible + # - second gives robust optimal x + self.assertEqual(results.iterations, 2) + # final objective is just sum of the parameter values + # that cause maximal violation of the inequality constraint + self.assertAlmostEqual(results.final_objective_value, 5.5, places=6) + self.assertAlmostEqual(m.x.value, 5.5, places=6) + + def test_solve_invalid_cartesian_product_set(self): + """ + Test PyROS solver cartesian product set validation failure. + """ + m = self.build_model() + cpset = CartesianProductSet( + [BoxSet([[0, 1]]), DiscreteScenarioSet([(0, 1, 3), (0, 0, 0)])] + ) + # exception raised during set validation due to involvement + # of discrete uncertainty in the cartesian product + disc_exc_str = r"CartesianProductSet.*entry.*with a discrete geometry" + with self.assertRaisesRegex(ValueError, disc_exc_str): + SolverFactory("pyros").solve( + model=m, + first_stage_variables=[m.x], + second_stage_variables=[], + uncertain_params=m.q, + uncertainty_set=cpset, + local_solver=SolverFactory("ipopt"), + global_solver=SolverFactory("ipopt"), + options={ + "objective_focus": ObjectiveType.worst_case, + "solve_master_globally": True, + }, + ) + + class TestPyROSRobustInfeasible(unittest.TestCase): @unittest.skipUnless(baron_available, "BARON is not available and licensed") def test_pyros_robust_infeasible(self): diff --git a/pyomo/contrib/pyros/tests/test_uncertainty_sets.py b/pyomo/contrib/pyros/tests/test_uncertainty_sets.py index 28e974598e3..5923a5c3ba5 100644 --- a/pyomo/contrib/pyros/tests/test_uncertainty_sets.py +++ b/pyomo/contrib/pyros/tests/test_uncertainty_sets.py @@ -32,6 +32,7 @@ AxisAlignedEllipsoidalSet, BoxSet, BudgetSet, + CartesianProductSet, CardinalitySet, DiscreteScenarioSet, EllipsoidalSet, @@ -3179,6 +3180,512 @@ def test_is_coordinate_fixed(self): ) +class TestCartesianProductSet(unittest.TestCase): + """ + Tests for CartesianProductSet. + """ + + def test_normal_construction_and_update(self): + """ + Test CartesianProductSet constructor works as expected. + """ + bset = BoxSet(bounds=[[-1, 1], [-1, 1]]) + aset = AxisAlignedEllipsoidalSet([0, 0, 0, 0], [1, 1, 1, 1]) + + cpset = CartesianProductSet([bset, aset]) + self.assertIs( + bset, + cpset._all_sets[0], + msg=( + "CartesianProductSet 'all_sets' attribute does not " + "contain expected BoxSet" + ), + ) + self.assertIs( + aset, + cpset._all_sets[1], + msg=( + "CartesianProductSet 'all_sets' attribute does not " + "contain expected AxisAlignedEllipsoidalSet" + ), + ) + # check defined attributes/methods inherited from base class + self.assertIs(cpset.geometry, Geometry.CONVEX_NONLINEAR) + self.assertEqual(cpset.type, "cartesian_product") + self.assertEqual(cpset.dim, 6) + + exc_str = ( + r"CartesianProductSet has an entry.*1 that is not of type UncertaintySet" + ) + with self.assertRaisesRegex(TypeError, exc_str): + CartesianProductSet([BoxSet([[0, 1]]), 1]) + + # iterable should be a sequence, and the constructor performs + # the iterable type check before doing anything else + iter_exc_str = r"`all_sets`.*Sequence.*but is of type set" + with self.assertRaisesRegex(TypeError, iter_exc_str): + CartesianProductSet({BoxSet([[0, 1]]), 1}) + + def test_set_as_constraint(self): + """ + Test method for setting up Cartesian product constraints + works correctly. + """ + m = ConcreteModel() + m.v = Var(range(8), initialize=0) + cpset = CartesianProductSet( + [ + BoxSet([(-0.5, 0.5)]), + FactorModelSet( + origin=[0, 1], number_of_factors=1, beta=0.75, psi_mat=[[1], [3]] + ), + CardinalitySet([-0.5, -0.5], [2, 2], 2), + AxisAlignedEllipsoidalSet([0, 0, 0], [0.25, 0.25, 0.25]), + ] + ) + + uq = cpset.set_as_constraint(uncertain_params=m.v, block=m) + + self.assertIs(uq.block, m) + self.assertEqual(uq.uncertain_param_vars, list(m.v.values())) + self.assertEqual(len(uq.auxiliary_vars), 3) + self.assertEqual(len(uq.uncertainty_cons), 8) + + # box constraint + assertExpressionsEqual( + self, + uq.uncertainty_cons[0].expr, + RangedExpression((np.float64(-0.5), m.v[0], np.float64(0.5)), False), + ) + + # factor model constraints + aux_vars = uq.auxiliary_vars + assertExpressionsEqual(self, uq.uncertainty_cons[1].expr, aux_vars[0] == m.v[1]) + assertExpressionsEqual( + self, uq.uncertainty_cons[2].expr, 1 + 3 * aux_vars[0] == m.v[2] + ) + assertExpressionsEqual( + self, + uq.uncertainty_cons[3].expr, + RangedExpression((-0.75, aux_vars[0], 0.75), False), + ) + self.assertEqual(aux_vars[0].bounds, (-1, 1)) + + # cardinality set constraints + assertExpressionsEqual( + self, uq.uncertainty_cons[4].expr, -0.5 + 2 * aux_vars[1] == m.v[3] + ) + assertExpressionsEqual( + self, uq.uncertainty_cons[5].expr, -0.5 + 2 * aux_vars[2] == m.v[4] + ) + assertExpressionsEqual( + self, uq.uncertainty_cons[6].expr, sum(aux_vars[1:3]) <= 2 + ) + self.assertEqual(aux_vars[1].bounds, (0, 1)) + self.assertEqual(aux_vars[2].bounds, (0, 1)) + + # axis-aligned ellipsoid constraint + assertExpressionsEqual( + self, + uq.uncertainty_cons[7].expr, + ( + m.v[5] ** 2 / np.float64(0.0625) + + m.v[6] ** 2 / np.float64(0.0625) + + m.v[7] ** 2 / np.float64(0.0625) + <= 1 + ), + ) + + def test_set_as_constraint_dim_mismatch(self): + """ + Check exception raised when writing Cartesian product constraints + if number of uncertain parameters does not match the dimension. + """ + m = ConcreteModel() + m.v1 = Var(initialize=0) + m.v2 = Var(initialize=0) + cpset = CartesianProductSet( + [BoxSet(bounds=[[1, 2], [3, 4]]), AxisAlignedEllipsoidalSet([0, 1], [5, 5])] + ) + with self.assertRaisesRegex(ValueError, ".*dimension"): + cpset.set_as_constraint(uncertain_params=[m.v1, m.v2], block=m) + + def test_set_as_constraint_type_mismatch(self): + """ + Check exception raised in Cartesian product set constraint + building if uncertain parameter variables are of invalid type. + """ + m = ConcreteModel() + m.p1 = Param([0, 1], initialize=0, mutable=True) + cpset = CartesianProductSet( + [BoxSet(bounds=[[1, 2], [3, 4]]), AxisAlignedEllipsoidalSet([0, 1], [5, 5])] + ) + with self.assertRaisesRegex(TypeError, ".*valid component type"): + cpset.set_as_constraint(uncertain_params=[m.p1[0], m.p1[1]], block=m) + + with self.assertRaisesRegex(TypeError, ".*valid component type"): + cpset.set_as_constraint(uncertain_params=m.p1, block=m) + + @unittest.skipUnless(baron_available, "BARON is not available.") + def test_compute_exact_parameter_bounds(self): + """ + Test Cartesian product set exact parameter bounds + computations give expected results. + """ + cpset = CartesianProductSet( + [ + BoxSet([(-0.5, 0.5)]), + FactorModelSet( + origin=[0, 0], + number_of_factors=2, + beta=0.75, + psi_mat=[[1, 1], [1, 2]], + ), + CardinalitySet([-0.5, -0.5], [2, 2], 2), + AxisAlignedEllipsoidalSet([0, 0, 1], [0.25, 0.8, 0.25]), + ] + ) + + computed_bounds = cpset._compute_exact_parameter_bounds(SolverFactory("baron")) + np.testing.assert_allclose( + computed_bounds, + [ + # box bounds + [-0.5, 0.5], + # factor model bounds + [-1.5, 1.5], + [-2.5, 2.5], + # cardinality bounds + [-0.5, 1.5], + [-0.5, 1.5], + # ellipsoid bounds + [-0.25, 0.25], + [-0.8, 0.8], + [0.75, 1.25], + ], + ) + # since all sets in the product allow for quick bounds, + # also check for parity + np.testing.assert_allclose(cpset.parameter_bounds, computed_bounds) + + # check that response to `index` argument is as expected + partial_index = [ + (True, False), + (False, True), + (False, False), + (True, True), + (False, False), + (True, False), + (False, True), + (True, True), + ] + partial_computed_bounds = cpset._compute_exact_parameter_bounds( + SolverFactory("baron"), index=partial_index + ) + partial_index_arr = np.array(partial_index) + for idx1, idx2 in zip(*np.where(partial_index_arr)): + self.assertAlmostEqual( + partial_computed_bounds[idx1][idx2], + computed_bounds[idx1][idx2], + msg=( + f"Bound for index {idx1, idx2} should be " + f"{computed_bounds[idx1][idx2]}, " + f"instead got {partial_computed_bounds[idx1][idx2]}" + ), + ) + for idx1, idx2 in zip(*np.where(~partial_index_arr)): + self.assertIsNone( + partial_computed_bounds[idx1][idx2], + msg=( + f"Bound for index {idx1, idx2} should be None, " + f"instead got {partial_computed_bounds[idx1][idx2]}" + ), + ) + + def test_parameter_bounds(self): + """ + Test CartesianProductSet `parameter_bounds` method works + as expected. + """ + cpset = CartesianProductSet( + [ + BoxSet([(-0.5, 0.5)]), + CardinalitySet([-0.5, -0.5], [2, 2], 2), + AxisAlignedEllipsoidalSet([0, 0, 1], [0.25, 0.8, 0.25]), + ] + ) + self.assertTrue(cpset._PARAMETER_BOUNDS_EXACT) + np.testing.assert_allclose( + cpset.parameter_bounds, + [ + # box bounds + [-0.5, 0.5], + # cardinality bounds + [-0.5, 1.5], + [-0.5, 1.5], + # ellipsoid bounds + [-0.25, 0.25], + [-0.8, 0.8], + [0.75, 1.25], + ], + ) + + # polyhedral set doesn't provide parameter bounds, + # so neither should cartesian product + cpset2 = CartesianProductSet( + [BoxSet([(0, 1)]), PolyhedralSet([[1, 1], [1, 1]], [1, 1])] + ) + self.assertFalse(cpset2.parameter_bounds) + self.assertFalse(cpset2._PARAMETER_BOUNDS_EXACT) + + # polyhedral set doesn't provide parameter bounds, + # so neither should cartesian product + cpset3 = CartesianProductSet( + [ + BoxSet([(0, 1)]), + IntersectionSet( + set1=BoxSet([[-1, 1], [-1, 1]]), + set2=AxisAlignedEllipsoidalSet([0, 0], [1, 1]), + ), + ] + ) + self.assertFalse(cpset3.parameter_bounds) + self.assertFalse(cpset3._PARAMETER_BOUNDS_EXACT) + + def test_point_in_set(self): + """ + Test Cartesian product set membership check. + """ + cpset = CartesianProductSet( + [BoxSet([(-0.5, 0.5)]), AxisAlignedEllipsoidalSet([0, 0], [0.25, 0.25])] + ) + + # in both sets + self.assertTrue(cpset.point_in_set([-0.5] + [0, -0.25])) + self.assertTrue(cpset.point_in_set([-0.5] + [0, 0.25])) + self.assertTrue(cpset.point_in_set([-0.5] + [-0.25, 0])) + self.assertTrue(cpset.point_in_set([-0.5] + [0.25, 0])) + self.assertTrue(cpset.point_in_set([-0.5] + [0, 0])) + self.assertTrue(cpset.point_in_set([0.5] + [0, -0.25])) + self.assertTrue(cpset.point_in_set([0.5] + [0, 0.25])) + self.assertTrue(cpset.point_in_set([0.5] + [-0.25, 0])) + self.assertTrue(cpset.point_in_set([0.5] + [0.25, 0])) + self.assertTrue(cpset.point_in_set([0.5] + [0, 0])) + + # in box set, outside ellipsoid + self.assertFalse(cpset.point_in_set([-0.5] + [-0.25, -0.25])) + self.assertFalse(cpset.point_in_set([-0.5] + [-0.25, 0.25])) + self.assertFalse(cpset.point_in_set([-0.5] + [0.25, -0.25])) + self.assertFalse(cpset.point_in_set([-0.5] + [0.25, 0.25])) + self.assertFalse(cpset.point_in_set([0.5] + [-0.25, -0.25])) + self.assertFalse(cpset.point_in_set([0.5] + [-0.25, 0.25])) + self.assertFalse(cpset.point_in_set([0.5] + [0.25, -0.25])) + self.assertFalse(cpset.point_in_set([0.5] + [0.25, 0.25])) + + # outside box, in ellipsoid + self.assertFalse(cpset.point_in_set([-0.6] + [0, -0.25])) + self.assertFalse(cpset.point_in_set([-0.6] + [0, 0.25])) + self.assertFalse(cpset.point_in_set([-0.6] + [-0.25, 0])) + self.assertFalse(cpset.point_in_set([-0.6] + [0.25, 0])) + self.assertFalse(cpset.point_in_set([-0.6] + [0, 0])) + + # outside both sets + self.assertFalse(cpset.point_in_set([-0.6] + [-0.25, -0.25])) + self.assertFalse(cpset.point_in_set([-0.6] + [-0.25, 0.25])) + self.assertFalse(cpset.point_in_set([-0.6] + [0.25, -0.25])) + self.assertFalse(cpset.point_in_set([-0.6] + [0.25, 0.25])) + + def test_add_bounds_on_uncertain_parameters(self): + m = ConcreteModel() + m.uncertain_param_vars = Var(range(6), initialize=0) + cpset = CartesianProductSet( + [ + BoxSet([(-0.5, 0.5)]), + CardinalitySet([-0.5, -0.5], [2, 2], 2), + AxisAlignedEllipsoidalSet([0, 0, 1], [0.25, 0.8, 0.25]), + ] + ) + cpset._add_bounds_on_uncertain_parameters( + uncertain_param_vars=m.uncertain_param_vars + ) + + # check bounds + np.testing.assert_allclose(m.uncertain_param_vars[0].bounds, (-0.5, 0.5)) + np.testing.assert_allclose(m.uncertain_param_vars[1].bounds, (-0.5, 1.5)) + np.testing.assert_allclose(m.uncertain_param_vars[2].bounds, (-0.5, 1.5)) + np.testing.assert_allclose(m.uncertain_param_vars[3].bounds, (-0.25, 0.25)) + np.testing.assert_allclose(m.uncertain_param_vars[4].bounds, (-0.8, 0.8)) + np.testing.assert_allclose(m.uncertain_param_vars[5].bounds, (0.75, 1.25)) + + # check exception raised if there is dimension mismatch + exc_str = r"Passed 2 VarData objects representing.* but.*of dimension 6." + with self.assertRaisesRegex(ValueError, exc_str): + cpset._add_bounds_on_uncertain_parameters( + uncertain_param_vars=[ + m.uncertain_param_vars[0], + m.uncertain_param_vars[1], + ] + ) + + def test_validate(self): + """ + Test Cartesian product set validation methods work as expected. + """ + CONFIG = pyros_config() + + # works if all sets are valid and none are discrete + bset = BoxSet(bounds=[[-1, 1]]) + aset = AxisAlignedEllipsoidalSet([0, 0, 0], [1, 1, 1]) + CartesianProductSet([bset, aset]).validate(CONFIG) + + # works if otherwise valid and nominal values provided + CONFIG.nominal_uncertain_param_vals = [0, 0.5, 0.5, 0.5] + CartesianProductSet([bset, aset]).validate(CONFIG) + # check that state of CONFIG is unchanged + self.assertEqual(CONFIG.nominal_uncertain_param_vals, [0, 0.5, 0.5, 0.5]) + + # allow repeated sets (set powers) + CONFIG.nominal_uncertain_param_vals = None + CartesianProductSet([bset, bset]).validate(CONFIG) + + # fails if a discrete set is involved in the product + disc_exc_str = r"CartesianProductSet.*entry.*with a discrete geometry" + with self.assertRaisesRegex(ValueError, disc_exc_str): + CartesianProductSet([bset, DiscreteScenarioSet([(0,), (1,)])]).validate( + CONFIG + ) + with self.assertRaisesRegex(ValueError, disc_exc_str): + CartesianProductSet( + [ + bset, + IntersectionSet( + set1=DiscreteScenarioSet([(0, 0), (0, 1)]), + set2=BoxSet([[0, 1]] * 2), + ), + ] + ).validate(CONFIG) + + # fails if at least one set is invalid + exc_str = "Lower bound.*exceeds upper bound" + with self.assertRaisesRegex(ValueError, exc_str): + # second box set invalid due to failed bounds + CartesianProductSet([bset, BoxSet([[1, 0], [0, 0]])]).validate(CONFIG) + + @unittest.skipUnless(baron_available, "BARON not available") + def test_is_coordinate_fixed(self): + """ + Test Cartesian product set method for checking whether + there are coordinates constrained to a single value. + """ + bset = BoxSet([[0, 0], [-1, 1]]) + aset = AxisAlignedEllipsoidalSet([0, 0], [1, 0]) + self.assertEqual( + CartesianProductSet([bset, aset])._is_coordinate_fixed( + # don't need a global solver, since exact bounds + # are given by the `parameter_bounds` method + config=Bunch() + ), + [True, False, False, True], + ) + + iset = IntersectionSet(set1=aset, set2=BoxSet([(0, 1), (0, 1)])) + self.assertEqual( + CartesianProductSet([bset, iset])._is_coordinate_fixed( + # need global solver to compute intersection set bounds + config=Bunch(global_solver=SolverFactory("baron")) + ), + [True, False, False, True], + ) + + @unittest.skipUnless(baron_available, "BARON not available") + def test_compute_auxiliary_param_vals(self): + """ + Test computations of Cartesian product set + auxiliary uncertain parameter values. + """ + # for case where all sets do not use auxiliary parameters, + # the return value should just be an empty 1D array + self.assertEqual( + CartesianProductSet( + [BoxSet([[0, 1]]), AxisAlignedEllipsoidalSet([0, 0], [1, 1])] + ) + .compute_auxiliary_uncertain_param_vals([0, 0, 0]) + .shape, + (0,), + ) + + # product of sets involving cardinality and factor model: + # should just reduce to concatenation of individual + # set calculations + cpset = CartesianProductSet( + [ + BoxSet([(-0.5, 0.5)]), + FactorModelSet( + origin=[0, 1], number_of_factors=1, beta=0.75, psi_mat=[[1], [4]] + ), + CardinalitySet([-0.5, -0.5], [2, 2], 1), + AxisAlignedEllipsoidalSet([0, 0, 0], [0.25, 0.25, 0.25]), + ] + ) + np.testing.assert_allclose( + cpset.compute_auxiliary_uncertain_param_vals( + [0.5] + [0, 1] + [-0.5, -0.5] + [0, 0, 0] + ), + [0, 0, 0], + ) + # deviations from factor model origin + np.testing.assert_allclose( + cpset.compute_auxiliary_uncertain_param_vals( + [0.5] + [0.75, 4] + [-0.5, -0.5] + [0, 0, 0] + ), + [0.75, 0, 0], + ) + np.testing.assert_allclose( + cpset.compute_auxiliary_uncertain_param_vals( + [0.5] + [-0.75, -2] + [-0.5, -0.5] + [0, 0, 0] + ), + [-0.75, 0, 0], + ) + # deviations from cardinality origin + np.testing.assert_allclose( + cpset.compute_auxiliary_uncertain_param_vals( + [0.5] + [0, 1] + [1.5, -0.5] + [0, 0, 0] + ), + [0, 1, 0], + ) + np.testing.assert_allclose( + cpset.compute_auxiliary_uncertain_param_vals( + [0.5] + [0, 1] + [-0.5, 1.5] + [0, 0, 0] + ), + [0, 0, 1], + ) + # deviations from cardinality and factor model origins + np.testing.assert_allclose( + cpset.compute_auxiliary_uncertain_param_vals( + [0.5] + [0.75, 4] + [-0.5, 1.5] + [0, 0, 0] + ), + [0.75, 0, 1], + ) + np.testing.assert_allclose( + cpset.compute_auxiliary_uncertain_param_vals( + [0.5] + [-0.75, -2] + [-0.5, 1.5] + [0, 0, 0] + ), + [-0.75, 0, 1], + ) + np.testing.assert_allclose( + cpset.compute_auxiliary_uncertain_param_vals( + [0.5] + [0.75, 4] + [1.5, -0.5] + [0, 0, 0] + ), + [0.75, 1, 0], + ) + np.testing.assert_allclose( + cpset.compute_auxiliary_uncertain_param_vals( + [0.5] + [-0.75, -2] + [1.5, -0.5] + [0, 0, 0] + ), + [-0.75, 1, 0], + ) + + class CustomUncertaintySet(UncertaintySet): """ Test simple custom uncertainty set subclass. diff --git a/pyomo/contrib/pyros/uncertainty_sets.py b/pyomo/contrib/pyros/uncertainty_sets.py index a925ab19396..6c6cdbc15c4 100644 --- a/pyomo/contrib/pyros/uncertainty_sets.py +++ b/pyomo/contrib/pyros/uncertainty_sets.py @@ -22,7 +22,7 @@ import itertools from numbers import Integral from collections import namedtuple -from collections.abc import Iterable, MutableSequence +from collections.abc import Iterable, MutableSequence, Sequence from enum import Enum from pyomo.common.dependencies import numpy as np, scipy as sp @@ -756,7 +756,7 @@ def point_in_set(self, point): def _compute_exact_parameter_bounds(self, solver, index=None): """ - Compute lower and upper coordinate value bounds + Compute specified tight lower and upper coordinate value bounds for every dimension of `self` by solving a bounding model. Parameters @@ -774,10 +774,12 @@ def _compute_exact_parameter_bounds(self, solver, index=None): Returns ------- - param_bounds : list of tuple of float - Each entry of the list is a 2-tuple - containing the lower and upper bound for - the corresponding dimension. + param_bounds : list of tuple + Every entry of the list is a 2-tuple, + each member of which is the corresponding dimension's + lower/upper bound + (if the corresponding entry of `index` is True) + or None (if the corresponding entry of `index` is False). Raises ------ @@ -971,7 +973,7 @@ def _is_coordinate_fixed(self, config, index=None): config : ConfigDict PyROS solver options. Should at least contain attribute `global_solver`. - index : iterable of int, optional + index : Sequence[int] | None Positional indices of the coordinates to check. If `None` is passed, then `index` is set to ``list(range(self.dim))``, so that all coordinates @@ -3849,3 +3851,300 @@ def validate(self, config): # check boundedness and nonemptiness of intersected set super().validate(config) + + +class CartesianProductSet(UncertaintySet): + """ + A Cartesian product of uncertainty sets. + + The order and identities of the uncertainty sets + involved in the Cartesian product are immutable, + and all sets in the product should be non-discrete. + + Parameters + ---------- + all_sets : Sequence[UncertaintySet] + Uncertainty sets of which the product is to be taken. + + Raises + ------ + TypeError + If any entry of ``all_sets`` is not of type `UncertaintySet`. + + Notes + ----- + Given uncertainty sets + :math:`\\mathcal{Q}_1 \\in \\mathbb{R}^{n_1}`, + :math:`\\mathcal{Q}_2 \\in \\mathbb{R}^{n_2}`, + :math:`\\dots`, + :math:`\\mathcal{Q}_m \\in \\mathbb{R}^{n_m}`, + collectively represented by the argument ``all_sets``, + the :math:`(n_1 + n_2 + \\dots + n_m)`-dimensional + Cartesian product set is defined by + + .. math:: + + \\mathcal{Q}_1 \\times \\mathcal{Q}_2 \\times \\cdots + \\times \\mathcal{Q}_m. + + Examples + -------- + Cartesian product of 1D box/interval and 2D + hypersphere (circle): + + >>> from pyomo.contrib.pyros import ( + ... BoxSet, AxisAlignedEllipsoidalSet, CartesianProductSet, + ... ) + >>> interval = BoxSet(bounds=[[-1.5, 1.5]]) + >>> circle = AxisAlignedEllipsoidalSet( + ... center=[0, 0], + ... half_lengths=[2, 2], + ... ) + >>> cartesian_product = CartesianProductSet([interval, circle]) + """ + + def __init__(self, all_sets): + """Initialize self (see class docstring).""" + if not isinstance(all_sets, Sequence): + raise TypeError( + f"Argument `all_sets` should be a {Sequence.__name__}-type " + f"iterable, but is of type {type(all_sets).__name__}." + ) + all_sets = tuple(all_sets) + for val in all_sets: + if not isinstance(val, UncertaintySet): + raise TypeError( + f"{type(self).__name__} has an entry of value {val!r} " + "that is not of type " + f"{UncertaintySet.__name__}. " + "Ensure that all entries are of type " + f"{UncertaintySet.__name__}." + ) + + # protect this attribute to make the Cartesian product set, + # and thus the set's dimension, effectively immutable, as + # instances of the other (pre-implemented) uncertainty set + # types are also of immutable dimension + self._all_sets = all_sets + + @property + def type(self): + """ + str : Brief description of the type of the uncertainty set. + """ + return "cartesian_product" + + @property + def dim(self): + """ + int : Dimension of the cartesian product set. + """ + return sum(uset.dim for uset in self._all_sets) + + @property + def geometry(self): + """ + Geometry : Geometry of the Cartesian product set, + assuming that there are no discrete sets. + See the `Geometry` class documentation. + """ + return Geometry(max(uset.geometry.value for uset in self._all_sets)) + + @property + def _PARAMETER_BOUNDS_EXACT(self): + """ + bool : True if the coordinate value bounds returned by + ``self.parameter_bounds`` are exact + (i.e., specify the minimum bounding box), + False otherwise. + + For the cartesian product set, parameter bounds are exact iff + the parameter bounds for each multiplicand are exact. + """ + return all(uset._PARAMETER_BOUNDS_EXACT for uset in self._all_sets) + + @property + def parameter_bounds(self): + """ + Bounds for the value of each uncertain parameter constrained + by the set (i.e. bounds for each set dimension). + + Returns + ------- + list[tuple[numbers.Real, numbers.Real]] + If the ``parameter_bounds`` method returns a nonempty + list for all sets involved in the Cartesian product, + then this list is of length ``self.dim`` and contain the + (lower, upper) bound pairs. Otherwise, the list is empty. + """ + parameter_bounds = [] + for uset in self._all_sets: + # NOTE: by assumption, the list of parameter bounds for + # `uset` is either empty or of length equal to + # ``uset.dim``. + uset_bounds = uset.parameter_bounds + if uset_bounds: + parameter_bounds.extend(uset_bounds) + else: + return [] + return parameter_bounds + + def _iterate_over_all_sets(self): + """ + Iterate over the sets contained in `self`. + + Yields + ------ + start_dim : int + Positional index for the first dimension of the + multiplicand set iterate. + stop_dim : int + One plus the positional index for the last dimension of the + multiplicand set iterate. + uset : UncertaintySet + The multiplicand set iterate. + """ + starting_dim = 0 + for uset in self._all_sets: + yield starting_dim, starting_dim + uset.dim, uset + starting_dim += uset.dim + + @copy_docstring(UncertaintySet.point_in_set) + def point_in_set(self, point): + for start_dim, stop_dim, uset in self._iterate_over_all_sets(): + in_uset = uset.point_in_set(point[start_dim:stop_dim]) + if not in_uset: + return False + return True + + @copy_docstring(UncertaintySet.compute_auxiliary_uncertain_param_vals) + def compute_auxiliary_uncertain_param_vals(self, point, solver=None): + validate_array( + arr=point, + arr_name="point", + dim=1, + valid_types=native_numeric_types, + valid_type_desc="numeric type", + required_shape=[self.dim], + required_shape_qual="to match the set dimension", + ) + + aux_vals = [] + for start_dim, stop_dim, uset in self._iterate_over_all_sets(): + uset_pt = point[start_dim:stop_dim] + aux_vals.extend( + uset.compute_auxiliary_uncertain_param_vals(uset_pt, solver=solver) + ) + + return np.array(aux_vals) + + @copy_docstring(UncertaintySet.set_as_constraint) + def set_as_constraint(self, uncertain_params=None, block=None): + block, param_var_data_list, uncertainty_conlist, aux_var_list = ( + _setup_standard_uncertainty_set_constraint_block( + block=block, + uncertain_param_vars=uncertain_params, + dim=self.dim, + num_auxiliary_vars=None, + ) + ) + + all_cons, all_aux_vars = [], [] + for idx, (start_dim, stop_dim, uset) in enumerate( + self._iterate_over_all_sets() + ): + sub_block = Block() + block.add_component( + unique_component_name(block, f"sub_block_{idx}"), sub_block + ) + set_quantification = uset.set_as_constraint( + block=sub_block, + uncertain_params=param_var_data_list[start_dim:stop_dim], + ) + all_cons.extend(set_quantification.uncertainty_cons) + all_aux_vars.extend(set_quantification.auxiliary_vars) + + return UncertaintyQuantification( + block=block, + uncertain_param_vars=param_var_data_list, + uncertainty_cons=all_cons, + auxiliary_vars=all_aux_vars, + ) + + def _compute_exact_parameter_bounds(self, solver, index=None): + """ + Compute specified tight lower and upper coordinate value bounds + for every dimension of `self` by solving a bounding model. + + Parameters + ---------- + solver : Pyomo solver type + Optimizer to invoke on the bounding problems. + index : list of 2-tuple of bool, optional + A list of tuples for each index of the coordinates for + which to compute bounds. A lower or upper bound is + computed for any value that is True, while False + indicates that the bound should be skipped. + If None is passed, then the argument is set to + ``[(True, True)]*self.dim``, so that the bounds + for all coordinates are computed. + + Returns + ------- + param_bounds : list of tuple + Every entry of the list is a 2-tuple, + each member of which is the corresponding dimension's + lower/upper bound + (if the corresponding entry of `index` is True) + or None (if the corresponding entry of `index` is False). + """ + if index is None: + index = [(True, True)] * self.dim + param_bounds = [(None, None)] * self.dim + for start_dim, stop_dim, uset in self._iterate_over_all_sets(): + param_bounds[start_dim:stop_dim] = uset._compute_exact_parameter_bounds( + solver, index=index[start_dim:stop_dim] + ) + return param_bounds + + def validate(self, config): + """ + Validate the Cartesian product set instance. + + Parameters + ---------- + config : ConfigDict + PyROS solver configuration. + + Raises + ------ + ValueError + If any set involved in the product has a discrete geometry. + """ + + full_nom_param_vals = config.nominal_uncertain_param_vals + for start_dim, stop_dim, uset in self._iterate_over_all_sets(): + # ensure there are no discrete sets + if uset.geometry == Geometry.DISCRETE_SCENARIOS: + raise ValueError( + f"{type(self).__name__} has an entry {uset!r} " + "with a discrete geometry. " + "Ensure that all entries do not have discrete geometries." + ) + + # instead of using the default validation method + # on `self` (generally slow), we are going to separately + # validate each set in the product (possibly fast). + # as the check for each set may require the nominal values + # of the set's corresponding uncertain parameters, we + # need to temporarily update the appropriate config attribute + if full_nom_param_vals: + config.nominal_uncertain_param_vals = full_nom_param_vals[ + start_dim:stop_dim + ] + + try: + uset.validate(config) + finally: + # ensure the config's state ultimately remains unchanged + config.nominal_uncertain_param_vals = full_nom_param_vals