~+;~o|32;~?!i&32768;~o|7       / run in silent mode if switch 15 is off
~/*****CONGTRAS: EMME MACRO FOR CONGESTED TRANSIT ASSIGNMENT (1.17)  ***
~/***  Copyright (C) INRO Consultants Inc., Montreal, Canada 1991-95 ***
~/***  Modified by Yolanda Noriega, CIRRELT, December 2007.          ***
~/***  Modified by Diane Larin, INRO, April 2009.                    ***
~/***                                                                ***
~/***  This Emme macro implements the congested transit assignment 
~/***  based on the Frank+Wolfe algorthm.  For the theoretical
~/***  details see section 5 in Spiess and Florian (1988).
~/***
~/***  Differences when compared to version 1.11:
~/***  - Uses only one scenario and extra attributes to save
~/***    temporary results
~/***  - Module 1.11 is used only at the end of the macro; the
~/***    rest of the computations are done with Module 2.41.
~/***  - The Muller method replaces the secant method in the step 
~/***    length computation.
~/***  - Changed approach to set the default values when preparing
~/***    the transit assignment (to deal with new spread factor).
~/***
~/***  Requirements:
~/***  - All transit time functions must be coded in the form
~/***    (<uncongested transit time>)*(1+us1)
~/***  - ms85 and ms86 contain the weight and exponent of the
~/***      congestion term (must be set by the user)
~/***  - ms83-ms84 & ms87-ms99 are reserved for the use of this macro
~/***      ms83: gradient at second point 
~/***      ms84: second approximation point
~/***      ms87: temporary scalar
~/***      ms88: total number of transit trips
~/***      ms89: current normalized gap
~/***      ms90: iteration counter
~/***      ms92: average transit impedance
~/***      ms93: average minimal transit impedance
~/***      ms94: current (optimal) step length
~/***      ms95: current gradient
~/***      ms96: gradient at first point
~/***      ms97: gradient at third point
~/***      ms98: first approximation point
~/***      ms99: third approximation point
~/***  - ui1, ul1, us1 and us2 attributes are used by the macro;
~/***      their values are temporarily stored in extra attributes 
~/***      and recovered at the end of the macro.
~/***  - Extra attributes @inboa, @fiali, @volax, @voltr, @timtr    
~/***      and @board are created and used during the computations.
~/***      They are deleted at the end of the macro.
~/***  - Segment attribute @ccost is created (if it doesn't exist)
~/***      and used to store the congestion costs.
~/***  - r1 to r12 registers are used by the macro.
~/***
~/***  Calling sequence:
~/***  ~<congtras <ftyp> <iter> <ngap> <gpq> [modes] [btim] [ww] [aw] [bw]
~/***  where
~/***    ftyp  type of congestion function to be used
~/***            (1=BPR, 2=Conical)
~/***    iter  no. of F+W iterations stopping criterion
~/***    ngap  normalized gap stopping criterion
~/***    gpq   transit demand matrix
~/***    modes active modes,          optional (default=all modes)
~/***    btim  boarding time,         optional (default=2.5)
~/***    ww    waiting time weight,   optional (default=2.0)
~/***    aw    auxiliary time weight, optional (default=2.0)
~/***    bw    boarding time weight,  optional (default=2.0)
~/***
~/***  Reference for the Muller method:
~/***    Burden, R. L., Faires, J. D., "Numerical Analysis", Seventh 
~/***    Edition, Books/Cole, Thomson Learning, 2001, p. 95-99. 
~/**********************************************************************
~x=%0%
~?x=0
~$stop
~y=3       /test if enough parameters are provided
~?x<4
~$error
~z=%1%
~%
~y=2       /test if congestion function is valid
~?z<1
~$error
~?z>2
~$error
~:start
off=2
~#set active modes to default if not specified
~+;~t1=%4%;~?t1=;~t0=%t0% *
~#set boarding time to default if not specified
~+;~t1=%5%;~?t1=;~t0=%t0% 2.5
~#set waiting time weight to default if not specified
~+;~t1=%6%;~?t1=;~t0=%t0% 2.0
~#set auxiliary time weight to default if not specified
~+;~t1=%7%;~?t1=;~t0=%t0% 2.0
~#set boarding time weight to default if not specified
~+;~t1=%8%;~?t1=;~t0=%t0% 2.0
c='START congtras %t0%'
~p=2005
~t1=%p%
~+;~p=2009;~r10=%%%p%%%;~p=2010;~r11=%%%p%%%;~p=2011;~r12=%%%p%%%
~/============ CONGESTED TRANSIT ASSIGNMENT =========== 
~/Initial time:         %r10%h%r11%m%r12%s  
~/Date:                 %d% 
~/Scenario:             %s% 
~/Max. iterations:      %1% 
~/Max. norm. gap:       %2% min 
~/Transit demand:       %3% 
~?z=1
~/Congestion function:  BPR 
~?z=2
~/Congestion function:  Conical 
~/Congestion weight:    ms85=%ms85% 
~/Congestion exponent:  ms86=%ms86% 
~?p=1  (UNIX)
~!rm -f congtras.jnk congtras.jnk
~?p=2  (DOS)
~+|~!del congtras.jnk|~!del congtras.rep
reports=congtras.jnk
~# ***** initialize the extra attributes *****
2.42
~+;2;1
~?e
~+|q|~/ *** No space for extra node attribute @ui1!|~$>stop
~+;@ui1;~?e
~+; ;3;@ui1;y;2;1;@ui1
~+;last value of ui1;
~+;2;1
~?e
~+|q|~/ *** No space for extra node attribute @inboa!|~$>stop
~+;@inboa;~?e
~+; ;3;@inboa;y;2;1;@inboa
~+;initial boardings at node i;
~+;2;1
~?e
~+|q|~/ *** No space for extra node attribute @fiali!|~$>stop
~+;@fiali;~?e
~+; ;3;@fiali;y;2;1;@fiali
~+;final alightings at node i;
~+;2;2
~?e
~+|q|~/ *** No space for extra link attribute @ul1!|~$>stop
~+;@ul1;~?e
~+; ;3;@ul1;y;2;2;@ul1
~+;last value of ul1;
~+;2;2
~?e
~+|q|~/ *** No space for extra link attribute @volax!|~$>stop
~+;@volax;~?e
~+; ;3;@volax;y;2;2;@volax
~+;aux. transit volume at link;
~+;2;4
~?e
~+|q|~/ *** No space for extra segment attribute @ccost!|~$>stop
~+;@ccost;~?e
~+; ;3;@ccost;y;2;4;@ccost
~+;congestion cost;
~+;2;4
~?e
~+|q|~/ *** No space for extra segment attribute @us1!|~$>stop
~+;@us1;~?e
~+; ;3;@us1;y;2;4;@us1
~+;last value of us1;
~+;2;4
~?e
~+|q|~/ *** No space for extra segment attribute @us2!|~$>stop
~+;@us2;~?e
~+; ;3;@us2;y;2;4;@us2
~+;last value of us2;
~+;2;4
~?e
~+|q|~/ *** No space for extra segment attribute @timtr!|~$>stop
~+;@timtr;~?e
~+; ;3;@timtr;y;2;4;@timtr
~+;auxiliary segment attribute;
~+;2;4
~?e
~+|q|~/ *** No space for extra segment attribute @voltr!|~$>stop
~+;@voltr;~?e
~+; ;3;@voltr;y;2;4;@voltr
~+;transit volume at segment;
~+;2;4
~?e
~+|q|~/ *** No space for extra segment attribute @board!|~$>stop
~+;@board;~?e
~+; ;3;@board;y;2;4;@board
~+;transit boardings at segment;
q
~# ***** store the original value of the attributes *****
2.41
~+;1;y;@ui1;n;ui1;;*;5;r
~+;1;y;@ul1;n;ul1;;*;5;r
~+;1;y;@us1;n;us1;;*;*;5;r
~+;1;y;@us2;n;us2;;*;*;5;r
~##### initialize US1 to zero
~+;1;y;us1;0;;*;*;5;r
q
3.21 /##### compute total number of trips
1
y
ms88
y
gtrtot
total number of transit trips
~?q=1
y
0
%3%


n
+
+
~?q=2
2
~/Total no. of trips:   %ms88% 
1    /set iteration counter to 0
y
ms90
y
iter
iteration counter
~?q=1
y
0
0


~?q=2
2
q
3.12 /##### change demand matrix to columnwise storage
3
%3%
4
2
q
5.11 /##### prepare initial transit assignment
2         / transit assignment
~?q=2
2
%3%       / demand matrix
ms93      / transit impedances
yes       / initialize
utr%ms90% / matrix name
'congested transit impedances iteration %ms90%' / matrix description
~?q=1
y
0         / default value
          / no in-vehicle time matrix
          / no auxiliary time matrix
          / no waiting time matrix
          / no first waiting time matrix
          / no boarding time matrix
          / no no. of boardings matrix
'%4%'     / active modes
1         / actual headways
1         / same boarding times
%5%       / minutes boarding time (default 2.5)
1         / same wait time factor
0.5       / wait time factor (default 0.5)
%6%       / wait time weight (default 2.0)
%7%       / aux time weight (default 2.0)
%8%       / boarding time weight (default 2.0)
~?q=1
n         / no additional options 
reports=^
reports=congtras.rep
5.31 /##### initial transit assignment
~?b=1
2         / report to printer
reports=^
reports=congtras.jnk
3.21 /##### compute average transit impedance
1
y
ms92
y
uav%ms90%
avg transit impedance iter %ms90%
~?q=1
y
0
ms93


~?q=2
2
q
2.41
1           /now compute the congestion penalties
y
@timtr
n
timtr

all
all
5
r
1           /now compute the congestion penalties
n
voltr*(timtr-dwt)/(1+us1)/ms(88)*ms(85)*
~?z=1
(voltr/(vcapt*60./hdwy))^ms(86)
~?z=2
(sqrt(put(ms(86)*(1-voltr/vcapt/60*hdwy))*get(1)+
~?z=2
put((ms(86)-.5)/(ms(86)-1))*get(2))+1-get(1)-get(2))

all
all
5
4
87
ccost
congestion costs
r
1
n        /0*hdwy+%ms92%+%ms87%
0*hdwy+ms(92)+ms(87)

ind=1

5
4
92


q
~/Avg congestion cost:  %ms87% min 
~/Avg trip impedance:   %ms92% min 
~t2=%z%
~p=2005
~z=%p%
~z-%t1%
~z/10
~t1=%p%
~/CPU Time:             %z% sec
~z=%t2%
~# ***** store the transit assignment results *****
2.41
~+;1;y;@fiali;n;fialii;;*;5;r
~+;1;y;@inboa;n;inboai;;*;5;r
~+;1;y;@volax;n;volax;;*;5;r
~+;1;y;@voltr;n;voltr;;*;*;5;r
~+;1;y;@board;n;board;;*;*;5;r
q
~:nextiter
~x=%ms90%
~?x=%1%
~$alldone
~x+1
~/===================== Iteration %x% ================== 
3.21 /##### increment iteration counter
1
y
ms90
n
ms90+1


~?q=2
2
q
2.41 /##### compute congestion penalties in US1
1           /first get current segment volumes into US2
y
us2
@voltr

all
all
5
r
1
n
((us2-60*vcapt/hdwy).max.0)*length
 
all
all
5
4
87
excess
excess passenger-kms
r
~/Excess pass-kms:      %ms87% 
1           /now compute the congestion penalties
y
us1
ms(85)*
~?z=1
(us2/(vcapt*60./hdwy))^ms(86)
~?z=2
(sqrt(put(ms(86)*(1-us2/vcapt/60*hdwy))*get(1)+
~?z=2
put((ms(86)-.5)/(ms(86)-1))*get(2))+1-get(1)-get(2))

all
all
5
q
5.11 /##### prepare transit assignment iteration %ms90%
2         / transit assignment
~?q=2
2
%3%       / demand matrix
ms93      / transit impedances
yes       / initialize
utr%ms90% / matrix name
'congested transit impedances iteration %ms90%' / matrix description
~?q=1
y
0         / default value
          / no in-vehicle time matrix
          / no auxiliary time matrix
          / no waiting time matrix
          / no first waiting time matrix
          / no boarding time matrix
          / no no. of boardings matrix
'%4%'     / active modes
1         / actual headways
1         / same boarding times
%5%       / minutes boarding time (default 2.5)
1         / same wait time factor
0.5       / wait time factor (default 0.5)
%6%       / wait time weight (default 2.0)
%7%       / aux time weight (default 2.0)
%8%       / boarding time weight (default 2.0)
~?q=1
n         / no additional options 
reports=^
reports=congtras.rep
5.31 /##### transit assignment iteration %ms90%
~?b=1
2         / iteration %ms90% uavg=%ms92% umin=%ms93% lambda=%ms94%
reports=^
reports=congtras.jnk
2.41 /##### compute the gap
~/Avg min. trip imp.:   %ms93% min 
1
n
0*hdwy+ms(93)-ms(92)

ind=1

5
4
96
grafp
gradient at first point
r
1           
n
(timtr-dwt)/(1+us1)*ms(85)*(
~?z=1
(voltr/(vcapt*60./hdwy))^ms(86)
~?z=1
-(us2/(vcapt*60./hdwy))^ms(86)
~?z=2
(sqrt(put(ms(86)*(1-voltr/vcapt/60*hdwy))*get(1)+
~?z=2
put((ms(86)-.5)/(ms(86)-1))*get(2))-get(1))
~?z=2
-(sqrt(put(ms(86)*(1-us2/vcapt/60*hdwy))*get(3)+
~?z=2
get(2)*get(2))-get(3))
)*(voltr-us2)/ms(88)

all
all
5
4
97
gradu
time difference
r
1
n
0*hdwy+ms(97)+ms(93)-ms(92)

ind=1

5
4
97
gratp
gradient at third point
r
1
n
0*hdwy-ms(96)

ind=1

5
4
89
gap%ms90%
normalized gap at iteration %ms90%
r
~/Normalized gap:       %ms89% min 
1
n
0*hdwy

ind=1

5
4
98
fapp
first approximation point
r
1
n
0*hdwy+1

ind=1

5
4
99
tapp
third approximation point
r
~# ***** third point for the Muller step *****
1
n
0*hdwy+0.5

ind=1

5
4
84
sapp
second approximation point
r
1           
n
(timtr-dwt)/(1+us1)*ms(85)*(
~?z=1
((us2+ms(84)*(voltr-us2))/(vcapt*60./hdwy))^ms(86)
~?z=1
-(us2/(vcapt*60./hdwy))^ms(86)
~?z=2
(sqrt(put(ms(86)*(1-(us2+ms(84)*(voltr-us2))/vcapt/60*hdwy))*get(1)+
~?z=2
put((ms(86)-.5)/(ms(86)-1))*get(2))-get(1))
~?z=2
-(sqrt(put(ms(86)*(1-us2/vcapt/60*hdwy))*get(3)+
~?z=2
get(2)*get(2))-get(3))
)*(voltr-us2)/ms(88)

all
all
5
4
83
grasp
gradient at second point
r
1
n
0*hdwy+ms(83)+ms(93)-ms(92)

ind=1

5
4
83
grasp
gradient at second point
r
~x=0
~:nextlambda
~x+1
~# Value of h1
~+;~r1=%ms84%;~r1-%ms98%
~# Value of h2
~+;~r2=%ms99%;~r2-%ms84%
~# Somme h1+h2
~+;~r3=%r1%;~r3+%r2%
~# Value of delta1
~+;~r4=%ms83%;~r4-%ms96%;~r4/%r1%
~# Value of delta2
~+;~r5=%ms97%;~r5-%ms83%;~r5/%r2%
~# Value of d
~+;~r6=%r5%;~r6-%r4%;~r6/%r3%
~# Value of b
~+;~r7=%r2%;~r7*%r6%;~r7+%r5%
~# Value of "4*a*c" for the radical
~+;~r8=%ms97%;~r8*%r6%;~r8*4
~# Value of "b^2" for the radical 
~+;~r9=%r7%;~r9*%r7%
1
n
0*hdwy+(%r9%>%r8%)*sqrt(%r9%-%r8%)

ind=1

5
4
87
D
auxiliary storage for D
r
1
n
0*hdwy+(abs(%r7%-ms(87)).lt.abs(%r7%+ms(87)))*(%r7%+ms(87))
+(abs(%r7%-ms(87)).ge.abs(%r7%+ms(87)))*(%r7%-ms(87))

ind=1

5
4
87
E
auxiliary storage for E
r
1
n
0*hdwy-2*ms(97)/ms(87)

ind=1

5
4
87
h
auxiliary storage for h
r
1
n
0*hdwy+ms(99)+ms(87)

ind=1

5
4
94
lambda
current step length
r
1
n
0*hdwy+abs(ms(87))*100000

ind=1

5
4
87
h
auxiliary storage for h
r
~y=%ms87%
~?y<100
~$stepfound
1           
n
(timtr-dwt)/(1+us1)*ms(85)*(
~?z=1
((us2+ms(94)*(voltr-us2))/(vcapt*60./hdwy))^ms(86)
~?z=1
-(us2/(vcapt*60./hdwy))^ms(86)
~?z=2
(sqrt(put(ms(86)*(1-(us2+ms(94)*(voltr-us2))/vcapt/60*hdwy))*get(1)+
~?z=2
put((ms(86)-.5)/(ms(86)-1))*get(2))-get(1))
~?z=2
-(sqrt(put(ms(86)*(1-us2/vcapt/60*hdwy))*get(3)+
~?z=2
get(2)*get(2))-get(3))
)*(voltr-us2)/ms(88)

all
all
5
4
95
grad
gradient at current lambda
r
1
n
0*hdwy+ms(95)+ms(93)-ms(92)

ind=1

5
4
95
grad
gradient at current lambda
r
~?i&32768
~/L%x% %ms94%/%ms95% %ms98%/%ms96% %ms99%/%ms97% %ms84%/%ms83%
1
n
0*hdwy+ms(84)

ind=1

5
4
98
fapp
first approximation point
r
1
n
0*hdwy+ms(99)

ind=1

5
4
84
sapp
second approximation point
r
1
n
0*hdwy+ms(94)

ind=1

5
4
99
tapp
third approximation point
r
1
n
0*hdwy+ms(83)

ind=1

5
4
96
grafp
gradient at first point
r
1
n
0*hdwy+ms(97)

ind=1

5
4
83
grasp
gradient at second point
r
1
n
0*hdwy+ms(95)

ind=1

5
4
97
gratp
gradient at third point
r
~?x<20
~$nextlambda
~:stepfound
1    / step on the range
n
ms(94) .min. 1 .max. 0*hdwy

ind=1

5
4
94
lambda
current step length
r
1    / compute complement of step size
n
0*hdwy+(1-ms(94))

ind=1

5
4
87
lbar
complement of step length
r
q
~/Optimal step length:  %ms94% after %x% Muller steps
~# ***** update the congested transit assignment results *****
2.41
~+;1;y;@inboa;n;inboai*ms(94)+@inboa*ms(87);;*;5;r
~+;1;y;@fiali;n;fialii*ms(94)+@fiali*ms(87);;*;5;r
~+;1;y;@volax;n;volax*ms(94)+@volax*ms(87);;*;5;r
~+;1;y;@voltr;n;voltr*ms(94)+@voltr*ms(87);;*;*;5;r
~+;1;y;@board;n;board*ms(94)+@board*ms(87);;*;*;5;r
q
2.41
1           
n
(timtr-dwt)/(1+us1)*ms(85)*(
~?z=1
((us2+ms(94)*(voltr-us2))/(vcapt*60./hdwy))^ms(86)
~?z=1
-(us2/(vcapt*60./hdwy))^ms(86)
~?z=2
(sqrt(put(ms(86)*(1-(us2+ms(94)*(voltr-us2))/vcapt/60*hdwy))*get(1)+
~?z=2
put((ms(86)-.5)/(ms(86)-1))*get(2))-get(1))
~?z=2
-(sqrt(put(ms(86)*(1-us2/vcapt/60*hdwy))*get(3)+
~?z=2
get(2)*get(2))-get(3))
)*(us2+ms(94)*(voltr-us2))/ms(88)

all
all
5
4
87
udiff
udiff
q
3.21 /##### update current linear bias
1
y
ms92
y
uav%ms90%
average transit impedance it %ms90%
n
ms(94)*ms(93)+(1-ms(94))*ms(92)+ms(87)


~?q=2
2
q
~/Avg trip impedance:   %ms92% min
2.41 /#####  test gap stopping criteria
1
n
(0*hdwy-ms(89)+%2%)*100000

ind=1

5
4
87
gdiff
normalized gap - stopping critera
q
~t2=%z%
~p=2005
~z=%p%
~z-%t1%
~t1=%p%
~z/10
~/CPU Time:             %z% sec
~z=%t2%
~x=%ms87%
~?x<0
~$nextiter
~:alldone
~# ***** store results in the appropiate attributes *****
~# ***** copy congestion costs into @ccost *****
~+;2.41;1;y;@ccost;n;timtr-@timtr;;*;*;4;q
~# copy results using module 1.11
~# find the internal number of the current scenario, and save it in x
1.11
~+;4;1;1,1,35,35;61,1,87,87;y;r
~x=%ms87%
~# find the max. number of the scenarios, and save it in y
~+;4;1;1,1,51,51;61,1,87,87;y;q
~y=%ms87%
~y+%x%
~# ***** copy @inboa into inboai *****
~+;2.41;1;y;ui1;@inboa;;*;4;q
~+;1.11;4;1;8,%x%;10,%x%;q
~# ***** copy @fiali into fialii *****
~+;2.41;1;y;ui1;@fiali;;*;4;q
~+;1.11;4;1;8,%x%;10,%y%;q
~# ***** copy @volax into volax *****
~+;2.41;1;y;ul1;@volax;;*;4;q
~+;1.11;4;1;16,%x%;19,%x%;q
~# ***** copy @timtr into timtr *****
~+;2.41;1;y;us1;@timtr;;*;*;4;q
~+;1.11;4;1;39,%x%;38,%x%;q
~# ***** copy @voltr into voltr *****
~+;2.41;1;y;us1;@voltr;;*;*;4;q
~+;1.11;4;1;39,%x%;36,%x%;q
~# ***** copy @board into board *****
~+;2.41;1;y;us1;@board;;*;*;4;q
~+;1.11;4;1;39,%x%;37,%x%;q
~/===================================================== 
~p=2004
~?p=1  (UNIX)
~!rm -f congtras.jnk
~?p=2  (DOS)
~!del congtras.jnk
~?i&32768
~$cleanup
~# ***** restore the values of the attributes *****
2.41
~+;1;y;ui1;@ui1;;*;5;r
~+;1;y;ul1;@ul1;;*;5;r
~+;1;y;us1;@us1;;*;*;5;r
~+;1;y;us2;@us2;;*;*;5;r
q
~# ***** delete auxiliary extra attributes *****
2.42
~+;3;@ui1;y
~+;3;@fiali;y
~+;3;@inboa;y
~+;3;@ul1;y
~+;3;@volax;y
~+;3;@us1;y
~+;3;@us2;y
~+;3;@timtr;y
~+;3;@voltr;y
~+;3;@board;y
q
3.12 /##### change demand matrix back to rowwise storage
3
%3%
4
1
q
~$cleanup
~:error
~/  
~?y=2
~/ Congestion function type must be 1 for BPR or 2 for Conical!
~?y=3
~/ Invalid number of macro parameters!
~$stop
~:cleanup
on=2
reports=^
c='end congtras after %ms90% iterations'
~+;~p=2009;~r10=%%%p%%%;~p=2010;~r11=%%%p%%%;~p=2011;~r12=%%%p%%%
~/ 
~/Final time:         %r10%h%r11%m%r12%s  
~/***************************************************************
~/`` ~<congtras %z% %1% %2% %3% %4% %5% %6% %7% %8%''
~/CONGESTED TRANSIT ASSIGNMENT MACRO CONGTRAS TERMINATED NORMALLY
~/CURRENT SCENARIO:      %s%
~/NO. OF TRANSIT TRIPS:  %ms88%
~/NO. OF ITERATIONS:     %ms90%
~/AVG MIN. IMPEDANCE:    %ms93% min
~/AVG TRIP IMPEDANCE:    %ms92% min
~/NORMALIZED GAP:        %ms89% min
~/***************************************************************
~:stop
~o=6