From 535da3f39dd1e7caa8f21b456deb5ba58f2afd98 Mon Sep 17 00:00:00 2001 From: Kevin Anderson Date: Fri, 6 Oct 2023 16:57:58 -0400 Subject: [PATCH 1/4] create tools._logwrightomega function --- pvlib/tests/test_tools.py | 16 ++++++++++++++++ pvlib/tools.py | 21 +++++++++++++++++++++ 2 files changed, 37 insertions(+) diff --git a/pvlib/tests/test_tools.py b/pvlib/tests/test_tools.py index 583141a726..c7229ad1bc 100644 --- a/pvlib/tests/test_tools.py +++ b/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) diff --git a/pvlib/tools.py b/pvlib/tools.py index fe1b79a5f1..c42e66f9b6 100644 --- a/pvlib/tools.py +++ b/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 From fb5abdaa4273cca08f35a36726ba3ef8f09cfef1 Mon Sep 17 00:00:00 2001 From: Kevin Anderson Date: Fri, 6 Oct 2023 17:04:32 -0400 Subject: [PATCH 2/4] create _logwright_i_from_v and v_from_i; integrate into _lambertw --- pvlib/singlediode.py | 90 ++++++++++++++++++++++++++++++++++++-------- 1 file changed, 75 insertions(+), 15 deletions(-) diff --git a/pvlib/singlediode.py b/pvlib/singlediode.py index 203b20cdae..6b7dc157ac 100644 --- a/pvlib/singlediode.py +++ b/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,7 +770,16 @@ 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, @@ -778,10 +787,10 @@ def _lambertw(photocurrent, saturation_current, 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 + From 09f140f7516f02052728db436bc2f67e51e758d3 Mon Sep 17 00:00:00 2001 From: Kevin Anderson Date: Fri, 6 Oct 2023 17:05:43 -0400 Subject: [PATCH 3/4] add method='logwright' in pvlib.pvsystem functions --- pvlib/pvsystem.py | 21 ++++++++++++++------- 1 file changed, 14 insertions(+), 7 deletions(-) diff --git a/pvlib/pvsystem.py b/pvlib/pvsystem.py index f0c1cd2cda..78345153ac 100644 --- a/pvlib/pvsystem.py +++ b/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 From e84c073a6820bd73057551b417d75a44e80694e9 Mon Sep 17 00:00:00 2001 From: Kevin Anderson Date: Fri, 6 Oct 2023 17:05:50 -0400 Subject: [PATCH 4/4] doc stub --- docs/sphinx/source/user_guide/singlediode.rst | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/sphinx/source/user_guide/singlediode.rst b/docs/sphinx/source/user_guide/singlediode.rst index 2fafa7d9f5..ca0b67a1b4 100644 --- a/docs/sphinx/source/user_guide/singlediode.rst +++ b/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]