add capability for reverse bias (#948)

* add capability for reverse bias

* lint, correct test_pvsyst_breakdown

* default value for breakdown_factor

* handle recombination in breakdown test, lint

* correct parameter parsing

* whatsnew

* Update pvlib/tests/test_singlediode.py

Co-Authored-By: Cameron Stark <CameronTStark@users.noreply.github.com>

* Update pvlib/singlediode.py

Co-Authored-By: Cameron Stark <CameronTStark@users.noreply.github.com>

* Update pvlib/singlediode.py

Co-Authored-By: Cameron Stark <CameronTStark@users.noreply.github.com>

* fix formatting

* improve docstring

* lint

Co-authored-by: Cameron Stark <CameronTStark@users.noreply.github.com>

Author: cwhanse

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

@@ -35,6 +35,10 @@ Enhancements
 * Add new module :py:mod:`pvlib.snow` to contain models related to snow coverage and effects on a PV system. (:pull:`764`)
 * Add snow coverage model :py:func:`pvlib.snow.coverage_nrel` and function to identify when modules are fully covered by snow :py:func:`pvlib.snow.fully_covered_nrel`. (:issue:`577`)
 * Add function :py:func:`pvlib.snow.dc_loss_nrel` for effect of snow coverage on DC output. (:pull:`764`)
+* Add capability to calculate current at reverse bias using an avalanche
+  breakdown model, affects :py:func:`pvlib.singlediode.bishop88`,
+  :py:func:`pvlib.singlediode.bishop88_i_from_v`, :py:func:`pvlib.singlediode.bishop88_v_from_i`,
+  :py:func:`pvlib.singlediode.bishop88_mpp`. (:pull:`948`)
 
 Bug fixes
 ~~~~~~~~~

pvlib/singlediode.py

@@ -69,42 +69,65 @@ def estimate_voc(photocurrent, saturation_current, nNsVth):
 
 def bishop88(diode_voltage, photocurrent, saturation_current,
              resistance_series, resistance_shunt, nNsVth, d2mutau=0,
-             NsVbi=np.Inf, gradients=False):
-    """
+             NsVbi=np.Inf, breakdown_factor=0., breakdown_voltage=-5.5,
+             breakdown_exp=3.28, gradients=False):
+    r"""
     Explicit calculation of points on the IV curve described by the single
     diode equation.  Values are calculated as described in [1]_.
 
+    The single diode equation with recombination current and reverse bias
+    breakdown is
+
+    .. math::
+
+        I = I_{L} - I_{0} (\exp \frac{V_{d}}{nNsVth} - 1)
+        - \frac{V_{d}}{R_{sh}}
+        - \frac{I_{L} \frac{d^{2}}{\mu \tau}{N_{s} V_{bi} - V_{d}}
+        - a \frac{V_{d}{R_{sh}} (1 - \frac{V_{d}}{V_{br}})^-m
+
+    The input `diode_voltage` must be :math:`V + I R_{s}`.
+
+
     .. warning::
+       * Usage of ``d2mutau`` is required with PVSyst
+         coefficients for cadmium-telluride (CdTe) and amorphous-silicon
+         (a:Si) PV modules only.
        * Do not use ``d2mutau`` with CEC coefficients.
-       * Usage of ``d2mutau`` with PVSyst coefficients is required for cadmium-
-         telluride (CdTe) and amorphous-silicon (a:Si) PV modules only.
 
     Parameters
     ----------
     diode_voltage : numeric
         diode voltages [V]
     photocurrent : numeric
-        photo-generated current [A]
+        photo-generated current :math:`I_{L}` [A]
     saturation_current : numeric
-        diode reverse saturation current [A]
+        diode reverse saturation current :math:`I_{0}` [A]
     resistance_series : numeric
-        series resistance [ohms]
+        series resistance :math:`R_{s}` [ohms]
     resistance_shunt: numeric
-        shunt resistance [ohms]
+        shunt resistance :math:`R_{sh}` [ohms]
     nNsVth : numeric
-        product of thermal voltage ``Vth`` [V], diode ideality factor ``n``,
-        and number of series cells ``Ns``
+        product of thermal voltage :math:`V_{th}` [V], diode ideality factor
+        ``n``, and number of series cells :math:`N_{s}`
     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]
+        :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].
+        cells :math:`N_{s}` and the builtin voltage :math:`V_{bi}` of the
+        intrinsic layer. [V].
+    breakdown_factor : numeric, default 0
+        fraction of ohmic current involved in avalanche breakdown :math:`a`.
+        Default of 0 excludes the reverse bias term from the model. [unitless]
+    breakdown_voltage : numeric, default -5.5
+        reverse breakdown voltage of the photovoltaic junction :math:`V_{br}`
+        [V]
+    breakdown_exp : numeric, default 3.28
+        avalanche breakdown exponent :math:`m` [unitless]
     gradients : bool
         False returns only I, V, and P. True also returns gradients
 
@@ -150,21 +173,39 @@ def bishop88(diode_voltage, photocurrent, saturation_current,
     # calculate temporary values to simplify calculations
     v_star = diode_voltage / nNsVth  # non-dimensional diode voltage
     g_sh = 1.0 / resistance_shunt  # conductance
-    i = (photocurrent - saturation_current * np.expm1(v_star)
-         - diode_voltage * g_sh - i_recomb)
+    if breakdown_factor > 0:  # reverse bias is considered
+        brk_term = 1 - diode_voltage / breakdown_voltage
+        brk_pwr = np.power(brk_term, -breakdown_exp)
+        i_breakdown = breakdown_factor * diode_voltage * g_sh * brk_pwr
+    else:
+        i_breakdown = 0.
+    i = (photocurrent - saturation_current * np.expm1(v_star)  # noqa: W503
+         - diode_voltage * g_sh - i_recomb - i_breakdown)   # noqa: W503
     v = diode_voltage - i * resistance_series
     retval = (i, v, i*v)
     if gradients:
         # calculate recombination loss current gradients where d2mutau > 0
         grad_i_recomb = np.where(is_recomb, i_recomb / v_recomb, 0)
         grad_2i_recomb = np.where(is_recomb, 2 * grad_i_recomb / v_recomb, 0)
         g_diode = saturation_current * np.exp(v_star) / nNsVth  # conductance
-        grad_i = -g_diode - g_sh - grad_i_recomb  # di/dvd
+        if breakdown_factor > 0:  # reverse bias is considered
+            brk_pwr_1 = np.power(brk_term, -breakdown_exp - 1)
+            brk_pwr_2 = np.power(brk_term, -breakdown_exp - 2)
+            brk_fctr = breakdown_factor * g_sh
+            grad_i_brk = brk_fctr * (brk_pwr + diode_voltage *
+                                     -breakdown_exp * brk_pwr_1)
+            grad2i_brk = (brk_fctr * -breakdown_exp        # noqa: W503
+                          * (2 * brk_pwr_1 + diode_voltage   # noqa: W503
+                             * (-breakdown_exp - 1) * brk_pwr_2))  # noqa: W503
+        else:
+            grad_i_brk = 0.
+            grad2i_brk = 0.
+        grad_i = -g_diode - g_sh - grad_i_recomb - grad_i_brk  # di/dvd
         grad_v = 1.0 - grad_i * resistance_series  # dv/dvd
         # dp/dv = d(iv)/dv = v * di/dv + i
         grad = grad_i / grad_v  # di/dv
         grad_p = v * grad + i  # dp/dv
-        grad2i = -g_diode / nNsVth - grad_2i_recomb  # d2i/dvd
+        grad2i = -g_diode / nNsVth - grad_2i_recomb - grad2i_brk  # d2i/dvd
         grad2v = -grad2i * resistance_series  # d2v/dvd
         grad2p = (
             grad_v * grad + v * (grad2i/grad_v - grad_i*grad2v/grad_v**2)
@@ -176,7 +217,9 @@ def bishop88(diode_voltage, photocurrent, saturation_current,
 
 def bishop88_i_from_v(voltage, photocurrent, saturation_current,
                       resistance_series, resistance_shunt, nNsVth,
-                      d2mutau=0, NsVbi=np.Inf, method='newton'):
+                      d2mutau=0, NsVbi=np.Inf, breakdown_factor=0.,
+                      breakdown_voltage=-5.5, breakdown_exp=3.28,
+                      method='newton'):
     """
     Find current given any voltage.
 
@@ -185,13 +228,13 @@ def bishop88_i_from_v(voltage, photocurrent, saturation_current,
     voltage : numeric
         voltage (V) in volts [V]
     photocurrent : numeric
-        photogenerated current (Iph or IL) in amperes [A]
+        photogenerated current (Iph or IL) [A]
     saturation_current : numeric
-        diode dark or saturation current (Io or Isat) in amperes [A]
+        diode dark or saturation current (Io or Isat) [A]
     resistance_series : numeric
-        series resistance (Rs) in ohms
+        series resistance (Rs) in [Ohm]
     resistance_shunt : numeric
-        shunt resistance (Rsh) in ohms
+        shunt resistance (Rsh) [Ohm]
     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]
@@ -206,18 +249,27 @@ def bishop88_i_from_v(voltage, photocurrent, saturation_current,
         (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.
+    breakdown_factor : numeric, default 0
+        fraction of ohmic current involved in avalanche breakdown :math:`a`.
+        Default of 0 excludes the reverse bias term from the model. [unitless]
+    breakdown_voltage : numeric, default -5.5
+        reverse breakdown voltage of the photovoltaic junction :math:`V_{br}`
+        [V]
+    breakdown_exp : numeric, default 3.28
+        avalanche breakdown exponent :math:`m` [unitless]
+    method : str, default 'newton'
+       Either ``'newton'`` or ``'brentq'``. ''method'' must be ``'newton'``
+       if ``breakdown_factor`` is not 0.
 
     Returns
     -------
     current : numeric
-        current (I) at the specified voltage (V) in amperes [A]
+        current (I) at the specified voltage (V). [A]
     """
     # collect args
     args = (photocurrent, saturation_current, resistance_series,
-            resistance_shunt, nNsVth, d2mutau, NsVbi)
+            resistance_shunt, nNsVth, d2mutau, NsVbi,
+            breakdown_factor, breakdown_voltage, breakdown_exp)
 
     def fv(x, v, *a):
         # calculate voltage residual given diode voltage "x"
@@ -230,9 +282,12 @@ 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, d2mutau, NsVbi):
+        def vd_from_brent(voc, v, iph, isat, rs, rsh, gamma, d2mutau, NsVbi,
+                          breakdown_factor, breakdown_voltage, breakdown_exp):
             return brentq(fv, 0.0, voc,
-                          args=(v, iph, isat, rs, rsh, gamma, d2mutau, NsVbi))
+                          args=(v, iph, isat, rs, rsh, gamma, d2mutau, NsVbi,
+                                breakdown_factor, breakdown_voltage,
+                                breakdown_exp))
 
         vd_from_brent_vectorized = np.vectorize(vd_from_brent)
         vd = vd_from_brent_vectorized(voc_est, voltage, *args)
@@ -250,7 +305,9 @@ def vd_from_brent(voc, v, iph, isat, rs, rsh, gamma, d2mutau, NsVbi):
 
 def bishop88_v_from_i(current, photocurrent, saturation_current,
                       resistance_series, resistance_shunt, nNsVth,
-                      d2mutau=0, NsVbi=np.Inf, method='newton'):
+                      d2mutau=0, NsVbi=np.Inf, breakdown_factor=0.,
+                      breakdown_voltage=-5.5, breakdown_exp=3.28,
+                      method='newton'):
     """
     Find voltage given any current.
 
@@ -259,13 +316,13 @@ def bishop88_v_from_i(current, photocurrent, saturation_current,
     current : numeric
         current (I) in amperes [A]
     photocurrent : numeric
-        photogenerated current (Iph or IL) in amperes [A]
+        photogenerated current (Iph or IL) [A]
     saturation_current : numeric
-        diode dark or saturation current (Io or Isat) in amperes [A]
+        diode dark or saturation current (Io or Isat) [A]
     resistance_series : numeric
-        series resistance (Rs) in ohms
+        series resistance (Rs) in [Ohm]
     resistance_shunt : numeric
-        shunt resistance (Rsh) in ohms
+        shunt resistance (Rsh) [Ohm]
     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]
@@ -280,9 +337,17 @@ def bishop88_v_from_i(current, photocurrent, saturation_current,
         (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.
+    breakdown_factor : numeric, default 0
+        fraction of ohmic current involved in avalanche breakdown :math:`a`.
+        Default of 0 excludes the reverse bias term from the model. [unitless]
+    breakdown_voltage : numeric, default -5.5
+        reverse breakdown voltage of the photovoltaic junction :math:`V_{br}`
+        [V]
+    breakdown_exp : numeric, default 3.28
+        avalanche breakdown exponent :math:`m` [unitless]
+    method : str, default 'newton'
+       Either ``'newton'`` or ``'brentq'``. ''method'' must be ``'newton'``
+       if ``breakdown_factor`` is not 0.
 
     Returns
     -------
@@ -291,7 +356,8 @@ def bishop88_v_from_i(current, photocurrent, saturation_current,
     """
     # collect args
     args = (photocurrent, saturation_current, resistance_series,
-            resistance_shunt, nNsVth, d2mutau, NsVbi)
+            resistance_shunt, nNsVth, d2mutau, NsVbi, breakdown_factor,
+            breakdown_voltage, breakdown_exp)
     # first bound the search using voc
     voc_est = estimate_voc(photocurrent, saturation_current, nNsVth)
 
@@ -303,9 +369,12 @@ 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, d2mutau, NsVbi):
+        def vd_from_brent(voc, i, iph, isat, rs, rsh, gamma, d2mutau, NsVbi,
+                          breakdown_factor, breakdown_voltage, breakdown_exp):
             return brentq(fi, 0.0, voc,
-                          args=(i, iph, isat, rs, rsh, gamma, d2mutau, NsVbi))
+                          args=(i, iph, isat, rs, rsh, gamma, d2mutau, NsVbi,
+                                breakdown_factor, breakdown_voltage,
+                                breakdown_exp))
 
         vd_from_brent_vectorized = np.vectorize(vd_from_brent)
         vd = vd_from_brent_vectorized(voc_est, current, *args)
@@ -323,20 +392,21 @@ def vd_from_brent(voc, i, iph, isat, rs, rsh, gamma, d2mutau, NsVbi):
 
 def bishop88_mpp(photocurrent, saturation_current, resistance_series,
                  resistance_shunt, nNsVth, d2mutau=0, NsVbi=np.Inf,
-                 method='newton'):
+                 breakdown_factor=0., breakdown_voltage=-5.5,
+                 breakdown_exp=3.28, method='newton'):
     """
     Find max power point.
 
     Parameters
     ----------
     photocurrent : numeric
-        photogenerated current (Iph or IL) in amperes [A]
+        photogenerated current (Iph or IL) [A]
     saturation_current : numeric
-        diode dark or saturation current (Io or Isat) in amperes [A]
+        diode dark or saturation current (Io or Isat) [A]
     resistance_series : numeric
-        series resistance (Rs) in ohms
+        series resistance (Rs) in [Ohm]
     resistance_shunt : numeric
-        shunt resistance (Rsh) in ohms
+        shunt resistance (Rsh) [Ohm]
     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]
@@ -351,9 +421,17 @@ def bishop88_mpp(photocurrent, saturation_current, resistance_series,
         (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.
+    breakdown_factor : numeric, default 0
+        fraction of ohmic current involved in avalanche breakdown :math:`a`.
+        Default of 0 excludes the reverse bias term from the model. [unitless]
+    breakdown_voltage : numeric, default -5.5
+        reverse breakdown voltage of the photovoltaic junction :math:`V_{br}`
+        [V]
+    breakdown_exp : numeric, default 3.28
+        avalanche breakdown exponent :math:`m` [unitless]
+    method : str, default 'newton'
+       Either ``'newton'`` or ``'brentq'``. ''method'' must be ``'newton'``
+       if ``breakdown_factor`` is not 0.
 
     Returns
     -------
@@ -363,7 +441,8 @@ def bishop88_mpp(photocurrent, saturation_current, resistance_series,
     """
     # collect args
     args = (photocurrent, saturation_current, resistance_series,
-            resistance_shunt, nNsVth, d2mutau, NsVbi)
+            resistance_shunt, nNsVth, d2mutau, NsVbi, breakdown_factor,
+            breakdown_voltage, breakdown_exp)
     # first bound the search using voc
     voc_est = estimate_voc(photocurrent, saturation_current, nNsVth)
 
@@ -373,9 +452,10 @@ 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, d2mutau, NsVbi:
-                brentq(fmpp, 0.0, voc,
-                       args=(iph, isat, rs, rsh, gamma, d2mutau, NsVbi))
+            lambda voc, iph, isat, rs, rsh, gamma, d2mutau, NsVbi, vbr_a, vbr,
+            vbr_exp: brentq(fmpp, 0.0, voc,
+                            args=(iph, isat, rs, rsh, gamma, d2mutau, NsVbi,
+                                  vbr_a, vbr, vbr_exp))
         )
         vd = vec_fun(voc_est, *args)
     elif method.lower() == 'newton':