Add lambertw_pvlib

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

@@ -23,9 +23,12 @@ Bug fixes
 
 Enhancements
 ~~~~~~~~~~~~
-* Use ``k`` and ``cap_adjustment`` from :py:func:`pvlib.pvsystem.Array.module_parameters` in :py:func:`pvlib.pvsystem.PVSystem.pvwatts_dc`
+* Use ``k`` and ``cap_adjustment`` from :py:func:`pvlib.pvsystem.Array.module_parameters`
+  in :py:func:`pvlib.pvsystem.PVSystem.pvwatts_dc`
   (:issue:`2714`, :pull:`2715`)
-
+* Add :py:func:`pvlib.tools.lambertw_pvlib` to speed up calculations when using
+  LambertW single diode equation methods.
+  (:discussion:`2720`, :pull:`2723`)  
 
 Documentation
 ~~~~~~~~~~~~~

pvlib/ivtools/utils.py

@@ -544,3 +544,60 @@ def astm_e1036(v, i, imax_limits=(0.75, 1.15), vmax_limits=(0.75, 1.15),
     result['mp_fit'] = mp_fit
 
     return result
+
+
+def _log_lambertw(logx):
+    r'''Computes W(x) starting from log(x).
+
+    Parameters
+    ----------
+    logx : numeric
+        Log(x) of
+
+    Returns
+    -------
+    numeric
+        Lambert's W(x)
+
+    '''
+    # handles overflow cases, but results in nan for x <= 1
+    w = logx - np.log(logx)  # initial guess, w = log(x) - log(log(x))
+
+    for _ in range(0, 3):
+        # Newton's. Halley's is not substantially faster or more accurate
+        # because f''(w) = -1 / (w**2) is small for large w
+        w = w * (1. - np.log(w) + logx) / (1. + w)
+    return w
+
+
+def _lambertw_pvlib(x):
+    r'''Lambert's W function principal branch, :math:`W_0(x)`, for :math:`x`
+    real valued.
+
+    Parameters
+    ----------
+    x : np.array
+        Must be real numbers.
+
+    Returns
+    -------
+    np.array
+
+    '''
+    w = np.full_like(x, np.nan)
+    small = x <= 10
+    # for large x, solve 0 = f(w) = w + log(w) - log(x) using Newton's
+    w[~small] = _log_lambertw(np.log(x[~small]))
+
+    # w will contain nan for these numbers due to log(w) = log(log(x))
+    # for small x, solve 0 = g(w) = w * exp(w) - x using Halley's method
+    if any(small):
+        z = x[small]
+        g = np.log(x[small] + 1) - np.log(np.log(x[small] + 1) + 1)
+        for _ in range(0, 3):
+            expg = np.exp(g)
+            g = g - (g*expg - z) * (g + 1) / \
+                (expg * (g + 1)**2 - 0.5*(g + 2)*(expg*g - z))
+        w[small] = g
+
+    return w

tests/ivtools/test_utils.py

@@ -3,6 +3,7 @@
 import pytest
 from pvlib.ivtools.utils import _numdiff, rectify_iv_curve, astm_e1036
 from pvlib.ivtools.utils import _schumaker_qspline
+from pvlib.ivtools.utils import _lambertw_pvlib
 
 from tests.conftest import TESTS_DATA_DIR
 
@@ -171,3 +172,29 @@ def test_astm_e1036_fit_points(v_array, i_array):
                 'ff': 0.7520255886236707}
     result.pop('mp_fit')
     assert result == pytest.approx(expected)
+
+
+def test_lambertw_pvlib():
+    test_exp = np.arange(-10., 300, step=10)
+    test_x = 10.**test_exp
+    # known solution from scipy.special.lambertw
+    # scipy 1.7.1, python 3.13.1, numpy 2.3.5
+    expected = np.array([
+        9.9999999989999997e-11, 5.6714329040978384e-01,
+        2.0028685413304952e+01, 4.2306755091738395e+01,
+        6.4904633770046118e+01, 8.7630277151947183e+01,
+        1.1042491882731335e+02, 1.3326278259180333e+02,
+        1.5613026581351718e+02, 1.7901931374150624e+02,
+        2.0192476320084489e+02, 2.2484310644511851e+02,
+        2.4777185185877809e+02, 2.7070916610249782e+02,
+        2.9365366103997610e+02, 3.1660426041503479e+02,
+        3.3956011295458728e+02, 3.6252053376149752e+02,
+        3.8548496362161768e+02, 4.0845294003314166e+02,
+        4.3142407612718210e+02, 4.5439804503371403e+02,
+        4.7737456808796901e+02, 5.0035340579834485e+02,
+        5.2333435083468805e+02, 5.4631722251791496e+02,
+        5.6930186244110166e+02, 5.9228813095427859e+02,
+        6.1527590431628334e+02, 6.3826507236734335e+02,
+        6.6125553661218726e+02])
+    result = _lambertw_pvlib(test_x)
+    assert np.allclose(result, expected, rtol=1e-14)