Add "logwright" SDE solver functions

docs/sphinx/source/user_guide/singlediode.rst

@@ -49,6 +49,8 @@ Then the module current can be solved using the Lambert W-function,
        - \frac{n Ns V_{th}}{R_s} W \left(z \right)
 
 
+TODO
+
 Bishop's Algorithm
 ------------------
 The function :func:`pvlib.singlediode.bishop88` uses an explicit solution [4]

pvlib/pvsystem.py

@@ -2405,7 +2405,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'``, ``'logwright'``, ``'newton'``, or
+        ``'brentq'``.
 
     Returns
     -------
@@ -2441,6 +2442,8 @@ def singlediode(photocurrent, saturation_current, resistance_series,
     implicit diode equation utilizes the Lambert W function to obtain an
     explicit function of :math:`V=f(I)` and :math:`I=f(V)` as shown in [2]_.
 
+    If the method is ``'logwright'`` TODO
+
     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.
@@ -2482,8 +2485,8 @@ def singlediode(photocurrent, saturation_current, resistance_series,
             resistance_shunt, nNsVth)  # collect args
     # Calculate points on the IV curve using the LambertW solution to the
     # single diode equation
-    if method.lower() == 'lambertw':
-        out = _singlediode._lambertw(*args, ivcurve_pnts)
+    if method.lower() in ['lambertw', 'logwright']:
+        out = _singlediode._lambertw(*args, ivcurve_pnts, how=method)
         points = out[:7]
         if ivcurve_pnts:
             ivcurve_i, ivcurve_v = out[7:]
@@ -2651,8 +2654,8 @@ def v_from_i(current, photocurrent, saturation_current, resistance_series,
         0 < nNsVth
 
     method : str
-        Method to use: ``'lambertw'``, ``'newton'``, or ``'brentq'``. *Note*:
-        ``'brentq'`` is limited to 1st quadrant only.
+        Method to use: ``'lambertw'``, ``'logwright'``, ``'newton'``,
+        or ``'brentq'``. *Note*: ``'brentq'`` is limited to 1st quadrant only.
 
     Returns
     -------
@@ -2668,6 +2671,8 @@ def v_from_i(current, photocurrent, saturation_current, resistance_series,
             resistance_series, resistance_shunt, nNsVth)
     if method.lower() == 'lambertw':
         return _singlediode._lambertw_v_from_i(*args)
+    elif method.lower() == 'logwright':
+        return _singlediode._logwright_v_from_i(*args)
     else:
         # Calculate points on the IV curve using either 'newton' or 'brentq'
         # methods. Voltages are determined by first solving the single diode
@@ -2733,8 +2738,8 @@ def i_from_v(voltage, photocurrent, saturation_current, resistance_series,
         0 < nNsVth
 
     method : str
-        Method to use: ``'lambertw'``, ``'newton'``, or ``'brentq'``. *Note*:
-        ``'brentq'`` is limited to 1st quadrant only.
+        Method to use: ``'lambertw'``, ``'logwright'``, ``'newton'``,
+        or ``'brentq'``. *Note*: ``'brentq'`` is limited to 1st quadrant only.
 
     Returns
     -------
@@ -2750,6 +2755,8 @@ def i_from_v(voltage, photocurrent, saturation_current, resistance_series,
             resistance_series, resistance_shunt, nNsVth)
     if method.lower() == 'lambertw':
         return _singlediode._lambertw_i_from_v(*args)
+    elif method.lower() == 'logwright':
+        return _singlediode._logwright_i_from_v(*args)
     else:
         # Calculate points on the IV curve using either 'newton' or 'brentq'
         # methods. Voltages are determined by first solving the single diode

pvlib/singlediode.py

@@ -3,7 +3,7 @@
 """
 
 import numpy as np
-from pvlib.tools import _golden_sect_DataFrame
+from pvlib.tools import _golden_sect_DataFrame, _logwrightomega
 
 from scipy.optimize import brentq, newton
 from scipy.special import lambertw
@@ -770,18 +770,27 @@ def _lambertw_i_from_v(voltage, photocurrent, saturation_current,
 
 
 def _lambertw(photocurrent, saturation_current, resistance_series,
-              resistance_shunt, nNsVth, ivcurve_pnts=None):
+              resistance_shunt, nNsVth, ivcurve_pnts=None, how='lambertw'):
+    if how == 'lambertw':
+        f_i_from_v = _lambertw_i_from_v
+        f_v_from_i = _lambertw_v_from_i
+    elif how == 'logwright':
+        f_i_from_v = _logwright_i_from_v
+        f_v_from_i = _logwright_v_from_i
+    else:
+        raise ValueError(f'invalid method: {how}')
+
     # collect args
     params = {'photocurrent': photocurrent,
               'saturation_current': saturation_current,
               'resistance_series': resistance_series,
               'resistance_shunt': resistance_shunt, 'nNsVth': nNsVth}
 
     # Compute short circuit current
-    i_sc = _lambertw_i_from_v(0., **params)
+    i_sc = f_i_from_v(0., **params)
 
     # Compute open circuit voltage
-    v_oc = _lambertw_v_from_i(0., **params)
+    v_oc = f_v_from_i(0., **params)
 
     # Set small elements <0 in v_oc to 0
     if isinstance(v_oc, np.ndarray):
@@ -792,15 +801,16 @@ def _lambertw(photocurrent, saturation_current, resistance_series,
 
     # 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)
+    p_mp, v_mp = _golden_sect_DataFrame(params, 0., v_oc * 1.14,
+                                        _make_pwr_optfcn(f_i_from_v))
 
     # Find Imp using Lambert W
-    i_mp = _lambertw_i_from_v(v_mp, **params)
+    i_mp = f_i_from_v(v_mp, **params)
 
     # Find Ix and Ixx using Lambert W
-    i_x = _lambertw_i_from_v(0.5 * v_oc, **params)
+    i_x = f_i_from_v(0.5 * v_oc, **params)
 
-    i_xx = _lambertw_i_from_v(0.5 * (v_oc + v_mp), **params)
+    i_xx = f_i_from_v(0.5 * (v_oc + v_mp), **params)
 
     out = (i_sc, v_oc, i_mp, v_mp, p_mp, i_x, i_xx)
 
@@ -809,21 +819,71 @@ def _lambertw(photocurrent, saturation_current, resistance_series,
         ivcurve_v = (np.asarray(v_oc)[..., np.newaxis] *
                      np.linspace(0, 1, ivcurve_pnts))
 
-        ivcurve_i = _lambertw_i_from_v(ivcurve_v.T, **params).T
+        ivcurve_i = f_i_from_v(ivcurve_v.T, **params).T
 
         out += (ivcurve_i, ivcurve_v)
 
     return out
 
 
-def _pwr_optfcn(df, loc):
+def _make_pwr_optfcn(i_from_v):
     '''
     Function to find power from ``i_from_v``.
     '''
+    def _pwr_optfcn(df, loc):
+        current = i_from_v(df[loc], df['photocurrent'],
+                           df['saturation_current'],
+                           df['resistance_series'],
+                           df['resistance_shunt'], df['nNsVth'])
+
+        return current * df[loc]
+    return _pwr_optfcn
+
+
+def _logwright_v_from_i(current, photocurrent, saturation_current,
+                        resistance_series, resistance_shunt, nNsVth):
+    # Record if inputs were all scalar
+    output_is_scalar = all(map(np.isscalar,
+                               (current, photocurrent, saturation_current,
+                                resistance_series, resistance_shunt, nNsVth)))
+
+    # Ensure that we are working with read-only views of numpy arrays
+    # Turns Series into arrays so that we don't have to worry about
+    #  multidimensional broadcasting failing
+    I, IL, I0, Rs, Rsh, a = \
+        np.broadcast_arrays(current, photocurrent, saturation_current,
+                            resistance_series, resistance_shunt, nNsVth)
+
+    log_I0_Rsh_a = np.log(I0 * Rsh / a)
+    x = log_I0_Rsh_a + Rsh * (-I + IL + I0) / a
+    V = a * (_logwrightomega(x) - log_I0_Rsh_a) - I * Rs
+
+    if output_is_scalar:
+        return V.item()
+    else:
+        return V
 
-    current = _lambertw_i_from_v(df[loc], df['photocurrent'],
-                                 df['saturation_current'],
-                                 df['resistance_series'],
-                                 df['resistance_shunt'], df['nNsVth'])
 
-    return current * df[loc]
+def _logwright_i_from_v(voltage, photocurrent, saturation_current,
+                        resistance_series, resistance_shunt, nNsVth):
+    # Record if inputs were all scalar
+    output_is_scalar = all(map(np.isscalar,
+                               (voltage, photocurrent, saturation_current,
+                                resistance_series, resistance_shunt, nNsVth)))
+
+    # Ensure that we are working with read-only views of numpy arrays
+    # Turns Series into arrays so that we don't have to worry about
+    #  multidimensional broadcasting failing
+    V, IL, I0, Rs, Rsh, a = \
+        np.broadcast_arrays(voltage, photocurrent, saturation_current,
+                            resistance_series, resistance_shunt, nNsVth)
+
+    log_term = np.log(I0 * Rs * Rsh / (a * (Rs + Rsh)))
+    x = log_term + (Rsh / (Rs + Rsh)) * (Rs * (IL + I0) + V) / a
+    I = (a * (_logwrightomega(x) - log_term) - V) / Rs
+
+    if output_is_scalar:
+        return I.item()
+    else:
+        return I
+

pvlib/tests/test_tools.py

@@ -3,6 +3,7 @@
 from pvlib import tools
 import numpy as np
 import pandas as pd
+import scipy
 
 
 @pytest.mark.parametrize('keys, input_dict, expected', [
@@ -120,3 +121,18 @@ def test_get_pandas_index(args, args_idx):
         assert index is None
     else:
         pd.testing.assert_index_equal(args[args_idx].index, index)
+
+
+def test__logwrightomega():
+    x = np.concatenate([-(10.**np.arange(100, -100, -0.1)),
+                        10.**np.arange(-100, 100, 0.1)])
+    with np.errstate(divide='ignore'):
+        expected = np.log(scipy.special.wrightomega(x))
+
+    expected[x < -750] = x[x < -750]
+
+    actual = tools._logwrightomega(x, n=2)
+    np.testing.assert_allclose(actual, expected, atol=2e-11, rtol=2e-11)
+
+    actual = tools._logwrightomega(x, n=3)
+    np.testing.assert_allclose(actual, expected, atol=2e-15, rtol=4e-16)

pvlib/tools.py

@@ -494,3 +494,24 @@ def get_pandas_index(*args):
         (a.index for a in args if isinstance(a, (pd.DataFrame, pd.Series))),
         None
     )
+
+
+def _logwrightomega(x, n=3):
+    # a faster alternative to np.log(scipy.special.wrightomega(x)).
+    # this calculation also more robust in that it avoids underflow
+    # problems that np.log(scipy.specia.wrightomega()) experiences
+    # for x < ~-745.
+
+    # TODO see ref X
+
+    y = np.where(x <= -np.e, x, 1)
+    y = np.where((-np.e < x) & (x < np.e), -np.e + (1 + np.e) * (x + np.e) / (2 * np.e), y)
+    np.log(x, out=y, where=x >= np.e)
+
+    for _ in range(n):
+        ey = np.exp(y)
+        ey_plus_one = 1 + ey
+        y_ey_x = y + ey - x
+        y = y - 2 * (y_ey_x) * (ey_plus_one) / ((2 * (ey_plus_one)**2 - (y_ey_x)*ey))
+
+    return y