Logical conditions

> How do I model
>   if u(i,j)=0 then x(i,j)=x(i,j-1)+1 & y(i,j)=0
>   else y(i,j)=y(i,j-1)+1 & x(i,j)=0

With big-M formulations:

• -M₁*u(i,j) ≤ x(i,j) − x(i,j−1) − 1 ≤ +M₁*u(i,j)
• y(i,j) ≤ M₂*u(i,j)
• -M₃*(1−u(i,j)) ≤ y(i,j)−y(i,j−1)−1 ≤ +M₃*(1−u(i,j))
• x(i,j) ≤ M₄*(1−u(i,j))

The constants M₁ through M₄ need to be chosen with care: as small as possible but large enough that they can make the corresponding constraint non-binding. Some conservative choices would be (assuming the integer variables have a lower-bound of zero):

• M₁=x.up+1
• M₂=y.up
• M₃=y.up+1
• M₄=x.up

where x.up and y.up are the upper-bounds on x and y (chosen as tight as possible).

With Cplex indicator constraints this can be modeled even simpler, as you don’t have to worry about big-M values. GAMS has limited but quite workable support for this through an option file; AMPL has more complete support for these directly in their language which is a cleaner approach IMHO. Implications are integral part of the model, and moving them to an option file makes the model less self-contained: you have to look in different spots to understand the model. I would prefer that options are reserved to solver options (tolerances, algorithm selection etc.) instead of modeling constructs.

The argument against augmenting the language is that this is a pure Cplex-only facility (although I believe SCIP also supports indicator constraints), that cannot be translated transparently to other MIP solvers. (In the long term one would expect constraint programming with their if-then-else constructs will get more accepted either separate or integrated into MIP solvers; modeling languages will need then to manage these constructs in a consistent way. AMPL had some efforts in this direction, although I don’t think their CP facilities are general available, and MS OML already supports both paradigms.)

The above could be modeled in AMPL as:

  e1{i in I, j in I}: u[i,j]=0 ==> y[i,j]=0 else x[i,j]=0;
  e2{i in I, j in I: j>1}: u[i,j]=0 ==> x[i,j]=x[i,j−1]+1 else y[i,j]=y[i,j−1]+1;

(see also this post).