diff --git a/Analysis.r b/Analysis.r index c07bd4f..947eb89 100644 --- a/Analysis.r +++ b/Analysis.r @@ -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)) diff --git a/Cost_Estimates.r b/Cost_Estimates.r new file mode 100644 index 0000000..ff66fe4 --- /dev/null +++ b/Cost_Estimates.r @@ -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]],] + + + diff --git a/Scripts/2_Compiled_Results_Data.r b/Scripts/2_Compiled_Results_Data.r index e757150..9e15094 100644 --- a/Scripts/2_Compiled_Results_Data.r +++ b/Scripts/2_Compiled_Results_Data.r @@ -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))) diff --git a/Scripts/5_All_Stages.r b/Scripts/5_All_Stages.r index 26c9b9e..b9b66c2 100644 --- a/Scripts/5_All_Stages.r +++ b/Scripts/5_All_Stages.r @@ -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")