Add recombination current params to all bishop88 functions (pvlib#763)

* Pass through extra recombination parameters to allow use in max_power_point().

* Add recombination current parameters to the other bishop88 functions.

* Update whatsnew.

* Add tests.

* Tests for both newton and brentq, and some formatting.

* Delete blank line.

* Add scipy requirement triggered by brentq tests.

* Add blank line.

* docstring edits

* redo docstring edits

* add parameterize

* remove blank lines

Author: adriesse

docs/sphinx/source/whatsnew/v0.7.0.rst

@@ -14,6 +14,8 @@ Enhancements
 ~~~~~~~~~~~~
 * Created two new incidence angle modifier functions: :py:func:`pvlib.pvsystem.iam_martin_ruiz`
   and :py:func:`pvlib.pvsystem.iam_interp`. (:issue:`751`)
+* Added recombination current parameters to bishop88 single-diode functions and also
+  to :py:func:`pvlib.pvsystem.max_power_point`. (:issue:`762`)
 
 Bug fixes
 ~~~~~~~~~

pvlib/pvsystem.py

@@ -2473,7 +2473,8 @@ def singlediode(photocurrent, saturation_current, resistance_series,
 
 
 def max_power_point(photocurrent, saturation_current, resistance_series,
-                    resistance_shunt, nNsVth, method='brentq'):
+                    resistance_shunt, nNsVth, d2mutau=0, NsVbi=np.Inf,
+                    method='brentq'):
     """
     Given the single diode equation coefficients, calculates the maximum power
     point (MPP).
@@ -2491,6 +2492,17 @@ def max_power_point(photocurrent, saturation_current, resistance_series,
     nNsVth : numeric
         product of thermal voltage ``Vth`` [V], diode ideality factor ``n``,
         and number of serices cells ``Ns``
+    d2mutau : numeric, default 0
+        PVsyst parameter for cadmium-telluride (CdTe) and amorphous-silicon
+        (a-Si) modules that accounts for recombination current in the
+        intrinsic layer. The value is the ratio of intrinsic layer thickness
+        squared :math:`d^2` to the diffusion length of charge carriers
+        :math:`\\mu \\tau`. [V]
+    NsVbi : numeric, default np.inf
+        PVsyst parameter for cadmium-telluride (CdTe) and amorphous-silicon
+        (a-Si) modules that is the product of the PV module number of series
+        cells ``Ns`` and the builtin voltage ``Vbi`` of the intrinsic layer.
+        [V].
     method : str
         either ``'newton'`` or ``'brentq'``
 
@@ -2508,7 +2520,8 @@ def max_power_point(photocurrent, saturation_current, resistance_series,
     """
     i_mp, v_mp, p_mp = _singlediode.bishop88_mpp(
         photocurrent, saturation_current, resistance_series,
-        resistance_shunt, nNsVth, method=method.lower()
+        resistance_shunt, nNsVth, d2mutau=0, NsVbi=np.Inf,
+        method=method.lower()
     )
     if isinstance(photocurrent, pd.Series):
         ivp = {'i_mp': i_mp, 'v_mp': v_mp, 'p_mp': p_mp}

pvlib/singlediode.py

@@ -94,14 +94,17 @@ def bishop88(diode_voltage, photocurrent, saturation_current,
     nNsVth : numeric
         product of thermal voltage ``Vth`` [V], diode ideality factor ``n``,
         and number of series cells ``Ns``
-    d2mutau : numeric
-        PVSyst thin-film recombination parameter that is the ratio of thickness
-        of the intrinsic layer squared :math:`d^2` and the diffusion length of
-        charge carriers :math:`\\mu \\tau`, in volts [V], defaults to 0[V]
-    NsVbi : numeric
-        PVSyst thin-film recombination parameter that is the product of the PV
-        module number of series cells ``Ns`` and the builtin voltage ``Vbi`` of
-        the intrinsic layer, in volts [V], defaults to ``np.inf``
+    d2mutau : numeric, default 0
+        PVsyst parameter for cadmium-telluride (CdTe) and amorphous-silicon
+        (a-Si) modules that accounts for recombination current in the
+        intrinsic layer. The value is the ratio of intrinsic layer thickness
+        squared :math:`d^2` to the diffusion length of charge carriers
+        :math:`\\mu \\tau`. [V]
+    NsVbi : numeric, default np.inf
+        PVsyst parameter for cadmium-telluride (CdTe) and amorphous-silicon
+        (a-Si) modules that is the product of the PV module number of series
+        cells ``Ns`` and the builtin voltage ``Vbi`` of the intrinsic layer.
+        [V].
     gradients : bool
         False returns only I, V, and P. True also returns gradients
 
@@ -116,8 +119,8 @@ def bishop88(diode_voltage, photocurrent, saturation_current,
     Notes
     -----
     The PVSyst thin-film recombination losses parameters ``d2mutau`` and
-    ``NsVbi`` are only applied to cadmium-telluride (CdTe) and amorphous-
-    silicon (a:Si) PV modules, [2]_, [3]_. The builtin voltage :math:`V_{bi}`
+    ``NsVbi`` should only be applied to cadmium-telluride (CdTe) and amorphous-
+    silicon (a-Si) PV modules, [2]_, [3]_. The builtin voltage :math:`V_{bi}`
     should account for all junctions. For example: tandem and triple junction
     cells would have builtin voltages of 1.8[V] and 2.7[V] respectively, based
     on the default of 0.9[V] for a single junction. The parameter ``NsVbi``
@@ -173,7 +176,7 @@ def bishop88(diode_voltage, photocurrent, saturation_current,
 
 def bishop88_i_from_v(voltage, photocurrent, saturation_current,
                       resistance_series, resistance_shunt, nNsVth,
-                      method='newton'):
+                      d2mutau=0, NsVbi=np.Inf, method='newton'):
     """
     Find current given any voltage.
 
@@ -192,6 +195,17 @@ def bishop88_i_from_v(voltage, photocurrent, saturation_current,
     nNsVth : numeric
         product of diode ideality factor (n), number of series cells (Ns), and
         thermal voltage (Vth = k_b * T / q_e) in volts [V]
+    d2mutau : numeric, default 0
+        PVsyst parameter for cadmium-telluride (CdTe) and amorphous-silicon
+        (a-Si) modules that accounts for recombination current in the
+        intrinsic layer. The value is the ratio of intrinsic layer thickness
+        squared :math:`d^2` to the diffusion length of charge carriers
+        :math:`\\mu \\tau`. [V]
+    NsVbi : numeric, default np.inf
+        PVsyst parameter for cadmium-telluride (CdTe) and amorphous-silicon
+        (a-Si) modules that is the product of the PV module number of series
+        cells ``Ns`` and the builtin voltage ``Vbi`` of the intrinsic layer.
+        [V].
     method : str
         one of two optional search methods: either ``'brentq'``, a reliable and
         bounded method or ``'newton'`` which is the default.
@@ -203,7 +217,7 @@ def bishop88_i_from_v(voltage, photocurrent, saturation_current,
     """
     # collect args
     args = (photocurrent, saturation_current, resistance_series,
-            resistance_shunt, nNsVth)
+            resistance_shunt, nNsVth, d2mutau, NsVbi)
 
     def fv(x, v, *a):
         # calculate voltage residual given diode voltage "x"
@@ -216,8 +230,9 @@ def fv(x, v, *a):
         # brentq only works with scalar inputs, so we need a set up function
         # and np.vectorize to repeatedly call the optimizer with the right
         # arguments for possible array input
-        def vd_from_brent(voc, v, iph, isat, rs, rsh, gamma):
-            return brentq(fv, 0.0, voc, args=(v, iph, isat, rs, rsh, gamma))
+        def vd_from_brent(voc, v, iph, isat, rs, rsh, gamma, d2mutau, NsVbi):
+            return brentq(fv, 0.0, voc,
+                          args=(v, iph, isat, rs, rsh, gamma, d2mutau, NsVbi))
 
         vd_from_brent_vectorized = np.vectorize(vd_from_brent)
         vd = vd_from_brent_vectorized(voc_est, voltage, *args)
@@ -235,7 +250,7 @@ def vd_from_brent(voc, v, iph, isat, rs, rsh, gamma):
 
 def bishop88_v_from_i(current, photocurrent, saturation_current,
                       resistance_series, resistance_shunt, nNsVth,
-                      method='newton'):
+                      d2mutau=0, NsVbi=np.Inf, method='newton'):
     """
     Find voltage given any current.
 
@@ -254,6 +269,17 @@ def bishop88_v_from_i(current, photocurrent, saturation_current,
     nNsVth : numeric
         product of diode ideality factor (n), number of series cells (Ns), and
         thermal voltage (Vth = k_b * T / q_e) in volts [V]
+    d2mutau : numeric, default 0
+        PVsyst parameter for cadmium-telluride (CdTe) and amorphous-silicon
+        (a-Si) modules that accounts for recombination current in the
+        intrinsic layer. The value is the ratio of intrinsic layer thickness
+        squared :math:`d^2` to the diffusion length of charge carriers
+        :math:`\\mu \\tau`. [V]
+    NsVbi : numeric, default np.inf
+        PVsyst parameter for cadmium-telluride (CdTe) and amorphous-silicon
+        (a-Si) modules that is the product of the PV module number of series
+        cells ``Ns`` and the builtin voltage ``Vbi`` of the intrinsic layer.
+        [V].
     method : str
         one of two optional search methods: either ``'brentq'``, a reliable and
         bounded method or ``'newton'`` which is the default.
@@ -265,7 +291,7 @@ def bishop88_v_from_i(current, photocurrent, saturation_current,
     """
     # collect args
     args = (photocurrent, saturation_current, resistance_series,
-            resistance_shunt, nNsVth)
+            resistance_shunt, nNsVth, d2mutau, NsVbi)
     # first bound the search using voc
     voc_est = estimate_voc(photocurrent, saturation_current, nNsVth)
 
@@ -277,8 +303,9 @@ def fi(x, i, *a):
         # brentq only works with scalar inputs, so we need a set up function
         # and np.vectorize to repeatedly call the optimizer with the right
         # arguments for possible array input
-        def vd_from_brent(voc, i, iph, isat, rs, rsh, gamma):
-            return brentq(fi, 0.0, voc, args=(i, iph, isat, rs, rsh, gamma))
+        def vd_from_brent(voc, i, iph, isat, rs, rsh, gamma, d2mutau, NsVbi):
+            return brentq(fi, 0.0, voc,
+                          args=(i, iph, isat, rs, rsh, gamma, d2mutau, NsVbi))
 
         vd_from_brent_vectorized = np.vectorize(vd_from_brent)
         vd = vd_from_brent_vectorized(voc_est, current, *args)
@@ -295,7 +322,8 @@ def vd_from_brent(voc, i, iph, isat, rs, rsh, gamma):
 
 
 def bishop88_mpp(photocurrent, saturation_current, resistance_series,
-                 resistance_shunt, nNsVth, method='newton'):
+                 resistance_shunt, nNsVth, d2mutau=0, NsVbi=np.Inf,
+                 method='newton'):
     """
     Find max power point.
 
@@ -312,6 +340,17 @@ def bishop88_mpp(photocurrent, saturation_current, resistance_series,
     nNsVth : numeric
         product of diode ideality factor (n), number of series cells (Ns), and
         thermal voltage (Vth = k_b * T / q_e) in volts [V]
+    d2mutau : numeric, default 0
+        PVsyst parameter for cadmium-telluride (CdTe) and amorphous-silicon
+        (a-Si) modules that accounts for recombination current in the
+        intrinsic layer. The value is the ratio of intrinsic layer thickness
+        squared :math:`d^2` to the diffusion length of charge carriers
+        :math:`\\mu \\tau`. [V]
+    NsVbi : numeric, default np.inf
+        PVsyst parameter for cadmium-telluride (CdTe) and amorphous-silicon
+        (a-Si) modules that is the product of the PV module number of series
+        cells ``Ns`` and the builtin voltage ``Vbi`` of the intrinsic layer.
+        [V].
     method : str
         one of two optional search methods: either ``'brentq'``, a reliable and
         bounded method or ``'newton'`` which is the default.
@@ -324,7 +363,7 @@ def bishop88_mpp(photocurrent, saturation_current, resistance_series,
     """
     # collect args
     args = (photocurrent, saturation_current, resistance_series,
-            resistance_shunt, nNsVth)
+            resistance_shunt, nNsVth, d2mutau, NsVbi)
     # first bound the search using voc
     voc_est = estimate_voc(photocurrent, saturation_current, nNsVth)
 
@@ -334,8 +373,9 @@ def fmpp(x, *a):
     if method.lower() == 'brentq':
         # break out arguments for numpy.vectorize to handle broadcasting
         vec_fun = np.vectorize(
-            lambda voc, iph, isat, rs, rsh, gamma:
-                brentq(fmpp, 0.0, voc, args=(iph, isat, rs, rsh, gamma))
+            lambda voc, iph, isat, rs, rsh, gamma, d2mutau, NsVbi:
+                brentq(fmpp, 0.0, voc,
+                       args=(iph, isat, rs, rsh, gamma, d2mutau, NsVbi))
         )
         vd = vec_fun(voc_est, *args)
     elif method.lower() == 'newton':

pvlib/test/test_singlediode.py

@@ -4,7 +4,8 @@
 
 import numpy as np
 from pvlib import pvsystem
-from pvlib.singlediode import bishop88, estimate_voc, VOLTAGE_BUILTIN
+from pvlib.singlediode import (bishop88_mpp, estimate_voc, VOLTAGE_BUILTIN,
+                               bishop88, bishop88_i_from_v, bishop88_v_from_i)
 import pytest
 from conftest import requires_scipy
 
@@ -153,9 +154,13 @@ def get_pvsyst_fs_495():
         'temp_ref': 25, 'irrad_ref': 1000, 'I_L_ref': 1.5743233463848496
     }
 
+# DeSoto @(888[W/m**2], 55[degC]) = {Pmp: 72.71, Isc: 1.402, Voc: 75.42)
 
+
+@requires_scipy
 @pytest.mark.parametrize(
     'poa, temp_cell, expected, tol', [
+        # reference conditions
         (
             get_pvsyst_fs_495()['irrad_ref'],
             get_pvsyst_fs_495()['temp_ref'],
@@ -167,9 +172,21 @@ def get_pvsyst_fs_495():
             },
             (5e-4, 0.04)
         ),
-        (POA, TCELL, {'pmp': 76.26, 'isc': 1.387, 'voc': 79.29}, (1e-3, 1e-3))]
-)  # DeSoto @(888[W/m**2], 55[degC]) = {Pmp: 72.71, Isc: 1.402, Voc: 75.42)
-def test_pvsyst_recombination_loss(poa, temp_cell, expected, tol):
+        # other conditions
+        (
+            POA,
+            TCELL,
+            {
+                'pmp': 76.262,
+                'isc': 1.3868,
+                'voc': 79.292
+            },
+            (1e-4, 1e-4)
+        )
+    ]
+)
+@pytest.mark.parametrize('method', ['newton', 'brentq'])
+def test_pvsyst_recombination_loss(method, poa, temp_cell, expected, tol):
     """test PVSst recombination loss"""
     pvsyst_fs_495 = get_pvsyst_fs_495()
     # first evaluate PVSyst model with thin-film recombination loss current
@@ -199,9 +216,30 @@ def test_pvsyst_recombination_loss(poa, temp_cell, expected, tol):
     )
     # test max power
     assert np.isclose(max(pvsyst[2]), expected['pmp'], *tol)
+
     # test short circuit current
     isc_pvsyst = np.interp(0, pvsyst[1], pvsyst[0])
     assert np.isclose(isc_pvsyst, expected['isc'], *tol)
-    # test open circuit current
+
+    # test open circuit voltage
     voc_pvsyst = np.interp(0, pvsyst[0][::-1], pvsyst[1][::-1])
     assert np.isclose(voc_pvsyst, expected['voc'], *tol)
+
+    # repeat tests as above with specialized bishop88 functions
+    y = dict(d2mutau=pvsyst_fs_495['d2mutau'],
+             NsVbi=VOLTAGE_BUILTIN*pvsyst_fs_495['cells_in_series'])
+
+    mpp_88 = bishop88_mpp(*x, **y, method=method)
+    assert np.isclose(mpp_88[2], expected['pmp'], *tol)
+
+    isc_88 = bishop88_i_from_v(0, *x, **y, method=method)
+    assert np.isclose(isc_88, expected['isc'], *tol)
+
+    voc_88 = bishop88_v_from_i(0, *x, **y, method=method)
+    assert np.isclose(voc_88, expected['voc'], *tol)
+
+    ioc_88 = bishop88_i_from_v(voc_88, *x, **y, method=method)
+    assert np.isclose(ioc_88, 0.0, *tol)
+
+    vsc_88 = bishop88_v_from_i(isc_88, *x, **y, method=method)
+    assert np.isclose(vsc_88, 0.0, *tol)