add capability for reverse bias

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

@@ -25,6 +25,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) -
+ diode_voltage * g_sh - i_recomb - i_breakdown)
  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 * (
+ 2 * brk_pwr_1 + diode_voltage * (-breakdown_exp - 1) *
+ brk_pwr_2)
+ 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'
+ must be either ``'newton'`` or ``'brentq'``. If ``'breakdown_factor'``
+ is not 0, must be ``'newton'``.
 
  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'
+ must be either ``'newton'`` or ``'brentq'``. If ``'breakdown_factor'``
+ is not 0, must be ``'newton'``.
 
  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'
+ must be either ``'newton'`` or ``'brentq'``. If ``'breakdown_factor'``
+ is not 0, must be ``'newton'``.
 
  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':