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;
\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
\[\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.00sError: Column is not of type VariableOption 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. Wanting proper feedback in the form of meaningful error messages is not asking too much.
References
- Unique values constraint in PuLP, https://or.stackexchange.com/questions/12004/unique-values-constraint-in-pulp
- PuLP mystery, https://yetanothermathprogrammingconsultant.blogspot.com/2020/10/pulp-mystery.html
No comments:
Post a Comment