Merge branch 'main' into indexed-components-as-numeric

Author: jsiirola

pyomo/contrib/fbbt/tests/test_fbbt.py

@@ -469,8 +469,8 @@ def test_pow4(self):
                 else:
                     xu = m.x.ub
                 _x = np.exp(np.log(y) / _exp_val)
-                self.assertTrue(np.all(xl <= _x))
-                self.assertTrue(np.all(xu >= _x))
+                self.assertTrue(np.all(xl - 1e-14 <= _x))
+                self.assertTrue(np.all(xu + 1e-14 >= _x))
 
     def test_sqrt(self):
         m = pyo.ConcreteModel()
@@ -640,8 +640,8 @@ def test_log(self):
             else:
                 xu = pyo.value(m.x.ub)
             x = np.exp(z)
-            self.assertTrue(np.all(xl <= x))
-            self.assertTrue(np.all(xu >= x))
+            self.assertTrue(np.all(xl - 1e-14 <= x))
+            self.assertTrue(np.all(xu + 1e-14 >= x))
 
     def test_log10(self):
         if not numpy_available:

pyomo/contrib/pyros/CHANGELOG.txt

@@ -2,6 +2,12 @@
 PyROS CHANGELOG
 ===============
 
+-------------------------------------------------------------------------------
+PyROS 1.1.2    31 May 2022
+-------------------------------------------------------------------------------
+- Fixes to PyROS ellipsoidal sets.
+
+
 -------------------------------------------------------------------------------
 PyROS 1.1.1    25 Apr 2022
 -------------------------------------------------------------------------------

pyomo/contrib/pyros/pyros.py

@@ -45,7 +45,7 @@
 from pyomo.contrib.pyros.uncertainty_sets import uncertainty_sets
 from pyomo.core.base import Constraint
 
-__version__ =  "1.1.1"
+__version__ =  "1.1.2"
 
 def NonNegIntOrMinusOne(obj):
     '''

pyomo/contrib/pyros/tests/test_grcs.py

@@ -673,6 +673,33 @@ def test_add_bounds_on_uncertain_parameters(self):
         self.assertNotEqual(m.util.uncertain_param_vars[1].ub, None,
                             "Bounds not added correctly for EllipsoidalSet")
 
+    def test_ellipsoidal_set_bounds(self):
+        """Check `EllipsoidalSet` parameter bounds method correct."""
+        cov = [[2, 1], [1, 2]]
+        scales=[0.5, 2]
+        mean = [1, 1]
+
+        for scale in scales:
+            ell = EllipsoidalSet(center=mean, shape_matrix=cov, scale=scale)
+            bounds = ell.parameter_bounds
+            actual_bounds = list()
+            for idx, val in enumerate(mean):
+                diff = (cov[idx][idx] * scale) ** 0.5
+                actual_bounds.append((val - diff, val + diff))
+            self.assertTrue(
+                np.allclose(
+                    np.array(bounds),
+                    np.array(actual_bounds),
+                ),
+                msg=(
+                    f"EllipsoidalSet bounds {bounds} do not match their actual"
+                    f" values {actual_bounds} (for scale {scale}"
+                    f" and shape matrix {cov})."
+                    " Check the `parameter_bounds`"
+                    " method for the EllipsoidalSet."
+                ),
+            )
+
 class testAxisAlignedEllipsoidalUncertaintySetClass(unittest.TestCase):
     '''
     Axis aligned ellipsoidal uncertainty sets. Required inputs are half-lengths, nominal point, and right-hand side.
@@ -747,6 +774,102 @@ def test_add_bounds_on_uncertain_parameters(self):
         self.assertNotEqual(m.util.uncertain_param_vars[1].lb, None, "Bounds not added correctly for AxisAlignedEllipsoidalSet")
         self.assertNotEqual(m.util.uncertain_param_vars[1].ub, None, "Bounds not added correctly for AxisAlignedEllipsoidalSet")
 
+    def test_set_with_zero_half_lengths(self):
+        # construct ellipsoid
+        half_lengths = [1, 0, 2, 0]
+        center = [1, 1, 1, 1]
+        ell = AxisAlignedEllipsoidalSet(center, half_lengths)
+
+        # construct model
+        m = ConcreteModel()
+        m.v1 = Var()
+        m.v2 = Var([1, 2])
+        m.v3 = Var()
+
+        # test constraints
+        conlist = ell.set_as_constraint([m.v1, m.v2, m.v3])
+        eq_cons = [con for con in conlist.values() if con.equality]
+
+        self.assertEqual(
+            len(conlist),
+            3,
+            msg=(
+                "Constraint list for this `AxisAlignedEllipsoidalSet` should"
+                f" be of length 3, but is of length {len(conlist)}"
+            ),
+        )
+        self.assertEqual(
+            len(eq_cons),
+            2,
+            msg=(
+                "Number of equality constraints for this"
+                "`AxisAlignedEllipsoidalSet` should be 2,"
+                f" there are {len(eq_cons)} such constraints"
+            ),
+        )
+
+    @unittest.skipUnless(SolverFactory('baron').license_is_valid(),
+                         "Global NLP solver is not available and licensed.")
+    def test_two_stg_mod_with_axis_aligned_set(self):
+        """
+        Test two-stage model with `AxisAlignedEllipsoidalSet`
+        as the uncertainty set.
+        """
+        # define model
+        m = ConcreteModel()
+        m.x1 = Var(initialize=0, bounds=(0, None))
+        m.x2 = Var(initialize=0, bounds=(0, None))
+        m.x3 = Var(initialize=0, bounds=(None, None))
+        m.u1 = Param(initialize=1.125, mutable=True)
+        m.u2 = Param(initialize=1, mutable=True)
+
+        m.con1 = Constraint(expr=m.x1 * m.u1**(0.5) - m.x2 * m.u1 <= 2)
+        m.con2 = Constraint(expr=m.x1 ** 2 - m.x2 ** 2 * m.u1 == m.x3)
+
+        m.obj = Objective(expr=(m.x1 - 4) ** 2 + (m.x2 - m.u2) ** 2)
+
+        # Define the uncertainty set
+        # we take the parameter `u2` to be 'fixed'
+        ellipsoid = AxisAlignedEllipsoidalSet(
+            center=[1.125, 1],
+            half_lengths=[1, 0],
+        )
+
+        # Instantiate the PyROS solver
+        pyros_solver = SolverFactory("pyros")
+
+        # Define subsolvers utilized in the algorithm
+        local_subsolver = SolverFactory('baron')
+        global_subsolver = SolverFactory("baron")
+
+        # Call the PyROS solver
+        results = pyros_solver.solve(
+            model=m,
+            first_stage_variables=[m.x1, m.x2],
+            second_stage_variables=[],
+            uncertain_params=[m.u1, m.u2],
+            uncertainty_set=ellipsoid,
+            local_solver=local_subsolver,
+            global_solver=global_subsolver,
+            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."
+        )
+        self.assertGreater(
+            results.iterations,
+            0,
+            msg="Robust infeasible model terminated in 0 iterations (nominal case)."
+        )
+
+
 class testPolyhedralUncertaintySetClass(unittest.TestCase):
     '''
     Polyhedral uncertainty sets. Required inputs are matrix A, right-hand-side b, and list of uncertain params.

pyomo/contrib/pyros/uncertainty_sets.py

@@ -734,8 +734,8 @@ def __init__(self, center, half_lengths):
             raise AttributeError("Vector of half-lengths must be real-valued and numeric.")
         if not all(isinstance(elem, (int, float)) for elem in center):
             raise AttributeError("Vector center must be real-valued and numeric.")
-        if any(elem <= 0 for elem in half_lengths):
-            raise AttributeError("Half length values must be > 0.")
+        if any(elem < 0 for elem in half_lengths):
+            raise AttributeError("Half length values must be nonnegative.")
         # === Valid variance dimensions
         if not len(center) == len(half_lengths):
             raise AttributeError("Half lengths and center of ellipsoid must have same dimensions.")
@@ -765,33 +765,52 @@ def parameter_bounds(self):
         parameter_bounds = [(nom_value[i] - half_length[i], nom_value[i] + half_length[i]) for i in range(len(nom_value))]
         return parameter_bounds
 
-    def set_as_constraint(self, uncertain_params, **kwargs):
+    def set_as_constraint(self, uncertain_params, model=None, config=None):
         """
-        Function to generate constraints for the AxisAlignedEllipsoid uncertainty set.
+        Generate constraint(s) for the `AxisAlignedEllipsoidSet`
+        class.
 
         Args:
-           uncertain_params: uncertain parameter objects for writing constraint objects
-        """
-        if len(uncertain_params) != len(self.center):
-               raise AttributeError("Center of ellipsoid must be same dimensions as vector of uncertain parameters.")
-        # square and invert half lengths
-        inverse_squared_half_lengths = list(1.0/(a**2) for a in self.half_lengths)
-        # Calculate row vector of differences
-        diff_squared = []
-        # === Assume VarList uncertain_param_vars
-        for idx, i in enumerate(uncertain_params):
-           if uncertain_params[idx].is_indexed():
-               for index in uncertain_params[idx]:
-                   diff_squared.append((uncertain_params[idx][index] - self.center[idx])**2)
-           else:
-               diff_squared.append((uncertain_params[idx] - self.center[idx])**2)
-
-        # Calculate inner product of difference vector and variance matrix
-        constraint = sum([x * y for x, y in zip(inverse_squared_half_lengths, diff_squared)])
-
+            uncertain_params: uncertain parameter objects for writing
+            constraint objects. Indexed parameters are accepted, and
+            are unpacked for constraint generation.
+        """
+        all_params = list()
+
+        # expand all uncertain parameters to a list.
+        # this accounts for the cases in which `uncertain_params`
+        # consists of indexed model components,
+        # or is itself a single indexed component
+        if not isinstance(uncertain_params, (tuple, list)):
+            uncertain_params = [uncertain_params]
+
+        all_params = []
+        for uparam in uncertain_params:
+            all_params.extend(uparam.values())
+
+        if len(all_params) != len(self.center):
+            raise AttributeError(
+                f"Center of ellipsoid is of dimension {len(self.center)},"
+                f" but vector of uncertain parameters is of dimension"
+                f" {len(all_params)}"
+            )
+
+        zip_all = zip(all_params, self.center, self.half_lengths)
+        diffs_squared = list()
+
+        # now construct the constraints
         conlist = ConstraintList()
         conlist.construct()
-        conlist.add(constraint <= 1)
+        for param, ctr, half_len in zip_all:
+            if half_len > 0:
+                diffs_squared.append((param - ctr) ** 2 / (half_len) ** 2)
+            else:
+                # equality constraints for parameters corresponding to
+                # half-lengths of zero
+                conlist.add(param == ctr)
+
+        conlist.add(sum(diffs_squared) <= 1)
+
         return conlist
 
 
@@ -802,13 +821,19 @@ class EllipsoidalSet(UncertaintySet):
 
     def __init__(self, center, shape_matrix, scale=1):
         """
-        EllipsoidalSet constructor
+        EllipsoidalSet constructor.
 
-        Args:
-            center: Vector (``list``) of uncertain parameter values around which deviations are restrained.
-            shape_matrix: Positive semi-definite matrix, effectively a covariance matrix for
-            constraint and bounds determination.
-            scale: Right-hand side value for the ellipsoid.
+        Parameters
+        ----------
+        center : (N,) array-like
+            Center of the ellipsoid.
+        shape_matrix : (N, N) array-like
+            A positive definite matrix characterizing the shape
+            and orientation of the ellipsoid.
+        scale : float
+            Square of the factor by which to scale the semi-axes
+            of the ellipsoid (i.e. the eigenvectors of the covariance
+            matrix).
         """
 
         # === Valid data in lists/matrixes
@@ -875,8 +900,8 @@ def parameter_bounds(self):
         scale = self.scale
         nom_value = self.center
         P = self.shape_matrix
-        parameter_bounds = [(nom_value[i] - np.power(P[i][i], 0.5) * scale,
-                             nom_value[i] + np.power(P[i][i], 0.5) * scale) for i in range(self.dim)]
+        parameter_bounds = [(nom_value[i] - np.power(P[i][i] * scale, 0.5),
+                             nom_value[i] + np.power(P[i][i] * scale, 0.5)) for i in range(self.dim)]
         return parameter_bounds
 
     def set_as_constraint(self, uncertain_params, **kwargs):