Source code for pvlib.singlediode

"""
Low-level functions for solving the single diode equation.
"""

import numpy as np
from pvlib.tools import _golden_sect_DataFrame

from scipy.optimize import brentq, newton
from scipy.special import lambertw

# newton method default parameters for this module
NEWTON_DEFAULT_PARAMS = {
    'tol': 1e-6,
    'maxiter': 100
}

# intrinsic voltage per cell junction for a:Si, CdTe, Mertens et al.
VOLTAGE_BUILTIN = 0.9  # [V]


[docs] def estimate_voc(photocurrent, saturation_current, nNsVth): """ Rough estimate of open circuit voltage useful for bounding searches for ``i`` of ``v`` when using :func:`~pvlib.pvsystem.singlediode`. Parameters ---------- photocurrent : numeric photo-generated current [A] saturation_current : numeric diode reverse saturation current [A] nNsVth : numeric product of thermal voltage ``Vth`` [V], diode ideality factor ``n``, and number of series cells ``Ns`` Returns ------- numeric rough estimate of open circuit voltage [V] Notes ----- Calculating the open circuit voltage, :math:`V_{oc}`, of an ideal device with infinite shunt resistance, :math:`R_{sh} \\to \\infty`, and zero series resistance, :math:`R_s = 0`, yields the following equation [1]. As an estimate of :math:`V_{oc}` it is useful as an upper bound for the bisection method. .. math:: V_{oc, est}=n Ns V_{th} \\log \\left( \\frac{I_L}{I_0} + 1 \\right) .. [1] http://www.pveducation.org/pvcdrom/open-circuit-voltage """ return nNsVth * np.log(np.asarray(photocurrent) / saturation_current + 1.0)
[docs] def bishop88(diode_voltage, photocurrent, saturation_current, resistance_series, resistance_shunt, nNsVth, d2mutau=0, NsVbi=np.inf, breakdown_factor=0., breakdown_voltage=-5.5, breakdown_exp=3.28, gradients=False): r""" Explicit calculation of points on the IV curve described by the single diode equation. Values are calculated as described in [1]_. The single diode equation with recombination current and reverse bias breakdown is .. math:: I = I_{L} - I_{0} \left (\exp \frac{V_{d}}{nN_{s}V_{th}} - 1 \right ) - \frac{V_{d}}{R_{sh}} - \frac{I_{L} \frac{d^{2}}{\mu \tau}}{N_{s} V_{bi} - V_{d}} - a \frac{V_{d}}{R_{sh}} \left (1 - \frac{V_{d}}{V_{br}} \right )^{-m} The input `diode_voltage` must be :math:`V + I R_{s}`. .. warning:: * Usage of ``d2mutau`` is required with PVSyst coefficients for cadmium-telluride (CdTe) and amorphous-silicon (a:Si) PV modules only. * Do not use ``d2mutau`` with CEC coefficients. Parameters ---------- diode_voltage : numeric diode voltage :math:`V_d` [V] photocurrent : numeric photo-generated current :math:`I_{L}` [A] saturation_current : numeric diode reverse saturation current :math:`I_{0}` [A] resistance_series : numeric series resistance :math:`R_{s}` [ohms] resistance_shunt: numeric shunt resistance :math:`R_{sh}` [ohms] nNsVth : numeric product of thermal voltage :math:`V_{th}` [V], diode ideality factor :math:`n`, and number of series cells :math:`N_{s}` [V] d2mutau : numeric, default 0 PVsyst parameter for cadmium-telluride (CdTe) and amorphous-silicon (a-Si) modules that accounts for recombination current in the intrinsic layer. The value is the ratio of intrinsic layer thickness squared :math:`d^2` to the diffusion length of charge carriers :math:`\mu \tau`. [V] NsVbi : numeric, default np.inf PVsyst parameter for cadmium-telluride (CdTe) and amorphous-silicon (a-Si) modules that is the product of the PV module number of series cells :math:`N_{s}` and the builtin voltage :math:`V_{bi}` of the intrinsic layer. [V]. breakdown_factor : float, default 0 fraction of ohmic current involved in avalanche breakdown :math:`a`. Default of 0 excludes the reverse bias term from the model. [unitless] breakdown_voltage : float, default -5.5 reverse breakdown voltage of the photovoltaic junction :math:`V_{br}` [V] breakdown_exp : float, default 3.28 avalanche breakdown exponent :math:`m` [unitless] gradients : bool False returns only I, V, and P. True also returns gradients Returns ------- tuple currents [A], voltages [V], power [W], and optionally :math:`\frac{dI}{dV_d}`, :math:`\frac{dV}{dV_d}`, :math:`\frac{dI}{dV}`, :math:`\frac{dP}{dV}`, and :math:`\frac{d^2 P}{dV dV_d}` Notes ----- The PVSyst thin-film recombination losses parameters ``d2mutau`` and ``NsVbi`` should only be applied to cadmium-telluride (CdTe) and amorphous- silicon (a-Si) PV modules, [2]_, [3]_. The builtin voltage :math:`V_{bi}` should account for all junctions. For example: tandem and triple junction cells would have builtin voltages of 1.8[V] and 2.7[V] respectively, based on the default of 0.9[V] for a single junction. The parameter ``NsVbi`` should only account for the number of series cells in a single parallel sub-string if the module has cells in parallel greater than 1. References ---------- .. [1] "Computer simulation of the effects of electrical mismatches in photovoltaic cell interconnection circuits" JW Bishop, Solar Cell (1988) :doi:`10.1016/0379-6787(88)90059-2` .. [2] "Improved equivalent circuit and Analytical Model for Amorphous Silicon Solar Cells and Modules." J. Mertens, et al., IEEE Transactions on Electron Devices, Vol 45, No 2, Feb 1998. :doi:`10.1109/16.658676` .. [3] "Performance assessment of a simulation model for PV modules of any available technology", André Mermoud and Thibault Lejeune, 25th EUPVSEC, 2010 :doi:`10.4229/25thEUPVSEC2010-4BV.1.114` """ # calculate recombination loss current where d2mutau > 0 is_recomb = d2mutau > 0 # True where there is thin-film recombination loss v_recomb = np.where(is_recomb, NsVbi - diode_voltage, np.inf) i_recomb = np.where(is_recomb, photocurrent * d2mutau / v_recomb, 0) # calculate temporary values to simplify calculations v_star = diode_voltage / nNsVth # non-dimensional diode voltage g_sh = 1.0 / resistance_shunt # conductance if breakdown_factor > 0: # reverse bias is considered brk_term = 1 - diode_voltage / breakdown_voltage brk_pwr = np.power(brk_term, -breakdown_exp) i_breakdown = breakdown_factor * diode_voltage * g_sh * brk_pwr else: i_breakdown = 0. i = (photocurrent - saturation_current * np.expm1(v_star) # noqa: W503 - diode_voltage * g_sh - i_recomb - i_breakdown) # noqa: W503 v = diode_voltage - i * resistance_series retval = (i, v, i*v) if gradients: # calculate recombination loss current gradients where d2mutau > 0 grad_i_recomb = np.where(is_recomb, i_recomb / v_recomb, 0) grad_2i_recomb = np.where(is_recomb, 2 * grad_i_recomb / v_recomb, 0) g_diode = saturation_current * np.exp(v_star) / nNsVth # conductance if breakdown_factor > 0: # reverse bias is considered brk_pwr_1 = np.power(brk_term, -breakdown_exp - 1) brk_pwr_2 = np.power(brk_term, -breakdown_exp - 2) brk_fctr = breakdown_factor * g_sh grad_i_brk = brk_fctr * (brk_pwr + diode_voltage * -breakdown_exp * brk_pwr_1) grad2i_brk = (brk_fctr * -breakdown_exp # noqa: W503 * (2 * brk_pwr_1 + diode_voltage # noqa: W503 * (-breakdown_exp - 1) * brk_pwr_2)) # noqa: W503 else: grad_i_brk = 0. grad2i_brk = 0. grad_i = -g_diode - g_sh - grad_i_recomb - grad_i_brk # di/dvd grad_v = 1.0 - grad_i * resistance_series # dv/dvd # dp/dv = d(iv)/dv = v * di/dv + i grad = grad_i / grad_v # di/dv grad_p = v * grad + i # dp/dv grad2i = -g_diode / nNsVth - grad_2i_recomb - grad2i_brk # d2i/dvd grad2v = -grad2i * resistance_series # d2v/dvd grad2p = ( grad_v * grad + v * (grad2i/grad_v - grad_i*grad2v/grad_v**2) + grad_i ) # d2p/dv/dvd retval += (grad_i, grad_v, grad, grad_p, grad2p) return retval
[docs] def bishop88_i_from_v(voltage, photocurrent, saturation_current, resistance_series, resistance_shunt, nNsVth, d2mutau=0, NsVbi=np.inf, breakdown_factor=0., breakdown_voltage=-5.5, breakdown_exp=3.28, method='newton', method_kwargs=None): """ Find current given any voltage. Parameters ---------- voltage : numeric voltage (V) in volts [V] photocurrent : numeric photogenerated current (Iph or IL) [A] saturation_current : numeric diode dark or saturation current (Io or Isat) [A] resistance_series : numeric series resistance (Rs) in [Ohm] resistance_shunt : numeric shunt resistance (Rsh) [Ohm] nNsVth : numeric product of diode ideality factor (n), number of series cells (Ns), and thermal voltage (Vth = k_b * T / q_e) in volts [V] d2mutau : numeric, default 0 PVsyst parameter for cadmium-telluride (CdTe) and amorphous-silicon (a-Si) modules that accounts for recombination current in the intrinsic layer. The value is the ratio of intrinsic layer thickness squared :math:`d^2` to the diffusion length of charge carriers :math:`\\mu \\tau`. [V] NsVbi : numeric, default np.inf PVsyst parameter for cadmium-telluride (CdTe) and amorphous-silicon (a-Si) modules that is the product of the PV module number of series cells ``Ns`` and the builtin voltage ``Vbi`` of the intrinsic layer. [V]. breakdown_factor : float, default 0 fraction of ohmic current involved in avalanche breakdown :math:`a`. Default of 0 excludes the reverse bias term from the model. [unitless] breakdown_voltage : float, default -5.5 reverse breakdown voltage of the photovoltaic junction :math:`V_{br}` [V] breakdown_exp : float, default 3.28 avalanche breakdown exponent :math:`m` [unitless] method : str, default 'newton' Either ``'newton'`` or ``'brentq'``. ''method'' must be ``'newton'`` if ``breakdown_factor`` is not 0. method_kwargs : dict, optional Keyword arguments passed to root finder method. See :py:func:`scipy:scipy.optimize.brentq` and :py:func:`scipy:scipy.optimize.newton` parameters. ``'full_output': True`` is allowed, and ``optimizer_output`` would be returned. See examples section. Returns ------- current : numeric current (I) at the specified voltage (V). [A] optimizer_output : tuple, optional, if specified in ``method_kwargs`` see root finder documentation for selected method. Found root is diode voltage in [1]_. Examples -------- Using the following arguments that may come from any `calcparams_.*` function in :py:mod:`pvlib.pvsystem`: >>> args = {'photocurrent': 1., 'saturation_current': 9e-10, 'nNsVth': 4., ... 'resistance_series': 4., 'resistance_shunt': 5000.0} Use default values: >>> i = bishop88_i_from_v(0.0, **args) Specify tolerances and maximum number of iterations: >>> i = bishop88_i_from_v(0.0, **args, method='newton', ... method_kwargs={'tol': 1e-3, 'rtol': 1e-3, 'maxiter': 20}) Retrieve full output from the root finder: >>> i, method_output = bishop88_i_from_v(0.0, **args, method='newton', ... method_kwargs={'full_output': True}) References ---------- .. [1] "Computer simulation of the effects of electrical mismatches in photovoltaic cell interconnection circuits" JW Bishop, Solar Cell (1988) :doi:`10.1016/0379-6787(88)90059-2` """ # collect args args = (photocurrent, saturation_current, resistance_series, resistance_shunt, nNsVth, d2mutau, NsVbi, breakdown_factor, breakdown_voltage, breakdown_exp) method = method.lower() # method_kwargs create dict if not provided # this pattern avoids bugs with Mutable Default Parameters if not method_kwargs: method_kwargs = {} def fv(x, v, *a): # calculate voltage residual given diode voltage "x" return bishop88(x, *a)[1] - v if method == 'brentq': # first bound the search using voc voc_est = estimate_voc(photocurrent, saturation_current, nNsVth) # start iteration slightly less than NsVbi when voc_est > NsVbi, to # avoid the asymptote at NsVbi xp = np.where(voc_est < NsVbi, voc_est, 0.9999*NsVbi) # brentq only works with scalar inputs, so we need a set up function # and np.vectorize to repeatedly call the optimizer with the right # arguments for possible array input def vd_from_brent(voc, v, iph, isat, rs, rsh, gamma, d2mutau, NsVbi, breakdown_factor, breakdown_voltage, breakdown_exp): return brentq(fv, 0.0, voc, args=(v, iph, isat, rs, rsh, gamma, d2mutau, NsVbi, breakdown_factor, breakdown_voltage, breakdown_exp), **method_kwargs) vd_from_brent_vectorized = np.vectorize(vd_from_brent) vd = vd_from_brent_vectorized(xp, voltage, *args) elif method == 'newton': x0, (voltage, *args), method_kwargs = \ _prepare_newton_inputs(voltage, (voltage, *args), method_kwargs) vd = newton(func=lambda x, *a: fv(x, voltage, *a), x0=x0, fprime=lambda x, *a: bishop88(x, *a, gradients=True)[4], args=args, **method_kwargs) else: raise NotImplementedError("Method '%s' isn't implemented" % method) # When 'full_output' parameter is specified, returned 'vd' is a tuple with # many elements, where the root is the first one. So we use it to output # the bishop88 result and return tuple(scalar, tuple with method results) if method_kwargs.get('full_output') is True: return (bishop88(vd[0], *args)[0], vd) else: return bishop88(vd, *args)[0]
[docs] def bishop88_v_from_i(current, photocurrent, saturation_current, resistance_series, resistance_shunt, nNsVth, d2mutau=0, NsVbi=np.inf, breakdown_factor=0., breakdown_voltage=-5.5, breakdown_exp=3.28, method='newton', method_kwargs=None): """ Find voltage given any current. Parameters ---------- current : numeric current (I) in amperes [A] photocurrent : numeric photogenerated current (Iph or IL) [A] saturation_current : numeric diode dark or saturation current (Io or Isat) [A] resistance_series : numeric series resistance (Rs) in [Ohm] resistance_shunt : numeric shunt resistance (Rsh) [Ohm] nNsVth : numeric product of diode ideality factor (n), number of series cells (Ns), and thermal voltage (Vth = k_b * T / q_e) in volts [V] d2mutau : numeric, default 0 PVsyst parameter for cadmium-telluride (CdTe) and amorphous-silicon (a-Si) modules that accounts for recombination current in the intrinsic layer. The value is the ratio of intrinsic layer thickness squared :math:`d^2` to the diffusion length of charge carriers :math:`\\mu \\tau`. [V] NsVbi : numeric, default np.inf PVsyst parameter for cadmium-telluride (CdTe) and amorphous-silicon (a-Si) modules that is the product of the PV module number of series cells ``Ns`` and the builtin voltage ``Vbi`` of the intrinsic layer. [V]. breakdown_factor : float, default 0 fraction of ohmic current involved in avalanche breakdown :math:`a`. Default of 0 excludes the reverse bias term from the model. [unitless] breakdown_voltage : float, default -5.5 reverse breakdown voltage of the photovoltaic junction :math:`V_{br}` [V] breakdown_exp : float, default 3.28 avalanche breakdown exponent :math:`m` [unitless] method : str, default 'newton' Either ``'newton'`` or ``'brentq'``. ''method'' must be ``'newton'`` if ``breakdown_factor`` is not 0. method_kwargs : dict, optional Keyword arguments passed to root finder method. See :py:func:`scipy:scipy.optimize.brentq` and :py:func:`scipy:scipy.optimize.newton` parameters. ``'full_output': True`` is allowed, and ``optimizer_output`` would be returned. See examples section. Returns ------- voltage : numeric voltage (V) at the specified current (I) in volts [V] optimizer_output : tuple, optional, if specified in ``method_kwargs`` see root finder documentation for selected method. Found root is diode voltage in [1]_. Examples -------- Using the following arguments that may come from any `calcparams_.*` function in :py:mod:`pvlib.pvsystem`: >>> args = {'photocurrent': 1., 'saturation_current': 9e-10, 'nNsVth': 4., ... 'resistance_series': 4., 'resistance_shunt': 5000.0} Use default values: >>> v = bishop88_v_from_i(0.0, **args) Specify tolerances and maximum number of iterations: >>> v = bishop88_v_from_i(0.0, **args, method='newton', ... method_kwargs={'tol': 1e-3, 'rtol': 1e-3, 'maxiter': 20}) Retrieve full output from the root finder: >>> v, method_output = bishop88_v_from_i(0.0, **args, method='newton', ... method_kwargs={'full_output': True}) References ---------- .. [1] "Computer simulation of the effects of electrical mismatches in photovoltaic cell interconnection circuits" JW Bishop, Solar Cell (1988) :doi:`10.1016/0379-6787(88)90059-2` """ # collect args args = (photocurrent, saturation_current, resistance_series, resistance_shunt, nNsVth, d2mutau, NsVbi, breakdown_factor, breakdown_voltage, breakdown_exp) method = method.lower() # method_kwargs create dict if not provided # this pattern avoids bugs with Mutable Default Parameters if not method_kwargs: method_kwargs = {} # first bound the search using voc voc_est = estimate_voc(photocurrent, saturation_current, nNsVth) # start iteration slightly less than NsVbi when voc_est > NsVbi, to avoid # the asymptote at NsVbi xp = np.where(voc_est < NsVbi, voc_est, 0.9999*NsVbi) def fi(x, i, *a): # calculate current residual given diode voltage "x" return bishop88(x, *a)[0] - i if method == 'brentq': # brentq only works with scalar inputs, so we need a set up function # and np.vectorize to repeatedly call the optimizer with the right # arguments for possible array input def vd_from_brent(voc, i, iph, isat, rs, rsh, gamma, d2mutau, NsVbi, breakdown_factor, breakdown_voltage, breakdown_exp): return brentq(fi, 0.0, voc, args=(i, iph, isat, rs, rsh, gamma, d2mutau, NsVbi, breakdown_factor, breakdown_voltage, breakdown_exp), **method_kwargs) vd_from_brent_vectorized = np.vectorize(vd_from_brent) vd = vd_from_brent_vectorized(xp, current, *args) elif method == 'newton': x0, (current, *args), method_kwargs = \ _prepare_newton_inputs(xp, (current, *args), method_kwargs) vd = newton(func=lambda x, *a: fi(x, current, *a), x0=x0, fprime=lambda x, *a: bishop88(x, *a, gradients=True)[3], args=args, **method_kwargs) else: raise NotImplementedError("Method '%s' isn't implemented" % method) # When 'full_output' parameter is specified, returned 'vd' is a tuple with # many elements, where the root is the first one. So we use it to output # the bishop88 result and return tuple(scalar, tuple with method results) if method_kwargs.get('full_output') is True: return (bishop88(vd[0], *args)[1], vd) else: return bishop88(vd, *args)[1]
[docs] def bishop88_mpp(photocurrent, saturation_current, resistance_series, resistance_shunt, nNsVth, d2mutau=0, NsVbi=np.inf, breakdown_factor=0., breakdown_voltage=-5.5, breakdown_exp=3.28, method='newton', method_kwargs=None): """ Find max power point. Parameters ---------- photocurrent : numeric photogenerated current (Iph or IL) [A] saturation_current : numeric diode dark or saturation current (Io or Isat) [A] resistance_series : numeric series resistance (Rs) in [Ohm] resistance_shunt : numeric shunt resistance (Rsh) [Ohm] nNsVth : numeric product of diode ideality factor (n), number of series cells (Ns), and thermal voltage (Vth = k_b * T / q_e) in volts [V] d2mutau : numeric, default 0 PVsyst parameter for cadmium-telluride (CdTe) and amorphous-silicon (a-Si) modules that accounts for recombination current in the intrinsic layer. The value is the ratio of intrinsic layer thickness squared :math:`d^2` to the diffusion length of charge carriers :math:`\\mu \\tau`. [V] NsVbi : numeric, default np.inf PVsyst parameter for cadmium-telluride (CdTe) and amorphous-silicon (a-Si) modules that is the product of the PV module number of series cells ``Ns`` and the builtin voltage ``Vbi`` of the intrinsic layer. [V]. breakdown_factor : numeric, default 0 fraction of ohmic current involved in avalanche breakdown :math:`a`. Default of 0 excludes the reverse bias term from the model. [unitless] breakdown_voltage : numeric, default -5.5 reverse breakdown voltage of the photovoltaic junction :math:`V_{br}` [V] breakdown_exp : numeric, default 3.28 avalanche breakdown exponent :math:`m` [unitless] method : str, default 'newton' Either ``'newton'`` or ``'brentq'``. ''method'' must be ``'newton'`` if ``breakdown_factor`` is not 0. method_kwargs : dict, optional Keyword arguments passed to root finder method. See :py:func:`scipy:scipy.optimize.brentq` and :py:func:`scipy:scipy.optimize.newton` parameters. ``'full_output': True`` is allowed, and ``optimizer_output`` would be returned. See examples section. Returns ------- tuple max power current ``i_mp`` [A], max power voltage ``v_mp`` [V], and max power ``p_mp`` [W] optimizer_output : tuple, optional, if specified in ``method_kwargs`` see root finder documentation for selected method. Found root is diode voltage in [1]_. Examples -------- Using the following arguments that may come from any `calcparams_.*` function in :py:mod:`pvlib.pvsystem`: >>> args = {'photocurrent': 1., 'saturation_current': 9e-10, 'nNsVth': 4., ... 'resistance_series': 4., 'resistance_shunt': 5000.0} Use default values: >>> i_mp, v_mp, p_mp = bishop88_mpp(**args) Specify tolerances and maximum number of iterations: >>> i_mp, v_mp, p_mp = bishop88_mpp(**args, method='newton', ... method_kwargs={'tol': 1e-3, 'rtol': 1e-3, 'maxiter': 20}) Retrieve full output from the root finder: >>> (i_mp, v_mp, p_mp), method_output = bishop88_mpp(**args, ... method='newton', method_kwargs={'full_output': True}) References ---------- .. [1] "Computer simulation of the effects of electrical mismatches in photovoltaic cell interconnection circuits" JW Bishop, Solar Cell (1988) :doi:`10.1016/0379-6787(88)90059-2` """ # collect args args = (photocurrent, saturation_current, resistance_series, resistance_shunt, nNsVth, d2mutau, NsVbi, breakdown_factor, breakdown_voltage, breakdown_exp) method = method.lower() # method_kwargs create dict if not provided # this pattern avoids bugs with Mutable Default Parameters if not method_kwargs: method_kwargs = {} # first bound the search using voc voc_est = estimate_voc(photocurrent, saturation_current, nNsVth) # start iteration slightly less than NsVbi when voc_est > NsVbi, to avoid # the asymptote at NsVbi xp = np.where(voc_est < NsVbi, voc_est, 0.9999*NsVbi) def fmpp(x, *a): return bishop88(x, *a, gradients=True)[6] if method == 'brentq': # break out arguments for numpy.vectorize to handle broadcasting vec_fun = np.vectorize( lambda voc, iph, isat, rs, rsh, gamma, d2mutau, NsVbi, vbr_a, vbr, vbr_exp: brentq(fmpp, 0.0, voc, args=(iph, isat, rs, rsh, gamma, d2mutau, NsVbi, vbr_a, vbr, vbr_exp), **method_kwargs) ) vd = vec_fun(xp, *args) elif method == 'newton': # make sure all args are numpy arrays if max size > 1 # if voc_est is an array, then make a copy to use for initial guess, v0 x0, args, method_kwargs = \ _prepare_newton_inputs(xp, args, method_kwargs) vd = newton(func=fmpp, x0=x0, fprime=lambda x, *a: bishop88(x, *a, gradients=True)[7], args=args, **method_kwargs) else: raise NotImplementedError("Method '%s' isn't implemented" % method) # When 'full_output' parameter is specified, returned 'vd' is a tuple with # many elements, where the root is the first one. So we use it to output # the bishop88 result and return # tuple(tuple with bishop88 solution, tuple with method results) if method_kwargs.get('full_output') is True: return (bishop88(vd[0], *args), vd) else: return bishop88(vd, *args)
def _shape_of_max_size(*args): return max(((np.size(a), np.shape(a)) for a in args), key=lambda t: t[0])[1] def _prepare_newton_inputs(x0, args, method_kwargs): """ Make inputs compatible with Scipy's newton by: - converting all arguments (`x0` and `args`) into numpy.ndarrays if any argument is not a scalar. - broadcasting the initial guess `x0` to the shape of the argument with the greatest size. Parameters ---------- x0: numeric Initial guess for newton. args: Iterable(numeric) Iterable of additional arguments to use in SciPy's newton. method_kwargs: dict Options to pass to newton. Returns ------- tuple The updated initial guess, arguments, and options for newton. """ if not (np.isscalar(x0) and all(map(np.isscalar, args))): args = tuple(map(np.asarray, args)) x0 = np.broadcast_to(x0, _shape_of_max_size(x0, *args)) # set abs tolerance and maxiter from method_kwargs if not provided # apply defaults, but giving priority to user-specified values method_kwargs = {**NEWTON_DEFAULT_PARAMS, **method_kwargs} return x0, args, method_kwargs def _lambertw_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))) # This transforms Gsh=1/Rsh, including ideal Rsh=np.inf into Gsh=0., which # is generally more numerically stable conductance_shunt = 1. / resistance_shunt # 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, Gsh, a = \ np.broadcast_arrays(current, photocurrent, saturation_current, resistance_series, conductance_shunt, nNsVth) # Intitalize output V (I might not be float64) V = np.full_like(I, np.nan, dtype=np.float64) # 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 # Explicit solutions where Gsh=0 if np.any(idx_z): V[idx_z] = a[idx_z] * np.log1p((IL[idx_z] - I[idx_z]) / I0[idx_z]) - \ I[idx_z] * Rs[idx_z] # Only compute using LambertW if there are cases with Gsh>0 if np.any(idx_p): # LambertW argument, cannot be float128, may overflow to np.inf # overflow is explicitly handled below, so ignore warnings here with np.errstate(over='ignore'): argW = (I0[idx_p] / (Gsh[idx_p] * a[idx_p]) * np.exp((-I[idx_p] + IL[idx_p] + I0[idx_p]) / (Gsh[idx_p] * a[idx_p]))) # lambertw typically returns complex value with zero imaginary part # may overflow to np.inf lambertwterm = lambertw(argW).real # Record indices where lambertw input overflowed output idx_inf = np.logical_not(np.isfinite(lambertwterm)) # Only re-compute LambertW if it overflowed if np.any(idx_inf): # Calculate using log(argW) in case argW is really big logargW = (np.log(I0[idx_p]) - np.log(Gsh[idx_p]) - np.log(a[idx_p]) + (-I[idx_p] + IL[idx_p] + I0[idx_p]) / (Gsh[idx_p] * a[idx_p]))[idx_inf] # Three 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, 3): w = w * (1. - np.log(w) + logargW) / (1. + w) lambertwterm[idx_inf] = w # Eqn. 3 in Jain and Kapoor, 2004 # V = -I*(Rs + Rsh) + IL*Rsh - a*lambertwterm + I0*Rsh # Recast in terms of Gsh=1/Rsh for better numerical stability. V[idx_p] = (IL[idx_p] + I0[idx_p] - I[idx_p]) / Gsh[idx_p] - \ I[idx_p] * Rs[idx_p] - a[idx_p] * lambertwterm if output_is_scalar: return V.item() else: return V def _lambertw_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))) # This transforms Gsh=1/Rsh, including ideal Rsh=np.inf into Gsh=0., which # is generally more numerically stable conductance_shunt = 1. / resistance_shunt # 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, Gsh, a = \ np.broadcast_arrays(voltage, photocurrent, saturation_current, resistance_series, conductance_shunt, nNsVth) # Intitalize output I (V might not be float64) I = np.full_like(V, np.nan, dtype=np.float64) # noqa: E741, N806 # Determine indices where 0 < Rs requires implicit model solution idx_p = 0. < Rs # Determine indices where 0 = Rs allows explicit model solution idx_z = 0. == Rs # Explicit solutions where Rs=0 if np.any(idx_z): I[idx_z] = IL[idx_z] - I0[idx_z] * np.expm1(V[idx_z] / a[idx_z]) - \ Gsh[idx_z] * V[idx_z] # Only compute using LambertW if there are cases with Rs>0 # Does NOT handle possibility of overflow, github issue 298 if np.any(idx_p): # LambertW argument, cannot be float128, may overflow to np.inf argW = Rs[idx_p] * I0[idx_p] / ( a[idx_p] * (Rs[idx_p] * Gsh[idx_p] + 1.)) * \ np.exp((Rs[idx_p] * (IL[idx_p] + I0[idx_p]) + V[idx_p]) / (a[idx_p] * (Rs[idx_p] * Gsh[idx_p] + 1.))) # lambertw typically returns complex value with zero imaginary part # may overflow to np.inf lambertwterm = lambertw(argW).real # Eqn. 2 in Jain and Kapoor, 2004 # I = -V/(Rs + Rsh) - (a/Rs)*lambertwterm + Rsh*(IL + I0)/(Rs + Rsh) # Recast in terms of Gsh=1/Rsh for better numerical stability. I[idx_p] = (IL[idx_p] + I0[idx_p] - V[idx_p] * Gsh[idx_p]) / \ (Rs[idx_p] * Gsh[idx_p] + 1.) - ( a[idx_p] / Rs[idx_p]) * lambertwterm if output_is_scalar: return I.item() else: return I def _lambertw(photocurrent, saturation_current, resistance_series, resistance_shunt, nNsVth, ivcurve_pnts=None): # 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) # Compute open circuit voltage v_oc = _lambertw_v_from_i(0., **params) # Set small elements <0 in v_oc to 0 if isinstance(v_oc, np.ndarray): v_oc[(v_oc < 0) & (v_oc > -1e-12)] = 0. elif isinstance(v_oc, (float, int)): if v_oc < 0 and v_oc > -1e-12: v_oc = 0. # 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) # Find Imp using Lambert W i_mp = _lambertw_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_xx = _lambertw_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) # create ivcurve if ivcurve_pnts: 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 out += (ivcurve_i, ivcurve_v) return out def _pwr_optfcn(df, loc): ''' Function to find power from ``i_from_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]