Sgouridis, Carbajales-Dale, Csala, Chiesa, Bardi
Use this Jupyter workbook to reproduce all results presented in the paper.
This is document has been created using Jupyter Notebook in the Anaconda distribution and it can be edited and run in active mode by clicking download in top right corner of this page. The code is partitioned into sections, called cells. When you are using this workbook in active mode, double-click on a cell to edit it and then run using Ctrl + Enter
. Hitting Shift + Enter
runs the code and steps into the next cell, while Alt + Enter
runs the code and adds a new, empty cell. If you are running this notebook on a presonal computer, you will need a machine with at least 1GB of memory (2GB recommended) and a processor of 1GHz.
Import dependencies
import pandas as pd, matplotlib.pyplot as plt, numpy as np
from matplotlib.colors import BoundaryNorm
from matplotlib.ticker import MaxNLocator
import matplotlib.patches as patches
import matplotlib.lines as lines
%matplotlib inline
import sys, string
from sympy import *
init_printing()
print('Running on Python',sys.version)
Running on Python 3.6.6 |Anaconda, Inc.| (default, Jun 28 2018, 11:27:44) [MSC v.1900 64 bit (AMD64)]
This section reproduces all the equations presented in the paper (supplementary material). This is just for easier visual guidance and for quick checking.
var('EROEI E_out E_in');
eq_S1=Eq(EROEI,E_out/E_in)
eq_S1
var('E_cap E_f E_OM');
eq_S1i=Eq(EROEI,E_out/(E_cap+E_OM+E_f))
eq_S1i
var('L s_OM');
eq_S2i=Eq(E_in,E_cap*(1+L*s_OM)+E_f)
eq_S2i
var('P L cf');
eq_S2ii=Eq(E_out,P*cf*L)
eq_S2ii
var('E_el');
eq_S2iii=Eq(E_out,E_el)
eq_S2iii
e=Symbol('eta')
var('E_th');
eq_S2iv=Eq(E_th,P*cf*L/e)
eq_S2iv
var('EROEI_th');
eq_S3i=Eq(E_f,E_th/EROEI_th)
eq_S3i
#define pretty variable substituion function
def sub(var,eq):
return [var,solve(eq,var)[0]]
#create multiple substituion method by chaining
def multisub(eq0,msub):
for i in range(len(msub)//2):
eq0=eq0.subs([sub(msub[i*2],msub[i*2+1])])
return eq0
var('EROEI_el');
eq_S3=Eq(EROEI_el,multisub(eq_S1,[E_out,eq_S2ii,E_in,eq_S2i,E_f,eq_S3i,E_th,eq_S2iv]).rhs)
eq_S3
var('CR f_op EROEI_CCS f_cap');
eq_S4=Eq(EROEI_CCS,E_el*(1-f_op(CR))/(E_cap*(1+f_cap)*(1+L*s_OM)+E_f))
eq_S4
var('CR f_op EROEI_CCS f_cap');
eq_S5=Eq(EROEI_CCS,(1-f_op(CR))*(E_cap*(1+L*s_OM)+E_f)*EROEI_el/(E_cap*(1+f_cap)*(1+L*s_OM)+E_f))
eq_S5
var('R');
eq_S6i=Eq(R,E_f/(E_cap*(1+L*s_OM)))
eq_S6i
eq_S6=Eq(EROEI_CCS,(1-f_op(CR))*(R+1)*EROEI_el/(R+1+f_cap))
eq_S6
We have made a major simplifaction step from Eq. S5 to reach Eq. S6. Here we prove by symbolic evaluation of the subtituted expressions that we were correct.
eq_S5.simplify()
eq_S6.subs(R,eq_S6i.rhs).simplify()
#eq_S6 is correct!
eq_S6.subs(R,eq_S6i.rhs).simplify()==eq_S5.simplify()
True
In this part we show you how to reproduce the graphs presented in the paper, which form the basis of our results discussion and conclusions. There is a similar section for the Renewables part. First, we illustrate a representative range for CCS operational parameters for a number of fossil fuel plants. We then take note of the min and max ranges and plot the values alongside those obtained from the Renewables part. The calculations are based on actual power plant data, stored in the Excel file plants.xlsx
, found in the project folder and presented in Table S1 in the supplementary material of the paper.
Using the formula in Eq.4 and the plant values from the Table S1 (reproduced from plants.xlsx
) calculate the $EROEI_{el}$ for all capacity factors $CF$s between 50 and 90. We then use these values as a proxy to create an a adjustment factor for the $EROEI_{el}$ due to the indirect effect of the change of the capacity factor. The same is true for the operational versus capital energy investment ratio, $R$. We therefore created capacity-factor adjusted values for $R$ as well.
In order to manipulate the Excel file on the fly, we use a python module called xlwings
. This will open an Excel application in the background. Be sure to either manually close it after you're done or run the cell just above the Auxiliary plots section.
import xlwings as xw #will open a MS Excel application, let it run in the background
wb = xw.Book('plants.xlsx')
We define a function eroei_r_getter
which takes as parameters the capacity factor $CF$ and the $EROEI_{fuel}$ for the two studies fossil fuels, Coal and Natural Gas. Then, this function sets the correct $CF$ and $EROEI_{fuel}$ values in the Excel workbook, which then recalculates the resulting $EROEI_{el}$ and $R$ values. Since we have 3 IGCC Coal, 2 Pulverized Coal and 1 NGCC Gas plant, we average over $EROEI_{el}$ and $R$ values within each group. The resulting pandas dataframe is returned by the function, which now contains 3 columns ['IGCC','PC','NGCC']
and 2 rows ['EROEI_el','R']
.
def eroei_r_getter(CF=100,EROEI_fuel_coal=58,EROEI_fuel_gas=87):
xw.Range('L7').value=CF/100.0 #hard-coded Excel workbook cell ID for capacity factor
xw.Range('M8').value=EROEI_fuel_coal #hard-coded Excel workbook cell ID for EROEI_fuel for Coal
xw.Range('O8').value=EROEI_fuel_gas #hard-coded Excel workbook cell ID for EROEI_fuel for Gas
wb.save() #write new values
df=pd.read_excel('plants.xlsx',usecols=range(10,24)).loc[[37,38]]
df=df[df.columns[1::2]]
df.columns=['IGCC 1','IGCC 3','IGCC 5','PC 9','PC 11','NGCC 13']
df.index=['EROEI_el','R']
df['IGCC']=(df['IGCC 1']+df['IGCC 3']+df['IGCC 5'])/3.0
df['PC']=(df['PC 9']+df['PC 11'])/2.0
df['NGCC']=df['NGCC 13']
return df[['IGCC','PC','NGCC']]
eroei_r_getter(95,58,87) #example to get the EROEI_el and R for capacity factor 95% and EROEI_fuels of 58 for Coal and 87 for Gas
IGCC | PC | NGCC | |
---|---|---|---|
EROEI_el | 12.0177 | 14.0583 | 31.9783 |
R | 1.5745 | 2.06751 | 2.91996 |
Create reasonable operating range for capacity factor $CF$ (%).
CFs=[i*5 for i in range (10,21)]
CFs
This cell is just an illustration of the pseudo-linear relationship between $CF$ and $EROEI_{el}$ (top row) and $R$ (bottom row), respectively. Note that the relationship is also nonlinear with regards to $EROEI_{fuel}$ (left vs. right column).
EROEI_fuel=[[58,87],[23,20]]
dfe=[]
dfr=[]
for j in range(len(EROEI_fuel)):
EROEI_fuel_coal=EROEI_fuel[j][0]
EROEI_fuel_gas=EROEI_fuel[j][1]
dfs1=[]
dfs2=[]
for i in range(len(CFs)):
df=eroei_r_getter(CFs[i],EROEI_fuel_coal,EROEI_fuel_gas)
df1=df.loc['EROEI_el']
df2=df.loc['R']
df1.name=CFs[i]
df2.name=CFs[i]
dfs1.append(df1)
dfs2.append(df2)
dfe.append(pd.concat(dfs1,axis=1))
dfr.append(pd.concat(dfs2,axis=1))
Figure S3. presents two $EROEI_{fuel}$ values for Coal and Gas (columns), respectively, and the dependence of $EROEI_{el}$ (top row) and the Fuel to capital ratio $R$ (bottom row) on the Capacity Factor $CF$.
#Plots init
fig,axes=plt.subplots(2,2,figsize=(8,6),subplot_kw=dict(facecolor='#EEEEEE',axisbelow=True))
for axi in axes:
for ax in axi:
ax.grid(color='white', linestyle='solid')
plt.subplots_adjust(wspace=0.15)
plt.subplots_adjust(hspace=0.1)
for j in range(len(EROEI_fuel)):
EROEI_fuel_coal=EROEI_fuel[j][0]
EROEI_fuel_gas=EROEI_fuel[j][1]
df=dfe[j]
ax=axes[0][j]
df.T.plot(ax=ax,grid=True,legend=False)
ax.set_title('$EROEI_{fuel}$ Coal$='+str(EROEI_fuel_coal)+'$\n$EROEI_{fuel}$ Gas$='+str(EROEI_fuel_gas)+'$')
df=dfr[j]
ax=axes[1][j]
df.T.plot(ax=ax,grid=True,legend=False)
for axi in range(len(axes[j])):
ax=axes[j][axi]
ax.text(0.04, 0.95, string.ascii_lowercase[axi*len(axes)+j],
horizontalalignment='left',
verticalalignment='top',
transform=ax.transAxes,size=13,alpha=0.7,
bbox=dict(boxstyle="round", fc="w", alpha=0.7))
axes[1][1].set_xlabel('Capacity Factor $CF$ (%)',fontsize=13,x=0)
axes[0][0].set_ylabel('$EROEI_{el}$',fontsize=13)
axes[1][0].set_ylabel('Fuel to capital ratio $R$',fontsize=13)
#axes[0][1].set_yticklabels('')
#axes[1][1].set_yticklabels('')
axes[0][0].set_xticklabels('')
axes[0][1].set_xticklabels('')
axes[0][0].legend(framealpha=0)
#plt.suptitle('Normalized Relationship of $EROEI_{el}$ and $R$ to Capacity Factor $CF$ ',fontsize=13,y=1.01)
plt.savefig('plot/png/figS3.png',dpi=300,bbox_inches = 'tight', pad_inches = 0.1)
plt.savefig('plot/eps/figS3.eps',format='eps',bbox_inches = 'tight', pad_inches = 0.1)
plt.savefig('plot/pdf/figS3.pdf',format='pdf',bbox_inches = 'tight', pad_inches = 0.1)
plt.show()
We need to create a energy penalty conversion ratio for the Capture Ratio $CR$. We use the Energy penalty figure Fig. S1 from the supplementary material. We then use this value to adjust the value of the operational energy penalty $f_{op}$. All this means that a higher Capture Ratio $CR$ lowers $EROEI_{CCS}$ more - and we would like to create a measure about how much. We use the middle solid black line, linear approximation from the plot.
CR | energy penalty |
---|---|
60 | 7.8 |
90 | 11.8 |
We will only test Capture Ratios $CR$ of $60$ and $90$.
CR=[90,60]
CR_ratio=[1, 7.8/11.8]
CR_ratio
Create a simple fuel name convereter function.
def fuel_name_converter(fuel):
if fuel=='NGCC\nGas': return 'NGCC'
elif fuel=='IGCC\nCoal': return 'IGCC'
elif fuel=='Pulverized\nCoal': return 'PC'
else: return None
Define $EROEI_{CCS}$ calculator function eroei_ccs_cycler
, including minimums and maximums over $f_{op}$ and $f_{cap}$. This function takes as input:
names
, which holds the fossil fuel type namest0s
in the format [min $f_{op}$,max $f_{op}$, min $f_{cap}$,max $f_{cap}$], with this array repreated for each name
The function will iterate over all inputs and return 4 multidimensional variables. First, uisng the eroei_r_getter
, the function will calculate the corresponding $EROEI_{el}$ and $R$ to the $CF$ and $EROEI_{fuel}$ input. Then it will calculate $EROEI_{CCS}$ using Eq. S6, for all variable ranges. The t0x
inner loop keeps track of the minimums and maximums. Finally, it returns the calculated 2-dimensional arrays $EROEI_{el}$ and $R$ in eroei_el_cf
and R_cf
, and multidimensional arrays eroei_ccs
, eroei_ccs_min
, eroei_ccs_max
, holdig the actual values, the minimums and the maximums for $EROEI_{CCS}$.
def eroei_ccs_cycler(names, eroei_fuels, t0s, CF, CR_ratio, fop, fcap, verbose=False):
eroei_ccs=np.zeros([len(names),len(CF),len(CR_ratio),len(fop),len(fcap)])
eroei_ccs_min=np.zeros([len(names),len(CF),len(CR_ratio)])+1000
eroei_ccs_max=np.zeros([len(names),len(CF),len(CR_ratio)])
eroei_el_cf=np.zeros([len(names),len(CF)])
R_cf=np.zeros([len(names),len(CF)])
for i in range(len(names)):
if verbose: print('Calculating '+names[i].replace('\n',' ')+'...')
for j in range(len(CF)):
eroei_r=eroei_r_getter(CF[j],eroei_fuels[i],eroei_fuels[i])[fuel_name_converter(names[i])]
#convert fuel EROEI to EROEI_el using the CF converstion table
eroei_el_cf[i][j]=eroei_r.loc['EROEI_el']
#calculate R using the CF converstion table
R_cf[i][j]=eroei_r.loc['R']
for r in range(len(CR_ratio)):
for k in range(len(fop)):
for l in range(len(fcap)):
#calculate eoris over maps using eq_S6
eroei_ccs[i][j][r][l,k]=(1-fop[k]*CR_ratio[r]/100.0)*\
((R_cf[i][j]+1)/(R_cf[i][j]+1+fcap[l]/100.0))*\
eroei_el_cf[i][j]
#calculate min-maxes
for t0x in range(len(t0s[i])//2):
for t0y in range(len(t0s[i])//2,len(t0s[i])):
x=t0s[i][t0x]
y=t0s[i][t0y]
if ((abs(fop[k]-x)<0.5)and(abs(fcap[l]-y)<1.5)):
eroei_ccs_max[i][j][r]=max(eroei_ccs_max[i][j][r],eroei_ccs[i][j][r][l,k])
eroei_ccs_min[i][j][r]=min(eroei_ccs_min[i][j][r],eroei_ccs[i][j][r][l,k])
if verbose:
print('CF='+str(CF[j])+', CR='+str(CR[r])+\
', EROEI_fuel='+str(eroei_fuels[i])+\
', EROEI_el_CF_'+str(CF[j])+'%='+str(np.round(eroei_el_cf[i][j],2))+\
', R_CF_'+str(CF[j])+'%='+str(np.round(R_cf[i][j],2))+\
', EROEI_CCS - min='+str(np.round(eroei_ccs_min[i][j][r],2))+\
', max='+str(np.round(eroei_ccs_max[i][j][r],2)))
return eroei_el_cf, R_cf, eroei_ccs, eroei_ccs_min, eroei_ccs_max
Having set up the eroei_ccs_cycler
function, we define our desired simulation ranges and run it. We also display min and max values for the $EROEI_{CCS}$. These are the bottom left and top right corners of the rectangles in Figure 2.
#Simulation options
fop=np.linspace(5,40,100) #in %
fcap=np.linspace(10,120,100) #in %
CF=[55,85] #in %
title=u'$EROEI_{CCS}$ comparison for Coal and Gas'
#eroei_els=[30,13,13] #deprecated, varies by name
#eroei_els=[32.4,14.3,12.25] #deprecated, 100% CF value from the CF converstion table
eroei_fuels=[87,58,58]
names=['NGCC\nGas','Pulverized\nCoal','IGCC\nCoal']
names_=['NGCC Gas','Pulverized Coal','IGCC Coal']
t0s=[[12,20,58,110],[17,31,45,95],[15,28,23,40]] # [fop min, fop max, fcap min fcap max] by each name, in %
eroei_el_cf, R_cf, eroei_ccs, eroei_ccs_min, eroei_ccs_max=eroei_ccs_cycler(names, eroei_fuels, t0s, CF, CR_ratio, fop, fcap, verbose=True)
Calculating NGCC Gas... CF=55, CR=90, EROEI_fuel=87, EROEI_el_CF_55%=26.97, R_CF_55%=1.69, EROEI_CCS - min=15.23, max=19.67 CF=55, CR=60, EROEI_fuel=87, EROEI_el_CF_55%=26.97, R_CF_55%=1.69, EROEI_CCS - min=16.54, max=20.56 CF=85, CR=90, EROEI_fuel=87, EROEI_el_CF_85%=31.05, R_CF_85%=2.61, EROEI_CCS - min=18.95, max=23.69 CF=85, CR=60, EROEI_fuel=87, EROEI_el_CF_85%=31.05, R_CF_85%=2.61, EROEI_CCS - min=20.57, max=24.76 Calculating Pulverized Coal... CF=55, CR=90, EROEI_fuel=58, EROEI_el_CF_55%=11.36, R_CF_55%=1.2, EROEI_CCS - min=5.45, max=7.87 CF=55, CR=60, EROEI_fuel=58, EROEI_el_CF_55%=11.36, R_CF_55%=1.2, EROEI_CCS - min=6.29, max=8.41 CF=85, CR=90, EROEI_fuel=58, EROEI_el_CF_85%=13.54, R_CF_85%=1.85, EROEI_CCS - min=6.98, max=9.76 CF=85, CR=60, EROEI_fuel=58, EROEI_el_CF_85%=13.54, R_CF_85%=1.85, EROEI_CCS - min=8.05, max=10.42 Calculating IGCC Coal... CF=55, CR=90, EROEI_fuel=58, EROEI_el_CF_55%=9.37, R_CF_55%=0.91, EROEI_CCS - min=5.52, max=7.17 CF=55, CR=60, EROEI_fuel=58, EROEI_el_CF_55%=9.37, R_CF_55%=0.91, EROEI_CCS - min=6.26, max=7.58 CF=85, CR=90, EROEI_fuel=58, EROEI_el_CF_85%=11.49, R_CF_85%=1.41, EROEI_CCS - min=7.03, max=8.99 CF=85, CR=60, EROEI_fuel=58, EROEI_el_CF_85%=11.49, R_CF_85%=1.41, EROEI_CCS - min=7.98, max=9.51
Generate Figure 2. This is a 2 by 2 multiplot, with each sub-plot being a filled contour area plot, sometimes also referred to as a colormap or shaded area plot. Here we have plotted the 4-dimensional variable $EROEI_{CCS}$ in the following manner:
#Plotting options
fig,axes=plt.subplots(2,2,figsize=(10,8))
plt.subplots_adjust(hspace=0.1)
plt.subplots_adjust(wspace=0.25)
eroei_ccs_scale_min2=1.2
eroei_ccs_scale_max2=2.3
ls=['-','--']
ls0='-'
lw0=0.5
t_x_offset=[6.5,12.5,6]
t_y_offset=[47,45,14]
breaker=['\n','\n','\n']
manual_locations=[[[(12,20),(20,20),(30,21),(35,20),(35,60)],\
[(7,70),(32,20),(35,60)]],\
[[(15,20),(23,19),(30,21),(35,20),(35,60)],\
[(7,70),(35,30),(35,70)]]]
levelsi=[[np.arange(11,24,1),np.arange(4,10.5,0.5)],[np.arange(14,27,1),np.arange(6,12.5,0.5)]]
for i in range(len(axes)):
for j in range(len(axes[i])):
eroei_el=eroei_el_cf[j][i]
ax=axes[i][j]
X, Y = np.meshgrid(fop, fcap)
z = eroei_ccs[j][i][0]#*100.0/eroei_el #third index is CR, keep at 0, its relative
levels = levelsi[i][j]
cmap = plt.get_cmap('coolwarm_r')
im = ax.contourf(X, Y, z, cmap=cmap, levels=levels,extend='both')
im2 = ax.contour(X, Y, z, levels=levels,extend='both',colors='k',linewidths=0.4)
ax.clabel(im2, fmt = '%.i',inline=1, fontsize=11, manual=manual_locations[i][j])
ax.set_xlim((fop.min(),fop.max()))
ax.set_ylim((fcap.min(),fcap.max()))
#fig.subplots_adjust(right=0.9)
cbar_ax = fig.add_axes([0.48+0.43*j, 0.12+0.395*(1-i), 0.015, 0.365])
fig.colorbar(im, cax=cbar_ax)
if j==1: cbar_ax.set_ylabel(u'$EROEI_{CCS}$',fontsize=14,labelpad=10)
for m in range(len(CR_ratio)):
n=0
for ti in range(len(t0s)):
if (((j==1)and(ti==2))or(ti==j)):
t0=t0s[ti]
t=[t0[0]*CR_ratio[m],t0[1]*CR_ratio[m],t0[2]+n,t0[3]+n]
rect = patches.Rectangle((t[0],t[2]),t[1]-t[0],t[3]-t[2],ls=ls[m],
edgecolor='k',facecolor='none')
ax.add_patch(rect)
'''
if m==1:
ax.text(t[0]+2,t[3]-14,'CR\n'+str(CR[m]),horizontalalignment='center')
else:
ax.text(t[0]+t_x_offset[ti],t[3]-t_y_offset[ti],'CR'+breaker[ti]+str(CR[m]),
horizontalalignment='center')
'''
if m==0:
ax.text(t[1]-0.5,t[3]-15,names[ti],horizontalalignment='right',weight='bold')
if (ti==1)and(m==1):
line = lines.Line2D((t0[0],t[1]),((t0[2]+t0[3])/2,(t[2]+t[3])/2),color='k',
linewidth=lw0,ls=ls0)
#ax.add_line(line)
'''
if j==0:
ax.set_title(names[0].replace('\n',' ')+' - $EROEI_{el}='+\
str(round(eroei_el_cf[0][i],1))+\
u'$\nCapacity factor $CF='+str(CF[i])+'$%')
else:
ax.set_title(names[1].replace('\n',' ')+' - $EROEI_{el}='+\
str(round(eroei_el_cf[1][i],1))+'$\n'+\
names[2].replace('\n',' ')+' - $EROEI_{el}='+\
str(round(eroei_el_cf[2][i],1))+\
u'$\nCapacity factor $CF='+str(CF[i])+'$%')
'''
ax.text(0.93, 0.95, string.ascii_lowercase[j*len(axes)+i],
horizontalalignment='left',
verticalalignment='top',
transform=ax.transAxes,size=13,alpha=0.7,
bbox=dict(boxstyle="round", fc="w", alpha=0.7))
axes[1][1].set_xlabel('Operational energy penalty $f_{op}$ (%)',fontsize=13,x=0)
axes[0][0].set_ylabel('Capital energy penalty $f_{cap}$ (%)',fontsize=13,y=0)
axes[0][1].set_yticklabels('')
axes[1][1].set_yticklabels('')
axes[0][0].set_xticklabels('')
axes[0][1].set_xticklabels('')
#plt.suptitle(title+'\nframes represent different capture ratios CR (%)',fontsize=13,y=1.04)
plt.savefig('plot/png/fig2.png',dpi=300,bbox_inches = 'tight', pad_inches = 0.1)
plt.savefig('plot/eps/fig2.eps',format='eps',bbox_inches = 'tight', pad_inches = 0.1)
plt.savefig('plot/pdf/fig2.pdf',format='pdf',bbox_inches = 'tight', pad_inches = 0.1)
plt.show()
l=Symbol('lambda')
epse=Symbol('epsilon')
var('ESOI D');
eq_S7=Eq(ESOI,l*e*D/epse)
eq_S7
epsp=Symbol('epsilon_p')
epss=Symbol('epsilon_s')
var('D P C');
eq_S8=Eq(ESOI,C*l*e*D/(P*epsp+C*epss))
eq_S8
f=Symbol('varphi_j')
Ei=Symbol('E_{in_i}')
Ej=Symbol('E_{in_j}')
var('EROEI_disp c j i I J');
eq_S9=Eq(EROEI_disp,(E_out-(1-e)*Sum(f*E_out,(j,1,J))-c*E_out)/(Sum(Ei,(i,1,I))+Sum(Ej,(j,1,J))))
eq_S9
a=Symbol('alpha_i')
ej=Symbol('eta_j')
f=Symbol('varphi_j')
Ei=Symbol('E_{in_i}')
Ej=Symbol('E_{in_j}')
var('EROEI_i u_j ESOI_j');
eq_S10=Eq(EROEI_disp,(E_out*(1-(1-e)*Sum(f,(j,1,J))-c))/(Sum(E_out*a/EROEI_i,(i,1,I))+Sum((E_out*f*ej)/(ESOI_j*u_j),(j,1,J))))
eq_S10
eq_S11=Eq(EROEI_disp,((1-(1-e)*Sum(f,(j,1,J))-c))/(Sum(a/EROEI_i,(i,1,I))+Sum((f*ej)/(ESOI_j*u_j),(j,1,J))))
eq_S11
Define $EROEI_{disp}$ calcuator function eroei_disp_cycler
, including minimums and maximums over Storage fraction $\varphi$. It takes as inputs:
ESOI
, of length how many storage scenarios we are considering, holding $ESOI$ values for eachphi
, of arbitrary length, with values between $0$ and $1$, corresponding to the simulated storage fraction levelseta
, holding the Roundtrip efficiency for each storage scenarioc
, holding the Curtailment for each storage scenarioeroei_el
, of arbitrary length, with values between $1$ and $100$, corresponding to the $EROEI_{el}$ levelsReturn one multidimensional array eroei_disp
, corresponding to the Dispatchable Scalable Renewable EROEI $EROEI_{disp}$ for each parameter combination.
def eroei_disp_cycler(ESOI, phi, eta, c, eroei_el):
#eroei_disp=np.zeros([len(ESOI),len(phi),len(eroei_el)])
eroei_disp=np.zeros([len(ESOI),len(eroei_el),len(phi)])
for i in range(len(ESOI)):
for j in range(len(phi)):
for k in range(len(eroei_el)):
#eq_S11
#Also Eq. 6 from http://pubs.rsc.org/en/content/articlepdf/2013/ee/c3ee41973h
eroei_disp[i][k,j]=(-c[i]-(+1-eta[i])*phi[j]+1)/\
((1.0/eroei_el[k])+(eta[i]*phi[j]/ESOI[i]))
return eroei_disp
Generate Figure 3. This is a 2 by 2 multiplot, with each sub-plot being a filled contour area plot, sometimes also referred to as a colormap or shaded area plot. Here we have plotted the 4-dimensional variable $EROEI_{disp}$ in the following manner:
segments=800
eroei_el=np.linspace(10,60,segments)
phi=np.linspace(0,35,segments)
ESOI=[11,24,11,249]
ESOI_tech=['Battery storage','Hydrogen storage','Battery + Curtailment','Pumped hydro storage']
eta=np.array([83,30,83,80])
c=[0,0,0.07,0]
eroei_disp=eroei_disp_cycler(ESOI, phi/100.0, eta/100.0, c, eroei_el)
fig,axes=plt.subplots(2,2,figsize=(10,8))
plt.subplots_adjust(wspace=0.1)
plt.subplots_adjust(hspace=0.15)
contour_min=10
contour_max=60
levels = MaxNLocator(nbins=10).tick_values(contour_min, contour_max)
manual_locations=[[(30,10),(7,20),(7,34),(7,50),(3,55)],\
[(30,10),(7,20),(7,30),(7,45),(7,55)],\
[(30,10),(7,25),(7,42),(7,55)],\
[(7,20),(7,30),(7,40),(7,50)]]
for i in range(len(ESOI)):
ax=axes[i//2][i%2]
z = eroei_disp[i]
cmap = plt.get_cmap('coolwarm_r')
X, Y = np.meshgrid(phi, eroei_el)
im = ax.contourf(X, Y, z, cmap=cmap, levels=levels,extend='both')
im2 = ax.contour(X, Y, z, levels=levels,extend='both',colors='k',linewidths=0.4)
ax.clabel(im2, fmt = '%i',inline=1, fontsize=11, manual=manual_locations[i])
if (i==0):
fig.subplots_adjust(right=0.9)
cbar_ax = fig.add_axes([0.93, 0.12, 0.015, 0.77])
fig.colorbar(im, cax=cbar_ax)
cbar_ax.set_ylabel(u'$EROEI_{disp}$',fontsize=13,labelpad=10)
ax.set_xlim((phi.min(),phi.max()))
ax.set_ylim((eroei_el.min(),eroei_el.max()))
ax.set_title(ESOI_tech[i],fontsize=13)
ax.text(0.93, 0.95, string.ascii_lowercase[(i%2)*2+i//2],
horizontalalignment='left',
verticalalignment='top',
transform=ax.transAxes,size=13,alpha=0.7,
bbox=dict(boxstyle="round", fc="w", alpha=0.7))
axes[1][1].set_xlabel('Storage fraction $\\varphi$ (%)',fontsize=13,x=0)
axes[0][0].set_ylabel('$EROEI_{el}$',fontsize=13,y=0)
axes[0][1].set_yticklabels('')
axes[1][1].set_yticklabels('')
axes[0][0].set_xticklabels('')
axes[0][1].set_xticklabels('')
#plt.suptitle(r'$EROEI_{disp}$ for Scalable Dispatchable Renewables - PV and Wind',fontsize=13,y=0.98)
plt.savefig('plot/png/fig3.png',dpi=150,bbox_inches = 'tight', pad_inches = 0.1)
plt.savefig('plot/eps/fig3.eps',format='eps',bbox_inches = 'tight', pad_inches = 0.1)
plt.savefig('plot/pdf/fig3.pdf',format='pdf',bbox_inches = 'tight', pad_inches = 0.1)
plt.show()
In this section we conduct a comparative sensitivity analysis over a range of parameters presented in the Fossils and Renewables sections, in order to be able to represent $EROEI_{CCS}$ and $EROEI_{disp}$ in one common axis, for a range of $EROEI_{el}$ values. In essence this means replicating the rectangles of Figure 2. subplots and taking their minimum and maximum ranges, over many starting $EROEI_{fuel}$ - and thus varying $EROEI_{el}$ and $R$ - values, to result in trapezoids. This is then compared with the 4 energy storage cases presented in Figure 3., with separate $EROEI_{el}$ range assumptions for Wind and PV.
First, we define a convex hull function to get outermost points of trapezoids in the sensitivity analysis area.
from scipy.spatial import ConvexHull
def convexhull(p):
p = np.array(p)
hull = ConvexHull(p)
return p[hull.vertices,:]
Define wieldy color functions.
def cbrew(i):
colors=['#984ea3',
'#e41a1c',
'#377eb8',
'#018571',
'#e66101']
return colors[i]
def cbrew2(i):
colors=['#984ea3',
'#e41a1c',
'#377eb8',
'#4daf4a',
'#ff7f00']
return colors[i]
This is main sensitivity analysis calculation. We calculate envelopes over a range of parameters for Fossils and Renewables.
Fossils parameters:
eroei_el_ns
with the number of steps indicated in segments
.#Fossil part
segments=800
eroei_el_ns=[np.linspace(20,87,segments),np.linspace(23,58,segments),np.linspace(23,58,segments)]
eroei_el_cf_N, eroei_ccs_N, eroei_ccs_min_N, eroei_ccs_max_N=[],[],[],[]
for i in range(len(names)):
print('calculating '+names[i].replace('\n',' ')+' area...')
eroei_el_cf_n, eroei_ccs_n, eroei_ccs_min_n, eroei_ccs_max_n=[],[],[],[]
for eroi in range(len(eroei_el_ns[i])):
eroei_el_cf, R_cf, eroei_ccs, eroei_ccs_min, eroei_ccs_max=\
eroei_ccs_cycler([names[i]], [eroei_el_ns[i][eroi]], [t0s[i]], CF, CR_ratio, fop, fcap)
eroei_el_cf_n.append(eroei_el_cf)
eroei_ccs_n.append(eroei_ccs)
eroei_ccs_min_n.append(eroei_ccs_min)
eroei_ccs_max_n.append(eroei_ccs_max)
eroei_el_cf_N.append(eroei_el_cf_n)
eroei_ccs_N.append(eroei_ccs_n)
eroei_ccs_min_N.append(eroei_ccs_min_n)
eroei_ccs_max_N.append(eroei_ccs_max_n)
calculating NGCC Gas area... calculating Pulverized Coal area... calculating IGCC Coal area...
OPTIONAL. Save generated data for easier loading later. No need to run these two cells if you have run the previous cell. GitHub does not hold large files, so you would have to generate your own data by running the sensitivity cell as least once.
#optional SAVE
np.save('eroei_el_cf_N',eroei_el_cf_N)
np.save('eroei_ccs_N',eroei_ccs_N)
np.save('eroei_ccs_min_N',eroei_ccs_min_N)
np.save('eroei_ccs_max_N',eroei_ccs_max_N)
#optional LOAD
eroei_el_cf_N=np.load('eroei_el_cf_N.npz')
eroei_ccs_N=np.load('eroei_ccs_N.npz')
eroei_ccs_min_N=np.load('eroei_ccs_min_N.npz')
eroei_ccs_max_N=np.load('eroei_ccs_max_N.npz')
Renewables parameters:
segments
number of steps (set to $800$ in the Fossils cell. This is stored in eroei_el_ren
.#RE part
eroei_el_ren=[np.linspace(9,30,segments),np.linspace(20,40,segments)]
eroei_el_pv=[[9,17],[12,23],[20,40]] #current medium and high PV EROEI_el ranges, future 2022 range
eroei_el_wind=[20,40] #current medium and high PV EROEI_el ranges, future 2022 range
eroei_el_eu=[[18,35],[0.33,0.67],21.9] #Table S3, current EU solar PV and wind EROEI_el, shares and resulting EROEI_disp
eroei_ren_names=['PV','Wind']
eroei_disp_n=[]
for i in range(len(eroei_ren_names)):
print('calculating '+eroei_ren_names[i]+' area...')
eroei_disp=eroei_disp_cycler(ESOI, phi/100.0, eta/100.0, c, eroei_el_ren[i])
eroei_disp_n.append(eroei_disp)
calculating PV area... calculating Wind area...
#optional SAVE
np.save('eroei_disp_n',eroei_disp_n)
#optional LOAD
eroei_disp_n=np.load('eroei_disp_n.npz')
We convert the sensitivity analysis mutlidimensional arrays into points on the canvas.
#Calculate points for convex hull
all_points={}
#Renewables
for i in range(len(eroei_disp_n)):
print('calculating '+eroei_ren_names[i]+' points...')
if eroei_ren_names[i] not in all_points:all_points[eroei_ren_names[i]]=[[] for ii in range(len(eroei_disp_n[i]))]
for j in range(len(eroei_disp_n[i])):
for eroi in range(len(eroei_disp_n[i][j])):
x0=eroei_el_ren[i][eroi]
y0min=min(eroei_disp_n[i][j][eroi].flatten())
y0max=max(eroei_disp_n[i][j][eroi].flatten())
all_points[eroei_ren_names[i]][j].append([x0,y0min])
all_points[eroei_ren_names[i]][j].append([x0,y0max])
#Fossil part
for i in range(len(names)):
print('calculating '+names[i].replace('\n',' ')+' points...')
if names[i] not in all_points:all_points[names[i]]=[]
for eroi in range(len(eroei_el_ns[i])):
for j in range(len(CF)):
x0=eroei_el_cf_N[i][eroi][0][j]
y0min=min(eroei_ccs_min_N[i][eroi][0][j])
y0max=max(eroei_ccs_max_N[i][eroi][0][j])
all_points[names[i]].append([x0,y0min])
all_points[names[i]].append([x0,y0max])
calculating PV points... calculating Wind points... calculating NGCC Gas points... calculating Pulverized Coal points... calculating IGCC Coal points...
Define a simple function for finding nearest values in an array, in order to be able to identify intersections between Fossil $EROEI_{CCS}$ and Renewable $EROEI_{disp}$ cases.
def find_nearest(array, value):
array = np.asarray(array)
idx = (np.abs(array - value)).argmin()
return idx
We define the $EROEI_{fuel}$ values, for each Fossil fuel source to calculate the intersections for. We reasonably use a Capacity Factor of $CF=85$% and a Capture Ratio of $CR=90$%. For the EROEI we select the highest value and the mid-point of the defined sensitvity range, for each fuel.
Once this is done, we cycle through all Fossil-Renewable combinations and take note of the matching $(EROEI_{el},\varphi)$ combinations, for all cases when $EROEI_{CCS}=EROEI_{disp}$.
eroei_ids=[-1,segments//2] #Fossil EROEI range IDs to calculate intersections for. Default: LAST and MID
eroei_ida=['high','mid'] #corresponding scenario names
CF_id=1 #85%
CR_id=0 #90%
all_intersections=[]
for eroei_idi in range(len(eroei_ids)):
eroei_id=eroei_ids[eroei_idi]
case_id=eroei_ida[eroei_idi]
for name_id in range(len(names)):
name=names[name_id].replace('\n',' ')
fs_x=eroei_el_cf_N[name_id][eroei_id][0][CF_id]
fs_y_max=eroei_ccs_max_N[name_id][eroei_id][0][CF_id][CR_id]
fs_y_min=eroei_ccs_min_N[name_id][eroei_id][0][CF_id][CR_id]
fs_y=np.round((fs_y_max+fs_y_min)/2.0,2)
for i in range(len(eroei_disp_n))[::-1]:
for j in range(len(eroei_disp_n[i])):
minimums=all_points[eroei_ren_names[i]][j][0::2] #bottom line on RE area plots
maximums=all_points[eroei_ren_names[i]][j][1::2] #top line on RE area plots
intersections=[]
for eroi in range(len(eroei_el_ren[i])):
if minimums[eroi][1]<fs_y<maximums[eroi][1]:
intersections.append(eroi)
if len(intersections)>0:
imin=min(intersections)
rmin=np.round(eroei_el_ren[i][imin],1)
emin=eroei_disp_n[i][j][imin]
phimin=np.round(phi[find_nearest(emin,fs_y)],1)
imax=max(intersections)
rmax=np.round(eroei_el_ren[i][imax],1)
emax=eroei_disp_n[i][j][imax]
phimax=np.round(phi[find_nearest(emax,fs_y)],1)
else:
imax=np.nan
emax=np.nan
rmax=np.nan
phimax=np.nan
if phimax>34:phimax=35.0
if phimin<0.31:phimin=0.0
all_intersections.append([case_id,name,eroei_ren_names[i],ESOI_tech[j],
np.round(eroei_el_ns[name_id][eroei_id],1),np.round(fs_x,1),np.round(fs_y,1),
rmin,phimin,rmax,phimax])
all_intersections=pd.DataFrame(all_intersections)
all_intersections.columns=['EROEI case','Fossil','Renewable','Storage',
'Fossil $EROEI_{fuel}$','Fossil $EROEI_{el}$','Fossil $EROEI_{CCS}$',
'RE $EROEI_{el}$ min',r'$\varphi$ min (%)',
'RE $EROEI_{el}$ max',r'$\varphi$ max (%)']
all_intersections['Case ID']=[i+1 for i in range(len(all_intersections))]
all_intersections2=all_intersections.dropna()
all_intersections2['Case ID']=[i+1 for i in range(len(all_intersections2))]
C:\Users\csala\AppData\Local\Continuum\anaconda2\envs\python3\lib\site-packages\ipykernel_launcher.py:52: SettingWithCopyWarning: A value is trying to be set on a copy of a slice from a DataFrame. Try using .loc[row_indexer,col_indexer] = value instead See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy
We present the obtained intersections in Table S4, sorting by EROEI case, Fossil fuel, Renewable type and Energy storage scenario.
all_intersections2.set_index(['EROEI case','Fossil','Renewable','Storage','Case ID'])
Fossil $EROEI_{fuel}$ | Fossil $EROEI_{el}$ | Fossil $EROEI_{CCS}$ | RE $EROEI_{el}$ min | $\varphi$ min (%) | RE $EROEI_{el}$ max | $\varphi$ max (%) | |||||
---|---|---|---|---|---|---|---|---|---|---|---|
EROEI case | Fossil | Renewable | Storage | Case ID | |||||||
high | NGCC Gas | Wind | Battery storage | 1 | 87.0 | 31.0 | 21.3 | 21.3 | 0.0 | 40.0 | 26.2 |
Hydrogen storage | 2 | 87.0 | 31.0 | 21.3 | 21.3 | 0.0 | 32.2 | 35.0 | |||
Battery + Curtailment | 3 | 87.0 | 31.0 | 21.3 | 22.9 | 0.0 | 40.0 | 22.3 | |||
Pumped hydro storage | 4 | 87.0 | 31.0 | 21.3 | 21.3 | 0.0 | 23.5 | 35.0 | |||
PV | Battery storage | 5 | 87.0 | 31.0 | 21.3 | 21.3 | 0.0 | 30.0 | 16.3 | ||
Hydrogen storage | 6 | 87.0 | 31.0 | 21.3 | 21.3 | 0.0 | 30.0 | 29.9 | |||
Battery + Curtailment | 7 | 87.0 | 31.0 | 21.3 | 22.9 | 0.0 | 30.0 | 12.4 | |||
Pumped hydro storage | 8 | 87.0 | 31.0 | 21.3 | 21.3 | 0.0 | 23.5 | 35.0 | |||
Pulverized Coal | PV | Battery storage | 9 | 58.0 | 13.5 | 8.4 | 9.0 | 8.7 | 11.6 | 35.0 | |
Hydrogen storage | 10 | 58.0 | 13.5 | 8.4 | 9.0 | 8.7 | 11.6 | 35.0 | |||
Battery + Curtailment | 11 | 58.0 | 13.5 | 8.4 | 9.0 | 0.4 | 12.9 | 35.0 | |||
Pumped hydro storage | 12 | 58.0 | 13.5 | 8.4 | 9.0 | 30.8 | 9.1 | 35.0 | |||
IGCC Coal | PV | Battery storage | 13 | 58.0 | 11.5 | 8.0 | 9.0 | 14.2 | 11.0 | 35.0 | |
Hydrogen storage | 14 | 58.0 | 11.5 | 8.0 | 9.0 | 13.8 | 11.1 | 35.0 | |||
Battery + Curtailment | 15 | 58.0 | 11.5 | 8.0 | 9.0 | 5.2 | 12.2 | 35.0 | |||
mid | NGCC Gas | Wind | Battery storage | 16 | 53.5 | 21.4 | 15.6 | 20.0 | 16.5 | 29.4 | 35.0 |
Hydrogen storage | 17 | 53.5 | 21.4 | 15.6 | 20.0 | 24.8 | 22.6 | 35.0 | |||
Battery + Curtailment | 18 | 53.5 | 21.4 | 15.6 | 20.0 | 11.3 | 33.8 | 35.0 | |||
PV | Battery storage | 19 | 53.5 | 21.4 | 15.6 | 15.6 | 0.0 | 29.4 | 35.0 | ||
Hydrogen storage | 20 | 53.5 | 21.4 | 15.6 | 15.6 | 0.0 | 22.6 | 35.0 | |||
Battery + Curtailment | 21 | 53.5 | 21.4 | 15.6 | 16.8 | 0.0 | 30.0 | 30.6 | |||
Pumped hydro storage | 22 | 53.5 | 21.4 | 15.6 | 15.6 | 0.0 | 17.0 | 35.0 | |||
Pulverized Coal | PV | Hydrogen storage | 23 | 40.5 | 10.6 | 6.8 | 9.0 | 31.0 | 9.4 | 35.0 | |
Battery + Curtailment | 24 | 40.5 | 10.6 | 6.8 | 9.0 | 25.4 | 9.8 | 35.0 | |||
IGCC Coal | PV | Hydrogen storage | 25 | 40.5 | 9.2 | 6.6 | 9.0 | 34.8 | 9.0 | 35.0 | |
Battery + Curtailment | 26 | 40.5 | 9.2 | 6.6 | 9.0 | 30.4 | 9.4 | 35.0 |
In this section we represent the sensitivity analysis results from both the Fossils and the Renewables in one graph.
For aesthetics, we define a custom shaded area plot, which fades towards the edges.
def alpharater(g,n,t):
if g<t*n:
return (g+1)*1.0/(t*n)
elif g+1>(1-t)*n:
return ((n-g))*1.0/(t*n)
else:
return 1
#Renewables plotter
def plot_renewables(axes,segment_wise=False,more_labels=False,edgeless=False,blow_up=False):
for i in range(len(eroei_disp_n))[::-1]:
print('plotting '+eroei_ren_names[i]+' area...')
for j in range(len(eroei_disp_n[i])):
points=all_points[eroei_ren_names[i]][j]
ax=axes[j//2][j%2]
colr=cbrew(4-i)
if segment_wise:
nfaces=400
if not blow_up:
tfaces=0.2
basealpha=0.15 #0.2
lws=0.95
else:
tfaces=0.2
basealpha=0.3 #0.4
lws=1.7
for g,p in enumerate(np.split(np.array(points),nfaces)):
psegment=convexhull(p)
poly = plt.Polygon(psegment, edgecolor=colr,
alpha=basealpha*alpharater(g,nfaces,tfaces),fill=False,lw=lws)
ax.add_patch(poly)
if not edgeless:
ax.plot(psegment.T[0][:2],psegment.T[1][:2],c=colr,alpha=0.9,lw=0.5)
ax.plot(psegment.T[0][2:],psegment.T[1][2:],c=colr,alpha=0.9,lw=0.5)
if g==0:
ax.plot(psegment.T[0][0::2],psegment.T[1][0::2],c=colr,alpha=0.9,lw=0.5)
if g==len(np.split(np.array(points),nfaces))-1:
ax.plot(psegment.T[0][1::2],psegment.T[1][1::2],c=colr,alpha=0.9,lw=0.5)
else:
poly = plt.Polygon(convexhull(points), edgecolor="none",facecolor=colr,alpha=0.4)
poly2 = plt.Polygon(convexhull(points), edgecolor=colr,alpha=0.8,fill=False,lw=1)
ax.add_patch(poly)
ax.add_patch(poly2)
legendpoly = plt.Polygon([[0,0],[0,-1],[-1,-1]], edgecolor=colr,facecolor=colr,alpha=legendalpha,lw=1,label=eroei_ren_names[i])
ax.add_patch(legendpoly)
if (i==0):
ax.set_title(ESOI_tech[j],fontsize=13)
ax.text(0.04, 0.12, string.ascii_lowercase[(j%2)*2+j//2],
horizontalalignment='left',
verticalalignment='top',
transform=ax.transAxes,size=13,alpha=0.7,
bbox=dict(boxstyle="round", fc="w", alpha=0.7))
ax.plot([0,100],[0,100],c='grey',lw=1,ls='--')
if more_labels:
if (j==2):
ax.text(0.22, 0.27, 'no storage',
horizontalalignment='center',
verticalalignment='center',
transform=ax.transAxes,size=11,rotation=40, color='grey')
ax.arrow(35, 35, 0, -1.5, head_width=0.5, head_length=0.5, fc='k', ec='k',ls='-',lw=0.5)
ax.arrow(35, 32, 0, -15, head_width=0.5, head_length=0.5, fc='k', ec='k',ls='-',lw=0.5)
ax.text(34.5, 34, 'curtailment',
horizontalalignment='right',
verticalalignment='center', rotation=0, color='k')
ax.text(34.5, 20, 'more\nstorage',
horizontalalignment='right',
verticalalignment='center', rotation=0, color='k')
if (j==3):
ax.arrow(31.5, 18.7, 0, 6, head_width=0.5, head_length=0.5, fc='k', ec='k',ls='-',lw=0.5)
ax.arrow(31.5, 18.7, -4, -4, head_width=0.5, head_length=0.5, fc='k', ec='k',ls='-',lw=0.5)
ax.text(32.5, 21.7, 'lower\ncapture\nratio',
horizontalalignment='left',
verticalalignment='center', rotation=0, color='k')
ax.text(30, 12.2, 'lower\ncapacity\nfactor',
horizontalalignment='left',
verticalalignment='center', rotation=0, color='k')
if (j==1):
arrowoffset=[1.5,0.75,1.5]
arrowtextoffset=[1.5,4.5,1]
arrowtext=['PV moderate ',' PV high','PV in 2022']
arrowtextoffsetx=[1,1,1]
arrowtextoffsety=[35.75,32.5,29.25]
#From main text + Breyer estimate
for a in range(len(eroei_el_pv)):
x1,y1=eroei_el_pv[a][0], arrowoffset[a]
w1=eroei_el_pv[a][1]-eroei_el_pv[a][0]
ax.arrow(x1,y1,w1,0,length_includes_head=True,
head_width=0.5, head_length=0.5, fc=colr, ec=colr,ls='-',lw=0.5)
ax.arrow(x1+w1,y1,-w1,0,length_includes_head=True,
head_width=0.5, head_length=0.5, fc=colr, ec=colr,ls='-',lw=0.5)
ax.text(arrowtextoffset[a]+x1+w1/2, 3, arrowtext[a],
horizontalalignment='center',
verticalalignment='center', rotation=0, color=colr)
ax.text(arrowtextoffsetx[a], arrowtextoffsety[a], arrowtext[a].strip()+\
': $'+str(eroei_el_pv[a][0])+'-'+str(eroei_el_pv[a][1])+'$',
horizontalalignment='left',
verticalalignment='top', rotation=0, color=colr)
x2=eroei_el_eu[0][0]*eroei_el_eu[1][0]+eroei_el_eu[0][1]*eroei_el_eu[1][1]
y2=eroei_el_eu[2]
ax.scatter([x2],[y2],color='k',zorder=100)
#From Table S2 in supplementary
ax.text(x2*1.04, y2*1.03, 'European\nsystem\nexample',
horizontalalignment='left',
verticalalignment='center', rotation=0, color='k')
#From Davidsson, S., Höök, M. & Wall, G. A review of life cycle assessments on wind energy systems
elif (j==1):
if more_labels:
xk=3
x1,y1=eroei_el_wind[0], 6
w1=eroei_el_wind[1]-eroei_el_wind[0]
ax.arrow(x1+w1-xk,y1,xk,0,length_includes_head=True,
head_width=0.5, head_length=0.5, fc=colr, ec=colr,ls='--',lw=0.5)
ax.arrow(x1+w1-xk,y1,-w1+xk,0,length_includes_head=True,
head_width=0.5, head_length=0.5, fc=colr, ec=colr,ls='-',lw=0.5)
ax.text(x1+w1/2, y1+1.5, 'Wind',
horizontalalignment='center',
verticalalignment='center', rotation=0, color=colr)
ax.text(1, 39,
'Wind: $EROEI_{el}='+str(eroei_el_wind[0])+'-'+str(eroei_el_wind[1])+'$+',
horizontalalignment='left',
verticalalignment='top', rotation=0, color=colr)
#Fossil plotter
def plot_fossils(axes,segment_wise=False,more_labels=False,edgeless=False,blow_up=False):
fossil_minmax=[[],[]]
for i in range(len(names)):
points=all_points[names[i]]
print('plotting '+names[i].replace('\n',' ')+' area...')
for axi in range(len(axes)):
for axj in range(len(axes[axi])):
ax=axes[axi][axj]
colr=cbrew(i)
if segment_wise:
nfaces=100
if blow_up:
tfaces=0.2
basealpha=0.15 #0.2
lws=1.05
else:
tfaces=0.24
basealpha=0.11 #0.15
lws=1
for g,p in enumerate(np.split(np.array(points),nfaces)):
psegment=convexhull(p)
poly = plt.Polygon(psegment, edgecolor=colr, facecolor=colr,
alpha=basealpha*alpharater(g,nfaces,tfaces),lw=lws)
ax.add_patch(poly)
if not edgeless: ax.plot([psegment.T[0][-1],psegment.T[0][-2]],[psegment.T[1][-1],psegment.T[1][-2]],
c=colr,alpha=0.9,lw=0.5)
if g==len(np.split(np.array(points),nfaces))-1:
if not edgeless: ax.plot([psegment.T[0][0],psegment.T[0][-1]],[psegment.T[1][0],psegment.T[1][-1]],
c=colr,alpha=0.9,lw=0.5)
if axi==0 and axj==0:
fossil_minmax[1].append(np.round(psegment.T[0][0],1))
psegment=convexhull(points)
if axi==0 and axj==0:
fossil_minmax[0].append(np.round(points[0][0],1))
for g in range(len(psegment.T[0])-2):
if not edgeless: ax.plot([psegment.T[0][g],psegment.T[0][g+1]],[psegment.T[1][g],psegment.T[1][g+1]],
c=colr,alpha=0.9,lw=0.5)
else:
poly = plt.Polygon(convexhull(points), edgecolor="k",facecolor=colr,alpha=0.4)
poly2 = plt.Polygon(convexhull(points), edgecolor=colr,alpha=0.8,fill=False,lw=1)
ax.add_patch(poly)
ax.add_patch(poly2)
legendpoly = plt.Polygon([[0,0],[0,-1],[-1,-1]], edgecolor=colr,facecolor=colr,alpha=legendalpha,lw=1,label=names[i].replace('\n',' '))
ax.add_patch(legendpoly)
if more_labels:
ax=axes[1][1]
arrowoffset=[0.75,2.25,1.5]
arrowtextoffsetx=[1,1,1]
arrowtextoffsety=[39,32.5,26]
arrowtextoffsetx2=[23,15,14.5]
arrowtextoffsety2=[2.25,5,2.25]
sep=['\t ','\t ','\t ']
x1,y1=fossil_minmax[0][i], arrowoffset[i]
w1=fossil_minmax[1][i]-fossil_minmax[0][i]
ax.arrow(x1,y1,w1,0,length_includes_head=True,
head_width=0.5, head_length=0.5, fc=colr, ec=colr,ls='-',lw=0.5)
ax.arrow(x1+w1,y1,-w1,0,length_includes_head=True,
head_width=0.5, head_length=0.5, fc=colr, ec=colr,ls='-',lw=0.5)
ax.text(arrowtextoffsetx[i], arrowtextoffsety[i],
'$EROEI_{th}='+str(int(min(eroei_el_ns[i])))+'-'+str(int(max(eroei_el_ns[i])))+'$'+\
'\n$EROEI_{el}='+str(fossil_minmax[0][i])+'-'+str(fossil_minmax[1][i])+'$',
horizontalalignment='left',
verticalalignment='top', rotation=0, color=colr)
ax.text(arrowtextoffsetx2[i], arrowtextoffsety2[i], names[i].replace('\n',''),
horizontalalignment='left', verticalalignment='center', rotation=0, color=colr)
#Intersections plotter
def plot_intersections(axes,blow_up=False):
print('plotting Intersections...')
for i in all_intersections2.index:
j=ESOI_tech.index(all_intersections2.loc[i]['Storage'])
z=names.index(all_intersections2.loc[i]['Fossil'].replace(' ','\n'))
w=eroei_ren_names.index(all_intersections2.loc[i]['Renewable'])
ax=axes[j//2][j%2]
x=all_intersections2.loc[i]['Fossil $EROEI_{el}$']
y=all_intersections2.loc[i]['Fossil $EROEI_{CCS}$']
if (i-1)%4==0:
for axi in axes:
for axj in axi:
axj.scatter([x],[y],c='k',s=20,zorder=100)
axj.scatter([x],[y],c=cbrew(z),s=12,zorder=100)
x1=float(all_intersections2.loc[i]['RE $EROEI_{el}$ min'])
x2=float(all_intersections2.loc[i]['RE $EROEI_{el}$ max'])
colr=cbrew(4-w)
if (not blow_up and x1>13) or (blow_up and x2<maxlabel):
ax.plot([x1,x2],[y,y],c='k',lw=2,alpha=0.3)
ax.plot([x1,x2],[y,y],c=colr,lw=1,alpha=0.5)
if blow_up: z=0.1
else: z=1
ax.text(x1+(x2-x1)/2, y+z+z*(w), all_intersections2.loc[i]['Case ID'],
horizontalalignment='center',
verticalalignment='bottom', fontsize=8,
bbox=dict(boxstyle="round", fc=colr, alpha=0.7))
We plot the Renewables and Fossils on a common axis, with $EROEI_{el}$ on the $x$ and $EROEI_{CCS}$ or $EROEI_{disp}$ on the $y$ axis. We split the figure into 4 sub-plots, split by enegeyr storage cases, exactly as in Figure 3. Each sub-plot will have 5 shaded trapezoid areas, for Wind, PV, NGCC Gas, Pulverized Coal and IGCC Coal, representing the extent ranges of the sensitivity analysis.
We produce a set of these area plots, the first of which is the explicit representation of the Intersection cases from Table S4.
#Plots init
fig,axes=plt.subplots(2,2,figsize=(9,8),subplot_kw=dict(facecolor='#EEEEEE',axisbelow=True))
for axi in axes:
for ax in axi:
ax.grid(color='white', linestyle='solid')
plt.subplots_adjust(wspace=0.1)
plt.subplots_adjust(hspace=0.15)
legendalpha=0.7
#Renewables part
plot_renewables(axes,segment_wise=True)
#Fossils part
plot_fossils(axes,segment_wise=True)
#Intersections
plot_intersections(axes)
#Format labels
axes[0][0].set_ylabel('$EROEI_{CCS}$ | $EROEI_{disp}$',fontsize=14,y=0)
axes[1][1].set_xlabel('$EROEI_{el}$',fontsize=14,x=0)
maxlabel=max(np.array(eroei_el_ren).flatten())
labelrange=[i if i%2==0 else '' for i in range(int(maxlabel)+1) ]
for axi in axes:
for ax in axi:
ax.set_xlim(0,maxlabel)
ax.set_ylim(0,maxlabel)
ax.set_xticks(range(int(maxlabel)+1))
ax.set_yticks(range(int(maxlabel)+1))
ax.set_yticks(range(int(maxlabel)+1))
ax.set_xticklabels(labelrange,fontsize=8)
ax.set_yticklabels(labelrange,fontsize=8)
axes[0][0].legend(loc=(0.03,0.52),framealpha=0,fontsize=11)
#plt.suptitle('Comparison of adjusted EROEIs for \nfossil fuels with carbon-capture and storage $EROEI_{CCS}$ and \ndispatchable renewables with energy storage $EROEI_{disp}$',y=1.03,fontsize=13)
plt.savefig('plot/extra/figS4.png',dpi=150,bbox_inches = 'tight', pad_inches = 0.1)
plt.show()
plotting Wind area... plotting PV area... plotting NGCC Gas area... plotting Pulverized Coal area... plotting IGCC Coal area... plotting Intersections...
Next, we produce a blow-up of the same plot of the $EROEI_{el}$ range of $3-14$ for better visuals.
#Plots init
fig,axes=plt.subplots(2,2,figsize=(9,8),subplot_kw=dict(facecolor='#EEEEEE',axisbelow=True))
for axi in axes:
for ax in axi:
ax.grid(color='white', linestyle='solid')
plt.subplots_adjust(wspace=0.1)
plt.subplots_adjust(hspace=0.15)
legendalpha=0.7
#Renewables part
plot_renewables(axes,segment_wise=True,blow_up=True)
#Fossil part
plot_fossils(axes,segment_wise=True,blow_up=True)
minlabel=3
maxlabel=14
#Intersections
plot_intersections(axes,blow_up=True)
#Format labels
axes[0][0].set_ylabel('$EROEI_{CCS}$ | $EROEI_{disp}$',fontsize=14,y=0)
axes[1][1].set_xlabel('$EROEI_{el}$',fontsize=14,x=0)
labelrange=[i for i in range(minlabel,int(maxlabel)+1) ]
for axi in axes:
for ax in axi:
ax.set_xlim(minlabel,maxlabel)
ax.set_ylim(minlabel,maxlabel)
ax.set_xticks(range(minlabel,int(maxlabel)+1))
ax.set_yticks(range(minlabel,int(maxlabel)+1))
ax.set_yticks(range(minlabel,int(maxlabel)+1))
ax.set_xticklabels(labelrange,fontsize=8)
ax.set_yticklabels(labelrange,fontsize=8)
axes[0][0].legend(loc=(0.03,0.52),framealpha=0,fontsize=11)
#plt.suptitle('Comparison of adjusted EROEIs for \nfossil fuels with carbon-capture and storage $EROEI_{CCS}$ and \ndispatchable renewables with energy storage $EROEI_{disp}$',y=1.03,fontsize=13)
plt.savefig('plot/extra/figS5.png',dpi=150,bbox_inches = 'tight', pad_inches = 0.1)
plt.show()
plotting Wind area... plotting PV area... plotting NGCC Gas area... plotting Pulverized Coal area... plotting IGCC Coal area... plotting Intersections...
And finally, we produce an annotated plots, with the $EROEI_{fuel}$s and $EROEI_{el}$ indicated as on-figure texts, as well as teh European system example presented in Table S3. This becomes Figure 4. in the paper.
#Plots init
fig,axes=plt.subplots(2,2,figsize=(9,8),subplot_kw=dict(facecolor='#FFFFFF',axisbelow=True))
for axi in axes:
for ax in axi:
ax.grid(color='white', linestyle='solid')
plt.subplots_adjust(wspace=0.1)
plt.subplots_adjust(hspace=0.15)
legendalpha=0.7
#Renewables part
plot_renewables(axes,segment_wise=True,more_labels=True)
#Fossil part
plot_fossils(axes,segment_wise=True,more_labels=True)
#Intersections
#plot_intersections(axes)
#Format labels
for axi in axes:
for ax in axi:
ax.set_xlim(0,max(np.array(eroei_el_ren).flatten()))
ax.set_ylim(0,max(np.array(eroei_el_ren).flatten()))
axes[0][0].set_ylabel('$EROEI_{CCS}$ | $EROEI_{disp}$',fontsize=14,y=0)
axes[1][1].set_xlabel('$EROEI_{el}$',fontsize=14,x=-0.05)
axes[0][1].set_yticklabels('')
axes[1][1].set_yticklabels('')
axes[0][0].set_xticklabels('')
axes[0][1].set_xticklabels('')
axes[0][0].legend(loc=(0.03,0.52),framealpha=0,fontsize=11)
#plt.suptitle('Comparison of adjusted EROEIs for \nfossil fuels with carbon-capture and storage $EROEI_{CCS}$ and \ndispatchable renewables with energy storage $EROEI_{disp}$',y=1.03,fontsize=13)
plt.savefig('plot/png/fig4.png',dpi=150,bbox_inches = 'tight', pad_inches = 0.1)
plt.savefig('plot/eps/fig4.eps',format='eps',bbox_inches = 'tight', pad_inches = 0.1)
plt.savefig('plot/pdf/fig4.pdf',format='pdf',bbox_inches = 'tight', pad_inches = 0.1)
plt.show()
plotting Wind area... plotting PV area... plotting NGCC Gas area... plotting Pulverized Coal area... plotting IGCC Coal area...
An alternative iteration of Figure 4. with the area edges removed.
#Plots init
fig,axes=plt.subplots(2,2,figsize=(9,8),subplot_kw=dict(facecolor='#FFFFFF',axisbelow=True))
for axi in axes:
for ax in axi:
ax.grid(color='white', linestyle='solid')
plt.subplots_adjust(wspace=0.1)
plt.subplots_adjust(hspace=0.15)
legendalpha=0.7
#Renewables part
plot_renewables(axes,segment_wise=True,more_labels=True,edgeless=True)
#Fossil part
plot_fossils(axes,segment_wise=True,more_labels=True,edgeless=True)
#Intersections
#plot_intersections(axes)
#Format labels
for axi in axes:
for ax in axi:
ax.set_xlim(0,max(np.array(eroei_el_ren).flatten()))
ax.set_ylim(0,max(np.array(eroei_el_ren).flatten()))
axes[0][0].set_ylabel('$EROEI_{CCS}$ | $EROEI_{disp}$',fontsize=14,y=0)
axes[1][1].set_xlabel('$EROEI_{el}$',fontsize=14,x=-0.05)
axes[0][1].set_yticklabels('')
axes[1][1].set_yticklabels('')
axes[0][0].set_xticklabels('')
axes[0][1].set_xticklabels('')
axes[0][0].legend(loc=(0.03,0.52),framealpha=0,fontsize=11)
#plt.suptitle('Comparison of adjusted EROEIs for \nfossil fuels with carbon-capture and storage $EROEI_{CCS}$ and \ndispatchable renewables with energy storage $EROEI_{disp}$',y=1.03,fontsize=13)
plt.savefig('plot/extra/fig4_edgeless.png',dpi=150,bbox_inches = 'tight', pad_inches = 0.1)
plt.show()
plotting Wind area... plotting PV area... plotting NGCC Gas area... plotting Pulverized Coal area... plotting IGCC Coal area...
Close Excel workbook.
wb.close()
Finally we create Figure S2. First, we read in the capacity factors for the years 2005 and 2015 of all Coal and Gas power plants in the US. We then produce 4 vertical histograms, for each fuel source and each year. We then calculate the mean of these distributions to yield the mean capacity factor per fuel, per year. Please note that due to the structure of the data, this is not a weighted value by plant capacity or yearly generated energy - all plants of all sizes have the same weight in the mean.
df=pd.read_excel('cfdata.xlsx')
fig,ax=plt.subplots(1,1,figsize=(6,5),subplot_kw=dict(facecolor='#FFFFFF',axisbelow=True))
ax.grid(color='white', linestyle='solid')
ax.barh(df.index,-df['PC 2005']/sum(df['PC 2005'])*100,0.4,color=cbrew(3),label='Coal 2005',alpha=legendalpha)
ax.barh(df.index,df['PC 2015']/sum(df['PC 2015'])*100,0.4,color=cbrew(2),label='Coal 2015',alpha=legendalpha)
ax.barh(df.index-0.4,-df['NG 2005']/sum(df['NG 2005'])*100,0.4,color=cbrew(1),label='Gas 2005',alpha=legendalpha)
ax.barh(df.index-0.4,df['NG 2015']/sum(df['NG 2015'])*100,0.4,color=cbrew(4),label='Gas 2015',alpha=legendalpha)
z=np.average(df.index+0.5,weights=df['NG 2005'])-1
ax.plot([-30,-25],[z,z],ls='--',color=cbrew(1))
ax.text(-29.7,z+0.2,str(np.round((z+1)*10,1))+ '%',color=cbrew(1))
z=np.average(df.index+0.5,weights=df['NG 2015'])-1
ax.plot([15,20],[z,z],ls='--',color=cbrew(4))
ax.text(15,z+0.15,str(np.round((z+1)*10,1))+ '%',color=cbrew(4))
z=np.average(df.index+0.5,weights=df['PC 2005'])-1
ax.plot([-30,-25],[z,z],ls='--',color=cbrew(3))
ax.text(-29.7,z+0.2,str(np.round((z+1)*10,1))+ '%',color=cbrew(3))
ax.text(-29.7,z-0.4,'mean',color=cbrew(3))
z=np.average(df.index+0.5,weights=df['PC 2015'])-1
#z=np.average(df.index+0.5)-1
ax.plot([15,20],[z,z],ls='--',color=cbrew(2))
ax.text(15,z-0.5,str(np.round((z+1)*10,1))+ '%',color=cbrew(2))
fig.canvas.draw()
labels = ['']+[item.get_text().replace(u'\u2212','') for item in ax.get_xticklabels()[2:-1]]+['']
ax.set_xticklabels(labels)
ax.set_yticks(df.index)
ax.set_xlim([-30,20])
ax.set_yticklabels(df['CF'][::-1])
ax.legend(bbox_to_anchor=(0.26, 0.22),framealpha=0.5,fontsize=9)
ax.set_ylabel('Plant capacity factor')
ax.set_xlabel('Percentage of plants (United States)\n source: US EIA')
#ax.set_title('Distribution of annual capacity factors of\ngas combined-cycle and coal steam plants')
plt.savefig('plot/png/figS2.png',dpi=150,bbox_inches = 'tight', pad_inches = 0.1)
plt.savefig('plot/eps/figS2.eps',format='eps',bbox_inches = 'tight', pad_inches = 0.1)
plt.savefig('plot/pdf/figS2.pdf',format='pdf',bbox_inches = 'tight', pad_inches = 0.1)
plt.show()
Denes Csala | MIT License | 2018