Speedup of singlediode._lambertw function

[toddkarin]

I have an application where I need to evaluate _lambertw many many times. Even a small speedup would be helpful. Currently _lambertw calculates IV curve points like I_x, I_xx that I don't need. What about adding an input flag that allows us to skip the evaluation of some parameters?

The version I wrote to do this is below, but there are other ways this could be done as well.

def _lambertw(photocurrent, saturation_current, resistance_series,
              resistance_shunt, nNsVth, ivcurve_pnts=None,
              calculate_all=False):
    if calculate_all:
        # Compute short circuit current
        i_sc = _lambertw_i_from_v(resistance_shunt, resistance_series, nNsVth,
                                  0.,
                                  saturation_current, photocurrent)

    # Compute open circuit voltage
    v_oc = _lambertw_v_from_i(resistance_shunt, resistance_series, nNsVth, 0.,
                              saturation_current, photocurrent)

    params = {'r_sh': resistance_shunt,
              'r_s': resistance_series,
              'nNsVth': nNsVth,
              'i_0': saturation_current,
              'i_l': photocurrent}

    # 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)

    # Find Imp using Lambert W
    i_mp = _lambertw_i_from_v(resistance_shunt, resistance_series, nNsVth,
                              v_mp, saturation_current, photocurrent)

    if calculate_all:
        # Find Ix and Ixx using Lambert W
        i_x = _lambertw_i_from_v(resistance_shunt, resistance_series, nNsVth,
                                 0.5 * v_oc, saturation_current, photocurrent)

        i_xx = _lambertw_i_from_v(resistance_shunt, resistance_series, nNsVth,
                                  0.5 * (v_oc + v_mp), saturation_current,
                                  photocurrent)

    if calculate_all:
        out = (i_sc, v_oc, i_mp, v_mp, p_mp, i_x, i_xx)
    else:
        out = (v_oc, i_mp, v_mp, p_mp)

    # create ivcurve
    if ivcurve_pnts:
        ivcurve_v = (np.asarray(v_oc)[..., np.newaxis] *
                     np.linspace(0, 1, ivcurve_pnts))

        ivcurve_i = _lambertw_i_from_v(resistance_shunt, resistance_series,
                                       nNsVth, ivcurve_v.T, saturation_current,
                                       photocurrent).T

        out += (ivcurve_i, ivcurve_v)

    return out

[cwhanse]

Not opposed but not excited about the specific proposal. What about a helper function _lambertw_mp that calculates the maximum power voltage and current? Then you could use _lambertw_v_from_i etc. to get Voc and Isc if desired.

[mikofski]

Or just use Newton, or the explicit method. Why do you need Lambert W? What is your use case?

• Trying to get a full of curve? You won't beat the performance of pvlib.singlediode.bishop88, you can also access this from pvsystem.singlediode by using method='newton'.
• Need just principal points? You can use the i_from_v and max power point functions again either from pvsystem with method='newton' or directly with the bishop88_ in singlediode

SciPy was specifically modified to allow us to use a vectorised version of newton to speed this up. The SciPy Brent method was also Cythonized but it hasn't been applied yet, but there's an example in this gist

See these issues

• singlediode is slow #408
• Single diode model helper function to wrap Lambert w and other methods #410

[toddkarin]

Thanks for the suggestions. I wrote a quick performance testing script and found that method='newton' is 2.4x faster than method='lambertw'. Is there a reason to prefer method='lambertw'?

import numpy as np
from time import time
from pvlib.pvsystem import singlediode

photocurrent = np.random.random(100000)*9
saturation_current = 1e-9 + np.zeros_like(photocurrent)
resistance_series = 0.5 + np.zeros_like(photocurrent)
resistance_shunt = 100 + np.zeros_like(photocurrent)
nNsVth = 1.8 + np.zeros_like(photocurrent)

start_time = time()
out = singlediode(photocurrent,saturation_current,resistance_series,resistance_shunt,
            nNsVth,method='lambertw')
print('Lambertw method elapsed time: {:.3f} s'.format( time() - start_time))

start_time = time()
out = singlediode(photocurrent,saturation_current,resistance_series,resistance_shunt,
            nNsVth,method='newton')
print('Newton method elapsed time: {:.3f} s'.format( time() - start_time))

# Lambertw method elapsed time: 0.685 s
# Newton method elapsed time: 0.281 s

[mikofski]

I believe it's because of the possibility that Newton could have trouble converging, but I also think that concern may have been overstated, b/c these curves are all monotonic. But I think this concern is also why we provided the Brent method in addition to Newton, b/c it's guaranteed to converge for monotonic curves to any point between the bounds of Isc and Voc. However the Brent method is still using numpy.vectorize (a broadcasting loop) instead of the Brent method from scipy.cython_optimize. I never followed up and completed the work, b/c newton seems to work well enough. It's possible that a Cythonized Brent could be even faster than the vectorized Newton method. In benchmarks of scipy.cython_optimize the Cython version was over 30x faster than looping. If there's an appetite to add the Brent method, I can look at making a PR, it shouldn't require too much effort, but I'm also okay if someone else wants to do it, or just benchmark it.

[cwhanse]

My personal testing has shown that if you are after the most accurate solutions, the LambertW method is more reliable than the Newton/Brentq methods, but, the difference in accuracy is typically in the 11th and lower significant digits. I doubt most users would notice, or care, about these differences.

[cwhanse]

And solving 10^5 diode equations in a fraction of a second isn't bad performance in my opinion :)

[toddkarin]

Thanks for this advice! I think the Newton method will be accurate enough for my application. We can then avoid making any changes to the library.