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)

$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

Keywords: mixed integer linear programming, relaxed mixed integer linear
          programming, general assignment problem, lagrangian relaxation, knapsack
$offText

$eolCom //

$sTitle Original Model Definition
Set
   i 'resources'
   j 'items';

Variable
   x(i,j) 'assignment of i to j'
   z      'total cost of assignment';

Binary Variable x;

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

Parameter
   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
Set
   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;

Parameter 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
Set
   id(i) 'dynamic version of i used to define a subset of i'
   iter  'subgradient iteration index' / iter1*iter20 /;

Alias (i,ii);

Parameter
   w(j)   'Lagrangian multipliers'
   improv 'has the Lagrangian bound improved over the previous iterations';

Variable zlrx 'relaxed objective';

Equation
   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 /;

Scalar
   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    /;

Parameter
   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;