diff --git a/mpisppy/extensions/sep_rho.py b/mpisppy/extensions/sep_rho.py index 4f8e731c..90b96326 100644 --- a/mpisppy/extensions/sep_rho.py +++ b/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): diff --git a/mpisppy/phbase.py b/mpisppy/phbase.py index 745e4bc8..bdd8a5a0 100644 --- a/mpisppy/phbase.py +++ b/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) * \ diff --git a/mpisppy/utils/prox_approx.py b/mpisppy/utils/prox_approx.py index 851663cb..b78a06a0 100644 --- a/mpisppy/utils/prox_approx.py +++ b/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=}") diff --git a/mpisppy/utils/sputils.py b/mpisppy/utils/sputils.py index 40558d6f..a3e22a88 100644 --- a/mpisppy/utils/sputils.py +++ b/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): """