diff --git a/pyomo/core/expr/compare.py b/pyomo/core/expr/compare.py index f75a5490e73..8adb16172f1 100644 --- a/pyomo/core/expr/compare.py +++ b/pyomo/core/expr/compare.py @@ -16,7 +16,7 @@ NPV_ProductExpression, NPV_DivisionExpression, NPV_PowExpression, NPV_SumExpression, NPV_NegationExpression, NPV_UnaryFunctionExpression, NPV_ExternalFunctionExpression, Expr_ifExpression, AbsExpression, - NPV_AbsExpression, NumericValue) + NPV_AbsExpression, NumericValue, QuadraticExpression) from pyomo.core.expr.logical_expr import ( RangedExpression, InequalityExpression, EqualityExpression ) @@ -32,6 +32,14 @@ def handle_linear_expression(node: LinearExpression, pn: List): return tuple() +def handle_quadratic_expression(node: QuadraticExpression, pn: List): + pn.append((type(node), 3*node.nargs())) + pn.extend(node.coefs) + pn.extend(node.vars_1) + pn.extend(node.vars_2) + return tuple() + + def handle_expression(node: ExpressionBase, pn: List): pn.append((type(node), node.nargs())) return node.args @@ -76,6 +84,7 @@ def handle_external_function_expression(node: ExternalFunctionExpression, pn: Li handler[RangedExpression] = handle_expression handler[InequalityExpression] = handle_expression handler[EqualityExpression] = handle_expression +handler[QuadraticExpression] = handle_quadratic_expression class PrefixVisitor(StreamBasedExpressionVisitor): diff --git a/pyomo/core/expr/current.py b/pyomo/core/expr/current.py index ae9f2407ce5..dfc74649726 100644 --- a/pyomo/core/expr/current.py +++ b/pyomo/core/expr/current.py @@ -66,6 +66,7 @@ class Mode(object): NPV_UnaryFunctionExpression, AbsExpression, NPV_AbsExpression, LinearExpression, + QuadraticExpression, _MutableLinearExpression, decompose_term, LinearDecompositionError, diff --git a/pyomo/core/expr/numeric_expr.py b/pyomo/core/expr/numeric_expr.py index cafd1f33cb2..d36c51841a1 100644 --- a/pyomo/core/expr/numeric_expr.py +++ b/pyomo/core/expr/numeric_expr.py @@ -1518,6 +1518,160 @@ def _combine_expr(self, etype, _other): return self +class QuadraticExpression(ExpressionBase): + """ + An expression object for quadratic expressions + + Args: + args: tuple or list + The children nodes + """ + __slots__ = ('_coefs', '_vars_1', '_vars_2', '_args_cache_') + + PRECEDENCE = 6 + + def __init__(self, args=None, coefs=None, vars_1=None, vars_2=None): + """ + A quadratic expression of the form sum_i(c_i*x_i*y_i). + + You can specify args OR (coefs, vars_1, and vars_2). If args is provided, it + should be a list or tuple that contains a series of + :py:class:`ProductExpression` objects with arguments + :py:class:`MonomialTermExpression` and :py:class:`_GeneralVarData`. + Alternatively, you can specify the lists of coefficients and variables + separately. + """ + if args: + if any(arg is not None for arg in (coefs, vars_1, vars_2)): + raise ValueError('Cannot specify both args and any of ' + '{coefs, vars_1, vars_2}') + self._args_ = args + else: + self._coefs = tuple(coefs) + self._vars_1 = tuple(vars_1) + self._vars_2 = tuple(vars_2) + self._args_cache_ = tuple() + + if len(self.vars_1) != self.nargs() or len(self.vars_2) != self.nargs(): + raise ValueError('The length of coefs, vars_1, and ' + 'vars_2 must be the same.') + + @property + def coefs(self): + return self._coefs + + @property + def vars_1(self): + return self._vars_1 + + @property + def vars_2(self): + return self._vars_2 + + def nargs(self): + return len(self.coefs) + + @property + def _args_(self): + nargs = self.nargs() + if len(self._args_cache_) != nargs: + tmp = map(MonomialTermExpression, zip(self.coefs, self.vars_1)) + tmp = map(ProductExpression, zip(tmp, self.vars_2)) + self._args_cache_ = tuple(tmp) + return self._args_cache_ + + @_args_.setter + def _args_(self, val): + self._args_cache_ = tuple(val) + coefs = list() + vars_1 = list() + vars_2 = list() + for term in val: + arg1, arg2 = term.args + if (type(arg1) is not MonomialTermExpression + or type(arg2) in native_numeric_types + or not arg2.is_variable_type()): + raise ValueError('When constructing a QuadraticExpression with args, ' + 'the args must all be ProductExpressions containing a ' + 'MonomialTermExpression for the first argument and a ' + 'variable for the second argument.') + c, v1 = arg1.args + v2 = arg2 + coefs.append(c) + vars_1.append(v1) + vars_2.append(v2) + self._coefs = tuple(coefs) + self._vars_1 = tuple(vars_1) + self._vars_2 = tuple(vars_2) + + def _precedence(self): + return QuadraticExpression.PRECEDENCE + + def create_node_with_local_data(self, args, classtype=None): + if classtype is None: + if not args: + classtype = self.__class__ + else: + correct_format = True + for arg in args: + if (type(arg) is not ProductExpression + or type(arg.args[0]) is not MonomialTermExpression + or type(arg.args[1]) in native_numeric_types + or not arg.args[1].is_variable_type()): + correct_format = False + break + if correct_format: + classtype = self.__class__ + else: + classtype = SumExpression + return classtype(args) + + def getname(self, *args, **kwds): + return 'sum' + + def _compute_polynomial_degree(self, result): + res = 0 + for v1, v2 in zip(self.vars_1, self.vars_2): + if v1.fixed: + if not v2.fixed: + res = max(res, 1) + else: + if v2.fixed: + res = max(res, 1) + else: + res = max(res, 2) + return res + + def is_constant(self): + return self.nargs() == 0 + + def _is_fixed(self, values=None): + return all(v.fixed for v in self.vars_1) and all(v.fixed for v in self.vars_2) + + def is_fixed(self): + return self._is_fixed() + + def _to_string(self, values, verbose, smap, compute_values): + if not values: + values = ['0'] + if verbose: + return "%s(%s)" % (self.getname(), ', '.join(values)) + + for i in range(1, len(values)): + term = values[i] + if term[0] not in '+-': + values[i] = '+ ' + term + elif term[1] != ' ': + values[i] = term[0] + ' ' + term[1:] + return ' '.join(values) + + def is_potentially_variable(self): + return self.nargs() > 0 + + def _apply_operation(self, result): + return sum(result) + + class _MutableLinearExpression(LinearExpression): __slots__ = () diff --git a/pyomo/core/tests/unit/test_numeric_expr.py b/pyomo/core/tests/unit/test_numeric_expr.py index a1ad6f7bca2..9e060cf7fd8 100644 --- a/pyomo/core/tests/unit/test_numeric_expr.py +++ b/pyomo/core/tests/unit/test_numeric_expr.py @@ -41,7 +41,7 @@ from pyomo.core.expr.numeric_expr import ( ExpressionBase, UnaryFunctionExpression, SumExpression, PowExpression, ProductExpression, NegationExpression, linear_expression, - MonomialTermExpression, LinearExpression, DivisionExpression, + MonomialTermExpression, LinearExpression, QuadraticExpression, DivisionExpression, NPV_NegationExpression, NPV_ProductExpression, NPV_PowExpression, NPV_DivisionExpression, decompose_term, clone_counter, nonlinear_expression, @@ -57,6 +57,7 @@ from pyomo.core.expr.template_expr import IndexTemplate from pyomo.core.expr import expr_common from pyomo.core.base.var import _GeneralVarData +from pyomo.core.expr.compare import compare_expressions from pyomo.repn import generate_standard_repn from pyomo.core.expr.numvalue import NumericValue @@ -4784,6 +4785,138 @@ def test_pow_other(self): self.assertIs(e.__class__, PowExpression) +class TestQuadraticExpression(unittest.TestCase): + def test_init_without_args(self): + m = ConcreteModel() + m.a = Set(initialize=[0,1,2]) + m.x = Var(m.a) + m.y = Var(m.a) + coefs = [0.1, 20, 3] + e = QuadraticExpression(coefs=coefs, vars_1=m.x.values(), vars_2=m.y.values()) + self.assertEqual(e._args_cache_, tuple()) + self.assertEqual(e.coefs, tuple(coefs)) + self.assertEqual(e.vars_1, (m.x[0], m.x[1], m.x[2])) + self.assertEqual(e.vars_2, (m.y[0], m.y[1], m.y[2])) + self.assertTrue(compare_expressions(e.args[0], 0.1*m.x[0]*m.y[0])) + self.assertTrue(compare_expressions(e.args[1], 20*m.x[1]*m.y[1])) + self.assertTrue(compare_expressions(e.args[2], 3*m.x[2]*m.y[2])) + self.assertEqual(len(e._args_cache_), 3) + + def test_init_with_args(self): + m = ConcreteModel() + m.a = Set(initialize=[0,1,2]) + m.x = Var(m.a) + m.y = Var(m.a) + coefs = [0.1, 20, 3] + e = QuadraticExpression(args=[coefs[i]*m.x[i]*m.y[i] for i in m.a]) + self.assertEqual(len(e._args_cache_), 3) + self.assertEqual(e.coefs, tuple(coefs)) + self.assertEqual(e.vars_1, (m.x[0], m.x[1], m.x[2])) + self.assertEqual(e.vars_2, (m.y[0], m.y[1], m.y[2])) + print(e.args[0]) + self.assertTrue(compare_expressions(e.args[0], 0.1*m.x[0]*m.y[0])) + self.assertTrue(compare_expressions(e.args[1], 20*m.x[1]*m.y[1])) + self.assertTrue(compare_expressions(e.args[2], 3*m.x[2]*m.y[2])) + + def test_init_errors(self): + m = ConcreteModel() + m.a = Set(initialize=[0,1,2]) + m.x = Var(m.a) + m.y = Var(m.a) + coefs = [0.1, 20, 3] + with self.assertRaisesRegex(ValueError, 'Cannot specify both args and any of {coefs, vars_1, vars_2}'): + e = QuadraticExpression(args=[coefs[i]*m.x[i]*m.y[i] for i in m.a], coefs=coefs) + with self.assertRaisesRegex(ValueError, 'Cannot specify both args and any of {coefs, vars_1, vars_2}'): + e = QuadraticExpression(args=[coefs[i]*m.x[i]*m.y[i] for i in m.a], vars_1=m.x.values()) + with self.assertRaisesRegex(ValueError, 'Cannot specify both args and any of {coefs, vars_1, vars_2}'): + e = QuadraticExpression(args=[coefs[i]*m.x[i]*m.y[i] for i in m.a], vars_2=m.x.values()) + with self.assertRaisesRegex(ValueError, 'The length of coefs, vars_1, and vars_2 must be the same.'): + e = QuadraticExpression(coefs=coefs, vars_1=[m.x[0]], vars_2=m.y.values()) + with self.assertRaisesRegex(ValueError, 'The length of coefs, vars_1, and vars_2 must be the same.'): + e = QuadraticExpression(coefs=coefs, vars_1=m.x.values(), vars_2=[m.y[0]]) + with self.assertRaisesRegex(ValueError, 'The length of coefs, vars_1, and vars_2 must be the same.'): + e = QuadraticExpression(coefs=coefs[0:2], vars_1=m.x.values(), vars_2=m.y.values()) + with self.assertRaisesRegex(ValueError, + 'When constructing a QuadraticExpression with ' + 'args, the args must all be ProductExpressions ' + 'containing a MonomialTermExpression for the first ' + 'argument and a variable for the second argument.'): + e = QuadraticExpression(args=[coefs[i]*(m.x[i]*m.y[i]) for i in m.a]) + + def test_to_string(self): + m = ConcreteModel() + m.a = Set(initialize=[0,1,2,3]) + m.x = Var(m.a) + m.y = Var(m.a) + coefs = [1, -1, 2, -2] + e = QuadraticExpression(coefs=[], vars_1=[], vars_2=[]) + self.assertEqual(e.to_string(), "0") + e = QuadraticExpression(coefs=coefs, vars_1=m.x.values(), vars_2=m.y.values()) + self.assertEqual(e.to_string(), "x[0]*y[0] - x[1]*y[1] + 2*x[2]*y[2] - 2*x[3]*y[3]") + + def test_polynomial_degree(self): + m = ConcreteModel() + m.a = Set(initialize=[0,1,2,3]) + m.x = Var(m.a) + m.y = Var(m.a) + coefs = [1, -1, 2, -2] + e = QuadraticExpression(coefs=coefs, vars_1=m.x.values(), vars_2=m.y.values()) + self.assertEqual(e.polynomial_degree(), 2) + m.x.fix(1) + m.y.fix(2) + self.assertEqual(e.polynomial_degree(), 0) + m.x.unfix() + self.assertEqual(e.polynomial_degree(), 1) + m.x.fix() + m.y.unfix() + self.assertEqual(e.polynomial_degree(), 1) + + def test_is_constant_and_fixed(self): + m = ConcreteModel() + m.a = Set(initialize=[0,1,2,3]) + m.x = Var(m.a) + m.y = Var(m.a) + coefs = [1, -1, 2, -2] + e = QuadraticExpression(coefs=coefs, vars_1=m.x.values(), vars_2=m.y.values()) + self.assertFalse(e.is_constant()) + self.assertFalse(e.is_fixed()) + self.assertTrue(e.is_potentially_variable()) + m.x.fix(1) + m.y.fix(2) + self.assertFalse(e.is_constant()) + self.assertTrue(e.is_fixed()) + self.assertTrue(e.is_potentially_variable()) + m.x.unfix() + self.assertFalse(e.is_constant()) + self.assertFalse(e.is_fixed()) + self.assertTrue(e.is_potentially_variable()) + m.x.fix() + m.y.unfix() + self.assertFalse(e.is_constant()) + self.assertFalse(e.is_fixed()) + self.assertTrue(e.is_potentially_variable()) + + e = QuadraticExpression(coefs=[], vars_1=[], vars_2=[]) + self.assertTrue(e.is_constant()) + self.assertTrue(e.is_fixed()) + self.assertFalse(e.is_potentially_variable()) + + def test_apply_operation(self): + m = ConcreteModel() + m.a = Set(initialize=[0,1,2]) + m.x = Var(m.a) + m.y = Var(m.a) + coefs = [0.1, 20, 3] + e = QuadraticExpression(coefs=coefs, vars_1=m.x.values(), vars_2=m.y.values()) + m.x[0].value = 1.2 + m.x[1].value = 2.3 + m.x[2].value = 3.4 + m.y[0].value = -1.5 + m.y[1].value = 2.6 + m.y[2].value = 3.8 + self.assertAlmostEqual(value(e), sum(coefs[i]*m.x[i].value*m.y[i].value for i in m.a)) + + class TestNonlinearExpression(unittest.TestCase): def test_sum_other(self): diff --git a/pyomo/core/tests/unit/test_visitor.py b/pyomo/core/tests/unit/test_visitor.py index 60805f61e7c..7ec91907dca 100644 --- a/pyomo/core/tests/unit/test_visitor.py +++ b/pyomo/core/tests/unit/test_visitor.py @@ -22,7 +22,7 @@ from pyomo.core.expr.numvalue import nonpyomo_leaf_types, NumericConstant from pyomo.core.expr.numeric_expr import ( SumExpression, ProductExpression, - MonomialTermExpression, LinearExpression, + MonomialTermExpression, LinearExpression, QuadraticExpression, NPV_SumExpression, NPV_ProductExpression, NegationExpression, NPV_NegationExpression, PowExpression, NPV_PowExpression, MaxExpression, NPV_MaxExpression, MinExpression, NPV_MinExpression, @@ -366,6 +366,30 @@ def test_replacement_walker0(self): self.assertTrue(compare_expressions(2*LinearExpression(linear_coefs=[i for i in M.z.values()], linear_vars=[i for i in M.w.values()]), f)) + def test_replacement_quadratic(self): + m = ConcreteModel() + m.x = Var(range(3)) + m.w = VarList() + e = QuadraticExpression(coefs=[0.1, 3], vars_1=[m.x[0], m.x[1]], + vars_2=[m.x[2], m.x[2]]) + walker = ReplacementWalkerTest1(m) + f = walker.dfs_postorder_stack(e) + self.assertTrue(compare_expressions( + QuadraticExpression(coefs=[0.1, 3], vars_1=[m.w[1], m.w[3]], + vars_2=[m.w[2], m.w[2]]), f)) + + del m.w + del m.w_index + m.w = VarList() + e = 2 * QuadraticExpression(coefs=[0.1, 3], vars_1=[m.x[0], m.x[1]], + vars_2=[m.x[2], m.x[2]]) + walker = ReplacementWalkerTest1(m) + f = walker.dfs_postorder_stack(e) + self.assertTrue(compare_expressions( + 2 * QuadraticExpression(coefs=[0.1, 3], vars_1=[m.w[1], m.w[3]], + vars_2=[m.w[2], m.w[2]]), f)) + + def test_replace_expressions_with_monomial_term(self): M = ConcreteModel() M.x = Var() @@ -529,6 +553,17 @@ def test_replacement_walker0(self): linear_vars=[i for i in M.x.values()]), e)) self.assertTrue(compare_expressions(2*(M.z[0]*(2*M.w[4]) + M.z[1]*(2*M.w[5]) + M.z[2]*(2*M.w[6])), f)) + def test_replacement_quadratic(self): + m = ConcreteModel() + m.x = Var(range(3)) + m.w = VarList() + e = QuadraticExpression(coefs=[0.1, 3], vars_1=[m.x[0], m.x[1]], + vars_2=[m.x[2], m.x[2]]) + walker = ReplacementWalkerTest2(m) + f = walker.dfs_postorder_stack(e) + self.assertTrue(compare_expressions(0.1*(2*m.w[1])*(2*m.w[2]) + 3*(2*m.w[3])*(2*m.w[2]), f)) + + # # Replace all mutable parameters with variables # @@ -672,6 +707,18 @@ def test_replacement_walker0(self): linear_vars=[i for i in M.x.values()]), e)) self.assertTrue(compare_expressions(2*(2*M.w[4]*M.x[0] + 2*M.w[5]*M.x[1] + 2*M.w[6]*M.x[2]), f)) + def test_replacement_quadratic(self): + m = ConcreteModel() + m.x = Var(range(3)) + m.w = VarList() + m.z = Param(range(2), mutable=True) + e = QuadraticExpression(coefs=[m.z[0], m.z[1]], vars_1=[m.x[0], m.x[1]], + vars_2=[m.x[2], m.x[2]]) + walker = ReplacementWalkerTest3(m) + f = walker.dfs_postorder_stack(e) + self.assertTrue(compare_expressions( + 2 * m.w[1] * m.x[0] * m.x[2] + 2 * m.w[2] * m.x[1] * m.x[2], f)) + # # Replace all mutable parameters with variables diff --git a/pyomo/repn/standard_repn.py b/pyomo/repn/standard_repn.py index ff28e6cb19a..f58ae90fb86 100644 --- a/pyomo/repn/standard_repn.py +++ b/pyomo/repn/standard_repn.py @@ -877,6 +877,7 @@ def _collect_external_fn(exp, multiplier, idMap, compute_values, verbose, quadra EXPR.AbsExpression : _collect_nonl, EXPR.NegationExpression : _collect_negation, EXPR.LinearExpression : _collect_linear, + EXPR.QuadraticExpression : _collect_sum, EXPR.InequalityExpression : _collect_comparison, EXPR.RangedExpression : _collect_comparison, EXPR.EqualityExpression : _collect_comparison, diff --git a/pyomo/repn/tests/test_standard.py b/pyomo/repn/tests/test_standard.py index 011e854d079..4c6707db4d4 100644 --- a/pyomo/repn/tests/test_standard.py +++ b/pyomo/repn/tests/test_standard.py @@ -24,6 +24,8 @@ from pyomo.environ import AbstractModel, ConcreteModel, Var, Param, Set, Expression, RangeSet, ExternalFunction, quicksum, cos, sin, summation, sum_product import pyomo.kernel from pyomo.core.expr.numvalue import native_numeric_types, as_numeric, value +from pyomo.core.expr.numeric_expr import QuadraticExpression +from pyomo.core.expr.compare import compare_expressions class frozendict(dict): @@ -2696,6 +2698,65 @@ def test_quadratic2(self): rep = pickle.loads(s) self.assertEqual(baseline1, repn_to_dict(rep)) + def test_quadratic3(self): + m = ConcreteModel() + m.a = Set(initialize=[0,1,2]) + m.x = Var(m.a) + m.y = Var(m.a) + coefs = [0.1, 20, 3] + e = QuadraticExpression(coefs=coefs, vars_1=m.x.values(), vars_2=m.y.values()) + + rep = generate_standard_repn(e) + self.assertFalse(rep.is_fixed()) + self.assertEqual(rep.polynomial_degree(), 2) + self.assertFalse(rep.is_constant()) + self.assertFalse(rep.is_linear()) + self.assertTrue(rep.is_quadratic()) + self.assertTrue(rep.is_nonlinear()) + self.assertEqual(rep.constant, 0) + self.assertEqual(len(rep.linear_vars), 0) + self.assertEqual(len(rep.linear_coefs), 0) + self.assertEqual(len(rep.quadratic_vars), 3) + self.assertEqual(len(rep.quadratic_coefs), 3) + self.assertIsNone(rep.nonlinear_expr) + self.assertEqual(len(rep.nonlinear_vars), 0) + self.assertEqual(rep.quadratic_coefs, tuple(coefs)) + expected = ((m.x[0], m.y[0]), (m.x[1], m.y[1]), (m.x[2], m.y[2])) + self.assertEqual(rep.quadratic_vars, expected) + + rep = generate_standard_repn(e, quadratic=False) + self.assertFalse(rep.is_fixed()) + self.assertEqual(rep.polynomial_degree(), None) + self.assertFalse(rep.is_constant()) + self.assertFalse(rep.is_linear()) + self.assertFalse(rep.is_quadratic()) + self.assertTrue(rep.is_nonlinear()) + self.assertEqual(rep.constant, 0) + self.assertEqual(len(rep.linear_vars), 0) + self.assertEqual(len(rep.linear_coefs), 0) + self.assertEqual(len(rep.quadratic_vars), 0) + self.assertEqual(len(rep.quadratic_coefs), 0) + self.assertEqual(rep.nonlinear_vars, + (m.y[0], m.y[1], m.y[2], m.x[0], m.x[1], m.x[2])) + self.assertTrue(compare_expressions(rep.nonlinear_expr, sum( + coefs[i] * m.x[i] * m.y[i] for i in m.a))) + + m.x.fix(2) + rep = generate_standard_repn(e) + self.assertFalse(rep.is_fixed()) + self.assertEqual(rep.polynomial_degree(), 1) + self.assertFalse(rep.is_constant()) + self.assertTrue(rep.is_linear()) + self.assertFalse(rep.is_quadratic()) + self.assertFalse(rep.is_nonlinear()) + self.assertEqual(rep.constant, 0) + self.assertEqual(rep.linear_vars, (m.y[0], m.y[1], m.y[2])) + self.assertEqual(rep.linear_coefs, (0.2, 40, 6)) + self.assertEqual(len(rep.quadratic_vars), 0) + self.assertEqual(len(rep.quadratic_coefs), 0) + self.assertIsNone(rep.nonlinear_expr) + self.assertEqual(len(rep.nonlinear_vars), 0) + def test_pow(self): # ^ # / \