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

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

@@ -18,7 +18,16 @@ Bug fixes
 
 Enhancements
 ~~~~~~~~~~~~
-
+* Added ``method='chandrupatla'`` (faster than ``brentq`` and slower than ``newton``,
+  but convergence is guaranteed) to
+  :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
 ~~~~~~~~~~~~~
@@ -44,3 +53,4 @@ Maintenance
 Contributors
 ~~~~~~~~~~~~
 * Elijah Passmore (:ghuser:`eljpsm`)
+* Kevin Anderson (:ghuser:`kandersolar`)

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 1st quadrant
+        only.
+
+        .. 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 1st quadrant
+        only.
+
+        .. note::
+           ``'chandrupatla'`` requires scipy 1.15 or greater.
+
 
     Returns
     -------

pvlib/singlediode.py

@@ -109,13 +109,13 @@
         (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 @@
     # 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 @@
         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 @@
     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 @@
     .. [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 @@
         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 @@
     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 @@
     .. [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 @@
         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 @@
     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 @@
     .. [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 @@
         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

@@ -4,6 +4,7 @@
 
 import pandas as pd
 import os
+import sys
 from packaging.version import Version
 import pytest
 from functools import wraps
@@ -194,6 +195,18 @@ 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.  That is only available for python 3.10 and above, so
+# we need to skip those tests on python 3.9.
+# TODO remove this when we drop support for python 3.9.
+chandrupatla_available = sys.version_info >= (3, 10)
+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)