diff --git a/.gitignore b/.gitignore index f26b794..bd7419f 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,9 @@ # ---> R +#Folders with large data not to be saved +Data/Cleaned_Data/ +Data/Results/ +#Do not save any swap files +*.swp # History files .Rhistory .Rapp.history diff --git a/Analysis.r b/Analysis.r index 26abda4..48f654c 100644 --- a/Analysis.r +++ b/Analysis.r @@ -1,132 +1,27 @@ library(tidyverse) -library(scales) -library(RcppRoll) -library(lpSolve) -RES <- readRDS("Results/Storage_Values_by_Facility_and_Variable_Discounts.Rds") -RES_REDUCED_FEE <- RES -RES_INCREASED_FEE <- RES - -CV1 <- 1.3074*(10607030-1060703) #Data from Texas Report, converted from 2018 to Dec 2025 dollars -CV2 <- 1.2874*(6984013-1117442) #Data from New Mexico Report, Converted from 2019 to Dec 2025 - -LENGTH <- length(RES) -SHIPPING_COST <- 1.2874*26000 #Inflation adjusted from New Mexico Report - -for(i in 1:LENGTH){ - RES[[i]] <- RES[[i]] %>% filter(!(Facility %in% c("Palo Verde","Vogtle"))) #The shipping cost is higher than the marginal value per SNF at these locations in the lat period of study 2083. - RES[[i]]$Revenue <- CV1*as.numeric(RES[[i]]$Revenue) - RES[[i]]$Year <- as.numeric(RES[[i]]$Year) - DISCOUNT <- as.numeric(gsub("Per_","",names(RES)[i]))/100 - RES[[i]] <- RES[[i]] %>% left_join(read_csv("Data/Raw_Data/Curie_Spent_Fuel_Site_Totals.csv")) %>% select(Year,Facility,Total_Tons,Revenue) -} - - -MAX_REV <- function(YEAR,TBL,CIFS_SIZE,COST=0){ - TBL <- TBL %>% filter(Year==YEAR) - REV <- TBL %>% pull(Revenue) - VOL <- TBL %>% pull(Total_Tons) - RES <- lp(direction = "max", objective.in = REV-COST, const.mat = matrix(VOL, nrow = 1),const.dir = "<=", const.rhs = CIFS_SIZE, all.bin = TRUE) - return(sum(RES$objective)) - } -FACILITY_REVENUE <-function(CAPACITY){OUTCOME <- do.call(cbind,lapply(1:length(RES),function(i){sapply(1960:2083,MAX_REV,TBL=RES[[i]],CIFS_SIZE=CAPACITY)})) - OUTCOME <- cbind(1960:2083,OUTCOME ) %>% as_tibble - colnames(OUTCOME) <- c("Year",parse_number(names(RES))/100) - OUTCOME <- OUTCOME %>% pivot_longer(-Year,names_to="Discount",values_to="Revenue") %>% mutate(Capacity=CAPACITY) %>% select(Year,Capacity,everything()) - return(OUTCOME) - } -REV_RES <- do.call(rbind,lapply(c(8680,8680+5000*19,5000,40000,10000),FACILITY_REVENUE)) -REV_RES %>% mutate(as.numeric(Discount)) -%>% mutate( -CIFS <- rbind(readRDS("Data/Cleaned_Data/Texas_CIFS_Costs.Rds"),readRDS("Data/Cleaned_Data/New_Mexico_CIFS_Costs.Rds")) -COSTS <- do.call(rbind,lapply(seq(0,1,by=0.025),function(DISCOUNT){CIFS %>% group_by(Location,Capacity,Cost_Assumption) %>% mutate(NPC=Total/(1+DISCOUNT)^Year) %>% summarize(Discount=DISCOUNT,Costs=sum(NPC))})) -RES <- REV_RES %>% left_join(COSTS) %>% mutate(Profit=Revenue-Costs) -PLOT_DATA <- RES %>% filter(Discount %in% c(0.03,0.05,0.07,0.1)) %>% filter(!is.na(Costs)) -PLOT_DATA -png("Revenue.png",height=6,width=11,units="in",res=900) - -ggplot(PLOT_DATA,aes(x=Year,y=Profit/10^6,color=Discount))+geom_point()+facet_wrap(~Capacity) - - -NPC <- CIFS %>% select(Year,Location,Phase,Capacity,Cost_Assumption,Total) %>% group_by(Year,Location,Cost_Assumption,Total,Phase,Capacity) %>% summarize(Total=mean(Total)) %>% arrange(Location,Phase,Year) %>% mutate(Cost_3=Total/(1+0.03)^Year,Cost_5=Total/(1+0.05)^Year,Cost_7=Total/(1+0.07)^Year,Cost_10=Total/(1+0.1)^Year) %>% ungroup %>% group_by(Location,Phase,Capacity,Cost_Assumption) %>% summarize('3%'=sum(Cost_3),'5%'=sum(Cost_5),"7%"=sum(Cost_7),"10%"=sum(Cost_10)) %>% pivot_longer(c(-Phase,-Location,-Cost_Assumption,-Capacity),names_to="Discount",values_to="CIFS_Cost") +Rev_No_Shipping <- readRDS("Data/Cleaned_Data/Model_Estimates.Rds") +Rev_No_Shipping_2018 <- readRDS("Data/Cleaned_Data/Model_Estimates_2018.Rds") +Rev_Shipping <- readRDS("Data/Cleaned_Data/Model_Estimates_With_Shipping_Costs.Rds") +Rev_Shipping_2018 <- readRDS("Data/Cleaned_Data/Model_Estimates_With_Shipping_Costs_2018.Rds") - PLOT_DATA <- OUTCOME %>% filter(Discount %in% c(0.03,0.05,0.1)) -png("Revenue.png",height=6,width=11,units="in",res=900) - ggplot(PLOT_DATA,aes(x=Year,y=Revenue,color=Discount))+geom_point() -dev.off() +TEST <- Rev_No_Shipping +TEST <- TEST %>% mutate(Current_Profit=Profit/((1+Discount)^(Year-2026))) %>% ungroup +TEST <- TEST %>% group_by(Year,Location,Discount) %>% filter(Current_Profit==max(Current_Profit)) %>% ungroup +TEST %>% select(Profit,Current_Profit,Year,Discount) %>% tail +TEST %>% group_by(Location,Discount) %>% filter(Current_Profit==max(Current_Profit)) %>% print(n=100) %>% +TEMP <- TEST %>% group_by(Location,Discount) %>% select(Year,Capacity,Profit,Current_Profit) %>% mutate(Delay_Profit=(lead(Current_Profit)-Current_Profit)/10^6) +TEMP %>% filter(Year>=2025) %>% summarize(Optimal_Year=sum(ifelse(Current_Profit==max(Current_Profit),Year,0)),Max_Profit=max(Current_Profit),Current_Profit=sum(ifelse(Year==2026,Current_Profit,0)),Yearly_Loss=(Max_Profit-Current_Profit)/(2026-Optimal_Year)/10^6,Next_Years_Loss=(Max_Profit-sum(ifelse(Year==(Optimal_Year+1),Current_Profit,0)))/10^6) +TEST %>% group_by(Location,Discount, +ggplot(TEST %>% filter(Discount==0.03) ,aes(x=Year,y=Current_Profit/10^9,color=Capacity))+geom_line()+facet_wrap(Discount~Location,ncol=1)+scale_x_continuous(breaks=seq(1960,2083,by=5)) -KEY_YEARS <- c(1986,2026,2066) -#Data for reduced and increased fee at 5% and 2026 -FEE_DATA <- rbind(KEY_DATA %>% filter(Year %in% 2026 ,Q>=500,Discount=="5%") %>% mutate(Fee="Current Fee"),KEY_DATA %>% filter(Year %in% 2026 ,Q>=500,Discount=="5%") %>% mutate(Marginal_Value=Marginal_Value/2,Fee="Reduced Fee"),KEY_DATA %>% filter(Year %in% 2026 ,Q>=500,Discount=="5%") %>% mutate(Marginal_Value=2*Marginal_Value,Fee="Increased Fee")) -FEE_DATA$Fee <- factor(FEE_DATA$Fee,levels=c("Reduced Fee","Current Fee","Increased Fee")) +ggplot(TEST %>% filter(Discount==0.05) ,aes(x=Year,y=Current_Profit/10^9,color=Capacity))+geom_line()+facet_wrap(Discount~Location,ncol=1)+scale_x_continuous(breaks=seq(1960,2083,by=5)) +ggplot(TEST %>% filter(Discount %in% c(0.03,0.05,0.07) ,Year>1970) ,aes(x=Year,y=Current_Profit/10^9,color=Capacity))+geom_line()+facet_wrap(Discount~Location,ncol=2)+scale_x_continuous(breaks=seq(1960,2083,by=5)) -DEMAND_CURVE_YEARS <- ggplot(KEY_DATA %>% filter(Year %in% KEY_YEARS ,Q>=500,Discount=="5%") ,aes(x=Q/1000,y=Marginal_Value/1000,group=Year,color=Year))+geom_step(linewidth=1,arrow = arrow(length = unit(0.25, "cm")))+theme_bw()+scale_color_binned(high= "#132B43", low= "#56B1F7",breaks = KEY_YEARS)+scale_y_continuous(breaks=seq(0,2000,by=50))+scale_x_continuous(breaks=seq(0,150,by=5))+ylab("Price ($1000 per ton)")+xlab("Quantity (Thousand Tons)")+theme(text = element_text(size = 16),legend.position = "top")+ guides(color = guide_legend(reverse = FALSE)) -DEMAND_CURVE_YEARS +ggplot(TEST %>% filter(Discount %in% c(0.05) ,Year>1970) ,aes(x=Year,y=Profit/10^9,color=Location))+geom_line()+scale_x_continuous(breaks=seq(1960,2083,by=5)) -DEMAND_CURVE_YEARS_WITH_SHIPPING <- ggplot(KEY_DATA_WITH_SHIPPING %>% filter(Year %in% KEY_YEARS ,Q>=500,Discount=="5%") ,aes(x=Q/1000,y=Marginal_Value/1000,group=Year,color=Year))+geom_step(linewidth=1,arrow = arrow(length = unit(0.25, "cm")))+theme_bw()+scale_color_binned(high= "#132B43", low= "#56B1F7",breaks = KEY_YEARS)+scale_y_continuous(breaks=seq(0,2000,by=50))+scale_x_continuous(breaks=seq(0,150,by=5))+ylab("Price ($1000 per ton)")+xlab("Quantity (Thousand Tons)")+theme(text = element_text(size = 16),legend.position = "top")+ guides(color = guide_legend(reverse = FALSE)) -DEMAND_CURVE_YEARS_WITH_SHIPPING +ggplot(TEST %>% filter(Discount %in% c(0.03,0.05,0.07),Year>1970,Year<2050) ,aes(x=Year,y=Marginal/10^6))+ geom_area(aes(y=ifelse(Marginal>0,Marginal/10^6,0)),fill="lightgreen",alpha=0.6)+ geom_area(aes(y=ifelse(Marginal<0,Marginal/10^6,0)),fill="firebrick2",alpha=0.5) +facet_wrap(Discount~Location,ncol=2,scales = "free_y")+scale_x_continuous(breaks=seq(1960,2083,by=5))+scale_y_continuous(breaks=seq(-300,50,by=50))+geom_hline(yintercept=0)+theme_bw()+ylab("Profit Change (Million Dollars)") - - -DEMAND_CURVE_YEARS_PLOT_FULL <-ggplot(KEY_DATA %>% filter(Year %in% KEY_YEARS,Discount=="5%" ) ,aes(x=Q/1000,y=log(Marginal_Value),group=Year,color=Year))+geom_step(linewidth=2,arrow = arrow(length = unit(0.3, "cm")))+theme_bw()+scale_color_binned(high= "#132B43", low= "#56B1F7",breaks = KEY_YEARS)+ guides(color = guide_legend(reverse = TRUE))+scale_y_continuous(breaks=c(seq(4,20,by=0.5)))+scale_x_continuous(breaks=seq(0,150,by=5))+ylab("Price (Log dollar per Ton)")+xlab("Quantity (Thousand Tons)")+theme(text = element_text(size = 16),legend.position = "top")+ guides(color = guide_legend(reverse = FALSE)) -DEMAND_CURVE_YEARS_PLOT_FULL - -DEMAND_CURVE_DISCOUNT <- ggplot(KEY_DATA %>% filter(Year %in% KEY_YEARS ,Q>=500,Year==2026) ,aes(x=Q/1000,y=Marginal_Value/1000,group=Discount,color=Discount))+geom_step(linewidth=1,arrow = arrow(length = unit(0.1, "cm")))+theme_bw()+scale_y_continuous(breaks=seq(0,2000,by=50))+scale_x_continuous(breaks=seq(0,150,by=5))+ylab("Price ($1000 per ton)")+xlab("Quantity (Thousand Tons)")+theme(text = element_text(size = 16),legend.position = "top")+ guides(color = guide_legend(reverse = FALSE))+scale_color_manual(values = c("tomato1", "tomato2", "tomato3","tomato4")) -DEMAND_CURVE_DISCOUNT - -DEMAND_CURVE_FEE <- ggplot(FEE_DATA ,aes(x=Q/1000,y=Marginal_Value/1000,group=Fee,color=Fee))+geom_step(linewidth=1,arrow = arrow(length = unit(0.25, "cm")))+theme_bw()+scale_y_continuous(breaks=seq(0,2000,by=50))+scale_x_continuous(breaks=seq(0,150,by=5))+ylab("Price ($1000 per ton)")+xlab("Quantity (Thousand Tons)")+theme(text = element_text(size = 16),legend.position = "top")+ guides(color = guide_legend(reverse = FALSE))+scale_color_manual(values = c("forestgreen", "darkblue", "red3")) -DEMAND_CURVE_FEE - - - - -DEMAND_CURVE_FACET <- ggplot(KEY_DATA %>% filter(Year %in% KEY_YEARS ,Q>=500) ,aes(x=Q/1000,y=Marginal_Value/1000,group=Year,color=Year))+geom_step(linewidth=1,arrow = arrow(length = unit(0.08, "cm")))+theme_bw()+scale_color_binned(high= "#132B43", low= "#56B1F7",breaks = KEY_YEARS)+scale_y_continuous(breaks=seq(0,2000,by=200))+scale_x_continuous(breaks=seq(0,150,by=10))+ylab("Price ($1000 per ton)")+xlab("Quantity (Thousand Tons)")+theme(text = element_text(size = 16),legend.position = "top")+ guides(color = guide_legend(reverse = FALSE))+facet_grid(Discount~Year) -DEMAND_CURVE_FACET - -DEMAND_CURVE_YEARS_PLOT_FULL -DEMAND_CURVE_YEARS -DEMAND_CURVE_YEARS_WITH_SHIPPING -DEMAND_CURVE_DISCOUNT -DEMAND_CURVE_FEE -DEMAND_CURVE_FACET - - - -### -CIFS <- rbind(readRDS("Data/Cleaned_Data/Texas_CIFS_Costs.Rds"),readRDS("Data/Cleaned_Data/New_Mexico_CIFS_Costs.Rds")) - -NPC <- CIFS %>% select(Year,Location,Phase,Capacity,Cost_Assumption,Total) %>% group_by(Year,Location,Cost_Assumption,Total,Phase,Capacity) %>% summarize(Total=mean(Total)) %>% arrange(Location,Phase,Year) %>% mutate(Cost_3=Total/(1+0.03)^Year,Cost_5=Total/(1+0.05)^Year,Cost_7=Total/(1+0.07)^Year,Cost_10=Total/(1+0.1)^Year) %>% ungroup %>% group_by(Location,Phase,Capacity,Cost_Assumption) %>% summarize('3%'=sum(Cost_3),'5%'=sum(Cost_5),"7%"=sum(Cost_7),"10%"=sum(Cost_10)) %>% pivot_longer(c(-Phase,-Location,-Cost_Assumption,-Capacity),names_to="Discount",values_to="CIFS_Cost") -TEMP1 <- NPC%>% filter(Location=='Texas') %>% select(-CIFS_Cost) %>% mutate(Phase="Extended",Capacity=10^5) %>% unique -TEMP2 <- NPC %>% filter(Location=='Texas') %>% group_by(Location,Discount,Cost_Assumption) %>% summarize(CAP_CHANGE=max(Capacity)-min(Capacity),ST_COST=min(CIFS_Cost),END_COST=max(CIFS_Cost),SLOPE=(END_COST-ST_COST)/CAP_CHANGE,CIFS_Cost=ST_COST+SLOPE*(10^5-5000)) %>% select(Location,Discount,CIFS_Cost,Cost_Assumption) -TEMP <- TEMP1 %>% left_join(TEMP2) -NPC <- rbind(NPC,TEMP ) -NPC <- NPC %>% filter(Cost_Assumption=="Average") %>% select(-Cost_Assumption) - -RETURN_DATA <- KEY_DATA %>% left_join(NPC) %>% filter(Q<=Capacity) %>% group_by(Year,Discount,Location,Phase,Capacity,CIFS_Cost) %>% summarize(Q=sum(Total_Tons),Revenue=sum(Marginal_Value*Total_Tons),P=Revenue/Q,Profit=sum(Revenue-CIFS_Cost)) %>% ungroup - -RETURN_DATA <- RETURN_DATA %>% group_by(Discount,Location,Phase,Capacity,CIFS_Cost) %>% arrange(Discount,Location,Phase,Capacity,CIFS_Cost,Year) %>% mutate(FOC=(lead(Profit)-Profit)-Profit*(parse_number(Discount)/100)) -####This is a plot of the supply curve in orange UNDER A COMPETITVE MARKET. Add in the monoplosty wating function next (value is gained by waiting another year. -#Added cost per unit -#M <- (3373397183-1401068722)/(40000-5000)/1000 -#Added cost for intial project -#C <- 1401068722/5000/1000 -#EXTENDED_COST <- (M*(100000-5000) - +C)*1000 -#RETURN_DATA - -#RETURN_DATA <- RETURN_DATA %>% mutate(size=ifelse(Phase=="Partial","5,000 MTU Capacity (Phase 1)","40,000 MTU Capacity (Full Build Out)")) -RETURN_DATA %>% filter(Phase=="Extended") -RETURN_DATA %>% pull(Capacity) %>% unique -RETURN_DATA %>% group_by(Year,Capacity,Discount) %>% filter(n()>1) %>% arrange(Year,Capacity) -png("Revenue.png",height=6,width=11,units="in",res=900) - ggplot(RETURN_DATA %>% filter(Discount=='5%'),aes(x=Year,y=Profit/10^9,color=Location))+geom_point()+scale_x_continuous(breaks=seq(1960,2085,by=5)) -dev.off() - - -FOC_PLOT <- ggplot(RETURN_DATA,aes(x=Year,y=FOC/10^6,group=Discount,color=Discount))+facet_wrap(~Capacity,ncol=1,scales="free")+scale_x_continuous(breaks=seq(1960,2083,by=10))+geom_point(size=0.5)+geom_step(linewidth=0.75)+ geom_hline(yintercept = 0, color = "black", linetype = "solid", size = 1)+theme(text = element_text(size = 16),legend.position = "top")+ylab("Million Dollars")png("FOC_PLOT.png",width=8,height=18, units="in",res=1800) -FOC_PLOT -dev.off() - -DEMAND_CURVE_YEARS+geom_segment(aes(x = 0, y = C, xend = 5, yend = C),color="orange")+geom_segment(aes(x = 5, y = C, xend = 5, yend = M),color="orange")+geom_segment(aes(x = 5, y = M, xend = 140, yend = M),color="orange",arrow = arrow(length = unit(0.4, "cm"))) diff --git a/Data_Proc.r b/Data_Proc.r deleted file mode 100644 index 417ad72..0000000 --- a/Data_Proc.r +++ /dev/null @@ -1,251 +0,0 @@ -#A script which attempts to pull in all data, and create a data frame with the maximum revenue values for each facility, year and discount rate. The output can then be used to make figures and graphs -library(tidyverse) -library(parallel) - NCORES <- detectCores()-1 -library(lpSolve) #For solving discrete value maximization for the power plants -####Manual inputs -#DISCOUNT_RATE_LIST <- seq(0,1,by=0.0025) -DISCOUNT_RATE_LIST <- c(0.03,0.0325,0.035,0.0375,0.04,0.045,0.0475,0.05,0.07,0.1) -SHIPPING_COST_PER_TON <- 1.2874*26000 #Inflation adjusted from New Mexico Report -CV1 <- 1.3074*(10607030-1060703) #Data from Texas Report, converted from 2018 to Dec 2025 dollars -#CV2 <- 1.2874*(6984013-1117442) #Data from New Mexico Report, Converted from 2019 to Dec 2025 -CV2 <- 1.2874*(6984013) #Data from New Mexico Report, Converted from 2019 to Dec 2025 -CV2_OP <- 1.2874*1117442 - -CV3 <- mean(CV1,CV2) #Average of the two - -###################################Cost results -#source("Cost_Data_Proc.r") -CIFS <- rbind(readRDS("Data/Cleaned_Data/Texas_CIFS_Costs.Rds"),readRDS("Data/Cleaned_Data/New_Mexico_CIFS_Costs.Rds")) -#Adjust for inflation -CIFS[,-1:-5] <- 1.2874*CIFS[,-1:-5] -COSTS <- do.call(rbind,lapply(DISCOUNT_RATE_LIST ,function(DISCOUNT){CIFS %>% group_by(Location,Capacity,Cost_Assumption) %>% mutate(NPC=Total/(1+DISCOUNT)^Year) %>% summarize(Discount=DISCOUNT,Costs=sum(NPC))})) %>% ungroup - -TEMP <- COSTS%>% group_by(Location,Cost_Assumption,Discount) %>% summarize(ST_COST=min(Costs),ST_CAPACITY=min(Capacity),M_COST=(max(Costs)-min(Costs))/(max(Capacity)-min(Capacity))) %>% ungroup -CAPACITY_INCREMENT <- 5000 -#rm(COST_DATA) -for(i in 1:nrow(TEMP)){ - LOC <- as.character(TEMP[i,"Location"]) - COST_LEVEL <- as.character(TEMP[i,"Cost_Assumption"]) - DISCOUNT <- as.numeric(TEMP[i,"Discount"]) - ST_CAP <- as.numeric(TEMP[i,"ST_CAPACITY"]) - ST_COST <- as.numeric(TEMP[i,"ST_COST"]) - COST_SLOPE <- as.numeric(TEMP[i,]$M_COST) - CAPACITY <- seq(ST_CAP,160000,by=5000) - COST <- ST_COST+COST_SLOPE*CAPACITY_INCREMENT*(0:(length(CAPACITY)-1)) - - C_RES <- cbind(CAPACITY,COST) %>% as_tibble %>% mutate(Location=LOC,Cost_Assumption=COST_LEVEL,Discount=DISCOUNT) %>% select(Location,Cost_Assumption,Discount,Capacity=CAPACITY,Cost=COST) - if(!exists("COST_DATA")){COST_DATA <- C_RES}else{COST_DATA<- rbind(COST_DATA,C_RES)} -} -saveRDS(COST_DATA ,"Data/Cleaned_Data/All_CIFS_Discounted_Costs.Rds") -#COST_DATA <- COST_DATA %>% filter(Cost_Assumption=='Average') %>% select(-Cost_Assumption) %>% unique -COST_DATA <- COST_DATA %>% filter(Cost_Assumption=='High') %>% select(-Cost_Assumption) %>% unique - - - -##All unique capacity levels that the revenues need to be calculated for -CAPACITY_LIST <- COST_DATA %>% pull(Capacity) %>% unique - -### -TOTAL <- read_csv("Data/Raw_Data/Curie_Spent_Fuel_Site_Totals.csv") %>% mutate(OP_YEAR=year(Op_Date_Min),CLOSE_YEAR=year(Close_Date_Max))%>% select(Facility,Total_Assemblies,Total_Tons,OP_YEAR,CLOSE_YEAR) -Duane Arnold -Nine Mile Point -Ginna -FACILITY_LIST <- TOTAL %>% pull(Facility) -#https://www.nrc.gov/reactors/operating/licensing/renewal/subsequent-license-renewal -SUBMITTED <-rbind(c(FACILITY_LIST[str_detect(FACILITY_LIST,"Duane*" )],2025), -c(FACILITY_LIST[str_detect(FACILITY_LIST,"Nine Mile*" )],2026), -c(FACILITY_LIST[str_detect(FACILITY_LIST,"Ginna*" )],2026), -c(FACILITY_LIST[str_detect(FACILITY_LIST,"Cooper*" )],2026), -c(FACILITY_LIST[str_detect(FACILITY_LIST,"Farley*" )],2027), -c(FACILITY_LIST[str_detect(FACILITY_LIST,"Prairie*" )],2027), -c(FACILITY_LIST[str_detect(FACILITY_LIST,"Brunswick*" )],2027), -c(FACILITY_LIST[str_detect(FACILITY_LIST,"Cook" )],2027), -c(FACILITY_LIST[str_detect(FACILITY_LIST,"Hope" )],2027), -c(FACILITY_LIST[str_detect(FACILITY_LIST,"Salem" )],2027), -c(FACILITY_LIST[str_detect(FACILITY_LIST,"Perry" )],2027), -c(FACILITY_LIST[str_detect(FACILITY_LIST,"Millstone" )],2028), -c(FACILITY_LIST[str_detect(FACILITY_LIST,"Palisades" )],2028), -c(FACILITY_LIST[str_detect(FACILITY_LIST,"Beaver" )],2028), -c(FACILITY_LIST[str_detect(FACILITY_LIST,"Callaway" )],2029), -c(FACILITY_LIST[str_detect(FACILITY_LIST,"Three Mile Island" )],2029), -c(FACILITY_LIST[str_detect(FACILITY_LIST,"Davis-Besse" )],2029), -c(FACILITY_LIST[str_detect(FACILITY_LIST,"Wolf" )],2030), -c(FACILITY_LIST[str_detect(FACILITY_LIST,"Lucie" )],2021), -c(FACILITY_LIST[str_detect(FACILITY_LIST,"Robinson" )],2025), -c(FACILITY_LIST[str_detect(FACILITY_LIST,"Hatch" )],2025)) -SUBMITTED <- SUBMITTED %>% as_tibble -colnames(SUBMITTED ) <- c("Facility","App_Date","Status") -SUBMITTED <- SUBMITTED%>% mutate(Status="Applied",App_Date=as.numeric(App_Date)) %>% select(Facility,Status,App_Date) - -#Issued -RENEWED <- rbind(c(FACILITY_LIST[str_detect(FACILITY_LIST,"Turkey" )],"Granted",2018,2033), -c(FACILITY_LIST[str_detect(FACILITY_LIST,"Peach" )],"Granted",2019,2034), -c(FACILITY_LIST[str_detect(FACILITY_LIST,"Surry" )],"Granted",2020,2033), -c(FACILITY_LIST[str_detect(FACILITY_LIST,"North" )],"Granted",2021,2040), -c(FACILITY_LIST[str_detect(FACILITY_LIST,"Monticello" )],"Granted",2024,2030), -c(FACILITY_LIST[str_detect(FACILITY_LIST,"Oconee" )],"Granted",2025,2034), -c(FACILITY_LIST[str_detect(FACILITY_LIST,"Summer" )],"Granted",2025,2042), -c(FACILITY_LIST[str_detect(FACILITY_LIST,"Beach" )],"Granted",2025,2033), -c(FACILITY_LIST[str_detect(FACILITY_LIST,"Browns" )],"Granted",2025,2036), -c(FACILITY_LIST[str_detect(FACILITY_LIST,"Dresden" )],"Granted",2025,2031)) -RENEWED <- RENEWED %>% as_tibble -colnames(RENEWED) <- c("Facility","Status","App_Date","Op_Date") -RENEWED <- RENEWED %>% mutate(Op_Date=as.numeric(Op_Date),App_Date=as.numeric(App_Date)) -AVG_LENGTH <- RENEWED %>% mutate(DIFF=Op_Date-App_Date) %>% pull(DIFF) %>% mean %>% round -SUBMITTED <- SUBMITTED %>% mutate(Op_Date=App_Date+AVG_LENGTH) -UPDATE <- rbind(RENEWED,SUBMITTED ) -UPDATE - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -str_detect() -TOTAL %>% print(n=50) -#!!!!!!!!!!!!#####TEST OF VARIABLITY IN CLOSING YEAR, CHANGE OR FORMALIZE LATER -#TOTAL <- TOTAL %>% mutate(CLOSE_YEAR=ifelse(CLOSE_YEAR>2018,CLOSE_YEAR+20,CLOSE_YEAR)) - -###########Series of functions to calculate the gross consumer surplus from a CIFS for each facility, in each year from 1960 to 2082. - -#Function to find the net present revenue of a facility,given a discount rate, and the current year, and year the facility will close. This is the sum of discounted costs that WOULD have taken place if the facility was not built -VALUE_ADD <- function(r,CURRENT_YEAR,CLOSE_YEAR){ - Years_Until_Close <- max(CLOSE_YEAR-CURRENT_YEAR+1,0) - VALUES <- (1+r)^-(1:10^4) - if(Years_Until_Close==0){return(sum(VALUES))} else{return(sum(VALUES[-1:-Years_Until_Close]))} -} -#A function to extend the revenues calculations to the closure date of all of the facilities. -VALUE_ADD_SINGLE_YEAR <- function(r,CURRENT_YEAR,CLOSE_YEARS){return(sapply(CLOSE_YEARS,function(x){VALUE_ADD(r,CURRENT_YEAR,x)}))} -#A function to extend the calculation of the net present revenues of each facility to all years of interest. That is what is the NPV of building the facility in each year, for each facility. -NPV_CALC <- function(r,DATA=TOTAL,YEARS=1960:2083){ - Facility <- pull(DATA,Facility) - RES <- cbind(Facility,do.call(cbind,lapply(YEARS,function(x){VALUE_ADD_SINGLE_YEAR(r,x,DATA$CLOSE_YEAR)}))) -colnames(RES) <- c("Facility",YEARS) - RES <- as_tibble(RES) %>% pivot_longer(cols=-Facility,names_to="Year",values_to="Revenue") %>% arrange(Year) %>% mutate(Year=parse_number(Year),Revenue=parse_number(Revenue),Discount=r) - return(RES) -} -#A function which returns a list, of net present revenue calculation tables (facility by year) for a range of possible discount rates. This allows for the results to be quickly looked up, when we want to adjust the time value of money. These results combine costs savings to calculate NPV -MULTI_DISCOUNT_RATE_NPV <- function(DISCOUNT_INCREMENT,DATA=TOTAL,YEARS=1960:2083,DOLLARS_SAVED_PER_YEAR){ - NCORES <- detectCores()-1 - RES <- mclapply(DISCOUNT_INCREMENT,NPV_CALC,mc.cores = NCORES) - RES <- do.call(rbind,RES) %>% mutate(Revenue=Revenue*DOLLARS_SAVED_PER_YEAR) - return(RES) -} - - -TOTAL_VALUE_METRICS <- MULTI_DISCOUNT_RATE_NPV(DISCOUNT_RATE_LIST ,DOLLARS_SAVED_PER_YEAR=CV2) -#TOTAL_VALUE_METRICS <- TOTAL_VALUE_METRICS %>% filter(!(Facility %in% c("Palo Verde","Vogtle")))#These facilities always have a negative NPV by the end date 2083 -TOTAL_VALUE_METRICS <- TOTAL_VALUE_METRICS %>% left_join(read_csv("Data/Raw_Data/Curie_Spent_Fuel_Site_Totals.csv")) %>% select(Year,Facility,Discount,Total_Tons,Revenue) - -#Function to find the maximum Revenue BUT it does not take dollars just relative savings -MAX_REV <- function(YEAR,TBL,DISCOUNT,CIFS_SIZE,SHIPPING_COST=0){ - TBL <- TBL %>% filter(Year==YEAR,Discount==DISCOUNT,Revenue/Total_Tons>SHIPPING_COST) - REV <- TBL %>% pull(Revenue) - VOL <- TBL %>% pull(Total_Tons) - RES <- lp(direction = "max", objective.in = REV, const.mat = matrix(VOL, nrow = 1),const.dir = "<=", const.rhs = CIFS_SIZE, all.bin = TRUE) - return(RES[[11]]) - } - -#Calculate all results for the start year to the end year. Discount rate, and capacities are fixed. -YEARLY_RESULTS <- function(DATA,DISCOUNT,CAPACITY,SHIPPING_COSTS=0){ - RES <- cbind(1960:2083,sapply(1960:2083, function(x){MAX_REV(YEAR=x,TBL=DATA,DISCOUNT,CAPACITY,SHIPPING_COSTS)})) %>% as_tibble - colnames(RES) <- c("Year","Revenue") - RES <- RES %>% mutate(Discount=DISCOUNT,Capacity=CAPACITY) %>% select(Year,Discount,Capacity,Revenue) %>% mutate(Shipping_Cost_Cuttoff=SHIPPING_COSTS) - return(RES) -} -#Calculate and individual facilities rate, for all discount rates, and all years - - Rev_No_Shipping <- do.call(rbind,mclapply(CAPACITY_LIST,function(CAPACITY){do.call(rbind,lapply(DISCOUNT_RATE_LIST,function(x){YEARLY_RESULTS(TOTAL_VALUE_METRICS,x,CAPACITY,0)}))},mc.cores=min(length(CAPACITY_LIST),NCORES))) - saveRDS(Rev_No_Shipping,"Data/Cleaned_Data/Model_Estimates.Rds") -TEST <- Rev_No_Shipping %>% left_join(COST_DATA) %>% mutate(Profit=Revenue-Cost) %>% group_by(Discount,Capacity,Location) %>% mutate(Time_Benefit=(1-Discount)*(lead(Profit)-Profit),Op_Cost=Discount*Profit,Marginal=Time_Benefit-Op_Cost) %>% unique -TEST <- TEST %>% mutate(Current_Profit=Profit/((1+Discount)^(Year-2026))) %>% ungroup -TEST <- TEST %>% group_by(Year,Location,Discount) %>% filter(Current_Profit==max(Current_Profit)) %>% ungroup -TEST %>% select(Profit,Current_Profit,Year,Discount) %>% tail -TEST %>% group_by(Location,Discount) %>% filter(Current_Profit==max(Current_Profit)) %>% print(n=100) %>% -TEMP <- TEST %>% group_by(Location,Discount) %>% select(Year,Capacity,Profit,Current_Profit) %>% mutate(Delay_Profit=(lead(Current_Profit)-Current_Profit)/10^6) -TEMP %>% filter(Year>=2025) %>% summarize(Optimal_Year=sum(ifelse(Current_Profit==max(Current_Profit),Year,0)),Max_Profit=max(Current_Profit),Current_Profit=sum(ifelse(Year==2026,Current_Profit,0)),Yearly_Loss=(Max_Profit-Current_Profit)/(2026-Optimal_Year)/10^6,Next_Years_Loss=(Max_Profit-sum(ifelse(Year==(Optimal_Year+1),Current_Profit,0)))/10^6) -TEST %>% group_by(Location,Discount, - -ggplot(TEST %>% filter(Discount==0.03) ,aes(x=Year,y=Current_Profit/10^9,color=Capacity))+geom_line()+facet_wrap(Discount~Location,ncol=1)+scale_x_continuous(breaks=seq(1960,2083,by=5)) - -ggplot(TEST %>% filter(Discount==0.05) ,aes(x=Year,y=Current_Profit/10^9,color=Capacity))+geom_line()+facet_wrap(Discount~Location,ncol=1)+scale_x_continuous(breaks=seq(1960,2083,by=5)) - -ggplot(TEST %>% filter(Discount %in% c(0.03,0.05,0.07) ,Year>1970) ,aes(x=Year,y=Current_Profit/10^9,color=Capacity))+geom_line()+facet_wrap(Discount~Location,ncol=2)+scale_x_continuous(breaks=seq(1960,2083,by=5)) - -ggplot(TEST %>% filter(Discount %in% c(0.05) ,Year>1970) ,aes(x=Year,y=Profit/10^9,color=Location))+geom_line()+scale_x_continuous(breaks=seq(1960,2083,by=5)) - -ggplot(TEST %>% filter(Discount %in% c(0.03,0.05,0.07),Year>1970,Year<2050) ,aes(x=Year,y=Marginal/10^6))+ geom_area(aes(y=ifelse(Marginal>0,Marginal/10^6,0)),fill="lightgreen",alpha=0.6)+ geom_area(aes(y=ifelse(Marginal<0,Marginal/10^6,0)),fill="firebrick2",alpha=0.5) +facet_wrap(Discount~Location,ncol=2,scales = "free_y")+scale_x_continuous(breaks=seq(1960,2083,by=5))+scale_y_continuous(breaks=seq(-300,50,by=50))+geom_hline(yintercept=0)+theme_bw()+ylab("Profit Change (Million Dollars)") - -+geom_line()+facet_wrap(Discount~Location,ncol=2,scales = "free_y")+scale_x_continuous(breaks=seq(1960,2083,by=5))+geom_hline(yintercept=0)+theme_bw()+geom_area(fill="green") - - -ggplot(TEST %>% filter(Year>2000) ,aes(x=Year,y=Marginal/10^9,color=Capacity))+geom_line()+facet_wrap(Discount~Location,ncol=2,scales = "free_y")+scale_x_continuous(breaks=seq(1960,2083,by=5)) - - - - -ggplot(TEST ,aes(x=Year,y=Current_Profit/10^9,color=Capacity))+geom_line()+facet_wrap(Location~Discount) - - - - -TEST2 <- TEST %>% filter(Discount==0.05,Location=='Texas') - -ggplot(TEST2,aes(x=Year,y=Marginal,group=as.factor(Capacity),color=as.factor(Capacity)))+geom_line() - - #Calculate the net present revenue of all facilities that can CURRENTLY be shipped to the site with a profit. Note that the total Revenue is more important because most of the sites will be willing to pay for the right to CIFS storage in the future even if the shipping costs are too high presently. This result is used to show what might be the current ideal facility size, even if future expansion is expected to maximize profit. - Rev_Shipping <- do.call(rbind,mclapply(CAPACITY_LIST,function(CAPACITY){do.call(rbind,lapply(DISCOUNT_RATE_LIST,function(x){YEARLY_RESULTS(TOTAL_VALUE_METRICS,x,CAPACITY,SHIPPING_COST_PER_TON )}))},mc.cores=min(length(CAPACITY_LIST),NCORES))) - -Revenue_Results <- rbind(Rev_No_Shipping,Rev_Shipping) %>% mutate(Type=ifelse(Shipping_Cost_Cuttoff==0,"Revenue","Revenue_Shipping")) %>% pivot_wider(values_from=Revenue,names_from=Type) %>% select(-Shipping_Cost_Cuttoff) %>% group_by(Year,Discount,Capacity) %>% summarize(Revenue=mean(Revenue,na.rm=TRUE),Revenue_Shipping=mean(Revenue_Shipping,na.rm=TRUE)) %>% ungroup -Revenue_Results -#################################################### -COSTS <- rbind(COSTS,EXTENDED_CAPACITY_NEW_MEXICO) -COST_DATA %>% pull(Cost_Assumption) %>% unique -saveRDS(COSTS,"Data/Cleaned_Data/All_CIFS_Discounted_Costs.Rds") -COSTS %>% group_by(Location,Capac -CIFS_Data <- Revenue_Results %>% left_join(COSTS) %>% mutate(Profit=Revenue-Costs,Profit_Shipping=Revenue_Shipping-Costs) %>% select(Year,Location,Capacity,Discount,Revenue,Costs,Profit,Revenue_Shipping,Profit_Shipping,everything()) -#CIFS_Data <- -CIFS_Data <- CIFS_Data %>% group_by(Location,Capacity,Discount) %>% arrange(Location,Capacity,Discount,Year) %>% mutate(Time_Benefit=lead(Profit)-Profit,Op_Cost=Profit*Discount,Marginal=Time_Benefit-Op_Cost) %>% ungroup -TEMP <- CIFS_Data %>% filter(Discount==0.05) %>% mutate(Time_Benefit=lead(Profit)-Profit,Op_Cost=Profit*Discount,Marginal=Time_Benefit-Op_Cost) - -CIFS_Data <- CIFS_Data %>% group_by(Location,Phase,Capacity,Discount) %>% mutate(Time_Benefit=(1-Discount)*(lead(Profit)-Profit),Op_Cost=Profit*Discount,Marginal=Time_Benefit-Op_Cost) -((-3504542954)-(-3530861857))-0.05*(-3530861857 ) -STARTING_YEARS <- CIFS_Data %>% group_by(Location,Phase,Capacity,Discount) %>% mutate(PROFITABLE=Profit>0 & Marginal<=0) %>% filter(PROFITABLE) %>% filter(Year==min(Year)) %>% select(Location,Phase,Capacity,Discount,Start_Year=Year,Profit) %>% ungroup -CIFS_Data %>% mutate(Current_Profit=Profit/((1+Discount)^(Year-2026))) %>% group_by(Location,Discount,Capacity) %>% filter(Current_Profit==max(Current_Profit)) %>% select(Start_Year=Year,Location,Phase,Capacity,Current_Profit) %>% ungroup %>% arrange(Discount,Location,desc(Current_Profit)) %>% select(Discount,Location,Phase,Start_Year,Current_Profit) - -############# -ggplot(CIFS_Data ,aes(x=Year,y=Marginal/10^6,group=Capacity,color=as.factor(Capacity)))+geom_line()+facet_wrap(~Discount,ncol=1)+scale_x_continuous(breaks=seq(1960,2083,by=5)) - -ggplot(CIFS_Data ,aes(x=Year,y=Marginal,group=Capacity,color=as.factor(Capacity)))+geom_line()+facet_grid(~Discount) - - - - -dir.create("./Results",showWarnings=FALSE) -saveRDS(TOTAL_VALUE_METRICS,"./Results/Storage_Values_by_Facility_and_Variable_Discounts.Rds") diff --git a/Run_All.sh b/Run_All.sh new file mode 100644 index 0000000..b48c7b3 --- /dev/null +++ b/Run_All.sh @@ -0,0 +1,2 @@ +Rscript ./Scripts/1_Raw_Cost_Data_Clean.r +Rscript ./Scripts/2_Compiled_Results_Data.r diff --git a/Cost_Data_Proc.r b/Scripts/1_Raw_Cost_Data_Clean.r similarity index 80% rename from Cost_Data_Proc.r rename to Scripts/1_Raw_Cost_Data_Clean.r index d4377c1..5fcebe2 100644 --- a/Cost_Data_Proc.r +++ b/Scripts/1_Raw_Cost_Data_Clean.r @@ -1,16 +1,21 @@ library(tidyverse) -######CIFS costs -dir.create("./Data/Cleaned_Data",recursive=TRUE,showWarnings=FALSE) -read_csv("./Data/Raw_Data/Cost_Tables/Texas/Table_C-3_Undiscounted_Cost_Estimates_Phase_1_Low.csv") +INFLATION_ADJUST <- 1.2874 #Adjust cost from 2019 to 2026 +#################CIFS costs +#Combine all the cost files from Texas NRC Environmental Report CIFS_TEXAS <- rbind( read_csv("./Data/Raw_Data/Cost_Tables/Texas/Table_C-3_Undiscounted_Cost_Estimates_Phase_1_Low.csv") %>% mutate(Phase='Partial',Cost_Assumption="Low"),read_csv("./Data/Raw_Data/Cost_Tables/Texas/Table_C-4_Undiscounted_Cost_Estimates_Phase_1_High.csv") %>% mutate(Phase='Partial',Cost_Assumption="High"), read_csv("./Data/Raw_Data/Cost_Tables/Texas/Table_C-5_Undiscounted_Cost_Estimates_Full_Low.csv") %>% mutate(Phase='Full',Cost_Assumption="Low"), read_csv("./Data/Raw_Data/Cost_Tables/Texas/Table_C-6_Undiscounted_Cost_Estimates_Full_High.csv") %>% mutate(Phase='Full',Cost_Assumption="High")) %>% mutate(Location="Texas",Capacity=ifelse(Phase=='Partial',5000,8*5000)) %>% select(Year,Location,Phase,Capacity,Cost_Assumption,Total,everything()) - +#Calcualte an average value CIFS_TEXAS <- rbind(CIFS_TEXAS,CIFS_TEXAS %>% group_by(Year,Location,Phase,Capacity) %>% summarize(Cost_Assumption="Average",Total=mean(Total),Construction=mean(Construction),Transportation_to_CISF=mean(Transportation_to_CISF),Operations=mean(Operations),Transportation_to_Repository=mean(Transportation_to_Repository),Decommissioning=mean(Decommissioning)) %>% ungroup) -#CIFS_TEXAS %>% group_by(Year,Phase, +CIFS_TEXAS[,-1:-5] <- INFLATION_ADJUST*CIFS_TEXAS[,-1:-5] #Adjust for inflation -saveRDS(CIFS_TEXAS,"Data/Cleaned_Data/Texas_CIFS_Costs.Rds") +#Combine all the cost files from New Mexico NRC Environmental Report CIFS_NEW_MEXICO <- rbind( read_csv("./Data/Raw_Data/Cost_Tables/New_Mexico/Table_C-3_Undiscounted_Cost_Estimates_Phase_1_Low.csv") %>% mutate(Phase='Partial',Cost_Assumption="Low"),read_csv("./Data/Raw_Data/Cost_Tables/New_Mexico/Table_C-4_Undiscounted_Cost_Estimates_Phase_1_High.csv") %>% mutate(Phase='Partial',Cost_Assumption="High"), read_csv("./Data/Raw_Data/Cost_Tables/New_Mexico/Table_C-5_Undiscounted_Cost_Estimates_Full_Low.csv") %>% mutate(Phase='Full',Cost_Assumption="Low"), read_csv("./Data/Raw_Data/Cost_Tables/New_Mexico/Table_C-6_Undiscounted_Cost_Estimates_Full_High.csv") %>% mutate(Phase='Full',Cost_Assumption="High"))%>% mutate(Location="New Mexico",Capacity=ifelse(Phase=='Partial',8680,8680+5000*19))%>% select(Year,Location,Phase,Capacity,Cost_Assumption,Total,everything()) CIFS_NEW_MEXICO <- rbind(CIFS_NEW_MEXICO,CIFS_NEW_MEXICO %>% group_by(Year,Location,Phase,Capacity) %>% summarize(Cost_Assumption="Average",Total=mean(Total),Construction=mean(Construction),Transportation_to_CISF=mean(Transportation_to_CISF),Operations=mean(Operations),Transportation_to_Repository=mean(Transportation_to_Repository),Decommissioning=mean(Decommissioning)) %>% ungroup) -saveRDS(CIFS_NEW_MEXICO,"Data/Cleaned_Data/New_Mexico_CIFS_Costs.Rds") +CIFS_NEW_MEXICO[,-1:-5] <- INFLATION_ADJUST*CIFS_NEW_MEXICO[,-1:-5] #Adjust for inflation +#Save results +dir.create("./Data/Cleaned_Data",recursive=TRUE,showWarnings=FALSE) + saveRDS(CIFS_TEXAS,"Data/Cleaned_Data/Texas_CIFS_Costs.Rds") + saveRDS(CIFS_NEW_MEXICO,"Data/Cleaned_Data/New_Mexico_CIFS_Costs.Rds") + diff --git a/Scripts/2_Compiled_Results_Data.r b/Scripts/2_Compiled_Results_Data.r new file mode 100644 index 0000000..c0f16c2 --- /dev/null +++ b/Scripts/2_Compiled_Results_Data.r @@ -0,0 +1,133 @@ +#A script which attempts to pull in all data, and create a data frame with the maximum revenue values for each facility, year and discount rate. The output can then be used to make figures and graphs +library(tidyverse) +library(parallel) + NCORES <- detectCores()-1 +library(lpSolve) #For solving discrete value maximization for the power plants +####Manual inputs +#DISCOUNT_RATE_LIST <- seq(0,1,by=0.0025) + #Range of discount rates to calculate in the model. Each facility will have each rate calculated, so more values slows the results but allows for more discount rates to be reported in the findings. +DISCOUNT_RATE_LIST <- c(0.03,0.0325,0.035,0.0375,0.04,0.045,0.0475,0.05,0.07,0.1) + #The cost per ton of shipping uranium, used to see what can be shipped on day one of the project. +SHIPPING_COST_PER_TON <- 1.2874*26000 #Inflation adjusted from New Mexico Report + #The savings per year of having a CIFS at a served reactor (cost to house the CIFS). +CV <- 1.2874*(6984013) #Data from New Mexico Report, Converted from 2019 to Dec 2025 + #Locations to save results +RES_DIR <- "./Data/Results/" +INTERMEDIATE_DIR <- "./Data/Results/Separate_Costs_and_Benefits_Data/" + #Directory where the individual CIFS project plan cost data can be found. + #This has data has already been combined with the original files (low,high) and inflation adjusted. +CIFS_INDIVIDUAL_COST_DATA_DIR <- 'Data/Cleaned_Data/' +#Create any need save locations + dir.create(RES_DIR,recursive=TRUE,showWarnings=FALSE) + dir.create(INTERMEDIATE_DIR,recursive=TRUE,showWarnings=FALSE) + +###################################Cost results +CIFS <- rbind(readRDS(paste0(CIFS_INDIVIDUAL_COST_DATA_DIR,"Texas_CIFS_Costs.Rds")),readRDS(paste0(CIFS_INDIVIDUAL_COST_DATA_DIR,"New_Mexico_CIFS_Costs.Rds"))) +#Adjust for inflation +#CIFS[,-1:-5] <- 1.2874*CIFS[,-1:-5] +COSTS <- do.call(rbind,lapply(DISCOUNT_RATE_LIST ,function(DISCOUNT){CIFS %>% group_by(Location,Capacity,Cost_Assumption) %>% mutate(NPC=Total/(1+DISCOUNT)^Year) %>% summarize(Discount=DISCOUNT,Costs=sum(NPC))})) %>% ungroup + +TEMP <- COSTS%>% group_by(Location,Cost_Assumption,Discount) %>% summarize(ST_COST=min(Costs),ST_CAPACITY=min(Capacity),M_COST=(max(Costs)-min(Costs))/(max(Capacity)-min(Capacity))) %>% ungroup +CAPACITY_INCREMENT <- 5000 +#rm(COST_DATA) +for(i in 1:nrow(TEMP)){ + LOC <- as.character(TEMP[i,"Location"]) + COST_LEVEL <- as.character(TEMP[i,"Cost_Assumption"]) + DISCOUNT <- as.numeric(TEMP[i,"Discount"]) + ST_CAP <- as.numeric(TEMP[i,"ST_CAPACITY"]) + ST_COST <- as.numeric(TEMP[i,"ST_COST"]) + COST_SLOPE <- as.numeric(TEMP[i,]$M_COST) + CAPACITY <- seq(ST_CAP,160000,by=5000) + COST <- ST_COST+COST_SLOPE*CAPACITY_INCREMENT*(0:(length(CAPACITY)-1)) + + C_RES <- cbind(CAPACITY,COST) %>% as_tibble %>% mutate(Location=LOC,Cost_Assumption=COST_LEVEL,Discount=DISCOUNT) %>% select(Location,Cost_Assumption,Discount,Capacity=CAPACITY,Cost=COST) + if(!exists("COST_DATA")){COST_DATA <- C_RES}else{COST_DATA<- rbind(COST_DATA,C_RES)} +} +saveRDS(COST_DATA ,paste0(INTERMEDIATE_DIR,"All_CIFS_Discounted_Costs.Rds")) +#COST_DATA <- COST_DATA %>% filter(Cost_Assumption=='Average') %>% select(-Cost_Assumption) %>% unique +COST_DATA <- COST_DATA %>% filter(Cost_Assumption=='High') %>% select(-Cost_Assumption) %>% unique + + +##All unique capacity levels that the revenues need to be calculated for +CAPACITY_LIST <- COST_DATA %>% pull(Capacity) %>% unique + +### +TOTAL <- read_csv("Data/Raw_Data/Curie_Spent_Fuel_Site_Totals.csv") %>% mutate(OP_YEAR=year(Op_Date_Min),CLOSE_YEAR=year(Close_Date_Max))%>% select(Facility,Total_Assemblies,Total_Tons,OP_YEAR,CLOSE_YEAR) +FACILITY_LIST <- TOTAL %>% pull(Facility) +#https://www.nrc.gov/reactors/operating/licensing/renewal/subsequent-license-renewal +SUBMITTED <-rbind(c(FACILITY_LIST[str_detect(FACILITY_LIST,"Duane*" )],2025), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Nine Mile*" )],2026), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Ginna*" )],2026), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Cooper*" )],2026), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Farley*" )],2027), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Prairie*" )],2027), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Brunswick*" )],2027), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Cook" )],2027), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Hope" )],2027), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Salem" )],2027), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Perry" )],2027), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Millstone" )],2028), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Palisades" )],2028), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Beaver" )],2028), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Callaway" )],2029), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Three Mile Island" )],2029), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Davis-Besse" )],2029), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Wolf" )],2030), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Lucie" )],2021), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Robinson" )],2025), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Hatch" )],2025)) +SUBMITTED <- SUBMITTED %>% as_tibble +colnames(SUBMITTED ) <- c("Facility","App_Date","Status") +SUBMITTED <- SUBMITTED%>% mutate(Status="Applied",App_Date=as.numeric(App_Date)) %>% select(Facility,Status,App_Date) + +#Issued +RENEWED <- rbind(c(FACILITY_LIST[str_detect(FACILITY_LIST,"Turkey" )],"Granted",2018,2033), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Peach" )],"Granted",2019,2034), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Surry" )],"Granted",2020,2033), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"North" )],"Granted",2021,2040), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Monticello" )],"Granted",2024,2030), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Oconee" )],"Granted",2025,2034), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Summer" )],"Granted",2025,2042), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Beach" )],"Granted",2025,2033), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Browns" )],"Granted",2025,2036), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Dresden" )],"Granted",2025,2031)) +RENEWED <- RENEWED %>% as_tibble +colnames(RENEWED) <- c("Facility","Status","App_Date","Op_Date") +RENEWED <- RENEWED %>% mutate(Op_Date=as.numeric(Op_Date),App_Date=as.numeric(App_Date)) +AVG_LENGTH <- RENEWED %>% mutate(DIFF=Op_Date-App_Date) %>% pull(DIFF) %>% mean %>% round +SUBMITTED <- SUBMITTED %>% mutate(Op_Date=App_Date+AVG_LENGTH) +UPDATE <- rbind(RENEWED,SUBMITTED ) +TOTAL_ORIG <- TOTAL +TOTAL <- TOTAL %>% left_join(UPDATE) %>% mutate(CLOSE_YEAR=ifelse(Op_Date>CLOSE_YEAR & !is.na(Status),Op_Date,CLOSE_YEAR)) %>% select(-Status,-App_Date,-Op_Date) +source("./Scripts/Functions/NPV_Functions.r") +#Calculate and individual facilities rate, for all discount rates, and all years + +TOTAL_VALUE_METRICS <- MULTI_DISCOUNT_RATE_NPV(DISCOUNT_RATE_LIST,TOTAL ,DOLLARS_SAVED_PER_YEAR=CV)%>% left_join(read_csv("Data/Raw_Data/Curie_Spent_Fuel_Site_Totals.csv")) %>% select(Year,Facility,Discount,Total_Tons,Revenue) +TOTAL_VALUE_METRICS_ORIG <- MULTI_DISCOUNT_RATE_NPV(DISCOUNT_RATE_LIST,TOTAL_ORIG ,DOLLARS_SAVED_PER_YEAR=CV) %>% left_join(read_csv("Data/Raw_Data/Curie_Spent_Fuel_Site_Totals.csv")) %>% select(Year,Facility,Discount,Total_Tons,Revenue) +####Find the results for each facility size + #Main Results + +CREATE_FULL_RESULTS <- function(REVENUE_RESULTS,COST_RESULTS){return(REVENUE_RESULTS %>% left_join(COST_COST_RESULTS) %>% mutate(Profit=Revenue-Cost) %>% group_by(Discount,Capacity,Location) %>% mutate(Time_Benefit=(1-Discount)*(lead(Profit)-Profit),Op_Cost=Discount*Profit,Marginal=Time_Benefit-Op_Cost) %>% unique)} + Rev_No_Shipping <- do.call(rbind,mclapply(CAPACITY_LIST,function(CAPACITY){do.call(rbind,lapply(DISCOUNT_RATE_LIST,function(x){YEARLY_RESULTS(TOTAL_VALUE_METRICS,x,CAPACITY,0)}))},mc.cores=min(length(CAPACITY_LIST),NCORES))) + saveRDS(Rev_No_Shipping,paste0(INTERMEDIATE_DIR ,"Revenue_Estimates.Rds")) + Rev_No_Shipping <- CREATE_FULL_RESULTS(Rev_No_Shipping,COST_DATA) + saveRDS(Rev_No_Shipping,paste0(RES_DIR,"Model_Estimates.Rds")) + + #Same as main results but the optimization is based on information known in 2018, using unadjusted CURIE map data + Rev_No_Shipping_2018 <- do.call(rbind,mclapply(CAPACITY_LIST,function(CAPACITY){do.call(rbind,lapply(DISCOUNT_RATE_LIST,function(x){YEARLY_RESULTS(TOTAL_VALUE_METRICS_ORIG,x,CAPACITY,0)}))},mc.cores=min(length(CAPACITY_LIST),NCORES))) + saveRDS(Rev_No_Shipping_2018,paste0(INTERMEDIATE_DIR ,"Revenue_Estimates_2018.Rds")) + Rev_No_Shipping_2018 <- CREATE_FULL_RESULTS(Rev_No_Shipping_2018,COST_DATA) + saveRDS(Rev_No_Shipping_2018,paste0(RES_DIR,"Model_Estimates_2018.Rds")) + + #Same as main results but only the reactors with a value ready to ship are included. That is to say other reactors may be willing to pay for future rights to hold the spent fuel but only those reactors with current values ready to ship are included, which set the bounds of the starting CIFS size. This is used for demonstration. + Rev_Shipping <- do.call(rbind,mclapply(CAPACITY_LIST,function(CAPACITY){do.call(rbind,lapply(DISCOUNT_RATE_LIST,function(x){YEARLY_RESULTS(TOTAL_VALUE_METRICS,x,CAPACITY,SHIPPING_COST_PER_TON )}))},mc.cores=min(length(CAPACITY_LIST),NCORES))) + saveRDS(Rev_Shipping,paste0(INTERMEDIATE_DIR ,"Revenue_Estimates_With_Shipping_Costs.Rds")) + Rev_Shipping <- CREATE_FULL_RESULTS(Rev_Shipping,COST_DATA) + saveRDS(Rev_Shipping,paste0(RES_DIR,"Model_Estimates_With_Shipping_Costs.Rds")) + + #Same as Rev_Shipping, but only information available in 2018 from the CURIE map is used (no updated operation dates). + Rev_Shipping_2018 <- do.call(rbind,mclapply(CAPACITY_LIST,function(CAPACITY){do.call(rbind,lapply(DISCOUNT_RATE_LIST,function(x){YEARLY_RESULTS(TOTAL_VALUE_METRICS_ORIG,x,CAPACITY,SHIPPING_COST_PER_TON )}))},mc.cores=min(length(CAPACITY_LIST),NCORES))) + saveRDS(Rev_Shipping_2018,paste0(INTERMEDIATE_DIR ,"Revenue_Estimates_With_Shipping_Costs_2018.Rds")) + Rev_Shipping_2018 <- CREATE_FULL_RESULTS(Rev_Shipping_2018,COST_DATA) + saveRDS(Rev_Shipping_2018,paste0(RES_DIR,"Model_Estimates_With_Shipping_Costs_2018.Rds")) + diff --git a/Scrape_Curie_Map_Totals.js b/Scripts/Functions/Java_Scripts_To_Collect_Raw_Data/Scrape_Curie_Map_Totals.js similarity index 100% rename from Scrape_Curie_Map_Totals.js rename to Scripts/Functions/Java_Scripts_To_Collect_Raw_Data/Scrape_Curie_Map_Totals.js diff --git a/Scrape_Curie_Map_Yearly_Values.js b/Scripts/Functions/Java_Scripts_To_Collect_Raw_Data/Scrape_Curie_Map_Yearly_Values.js similarity index 100% rename from Scrape_Curie_Map_Yearly_Values.js rename to Scripts/Functions/Java_Scripts_To_Collect_Raw_Data/Scrape_Curie_Map_Yearly_Values.js diff --git a/Scripts/Functions/NPV_Functions.r b/Scripts/Functions/NPV_Functions.r new file mode 100644 index 0000000..1c30462 --- /dev/null +++ b/Scripts/Functions/NPV_Functions.r @@ -0,0 +1,44 @@ +###########Series of functions to calculate the gross consumer surplus from a CIFS for each facility, in each year from 1960 to 2082. + +#Function to find the net present revenue of a facility,given a discount rate, and the current year, and year the facility will close. This is the sum of discounted costs that WOULD have taken place if the facility was not built +VALUE_ADD <- function(r,CURRENT_YEAR,CLOSE_YEAR){ + Years_Until_Close <- max(CLOSE_YEAR-CURRENT_YEAR+1,0) + VALUES <- (1+r)^-(1:10^4) + if(Years_Until_Close==0){return(sum(VALUES))}else if(Years_Until_Close>40){return(0)} else{return(sum(VALUES[-1:-Years_Until_Close]))} +} +#A function to extend the revenues calculations to the closure date of all of the facilities. +VALUE_ADD_SINGLE_YEAR <- function(r,CURRENT_YEAR,CLOSE_YEARS){return(sapply(CLOSE_YEARS,function(x){VALUE_ADD(r,CURRENT_YEAR,x)}))} +#A function to extend the calculation of the net present revenues of each facility to all years of interest. That is what is the NPV of building the facility in each year, for each facility. +NPV_CALC <- function(r,DATA,YEARS=1960:2083){ + Facility <- pull(DATA,Facility) + RES <- cbind(Facility,do.call(cbind,lapply(YEARS,function(x){VALUE_ADD_SINGLE_YEAR(r,x,DATA$CLOSE_YEAR)}))) +colnames(RES) <- c("Facility",YEARS) + RES <- as_tibble(RES) %>% pivot_longer(cols=-Facility,names_to="Year",values_to="Revenue") %>% arrange(Year) %>% mutate(Year=parse_number(Year),Revenue=parse_number(Revenue),Discount=r) + return(RES) +} +#A function which returns a list, of net present revenue calculation tables (facility by year) for a range of possible discount rates. This allows for the results to be quickly looked up, when we want to adjust the time value of money. These results combine costs savings to calculate NPV +MULTI_DISCOUNT_RATE_NPV <- function(DISCOUNT_INCREMENT,DATA,YEARS=1960:2083,DOLLARS_SAVED_PER_YEAR){ + NCORES <- detectCores()-1 + RES <- mclapply(DISCOUNT_INCREMENT,function(r,TABLE=DATA){NPV_CALC(r,TABLE)},mc.cores = NCORES) + RES <- do.call(rbind,RES) %>% mutate(Revenue=Revenue*DOLLARS_SAVED_PER_YEAR) + return(RES) +} + + +#Function to find the maximum Revenue BUT it does not take dollars just relative savings +MAX_REV <- function(YEAR,TBL,DISCOUNT,CIFS_SIZE,SHIPPING_COST=0){ + TBL <- TBL %>% filter(Year==YEAR,Discount==DISCOUNT,Revenue/Total_Tons>SHIPPING_COST) + REV <- TBL %>% pull(Revenue) + VOL <- TBL %>% pull(Total_Tons) + RES <- lp(direction = "max", objective.in = REV, const.mat = matrix(VOL, nrow = 1),const.dir = "<=", const.rhs = CIFS_SIZE, all.bin = TRUE) + return(RES[[11]]) + } + +#Calculate all results for the start year to the end year. Discount rate, and capacities are fixed. +#DATA <- TOTAL_VALUE_METRICS;DISCOUNT<-0.05;SHIPPING_COSTS <-0 +YEARLY_RESULTS <- function(DATA,DISCOUNT,CAPACITY,SHIPPING_COSTS=0){ + RES <- cbind(1960:2083,sapply(1960:2083, function(x){MAX_REV(YEAR=x,TBL=DATA,DISCOUNT,CAPACITY,SHIPPING_COSTS)})) %>% as_tibble + colnames(RES) <- c("Year","Revenue") + RES <- RES %>% mutate(Discount=DISCOUNT,Capacity=CAPACITY) %>% select(Year,Discount,Capacity,Revenue) %>% mutate(Shipping_Cost_Cuttoff=SHIPPING_COSTS) + return(RES) +} diff --git a/Single_Price_Data_Proc.r b/Single_Price_Data_Proc.r deleted file mode 100644 index 39c4181..0000000 --- a/Single_Price_Data_Proc.r +++ /dev/null @@ -1,168 +0,0 @@ -#A script which attempts to pull in all data, and create a data frame with the maximum revenue values for each facility, year and discount rate. The output can then be used to make figures and graphs -library(tidyverse) -library(parallel) - NCORES <- detectCores()-1 -library(lpSolve) #For solving discrete value maximization for the power plants -####Manual inputs -#DISCOUNT_RATE_LIST <- seq(0,1,by=0.0025) -DISCOUNT_RATE_LIST <- c(0.03,0.05,0.07,0.1) -SHIPPING_COST_PER_TON <- 1.2874*26000 #Inflation adjusted from New Mexico Report -CV1 <- 1.3074*(10607030-1060703) #Data from Texas Report, converted from 2018 to Dec 2025 dollars -CV2 <- 1.2874*(6984013-1117442) #Data from New Mexico Report, Converted from 2019 to Dec 2025 -CV3 <- mean(CV1,CV2) #Average of the two - -###################################Cost results -#source("Cost_Data_Proc.r") -CIFS <- rbind(readRDS("Data/Cleaned_Data/Texas_CIFS_Costs.Rds"),readRDS("Data/Cleaned_Data/New_Mexico_CIFS_Costs.Rds")) -#Adjust for inflation -CIFS[,-1:-5] <- 1.2874*CIFS[,-1:-5] -COSTS <- do.call(rbind,lapply(DISCOUNT_RATE_LIST ,function(DISCOUNT){CIFS %>% group_by(Location,Capacity,Cost_Assumption) %>% mutate(NPC=Total/(1+DISCOUNT)^Year) %>% summarize(Discount=DISCOUNT,Costs=sum(NPC))})) %>% ungroup - -TEMP <- COSTS %>% group_by(Location,Cost_Assumption,Discount) %>% summarize(ST_COST=min(Costs),ST_CAPACITY=min(Capacity),M_COST=(max(Costs)-min(Costs))/(max(Capacity)-min(Capacity))) %>% ungroup -CAPACITY_INCREMENT <- 5000 -#rm(COST_DATA) -for(i in 1:nrow(TEMP)){ - LOC <- as.character(TEMP[i,"Location"]) - COST_LEVEL <- as.character(TEMP[i,"Cost_Assumption"]) - DISCOUNT <- as.numeric(TEMP[i,"Discount"]) - ST_CAP <- as.numeric(TEMP[i,"ST_CAPACITY"]) - ST_COST <- as.numeric(TEMP[i,"ST_COST"]) - COST_SLOPE <- as.numeric(TEMP[i,]$M_COST) - CAPACITY <- seq(ST_CAP,160000,by=5000) - COST <- ST_COST+COST_SLOPE*CAPACITY_INCREMENT*(0:(length(CAPACITY)-1)) - - C_RES <- cbind(CAPACITY,COST) %>% as_tibble %>% mutate(Location=LOC,Cost_Assumption=COST_LEVEL,Discount=DISCOUNT) %>% select(Location,Cost_Assumption,Discount,Capacity=CAPACITY,Cost=COST) - if(!exists("COST_DATA")){COST_DATA <- C_RES}else{COST_DATA<- rbind(COST_DATA,C_RES)} -} -saveRDS(COST_DATA ,"Data/Cleaned_Data/All_CIFS_Discounted_Costs.Rds") -#COST_DATA <- COST_DATA %>% filter(Cost_Assumption=='Average') %>% select(-Cost_Assumption) %>% unique -COST_DATA <- COST_DATA %>% filter(Cost_Assumption=='High') %>% select(-Cost_Assumption) %>% unique - - - -##All unique capacity levels that the revenues need to be calculated for -CAPACITY_LIST <- COST_DATA %>% pull(Capacity) %>% unique - -### -TOTAL <- read_csv("Data/Raw_Data/Curie_Spent_Fuel_Site_Totals.csv") %>% mutate(OP_YEAR=year(Op_Date_Min),CLOSE_YEAR=year(Close_Date_Max))%>% select(Facility,Total_Assemblies,Total_Tons,OP_YEAR,CLOSE_YEAR) -#!!!!!!!!!!!!#####TEST OF VARIABLITY IN CLOSING YEAR, CHANGE OR FORMALIZE LATER -#TOTAL <- TOTAL %>% mutate(CLOSE_YEAR=ifelse(CLOSE_YEAR>2018,CLOSE_YEAR+20,CLOSE_YEAR)) - -###########Series of functions to calculate the gross consumer surplus from a CIFS for each facility, in each year from 1960 to 2082. - -#Function to find the net present revenue of a facility,given a discount rate, and the current year, and year the facility will close. This is the sum of discounted costs that WOULD have taken place if the facility was not built -VALUE_ADD <- function(r,CURRENT_YEAR,CLOSE_YEAR){ - Years_Until_Close <- max(CLOSE_YEAR-CURRENT_YEAR+1,0) - VALUES <- (1+r)^-(1:10^4) - if(Years_Until_Close==0){return(sum(VALUES))} else{return(sum(VALUES[-1:-Years_Until_Close]))} -} -#A function to extend the revenues calculations to the closure date of all of the facilities. -VALUE_ADD_SINGLE_YEAR <- function(r,CURRENT_YEAR,CLOSE_YEARS){return(sapply(CLOSE_YEARS,function(x){VALUE_ADD(r,CURRENT_YEAR,x)}))} -#A function to extend the calculation of the net present revenues of each facility to all years of interest. That is what is the NPV of building the facility in each year, for each facility. -NPV_CALC <- function(r,DATA=TOTAL,YEARS=1960:2083){ - Facility <- pull(DATA,Facility) - RES <- cbind(Facility,do.call(cbind,lapply(YEARS,function(x){VALUE_ADD_SINGLE_YEAR(r,x,DATA$CLOSE_YEAR)}))) -colnames(RES) <- c("Facility",YEARS) - RES <- as_tibble(RES) %>% pivot_longer(cols=-Facility,names_to="Year",values_to="Revenue") %>% arrange(Year) %>% mutate(Year=parse_number(Year),Revenue=parse_number(Revenue),Discount=r) - return(RES) -} -#A function which returns a list, of net present revenue calculation tables (facility by year) for a range of possible discount rates. This allows for the results to be quickly looked up, when we want to adjust the time value of money. These results combine costs savings to calculate NPV -MULTI_DISCOUNT_RATE_NPV <- function(DISCOUNT_INCREMENT,DATA=TOTAL,YEARS=1960:2083,DOLLARS_SAVED_PER_YEAR){ - NCORES <- detectCores()-1 - RES <- mclapply(DISCOUNT_INCREMENT,NPV_CALC,mc.cores = NCORES) - RES <- do.call(rbind,RES) %>% mutate(Revenue=Revenue*DOLLARS_SAVED_PER_YEAR) - return(RES) -} - - -TOTAL_VALUE_METRICS <- MULTI_DISCOUNT_RATE_NPV(DISCOUNT_RATE_LIST ,DOLLARS_SAVED_PER_YEAR=CV2) -#TOTAL_VALUE_METRICS <- TOTAL_VALUE_METRICS %>% filter(!(Facility %in% c("Palo Verde","Vogtle")))#These facilities always have a negative NPV by the end date 2083 -TOTAL_VALUE_METRICS <- TOTAL_VALUE_METRICS %>% left_join(read_csv("Data/Raw_Data/Curie_Spent_Fuel_Site_Totals.csv")) %>% select(Year,Facility,Discount,Total_Tons,Revenue) - -#Function to find the maximum Revenue BUT it does not take dollars just relative savings -TBL <- TOTAL_VALUE_METRICS -YEAR <- 2020 -CIFS_SIZE=40000 -SHIPPING_COST <- 0 - -MAX_REV <- function(YEAR,TBL,DISCOUNT,CIFS_SIZE,SHIPPING_COST=0){ - TBL <- TBL %>% filter(Year==YEAR,Discount==DISCOUNT,Revenue/Total_Tons>SHIPPING_COST)%>% mutate(Marginal=Revenue/Total_Tons) - VALUES <- TBL %>% pull(Marginal) %>% unique - RESULT <- 0 - for(i in 1:length(VALUES)){ - VOL <- pull(TBL[TBL$Marginal<=VALUES[i],],Total_Tons) - REV <- VALUES[i]*VOL - C_MAX_VALUE <- as.numeric(lp(direction = "max", objective.in = REV, const.mat = matrix(VOL, nrow = 1),const.dir = "<=", const.rhs = CIFS_SIZE, all.bin = TRUE)[11]) - C_MAX_VALUE - RESULT <- max(RESULT,C_MAX_VALUE) - } - return(RESULT) - } - -#Calculate all results for the start year to the end year. Discount rate, and capacities are fixed. -YEARLY_RESULTS <- function(DATA,DISCOUNT,CAPACITY,SHIPPING_COSTS=0){ - RES <- cbind(1960:2083,sapply(1960:2083, function(x){MAX_REV(YEAR=x,TBL=DATA,DISCOUNT,CAPACITY,SHIPPING_COSTS)})) %>% as_tibble - colnames(RES) <- c("Year","Revenue") - RES <- RES %>% mutate(Discount=DISCOUNT,Capacity=CAPACITY) %>% select(Year,Discount,Capacity,Revenue) %>% mutate(Shipping_Cost_Cuttoff=SHIPPING_COSTS) - return(RES) -} -#Calculate and individual facilities rate, for all discount rates, and all years -TOTAL_VALUE_METRICS - Rev_No_Shipping <- do.call(rbind,mclapply(CAPACITY_LIST,function(CAPACITY){do.call(rbind,lapply(DISCOUNT_RATE_LIST,function(x){YEARLY_RESULTS(TOTAL_VALUE_METRICS,x,CAPACITY,0)}))},mc.cores=min(length(CAPACITY_LIST),NCORES))) -Rev_No_Shipping -TEST <- Rev_No_Shipping %>% left_join(COST_DATA) %>% mutate(Profit=Revenue-Cost) %>% group_by(Discount,Capacity,Location) %>% mutate(Time_Benefit=(1-Discount)*(lead(Profit)-Profit),Op_Cost=Discount*Profit,Marginal=Time_Benefit-Op_Cost) %>% unique -TEST -TEST <- TEST %>% mutate(Current_Profit=Profit/((1+Discount)^(Year-2026))) %>% ungroup -TEST <- TEST %>% group_by(Year,Location,Discount) %>% filter(Current_Profit==max(Current_Profit)) %>% ungroup -TEST %>% select(Profit,Current_Profit,Year,Discount) %>% tail -TEST %>% group_by(Location,Discount) %>% filter(Current_Profit==max(Current_Profit)) %>% print(n=100) -TEST %>% group_by(Location,Discount) %>% filter(Discount==0.04) %>% select(Year,Capacity,Profit,Current_Profit) %>% mutate(Delay_Profit=(lead(Current_Profit)-Current_Profit)/10^6) %>% print(n=70) - -ggplot(TEST %>% filter(Discount==0.03) ,aes(x=Year,y=Current_Profit/10^9,color=Capacity))+geom_line()+facet_wrap(Discount~Location,ncol=1)+scale_x_continuous(breaks=seq(1960,2083,by=5)) - -ggplot(TEST %>% filter(Discount==0.05) ,aes(x=Year,y=Current_Profit/10^9,color=Capacity))+geom_line()+facet_wrap(Discount~Location,ncol=1)+scale_x_continuous(breaks=seq(1960,2083,by=5)) - -ggplot(TEST %>% filter(Year>1970) ,aes(x=Year,y=Current_Profit/10^9,color=Capacity))+geom_line()+facet_wrap(Discount~Location,ncol=2,scales = "free_y")+scale_x_continuous(breaks=seq(1960,2083,by=5)) - - - - -ggplot(TEST ,aes(x=Year,y=Current_Profit/10^9,color=Capacity))+geom_line()+facet_wrap(Discount~Location) - - - - -TEST2 <- TEST %>% filter(Discount==0.05,Location=='Texas') - -ggplot(TEST2,aes(x=Year,y=Marginal,group=as.factor(Capacity),color=as.factor(Capacity)))+geom_line() - - #Calculate the net present revenue of all facilities that can CURRENTLY be shipped to the site with a profit. Note that the total Revenue is more important because most of the sites will be willing to pay for the right to CIFS storage in the future even if the shipping costs are too high presently. This result is used to show what might be the current ideal facility size, even if future expansion is expected to maximize profit. - Rev_Shipping <- do.call(rbind,mclapply(CAPACITY_LIST,function(CAPACITY){do.call(rbind,lapply(DISCOUNT_RATE_LIST,function(x){YEARLY_RESULTS(TOTAL_VALUE_METRICS,x,CAPACITY,SHIPPING_COST_PER_TON )}))},mc.cores=min(length(CAPACITY_LIST),NCORES))) - -Revenue_Results <- rbind(Rev_No_Shipping,Rev_Shipping) %>% mutate(Type=ifelse(Shipping_Cost_Cuttoff==0,"Revenue","Revenue_Shipping")) %>% pivot_wider(values_from=Revenue,names_from=Type) %>% select(-Shipping_Cost_Cuttoff) %>% group_by(Year,Discount,Capacity) %>% summarize(Revenue=mean(Revenue,na.rm=TRUE),Revenue_Shipping=mean(Revenue_Shipping,na.rm=TRUE)) %>% ungroup -Revenue_Results -#################################################### -COSTS <- rbind(COSTS,EXTENDED_CAPACITY_NEW_MEXICO) -COST_DATA %>% pull(Cost_Assumption) %>% unique -saveRDS(COSTS,"Data/Cleaned_Data/All_CIFS_Discounted_Costs.Rds") -COSTS %>% group_by(Location,Capac -CIFS_Data <- Revenue_Results %>% left_join(COSTS) %>% mutate(Profit=Revenue-Costs,Profit_Shipping=Revenue_Shipping-Costs) %>% select(Year,Location,Capacity,Discount,Revenue,Costs,Profit,Revenue_Shipping,Profit_Shipping,everything()) -#CIFS_Data <- -CIFS_Data <- CIFS_Data %>% group_by(Location,Capacity,Discount) %>% arrange(Location,Capacity,Discount,Year) %>% mutate(Time_Benefit=lead(Profit)-Profit,Op_Cost=Profit*Discount,Marginal=Time_Benefit-Op_Cost) %>% ungroup -TEMP <- CIFS_Data %>% filter(Discount==0.05) %>% mutate(Time_Benefit=lead(Profit)-Profit,Op_Cost=Profit*Discount,Marginal=Time_Benefit-Op_Cost) - -CIFS_Data <- CIFS_Data %>% group_by(Location,Phase,Capacity,Discount) %>% mutate(Time_Benefit=(1-Discount)*(lead(Profit)-Profit),Op_Cost=Profit*Discount,Marginal=Time_Benefit-Op_Cost) -((-3504542954)-(-3530861857))-0.05*(-3530861857 ) -STARTING_YEARS <- CIFS_Data %>% group_by(Location,Phase,Capacity,Discount) %>% mutate(PROFITABLE=Profit>0 & Marginal<=0) %>% filter(PROFITABLE) %>% filter(Year==min(Year)) %>% select(Location,Phase,Capacity,Discount,Start_Year=Year,Profit) %>% ungroup -CIFS_Data %>% mutate(Current_Profit=Profit/((1+Discount)^(Year-2026))) %>% group_by(Location,Discount,Capacity) %>% filter(Current_Profit==max(Current_Profit)) %>% select(Start_Year=Year,Location,Phase,Capacity,Current_Profit) %>% ungroup %>% arrange(Discount,Location,desc(Current_Profit)) %>% select(Discount,Location,Phase,Start_Year,Current_Profit) - -############# -ggplot(CIFS_Data ,aes(x=Year,y=Marginal/10^6,group=Capacity,color=as.factor(Capacity)))+geom_line()+facet_wrap(~Discount,ncol=1)+scale_x_continuous(breaks=seq(1960,2083,by=5)) - -ggplot(CIFS_Data ,aes(x=Year,y=Marginal,group=Capacity,color=as.factor(Capacity)))+geom_line()+facet_grid(~Discount) - - - - -dir.create("./Results",showWarnings=FALSE) -saveRDS(TOTAL_VALUE_METRICS,"./Results/Storage_Values_by_Facility_and_Variable_Discounts.Rds")