Use chandrupatla to find MPP in lambertw method (pvlib#2571)

* restart from PR2567

* left interval to neg., calculate i_mp

* relax tolerance for i_mp

* docstring work

* indent

* Apply suggestions from code review

Co-authored-by: Kevin Anderson <kevin.anderso@gmail.com>

* better inspection of values

* shorten string in test system

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

Co-authored-by: Kevin Anderson <kevin.anderso@gmail.com>

---------

Co-authored-by: Kevin Anderson <kevin.anderso@gmail.com>

Author: cwhanse

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

@@ -30,7 +30,8 @@ Enhancements
   :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`)
-
+* Accelerate :py:func:`~pvlib.pvsystem.singlediode` when scipy>=1.15 is
+  installed. (:issue:`2497`, :pull:`2571`)
 
 
 Documentation
@@ -56,4 +57,4 @@ Maintenance
 
 Contributors
 ~~~~~~~~~~~~
-
+* Cliff Hansen (:ghuser:`cwhanse`)

pvlib/pvsystem.py

@@ -2534,28 +2534,35 @@ def singlediode(photocurrent, saturation_current, resistance_series,
     explicit function of :math:`V=f(I)` and :math:`I=f(V)` as shown in [2]_.
 
     If the method is ``'newton'`` then the root-finding Newton-Raphson method
-    is used. It should be safe for well behaved IV-curves, but the ``'brentq'``
-    method is recommended for reliability.
+    is used. It should be safe for well-behaved IV curves, otherwise the
+    ``'chandrupatla``` or ``'brentq'`` methods are recommended for reliability.
 
     If the method is ``'brentq'`` then Brent's bisection search method is used
     that guarantees convergence by bounding the voltage between zero and
-    open-circuit.
+    open-circuit. ``'brentq'`` is generally slower than the other options.
+
+    If the method is ``'chandrupatla'`` then Chandrupatla's method is used
+    that guarantees convergence.
 
     References
     ----------
-    .. [1] S.R. Wenham, M.A. Green, M.E. Watt, "Applied Photovoltaics" ISBN
-       0 86758 909 4
+    .. [1] S. R. Wenham, M. A. Green, M. E. Watt, "Applied Photovoltaics",
+       Centre for Photovoltaic Devices and Systems, 1995. ISBN
+       0867589094
 
     .. [2] A. Jain, A. Kapoor, "Exact analytical solutions of the
        parameters of real solar cells using Lambert W-function", Solar
-       Energy Materials and Solar Cells, 81 (2004) 269-277.
+       Energy Materials and Solar Cells, vol. 81 no. 2, pp. 269-277, Feb. 2004.
+       :doi:`10.1016/j.solmat.2003.11.018`.
 
-    .. [3] D. King et al, "Sandia Photovoltaic Array Performance Model",
-       SAND2004-3535, Sandia National Laboratories, Albuquerque, NM
+    .. [3] D. L. King, E. E. Boyson and J. A. Kratochvil "Photovoltaic Array
+       Performance Model", Sandia National Laboratories, Albuquerque, NM, USA.
+       Report SAND2004-3535, 2004.
 
-    .. [4] "Computer simulation of the effects of electrical mismatches in
-       photovoltaic cell interconnection circuits" JW Bishop, Solar Cell (1988)
-       https://doi.org/10.1016/0379-6787(88)90059-2
+    .. [4] J.W. Bishop, "Computer simulation of the effects of electrical
+       mismatches in photovoltaic cell interconnection circuits" Solar Cells,
+       vol. 25 no. 1, pp. 73-89, Oct. 1988.
+       :doi:`doi.org/10.1016/0379-6787(88)90059-2`
     """
     args = (photocurrent, saturation_current, resistance_series,
             resistance_shunt, nNsVth)  # collect args
@@ -2565,8 +2572,9 @@ def singlediode(photocurrent, saturation_current, resistance_series,
         out = _singlediode._lambertw(*args)
         points = out[:7]
     else:
-        # Calculate points on the IV curve using either 'newton' or 'brentq'
-        # methods. Voltages are determined by first solving the single diode
+        # Calculate points on the IV curve using Bishop's algorithm and solving
+        # with 'newton', 'brentq' or 'chandrupatla' method.
+        # Voltages are determined by first solving the single diode
         # equation for the diode voltage V_d then backing out voltage
         v_oc = _singlediode.bishop88_v_from_i(
             0.0, *args, method=method.lower()

pvlib/singlediode.py

@@ -141,18 +141,20 @@ def bishop88(diode_voltage, photocurrent, saturation_current,
 
     References
     ----------
-    .. [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`
-
-    .. [2] "Improved equivalent circuit and Analytical Model for Amorphous
-       Silicon Solar Cells and Modules." J. Mertens, et al., IEEE Transactions
-       on Electron Devices, Vol 45, No 2, Feb 1998.
+    .. [1] J.W. Bishop, "Computer simulation of the effects of electrical
+       mismatches in photovoltaic cell interconnection circuits" Solar Cells,
+       vol. 25 no. 1, pp. 73-89, Oct. 1988.
+       :doi:`doi.org/10.1016/0379-6787(88)90059-2`
+
+    .. [2] J. Merten, J. M. Asensi, C. Voz, A. V. Shah, R. Platz and J. Andreu,
+       "Improved equivalent circuit and Analytical Model for Amorphous
+       Silicon Solar Cells and Modules." , IEEE Transactions
+       on Electron Devices, vol. 45, no. 2, pp. 423-429, Feb 1998.
        :doi:`10.1109/16.658676`
 
-    .. [3] "Performance assessment of a simulation model for PV modules of any
-       available technology", André Mermoud and Thibault Lejeune, 25th EUPVSEC,
-       2010
+    .. [3] A. Mermoud and T. Lejeune, "Performance assessment of a simulation
+       model for PV modules of any available technology", In Proc. of the 25th
+       European PVSEC, Valencia, ES, 2010.
        :doi:`10.4229/25thEUPVSEC2010-4BV.1.114`
     """
     # calculate recombination loss current where d2mutau > 0
@@ -913,10 +915,25 @@ def _lambertw(photocurrent, saturation_current, resistance_series,
             v_oc = 0.
 
     # Find the voltage, v_mp, where the power is maximized.
-    # Start the golden section search at v_oc * 1.14
-    p_mp, v_mp = _golden_sect_DataFrame(params, 0., v_oc * 1.14, _pwr_optfcn)
+    # use scipy.elementwise if available
+    # remove try/except when scipy>=1.15, and golden mean is retired
+    try:
+        from scipy.optimize.elementwise import find_minimum
+        # left negative to insure strict inequality
+        init = (-1., 0.8*v_oc, v_oc)
+        res = find_minimum(_vmp_opt, init,
+                           args=(params['photocurrent'],
+                                 params['saturation_current'],
+                                 params['resistance_series'],
+                                 params['resistance_shunt'],
+                                 params['nNsVth'],))
+        v_mp = res.x
+        p_mp = -1.*res.f_x
+    except ModuleNotFoundError:
+        # switch to old golden section method
+        p_mp, v_mp = _golden_sect_DataFrame(params, 0., v_oc * 1.14,
+                                            _pwr_optfcn)
 
-    # Find Imp using Lambert W
     i_mp = _lambertw_i_from_v(v_mp, **params)
 
     # Find Ix and Ixx using Lambert W
@@ -938,6 +955,15 @@ def _lambertw(photocurrent, saturation_current, resistance_series,
     return out
 
 
+def _vmp_opt(v, iph, io, rs, rsh, nNsVth):
+    '''
+    Function to find negative of power from ``i_from_v``.
+    '''
+    current = _lambertw_i_from_v(v, iph, io, rs, rsh, nNsVth)
+
+    return -v * current
+
+
 def _pwr_optfcn(df, loc):
     '''
     Function to find power from ``i_from_v``.

tests/test_modelchain.py

@@ -140,7 +140,7 @@ def cec_dc_adr_ac_system(sam_data, cec_module_cs5p_220m,
                       module_parameters=module_parameters,
                       temperature_model_parameters=temp_model_params,
                       inverter_parameters=inverter,
-                      modules_per_string=14)
+                      modules_per_string=13)
     return system
 
 
@@ -1383,6 +1383,12 @@ def test_ac_models(sapm_dc_snl_ac_system, cec_dc_adr_ac_system,
         'adr': 'adr',
         'pvwatts': 'pvwatts',
         'pvwatts_multi': 'pvwatts'}
+    inverter_to_ac_model_param = {
+        'sandia': 'Paco',
+        'sandia_multi': 'Paco',
+        'adr': 'Pnom',
+        'pvwatts': 'pdc0',
+        'pvwatts_multi': 'pdc0'}
     ac_model = inverter_to_ac_model[inverter_model]
     system = ac_systems[inverter_model]
 
@@ -1399,7 +1405,12 @@ def test_ac_models(sapm_dc_snl_ac_system, cec_dc_adr_ac_system,
     assert m.call_count == 1
     assert isinstance(mc.results.ac, pd.Series)
     assert not mc.results.ac.empty
-    assert mc.results.ac.iloc[1] < 1
+    # irradiance 800 W/m2 at 1st timestamp
+    inv_param = mc.system.inverter_parameters[
+        inverter_to_ac_model_param[inverter_model]]
+    assert (mc.results.ac.iloc[0] > inv_param / 2.)
+    # irradiance 0 W/m2 at the 2nd timestamp
+    assert (np.isnan(mc.results.ac.iloc[1]) or (mc.results.ac.iloc[1] < 1))
 
 
 def test_ac_model_user_func(pvwatts_dc_pvwatts_ac_system, location, weather,

tests/test_singlediode.py

@@ -159,7 +159,7 @@ def test_singlediode_precision(method, precise_iv_curves):
 
     assert np.allclose(pc['i_sc'], outs['i_sc'], atol=1e-10, rtol=0)
     assert np.allclose(pc['v_oc'], outs['v_oc'], atol=1e-10, rtol=0)
-    assert np.allclose(pc['i_mp'], outs['i_mp'], atol=7e-8, rtol=0)
+    assert np.allclose(pc['i_mp'], outs['i_mp'], atol=1e-7, rtol=0)
     assert np.allclose(pc['v_mp'], outs['v_mp'], atol=1e-6, rtol=0)
     assert np.allclose(pc['p_mp'], outs['p_mp'], atol=1e-10, rtol=0)
     assert np.allclose(pc['i_x'], outs['i_x'], atol=1e-10, rtol=0)