Working on new methods using costs supplied

This commit is contained in:
Alex Gebben Work 2026-02-19 17:06:34 -07:00
parent 83061c02c8
commit 7457ab381d
4 changed files with 118 additions and 7 deletions

View File

@ -1,11 +1,12 @@
library(tidyverse)
#Rev_No_Shipping <- readRDS("Data/Results/Model_Estimates.Rds")
Rev_No_Shipping <- readRDS("Data/Results/Model_Estimates.Rds")
#Rev_No_Shipping %>% pull(Discount) %>% unique
#Rev_No_Shipping <- Rev_No_Shipping %>% mutate(Discount=round(Discount,5))
#Rev_No_Shipping_2018 <- readRDS("Data/Results/Model_Estimates_2018.Rds")
Rev_Shipping <- readRDS("Data/Results/Model_Estimates_With_Shipping_Costs.Rds")
Rev_Shipping
#Rev_Shipping_2018 <- readRDS("Data/Results/Model_Estimates_With_Shipping_Costs_2018.Rds")
#Rev_No_Shipping_2018 <- readRDS("Data/Results/Model_Estimates_2018.Rds")
@ -14,9 +15,9 @@ Rev_Shipping <- readRDS("Data/Results/Model_Estimates_With_Shipping_Costs.Rds")
#TEST <- Rev_No_Shipping
TEST <- Rev_No_Shipping
#TEST <- Rev_No_Shipping_2018
TEST <- Rev_Old_Sites
#TEST <- Rev_Old_Sites
TEST <- TEST %>% mutate(Discount=round(Discount,5))

103
Cost_Estimates.r Normal file
View File

@ -0,0 +1,103 @@
library(tidyverse)
##Value from table C-2 of Holtec report, cost categories
#"Initial Annual Estimated Costs and 2019 Constant Dollar Values for the Various Activities for the Proposed CISF and the No-Action Alternative"
#########
#CURRENT_YEAR <- 30
#END_YEAR <- 40
#PHASE_SIZE <- 5000
#ST_CAP <- 8680
#TRANSPORT_COST_RATIO <- 155462880/5000 # or 269883561/8680 since transportation cost scale linearly
#VOLUME <- 10^5/2+480
#TRANS_CONSTRAINT <- (ST_CAP+5000*19)/19
#PHASES_ADDED <- max(ceiling((VOLUME-ST_CAP)/PHASE_SIZE),0) #Make sure added phases are not negative
OPTIMAL_COST <- function(PHASES_ADDED,VOLUME,CURRENT_YEAR,DISCOUNT=0.05,STARTING_VOLUME=8680,PHASE_SIZE=5000,END_YEAR=40,PHASE_CONST_COST=103399272,TRANSPORT_COST_RATIO=155462880/5000,TRANS_CONSTRAINT =(8680+5000*19)/19,DECOM_COST_PER_TON =24822656/5000,INFLATION_ADJUST=1.2874){
ADDED_VOLUME=max(VOLUME-STARTING_VOLUME,0)
if(ADDED_VOLUME/PHASE_SIZE>PHASES_ADDED){stop("Not enough capacity for the requested volume")}
###Construction Cost to cover supplied SNF volume
CONSTRUCT_COST <- PHASES_ADDED*PHASE_CONST_COST
#Cost to transport SNF from reactors to the CIFS in the current year and to repository at end of project period
SHIPPING_TIME<- VOLUME/TRANS_CONSTRAINT
SHIPPING_YEARS <- ceiling(SHIPPING_TIME )
if(SHIPPING_YEARS>END_YEAR-CURRENT_YEAR){stop("Not enough time to ship the requested SNF volume")}
AT_CAPACITY_VOLUME <- TRANS_CONSTRAINT*floor(SHIPPING_TIME )
UNDER_CAPACITY_VOLUME <- VOLUME-AT_CAPACITY_VOLUME
SHIPPING_SCHEDULE_OUT <- rep(AT_CAPACITY_VOLUME,SHIPPING_YEARS)
if(UNDER_CAPACITY_VOLUME!=0){SHIPPING_SCHEDULE_OUT[1] <- UNDER_CAPACITY_VOLUME}
SHIPPING_SCHEDULE_OUT <- TRANSPORT_COST_RATIO*SHIPPING_SCHEDULE_OUT
SHIPPING_SCHEDULE_IN <- rev(SHIPPING_SCHEDULE_OUT)
SHIPPING_DISCOUNT <- (1/((1+DISCOUNT)^(1:SHIPPING_YEARS)))
SHIPPING_IN_TOTAL_COST <- sum(SHIPPING_SCHEDULE_IN*SHIPPING_DISCOUNT)
OUT_DISCOUNT <- 1/((1+DISCOUNT)^(END_YEAR-CURRENT_YEAR-SHIPPING_YEARS))
SHIPPING_OUT_TOTAL_COST <- OUT_DISCOUNT*sum(SHIPPING_SCHEDULE_OUT*SHIPPING_DISCOUNT)
SHIPPING_IN_TOTAL_COST +SHIPPING_OUT_TOTAL_COST +CONSTRUCT_COST
##Decommsioning cost: See Holtec Report, the NRC applies a adjustment factor for the larger capacity which we apply also to the small capacity this is in Section C-2 (just after the table and above section C-3)
YEARS_UNTIL_DECOM <- END_YEAR-CURRENT_YEAR+1
#Assume decom is paid by phase not average volume
DECOM_TOTAL_COST <-(DECOM_COST_PER_TON*PHASES_ADDED*PHASE_SIZE )/((1+DISCOUNT)^YEARS_UNTIL_DECOM)
TOTAL_COST <- CONSTRUCT_COST+SHIPPING_IN_TOTAL_COST+SHIPPING_OUT_TOTAL_COST+ DECOM_TOTAL_COST
TOTAL_COST <- TOTAL_COST*INFLATION_ADJUST
return(TOTAL_COST)
}
OPTIMAL_COST(20,10^4,20)
CHECK_FEASIBLE_SHIPPING <- function(VOLUME,ST_YEAR){
RESULT <- try(OPTIMAL_COST(10^5,VOLUME,ST_YEAR),silent=TRUE)
return(class(RESULT)!="try-error")
}
FIND_FEASIBLE_LIMIT <- function(STARTING_TIME){
ST_BOUND <-0
END_BOUND <-140*10^3
MID <- ceiling((END_BOUND+ST_BOUND)/2)
while(END_BOUND-ST_BOUND>10){
CHECK <-CHECK_FEASIBLE_SHIPPING(MID,STARTING_TIME)
if(CHECK){ST_BOUND <- MID} else{END_BOUND <- MID}
MID <- ceiling((END_BOUND+ST_BOUND)/2)
}
MAX_CAPACITY <- ST_BOUND+max(which(sapply(ST_BOUND:END_BOUND,CHECK_FEASIBLE_SHIPPING,ST_YEAR=STARTING_TIME)))
return(MAX_CAPACITY)
}
SHIPPING_CAPACITY_LIMITS <- cbind(1:40,sapply(1:40,FIND_FEASIBLE_LIMIT)) %>% as_tibble %>% rename(Year=V1,Max_Capacity=V2)
SHIPPING_CAPACITY_LIMITS %>% print(n=100)
REACTOR_VALUES <- readRDS("Data/Cleaned_Data/Reactor_Values.Rds")
C_VALUES <- REACTOR_VALUES %>% filter(Year==2026+40,Discount==0.05)
C_VALUES
TBL <- CURRENT
CIFS_SIZE <- ST_CAP
MAX_REV <- function(TBL,CIFS_SIZE){
# 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,byrow=TRUE),const.dir = "<=", const.rhs = CIFS_SIZE, all.bin = TRUE)
return(RES)
}
ST_CAP <- 8680
PHASE_SIZE <- 5000
RES <-
names(RES)
PROFIT_EST <- function(ADDED_PHASES,ST_YEAR,YEARS_AHEAD,DATA=REACTOR_VALUES,DISCOUNT_RATE=0.05,ST_CAP=8680,PHASE_SIZE=5000){
CURRENT <- DATA%>% filter(Year==ST_YEAR+YEARS_AHEAD,Discount==DISCOUNT_RATE)
RES <- MAX_REV(CURRENT,ST_CAP+PHASE_SIZE*ADDED_PHASES)
REVENUE <- RES$objval
TONS_STORED <- sum(CURRENT$Total_Tons*RES$solution)
PROFIT <- rbind(REVENUE,OPTIMAL_COST(ADDED_PHASES,TONS_STORED,YEARS_AHEAD))
return(PROFIT)
}
MAX_REV(CURRENT,ST_CAP+PHASE_SIZE*1)
OPTIMAL_COST(0,TONS_STORED,YEARS_AHEAD)
t(sapply(0:3,PROFIT_EST,ST_YEAR=2026,YEARS_AHEAD=15)/10^6) %>% as_tibble %>% rename(Rev=V1,Cost=V2) %>% mutate(Profit=Rev-Cost)
RES$solution
names(RES)
length(RES$objective)
CURRENT[RES[[9]],]
CURRENT[-RES[[9]],]
REACTOR_VALUES %>% filter(Year=[RES[[9]],]

View File

@ -105,7 +105,8 @@ source("./Scripts/Functions/NPV_Functions.r")
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
saveRDS(TOTAL_VALUE_METRICS ,"Data/Cleaned_Data/Reactor_Values.Rds")
####Find the results for each facility size
#Main Results
CREATE_FULL_RESULTS <- function(REVENUE_RESULTS,COST_RESULTS){return(REVENUE_RESULTS %>% left_join(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)))

View File

@ -6,7 +6,12 @@ library(ompr.roi)
library(ROI)
library(ROI.plugin.glpk)
source("./Scripts/Data_Setup.r")
CURRENT_YEAR <- 2026
TOTAL
r <- 0.05
TEMP <- TOTAL %>% mutate(TIME_VALUE=CV/((1+r)^(CLOSE_YEAR-CURRENT_YEAR))/Total_Tons,Marginal=CV/Total_Tons) %>% mutate(ORDERED=rank(-Marginal),ORDERED_TIME=rank(-TIME_VALUE) ) %>% arrange(desc(Marginal)) %>% mutate(MARGINAL_TONS=cumsum(Total_Tons),MARGINAL_INCLUDE=ifelse(MARGINAL_TONS<=10^5+5000,1,0)) %>% arrange(ORDERED_TIME) %>% mutate(TIME_TONS=cumsum(Total_Tons),TIME_INCLUDE=ifelse(TIME_TONS<=10^5+5000,1,0))
TEMP %>% filter(MARGINAL_INCLUDE==0, TIME_INCLUDE==0)
factorial(30)
# Solve fixed-charge multi-knapsack (Generalized Assignment with knapsack costs)
# profit: n x m matrix, profit[i, k] = value of assigning item i to knapsack k
# weights: length-n vector (default) OR n x m matrix for knapsack-specific weights
@ -185,12 +190,13 @@ profit <- sapply(1:30,CURRENT_VALUE,YEARS_UNTIL_CLOSE_LIST=YEARS_UNTIL_CLOSE)
# profit <- DISCOUNT_MAT*PROFIT_MAT
res <- solve_fixed_charge_multi_knapsack( profit = profit, weights = weights, caps = caps, costs = costs, force_k1_cost = TRUE # pay K1 cost even if unused
)
gc()
return(res )
}
YEARS <- 1970:2056
YEARS <- 1970:2050
#ADJUST_FACTOR <- 1/((1+0.05)^(YEARS-2026))
TEST <- mclapply(YEARS,OPTIMIZE_CIFS,mc.cores=NCORES)
saveRDS(TEST,"../GRADATED_CAPACITY_TEST_RESULTS.Rds")
RES <- cbind(YEARS,(sapply(1:length(TEST),function(x){TEST[[x]]$total_profit})/10^6)) %>% as_tibble
colnames(RES) <- c("Year","Profit")
plot(x=RES$Year,y=RES$Profit)
saveRDS(TEST,"../GRADATED_CAPACITY_TEST_RESULTS.Rds")