gapmin.gms : Lagrangian Relaxation of Assignment Problem

Description

A general assignment problem is solved via Lagrangian Relaxation
by dualizing the multiple choice constraints and solving
the remaining knapsack subproblems.

The data for this problem are taken from Martello.
The optimal value is 223 and the optimal solution is:
    1 1 4 2 3 5 1 4 3 5, where
in columns 1 and 2, the variable in the first row is equal to 1,
in column 3, the variable in the fourth row is equal to 1, etc...


References

  • Martello, S, and Toth, P, Knapsack Problems: Algorithms and Computer Implementations. John Wiley and Sons, Chichester, 1990.
  • Guignard, M, and Rosenwein, M, An Improved Dual-Based Algorthm for the Generalized Assignment Problem. Operations Research 37 (1989), 658-663.

Small Model of Types : MIP rmip


Category : GAMS Model library


Main file : gapmin.gms

$title Lagrangian Relaxation for Generalized Assignment (GAPMIN,SEQ=182)
$stitle original model definition
$eolcom !
$Ontext
  A general assignment problem is solved via Lagrangian Relaxation
  by dualizing the multiple choice constraints and solving
  the remaining knapsack subproblems.

  The data for this problem are taken from Martello.
  The optimal value is 223 and the optimal solution is:
      1 1 4 2 3 5 1 4 3 5, where
  in columns 1 and 2, the variable in the first row is equal to 1,
  in column 3, the variable in the fourth row is equal to 1, etc...


  Martello, S, and Toth, P, Knapsack Problems: Algorithms and Computer
  Implementations. John Wiley and Sons, Chichester, 1990.

  Guignard, M, and Rosenwein, M, An Improved Dual-Based Algorithm for
  the Generalized Assignment Problem. Operations Research 37 (1989), 658-663.


 --- original model definition

$Offtext

sets  i   resources
      j   items

binary variable x(i,j) assignment of i to j
       variable z      total cost of assignment

equations capacity(i)  resource availability
          choice(j)    assignment constraint.. one resource per item
          defz         definition of total cost;

parameters  a(i,j)  utilization of resource i by item j
            f(i,j)  cost of assigning item j to resource i
            b(i)    available resources;

capacity(i)..  sum(j, a(i,j)*x(i,j)) =l= b(i);

choice(j)..    sum(i, x(i,j)) =e= 1;

defz..         z =e= sum((i,j), f(i,j)*x(i,j));

model assign original assignment model / capacity, choice, defz /;

* --- data for Martello model
sets  i   resources / r1 *r5   /
      j   items     / i1 * i10 /

      xopt(i,j) optimal assignment
      / r1.(i1,i2,i7),r2.i4, r3.(i5,i9), r4.(i3,i8), r5.(i6,i10) /;

table  a(i,j)  utilization of resource i by item j

       i1  i2  i3  i4  i5  i6  i7  i8  i9 i10
  r1   12   8  25  17  19  22   6  22  20  25
  r2    5  15  15  14   7  11  14  16  17  15
  r3   21  24  13  24  12  16  23  20  15   5
  r4   23  17  10   6  24  20  15  10  19   9
  r5   17  20  15  16   5  13   7  16   8   5

table  f(i,j) cost of assigning item j to resource i

       i1  i2  i3  i4  i5  i6  i7  i8  i9 i10
  r1   16  26  30  47  18  19  33  37  42  31
  r2   38  42  15  21  26  11  11  50  24  19
  r3   48  17  14  22  14  18  47  32  17  42
  r4   22  32  28  39  37  23  25  12  44  17
  r5   31  42  31  40  16  15  29  31  44  41

parameters b(i) available resources  / r1 28, r2 20, r3 27, r4 24, r5 19/


$ontext  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

          !!  if one wants to check the data, one can
          !!  solve the MIP problem, this is just a check

assign.optcr=0; solve assign minimizing z using mip;

if( sum(xopt, x.l(xopt) <> 1),
  abort '*** Something wrong with this solution', x.l, xopt );

$offtext   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

$stitle Relaxed Problem Definition and Subgradient Optimization
* Lagrangian subproblem definition
* uses dynamic set to define WHICH knapsack to solve

sets id(i)  dynamic version of i   used to define a subset of i
     iter   subgradient iteration index / iter1 * iter20  /
alias (i,ii)

parameters  w(j)   Lagrangian multipliers
            improv has the Lagrangian bound improved over the previous iterations

variable  zlrx    relaxed objective

equations knapsack(i) capacity with dynamic sets
          defzlrx     definition of zlrx;

knapsack(id)..  sum(j, a(id,j)*x(id,j)) =l= b(id);

defzlrx..       zlrx =e= sum((id,j), (f(id,j)-w(j))*x(id,j));

model pknap / knapsack, defzlrx /;



scalars target target objective function value
        alpha  step adjuster / 1 /
        norm   norm of slacks
        step   step size for subgradient / na /
        zfeas  value for best known solution or valid upper bound
        zlr    Lagrangian objective value
        zl     Lagrangian objective value
        zlbest current best Lagrangian lower bound
        count  count of iterations without improvement
        reset  reset count counter / 5 /
        tol    termination tolerance / 1e-5 /
        status outer loop status /0/

parameters
   s(j)             slack variable
   report(iter,*)   iteration log
   xrep(j,i,*)      x iteration report
   srep(iter,j)     slack report
   wrep(iter,j)     w iteration report ;

* --- calculate initial Lagrangian multipliers
*     There are many possible ways to find initial multipliers.
*     The choice of initial multipliers is very important for the
*     overall performance. The marginals of the relaxed problem
*     are often used to initialize the multipliers. Another choice
*     is simply to start with zero multipliers.


* replace 'default' with solver of your choice.
option   mip = default
        rmip = default;

file results writes iteration report / solution /;

put results 'solvers used: RMIP = ' system.rmip /
            '               MIP = ' system.mip  /;

* --- solve relaxed problem to get initial multipliers
*     Note that different solvers get different dual solutions
*     which are not as good as a zero set of initial multipliers.

solve assign minimizing z using rmip;
put / 'RMIP objective value = ', z.l:12:6/;

if(assign.modelstat = %modelstat.Optimal%,
   status = %modelstat.Optimal%                        ! everything ok
else
   abort '*** relaxed MIP not optimal',
         '    no subgradient iterations', x.l );
xrep(j,i,'initial') = x.l(i,j);
xrep(j,i,'optimal') = 1$xopt(i,j);

parameter wopt(j) an optimal set of multipliers
          / i1 35,i2 40, i3 60, i4 69, i5 21, i6 49, i7 42, i8 47, i9 64, i10 46 /;

zlbest = z.l;

* --- use RMIP duals
w(j)   = choice.m(j);

* --- use optimal duals
*w(j) = wopt(j);

* --- use zero starting point
*w(j) = 0;
*zlbest=0;

put // 'zlbest                    objective value  = ', zlbest:12:6;
put // "Dual values on assignment constraint"/ ;
loop(j, put / "w('",j.tl,"') =  ", w(j):16:6 ";" );

*  one needs a value for zfeas
*  one can compute a valid upper bound as follows:

zfeas  = sum(j, smax(i, f(i,j)));
put // 'zfeas quick and dirty bound obj value      = ', zfeas:12:6;
display 'a priori upper bound',zfeas;

$ontext
*  another alternative to compute a value for zfeas is
*  to solve gapmin by B-B and stop
*  at first 0-1 feasible solution found
*  using gapmin.optcr = 1, as follows

assign.optcr=1;assign.solprint=%solprint.Quiet%;

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

solve assign minimizing z using mip;
zfeas=min(zfeas,z.l);
display 'final zfeas',zfeas;
display 'heuristic solution by B-B ',x.l,z.l;
put  / 'zfeas IP solution bound objective value    = ', zfeas.l:12:6;
$offtext

put /// 'Iteration         New Bound   Previous Bound            norm      abs(zl-zf)'/;


*  then keep the smaller of the two values as zfeas

pknap.optcr = 0;                    ! ask for global solution
pknap.solprint = %solprint.Quiet%;  ! turn off all solution output

*============================================================================*
*                                                                            *
*  beginning of subgradient loop                                             *
*                                                                            *
*============================================================================*

id(i)  = no;  ! initially empty
count  = 1;
alpha  = 1;

display status;

loop(iter$(status = 1),    ! i.e., repeat while status is 1

*  solve Lagrangian subproblems by solving nonoverlapping knapsack
*  problems. Note the use of the dynamic set id(i) which will
*  contain the current knapsack descriptor.

   zlr = 0;
   loop(ii,
      id(ii) = yes;                        ! assume id was empty
      solve pknap using mip minimizing zlrx;
      zlr = zlr + zlrx.l;
      id(ii) = no );                       ! make set empty again
   improv = 0;
   zl     = zlr + sum(j, w(j));
   improv$(zl > zlbest) = 1;               ! is zl better than zlbest?
   zlbest = max(zlbest,zl);
   s(j)   = 1 - sum(i, x.l(i,j));          ! subgradient
   norm   = sum(j, sqr(s(j)));

   status$(norm < tol)                             = 2;
   status$(abs(zlbest-zfeas) < 1e-4)               = 3;
   status$(pknap.modelstat <> %modelstat.Optimal%) = 4;
   put results / iter.tl ,zl:16:6,zlbest:16:6,norm:16:6,abs(zlbest-zfeas):16:6;
   if((status = 2),
      put //"subgr. method has converged, status = ",status:5:0//;
      put //"last solution found is optimal for IP problem"//;
     );    ! end if
   if((status = 3),
      put //"subgr. method has converged, status = ",status:5:0//;
      put //"no duality gap, best sol. found is optimal "//;
    );     ! end if
   if ((status = 4),
      put //"something wrong with last Lag. subproblem"//;
      put //"status = ",status:5:0//;
    );     ! end if

   report(iter,'zlr')    = zlr;
   report(iter,'zl')     = zl;
   report(iter,'zlbest') = zlbest;
   report(iter,'norm')   = norm;
   report(iter,'step')   = step;

   wrep(iter,j) = w(j);
   srep(iter,j) = s(j);
   xrep(j,i,iter) = x.l(i,j);

   if(status=1,
      target = (zlbest+zfeas)/2;
      step   = (alpha*(target-zl)/norm)$(norm > tol);
      w(j)   = w(j)+step*s(j);
      if(count>reset,          ! too many iterations w/o improvement
         alpha = alpha/2;
         count = 1
      else if(improv,          ! reset count if improvement
              count = 1
           else
              count = count +1 ! update count if no improvement
             )
        )
     )
);                                         ! end loop iter

display report, wrep, srep, xrep;
put results // "Dual values on assignment constraint" /;
loop(j, put /  "w('",j.tl,"') =  ", w(j):16:6  ";" )
put //"best Lagrangian bound   =   ",zlbest:10:5;