alex/CIFS_Storage_Demand_Model
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Major cleanup of code

Browse Source
This commit is contained in:
Alex Gebben Work 2026-02-05 11:59:59 -07:00
parent 9841ac00e7
commit aa6c892ced
10 changed files with 213 additions and 548 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
.gitignore vendored
View File

@ -1,4 +1,9 @@
# ---> R # ---> R
#Folders with large data not to be saved
Data/Cleaned_Data/
Data/Results/
#Do not save any swap files
*.swp
# History files # History files
.Rhistory .Rhistory
.Rapp.history .Rapp.history

139
Analysis.r
View File

@ -1,132 +1,27 @@
library(tidyverse) library(tidyverse)
library(scales) Rev_No_Shipping <- readRDS("Data/Cleaned_Data/Model_Estimates.Rds")
library(RcppRoll) Rev_No_Shipping_2018 <- readRDS("Data/Cleaned_Data/Model_Estimates_2018.Rds")
library(lpSolve) Rev_Shipping <- readRDS("Data/Cleaned_Data/Model_Estimates_With_Shipping_Costs.Rds")
RES <- readRDS("Results/Storage_Values_by_Facility_and_Variable_Discounts.Rds") Rev_Shipping_2018 <- readRDS("Data/Cleaned_Data/Model_Estimates_With_Shipping_Costs_2018.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")
PLOT_DATA <- OUTCOME %>% filter(Discount %in% c(0.03,0.05,0.1)) TEST <- Rev_No_Shipping
png("Revenue.png",height=6,width=11,units="in",res=900) TEST <- TEST %>% mutate(Current_Profit=Profit/((1+Discount)^(Year-2026))) %>% ungroup
ggplot(PLOT_DATA,aes(x=Year,y=Revenue,color=Discount))+geom_point() TEST <- TEST %>% group_by(Year,Location,Discount) %>% filter(Current_Profit==max(Current_Profit)) %>% ungroup
dev.off() 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) 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))
#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 %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)) 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
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)) 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_WITH_SHIPPING
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")))

251
Data_Proc.r
View File

@ -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")

2
Run_All.sh Normal file
View File

@ -0,0 +1,2 @@
Rscript ./Scripts/1_Raw_Cost_Data_Clean.r
Rscript ./Scripts/2_Compiled_Results_Data.r

19
Cost_Data_Proc.r → Scripts/1_Raw_Cost_Data_Clean.r
View File

@ -1,16 +1,21 @@
library(tidyverse) library(tidyverse)
######CIFS costs INFLATION_ADJUST <- 1.2874 #Adjust cost from 2019 to 2026
dir.create("./Data/Cleaned_Data",recursive=TRUE,showWarnings=FALSE) #################CIFS costs
read_csv("./Data/Raw_Data/Cost_Tables/Texas/Table_C-3_Undiscounted_Cost_Estimates_Phase_1_Low.csv") #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()) 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 <- 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( 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) 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")

133
Scripts/2_Compiled_Results_Data.r Normal file
View File

@ -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"))

0
Scrape_Curie_Map_Totals.js → Scripts/Functions/Java_Scripts_To_Collect_Raw_Data/Scrape_Curie_Map_Totals.js
View File

0
Scrape_Curie_Map_Yearly_Values.js → Scripts/Functions/Java_Scripts_To_Collect_Raw_Data/Scrape_Curie_Map_Yearly_Values.js
View File

44
Scripts/Functions/NPV_Functions.r Normal file
View File

@ -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)
}

168
Single_Price_Data_Proc.r
View File

@ -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")