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.


References

  • 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.

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.

$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

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

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

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

Equations
    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

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