From 3a40664bc73892c54955ce1f2c11e6596bdf07f9 Mon Sep 17 00:00:00 2001 From: Michael Bynum Date: Fri, 7 Jul 2023 22:38:10 -0600 Subject: [PATCH 1/9] initial draft of appsi interface to wntr --- pyomo/contrib/appsi/base.py | 2 +- pyomo/contrib/appsi/solvers/__init__.py | 1 + pyomo/contrib/appsi/solvers/wntr.py | 471 ++++++++++++++++++++++++ 3 files changed, 473 insertions(+), 1 deletion(-) create mode 100644 pyomo/contrib/appsi/solvers/wntr.py diff --git a/pyomo/contrib/appsi/base.py b/pyomo/contrib/appsi/base.py index ca7255d5628..00c945235e5 100644 --- a/pyomo/contrib/appsi/base.py +++ b/pyomo/contrib/appsi/base.py @@ -933,7 +933,7 @@ def update_config(self, val: UpdateConfig): def set_instance(self, model): saved_update_config = self.update_config - self.__init__() + self.__init__(only_child_vars=self._only_child_vars) self.update_config = saved_update_config self._model = model if self.use_extensions and cmodel_available: diff --git a/pyomo/contrib/appsi/solvers/__init__.py b/pyomo/contrib/appsi/solvers/__init__.py index df58a0cb245..20755d1eb07 100644 --- a/pyomo/contrib/appsi/solvers/__init__.py +++ b/pyomo/contrib/appsi/solvers/__init__.py @@ -3,3 +3,4 @@ from .cbc import Cbc from .cplex import Cplex from .highs import Highs +from .wntr import Wntr, WntrResults diff --git a/pyomo/contrib/appsi/solvers/wntr.py b/pyomo/contrib/appsi/solvers/wntr.py new file mode 100644 index 00000000000..655a620077c --- /dev/null +++ b/pyomo/contrib/appsi/solvers/wntr.py @@ -0,0 +1,471 @@ +from pyomo.contrib.appsi.base import ( + PersistentBase, + PersistentSolver, + SolverConfig, + Results, + TerminationCondition, + PersistentSolutionLoader +) +from pyomo.core.expr.numeric_expr import ( + ProductExpression, + DivisionExpression, + PowExpression, + SumExpression, + MonomialTermExpression, + NegationExpression, + UnaryFunctionExpression, + LinearExpression, + AbsExpression, + NPV_ProductExpression, + NPV_DivisionExpression, + NPV_PowExpression, + NPV_SumExpression, + NPV_NegationExpression, + NPV_UnaryFunctionExpression, + NPV_AbsExpression, +) +from pyomo.common.errors import PyomoException +from pyomo.common.collections import ComponentMap +from pyomo.core.expr.numvalue import native_numeric_types +from typing import Dict, Optional, List +from pyomo.core.base.block import _BlockData +from pyomo.core.base.var import _GeneralVarData +from pyomo.core.base.param import _ParamData +from pyomo.core.base.constraint import _GeneralConstraintData +from pyomo.common.timing import HierarchicalTimer +from pyomo.core.base import SymbolMap, NumericLabeler, TextLabeler +from pyomo.common.dependencies import attempt_import +from pyomo.core.staleflag import StaleFlagManager +from pyomo.contrib.appsi.cmodel import cmodel, cmodel_available +wntr, wntr_available = attempt_import('wntr') +import wntr +import logging +import time +from pyomo.core.expr.visitor import ExpressionValueVisitor + + +logger = logging.getLogger(__name__) + + +class WntrConfig(SolverConfig): + def __init__( + self, + description=None, + doc=None, + implicit=False, + implicit_domain=None, + visibility=0, + ): + super().__init__( + description=description, + doc=doc, + implicit=implicit, + implicit_domain=implicit_domain, + visibility=visibility, + ) + + +class WntrResults(Results): + def __init__(self, solver): + super().__init__() + self.wallclock_time = None + self.solution_loader = PersistentSolutionLoader(solver=solver) + + +class Wntr(PersistentBase, PersistentSolver): + def __init__(self, only_child_vars=True): + super().__init__(only_child_vars=only_child_vars) + self._config = WntrConfig() + self._solver_options = dict() + self._solver_model = None + self._symbol_map = SymbolMap() + self._labeler = None + self._pyomo_var_to_solver_var_map = dict() + self._pyomo_con_to_solver_con_map = dict() + self._pyomo_param_to_solver_param_map = dict() + self._needs_updated = True + self._last_results_object: Optional[WntrResults] = None + self._pyomo_to_wntr_visitor = PyomoToWntrVisitor( + self._pyomo_var_to_solver_var_map, + self._pyomo_param_to_solver_param_map + ) + + def available(self): + if wntr_available: + return self.Availability.FullLicense + else: + return self.Availability.NotFound + + def version(self): + return tuple(int(i) for i in wntr.__version__.split('.')) + + @property + def config(self) -> WntrConfig: + return self._config + + @config.setter + def config(self, val: WntrConfig): + self._config = val + + @property + def wntr_options(self): + return self._solver_options + + @wntr_options.setter + def wntr_options(self, val: Dict): + self._solver_options = val + + @property + def symbol_map(self): + return self._symbol_map + + def _solve(self, timer: HierarchicalTimer): + t0 = time.time() + opt = wntr.sim.solvers.NewtonSolver(self.wntr_options) + if self._needs_updated: + timer.start('set_structure') + self._solver_model.set_structure() + timer.stop('set_structure') + self._needs_updated = False + timer.start('newton solve') + status, msg, num_iter = opt.solve(self._solver_model) + timer.stop('newton solve') + tf = time.time() + + results = WntrResults(self) + results.wallclock_time = tf - t0 + if status == wntr.sim.solvers.SolverStatus.converged: + results.termination_condition = TerminationCondition.optimal + else: + results.termination_condition = TerminationCondition.error + results.best_feasible_objective = None + results.best_objective_bound = None + + if self.config.load_solution: + if status == wntr.sim.solvers.SolverStatus.converged: + timer.start('load solution') + self.load_vars() + timer.stop('load solution') + else: + raise RuntimeError( + 'A feasible solution was not found, so no solution can be loaded.' + 'Please set opt.config.load_solution=False and check ' + 'results.termination_condition and ' + 'results.best_feasible_objective before loading a solution.' + ) + return results + + def solve(self, model: _BlockData, timer: HierarchicalTimer = None) -> Results: + StaleFlagManager.mark_all_as_stale() + if self._last_results_object is not None: + self._last_results_object.solution_loader.invalidate() + if timer is None: + timer = HierarchicalTimer() + if model is not self._model: + timer.start('set_instance') + self.set_instance(model) + timer.stop('set_instance') + else: + timer.start('update') + self.update(timer=timer) + timer.stop('update') + res = self._solve(timer) + self._last_results_object = res + if self.config.report_timing: + logger.info('\n' + str(timer)) + return res + + def _reinit(self): + saved_config = self.config + saved_options = self.wntr_options + saved_update_config = self.update_config + self.__init__(only_child_vars=self._only_child_vars) + self.config = saved_config + self.wntr_options = saved_options + self.update_config = saved_update_config + + def set_instance(self, model): + if self._last_results_object is not None: + self._last_results_object.solution_loader.invalidate() + if not self.available(): + c = self.__class__ + raise PyomoException( + f'Solver {c.__module__}.{c.__qualname__} is not available ' + f'({self.available()}).' + ) + self._reinit() + self._model = model + if self.use_extensions and cmodel_available: + self._expr_types = cmodel.PyomoExprTypes() + + if self.config.symbolic_solver_labels: + self._labeler = TextLabeler() + else: + self._labeler = NumericLabeler('x') + + self._solver_model = wntr.sim.aml.aml.Model() + + self.add_block(model) + + def _add_variables(self, variables: List[_GeneralVarData]): + aml = wntr.sim.aml.aml + for var in variables: + varname = self._symbol_map.getSymbol(var, self._labeler) + _v, _lb, _ub, _fixed, _domain_interval, _value = self._vars[id(var)] + lb, ub, step = _domain_interval + if ( + _lb is not None + or _ub is not None + or lb is not None + or ub is not None + or step != 0 + ): + raise ValueError(f"WNTR's newton solver only supports continuous variables without bounds: {var.name}") + if _value is None: + _value = 0 + wntr_var = aml.Var(_value) + setattr(self._solver_model, varname, wntr_var) + self._pyomo_var_to_solver_var_map[id(var)] = wntr_var + if _fixed: + self._solver_model._wntr_fixed_var_params[id(var)] = param = aml.Param(_value) + wntr_expr = self._pyomo_to_wntr_visitor.dfs_postorder_stack(var - param) + self._solver_model._wntr_fixed_var_cons[id(var)] = aml.Constraint(wntr_expr) + self._needs_updated = True + + def _add_params(self, params: List[_ParamData]): + aml = wntr.sim.aml.aml + for p in params: + pname = self._symbol_map.getSymbol(p, self._labeler) + wntr_p = aml.Param(p.value) + setattr(self._solver_model, pname, wntr_p) + self._pyomo_param_to_solver_param_map[id(p)] = wntr_p + + def _add_constraints(self, cons: List[_GeneralConstraintData]): + aml = wntr.sim.aml.aml + for con in cons: + if not con.equality: + raise ValueError(f"WNTR's newtwon solver only supports equality constraints: {con.name}") + conname = self._symbol_map.getSymbol(con, self._labeler) + wntr_expr = self._pyomo_to_wntr_visitor.dfs_postorder_stack(con.body - con.upper) + wntr_con = aml.Constraint(wntr_expr) + setattr(self._solver_model, conname, wntr_con) + self._pyomo_con_to_solver_con_map[con] = wntr_con + self._needs_updated = True + + def _remove_constraints(self, cons: List[_GeneralConstraintData]): + for con in cons: + solver_con = self._pyomo_con_to_solver_con_map[con] + delattr(self._solver_model, solver_con.name) + self._symbol_map.removeSymbol(con) + del self._pyomo_con_to_solver_con_map[con] + self._needs_updated = True + + def _remove_variables(self, variables: List[_GeneralVarData]): + for var in variables: + v_id = id(var) + solver_var = self._pyomo_var_to_solver_var_map[v_id] + delattr(self._solver_model, solver_var.name) + self._symbol_map.removeSymbol(var) + del self._pyomo_var_to_solver_var_map[v_id] + if v_id in self._solver_model._wntr_fixed_var_params: + del self._solver_model._wntr_fixed_var_params[v_id] + del self._solver_model._wntr_fixed_var_cons[v_id] + self._needs_updated = True + + def _remove_params(self, params: List[_ParamData]): + for p in params: + p_id = id(p) + solver_param = self._pyomo_param_to_solver_param_map[p_id] + delattr(self._solver_model, solver_param.name) + self._symbol_map.removeSymbol(p) + del self._pyomo_param_to_solver_param_map[p_id] + + def _update_variables(self, variables: List[_GeneralVarData]): + for var in variables: + v_id = id(var) + solver_var = self._pyomo_var_to_solver_var_map[v_id] + _v, _lb, _ub, _fixed, _domain_interval, _value = self._vars[v_id] + lb, ub, step = _domain_interval + if ( + _lb is not None + or _ub is not None + or lb is not None + or ub is not None + or step != 0 + ): + raise ValueError(f"WNTR's newton solver only supports continuous variables without bounds: {var.name}") + if _value is None: + _value = 0 + solver_var.value = _value + if _fixed: + if v_id not in self._solver_model._wntr_fixed_var_params: + self._solver_model._wntr_fixed_var_params[v_id] = param = aml.Param(_value) + wntr_expr = self._pyomo_to_wntr_visitor.dfs_postorder_stack(var - param) + self._solver_model._wntr_fixed_var_cons[v_id] = aml.Constraint(wntr_expr) + self._needs_updated = True + else: + self._solver_model._wntr_fixed_var_params[v_id].value = _value + else: + if v_id in self._solver_model._wntr_fixed_var_params: + del self._solver_model._wntr_fixed_var_params[v_id] + del self._solver_model._wntr_fixed_var_cons[v_id] + self._needs_updated = True + + def update_params(self): + for p_id, solver_p in self._pyomo_param_to_solver_param_map.items(): + p = self._params[p_id] + solver_p.value = p.value + + def _set_objective(self, obj): + raise NotImplementedError(f"WNTR's newton solver can only solve square problems") + + def load_vars(self, vars_to_load=None): + if vars_to_load is None: + vars_to_load = [i[0] for i in self._vars.values()] + for v in vars_to_load: + v_id = id(v) + solver_v = self._pyomo_var_to_solver_var_map[v_id] + v.value = solver_v.value + + def get_primals(self, vars_to_load=None): + if vars_to_load is None: + vars_to_load = [i[0] for i in self._vars.values()] + res = ComponentMap() + for v in vars_to_load: + v_id = id(v) + solver_v = self._pyomo_var_to_solver_var_map[v_id] + res[v] = solver_v.value + + def _add_sos_constraints(self, cons): + if len(cons) > 0: + raise NotImplementedError(f"WNTR's newton solver does not support SOS constraints") + + def _remove_sos_constraints(self, cons): + if len(cons) > 0: + raise NotImplementedError(f"WNTR's newton solver does not support SOS constraints") + + +def _handle_product_expression(node, values): + arg1, arg2 = values + return arg1 * arg2 + + +def _handle_sum_expression(node, values): + return sum(values) + + +def _handle_division_expression(node, values): + arg1, arg2 = values + return arg1 / arg2 + + +def _handle_pow_expression(node, values): + arg1, arg2 = values + return arg1 ** arg2 + + +def _handle_negation_expression(node, values): + return -values[0] + + +def _handle_exp_expression(node, values): + return wntr.sim.aml.exp(values[0]) + + +def _handle_log_expression(node, values): + return wntr.sim.aml.log(values[0]) + + +def _handle_sin_expression(node, values): + return wntr.sim.aml.sin(values[0]) + + +def _handle_cos_expression(node, values): + return wntr.sim.aml.cos(values[0]) + + +def _handle_tan_expression(node, values): + return wntr.sim.aml.tan(values[0]) + + +def _handle_asin_expression(node, values): + return wntr.sim.aml.asin(values[0]) + + +def _handle_acos_expression(node, values): + return wntr.sim.aml.acos(values[0]) + + +def _handle_atan_expression(node, values): + return wntr.sim.aml.atan(values[0]) + + +def _handle_sqrt_expression(node, values): + return (values[0])**0.5 + + +def _handle_abs_expression(node, values): + return wntr.sim.aml.abs(values[0]) + + +_unary_handler_map = dict() +_unary_handler_map['exp'] = _handle_exp_expression +_unary_handler_map['log'] = _handle_log_expression +_unary_handler_map['sin'] = _handle_sin_expression +_unary_handler_map['cos'] = _handle_cos_expression +_unary_handler_map['tan'] = _handle_tan_expression +_unary_handler_map['asin'] = _handle_asin_expression +_unary_handler_map['acos'] = _handle_acos_expression +_unary_handler_map['atan'] = _handle_atan_expression +_unary_handler_map['sqrt'] = _handle_sqrt_expression +_unary_handler_map['abs'] = _handle_abs_expression + + +def _handle_unary_function_expression(node, values): + if node.getname() in _unary_handler_map: + return _unary_handler_map[node.getname()](node, values) + else: + raise NotImplementedError(f'Unrecognized unary function expression: {node.getname()}') + + +_handler_map = dict() +_handler_map[ProductExpression] = _handle_product_expression +_handler_map[DivisionExpression] = _handle_division_expression +_handler_map[PowExpression] = _handle_pow_expression +_handler_map[SumExpression] = _handle_sum_expression +_handler_map[MonomialTermExpression] = _handle_product_expression +_handler_map[NegationExpression] = _handle_negation_expression +_handler_map[UnaryFunctionExpression] = _handle_unary_function_expression +_handler_map[LinearExpression] = _handle_sum_expression +_handler_map[AbsExpression] = _handle_abs_expression +_handler_map[NPV_ProductExpression] = _handle_product_expression +_handler_map[NPV_DivisionExpression] = _handle_division_expression +_handler_map[NPV_PowExpression] = _handle_pow_expression +_handler_map[NPV_SumExpression] = _handle_sum_expression +_handler_map[NPV_NegationExpression] = _handle_negation_expression +_handler_map[NPV_UnaryFunctionExpression] = _handle_unary_function_expression +_handler_map[NPV_AbsExpression] = _handle_abs_expression + + +class PyomoToWntrVisitor(ExpressionValueVisitor): + def __init__(self, var_map, param_map): + self.var_map = var_map + self.param_map = param_map + + def visit(self, node, values): + if node.__class__ in _handler_map: + return _handler_map[node.__class__](node, values) + else: + raise NotImplementedError(f'Unrecognized expression type: {node.__class__}') + + def visiting_potential_leaf(self, node): + if node.__class__ in native_numeric_types: + return True, node + + if node.is_variable_type(): + return True, self.var_map[id(node)] + + if node.is_parameter_type(): + return True, self.param_map[id(node)] + + return False, None From a5b1eb35c52982161400445e1875279e2d2a3c7b Mon Sep 17 00:00:00 2001 From: Michael Bynum Date: Sun, 9 Jul 2023 20:19:05 -0600 Subject: [PATCH 2/9] tests for persistent interface to wntr --- .../solvers/tests/test_wntr_persistent.py | 184 ++++++++++++++++++ pyomo/contrib/appsi/solvers/wntr.py | 30 ++- 2 files changed, 210 insertions(+), 4 deletions(-) create mode 100644 pyomo/contrib/appsi/solvers/tests/test_wntr_persistent.py diff --git a/pyomo/contrib/appsi/solvers/tests/test_wntr_persistent.py b/pyomo/contrib/appsi/solvers/tests/test_wntr_persistent.py new file mode 100644 index 00000000000..b172a9204e8 --- /dev/null +++ b/pyomo/contrib/appsi/solvers/tests/test_wntr_persistent.py @@ -0,0 +1,184 @@ +import pyomo.environ as pe +import pyomo.common.unittest as unittest +from pyomo.contrib.appsi.base import TerminationCondition, Results, PersistentSolver +from pyomo.contrib.appsi.solvers.wntr import Wntr, wntr_available +import math + + +_default_wntr_options = dict( + TOL=1e-8, +) + + +@unittest.skipUnless(wntr_available, 'wntr is not available') +class TestWntrPersistent(unittest.TestCase): + def test_param_updates(self): + m = pe.ConcreteModel() + m.x = pe.Var() + m.p = pe.Param(initialize=1, mutable=True) + m.c = pe.Constraint(expr=m.x == m.p) + opt = Wntr() + opt.wntr_options.update(_default_wntr_options) + res = opt.solve(m) + self.assertEqual(res.termination_condition, TerminationCondition.optimal) + self.assertAlmostEqual(m.x.value, 1) + + m.p.value = 2 + res = opt.solve(m) + self.assertEqual(res.termination_condition, TerminationCondition.optimal) + self.assertAlmostEqual(m.x.value, 2) + + def test_remove_add_constraint(self): + m = pe.ConcreteModel() + m.x = pe.Var() + m.y = pe.Var() + m.c1 = pe.Constraint(expr=m.y == (m.x - 1)**2) + m.c2 = pe.Constraint(expr=m.y == pe.exp(m.x)) + opt = Wntr() + opt.config.symbolic_solver_labels = True + opt.wntr_options.update(_default_wntr_options) + res = opt.solve(m) + self.assertEqual(res.termination_condition, TerminationCondition.optimal) + self.assertAlmostEqual(m.x.value, 0) + self.assertAlmostEqual(m.y.value, 1) + + del m.c2 + m.c2 = pe.Constraint(expr=m.y == pe.log(m.x)) + m.x.value = 0.5 + m.y.value = 0.5 + res = opt.solve(m) + self.assertEqual(res.termination_condition, TerminationCondition.optimal) + self.assertAlmostEqual(m.x.value, 1) + self.assertAlmostEqual(m.y.value, 0) + + def test_fixed_var(self): + m = pe.ConcreteModel() + m.x = pe.Var() + m.y = pe.Var() + m.c1 = pe.Constraint(expr=m.y == (m.x - 1)**2) + m.x.fix(0.5) + opt = Wntr() + opt.wntr_options.update(_default_wntr_options) + res = opt.solve(m) + self.assertEqual(res.termination_condition, TerminationCondition.optimal) + self.assertAlmostEqual(m.x.value, 0.5) + self.assertAlmostEqual(m.y.value, 0.25) + + m.x.unfix() + m.c2 = pe.Constraint(expr=m.y == pe.exp(m.x)) + res = opt.solve(m) + self.assertEqual(res.termination_condition, TerminationCondition.optimal) + self.assertAlmostEqual(m.x.value, 0) + self.assertAlmostEqual(m.y.value, 1) + + m.x.fix(0.5) + del m.c2 + res = opt.solve(m) + self.assertEqual(res.termination_condition, TerminationCondition.optimal) + self.assertAlmostEqual(m.x.value, 0.5) + self.assertAlmostEqual(m.y.value, 0.25) + + def test_remove_variables_params(self): + m = pe.ConcreteModel() + m.x = pe.Var() + m.y = pe.Var() + m.z = pe.Var() + m.z.fix(0) + m.px = pe.Param(mutable=True, initialize=1) + m.py = pe.Param(mutable=True, initialize=1) + m.c1 = pe.Constraint(expr=m.x == m.px) + m.c2 = pe.Constraint(expr=m.y == m.py) + opt = Wntr() + opt.wntr_options.update(_default_wntr_options) + res = opt.solve(m) + self.assertEqual(res.termination_condition, TerminationCondition.optimal) + self.assertAlmostEqual(m.x.value, 1) + self.assertAlmostEqual(m.y.value, 1) + self.assertAlmostEqual(m.z.value, 0) + + del m.c2 + del m.y + del m.py + m.z.value = 2 + m.px.value = 2 + res = opt.solve(m) + self.assertEqual(res.termination_condition, TerminationCondition.optimal) + self.assertAlmostEqual(m.x.value, 2) + self.assertAlmostEqual(m.z.value, 2) + + del m.z + m.px.value = 3 + res = opt.solve(m) + self.assertEqual(res.termination_condition, TerminationCondition.optimal) + self.assertAlmostEqual(m.x.value, 3) + + def test_get_primals(self): + m = pe.ConcreteModel() + m.x = pe.Var() + m.y = pe.Var() + m.c1 = pe.Constraint(expr=m.y == (m.x - 1)**2) + m.c2 = pe.Constraint(expr=m.y == pe.exp(m.x)) + opt = Wntr() + opt.config.load_solution = False + opt.wntr_options.update(_default_wntr_options) + res = opt.solve(m) + self.assertEqual(res.termination_condition, TerminationCondition.optimal) + self.assertAlmostEqual(m.x.value, None) + self.assertAlmostEqual(m.y.value, None) + primals = opt.get_primals() + self.assertAlmostEqual(primals[m.x], 0) + self.assertAlmostEqual(primals[m.y], 1) + + def test_operators(self): + m = pe.ConcreteModel() + m.x = pe.Var(initialize=1) + m.c1 = pe.Constraint(expr=2/m.x == 1) + opt = Wntr() + opt.wntr_options.update(_default_wntr_options) + res = opt.solve(m) + self.assertEqual(res.termination_condition, TerminationCondition.optimal) + self.assertAlmostEqual(m.x.value, 2) + + del m.c1 + m.x.value = 0 + m.c1 = pe.Constraint(expr=pe.sin(m.x) == math.sin(math.pi/4)) + res = opt.solve(m) + self.assertEqual(res.termination_condition, TerminationCondition.optimal) + self.assertAlmostEqual(m.x.value, math.pi/4) + + del m.c1 + m.c1 = pe.Constraint(expr=pe.cos(m.x) == 0) + res = opt.solve(m) + self.assertEqual(res.termination_condition, TerminationCondition.optimal) + self.assertAlmostEqual(m.x.value, math.pi/2) + + del m.c1 + m.c1 = pe.Constraint(expr=pe.tan(m.x) == 1) + m.x.value = 0 + res = opt.solve(m) + self.assertEqual(res.termination_condition, TerminationCondition.optimal) + self.assertAlmostEqual(m.x.value, math.pi/4) + + del m.c1 + m.c1 = pe.Constraint(expr=pe.asin(m.x) == math.asin(0.5)) + res = opt.solve(m) + self.assertEqual(res.termination_condition, TerminationCondition.optimal) + self.assertAlmostEqual(m.x.value, 0.5) + + del m.c1 + m.c1 = pe.Constraint(expr=pe.acos(m.x) == math.acos(0.6)) + res = opt.solve(m) + self.assertEqual(res.termination_condition, TerminationCondition.optimal) + self.assertAlmostEqual(m.x.value, 0.6) + + del m.c1 + m.c1 = pe.Constraint(expr=pe.atan(m.x) == math.atan(0.5)) + res = opt.solve(m) + self.assertEqual(res.termination_condition, TerminationCondition.optimal) + self.assertAlmostEqual(m.x.value, 0.5) + + del m.c1 + m.c1 = pe.Constraint(expr=pe.sqrt(m.x) == math.sqrt(0.6)) + res = opt.solve(m) + self.assertEqual(res.termination_condition, TerminationCondition.optimal) + self.assertAlmostEqual(m.x.value, 0.6) diff --git a/pyomo/contrib/appsi/solvers/wntr.py b/pyomo/contrib/appsi/solvers/wntr.py index 655a620077c..1e0952cd867 100644 --- a/pyomo/contrib/appsi/solvers/wntr.py +++ b/pyomo/contrib/appsi/solvers/wntr.py @@ -41,6 +41,7 @@ import wntr import logging import time +import sys from pyomo.core.expr.visitor import ExpressionValueVisitor @@ -120,15 +121,25 @@ def symbol_map(self): return self._symbol_map def _solve(self, timer: HierarchicalTimer): + options = dict() + if self.config.time_limit is not None: + options['TIME_LIMIT'] = self.config.time_limit + options.update(self.wntr_options) + opt = wntr.sim.solvers.NewtonSolver(options) + + if self.config.stream_solver: + ostream = sys.stdout + else: + ostream = None + t0 = time.time() - opt = wntr.sim.solvers.NewtonSolver(self.wntr_options) if self._needs_updated: timer.start('set_structure') self._solver_model.set_structure() timer.stop('set_structure') self._needs_updated = False timer.start('newton solve') - status, msg, num_iter = opt.solve(self._solver_model) + status, msg, num_iter = opt.solve(self._solver_model, ostream) timer.stop('newton solve') tf = time.time() @@ -168,6 +179,13 @@ def solve(self, model: _BlockData, timer: HierarchicalTimer = None) -> Results: else: timer.start('update') self.update(timer=timer) + timer.start('initial values') + for v_id, solver_v in self._pyomo_var_to_solver_var_map.items(): + pyomo_v = self._vars[v_id][0] + val = pyomo_v.value + if val is not None: + solver_v.value = val + timer.stop('initial values') timer.stop('update') res = self._solve(timer) self._last_results_object = res @@ -204,6 +222,8 @@ def set_instance(self, model): self._labeler = NumericLabeler('x') self._solver_model = wntr.sim.aml.aml.Model() + self._solver_model._wntr_fixed_var_params = wntr.sim.aml.aml.ParamDict() + self._solver_model._wntr_fixed_var_cons = wntr.sim.aml.aml.ConstraintDict() self.add_block(model) @@ -228,7 +248,7 @@ def _add_variables(self, variables: List[_GeneralVarData]): self._pyomo_var_to_solver_var_map[id(var)] = wntr_var if _fixed: self._solver_model._wntr_fixed_var_params[id(var)] = param = aml.Param(_value) - wntr_expr = self._pyomo_to_wntr_visitor.dfs_postorder_stack(var - param) + wntr_expr = wntr_var - param self._solver_model._wntr_fixed_var_cons[id(var)] = aml.Constraint(wntr_expr) self._needs_updated = True @@ -281,6 +301,7 @@ def _remove_params(self, params: List[_ParamData]): del self._pyomo_param_to_solver_param_map[p_id] def _update_variables(self, variables: List[_GeneralVarData]): + aml = wntr.sim.aml.aml for var in variables: v_id = id(var) solver_var = self._pyomo_var_to_solver_var_map[v_id] @@ -300,7 +321,7 @@ def _update_variables(self, variables: List[_GeneralVarData]): if _fixed: if v_id not in self._solver_model._wntr_fixed_var_params: self._solver_model._wntr_fixed_var_params[v_id] = param = aml.Param(_value) - wntr_expr = self._pyomo_to_wntr_visitor.dfs_postorder_stack(var - param) + wntr_expr = solver_var - param self._solver_model._wntr_fixed_var_cons[v_id] = aml.Constraint(wntr_expr) self._needs_updated = True else: @@ -335,6 +356,7 @@ def get_primals(self, vars_to_load=None): v_id = id(v) solver_v = self._pyomo_var_to_solver_var_map[v_id] res[v] = solver_v.value + return res def _add_sos_constraints(self, cons): if len(cons) > 0: From 797c918e9c01ae07b2ad239673200df40cc61e95 Mon Sep 17 00:00:00 2001 From: Michael Bynum Date: Sun, 9 Jul 2023 20:19:53 -0600 Subject: [PATCH 3/9] tests for persistent interface to wntr --- pyomo/contrib/appsi/solvers/wntr.py | 1 - 1 file changed, 1 deletion(-) diff --git a/pyomo/contrib/appsi/solvers/wntr.py b/pyomo/contrib/appsi/solvers/wntr.py index 1e0952cd867..09fe53dbb98 100644 --- a/pyomo/contrib/appsi/solvers/wntr.py +++ b/pyomo/contrib/appsi/solvers/wntr.py @@ -38,7 +38,6 @@ from pyomo.core.staleflag import StaleFlagManager from pyomo.contrib.appsi.cmodel import cmodel, cmodel_available wntr, wntr_available = attempt_import('wntr') -import wntr import logging import time import sys From 7065274ec8c0265b198395e0e0bb9e9b6f7b6162 Mon Sep 17 00:00:00 2001 From: Michael Bynum Date: Mon, 10 Jul 2023 05:47:48 -0600 Subject: [PATCH 4/9] add wntr to GHA --- .github/workflows/test_branches.yml | 2 ++ .github/workflows/test_pr_and_main.yml | 2 ++ 2 files changed, 4 insertions(+) diff --git a/.github/workflows/test_branches.yml b/.github/workflows/test_branches.yml index 48feb7a4465..d6de4022f4c 100644 --- a/.github/workflows/test_branches.yml +++ b/.github/workflows/test_branches.yml @@ -280,6 +280,8 @@ jobs: || echo "WARNING: Gurobi is not available" python -m pip install --cache-dir cache/pip xpress \ || echo "WARNING: Xpress Community Edition is not available" + python -m pip install --cache-dir cache/pip wntr \ + || echo "WARNING: WNTR is not available" fi python -c 'import sys; print("PYTHON_EXE=%s" \ % (sys.executable,))' >> $GITHUB_ENV diff --git a/.github/workflows/test_pr_and_main.yml b/.github/workflows/test_pr_and_main.yml index fd52c9610a8..5d3503f17c7 100644 --- a/.github/workflows/test_pr_and_main.yml +++ b/.github/workflows/test_pr_and_main.yml @@ -298,6 +298,8 @@ jobs: || echo "WARNING: Gurobi is not available" python -m pip install --cache-dir cache/pip xpress \ || echo "WARNING: Xpress Community Edition is not available" + python -m pip install --cache-dir cache/pip wntr \ + || echo "WARNING: WNTR is not available" fi python -c 'import sys; print("PYTHON_EXE=%s" \ % (sys.executable,))' >> $GITHUB_ENV From b7c54f3ebd9cdd2a431e59f37024a7d39888f9fc Mon Sep 17 00:00:00 2001 From: Michael Bynum Date: Mon, 10 Jul 2023 05:57:03 -0600 Subject: [PATCH 5/9] run black --- .../solvers/tests/test_wntr_persistent.py | 20 +++-- pyomo/contrib/appsi/solvers/wntr.py | 78 ++++++++++++------- 2 files changed, 60 insertions(+), 38 deletions(-) diff --git a/pyomo/contrib/appsi/solvers/tests/test_wntr_persistent.py b/pyomo/contrib/appsi/solvers/tests/test_wntr_persistent.py index b172a9204e8..d250923f104 100644 --- a/pyomo/contrib/appsi/solvers/tests/test_wntr_persistent.py +++ b/pyomo/contrib/appsi/solvers/tests/test_wntr_persistent.py @@ -5,9 +5,7 @@ import math -_default_wntr_options = dict( - TOL=1e-8, -) +_default_wntr_options = dict(TOL=1e-8) @unittest.skipUnless(wntr_available, 'wntr is not available') @@ -32,7 +30,7 @@ def test_remove_add_constraint(self): m = pe.ConcreteModel() m.x = pe.Var() m.y = pe.Var() - m.c1 = pe.Constraint(expr=m.y == (m.x - 1)**2) + m.c1 = pe.Constraint(expr=m.y == (m.x - 1) ** 2) m.c2 = pe.Constraint(expr=m.y == pe.exp(m.x)) opt = Wntr() opt.config.symbolic_solver_labels = True @@ -55,7 +53,7 @@ def test_fixed_var(self): m = pe.ConcreteModel() m.x = pe.Var() m.y = pe.Var() - m.c1 = pe.Constraint(expr=m.y == (m.x - 1)**2) + m.c1 = pe.Constraint(expr=m.y == (m.x - 1) ** 2) m.x.fix(0.5) opt = Wntr() opt.wntr_options.update(_default_wntr_options) @@ -116,7 +114,7 @@ def test_get_primals(self): m = pe.ConcreteModel() m.x = pe.Var() m.y = pe.Var() - m.c1 = pe.Constraint(expr=m.y == (m.x - 1)**2) + m.c1 = pe.Constraint(expr=m.y == (m.x - 1) ** 2) m.c2 = pe.Constraint(expr=m.y == pe.exp(m.x)) opt = Wntr() opt.config.load_solution = False @@ -132,7 +130,7 @@ def test_get_primals(self): def test_operators(self): m = pe.ConcreteModel() m.x = pe.Var(initialize=1) - m.c1 = pe.Constraint(expr=2/m.x == 1) + m.c1 = pe.Constraint(expr=2 / m.x == 1) opt = Wntr() opt.wntr_options.update(_default_wntr_options) res = opt.solve(m) @@ -141,23 +139,23 @@ def test_operators(self): del m.c1 m.x.value = 0 - m.c1 = pe.Constraint(expr=pe.sin(m.x) == math.sin(math.pi/4)) + m.c1 = pe.Constraint(expr=pe.sin(m.x) == math.sin(math.pi / 4)) res = opt.solve(m) self.assertEqual(res.termination_condition, TerminationCondition.optimal) - self.assertAlmostEqual(m.x.value, math.pi/4) + self.assertAlmostEqual(m.x.value, math.pi / 4) del m.c1 m.c1 = pe.Constraint(expr=pe.cos(m.x) == 0) res = opt.solve(m) self.assertEqual(res.termination_condition, TerminationCondition.optimal) - self.assertAlmostEqual(m.x.value, math.pi/2) + self.assertAlmostEqual(m.x.value, math.pi / 2) del m.c1 m.c1 = pe.Constraint(expr=pe.tan(m.x) == 1) m.x.value = 0 res = opt.solve(m) self.assertEqual(res.termination_condition, TerminationCondition.optimal) - self.assertAlmostEqual(m.x.value, math.pi/4) + self.assertAlmostEqual(m.x.value, math.pi / 4) del m.c1 m.c1 = pe.Constraint(expr=pe.asin(m.x) == math.asin(0.5)) diff --git a/pyomo/contrib/appsi/solvers/wntr.py b/pyomo/contrib/appsi/solvers/wntr.py index 09fe53dbb98..0a358c6aedf 100644 --- a/pyomo/contrib/appsi/solvers/wntr.py +++ b/pyomo/contrib/appsi/solvers/wntr.py @@ -4,7 +4,7 @@ SolverConfig, Results, TerminationCondition, - PersistentSolutionLoader + PersistentSolutionLoader, ) from pyomo.core.expr.numeric_expr import ( ProductExpression, @@ -37,6 +37,7 @@ from pyomo.common.dependencies import attempt_import from pyomo.core.staleflag import StaleFlagManager from pyomo.contrib.appsi.cmodel import cmodel, cmodel_available + wntr, wntr_available = attempt_import('wntr') import logging import time @@ -86,8 +87,7 @@ def __init__(self, only_child_vars=True): self._needs_updated = True self._last_results_object: Optional[WntrResults] = None self._pyomo_to_wntr_visitor = PyomoToWntrVisitor( - self._pyomo_var_to_solver_var_map, - self._pyomo_param_to_solver_param_map + self._pyomo_var_to_solver_var_map, self._pyomo_param_to_solver_param_map ) def available(self): @@ -233,22 +233,28 @@ def _add_variables(self, variables: List[_GeneralVarData]): _v, _lb, _ub, _fixed, _domain_interval, _value = self._vars[id(var)] lb, ub, step = _domain_interval if ( - _lb is not None - or _ub is not None - or lb is not None - or ub is not None - or step != 0 + _lb is not None + or _ub is not None + or lb is not None + or ub is not None + or step != 0 ): - raise ValueError(f"WNTR's newton solver only supports continuous variables without bounds: {var.name}") + raise ValueError( + f"WNTR's newton solver only supports continuous variables without bounds: {var.name}" + ) if _value is None: _value = 0 wntr_var = aml.Var(_value) setattr(self._solver_model, varname, wntr_var) self._pyomo_var_to_solver_var_map[id(var)] = wntr_var if _fixed: - self._solver_model._wntr_fixed_var_params[id(var)] = param = aml.Param(_value) + self._solver_model._wntr_fixed_var_params[id(var)] = param = aml.Param( + _value + ) wntr_expr = wntr_var - param - self._solver_model._wntr_fixed_var_cons[id(var)] = aml.Constraint(wntr_expr) + self._solver_model._wntr_fixed_var_cons[id(var)] = aml.Constraint( + wntr_expr + ) self._needs_updated = True def _add_params(self, params: List[_ParamData]): @@ -263,9 +269,13 @@ def _add_constraints(self, cons: List[_GeneralConstraintData]): aml = wntr.sim.aml.aml for con in cons: if not con.equality: - raise ValueError(f"WNTR's newtwon solver only supports equality constraints: {con.name}") + raise ValueError( + f"WNTR's newtwon solver only supports equality constraints: {con.name}" + ) conname = self._symbol_map.getSymbol(con, self._labeler) - wntr_expr = self._pyomo_to_wntr_visitor.dfs_postorder_stack(con.body - con.upper) + wntr_expr = self._pyomo_to_wntr_visitor.dfs_postorder_stack( + con.body - con.upper + ) wntr_con = aml.Constraint(wntr_expr) setattr(self._solver_model, conname, wntr_con) self._pyomo_con_to_solver_con_map[con] = wntr_con @@ -307,21 +317,27 @@ def _update_variables(self, variables: List[_GeneralVarData]): _v, _lb, _ub, _fixed, _domain_interval, _value = self._vars[v_id] lb, ub, step = _domain_interval if ( - _lb is not None - or _ub is not None - or lb is not None - or ub is not None - or step != 0 + _lb is not None + or _ub is not None + or lb is not None + or ub is not None + or step != 0 ): - raise ValueError(f"WNTR's newton solver only supports continuous variables without bounds: {var.name}") + raise ValueError( + f"WNTR's newton solver only supports continuous variables without bounds: {var.name}" + ) if _value is None: _value = 0 solver_var.value = _value if _fixed: if v_id not in self._solver_model._wntr_fixed_var_params: - self._solver_model._wntr_fixed_var_params[v_id] = param = aml.Param(_value) + self._solver_model._wntr_fixed_var_params[v_id] = param = aml.Param( + _value + ) wntr_expr = solver_var - param - self._solver_model._wntr_fixed_var_cons[v_id] = aml.Constraint(wntr_expr) + self._solver_model._wntr_fixed_var_cons[v_id] = aml.Constraint( + wntr_expr + ) self._needs_updated = True else: self._solver_model._wntr_fixed_var_params[v_id].value = _value @@ -337,7 +353,9 @@ def update_params(self): solver_p.value = p.value def _set_objective(self, obj): - raise NotImplementedError(f"WNTR's newton solver can only solve square problems") + raise NotImplementedError( + f"WNTR's newton solver can only solve square problems" + ) def load_vars(self, vars_to_load=None): if vars_to_load is None: @@ -359,11 +377,15 @@ def get_primals(self, vars_to_load=None): def _add_sos_constraints(self, cons): if len(cons) > 0: - raise NotImplementedError(f"WNTR's newton solver does not support SOS constraints") + raise NotImplementedError( + f"WNTR's newton solver does not support SOS constraints" + ) def _remove_sos_constraints(self, cons): if len(cons) > 0: - raise NotImplementedError(f"WNTR's newton solver does not support SOS constraints") + raise NotImplementedError( + f"WNTR's newton solver does not support SOS constraints" + ) def _handle_product_expression(node, values): @@ -382,7 +404,7 @@ def _handle_division_expression(node, values): def _handle_pow_expression(node, values): arg1, arg2 = values - return arg1 ** arg2 + return arg1**arg2 def _handle_negation_expression(node, values): @@ -422,7 +444,7 @@ def _handle_atan_expression(node, values): def _handle_sqrt_expression(node, values): - return (values[0])**0.5 + return (values[0]) ** 0.5 def _handle_abs_expression(node, values): @@ -446,7 +468,9 @@ def _handle_unary_function_expression(node, values): if node.getname() in _unary_handler_map: return _unary_handler_map[node.getname()](node, values) else: - raise NotImplementedError(f'Unrecognized unary function expression: {node.getname()}') + raise NotImplementedError( + f'Unrecognized unary function expression: {node.getname()}' + ) _handler_map = dict() From 2451c444c3f923d83131974e6048d1188f85a894 Mon Sep 17 00:00:00 2001 From: Michael Bynum Date: Mon, 10 Jul 2023 07:19:48 -0600 Subject: [PATCH 6/9] add wntr to GHA --- .github/workflows/test_branches.yml | 2 +- .github/workflows/test_pr_and_main.yml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/test_branches.yml b/.github/workflows/test_branches.yml index d6de4022f4c..69028e55a17 100644 --- a/.github/workflows/test_branches.yml +++ b/.github/workflows/test_branches.yml @@ -280,7 +280,7 @@ jobs: || echo "WARNING: Gurobi is not available" python -m pip install --cache-dir cache/pip xpress \ || echo "WARNING: Xpress Community Edition is not available" - python -m pip install --cache-dir cache/pip wntr \ + python -m pip install git+https://github.com/michaelbynum/wntr.git@working \ || echo "WARNING: WNTR is not available" fi python -c 'import sys; print("PYTHON_EXE=%s" \ diff --git a/.github/workflows/test_pr_and_main.yml b/.github/workflows/test_pr_and_main.yml index 5d3503f17c7..357a2fe866e 100644 --- a/.github/workflows/test_pr_and_main.yml +++ b/.github/workflows/test_pr_and_main.yml @@ -298,7 +298,7 @@ jobs: || echo "WARNING: Gurobi is not available" python -m pip install --cache-dir cache/pip xpress \ || echo "WARNING: Xpress Community Edition is not available" - python -m pip install --cache-dir cache/pip wntr \ + python -m pip install git+https://github.com/michaelbynum/wntr.git@working \ || echo "WARNING: WNTR is not available" fi python -c 'import sys; print("PYTHON_EXE=%s" \ From 26375cb1e5363daddd4dc1d896ba6c7f54854912 Mon Sep 17 00:00:00 2001 From: Michael Bynum Date: Fri, 18 Aug 2023 08:33:56 -0600 Subject: [PATCH 7/9] update workflows --- .github/workflows/test_branches.yml | 2 +- .github/workflows/test_pr_and_main.yml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/test_branches.yml b/.github/workflows/test_branches.yml index 69028e55a17..1c865462f60 100644 --- a/.github/workflows/test_branches.yml +++ b/.github/workflows/test_branches.yml @@ -280,7 +280,7 @@ jobs: || echo "WARNING: Gurobi is not available" python -m pip install --cache-dir cache/pip xpress \ || echo "WARNING: Xpress Community Edition is not available" - python -m pip install git+https://github.com/michaelbynum/wntr.git@working \ + python -m pip install git+https://github.com/usepa/wntr.git@main \ || echo "WARNING: WNTR is not available" fi python -c 'import sys; print("PYTHON_EXE=%s" \ diff --git a/.github/workflows/test_pr_and_main.yml b/.github/workflows/test_pr_and_main.yml index 357a2fe866e..547e7972aa6 100644 --- a/.github/workflows/test_pr_and_main.yml +++ b/.github/workflows/test_pr_and_main.yml @@ -298,7 +298,7 @@ jobs: || echo "WARNING: Gurobi is not available" python -m pip install --cache-dir cache/pip xpress \ || echo "WARNING: Xpress Community Edition is not available" - python -m pip install git+https://github.com/michaelbynum/wntr.git@working \ + python -m pip install git+https://github.com/usepa/wntr.git@main \ || echo "WARNING: WNTR is not available" fi python -c 'import sys; print("PYTHON_EXE=%s" \ From 15e85a4df798cc5665425f864072f28c5ffdbec5 Mon Sep 17 00:00:00 2001 From: Michael Bynum Date: Mon, 27 Nov 2023 09:03:24 -0700 Subject: [PATCH 8/9] update wntr install for tests --- .github/workflows/test_branches.yml | 2 +- .github/workflows/test_pr_and_main.yml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/test_branches.yml b/.github/workflows/test_branches.yml index 69c835f7c34..6faa7f167c9 100644 --- a/.github/workflows/test_branches.yml +++ b/.github/workflows/test_branches.yml @@ -258,7 +258,7 @@ jobs: || echo "WARNING: Gurobi is not available" python -m pip install --cache-dir cache/pip xpress \ || echo "WARNING: Xpress Community Edition is not available" - python -m pip install git+https://github.com/usepa/wntr.git@main \ + python -m pip install wntr \ || echo "WARNING: WNTR is not available" fi python -c 'import sys; print("PYTHON_EXE=%s" \ diff --git a/.github/workflows/test_pr_and_main.yml b/.github/workflows/test_pr_and_main.yml index b771c314ac2..3e5ca2b3110 100644 --- a/.github/workflows/test_pr_and_main.yml +++ b/.github/workflows/test_pr_and_main.yml @@ -288,7 +288,7 @@ jobs: || echo "WARNING: Gurobi is not available" python -m pip install --cache-dir cache/pip xpress \ || echo "WARNING: Xpress Community Edition is not available" - python -m pip install git+https://github.com/usepa/wntr.git@main \ + python -m pip install wntr \ || echo "WARNING: WNTR is not available" fi python -c 'import sys; print("PYTHON_EXE=%s" \ From b474f2f9f1c22fd152bebe0fe6130f095dbc4fb3 Mon Sep 17 00:00:00 2001 From: Michael Bynum Date: Mon, 27 Nov 2023 12:02:45 -0700 Subject: [PATCH 9/9] update workflows --- .github/workflows/test_branches.yml | 8 ++++++-- .github/workflows/test_pr_and_main.yml | 8 ++++++-- 2 files changed, 12 insertions(+), 4 deletions(-) diff --git a/.github/workflows/test_branches.yml b/.github/workflows/test_branches.yml index 6faa7f167c9..529299bc73f 100644 --- a/.github/workflows/test_branches.yml +++ b/.github/workflows/test_branches.yml @@ -258,8 +258,12 @@ jobs: || echo "WARNING: Gurobi is not available" python -m pip install --cache-dir cache/pip xpress \ || echo "WARNING: Xpress Community Edition is not available" - python -m pip install wntr \ - || echo "WARNING: WNTR is not available" + if [[ ${{matrix.python}} == pypy* ]]; then + echo "skipping wntr for pypy" + else + python -m pip install wntr \ + || echo "WARNING: WNTR is not available" + fi fi python -c 'import sys; print("PYTHON_EXE=%s" \ % (sys.executable,))' >> $GITHUB_ENV diff --git a/.github/workflows/test_pr_and_main.yml b/.github/workflows/test_pr_and_main.yml index 3e5ca2b3110..36c9c45c6a4 100644 --- a/.github/workflows/test_pr_and_main.yml +++ b/.github/workflows/test_pr_and_main.yml @@ -288,8 +288,12 @@ jobs: || echo "WARNING: Gurobi is not available" python -m pip install --cache-dir cache/pip xpress \ || echo "WARNING: Xpress Community Edition is not available" - python -m pip install wntr \ - || echo "WARNING: WNTR is not available" + if [[ ${{matrix.python}} == pypy* ]]; then + echo "skipping wntr for pypy" + else + python -m pip install wntr \ + || echo "WARNING: WNTR is not available" + fi fi python -c 'import sys; print("PYTHON_EXE=%s" \ % (sys.executable,))' >> $GITHUB_ENV