! Using Lagrangean relaxation to try to get
better bounds, and perhaps solve a 
Task to Machine assignment problem;
!Reference:
   Geoffrion, A. (1974),"Lagrangean Relaxation for Integer Programming",
Math. Programming Study 2, pp.82-114;
! Keywords: Lagrangean relaxation, Task Assignment;
SETS:
 TASK: LAMBDA;
 MACH: TIMEA;
 TXM( TASK, MACH): VAL, USAGE, Y, YCUR;
ENDSETS
DATA:
! A set of tasks...;
 TASK = T01  T02  T03  T04  T05  T06  T07;
! and machines to which tasks can be assigned;
 MACH =   M01  M02  M03;
! Time capacity available on each machine;
 TIMEA =  40  40  60;
! Value of assigning task i to machine j;
 VAL =    24  43  41
          17  52  28
          81  73  74
          39  48  47
          55  62  58
          73  71  68
          29  38  36;
! Usage by task i of machine j capacity;
 USAGE =  24  28  19
          32  31  29
          21  18  15
          17  14  21
          11  13  19
          17  12  17
          15  15  17;
ENDDATA

SUBMODEL REAL_PROB:
 ! The original, full model is...;
  MAX = OBJ;
    OBJ = @SUM( TASK(i): 
          @SUM( MACH(j): VAL(i,j)*Y(i,j)));

! Cannot assign more work to machine j than capacity available;
   @FOR( MACH(j):
   [CAP]  @SUM( TASK(i): USAGE(i, j)*Y(i,j)) <= TIMEA(j)
       );

! Can assign task i to at most one machine.
  These constraints will be relaxed vial Lagrange multipliers;
   @FOR( TASK(i):
    [ATMOST1] @SUM( MACH(j): Y(i,j)) <= 1; ! Apply muliplier LAMBD(i);
       );

!  The Y(i,j) must be 0 or 1. Can turn this off by
  setting NORELAX = 0;
   @FOR( TXM( i,j) | NORELAX : @BIN( Y(i,j)));
! End of original full model;
ENDSUBMODEL


SUBMODEL SOLVE_MACH:
! Relaxed problem for machine jnow if the ATMOST1 constraints
 are relaxed.  The problem decomposes into a separate problem
 for each machine;
 ! Inputs:
    jnow = machine for which we are solving knapsack/loading problem;
 ! Maximize value minus lagrange multiplier charge of tasks assigned to machine jnow;
  MAX = OBJ;
   OBJ = @SUM( TASK(i): VAL(i,jnow)*Y(i,jnow)) 
       - @SUM( TASK(i): LAMBDA(i)*Y(i,jnow));
 ! Cannot assign more work to machine jnow than time available;
  [CAPJ] @SUM( TASK(i): USAGE(i, jnow)*Y(i,jnow)) <= TIMEA(jnow);
 ! The Y(i,j) must be 0 or 1, i.e., binary;
  @FOR( TASK(i): @BIN( Y(i,jnow)));
ENDSUBMODEL


CALC:
  @SET('TERSEO',2); ! Output level (0:verb, 1:terse, 2:errors, 3:none);
  @SET('OROUTE',1); ! Route every line immediately to the window;

! If we want to cheat and solve the orginal problem, turn this on;
!  NORELAX = 0; !  @SOLVE( REAL_PROB);
!  @WRITE(' True optimal value= ', OBJ, @NEWLINE(1));

 ! Guess initial multipliers;
  @FOR( TASK(i):
    LAMBDA(i) = 0;
     );
 ! Initial multiplier adjustment factor;
  DELTA = 10;
 ! Max iteration limit;
  ITERMX = 100;
 ! Current interation number;
  iter = 0;

  @WRITE('  Iteration    Upper bound   Lower bound   Step size', @NEWLINE(1));

 ! Main loop, adjusting Lagrange multipliers each pass.
    One may wish to add additional termination conditions
    besides a max iterations limit;
  @WHILE( iter #LT# ITERMX:
   iter = iter + 1;

 ! For given set of Lagrange Multipliers, solve current relaxed problem.
   The relaxation separates into independent problems for each machine;
  UBOUND = @SUM( TASK(i): LAMBDA(i)*1);
  @FOR( MACH( j):
    jnow = j;
    @SOLVE( SOLVE_MACH);
    UBOUND = UBOUND + OBJ;
! Save the current Y values;
    @FOR( TASK(i):
      YCUR(i,jnow) = Y(i,jnow);
       );
     );

! Begin code to adjust the multipliers;
 @FOR( TASK(i):
! Compute the LHS of the relaxed constraint;
   LHS = @SUM( MACH(j) | YCUR(i,j) #GT# 0.5: 1);
 ! Begin code to adjust multiplier based on the amount of constraint violation.
   Increase LAMBDA(i) if violated, decrease if there is slack.
    There may be better ways to do this;
   LAMBDA(i) = @SMAX(0, LAMBDA(i) + DELTA*(LHS - 1));
!   @WRITE(i,' LHS=', LHS,' Adjust multiplier to ', LAMBDA(i),@NEWLINE(1));  
      );
    
 ! Multiplier adjustment step size adjustment.
     The stepsize must decrease towards 0 to get convergence.
     There may be better ways to this;
  DELTA = DELTA*(iter+1)/(iter+2);
 ! End code to adjust multipliers and stepsize;

 ! Begin code to try to find a feasible solution from
    the relaxed solution.
    There may be better ways of guessing a feasible solution;
 @FOR( TASK(i):
! Compute the LHS of the relaxed constraint;
   LHS = @SUM( MACH(j) | YCUR(i,j) #GT# 0.5: 1);
 ! If LHS > 1, then decrease LHS by setting some of the YCUR(i,j)
   = 0 until constraint is feasible;
     @FOR( MACH(j) | YCUR(i,j) #GE# 0.5 #AND# LHS #GT# 1:
         YCUR(i,j) = 0;
         LHS = LHS - 1;
       );
     );

 !  Compute the objective value of this feasible solution;
   LBOUND = @SUM( TXM(i,j): VAL(i,j)*YCUR(i,j));
 ! End code to try to find a feasible solution;

 ! If the original problem is nonconvex, e.g., has integer variables,
  then it is generally impossible to reduce the gap between
  UBOUND and LBOUND to 0;
 @WRITE(@FORMAT(iter,'5.0f'),'      ', @FORMAT( UBOUND,'12.2f'),
        '    ', @FORMAT(LBOUND,'12.2f'),'  ',@FORMAT(DELTA,'10.3f'),@NEWLINE(1));

   );  ! @WHILE( iter...;
    
ENDCALC