lop.gms : Line Optimization

Description

The problem finds line plans for a given rail network and origin
destination demand data. Models for minimum cost and direct traveler
objectives are given. The set of possible lines is defined by the
shortest paths in the rail network.


Large Model of Type : MIP


Category : GAMS Model library


Main file : lop.gms

$title Line Optimization (LOP,SEQ=221)

$onText
The problem finds line plans for a given rail network and origin
destination demand data. Models for minimum cost and direct traveler
objectives are given. The set of possible lines is defined by the
shortest paths in the rail network.


Bussieck, M R, Optimal Lines in Public Rail Transport. PhD thesis,
TU Braunschweig, 1998.

Bussieck, M R, Kreuzer, P, and Zimmermann, U T, Optimal Lines for
Railway Systems. European Journal of Operation Research 96, 1 (1996),
54-63.

Claessens, M T, van Dijk, N M, and Zwaneveld, P J, Cost Optimal
Allocation of Rail Passenger Lines. European Journal Operation
Research 110, 3 (1998), 474-489.

Keywords: linear programming, mixed integer linear programming, passenger railway optimization,
          shortest path, dutch railway, public rail transport, network optimization
$offText

$eolCom //

Set s 'stations'
      / Ah   'Arnhem',             Apd  'Apeldoorn',        Asd 'Amsterdam CS'
        Asdz 'Amsterdam Zuid WTC', Asn  'Assen',            Bd  'Breda'
        Ehv  'Eindhoven',          Gn   'Groningen',        Gv  'Den Haag HS'
        Gvc  'Den Haag CS',        Hgl  'Hengelo',          Hr  'Heerenveen'
        Lls  'Lelystad Centrum',   Lw   'Leeuwarden',       Mt  'Maastricht'
        Odzg 'Oldenzaal Grens',    Rsdg 'Roosendaal Grens', Rtd 'Rotterdam CS'
        Shl  'Schiphol',           Std  'Sittard',          Ut  'Utrecht CS'
        Zl   'Zwolle',             Zvg  'Zevenaar Grens'                      /;

Table rt(s,s) 'running time (defines edges of rail network)'
        Hr Asn Zl Hgl Ut Shl Asdz Asd Gv Gvc Rtd Bd Ehv Std Mt Lls Rsdg Zvg Odzg
     Lw 29
     Gn     28
     Hr        66
    Asn        78
     Zl               85                                        50
    Apd            69 64           89
    Hgl                                                                       18
     Ah               58                                                 19
     Ut                        34  39     61  57 92  81
    Shl                         9  19 43  43
   Asdz                    9                                    56
    Asd                                                         54
     Gv                                    1  23
    Rtd                                          49                  18
    Ehv                                                  78
    Std                                                     21                 ;

Parameter tt(s) 'turn-a-round time'
                / (Ehv,Gn,Gvc,Hgl,Lls,Rsdg,Std)       5.0
                  (Apd,Asn,Gv,Hr,Mt,Shl,Zvg)         13.8
                  (Ah,Asd,Asdz,Bd,Lw,Odzg,Rtd,Ut,Zl) 14.1 /;

Table lfr(s,s) 'line frequency requirement'
        Hr Asn Zl Hgl Ut Shl Asdz Asd Gv Gvc Rtd Bd Ehv Std Mt Lls Rsdg Zvg Odzg
     Lw  1
     Gn      1
     Hr         1
    Asn         1
     Zl                1                                         1
    Apd             1  1            1
    Hgl                                                                        1
     Ah                2                                                  1
     Ut                         2   2      2   1  1   2
    Shl                         1   3  2   1
   Asdz                    1                                     1
    Asd                                                          1
     Gv                                    3   1
    Rtd                                           2                   1
    Ehv                                                   1
    Std                                                      1                 ;

Table od(s,s) 'origin destination matrix'
         Hr  Asn  Zl Hgl Ah   Ut Shl Asdz  Asd   Gv  Gvc  Rtd   Bd  Ehv Std  Mt  Lls Rsdg Zvg Odzg
     Lw 478      380     13  145  20   21   90    6   26   36   14    9   9   4   77    7  14
     Gn     1720 720         331  48   88  205   12   73   75   34   28  29  13  200   33  14
     Hr          511     11  209  20   16  115   10   48   58   16   11   8   4   77   10  19
    Asn          854     16  502  32   58  235   13  117  125   42   33  28  14  152   48  19
     Zl                  56 1112  64  171  400   33  163  182   79   47  46  21  390  100  32
    Apd              468    1160  32   76  917   21  202  143   57   62  10   5        47  83   71
    Hgl                      422  11   24  287   20   81   52   39   28  20  12        24       75
     Ah                     4244  60  721  726  109  741  180  136  101            8  320 602
     Ut                          278 5826 4919  225 3138 2260 1165 3109 720 359   89  325 996   21
    Shl                              1456 6469 1339 1503  509    7   99  44  29  103  164
   Asdz                                         461  207  369  138  542 203 149  819    6 155
    Asd                                         730 2540 1756  154  437 155  37 2783 2258 489   22
     Gv                                              785 4586  531   35  22   8   29  890
    Gvc                                                  2829  228  335 104  41   31    3 229    7
    Rtd                                                       1829  569 179  73   46 1077 157   11
     Bd                                                             950 157  79    6  329  14    5
    Ehv                                                                 936 404    8   75  11    3
    Std                                                                     863    2   19
     Mt                                                                            1   22
    Lls                                                                                15         ;

Set
   lf 'line frequencies'          / 1 'every 60 minutes', 2 'every 30 minutes' /
   ac 'additional cars per train' / 0*9 /;

Scalar
   mincars 'minimum number of cars per train' /      3 /
   ccap    'passenger capacity per car'       /    467 /
   cfx     'fix cost per car'                 / 353100 /
   crm     'cost per ride minute per car'     /   5803 /
   trm     'cost per ride minute per train'   /  44959 /
   cmp     'additional cost multiplier'       /     90 /
   maxtcap 'maximum train capacity';

maxtcap = (mincars + card(ac) - 1)*ccap;

Alias (s,s1,s2,s3);

* Some data checks
Set error01(s,s);
error01(s1,s2) = rt(s1,s2) and not lfr(s1,s2) or not rt(s1,s2) and lfr(s1,s2);
abort$card(error01) 'Inconsistent Edge data', error01;

$onText
Generate Shortest Path from all nodes to all nodes. Instead of solving s
network problems we solve all shortest path problems simultaneously.
$offText

Set d(s,s) 'directed version of the network';
d(s1,s2) = rt(s1,s2) or rt(s2,s1);

Variable
   f(s,s1,s2) 'flow from s on edge s1s2'
   spobj      'objective variable';

Positive Variable f;

Equation
   balance(s,s1) 'flow balance constraint for flow from s in s1'
   defspobj      'definition of the combined min cost flow objective';

balance(s,s1)..
   sum(d(s1,s2), f(s,d)) =e= sum(d(s2,s1), f(s,d)) + sameas(s,s1)*card(s) - 1;

defspobj..
   spobj =e= sum((s,d(s1,s2)), f(s,d)*max(rt(s1,s2),rt(s2,s1)));

Model sp 'shortest path model' / balance, defspobj /;

solve sp minimzing spobj using lp;

Set tree(s,s1,s2) 'shortest path tree from s';
tree(s,s1,s2) = f.l(s,s1,s2);
abort$(card(tree) <> card(s)*(card(s) - 1)) 'wrong tree computation';

* From the trees we generate all shortest path lines. We perfom a BFS.
Set
   r             'rank' / 1*100 /
   k(s,s)        'arcs from root to a node'
   v(s,r)        'nodes with rank from root to a node'
   unvisit(s)    'unvisited nodes'
   visit(s)      'visited nodes'
   from(s)       'from nodes'
   to(s)         'to nodes'
   l(s,s1,s2,s3) 'line from s to s1 with edge s2s3'
   lr(s,s1,s2,r) 'rank of s2 in line from s to s1';

Alias (s,root), (r,r1);

l(root,s,s1,s2) = no;
lr(root,s,s1,r) = no;

loop(root,
   from(root) = yes;  // We start with the root node
   unvisit(s) = yes;  // All nodes are unvisited
   visit(s)   =  no;  // No node is visited
   loop(r$(ord(r) > 1 and card(unvisit)),
      unvisit(from) =  no;
      visit(from)   = yes;
      // nodes that can be reach from a node in from
      to(unvisit) = sum(tree(root,from,unvisit), yes);
      loop(from,
         k(s2,s3)$l(root,from,s2,s3)  = yes;  // arcs of line root-from
         v(s2,r1)$lr(root,from,s2,r1) = yes;  // nodes of line root-from
         v(from,"1")$(card(k) = 0)    = yes;  // this is the root node

         // the line root-to has the arcs/nodes of line root-from
         l(root,to,k)$tree(root,from,to)       = yes;
         lr(root,to,v)$tree(root,from,to)      = yes;
         // plus the arc from-to
         l(root,to,from,to)$tree(root,from,to) = yes;
         lr(root,to,to,r)$tree(root,from,to)   = yes;

         k(s2,s3) = no;
         v(s2,r1) = no;
      );
      from(s)  =  no;
      from(to) = yes;  // move one layer down
      to(s)    =  no;
   );
   from(s) = no;
);

Set error02(s1,s2) 'arcs not covered by shortest path lines';
error02(s1,s2) = lfr(s1,s2) and sum(l(root,s,s1,s2), 1) = 0;
abort$card(error02) error02;

* Lines are symetric, so delete one half of them
Set ll(s,s) 'station pair represening a line';
ll(s1,s2) = ord(s1) < ord(s2);

l(root,s,s1,s2)$(not ll(root,s)) = no;
lr(root,s,s1,r)$(not ll(root,s)) = no;

* and order the edges in the lines in the way we stored them
l(root,s,s1,s2)$(l(root,s,s2,s1) and rt(s1,s2)) = yes;
l(root,s,s1,s2)$(not rt(s1,s2)) = no;

Parameter
   rp(s,s,s)   'rank of node'
   lastrp(s,s) 'rank of the last node in line';

rp(ll,s)   = sum(r$lr(ll,s,r), ord(r));
lastrp(ll) = smax(s,rp(ll,s));

Parameter load(s1,s2) 'passenger load of an edge';
load(s1,s2)$rt(s1,s2) = sum(l(root,s,s1,s2)$od(root,s), od(root,s));

$onText
Model dtlop:
   Determines a line plan with a maximizing the number of direct
   travelers. The number of direct traveler represents an upper
   bound because the capcity constraint is relaxed.
$offText

Variable
   dt(s1,s2)   'direct traveler between s1 and s2'
   freq(s1,s2) 'frequency on arc s1s2'
   phi(s1,s2)  'frequency of line between s1 and s2'
   obj         'objective variable';

Integer Variable phi;

Equation
   deffreqlop(s1,s2) 'definition of the frequency for each edge'
   dtlimit(s1,s2)    'limit the direct travelers'
   defobjdtlop       'objective function';

deffreqlop(s1,s2)$rt(s1,s2)..
   freq(s1,s2) =e= sum(l(ll,s1,s2), phi(ll));

dtlimit(s1,s2)$od(s1,s2)..
   dt(s1,s2) =l= min(od(s1,s2),maxtcap)*sum(ll$(rp(ll,s1) and rp(ll,s2)), phi(ll));

defobjdtlop..
   obj =e= sum((s1,s2)$od(s1,s2), dt(s1,s2));

Model lopdt / deffreqlop, dtlimit, defobjdtlop /;

freq.lo(s1,s2)$rt(s1,s2) = max(lfr(s1,s2),ceil(load(s1,s2)/maxtcap));
freq.up(s1,s2)$rt(s1,s2) = freq.lo(s1,s2);
dt.up(s1,s2)$od(s1,s2)   = od(s1,s2);

solve lopdt maximizing obj using mip;

* Store the solution for further reporting
Parameter solrep, solsum;
solrep('DT',ll,'freq') = phi.l(ll);
solrep('DT',ll,'cars')$phi.l(ll) = mincars + card(ac) - 1;

$onText
Model ILP:
   Determines a line plan of minimum cost.
$offText

Parameter
   xcost(root,s,lf) 'operating and capcital cost for line with mincars cars'
   ycost(root,s,lf) 'operating and capcital cost for additional cars'
   len(s,s)         'length of line'
   sigma(s,s)       'line circulation factor';

len(ll)      =  sum(l(ll,s1,s2), rt(s1,s2));
sigma(ll)    = (len(ll) + sum(s$lr(ll,s,"1"), tt(s)) + sum(s$(rp(ll,s) = lastrp(ll)), tt(s)))/60;
xcost(ll,lf) =  ord(lf)*len(ll)*(trm + mincars*crm) + mincars*ceil(sigma(ll)*ord(lf))*cfx;
ycost(ll,lf) =  ord(lf)*len(ll)*crm + ceil(sigma(ll)*ord(lf))*cfx;

Variable
   x(s1,s2,lf) 'line frequency indicator of line s1-s2'
   y(s1,s2,lf) 'additional cars on line s1-s2 with frequency lf';

Integer Variable y;
Binary  Variable x;

Equation
   deffreqilp(s,s)  'definition of the frequency for each edge'
   defloadilp(s,s)  'capacity of lines fulfill the demand'
   oneilp(s,s)      'only one frequency per line'
   couplexy(s,s,lf) 'coupling constraints'
   defobjilp        'definition of the objective';

deffreqilp(s1,s2)$rt(s1,s2)..
   freq(s1,s2) =e= sum((l(ll,s1,s2),lf), ord(lf)*x(ll,lf));

defloadilp(s1,s2)$rt(s1,s2).. ceil(load(s1,s2)/ccap) =l=
   sum((l(ll,s1,s2),lf), ord(lf)*(mincars*x(ll,lf) + y(ll,lf)));

oneilp(ll)..
   sum(lf, x(ll,lf)) =l= 1;

couplexy(ll,lf)..
   y(ll,lf) =l= y.up(ll,lf)*x(ll,lf);

defobjilp..
   obj =e= sum((ll,lf), xcost(ll,lf)*x(ll,lf) + ycost(ll,lf)*y(ll,lf));

Model ilp / defobjilp, deffreqilp, defloadilp, oneilp, couplexy /;

y.up(ll,lf) = card(ac) - 1;
freq.up(s1,s2)$rt(s1,s2) = 100;

ilp.optCr  = 0;
ilp.resLim = 100;

solve ilp minimizing obj using mip;

solrep('ILP',ll,'freq') = sum(lf$x.l(ll,lf), ord(lf));
solrep('ILP',ll,'cars') = sum(lf$x.l(ll,lf), mincars + y.l(ll,lf));
solsum('ILP','cost')    = obj.l;

* We have now two line plans. Lets make a comparison: Cost and Direct
* Travelers. The Direct Traveler Evaluation requires another more
* detailed model. Model lopdt gave just an upper bound.

$onText
Model EvalDT:
   Determines for a given line plan (routes, frequency, capacity) the
   maximum number of direct travelers.
$offText

Parameter cap(s,s) 'the capacity of a line';
Set       sol(s,s) 'the actual lines in a line plan';

Positive Variable dtr(s,s,s,s) 'direct travelers of OD pair u v in line on route s s';

Equation
   dtllimit(s,s1,s2,s3) 'limit direct travelers in line s-s1 on edge s2-s3'
   sumbound(s,s)        'sum of direct travels <= total number of travelers';

dtllimit(l(sol,s,s1))..
   sum((s2,s3)$(od(s2,s3) and rp(sol,s2) and rp(sol,s3) and
       (min(rp(sol,s),rp(sol,s1)) >= rp(sol,s2) and // s and s1 must be
        max(rp(sol,s),rp(sol,s1)) <= rp(sol,s3)  or // between the nodes of the
        min(rp(sol,s),rp(sol,s1)) >= rp(sol,s3) and // origin destination pair
        max(rp(sol,s),rp(sol,s1)) <= rp(sol,s2))),  // s2-s3 in order to
   dtr(sol,s2,s3)) =l= cap(sol);                    // occupy capacity of s s1

sumbound(s2,s3)$od(s2,s3)..
   sum(sol$(rp(sol,s2) and rp(sol,s3)), dtr(sol,s2,s3)) =e= dt(s2,s3);

Model evaldt / dtllimit, sumbound, defobjdtlop /;

* Evaluate direct travelers for DT line plan
sol(ll)  = solrep('DT',ll,'freq');
cap(sol) = solrep('DT',sol,'freq')*solrep('DT',sol,'cars')*ccap;

solve evaldt maximizing obj using lp;

solsum('DT','dtrav') = obj.l;
solsum('DT','cost')  = sum(sol,  solrep('DT',sol,'freq')*len(sol)*trm
                              + (solrep('DT',sol,'freq')*len(sol)*crm
                              +  ceil(sigma(sol)*solrep('DT',sol,'freq'))*cfx)
                              *  solrep('DT',sol,'cars'));

* Evaluate DT for ILP line plan
sol(ll)  = solrep('ILP',ll,'freq');
cap(sol) = solrep('DT',sol,'freq')*solrep('DT',sol,'cars')*ccap;

solve evaldt maximizing obj using lp;
solsum('ILP','dtrav') = obj.l;

display solrep, solsum;