From 83061c02c8524af89d393e56a4ec7e2de6d90aec Mon Sep 17 00:00:00 2001 From: Alex Gebben Work Date: Tue, 10 Feb 2026 15:18:34 -0700 Subject: [PATCH] Working to add CIFS stages to optimization --- Analysis.r | 25 ++++-- Scripts/5_All_Stages.r | 196 +++++++++++++++++++++++++++++++++++++++++ Scripts/Data_Setup.r | 102 +++++++++++++++++++++ 3 files changed, 316 insertions(+), 7 deletions(-) create mode 100644 Scripts/5_All_Stages.r create mode 100644 Scripts/Data_Setup.r diff --git a/Analysis.r b/Analysis.r index ca5e0bc..c07bd4f 100644 --- a/Analysis.r +++ b/Analysis.r @@ -1,24 +1,35 @@ library(tidyverse) -Rev_No_Shipping <- readRDS("Data/Results/Model_Estimates.Rds") -Rev_No_Shipping_2018 <- readRDS("Data/Results/Model_Estimates_2018.Rds") -#Rev_Shipping <- readRDS("Data/Results/Model_Estimates_With_Shipping_Costs.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_2018 <- readRDS("Data/Results/Model_Estimates_With_Shipping_Costs_2018.Rds") +#Rev_No_Shipping_2018 <- readRDS("Data/Results/Model_Estimates_2018.Rds") +#Rev_Old_Sites <- readRDS("Data/Results/Model_Estimates_Old_Sites_Cost_Addition.Rds") -TEST <- Rev_No_Shipping + + +#TEST <- Rev_No_Shipping #TEST <- Rev_No_Shipping_2018 +TEST <- Rev_Old_Sites +TEST <- TEST %>% mutate(Discount=round(Discount,5)) + TEST <- TEST %>% mutate(Current_Profit=Profit/((1+Discount)^(Year-2026))) %>% ungroup TEST <- TEST %>% group_by(Year,Location,Discount) %>% filter(Current_Profit==max(Current_Profit)) %>% ungroup -ggplot(TEST %>% filter(Discount==0.03) ,aes(x=Year,y=Current_Profit/10^9,color=Capacity))+geom_line()+facet_wrap(Discount~Location,ncol=1)+scale_x_continuous(breaks=seq(1960,2083,by=5)) +ggplot(TEST %>% filter(Discount %in% c(0.03,0.045,0.05,0.07,0.1)) ,aes(x=Year,y=Current_Profit/10^9,color=Capacity))+geom_line()+facet_wrap(Discount~Location,ncol=1)+scale_x_continuous(breaks=seq(1960,2083,by=5)) -ggplot(TEST %>% filter(Discount==0.05) ,aes(x=Year,y=Current_Profit/10^9,color=Capacity))+geom_line()+facet_wrap(Discount~Location,ncol=1)+scale_x_continuous(breaks=seq(1960,2083,by=5)) +ggplot(TEST %>% filter(Discount%in% c(0.03,0.045)) ,aes(x=Year,y=Current_Profit/10^9,color=Capacity))+geom_line()+facet_wrap(Discount~Location,ncol=1)+scale_x_continuous(breaks=seq(1960,2083,by=5)) ggplot(TEST %>% filter(Discount %in% c(0.03,0.05,0.07) ,Year>1970) ,aes(x=Year,y=Current_Profit/10^9,color=Capacity))+geom_line()+facet_wrap(Discount~Location,ncol=2)+scale_x_continuous(breaks=seq(1960,2083,by=5)) ggplot(TEST %>% filter(Discount %in% c(0.05) ,Year>1970) ,aes(x=Year,y=Profit/10^9,color=Location))+geom_line()+scale_x_continuous(breaks=seq(1960,2083,by=5)) -ggplot(TEST %>% filter(Discount %in% c(0.03,0.05,0.07),Year>1970,Year<2050) ,aes(x=Year,y=Marginal/10^6))+ geom_area(aes(y=ifelse(Marginal>0,Marginal/10^6,0)),fill="lightgreen",alpha=0.6)+ geom_area(aes(y=ifelse(Marginal<0,Marginal/10^6,0)),fill="firebrick2",alpha=0.5) +facet_wrap(Discount~Location,ncol=2,scales = "free_y")+scale_x_continuous(breaks=seq(1960,2083,by=5))+scale_y_continuous(breaks=seq(-300,300,by=50))+geom_hline(yintercept=0)+theme_bw()+ylab("Profit Change (Million Dollars)") +ggplot(TEST %>% filter(Discount %in% c(0.03,0.045,0.05,0.07,0.15),Year>1970,Year<2045) ,aes(x=Year,y=Marginal/10^6))+ geom_area(aes(y=ifelse(Marginal>=0,Marginal/10^6,0)),fill="lightgreen",alpha=0.6)+ geom_area(aes(y=ifelse(Marginal<0,Marginal/10^6,0)),fill="firebrick2",alpha=0.5) +facet_wrap(Discount~Location,ncol=2,scales = "free_y")+scale_x_continuous(breaks=seq(1960,2083,by=5))+scale_y_continuous(breaks=seq(-300,300,by=50))+geom_hline(yintercept=0)+theme_bw()+ylab("Profit Change (Million Dollars)") diff --git a/Scripts/5_All_Stages.r b/Scripts/5_All_Stages.r new file mode 100644 index 0000000..26c9b9e --- /dev/null +++ b/Scripts/5_All_Stages.r @@ -0,0 +1,196 @@ +#A script which attempts to pull in all data, and create a data frame with the maximum revenue values for each facility, year and discount rate. The output can then be used to make figures and graphs +library(tidyverse) +library(lpSolve) #For solving discrete value maximization for the power plants +library(ompr) +library(ompr.roi) +library(ROI) +library(ROI.plugin.glpk) +source("./Scripts/Data_Setup.r") + +# 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 +# caps: length-m vector of knapsack capacities +# costs: length-m vector of fixed costs for using each knapsack +# force_k1_cost: if TRUE, forces paying cost of knapsack 1 even if unused (y1 = 1) +# bigM: "auto" (default) or length-m numeric vector for linking constraints +# returns: list with assignment, activated knapsacks, objective, capacity usage, etc. +solve_fixed_charge_multi_knapsack <- function(profit, weights, caps, costs, + force_k1_cost = TRUE, + bigM = c("auto")) { + stopifnot(is.matrix(profit)) + n <- nrow(profit) + m <- ncol(profit) + stopifnot(length(caps) == m) + stopifnot(length(costs) == m) + + # Normalize weights to an n x m matrix + if (is.vector(weights)) { + stopifnot(length(weights) == n) + W <- matrix(rep(weights, m), nrow = n, ncol = m) + } else if (is.matrix(weights)) { + stopifnot(nrow(weights) == n, ncol(weights) == m) + W <- weights + } else { + stop("`weights` must be a length-n vector or an n x m matrix.") + } + + # Disallow item-knapsack placements where profit is NA + allowed <- !is.na(profit) + profit[!allowed] <- 0 # values won't count; we'll also fix x=0 via constraints + + # Compute Big-M for linking constraints: sum_i x[i,k] <= M[k] * y[k] + if (length(bigM) == 1 && bigM[1] == "auto") { + # A tighter M is min(n, floor(cap_k / min_weight_k)) + min_w_by_k <- apply(W, 2, function(col) { + mw <- min(col, na.rm = TRUE) + if (is.infinite(mw) || is.na(mw) || mw <= 0) NA_real_ else mw + }) + M_auto <- rep(n, m) # fallback + for (k in 1:m) { + if (!is.na(min_w_by_k[k]) && is.finite(caps[k]) && caps[k] > 0) { + M_auto[k] <- min(n, floor(caps[k] / min_w_by_k[k])) + if (!is.finite(M_auto[k]) || M_auto[k] < 1) M_auto[k] <- n + } + } + M <- M_auto + } else { + stopifnot(length(bigM) == m, all(bigM >= 1)) + M <- bigM + } + + model <- MIPModel() %>% + # Decision vars: + add_variable(x[i, k], i = 1:n, k = 1:m, type = "binary") %>% # assignment + add_variable(y[k], k = 1:m, type = "binary") %>% # knapsack usage + # Objective: maximize total profit - fixed costs of used knapsacks + set_objective( + sum_expr(profit[i, k] * x[i, k], i = 1:n, k = 1:m) - + sum_expr(costs[k] * y[k], k = 1:m), + "max" + ) %>% + # Capacity per knapsack + add_constraint(sum_expr(W[i, k] * x[i, k], i = 1:n) <= caps[k], k = 1:m) %>% + # Each item assigned to at most one knapsack + add_constraint(sum_expr(x[i, k], k = 1:m) <= 1, i = 1:n) %>% + # Linking: using knapsack k (y[k]=1) if and only if at least one item is assigned. + # We need only the "if assigned then y=1" direction: + # sum_i x[i,k] <= M[k] * y[k] + add_constraint(sum_expr(x[i, k], i = 1:n) <= M[k] * y[k], k = 1:m) + + # Force paying cost for knapsack 1 even if unused + if (isTRUE(force_k1_cost)) { + model <- model %>% add_constraint(y[1] == 1) + } + + # Disallow NA placements explicitly + disallowed_idx <- which(!allowed, arr.ind = TRUE) + if (nrow(disallowed_idx) > 0) { + for (r in seq_len(nrow(disallowed_idx))) { + i <- disallowed_idx[r, 1] + k <- disallowed_idx[r, 2] + model <- model %>% add_constraint(x[i, k] == 0) + } + } + + result <- solve_model(model, with_ROI(solver = "glpk")) + + # Extract solutions + sol_x <- get_solution(result, x[i, k]) + sol_y <- get_solution(result, y[k]) + + # Assigned pairs + chosen <- subset(sol_x, value > 0.5) + assignment <- rep(NA_integer_, n) + if (nrow(chosen) > 0) assignment[chosen$i] <- chosen$k + + # Activated knapsacks (paying costs) + y_on <- rep(FALSE, m) + if (nrow(sol_y) > 0) y_on[sol_y$k] <- sol_y$value > 0.5 + + # Compute totals + total_profit <- if (nrow(chosen) > 0) sum(profit[cbind(chosen$i, chosen$k)]) else 0 + total_cost <- sum(costs[y_on]) + objective <- total_profit - total_cost + + # Capacity usage + cap_used <- numeric(m) + if (nrow(chosen) > 0) { + for (k in 1:m) { + idx <- chosen$i[chosen$k == k] + if (length(idx) > 0) { + cap_used[k] <- sum(W[cbind(idx, rep(k, length(idx)))]) + } + } + } + + list( + status = result$status, + objective_value = objective, + total_profit = total_profit, + total_fixed_cost = total_cost, + assignment = assignment, # per item: knapsack index (NA if unassigned) + knapsack_used = y_on, # logical vector of length m + capacity_used = cap_used, + capacity_remaining = caps - cap_used, + chosen_pairs = chosen[, c("i", "k", "value")], + profit_matrix = profit, + weight_matrix = W, + costs = costs, + bigM = M + ) +} +############################################### + +COST_DATA <- COST_DATA %>% group_by(Location,Discount) %>% arrange(Location,Discount,Capacity) %>% mutate(Cost_Change=lead(Cost)-Cost) %>% filter(Capacity==min(Capacity)) %>% ungroup + + + +#r <- 0.05 +#FACILITY <- 'New Mexico' +#CURRENT_YEAR <- 2026 +#CIFS_COST_DATA <- COST_DATA +#REACTOR_DATA <- TOTAL +#SAVED_PER_YEAR=CV + +OPTIMIZE_CIFS <- function(CURRENT_YEAR,r=0.05,FACILITY='New Mexico',CIFS_COST_DATA=COST_DATA,REACTOR_DATA=TOTAL,SAVED_PER_YEAR=CV){ + REACTOR_DATA <- REACTOR_DATA %>% filter(CLOSE_YEAR-CURRENT_YEAR<30) + CURRENT_VALUES <- CIFS_COST_DATA %>% filter(Discount==r,Location==FACILITY) + DISCOUNTS <- 1/((1+r)^(0:29)) + #Cost per addtion + costs <- DISCOUNTS*c(pull(CURRENT_VALUES ,Cost),rep(pull(CURRENT_VALUES ,Cost_Change),29)) + # Capacities + caps <- c(pull(CURRENT_VALUES ,Capacity),rep(5000,29)) + # Weights per item (volume of SNF per facility) + weights <- REACTOR_DATA %>% pull(Total_Tons) + #Revenue added per reactor added + YEARS_UNTIL_CLOSE <- ifelse(REACTOR_DATA$CLOSE_YEAR<=CURRENT_YEAR,0,REACTOR_DATA$CLOSE_YEAR-CURRENT_YEAR) + YEARS_UNTIL_CLOSE <- ifelse(YEARS_UNTIL_CLOSE >30,30,YEARS_UNTIL_CLOSE ) + +YEARS_UNTIL_CLOSE +CURRENT_VALUE <- function(YEARS_UNTIL_CLOSE_LIST,YEARS_AHEAD){ + YEARS_UNTIL_CLOSE_LIST <- ifelse(YEARS_UNTIL_CLOSE>YEARS_AHEAD,YEARS_UNTIL_CLOSE-YEARS_AHEAD,0) + TOTAL_VALUE <- sum(1/(1+r)^1:1000) + PROFIT_LIST <- c() + for(i in YEARS_UNTIL_CLOSE_LIST){ + PROFIT_LIST[length(PROFIT_LIST)+1] <- SAVED_PER_YEAR*(TOTAL_VALUE-ifelse(i>0,sum(1/((1+r)^1:i)),0)) + } + PROFIT_LIST <- PROFIT_LIST/((1+r)^YEARS_AHEAD) + return(PROFIT_LIST) +} +profit <- sapply(1:30,CURRENT_VALUE,YEARS_UNTIL_CLOSE_LIST=YEARS_UNTIL_CLOSE) +#PROFIT_LIST <- SAVED_PER_YEAR*PROFIT_LIST +# PROFIT_MAT <- matrix(sapply(YEARS_UNTIL_CLOSE, function(x){c(rep(0,x),rep(CV,30-x))}),ncol=30,byrow=TRUE) +# DISCOUNT_MAT <- matrix(DISCOUNTS,ncol=30,nrow=nrow(PROFIT_MAT),byrow=TRUE) +# 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 +) + return(res ) +} +YEARS <- 1970:2056 +#ADJUST_FACTOR <- 1/((1+0.05)^(YEARS-2026)) +TEST <- mclapply(YEARS,OPTIMIZE_CIFS,mc.cores=NCORES) +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") diff --git a/Scripts/Data_Setup.r b/Scripts/Data_Setup.r new file mode 100644 index 0000000..c30fef1 --- /dev/null +++ b/Scripts/Data_Setup.r @@ -0,0 +1,102 @@ +#A script which attempts to pull in all data, and create a data frame with the maximum revenue values for each facility, year and discount rate. The output can then be used to make figures and graphs +library(tidyverse) +library(parallel) + NCORES <- detectCores()-1 +library(lpSolve) #For solving discrete value maximization for the power plants +####Manual inputs + #Range of discount rates to calculate in the model. Each facility will have each rate calculated, so more values slows the results but allows for more discount rates to be reported in the findings. +#DISCOUNT_RATE_LIST <- c(0.03,0.05,0.07) +#DISCOUNT_RATE_LIST <- c(0.03,0.0325,0.035,0.0375,0.04,0.045,0.0475,0.05,0.07,0.1) +DISCOUNT_RATE_LIST <- seq(0.01,0.15,by=0.0025) + #The cost per ton of shipping uranium, used to see what can be shipped on day one of the project. +SHIPPING_COST_PER_TON <- 1.2874*26000 #Inflation adjusted from New Mexico Report + #The savings per year of having a CIFS at a served reactor (cost to house the CIFS). +CV <- 1.2874*(6984013) #Data from New Mexico Report, Converted from 2019 to Dec 2025 + #Locations to save results +RES_DIR <- "./Data/Results/" +INTERMEDIATE_DIR <- "./Data/Results/Separate_Costs_and_Benefits_Data/" + #Directory where the individual CIFS project plan cost data can be found. + #This has data has already been combined with the original files (low,high) and inflation adjusted. +CIFS_INDIVIDUAL_COST_DATA_DIR <- 'Data/Cleaned_Data/' +#Create any need save locations + dir.create(RES_DIR,recursive=TRUE,showWarnings=FALSE) + dir.create(INTERMEDIATE_DIR,recursive=TRUE,showWarnings=FALSE) + +###################################Cost results +CIFS <- rbind(readRDS(paste0(CIFS_INDIVIDUAL_COST_DATA_DIR,"Texas_CIFS_Costs.Rds")),readRDS(paste0(CIFS_INDIVIDUAL_COST_DATA_DIR,"New_Mexico_CIFS_Costs.Rds"))) +#Adjust for inflation +#CIFS[,-1:-5] <- 1.2874*CIFS[,-1:-5] +COSTS <- do.call(rbind,lapply(DISCOUNT_RATE_LIST ,function(DISCOUNT){CIFS %>% group_by(Location,Capacity,Cost_Assumption) %>% mutate(NPC=Total/(1+DISCOUNT)^Year) %>% summarize(Discount=DISCOUNT,Costs=sum(NPC))})) %>% ungroup + +TEMP <- COSTS%>% group_by(Location,Cost_Assumption,Discount) %>% summarize(ST_COST=min(Costs),ST_CAPACITY=min(Capacity),M_COST=(max(Costs)-min(Costs))/(max(Capacity)-min(Capacity))) %>% ungroup +CAPACITY_INCREMENT <- 5000 +#rm(COST_DATA) +for(i in 1:nrow(TEMP)){ + LOC <- as.character(TEMP[i,"Location"]) + COST_LEVEL <- as.character(TEMP[i,"Cost_Assumption"]) + DISCOUNT <- as.numeric(TEMP[i,"Discount"]) + ST_CAP <- as.numeric(TEMP[i,"ST_CAPACITY"]) + ST_COST <- as.numeric(TEMP[i,"ST_COST"]) + COST_SLOPE <- as.numeric(TEMP[i,]$M_COST) + CAPACITY <- seq(ST_CAP,160000,by=5000) + COST <- ST_COST+COST_SLOPE*CAPACITY_INCREMENT*(0:(length(CAPACITY)-1)) + + C_RES <- cbind(CAPACITY,COST) %>% as_tibble %>% mutate(Location=LOC,Cost_Assumption=COST_LEVEL,Discount=DISCOUNT) %>% select(Location,Cost_Assumption,Discount,Capacity=CAPACITY,Cost=COST) + if(!exists("COST_DATA")){COST_DATA <- C_RES}else{COST_DATA<- rbind(COST_DATA,C_RES)} +} +saveRDS(COST_DATA ,paste0(INTERMEDIATE_DIR,"All_CIFS_Discounted_Costs.Rds")) +#COST_DATA <- COST_DATA %>% filter(Cost_Assumption=='Average') %>% select(-Cost_Assumption) %>% unique +COST_DATA <- COST_DATA %>% filter(Cost_Assumption=='High') %>% select(-Cost_Assumption) %>% unique + + +##All unique capacity levels that the revenues need to be calculated for +CAPACITY_LIST <- COST_DATA %>% pull(Capacity) %>% unique + +### +TOTAL <- read_csv("Data/Raw_Data/Curie_Spent_Fuel_Site_Totals.csv") %>% mutate(OP_YEAR=year(Op_Date_Min),CLOSE_YEAR=year(Close_Date_Max))%>% select(Facility,Total_Assemblies,Total_Tons,OP_YEAR,CLOSE_YEAR) +FACILITY_LIST <- TOTAL %>% pull(Facility) +#https://www.nrc.gov/reactors/operating/licensing/renewal/subsequent-license-renewal +SUBMITTED <-rbind(c(FACILITY_LIST[str_detect(FACILITY_LIST,"Duane*" )],2025), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Nine Mile*" )],2026), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Ginna*" )],2026), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Cooper*" )],2026), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Farley*" )],2027), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Prairie*" )],2027), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Brunswick*" )],2027), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Cook" )],2027), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Hope" )],2027), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Salem" )],2027), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Perry" )],2027), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Millstone" )],2028), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Palisades" )],2028), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Beaver" )],2028), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Callaway" )],2029), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Three Mile Island" )],2029), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Davis-Besse" )],2029), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Wolf" )],2030), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Lucie" )],2021), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Robinson" )],2025), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Hatch" )],2025)) +SUBMITTED <- SUBMITTED %>% as_tibble +colnames(SUBMITTED ) <- c("Facility","App_Date","Status") +SUBMITTED <- SUBMITTED%>% mutate(Status="Applied",App_Date=as.numeric(App_Date)) %>% select(Facility,Status,App_Date) + +#Issued +RENEWED <- rbind(c(FACILITY_LIST[str_detect(FACILITY_LIST,"Turkey" )],"Granted",2018,2033), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Peach" )],"Granted",2019,2034), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Surry" )],"Granted",2020,2033), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"North" )],"Granted",2021,2040), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Monticello" )],"Granted",2024,2030), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Oconee" )],"Granted",2025,2034), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Summer" )],"Granted",2025,2042), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Beach" )],"Granted",2025,2033), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Browns" )],"Granted",2025,2036), +c(FACILITY_LIST[str_detect(FACILITY_LIST,"Dresden" )],"Granted",2025,2031)) +RENEWED <- RENEWED %>% as_tibble +colnames(RENEWED) <- c("Facility","Status","App_Date","Op_Date") +RENEWED <- RENEWED %>% mutate(Op_Date=as.numeric(Op_Date),App_Date=as.numeric(App_Date)) +AVG_LENGTH <- RENEWED %>% mutate(DIFF=Op_Date-App_Date) %>% pull(DIFF) %>% mean %>% round +SUBMITTED <- SUBMITTED %>% mutate(Op_Date=App_Date+AVG_LENGTH) +UPDATE <- rbind(RENEWED,SUBMITTED ) +TOTAL_ORIG <- TOTAL +TOTAL <- TOTAL %>% left_join(UPDATE) %>% mutate(CLOSE_YEAR=ifelse(Op_Date>CLOSE_YEAR & !is.na(Status),Op_Date,CLOSE_YEAR)) %>% select(-Status,-App_Date,-Op_Date)