Most LP/MIP modeling tools stay close to what is the underpinning of a MixedInteger Programming problem: a system of linear equations (equalities and inequalities) plus a linear objective. Examples are Pulp and GAMS. In Constraint Programming (CP), modeling tools need to provide access to higherlevel global constraints. Without these global constraints (such as the famous alldifferent constraint), CP solvers would not perform very well.
Problem
Highlevel Constraint 

\[\begin{align}& \color{darkred}x_i = \left \{ j: \color{darkred}x_j=i\} \right \\ & \color{darkred}x_i \in \{0,1,2,\dots\} \end{align}\] 
x = [2,1,2,0,0]
Traditional MIP formulation
Linear MIP constraints 

\[\begin{align}& \sum_j j\cdot\color{darkred}y_{i,j} = \sum_j \color{darkred}y_{j,i} && \forall i \\ & \sum_j \color{darkred}y_{i,j} = 1 && \forall i \\ & \color{darkred}y_{i,j}\in \{0,1\} \end{align}\] 
Constraint Programming Environments
Models have two different audiences. The first is obviously the modeling tool and solver. The second, at least as important, is the human reader. I think the mathematical notation using the sum is way easier to understand. Obviously, the solver could analyze the model and replace the summation with the cardinality constraint. We see this more and more in highend tools and MIP solvers. The solver tries to understand the model structure and may apply automatic reformulations. Another argument for this is that there are just way too many global constraints: it is left to the modeler to remember them all. I count 121 different global constraints listed in the Minizinc documentation.
Automatic reformulations
Nerdy side note: if you use sum anyway, you will get a syntax error with the message "syntax error." Such poor error messages are sometimes the result of using certain parser generators. Many compilers choose to implement hardcrafted parsers to improve error messaging.
If we use a MIP solver, AMPL cannot just pass this on to the solver. Rather it will need to apply some reformulation. Interestingly, what AMPL generates, depends on the exact declaration of the integer variable. If we specify \(\color{darkred}x_i \in \{0,1,\dots\}\), then AMPL will generate something like: \[\begin{align}&\color{darkred}x_i = \sum_j \color{darkred}\delta^=_{i,j} \\& \color{darkred}\delta^{\lt}_{i,j} + \color{darkred}\delta^=_{i,j} + \color{darkred}\delta^{\gt}_{i,j} \ge 1 \\ & \color{darkred}\delta^{\lt}_{i,j}=1 \Rightarrow \color{darkred}x_j \le i1 \\& \color{darkred}\delta^{=}_{i,j}=1 \Rightarrow \color{darkred}x_j = i \\ & \color{darkred}\delta^{\gt}_{i,j}=1 \Rightarrow \color{darkred}x_j \ge i+1 \\ & \color{darkred}\delta^{\lt}_{i,j},\color{darkred}\delta^=_{i,j},\color{darkred}\delta^{\gt}_{i,j} \in \{0,1\}\end{align}\] I think I would use: \( \color{darkred}\delta^{\lt}_{i,j} + \color{darkred}\delta^=_{i,j} + \color{darkred}\delta^{\gt}_{i,j} = 1\) but the \(\ge 1\) is also fine as the indicator constraints make sure only one of them can be 1.
However, when we declare: \(\color{darkred}x_i \in \{0,\dots,n\}\), then AMPL will generate a similar thing as our linear MIP constraint. That is actually pretty cool. In this case, the difference in the generated MIP model does not really matter (Gurobi solves both cases in the presolve). But, in general, it may be better to be as tight as possible. So instead of
var x{I} integer >= 0; (or even var x{I} integer;)
use
var x{I} in I;
The generated MIP model will be different and most likely more efficient if we use the latter declaration.
This is a good example where it can help to look a bit under the hood. Of course, this leads to the observation: these transformations are somewhat of a black box, and it is not always easy to assess that the modeler provided the best information. A slight difference in some detail in the AMPL code, may cause very different models to be generated. Without letting the modeler know. When specifying the linear MIP model with binary variables myself, I know better what the solver will receive.
One other remark. I noticed that the MIP model generated by AMPL relies on further presolving by the solver. AMPL also has its own presolver. I think the translation is done after AMPLs own presolving.
Not all solvers support indicator constraints. I have not checked what AMPL does in the case of missing support for that. I assume some bigM constraints will be generated. It would be interesting to see what values for the bigM constants are used. (Solvers are very sensitive to that).
To see what AMPL generated, I let Gurobi write an LP file. I don't know if there is a better way to see what AMPL did. Digesting the LP file is a bit difficult. I don't know if there is a better way to get some feedback.
Conclusion
 Improve readability,
 Novice modelers can implement constructs that were otherwise beyond their reach,
 Experienced modelers may use it for fast prototyping.
References
 A strange series, http://yetanothermathprogrammingconsultant.blogspot.com/2022/11/astrangeseries.html
 Global_cardinality, https://www.minizinc.org/doc2.6.4/en/libglobalscounting.html#mznrefglobalscountingglobalcardinality
Appendix: models
OPL CP Optimizer
using CP;
int n = 100; range I = 0..n; // dvar int x[I]; dvar int x[I] in I; // slight performance difference
constraints { forall (i in I) x[i] == sum(j in I) (x[j]==i); }
Minizinc/Gecode
int: n = 4; set of int: I = 0..n;
array[I] of var I: x;
% Constraint variant 1: Individual sums for each value constraint forall (i in I) (x[i] = sum(j in I) (x[j] = i));
%Constraint variant 2: Global cardinality constraint %include "globals.mzn"; %array[I] of I: index = array1d(I, [i  i in I]); %constraint global_cardinality(x, index, x);
solve satisfy;
output [show(x)];
AMPL
param n integer = 4; set I = {0..n};
var x{I} integer; #var x{I} in I;
subject to con{i in I}: x[i] = count{j in I} (x[j] = i);
option solver gurobi; option gurobi_options 'writeprob=x.lp outlev=1'; solve;

CVX, and even CVXPY, are "primitive" compared to what they could be. http://ask.cvxr.com/t/phdthesisideanextgenerationofcvxorcvxpywhichautomaticallyfiguesoutreformulationsusingcomputeralgebratechniques/10515
ReplyDeleteInteresting post! Coming from a CP background, I generally have a much easier timer reading models that use global constraints. While there are many different ones proposed in the literature, in most cases there is only a handful of global constraints that are used in most models.
ReplyDeleteI typically compare the use of global constraints in modelling to programming with a standard library. Sure, I could implement a sorting routine every time a write a program that needs to sort values. But it is easier to read the code and be confident that it does what I think it does if the most common operations are named entities in the library.
Sometimes the library code will be good and sometimes it will be bad. The same holds true for a global constraint, especially if it is decomposed for a solver that does not support it. I prefer using the global constraint first as a convenient shorthand, and only when a model needs to be fast will I look into if there is anything smarter that could be done. I think that systems like MiniZinc with the focus on translating into lowlevel constraints have a lot of potential in making this a structured process, such as automating trying out different decompositions and comparing for different solving technologies.
For the problem in question, when I read the original description it naturally read as a global cardinality problem to me (counting occurrences of multiple values). Reading your original model, the best way for me to understand the constraints was to mentally map it to the corresponding global cardinality, with the slightly nonstandard construction of the same variables doing both the counting and being counted. Given that, I really just wanted to try it out in MiniZinc to see if it made any difference in solving performance for a CP solver.
I find it interesting just how much the modelling tradition one is used to influences how we think about models!