# oba_emissions_impact_approximation.R
#
# This script uses IPM output from the mass-based scenario to approximate the 
# impact of implementing output based allocation (OBA) as proposed in the 
# federal plan. The algorithm considers each IPM region with existing NGCC and 
# fossil steam capacity and new NGCC generation. Each region is considered as an
# island with no interaction with other regions. For each qualifying region 
# model plants are sorted based on their marginal cost of generation including
# the CO2 allowance price for the state in which the model plant is located. The
# marginal cost of existing NGCC units is then adjusted downwards based on the 
# value of allowances from OBA. Given the new marginal generation costs, the 
# model plants in each region are resorted. For each case where an existing NGCC
# model plant passess an existing fossil steam model plant on this dispatch 
# curve approximation a change in the plants' generation is considered. The 
# existing NGCC model plant is ramped up to a capacity factor of 87% with the 
# model plant earning OBA allowances after a capacity factor of 50% at an 
# allocation rate of 1,030 lbs/MWh. Based on the emissions rate of the existing
# NGCC plant relative to the fossil steam plant, the generation of the fossil
# steam plant that could be displaced and keep the covered sources aggregate 
# emissions constant for the region is calculated. If the average cost per MWh
# of the additional generation from the existing NGCC model plant under OBA
# is less than for the existing fossil steam generation and new NGCC plant
# then the generation shift is made. The additional generation from the existing 
# NGCC plant over what is displaced at the fossil steam plant is subtracted from
# the new NGCC model plant in the region. Bounds are considered to ensure that 
# regional generation remains equal to the pre-OBA case and that no fossil steam
# plant goes below a capacity factor of 0%.
#
#
# Directions:
#
#   1.) Open the IPM RPE replcement file in Excel. Go to 'File' and 'Save As'. 
#       Select the save as tpye 'CSV (Comma delimited) (*.csv)' and save with 
#       the file name 'RPE Replacement File.csv'.
#
#   2.) Open the IPM summary file in Excel and click on the sheet 
#       'All Constratints'. Go to 'File' and 'Save As'. Select the save as tpye 
#       'CSV (Comma delimited) (*.csv)' and save with the file name
#       'All Constraints.csv'.
#
#   3.) Run this script in R. The ouput will be written to the file
#       'oba approximation results 2030.csv'.
#
# 
# R Interpreter:
#
#   The R interpreter to run this script is free software avialble at 
#  https://www.r-project.org/
#
#
# Output:
#
#   The results are stroed in a CSV file 'oba approximation results 2030.csv', 
#   which contains the following columns for each IPM region that has existing
#   and new NGCC units and existing fossil steam units.
#
#   Region                            IPM regiona abbreviation
#   Steam.Capacity                    Existing fossil steam capacity [MW]
#   Exist.NGCC.Capacity               Existing NGCC capacity [MW]
#   New.NGCC.Capacity                 New NGCC capacity [MW]
#   Steam.Generation                  Existing fossil steam generation [MWh]
#   Exist.NGCC.Generation             Existing NGCC generation [MWh]
#   New.NGCC.Generation               New NGCC generation [MWh]
#   Steam.Capacity.Factor             Existing fossil steam capacity factor for region
#   Exist.NGCC.Capacity.Factor        Existing NGCC capacity factor for region
#   New.NGCC.Capacity.Factor          New NGCC capacity factor for region
#   Steam.Generation.OBA              Existing fossil steam generation with OBA [MWh]
#   Exist.NGCC.Generation.OBA         Existing NGCC generation with OBA [MWh]
#   New.NGCC.Generation.OBA           New NGCC generation with OBA [MWh]
#   Steam.Capacity.Factor.OBA         Existing fossil steam capacity factor for region with OBA
#   Exist.NGCC.Capacity.Factor.OBA    Existing NGCC capacity factor for region with OBA
#   New.NGCC.Capacity.Factor.OBA      New NGCC capacity factor for region with OBA
#   New.NGCC.Generation.Change        Percentage change in new NGCC generation with OBA [%]
#   Emissions                         CO2 emissions without OBA [million tons]
#   Emissions.Change                  Percentage change in CO2 emissions with OBA [%]
#   Emissions.Diff                    Change in CO2 emissions with OBA [million tons]
#
#



# Parameters for analysis
# ------------------------------------------------------------------------------

# Year of analysis
year = 2030

# Minimum capacity for applicability requirements [MW]
min.capacity = 25

# Maximum capacity factor that existing NGCC can reach
max.ngcc.capacity.factor = 0.87

# Allocation rate for OBA [lbs/MWh]
allocation.rate = 1030

# Allocation capacity factor (above this OBA is triggered)
allocation.cf = .50

# Existing NGCC capacity reporting type
exist.ngcc.label = '00 Exist Combined Cycle'

# New NGCC capacity reporting type
new.ngcc.label = '00 New Combined Cycle'

# Existing fossil steam capacity reporting types
# This does not include other fossil ('00 Exist Fossil_Other') which covers waste
fossil.steam.label = c('00 Exist Coal Steam_Uncontrolled',
                       '00 Exist Coal Steam_DryFGD_SCR',
                       '00 Exist Coal Steam_DryFGD_SNCR',
                       '00 Exist Coal Steam_DryFGD',
                       '00 Exist Oil/Gas Steam_SNCR',
                       '00 Exist FBC_DryFGD',
                       '00 Exist FBC_Uncontrolled',
                       '00 Exist FBC_SNCR',
                       '00 Exist Oil/Gas Steam',
                       '00 Exist Oil/Gas Steam_SCR',
                       '00 Exist Coal Steam_WetFGD_SCR_ACI',
                       '00 Exist Coal Steam_WetFGD_SNCR_ACI',
                       '00 Exist Coal Steam_WetFGD_ACI',
                       '00 Exist Coal Steam_DryFGD_SCR_ACI',
                       '00 Exist Coal Steam_DryFGD_SNCR_ACI',
                       '00 Exist Coal Steam_DryFGD_ACI',
                       '00 Exist Coal Steam_ACI',
                       '00 Exist FBC_SNCR_ACI',
                       '00 Exist FBC_ACI',
                       '00 Exist IGCC',
                       '00 Exist FBC_DryFGD_SNCR',
                       '00 Exist Coal Steam_SNCR',
                       '00 Exist FBC_DryFGD_SNCR_ACI',
                       '00 Exist Coal Steam_ACI_DSI',
                       '00 Exist FBC_WetFGD_SNCR_ACI',
                       '00 Exist Coal Steam_WetFGD',
                       '00 Exist Coal Steam_WetFGD_SCR_DSI',
                       '00 Exist Coal Steam_WetFGD_SCR_ACI_DSI',
                       '00 Exist Coal Steam_WetFGD_ACI_DSI',
                       '00 Exist Coal Steam_DryFGD_SCR_ACI_DSI',
                       '00 Exist Coal Steam_WetFGD_SNCR_ACI_DSI',
                       '00 Exist Coal Steam_WetFGD_SCR',
                       '00 Exist Coal Steam_WetFGD_SNCR',
                       '00 Exist Coal Steam_SCR_ACI',
                       '00 Ret.DryFGD',
                       '00 Ret.HRI',
                       '00 Ret.SCR',
                       '00 Ret.ACI',
                       '00 Ret.C2G',
                       '00 Ret.WetFGD',
                       '00 Ret.DSI',
                       '00 Ret.ACI & WetFGD',
                       '00 Ret.ExistDryFGD & ACI',
                       '00 Ret.FB & ACI',
                       '00 Ret.Oil/Gas Steam SCR',
                       '00 Ret.ExistSNCR & ACI & WetFGD',
                       '00 Ret.ExistACI_DSI & DryFGD',
                       '00 Ret.ExistACI_DSI & WetFGD',
                       '00 Ret.ExistWetFGD & ACI & SCR',
                       '00 Ret.ExistSNCR & ACI & DSI',
                       '00 Ret.ExistACI_DSI & SCR',
                       '00 Ret.C2G & SCR',
                       '00 Ret.ACI & DSI',
                       '00 Ret.ExistSCR_DryFGD & HRI',
                       '00 Ret.ExistSCR_DryFGD & ACI',
                       '00 Ret.ExistSNCR_DryFGD & HRI',
                       '00 Ret.ExistDryFGD & HRI',
                       '00 Ret.ExistDryFGD & SCR',
                       '00 Ret.FB & HRI',
                       '00 Ret.FB-ExistSNCR & ACI',
                       '00 Ret.FB-ExistSNCR & HRI',
                       '00 Ret.ExistSCR_WetFGD_ACI & HRI',
                       '00 Ret.ExistSNCR_WetFGD_ACI & HRI',
                       '00 Ret.ExistWetFGD_ACI & HRI',
                       '00 Ret.ExistWetFGD_ACI & SCR',
                       '00 Ret.ExistSCR_DryFGD_ACI & HRI',
                       '00 Ret.ExistDryFGD_ACI & SCR',
                       '00 Ret.ExistSNCR_DryFGD_ACI & HRI',
                       '00 Ret.ExistDryFGD_ACI & HRI',
                       '00 Ret.ExistSNCR_ACI & DryFGD',
                       '00 Ret.ExistACI & DSI',
                       '00 Ret.ExistACI & SCR & WetFGD',
                       '00 Ret.ExistSCR & ACI & DSI',
                       '00 Ret.ExistSCR & WetFGD',
                       '00 Ret.FB-ExistDryFGD_SNCR & HRI',
                       '00 Ret.ExistSNCR & HRI',
                       '00 Ret.FB-ExistDryFGD_SNCR_ACI & HRI',
                       '00 Ret.ExistACI_DSI & HRI',
                       '00 Ret.ExistWetFGD & ACI',
                       '00 Ret.ExistWetFGD & SCR',
                       '00 Ret.ExistWetFGD & HRI',
                       '00 Ret.ExistSCR_WetFGD_ACI_DSI & HRI',
                       '00 Ret.ExistWetFGD_ACI_DSI & HRI',
                       '00 Ret.ExistSCR_DryFGD_ACI_DSI & HRI',
                       '00 Ret.ExistSCR_WetFGD & ACI',
                       '00 Ret.ExistSCR_WetFGD & HRI',
                       '00 Ret.ExistSNCR_WetFGD & HRI',
                       '00 Ret.ExistSNCR_WetFGD & ACI',
                       '00 Ret.ACI & DSI & HRI',
                       '00 Ret.ACI & DryFGD & HRI',
                       '00 Ret.WetFGD & HRI',
                       '00 Ret.ACI & HRI',
                       '00 Ret.ACI & DryFGD',
                       '00 Ret.ACI & WetFGD & HRI',
                       '00 Ret.ExistDryFGD_SCR & ACI & HRI',
                       '00 Ret.ExistDryFGD & ACI & HRI',
                       '00 Ret.ExistDryFGD & ACI & SCR',
                       '00 Ret.FB & ACI & HRI',
                       '00 Ret.ExistWetFGD_ACI & SCR & HRI',
                       '00 Ret.ExistDryFGD_ACI & SCR & HRI',
                       '00 Ret.ExistSCR_ACI & DryFGD & HRI',
                       '00 Ret.ExistSNCR_ACI & WetFGD & HRI',
                       '00 Ret.ExistACI & DSI & HRI',
                       '00 Ret.ExistSCR & ACI & DSI & HRI',
                       '00 Ret.ExistSNCR & ACI & DSI & HRI',
                       '00 Ret.SCR & HRI',
                       '00 Ret.ExistSNCR & ACI & WetFGD & HRI',
                       '00 Ret.ExistACI_DSI & DryFGD & HRI',
                       '00 Ret.ExistACI_DSI & WetFGD & HRI',
                       '00 Ret.ExistWetFGD & SCR & HRI',
                       '00 Ret.ExistWetFGD & ACI & HRI',
                       '00 Ret.ExistWetFGD_SCR & ACI & HRI',
                       '00 Ret.ExistWetFGD_SNCR & ACI & HRI',
                       '00 Ret.ExistACI_DSI & SCR & WetFGD & HRI',
                       '00 Ret.ExistWetFGD & ACI & SCR & HRI')



# Load in the IPM output
# ------------------------------------------------------------------------------

# Load in the RPE file for the baseline
rpe = read.csv('RPE Replacement File.csv',stringsAsFactors=FALSE)

# Load in the constraint data
const = read.csv('All Constraints.csv',stringsAsFactors=FALSE)



# Extract the shadow prices on the state mass-based constraints
# ------------------------------------------------------------------------------

# Location of NSPS shadow prices
ind = substr(const[,1],1,26)=='#NSPS State CO2 Constraint' & const$Data.Field=='Constraint Shadow Price'

# To avoid duplication with Utah abbreviated UTE as UE
const[substr(const[,1],29,31)=='UTE',1] = '#NSPS State CO2 Constraint -UE'

# Extract the allowance prices
allowance.price = matrix(NA,nrow=50,ncol=7)
for (i in 1:dim(allowance.price)[2])
  allowance.price[,i] = as.numeric(const[ind,9+i])
rownames(allowance.price) = substr(const[ind,1],29,30)
colnames(allowance.price) = c(2016,2018,2020,2025,2030,2040,2050)



# Remove unnecessary information from the RPE data frame
# ------------------------------------------------------------------------------

# Pull out the data for the specific year
rpe = rpe[rpe$Year==year,]

# Keep only applicable fossil units and NGCC
rpe = rpe[rpe$Capacity.Reporting.Type %in% c(exist.ngcc.label,new.ngcc.label,
                                             fossil.steam.label),]
rpe = rpe[rpe$Dispatchable.Capacity.MW>=min.capacity,]

# Drop Canada
rpe = rpe[substr(rpe$Region.Group.Acronym,1,2)!='CN',]

# Remove biomass co-firing
rpe = rpe[rpe$Fuel.Type!='Biomass',]



# Add the allowance prices and total allowance cost to the RPE data frame
# ------------------------------------------------------------------------------

# States abbreviations
states = rownames(allowance.price)

# Add state abbreviations to the RPE data
rpe$state = substr(rpe$Unit.Long.Name,5,6)

# Add the allowance prices for the units except where new gas doesn't have a state
for (i in 1:dim(rpe)[1]) {
  if (rpe$state[i]!="  ")
    rpe$allowance.price[i] = allowance.price[states==rpe$state[i],as.character(year)]
}

# New NGCC isn't covered
rpe$allowance.price[rpe$Capacity.Reporting.Type==new.ngcc.label] = 0

# Add in the total CO2 emissions cost
rpe$allowance.cost = rpe$CO2.Emissions.Total.Thousand.Tons*1000*rpe$allowance.price



# Create the dispatch curve of model plants by region
# ------------------------------------------------------------------------------

# Get the list of regions
regions = unique(rpe$Region.Group.Acronym)

# New container for regional outptu
data = list()

# Generation by region
generation = NULL

# Loop through each region and create the dispatch curve
for (r in regions) {
  
  # Set of fossil steam and existing NGCC units in the region and year
  ind = rpe$Region.Group.Acronym==r & rpe$Capacity.Reporting.Type %in% c(fossil.steam.label,exist.ngcc.label)
  
  # Set of new NGCC units in the region and year
  ind.new = rpe$Region.Group.Acronym==r & rpe$Capacity.Reporting.Type==new.ngcc.label
  
  if (any(ind) & any(ind.new)) {
    
    # New NGCC generation in the region
    new.g = sum(rpe$Generation.Total.GWh[ind.new]*1000)
    
    # "Total" generation in the region 
    total.generation = new.g
    
    # "Variable" cost of new NGCC generation including FOM and capital costs 
    # since these would not be incurred if the "unit" was not constructed due to
    # different incentives under OBA
    c = new.g*sum(rpe$VOM...MWh[ind.new])+sum(rpe$Fuel...MMBtu[ind.new]*rpe$Fuel.Consumption.Total.TBtu[ind.new]*1000000)+sum((rpe$FOM...kW.year[ind.new]+rpe$Capital.Cost...kw.year[ind.new])*rpe$Dispatchable.Capacity.MW[ind.new]*1000)
    
    # Marginal cost of new NGCC generation
    new.mc = c/new.g
    
    # Emissions rate of new NGCC
    e = sum(rpe$CO2.Emissions.Total.Thousand.Tons[ind.new]*1000)/new.g
    
    # Existing fossil Units in the region
    units = rpe[ind,]
    
    # Initialize the output data frame for the region's dispatch curve
    data[[r]] = data.frame(type=NA,unit.id=NA,state=NA,mc=NA,mc.oba=NA,capacity=NA,generation=NA,max.generation=NA,emissions.rate=NA)

    # Fill in the output for new NGCC in the region
    data[[r]][1,] = c('new.ngcc',NA,rpe$state[ind.new],new.mc,new.mc,sum(rpe$Dispatchable.Capacity.MW[ind.new]),new.g,new.g,e)
   
    for (i in 1:sum(ind)) {
      
      # Unit's generation [MWh]
      g = units$Generation.Total.GWh[i]*1000
      
      # VOM plus fuel costs plus emissions costs
      c = g*units$VOM...MWh[i]+units$Fuel...MMBtu[i]*units$Fuel.Consumption.Total.TBtu[i]*1000000+units$allowance.cost[i]
      
      # Marginal cost
      mc = c/g

      # Total regional generation
      total.generation = total.generation+g
      
      # Emissions rate
      e = units$CO2.Emissions.Total.Thousand.Tons[i]*1000/g
      
      # If the model plant is an NGCC unit calculate the maximum dispatchable 
      # generation under OBA and adjust the MC for the OBA case
      if (units$Capacity.Reporting.Type[i]==exist.ngcc.label){
        type = 'ngcc'
        max.g = max.ngcc.capacity.factor*units$Dispatchable.Capacity.MW[i]*365*24
        mc.oba = mc-allocation.rate/2000*units$allowance.price[i]
        
      # If the model plant is a fossil steam unit then nothing changes in the
      # OBA case
      } else {
        type = 'steam'
        mc.oba = mc
        max.g = g
      }

      # Add the model plant to the dispatch curve data frame
      data[[r]][i+1,] = c(type,units$UnitId[i],units$state[i],mc,mc.oba,units$Dispatchable.Capacity.MW[i],g,max.g,e)
      
    }
        
    # Total generation for the region
    generation[r] = total.generation
     
  }

}

# Sort each regions affected sources by the marginal cost in the original case
data = lapply(data,function(x) x[c(1,sort(as.numeric(x[-1,'mc']),index.return=TRUE)$ix+1),])

# Add unique identifiers for the units in each region
data = lapply(data,function(x) cbind(id=1:dim(x)[1],x))



# Approximate the new dispatch curve under the case of OBA
# ------------------------------------------------------------------------------

# Sort each regions affected sources by the marginal cost in the OBA case
data = lapply(data,function(x) x[c(1,sort(as.numeric(x[-1,'mc.oba']),index.return=TRUE)$ix+1),])

# Initialize the national additional emissions reductions from OBA
emissions.reductions = 0

# Loop through each region and build up the new dispatch curve
for (r in regions) {
  
  # If the region has no generation from applicable sources, skip it
  if (is.na(generation[r]))
    next
  
  # If the region has no fossil steam units, skip it
  if (!any(data[[r]][,'type']=='steam'))
    next
  
  # Initialize model plant generation in the OBA case
  data[[r]][,'generation.oba'] = data[[r]][,'generation']
  
  # Initialize model OBA allowances received to zero
  data[[r]][,'oba.allowances'] = 0
  
  # Locate the fossil steam units in the original and OBA dispatch curves
  steam.new = which(data[[r]][,'type']=='steam')
  
  # Loop through each fossil steam unit and check if it has been over taken by
  # an NGCC unit in the dispatch curve under OBA
  for (i in 1:length(steam.new)) {
    
    # Loop through each unit that is above the current fossil steam unit in
    # the OBA case but was not in the original case
    for (j in 2:steam.new[i]) {
      
      # If the unit that moved up is an NGCC unit adjust its generation in 
      # response
      if (data[[r]][j,'id']>data[[r]][steam.new[i],'id'] & data[[r]][j,'type']=='ngcc') {
        
        # Potential increase in the NGCC unit's generation
        potential = as.numeric(data[[r]][j,'max.generation'])-as.numeric(data[[r]][j,'generation.oba'])
        
        # Make sure we don't exceed total demand in the region
        potential = min(potential,as.numeric(generation[r])-sum(as.numeric(data[[r]][-1,'generation.oba'])))
        
        # Emissions of the NGCC unit at the potential generation level
        e = potential*as.numeric(data[[r]][j,'emissions.rate'])
        
        # Amount of the fossil steam unit's generation that could be displaced
        # without emissions exceeding the mass-based standard
        displaced.gen = e/as.numeric(data[[r]][steam.new[i],'emissions.rate'])
        
        # Make sure we are not displacing more generation than remains at the
        # current fossil steam unit
        displaced.gen = min(displaced.gen,as.numeric(data[[r]][steam.new[i],'generation.oba']))
        
        # Make sure that the potential increase in NGCC generation is adjusted
        # if the displaced generation was curtailed to not exceed its
        # remaining generation
        potential = displaced.gen*as.numeric(data[[r]][steam.new[i],'emissions.rate'])/as.numeric(data[[r]][j,'emissions.rate'])
        
        # Compute the cost of the additional NGCC generation given that it 
        # will only earn OBA above the allocation capacity factor
        original.mc.gen = min(potential,max(0,allocation.cf*365*24*as.numeric(data[[r]][j,'capacity'])-as.numeric(data[[r]][j,'generation'])))
        ngcc.cost = original.mc.gen*as.numeric(data[[r]][j,'mc'])+max(0,potential-original.mc.gen)*as.numeric(data[[r]][j,'mc.oba'])
        
        # Average cost of the fossil steam and new NGCC generation being 
        # displaced by the existing NGCC unit
        displaced.mc = displaced.gen/potential*as.numeric(data[[r]][steam.new[i],'mc'])+(1-displaced.gen/potential)*as.numeric(data[[r]][1,'mc'])
        
        # If the average MWh generation cost for NGCC with OBA is less than 
        # for the fossil steam and new NGCC then don't displace any generation
        if (potential>0 & (ngcc.cost/potential)>displaced.mc) {
          potential = 0
          displaced.gen = 0
        }
        
        # Increase the NGCC generation in the OBA case
        data[[r]][j,'generation.oba'] = as.numeric(data[[r]][j,'generation.oba'])+potential
        
        # Reduce the fossil steam unit's generation
        data[[r]][steam.new[i],'generation.oba'] = as.numeric(data[[r]][steam.new[i],'generation.oba'])-displaced.gen
        
        # Track the number of OBA allowances allocated due to the displacement
        data[[r]][j,'oba.allowances'] = data[[r]][j,'oba.allowances']+max(0,potential-original.mc.gen)*allocation.rate/2000
        
      }
      
    }
    
  }
  
  # Decrease the generation of new NGCC by the increase that occurs under the 
  # mass-based standard
  data[[r]][1,'generation.oba'] = generation[r]-sum(as.numeric(data[[r]][-1,'generation.oba']))
  
  # Calculate the model plant's total emissions in the original and OBA cases
  data[[r]][,'emissions'] = as.numeric(data[[r]][,'emissions.rate'])*as.numeric(data[[r]][,'generation'])
  data[[r]][,'emissions.oba'] = as.numeric(data[[r]][,'emissions.rate'])*as.numeric(data[[r]][,'generation.oba'])
  
  # Add this region's additional OBA emissions reductions to the national total
  emissions.reductions = emissions.reductions+sum(as.numeric(data[[r]][,'emissions'])-as.numeric(data[[r]][,'emissions.oba']))
  
}



# Perform validity checks on the OBA dispatch curve approximation
# ------------------------------------------------------------------------------

for (r in regions) {
  
  # If the region has no generation from applicable sources, skip it
  if (is.na(generation[r]))
    next
  
  # If the region has no fossil steam units, skip it
  if (!any(data[[r]][,'type']=='steam'))
    next
  
  # Check to make sure total demand is met
  if (abs(as.numeric(generation[r])-sum(as.numeric(data[[r]][,'generation.oba'])))>100)
    warning(paste('demand constraint violdated in',r))
  
  # Check to make sure no ngcc max capacity factors are violated
  if (any(as.numeric(data[[r]][,'generation.oba'])>as.numeric(data[[r]][,'max.generation'])))
    warning(paste('ngcc max cf constraint violdated in',r))
    
  # Check to make sure all generation is weakly positive
  if (any(as.numeric(data[[r]][,'generation.oba'])<0))
      warning(paste('weakly positive generation violdated in',r))
  
}



# Output results
# ------------------------------------------------------------------------------

# Percentage emissions change in each region
emissions.change = lapply(data,function(x) if ('emissions.oba'%in%names(x)) (sum(as.numeric(x[,'emissions.oba']))/sum(as.numeric(x[,'emissions']))-1)*100)

# Percentage change in new NGCC generation in each region
new.ngcc.change = lapply(data,function(x) (as.numeric(x[1,'generation.oba'])/as.numeric(x[1,'generation'])-1)*100)

# Initialize space for regional statistics
stats = matrix(NA,nrow=length(regions),ncol=19)

# Loop through each region and compile the OBA statistics
for (i in 1:length(regions)) {
  
  # Fossil steam capacity
  ind = rpe$Region.Group.Acronym==regions[i] & rpe$Capacity.Reporting.Type %in% fossil.steam.label
  stats[i,1] = sum(rpe$Dispatchable.Capacity.MW[ind])
  
  # Existing NGCC capacity
  ind = rpe$Region.Group.Acronym==regions[i] & rpe$Capacity.Reporting.Type %in% exist.ngcc.label
  stats[i,2] = sum(rpe$Dispatchable.Capacity.MW[ind])
  
  # New NGCC capacity
  ind = rpe$Region.Group.Acronym==regions[i] & rpe$Capacity.Reporting.Type %in% new.ngcc.label
  stats[i,3] = sum(rpe$Dispatchable.Capacity.MW[ind])
  
  # Generation in original case
  stats[i,4] = sum(as.numeric(data[[regions[i]]][data[[regions[i]]]$type=='steam','generation']))
  stats[i,5] = sum(as.numeric(data[[regions[i]]][data[[regions[i]]]$type=='ngcc','generation']))
  stats[i,6] = sum(as.numeric(data[[regions[i]]][data[[regions[i]]]$type=='new.ngcc','generation']))
  
  # Capacity factors in original case
  stats[i,7] = stats[i,4]/(stats[i,1]*365*24)
  stats[i,8] = stats[i,5]/(stats[i,2]*365*24)
  stats[i,9] = stats[i,6]/(stats[i,3]*365*24)
  
  # Generation in OBA case
  stats[i,10] = sum(as.numeric(data[[regions[i]]][data[[regions[i]]]$type=='steam','generation.oba']))
  stats[i,11] = sum(as.numeric(data[[regions[i]]][data[[regions[i]]]$type=='ngcc','generation.oba']))
  stats[i,12] = sum(as.numeric(data[[regions[i]]][data[[regions[i]]]$type=='new.ngcc','generation.oba']))
  
  # Capacity factors in original case
  stats[i,13] = stats[i,10]/(stats[i,1]*365*24)
  stats[i,14] = stats[i,11]/(stats[i,2]*365*24)
  stats[i,15] = stats[i,12]/(stats[i,3]*365*24)
  
  # Percentage change in new NGCC generation
  if (regions[i] %in% names(new.ngcc.change))
    stats[i,16] = unlist(as.numeric(new.ngcc.change[regions[i]]))

  # Emissions
  ind = rpe$Region.Group.Acronym==regions[i]
  stats[i,17] = sum(rpe$CO2.Emissions.Total.Thousand.Tons[ind]/1000)
  
  # Percentage change in emissions
  if (regions[i] %in% names(emissions.change)) {
    temp = as.numeric(unlist(emissions.change[regions[i]]))
    if (length(temp)!=0)
      stats[i,18] = temp
    else
      stats[i,18] = 0
  } else {
    stats[i,18] = 0}
  
  # Emissions change in million tons
  stats[i,19] = stats[i,17]*stats[i,18]/100
  
}

# Add row and column names to the OBA statistics
rownames(stats) = regions
colnames(stats) = c('Steam.Capacity','Exist.NGCC.Capacity','New.NGCC.Capacity',
                    'Steam.Generation','Exist.NGCC.Generation','New.NGCC.Generation',
                    'Steam.Capacity.Factor','Exist.NGCC.Capacity.Factor',
                    'New.NGCC.Capacity.Factor','Steam.Generation.OBA',
                    'Exist.NGCC.Generation.OBA','New.NGCC.Generation.OBA',
                    'Steam.Capacity.Factor.OBA','Exist.NGCC.Capacity.Factor.OBA',
                    'New.NGCC.Capacity.Factor.OBA','New.NGCC.Generation.Change',
                    'Emissions','Emissions.Change','Emissions.Diff')

# Remove any rows for regions where OBA is not applicable in this analysis
stats = stats[stats[,4]!=0 & stats[,5]!=0,]

# Write the results to file
write.csv(stats,file=paste0('oba approximation results ',year,'.csv'))