~?!i&32768 ~o=263 ~/*****CONGTRAS: EMME/2 MACRO FOR CONGESTED TRANSIT ASSIGNMENT (1.11)*** ~/*** Copyright (C) INRO Consultants Inc., Montreal, Canada 1991-95 *** ~/*** *** ~/*** This EMME/2 macro implements the congested transit assignment ~/*** based on the Frank+Wolfe algorthm. For the theoretical ~/*** details see section 5 in Spiess and Florian (1988). ~/*** ~/*** Requirements: ~/*** - All transit time functions must be coded in the form ~/*** ()*(1+us1) ~/*** - ms85 and ms86 contain the weight and exponent of the ~/*** congestion term (must set by the user) ~/*** - us1, us2, ms87-ms99 are reserved for the of of this macro ~/*** us1: congestion penalty for segment ~/*** ms87: temporary scalar ~/*** ms88: total number of transit trips ~/*** ms89: current normalized gap ~/*** ms90: iteration counter ~/*** ms91: macro containing scenario pointers ~/*** ms92: avg transit impedance ~/*** ms93: avg minimal transit impedance ~/*** ms94: current (optimal) step length ~/*** ms95: current gradient ~/*** ms96: gradient at lower limit ~/*** ms97: gradient at upper limit ~/*** ms98: lower limit of lambda interval ~/*** ms99: upper limit of lambda interval ~/*** ~/*** Calling sequence: ~/*** ~ ~/*** where ~/*** type of congestion function to be used ~/*** (1=BPR, 2=Conical) ~/*** no. of F+W iterations stopping criterion ~/*** normalized gap stopping criterion ~/*** transit demand matrix ~/*** active modes ~/*** boarding time ~/*** waiting time weight ~/*** auxiliary time weight ~/*** boarding time weight ~/*** ~/********************************************************************** ~x=%0% ~?x=0 ~$stop ~y=3 /test if enough parameters are provided ~?x<5 ~$error ~z=%1% ~% ~y=1 /test if scenario number is 4 digits or less ~?s>9999 ~$error ~y=2 /test if congestion function is valid ~?z<1 ~$error ~?z>2 ~$error ~x=%s% ~x+90000 s=%x% / test if scenario %x% exists already ~?e ~$start ~x-90000 s=%x% ~y=4 ~$error ~:start off=2 c='START congtras %1% %2% %3% %4% %5% %6%' ~p=2005 ~t1=%p% ~/============ CONGESTED TRANSIT ASSIGNMENT =========== ~/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 2.41 /##### initialize US1 to zero 1 y us1 0 all all 5 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% 2.5 / minutes boarding time (default 2.5) 1 / same wait time factor 0.5 / wait time factor (default 0.5) %6% 2.0 / wait time weight (default 2.0) %7% 2.0 / aux time weight (default 2.0) %8% 2.0 / 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 avg transit impedance 1 y ms92 y uavg%ms90% avg transit impedance iter %ms90% ~?q=1 y 0 ms93 ~?q=2 2 q 2.41 1 /now compute the congestion penalties y us2 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 6 91 scenp scenario pointers 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% 1.22 /##### create a temporary copy copy of scenario 3 %s% ~x=%s% ~x+90000 %x% ~?e ~+|q|~/No space available for creating temporary scenario %x%!|~$>stop temp. scen. for congested transit assignment y q 1.11 /##### get scenario pointers into ms91 4 1 1,1,35,35 61,1,87,87 y 1 1,1,51,51 61,1,91,91 y 3 61,1,87,87 61,1,91,91 100,1,0 q ~x=%s% ~x-90000 s=%x% 1.11 /##### get scenario pointers 4 1 1,1,35,35 61,1,87,87 y 3 61,1,87,87 61,1,91,91 10000,1,0 q ~x+90000 s=%x% /Scenario pointers: %ms91% ~: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 2 /first get current segment volumes into US2 ~x=%s% ~x-90000 %x% voltr us2 all all 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% 2.5 / minutes boarding time (default 2.5) 1 / same wait time factor 0.5 / wait time factor (default 0.5) %6% 2.0 / wait time weight (default 2.0) %7% 2.0 / aux time weight (default 2.0) %8% 2.0 / 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 gradl gradient at lower limit 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 gradu gradient at upper limit 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 lambdl lower bound of step length r 1 n 0*hdwy+1 ind=1 5 4 99 lambdu upper bound of step length r ~x=0 ~:nextlambda 1 n / alternate between secant and bisection steps ~x+1 ms(98)-(ms(99)-ms(98))/(ms(97)-ms(96))*ms(96) .min. 1 .max. 0*hdwy ind=1 5 4 94 lambda current step length r 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% 1 n 0*hdwy+((ms(94)-ms(98)).min.(ms(99)-ms(94)))*100000 ind=1 5 4 87 ldiff current step length delta r ~y=%ms87% ~?y<100 ~$stepfound 1 n 0*hdwy+ms(98)+(ms(95).le.0.)*(ms(94)-ms(98)) ind=1 5 4 98 lambdl lower bound of step length r 1 n 0*hdwy+ms(99)+(ms(95).ge.0.)*(ms(94)-ms(99)) ind=1 5 4 99 lambdu upper bound of step length r 1 n 0*hdwy+ms(96)+(ms(95).le.0.)*(ms(95)-ms(96)) ind=1 5 4 96 grad1 gradient at lower bound r 1 n 0*hdwy+ms(97)+(ms(95).ge.0.)*(ms(95)-ms(97)) ind=1 5 4 97 grad2 gradient at upper bound r ~?x<20 ~$nextlambda ~:stepfound 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% secant steps ~x=%s% ~x-90000 1.11 4 / update internal files 10,19,36,37 3 ~y=%ms91% ~y/100 ~y%100 10,%y% ~y=%ms91% ~y/10000 10,%y% %ms94%,%ms87%,0 3 ~y=%ms91% ~y/100 ~y+%ms91% ~y%100 10,%y% ~y=%ms91% ~y/10000 ~y+%ms91% ~y%100 10,%y% %ms94%,%ms87%,0 3 19,-%s% 19,-%x% %ms94%,%ms87%,0 3 36,-%s% 36,-%x% %ms94%,%ms87%,0 3 37,-%s% 37,-%x% %ms94%,%ms87%,0 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 uavg%ms90% average transit impedance iter %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 ~x=%s% ~?x>90000 ~x-90000 s=%x% 2.41 1 /compute the congestion costs and store into US1 y us1 timtr-us2 all all 5 q 1.11 4 / copy uncongested times back from US1 back into TIMTR 1 ~y=%ms91% ~y/10000 ~y+%ms91% ~y%100 39,%y% /source is US2 ~y=%ms91% ~y/10000 /destination is TIMTR 38,%y% q ~/===================================================== ~p=2004 ~?p=1 (UNIX) ~!rm -f congtras.jnk ~?p=2 (DOS) ~!del congtras.jnk ~x+90000 ~?i&32768 ~$cleanup 1.22 /##### delete the temporary scenario %x% 2 %x% y q 3.12 /##### change demand matrix back to rowwise storage 3 %3% 4 1 q ~$cleanup ~:error ~/ ~?y=1 ~/ Scenario number must not exceed 9999! ~?y=2 ~/ Congestion function type must be 1 for BPR or 2 for Conical! ~?y=3 ~/ Invalid number of macro parameters! ~x=%s% ~x+90000 ~?y=4 ~/ Scenario %x% exists already! ~$stop ~:cleanup on=2 reports=^ c='end congtras after %ms90% iterations' ~/ ~/*************************************************************** ~/`` ~