133 lines
11 KiB
R
133 lines
11 KiB
R
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")
|
|
|
|
|
|
|
|
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()
|
|
|
|
|
|
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"))
|
|
|
|
|
|
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
|
|
|
|
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
|
|
|
|
|
|
|
|
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")))
|