diff --git a/docs/examples/irradiance-transposition/plot_rtranpose_limitations.py b/docs/examples/irradiance-transposition/plot_rtranpose_limitations.py new file mode 100644 index 0000000000..8df8339ad4 --- /dev/null +++ b/docs/examples/irradiance-transposition/plot_rtranpose_limitations.py @@ -0,0 +1,181 @@ +""" +Reverse transposition limitations +==================================== + +Unfortunately, sometimes there is not a unique solution. + +Author: Anton Driesse + +""" + +# %% +# +# Introduction +# ------------ +# When irradiance is measured on a tilted plane, it is useful to be able to +# estimate the GHI that produces the POA irradiance. +# The estimation requires inverting a GHI-to-POA irradiance model, +# which involves two parts: +# a decomposition of GHI into direct and diffuse components, +# and a transposition model that calculates the direct and diffuse irradiance +# on the tilted plane. +# Recovering GHI from POA irradiance is termed "reverse transposition." +# +# Unfortunately, for a given POA irradiance value, sometimes there is not a +# unique solution for GHI. +# Different GHI values can produce different combinations of direct and +# diffuse irradiance that sum to the same POA irradiance value. +# +# In this example we look at a single point in time and consider a full range +# of possible GHI and POA global values as shown in figures 3 and 4 of [1]_. +# Then we use :py:meth:`pvlib.irradiance.ghi_from_poa_driesse_2023` to estimate +# the original GHI from POA global. +# +# References +# ---------- +# .. [1] Driesse, A., Jensen, A., Perez, R., 2024. A Continuous form of the +# Perez diffuse sky model for forward and reverse transposition. +# Solar Energy vol. 267. :doi:`10.1016/j.solener.2023.112093` +# + +import numpy as np + +import matplotlib +import matplotlib.pyplot as plt + +from pvlib.irradiance import (erbs_driesse, + get_total_irradiance, + ghi_from_poa_driesse_2023, + ) + +matplotlib.rcParams['axes.grid'] = True + +# %% +# +# Define the conditions that were used for figure 3 in [1]_. +# + +dni_extra = 1366.1 +albedo = 0.25 +surface_tilt = 40 +surface_azimuth = 180 + +solar_azimuth = 82 +solar_zenith = 75 + +# %% +# +# Define a range of possible GHI values and calculate the corresponding +# POA global. First estimate DNI and DHI using the Erbs-Driesse model, then +# transpose using the Perez-Driesse model. +# + +ghi = np.linspace(0, 500, 100+1) + +erbsout = erbs_driesse(ghi, solar_zenith, dni_extra=dni_extra) + +dni = erbsout['dni'] +dhi = erbsout['dhi'] + +irrads = get_total_irradiance(surface_tilt, surface_azimuth, + solar_zenith, solar_azimuth, + dni, ghi, dhi, + dni_extra, + model='perez-driesse') + +poa_global = irrads['poa_global'] + +# %% +# +# Suppose you measure that POA global is 200 W/m2. What would GHI be? +# + +poa_test = 200 + +ghi_hat = ghi_from_poa_driesse_2023(surface_tilt, surface_azimuth, + solar_zenith, solar_azimuth, + poa_test, + dni_extra, + full_output=False) + +print('Estimated GHI: %.2f W/m².' % ghi_hat) + +# %% +# +# Show this result on the graph of all possible combinations of GHI and POA. +# + +plt.figure() +plt.plot(ghi, poa_global, 'k-') +plt.axvline(ghi_hat, color='g', lw=1) +plt.axhline(poa_test, color='g', lw=1) +plt.plot(ghi_hat, poa_test, 'gs') +plt.annotate('GHI=%.2f' % (ghi_hat), + xy=(ghi_hat-2, 200+2), + xytext=(ghi_hat-20, 200+20), + ha='right', + arrowprops={'arrowstyle': 'simple'}) +plt.xlim(0, 500) +plt.ylim(0, 250) +plt.xlabel('GHI [W/m²]') +plt.ylabel('POA [W/m²]') +plt.show() + +# %% +# +# Now change the solar azimuth to match the conditions for figure 4 in [1]_. +# + +solar_azimuth = 76 + +# %% +# +# Again, estimate DNI and DHI using the Erbs-Driesse model, then +# transpose using the Perez-Driesse model. +# + +erbsout = erbs_driesse(ghi, solar_zenith, dni_extra=dni_extra) + +dni = erbsout['dni'] +dhi = erbsout['dhi'] + +irrads = get_total_irradiance(surface_tilt, surface_azimuth, + solar_zenith, solar_azimuth, + dni, ghi, dhi, + dni_extra, + model='perez-driesse') + +poa_global = irrads['poa_global'] + +# %% +# +# Now reverse transpose all the POA values and observe that the original +# GHI cannot always be found. There is a range of POA values that +# maps to three possible GHI values, and there is not enough information +# to choose one of them. Sometimes we get lucky and the right one comes +# out, other times not. +# + +result = ghi_from_poa_driesse_2023(surface_tilt, surface_azimuth, + solar_zenith, solar_azimuth, + poa_global, + dni_extra, + full_output=True, + ) + +ghi_hat, conv, niter = result +correct = np.isclose(ghi, ghi_hat, atol=0.01) + +plt.figure() +plt.plot(np.where(correct, ghi, np.nan), np.where(correct, poa_global, np.nan), + 'g.', label='correct GHI found') +plt.plot(ghi[~correct], poa_global[~correct], 'r.', label='unreachable GHI') +plt.plot(ghi[~conv], poa_global[~conv], 'm.', label='out of range (kt > 1.25)') +plt.axhspan(88, 103, color='y', alpha=0.25, label='problem region') + +plt.xlim(0, 500) +plt.ylim(0, 250) +plt.xlabel('GHI [W/m²]') +plt.ylabel('POA [W/m²]') +plt.legend() +plt.show() diff --git a/docs/examples/irradiance-transposition/plot_rtranpose_year.py b/docs/examples/irradiance-transposition/plot_rtranpose_year.py new file mode 100644 index 0000000000..c0b9860bba --- /dev/null +++ b/docs/examples/irradiance-transposition/plot_rtranpose_year.py @@ -0,0 +1,152 @@ +""" +Reverse transposition using one year of hourly data +=================================================== + +With a brief look at accuracy and speed. + +Author: Anton Driesse + +""" +# %% +# +# Introduction +# ------------ +# When irradiance is measured on a tilted plane, it is useful to be able to +# estimate the GHI that produces the POA irradiance. +# The estimation requires inverting a GHI-to-POA irradiance model, +# which involves two parts: +# a decomposition of GHI into direct and diffuse components, +# and a transposition model that calculates the direct and diffuse +# irradiance on the tilted plane. +# Recovering GHI from POA irradiance is termed "reverse transposition." +# +# In this example we start with a TMY file and calculate POA global irradiance. +# Then we use :py:meth:`pvlib.irradiance.ghi_from_poa_driesse_2023` to estimate +# the original GHI from POA global. Details of the method found in [1]_. +# +# Another method for reverse tranposition called GTI-DIRINT is also +# available in pvlib python (:py:meth:`pvlib.irradiance.gti_dirint`). +# More information is available in [2]_. +# +# References +# ---------- +# .. [1] Driesse, A., Jensen, A., Perez, R., 2024. A Continuous form of the +# Perez diffuse sky model for forward and reverse transposition. +# Solar Energy vol. 267. :doi:`10.1016/j.solener.2023.112093` +# +# .. [2] B. Marion, A model for deriving the direct normal and +# diffuse horizontal irradiance from the global tilted +# irradiance, Solar Energy 122, 1037-1046. +# :doi:`10.1016/j.solener.2015.10.024` + +import os +import time +import pandas as pd + +import matplotlib.pyplot as plt + +import pvlib +from pvlib import iotools, location +from pvlib.irradiance import (get_extra_radiation, + get_total_irradiance, + ghi_from_poa_driesse_2023, + aoi, + ) + +# %% +# +# Read a TMY3 file containing weather data and select needed columns. +# + +PVLIB_DIR = pvlib.__path__[0] +DATA_FILE = os.path.join(PVLIB_DIR, 'data', '723170TYA.CSV') + +tmy, metadata = iotools.read_tmy3(DATA_FILE, coerce_year=1990, + map_variables=True) + +df = pd.DataFrame({'ghi': tmy['ghi'], 'dhi': tmy['dhi'], 'dni': tmy['dni'], + 'temp_air': tmy['temp_air'], + 'wind_speed': tmy['wind_speed'], + }) + +# %% +# +# Shift the timestamps to the middle of the hour and calculate sun positions. +# + +df.index = df.index - pd.Timedelta(minutes=30) + +loc = location.Location.from_tmy(metadata) +solpos = loc.get_solarposition(df.index) + +# %% +# +# Estimate global irradiance on a fixed-tilt array (forward transposition). +# The array is tilted 30 degrees and oriented 30 degrees east of south. +# + +TILT = 30 +ORIENT = 150 + +df['dni_extra'] = get_extra_radiation(df.index) + +total_irrad = get_total_irradiance(TILT, ORIENT, + solpos.apparent_zenith, + solpos.azimuth, + df.dni, df.ghi, df.dhi, + dni_extra=df.dni_extra, + model='perez-driesse') + +df['poa_global'] = total_irrad.poa_global +df['aoi'] = aoi(TILT, ORIENT, solpos.apparent_zenith, solpos.azimuth) + +# %% +# +# Now estimate ghi from poa_global using reverse transposition. +# The algorithm uses a simple bisection search, which is quite slow +# because scipy doesn't offer a vectorized version (yet). +# For this reason we'll process a random sample of 1000 timestamps +# rather than the whole year. +# + +df = df[df.ghi > 0].sample(n=1000) +solpos = solpos.reindex(df.index) + +start = time.process_time() + +df['ghi_rev'] = ghi_from_poa_driesse_2023(TILT, ORIENT, + solpos.apparent_zenith, + solpos.azimuth, + df.poa_global, + dni_extra=df.dni_extra) +finish = time.process_time() + +print('Elapsed time for reverse transposition: %.1f s' % (finish - start)) + +# %% +# +# This graph shows the reverse transposed values vs. the original values. +# The markers are color-coded by angle-of-incidence to show that +# errors occur primarily with incidence angle approaching 90° and beyond. +# +# Note that the results look particularly good because the POA values +# were calculated using the same models as used in reverse transposition. +# This isn't cheating though. It's a way of ensuring that the errors +# we see are really due to the reverse transposition algorithm. +# Expect to see larger errors with real-word POA measurements +# because errors from forward and reverse transposition will both be present. +# + +df = df.sort_values('aoi') + +plt.figure() +plt.gca().grid(True, alpha=.5) +pc = plt.scatter(df['ghi'], df['ghi_rev'], c=df['aoi'], s=15, + cmap='jet', vmin=60, vmax=120) +plt.colorbar(label='AOI [°]') +pc.set_alpha(0.5) + +plt.xlabel('GHI original [W/m²]') +plt.ylabel('GHI from POA [W/m²]') + +plt.show() diff --git a/docs/sphinx/source/reference/irradiance/transposition.rst b/docs/sphinx/source/reference/irradiance/transposition.rst index 7b3624e692..22136f0c58 100644 --- a/docs/sphinx/source/reference/irradiance/transposition.rst +++ b/docs/sphinx/source/reference/irradiance/transposition.rst @@ -15,3 +15,4 @@ Transposition models irradiance.klucher irradiance.reindl irradiance.king + irradiance.ghi_from_poa_driesse_2023 diff --git a/docs/sphinx/source/whatsnew/v0.10.3.rst b/docs/sphinx/source/whatsnew/v0.10.3.rst index 007eb8d34b..63f7b90b84 100644 --- a/docs/sphinx/source/whatsnew/v0.10.3.rst +++ b/docs/sphinx/source/whatsnew/v0.10.3.rst @@ -9,6 +9,9 @@ Enhancements ~~~~~~~~~~~~ * Added the continuous Perez-Driesse transposition model. :py:func:`pvlib.irradiance.perez_driesse` (:issue:`1841`, :pull:`1876`) +* Added a reverse transposition algorithm using the Perez-Driesse model. + :py:func:`pvlib.irradiance.ghi_from_poa_driesse_2023` + (:issue:`1901`, :pull:`1907`) * :py:func:`pvlib.bifacial.infinite_sheds.get_irradiance` and :py:func:`pvlib.bifacial.infinite_sheds.get_irradiance_poa` now include shaded fraction in returned variables. (:pull:`1871`) @@ -27,8 +30,10 @@ Documentation ~~~~~~~~~~~~~ * Create :ref:`weatherdata` User's Guide page. (:pull:`1754`) * Fixed a plotting issue in the IV curve gallery example (:pull:`1895`) +* Added two examples to demonstrate reverse transposition (:pull:`1907`) * Fixed `detect_clearsky` example in `clearsky.rst` (:issue:`1914`) + Requirements ~~~~~~~~~~~~ * Minimum version of scipy advanced from 1.4.0 to 1.5.0. (:issue:`1918`, :pull:`1919`) diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py index efae7f6236..f5e9a44b7a 100644 --- a/pvlib/irradiance.py +++ b/pvlib/irradiance.py @@ -11,6 +11,7 @@ import numpy as np import pandas as pd from scipy.interpolate import splev +from scipy.optimize import bisect from pvlib import atmosphere, solarposition, tools import pvlib # used to avoid dni name collision in complete_irradiance @@ -1380,9 +1381,9 @@ def perez_driesse(surface_tilt, surface_azimuth, dhi, dni, dni_extra, References ---------- - .. [1] A. Driesse, A. Jensen, R. Perez, A Continuous Form of the Perez - Diffuse Sky Model for Forward and Reverse Transposition, accepted - for publication in the Solar Energy Journal. + .. [1] Driesse, A., Jensen, A., Perez, R., 2024. A Continuous form of the + Perez diffuse sky model for forward and reverse transposition. + Solar Energy vol. 267. :doi:`10.1016/j.solener.2023.112093` .. [2] Perez, R., Ineichen, P., Seals, R., Michalsky, J., Stewart, R., 1990. Modeling daylight availability and irradiance components from @@ -1444,6 +1445,170 @@ def perez_driesse(surface_tilt, surface_azimuth, dhi, dni, dni_extra, return sky_diffuse +def _poa_from_ghi(surface_tilt, surface_azimuth, + solar_zenith, solar_azimuth, + ghi, + dni_extra, airmass, albedo): + ''' + Transposition function that includes decomposition of GHI using the + continuous Erbs-Driesse model. + + Helper function for ghi_from_poa_driesse_2023. + ''' + # Contributed by Anton Driesse (@adriesse), PV Performance Labs. Nov., 2023 + + erbsout = erbs_driesse(ghi, solar_zenith, dni_extra=dni_extra) + + dni = erbsout['dni'] + dhi = erbsout['dhi'] + + irrads = get_total_irradiance(surface_tilt, surface_azimuth, + solar_zenith, solar_azimuth, + dni, ghi, dhi, + dni_extra, airmass, albedo, + model='perez-driesse') + + return irrads['poa_global'] + + +def _ghi_from_poa(surface_tilt, surface_azimuth, + solar_zenith, solar_azimuth, + poa_global, + dni_extra, airmass, albedo, + xtol=0.01): + ''' + Reverse transposition function that uses the scalar bisection from scipy. + + Helper function for ghi_from_poa_driesse_2023. + ''' + # Contributed by Anton Driesse (@adriesse), PV Performance Labs. Nov., 2023 + + # propagate nans and zeros quickly + if np.isnan(poa_global): + return np.nan, False, 0 + if poa_global <= 0: + return 0.0, True, 0 + + # function whose root needs to be found + def poa_error(ghi): + poa_hat = _poa_from_ghi(surface_tilt, surface_azimuth, + solar_zenith, solar_azimuth, + ghi, + dni_extra, airmass, albedo) + return poa_hat - poa_global + + # calculate an upper bound for ghi using clearness index 1.25 + ghi_clear = dni_extra * tools.cosd(solar_zenith) + ghi_high = np.maximum(10, 1.25 * ghi_clear) + + try: + result = bisect(poa_error, + a=0, + b=ghi_high, + xtol=xtol, + maxiter=25, + full_output=True, + disp=False, + ) + except ValueError: + # this occurs when poa_error has the same sign at both end points + ghi = np.nan + conv = False + niter = -1 + else: + ghi = result[0] + conv = result[1].converged + niter = result[1].iterations + + return ghi, conv, niter + + +def ghi_from_poa_driesse_2023(surface_tilt, surface_azimuth, + solar_zenith, solar_azimuth, + poa_global, + dni_extra=None, airmass=None, albedo=0.25, + xtol=0.01, + full_output=False): + ''' + Estimate global horizontal irradiance (GHI) from global plane-of-array + (POA) irradiance. This reverse transposition algorithm uses a bisection + search together with the continuous Perez-Driesse transposition and + continuous Erbs-Driesse decomposition models, as described in [1]_. + + Parameters + ---------- + surface_tilt : numeric + Panel tilt from horizontal. [degree] + surface_azimuth : numeric + Panel azimuth from north. [degree] + solar_zenith : numeric + Solar zenith angle. [degree] + solar_azimuth : numeric + Solar azimuth angle. [degree] + poa_global : numeric + Plane-of-array global irradiance, aka global tilted irradiance. [W/m^2] + dni_extra : None or numeric, default None + Extraterrestrial direct normal irradiance. [W/m^2] + airmass : None or numeric, default None + Relative airmass (not adjusted for pressure). [unitless] + albedo : numeric, default 0.25 + Ground surface albedo. [unitless] + xtol : numeric, default 0.01 + Convergence criterion. The estimated GHI will be within xtol of the + true value. [W/m^2] + full_output : boolean, default False + If full_output is False, only ghi is returned, otherwise the return + value is (ghi, converged, niter). (see Returns section for details). + + Returns + ------- + ghi : numeric + Estimated GHI. [W/m^2] + converged : boolean, optional + Present if full_output=True. Indicates which elements converged + successfully. + niter : integer, optional + Present if full_output=True. Indicates how many bisection iterations + were done. + + Notes + ----- + Since :py:func:`scipy.optimize.bisect` is not vectorized, high-resolution + time series can be quite slow to process. + + References + ---------- + .. [1] Driesse, A., Jensen, A., Perez, R., 2024. A Continuous form of the + Perez diffuse sky model for forward and reverse transposition. + Solar Energy vol. 267. :doi:`10.1016/j.solener.2023.112093` + + See also + -------- + perez_driesse + erbs_driesse + gti_dirint + ''' + # Contributed by Anton Driesse (@adriesse), PV Performance Labs. Nov., 2023 + + ghi_from_poa_array = np.vectorize(_ghi_from_poa) + + ghi, conv, niter = ghi_from_poa_array(surface_tilt, surface_azimuth, + solar_zenith, solar_azimuth, + poa_global, + dni_extra, airmass, albedo, + xtol=0.01) + + if isinstance(poa_global, pd.Series): + ghi = pd.Series(ghi, poa_global.index) + conv = pd.Series(conv, poa_global.index) + niter = pd.Series(niter, poa_global.index) + + if full_output: + return ghi, conv, niter + else: + return ghi + + def clearsky_index(ghi, clearsky_ghi, max_clearsky_index=2.0): """ Calculate the clearsky index. @@ -2567,8 +2732,9 @@ def erbs_driesse(ghi, zenith, datetime_or_doy=None, dni_extra=None, References ---------- - .. [1] A. Driesse, A. Jensen, R. Perez, A Continuous Form of the Perez - Diffuse Sky Model for Forward and Reverse Transposition, forthcoming. + .. [1] Driesse, A., Jensen, A., Perez, R., 2024. A Continuous form of the + Perez diffuse sky model for forward and reverse transposition. + Solar Energy vol. 267. :doi:`10.1016/j.solener.2023.112093` .. [2] D. G. Erbs, S. A. Klein and J. A. Duffie, Estimation of the diffuse radiation fraction for hourly, daily and monthly-average diff --git a/pvlib/tests/test_irradiance.py b/pvlib/tests/test_irradiance.py index ba66f4dc36..55fef490ba 100644 --- a/pvlib/tests/test_irradiance.py +++ b/pvlib/tests/test_irradiance.py @@ -782,6 +782,50 @@ def test_dirint_min_cos_zenith_max_zenith(): assert_series_equal(out, expected, check_less_precise=True) +def test_ghi_from_poa_driesse(): + # inputs copied from test_gti_dirint + times = pd.DatetimeIndex( + ['2014-06-24T06-0700', '2014-06-24T09-0700', '2014-06-24T12-0700']) + poa_global = np.array([20, 300, 1000]) + zenith = np.array([80, 45, 20]) + azimuth = np.array([90, 135, 180]) + surface_tilt = 30 + surface_azimuth = 180 + + # test core function + output = irradiance.ghi_from_poa_driesse_2023( + surface_tilt, surface_azimuth, zenith, azimuth, + poa_global, dni_extra=1366.1) + + expected = [22.089, 304.088, 931.143] + assert_allclose(expected, output, atol=0.001) + + # test series output + poa_global = pd.Series([20, 300, 1000], index=times) + + output = irradiance.ghi_from_poa_driesse_2023( + surface_tilt, surface_azimuth, zenith, azimuth, + poa_global, dni_extra=1366.1) + + assert isinstance(output, pd.Series) + + # test full_output option and special cases + poa_global = np.array([0, 1500, np.nan]) + + ghi, conv, niter = irradiance.ghi_from_poa_driesse_2023( + surface_tilt, surface_azimuth, zenith, azimuth, + poa_global, dni_extra=1366.1, full_output=True) + + expected = [0, np.nan, np.nan] + assert_allclose(expected, ghi, atol=0.001) + + expected = [True, False, False] + assert_allclose(expected, conv) + + expected = [0, -1, 0] + assert_allclose(expected, niter) + + def test_gti_dirint(): times = pd.DatetimeIndex( ['2014-06-24T06-0700', '2014-06-24T09-0700', '2014-06-24T12-0700'])