WIP: Speed up singlediode._lambertw

benchmarks/benchmarks/singlediode.py

@@ -0,0 +1,41 @@
+"""
+ASV benchmarks for singlediode.py
+"""
+# slower MT implementation in numpy<=1.16
+from numpy.random import RandomState
+# replace with below when ASV migrates to numpy>=1.17
+# and replace 'rng.rang()' with 'rng()'
+# from numpy.random import Generator, MT19937
+from pvlib import singlediode as _singlediode
+
+
+def b88(params):
+    # for a fair comparison, need to also compute isc, voc, i_x and i_xx
+    isc = _singlediode.bishop88_i_from_v(0., *params)
+    voc = _singlediode.bishop88_v_from_i(0., *params)
+    imp, vmp, pmp = _singlediode.bishop88_mpp(*params)
+    ix = _singlediode.bishop88_i_from_v(vmp/2., *params)
+    ixx = _singlediode.bishop88_i_from_v((voc + vmp)/2., *params)
+    return imp, vmp, pmp, isc, voc, ix, ixx
+
+
+class SingleDiode:
+
+    def setup(self):
+        seed = 11471
+        rng = RandomState(seed)
+        nsamples = 10000
+        il = 9. * rng.rand(nsamples) + 1.  # 1.- 10. A
+        io = 10**(-9 + 3. * rng.rand(nsamples))  # 1e-9 to 1e-6 A
+        rs = 5. * rng.rand(nsamples) + 0.05  # 0.05 to 5.05 Ohm
+        rsh = 10**(2 + 2 * rng.rand(nsamples))  # 100 to 10000 Ohm
+        n = 1 + 0.7 * rng.rand(nsamples)  # 1.0 to 1.7
+        # 72 cells in series, roughly 25C Tcell
+        nNsVth = 72 * n * 0.025
+        self.params = (il, io, rs, rsh, nNsVth)
+
+    def time_bishop88(self):
+        b88(self.params)
+
+    def time_lambertw(self):
+        _singlediode._lambertw(*self.params)

docs/sphinx/source/user_guide/singlediode.rst

@@ -6,57 +6,116 @@ Single Diode Equation
 This section reviews the solutions to the single diode equation used in
 pvlib-python to generate an IV curve of a PV module.
 
-pvlib-python supports two ways to solve the single diode equation:
-
-1. Lambert W-Function
-2. Bishop's Algorithm
-
-The :func:`pvlib.pvsystem.singlediode` function allows the user to choose the
-method using the ``method`` keyword.
-
-Lambert W-Function
-------------------
-When ``method='lambertw'``, the Lambert W-function is used as previously shown
-by Jain, Kapoor [1, 2] and Hansen [3]. The following algorithm can be found on
-`Wikipedia: Theory of Solar Cells
-<https://en.wikipedia.org/wiki/Theory_of_solar_cells>`_, given the basic single
-diode model equation.
+The single diode equation describes the current-voltage characteristic of
+an ideal single-junction PV device:
 
 .. math::
 
    I = I_L - I_0 \left(\exp \left(\frac{V + I R_s}{n Ns V_{th}} \right) - 1 \right)
        - \frac{V + I R_s}{R_{sh}}
 
-Lambert W-function is the inverse of the function
-:math:`f \left( w \right) = w \exp \left( w \right)` or
-:math:`w = f^{-1} \left( w \exp \left( w \right) \right)` also given as
-:math:`w = W \left( w \exp \left( w \right) \right)`. Defining the following
-parameter, :math:`z`, is necessary to transform the single diode equation into
-a form that can be expressed as a Lambert W-function.
+where
+
+* :math:`I` is electrical current in Amps
+* :math:`V` is voltage in Amps
+* :math:`I_L` is the photocurrent in Amps
+* :math:`I_0` is the dark, or saturation current, in Amps
+* :math:`R_{sh}` is the shunt resistance in Ohm
+* :math:`R_s` is the series resistance in Ohm
+* :math:`n` is the diode (ideality) factor (unitless)
+* :math:`Ns` is the number of cells in series. Cells are assumed to be identical.
+* :math:`V_{th}` is the thermal voltage at each cell's junction, given by :math:`V_{th} = \frac{k}{q} T_K`,
+  where :math:`k` is the Boltzmann constant (J/K), :math:`q` is the elementary charge (Couloumb) and :math:`T_k`
+  is the cell temperature in K.
+
+pvlib-python supports two ways to solve the single diode equation:
+
+1. Using the Lambert W function
+2. Bishop's algorithm
+
+The :func:`pvlib.pvsystem.singlediode` function's ``method`` keyword allows the user to choose the solution method.
+
+Lambert W-Function
+------------------
+When ``method='lambertw'``, the Lambert W function is used to find solutions :math:`(V, I)`.
+The Lambert W function is the converse relation of the function :math:`f \left( w \right) = w \exp \left( w \right)`,
+that is, if :math:`y \exp \left( y \right) = x`, then :math:`y = W(x)`.
+As previously shown by Jain, Kapoor [1, 2] and Hansen [3], solutions to the single diode equation
+may be written in terms of :math:`W`. Define a variable :math:`\theta` as 
 
 .. math::
 
-   z = \frac{R_s I_0}{n Ns V_{th} \left(1 + \frac{R_s}{R_{sh}} \right)} \exp \left(
+   \theta = \frac{R_s I_0}{n Ns V_{th} \left(1 + \frac{R_s}{R_{sh}} \right)} \exp \left(
        \frac{R_s \left( I_L + I_0 \right) + V}{n Ns V_{th} \left(1 + \frac{R_s}{R_{sh}}\right)}
        \right)
 
-Then the module current can be solved using the Lambert W-function,
+Then the module current can be written as a function of voltage, using the Lambert W-function,
 :math:`W \left(z \right)`.
 
 .. math::
 
    I = \frac{I_L + I_0 - \frac{V}{R_{sh}}}{1 + \frac{R_s}{R_{sh}}}
-       - \frac{n Ns V_{th}}{R_s} W \left(z \right)
+       - \frac{n Ns V_{th}}{R_s} W \left(\theta \right)
+
+
+Similarly, the voltage can be written as a function of current by defining a variable :math:`\psi` as
+
+.. math::
+
+   \psi = \frac{I_0 R_{sh}}{n Ns V_{th}} \exp \left(\frac{\left(I_L + I_0 - I\right) R_{sh}}{n Ns V_{th}} \right)
+
+Then
+
+.. math::
+
+   V = \left(I_L + I_0 - I\right) R_sh - I R_s - n Ns V_th W\left( \psi \right)
+
+When computing :math:`V = V\left( I \right)`, care must be taken to avoid overflow errors because the argument
+of the exponential function in :math:`\psi` can become large.
+
+The pvlib function :func:`pvlib.pvsystem.singlediode` uses these expressions :math:`I = I\left(V\right)` and
+:math:`V = V\left( I \right)` to calculate :math:`I_{sc}` and :math:`V_{oc}` respectively.
+
+To calculate the maximum power point :math:`\left( V_{mp}, I_{mp} \right)` a root-finding method is used. At the
+maximum power point, the derivative of power with respect to current (or voltage) is zero. Differentiating
+the equation :math:`P = V I` and evaluating at the maximum power current
+
+.. math::
+
+   0 = \frac{dP}{dI} \Bigr|_{I=I_{mp}} = \left(V + \frac{dV}{dI}\right) \Bigr|_{I=I_{mp}}
+
+obtains
+
+.. math::
+
+   \frac{dV}{dI}\Bigr|_{I=I_{mp}} = \frac{-V_{mp}}{I_{mp}}
+
+Differentiating :math:`V = V(I)` with respect to current, and applying the identity
+:math:`\frac{dW\left( x \right)}{dx} = \frac{W\left( x \right)}{x \left( 1 + W \left( x \right) \right)}` obtains
+
+.. math::
+
+   \frac{dV}{dI}\Bigr|_{I=I_{mp}} = -\left(R_s + \frac{R_{sh}}{1 + W\left( \psi \right)} \right)\Bigr|_{I=I_{mp}}
+
+Equating the two expressions for :math:`\frac{dV}{dI}\Bigr|_{I=I_{mp}}` and rearranging yields
+
+.. math::
+
+   \frac{\left(I_L + I_0 - I\right) R_{sh} - I R_s - n Ns V_{th} W\left( \psi \right)}{R_s + \frac{R_{sh}}{1 + W\left( \psi \right)}}\Bigr|_{I=I_{mp}} - I_{mp} = 0.
+
+The above equation is solved for :math:`I_{mp}` using Newton's method, and then :math:`V_{mp} = V \left( I_{mp} \right)` is computed.
 
 
 Bishop's Algorithm
 ------------------
 The function :func:`pvlib.singlediode.bishop88` uses an explicit solution [4]
 that finds points on the IV curve by first solving for pairs :math:`(V_d, I)`
-where :math:`V_d` is the diode voltage :math:`V_d = V + I*Rs`. Then the voltage
-is backed out from :math:`V_d`. Points with specific voltage, such as open
-circuit, are located using the bisection search method, ``brentq``, bounded
-by a zero diode voltage and an estimate of open circuit voltage given by
+where :math:`V_d` is the diode voltage :math:`V_d = V + I R_s`. Then the voltage
+is backed out from :math:`V_d`. Points with specific voltage or current are located
+using either Newton's or Brent's method, ``method='newton'`` or ``method='brentq'``,
+respectvely.
+
+For example, to find the open circuit voltage, we start with an estimate given by
 
 .. math::
 

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

@@ -18,6 +18,8 @@ Enhancements
 ~~~~~~~~~~~~
 * Added the optional `string_factor` parameter to
   :py:func:`pvlib.snow.loss_townsend` (:issue:`1636`, :pull:`1653`)
+* Improved performance of :py:func:`pvlib.pvsystem.singlediode` when using
+  ``method='lambertw'`` (:issue:`1649`, :pull:`1661`)
 
 Bug fixes
 ~~~~~~~~~

pvlib/singlediode.py

@@ -4,7 +4,6 @@
 
 from functools import partial
 import numpy as np
-from pvlib.tools import _golden_sect_DataFrame
 
 from scipy.optimize import brentq, newton
 from scipy.special import lambertw
@@ -640,16 +639,20 @@ def _lambertw(photocurrent, saturation_current, resistance_series,
     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)
+    conductance_shunt = 1. / resistance_shunt
+    il, io, rs, gsh, a = \
+        np.broadcast_arrays(photocurrent, saturation_current,
+                            resistance_series, conductance_shunt, nNsVth)
+
+    # Compute maximum power point quantities
+    params = (il, io, rs, gsh, a)
+    imp_est = 0.8 * il
+    args, i0 = _prepare_newton_inputs((), params, imp_est)
+    i_mp = newton(func=_imp_zero, x0=i0, fprime=_imp_zero_prime,
+                  args=args)
+    v_mp = _lambertw_v_from_i(resistance_shunt, resistance_series, nNsVth,
+                              i_mp, saturation_current, photocurrent)
+    p_mp = i_mp * v_mp
 
     # Find Imp using Lambert W
     i_mp = _lambertw_i_from_v(resistance_shunt, resistance_series, nNsVth,
@@ -679,6 +682,105 @@ def _lambertw(photocurrent, saturation_current, resistance_series,
     return out
 
 
+def _w_psi(i, il, io, gsh, a):
+    ''' Computes W(psi), where psi = io * rsh / a exp((il + io - i) * rsh / a)
+    This term is part of the equation V=V(I) solving the single diode equation
+    V = (il + io - i)*rsh - i*rs - a W(psi)
+
+    Parameters
+    ----------
+    i : numeric
+        Current (A)
+    il : numeric
+        Photocurrent (A)
+    io : numeric
+        Saturation current (A)
+    gsh : numeric
+        Shunt conductance (1/Ohm)
+    a : numeric
+        The product n*Ns*Vth (V).
+
+    Returns
+    -------
+    lambertwterm : numeric
+        The value of W(psi)
+
+    '''
+    with np.errstate(over='ignore'):
+        argW = (io / (gsh * a) *
+                np.exp((-i + il + io) /
+                       (gsh * a)))
+    lambertwterm = np.array(lambertw(argW).real)
+
+    idx_inf = np.logical_not(np.isfinite(lambertwterm))
+    if np.any(idx_inf):
+        # Calculate using log(argW) in case argW is really big
+        logargW = (np.log(io) - np.log(gsh) -
+                   np.log(a) +
+                   (-i + il + io) /
+                   (gsh * a))[idx_inf]
+
+        # Six iterations of Newton-Raphson method to solve
+        #  w+log(w)=logargW. The initial guess is w=logargW. Where direct
+        #  evaluation (above) results in NaN from overflow, 3 iterations
+        #  of Newton's method gives approximately 8 digits of precision.
+        w = logargW
+        for _ in range(0, 6):
+            w = w * (1. - np.log(w) + logargW) / (1. + w)
+        lambertwterm[idx_inf] = w
+    return lambertwterm
+
+
+def _imp_est(i, il, io, rs, gsh, a):
+    wma = _w_psi(i, il, io, gsh, a)
+    f = (il + io - i) / gsh - i * rs - a * wma
+    fprime = -rs - 1. / (gsh * (1. + wma))
+    return -f / fprime
+
+
+def _split_on_gsh(gsh):
+    # Determine indices where 0 < Gsh requires implicit model solution
+    idx_p = 0. < gsh
+    # Determine indices where 0 = Gsh allows explicit model solution
+    idx_z = 0. == gsh
+    return idx_p, idx_z
+
+
+def _imp_zero(i, il, io, rs, gsh, a):
+    ''' Root of this function is at imp
+    '''
+    idx_p, idx_z = _split_on_gsh(gsh)
+    res = np.full_like(i, np.nan, dtype=np.float64)
+
+    if np.any(idx_z):
+        # explicit solution for gsh=0
+        t = il[idx_z] + io[idx_z] - i[idx_z]
+        res[idx_z] = i[idx_z] / t - np.log(t / io[idx_z])
+    if np.any(idx_p):
+        res[idx_p] = _imp_est(i[idx_p], il[idx_p], io[idx_p], rs[idx_p],
+                              gsh[idx_p], a[idx_p]) - i[idx_p]
+    return res
+
+
+def _imp_zero_prime(i, il, io, rs, gsh, a):
+    ''' Derivative of _imp_zero with respect to current i
+    '''
+    idx_p, idx_z = _split_on_gsh(gsh)
+    res = np.full_like(i, np.nan, dtype=np.float64)
+    if np.any(idx_z):
+        # explicit solution for gsh=0
+        t = il[idx_z] + io[idx_z] - i[idx_z]
+        res[idx_z] = 2. / t + i[idx_z] / t**2.
+    if np.any(idx_p):
+        wma = _w_psi(i[idx_p], il[idx_p], io[idx_p], gsh[idx_p], a[idx_p])
+        f = (il[idx_p] + io[idx_p] - i[idx_p]) / gsh[idx_p] - \
+            i[idx_p] * rs[idx_p] - a[idx_p] * wma
+        fprime = -rs[idx_p] - 1. / (gsh[idx_p] * (1. + wma))
+        fprime2 = -1. / (gsh[idx_p]**2. * a[idx_p]) * wma / (1 + wma)**3.
+        res[idx_p] = f / fprime**2. * fprime2 - 2.
+    return res
+
+
 def _pwr_optfcn(df, loc):
     '''
     Function to find power from ``i_from_v``.