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