Postprocessor.py

def min_LC(write_flag, GT_flag):
	from os import getcwd
	from os import chdir
	from os import mkdir
	from csv import reader
	from math import exp
	from decimal import Decimal
	from Preprocessor import data_generator
	from Preprocessor import wind_generator
	from Preprocessor import wave_generator
	from Preprocessor import CF_natural_up
	from Preprocessor import CF_OCAES_up
	from Preprocessor import CF_interval_noOCAES
	from Preprocessor import CF_interval_OCAES
	
	dir_default	= getcwd()
	
	data = data_generator()
	wind = wind_generator()
	wave = wave_generator()

	Xw1 	= data['Xw1']
	etaT1 	= data['etaT1']
	C1 		= data['C1']
	F1 		= data['F1']
	V1 		= data['V1']
	Xw2 	= data['Xw2']
	etaT2 	= data['etaT2']
	C2 		= data['C2']
	F2 		= data['F2']
	V2 		= data['V2']
	Xc 		= data['Xc']
	Xe 		= data['Xe']
	etaC 	= data['etaC']
	etaE 	= data['etaE']
	Cl		= data['Cl']/1.0
	Cs	 	= data['Cs']/1.0
	Xt0 	= data['Xt0']
	etaT0 	= data['etaT0']
	T 		= data['T']
	r	 	= data['r']
	L	 	= data['L']
	Ls	 	= data['Ls']
	s0 		= data['s0']

	CCR		= r*(1 + r)**L/((1 + r)**L - 1)
	T		= int(T)
	Xt0		= int(Xt0)
	steps	= int(200/Xt0)
		
	w = list()	
	for i in range(0, T):
		w.append(WindMode*wind[i] + WaveMode*wave[i])

	# The CF limit when no OCAES incorporated
	CF_up_1					= CF_natural_up(w, Xt0)
	SCF_interval_noOCAES	= CF_interval_noOCAES(CF_up_1)

	# The system feasible CF when OCAES is involved
	CF_up_2					= CF_OCAES_up(w, Xt0, etaC, etaE)
	SCF_interval 			= CF_interval_OCAES(CF_up_1, CF_up_2)

	M = float('inf')

	################################################################################
	# Input result from result file generated by optimization script
	if GT_flag is True:
		dir_result = 'Result_' + str(T) + 'h' + '_' + str(steps) + 'x' + str(Xt0) + 'MW' + '_GT'
		
		chdir(dir_result)
		# Initialize the minimum container.
		# [n0, Target CF, CF, Xgt, LC ($/MWh), Time (s)]
		# Attention: the key to minimum is str(CF x 100), not the exact str(CF).
		minimum = {str(i): [M, M, M, M, M, M] for i in range(CF_up_1[str(steps - 1)], CF_up_2['1'] + 1)}
		result = dict()
		
		# Read optimization results from result log files
		for n0 in range(1, steps):
			csvname = 'n0_' + str(n0) + '_' + str(T) + 'h' + '.csv'
			with open(csvname, 'rb') as csvfile:
				csv_reader = reader(csvfile, dialect = 'excel')
				for row in csv_reader:
					try:
						result[(int(row[0]), int(100*Decimal(row[1])))] = [float(i) for i in row]
					except ValueError:
						continue

		# The minimum levelized cost when there is no OCAES involved
		for CF in [0.01*i for i in range(CF_up_1[str(steps - 1)], CF_up_1['1'] + 1)]:
			LC = M
			key = Decimal(str(CF))
			key = str(int(key*100))
			minimum[key] = [0, CF, CF, 0, M, M]
			for n0 in range(1, steps):
				if CF not in SCF_interval_noOCAES[str(n0)]:
					break
				C0 = 16.59*(1.971E6 + 0.209E6*exp(1.66*n0*Xt0*1E6/1E8))*1.6*0.15 +\
					 16.59*2400*1600*0.15 +\
					 0.03327*(n0*Xt0)**0.7513*1E6*1.35
				AnnualCapital		= CCR*(C0 + C1 + C2)
				AnnualFixed			= F1 + F2
				# AnnualVariable		= 8760.0/T*(V1*sum(wind[0: T]) + V2*sum(wave[0: T]))
				AnnualVariable		= 0
				AnnualProduced		= 8760.0/T*CF*n0*Xt0*T
				LC = (AnnualCapital + AnnualFixed + AnnualVariable)/AnnualProduced
				if minimum[key][4] >= LC:
					minimum[key] = [n0, CF, CF, 0, LC, M]
		
		chdir(dir_default)

		# The minimum levelized cost over the whole interval we are interested
		for key in result:
			if result[key][4] < minimum[str(key[1])][4]:
				minimum[str(key[1])] = result[key]
			else:
				continue
	else:
		dir_result = 'Result_' + str(T) + 'h' + '_' + str(steps) + 'x' + str(Xt0) + 'MW'
	
		chdir(dir_result)
		# Initialize the minimum container.
		# [n0, Target CF, CF, nl, Xs, LC ($/MWh), Time (s)]
		# Attention: the key to minimum is str(CF x 100), not the exact str(CF).
		minimum = {str(i): [M, M, M, M, M, M, M] for i in range(CF_up_1[str(steps - 1)], CF_up_2['1'] + 1)}
		result = dict()
		
		# Read optimization results from result log files
		for n0 in range(1, steps):
			csvname = 'n0_' + str(n0) + '_' + str(T) + 'h' + '.csv'
			with open(csvname, 'rb') as csvfile:
				csv_reader = reader(csvfile, dialect = 'excel')
				for row in csv_reader:
					try:
						result[(int(row[0]), int(100*Decimal(row[1])))] = [float(i) for i in row]
					except ValueError:
						continue

		# The minimum levelized cost when there is no OCAES involved
		for CF in [0.01*i for i in range(CF_up_1[str(steps - 1)], CF_up_1['1'] + 1)]:
			LC = M
			key = Decimal(str(CF))
			key = str(int(key*100))
			minimum[key] = [0, CF, CF, 0, 0, M, M]
			for n0 in range(1, steps):
				if CF not in SCF_interval_noOCAES[str(n0)]:
					break
				C0 = 16.59*(1.971E6 + 0.209E6*exp(1.66*n0*Xt0*1E6/1E8))*1.6*0.15 +\
					 16.59*2400*1600*0.15 +\
					 0.03327*(n0*Xt0)**0.7513*1E6*1.35
				AnnualCapital		= CCR*(C0 + C1 + C2)
				AnnualFixed			= F1 + F2
				# AnnualVariable		= 8760.0/T*(V1*sum(wind[0: T]) + V2*sum(wave[0: T]))
				AnnualVariable		= 0
				AnnualProduced		= 8760.0/T*CF*n0*Xt0*T
				LC = (AnnualCapital + AnnualFixed + AnnualVariable)/AnnualProduced
				if minimum[key][5] >= LC:
					minimum[key] = [n0, CF, CF, 0, 0, LC, M]
		
		chdir(dir_default)

		# The minimum levelized cost over the whole interval we are interested
		for key in result:
			if result[key][5] < minimum[str(key[1])][5]:
				minimum[str(key[1])] = result[key]
			else:
				continue
	# II()
	if write_flag is True:
		from csv import writer
		if GT_flag is True:
			csvname = 'minimum_' + str(T) + '_' + str(steps) + 'x' + str(Xt0) + 'MW_GT.csv'
		else:
			csvname = 'minimum_' + str(T) + '_' + str(steps) + 'x' + str(Xt0) + 'MW.csv'
		
		with open(csvname, 'wb') as f:
			csv_writer = writer(f)
			for CF in [0.01*i for i in range(CF_up_1[str(steps - 1)], CF_up_2['1'] + 1)]:
				key = Decimal(str(CF))
				key = str(int(key*100))
				csv_writer.writerow(minimum[key])
	
	return minimum

# Return the hours of simultaneous store and dispatch in a list.
def hrs_simultaneity(list_1, list_2, tolerance):
	hours = list()
	for i in range(0, len(list_1)):
		if abs(list_1[i]) <= tolerance or abs(list_2[i]) <= tolerance:
			continue
		if list_1[i]*list_2[i] == 0:
			continue
		hours.append(i + 1)
	return hours

# Return a dictionary contains list of each column of the result file.
# May not work for the n0xXt0 folder.
def read_result(T, n0, SCF, GT_flag):
	from os import getcwd
	from os import chdir
	from csv import reader
	
	dir_default	= getcwd()
	csvname		= 'result' + '_' + str(n0) + '_' + str(int(100*SCF)) + '.csv'
	result		= dict()
	t			= list()
	w			= list()
	w_s			= list()
	drop		= list()
	
	if GT_flag is True:
		GT_output	= list()
		dir_result	= 'Result_' + str(T) + 'h' + '_' + str(steps) + 'x' + str(Xt0) + 'MW' + '_GT'
		chdir(dir_result)
		try:
			with open(csvname, 'rb') as csvfile:
				csv_reader = reader(csvfile, dialect = 'excel')
				for row in csv_reader:
					try:
						t.append(int(row[0]))
						w.append(float(row[1]))
						w_s.append(float(row[2]))
						drop.append(float(row[3]))
						GT_output.append(float(row[4]))
					except ValueError:
						continue
				result['t']			= t
				result['w']			= w
				result['w_s']		= w_s
				result['drop']		= drop
				result['GT_output']	= GT_output
		except IOError:
			print '%s\t%s' % (csvname, ' does not exist.')
			chdir(dir_default)
			return None
	else:
		dispatch	= list()
		store		= list()
		dir_result	= 'Result_' + str(T) + 'h' + '_' + str(steps) + 'x' + str(Xt0) + 'MW'
		chdir(dir_result)
		try:
			with open(csvname, 'rb') as csvfile:
				csv_reader = reader(csvfile, dialect = 'excel')
				for row in csv_reader:
					try:
						t.append(int(row[0]))
						w.append(float(row[1]))
						w_s.append(float(row[2]))
						store.append(float(row[3]))
						dispatch.append(float(row[4]))
						drop.append(float(row[5]))
					except ValueError:
						continue
				result['t']			= t
				result['w']			= w
				result['w_s']		= w_s
				result['store']		= store
				result['dispatch']	= dispatch
				result['drop']		= drop
		except IOError:
			print '%s\t%s' % (csvname, ' does not exist.')
			chdir(dir_default)
			return None
	chdir(dir_default)
	return result

#This function is used by OCAES scenario only.
def equivalency(result, etaC, etaE):
	from copy import copy
	e_result	= copy(result)
	x			= copy(result['store'])
	y			= copy(result['dispatch'])
	z			= copy(result['drop'])
	for i in range(0, len(x)):
		if etaC*etaE*x[i] >= y[i]:
			e_result['store'][i]	= x[i] - y[i]/(etaC*etaE)
			e_result['dispatch'][i]	= 0
			e_result['drop'][i]		= z[i] - y[i] + y[i]/(etaC*etaE)
		else:
			e_result['store'][i]	= 0
			e_result['dispatch'][i]	= y[i] - etaC*etaE*x[i]
			e_result['drop'][i]		= z[i] + x[i] - etaC*etaE*x[i]
	return e_result

# If there is simultaneous store/dispatch, print the result file's name.
# ATTENTION: may not work for the gas turbine scenario and the 'n0xXt0' folder.
def scan_for_simul():
	from csv import writer
	from Preprocessor import data_generator
	from Preprocessor import wind_generator
	from Preprocessor import wave_generator
	from Preprocessor import CF_natural_up
	from Preprocessor import CF_OCAES_up
	from Preprocessor import CF_interval_noOCAES
	from Preprocessor import CF_interval_OCAES
	
	wind = wind_generator()
	wave = wave_generator()
	data = data_generator()
	
	T		= data['T']
	Xt0		= data['Xt0']
	etaC	= data['etaC']
	etaE	= data['etaE']
	T		= int(T)
	
	# WindMode = True
	# WaveMode = False
	
	w = list()	
	for i in range(0, T):
		w.append(wind[i] + wave[i])
	
	CF_up_1			= CF_natural_up(w, Xt0)
	CF_up_2			= CF_OCAES_up(w, Xt0, etaC, etaE)
	SCF_interval	= CF_interval_OCAES(CF_up_1, CF_up_2)
	
	tolerance	= 1.0E-5
	
	for n0 in range(1, 10):
		for SCF in SCF_interval[str(n0)]:
			result	= read_result(n0, SCF)
			if result == None:
				continue
			else:
				simu_sto_dis	= hrs_simultaneity(result['store'], result['dispatch'], tolerance)
				simu_sto_drp	= hrs_simultaneity(result['store'], result['drop'], tolerance)
				simu_dis_drp	= hrs_simultaneity(result['dispatch'], result['drop'], tolerance)
				if len(simu_sto_dis) != 0:
					csvname		= 'result' + '_' + str(n0) + '_' + str(int(100*SCF)) + '.csv'
					print '%s%s%i' % (csvname, ': ', len(simu_sto_dis))
	
				# e_result	= equivalency(result, etaC, etaE)
				# csvname		= 'e_result' + '_' + str(n0) + '_' + str(int(100*SCF)) + '.csv'
				# with open(csvname, 'wb') as f:
					# csv_writer = writer(f)
					# zipped = zip(e_result['t'], e_result['w'], e_result['w_s'], e_result['store'], e_result['dispatch'], e_result['drop'])
					# zipped[:0] = [('Time', 'Input Power', 'Power to grid', 'Store', 'Dispatch', 'Drop')]
					# csv_writer.writerows(zipped)

def plot_minimum(minimum, GT_flag):
	import matplotlib.pyplot as plt
	if GT_flag is True:
		x		= [int(i) for i in minimum.keys()]
		x.sort()
		y_LC	= [minimum[str(i)][4] for i in x]
		y_Xgt	= [minimum[str(i)][3] for i in x]
		
		figure	= plt.figure()
		f_LC	= plt.figure(1)
		plt.plot(x, y_LC, '-*b')
		plt.xlabel('Capacity Factor (%)')
		plt.ylabel('Levelized cost ($/MWh)')
		plt.title('Levelized cost')

		f_nl	= plt.figure(2)
		plt.plot(x, y_Xgt, '-b*')
		plt.xlabel('Capacity Factor (%)')
		plt.ylabel('Capacity of GT (MW)')
		plt.title('Capacity of Gas Turbine')
	else:
		x		= [int(i) for i in minimum.keys()]
		x.sort()
		y_LC	= [minimum[str(i)][5] for i in x]
		y_nl	= [minimum[str(i)][3] for i in x]
		y_Xs	= [minimum[str(i)][4] for i in x]
		
		figure	= plt.figure()
		f_LC	= plt.figure(1)
		plt.plot(x, y_LC, '-*b')
		plt.xlabel('Capacity Factor (%)')
		plt.ylabel('Levelized cost ($/MWh)')
		plt.title('Levelized cost')

		f_nl	= plt.figure(2)
		plt.plot(x, y_nl, '-b*')
		plt.xlabel('Capacity Factor (%)')
		plt.ylabel('# of LP')
		plt.title('# of Liquid Pistons')

		f_Xs	= plt.figure(3)
		plt.plot(x, y_Xs, '-b*')
		plt.xlabel('Capacity Factor (%)')
		plt.ylabel('Capacity (MWh)')
		plt.title('Capacity of storage')
	plt.show()

def plot_all():
	from os import getcwd
	from os import chdir
	from csv import reader
	from decimal import Decimal
	from Preprocessor import data_generator
	import matplotlib.pyplot as plt
	
	data	= data_generator()
	T		= data['T']
	Xt0		= data['Xt0']
	
	T		= int(T)
	Xt0		= int(Xt0)
	steps	= int(200/Xt0)
	
	dir_default	= getcwd()
	dir_OCAES	= '.\\Result_' + str(T) + 'h' + '_' + str(steps) + 'x' + str(Xt0) + 'MW'
	dir_GT		= '.\\Result_' + str(T) + 'h' + '_' + str(steps) + 'x' + str(Xt0) + 'MW' + '_GT'
	
	CF			= list()
	# n0_OCAES	= list()
	# n0_GT		= list()
	LC_OCAES	= list()
	LC_GT		= list()
	
	chdir(dir_OCAES)
	for n0 in range(1, steps):
		CF.append([])
		LC_OCAES.append([])
		csvname = 'n0_' + str(n0) + '_' + str(T) + 'h' + '.csv'
		with open(csvname, 'rb') as f:
			csv_reader	= reader(f, dialect = 'excel')
			for row in csv_reader:
				try:
					CF[n0 - 1].append(float(row[1]))
					LC_OCAES[n0 - 1].append(float(row[5]))
				except ValueError:
					continue
	chdir(dir_default)
	
	chdir(dir_GT)
	for n0 in range(1, steps):
		LC_GT.append([])
		csvname = 'n0_' + str(n0) + '_' + str(T) + 'h' + '.csv'
		with open(csvname, 'rb') as f:
			csv_reader	= reader(f, dialect = 'excel')
			for row in csv_reader:
				try:
					LC_GT[n0 - 1].append(float(row[4]))
				except ValueError:
					continue
	chdir(dir_default)
	
	figure1	= plt.figure('OCAES')
	figure1.hold(True)
	figure2	= plt.figure('GT')
	figure2.hold(True)
	figure3	= plt.figure('Comparison')
	figure3.hold(True)
	for n0 in range(1, steps):
		plt.figure('OCAES')
		plt.plot(CF[n0 - 1], LC_OCAES[n0 - 1], '-b*')
		plt.xlabel('CF')
		plt.ylabel('LC ($/MWh)')
		
		plt.figure('GT')
		plt.plot(CF[n0 - 1], LC_GT[n0 - 1], '-ro')
		plt.xlabel('CF')
		plt.ylabel('LC ($/MWh)')
		
		plt.figure('Comparison')
		plt.plot(CF[n0 - 1], LC_OCAES[n0 - 1], '-b*', CF[n0 - 1], LC_GT[n0 - 1], '-ro')
		plt.xlabel('CF')
		plt.ylabel('LC ($/MWh)')		
	plt.show()

if __name__ == '__main__':
	from IPython import embed as II
	from csv import writer
	from sys import argv
	
	WindMode	= True
	WaveMode	= False
	
	if argv[1] == 'GT':
		GT_flag	= True
	else:
		GT_flag = False
		
	if argv[2] == 'W':
		write_flag = True # If true, minimum value will be output as csv file
	else:
		write_flag = False
	####################################################################
	minimum = min_LC(write_flag, GT_flag)
	plot_minimum(minimum, GT_flag)
	####################################################################
	
	# scan_for_simul()
	
	####################################################################
	# plot_all()
	####################################################################
	# while True:
		# try:
			# T			= input('Input T: ')
			# n0			= input('Input n0: ')
			# SCF			= input('Input CF: ')
			# tolerance	= input('Input tolerance: ')
			# result		= read_result(T, n0, SCF, GT_flag)
			# # print hrs_simultaneity(result['store'], result['dispatch'], tolerance)
			# hours		= hrs_simultaneity(result['drop'], result['GT_output'], tolerance)
			# for i in hours:
				# print result['w_s'][i - 1]
		# except TypeError:
			# break
	
	# II()