Merge pull request #44 from michaelbynum/bounds_push

Relax Bounds

parapint/algorithms/interior_point.py

@@ -8,7 +8,7 @@
 from parapint.interfaces.interface import BaseInteriorPointInterface
 from parapint.linalg.base_linear_solver_interface import LinearSolverInterface
 from typing import Optional
-from pyomo.common.config import ConfigDict, ConfigValue, PositiveFloat, NonNegativeInt
+from pyomo.common.config import ConfigDict, ConfigValue, PositiveFloat, NonNegativeInt, NonNegativeFloat
 
 
 """
@@ -54,7 +54,7 @@ def __init__(self,
         self.declare('factor_decrease', ConfigValue(domain=PositiveFloat))
         self.declare('max_coef', ConfigValue(domain=PositiveFloat))
 
-        self.init_coef = 1e-4
+        self.init_coef = 1e-8
         self.factor_increase = 10
         self.factor_decrease = 1/3
         self.max_coef = 1e9
@@ -154,11 +154,12 @@ def __init__(self,
         self.declare('line_search', LineSearchOptions())
         self.declare('unified_step', ConfigValue(domain=bool))
         self.declare('error_scaling', ConfigValue(domain=PositiveFloat))
+        self.declare('bounds_relaxation_factor', ConfigValue(domain=NonNegativeFloat))
 
-        self.max_iter = 100
+        self.max_iter = 1000
         self.tol = 1e-8
         self.init_barrier_parameter = 0.1
-        self.minimum_barrier_parameter = 1e-12
+        self.minimum_barrier_parameter = 1e-9
         self.barrier_decrease = 10
         self.report_timing = False
         self.use_inertia_correction = True
@@ -167,6 +168,7 @@ def __init__(self,
         self.line_search: LineSearchOptions = LineSearchOptions()
         self.unified_step: bool = False
         self.error_scaling: float = 100
+        self.bounds_relaxation_factor: float = 1e-8
 
 
 def check_convergence(interface, barrier, error_scaling: float, timer=None):
@@ -421,6 +423,8 @@ def ip_solve(interface: BaseInteriorPointInterface,
     timer.start('IP solve')
     timer.start('init')
 
+    interface.set_bounds_relaxation_factor(options.bounds_relaxation_factor)
+
     barrier_parameter = options.init_barrier_parameter
     inertia_coef = options.inertia_correction.init_coef
     used_inertia_coef = 0

parapint/algorithms/tests/test_reg.py

@@ -58,7 +58,7 @@ def _test_regularization(self, linear_solver):
                                          inertia_coef=options.inertia_correction.init_coef)
 
         # Expected regularization coefficient:
-        self.assertAlmostEqual(reg_coef, 1e-4)
+        self.assertGreaterEqual(reg_coef, 1e-8)
 
         desired_n_neg_evals = (interface.n_eq_constraints() +
                                interface.n_ineq_constraints())

parapint/examples/tests/test_examples.py

@@ -28,9 +28,9 @@ def test_stochastic(self):
         interface = examples.stochastic.main(farmer=farmer,
                                              subproblem_solver_class=parapint.linalg.ScipyInterface,
                                              subproblem_solver_options={'compute_inertia': True})
-        self.assertAlmostEqual(interface.pyomo_model(farmer.scenarios[rank]).devoted_acreage['CORN'].value, 80)
-        self.assertAlmostEqual(interface.pyomo_model(farmer.scenarios[rank]).devoted_acreage['SUGAR_BEETS'].value, 250)
-        self.assertAlmostEqual(interface.pyomo_model(farmer.scenarios[rank]).devoted_acreage['WHEAT'].value, 170)
+        self.assertAlmostEqual(interface.pyomo_model(farmer.scenarios[rank]).devoted_acreage['CORN'].value, 80, 5)
+        self.assertAlmostEqual(interface.pyomo_model(farmer.scenarios[rank]).devoted_acreage['SUGAR_BEETS'].value, 250, 5)
+        self.assertAlmostEqual(interface.pyomo_model(farmer.scenarios[rank]).devoted_acreage['WHEAT'].value, 170, 5)
 
     @pytest.mark.parallel
     @pytest.mark.medium

parapint/interfaces/interface.py

@@ -12,6 +12,14 @@ class BaseInteriorPointInterface(ABC):
     responsible for function evaluations and for building the KKT system (matrix and rhs).
     """
 
+    @abstractmethod
+    def get_bounds_relaxation_factor(self) -> float:
+        pass
+
+    @abstractmethod
+    def set_bounds_relaxation_factor(self, val: float):
+        pass
+
     @abstractmethod
     def n_primals(self):
         pass
@@ -246,14 +254,17 @@ def __init__(self, pyomo_model):
             self._nlp = ampl_nlp.AmplNLP(pyomo_model)
         else:
             self._nlp = pyomo_nlp.PyomoNLP(pyomo_model, nl_file_options={'skip_trivial_constraints': True})
+
+        self.bounds_relaxation_factor = 0
+
         self._slacks = self.init_slacks()
 
         # set the init_duals_primals_lb/ub from ipopt_zL_out, ipopt_zU_out if available
         # need to compress them as well and initialize the duals_primals_lb/ub 
         self._init_duals_primals_lb, self._init_duals_primals_ub =\
             self._get_full_duals_primals_bounds()
-        self._init_duals_primals_lb[np.isneginf(self._nlp.primals_lb())] = 0
-        self._init_duals_primals_ub[np.isinf(self._nlp.primals_ub())] = 0
+        self._init_duals_primals_lb[np.isneginf(self.primals_lb())] = 0
+        self._init_duals_primals_ub[np.isinf(self.primals_ub())] = 0
         self._duals_primals_lb = self._init_duals_primals_lb.copy()
         self._duals_primals_ub = self._init_duals_primals_ub.copy()
 
@@ -276,6 +287,12 @@ def __init__(self, pyomo_model):
         self._delta_duals_ineq = None
         self._barrier = None
 
+    def get_bounds_relaxation_factor(self) -> float:
+        return self.bounds_relaxation_factor
+
+    def set_bounds_relaxation_factor(self, val: float):
+        self.bounds_relaxation_factor = val
+
     def n_primals(self):
         return self._nlp.n_primals()
 
@@ -375,16 +392,36 @@ def get_duals_slacks_ub(self):
         return self._duals_slacks_ub
 
     def primals_lb(self):
-        return self._nlp.primals_lb()
+        lbs = self._nlp.primals_lb()
+        if self.bounds_relaxation_factor == 0:
+            return lbs
+        eye = np.ones(lbs.size)
+        lbs_mod = lbs - self.bounds_relaxation_factor * np.max(np.array([eye, np.abs(lbs)]), axis=0)
+        return lbs_mod
 
     def primals_ub(self):
-        return self._nlp.primals_ub()
+        ubs = self._nlp.primals_ub()
+        if self.bounds_relaxation_factor == 0:
+            return ubs
+        eye = np.ones(ubs.size)
+        ubs_mod = ubs + self.bounds_relaxation_factor * np.max(np.array([eye, np.abs(ubs)]), axis=0)
+        return ubs_mod
 
     def ineq_lb(self):
-        return self._nlp.ineq_lb()
+        lbs = self._nlp.ineq_lb()
+        if self.bounds_relaxation_factor == 0:
+            return lbs
+        eye = np.ones(lbs.size)
+        lbs_mod = lbs - self.bounds_relaxation_factor * np.max(np.array([eye, np.abs(lbs)]), axis=0)
+        return lbs_mod
 
     def ineq_ub(self):
-        return self._nlp.ineq_ub()
+        ubs = self._nlp.ineq_ub()
+        if self.bounds_relaxation_factor == 0:
+            return ubs
+        eye = np.ones(ubs.size)
+        ubs_mod = ubs + self.bounds_relaxation_factor * np.max(np.array([eye, np.abs(ubs)]), axis=0)
+        return ubs_mod
 
     def set_barrier_parameter(self, barrier):
         self._barrier = barrier
@@ -410,8 +447,8 @@ def evaluate_primal_dual_kkt_matrix(self, timer=None):
         primals = self._nlp.get_primals()
 
         timer.start('hess block')
-        data = (duals_primals_lb/(primals - self._nlp.primals_lb()) +
-                duals_primals_ub/(self._nlp.primals_ub() - primals))
+        data = (duals_primals_lb/(primals - self.primals_lb()) +
+                duals_primals_ub/(self.primals_ub() - primals))
         n = self._nlp.n_primals()
         indices = np.arange(n)
         hess_block.row = np.concatenate([hess_block.row, indices])
@@ -420,8 +457,8 @@ def evaluate_primal_dual_kkt_matrix(self, timer=None):
         timer.stop('hess block')
 
         timer.start('slack block')
-        data = (duals_slacks_lb/(self._slacks - self._nlp.ineq_lb()) +
-                duals_slacks_ub/(self._nlp.ineq_ub() - self._slacks))
+        data = (duals_slacks_lb/(self._slacks - self.ineq_lb()) +
+                duals_slacks_ub/(self.ineq_ub() - self._slacks))
         n = self._nlp.n_ineq_constraints()
         indices = np.arange(n)
         slack_block = scipy.sparse.coo_matrix((data, (indices, indices)), shape=(n, n))
@@ -472,14 +509,14 @@ def evaluate_primal_dual_kkt_rhs(self, timer=None):
         grad_lag_primals = (grad_obj +
                             jac_eq.transpose() * self._nlp.get_duals_eq() +
                             jac_ineq.transpose() * self._nlp.get_duals_ineq() -
-                            self._barrier / (self._nlp.get_primals() - self._nlp.primals_lb()) +
-                            self._barrier / (self._nlp.primals_ub() - self._nlp.get_primals()))
+                            self._barrier / (self._nlp.get_primals() - self.primals_lb()) +
+                            self._barrier / (self.primals_ub() - self._nlp.get_primals()))
         timer.stop('grad_lag_primals')
 
         timer.start('grad_lag_slacks')
         grad_lag_slacks = (-self._nlp.get_duals_ineq() -
-                           self._barrier / (self._slacks - self._nlp.ineq_lb()) +
-                           self._barrier / (self._nlp.ineq_ub() - self._slacks))
+                           self._barrier / (self._slacks - self.ineq_lb()) +
+                           self._barrier / (self.ineq_ub() - self._slacks))
         timer.stop('grad_lag_slacks')
 
         rhs = BlockVector(4)
@@ -510,25 +547,25 @@ def get_delta_duals_ineq(self):
 
     def get_delta_duals_primals_lb(self):
         res = (((self._barrier - self._duals_primals_lb * self._delta_primals) /
-                (self._nlp.get_primals() - self._nlp.primals_lb())) -
+                (self._nlp.get_primals() - self.primals_lb())) -
                self._duals_primals_lb)
         return res
 
     def get_delta_duals_primals_ub(self):
         res = (((self._barrier + self._duals_primals_ub * self._delta_primals) /
-                (self._nlp.primals_ub() - self._nlp.get_primals())) -
+                (self.primals_ub() - self._nlp.get_primals())) -
                self._duals_primals_ub)
         return res
 
     def get_delta_duals_slacks_lb(self):
         res = (((self._barrier - self._duals_slacks_lb * self._delta_slacks) /
-                (self._slacks - self._nlp.ineq_lb())) -
+                (self._slacks - self.ineq_lb())) -
                self._duals_slacks_lb)
         return res
 
     def get_delta_duals_slacks_ub(self):
         res = (((self._barrier + self._duals_slacks_ub * self._delta_slacks) /
-                (self._nlp.ineq_ub() - self._slacks)) -
+                (self.ineq_ub() - self._slacks)) -
                self._duals_slacks_ub)
         return res
 

parapint/interfaces/schur_complement/mpi_sc_ip_interface.py

@@ -135,6 +135,9 @@ def __init__(self, start_t: float, end_t: float, num_time_blocks: int, comm: MPI
         self._setup_jacs()
         self._setup_kkt_and_rhs_structure()
 
+        self._bounds_relaxation_factor = 0
+        self.set_bounds_relaxation_factor(self._bounds_relaxation_factor)
+
     def _build_mpi_block_matrix(self, extra_row: bool = False, extra_col: bool = False) -> MPIBlockMatrix:
         if extra_row:
             nrows = self._num_time_blocks + 1
@@ -383,6 +386,9 @@ def __init__(self,
         self._setup_jacs()
         self._setup_kkt_and_rhs_structure()
 
+        self._bounds_relaxation_factor = 0
+        self.set_bounds_relaxation_factor(self._bounds_relaxation_factor)
+
     def _build_mpi_block_matrix(self, extra_row: bool = False, extra_col: bool = False) -> MPIBlockMatrix:
         if extra_row:
             nrows = self._num_scenarios + 1

parapint/interfaces/schur_complement/sc_ip_interface.py

@@ -102,6 +102,9 @@ def __init__(self, start_t: float, end_t: float, num_time_blocks: int):
         self._setup_jacs()
         self._setup_kkt_and_rhs_structure()
 
+        self._bounds_relaxation_factor = 0
+        self.set_bounds_relaxation_factor(self._bounds_relaxation_factor)
+
     @abstractmethod
     def build_model_for_time_block(self,
                                    ndx: int,
@@ -487,6 +490,14 @@ def n_primals(self) -> int:
     def nnz_hessian_lag(self) -> int:
         raise NotImplementedError('This is not done yet')
 
+    def get_bounds_relaxation_factor(self) -> float:
+        return self._bounds_relaxation_factor
+
+    def set_bounds_relaxation_factor(self, val: float):
+        self._bounds_relaxation_factor = val
+        for nlp in self._nlps.values():
+            nlp.set_bounds_relaxation_factor(val)
+
     def primals_lb(self) -> BlockVector:
         """
         Returns
@@ -495,6 +506,8 @@ def primals_lb(self) -> BlockVector:
             The lower bounds for each primal variable. This BlockVector has one block for every time block
             and one block for the coupling variables.
         """
+        for ndx, nlp in self._nlps.items():
+            self._primals_lb.set_block(ndx, nlp.primals_lb())
         return self._primals_lb
 
     def primals_ub(self) -> BlockVector:
@@ -505,6 +518,8 @@ def primals_ub(self) -> BlockVector:
             The upper bounds for each primal variable. This BlockVector has one block for every time block
             and one block for the coupling variables.
         """
+        for ndx, nlp in self._nlps.items():
+            self._primals_ub.set_block(ndx, nlp.primals_ub())
         return self._primals_ub
 
     def init_primals(self) -> BlockVector:
@@ -601,6 +616,8 @@ def ineq_lb(self) -> BlockVector:
         ineq_lb: BlockVector
             The lower bounds for each inequality constraint. This BlockVector has one block for every time block.
         """
+        for ndx, nlp in self._nlps.items():
+            self._ineq_lb.set_block(ndx, nlp.ineq_lb())
         return self._ineq_lb
 
     def ineq_ub(self) -> BlockVector:
@@ -610,6 +627,8 @@ def ineq_ub(self) -> BlockVector:
         ineq_lb: BlockVector
             The lower bounds for each inequality constraint. This BlockVector has one block for every time block.
         """
+        for ndx, nlp in self._nlps.items():
+            self._ineq_ub.set_block(ndx, nlp.ineq_ub())
         return self._ineq_ub
 
     def init_duals_eq(self) -> BlockVector:
@@ -1097,6 +1116,9 @@ def __init__(self, scenarios: Sequence, nonanticipative_var_identifiers: Sequenc
         self._setup_jacs()
         self._setup_kkt_and_rhs_structure()
 
+        self._bounds_relaxation_factor = 0
+        self.set_bounds_relaxation_factor(self._bounds_relaxation_factor)
+
     @abstractmethod
     def build_model_for_scenario(self,
                                  scenario_identifier: Any) -> Tuple[_BlockData, Dict[Any, _GeneralVarData]]:
@@ -1310,6 +1332,14 @@ def n_primals(self) -> int:
     def nnz_hessian_lag(self) -> int:
         raise NotImplementedError('This is not done yet')
 
+    def get_bounds_relaxation_factor(self) -> float:
+        return self._bounds_relaxation_factor
+
+    def set_bounds_relaxation_factor(self, val: float):
+        self._bounds_relaxation_factor = val
+        for nlp in self._nlps.values():
+            nlp.set_bounds_relaxation_factor(val)
+
     def primals_lb(self) -> BlockVector:
         """
         Returns
@@ -1318,6 +1348,8 @@ def primals_lb(self) -> BlockVector:
             The lower bounds for each primal variable. This BlockVector has one block for every time block
             and one block for the coupling variables.
         """
+        for ndx, nlp in self._nlps.items():
+            self._primals_lb.set_block(ndx, nlp.primals_lb())
         return self._primals_lb
 
     def primals_ub(self) -> BlockVector:
@@ -1328,6 +1360,8 @@ def primals_ub(self) -> BlockVector:
             The upper bounds for each primal variable. This BlockVector has one block for every time block
             and one block for the coupling variables.
         """
+        for ndx, nlp in self._nlps.items():
+            self._primals_ub.set_block(ndx, nlp.primals_ub())
         return self._primals_ub
 
     def init_primals(self) -> BlockVector:
@@ -1424,6 +1458,8 @@ def ineq_lb(self) -> BlockVector:
         ineq_lb: BlockVector
             The lower bounds for each inequality constraint. This BlockVector has one block for every time block.
         """
+        for ndx, nlp in self._nlps.items():
+            self._ineq_lb.set_block(ndx, nlp.ineq_lb())
         return self._ineq_lb
 
     def ineq_ub(self) -> BlockVector:
@@ -1433,6 +1469,8 @@ def ineq_ub(self) -> BlockVector:
         ineq_lb: BlockVector
             The lower bounds for each inequality constraint. This BlockVector has one block for every time block.
         """
+        for ndx, nlp in self._nlps.items():
+            self._ineq_ub.set_block(ndx, nlp.ineq_ub())
         return self._ineq_ub
 
     def init_duals_eq(self) -> BlockVector: