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

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

@@ -18,7 +18,14 @@ Bug fixes
 
 Enhancements
 ~~~~~~~~~~~~
-
+* Added ``method='chandrupatla'`` 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`,
+  :py:func:`~pvlib.singlediode.bishop88_v_from_i`, and
+  :py:func:`~pvlib.singlediode.bishop88_i_from_v`. (:issue:`2497`, :pull:`2498`)
+
 
 Documentation
 ~~~~~~~~~~~~~
@@ -44,3 +51,4 @@ Maintenance
 Contributors
 ~~~~~~~~~~~~
 * Elijah Passmore (:ghuser:`eljpsm`)
+* Kevin Anderson (:ghuser:`kandersolar`)

pvlib/pvsystem.py

@@ -2498,7 +2498,8 @@ 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'`` (requires scipy 1.15 or greater).
 
     Returns
     -------
@@ -2630,7 +2631,8 @@ 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'`` (requires
+        scipy 1.15 or greater)
 
     Returns
     -------
@@ -2713,7 +2715,8 @@ def v_from_i(current, photocurrent, saturation_current, resistance_series,
         0 < nNsVth
 
     method : str
-        Method to use: ``'lambertw'``, ``'newton'``, or ``'brentq'``. *Note*:
+        Method to use: ``'lambertw'``, ``'newton'``, ``'brentq'``, or
+        ``'chandrupatla'`` (requires scipy 1.15 or greater). *Note*:
         ``'brentq'`` is limited to 1st quadrant only.
 
     Returns
@@ -2795,7 +2798,8 @@ def i_from_v(voltage, photocurrent, saturation_current, resistance_series,
         0 < nNsVth
 
     method : str
-        Method to use: ``'lambertw'``, ``'newton'``, or ``'brentq'``. *Note*:
+        Method to use: ``'lambertw'``, ``'newton'``, ``'brentq'``, or
+        ``'chandrupatla'`` (requires scipy 1.15 or greater). *Note*:
         ``'brentq'`` is limited to 1st quadrant only.
 
     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,14 @@
     breakdown_exp : float, default 3.28
         avalanche breakdown exponent :math:`m` [unitless]
     method : str, default 'newton'
-       Either ``'newton'`` or ``'brentq'``. ''method'' must be ``'newton'``
+       Either ``'newton'``, ``'brentq'``, or ``'chandrupatla'`` (requires
+       scipy 1.15 or greater). ''method'' must be ``'newton'``
        if ``breakdown_factor`` is not 0.
     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.
+        :py:func:`scipy:scipy.optimize.brentq`,
+        :py:func:`scipy:scipy.optimize.newton`, and
+        :py:func:`scipy:scipy.optimize.elementwise.find_root` for parameters.
         ``'full_output': True`` is allowed, and ``optimizer_output`` would be
         returned. See examples section.
 
@@ -333,6 +330,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. "
+                "Select another method, or update your version of scipy. "
+                f"({str(e)})"
+            )
+            raise ImportError(msg)
+
+        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,7 +409,8 @@
     breakdown_exp : float, default 3.28
         avalanche breakdown exponent :math:`m` [unitless]
     method : str, default 'newton'
-       Either ``'newton'`` or ``'brentq'``. ''method'' must be ``'newton'``
+       Either ``'newton'``, ``'brentq'``, or ``'chandrupatla'`` (requires
+       scipy 1.15 or greater). ''method'' must be ``'newton'``
        if ``breakdown_factor`` is not 0.
     method_kwargs : dict, optional
         Keyword arguments passed to root finder method. See
@@ -474,6 +496,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. "
+                "Select another method, or update your version of scipy. "
+                f"({str(e)})"
+            )
+            raise ImportError(msg)
+
+        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,7 +571,8 @@
     breakdown_exp : numeric, default 3.28
         avalanche breakdown exponent :math:`m` [unitless]
     method : str, default 'newton'
-       Either ``'newton'`` or ``'brentq'``. ''method'' must be ``'newton'``
+       Either ``'newton'``, ``'brentq'``, or ``'chandrupatla'`` (requires
+       scipy 1.15 or greater). ''method'' must be ``'newton'``
        if ``breakdown_factor`` is not 0.
     method_kwargs : dict, optional
         Keyword arguments passed to root finder method. See
@@ -611,6 +657,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. "
+                "Select another method, or update your version of scipy. "
+                f"({str(e)})"
+            )
+            raise ImportError(msg)
+
+        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/test_pvsystem.py

@@ -1371,7 +1371,8 @@ def fixture_i_from_v(request):
 
 
 @pytest.mark.parametrize(
-    'method, atol', [('lambertw', 1e-11), ('brentq', 1e-11), ('newton', 1e-11)]
+    'method, atol', [('lambertw', 1e-11), ('brentq', 1e-11), ('newton', 1e-11),
+                     ('chandrupatla', 1e-11)]
 )
 def test_i_from_v(fixture_i_from_v, method, atol):
     # Solution set loaded from fixture
@@ -1406,6 +1407,9 @@ def test_i_from_v_size():
     with pytest.raises(ValueError):
         pvsystem.i_from_v([7.5] * 3, 7., 6e-7, [0.1] * 2, 20, 0.5,
                           method='brentq')
+    with pytest.raises(ValueError):
+        pvsystem.i_from_v([7.5] * 3, 7., 6e-7, [0.1] * 2, 20, 0.5,
+                          method='chandrupatla')
     with pytest.raises(ValueError):
         pvsystem.i_from_v([7.5] * 3, np.array([7., 7.]), 6e-7, 0.1, 20, 0.5,
                           method='newton')
@@ -1417,6 +1421,9 @@ def test_v_from_i_size():
     with pytest.raises(ValueError):
         pvsystem.v_from_i([3.] * 3, 7., 6e-7, [0.1] * 2, 20, 0.5,
                           method='brentq')
+    with pytest.raises(ValueError):
+        pvsystem.v_from_i([3.] * 3, 7., 6e-7, [0.1] * 2, 20, 0.5,
+                          method='chandrupatla')
     with pytest.raises(ValueError):
         pvsystem.v_from_i([3.] * 3, np.array([7., 7.]), 6e-7, [0.1], 20, 0.5,
                           method='newton')
@@ -1437,7 +1444,8 @@ def test_mpp_floats():
         assert np.isclose(v, expected[k])
 
 
-def test_mpp_recombination():
+@pytest.mark.parametrize('method', ['brentq', 'newton', 'chandrupatla'])
+def test_mpp_recombination(method):
     """test max_power_point"""
     pvsyst_fs_495 = get_pvsyst_fs_495()
     IL, I0, Rs, Rsh, nNsVth = pvsystem.calcparams_pvsyst(
@@ -1455,7 +1463,7 @@ def test_mpp_recombination():
         IL, I0, Rs, Rsh, nNsVth,
         d2mutau=pvsyst_fs_495['d2mutau'],
         NsVbi=VOLTAGE_BUILTIN*pvsyst_fs_495['cells_in_series'],
-        method='brentq')
+        method=method)
     expected_imp = pvsyst_fs_495['I_mp_ref']
     expected_vmp = pvsyst_fs_495['V_mp_ref']
     expected_pmp = expected_imp*expected_vmp
@@ -1465,46 +1473,35 @@ def test_mpp_recombination():
     assert isinstance(out, dict)
     for k, v in out.items():
         assert np.isclose(v, expected[k], 0.01)
-    out = pvsystem.max_power_point(
-        IL, I0, Rs, Rsh, nNsVth,
-        d2mutau=pvsyst_fs_495['d2mutau'],
-        NsVbi=VOLTAGE_BUILTIN*pvsyst_fs_495['cells_in_series'],
-        method='newton')
-    for k, v in out.items():
-        assert np.isclose(v, expected[k], 0.01)
 
 
-def test_mpp_array():
+@pytest.mark.parametrize('method', ['brentq', 'newton', 'chandrupatla'])
+def test_mpp_array(method):
     """test max_power_point"""
     IL, I0, Rs, Rsh, nNsVth = (np.array([7, 7]), 6e-7, .1, 20, .5)
-    out = pvsystem.max_power_point(IL, I0, Rs, Rsh, nNsVth, method='brentq')
+    out = pvsystem.max_power_point(IL, I0, Rs, Rsh, nNsVth, method=method)
     expected = {'i_mp': [6.1362673597376753] * 2,
                 'v_mp': [6.2243393757884284] * 2,
                 'p_mp': [38.194210547580511] * 2}
     assert isinstance(out, dict)
     for k, v in out.items():
         assert np.allclose(v, expected[k])
-    out = pvsystem.max_power_point(IL, I0, Rs, Rsh, nNsVth, method='newton')
-    for k, v in out.items():
-        assert np.allclose(v, expected[k])
 
 
-def test_mpp_series():
+@pytest.mark.parametrize('method', ['brentq', 'newton', 'chandrupatla'])
+def test_mpp_series(method):
     """test max_power_point"""
     idx = ['2008-02-17T11:30:00-0800', '2008-02-17T12:30:00-0800']
     IL, I0, Rs, Rsh, nNsVth = (np.array([7, 7]), 6e-7, .1, 20, .5)
     IL = pd.Series(IL, index=idx)
-    out = pvsystem.max_power_point(IL, I0, Rs, Rsh, nNsVth, method='brentq')
+    out = pvsystem.max_power_point(IL, I0, Rs, Rsh, nNsVth, method=method)
     expected = pd.DataFrame({'i_mp': [6.1362673597376753] * 2,
                              'v_mp': [6.2243393757884284] * 2,
                              'p_mp': [38.194210547580511] * 2},
                             index=idx)
     assert isinstance(out, pd.DataFrame)
     for k, v in out.items():
         assert np.allclose(v, expected[k])
-    out = pvsystem.max_power_point(IL, I0, Rs, Rsh, nNsVth, method='newton')
-    for k, v in out.items():
-        assert np.allclose(v, expected[k])
 
 
 def test_singlediode_series(cec_module_params):