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

kandersolar · echedey-ls · cwhanse · web-flow · commit 74ded020440f · 2025-10-03T09:05:24.000-04:00
* add method='chandrupatla' for bishop88 functions

* update tests

* whatsnew

* doc tweaks

* test tweak

* try out skipping chandrupatla on py3.9

* fix py3.9 skips

* Apply suggestions from code review

Co-authored-by: Echedey Luis &lt;80125792+echedey-ls@users.noreply.github.com&gt;

* more edits from review

* docs fixes

* Apply suggestions from code review

Co-authored-by: Cliff Hansen &lt;cwhanse@sandia.gov&gt;

* move whatsnew entry to v0.13.2

* fix test and lint issues

---------

Co-authored-by: Echedey Luis &lt;80125792+echedey-ls@users.noreply.github.com&gt;
Co-authored-by: Cliff Hansen &lt;cwhanse@sandia.gov&gt;
diff --git a/docs/sphinx/source/whatsnew/v0.13.2.rst b/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
diff --git a/pvlib/pvsystem.py b/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
     -------
diff --git a/pvlib/singlediode.py b/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)
 
diff --git a/tests/conftest.py b/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)
diff --git a/tests/test_pvsystem.py b/tests/test_pvsystem.py
diff --git a/tests/test_singlediode.py b/tests/test_singlediode.py
