diff --git a/docs/sphinx/source/whatsnew/v0.15.1.rst b/docs/sphinx/source/whatsnew/v0.15.1.rst index f86fed987b..504efb588b 100644 --- a/docs/sphinx/source/whatsnew/v0.15.1.rst +++ b/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. + (:pull:`2723`) Documentation ~~~~~~~~~~~~~ diff --git a/pvlib/tools.py b/pvlib/tools.py index 6cb631f852..c0aecc6cca 100644 --- a/pvlib/tools.py +++ b/pvlib/tools.py @@ -586,3 +586,60 @@ def _file_context_manager(filename_or_object, mode='r', encoding=None): # otherwise, assume a filename or path context = open(str(filename_or_object), mode=mode, encoding=encoding) return context + + +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. + + Parameters + ---------- + x : numeric + Must be real numbers. + + Returns + ------- + numeric + Lambert's W(x). Principal branch only. + + ''' + 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 diff --git a/tests/test_tools.py b/tests/test_tools.py index 4b733ad711..d053b24c08 100644 --- a/tests/test_tools.py +++ b/tests/test_tools.py @@ -283,3 +283,28 @@ def test__file_context_manager(): buffer = StringIO("test content") with tools._file_context_manager(buffer) as obj: assert obj.read() == "test content" + + +def test_lambertw_pvlib(): + test_exp = np.arange(-10., 300, step=10) + test_x = 10.**test_exp + # known solution from scipy.special.lambertw + 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 = tools.lambertw_pvlib(test_x) + assert np.allclose(result, expected, rtol=1e-14)