$title Data Envelopment Analysis - DEA (extended,GUSS)
$ontext
Data Envelopment Analysis (DEA) is a technique for measuring the relative
performance of organizational units where presence of multiple inputs and
outputs makes comparison difficult.

            efficiency = weighted sum of output / weighted sum of input

Find weights that maximize the efficiency for one unit while ensuring
that no other units has an efficiency < 1 using these weights. A primal
and dual formulation is presented.


Dyson, Thanassoulis, and Boussofiane, A DEA Tutorial.
Warwick Business School

$offtext

sets  i  units / Depot1*Depot20 /
      j     inputs and outputs / stock, wages, issues, receipts, reqs /
      ji(j) inputs             / stock, wages                         /
      jo(j)            outputs /               issues, receipts, reqs /;
alias(i,k);

Table data(i,j)
         stock  wages  issues  receipts  reqs
Depot1     3      5      40       55      30
Depot2     2.5    4.5    45       50      40
Depot3     4      6      55       45      30
Depot4     6      7      48       20      60
Depot5     2.3    3.5    28       50      25
Depot6     4      6.5    48       20      65
Depot7     7     10      80       65      57
Depot8     4.4    6.4    25       48      30
Depot9     3      5      45       64      42
Depot10    5      7      70       65      48
Depot11    5      7      45       65      40
Depot12    2      4      45       40      44
Depot13    5      7      65       25      35
Depot14    4      4      38       18      64
Depot15    2      3      20       50      15
Depot16    3      6      38       20      60
Depot17    7     11      68       64      54
Depot18    4      6      25       38      20
Depot19    3      4      45       67      32
Depot20    3      6      57       60      40
;

Parameter vlo       v lower bound
          ulo       u lower bound
          norm      normalizing constant
          slice(j);

Variables v(ji) input weights
          u(jo) output weights
          eff   efficiency
          var   dual convexity

          lam(i) dual weights
          vs(ji) input duals
          us(jo) output duals
          z;

positive variables u,v,vs,us,lam;

Equations defe  efficiency definition - weighted output
          denom weighted input
          lime(i)  'output / input < 1'
          dii(ji) input duals
          dio(jo) output dual
          defvar    variable return to scale
          dobj      dual objective;

*  primal model

defe..   eff =e= sum(jo, u(jo)*slice(jo)) - 1*var;

denom..  sum(ji, v(ji)*slice(ji)) =e= norm;

lime(i)..    sum(jo, u(jo)*data(i,jo)) =l= sum(ji, v(ji)*data(i,ji)) + var;

*  dual model

dii(ji).. sum(i, lam(i)*data(i,ji)) + vs(ji) =e= z*slice(ji);

dio(jo).. sum(i, lam(i)*data(i,jo)) - us(jo) =e=     slice(jo);

defvar..     sum(i, lam(i)) =e= 1;

dobj.. eff =e= norm*z - vlo*sum(ji, vs(ji)) - ulo*sum(jo, us(jo));


model deap  primal / defe, denom,lime /
      deadv dual with VRS / dobj, dii, dio, defvar  /;

* to run VRS with the primal model
var.lo = -inf; var.up = +inf;
norm=100;

vlo=1e-4; ulo=1e-4;
v.lo(ji) = vlo; u.lo(jo) = ulo;

set headers report / modelstat, solvestat, objval /;
parameter scenrep(k,headers) solution report summary
          scopt / SkipBaseCase 1 /
          eff_k(k) 'efficiency report parameter';

set dict / k.     scenario.''
           scopt. opt.     scenrep
           slice. param.   data
           eff.   level.   eff_k /;

slice(j) = 0; option lp=gurobi;
solve deap using lp max eff scenario dict;

display "Solution primal problem:",scenrep;

slice(j) = 0; option lp=gurobi;
solve deadv using lp min eff scenario dict;

display "Solution dual problem:",scenrep,eff_k;

$if not set runtraditional $exit

parameter rep summary report;

loop(k,
   slice(j) = data(k,j);
   solve deap using lp max eff;
   rep(i,k) =  sum(jo, u.l(jo)*data(i,jo))/sum(ji, v.l(ji)*data(i,ji));
   rep('MStat-p',k) = deap.modelstat;
   solve deadv using lp min eff ;
   rep('MStat-d',k) = deadv.modelstat;
   rep('obj-check',k) = deadv.objval - deap.objval;
);

rep(i,'Min') = smin(k, rep(i,k));
rep(i,'Max') = smax(k, rep(i,k));
rep(i,'Avg') =  sum(k, rep(i,k))/card(k);

display rep;