Complex numbers in expression trees

[jsiirola]

Working on a PR #872 highlighted a rather nefarious situation in the Pyomo expression system. Python is rather inconsistent in how it treats complex numbers. For example, in Python 2.x:

>>> import math
>>> import cmath
>>> cmath.sqrt(-1)
1j
>>> math.sqrt(-1)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ValueError: math domain error
>>> (-1)**(.5)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ValueError: negative number cannot be raised to a fractional power

but in Python 3, we see:

>>> import math
>>> import cmath
>>> cmath.sqrt(-1)
1j
>>> math.sqrt(-1)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ValueError: math domain error
>>> (-1)**(.5)
(6.123233995736766e-17+1j)

This is problematic in a number of ways. First, in some versions of Python, the pow() operator will raise an exception and in others it will happily return a complex number:

>>> m = ConcreteModel()
>>> m.p = Param(initialize=-1, mutable=True)
>>> m.e = Expression(expr = m.p**0.5)

# Python2.x:
>>> value(m.e)
ERROR: evaluating object as numeric value: e
        (object: <class 'pyomo.core.base.expression.SimpleExpression'>)
    negative number cannot be raised to a fractional power

# Python 3.x:
>>> value(m.e)
(6.123233995736766e-17+1j)

This is further confounded by our use of the math library to evaluate unary operators. For ALL versions of Python, the math library operations will raise exceptions when they encounter complex numbers, so the following will raise exceptions in all Python versions; however, the exception that is raised will be different:

>>> m.f = Expression(expr=exp(m.e))

# Python 2:
>>> value(m.f)
ERROR: evaluating object as numeric value: f
        (object: <class 'pyomo.core.base.expression.SimpleExpression'>)
    negative number cannot be raised to a fractional power
    [...]
ValueError: negative number cannot be raised to a fractional power

# Python 3:
>>> value(m.f)
ERROR: evaluating object as numeric value: f
        (object: <class 'pyomo.core.base.expression.SimpleExpression'>)
    can't convert complex to float
    [...]
TypeError: can't convert complex to float

Finally, in all versions of Python, complex numbers are not comparable, and and relative comparison (<,<=) will raise a TypeError.

So, what should we do? I can see several options:

1. Disallow all complex numbers in Pyomo. We can trap for complex numbers at key locations (i.e., in pow(), and Var/Param set_value()) and raise an exception immediately. For consistency with Python2, we should probably raise a ValueError.

   • Pro: No (current) Pyomo solver supports complex variables, so we shouldn't incur the pain / technical debt of them either.
   • Pro: Minimal change to the current codebase.
   • Con: this will explicitly preclude support for complex numbers in Pyomo (and any future solver that can support complex numbers)

2. We could allow (more) use of complex numbers in the expression system. All of our unary operations have equivalents in the cmath library that we could (quietly) substitute in, except for floor and ceil. Evaluating expressions that had invalid operators in them (inequality, floor, ceil) would continue to raise a TypeError.

   • Pro: this would make Pyomo more flexible, and would make our expression system behave more like Python's.
   • Pro: this would allow for future solvers that could explicitly deal with complex numbers.
   • Pro: this is a simple change (just changing where we import the unary functions from); BUT:
   • Con: This could open up a testing nightmare (although I would argue that it already exists and we have just not triggered it before now). The bulk of our codebase assumes that all variables / parameters are real values. We could be opening ourselves to another lifetime of tracking down exceptions where complex numbers appeared where we had never thought they might.

[michaelbynum]

This is indeed a complex issue. If we went with option 1, and later decided we needed to support complex numbers, would that be a problem? It seems much easier to add support for complex numbers later than to remove support for complex numbers after attempting to support them.

[carldlaird]

I believe we should go with option 1. In most cases, obtaining a complex number during evaluation is the result of a numerical mistake (and I would want to know as soon as possible, rather than carrying around a complex number to cause an error of a different sort much later). Also, since this hasn't been brought up before, I don't believe we interface to any solvers that support this (or have any users of such a solver if we do). I say we throw an error early and add support later if a good reason emerges.

[whipstein]

Actually, this is posing a big problem for me currently. It may just be my lack of experience with using Pyomo. I did just start using it after having issues getting satisfactory results with other optimization libraries, such as scipy.

I am attempting to perform some complex number linear algebra to combine Y & Z parameter matrices of circuit elements to try to extract model values of devices. Because converting between Y & Z parameters involve matrix inversion, it is anything but straightforward to pre-compute the real and imaginary components. Since Pyomo does not support this, it puts me in a big bind. I have to figure out a way to precompute some very complex matrices. Also, I cannot just define a complete matrix (model) and pull the value(s) I need to send to the optimizer, since these optimizations involve multiple steps.

[mrmundt]

@jsiirola - At least part of this is resolved (in that we don't support Python 2 anymore), but do you know if this issue is otherwise still relevant/valid?