diff --git a/benchmarks/benchmarks/singlediode.py b/benchmarks/benchmarks/singlediode.py new file mode 100644 index 0000000000..525460849e --- /dev/null +++ b/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) diff --git a/docs/sphinx/source/user_guide/singlediode.rst b/docs/sphinx/source/user_guide/singlediode.rst index 2fafa7d9f5..f48f66494b 100644 --- a/docs/sphinx/source/user_guide/singlediode.rst +++ b/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 -`_, 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 + I \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:: diff --git a/docs/sphinx/source/whatsnew/v0.9.5.rst b/docs/sphinx/source/whatsnew/v0.9.5.rst index 6c50f32642..87e0f4e2ac 100644 --- a/docs/sphinx/source/whatsnew/v0.9.5.rst +++ b/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 ~~~~~~~~~ diff --git a/pvlib/singlediode.py b/pvlib/singlediode.py index 9f5fd336ef..0b7093711a 100644 --- 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,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``.