switch lambertw mpp to newton technique

pvlib · cwhanse · Feb 13, 2023 · Feb 13, 2023 · Feb 13, 2023 · Feb 13, 2023
commit c294e40213dc10979feea9729b73c5256b54f267
diff --git a/docs/sphinx/source/user_guide/singlediode.rst b/docs/sphinx/source/user_guide/singlediode.rst
@@ -6,54 +6,117 @@ 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.
 
+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}}
+
+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. Lambert W-Function
-2. Bishop's Algorithm
+1. Lamberts 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.
-
-.. 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.
+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 `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) R_{sh}}{n Ns V_{th}} \right)
+
+Then
+
+.. math::
 
+   V = \left(I_L + I_0 - I) 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:`phi` 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)` above 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}}
+
+Combining the two expressions for \frac{dV}{dI}\Bigr|_{I=I_{mp}} and rearranging yields
+
+.. math::
+
+   \frac{\left(I_L + I_0 - I) R_sh} - I R_s - n Ns V_th W\left( \psi \right)}{\left(R_s + \frac{R_{sh}}{1 + W\left( psi \right)} \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.
+
+
+`Wikipedia: Theory of Solar Cells
+<https://en.wikipedia.org/wiki/Theory_of_solar_cells>`_, 
 
 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
+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

diff --git a/pvlib/singlediode.py b/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,17 +639,18 @@ 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)
-
+    params = (photocurrent, saturation_current, resistance_series,
+              resistance_shunt, nNsVth)
+
+    # Compute maximum power point quantities
+    imp_est = 0.8 * photocurrent
+    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,
                               v_mp, saturation_current, photocurrent)
@@ -679,6 +679,75 @@ def _lambertw(photocurrent, saturation_current, resistance_series,
     return out
 
 
+def _w_psi(i, il, io, rsh, 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)
+    rsh : numeric
+        Shunt resistance (Ohm)
+    a : numeric
+        The product n*Ns*Vth (V).
+
+    Returns
+    -------
+    lambertwterm : numeric
+        The value of W(psi)
+
+    '''
+    gsh = 1. / rsh
+    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, rsh, a):
+    wma = _w_psi(i, il, io, rsh, a)
+    f = (il + io - i) * rsh - i * rs - a * wma
+    fprime = -rs - rsh / (1 + wma)
+    return -f / fprime
+
+
+def _imp_zero(i, il, io, rs, rsh, a):
+    return _imp_est(i, il, io, rs, rsh, a) - i
+
+
+def _imp_zero_prime(i, il, io, rs, rsh, a):
+    wma = _w_psi(i, il, io, rsh, a)
+    f = (il + io - i) * rsh - i * rs - a * wma
+    fprime = -rs - rsh / (1 + wma)
+    fprime2 = -rsh**2. / a * wma / (1 + wma)**3.
+    return f / fprime**2. * fprime2 - 2.
+
+
 def _pwr_optfcn(df, loc):
     '''
     Function to find power from ``i_from_v``.