.Rhistory

rename_columns <- function(df, suffix) {
original_names <- names(df)
# Keep the first two names unchanged, rename the rest
new_names <- c(original_names[1:2], paste(original_names[-c(1,2)], suffix, sep = "."))
names(df) <- new_names
return(df)
}
# These information are saved in the RDS Files folder including the scenarios
folder_path <- "/Users/amirgazar/Documents/GitHub/EPA_Debarbonization/R Files/RDS Files/"
rds_files <- list.files(path = folder_path, pattern = "\\.rds$", full.names = TRUE)
for (file_path in rds_files) {
variable_name <- tools::file_path_sans_ext(basename(file_path))
dataset <- readRDS(file_path)
assign(variable_name, dataset)
}
# Load Nuclear, hydro and biomass (MW)
# Solar, onshore and offshore wind Name plate capacity based on NE ISO
path <- "/Users/amirgazar/Documents/GitHub/EPA_Debarbonization/Capacity"
setwd(path)
files <- list.files(path = path, pattern = "\\.csv$", full.names = TRUE)
datasets <- list()
for (file in files) {
dataset_name <- gsub(".*/|\\.csv$", "", file)
datasets[[dataset_name]] <- read.csv(file, check.names = FALSE)
}
list2env(datasets, envir = .GlobalEnv) # Releasing the list
dates <- seq(as.Date("2024-01-01"), as.Date("2050-12-31"), by="day")
hours <- 1:24
date_hour_combos <- CJ(Date=dates, Hour=hours)
date_hour_combos[, DayLabel := as.integer(format(Date, "%j"))]
date_hour_combos[, Year := year(Date)]
# If we miss a value for the leap year we use Feb 28 value
date_hour_combos[, DayLabel := ifelse(Year %% 4 == 0 & DayLabel > 59, DayLabel - 1, DayLabel), by = Year]
# Hourly Solar and wind generation P50 and ISONE location for 2024-2050
# Apply the functions to each data set
Solar_CF <- rename_columns(solar_CF, "solar")
Onwind_CF <- rename_columns(onwind_CF, "onwind")
Offwind_CF <- rename_columns(offwind_CF, "offwind")
setDT(Solar_CF)
setDT(Onwind_CF)
setDT(Offwind_CF)
setDT(Solar_NPC)
setDT(Onwind_NPC)
setDT(Offwind_NPC)
# Solar
solar_cf_long <- melt(Solar_CF, id.vars = c("DayLabel", "Hour"),
variable.name = "Scenario", value.name = "CF_Value")
solar_gen <- merge(date_hour_combos, Solar_NPC, by="Year")
solar_gen <- solar_gen[Solar_CF, on = .(DayLabel, Hour), .(Date, Hour, Solar = ISONE * i.P50.solar)]
# Off wind
offwind_cf_long <- melt(Offwind_CF, id.vars = c("DayLabel", "Hour"),
variable.name = "Scenario", value.name = "CF_Value")
offwind_gen <- merge(date_hour_combos, Offwind_NPC, by="Year")
offwind_gen <- offwind_gen[Offwind_CF, on = .(DayLabel, Hour), .(Date, Hour, Offwind = ISONE * i.P50.offwind)]
# on wind
onwind_cf_long <- melt(Onwind_CF, id.vars = c("DayLabel", "Hour"),
variable.name = "Scenario", value.name = "CF_Value")
onwind_gen <- merge(date_hour_combos, Onwind_NPC, by="Year")
onwind_gen <- onwind_gen[Onwind_CF, on = .(DayLabel, Hour), .(Date, Hour, Onwind = ISONE * i.P50.onwind)]
# Hourly Import/export and storage capacity
# Import/export
setDT(Import_export_capacity)
Import_export_capacity[, Date := as.Date(paste(Year, Month, Day, sep="-"), format="%Y-%B-%d")]
Import_export_capacity[, DayLabel := yday(Date)]
Import_export_capacity[, DayLabel := ifelse(Year %% 4 == 0 & DayLabel > 59, DayLabel - 1, DayLabel), by = Year]
Import_export_capacity <- Import_export_capacity[, .(Date, DayLabel, NYISO_Imports, NYISO_Exports, NB_Imports, NB_Exports, HQ_Imports, HQ_Exports)]
dt <- data.table(Date = rep(dates, each = 24))
dt[, Hour := rep(1:24, times = length(dates))]
dt[, DayLabel := yday(Date)]
dt[, Year := year(Date)]
dt[, DayLabel := ifelse((Year %% 4 == 0 & Year %% 100 != 0 | Year %% 400 == 0) & DayLabel > 59, DayLabel - 1, DayLabel), by = Year]
merged_data <- merge(dt, Import_export_capacity, by = "DayLabel", all = TRUE)
Import_export_capacity <- merged_data[, .(Date = Date.x, Hour, NYISO_Imports, NYISO_Exports, NB_Imports, NB_Exports, HQ_Imports, HQ_Exports)]
base_year <- 2023
Import_export_capacity[, YearsSinceBase := as.integer(format(Date, "%Y")) - base_year]
# Storage
setDT(Storage_capacity)
merged_data <- merge(dt, Storage_capacity, by = "Year", all = TRUE)
Storage_capacity <- merged_data[, .(Date, Hour, Storage = ISONE)]
# Nuclear, Bio and Hydro
setDT(Nuclear_hydro_bio_NPC)
merged_data <- merge(dt, Nuclear_hydro_bio_NPC, by = "Year", all = TRUE)
Nuclear_hydro_bio <- merged_data[, .(Date, Hour, Nuclear, Hydro, Bio)]
# Scenario 2, no transmission increases
all_data <- Reduce(function(x, y) merge(x, y, by=c("Date", "Hour"), all=TRUE),
list(Nuclear_hydro_bio, solar_gen,
onwind_gen, offwind_gen, Import_export_capacity, Storage_capacity))
all_data[is.na(all_data)] <- 0
scenario_2 <- all_data
scenario_2$Date <- as.Date(scenario_2$Date)
setDT(scenario_2)
saveRDS(scenario_2, "/Users/amirgazar/Documents/GitHub/EPA_Debarbonization/R Files/RDS Files Scenarios/scenario_2.rds")
# increase the Import_export_capacity based on all options pathway
cols_to_adjust <- c("HQ_Imports", "HQ_Exports") # We only increase HQ Capacity, if we want to increase NY OR NB "NYISO_Imports", "NYISO_Exports", "NB_Imports", "NB_Exports"
Import_export_capacity[, YearsSinceBase := as.numeric(YearsSinceBase)]
Import_export_capacity[, (cols_to_adjust) := lapply(.SD, function(x) as.numeric(x)), .SDcols = cols_to_adjust]
Import_export_capacity[, (cols_to_adjust) := lapply(.SD, function(x) x * (1 + YearsSinceBase/15)), .SDcols = cols_to_adjust]
cols_to_adjust <- c("NYISO_Imports", "NYISO_Exports", "NB_Imports", "NB_Exports") # We only increase HQ Capacity, if we want to increase NY OR NB
Import_export_capacity[, (cols_to_adjust) := lapply(.SD, function(x) as.numeric(x)), .SDcols = cols_to_adjust]
Import_export_capacity[, (cols_to_adjust) := lapply(.SD, function(x) x * (1 + YearsSinceBase/50)), .SDcols = cols_to_adjust]
Import_export_capacity <- Import_export_capacity[, .(Date, Hour, NYISO_Imports, NYISO_Exports, NB_Imports, NB_Exports, HQ_Imports, HQ_Exports)]
# Scenario 1, everything is average i.e P50
all_data <- Reduce(function(x, y) merge(x, y, by=c("Date", "Hour"), all=TRUE),
list(Nuclear_hydro_bio, solar_gen,
onwind_gen, offwind_gen, Import_export_capacity, Storage_capacity))
all_data[is.na(all_data)] <- 0
scenario_1 <- all_data
scenario_1$Date <- as.Date(scenario_1$Date)
setDT(scenario_1)
saveRDS(scenario_1, "/Users/amirgazar/Documents/GitHub/EPA_Debarbonization/R Files/RDS Files Scenarios/scenario_1.rds")
# Scenario 3, high transmission increases (2x transmission, limited offshore)
# 10 GW Northeast Cap
# Off wind
offwind_cf_long <- melt(Offwind_CF, id.vars = c("DayLabel", "Hour"),
variable.name = "Scenario", value.name = "CF_Value")
offwind_gen <- merge(date_hour_combos, Offwind_NPC, by="Year")
offwind_gen <- offwind_gen[Offwind_CF, on = .(DayLabel, Hour), .(Date, Hour, Offwind = ISONE * i.P50.offwind/2)]
base_year <- 2023
Import_export_capacity[, YearsSinceBase := as.integer(format(Date, "%Y")) - base_year]
cols_to_adjust <- c("HQ_Imports", "HQ_Exports") # We only increase HQ Capacity, if we want to increase NY OR NB "NYISO_Imports", "NYISO_Exports", "NB_Imports", "NB_Exports"
Import_export_capacity[, YearsSinceBase := as.numeric(YearsSinceBase)]
Import_export_capacity[, (cols_to_adjust) := lapply(.SD, function(x) as.numeric(x)), .SDcols = cols_to_adjust]
Import_export_capacity[, (cols_to_adjust) := lapply(.SD, function(x) x * (1 + YearsSinceBase/20)), .SDcols = cols_to_adjust]
all_data <- Reduce(function(x, y) merge(x, y, by=c("Date", "Hour"), all=TRUE),
list(Nuclear_hydro_bio, solar_gen,
onwind_gen, offwind_gen, Import_export_capacity, Storage_capacity))
all_data[is.na(all_data)] <- 0
scenario_3 <- all_data
scenario_3$Date <- as.Date(scenario_3$Date)
setDT(scenario_3)
saveRDS(scenario_3, "/Users/amirgazar/Documents/GitHub/EPA_Debarbonization/R Files/RDS Files Scenarios/scenario_3.rds")
plot(scenario_3$HQ_Imports)
# Load libraries
library(httr)
library(htmltools)
library(jsonlite)
library(data.table)
library(stringr)
library(dplyr)
library(readxl)
library(ggplot2)
library(lubridate)
library(zoo)
# Defining the scenario values, NOTE: Summer:June-Sept and Winter:Oct-May
# Function to rename columns, skipping the first two
rename_columns <- function(df, suffix) {
original_names <- names(df)
# Keep the first two names unchanged, rename the rest
new_names <- c(original_names[1:2], paste(original_names[-c(1,2)], suffix, sep = "."))
names(df) <- new_names
return(df)
}
# These information are saved in the RDS Files folder including the scenarios
folder_path <- "/Users/amirgazar/Documents/GitHub/EPA_Debarbonization/R Files/RDS Files/"
rds_files <- list.files(path = folder_path, pattern = "\\.rds$", full.names = TRUE)
for (file_path in rds_files) {
variable_name <- tools::file_path_sans_ext(basename(file_path))
dataset <- readRDS(file_path)
assign(variable_name, dataset)
}
# Load Nuclear, hydro and biomass (MW)
# Solar, onshore and offshore wind Name plate capacity based on NE ISO
path <- "/Users/amirgazar/Documents/GitHub/EPA_Debarbonization/Capacity"
setwd(path)
files <- list.files(path = path, pattern = "\\.csv$", full.names = TRUE)
datasets <- list()
for (file in files) {
dataset_name <- gsub(".*/|\\.csv$", "", file)
datasets[[dataset_name]] <- read.csv(file, check.names = FALSE)
}
list2env(datasets, envir = .GlobalEnv) # Releasing the list
dates <- seq(as.Date("2024-01-01"), as.Date("2050-12-31"), by="day")
hours <- 1:24
date_hour_combos <- CJ(Date=dates, Hour=hours)
date_hour_combos[, DayLabel := as.integer(format(Date, "%j"))]
date_hour_combos[, Year := year(Date)]
# If we miss a value for the leap year we use Feb 28 value
date_hour_combos[, DayLabel := ifelse(Year %% 4 == 0 & DayLabel > 59, DayLabel - 1, DayLabel), by = Year]
# Hourly Solar and wind generation P50 and ISONE location for 2024-2050
# Apply the functions to each data set
Solar_CF <- rename_columns(solar_CF, "solar")
Onwind_CF <- rename_columns(onwind_CF, "onwind")
Offwind_CF <- rename_columns(offwind_CF, "offwind")
setDT(Solar_CF)
setDT(Onwind_CF)
setDT(Offwind_CF)
setDT(Solar_NPC)
setDT(Onwind_NPC)
setDT(Offwind_NPC)
# Solar
solar_cf_long <- melt(Solar_CF, id.vars = c("DayLabel", "Hour"),
variable.name = "Scenario", value.name = "CF_Value")
solar_gen <- merge(date_hour_combos, Solar_NPC, by="Year")
solar_gen <- solar_gen[Solar_CF, on = .(DayLabel, Hour), .(Date, Hour, Solar = ISONE * i.P50.solar)]
# Off wind
offwind_cf_long <- melt(Offwind_CF, id.vars = c("DayLabel", "Hour"),
variable.name = "Scenario", value.name = "CF_Value")
offwind_gen <- merge(date_hour_combos, Offwind_NPC, by="Year")
offwind_gen <- offwind_gen[Offwind_CF, on = .(DayLabel, Hour), .(Date, Hour, Offwind = ISONE * i.P50.offwind)]
# on wind
onwind_cf_long <- melt(Onwind_CF, id.vars = c("DayLabel", "Hour"),
variable.name = "Scenario", value.name = "CF_Value")
onwind_gen <- merge(date_hour_combos, Onwind_NPC, by="Year")
onwind_gen <- onwind_gen[Onwind_CF, on = .(DayLabel, Hour), .(Date, Hour, Onwind = ISONE * i.P50.onwind)]
# Hourly Import/export and storage capacity
# Import/export
setDT(Import_export_capacity)
Import_export_capacity[, Date := as.Date(paste(Year, Month, Day, sep="-"), format="%Y-%B-%d")]
Import_export_capacity[, DayLabel := yday(Date)]
Import_export_capacity[, DayLabel := ifelse(Year %% 4 == 0 & DayLabel > 59, DayLabel - 1, DayLabel), by = Year]
Import_export_capacity <- Import_export_capacity[, .(Date, DayLabel, NYISO_Imports, NYISO_Exports, NB_Imports, NB_Exports, HQ_Imports, HQ_Exports)]
dt <- data.table(Date = rep(dates, each = 24))
dt[, Hour := rep(1:24, times = length(dates))]
dt[, DayLabel := yday(Date)]
dt[, Year := year(Date)]
dt[, DayLabel := ifelse((Year %% 4 == 0 & Year %% 100 != 0 | Year %% 400 == 0) & DayLabel > 59, DayLabel - 1, DayLabel), by = Year]
merged_data <- merge(dt, Import_export_capacity, by = "DayLabel", all = TRUE)
Import_export_capacity <- merged_data[, .(Date = Date.x, Hour, NYISO_Imports, NYISO_Exports, NB_Imports, NB_Exports, HQ_Imports, HQ_Exports)]
base_year <- 2023
Import_export_capacity[, YearsSinceBase := as.integer(format(Date, "%Y")) - base_year]
# Storage
setDT(Storage_capacity)
merged_data <- merge(dt, Storage_capacity, by = "Year", all = TRUE)
Storage_capacity <- merged_data[, .(Date, Hour, Storage = ISONE)]
# Nuclear, Bio and Hydro
setDT(Nuclear_hydro_bio_NPC)
merged_data <- merge(dt, Nuclear_hydro_bio_NPC, by = "Year", all = TRUE)
Nuclear_hydro_bio <- merged_data[, .(Date, Hour, Nuclear, Hydro, Bio)]
# Scenario 2, no transmission increases
all_data <- Reduce(function(x, y) merge(x, y, by=c("Date", "Hour"), all=TRUE),
list(Nuclear_hydro_bio, solar_gen,
onwind_gen, offwind_gen, Import_export_capacity, Storage_capacity))
all_data[is.na(all_data)] <- 0
scenario_2 <- all_data
scenario_2$Date <- as.Date(scenario_2$Date)
setDT(scenario_2)
saveRDS(scenario_2, "/Users/amirgazar/Documents/GitHub/EPA_Debarbonization/R Files/RDS Files Scenarios/scenario_2.rds")
# increase the Import_export_capacity based on all options pathway
cols_to_adjust <- c("HQ_Imports", "HQ_Exports") # We only increase HQ Capacity, if we want to increase NY OR NB "NYISO_Imports", "NYISO_Exports", "NB_Imports", "NB_Exports"
Import_export_capacity[, YearsSinceBase := as.numeric(YearsSinceBase)]
Import_export_capacity[, (cols_to_adjust) := lapply(.SD, function(x) as.numeric(x)), .SDcols = cols_to_adjust]
Import_export_capacity[, (cols_to_adjust) := lapply(.SD, function(x) x * (1 + YearsSinceBase/15)), .SDcols = cols_to_adjust]
cols_to_adjust <- c("NYISO_Imports", "NYISO_Exports", "NB_Imports", "NB_Exports") # We only increase HQ Capacity, if we want to increase NY OR NB
Import_export_capacity[, (cols_to_adjust) := lapply(.SD, function(x) as.numeric(x)), .SDcols = cols_to_adjust]
Import_export_capacity[, (cols_to_adjust) := lapply(.SD, function(x) x * (1 + YearsSinceBase/50)), .SDcols = cols_to_adjust]
Import_export_capacity <- Import_export_capacity[, .(Date, Hour, NYISO_Imports, NYISO_Exports, NB_Imports, NB_Exports, HQ_Imports, HQ_Exports)]
# Scenario 1, everything is average i.e P50
all_data <- Reduce(function(x, y) merge(x, y, by=c("Date", "Hour"), all=TRUE),
list(Nuclear_hydro_bio, solar_gen,
onwind_gen, offwind_gen, Import_export_capacity, Storage_capacity))
all_data[is.na(all_data)] <- 0
scenario_1 <- all_data
scenario_1$Date <- as.Date(scenario_1$Date)
setDT(scenario_1)
saveRDS(scenario_1, "/Users/amirgazar/Documents/GitHub/EPA_Debarbonization/R Files/RDS Files Scenarios/scenario_1.rds")
# Scenario 3, high transmission increases (2x transmission, limited offshore)
# 10 GW Northeast Cap
# Off wind
offwind_cf_long <- melt(Offwind_CF, id.vars = c("DayLabel", "Hour"),
variable.name = "Scenario", value.name = "CF_Value")
offwind_gen <- merge(date_hour_combos, Offwind_NPC, by="Year")
offwind_gen <- offwind_gen[Offwind_CF, on = .(DayLabel, Hour), .(Date, Hour, Offwind = ISONE * i.P50.offwind/2)]
base_year <- 2023
Import_export_capacity[, YearsSinceBase := as.integer(format(Date, "%Y")) - base_year]
cols_to_adjust <- c("HQ_Imports", "HQ_Exports") # We only increase HQ Capacity, if we want to increase NY OR NB "NYISO_Imports", "NYISO_Exports", "NB_Imports", "NB_Exports"
Import_export_capacity[, YearsSinceBase := as.numeric(YearsSinceBase)]
Import_export_capacity[, (cols_to_adjust) := lapply(.SD, function(x) as.numeric(x)), .SDcols = cols_to_adjust]
Import_export_capacity[, (cols_to_adjust) := lapply(.SD, function(x) x * (1 + YearsSinceBase/25)), .SDcols = cols_to_adjust]
all_data <- Reduce(function(x, y) merge(x, y, by=c("Date", "Hour"), all=TRUE),
list(Nuclear_hydro_bio, solar_gen,
onwind_gen, offwind_gen, Import_export_capacity, Storage_capacity))
all_data[is.na(all_data)] <- 0
scenario_3 <- all_data
scenario_3$Date <- as.Date(scenario_3$Date)
setDT(scenario_3)
saveRDS(scenario_3, "/Users/amirgazar/Documents/GitHub/EPA_Debarbonization/R Files/RDS Files Scenarios/scenario_3.rds")
plot(scenario_3$HQ_Imports)
View(scenario_3)
mean((scenario_3$HQ_Imports))
start_time <- Sys.time()
# Dispatch curve, offsetting demand with generation
# Load libraries
library(httr)
library(htmltools)
library(jsonlite)
library(data.table)
library(stringr)
library(dplyr)
library(readxl)
library(ggplot2)
library(lubridate)
library(zoo)
# Load the predicted hourly demand
# Load wind and solar simulations
# Load oil and gas simulations
# Load imports
# These information are saved in the RDS Files folder including the scenarios
folder_path <- "/Users/amirgazar/Documents/GitHub/EPA_Debarbonization/R Files/RDS Files Scenarios/"
rds_files <- list.files(path = folder_path, pattern = "\\.rds$", full.names = TRUE)
for (file_path in rds_files) {
variable_name <- tools::file_path_sans_ext(basename(file_path))
dataset <- readRDS(file_path)
assign(variable_name, dataset)
}
setDT(demand_data)
# Grid includes Thermal Energy
# Functions:
# Displaces x MW each hour based on each scenario given and returns the emissions and loads removed
eliminate_load_and_emissions <- function(scenario) {
results <- list()
# Setting the dates and hours to map the data
unique_dates <- as.Date(sort(unique(demand_data$Date)), origin = "1970-01-01")
unique_hours <- sort(unique(demand_data$Hour))
storage_status <- 0
storage_hr <- 0
storage_hr_max <- 11
import_sources <- c("HQ_Imports", "NB_Imports", "NYISO_Imports")
for (date in unique_dates) {
for (hour in unique_hours) {
date_selected <- as.Date(date, origin = "1970-01-01")
demand_hour <- demand_data[demand_data$Date == date_selected & demand_data$Hour == hour, "Demand"]
demand_hour <- demand_hour$Demand
gen_hour_clean <-  scenario[Date == date_selected & Hour == hour]
gen_clean <- sum(gen_hour_clean$Nuclear, gen_hour_clean$Hydro, gen_hour_clean$Bio, gen_hour_clean$Solar, gen_hour_clean$Onwind, gen_hour_clean$Offwind)
# Save clean energy gen
solar <- scenario[Date == date_selected & Hour == hour, "Solar"]
onwind <- scenario[Date == date_selected & Hour == hour, "Onwind"]
offwind <- scenario[Date == date_selected & Hour == hour, "Offwind"]
nuclear <- scenario[Date == date_selected & Hour == hour, "Nuclear"]
hydro <- scenario[Date == date_selected & Hour == hour, "Hydro"]
bio <- scenario[Date == date_selected & Hour == hour, "Bio"]
# If clean gen exceeds demand
if (demand_hour < gen_clean) {
demand_hour_fossil_fuels <- 0  # No fossil fuels needed
excess_energy <- gen_clean - demand_hour
storage_status_new <- min(storage_status + excess_energy, gen_hour_clean$Storage)
} else {
demand_hour_fossil_fuels <- max(0, demand_hour - gen_clean - storage_status)
storage_status <- 0
}
# Initializing variables to track used imports
used_imports <- list(HQ_Imports=0, NB_Imports=0, NYISO_Imports=0)
# Sequentially reduce demand using imports, updating tracking variables
if (demand_hour_fossil_fuels > 0) {
for (source in import_sources) {
available_import <- -gen_hour_clean[[source]]
import_needed <- min(demand_hour_fossil_fuels, available_import)
demand_hour_fossil_fuels <- demand_hour_fossil_fuels - import_needed
# Update used imports directly using the source variable
used_imports[[source]] <- used_imports[[source]] + import_needed
}
}
# Fossil fuel loading
fossil_fuels <- Fossil_Fuels_data[Date == date_selected & Hour == hour]
setorder(fossil_fuels, -CF)
# Calculate Gen CO2 and NOx
fossil_fuels[, Hourly_Gen := CF * Nameplate_capacity]
fossil_fuels[, CO2_total := CO2 * Hourly_Gen]
fossil_fuels[, NOx_total := NOx * Hourly_Gen]
sum_gen_hour_fossil_fuels <- sum(fossil_fuels$Hourly_Gen)
gen_hourly <- 0
sorted_loads <- fossil_fuels$Hourly_Gen
sorted_facilities <- fossil_fuels$Facility_Unit.ID
if (sum_gen_hour_fossil_fuels > demand_hour_fossil_fuels) {
displaced_load <- sum_gen_hour_fossil_fuels - demand_hour_fossil_fuels
gen_shortage <- 0
gen_hourly <- demand_hour
fossil <- demand_hour_fossil_fuels
} else {
displaced_load <- 0
gen_shortage <- demand_hour_fossil_fuels - sum_gen_hour_fossil_fuels
gen_hourly <- gen_clean + sum_gen_hour_fossil_fuels
fossil <- sum_gen_hour_fossil_fuels
}
if (gen_shortage <= 0) {
total_consumed <- 0
consumed_loads <- 0
consumed_facilities <- c()
eliminated_facilities <- c()
for (i in 1:length(sorted_loads)) {
if (sorted_loads[i] != 0) {
if (total_consumed + sorted_loads[i] <= demand_hour_fossil_fuels) {
total_consumed <- total_consumed + sorted_loads[i]
consumed_loads <- c(consumed_loads, sorted_loads[i])
consumed_facilities <- c(consumed_facilities, sorted_facilities[i])
} else {
remaining_to_eliminate <- demand_hour_fossil_fuels - total_consumed
total_consumed <- demand_hour_fossil_fuels
consumed_loads <- c(consumed_loads, remaining_to_eliminate)
eliminated_facilities <- c(eliminated_facilities, sorted_facilities[i])
}
}
}
eliminated_facilities_hour <- data.table(Facility = eliminated_facilities)
eliminated_facility_ids <- eliminated_facilities_hour$Facility
loads_for_eliminated_facilities <- fossil_fuels[eliminated_facility_ids, on = .(Facility_Unit.ID), nomatch = 0, .(Facility_Unit.ID, Hourly_Gen)]
sum_load_before_last <- sum(loads_for_eliminated_facilities$Hourly_Gen, na.rm = TRUE)
remaining_load <-  sum_load_before_last - displaced_load
last_facility_load <- loads_for_eliminated_facilities$Hourly_Gen[[1]]
percentage_reduction <- 1- (remaining_load) / last_facility_load
co2_for_eliminated_facilities <- fossil_fuels[eliminated_facility_ids, on = .(Facility_Unit.ID), nomatch = 0, .(Facility_Unit.ID, CO2_total)]
nox_for_eliminated_facilities <- fossil_fuels[eliminated_facility_ids, on = .(Facility_Unit.ID), nomatch = 0, .(Facility_Unit.ID, NOx_total)]
co2_for_eliminated_facilities$CO2[[1]] <- co2_for_eliminated_facilities$CO2[[1]] * (percentage_reduction)
nox_for_eliminated_facilities$NOx[[1]] <- nox_for_eliminated_facilities$NOx[[1]] * (percentage_reduction)
eliminated_co2 <- sum(co2_for_eliminated_facilities$CO2, na.rm = TRUE)
eliminated_nox <- sum(nox_for_eliminated_facilities$NOx, na.rm = TRUE)
loads_for_eliminated_facilities$Hourly_Gen[[1]] <- loads_for_eliminated_facilities$Hourly_Gen[[1]] * (percentage_reduction)
co2_vector <- unlist(co2_for_eliminated_facilities$CO2, use.names = FALSE)
co2_vector[is.na(co2_vector)] <- 0
nox_vector <- unlist(nox_for_eliminated_facilities$NOx, use.names = FALSE)
nox_vector[is.na(nox_vector)] <- 0
load_vector <- unlist(loads_for_eliminated_facilities$Hourly_Gen, use.names = FALSE)
load_vector[is.na(load_vector)] <- 0
} else {
eliminated_facilities_hour <- data.table(Facility = NA)
co2_vector <- NA
nox_vector <- NA
load_vector <- NA
eliminated_co2 <- 0
eliminated_nox <- 0
}
result <- data.table(Date = date_selected, Hour = hour,
Facility = eliminated_facilities_hour$Facility,
CO2 = co2_vector, NOx = nox_vector, Displaced.Load = load_vector,
CO2.total_hr = eliminated_co2, NOx.total_hr = eliminated_nox,
Displaced.Load.total_hr = displaced_load, Gen.total_hr = gen_hourly,
Demand.total_hr = demand_hour, fossil,
HQ_imports = used_imports$HQ_Imports,
NB_imports = used_imports$NB_Imports,
NYISO_imports = used_imports$NYISO_Imports,
storage_status, solar, onwind, offwind, nuclear, hydro, bio)
results <- append(results, list(result))
}
}
return(results)
}
# Annual function
eliminate_load_and_emissions_annually <- function(scenario) {
results <- list()
Fossil_Fuels_data$Year <- as.numeric(format(as.Date(Fossil_Fuels_data$Date), "%Y"))
unique_years <- unique(as.numeric(format(as.Date(demand_data$Date), "%Y")))
for (year in unique_years) {
annual_demand <- sum(demand_data[format(as.Date(demand_data$Date), "%Y") == year, "Demand"])
scenario_yearly <- scenario[format(as.Date(scenario$Date), "%Y") == year,]
gen_yearly_clean <- sum(scenario_yearly$Nuclear + scenario_yearly$Hydro + scenario_yearly$Bio +
scenario_yearly$Solar + scenario_yearly$Onwind + scenario_yearly$Offwind)
imports <- -sum(scenario_yearly$HQ_Imports + scenario_yearly$NB_Imports + scenario_yearly$NYISO_Imports)
demand_fossil_fuels <- max((annual_demand - gen_yearly_clean - imports), 0)
annual_fossil_fuels_gen <- sum(Fossil_Fuels_data[Fossil_Fuels_data$Year == year, "Hourly_Gen"])
annual_CO2_total <- sum(Fossil_Fuels_data[Fossil_Fuels_data$Year == year, "CO2_total"], na.rm = TRUE)
annual_NOx_total <- sum(Fossil_Fuels_data[Fossil_Fuels_data$Year == year, "NOx_total"], na.rm = TRUE)
fossil_fuels_used <- min(demand_fossil_fuels, annual_fossil_fuels_gen)
annual_CO2_total <- annual_CO2_total * fossil_fuels_used/annual_fossil_fuels_gen
annual_NOx_total <- annual_NOx_total * fossil_fuels_used/annual_fossil_fuels_gen
Nuclear <- sum(scenario_yearly$Nuclear)/1000000
Hydro = sum(scenario_yearly$Hydro)/1000000
Bio = sum(scenario_yearly$Bio)/1000000
Solar = sum(scenario_yearly$Solar)/1000000
Onwind = sum(scenario_yearly$Onwind)/1000000
Offwind = sum(scenario_yearly$Offwind)/1000000
HQ_Imports = -sum(scenario_yearly$HQ_Imports)/1000000
NB_Imports = -sum(scenario_yearly$NB_Imports)/1000000
NYISO_Imports = -sum(scenario_yearly$NYISO_Imports)/1000000
result <- data.table(Year = year,
Demand = annual_demand/1000000,
Clean_Generation = gen_yearly_clean/1000000,
Fossil_Fuels_Used = fossil_fuels_used/1000000,
Total_CO2_Emissions = annual_CO2_total/1000000,
Total_NOx_Emissions = annual_NOx_total/1000000,
Nuclear, Hydro, Bio, Solar, Onwind, Offwind, HQ_Imports,
NB_Imports, NYISO_Imports)
results <- append(results, list(result))
}
return(results)
}
# Running Annual Scenario , Refer to Scenarios file to modify scenario specifications.
scenario <- scenario_1
scenario_1_results_annual <- rbindlist(results <- eliminate_load_and_emissions_annually(scenario))
scenario <- scenario_2
scenario_2_results_annual <- rbindlist(results <- eliminate_load_and_emissions_annually(scenario))
scenario <- scenario_3
scenario_3_results_annual <- rbindlist(results <- eliminate_load_and_emissions_annually(scenario))
saveRDS(scenario_1_results_annual, "/Users/amirgazar/Documents/GitHub/EPA_Debarbonization/R Files/RDS Results/scenario_1_results_annual.rds")
saveRDS(scenario_2_results_annual, "/Users/amirgazar/Documents/GitHub/EPA_Debarbonization/R Files/RDS Results/scenario_2_results_annual.rds")
saveRDS(scenario_3_results_annual, "/Users/amirgazar/Documents/GitHub/EPA_Debarbonization/R Files/RDS Results/scenario_3_results_annual.rds")
# Do electrical interties stimulate Canadian hydroelectric development? Using causal inference to scope environmental impact assessment in evolving sociotechnical systems
# Amir M. Gazar1,2,*, Mark E. Borsuk3, Ryan S.D. Calder1,2,3,4,5
# 1 Department of Population Health Sciences, Virginia Tech, Blacksburg, VA, 24061, USA
# 2 Global Change Center, Virginia Tech, Blacksburg, VA, 24061, USA
# 3 Department of Civil and Environmental Engineering, Duke University, Durham, NC, 277015, USA
# 4 Faculty of Health Sciences, Virginia Tech, Roanoke, VA, 24016, USA
# 5 Department of Civil and Environmental Engineering, Virginia Tech, Blacksburg, VA, 24061, USA
# *Contact: amirgazar@vt.edu.
# All rights reserved under Creative Commons 4.0
# Install Rgraphviz from Bioconductor
#if (!requireNamespace("BiocManager", quietly = TRUE))
#  install.packages("BiocManager")
#BiocManager::install("Rgraphviz")
# List of other packages to install
#packages <- c("bnlearn", "gRain", "visNetwork", "ggplot2", "zoo", "scales", "gridExtra", "dplyr", "MASS","svglite", "tidyverse")
# Install packages
#install.packages(packages)
# Load all the required packages
invisible(lapply(c("Rgraphviz", "bnlearn", "gRain", "visNetwork", "ggplot2",
"stats", "zoo", "scales", "gridExtra", "dplyr", "MASS","svglite","tidyverse"), library, character.only = TRUE))
# Set Working Directory (change for your setup)
setwd("/Users/amirgazar/Documents/Hydro EIA Code/data")
# Set Working Directory (change for your setup)
setwd("/Users/amirgazar/Documents/Hydro EIA Code/data")
setwd("~/Documents/GitHub/Hydro EIA Code/data")
# Set Working Directory (change for your setup)
setwd("/Users/amirgazar/Documents/Hydro EIA Code/data")