This tutorial will walk through the process of going from TMY data to AC power using the SAPM.
Table of contents:
This tutorial has been tested against the following package versions:
It should work with other Python and Pandas versions. It requires pvlib >= 0.3.0 and IPython >= 3.0.
Authors:
These are just your standard interactive scientific python imports that you'll get very used to using.
# built-in python modules
import os
import inspect
# scientific python add-ons
import numpy as np
import pandas as pd
# plotting stuff
# first line makes the plots appear in the notebook
%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib as mpl
# seaborn makes your plots look better
try:
import seaborn as sns
sns.set(rc={"figure.figsize": (12, 6)})
sns.set_color_codes()
except ImportError:
print('We suggest you install seaborn using conda or pip and rerun this cell')
# finally, we import the pvlib library
import pvlib
pvlib comes with a couple of TMY files, and we'll use one of them for simplicity. You could also load a file from disk, or specify a url. See this NREL website for a list of TMY files:
http://rredc.nrel.gov/solar/old_data/nsrdb/1991-2005/tmy3/by_state_and_city.html
# Find the absolute file path to your pvlib installation
pvlib_abspath = os.path.dirname(os.path.abspath(inspect.getfile(pvlib)))
# absolute path to a data file
datapath = os.path.join(pvlib_abspath, 'data', '703165TY.csv')
# read tmy data with year values coerced to a single year
tmy_data, meta = pvlib.tmy.readtmy3(datapath, coerce_year=2015)
tmy_data.index.name = 'Time'
# TMY data seems to be given as hourly data with time stamp at the end
# shift the index 30 Minutes back for calculation of sun positions
tmy_data = tmy_data.shift(freq='-30Min')
The file handling above looks complicated because we're trying to account for the many different ways that people will run this notebook on their systems. You can just put a simple string path into the readtmy3
function if you know where the file is.
Let's look at the imported version of the TMY file.
tmy_data.head()
ETR | ETRN | GHI | GHISource | GHIUncertainty | DNI | DNISource | DNIUncertainty | DHI | DHISource | ... | AOD | AODSource | AODUncertainty | Alb | AlbSource | AlbUncertainty | Lprecipdepth | Lprecipquantity | LprecipSource | LprecipUncertainty | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
Time | |||||||||||||||||||||
2015-01-01 00:30:00-09:00 | 0 | 0 | 0 | 1 | 0 | 0 | 1 | 0 | 0 | 1 | ... | 0.051 | F | 8 | 0.24 | F | 8 | -9900 | -9900 | ? | 0 |
2015-01-01 01:30:00-09:00 | 0 | 0 | 0 | 1 | 0 | 0 | 1 | 0 | 0 | 1 | ... | 0.051 | F | 8 | 0.24 | F | 8 | -9900 | -9900 | ? | 0 |
2015-01-01 02:30:00-09:00 | 0 | 0 | 0 | 1 | 0 | 0 | 1 | 0 | 0 | 1 | ... | 0.051 | F | 8 | 0.24 | F | 8 | -9900 | -9900 | ? | 0 |
2015-01-01 03:30:00-09:00 | 0 | 0 | 0 | 1 | 0 | 0 | 1 | 0 | 0 | 1 | ... | 0.051 | F | 8 | 0.24 | F | 8 | -9900 | -9900 | ? | 0 |
2015-01-01 04:30:00-09:00 | 0 | 0 | 0 | 1 | 0 | 0 | 1 | 0 | 0 | 1 | ... | 0.051 | F | 8 | 0.24 | F | 8 | -9900 | -9900 | ? | 0 |
5 rows × 66 columns
This is a pandas DataFrame
object. It has a lot of great properties that are beyond the scope of our tutorials.
Plot the GHI data from the TMY file
tmy_data['GHI'].plot()
plt.ylabel('Irradiance (W/m**2)')
<matplotlib.text.Text at 0xc4d0b53390>
Before we can calculate power for all times in the TMY file, we will need to calculate:
First, define some PV system parameters.
surface_tilt = 30
surface_azimuth = 180 # pvlib uses 0=North, 90=East, 180=South, 270=West convention
albedo = 0.2
# create pvlib Location object based on meta data
sand_point = pvlib.location.Location(meta['latitude'], meta['longitude'], tz='US/Alaska',
altitude=meta['altitude'], name=meta['Name'].replace('"',''))
print(sand_point)
SAND POINT: latitude=55.317, longitude=-160.517, tz=US/Alaska, altitude=7.0
Calculate the solar position for all times in the TMY file.
The default solar position algorithm is based on Reda and Andreas (2004). Our implementation is pretty fast, but you can make it even faster if you install numba
and use add method='nrel_numba'
to the function call below.
solpos = pvlib.solarposition.get_solarposition(tmy_data.index, sand_point.latitude, sand_point.longitude)
solpos.plot()
<matplotlib.axes._subplots.AxesSubplot at 0xc4cae78ba8>
The funny looking jump in the azimuth is just due to the coarse time sampling in the TMY file.
Calculate extra terrestrial radiation. This is needed for many plane of array diffuse irradiance models.
# the extraradiation function returns a simple numpy array
# instead of a nice pandas series. We will change this
# in a future version
dni_extra = pvlib.irradiance.extraradiation(tmy_data.index)
dni_extra = pd.Series(dni_extra, index=tmy_data.index)
dni_extra.plot()
plt.ylabel('Extra terrestrial radiation (W/m**2)')
<matplotlib.text.Text at 0xc4d10d5470>
Calculate airmass. Lots of model options here, see the atmosphere
module tutorial for more details.
airmass = pvlib.atmosphere.relativeairmass(solpos['apparent_zenith'])
airmass.plot()
plt.ylabel('Airmass')
<matplotlib.text.Text at 0xc4d20d9400>
The funny appearance is due to aliasing and setting invalid numbers equal to NaN
. Replot just a day or two and you'll see that the numbers are right.
Use the Hay Davies model to calculate the plane of array diffuse sky radiation. See the irradiance
module tutorial for comparisons of different models.
poa_sky_diffuse = pvlib.irradiance.haydavies(surface_tilt, surface_azimuth,
tmy_data['DHI'], tmy_data['DNI'], dni_extra,
solpos['apparent_zenith'], solpos['azimuth'])
poa_sky_diffuse.plot()
plt.ylabel('Irradiance (W/m**2)')
<matplotlib.text.Text at 0xc4d1666c18>
Calculate ground diffuse. We specified the albedo above. You could have also provided a string to the surface_type
keyword argument.
poa_ground_diffuse = pvlib.irradiance.grounddiffuse(surface_tilt, tmy_data['GHI'], albedo=albedo)
poa_ground_diffuse.plot()
plt.ylabel('Irradiance (W/m**2)')
<matplotlib.text.Text at 0xc4d363deb8>
Calculate AOI
aoi = pvlib.irradiance.aoi(surface_tilt, surface_azimuth, solpos['apparent_zenith'], solpos['azimuth'])
aoi.plot()
plt.ylabel('Angle of incidence (deg)')
<matplotlib.text.Text at 0xc4d3a3f5f8>
Note that AOI has values greater than 90 deg. This is ok.
Calculate POA irradiance
poa_irrad = pvlib.irradiance.globalinplane(aoi, tmy_data['DNI'], poa_sky_diffuse, poa_ground_diffuse)
poa_irrad.plot()
plt.ylabel('Irradiance (W/m**2)')
plt.title('POA Irradiance')
<matplotlib.text.Text at 0xc4d397f908>
Calculate pv cell and module temperature
pvtemps = pvlib.pvsystem.sapm_celltemp(poa_irrad['poa_global'], tmy_data['Wspd'], tmy_data['DryBulb'])
pvtemps.plot()
plt.ylabel('Temperature (C)')
<matplotlib.text.Text at 0xc4d3db4828>
Get module data from the web.
sandia_modules = pvlib.pvsystem.retrieve_sam(name='SandiaMod')
Choose a particular module
sandia_module = sandia_modules.Canadian_Solar_CS5P_220M___2009_
sandia_module
Vintage 2009 Area 1.701 Material c-Si Cells_in_Series 96 Parallel_Strings 1 Isco 5.09115 Voco 59.2608 Impo 4.54629 Vmpo 48.3156 Aisc 0.000397 Aimp 0.000181 C0 1.01284 C1 -0.0128398 Bvoco -0.21696 Mbvoc 0 Bvmpo -0.235488 Mbvmp 0 N 1.4032 C2 0.279317 C3 -7.24463 A0 0.928385 A1 0.068093 A2 -0.0157738 A3 0.0016606 A4 -6.93e-05 B0 1 B1 -0.002438 B2 0.0003103 B3 -1.246e-05 B4 2.11e-07 B5 -1.36e-09 DTC 3 FD 1 A -3.40641 B -0.0842075 C4 0.996446 C5 0.003554 IXO 4.97599 IXXO 3.18803 C6 1.15535 C7 -0.155353 Notes Source: Sandia National Laboratories Updated 9... Name: Canadian_Solar_CS5P_220M___2009_, dtype: object
Calculate the effective irradiance
effective_irradiance = pvlib.pvsystem.sapm_effective_irradiance(poa_irrad.poa_direct, poa_irrad.poa_diffuse, airmass, aoi, sandia_module)
Run the SAPM using the parameters we calculated above.
sapm_out = pvlib.pvsystem.sapm(effective_irradiance, pvtemps.temp_cell, sandia_module)
print(sapm_out.head())
sapm_out[['p_mp']].plot()
plt.ylabel('DC Power (W)')
i_sc i_mp v_oc v_mp p_mp i_x i_xx Time 2015-01-01 00:30:00-09:00 0.0 0.0 0.0 0.0 0.0 0.0 0.0 2015-01-01 01:30:00-09:00 0.0 0.0 0.0 0.0 0.0 0.0 0.0 2015-01-01 02:30:00-09:00 0.0 0.0 0.0 0.0 0.0 0.0 0.0 2015-01-01 03:30:00-09:00 0.0 0.0 0.0 0.0 0.0 0.0 0.0 2015-01-01 04:30:00-09:00 0.0 0.0 0.0 0.0 0.0 0.0 0.0
<matplotlib.text.Text at 0xc4a4e19cc0>
cec_modules = pvlib.pvsystem.retrieve_sam(name='CECMod')
cec_module = cec_modules.Canadian_Solar_CS5P_220M
photocurrent, saturation_current, resistance_series, resistance_shunt, nNsVth = (
pvlib.pvsystem.calcparams_desoto(poa_irrad.poa_global,
temp_cell=pvtemps['temp_cell'],
alpha_isc=cec_module['alpha_sc'],
module_parameters=cec_module,
EgRef=1.121,
dEgdT=-0.0002677) )
single_diode_out = pvlib.pvsystem.singlediode(photocurrent, saturation_current,
resistance_series, resistance_shunt, nNsVth)
single_diode_out[['p_mp']].plot()
plt.ylabel('DC Power (W)')
<matplotlib.text.Text at 0xc4a50b1da0>
Get the inverter database from the web
sapm_inverters = pvlib.pvsystem.retrieve_sam('sandiainverter')
Choose a particular inverter
sapm_inverter = sapm_inverters['ABB__MICRO_0_25_I_OUTD_US_208_208V__CEC_2014_']
sapm_inverter
Vac 208.000000 Paco 250.000000 Pdco 259.522050 Vdco 40.242603 Pso 1.771614 C0 -0.000025 C1 -0.000090 C2 0.000669 C3 -0.018900 Pnt 0.020000 Vdcmax 65.000000 Idcmax 10.000000 Mppt_low 20.000000 Mppt_high 50.000000 Name: ABB__MICRO_0_25_I_OUTD_US_208_208V__CEC_2014_, dtype: float64
p_acs = pd.DataFrame()
p_acs['sapm'] = pvlib.pvsystem.snlinverter(sapm_out.v_mp, sapm_out.p_mp, sapm_inverter)
p_acs['sd'] = pvlib.pvsystem.snlinverter(single_diode_out.v_mp, single_diode_out.p_mp, sapm_inverter)
p_acs.plot()
plt.ylabel('AC Power (W)')
<matplotlib.text.Text at 0xc4a5d3cd68>
diff = p_acs['sapm'] - p_acs['sd']
diff.plot()
plt.ylabel('SAPM - SD Power (W)')
<matplotlib.text.Text at 0xc4a5768160>
Plot just a few days.
p_acs['2015-07-05':'2015-07-06'].plot()
<matplotlib.axes._subplots.AxesSubplot at 0xc4a576d588>
Some statistics on the AC power
p_acs.describe()
C:\Anaconda3\lib\site-packages\numpy\lib\function_base.py:3823: RuntimeWarning: Invalid value encountered in percentile RuntimeWarning)
sapm | sd | |
---|---|---|
count | 8760.000000 | 4620.000000 |
mean | 23.108233 | 46.013726 |
std | 42.330294 | 49.976381 |
min | -0.020000 | -0.020000 |
25% | -0.020000 | NaN |
50% | -0.020000 | NaN |
75% | 27.946397 | NaN |
max | 211.290155 | 213.051148 |
p_acs.sum()
sapm 202428.122770 sd 212583.414811 dtype: float64
# create data for a y=x line
p_ac_max = p_acs.max().max()
yxline = np.arange(0, p_ac_max)
fig = plt.figure(figsize=(12,12))
ax = fig.add_subplot(111, aspect='equal')
sc = ax.scatter(p_acs['sd'], p_acs['sapm'], c=poa_irrad.poa_global, alpha=1, cmap=mpl.cm.YlGnBu_r)
ax.plot(yxline, yxline, 'r', linewidth=3)
ax.set_xlim(0, None)
ax.set_ylim(0, None)
ax.set_xlabel('Single Diode model')
ax.set_ylabel('Sandia model')
fig.colorbar(sc, label='POA Global (W/m**2)')
<matplotlib.colorbar.Colorbar at 0xc4a563f668>
We can change the value of color value c
to see the sensitivity of model accuracy to measured meterological conditions. It can be useful to define a simple plotting function for this kind of exploratory analysis.
def sapm_sd_scatter(c_data, label=None, **kwargs):
"""Display a scatter plot of SAPM p_ac vs. single diode p_ac.
You need to re-execute this cell if you re-run the p_ac calculation.
Parameters
----------
c_data : array-like
Determines the color of each point on the scatter plot.
Must be same length as p_acs.
kwargs passed to ``scatter``.
Returns
-------
tuple of fig, ax objects
"""
fig = plt.figure(figsize=(12,12))
ax = fig.add_subplot(111, aspect='equal')
sc = ax.scatter(p_acs['sd'], p_acs['sapm'], c=c_data, alpha=1, cmap=mpl.cm.YlGnBu_r, **kwargs)
ax.plot(yxline, yxline, 'r', linewidth=3)
ax.set_xlim(0, None)
ax.set_ylim(0, None)
ax.set_xlabel('Single diode model power (W)')
ax.set_ylabel('Sandia model power (W)')
fig.colorbar(sc, label='{}'.format(label), shrink=0.75)
return fig, ax
sapm_sd_scatter(tmy_data.DryBulb, label='Temperature (deg C)')
(<matplotlib.figure.Figure at 0xc4a5db34e0>, <matplotlib.axes._subplots.AxesSubplot at 0xc4a5dba470>)
sapm_sd_scatter(tmy_data.DNI, label='DNI (W/m**2)')
(<matplotlib.figure.Figure at 0xc4a5e50ba8>, <matplotlib.axes._subplots.AxesSubplot at 0xc4a5e5eeb8>)
sapm_sd_scatter(tmy_data.AOD, label='AOD')
(<matplotlib.figure.Figure at 0xc4a5ec34e0>, <matplotlib.axes._subplots.AxesSubplot at 0xc4a5dbacf8>)
sapm_sd_scatter(tmy_data.Wspd, label='Wind speed', vmax=10)
(<matplotlib.figure.Figure at 0xc4a8311080>, <matplotlib.axes._subplots.AxesSubplot at 0xc4a830eda0>)
Notice the use of the vmax
keyword argument in the above example. The **kwargs
pattern allows us to easily pass non-specified arguments to nested functions.
def sapm_other_scatter(c_data, x_data, clabel=None, xlabel=None, aspect_equal=False, **kwargs):
"""Display a scatter plot of SAPM p_ac vs. something else.
You need to re-execute this cell if you re-run the p_ac calculation.
Parameters
----------
c_data : array-like
Determines the color of each point on the scatter plot.
Must be same length as p_acs.
x_data : array-like
kwargs passed to ``scatter``.
Returns
-------
tuple of fig, ax objects
"""
fig = plt.figure(figsize=(12,12))
if aspect_equal:
ax = fig.add_subplot(111, aspect='equal')
else:
ax = fig.add_subplot(111)
sc = ax.scatter(x_data, p_acs['sapm'], c=c_data, alpha=1, cmap=mpl.cm.YlGnBu_r, **kwargs)
ax.set_xlim(0, None)
ax.set_ylim(0, None)
ax.set_xlabel('{}'.format(xlabel))
ax.set_ylabel('Sandia model power (W)')
fig.colorbar(sc, label='{}'.format(clabel), shrink=0.75)
return fig, ax
sapm_other_scatter(tmy_data.DryBulb, tmy_data.GHI, clabel='Temperature (deg C)', xlabel='GHI (W/m**2)')
(<matplotlib.figure.Figure at 0x11dce1a90>, <matplotlib.axes._subplots.AxesSubplot at 0x11dceb278>)
Next, we will assume that the SAPM model is representative of the real world performance so that we can use scipy's optimization routine to derive simulated PVUSA coefficients. You will need to install scipy to run these functions.
Here's one PVUSA reference:
def pvusa(pvusa_data, a, b, c, d):
"""
Calculates system power according to the PVUSA equation
P = I * (a + b*I + c*W + d*T)
where
P is the output power,
I is the plane of array irradiance,
W is the wind speed, and
T is the temperature
Parameters
----------
pvusa_data : pd.DataFrame
Must contain the columns 'I', 'W', and 'T'
a : float
I coefficient
b : float
I*I coefficient
c : float
I*W coefficient
d : float
I*T coefficient
Returns
-------
power : pd.Series
Power calculated using the PVUSA model.
"""
return pvusa_data['I'] * (a + b*pvusa_data['I'] + c*pvusa_data['W'] + d*pvusa_data['T'])
from scipy import optimize
pvusa_data = pd.DataFrame()
pvusa_data['I'] = poa_irrad.poa_global
pvusa_data['W'] = tmy_data.Wspd
pvusa_data['T'] = tmy_data.DryBulb
popt, pcov = optimize.curve_fit(pvusa, pvusa_data.dropna(), p_acs.sapm.values, p0=(.0001,0.0001,.001,.001))
print('optimized coefs:\n{}'.format(popt))
print('covariances:\n{}'.format(pcov))
optimized coefs: [ 2.09753481e-01 5.39374811e-06 9.34193018e-04 -1.37476822e-03] covariances: [[ 1.31102469e-07 -9.53286064e-11 -7.65512720e-09 -2.26637729e-09] [ -9.53286064e-11 1.83594105e-13 -7.52538706e-13 -2.45362115e-12] [ -7.65512720e-09 -7.52538706e-13 1.22446507e-09 1.87304284e-10] [ -2.26637729e-09 -2.45362115e-12 1.87304284e-10 3.68639894e-10]]
power_pvusa = pvusa(pvusa_data, *popt)
fig, ax = sapm_other_scatter(tmy_data.DryBulb, power_pvusa, clabel='Temperature (deg C)',
aspect_equal=True, xlabel='PVUSA (W)')
maxmax = max(ax.get_xlim()[1], ax.get_ylim()[1])
ax.set_ylim(None, maxmax)
ax.set_xlim(None, maxmax)
ax.plot(np.arange(maxmax), np.arange(maxmax), 'r')
[<matplotlib.lines.Line2D at 0xc4a84aea90>]