# 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.gener

ation=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+un

its$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(da

ta[[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/potentia

l)*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]][,'emissi

ons.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'))