" Second By Second Emissions Generator -- This code simply takes information from a MySQL Databse through query, uploads the data into the R environment, uses MOVES created by the EPA to generate emissions rates for vehicles based on speed and acceleration second by second for all specified vehicles." library(RMySQL) rm(list = ls()) #Choose accordingly outpath <- "C:/ProgramData/MySQL/MySQL Server 8.0/Data/" user <<- 'xxx' password <<- 'yyy' killDbConnections <- function () { all_cons <- dbListConnections(MySQL()) print(all_cons) for(con in all_cons) + dbDisconnect(con) print(paste(length(all_cons), " connections killed.")) } GetYearIndex <- function(modelYear) { return (2018 - modelYear +1) } operatingmodes_min = 1 operatingmodes_max = 41 fueltypes = 2 # 1 = petrol, 2=Diesel #source, Fuel, modelyear (back to 1988), opmode emissionsmatrix <<- array(data=-100, dim=c(62,2,31,41)) CreateEmissionsHashTable <- function() { killDbConnections() SQL = "SELECT sourceTypeID, fuelTypeID, modelYearID, OperatingModeID, emissionrateGramperVehHour FROM operatingmoderates WHERE pollutantID = 100" mydb <- dbConnect(MySQL(), user=user, password=password, dbname='yourdatabase', host='yourhost') res <- dbSendQuery(mydb, SQL) emissions <<- fetch(res, n=-1) #now we have the data in a dataframe - put it in a hash table for (i in 1:NROW(emissions)) { src_i = emissions[i,1] f_i = emissions[i,2] #1 is petrol 2 is diesel yr_i = GetYearIndex(emissions[i,3]) op_i = emissions[i,4] + 1 #Opmode bins start at 0 not 1! R indices start at 1 data = emissions[i,5] #rows = NROW(emissionsmatrix) #print(rows) if (src_i <= NROW(emissionsmatrix) && f_i <= NROW(emissionsmatrix) && op_i <= NROW(emissionsmatrix) && yr_i <= NROW(emissionsmatrix)) { emissionsmatrix[src_i, f_i, yr_i, op_i] <<- data } else { print (paste("warning in CreateEmissions ", src_i)) print (f_i) print (yr_i) print (op_i) } } } # VSP=(Av+Bv^2+Cv^3+mva)/m Equation 7 # Where, # VSP vehicle specific power for a specific MOVES vehicle type, # v = instantaneous speed (derived from the simulation data) # a = instantaneous acceleration (derived from the simulation data) # A = rolling resistance term (defined by MOVES for a specific vehicle type) # B= rotating resistance term (defined by MOVES for a specific vehicle type) # C = aerodynamic drag term (defined by MOVES for a specific vehicle type) # M = mass of the vehicle (defined by MOVES) CalculateVSP = function(speed, acceleration, sourceid) { A = physicsmatrix[sourceid, 1];#lookup from parameter set B = physicsmatrix[sourceid, 2]; C = physicsmatrix[sourceid, 3]; M = physicsmatrix[sourceid, 4]; Fs = physicsmatrix[sourceid, 5]; vsp <<- 0.0; vsp <<- (A*speed+B*(speed^2)+C*(speed^3)+M*speed*acceleration) vsp <<- vsp/M; #now calculate STP WHICH IS OK BECAUSE LIGHT DUTY VEHICLES CANCEL OUT stp<<- vsp*M stp<<-stp/Fs print(paste("VSP=", vsp, sep="")) print(speed) print(acceleration) print (A) print (B) print (C) print (M) print (paste("VSP", vsp, sep=" ")) return (stp) } CalculateOpModeBin <- function (vsp, speed, acceleration) { Opmode = 0; if (speed < 1) { Opmode = 1; #Idling } else if (speed < 25) { if (vsp < 0) { Opmode = 11; } else if(vsp > 0 && vsp <= 3) { Opmode = 12; } else if(vsp > 3 && vsp <= 6) { Opmode = 13; } else if (vsp > 6 && vsp <=9) { Opmode = 14; } else if (vsp > 12) { Opmode = 15; } else { Opmode = 16; } } else if (speed <= 50) { if (vsp < 0) { Opmode = 21; } else if(vsp > 0 && vsp <= 3) { Opmode = 22; } else if(vsp > 3 && vsp <= 6) { Opmode = 23; } else if (vsp > 6 && vsp <= 9) { Opmode = 24; } else if (vsp > 9 && vsp <= 12) { Opmode = 25; } else if (vsp > 12 && vsp <= 18) { Opmode = 27; } else if (vsp > 18 && vsp <= 24) { Opmode = 28; } else if (vsp > 24 && vsp <= 30) { Opmode = 29; } else if (vspindex > 30) { Opmode = 30; } } else if (speed >= 50) { if (vspindex < 6) { Opmode = 33; } else if(vsp > 6 && vsp <= 12) { Opmode = 35; } else if(vsp > 12 && vsp <= 18) { Opmode = 37; } else if (vsp > 18 && vsp <= 24) { Opmode = 38; } else if (vsp > 24 && vsp <= 30) { Opmode = 39; } else if (vsp > 30 ) { Opmode = 40; } } return (Opmode) } #todo load up paranmeters for VSP by vehicle type physicsmatrix <<- array(data=0, dim=c(62,5)) CreateVSPParameterTable <- function() { killDbConnections() mydb <- dbConnect(MySQL(), user=user, password=password, dbname='yourdatabase', host='yourhost') sql <- paste("select * from sourceusetypephysics WHERE endModelYearID = 2050 order by sourceTypeID", sep="") res <- dbSendQuery(mydb, sql) tbl <<- fetch(res, n=-1) for (i in 1:NROW(tbl)) { srci = tbl$sourceTypeID[i] physicsmatrix[srci, 1] <<- tbl$rollingTermA[i] physicsmatrix[srci, 2] <<- tbl$rotatingTermB[i] physicsmatrix[srci, 3] <<- tbl$dragTermC[i] physicsmatrix[srci, 4] <<- tbl$sourceMass[i] physicsmatrix[srci, 5] <<- tbl$fixedMassFactor[i] } } #Selects specific car trips from the database GetDistinctTrips <- function() { killDbConnections() mydb <- dbConnect(MySQL(), user=user, password=password, dbname='yourdatabase', host='yourhost') sql <- paste("SELECT distinct(VEHICLEID) FROM yourdatabase.maindata where ID=1;", sep="") res <- dbSendQuery(mydb, sql) alltrips <<- fetch(res, n=-1) } #to do - get each trip GetTripPoints <- function(index) { #query the db for each trip. Store the points in an array mydb <- dbConnect(MySQL(), user=user, password=password, dbname='yourdatabase', host='yourhost') sql <- paste("SELECT * FROM yourdatabase.maindata where ID=1 and VEHICLEID='", alltrips[index,1], "' order by TIME_S", sep="") print(sql) res <- dbSendQuery(mydb, sql) trip <<- fetch(res, n=-1) } #todo - query each trip - find vehicle type, calculate emissions at each second CalculateEmissionsByTrip <- function(sourceID, fueltype, modelyear) { acceleration = 0; rows = NROW(trip)-1 #create some temporary storage for emissions results <<- data.frame( c(1:rows),#ID vector("character", length = rows),#TRIPID c(1:rows),#Time vector("numeric", length = rows),#SUT, c(1:rows),#YEAR c(1:rows),#FUEL vector("numeric", length = rows),#SPEED_MS vector("numeric", length = rows),#SPEED_MPH vector("numeric", length = rows),#Acceleration vector("numeric", length = rows),#VSP, vector("numeric", length = rows),#OpMode, vector("numeric", length = rows),#Pollutant vector("numeric", length = rows),#XPos vector("numeric", length = rows),#YPOS vector("numeric", length = rows),#Pollutant SUMO stringsAsFactors = FALSE )#end create dataframe names=c("SIMID", "TRIPID", "TIME", "SUT", "MODELYR", "FUEL", "SPEED_MS", "SPEED_MPH", "ACCELERATION", "VSP", "OPMODE", "POLLUTANT", "XPOS", "YPOS") names(results) <<- names; for (i in 1:rows) { speed <- as.numeric(trip[i, 7])*0.44704 #MPH to M/S speed_next <-as.numeric(trip[i+1, 7])*0.44704 # MPH to M/S acceleration <- speed_next-speed; #print(paste(i, speed, speed_next, acceleration, sep=", ")) VSP <- CalculateVSP(speed, acceleration, sourceID) OpModeBin <- CalculateOpModeBin(VSP, speed, acceleration) # Get the raw Opmode bin - 0-40 emissions <- CalculateSecBySecEmissions(OpModeBin, sourceID, fueltype, modelyear) results[i,1] <<- trip[i,1] results[i,2] <<- trip[i,3] results[i,3] <<- trip[i,2] results[i,4] <<- sourceID results[i,5] <<- modelyear results[i,6] <<- fueltype results[i,7] <<- speed results[i,8] <<- speed*(1/0.44704) #M/S to MPH results[i,9] <<- acceleration results[i,10] <<- VSP results[i,11] <<- OpModeBin results[i,12] <<- emissions results[i,13] <<- trip[i,5] #XPos results[i,14] <<- trip[i,6] #YPOS results[i,15] <<- trip[i,9] #SUMO pollutant if (emissions < 0) { print(sourceID) print(modelyear) print(OpModeBin) } } #Store results table filepath <- paste(outpath, "temp1.csv", sep="") write.csv(results, file = filepath, row.names=FALSE, quote = FALSE) #finally reload into db table mydb <- dbConnect(MySQL(), user=user, password=password, dbname='yourdatabase', host='yourhost') sql <- paste("LOAD DATA LOCAL INFILE '", filepath , "' INTO TABLE emissions FIELDS TERMINATED BY ',' LINES TERMINATED BY '\r\n' IGNORE 1 LINES;", sep="") print(sql) res <- dbSendQuery(mydb, sql) dbDisconnect(mydb) } CalculateSecBySecEmissions <- function(OpModeBin, sourceID, fueltype, modelyear) { em = 0; print(paste("OpmodeBin ", OpModeBin)) #now look up emissions if (OpModeBin < 0) { print(OpModeBin) } modelyearindex = GetYearIndex(modelyear) em <- emissionsmatrix[sourceID, fueltype, modelyearindex, OpModeBin+1] #OPmode zero indexed, R is 1 indexed # Turning back to second by second em<- em/(60*60) return (em) } GetAgeDistribution <- function() { #Sets up the age distribution matrix from the db killDbConnections() SQL = "SELECT * FROM sourcetypeagedistribution" mydb <- dbConnect(MySQL(), user=user, password=password, dbname='yourdatabase', host='yourhost') res <- dbSendQuery(mydb, SQL) sourcetypedist <<- fetch(res, n=-1) } GetVMTMix <- function() { killDbConnections() SQL = "SELECT * FROM yourdatabase.vmtmix where Daytype = 'Weekday' and period = 'AM' and VMX_RDcode = 2;" mydb <- dbConnect(MySQL(), user=user, password=password, dbname='yourdatabase', host='yourhost') res <- dbSendQuery(mydb, SQL) VMTmix <<- fetch(res, n=-1) } GetRandomVehicle <- function() { sut = 0; year = 0; fuel = 0; #first data vehicle type and fuel #declare a array of indexes rows = NROW(VMTmix) indices = seq(1:rows) index = sample(x = indices, 1, replace = T, prob = VMTmix$VMTmix) sut = VMTmix$MOVES_STcode[index] sut_mod<-sut; #Translates VMTmixSUTTypes to 5 Model types which are : 11 21 31 61 62 if (sut > 31 & sut < 61) { sut_mod<-61 } fuel = VMTmix$MOVES_FTcode[index] #now draw age ages <<- sourcetypedist[sourcetypedist$sourceTypeID==sut,] year = sample(x = ages$ageID, 1, replace = T, prob = ages$ageFraction) year = 2018-year newlist <- list("source" = sut_mod, "year" = year, "fuel" = fuel) return (newlist) } MainFunction <- function() { GetAgeDistribution() GetVMTMix() #First create Hash table CreateEmissionsHashTable() #Now create Parameters for VSP and STP CreateVSPParameterTable() #GetDistinct trips GetDistinctTrips() #now for each trip, calculate emissions and put them back in the db for (i in 1:NROW(alltrips)) { veh <<- GetRandomVehicle() print (veh) GetTripPoints(i) CalculateEmissionsByTrip(veh$source, veh$fuel, veh$year) killDbConnections() } } MainFunction()