Merge pull request #467 from bknueven/prox_approx_6

mpisppy/extensions/sep_rho.py

@@ -140,7 +140,10 @@ def _compute_xmin(ph):
 
     def nonant_cost_coeffs(self, s):
         if s not in self._nonant_cost_coeffs:
-            self._nonant_cost_coeffs[s] = nonant_cost_coeffs(s)
+            if hasattr(s._mpisppy_data, "nonant_cost_coeffs"):
+                self._nonant_cost_coeffs[s] = s._mpisppy_model.nonant_cost_coeffs
+            else:
+                self._nonant_cost_coeffs[s] = nonant_cost_coeffs(s)
         return self._nonant_cost_coeffs[s]
 
     def _compute_and_update_rho(self):

mpisppy/phbase.py

@@ -728,6 +728,10 @@ def attach_PH_to_objective(self, add_duals, add_prox, add_smooth=0):
                 scenario._mpisppy_model.xsqvar = pyo.Var(scenario._mpisppy_data.nonant_indices, dense=False, bounds=(0, None))
                 scenario._mpisppy_model.xsqvar_cuts = pyo.Constraint(scenario._mpisppy_data.nonant_indices, pyo.Integers)
                 scenario._mpisppy_data.xsqvar_prox_approx = {}
+                try:
+                    scenario._mpisppy_data.nonant_cost_coeffs = sputils.nonant_cost_coeffs(scenario)
+                except sputils.NonLinearProblemFound:
+                    raise RuntimeError("The proximal term approximation can only be used with a linear objective function")
             else:
                 scenario._mpisppy_model.xsqvar = None
                 scenario._mpisppy_data.xsqvar_prox_approx = False
@@ -753,7 +757,7 @@ def attach_PH_to_objective(self, add_duals, add_prox, add_smooth=0):
                     elif self._prox_approx:
                         xvarsqrd = scenario._mpisppy_model.xsqvar[ndn_i]
                         scenario._mpisppy_data.xsqvar_prox_approx[ndn_i] = \
-                                ProxApproxManager(scenario._mpisppy_model, xvar, ndn_i)
+                                ProxApproxManager(scenario, xvar, ndn_i)
                     else:
                         xvarsqrd = xvar**2
                     prox_expr += (scenario._mpisppy_model.rho[ndn_i] / 2.0) * \

mpisppy/utils/prox_approx.py

@@ -34,13 +34,14 @@ class _ProxApproxManager:
     '''
     __slots__ = ()
 
-    def __init__(self, mpisppy_model, xvar, ndn_i):
+    def __init__(self, scenario, xvar, ndn_i):
         self.xvar = xvar
-        self.xvarsqrd = mpisppy_model.xsqvar[ndn_i]
-        self.cuts = mpisppy_model.xsqvar_cuts
-        self.xbar = mpisppy_model.xbars[ndn_i]
-        self.rho = mpisppy_model.rho[ndn_i]
-        self.W = mpisppy_model.W[ndn_i]
+        self.xvarsqrd = scenario._mpisppy_model.xsqvar[ndn_i]
+        self.cuts = scenario._mpisppy_model.xsqvar_cuts
+        self.xbar = scenario._mpisppy_model.xbars[ndn_i]
+        self.rho = scenario._mpisppy_model.rho[ndn_i]
+        self.W = scenario._mpisppy_model.W[ndn_i]
+        self.c = scenario._mpisppy_data.nonant_cost_coeffs[ndn_i]
         self.var_index = ndn_i
         self.cut_index = 0
         self.cut_values = array("d")
@@ -88,6 +89,68 @@ def check_and_add_value(self, val, tolerance):
         self.cut_values.insert(idx, val)
         return True
 
+    def add_initial_cuts_var(self, x_val, other_bound, sign, tolerance, persistent_solver):
+        if other_bound is not None:
+            min_approx_error_pnt = sign * (2.0 / 3.0) * (other_bound - x_val)
+        else:
+            min_approx_error_pnt = None
+        if self.c > 0:
+            cost_pnt = self.c / self.rho.value
+        elif self.c < 0:
+            cost_pnt = -self.c / self.rho.value
+        else: # c == 0
+            cost_pnt = None
+        num_cuts = 0
+        if cost_pnt is None and min_approx_error_pnt is None:
+            # take a wild guess
+            step = max(1, tolerance)
+            num_cuts += self.add_cut(x_val + sign*step, tolerance, persistent_solver)
+            num_cuts += self.add_cut(x_val + sign*10*step, tolerance, persistent_solver)
+            return num_cuts
+        # no upper bound, or weak upper bound compared to cost_pnt
+        if (min_approx_error_pnt is None) or (cost_pnt is not None and min_approx_error_pnt > 10*cost_pnt):
+            step = max(cost_pnt, tolerance)
+            num_cuts += self.add_cut(x_val + sign*step, tolerance, persistent_solver)
+            num_cuts += self.add_cut(x_val + sign*10*step, tolerance, persistent_solver)
+            return num_cuts
+        # no objective, or the objective is "large" compared to min_approx_error_pnt
+        if (cost_pnt is None) or (cost_pnt >= min_approx_error_pnt / 2.0):
+            step = max(min_approx_error_pnt, tolerance)
+            num_cuts += self.add_cut(x_val + step, tolerance, persistent_solver)
+            # guard against horrible bounds
+            if min_approx_error_pnt / 2.0 > max(1, tolerance):
+                step = max(1, tolerance)
+                num_cuts += self.add_cut(x_val + sign*step, tolerance, persistent_solver)
+            else:
+                step = max(min_approx_error_pnt/2.0, tolerance)
+                num_cuts += self.add_cut(x_val + sign*step, tolerance, persistent_solver)
+            return num_cuts
+        # cost_pnt and min_approx_error_pnt are *not* None
+        # and (cost_pnt < min_approx_error_pnt / 2.0 <= min_approx_error_pnt)
+        step = max(cost_pnt, tolerance)
+        num_cuts += self.add_cut(x_val + sign*step, tolerance, persistent_solver)
+        step = max(min_approx_error_pnt, tolerance)
+        num_cuts += self.add_cut(x_val + sign*step, tolerance, persistent_solver)
+        return num_cuts
+
+
+    def add_initial_cuts_var_right(self, x_val, bounds, tolerance, persistent_solver):
+        return self.add_initial_cuts_var(x_val, bounds[1], 1, tolerance, persistent_solver)
+
+    def add_initial_cuts_var_left(self, x_val, bounds, tolerance, persistent_solver):
+        return self.add_initial_cuts_var(x_val, bounds[0], -1, tolerance, persistent_solver)
+
+    def add_initial_cuts(self, x_val, tolerance, persistent_solver):
+        num_cuts = 0
+        bounds = self.xvar.bounds
+        # no lower bound **or** sufficiently inside lower bound
+        if bounds[0] is None or (x_val > bounds[0] + 3*tolerance):
+            num_cuts += self.add_initial_cuts_var_left(x_val, bounds, tolerance, persistent_solver)
+        # no upper bound **or** sufficiently inside upper bound
+        if bounds[1] is None or (x_val < bounds[1] - 3*tolerance):
+            num_cuts += self.add_initial_cuts_var_right(x_val, bounds, tolerance, persistent_solver)
+        return num_cuts
+
     def add_cuts(self, x_val, tolerance, persistent_solver):
         x_bar = self.xbar.value
         rho = self.rho.value
@@ -112,17 +175,7 @@ def add_cuts(self, x_val, tolerance, persistent_solver):
         # to capture something of the proximal term. This can happen
         # when x_bar == x_val and W == 0.
         if self.cut_index <= 1:
-            lb, ub = self.xvar.bounds
-            upval = 1
-            dnval = 1
-            if ub is not None:
-                # after a lot of calculus, you can show that
-                # this point minimizes the error in the approximation
-                upval = 2.0 * (ub - x_val) / 3.0
-            if lb is not None:
-                dnval = 2.0 * (x_val - lb) / 3.0
-            num_cuts += self.add_cut(x_val + max(upval, tolerance+1e-06), tolerance, persistent_solver)
-            num_cuts += self.add_cut(x_val - max(dnval, tolerance+1e-06), tolerance, persistent_solver)
+            num_cuts += self.add_initial_cuts(x_val, tolerance, persistent_solver)
         # print(f"{x_val=}, {x_bar=}, {W=}")
         # print(f"{self.cut_values=}")
         # print(f"{self.cut_index=}")

mpisppy/utils/sputils.py

@@ -1004,6 +1004,8 @@ def find_objective(pyomomodel, active=False):
     raise RuntimeError("Could not identify exactly one objective for model "
                        f"{pyomomodel.name} (found {len(obj)} objectives)")
 
+class NonLinearProblemFound(RuntimeError):
+    pass
 
 def nonant_cost_coeffs(s):
     """