ADR PV efficiency model from pvpltools-python (take 2) (#1602)

* Add ADR model and fitting function.

* Add tests.

* Add examples.

* Register API changes.

* Fix bullet lists in example.

* Move new model to new module pvarray.py.

* Add another gallery example and tweak fitting function.

* Adjust expected values in tests and optimistically tighten the tolerances, plus make stickler happy again.

* Roll back my optimism on those test tolerances again.

* Back off the tolelerance one more time to get the last three environments to pass.

* Apply suggestions from code review

Kevin's bug observations and identifications.  :)

Co-authored-by: Kevin Anderson <[email protected]>

* Improvements after first review.

* Apply suggestions from code review

Co-authored-by: Cliff Hansen <[email protected]>

* Improvements after second review.

* Stickler, and various other small mostly appearance-related changes.

* Changes after documentation review.

Co-authored-by: Kevin Anderson <[email protected]>
Co-authored-by: Cliff Hansen <[email protected]>

Author: 3 people

docs/examples/adr-pvarray/README.rst

@@ -0,0 +1,3 @@
+ADR Model for PV Module Efficiency
+----------------------------------
+

docs/examples/adr-pvarray/plot_fit_to_matrix.py

@@ -0,0 +1,108 @@
+"""
+Obtaining ADR model parameters from IEC 61853 matrix measurements
+=================================================================
+
+There's a fitting function provided in pvlib to do exactly that.
+
+Since PV module efficiency varies with irradiance and temperature
+what better way to train a model than using efficiency measurement
+over a broad range of temperature and irradiance levels?
+The standard IEC 61853-1 defines a standard matrix of conditions
+for such measurements and this example shows how the ADR model
+parameters can be determined with just a few lines of code using
+functions in pvlib-python.
+
+Author: Anton Driesse
+"""
+
+from io import StringIO
+import pandas as pd
+import matplotlib.pyplot as plt
+
+from pvlib.pvarray import pvefficiency_adr, fit_pvefficiency_adr
+
+# %% The text on this line is not displayed
+#
+# Here are some matrix measurements:
+#
+
+iec61853data = '''
+    irradiance  temperature     p_mp
+0          100         15.0   30.159
+1          200         15.0   63.057
+2          400         15.0  129.849
+3          600         15.0  197.744
+4          800         15.0  264.825
+5         1000         15.0  330.862
+6          100         25.0   29.250
+7          200         25.0   61.137
+8          400         25.0  126.445
+9          600         25.0  192.278
+10         800         25.0  257.561
+11        1000         25.0  322.305
+12        1100         25.0  354.174
+15         400         50.0  117.062
+16         600         50.0  177.959
+17         800         50.0  238.626
+18        1000         50.0  298.954
+19        1100         50.0  328.413
+23         600         75.0  162.966
+24         800         75.0  218.585
+25        1000         75.0  273.651
+26        1100         75.0  301.013
+'''
+df = pd.read_csv(StringIO(iec61853data), delim_whitespace=True)
+
+# %%
+#
+# Now calculate the normalized or relative efficiency values
+# and use the fitting function to determine the parameters.
+# The parameters (shown below) can now be used to
+# simulate the module operating in a PV system.
+#
+
+P_REF = 322.305   # (W) STC value from the table above
+G_REF = 1000.     # (W/m2)
+
+df['eta_rel'] = (df['p_mp'] / P_REF) / (df['irradiance'] / G_REF)
+
+adr_params = fit_pvefficiency_adr(df['irradiance'], df['temperature'],
+                                  df['eta_rel'])
+
+for k, v in adr_params.items():
+    print('%-5s = %8.5f' % (k, v))
+
+# %%
+#
+# Compare the model output to the original measurements.
+# The chart below shows minor differences but due to their random nature
+# they are most likely evidence of the limitations of measurement accuracy.
+#
+
+eta_rel_adr = pvefficiency_adr(df['irradiance'],
+                               df['temperature'], **adr_params)
+
+plt.figure()
+plt.plot(df['irradiance'], df['eta_rel'], 'oc', ms=8)
+plt.plot(df['irradiance'], eta_rel_adr, '.k')
+plt.legend(['Lab measurements', 'ADR model fit'], loc='lower right')
+plt.xlabel('Irradiance [W/m²]')
+plt.ylabel('Relative efficiency [-]')
+plt.grid(alpha=0.5)
+plt.xlim(0, 1200)
+plt.ylim(0.7, 1.1)
+plt.show()
+
+# %%
+#
+# References
+# ----------
+# .. [1] A. Driesse and J. S. Stein, "From IEC 61853 power measurements
+#    to PV system simulations", Sandia Report No. SAND2020-3877, 2020.
+#    :doi:`10.2172/1615179`
+#
+# .. [2] A. Driesse, M. Theristis and J. S. Stein, "A New Photovoltaic Module
+#    Efficiency Model for Energy Prediction and Rating," in IEEE Journal
+#    of Photovoltaics, vol. 11, no. 2, pp. 527-534, March 2021.
+#    :doi:`10.1109/JPHOTOV.2020.3045677`
+#

docs/examples/adr-pvarray/plot_simulate_fast.py

@@ -0,0 +1,176 @@
+"""
+Fast simulation using the ADR efficiency model starting from PVsyst parameters
+==============================================================================
+
+Would you like to increase simulation speed by a factor of 4000+?
+
+Simulation using single-diode models can be slow because the maximum
+power point is usually found by an iterative search.
+In this example we use the PVsyst single diode model to generate
+a matrix of efficiency values, then determine the ADR model
+parameters to approximate the behavior of the PVsyst model.
+This way both PVsyst and ADR models can simulate the same PV module type.
+
+To compare simulation speed, we run them using ``timeit``.
+
+Author: Anton Driesse
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+from pvlib.pvsystem import calcparams_pvsyst, max_power_point
+from pvlib.pvarray import fit_pvefficiency_adr, pvefficiency_adr
+
+from timeit import timeit
+
+# %% The text on this line is not displayed
+#
+# Generate a matrix of power values
+#
+
+pvsyst_params = {'alpha_sc': 0.0015,
+                 'gamma_ref': 1.20585,
+                 'mu_gamma': -9.41066e-05,
+                 'I_L_ref': 5.9301,
+                 'I_o_ref': 2.9691e-10,
+                 'R_sh_ref': 1144,
+                 'R_sh_0': 3850,
+                 'R_s': 0.6,
+                 'cells_in_series': 96,
+                 'R_sh_exp': 5.5,
+                 'EgRef': 1.12,
+                 }
+
+G_REF = 1000
+T_REF = 25
+
+params_stc = calcparams_pvsyst(G_REF, T_REF, **pvsyst_params)
+mpp_stc = max_power_point(*params_stc)
+
+P_REF = mpp_stc['p_mp']
+
+g, t = np.meshgrid(np.linspace(100, 1100, 11),
+                   np.linspace(0, 75, 4))
+
+adjusted_params = calcparams_pvsyst(g, t, **pvsyst_params)
+mpp = max_power_point(*adjusted_params)
+p_mp = mpp['p_mp']
+
+print('irradiance')
+print(g[:1].round(0))
+
+print('maximum power')
+print(p_mp.round(1))
+
+# %%
+#
+# Convert power matrix to efficiency and fit the ADR model to all the points
+#
+
+eta_rel_pvs = (p_mp / P_REF) / (g / G_REF)
+
+adr_params = fit_pvefficiency_adr(g, t, eta_rel_pvs, dict_output=True)
+
+for k, v in adr_params.items():
+    print('%-5s = %8.5f' % (k, v))
+
+# %%
+#
+# Compare the ADR model output to the PVsyst model output
+#
+
+eta_rel_adr = pvefficiency_adr(g, t, **adr_params)
+mbe = np.mean(eta_rel_adr - eta_rel_pvs)
+rmse = np.sqrt(np.mean(np.square(eta_rel_adr - eta_rel_pvs)))
+
+plt.figure()
+plt.plot(g.flat, eta_rel_pvs.flat, 'oc', ms=8)
+plt.plot(g.flat, eta_rel_adr.flat, '.k')
+plt.grid(alpha=0.5)
+plt.xlim(0, 1200)
+plt.ylim(0.7, 1.1)
+
+plt.xlabel('Irradiance [W/m²]')
+plt.ylabel('Relative efficiency [-]')
+plt.legend(['PVsyst model output', 'ADR model fit'], loc='lower right')
+plt.title('Differences: mean %.5f, RMS %.5f' % (mbe, rmse))
+plt.show()
+
+# %%
+#
+# Generate some random irradiance and temperature data
+#
+
+g = np.random.uniform(0, 1200, 8760)
+t = np.random.uniform(20,  80, 8760)
+
+
+def run_adr():
+    eta_rel = pvefficiency_adr(g, t, **adr_params)
+    p_adr = P_REF * eta_rel * (g / G_REF)
+    return p_adr
+
+
+def run_pvsyst():
+    adjusted_params = calcparams_pvsyst(g, t, **pvsyst_params)
+    mpp = max_power_point(*adjusted_params)
+    p_pvs = mpp['p_mp']
+    return p_pvs
+
+
+elapsed_adr = timeit('run_adr()', number=1, globals=globals())
+elapsed_pvs = timeit('run_pvsyst()', number=1, globals=globals())
+
+print('Elapsed time for the PVsyst model: %9.6f s' % elapsed_pvs)
+print('Elapsed time for the ADR    model: %9.6f s' % elapsed_adr)
+print('ADR acceleration ratio:           %9.0f x' % (elapsed_pvs/elapsed_adr))
+
+# %%
+#
+# That's fast, but is it accurate?
+# Run them again to compare the simulated power values
+#
+
+p_pvs = run_pvsyst()
+p_adr = run_adr()
+
+mbe = np.mean(p_adr - p_pvs)
+rmse = np.sqrt(np.mean(np.square(p_adr - p_pvs)))
+
+# sphinx_gallery_thumbnail_number = 2
+plt.figure()
+pc = plt.scatter(p_pvs, p_adr-p_pvs, c=t, cmap='jet')
+plt.colorbar()
+pc.set_alpha(0.25)
+plt.ylim(-1.4, 1.4)
+plt.grid(alpha=0.5)
+
+plt.xlabel('Power calculated using the PVsyst model [W]')
+plt.ylabel('ADR model power - PVsyst model power [W]')
+plt.title('Differences: mean %.2f W, RMS %.2f W' % (mbe, rmse))
+plt.show()
+
+# %%
+#
+# There are some small systematic differences between the original
+# PVsyst model output and the ADR fit.  But these differences are
+# much smaller than the typical uncertainty in measured output
+# of modules of this type. The PVsyst model and the parameters
+# we started with are of course also only approximations of the
+# true module behavior.
+#
+
+# %%
+#
+# References
+# ----------
+# .. [1] A. Driesse and J. S. Stein, "From IEC 61853 power measurements
+#    to PV system simulations", Sandia Report No. SAND2020-3877, 2020.
+#    :doi:`10.2172/1615179`
+#
+# .. [2] A. Driesse, M. Theristis and J. S. Stein, "A New Photovoltaic Module
+#    Efficiency Model for Energy Prediction and Rating," in IEEE Journal
+#    of Photovoltaics, vol. 11, no. 2, pp. 527-534, March 2021.
+#    :doi:`10.1109/JPHOTOV.2020.3045677`
+#