! Job Load Leveling  (LoadLevel).
  Given a
    Set of jobs, each with an
     earliest begin time,
     latest finish time,
     processing time, and
     load,  
   Decide 
     when each job should be started
     so that the total load in process at
     any instant is minimized.

   The loads could be amount of electricity used by
 various appliances, e.g., dishwasher, hot water heater, iron.
 We want to find when to run each appliance so maximum electricity
 used at any instant is minimized;

! Keywords: Load leveling, Sequencing, Scheduling, Due dates, Ready time,
            Disjunctive formulation;
SETS:
  job: r, p, d, v, s, load;
  jxj(job,job) | &1 #LT# &2: y1, y2, y3, y4, 
                       ds1, ds2, ds3, ds4;
ENDSETS
DATA:
 !      1  2  3  4  5  6  7  8  9 10 11
 ! The processing times;
   p =  1  3 18  5 24  4  8 13 12  2  6;
 ! Ready times;
   r =  5 12  0 20  0  2  2  3  5  6 20;
 ! Due date of each job;
   d = 10 16 22 32 25  8 11 29 25 14 29;
 ! Load associated with each job;
   v = 5   4  2  1  6  1  2  3  7  4  3;
ENDDATA
! Variables:
    s(j) = start time for job j,
    load(j) = load in process at instant just after 
              job j starts,
  for i < j: ( 4 overlap configurations for i and j)
    y1(i,j) = 1 if i finishes before j starts,
      ds1(i,j) = s(j) - s(i), 0 otherwise,
    y2(i,j) = 1 if i starts before j starts but
                        finishes after j starts,
      ds2(i,j) = s(j) - s(i), 0 otherwise,
    y3(i,j) = 1 if i starts after j starts but
                         not before j finishes,
      ds3(i,j) = s(i) - s(j), 0 otherwise,
    y4(i,j) = 1 is i starts after j finishes,
      ds4(i,j) = s(i) - s(j), 0 otherwise;

! Minimize the maximum load at any instant;
 MIN = mxload;

! Cannot start job i too early or too late;
 @FOR( job(i):
     s(i) >= r(i);
     s(i) + p(i) <= d(i);
     );

! The sequencing constraints. Define the disjunctions/scenarios
   of the 4 possible overlap types for i and j;
   @FOR( jxj( i,j) :
  ! Case 1, i finishes before j starts;
         ds1(i,j) >= p(i)*y1(i,j);
         ds1(i,j) <= (d(j)-p(j)-r(i))*y1(i,j);
  ! Case 2, i starts before j starts, but finishes after j starts;
         ds2(i,j) <= p(i)*y2(i,j);
  ! Case 3, i starts after j starts but before j finishes;
         ds3(i,j) <= p(j)*y3(i,j);
  ! Case 4, i starts after j finishes;
         ds4(i,j) >= p(j)*y4(i,j);
         ds4(i,j) <= (d(i)-p(i)-r(j))*y4(i,j);

  ! Tie the cases together;
  ! Jobs i and j must be performed in one of the 4 configurations;
        y1(i,j)  + y2(i,j)  + y3(i,j)  + y4(i,j)  = 1;
  ! Difference in start times of i and j;
        ds1(i,j) + ds2(i,j) - ds3(i,j) - ds4(i,j) = s(j) - s(i) ;
  ! Must be binary, 0 or 1;
    @BIN(y1(i,j)); @BIN(y2(i,j)); @BIN(y3(i,j)); @BIN(y4(i,j));
            );

  ! Compute load at various instants.
    We need look only at job start times.
    The load just after start of job k is
     v(k) + load of jobs i or j starting before k but
          not yet finished when k starts;
   @FOR( job(k):      
      load(k) = v(k) +  
       @SUM( jxj(i,k): y2(i,k)*v(i)) +  ! i < k;
       @SUM( jxj(k,j): y3(k,j)*v(j));   ! j > k;
      mxload >= load;   ! Force the maximum;
       );