Merge pull request #2620 from shermanjasonaf/fix-factor-model-bounds

Generalize and Fix Uncertain Parameter Bounds Evaluation for PyROS `FactorModelSet`

Author: blnicho

pyomo/contrib/pyros/CHANGELOG.txt

@@ -2,6 +2,14 @@
 PyROS CHANGELOG
 ===============
 
+-------------------------------------------------------------------------------
+PyROS 1.2.3    22 Nov 2022
+-------------------------------------------------------------------------------
+- Generalize FactorModelSet.
+- Resolve issues with FactorModelSet parameter bounds.
+- Modularize construction of uncertainty set bounding problems.
+
+
 -------------------------------------------------------------------------------
 PyROS 1.2.2    09 Nov 2022
 -------------------------------------------------------------------------------

pyomo/contrib/pyros/pyros.py

@@ -46,7 +46,7 @@
 from pyomo.core.base import Constraint
 
 
-__version__ = "1.2.2"
+__version__ = "1.2.3"
 
 
 def NonNegIntOrMinusOne(obj):

pyomo/contrib/pyros/tests/test_grcs.py

@@ -2460,6 +2460,93 @@ def test_error_on_invalid_beta(self):
         with self.assertRaisesRegex(ValueError, big_exc_str):
             fset.beta = big_beta
 
+    @unittest.skipUnless(
+        SolverFactory("cbc").available(exception_flag=False),
+        "LP solver CBC not available",
+    )
+    def test_factor_model_parameter_bounds_correct(self):
+        """
+        If LP solver is available, test parameter bounds method
+        for factor model set is correct (check against
+        results from an LP solver).
+        """
+        def eval_parameter_bounds(uncertainty_set, solver):
+            """
+            Evaluate parameter bounds of uncertainty set by solving
+            bounding problems (as opposed to via the `parameter_bounds`
+            method).
+            """
+            bounding_mdl = uncertainty_set.bounding_model()
+
+            param_bounds = []
+            for idx, obj in bounding_mdl.param_var_objectives.items():
+                # activate objective for corresponding dimension
+                obj.activate()
+                bounds = []
+
+                # solve for lower bound, then upper bound
+                # solve should be successful
+                for sense in (minimize, maximize):
+                    obj.sense = sense
+                    solver.solve(bounding_mdl)
+                    bounds.append(value(obj))
+
+                # add parameter bounds for current dimension
+                param_bounds.append(tuple(bounds))
+
+                # ensure sense is minimize when done, deactivate
+                obj.sense = minimize
+                obj.deactivate()
+
+            return param_bounds
+
+        solver = SolverFactory("cbc")
+
+        # four cases where prior parameter bounds
+        # approximations were probably too tight
+        fset1 = FactorModelSet(
+            origin=[0, 0],
+            number_of_factors=3,
+            psi_mat=[[1, -1, 1], [1, 0.1, 1]],
+            beta=1/6,
+        )
+        fset2 = FactorModelSet(
+            origin=[0],
+            number_of_factors=3,
+            psi_mat=[[1, 6, 8]],
+            beta=1/2,
+        )
+        fset3 = FactorModelSet(
+            origin=[1],
+            number_of_factors=2,
+            psi_mat=[[1, 2]],
+            beta=1/4,
+        )
+        fset4 = FactorModelSet(
+            origin=[1],
+            number_of_factors=3,
+            psi_mat=[[-1, -6, -8]],
+            beta=1/2,
+        )
+
+        # check parameter bounds matches LP results
+        # exactly for each case
+        for fset in [fset1, fset2, fset3, fset4]:
+            param_bounds = fset.parameter_bounds
+            lp_param_bounds = eval_parameter_bounds(fset, solver)
+
+            self.assertTrue(
+                np.allclose(param_bounds, lp_param_bounds),
+                msg=(
+                    "Parameter bounds not consistent with LP values for "
+                    "FactorModelSet with parameterization:\n"
+                    f"F={fset.number_of_factors},\n"
+                    f"beta={fset.beta},\n"
+                    f"psi_mat={fset.psi_mat},\n"
+                    f"origin={fset.origin}."
+                ),
+            )
+
     @unittest.skipIf(not numpy_available, 'Numpy is not available.')
     def test_uncertainty_set_with_correct_params(self):
         '''

pyomo/contrib/pyros/uncertainty_sets.py

@@ -348,6 +348,38 @@ def parameter_bounds(self):
         """
         raise NotImplementedError
 
+    def bounding_model(self):
+        """
+        Make uncertain parameter value bounding problems (optimize
+        value of each uncertain parameter subject to constraints on the
+        uncertain parameters).
+
+        Returns
+        -------
+        model : ConcreteModel
+            Bounding problem, with all Objectives deactivated.
+        """
+        model = ConcreteModel()
+        model.util = Block()
+
+        # construct param vars, initialize to nominal point
+        model.param_vars = Var(range(self.dim))
+
+        # add constraints
+        model.cons = self.set_as_constraint(
+            uncertain_params=model.param_vars,
+            model=model,
+        )
+
+        @model.Objective(range(self.dim))
+        def param_var_objectives(self, idx):
+            return model.param_vars[idx]
+
+        # deactivate all objectives
+        model.param_var_objectives.deactivate()
+
+        return model
+
     def is_bounded(self, config):
         """
         Determine whether the uncertainty set is bounded.
@@ -374,35 +406,36 @@ def is_bounded(self, config):
         This method is invoked during the validation step of a PyROS
         solver call.
         """
-        # === Determine bounds on all uncertain params
-        bounding_model = ConcreteModel()
-        bounding_model.util = Block() # So that boundedness checks work for Cardinality and FactorModel sets
-        bounding_model.uncertain_param_vars = IndexedVar(range(len(config.uncertain_params)), initialize=1)
-        for idx, param in enumerate(config.uncertain_params):
-            bounding_model.uncertain_param_vars[idx].set_value(
-                param.value, skip_validation=True)
-
-        bounding_model.add_component("uncertainty_set_constraint",
-                                     config.uncertainty_set.set_as_constraint(
-                                         uncertain_params=bounding_model.uncertain_param_vars,
-                                         model=bounding_model,
-                                         config=config
-                                     ))
-
-        for idx, param in enumerate(list(bounding_model.uncertain_param_vars.values())):
-            bounding_model.add_component("lb_obj_" + str(idx), Objective(expr=param, sense=minimize))
-            bounding_model.add_component("ub_obj_" + str(idx), Objective(expr=param, sense=maximize))
-
-        for o in bounding_model.component_data_objects(Objective):
-            o.deactivate()
-
-        for i in range(len(bounding_model.uncertain_param_vars)):
-            for limit in ("lb", "ub"):
-                getattr(bounding_model, limit + "_obj_" + str(i)).activate()
-                res = config.global_solver.solve(bounding_model, tee=False)
-                getattr(bounding_model, limit + "_obj_" + str(i)).deactivate()
+        bounding_model = self.bounding_model()
+        solver = config.global_solver
+
+        # initialize uncertain parameter variables
+        for param, param_var in zip(
+                config.uncertain_params,
+                bounding_model.param_vars.values(),
+                ):
+            param_var.set_value(param.value, skip_validation=True)
+
+        for idx, obj in bounding_model.param_var_objectives.items():
+            # activate objective for corresponding dimension
+            obj.activate()
+
+            # solve for lower bound, then upper bound
+            for sense in (minimize, maximize):
+                obj.sense = sense
+                res = solver.solve(
+                    bounding_model,
+                    load_solutions=False,
+                    tee=False,
+                )
+
                 if not check_optimal_termination(res):
                     return False
+
+            # ensure sense is minimize when done, deactivate
+            obj.sense = minimize
+            obj.deactivate()
+
         return True
 
     def is_nonempty(self, config):
@@ -1572,10 +1605,9 @@ class FactorModelSet(UncertaintySet):
         Natural number representing the dimensionality of the
         space to which the set projects.
     psi_mat : (N, `number_of_factors`) array_like
-        Matrix with nonnegative entries designating each
-        uncertain parameter's contribution to each factor.
-        Each row is associated with a separate uncertain parameter.
-        Each column is associated with a separate factor.
+        Matrix designating each uncertain parameter's contribution to
+        each factor.  Each row is associated with a separate uncertain
+        parameter.  Each column is associated with a separate factor.
     beta : numeric type
         Real value between 0 and 1 specifying the fraction of the
         independent factors that can simultaneously attain
@@ -1661,8 +1693,7 @@ def psi_mat(self):
         (N, `number_of_factors`) numpy.ndarray : Matrix designating each
         uncertain parameter's contribution to each factor. Each row is
         associated with a separate uncertain parameter. Each column with
-        a separate factor.  Every entry of the matrix must be
-        nonnegative.
+        a separate factor.
         """
         return self._psi_mat
 
@@ -1697,13 +1728,6 @@ def psi_mat(self, val):
                     "one nonzero entry"
                 )
 
-            for entry in column:
-                if entry < 0:
-                    raise ValueError(
-                        f"Entry {entry} of attribute 'psi_mat' is negative. "
-                        "Ensure all entries are nonnegative"
-                    )
-
         self._psi_mat = psi_mat_arr
 
     @property
@@ -1758,27 +1782,51 @@ def parameter_bounds(self):
             the uncertain parameter bounds for the corresponding set
             dimension.
         """
-        nom_val = self.origin
+        F = self.number_of_factors
         psi_mat = self.psi_mat
 
-        F = self.number_of_factors
-        beta_F = self.beta * F
-        floor_beta_F = math.floor(beta_F)
+        # evaluate some important quantities
+        beta_F = self.beta * self.number_of_factors
+        crit_pt_type = int((beta_F + F) / 2)
+        beta_F_fill_in = (beta_F + F) - 2 * crit_pt_type - 1
+
+        # argsort rows of psi_mat in descending order
+        row_wise_args = np.argsort(-psi_mat, axis=1)
+
         parameter_bounds = []
-        for i in range(len(nom_val)):
-            non_decreasing_factor_row = sorted(psi_mat[i], reverse=True)
-            # deviation = sum_j=1^floor(beta F) {psi_if_j} + (beta F - floor(beta F)) psi_{if_{betaF +1}}
-            # because indexing starts at 0, we adjust the limit on the sum and the final factor contribution
-            if beta_F - floor_beta_F == 0:
-                deviation = sum(non_decreasing_factor_row[j] for j in range(floor_beta_F - 1))
+        for idx, orig_val in enumerate(self.origin):
+            # number nonnegative values in row
+            M = len(psi_mat[idx][psi_mat[idx] >= 0])
+
+            # argsort psi matrix row in descending order
+            sorted_psi_row_args = row_wise_args[idx]
+            sorted_psi_row = psi_mat[idx, sorted_psi_row_args]
+
+            # now evaluate max deviation from origin
+            # (depends on number nonneg entries and critical point type)
+            if M > crit_pt_type:
+                max_deviation = (
+                    sorted_psi_row[:crit_pt_type].sum()
+                    + beta_F_fill_in * sorted_psi_row[crit_pt_type]
+                    - sorted_psi_row[crit_pt_type + 1:].sum()
+                )
+            elif M < F - crit_pt_type:
+                max_deviation = (
+                    sorted_psi_row[:F - crit_pt_type - 1].sum()
+                    - beta_F_fill_in * sorted_psi_row[F - crit_pt_type - 1]
+                    - sorted_psi_row[F - crit_pt_type:].sum()
+                )
             else:
-                deviation = sum(non_decreasing_factor_row[j] for j in range(floor_beta_F - 1)) + (
-                            beta_F - floor_beta_F) * psi_mat[i][floor_beta_F]
-            lb = nom_val[i] - deviation
-            ub = nom_val[i] + deviation
-            if lb > ub:
-                raise AttributeError("The computed lower bound on uncertain parameters must be less than or equal to the upper bound.")
-            parameter_bounds.append((lb, ub))
+                max_deviation = (
+                    sorted_psi_row[:M].sum()
+                    - sorted_psi_row[M:].sum()
+                )
+
+            # finally, evaluate the bounds for this dimension
+            parameter_bounds.append(
+                (orig_val - max_deviation, orig_val + max_deviation),
+            )
+
         return parameter_bounds
 
     def set_as_constraint(self, uncertain_params, **kwargs):