-
Notifications
You must be signed in to change notification settings - Fork 0
/
TD_Calibration_Model_OntarioOnly_HeavyOnly.R
174 lines (145 loc) · 6.23 KB
/
TD_Calibration_Model_OntarioOnly_HeavyOnly.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
library(dplyr)
library(reshape2)
library(tidyr)
library(foreign)
###############################################################################
#' Title: "Calibrating a Doubly Constrained Gravity Model - Ontario Only"
#' Programmer: "Mausam Duggal, Systems Analysis Group, WSP|PB"
#' Date: "January 11, 2017"
###############################################################################
#' This code was written to calibrate the beta value in a doubly constrained
#' for the empty truck model. The data used was the 2012 CVS data. Bryce
#' prepared the empty truck matrices for both heavies and medium trucks.
#' The trucks have been currently combined because the commodity flow model is
#' not currently distinguishing between heavy and medium truck tonnage. Further,
#' while Rolf's Empty Truck model can be developed for different truck types it
#' was not deemed necessary in this version of the model.
###############################################################################
########################### INPUTS ############################################
#' read observed truck trip matrix and combine
obs_hvy <- read.csv("C:/Projects/Province-wide/Empties/empty_heavy_truck_matrix.csv",
stringsAsFactors = FALSE)
#' Combine all empties and also cast it in a format desired and getting rid of
#' intraCSD trips.
all_obs <- transform(obs_hvy, Flag = paste0(Orig,Dest)) %>%
subset(., Orig != Dest)
#' wide format and set NAs to zero
all_obs_cast <- dcast(all_obs, Orig ~ Dest, value.var = "Trips")
rownames(all_obs_cast) <- all_obs_cast[, 1]
all_obs_cast <- all_obs_cast[, 2:ncol(all_obs_cast)]
all_obs_cast[is.na(all_obs_cast)] <- 0
#' get total trips
totals <- sum(apply(all_obs_cast, 1, function(x) sum(x)))
write.csv(all_obs_cast, "C:/Projects/Province-wide/Empties/dome_hv.csv", row.names = FALSE)
#' batch in distance skims and cast it in a format desired. There is also a need
#' to only keep distance skims for those O-D pairs that are in the combined trips.
#' The diagonal value of the distance matrix is set to 0 to avoid making that
#' cell receive any values during balancing
dist <- read.csv("C:/Projects/Province-wide/Empties/od_road_distances.csv",
stringsAsFactors = FALSE) %>%
transform(., Flag = paste0(Orig, Dest))
#' only keep those o-D pairs that match the all_obs truck trip table
dist1 <- subset(dist, dist$Flag %in% all_obs$Flag)
#' wide format and set NAs to zero
dist1_cast <- dcast(dist1, Orig ~ Dest, value.var = "Dist")
rownames(dist1_cast) <- dist1_cast[, 1]
dist1_cast <- dist1_cast[, 2:ncol(dist1_cast)]
dist1_cast[is.na(dist1_cast)] <- 0
#' set diagonal to global value of 10
#diag(dist1_cast) <- 0
####################### COST MATRIX ###########################################
#' generate a cost matrix using distance between CSDs and set the beta to -0.01
#' Rolf uses a very small beta of -0.001. This will have a tendency to require
#' longer to calibrate and also it does not create enough variation in the
#' value of the various interchanges, especially those that are far apart.
cost <- dist1_cast
initbeta <- -0.008
################ COST FUNCTION MATRIX #########################################
#' Cost Function Matrix. We are using an exponential function here. Other formulations
#' can be tested as well, if desired. All those cells that had a value of 0
#' would automatically get set to 1 after taking the exponent. This would
#' incorrectly make these interchanges a candidate to receive trips thereby
#' making it difficult to calibrate and converge the betas. These cells are
#' forced to zero.
cfunc <- exp(initbeta*cost)
cfunc[cfunc == 1] <- 0
#' Estimate observed average trip lengths
tlfd_obs <- sum(dist1_cast * all_obs_cast)/totals
print(tlfd_obs)
############################ FURNESSING ##############################
#' This is the furnessing function, which is the backbone of the calibration
#' procedure.
furness <- function(cmat){
"
This function conducts furnessing or two dimensional balancing
Inputs: observed truck matrix
Arguments: cost function
Results: balanced matrix
"
#' row sum of cost function matrix
skim_trow <- as.data.frame(t(colSums(cmat)))
#' row sum of observed matrix
obs_tatt <- as.data.frame(t(colSums(all_obs_cast)))
#' row balancing factors
rowbal <- obs_tatt/skim_trow
rowbal1 <- do.call("rbind", replicate(nrow(all_obs_cast), rowbal,
simplify = FALSE))
#first iteration
cfunc1 <- cmat*rowbal1
#' col sum of cost function matrix
skim_tcol <- as.data.frame(rowSums(cfunc1))
#' col sum of observed matrix
obs_tprod <- as.data.frame(rowSums(all_obs_cast))
#' column balancing factors
colbal <- obs_tprod/skim_tcol
colbal1 <- do.call("cbind", replicate(ncol(all_obs_cast), colbal,
simplify = FALSE))
# next iteration
cfunc1 <- cfunc1*colbal1
return(cfunc1)
}
############################### CALIBRATION PROCESS ###########################
#' the outer loop is to calculate revised beta values. The user can set the
#' number of outer loops one wants to evaluate.
outer_iter <- 20
#' Start outer loop
for (j in 1:outer_iter) {
################## BETA COEFFICIENTS ########################################
#' set initial beta coefficient and cost matrix
if(j == 1){
beta <- initbeta
print(beta)
#' compute cfunc again
cfunc <- exp(initbeta*cost)
cfunc[cfunc == 1] <- 0
} else {
beta <- beta1
print(beta1)
#' compute cfunc again
cfunc <- exp(beta*cost)
cfunc[cfunc == 1] <- 0
}
#' run the loop with user defined iterations
maxiter <- 15
i <- 1
#' run the first iteration of furnessing
for (i in 1:1){
t <- furness(cfunc)
}
#' now loop furnessing through till the maxiterations
for (i in 2:maxiter){
t1 <- furness(t)
t <- t1
#print(i)
}
###############################################################################
#' Estimate simulated average trip lengths
tlfd_sim <- sum(dist1_cast * t1)/totals
#' calculate revised beta value to feed back to outer loop
if (j == 1){
beta1 <- initbeta * (tlfd_sim/tlfd_obs)
} else {
beta1 <- beta1 * (tlfd_sim/tlfd_obs)
}
}
write.csv(t, "/Projects/Province-wide/Empties/t.csv")