Add scipy.optimize.elementwise.find_root (method='chandrupatla') to bishop88 functions

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

@@ -18,6 +18,16 @@ Bug fixes
 
 Enhancements
 ~~~~~~~~~~~~
+* Add ``method='chandrupatla'`` (faster than ``brentq`` and slower than ``newton``,
+ but convergence is guaranteed) as an option for
+ :py:func:`pvlib.pvsystem.singlediode`,
+ :py:func:`~pvlib.pvsystem.i_from_v`,
+ :py:func:`~pvlib.pvsystem.v_from_i`,
+ :py:func:`~pvlib.pvsystem.max_power_point`,
+ :py:func:`~pvlib.singlediode.bishop88_mpp`,
+ :py:func:`~pvlib.singlediode.bishop88_v_from_i`, and
+ :py:func:`~pvlib.singlediode.bishop88_i_from_v`. (:issue:`2497`, :pull:`2498`)
+
 
 
 Documentation

pvlib/pvsystem.py

@@ -2498,7 +2498,11 @@ def singlediode(photocurrent, saturation_current, resistance_series,
 
  method : str, default 'lambertw'
  Determines the method used to calculate points on the IV curve. The
- options are ``'lambertw'``, ``'newton'``, or ``'brentq'``.
+ options are ``'lambertw'``, ``'newton'``, ``'brentq'``, or
+ ``'chandrupatla'``.
+
+ .. note::
+ ``'chandrupatla'`` requires scipy 1.15 or greater.
 
  Returns
  -------
@@ -2630,7 +2634,11 @@ def max_power_point(photocurrent, saturation_current, resistance_series,
  cells ``Ns`` and the builtin voltage ``Vbi`` of the intrinsic layer.
  [V].
  method : str
- either ``'newton'`` or ``'brentq'``
+ either ``'newton'``, ``'brentq'``, or ``'chandrupatla'``.
+
+ .. note::
+ ``'chandrupatla'`` requires scipy 1.15 or greater.
+
 
  Returns
  -------
@@ -2713,8 +2721,13 @@ def v_from_i(current, photocurrent, saturation_current, resistance_series,
  0 < nNsVth
 
  method : str
- Method to use: ``'lambertw'``, ``'newton'``, or ``'brentq'``. *Note*:
- ``'brentq'`` is limited to 1st quadrant only.
+ Method to use: ``'lambertw'``, ``'newton'``, ``'brentq'``, or
+ ``'chandrupatla'``. *Note*: ``'brentq'`` is limited to
+ non-negative current.
+
+ .. note::
+ ``'chandrupatla'`` requires scipy 1.15 or greater.
+
 
  Returns
  -------
@@ -2795,8 +2808,13 @@ def i_from_v(voltage, photocurrent, saturation_current, resistance_series,
  0 < nNsVth
 
  method : str
- Method to use: ``'lambertw'``, ``'newton'``, or ``'brentq'``. *Note*:
- ``'brentq'`` is limited to 1st quadrant only.
+ Method to use: ``'lambertw'``, ``'newton'``, ``'brentq'``, or
+ ``'chandrupatla'``. *Note*: ``'brentq'`` is limited to
+ non-negative current.
+
+ .. note::
+ ``'chandrupatla'`` requires scipy 1.15 or greater.
+
 
  Returns
  -------

pvlib/singlediode.py

@@ -109,13 +109,13 @@ def bishop88(diode_voltage, photocurrent, saturation_current,
  (a-Si) modules that is the product of the PV module number of series
  cells :math:`N_{s}` and the builtin voltage :math:`V_{bi}` of the
  intrinsic layer. [V].
- breakdown_factor : float, default 0
+ 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 : float, default -5.5
+ breakdown_voltage : numeric, default -5.5
  reverse breakdown voltage of the photovoltaic junction :math:`V_{br}`
  [V]
- breakdown_exp : float, default 3.28
+ 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
@@ -162,12 +162,11 @@ 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
- 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.
+
+ 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
+
  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
@@ -177,18 +176,14 @@ def bishop88(diode_voltage, photocurrent, saturation_current,
  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
- 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.
+ 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
  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
@@ -247,12 +242,19 @@ def bishop88_i_from_v(voltage, photocurrent, saturation_current,
  breakdown_exp : float, 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.
+ Either ``'newton'``, ``'brentq'``, or ``'chandrupatla'``.
+ ''method'' must be ``'newton'`` if ``breakdown_factor`` is not 0.
+
+ .. note::
+ ``'chandrupatla'`` requires scipy 1.15 or greater.
+
  method_kwargs : dict, optional
- Keyword arguments passed to root finder method. See
- :py:func:`scipy:scipy.optimize.brentq` and
- :py:func:`scipy:scipy.optimize.newton` parameters.
+ Keyword arguments passed to the root finder. For options, see:
+
+ * ``method='brentq'``: :py:func:`scipy:scipy.optimize.brentq`
+ * ``method='newton'``: :py:func:`scipy:scipy.optimize.newton`
+ * ``method='chandrupatla'``: :py:func:`scipy:scipy.optimize.elementwise.find_root`
+
  ``'full_output': True`` is allowed, and ``optimizer_output`` would be
  returned. See examples section.
 
@@ -291,7 +293,7 @@ def bishop88_i_from_v(voltage, photocurrent, saturation_current,
  .. [1] "Computer simulation of the effects of electrical mismatches in
  photovoltaic cell interconnection circuits" JW Bishop, Solar Cell (1988)
  :doi:`10.1016/0379-6787(88)90059-2`
- """
+ """ # noqa: E501
  # collect args
  args = (photocurrent, saturation_current,
  resistance_series, resistance_shunt, nNsVth, d2mutau, NsVbi,
@@ -333,6 +335,30 @@ def vd_from_brent(voc, v, iph, isat, rs, rsh, gamma, d2mutau, NsVbi,
  vd = newton(func=lambda x, *a: fv(x, voltage, *a), x0=x0,
  fprime=lambda x, *a: bishop88(x, *a, gradients=True)[4],
  args=args, **method_kwargs)
+ elif method == 'chandrupatla':
+ try:
+ from scipy.optimize.elementwise import find_root
+ except ModuleNotFoundError as e:
+ # TODO remove this when our minimum scipy version is >=1.15
+ msg = (
+ "method='chandrupatla' requires scipy v1.15 or greater "
+ "(available for Python 3.10+). "
+ "Select another method, or update your version of scipy."
+ )
+ raise ImportError(msg) from e
+
+ voc_est = estimate_voc(photocurrent, saturation_current, nNsVth)
+ shape = _shape_of_max_size(voltage, voc_est)
+ vlo = np.zeros(shape)
+ vhi = np.full(shape, voc_est)
+ bounds = (vlo, vhi)
+ kwargs_trimmed = method_kwargs.copy()
+ kwargs_trimmed.pop("full_output", None) # not valid for find_root
+
+ result = find_root(fv, bounds, args=(voltage, *args), **kwargs_trimmed)
+ vd = result.x
+ if method_kwargs.get('full_output'):
+ vd = (vd, result) # mimic the other methods
  else:
  raise NotImplementedError("Method '%s' isn't implemented" % method)
 
@@ -388,12 +414,19 @@ def bishop88_v_from_i(current, photocurrent, saturation_current,
  breakdown_exp : float, 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.
+ Either ``'newton'``, ``'brentq'``, or ``'chandrupatla'``.
+ ''method'' must be ``'newton'`` if ``breakdown_factor`` is not 0.
+
+ .. note::
+ ``'chandrupatla'`` requires scipy 1.15 or greater.
+
  method_kwargs : dict, optional
- Keyword arguments passed to root finder method. See
- :py:func:`scipy:scipy.optimize.brentq` and
- :py:func:`scipy:scipy.optimize.newton` parameters.
+ Keyword arguments passed to the root finder. For options, see:
+
+ * ``method='brentq'``: :py:func:`scipy:scipy.optimize.brentq`
+ * ``method='newton'``: :py:func:`scipy:scipy.optimize.newton`
+ * ``method='chandrupatla'``: :py:func:`scipy:scipy.optimize.elementwise.find_root`
+
  ``'full_output': True`` is allowed, and ``optimizer_output`` would be
  returned. See examples section.
 
@@ -432,7 +465,7 @@ def bishop88_v_from_i(current, photocurrent, saturation_current,
  .. [1] "Computer simulation of the effects of electrical mismatches in
  photovoltaic cell interconnection circuits" JW Bishop, Solar Cell (1988)
  :doi:`10.1016/0379-6787(88)90059-2`
- """
+ """ # noqa: E501
  # collect args
  args = (photocurrent, saturation_current,
  resistance_series, resistance_shunt, nNsVth, d2mutau, NsVbi,
@@ -474,6 +507,29 @@ def vd_from_brent(voc, i, iph, isat, rs, rsh, gamma, d2mutau, NsVbi,
  vd = newton(func=lambda x, *a: fi(x, current, *a), x0=x0,
  fprime=lambda x, *a: bishop88(x, *a, gradients=True)[3],
  args=args, **method_kwargs)
+ elif method == 'chandrupatla':
+ try:
+ from scipy.optimize.elementwise import find_root
+ except ModuleNotFoundError as e:
+ # TODO remove this when our minimum scipy version is >=1.15
+ msg = (
+ "method='chandrupatla' requires scipy v1.15 or greater "
+ "(available for Python 3.10+). "
+ "Select another method, or update your version of scipy."
+ )
+ raise ImportError(msg) from e
+
+ shape = _shape_of_max_size(current, voc_est)
+ vlo = np.zeros(shape)
+ vhi = np.full(shape, voc_est)
+ bounds = (vlo, vhi)
+ kwargs_trimmed = method_kwargs.copy()
+ kwargs_trimmed.pop("full_output", None) # not valid for find_root
+
+ result = find_root(fi, bounds, args=(current, *args), **kwargs_trimmed)
+ vd = result.x
+ if method_kwargs.get('full_output'):
+ vd = (vd, result) # mimic the other methods
  else:
  raise NotImplementedError("Method '%s' isn't implemented" % method)
 
@@ -526,12 +582,19 @@ def bishop88_mpp(photocurrent, saturation_current, resistance_series,
  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.
+ Either ``'newton'``, ``'brentq'``, or ``'chandrupatla'``.
+ ''method'' must be ``'newton'`` if ``breakdown_factor`` is not 0.
+
+ .. note::
+ ``'chandrupatla'`` requires scipy 1.15 or greater.
+
  method_kwargs : dict, optional
- Keyword arguments passed to root finder method. See
- :py:func:`scipy:scipy.optimize.brentq` and
- :py:func:`scipy:scipy.optimize.newton` parameters.
+ Keyword arguments passed to the root finder. For options, see:
+
+ * ``method='brentq'``: :py:func:`scipy:scipy.optimize.brentq`
+ * ``method='newton'``: :py:func:`scipy:scipy.optimize.newton`
+ * ``method='chandrupatla'``: :py:func:`scipy:scipy.optimize.elementwise.find_root`
+
  ``'full_output': True`` is allowed, and ``optimizer_output`` would be
  returned. See examples section.
 
@@ -571,7 +634,7 @@ def bishop88_mpp(photocurrent, saturation_current, resistance_series,
  .. [1] "Computer simulation of the effects of electrical mismatches in
  photovoltaic cell interconnection circuits" JW Bishop, Solar Cell (1988)
  :doi:`10.1016/0379-6787(88)90059-2`
- """
+ """ # noqa: E501
  # collect args
  args = (photocurrent, saturation_current,
  resistance_series, resistance_shunt, nNsVth, d2mutau, NsVbi,
@@ -611,6 +674,31 @@ def fmpp(x, *a):
  vd = newton(func=fmpp, x0=x0,
  fprime=lambda x, *a: bishop88(x, *a, gradients=True)[7],
  args=args, **method_kwargs)
+ elif method == 'chandrupatla':
+ try:
+ from scipy.optimize.elementwise import find_root
+ except ModuleNotFoundError as e:
+ # TODO remove this when our minimum scipy version is >=1.15
+ msg = (
+ "method='chandrupatla' requires scipy v1.15 or greater "
+ "(available for Python 3.10+). "
+ "Select another method, or update your version of scipy."
+ )
+ raise ImportError(msg) from e
+
+ vlo = np.zeros_like(photocurrent)
+ vhi = np.full_like(photocurrent, voc_est)
+ kwargs_trimmed = method_kwargs.copy()
+ kwargs_trimmed.pop("full_output", None) # not valid for find_root
+
+ result = find_root(fmpp,
+ (vlo, vhi),
+ args=args,
+ **kwargs_trimmed)
+ vd = result.x
+ if method_kwargs.get('full_output'):
+ vd = (vd, result) # mimic the other methods
+
  else:
  raise NotImplementedError("Method '%s' isn't implemented" % method)
 

tests/conftest.py

@@ -3,6 +3,7 @@
 import warnings
 
 import pandas as pd
+import scipy
 import os
 from packaging.version import Version
 import pytest
@@ -194,6 +195,17 @@ def has_spa_c():
  reason="requires pandas>=2.0.0")
 
 
+# single-diode equation functions have method=='chandrupatla', which relies
+# on scipy.optimize.elementwise.find_root, which is only available in
+# scipy>=1.15.
+# TODO remove this when our minimum scipy is >=1.15
+chandrupatla_available = Version(scipy.__version__) >= Version("1.15.0")
+chandrupatla = pytest.param(
+ "chandrupatla", marks=pytest.mark.skipif(not chandrupatla_available,
+ reason="needs scipy 1.15")
+)
+
+
 @pytest.fixture()
 def golden():
  return Location(39.742476, -105.1786, 'America/Denver', 1830.14)