From 520f79ddd7a27422b8f6e467e6987ab226d5965d Mon Sep 17 00:00:00 2001 From: Alex Gebben Work Date: Wed, 25 Feb 2026 17:00:54 -0700 Subject: [PATCH] Testing new function for cost --- Update.r | 237 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 237 insertions(+) create mode 100644 Update.r diff --git a/Update.r b/Update.r new file mode 100644 index 0000000..dce9616 --- /dev/null +++ b/Update.r @@ -0,0 +1,237 @@ +library(tidyverse) +library(lpSolve) +##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" + +######### +NUM_PHASES1 <- 1 +NUM_PHASES2 <- 5 +STAGE1_YEAR <-10 +STAGE2_YEAR <- 20 +ST_VOLUME=8680 +VOLUME1<- 5000 +VOLUME2 <-1000 +DISCOUNT <- 0.05 +STARTING_CAPACITY=8680 +PHASE_SIZE <- 5000 +END_TIME <- 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 + + +################################################################ +OPTIMAL_COST <- function(NUM_PHASES1,NUM_PHASES2,STAGE1_YEAR,STAGE2_YEAR,ST_VOLUME,VOLUME1,VOLUME2,DISCOUNT=0.05,STARTING_CAPACITY=8680,PHASE_SIZE=5000,END_TIME=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){ + #Total volume not being analyzed, Other capacity is the capacity known in the future no in these two phases + + VOLUME <- ST_VOLUME+VOLUME1+VOLUME2 + OTHER_CAPACITY <- NUM_PHASES2*PHASE_SIZE+STARTING_CAPACITY + TIME_BETWEEN_PHASES <- STAGE2_YEAR-STAGE1_YEAR + TIME_LEFT <- END_TIME-STAGE1_YEAR +####Check if there is enough capacity to store the requested SNF + TOTAL_CAPACITY <- STARTING_CAPACITY+PHASE_SIZE*(NUM_PHASES1+NUM_PHASES2) +if(TOTAL_CAPACITYST_VOLUME){stop("Not enough capacity for the requested volume in the starting phase")} +if(NUM_PHASES1*PHASE_SIZETIME_BETWEEN_PHASES){stop("Not enough time to ship SNF to the CIFS")} +####Check if there is enough time to ship all SNF for disposal. + SHIPPING_YEARS <- ceiling(VOLUME/TRANS_CONSTRAINT) + if(TIME_LEFT% filter(Discount==0.05,Year==2046) +MAX_REV2 <- function(TBL,TBL2,CIFS_SIZE,LATER_PHASE_SIZE,LATER_PHASE_YEARS_AHEAD,DISCOUNT){ + REV1 <- TBL %>% pull(Revenue) + REV2 <- TBL2 %>% pull(Revenue) + REV2 <- REV2/((1+DISCOUNT)^LATER_PHASE_YEARS_AHEAD) + REV_ALL <- c(REV1,REV2) + VOL <- TBL %>% pull(Total_Tons) +VOL_CONST <- c(VOL,rep(0,length(VOL))) +ONLY_ONE_CONST <- cbind(diag(length(VOL)),diag(length(VOL))) +CONSTRAINTS <- rbind(VOL_CONST,rev(VOL_CONST),ONLY_ONE_CONST) +MAX_SIZE <- c(CIFS_SIZE,LATER_PHASE_SIZE,rep(1,length(VOL))) + RES <- lp(direction = "max", objective.in = REV_ALL, const.mat =CONSTRAINTS,const.dir = "<=", const.rhs = MAX_SIZE, all.bin = TRUE) +#TBL[RES$solution[1:(length(RES$solution)/2)]==1,] +#RES$solution[(length(RES$solution)/2+1):length(RES$solution)] + return(RES) + } +PROFIT_EST2 <- function(ADDED_PHASES,ST_YEAR,YEARS_AHEAD,NEXT_PHASE_YEARS_AHEAD,DATA=REACTOR_VALUES,DISCOUNT_RATE=0.05,ST_CAP=8680,PHASE_SIZE=5000,NEXT_PHASE_AHEAD_NUM=1){ + CURRENT <- DATA%>% filter(Year==ST_YEAR+YEARS_AHEAD,Discount==DISCOUNT_RATE) + LATER <- DATA%>% filter(Year==ST_YEAR+NEXT_PHASE_YEARS_AHEAD,Discount==DISCOUNT_RATE) + + RES <- MAX_REV2(CURRENT,LATER,ST_CAP+PHASE_SIZE*ADDED_PHASES,PHASE_SIZE*NEXT_PHASE_AHEAD_NUM,NEXT_PHASE_YEARS_AHEAD,DISCOUNT_RATE) + return(RES) +} +CURRENT <- REACTOR_VALUES %>% filter(Year==2042,Discount==0.05) +LATER <- REACTOR_VALUES %>% filter(Year==2043,Discount==0.05) +function(TBL,TBL2,CIFS_SIZE,LATER_PHASE_SIZE,LATER_PHASE_YEARS_AHEAD,DISCOUNT){ +MAX_REV2(CURRENT,LATER,8680+10000,5000,1,0.05) + +##################### +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) +} + + +################################ +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) +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) + } + +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) +} + +REACTOR_VALUES <- readRDS("Data/Cleaned_Data/Reactor_Values.Rds") +#ADDED <- 1 +#ADD_RES <- MAX_REV(CURRENT,8680+5000*ADDED) +#REACTOR_DATA <- CURRENT +#STARTING_YEAR <- 2026 +#YEARS_AHEAD=20 +#STARTING_CAP=8680 +#SINGLE_PHASE_CAP=5000 +#ADDED_UNITS <- 1 +#Find the optimal profit and cost, plus if the capacity constraint of an addtion in binding. +ADDITION_CHECK <- function(REACTOR_DATA,ADDED_UNITS,YEARS_AHEAD,Discount_Rate=0.05,STARTING_CAP=8680,SINGLE_PHASE_CAP=5000){ + OPTIM_GUESS <- MAX_REV(REACTOR_DATA,STARTING_CAP+SINGLE_PHASE_CAP*ADDED_UNITS) + ADD_RES <- MAX_REV(CURRENT,STARTING_CAP+5000*ADDED_UNITS) + SELECTED_REACTORS <- REACTOR_DATA[which(ADD_RES$solution==1),]%>% mutate(MARGINAL_VALUE=Revenue/Total_Tons) + FOUND_VOLUME <- OPTIM_GUESS$constraint[76] + LOWEST_VALUE_REACTOR <- SELECTED_REACTORS[SELECTED_REACTORS$MARGINAL_VALUE== min(SELECTED_REACTORS$MARGINAL_VALUE),] + MARGINAL_VALUE <- LOWEST_VALUE_REACTOR$MARGINAL_VALUE + FULL_COST_AT_CAPACITY <- OPTIMAL_COST(ADDED_UNITS,FOUND_VOLUME,YEARS_AHEAD) + MARGINAL_COST <- FULL_COST_AT_CAPACITY - OPTIMAL_COST(ADDED_UNITS,FOUND_VOLUME-1,YEARS_AHEAD) + BOUNDED <- MARGINAL_VALUE>MARGINAL_COST +if(BOUNDED){ + FOUND_VOLUME <- FOUND_VOLUME-STARTING_CAP + State <- "Capacity Constrainted" + Optimal_Rev <- OPTIM_GUESS$objval + Optimal_Cost <- FULL_COST_AT_CAPACITY +} +else { + HIGHEST_VALUE_REACTOR <-SELECTED_REACTORS[SELECTED_REACTORS$MARGINAL_VALUE== max(SELECTED_REACTORS$MARGINAL_VALUE),] + MARGINAL_VALUE <- unique(HIGHEST_VALUE_REACTOR$MARGINAL_VALUE) + if(MARGINAL_VALUE>=MARGINAL_COST){ + State <- "Not a Profitable Phase" + FOUND_VOLUME <- 0 + OPTIM_STARTING <- MAX_REV(REACTOR_DATA,STARTING_CAP) + Optimal_Rev <- OPTIM_STARTING$objval + Optimal_Cost <- OPTIMAL_COST(0,STARTING_CAP,YEARS_AHEAD) + }else{ + State <- "No Binding Constraints" + FOUND_VOLUME <- NA + Optimal_Rev <- NA + Optimal_Cost <- NA + + } + +} + Profit <- Optimal_Rev-Optimal_Cost + return(c(ADDED_UNITS,State,FOUND_VOLUME,Profit,Optimal_Rev,Optimal_Cost)) +} +#Note for self: By running addition's from 1 to 20 (Roughly) at the same number of years ahead the number of 5000 unit addtions which maximizes profit in that year can be found. It looks like at least 22 units will be built which is enough for the whole US, but the timing of addtions needs to be worked out by backwards induction using the years. +ST_YEAR <- 2026 +CONSTRAINED <- matrix(NA,nrow=30,ncol=40) +PROFIT <- matrix(NA,nrow=30,ncol=40) + +for(i in ST_YEAR:(ST_YEAR+40)){ + COL <- i-ST_YEAR + CURRENT <- REACTOR_VALUES %>% filter(Year==i,Discount==0.05) +for(n in 1:30){ +C_RES <- try(ADDITION_CHECK(CURRENT,n,COL)) + + CONSTRAINED[n,COL] <- C_RES[2] + PROFIT[n,COL] <- C_RES[6] + + } +} +CONSTRAINED <- CONSTRAINED %>% as_tibble +PROFIT <- PROFIT %>% as_tibble +colnames(CONSTRAINED) <- ST_YEAR:(ST_YEAR+40) +colnames(PROFIT) <- ST_YEAR:(ST_YEAR+40) + +CONSTRAINED$Size <- 1:30 +PROFIT$Size <- 1:30 + + +CONSTRAINED <- CONSTRAINED %>% select(Size,everything()) %>% pivot_longer(-Size,names_to="Year",values_to="Status") +PROFIT <- PROFIT%>% select(Size,everything()) %>% pivot_longer(-Size,names_to="Year",values_to="Profit") +PROFIT_RES <- PROFIT %>% filter(!is.na(Profit)) %>% group_by(Year) %>% mutate(Profit=as.numeric(Profit)) %>% filter(Profit==max(Profit)) %>% print(n=100) +CONSTRAINED$Status <- ifelse(grepl("Not enough time to",CONSTRAINED$Status),"Time Limited",CONSTRAINED$Status) + +CONSTRAINED$Status %>% unique +MAX_VOLUME <- CONSTRAINED %>% group_by(Year,Status) %>% summarize(Size=max(Size)) %>% filter(Status=="Capacity Constrainted") +ggplot(MAX_VOLUME,aes(x=Year,y=Size))+geom_point()+scale_y_continuous(breaks=1:30) +ggplot(PROFIT_RES,aes(x=Year,y=Size))+geom_point()+scale_y_continuous(breaks=1:30) + +9*5000+8680 +