Working to add CIFS stages to optimization

Analysis.r

@@ -1,24 +1,35 @@
 library(tidyverse)
-Rev_No_Shipping <- readRDS("Data/Results/Model_Estimates.Rds")
+#Rev_No_Shipping <- readRDS("Data/Results/Model_Estimates.Rds")
-Rev_No_Shipping_2018 <- readRDS("Data/Results/Model_Estimates_2018.Rds")
+#Rev_No_Shipping  %>% pull(Discount) %>% unique
-#Rev_Shipping <- readRDS("Data/Results/Model_Estimates_With_Shipping_Costs.Rds")
+
+#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)")  
 

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")

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)