~o|16 ~?!i&32768 ~o=39 ~$>end_of_copyright ~/!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ~/Copyright (C) Heinz Spiess, CH-2558 Aegerten, 1989-99. All rights reserved. ~/ ~/The right to use this macro is granted to all EMME/2 users, provided the ~/following conditions are met: ~/ 1) The macro cannot be sold for a fee (but it can be used and distributed ~/ without charge within consulting projects). ~/ 2) The user is aware that this macro is not a part of the EMME/2 software ~/ licence and there is no explicit or implied warranty or support ~/ provided with this macro. ~/ 3) This macro must not be modified without also changing its name. ~/ 4) The comments in this macros must not be removed and any additions or ~/ modification must be appropriately identified as such and give at least ~/ the date, the name and the reason of the modification. ~/!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ~:end_of_copyright ~# ~# Revision history: ~# 99-07-14 Version 2.2: ~# - Correction of a small bug which slipped into version 2.1 and made ~# the macro crash when being used without turn counts. (This bug ~# was found and reported by Zvi Leve - Thank you, Zvi!) ~# - Correction of two minor bugs regarding the case when a turn attribute ~# is specified which does not contain any counts. ~# ~# 99-07-13 Version 2.1: ~# - Support for counts on turns added. In order to keep the calling ~# sequence compatible with earlier versions, the turn attribute ~# containing the counts on turns is also specified in the ~# argument, separated from the link count attribute by a "/". ~# - R^2 is now reported on screen and in the summary file. Note that the ~# formula for R^2 is the one without correction for the mean (no constant ~# term!), thus it is *not* the same as the R^2 given in module 2.43. ~# - The summary file is now named "demadj.sum". It no longer uses command ~# espace commands ~!..., but the new write file feature. ~# - The number of link and turn counts is displayed and a tests were added ~# to detect negative values in the count data. ~# - Customization information can now be passed to DEMADJ in the registers ~# t3-t7. This avoids the need to modify the macro itself - instead it ~# is simply called via a wrapper macro which set the registers accordingly ~# before calling macro DEMADJ. ~# - Scalars ms91 - ms94 are no longer needed, floating point registers ~# r91 - r94 are now used instead. ~# - The macro will no longer work with releases older than Rel.9. ~# ~/***************************************************************************** ~/********* DEMADJ 2.2 - DEMAND ADJUSTMENT USING GRADIENT METHOD *********** ~/********* Copyright (C) Heinz Spiess, CH-2558 Aegerten, 1989-99 *********** ~/ ~/ This macro implements the gradient method for the O-D matrix adjustment ~/ problem. Using observed volumes on a subset of links and a starting O-D ~/ matrix, the model will adjust the matrix to better fit with the observed ~/ volumes. The "steepest descent" property of the gradient approach ensures ~/ that the O-D matrix is not changed more than necessary. In a nut shell, ~/ here is a summary of the computations performed by this macro at each ~/ gradient step: ~/ - auto assignment to get the link volumes 5.11/5.21 ~/ - computation link derivatives and objective funct. 2.41 ~/ - generate scattergram observed vs predicted volumes 2.43 ~/ - addl options assignment to compute gradient matrix 5.11/5.21 ~/ - multiply gradient matrix with weights (optional) 3.21 ~/ - addl options assignment to assign gradient matrix 5.11/5.21 ~/ - compute maximal gradient 3.21 ~/ - compute optimal step length 2.41 ~/ - update demand matrix and step counter 3.21 ~/ ~/ This method has shown to work very well and, using this macro, is very ~/ simple to apply. However, as all matrix adjustment methods, it is extremely ~/ important to be sure that the differences between observed and assigned ~/ volumes are indeed due to the demand matrix and not a result of a coding ~/ errors in the network, badly calibrated volume delay functions or wrong ~/ counts. Very careful analysis of the data before and after the adjustment ~/ is of utmost importance!!! ~/ ~/ Reference: Spiess, H., "A GRADIENT APPROACH FOR THE O-D MATRIX ADJUSTMENT ~/ PROBLEM", Publication 693, CRT, University of Montreal, 1990. ~/ ~/ Data requirements: ~/ @ltmp, @ptmp: The extra attributes @ltmp and @ptmp are used for ~/ temporary storage of the link, resp. turn volume differences. ~/ If these attributes already exist, their contents will be lost. ~/ ms90: This scalar contains the current adjustment step number. If it ~/ exists, it must be initialized to 0 before a new adjustment ~/ is started. This scalar will be incremented after each adjustment ~/ step. It is compared with the stopping criterion at each ~/ gradient step. ~/ ~/ Other requirements: ~/ 1)This macro must be started at the main menu. ~/ 2)This macro needs a release 9.0 or later of the EMME/2 system. ~/ 3)This macro produces for each gradient step a plot file named "scatNN.plt" ~/ containing the observed vs. assigned scattergram. ~/ 4)For each gradient step a one-line summary is written to the file ~/ demadj.sum containing the value of the objective function Z, the maximum ~/ gradient M, the optimal step length L and the R^2 measure of fit. ~/ 5)The demand adjustments can be restricted to a subset of O-D pairs by ~/ specifying an adjustability matrix in register t4 which contains a 1 ~/ if the demand of an OD-pair can be adjusted and 0 if it should not be ~/ changed. (Value of t4 can be defined in an upper level wrapper macro.) ~/ 6)If needed, optional weight factors for the counts can be specified ~/ as an arbitrary link expression in text registers t6/t7. ~/ (Values of t6/t7 can be defined in an upper level wrapper macro.) ~/ 7)The R^2 measure reported in the progress report and the summary file is ~/ based on a linear regression *without* a constant term. ~/ ~/ Usage: ~ [] ~/ ~/ Macro Parameters: ~/ : Link attribute containing link counts (0=no count), e.g. "ul1". ~/ An optional turn attribute containing counts on turns can be ~/ specified optionally after "/" (no blanks!), e.g. "ul2/up2". ~/ If counts are on turns only, use "0" for the link attribute. ~/ : Demand matrix to be adjusted (e.g. "mf6"). Since the initial ~/ contents of this matrix will be overwritten with the resulting ~/ adjusted matrix, it is recommended to first do a copy. ~/ This matrix must not be write protected. ~/ : Temporary matrix of type mf which is used to hold the gradient ~/ direction matrix. Must not be write protected. If this matrix ~/ already exists, its original contents will be lost. ~/ : Number of iteration stopping criteria used in assignment. ~/ : Relative gap stopping criteria used in assignment. ~/ : Normalized gap stopping criteria used in assignment. ~/ : Number of gradient steps to be performed (adjustment ~/ iterations). ~/ : Optional argument. If set to 1, the first "half step" (reg. ~/ assignment to compute volumes) is skipped, assuming that the ~/ assignment was already performed. Use only when continuing ~/ with more step on the same adjustment (i.e. higher ~/ value). ~/ ~x=%0% ~?x<7 ~$STOP ~p=50 ~?p<9 ~+|~/EMME/2 Release %p% - this version of DEMADJ needs Release 9 or later!|~$STOP ~?!m=000 ~+|~/Macro DEMADJ must be started from the main menu of EMME/2!|~$STOP ~?!q=0 ~+|~/Macro DEMADJ must be started from the main menu of EMME/2!|~$STOP ~t1=%1% ~z=0 ~:SPLIT-NEXT ~+|~t2=%%%t1.-%z%%%%|~t2=%%%t2.1%%%|~?t2=|~$SPLIT-DONE|~z+1|~?!t2=/|~$SPLIT-NEXT ~+|~t2=%%%t1.-%z%%%%|~z-1 ~+|~t1=%%%t1.%z%%%%|~?z=0|~t1=| ~:SPLIT-DONE ~?t1= ~+|~/No link counts specified!|~$STOP ~z=1 ~+|~t9=%%%t1.-%z% ~# Registers t3 to t7 contain various configuration parameters. Usually, ~# they do not need to be changed. Should a change ever become necessary, ~# please do not change the default in the macro itself, rather set the ~# corresponding registers in a wrapper macro which then calls DEMADJ, ~# e.g.: ~# ~t4=mf4.ne.0 ~# ~t6=@lwght ~# ~t7=@pwght ~# ~||" ~# ~# Set assignment type to default if it is not defined in calling macro: ~?t5= ~t5=1 ~# Registers t6/t7 contain optional link/turn weight factors for count posts. ~# By default these registers are empty, which implies that all count posts ~# have an implicit weight of 1. If needed t6 (t7) can be set to an arbitrary ~# link (turn) expression which provides positive weight factors for all count ~# post links (e.g. "t6=@cwght" if the weights are stored in the extra link ~# attribute @cwght). ~# ~# Set link weights to default value if they are not defined in calling macro: ~?t6= ~t6= ~# Set turn weights to default value if they are not defined in calling macro: ~?t7= ~t7= ~% ~/***************************************************************************** ~+|2.41|1|n|put(get(1)+((%t1%)>0))+put(get(2)+((%t1%)<0))| ~?e ~+|~/Invalid link counts (%t1%)!||q|~$STOP ~+|all|5|~r10=%%%g1%%%|~r11=%%%g2%%%|q ~/Link counts: %t1% (contains counts on %r10% links) ~?r11>0 ~+|~/%r11% negative link counts detected!|~$STOP ~?t2= ~$>NO TURN COUNTS ~+|2.41|1|n|0*tpf+put(get(1)+((%t2%)>0))+put(get(2)+((%t1%)<0))| ~?e ~+|~/Invalid turn counts (%t2%)!||q|~$STOP ~+|all|all|5|~r12=%%%g1%%%|~r13=%%%g2%%%|q ~/Turn counts: %t2% (contains counts on %r12% turns) ~?r13>0 ~+|~/%r13% negative turn counts detected!|~$STOP ~:NO TURN COUNTS ~x=%r10% ~x+%r12% ~/Current demand: %1% ~/Demand direction: %2% ~?!t4= ~/Adjustability factor: %t4% ~?!t6= ~/Countpost weight factors: %t6% (on links) ~?!t7= ~/ %t7% (on turn) ~/Assignment: %3% iterations, %4% % rel.gap, %5% norm.gap ~/Gradient steps: %ms90% - %6% (continue: %7%) c=demadj %t1%/%t2% %1% %2% %3% %4% %5% %6% %7% ~x=%7% 0 ~?x>0 ~$CONTINUE STEP 3.12 /##### check if matrices exist and create them if necessary ############## 1 /create scalar ms90, if it does not already exist ms90 ~?e ~$>CHECK DPQ step matrix adjustment step counter 0 1 /create full matrix dpq, if it does not already exist ~:CHECK DPQ %2% ~?e ~$>MATRIX CHECKS DONE dpq-0 demand gradient 0 ~:MATRIX CHECKS DONE ~?q=0 q 2.42 /##### check if extra attribute @ltmp exists 3 / remove @ltmp if it exists @ltmp ~?e ~?q=1 yes 2 2 ~?e ~+|q|~/No space for extra link attribute @ltmp!|~$STOP @ltmp temporary link attribute (DEMADJ.MAC) 0 ~?t2= ~$>NO PTMP 3 / remove @ptmp if it exists @ptmp ~?e ~?q=1 yes 2 5 ~?e ~+|q|~/No space for extra turn attribute @ptmp!|~$STOP @ptmp temporary turn attribute (DEMADJ.MAC) 0 ~:NO PTMP q ~:NEXT STEP ~/ ~/Gradient step %ms90%: 5.11 /##### prepare auto assignment to compute auto volumes ################### ~/.... Assign demand to get link and turn volumes 1 / fixed demand auto assignment ~?q=2 2 / new assignment %t5% / type of assignment (time based or generalized cost) 1 / no additional volumes %1% / demand / no auto occupancy matrix / no addl demand / no auto times saved ~?q=1 y / initialize auto volumes %3% / max iterations %4% / % gap %5% / time diff. 5.21 /##### perform auto assignment to compute auto volumes ################### ~?b=1 2 2.41 /##### compute R^2 observed vs predicted volumes ######################### 1 / network calculation - compute R^2 n put(get(1)+%t3%*(%t1%>0)*(%t1%)*volau)+ put(get(2)+%t3%*(%t1%>0)*(%t1%)*(%t1%))+ put(get(3)+%t3%*(%t1%>0)*volau*volau) all 4 / none ~?t2= ~$>NO TURN R2 1 / network calculation - compute turn gradient n put(get(1)+%t3%*(%t2%>0)*(%t2%)*pvolau)+ put(get(2)+%t3%*(%t2%>0)*(%t2%)*(%t2%))+ put(get(3)+%t3%*(%t2%>0)*pvolau*pvolau) all all 4 / none ~:NO TURN R2 ~# g1=%g1% g2=%g2% g3=%g3% ~r2=%g1% ~r2*%g1% ~r2/%g2% ~r2/%g3% ~/.... R^2 observed vs predicted volumes is %r2.5% q 2.41 /##### compute gradient and objective function ########################### 1 / network calculation - compute link gradient y @ltmp y link gradient - macro demadj (it %ms90%) %t3%*(%t1%>0)*(%t1%-volau) ~?!t6= *(%t6%) all 4 / none ~?t2= ~$>NO TURN GRADIENT 1 / network calculation - compute turn gradient y @ptmp y turn gradient - macro demadj (it %ms90%) %t3%*(%t2%>0)*(%t2%-pvolau) ~?!t7= *(%t7%) all all 4 / none ~:NO TURN GRADIENT 1 / network calculation - compute objective function n put(get(puti(91))+(%t1%>0)*(%t1%-volau)*(%t1%-volau)*%t3% ~?!t6= *(%t6%) ) all 4 / no report ~r3=%g91% ~?t2= ~$>NO TURN OBJECTIVE 1 / network calculation - compute objective function n put(get(puti(91))+ (%t2%>0)*(%t2%-pvolau)*(%t2%-pvolau)*%t3% ~?!t7= *(%t7%) ) all all 4 / no report ~:NO TURN OBJECTIVE ~r91=%g91% q ~r4=%r91% ~r4-%r3% ~/.... Objective function at step %ms90% is %r91% ~?!t2= ~/ (%r3% on links, %r4% on turns) on=1 plots=scat%ms90%.plt 2.43 /##### plot scattergram observed vs predicted volumes #################### 2 / link scattergram %t1% / observed volumes vs ~?e ~+||~$>NO LINK SCATTERGRAM volau / predicted volumes at step %ms90% n ~?q=1 n / no symbol index ~?q=1 n / no color index i=%ms90% and 0 not %t1%=0 y / compute regression ~?q>0 ~$>NO LINK SCATTERGRAM 0 / default x range maximum 0 / default y range maximum ~?q=2 2 ~:NO LINK SCATTERGRAM ~?t2= ~$>NO TURN SCATTERGRAM 5 / turn scattergram %t2% / observed volumes vs ~?e ~+||~$>NO TURN SCATTERGRAM pvolau / predicted volumes at step %ms90% n ~?q=1 n / no symbol index ~?q=1 n / no color index i=%ms90% or not 0 i=%ms90% or not 0 .00001,99999999 y / compute regression ~?q>0 ~$>NO TURN SCATTERGRAM 0 / default x range maximum 0 / default y range maximum ~?q=2 2 ~:NO TURN SCATTERGRAM q off=1 plots=^ ~:CONTINUE STEP ~x=%6% ~?x<%ms90% ~$STOP 5.11 /##### prepare add option assignment to compute gradient matrix ########## ~/.... Compute gradient matrix 1 / fixed demand auto assignment ~?q=2 2 / new assignment %t5% / type of assignment (time based or generalized cost) 5 / additional options %1% / demand matrix / no auto occupancy / no add. demand = auto demand / no auto times ~?q=1 y / initialize volumes ~?t2= ~+|6|@ltmp / addl path attribute just on links ~?!t2= ~+|7|@ltmp|@ptmp / addl path attribute on links and on turns + / addl path operator -99999,99999 / addl path threshold %2% / addl attribute matrix (4.3 and later) y D-%ms90% / matrix name demand gradient step %ms90% ~?q=1 y / initialize matrix 0 / to zero 3 / save weighted addl path attributes in matrix ~?q=1 y / initialize volumes %3% / number of iterations %4% / relative gap %5% / normalized gap 5.21 /##### perform add. opt. assignment to obtain gradient matrix ############ ~?q=2 2 ~?t4= ~$>No adjustability factor specified in register t4 3.21 /##### multiply gradient matrix with adjustability factor - if applicable 1 / matrix calculations y / save result %2% / gradient matrix n / don't change header %2%*(%t4%) / no constraint matrix n / no submatrix ~?q=2 2 / report to print file q ~:No adjustability factor specified in register t4 5.11 /##### prepare addl options assignment to get volume direction ########### ~/.... Assign gradient matrix to obtain volume direction 1 / fixed demand auto assignment 2 / new assignment %t5% / type of assignment (time based or generalized cost) 5 / addl options %1% / demand matrix / no auto occupancy %2% / demand gradient matrix / no auto times ~?q=1 y / initialize auto volumes 5 / no addl link attribute %3% / number of iterations %4% / relative gap %5% / normalized gap 5.21 /##### perform addl options assignment to get volume direction ########### ~?q=2 2 3.21 /##### compute maximum gradient ########################################## 1 / matrix calculation n / save result put(get(puti(92)).max.abs(%2%/%1%)) %1% 0,0,ex n .max. .max. ~?q=2 2 ~r92=%g92% q ~/.... Maximum gradient is %r92% 2.41 /##### compute denominator of optimal step length ################## 1 / compute denominator of optimal step (link part) n put(get(puti(94))+(%t1%>0)*volad*volad/%r92%) all 4 ~?t2= ~$>NO TURNS 1 1 / compute denominator of optimal step (turn part) n put(get(puti(94))+(%t2%>0)*pvolad*pvolad/%r92%) all all 4 ~:NO TURNS 1 1 / compute optimal step length now n put(get(puti(93))+(%t1%>0)*(%t1%-volau)*volad/get(94) ~?!t6= *(%t6%) ) all 4 ~?t2= ~$>NO TURNS 2 1 / compute optimal step length now n put(get(puti(93))+ (%t2%>0)*(%t2%-pvolau)*pvolad/get(94) ~?!t7= *(%t7%) ) all all 4 ~:NO TURNS 2 ~r93=%g93% q ~/.... Optimal step length is %r93% 3.21 /##### update demand matrix ############################################## 1 / matrix calculations y / save result %1% / in y / modify header gpq%ms90% Demand at step %ms90% ~?q=1 n %1%+(%r93%.min.1)*%2%/%r92% ~/.... Update demand: %1% := %1%+(%r93%.min.1)*%2%/%r92% n ~?q=2 2 ~?m=321 q 3.21 /##### increment step counter ms90 ####################################### 1 / matrix calculations ~>>demadj.sum ~"I%ms90%: Z%r91% M%r92% L%r93% (R^2=%r2%) ~> y / save result ms90 / in scalar 90 n / don't change header ms90+1 ~?q=2 2 ~?m=321 q c= I%ms90%: Z%r91% M%r92% L%r93% ~$NEXT STEP ~:STOP ~/***************************************************************************** ~o=6 ~?b=2 q