Monday, March 23, 2009

Indexing by a variable

> I've:
>
> param N_MAX = 1000;
> var N >= 1, <= N_MAX;
>
> var x {0 .. N_MAX} >= 0;
>
> I want to do *something like*:
>
> subject to {i in  {(N/2)-1, (N/2)+1} }: #calculation of indices
> simplified for sake of discussion
>     x[i] = 1000; #RHS simplified for sake of discussion
>
> But of course I can't use variables for indexing. Neither can I do it
> explicitly as there are too many possible values that N can take. So,
> I'm looking for suggestions.
>
> Sincerely,

param N_MAX := 10;

set I := {1..N_MAX};
set J := {0..N_MAX};
set IJ{i in I} := setof{j in J: ((j >= i/2-1) and (j <= i/2+1))}(j);
display IJ;

var N >= 1, <= N_MAX;
param xup := 10000;
var x {J} >= 0, <= xup;
var b {I} binary;

minimize obj: sum{j in J} x[j];
s.t. e1: sum{i in I} i*b[i] = N;
s.t. e2: sum{i in I} b[i] = 1;
s.t. e3{i in I, j in IJ[i]}:  x[j] <= 1000*b[i] + xup*(1-b[i]);
s.t. e4{i in I, j in IJ[i]}:  x[j] >= 1000*b[i];

solve;
display N,b,x;
end;

Usually it is best to introduce a binary variable b[i] with b[i]=1 <=> i=N. Then you can introduce constraints based on b[i]. xup is here an upper bound on x[j]. Choose this number as tight as possible: it is a big-M constant. As usual it is a good idea to move as much as possible complexity away from the constraints. In this case we use an intermediate set IJ to handle a more difficult concept. Such a set can be tested and debugged before a complete model is ready to run.