Modeling surprises

Here is an example where the PuLP modeling tool goes berserk.

In standard linear programming, only \(\ge\), \(=\) and \(\le\) constraints are supported. Some tools also allow \(\ne\), which for MIP models needs to be reformulated into a disjunctive constraint. Here is an attempt to do this in PuLP [1]. PuLP does not support this relational operator in its constraints, so we would expect a meaningful error message.

#-------------------------------------------
# funny behavior when using a != constraint
# in a PuLP model
#-------------------------------------------

import pulp 

model = pulp.LpProblem("notequal")
x = pulp.LpVariable("x",lowBound=0,cat=pulp.LpInteger)
y = pulp.LpVariable("y",lowBound=0,cat=pulp.LpInteger)

model += x+y     # obj
model += pulp.LpConstraint(x != y) 

print(model)
# this is an alternative debugging tool: 
#    model.writeLP("/tmp/notequal.lp")

# note that
#    model += x != y 
# produces a different model!!

This prints:

notequal:
MINIMIZE
1*x + 1*y + 0
SUBJECT TO
_C1:0 = -1

VARIABLES
0 <= x Integer
0 <= y Integer

Obviously, this is a surprise: the generated LP is just nonsensical. Surprises are bad when it comes to programming. We want reliability and predictability. It goes without saying that we are not getting that here.

The underlying issue is that operator overloading (redefining what an operator does for certain arguments) is very difficult to make bulletproof. Another example of this problem is shown in [2]. 

If a modeling tool wants to implement \(\color{darkred}x\ne \color{darkred}y\) for integer variables  \(\color{darkred}x,\color{darkred}y\) , it would need to generate something like : \(\color{darkred}x \le \color{darkred}y-1 \mathbf{\ or\ } \color{darkred}x \ge \color{darkred}y+1\). In a MIP model that could look like \[\begin{align}&\color{darkred}x\le\color{darkred}y-1+\color{darkblue}M_1\cdot\color{darkred}\delta\\&\color{darkred}x\ge\color{darkred}y+1-\color{darkblue}M_2\cdot(1-\color{darkred}\delta)\\&\color{darkred}\delta\in\{0,1\} \end{align}\] where \(\color{darkblue}M_1\), \(\color{darkblue}M_2\) are judiciously chosen large constants. If \(\color{darkred}x \in \{0,1,\dots,\color{darkblue}U_x\}\), \(\color{darkred}y \in \{0,1,\dots,\color{darkblue}U_y\}\), then the tightest values would be \(\color{darkblue}M_1 := 1+\color{darkblue}U_x\), \(\color{darkblue}M_2 := 1+\color{darkblue}U_y\). If no good upperbounds \(\color{darkblue}U\) are available, we can use either indicator constraints \[\begin{align}&\color{darkred}\delta=0\implies\color{darkred}x\le\color{darkred}y-1\\&\color{darkred}\delta=1\implies \color{darkred}x\ge \color{darkred}y+1\\&\color{darkred}\delta \in \{0,1\}\end{align}\] or SOS1 variables: \[\begin{align}&\color{darkred}x\le\color{darkred}y-1+\color{darkred}s_1\\& \color{darkred}x\ge \color{darkred}y+1 - \color{darkred}s_2 \\ & \color{darkred}s_1,\color{darkred}s_2\ge 0 \\& \{\color{darkred}s_1,\color{darkred}s_2\} \in \mathbf{SOS1} \end{align}\] depending on the capabilities of the modeling tools and solvers that are used.

Example 1: Minizinc

The modeling language Minizinc can do the binary-variable transformation automatically when bounds are available. My example model looks like:

var 0..5: x;
var 0..7: y;

constraint x != y;

solve minimize x+y;

It will generate a MIP model like:

\ENCODING=ISO-8859-1
\Problem name: MIPCplexWrapper

Minimize
 obj1: X_INTRODUCED_0_
Subject To
 p_lin_0: x - y + 6 X_INTRODUCED_7_ <= 5
 p_lin_1: - x + y - 8 X_INTRODUCED_7_ <= -1
 p_lin_2: x + y - X_INTRODUCED_0_  = 0
Bounds
 0 <= x <= 5
 0 <= y <= 7
 0 <= X_INTRODUCED_0_ <= 12
 0 <= X_INTRODUCED_7_ <= 1
Generals
 x  y  X_INTRODUCED_0_  X_INTRODUCED_7_
End

Let's rewrite this a little bit to make it more readable. We can recognize:

\[\begin{align}\min\>&\color{darkred}z \\ & \color{darkred}x \le \color{darkred}y - 1 + \color{darkblue}M_1 \cdot (1-\color{darkred}\delta) \\ &\color{darkred}x \ge \color{darkred}y + 1 - \color{darkblue}M_2 \cdot \color{darkred}\delta \\ & \color{darkred}z = \color{darkred}x+\color{darkred}y\end{align}\]  

Notice that the values of the big-M's are the same as we calculated before: \(\color{darkblue}M_1=6\) and \(\color{darkblue}M_2=8\).

Minizinc does not employ my suggestions when no upper bounds are available. In that case, it stops with an error message:

MiniZinc: assertion failed: Variable X_INTRODUCED_10_ needs finite upper bound for a big-M constraint, current domain -infinity..infinity

The message is a bit problematic: we don't have a X_INTRODUCED_10_ in our model. It would be better to point back to the cause of the problem: the variables x and y.

Example 2: GAMS

Even expensive commercial tools can show rather embarrassing behavior. Here, we use the indicator constraint approach discussed above. Indicator constraints are very poorly implemented in GAMS. They have to be specified in a solver option file. Obviously, that is very, very wrong. Using an option file should not change the meaning of a model. There are more problems, however. This is a good example of where things go wrong with an error message that is likely confusing.

$onText     Illustrates a problem with indicator constraints in GAMS   $offText   integer variable x,y; variable z 'objective'; binary variable delta;   Equations    Obj    e1    e2 ;   Obj.. z =e= x+y; e1..  x =l= y - 1; e2..  x =g= y + 1;   $onEcho > cplex.opt indic e1$delta 1 indic e2$delta 0 $offEcho   model m /all/; * activate option file m.optfile=1; solve m minimizing z using mip;

 

The error message looks like:

--- GMO setup time: 0.00s

Error: Column is not of type Variable

Option record: indic e1$delta 1

*** GMO ERROR: Failed calling GMO: get indicator map

*** Contact support@gams.com for help.

The underlying problem is that the variable delta is not part of the generated model. The workaround is to add a dummy constraint that contains delta (e.g. \(\color{darkred}\delta \ge 0\)).

Conclusion

In some sense, we are still in the stone age of computing. We want predictability in the behavior of optimization tools and meaningful error messages if we do something wrong.

References

1. Unique values constraint in PuLP, https://or.stackexchange.com/questions/12004/unique-values-constraint-in-pulp
2. PuLP mystery, https://yetanothermathprogrammingconsultant.blogspot.com/2020/10/pulp-mystery.html