PuLP mystery

PuLP is a popular Python-based modeling tool for LP and MIP models. In [1], a user asked a (somewhat trivial) question about PuLP syntax. But there was an interesting wrinkle in the model that was unfamiliar to me. The basic issue is that sometimes PuLP accepts a NumPy array in expressions:

 

import pulp as p
import numpy as np

a=np.array([1,2,3])

x=p.LpVariable('x')
y=p.LpVariable('y')

prob = p.LpProblem("test",p.LpMaximize)
prob += x+(a-y)
prob

This fragment creates an objective. But we have a strange term here. a-y is funny because a is a NumPy array (vector) and y is a scalar PuLP variable. Usually, adding a scalar to a vector is interpreted as elementwise addition: \[\begin{pmatrix} a_0 - y \\ a_1 - y \\ a_2-y \end{pmatrix}\] How one would interpret an objective like: \[ x + \begin{pmatrix} a_0 - y \\ a_1 - y \\ a_2-y \end{pmatrix}\] is not clear to me. I would say this is an error. However, PuLP accepts this, and prints the model as:

test:
MAXIMIZE
1*x + -3*y + 6
VARIABLES
x free Continuous
y free Continuous

So, apparently the interpretation is: \[ x + \sum_i (a_i-y)=x + \sum_i a_i - n\cdot y\] where \(n\) is the length of vector \(a\). But then again, we would expect PuLP to accept

 

  prob += (a-y)

as objective. However, this produces the error message: 

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-6-bf3a3a41e9da> in <module>()
      8 
      9 prob = p.LpProblem("test",p.LpMaximize)
---> 10 prob += (a-y)
     11 prob

/usr/local/lib/python3.6/dist-packages/pulp/pulp.py in __iadd__(self, other)
   1528             self.objective.name = name
   1529         else:
-> 1530             raise TypeError("Can only add LpConstraintVar, LpConstraint, LpAffineExpression or True objects")
   1531         return self
   1532 

TypeError: Can only add LpConstraintVar, LpConstraint, LpAffineExpression or True objects

I have no idea what is going on here. It could be all just a bug, but even that can not easily explain the behavior of sometimes accepting (a-y) and sometimes refusing it.

The underlying problem can actually be traced back to operator overloading as shown by [2]. Who owns the + or -? The following experiment demonstrates the issue in stark terms:

>>> a+x
array([1*x + 1, 1*x + 2, 1*x + 3], dtype=object)
>>> x+a
1*x + 6

In the first case, NumPy handles the +, leading to the interpretation: \[\begin{pmatrix} a_0  \\ a_1  \\ a_2 \end{pmatrix} + x = \begin{pmatrix} x+a_0  \\ x+a_1  \\x+a_2  \end{pmatrix}\] In the second case, the + is dealt with by PuLP. This gives \[x+\begin{pmatrix} a_0  \\ a_1  \\ a_2 \end{pmatrix}=x+\sum_i a_i\]  

Conclusions

This is mind-bending (a.k.a. crazy). It would be better if things were a bit more predictable ("orthogonal" is the term used in programming language design). Preferably, we don't want to spend time on figuring out how to interpret a + or a -. 

Advise: it is better not to mix PuLP with NumPy.

References

1. 'LpAffineExpression' object has no attribute 'solve' in pulp Lp problem- Optimization, https://stackoverflow.com/questions/64405352/lpaffineexpression-object-has-no-attribute-solve-in-pulp-lp-problem-optimiz
2. Ayhan Fuat Çelik, PuLP and NumPy, https://gist.github.com/ayhanfuat/5c72f6cec33fa6d50ae238d8f4345f44
3. The picture is from https://pulpcovers.com/tag/thrillingmystery/