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"

#########
#CURRENT_YEAR <- 30 
#END_YEAR <- 40
#PHASE_SIZE <- 5000
#ST_CAP <- 8680
#TRANSPORT_COST_RATIO <- 155462880/5000 # or 269883561/8680 since transportation cost scale linearly
#VOLUME <- 10 
#TRANS_CONSTRAINT <- (ST_CAP+5000*19)/19
#PHASES_ADDED <- max(ceiling((VOLUME-ST_CAP)/PHASE_SIZE),0) #Make sure added phases are not negative
#STARTING_VOLUME=8680
#PHASE_CONST_COST=103399272
#DISCOUNT <- 0.05
#VOLUME <- 8660+5000
#PHASES_ADDED <- 1
#DECOM_COST_PER_TON =24822656/5000
#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
	#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 
	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)
}

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)
	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.
ADDITION_CHECK(CURRENT,22,1)