alex/CIFS_Storage_Demand_Model
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Testing new function for cost

Browse Source
This commit is contained in:
Alex Gebben Work 2026-02-25 17:00:54 -07:00
parent df25529980
commit 520f79ddd7
1 changed files with 237 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

237
Update.r Normal file
View File

@ -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_CAPACITY<VOLUME){stop("Not enough capacity for the requested volume")}
if(STARTING_CAPACITY>ST_VOLUME){stop("Not enough capacity for the requested volume in the starting phase")}
if(NUM_PHASES1*PHASE_SIZE<VOLUME1){stop("Not enough capacity for the requested volume in phase 1")}
if(NUM_PHASES2*PHASE_SIZE<VOLUME2){stop("Not enough capacity for the requested volume in phase 2")}
####Check if there is enough time to ship SNF into phase 1 before phase 2 starts
if(ceiling(VOLUME1/TRANS_CONSTRAINT)>TIME_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<SHIPPING_YEARS){stop("Not enough time to ship the requested SNF volume for disposal")}
###Construction Cost to cover supplied SNF volume
CONSTRUCT_COST <- NUM_PHASES1*PHASE_CONST_COST+NUM_PHASES2*PHASE_CONST_COST/((1+DISCOUNT)^TIME_BETWEEN_PHASES)
#Cost to transport SNF from reactors to the CIFS in the current year and to repository at end of project period
UNDER_CAPACITY_VOLUME <- VOLUME-(SHIPPING_YEARS-1)*TRANS_CONSTRAINT
SHIPPING_SCHEDULE_OUT <- rep(TRANS_CONSTRAINT,SHIPPING_YEARS)
if(UNDER_CAPACITY_VOLUME!=0){SHIPPING_SCHEDULE_OUT[1] <- UNDER_CAPACITY_VOLUME}
SHIPPING_SCHEDULE_OUT <- TRANSPORT_COST_RATIO*SHIPPING_SCHEDULE_OUT
SHIPPING_OUT_COST <- sum(SHIPPING_SCHEDULE_OUT*(1/((1+DISCOUNT)^(1:SHIPPING_YEARS))))*(1/((1+DISCOUNT)^(END_TIME-STAGE1_YEAR)))
SHIPPING_YEARS1 <- ceiling(VOLUME1/TRANS_CONSTRAINT)
SHIPPING_YEARS2 <- ceiling(VOLUME2/TRANS_CONSTRAINT)
UNDER_CAPACITY_VOLUME_PHASE1 <- VOLUME1-(SHIPPING_YEARS1-1)*TRANS_CONSTRAINT
UNDER_CAPACITY_VOLUME_PHASE2 <- VOLUME2-(SHIPPING_YEARS1-2)*TRANS_CONSTRAINT
PHASE_ONE_IN <- rep(TRANS_CONSTRAINT,SHIPPING_YEARS1)
PHASE_TWO_IN <- rep(TRANS_CONSTRAINT,SHIPPING_YEARS2)
UNDER_CAPACITY_VOLUME_PHASE1
if(UNDER_CAPACITY_VOLUME_PHASE1!=0){PHASE_ONE_IN[SHIPPING_YEARS1] <- UNDER_CAPACITY_VOLUME_PHASE1}
if(UNDER_CAPACITY_VOLUME_PHASE2 !=0){PHASE_TWO_IN[SHIPPING_YEARS2] <- UNDER_CAPACITY_VOLUME_PHASE2 }
SHIPPING_COST_IN1 <- sum(TRANSPORT_COST_RATIO*PHASE_ONE_IN*(1/((1+DISCOUNT)^(1:SHIPPING_YEARS1))))
SHIPPING_COST_IN2 <- sum((TRANSPORT_COST_RATIO*PHASE_TWO_IN*(1/((1+DISCOUNT)^(1:SHIPPING_YEARS2))))/((1+DISCOUNT)^(TIME_BETWEEN_PHASES )))
# SHIPPING_COST_IN1+SHIPPING_COST_IN2+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_TIME-STAGE1_YEAR+1
#Assume decom is paid by phase not average volume
DECOM_TOTAL_COST <-(DECOM_COST_PER_TON*(NUM_PHASES1+NUM_PHASES2)*PHASE_SIZE )/((1+DISCOUNT)^YEARS_UNTIL_DECOM)
TOTAL_COST <- SHIPPING_COST_IN1+SHIPPING_COST_IN2+CONSTRUCT_COST + DECOM_TOTAL_COST
TOTAL_COST <- TOTAL_COST*INFLATION_ADJUST
return(TOTAL_COST)
}
####
TBL2 <- REACTOR_VALUES %>% 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