diff --git a/.github/workflows/ubuntu_python_matrix_test.yml b/.github/workflows/ubuntu_python_matrix_test.yml index 1591597..fbd1286 100644 --- a/.github/workflows/ubuntu_python_matrix_test.yml +++ b/.github/workflows/ubuntu_python_matrix_test.yml @@ -7,6 +7,8 @@ on: pull_request: branches: - master + schedule: + - cron: "0 2 * * 1-5" jobs: linux: @@ -32,10 +34,10 @@ jobs: conda init bash conda env list - - uses: actions/checkout@v1 + - uses: actions/checkout@v3 - name: Set up Python ${{ matrix.python-version }} - uses: actions/setup-python@v2 + uses: actions/setup-python@v4 with: python-version: ${{ matrix.python-version }} @@ -46,11 +48,12 @@ jobs: conda activate test conda config --set always_yes yes python -m pip install --upgrade pip + # scipy >= 1.8 does not support Python 3.7 + pip install nose coverage codecov pyyaml scipy==1.7.3 sphinx sphinx_rtd_theme gurobipy coramin pip install --quiet git+https://github.com/PyUtilib/pyutilib pip install --quiet git+https://github.com/Pyomo/pyomo python setup.py develop pip install wheel - pip install nose coverage codecov pyyaml numpy scipy sphinx sphinx_rtd_theme conda install -c conda-forge coincbc glpk ipopt conda list diff --git a/doc/readthedocs/solvers.rst b/doc/readthedocs/solvers.rst index fc7b8be..c5d0d30 100644 --- a/doc/readthedocs/solvers.rst +++ b/doc/readthedocs/solvers.rst @@ -287,7 +287,8 @@ termination condition indicates that an optimal solution was found. For example >>> nlp = pao.Solver('ipopt', print_level=3) >>> opt = pao.Solver('pao.pyomo.REG', nlp_solver=nlp) - >>> results = opt.solve(M) + >>> results = opt.solve(M) #doctest:+ELLIPSIS + W... >>> print(results.solver.termination_condition) TerminationCondition.optimal >>> results.check_optimal_termination() @@ -329,7 +330,8 @@ Pyomo solver is initially created: >>> nlp = pao.Solver('ipopt', print_level=3) >>> opt = pao.Solver('pao.pyomo.REG', nlp_solver=nlp) - >>> results = opt.solve(M) + >>> results = opt.solve(M) #doctest:+ELLIPSIS + W... When executing locally, the :keyword:`executable` option can be used to explicitly specify the path to the executable that is used by this solver. diff --git a/pao/duality/dual.py b/pao/duality/dual.py new file mode 100644 index 0000000..db8b8cf --- /dev/null +++ b/pao/duality/dual.py @@ -0,0 +1,281 @@ +import pyomo.environ as pe +from pyomo.core.base.block import _BlockData, ScalarBlock +from pyomo.core.base.var import _GeneralVarData, ScalarVar, IndexedVar +from typing import Sequence, Optional, Iterable, List +import math +from pyomo.common.collections import ComponentSet +from pyomo.core.expr.visitor import identify_variables, replace_expressions +from pyomo.core.expr.calculus.diff_with_pyomo import reverse_sd +from pyomo.contrib.appsi.base import PersistentBase +from pyomo.core.base.param import _ParamData +from pyomo.core.base.constraint import _GeneralConstraintData, IndexedConstraint +from .utils import ComponentHasher +from pyomo.core.base.sos import _SOSConstraintData +from pyomo.core.base.objective import _GeneralObjectiveData +from pyomo.repn.standard_repn import generate_standard_repn +from pyomo.contrib.appsi.cmodel import cmodel, cmodel_available + + +class Dual(PersistentBase): + def __init__(self, fixed_vars: Optional[Iterable[_GeneralVarData]] = None): + super().__init__(only_child_vars=False) + self._dual = ScalarBlock(concrete=True) + self._setup_dual() + self._grad_lag_map = pe.ComponentMap() + self._lagrangian_terms = dict() + if fixed_vars is None: + self._fixed_vars = ComponentSet() + else: + self._fixed_vars = ComponentSet(fixed_vars) + self._old_obj = None + self._old_obj_vars = list() + self._linear_primals = ComponentSet() + + def _setup_dual(self): + self._dual.grad_lag_set = pe.Set(dimen=1) + self._dual.eq_dual_set = pe.Set(dimen=1) + self._dual.ineq_dual_set = pe.Set(dimen=1) + + self._dual.eq_duals = IndexedVar(self._dual.eq_dual_set) + self._dual.ineq_duals = IndexedVar(self._dual.ineq_dual_set, bounds=(0, None)) + + self._dual.grad_lag = IndexedConstraint(self._dual.grad_lag_set) + + self._dual.objective = pe.Objective(expr=0, sense=pe.maximize) + + def _set_dual_obj(self): + new_obj = sum(self._lagrangian_terms.values()) + sub_map = dict() + for v in self._linear_primals: + sub_map[id(v)] = 0 + new_obj = replace_expressions(new_obj, substitution_map=sub_map) + self._dual.objective.expr = new_obj + + def set_instance(self, model): + saved_update_config = self.update_config + self.__init__(fixed_vars=self._fixed_vars) + self.update_config = saved_update_config + self._model = model + if self.use_extensions and cmodel_available: + self._expr_types = cmodel.PyomoExprTypes() + self.add_block(model) + if self._objective is None: + self.set_objective(None) + + def dual(self, model: _BlockData): + if model is not self._model: + self._dual = ScalarBlock(concrete=True) + self._setup_dual() + self.set_instance(model) + else: + self.update() + self._set_dual_obj() + return self._dual + + def _update_var_bounds(self, v): + v_lb, v_ub = v.bounds + lb_hasher = ComponentHasher(v, 'lb') + ub_hasher = ComponentHasher(v, 'ub') + + for v_bound, hasher in [(v_lb, lb_hasher), (v_ub, ub_hasher)]: + if v_bound is None: + if hasher in self._grad_lag_map[v]: + self._grad_lag_map[v].pop(hasher) + del self._dual.ineq_duals[hasher] + self._dual.ineq_dual_set.remove(hasher) + self._lagrangian_terms.pop(hasher) + else: + if hasher not in self._grad_lag_map[v]: + self._dual.ineq_dual_set.add(hasher) + if hasher.bound == 'lb': + e = self._dual.ineq_duals[hasher] * (v_lb - v) + else: + e = self._dual.ineq_duals[hasher] * (v - v_ub) + if hasher in self._grad_lag_map[v]: + self._grad_lag_map[v][hasher] = (reverse_sd(e)[v], True) + self._lagrangian_terms[hasher] = e + else: + self._process_lagrangian_term(e, [v], hasher) + + def _add_variables(self, variables: List[_GeneralVarData]): + for v in variables: + if v in self._fixed_vars: + continue + if v.fixed: + continue + self._grad_lag_map[v] = dict() + v_hasher = ComponentHasher(v, None) + self._dual.grad_lag_set.add(v_hasher) + self._dual.grad_lag[v_hasher] = (0, 0) + self._update_var_bounds(v) + self._linear_primals.add(v) + + def _remove_variables(self, variables: List[_GeneralVarData]): + for v in variables: + if v in self._fixed_vars: + continue + + self._linear_primals.discard(v) + + v_hasher = ComponentHasher(v, None) + lb_hasher = ComponentHasher(v, 'lb') + ub_hasher = ComponentHasher(v, 'ub') + + if lb_hasher in self._dual.ineq_dual_set: + del self._dual.ineq_duals[lb_hasher] + self._dual.ineq_dual_set.remove(lb_hasher) + self._lagrangian_terms.pop(lb_hasher) + + if ub_hasher in self._dual.ineq_dual_set: + del self._dual.ineq_duals[ub_hasher] + self._dual.ineq_dual_set.remove(ub_hasher) + self._lagrangian_terms.pop(ub_hasher) + + if v_hasher in self._dual.grad_lag_set: + # the variable may have been removed from these already if it was fixed + # sometime since calling add_variables + del self._dual.grad_lag[v_hasher] + self._dual.grad_lag_set.remove(v_hasher) + self._grad_lag_map.pop(v) + + def _regenerate_grad_lag_for_var(self, v): + v_hasher = ComponentHasher(v, None) + new_body = 0 + self._linear_primals.add(v) + for c, (der, linear) in self._grad_lag_map[v].items(): + new_body += der + if not linear: + self._linear_primals.discard(v) + self._dual.grad_lag[v_hasher] = (new_body, 0) + + def _add_params(self, params: List[_ParamData]): + pass + + def _remove_params(self, params: List[_ParamData]): + pass + + def _process_lagrangian_term(self, expr, variables, c_hasher): + self._lagrangian_terms[c_hasher] = expr + ders = reverse_sd(expr) + + vars_in_expr = identify_variables(expr, include_fixed=False) + unfixed_vars_in_expr = [v for v in vars_in_expr if not v.fixed] + orig_values = [v.value for v in unfixed_vars_in_expr] + unfixed_var_set = ComponentSet(unfixed_vars_in_expr) + + for v in unfixed_vars_in_expr: + v.fix(0) + + for v in variables: + if v in self._fixed_vars: + continue + if v not in unfixed_var_set: + continue + v.unfix() + v_hasher = ComponentHasher(v, None) + orig_body = self._dual.grad_lag[v_hasher].body + new_body = orig_body + ders[v] + self._dual.grad_lag[v_hasher] = (new_body, 0) + repn = generate_standard_repn(expr, quadratic=False) + linear = repn.is_linear() + if not linear: + self._linear_primals.discard(v) + self._grad_lag_map[v][c_hasher] = (ders[v], linear) + v.fix() + + for v, val in zip(unfixed_vars_in_expr, orig_values): + v.unfix() + v.value = val + + def _add_constraints(self, cons: List[_GeneralConstraintData]): + for c in cons: + if c.equality or (c.lb is not None and c.ub is not None and c.lb == c.ub): + c_hasher = ComponentHasher(c, None) + self._dual.eq_dual_set.add(c_hasher) + e = self._dual.eq_duals[c_hasher] * (c.body - c.lb) + self._process_lagrangian_term(e, self._vars_referenced_by_con[c], c_hasher) + else: + if c.lb is not None: + c_hasher = ComponentHasher(c, 'lb') + self._dual.ineq_dual_set.add(c_hasher) + e = self._dual.ineq_duals[c_hasher] * (c.lb - c.body) + self._process_lagrangian_term(e, self._vars_referenced_by_con[c], c_hasher) + if c.ub is not None: + c_hasher = ComponentHasher(c, 'ub') + self._dual.ineq_dual_set.add(c_hasher) + e = self._dual.ineq_duals[c_hasher] * (c.body - c.ub) + self._process_lagrangian_term(e, self._vars_referenced_by_con[c], c_hasher) + + def _removal_helper(self, c_hasher, variables): + if c_hasher in self._dual.eq_dual_set or c_hasher in self._dual.ineq_dual_set: + self._lagrangian_terms.pop(c_hasher) + for v in variables: + if v in self._grad_lag_map: + self._grad_lag_map[v].pop(c_hasher, None) + if c_hasher in self._dual.eq_dual_set: + del self._dual.eq_duals[c_hasher] + self._dual.eq_dual_set.remove(c_hasher) + else: + del self._dual.ineq_duals[c_hasher] + self._dual.ineq_dual_set.remove(c_hasher) + + def _remove_constraints(self, cons: List[_GeneralConstraintData]): + affected_variables = ComponentSet() + for c in cons: + for c_hasher in [ + ComponentHasher(c, None), + ComponentHasher(c, 'lb'), + ComponentHasher(c, 'ub') + ]: + variables = self._vars_referenced_by_con[c] + self._removal_helper(c_hasher, variables) + affected_variables.update(variables) + + for v in affected_variables: + if v in self._grad_lag_map: + self._regenerate_grad_lag_for_var(v) + + def _add_sos_constraints(self, cons: List[_SOSConstraintData]): + if len(cons) > 0: + raise NotImplementedError('Dual does not support SOS constraints') + + def _remove_sos_constraints(self, cons: List[_SOSConstraintData]): + if len(cons) > 0: + raise NotImplementedError('Dual does not support SOS constraints') + + def _set_objective(self, obj: _GeneralObjectiveData): + if obj.sense != pe.minimize: + raise NotImplementedError('Dual does not support maximization problems yet') + if self._old_obj is not None: + old_hasher = ComponentHasher(self._old_obj, None) + self._lagrangian_terms.pop(old_hasher) + for v in self._old_obj_vars: + if v in self._grad_lag_map: + self._grad_lag_map[v].pop(ComponentHasher(self._old_obj, None), None) + self._regenerate_grad_lag_for_var(v) + + hasher = ComponentHasher(obj, None) + e = obj.expr + self._lagrangian_terms[hasher] = e + self._process_lagrangian_term(e, self._vars_referenced_by_obj, hasher) + self._old_obj = obj + self._old_obj_vars = self._vars_referenced_by_obj + + def _update_variables(self, variables: List[_GeneralVarData]): + for v in variables: + if v in self._fixed_vars: + continue + if v.fixed: + self._remove_variables([v]) + else: + if v in self._grad_lag_map: + self._update_var_bounds(v) + self._regenerate_grad_lag_for_var(v) + else: + # the variable was fixed but is now unfixed + # the persistent base class will remove and add any + # constraints that involve this variable + self._add_variables([v]) + + def update_params(self): + pass diff --git a/pao/duality/kkt.py b/pao/duality/kkt.py new file mode 100644 index 0000000..0f00c29 --- /dev/null +++ b/pao/duality/kkt.py @@ -0,0 +1,195 @@ +import pyomo.environ as pe +from pyomo.core.base.block import _BlockData, ScalarBlock +from pyomo.core.base.var import _GeneralVarData, ScalarVar, IndexedVar +from typing import Sequence, Optional, Iterable +import math +from pyomo.common.collections import ComponentSet +from pyomo.core.expr.visitor import identify_variables +from pyomo.core.expr.calculus.diff_with_pyomo import reverse_sd +from .utils import ComponentHasher + + +def convert_integers_to_binaries( + model: _BlockData, + integer_vars: Iterable[_GeneralVarData], +) -> _BlockData: + parent = ScalarBlock(concrete=True) + parent.int_set = pe.Set(dimen=1) + parent.bin_set = pe.Set(dimen=2) + + parent.bins = IndexedVar(parent.bin_set, domain=pe.Binary, bounds=(0, 1)) + parent.int_cons = pe.Constraint(parent.int_set) + + for iv in integer_vars: + if not iv.is_integer(): + raise ValueError(f'{iv} is not an integer variable!') + iv_hasher = ComponentHasher(iv, None) + lb, ub = iv.bounds + iv.domain = pe.Reals + iv.setlb(lb) + iv.setub(ub) + parent.int_set.add(iv_hasher) + n = math.ceil(math.log2(max(abs(lb), abs(ub)) + 1) + 1) + expr = 0 + for ndx in range(n-1): + parent.bin_set.add((iv_hasher, ndx)) + bv = parent.bins[iv_hasher, ndx] + expr += 2**ndx * bv + if lb < 0: + parent.bin_set.add((iv_hasher, n-1)) + bv = parent.bins[iv_hasher, n-1] + expr -= 2**(n-1) * bv + parent.int_cons[iv_hasher] = iv == expr + + parent.block = model + + return parent + + +def convert_binary_domain_to_constraint( + model: _BlockData, + binary_vars: Iterable[_GeneralVarData], +) -> _BlockData: + parent = ScalarBlock(concrete=True) + parent.bin_set = pe.Set(dimen=1) + parent.bin_cons = pe.Constraint(parent.bin_set) + + for bv in binary_vars: + if not bv.is_binary(): + raise ValueError(f'{bv} is not binary variable!') + bv_hasher = ComponentHasher(bv, None) + bv.domain = pe.Reals + bv.setlb(0) + bv.setub(1) + parent.bin_set.add(bv_hasher) + parent.bin_cons[bv_hasher] = bv - bv**2 <= 0 + + parent.block = model + + return parent + + +def _get_vars_in_expr(expr, all_vars, bin_vars, int_vars): + for v in identify_variables(expr, include_fixed=False): + all_vars.add(v) + if v.is_binary(): + bin_vars.add(v) + elif v.is_integer(): + int_vars.add(v) + + +def construct_kkt( + model: _BlockData, + variables: Optional[Sequence[_GeneralVarData]] = None +) -> _BlockData: + all_vars = ComponentSet() + int_vars = ComponentSet() + bin_vars = ComponentSet() + + all_cons = ComponentSet() + obj = None + + for _obj in model.component_data_objects( + pe.Objective, active=True, descend_into=True + ): + if obj is not None: + raise ValueError('found multiple active objectives') + obj = _obj + if variables is None: + _get_vars_in_expr(obj.expr, all_vars, bin_vars, int_vars) + + for con in model.component_data_objects( + pe.Constraint, active=True, descend_into=True + ): + if con.lb is not None or con.ub is not None: + all_cons.add(con) + if variables is None: + _get_vars_in_expr(con.body, all_vars, bin_vars, int_vars) + + if variables is not None: + for v in variables: + all_vars.add(v) + if v.is_binary(): + bin_vars.add(v) + elif v.is_integer(): + int_vars.add(v) + + if len(int_vars) > 0: + model = convert_integers_to_binaries(model, int_vars) + for v in model.bins.values(): + bin_vars.add(v) + all_vars.add(v) + for c in model.int_cons.values(): + all_cons.add(c) + + if len(bin_vars) > 0: + model = convert_binary_domain_to_constraint(model, bin_vars) + for c in model.bin_cons.values(): + all_cons.add(c) + + parent = ScalarBlock(concrete=True) + parent.block = model + parent.eq_dual_set = pe.Set(dimen=1) + parent.ineq_dual_set = pe.Set(dimen=1) + parent.eq_duals = pe.Var(parent.eq_dual_set) + parent.ineq_duals = pe.Var(parent.ineq_dual_set, bounds=(0, None)) + + if obj is None: + raise ValueError('no active objective was found') + if obj.sense == pe.minimize: + lagrangian = obj.expr + else: + lagrangian = -obj.expr + + for c in all_cons: + if c.equality or (c.lb is not None and c.ub is not None and c.lb == c.ub): + c_hasher = ComponentHasher(c, None) + parent.eq_dual_set.add(c_hasher) + lagrangian += parent.eq_duals[c_hasher] * (c.body - c.lb) + else: + if c.lb is not None: + c_hasher = ComponentHasher(c, 'lb') + parent.ineq_dual_set.add(c_hasher) + lagrangian += parent.ineq_duals[c_hasher] * (c.lb - c.body) + if c.ub is not None: + c_hasher = ComponentHasher(c, 'ub') + parent.ineq_dual_set.add(c_hasher) + lagrangian += parent.ineq_duals[c_hasher] * (c.body - c.ub) + + for v in all_vars: + v_lb, v_ub = v.bounds + if v_lb is not None: + v_hasher = ComponentHasher(v, 'lb') + parent.ineq_dual_set.add(v_hasher) + lagrangian += parent.ineq_duals[v_hasher] * (v_lb - v) + if v_ub is not None: + v_hasher = ComponentHasher(v, 'ub') + parent.ineq_dual_set.add(v_hasher) + lagrangian += parent.ineq_duals[v_hasher] * (v - v_ub) + + parent.grad_lag_set = pe.Set(dimen=1) + parent.grad_lag = pe.Constraint(parent.grad_lag_set) + ders = reverse_sd(lagrangian) + for v in all_vars: + v_hasher = ComponentHasher(v, None) + parent.grad_lag_set.add(v_hasher) + parent.grad_lag[v_hasher] = ders[v] == 0 + + parent.complimentarity = pe.Constraint(parent.ineq_dual_set) + for c_hasher in parent.ineq_dual_set: + c = c_hasher.comp + bound = c_hasher.bound + if bound == 'lb': + if isinstance(c, _GeneralVarData): + expr = c.lb - c + else: + expr = c.lb - c.body + else: + assert bound == 'ub' + if isinstance(c, _GeneralVarData): + expr = c - c.ub + else: + expr = c.body - c.ub + parent.complimentarity[c_hasher] = parent.ineq_duals[c_hasher] * expr == 0 + + return parent diff --git a/pao/duality/tests/test_kkt.py b/pao/duality/tests/test_kkt.py new file mode 100644 index 0000000..3a986df --- /dev/null +++ b/pao/duality/tests/test_kkt.py @@ -0,0 +1,115 @@ +import unittest +import pyomo.environ as pe +from pyomo.contrib import appsi +from pao.duality.kkt import construct_kkt +try: + import coramin + coramin_available = True +except ImportError: + coramin_available = False + + +class TestKKT(unittest.TestCase): + def test_linear(self): + m = pe.ConcreteModel() + m.x = pe.Var(bounds=(-10, 10)) + m.y = pe.Var(bounds=(-10, 10)) + m.obj = pe.Objective(expr=m.y) + m.c1 = pe.Constraint(expr=m.y >= m.x + 2) + m.c2 = pe.Constraint(expr=m.y >= -m.x) + + opt1 = appsi.solvers.Gurobi() + res1 = opt1.solve(m) + self.assertEqual( + res1.termination_condition, + appsi.base.TerminationCondition.optimal + ) + + d = construct_kkt(m) + opt2 = appsi.solvers.Gurobi() + opt2.gurobi_options['nonconvex'] = 2 + m.obj.sense = pe.maximize + res2 = opt2.solve(d) + self.assertEqual( + res2.termination_condition, + appsi.base.TerminationCondition.optimal + ) + + res1 = res1.solution_loader.get_primals() + res2 = res2.solution_loader.get_primals() + + for k, v in res1.items(): + self.assertAlmostEqual(v, res2[k], 5) + + @unittest.skipIf(not coramin_available, "coramin is not available") + def test_milp(self): + m = pe.ConcreteModel() + m.x = pe.Var(bounds=(0, 2), domain=pe.Integers) + m.y = pe.Var(bounds=(0, 2), domain=pe.Integers) + + m.obj = pe.Objective(expr=-2*m.x - 3*m.y) + m.c1 = pe.Constraint(expr=m.x + m.y <= 2) + m.c2 = pe.Constraint(expr=2*m.x + 4*m.y <= 6) + + opt1 = appsi.solvers.Gurobi() + res1 = opt1.solve(m) + self.assertEqual( + res1.termination_condition, + appsi.base.TerminationCondition.optimal + ) + + d = construct_kkt(m) + self.assertEqual(len(d.grad_lag.index_set()), 6) + + coramin.relaxations.relax(d, in_place=True, use_alpha_bb=False) + for b in coramin.relaxations.relaxation_data_objects(d, active=True): + b.rebuild(build_nonlinear_constraint=True) + + opt2 = appsi.solvers.Gurobi() + opt2.gurobi_options['nonconvex'] = 2 + res2 = opt2.solve(d) + self.assertEqual( + res2.termination_condition, + appsi.base.TerminationCondition.optimal + ) + + res1 = res1.solution_loader.get_primals() + res2 = res2.solution_loader.get_primals() + + for k, v in res1.items(): + self.assertAlmostEqual(v, res2[k], 5) + + def test_nonlinear(self): + m = pe.ConcreteModel() + m.x = pe.Var() + m.y = pe.Var() + m.obj = pe.Objective(expr=m.x**2 + m.y**2) + m.c1 = pe.Constraint(expr=m.y >= pe.exp(m.x)) + m.c2 = pe.Constraint(expr=m.y >= (m.x - 1)**2) + + opt1 = appsi.solvers.Ipopt() + res1 = opt1.solve(m) + + d = construct_kkt(m) + opt2 = appsi.solvers.Ipopt() + m.obj.sense = pe.maximize + res2 = opt2.solve(d) + + self.assertEqual( + res1.termination_condition, + appsi.base.TerminationCondition.optimal + ) + self.assertEqual( + res2.termination_condition, + appsi.base.TerminationCondition.optimal + ) + + res1 = res1.solution_loader.get_primals() + res2 = res2.solution_loader.get_primals() + + for k, v in res1.items(): + self.assertAlmostEqual(v, res2[k]) + + +if __name__ == '__main__': + unittest.main() diff --git a/pao/duality/tests/test_nonlinear_dual.py b/pao/duality/tests/test_nonlinear_dual.py new file mode 100644 index 0000000..91ed0a7 --- /dev/null +++ b/pao/duality/tests/test_nonlinear_dual.py @@ -0,0 +1,216 @@ +import unittest +import pyomo.environ as pe +from pyomo.contrib import appsi +from pao.duality.dual import Dual +from itertools import product +from pyomo.core.base.var import ScalarVar + + +class TestNonlinearDual(unittest.TestCase): + def test_linear(self): + m = pe.ConcreteModel() + m.x = pe.Var(bounds=(-10, 10)) + m.y = pe.Var(bounds=(-10, 10)) + m.obj = pe.Objective(expr=m.y) + m.c1 = pe.Constraint(expr=m.y >= m.x + 2) + m.c2 = pe.Constraint(expr=m.y >= -m.x) + + opt = appsi.solvers.Gurobi() + res1 = opt.solve(m) + + d = Dual().dual(m) + res2 = opt.solve(d) + + self.assertAlmostEqual( + res1.best_feasible_objective, + res2.best_feasible_objective + ) + + def test_nonlinear(self): + m = pe.ConcreteModel() + m.x = pe.Var() + m.y = pe.Var() + m.obj = pe.Objective(expr=m.x**2 + m.y**2) + m.c1 = pe.Constraint(expr=m.y >= pe.exp(m.x)) + m.c2 = pe.Constraint(expr=m.y >= (m.x - 1)**2) + + opt = appsi.solvers.Ipopt(only_child_vars=False) + res1 = opt.solve(m) + + d = Dual().dual(m) + res2 = opt.solve(d) + + self.assertAlmostEqual( + res1.best_feasible_objective, + res2.best_feasible_objective + ) + + def test_equality(self): + m = pe.ConcreteModel() + m.x = pe.Var(bounds=(0, 2)) + m.y = pe.Var(bounds=(0, 2)) + + m.obj = pe.Objective(expr=-2*m.x - 3*m.y) + m.c1 = pe.Constraint(expr=m.x + m.y == 2) + m.c2 = pe.Constraint(expr=-2*m.x - 4*m.y >= -6) + + opt1 = appsi.solvers.Gurobi(only_child_vars=False) + opt2 = appsi.solvers.Gurobi(only_child_vars=False) + + pd = Dual() + d = pd.dual(m) + + res1 = opt1.solve(m) + res2 = opt2.solve(d) + + self.assertAlmostEqual( + res1.best_feasible_objective, + res2.best_feasible_objective + ) + + del m.c1 + m.c1 = pe.Constraint(expr=m.x + m.y == 1) + d = pd.dual(m) + res1 = opt1.solve(m) + res2 = opt2.solve(d) + self.assertAlmostEqual( + res1.best_feasible_objective, + res2.best_feasible_objective + ) + + def test_persistent_binary_branching(self): + m = pe.ConcreteModel() + m.x = pe.Var(bounds=(0, 2)) + m.y = pe.Var(bounds=(0, 2)) + + m.obj = pe.Objective(expr=-2*m.x - 3*m.y) + m.c1 = pe.Constraint(expr=m.x + m.y <= 2) + m.c2 = pe.Constraint(expr=2*m.x + 4*m.y <= 6) + + opt1 = appsi.solvers.Gurobi(only_child_vars=False) + opt2 = appsi.solvers.Gurobi(only_child_vars=False) + + opt1.gurobi_options['DualReductions'] = 0 + opt2.gurobi_options['DualReductions'] = 0 + + pd = Dual() + d = pd.dual(m) + + res1 = opt1.solve(m) + res2 = opt2.solve(d) + + self.assertAlmostEqual( + res1.best_feasible_objective, + res2.best_feasible_objective + ) + + opt1.config.load_solution = False + opt2.config.load_solution = False + + for xval, yval in product([None, 0, 1, 2], [None, 0, 1, 2]): + if xval is None: + m.x.unfix() + else: + m.x.fix(xval) + + if yval is None: + m.y.unfix() + else: + m.y.fix(yval) + + res1 = opt1.solve(m) + d = pd.dual(m) + res2 = opt2.solve(d) + + if res1.termination_condition == appsi.base.TerminationCondition.infeasible: + self.assertEqual( + res2.termination_condition, + appsi.base.TerminationCondition.unbounded + ) + else: + self.assertIsNotNone(res1.best_feasible_objective) + self.assertAlmostEqual( + res1.best_feasible_objective, + res2.best_feasible_objective + ) + + m.x.unfix() + m.y.unfix() + + res1 = opt1.solve(m) + d = pd.dual(m) + res2 = opt2.solve(d) + + self.assertIsNotNone(res1.best_feasible_objective) + self.assertAlmostEqual( + res1.best_feasible_objective, + res2.best_feasible_objective + ) + + def test_persistent_changing_bounds(self): + m = pe.ConcreteModel() + m.x = pe.Var() + m.y = pe.Var() + m.obj = pe.Objective(expr=m.x**2 + m.y**2) + m.c1 = pe.Constraint(expr=m.y >= pe.exp(m.x)) + m.c2 = pe.Constraint(expr=m.y >= (m.x - 1)**2) + + opt1 = appsi.solvers.Ipopt(only_child_vars=False) + opt2 = appsi.solvers.Ipopt(only_child_vars=False) + + pd = Dual() + + def check(): + d = pd.dual(m) + res1 = opt1.solve(m) + res2 = opt2.solve(d) + self.assertAlmostEqual( + res1.best_feasible_objective, + res2.best_feasible_objective + ) + + check() + + m.c1.deactivate() + check() + + m.x.setlb(1) + check() + + m.x.setlb(0.1) + check() + + m.x.setlb(None) + m.x.setub(2) + check() + + m.x.setub(0.1) + check() + + m.x.setub(None) + check() + + m.c1.activate() + check() + + def test_fixed_vars(self): + m = pe.ConcreteModel() + m.x = ScalarVar(bounds=(0, 2)) + m.y = ScalarVar(bounds=(0, 2)) + m.c = ScalarVar(bounds=(-3, 3)) + + m.obj = pe.Objective(expr=-2*m.x - m.c*m.y) + m.c1 = pe.Constraint(expr=m.x + m.y <= 2) + m.c2 = pe.Constraint(expr=-2*m.c*m.x + 4*m.y <= 6) + + opt = appsi.solvers.Gurobi(only_child_vars=False) + opt.gurobi_options['nonconvex'] = 2 + + pd = Dual(fixed_vars=[m.c]) + d = pd.dual(m) + + res = opt.solve(d) + self.assertAlmostEqual(res.best_feasible_objective, -2) + self.assertEqual(m.x.value, None) + self.assertEqual(m.y.value, None) + self.assertAlmostEqual(m.c.value, -3) diff --git a/pao/duality/utils.py b/pao/duality/utils.py new file mode 100644 index 0000000..e8baa99 --- /dev/null +++ b/pao/duality/utils.py @@ -0,0 +1,29 @@ + +class ComponentHasher(object): + def __init__(self, var, bound): + self.comp = var + self.bound = bound + + def __eq__(self, other): + if not isinstance(other, ComponentHasher): + return False + if self.comp is other.comp and self.bound is other.bound: + return True + else: + return False + + def __hash__(self): + return hash((id(self.comp), self.bound)) + + def __repr__(self): + first = str(self.comp) + if self.bound is None: + return first + else: + second = str(self.bound) + res = str((first, second)) + res = res.replace("'", "") + return res + + def __str__(self): + return self.__repr__()