> 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
_{1}*u(i,j) ≤ x(i,j) − x(i,j−1) − 1 ≤ +M_{1}*u(i,j) - y(i,j) ≤ M
_{2}*u(i,j) - -M
_{3}*(1−u(i,j)) ≤ y(i,j)−y(i,j−1)−1 ≤ +M_{3}*(1−u(i,j)) - x(i,j) ≤ M
_{4}*(1−u(i,j))

The constants M_{1} through M_{4} 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
_{1}=x.up+1 - M
_{2}=y.up - M
_{3}=y.up+1 - M
_{4}=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).