From c294e40213dc10979feea9729b73c5256b54f267 Mon Sep 17 00:00:00 2001 From: Cliff Hansen Date: Mon, 13 Feb 2023 11:17:53 -0700 Subject: [PATCH 01/14] switch lambertw mpp to newton technique --- docs/sphinx/source/user_guide/singlediode.rst | 109 ++++++++++++++---- pvlib/singlediode.py | 93 +++++++++++++-- 2 files changed, 167 insertions(+), 35 deletions(-) diff --git a/docs/sphinx/source/user_guide/singlediode.rst b/docs/sphinx/source/user_guide/singlediode.rst index 2fafa7d9f5..0e4cc038b8 100644 --- 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 -`_, 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 +`_, 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 index 9f5fd336ef..792dc64874 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,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``. From 62a8996f4786d6bee25cee05db283a77d799b289 Mon Sep 17 00:00:00 2001 From: Cliff Hansen Date: Mon, 13 Feb 2023 14:19:56 -0700 Subject: [PATCH 02/14] handle Rsh=inf --- pvlib/singlediode.py | 72 +++++++++++++++++++++++++++++++------------- 1 file changed, 51 insertions(+), 21 deletions(-) diff --git a/pvlib/singlediode.py b/pvlib/singlediode.py index 792dc64874..a214166552 100644 --- a/pvlib/singlediode.py +++ b/pvlib/singlediode.py @@ -639,18 +639,21 @@ def _lambertw(photocurrent, saturation_current, resistance_series, v_oc = _lambertw_v_from_i(resistance_shunt, resistance_series, nNsVth, 0., saturation_current, photocurrent) - params = (photocurrent, saturation_current, resistance_series, - resistance_shunt, nNsVth) + 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 - imp_est = 0.8 * photocurrent + 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) + 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,7 +682,7 @@ def _lambertw(photocurrent, saturation_current, resistance_series, return out -def _w_psi(i, il, io, rsh, a): +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) @@ -692,8 +695,8 @@ def _w_psi(i, il, io, rsh, a): Photocurrent (A) io : numeric Saturation current (A) - rsh : numeric - Shunt resistance (Ohm) + gsh : numeric + Shunt conductance (1/Ohm) a : numeric The product n*Ns*Vth (V). @@ -703,7 +706,6 @@ def _w_psi(i, il, io, rsh, a): The value of W(psi) ''' - gsh = 1. / rsh with np.errstate(over='ignore'): argW = (io / (gsh * a) * np.exp((-i + il + io) / @@ -729,23 +731,51 @@ def _w_psi(i, il, io, rsh, a): 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) +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 _imp_zero(i, il, io, rs, rsh, a): - return _imp_est(i, il, io, rs, rsh, a) - i +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) -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. + 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): + 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): From a97b383f0c72916cf01036c485f5c992a3bc1bfa Mon Sep 17 00:00:00 2001 From: Cliff Hansen Date: Mon, 13 Feb 2023 14:20:27 -0700 Subject: [PATCH 03/14] document edits --- docs/sphinx/source/user_guide/singlediode.rst | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/docs/sphinx/source/user_guide/singlediode.rst b/docs/sphinx/source/user_guide/singlediode.rst index 0e4cc038b8..482f33093b 100644 --- a/docs/sphinx/source/user_guide/singlediode.rst +++ b/docs/sphinx/source/user_guide/singlediode.rst @@ -25,8 +25,8 @@ where * :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. + 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: @@ -64,19 +64,19 @@ Similarly, the voltage can be written as a function of current by defining a var .. 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) + \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) \right) Then .. math:: - V = \left(I_L + I_0 - I) R_sh} - I R_s - n Ns V_th W\left( \psi \right) + 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:`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. +: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. @@ -93,18 +93,18 @@ obtains \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 +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}} -Combining the two expressions for \frac{dV}{dI}\Bigr|_{I=I_{mp}} and rearranging yields +Combining the two expressions for :math:`\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. + \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. From d84c83c3c73f481348699912f5d5f21d6daece45 Mon Sep 17 00:00:00 2001 From: Cliff Hansen Date: Mon, 13 Feb 2023 14:46:11 -0700 Subject: [PATCH 04/14] pacify stickler --- pvlib/singlediode.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/pvlib/singlediode.py b/pvlib/singlediode.py index a214166552..0b7093711a 100644 --- a/pvlib/singlediode.py +++ b/pvlib/singlediode.py @@ -745,7 +745,7 @@ def _split_on_gsh(gsh): 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 ''' @@ -763,6 +763,8 @@ def _imp_zero(i, il, io, rs, gsh, a): 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): @@ -771,7 +773,8 @@ def _imp_zero_prime(i, il, io, rs, gsh, a): 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 + 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. From e6d0a33dde4b8ab00533c8d1459801edfd2e4ad6 Mon Sep 17 00:00:00 2001 From: Cliff Hansen Date: Mon, 13 Feb 2023 14:49:00 -0700 Subject: [PATCH 05/14] doc fixes --- docs/sphinx/source/user_guide/singlediode.rst | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/docs/sphinx/source/user_guide/singlediode.rst b/docs/sphinx/source/user_guide/singlediode.rst index 482f33093b..599fe237c1 100644 --- a/docs/sphinx/source/user_guide/singlediode.rst +++ b/docs/sphinx/source/user_guide/singlediode.rst @@ -25,8 +25,8 @@ where * :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. + 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: @@ -64,7 +64,7 @@ Similarly, the voltage can be written as a function of current by defining a var .. 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) \right) + \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 @@ -104,7 +104,7 @@ Combining the two expressions for :math:`\frac{dV}{dI}\Bigr|_{I=I_{mp}}` and rea .. 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. + \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. From e1aff50078e5c301f10c179b9017bffb2e9c1653 Mon Sep 17 00:00:00 2001 From: Cliff Hansen Date: Mon, 13 Feb 2023 14:51:02 -0700 Subject: [PATCH 06/14] whatsnew --- docs/sphinx/source/whatsnew/v0.9.5.rst | 2 ++ 1 file changed, 2 insertions(+) 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 ~~~~~~~~~ From d362d7bbd562058698a98ed82deef5e6146c14f6 Mon Sep 17 00:00:00 2001 From: Cliff Hansen Date: Mon, 13 Feb 2023 15:25:11 -0700 Subject: [PATCH 07/14] add a benchmark --- benchmarks/benchmarks/singlediode.py | 39 ++++++++++++++++++++++++++++ 1 file changed, 39 insertions(+) create mode 100644 benchmarks/benchmarks/singlediode.py diff --git a/benchmarks/benchmarks/singlediode.py b/benchmarks/benchmarks/singlediode.py new file mode 100644 index 0000000000..154b9195f5 --- /dev/null +++ b/benchmarks/benchmarks/singlediode.py @@ -0,0 +1,39 @@ +""" +ASV benchmarks for singlediode.py +""" +from numpy.random import Generator, MT19937 +from pvlib import singlediode as _singlediode + +seed = 11471 + +rng = Generator(MT19937(seed)) +base_params = (1., 5.e-9, 0.5, 2000., 72 * 1.1 * 0.025) +nsamples = 10000 + + +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, base_params, nsamples): + self.il = 9 * rng(nsamples) + 1. # 1.- 10. A + self.io = 10**(-9 + 3. * rng(nsamples)) # 1e-9 to 1e-6 A + self.rs = 5 * rng(nsamples) + 0.05 # 0.05 to 5.05 Ohm + self.rsh = 10**(2 + 2 * rng(nsamples)) # 100 to 10000 Ohm + self.n = 1 + 0.7 * rng(nsamples) # 1.0 to 1.7 + self.nNsVth = 72 * self.n * 0.025 # 72 cells in series, roughly 25C Tcell + self.params = (self.il, self.io, self.rs, self.rsh, self.nNsVth) + + def bishop88(self): + b88(*self.params) + + def lambertw(self): + _singlediode.lambertw(*self.params) From 063abf80e6a18316ccfd4a61cb4cf6c94cb8cc7f Mon Sep 17 00:00:00 2001 From: Cliff Hansen Date: Mon, 13 Feb 2023 15:51:39 -0700 Subject: [PATCH 08/14] stickler, doc edits --- benchmarks/benchmarks/singlediode.py | 5 +++-- docs/sphinx/source/user_guide/singlediode.rst | 11 ++++------- 2 files changed, 7 insertions(+), 9 deletions(-) diff --git a/benchmarks/benchmarks/singlediode.py b/benchmarks/benchmarks/singlediode.py index 154b9195f5..16681cadb3 100644 --- a/benchmarks/benchmarks/singlediode.py +++ b/benchmarks/benchmarks/singlediode.py @@ -28,8 +28,9 @@ def setup(self, base_params, nsamples): self.io = 10**(-9 + 3. * rng(nsamples)) # 1e-9 to 1e-6 A self.rs = 5 * rng(nsamples) + 0.05 # 0.05 to 5.05 Ohm self.rsh = 10**(2 + 2 * rng(nsamples)) # 100 to 10000 Ohm - self.n = 1 + 0.7 * rng(nsamples) # 1.0 to 1.7 - self.nNsVth = 72 * self.n * 0.025 # 72 cells in series, roughly 25C Tcell + self.n = 1 + 0.7 * rng(nsamples) # 1.0 to 1.7 + # 72 cells in series, roughly 25C Tcell + self.nNsVth = 72 * self.n * 0.025 self.params = (self.il, self.io, self.rs, self.rsh, self.nNsVth) def bishop88(self): diff --git a/docs/sphinx/source/user_guide/singlediode.rst b/docs/sphinx/source/user_guide/singlediode.rst index 599fe237c1..fb22e47500 100644 --- a/docs/sphinx/source/user_guide/singlediode.rst +++ b/docs/sphinx/source/user_guide/singlediode.rst @@ -25,8 +25,8 @@ where * :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. + 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: @@ -98,20 +98,17 @@ Differentiating :math:`V = V(I)` with respect to current, and applying the ident .. math:: - \frac{dV}{dI}\Bigr|_{I=I_{mp}} = -\left(R_s + \frac{R_{sh}}{1 + W\left( psi \right)} \right)\Bigr|_{I=I_{mp}} + \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 :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. + \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. -`Wikipedia: Theory of Solar Cells -`_, - Bishop's Algorithm ------------------ The function :func:`pvlib.singlediode.bishop88` uses an explicit solution [4] From 5a0c80ba385f614316c6ea9da26de7ca7d0812d9 Mon Sep 17 00:00:00 2001 From: Cliff Hansen Date: Tue, 14 Feb 2023 10:48:39 -0700 Subject: [PATCH 09/14] revert to numpy 1.16 --- benchmarks/benchmarks/singlediode.py | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/benchmarks/benchmarks/singlediode.py b/benchmarks/benchmarks/singlediode.py index 16681cadb3..d08d528a17 100644 --- a/benchmarks/benchmarks/singlediode.py +++ b/benchmarks/benchmarks/singlediode.py @@ -1,12 +1,16 @@ """ ASV benchmarks for singlediode.py """ -from numpy.random import Generator, MT19937 +# 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 seed = 11471 -rng = Generator(MT19937(seed)) +rng = RandomState(seed) base_params = (1., 5.e-9, 0.5, 2000., 72 * 1.1 * 0.025) nsamples = 10000 @@ -24,11 +28,11 @@ def b88(params): class SingleDiode: def setup(self, base_params, nsamples): - self.il = 9 * rng(nsamples) + 1. # 1.- 10. A - self.io = 10**(-9 + 3. * rng(nsamples)) # 1e-9 to 1e-6 A - self.rs = 5 * rng(nsamples) + 0.05 # 0.05 to 5.05 Ohm - self.rsh = 10**(2 + 2 * rng(nsamples)) # 100 to 10000 Ohm - self.n = 1 + 0.7 * rng(nsamples) # 1.0 to 1.7 + self.il = 9. * rng.rand(nsamples) + 1. # 1.- 10. A + self.io = 10**(-9 + 3. * rng.rand(nsamples)) # 1e-9 to 1e-6 A + self.rs = 5. * rng.rand(nsamples) + 0.05 # 0.05 to 5.05 Ohm + self.rsh = 10**(2 + 2 * rng.rand(nsamples)) # 100 to 10000 Ohm + self.n = 1 + 0.7 * rng.rand(nsamples) # 1.0 to 1.7 # 72 cells in series, roughly 25C Tcell self.nNsVth = 72 * self.n * 0.025 self.params = (self.il, self.io, self.rs, self.rsh, self.nNsVth) From 9e88b919d29a22ba61949031ec5c3cc5c5930b65 Mon Sep 17 00:00:00 2001 From: Cliff Hansen Date: Tue, 14 Feb 2023 14:08:15 -0700 Subject: [PATCH 10/14] rename benchmark functions --- benchmarks/benchmarks/singlediode.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/benchmarks/benchmarks/singlediode.py b/benchmarks/benchmarks/singlediode.py index d08d528a17..c37f56fb9a 100644 --- a/benchmarks/benchmarks/singlediode.py +++ b/benchmarks/benchmarks/singlediode.py @@ -37,8 +37,8 @@ def setup(self, base_params, nsamples): self.nNsVth = 72 * self.n * 0.025 self.params = (self.il, self.io, self.rs, self.rsh, self.nNsVth) - def bishop88(self): + def time_bishop88(self): b88(*self.params) - def lambertw(self): + def time_lambertw(self): _singlediode.lambertw(*self.params) From 91ecf1c7162120f009564c944d68e84305a0c5ab Mon Sep 17 00:00:00 2001 From: Cliff Hansen Date: Tue, 14 Feb 2023 14:28:04 -0700 Subject: [PATCH 11/14] more doc edits --- docs/sphinx/source/user_guide/singlediode.rst | 28 +++++++++---------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/docs/sphinx/source/user_guide/singlediode.rst b/docs/sphinx/source/user_guide/singlediode.rst index fb22e47500..508e5f2ede 100644 --- a/docs/sphinx/source/user_guide/singlediode.rst +++ b/docs/sphinx/source/user_guide/singlediode.rst @@ -25,9 +25,8 @@ where * :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. - + 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: @@ -41,7 +40,7 @@ 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)`. +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 @@ -73,14 +72,13 @@ Then 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:`phi` can become large. +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 +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:: @@ -100,11 +98,11 @@ Differentiating :math:`V = V(I)` with respect to current, and applying the ident \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 :math:`\frac{dV}{dI}\Bigr|_{I=I_{mp}}` and rearranging yields +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. + \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. @@ -113,10 +111,12 @@ 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:: From 60988858642d6eecbc1f0a404dcff0af8f976c6f Mon Sep 17 00:00:00 2001 From: Cliff Hansen Date: Tue, 14 Feb 2023 15:02:07 -0700 Subject: [PATCH 12/14] fix stuff --- benchmarks/benchmarks/singlediode.py | 25 ++++++++----------- docs/sphinx/source/user_guide/singlediode.rst | 5 ++-- 2 files changed, 13 insertions(+), 17 deletions(-) diff --git a/benchmarks/benchmarks/singlediode.py b/benchmarks/benchmarks/singlediode.py index c37f56fb9a..32806e281a 100644 --- a/benchmarks/benchmarks/singlediode.py +++ b/benchmarks/benchmarks/singlediode.py @@ -8,12 +8,6 @@ # from numpy.random import Generator, MT19937 from pvlib import singlediode as _singlediode -seed = 11471 - -rng = RandomState(seed) -base_params = (1., 5.e-9, 0.5, 2000., 72 * 1.1 * 0.025) -nsamples = 10000 - def b88(params): # for a fair comparison, need to also compute isc, voc, i_x and i_xx @@ -27,15 +21,18 @@ def b88(params): class SingleDiode: - def setup(self, base_params, nsamples): - self.il = 9. * rng.rand(nsamples) + 1. # 1.- 10. A - self.io = 10**(-9 + 3. * rng.rand(nsamples)) # 1e-9 to 1e-6 A - self.rs = 5. * rng.rand(nsamples) + 0.05 # 0.05 to 5.05 Ohm - self.rsh = 10**(2 + 2 * rng.rand(nsamples)) # 100 to 10000 Ohm - self.n = 1 + 0.7 * rng.rand(nsamples) # 1.0 to 1.7 + 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 - self.nNsVth = 72 * self.n * 0.025 - self.params = (self.il, self.io, self.rs, self.rsh, self.nNsVth) + nNsVth = 72 * n * 0.025 + self.params = (il, io, rs, rsh, nNsVth) def time_bishop88(self): b88(*self.params) diff --git a/docs/sphinx/source/user_guide/singlediode.rst b/docs/sphinx/source/user_guide/singlediode.rst index 508e5f2ede..c8242dcfc8 100644 --- a/docs/sphinx/source/user_guide/singlediode.rst +++ b/docs/sphinx/source/user_guide/singlediode.rst @@ -30,11 +30,10 @@ where pvlib-python supports two ways to solve the single diode equation: -1. Lamberts W function +1. Using the Lambert W function 2. Bishop's algorithm -The :func:`pvlib.pvsystem.singlediode` function allows the user to choose the -method using the ``method`` keyword. +The :func:`pvlib.pvsystem.singlediode` function's ``method`` keyword allows the user to choose the solution method. Lambert W-Function ------------------ From 0ca692c8c53f8cf0ca0089128da1ac9aa8de4a33 Mon Sep 17 00:00:00 2001 From: Cliff Hansen Date: Tue, 14 Feb 2023 16:03:15 -0700 Subject: [PATCH 13/14] fix errors --- benchmarks/benchmarks/singlediode.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/benchmarks/benchmarks/singlediode.py b/benchmarks/benchmarks/singlediode.py index 32806e281a..525460849e 100644 --- a/benchmarks/benchmarks/singlediode.py +++ b/benchmarks/benchmarks/singlediode.py @@ -35,7 +35,7 @@ def setup(self): self.params = (il, io, rs, rsh, nNsVth) def time_bishop88(self): - b88(*self.params) + b88(self.params) def time_lambertw(self): - _singlediode.lambertw(*self.params) + _singlediode._lambertw(*self.params) From 4da7d2fba9ab086296c337be82f81034fdcd965e Mon Sep 17 00:00:00 2001 From: Cliff Hansen Date: Wed, 15 Feb 2023 14:36:48 -0700 Subject: [PATCH 14/14] Apply suggestions from code review Co-authored-by: Kevin Anderson --- docs/sphinx/source/user_guide/singlediode.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/sphinx/source/user_guide/singlediode.rst b/docs/sphinx/source/user_guide/singlediode.rst index c8242dcfc8..f48f66494b 100644 --- a/docs/sphinx/source/user_guide/singlediode.rst +++ b/docs/sphinx/source/user_guide/singlediode.rst @@ -68,7 +68,7 @@ Then .. math:: - V = \left(I_L + I_0 - I\right) R_sh - I R_s - n Ns V_th W\left( \psi \right) + 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. @@ -82,7 +82,7 @@ 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}} + 0 = \frac{dP}{dI} \Bigr|_{I=I_{mp}} = \left(V + I \frac{dV}{dI}\right) \Bigr|_{I=I_{mp}} obtains