Skip to content

Add "logwright" SDE solver functions #1884

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
wants to merge 4 commits into
base: main
Choose a base branch
from
Draft
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions docs/sphinx/source/user_guide/singlediode.rst
Original file line number Diff line number Diff line change
@@ -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]
21 changes: 14 additions & 7 deletions pvlib/pvsystem.py
Original file line number Diff line number Diff line change
@@ -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
90 changes: 75 additions & 15 deletions pvlib/singlediode.py
Original file line number Diff line number Diff line change
@@ -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,18 +770,27 @@ 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,
'resistance_series': 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

16 changes: 16 additions & 0 deletions pvlib/tests/test_tools.py
Original file line number Diff line number Diff line change
@@ -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)
21 changes: 21 additions & 0 deletions pvlib/tools.py
Original file line number Diff line number Diff line change
@@ -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