# Solution to sheet 4, exercise 2.
# Lecture: Integer Programming in Public Transport  (TU Berlin, winter term 2006/07).
# Ralf Borndoerfer, Christian Liebchen, and Marc Pfetsch
#
# Cost minimizing model for line planning.
#
# Marc Pfetsch
# 11/2006
# ------------------------------------------------------------------------------
#
# The data are:
# - edges   edges of the graph. They are directed and the data should contain
#           forward and backward directions if necessary.
# - load    passenger load on each edge
# - odnodes contains a list of all OD-nodes, which are a subset of the nodes of
#           the graph. We need this for reading and accessing the OD-matrix.
# - od      OD-matrix
# - lines   edges on lines. The lines are taken as undirected, i.e., the file
#           only contains one direction and the model has to take care to check
#           both directions.
## - times  travel times on each edge
# ------------------------------------------------------------------------------
# sets:

# OD-Nodes
# The OD-nodes could be inferred from the edges "E", but
# for reading the OD-matrix we need the correct order.
set O := { read "odnodes.dat" as "<1s>" comment "#" };

# edges (directed!)
set E := { read "edges.dat" as "<1s,2s>" comment "#" };

# line set
set line := { read "lines.dat" as "<1s,2s,3s>" comment "#" };

# line names
set L := proj(line,<1>);

# nodes on lines
set N := proj(line, <1,2>) union proj(line, <1,3>);

do print "Number of OD-nodes:";
do print card(O);

do print "Number of edges:";
do print card(E);

do print "Number of lines:";
do print card(L);

# ------------------------------------------------------------------------------
# parameters:

# load on each edge
param load[E] := read "load.dat" as "<1s,2s> 3n" comment "#";

# OD-matrix
param d[O*O] := read "od.dat" as "n+" comment "#";

# travel times
param time[E] := read "times.dat" as "<1s,2s> 3n" comment "#";

# Parameters for computing the capacity
param mincars := 3;
param carcap := 467;
param maxcars := 5;
param traincap := maxcars * carcap;
param maxfreq := 10;

# costs
param carfixedcost := 353100;
param trainfixedcost := carfixedcost;
param caropercost := 5803;
param trainopercost := 44959;

# compute lower bound on frequency
param lbf[<u,v> in E] := ceil(load[u,v]/traincap);

# turn-around-times:
param tt[O] := <"Ehv"> 5, <"Gn"> 5, <"Gvc"> 5, <"Hgl"> 5, <"Lls"> 5, <"Rsdg"> 5, <"Std"> 5,
	<"Apd"> 13.8, <"Asn"> 13.8, <"Gv"> 13.8, <"Hr"> 13.8, <"Mt"> 13.8, <"Shl"> 13.8, <"Zvg"> 13.8,
	<"Ah"> 14.1, <"Asd"> 14.1, <"Asdz"> 14.1, <"Bd"> 14.1, <"Lw"> 14.1, <"Odzg"> 14.1, <"Rtd"> 14.1, 
	<"Ut"> 14.1, <"Zl"> 14.1;

# compute total traveling time including turn-around-times:
param TT[<l> in L] := sum <l,u,v> in line: time[u,v] + 
	tt[ord({<ll,u,v> in line with ll == l}, 1, 2)]
	+ tt[ord({<ll,u,v> in line with ll == l}, card({<ll,u,v> in line with ll == l}), 3)];

# ------------------------------------------------------------------------------
# sets once again:

# demand pairs
set D := { <s,t> in O * O with d[s,t] > 0 };

# frequencies
set F := {1,2};

# carriages
set C := {mincars to maxcars};

# possible combinations for lines, frequencies, and number of carriages
set R := { <l,f,c> in L*F*C };

# ------------------------------------------------------------------------------
# checks
do forall <l,u,v> in line do check <u,v> in E and <v,u> in E;

# ------------------------------------------------------------------------------
# variables:

# y[l,f,c] whether combination of line l, frequency f, and additional number 
# of carriages c is taken
var y[R] binary;

# ------------------------------------------------------------------------------
# objective

maximize obj: sum <l,f,c> in R do 
		(ceil(f*TT[l]) * (trainfixedcost + carfixedcost * c) +
		(sum <ll,u,v> in line do time[u,v]) * (trainopercost + f * caropercost)) * y[l,f,c];

# ------------------------------------------------------------------------------
# constraints
subto lb: 
	forall <u,v> in E:
	sum <l,f,c> in R with <l,u,v> in line or <l,v,u> in line: f * y[l,f,c] >= lbf[u,v];

subto ub: 
	forall <u,v> in E:
	sum <l,f,c> in R with <l,u,v> in line or <l,v,u> in line: f * y[l,f,c] <= maxfreq;

subto capacity:
	forall <u,v> in E:
	sum <l,f,c> in R with <l,u,v> in line or <l,v,u> in line: carcap * f * c * y[l,f,c] >= load[u,v];

subto partitioning:
	forall <l> in L:
	sum <l,f,c> in R: y[l,f,c] <= 1;