Added non-linear results. Working to clean code

This commit is contained in:
Alex 2026-01-30 17:41:40 -07:00
parent 2d2de480df
commit d55898ce51
2 changed files with 88 additions and 30 deletions

View File

@ -1,6 +1,7 @@
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
@ -9,36 +10,52 @@ CV1 <- 1.3074*(10607030-1060703) #Data from Texas Report, converted from 2018 to
CV2 <- 1.2874*(6984013-1117442) #Data from New Mexico Report, Converted from 2019 to Dec 2025
LENGTH <- length(RES)
SUPPLY_RESULTS <- list()
SUPPLY_RESULTS_WITH_SHIPPING <- list() #This shows the amount availble to ship today, the SUPPLY_RESULTS shows the demand to have a CIFS available sometime in the future (present value brought forward with a discount). This SUPPLY_RESULTS_WITH_SHIPPING shows only SNF that has a current Net Present Cost high enough to exceed the shipping costs.
SHIPPING_COST <- 1.2874*26000 #Inflation adjusted from New Mexico Report
TEST <- RES$Per_5%>% left_join(read_csv("Data/Raw_Data/Curie_Spent_Fuel_Site_Totals.csv")) %>% mutate(Marginal_Value=Revenue/Total_Tons) %>% select(Year,Facility,Marginal_Value,Total_Tons) %>% group_by(Year) %>% arrange(Year,desc(Marginal_Value)) %>% mutate(Q=cumsum(Total_Tons)) %>% mutate(IN=Q<10000) %>% mutate(MIN=ZZMIN
TEST %>% filter(IN) %>% filter(Marginal_Value==min(Marginal_Value)) %>% print(n=100)
TEST %>% print(n=30)
TEST %>% arrange(Year,Marginal_Value)
TEST
TEST %>% filter(IN) %>% filter(Marginal_Value )
TEST <- RES[[1]]%>% left_join(read_csv("Data/Raw_Data/Curie_Spent_Fuel_Site_Totals.csv")) %>% select(Facility,Year,Revenue,Total_Tons) %>% mutate(Marginal=Revenue/Total_Tons) %>% group_by(Year) %>% arrange(Marginal) %>% mutate(IN=cumsum(Total_Tons)<10000)
TEST
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)
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
SUPPLY_RESULTS[[i]] <- RES[[i]] %>% left_join(read_csv("Data/Raw_Data/Curie_Spent_Fuel_Site_Totals.csv")) %>% mutate(Marginal_Value=Revenue/Total_Tons) %>% select(Year,Marginal_Value,Total_Tons) %>% group_by(Year) %>% arrange(Year,desc(Marginal_Value)) %>% mutate(Q=cumsum(Total_Tons),Discount=factor(percent(DISCOUNT,1),levels=c("3%","5%","7%","10%")))
SUPPLY_RESULTS_WITH_SHIPPING[[i]] <- RES[[i]] %>% left_join(read_csv("Data/Raw_Data/Curie_Spent_Fuel_Site_Totals.csv")) %>% mutate(Marginal_Value=Revenue/Total_Tons) %>% filter(Marginal_Value>=SHIPPING_COST) %>% select(Year,Marginal_Value,Total_Tons) %>% group_by(Year) %>% arrange(Year,desc(Marginal_Value)) %>% mutate(Q=cumsum(Total_Tons),Discount=factor(percent(DISCOUNT,1),levels=c("3%","5%","7%","10%")))
RES[[i]] <- RES[[i]] %>% left_join(read_csv("Data/Raw_Data/Curie_Spent_Fuel_Site_Totals.csv")) %>% select(Year,Facility,Total_Tons,Revenue)
}
names(SUPPLY_RESULTS) <- names(RES)
names(SUPPLY_RESULTS_WITH_SHIPPING) <- names(RES)
KEY_DATA <- do.call(rbind,SUPPLY_RESULTS[c("Per_3","Per_5","Per_7","Per_10")]) %>% ungroup
SUPPLY_RESULTS_WITH_SHIPPING
KEY_DATA_WITH_SHIPPING <- do.call(rbind,SUPPLY_RESULTS_WITH_SHIPPING[c("Per_3","Per_5","Per_7","Per_10")]) %>% ungroup
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

View File

@ -1,14 +1,18 @@
#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)
library(lpSolve) #For solving discrete value maximization for the power plants
TS <- read_csv("Data/Raw_Data//Curie_Spent_Fuel_Timeline.csv")
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)
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
###########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^6)
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.
@ -18,17 +22,54 @@ 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))
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 will need to be combined with costs to calculate NPV
MULTI_DISCOUNT_RATE_NPV <- function(INCREMENT=0.005,DATA=TOTAL,YEARS=1960:2083){
#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(INCREMENT=0.005,DATA=TOTAL,YEARS=1960:2083,DOLLARS_SAVED_PER_YEAR){
NCORES <- detectCores()-1
RES <- mclapply(seq(0,1,by=INCREMENT),NPV_CALC,mc.cores = NCORES)
names(RES) <- paste0("Per_",100*seq(0,1,by=INCREMENT)) #Name with the discount rate of the given table
RES <- do.call(rbind,RES) %>% mutate(Revenue=Revenue*DOLLARS_SAVED_PER_YEAR)
return(RES)
}
TOTAL_VALUE_METRICS <- MULTI_DISCOUNT_RATE_NPV(INCREMENT=0.0025)
TOTAL_VALUE_METRICS <- MULTI_DISCOUNT_RATE_NPV(INCREMENT=0.1,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)
#TBL <- TOTAL_VALUE_METRICS
#DICOUNT=0.1
#CIFS_SIZE=1
#SHIPPING_COST <- 0
#YEAR=2026
#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)
VOL
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]])
}
SHIPPING_COST_PER_TON <- 1.2874*26000 #Inflation adjusted from New Mexico Report
#Unique values of discount rates used
DISCOUNTS <- TOTAL_VALUE_METRICS$Discount %>% unique
#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)
return(RES)
}
#Calculate and individual facilities rate, for all discount rates, and all years
FACILITY_RESULTS <- function(CAPACITY){do.call(rbind,lapply(DISCOUNTS,function(x){YEARLY_RESULTS(TOTAL_VALUE_METRICS,x,CAPACITY)}))}
lapply(FACILITY_RESULTS()
FACILITY_RESULTS(10^5)
MAX_REV(YEAR=2020,TBL=TOTAL_VALUE_METRICS,DISCOUNT=0.1,CIFS_SIZE=100000,SHIPPING_COST_PER_TON )
dir.create("./Results",showWarnings=FALSE)
saveRDS(TOTAL_VALUE_METRICS,"./Results/Storage_Values_by_Facility_and_Variable_Discounts.Rds")