From df25529980e7348422553fd87216fbd8d5d5ba60 Mon Sep 17 00:00:00 2001 From: Alex Gebben Work Date: Mon, 23 Feb 2026 16:29:58 -0700 Subject: [PATCH] Mid way through getting two phase optimization --- Optimal_Deployments.r | 118 +++++++++++++++++++++++++++++++++++++++--- 1 file changed, 112 insertions(+), 6 deletions(-) diff --git a/Optimal_Deployments.r b/Optimal_Deployments.r index ff0111e..6dae785 100644 --- a/Optimal_Deployments.r +++ b/Optimal_Deployments.r @@ -21,9 +21,6 @@ library(lpSolve) #155462880/10^9 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) -ADDED_VOLUME -STARTING_VOLUME -VOLUME 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 @@ -52,6 +49,81 @@ if(ADDED_VOLUME/PHASE_SIZE>PHASES_ADDED){stop("Not enough capacity for the reque return(TOTAL_COST) } +################################################################ +OPTIMAL_COST2 <- function(PHASES_ADDED,VOLUME,CURRENT_YEAR,OTHER_PHASES=1,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){ +TOTAL_VOLUME_CONSTRAINT <- TRANS_CONSTRAINT*(END_YEAR-CURRENT_YEAR) + #Adjust volume constraint to include the SNF that is sent in the later phases. Find the percentage of the total possible volume that is already acounted for by the later SNF + TRANS_CONSTRAINT <- TRANS_CONSTRAINT*((TOTAL_VOLUME_CONSTRAINT-OTHER_PHASES*PHASES_ADDED)/TOTAL_VOLUME_CONSTRAINT) + 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 + UNDER_CAPACITY_VOLUME <- VOLUME-(SHIPPING_YEARS-1)*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) +} + +#### +CURRENT +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") @@ -99,12 +171,13 @@ REACTOR_VALUES <- readRDS("Data/Cleaned_Data/Reactor_Values.Rds") #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),] + 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) + 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 @@ -134,6 +207,39 @@ else { 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. -ADDITION_CHECK(CURRENT,22,1) +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 +