+ Show code- Hide code
library(OpasnetUtils)
library(ggplot2)
library(reshape2)
objects.latest("Op_en5136", code_name = "initiate")
# Trip aggregator, sampled passenger data as input
n.intervals.per.h <- 5
trips.next <- data.frame()
trips.left <- data.frame()
trips.out <- data.frame()
trips.secondary <- data.frame()
times <- seq(7, 9, 1 / n.intervals.per.h)
times[length(times)] <- 1
speed <- 40 # kmh^-1
tdelay <- distance[distance$Mode == "Car" & distance$Time == "1", c("From", "To", "Result")]
tdelay$Result <- ceiling((tdelay$Result / speed) * n.intervals.per.h) / n.intervals.per.h
colnames(tdelay)[colnames(tdelay) %in% "Result"] <- "Delay"
for (i in 1:(length(times) - 2)) {
if(i == 1) {
trips.sample.1 <- trips[trips$Time == times[1],]
trips.sample.1$Secondary <- 0
} else {
trips.sample.1 <- trips.sample.2
trips.sample.1 <- merge(trips.sample.1, trips.secondary, all.x = TRUE)
trips.sample.1$Secondary[is.na(trips.sample.1$Secondary)] <- 0
}
trips.sample.2 <- trips[trips$Time == times[i + 1],]
# Optimizer main code
optimal.d.trips <- 0
sub.optimal.d.trips <- 0
optimal.d.trips <- (trips.sample.1$Result + trips.sample.1$Secondary) %/% 4 * 4
sub.optimal.d.trips <- ifelse(trips.sample.1$Secondary - optimal.d.trips > 0, trips.sample.1$Result + trips.sample.1$Secondary - optimal.d.trips, 0)
busiest <- tapply(trips.sample.2$Result, trips.sample.2$From, sum)
busiest <- sort(busiest, decreasing = TRUE)
condition <- trips.sample.1$Result + trips.sample.1$Secondary - optimal.d.trips - sub.optimal.d.trips > 0
trips.next <- merge(trips.sample.1[condition, c("From","To")], roads[,c("From","To","Through")], all.x = TRUE)
trips.next$Through <- match(trips.next$Through, names(busiest))
checkpoints <- tapply(trips.next$Through, trips.next[,c("From", "To")], min)
trips.sample.1 <- merge(trips.sample.1, as.data.frame(as.table(checkpoints)), all.x = TRUE)
colnames(trips.sample.1)[colnames(trips.sample.1) == "Freq"] <- "Checkpoint"
trips.sample.1$Checkpoint <- names(busiest)[trips.sample.1$Checkpoint]
# Take into account those that don't have a checkpoint
condition2 <- is.na(trips.sample.1$Checkpoint)
condition3 <- trips.sample.1$Result + trips.sample.1$Secondary - optimal.d.trips - sub.optimal.d.trips > 0
no.transfer <- ifelse(condition2 & condition3, (trips.sample.1$Result + trips.sample.1$Secondary -
optimal.d.trips - sub.optimal.d.trips)[condition2 & condition3], 0)
trips.sample.1$Optim.d.trips <- optimal.d.trips
trips.sample.1$Sub.optim.d.trips <- sub.optimal.d.trips
trips.sample.1$No.transfer <- no.transfer
# Transfers
trips.left.trans <- data.frame(trips.sample.1[!condition2 & condition3, c("From", "Checkpoint", "To", "Time")],
Transferred = (trips.sample.1$Result + trips.sample.1$Secondary - optimal.d.trips - sub.optimal.d.trips)[!condition2 & condition3])
colnames(trips.left.trans)[1] <- "From"
colnames(trips.left.trans)[3] <- "Destination"
colnames(trips.left.trans)[2] <- "To"
# print(i)
# if(i >9) print(head(trips.left.trans))
trips.sample.1 <- merge(trips.sample.1, as.data.frame(as.table(tapply(trips.left.trans$Transferred, trips.left.trans[,c("From","To")], sum))),
all.x = TRUE)
colnames(trips.sample.1)[colnames(trips.sample.1) %in% "Freq"] <- "Transferred"
trips.sample.1$Transferred[is.na(trips.sample.1$Transferred)] <- 0
# Now divide passengers to cars
n.full.8.cars <- (trips.sample.1$Optim.d.trip + trips.sample.1$Sub.optim.d.trip + trips.sample.1$No.transfer + trips.sample.1$Transferred) %/% 8
n.full.4.cars <- (trips.sample.1$Optim.d.trip + trips.sample.1$Sub.optim.d.trip + trips.sample.1$No.transfer + trips.sample.1$Transferred -
n.full.8.cars * 8) %/% 4
n.4.cars.3.pas <- (trips.sample.1$Optim.d.trip + trips.sample.1$Sub.optim.d.trip + trips.sample.1$No.transfer + trips.sample.1$Transferred -
n.full.8.cars * 8 - n.full.4.cars * 4) %/% 3
n.4.cars.2.pas <- (trips.sample.1$Optim.d.trip + trips.sample.1$Sub.optim.d.trip + trips.sample.1$No.transfer + trips.sample.1$Transferred -
n.full.8.cars * 8 - n.full.4.cars * 4 - n.4.cars.3.pas * 3) %/% 2
n.4.cars.1.pas <- trips.sample.1$Optim.d.trip + trips.sample.1$Sub.optim.d.trip + trips.sample.1$No.transfer + trips.sample.1$Transferred -
n.full.8.cars * 8 - n.full.4.cars * 4 - n.4.cars.3.pas * 3 - n.4.cars.2.pas * 2
d8 <- ifelse(trips.sample.1$Optim.d.trip + trips.sample.1$Sub.optim.d.trip + trips.sample.1$No.transfer < 8 * n.full.8.cars,
trips.sample.1$Optim.d.trip + trips.sample.1$Sub.optim.d.trip + trips.sample.1$No.transfer, 8 * n.full.8.cars)
d4 <- ifelse(trips.sample.1$Optim.d.trip + trips.sample.1$Sub.optim.d.trip + trips.sample.1$No.transfer - d8 < 4 * n.full.4.cars,
trips.sample.1$Optim.d.trip + trips.sample.1$Sub.optim.d.trip + trips.sample.1$No.transfer - d8, 4 * n.full.4.cars)
c8 <- 8 * n.full.8.cars - d8
c4 <- 4 * n.full.4.cars - d4
c3 <- n.4.cars.3.pas * (trips.sample.1$Transferred - c8 - c4) # Note: there will be only 1 partially filled car
d3 <- 3 * n.4.cars.3.pas - c3
c2 <- n.4.cars.2.pas * (trips.sample.1$Transferred - c8 - c4)
d2 <- 2 * n.4.cars.2.pas - c2
c1 <- n.4.cars.1.pas * (trips.sample.1$Transferred - c8 - c4)
d1 <- n.4.cars.1.pas - c1
# Add transferred passengers to the next time slot as secondary passengers
trips.left.trans <- merge(trips.left.trans, tdelay)
colnames(trips.left.trans)[5] <- "Secondary"
colnames(trips.left.trans)[1] <- "Origin"
colnames(trips.left.trans)[2] <- "From"
colnames(trips.left.trans)[3] <- "To"
trips.left.trans$Time <- as.character(as.numeric(as.character(trips.left.trans$Time)) + trips.left.trans$Delay)
trips.left.trans <- as.data.frame(as.table(tapply(trips.left.trans$Secondary, trips.left.trans[,
c("From","To","Time")], sum)))
colnames(trips.left.trans)[4] <- "Secondary"
trips.left.trans <- trips.left.trans[!is.na(trips.left.trans$Secondary),]
trips.secondary <- rbind(trips.secondary, trips.left.trans)
trips.out <- rbind(trips.out, data.frame(trips.sample.1[, c("From", "To", "Time")], d8, d4, d3, d2, d1, c8, c4, c3, c2, c1))
}
# Summary
d8 <- tapply(trips.out$d8, trips.out$Time, sum)
d4 <- tapply(trips.out$d4, trips.out$Time, sum)
d3 <- tapply(trips.out$d3, trips.out$Time, sum)
d2 <- tapply(trips.out$d2, trips.out$Time, sum)
d1 <- tapply(trips.out$d1, trips.out$Time, sum)
c8 <- tapply(trips.out$c8, trips.out$Time, sum)
c4 <- tapply(trips.out$c4, trips.out$Time, sum)
c3 <- tapply(trips.out$c3, trips.out$Time, sum)
c2 <- tapply(trips.out$c2, trips.out$Time, sum)
c1 <- tapply(trips.out$c1, trips.out$Time, sum)
test <- data.frame(Time = names(d8), Type = rep(colnames(trips.out)[4:13], each = length(d8)),
Result = c(d8, d4, d3, d2, d1, c8, c4, c3, c2, c1))
test$Time <- as.numeric(as.character(test$Time))
ggplot(test, aes(x = Time, y = Result, fill = Type)) + geom_area()
# Vehicle kilometers
d8 <- tapply(trips.out$d8, trips.out[,c("From","To")], sum)
d4 <- tapply(trips.out$d4, trips.out[,c("From","To")], sum)
d3 <- tapply(trips.out$d3, trips.out[,c("From","To")], sum)
d2 <- tapply(trips.out$d2, trips.out[,c("From","To")], sum)
d1 <- tapply(trips.out$d1, trips.out[,c("From","To")], sum)
c8 <- tapply(trips.out$c8, trips.out[,c("From","To")], sum)
c4 <- tapply(trips.out$c4, trips.out[,c("From","To")], sum)
c3 <- tapply(trips.out$c3, trips.out[,c("From","To")], sum)
c2 <- tapply(trips.out$c2, trips.out[,c("From","To")], sum)
c1 <- tapply(trips.out$c1, trips.out[,c("From","To")], sum)
n.full.8.cars <- as.data.frame(as.table((d8 + c8) / 8))
n.full.4.cars <- as.data.frame(as.table((d4 + c4) / 4))
n.4.cars.3.pas <- as.data.frame(as.table((d3 + c3) / 3))
n.4.cars.2.pas <- as.data.frame(as.table((d2 + c2) / 2))
n.4.cars.1.pas <- as.data.frame(as.table(d1 + c1))
colnames(n.full.8.cars)[3] <- "n.full.8.cars"
colnames(n.full.4.cars)[3] <- "n.full.4.cars"
colnames(n.4.cars.3.pas)[3] <- "n.4.cars.3.pas"
colnames(n.4.cars.2.pas)[3] <- "n.4.cars.2.pas"
colnames(n.4.cars.1.pas)[3] <- "n.4.cars.1.pas"
test2 <- merge(n.full.8.cars, n.full.4.cars)
test2 <- merge(test2, n.4.cars.3.pas)
test2 <- merge(test2, n.4.cars.2.pas)
test2 <- merge(test2, n.4.cars.1.pas)
test2 <- melt(test2, id.vars = c("From","To"), variable_name = "Type")
test2 <- merge(test2, distance[,c("From", "To", "Result")])
colnames(test2)[5] <- "Distance"
test2$vehicle.kilometers <- test2$value * test2$Distance
# Costs
oprint(head(test2))
test3 <- as.data.frame(as.table(tapply(test2$vehicle.kilometers, test2$Type, sum)))
colnames(test3) <- c("Type", "vehicle.kilometers")
fuel.cons.emis.CO2 <- data.frame(Type = c("n.full.8.cars", "n.full.4.cars", "n.4.cars.3.pas", "n.4.cars.2.pas", "n.4.cars.1.pas"),
Consumption = c(8.7, 5.7, 5.7, 5.7, 5.7) / 100, Emis.factor.CO2 = c(232, 153, 153, 153, 153))
emis.factor.PM <- 0.1 # gkm^-1
fuel.price <- 1.374 # 4.8. average cost of a liter of diesel fuel in Finland
maint.price <- 0.0582 # €km^-1
driver.salary <- 2313 / 160 * 1.35 # €h^-1
test3 <- merge(test3, fuel.cons.emis.CO2)
test3$Fuel.cost <- test3$vehicle.kilometers * test3$Consumption * fuel.price
test3$Maint.cost <- test3$vehicle.kilometers * maint.price
test3$Emis.CO2 <- test3$vehicle.kilometers * test3$Emis.factor.CO2
test3$Emis.PM <- test3$vehicle.kilometers * emis.factor.PM
test3$Driver.cost <- test3$vehicle.kilometers / speed * driver.salary
PM.lethality <- c(-7.223e-004, 5.640e-006, 4.228e-005, 5.987e-005, 8.013e-005, 1.150e-004, 2.037e-004, 2.939e-004, 3.598e-004, 4.132e-004,
4.640e-004, 5.139e-004, 5.662e-004, 6.233e-004, 6.854e-004, 7.577e-004, 8.441e-004, 9.519e-004, 1.093e-003, 1.314e-003, 2.805e-003)
PM.lethality <- median(PM.lethality) # deaths / kg
value.of.life <- (2e6 + 0.98e6) / 2 # runif(n, 0.98e6, 2e6) # € / death
emis.price.PM <- PM.lethality * value.of.life # 201.879 # €kg^-1
emis.price.CO2 <- 10e-3 # runif(n, 5e-3, 40e-3) # €kg^-1, the price of CO2 allowances fluctuates quite a bit, in 8.8.2011 they
# were going for 10.74 € per ton of CO2 equivalent
test3$CO2.cost <- test3$Emis.CO2 * emis.price.CO2
test3$PM.cost <- test3$Emis.PM * emis.price.PM
oprint(test3)
| |