Add Batzelis 2017 simple nonlinear PV model (#2563)

* add function

* docs

* tests

* add SDM parameter estimation and key point functions

* lint

* tests for the two additional functions

* whatsnew

* lint

* Apply suggestions from code review

Co-authored-by: Cliff Hansen <[email protected]>

* lint

* better handling of Rsh=np.inf

* rename parameters

* switch to non-normalized temp coeffs

* fix whatsnew

* changes from review

---------

Co-authored-by: Cliff Hansen <[email protected]>

Author: kandersolar

docs/sphinx/source/reference/pv_modeling/sdm.rst

@@ -17,6 +17,7 @@ Functions relevant for single diode models.
    pvsystem.v_from_i
    pvsystem.max_power_point
    ivtools.sdm.pvsyst_temperature_coeff
+   singlediode.batzelis
 
 Low-level functions for solving the single diode equation.
 
@@ -37,3 +38,4 @@ Functions for fitting diode models
     ivtools.sde.fit_sandia_simple
     ivtools.sdm.fit_cec_sam
     ivtools.sdm.fit_desoto
+    ivtools.sdm.fit_desoto_batzelis

docs/sphinx/source/reference/pv_modeling/system_models.rst

@@ -55,3 +55,11 @@ PVGIS model
    :toctree: ../generated/
 
     pvarray.huld
+
+Other
+^^^^^
+
+.. autosummary::
+   :toctree: ../generated/
+
+    pvarray.batzelis

docs/sphinx/source/whatsnew/v0.13.2.rst

@@ -21,6 +21,13 @@ Bug fixes
 
 Enhancements
 ~~~~~~~~~~~~
+* Add :py:func:`~pvlib.ivtools.sdm.fit_desoto_batzelis`, a function to estimate
+  parameters for the De Soto single-diode model from datasheet values. (:pull:`2563`)
+* Add :py:func:`~pvlib.singlediode.batzelis`, a function to estimate
+  maximum power, open circuit, and short circuit points using parameters for
+  the single-diode equation. (:pull:`2563`)
+* Add :py:func:`~pvlib.pvarray.batzelis`, a function to estimate maximum power
+  open circuit, and short circuit points from datasheet values. (:pull:`2563`)
 * Add ``method='chandrupatla'`` (faster than ``brentq`` and slower than ``newton``,
   but convergence is guaranteed) as an option for
   :py:func:`pvlib.pvsystem.singlediode`,

pvlib/ivtools/sdm/__init__.py

@@ -10,7 +10,8 @@
 
 from pvlib.ivtools.sdm.desoto import (  # noqa: F401
     fit_desoto,
-    fit_desoto_sandia
+    fit_desoto_batzelis,
+    fit_desoto_sandia,
 )
 
 from pvlib.ivtools.sdm.pvsyst import (  # noqa: F401

pvlib/ivtools/sdm/desoto.py

@@ -2,6 +2,7 @@
 
 from scipy import constants
 from scipy import optimize
+from scipy.special import lambertw
 
 from pvlib.ivtools.utils import rectify_iv_curve
 from pvlib.ivtools.sde import _fit_sandia_cocontent
@@ -399,3 +400,74 @@ def _fit_desoto_sandia_diode(ee, voc, vth, tc, specs, const):
     new_x = sm.add_constant(x)
     res = sm.RLM(y, new_x).fit()
     return np.array(res.params)[1]
+
+
+def fit_desoto_batzelis(v_mp, i_mp, v_oc, i_sc, alpha_sc, beta_voc):
+    """
+    Determine De Soto single-diode model parameters from datasheet values
+    using Batzelis's method.
+
+    This method is described in Section II.C of [1]_ and fully documented
+    in [2]_.
+
+    Parameters
+    ----------
+    v_mp : float
+        Maximum power point voltage at STC. [V]
+    i_mp : float
+        Maximum power point current at STC. [A]
+    v_oc : float
+        Open-circuit voltage at STC. [V]
+    i_sc : float
+        Short-circuit current at STC. [A]
+    alpha_sc : float
+        Short-circuit current temperature coefficient at STC. [A/K]
+    beta_voc : float
+        Open-circuit voltage temperature coefficient at STC. [V/K]
+
+    Returns
+    -------
+    dict
+        The returned dict contains the keys:
+
+        * ``alpha_sc`` [A/K]
+        * ``a_ref`` [V]
+        * ``I_L_ref`` [A]
+        * ``I_o_ref`` [A]
+        * ``R_sh_ref`` [Ohm]
+        * ``R_s`` [Ohm]
+
+    References
+    ----------
+    .. [1] E. I. Batzelis, "Simple PV Performance Equations Theoretically Well
+       Founded on the Single-Diode Model," Journal of Photovoltaics vol. 7,
+       no. 5, pp. 1400-1409, Sep 2017, :doi:`10.1109/JPHOTOV.2017.2711431`
+    .. [2] E. I. Batzelis and S. A. Papathanassiou, "A method for the
+       analytical extraction of the single-diode PV model parameters,"
+       IEEE Trans. Sustain. Energy, vol. 7, no. 2, pp. 504-512, Apr 2016.
+       :doi:`10.1109/TSTE.2015.2503435`
+    """
+    # convert temp coeffs from A/K and V/K to 1/K
+    alpha_sc = alpha_sc / i_sc
+    beta_voc = beta_voc / v_oc
+
+    # Equation numbers refer to [1]
+    t0 = 298.15  # K
+    del0 = (1 - beta_voc * t0) / (50.1 - alpha_sc * t0)  # Eq 9
+    w0 = np.real(lambertw(np.exp(1/del0 + 1)))
+
+    # Eqs 11-15
+    a0 = del0 * v_oc
+    Rs0 = (a0 * (w0 - 1) - v_mp) / i_mp
+    Rsh0 = a0 * (w0 - 1) / (i_sc * (1 - 1/w0) - i_mp)
+    Iph0 = (1 + Rs0 / Rsh0) * i_sc
+    Isat0 = Iph0 * np.exp(-1/del0)
+
+    return {
+        'alpha_sc': alpha_sc * i_sc,  # convert 1/K to A/K
+        'a_ref': a0,
+        'I_L_ref': Iph0,
+        'I_o_ref': Isat0,
+        'R_sh_ref': Rsh0,
+        'R_s': Rs0,
+    }

pvlib/pvarray.py

@@ -9,8 +9,9 @@
 """
 
 import numpy as np
+import pandas as pd
 from scipy.optimize import curve_fit
-from scipy.special import exp10
+from scipy.special import exp10, lambertw
 
 
 def pvefficiency_adr(effective_irradiance, temp_cell,
@@ -394,3 +395,131 @@ def huld(effective_irradiance, temp_mod, pdc0, k=None, cell_type=None,
         k[5] * tprime**2
     )
     return pdc
+
+
+def batzelis(effective_irradiance, temp_cell,
+             v_mp, i_mp, v_oc, i_sc, alpha_sc, beta_voc):
+    """
+    Compute maximum power point, open circuit, and short circuit
+    values using Batzelis's method.
+
+    Batzelis's method (described in Section III of [1]_) is a fast method
+    of computing the maximum power current and voltage.  The calculations
+    are rooted in the De Soto single-diode model, but require only typical
+    datasheet information.
+
+    Parameters
+    ----------
+    effective_irradiance : numeric, non-negative
+        Effective irradiance incident on the PV module. [Wm⁻²]
+    temp_cell : numeric
+        PV module operating temperature. [°C]
+    v_mp : float
+        Maximum power point voltage at STC. [V]
+    i_mp : float
+        Maximum power point current at STC. [A]
+    v_oc : float
+        Open-circuit voltage at STC. [V]
+    i_sc : float
+        Short-circuit current at STC. [A]
+    alpha_sc : float
+        Short-circuit current temperature coefficient at STC. [A/K]
+    beta_voc : float
+        Open-circuit voltage temperature coefficient at STC. [V/K]
+
+    Returns
+    -------
+    dict
+        The returned dict-like object contains the keys/columns:
+
+        * ``p_mp`` - power at maximum power point. [W]
+        * ``i_mp`` - current at maximum power point. [A]
+        * ``v_mp`` - voltage at maximum power point. [V]
+        * ``i_sc`` - short circuit current. [A]
+        * ``v_oc`` - open circuit voltage. [V]
+
+    Notes
+    -----
+    This method is the combination of three sub-methods for:
+
+    1. estimating single-diode model parameters from datasheet information
+    2. translating SDM parameters from STC to operating conditions
+       (taken from the De Soto model)
+    3. estimating the MPP, OC, and SC points on the resulting I-V curve.
+
+    At extremely low irradiance (e.g. 1e-10 Wm⁻²), this model can produce
+    negative voltages.  This function clips any negative voltages to zero.
+
+    References
+    ----------
+    .. [1] E. I. Batzelis, "Simple PV Performance Equations Theoretically Well
+       Founded on the Single-Diode Model," Journal of Photovoltaics vol. 7,
+       no. 5, pp. 1400-1409, Sep 2017, :doi:`10.1109/JPHOTOV.2017.2711431`
+
+    Examples
+    --------
+    >>> params = {'i_sc': 15.98, 'v_oc': 50.26, 'i_mp': 15.27, 'v_mp': 42.57,
+    ...           'alpha_sc': 0.007351, 'beta_voc': -0.120624}
+    >>> batzelis(np.array([1000, 800]), np.array([25, 30]), **params)
+    {'p_mp': array([650.0439    , 512.99199048]),
+     'i_mp': array([15.27      , 12.23049303]),
+     'v_mp': array([42.57      , 41.94368856]),
+     'i_sc': array([15.98    , 12.813404]),
+     'v_oc': array([50.26      , 49.26532902])}
+    """
+    # convert temp coeffs from A/K and V/K to 1/K
+    alpha_sc = alpha_sc / i_sc
+    beta_voc = beta_voc / v_oc
+
+    t0 = 298.15
+    delT = temp_cell - (t0 - 273.15)
+    lamT = (temp_cell + 273.15) / t0
+    g = effective_irradiance / 1000
+    # for zero/negative irradiance, use lnG=large negative number so that
+    # computed voltages are negative and then clipped to zero
+    with np.errstate(divide='ignore'):  # needed for pandas for some reason
+        lnG = np.log(g, out=np.full_like(g, -9e9), where=(g > 0))
+        lnG = np.where(np.isfinite(g), lnG, np.nan)  # also preserve nans
+
+    # Eq 9-10
+    del0 = (1 - beta_voc * t0) / (50.1 - alpha_sc * t0)
+    w0 = np.real(lambertw(np.exp(1/del0 + 1)))
+
+    # Eqs 27-28
+    alpha_imp = alpha_sc + (beta_voc - 1/t0) / (w0 - 1)
+    beta_vmp = (v_oc / v_mp) * (
+        beta_voc / (1 + del0) +
+        (del0 * (w0 - 1) - 1/(1 + del0)) / t0
+    )
+
+    # Eq 26
+    eps0 = (del0 / (1 + del0)) * (v_oc / v_mp)
+    eps1 = del0 * (w0 - 1) * (v_oc / v_mp) - 1
+
+    # Eqs 22-25
+    isc = g * i_sc * (1 + alpha_sc * delT)
+    voc = v_oc * (1 + del0 * lamT * lnG + beta_voc * delT)
+    imp = g * i_mp * (1 + alpha_imp * delT)
+    vmp = v_mp * (1 + eps0 * lamT * lnG + eps1 * (1 - g) + beta_vmp * delT)
+
+    # handle negative voltages from zero and extremely small irradiance
+    vmp = np.clip(vmp, a_min=0, a_max=None)
+    voc = np.clip(voc, a_min=0, a_max=None)
+
+    out = {
+        'p_mp': vmp * imp,
+        'i_mp': imp,
+        'v_mp': vmp,
+        'i_sc': isc,
+        'v_oc': voc,
+    }
+
+    # if pandas in, ensure pandas out
+    pandas_inputs = [
+        x for x in [effective_irradiance, temp_cell]
+        if isinstance(x, pd.Series)
+    ]
+    if pandas_inputs:
+        out = pd.DataFrame(out, index=pandas_inputs[0].index)
+
+    return out

pvlib/singlediode.py

@@ -3,6 +3,7 @@
 """
 
 import numpy as np
+import pandas as pd
 from pvlib.tools import _golden_sect_DataFrame
 
 from scipy.optimize import brentq, newton
@@ -975,3 +976,85 @@ def _pwr_optfcn(df, loc):
                                  df['resistance_shunt'], df['nNsVth'])
 
     return current * df[loc]
+
+
+def batzelis(photocurrent, saturation_current, resistance_series,
+             resistance_shunt, nNsVth):
+    """
+    Estimate maximum power, open-circuit, and short-circuit points from
+    single-diode equation parameters using Batzelis's method.
+
+    This method is described in Section II.B of [1]_.
+
+    Parameters
+    ----------
+    photocurrent : numeric
+        Light-generated current. [A]
+    saturation_current : numeric
+        Diode saturation current. [A]
+    resistance_series : numeric
+        Series resistance. [Ohm]
+    resistance_shunt : numeric
+        Shunt resistance. [Ohm]
+    nNsVth : numeric
+        The product of the usual diode ideality factor (n, unitless),
+        number of cells in series (Ns), and cell thermal voltage at
+        specified effective irradiance and cell temperature. [V]
+
+    Returns
+    -------
+    dict
+        The returned dict-like object contains the keys/columns:
+
+        * ``p_mp`` - power at maximum power point. [W]
+        * ``i_mp`` - current at maximum power point. [A]
+        * ``v_mp`` - voltage at maximum power point. [V]
+        * ``i_sc`` - short circuit current. [A]
+        * ``v_oc`` - open circuit voltage. [V]
+
+    References
+    ----------
+    .. [1] E. I. Batzelis, "Simple PV Performance Equations Theoretically Well
+       Founded on the Single-Diode Model," Journal of Photovoltaics vol. 7,
+       no. 5, pp. 1400-1409, Sep 2017, :doi:`10.1109/JPHOTOV.2017.2711431`
+    """
+    # convenience variables
+    Iph = photocurrent
+    Is = saturation_current
+    Rsh = resistance_shunt
+    Rs = resistance_series
+    a = nNsVth
+
+    # Eqs 3-4
+    isc = Iph / (Rs / Rsh + 1)  # manipulated to handle Rsh=np.inf correctly
+    with np.errstate(divide='ignore'):  # zero Iph
+        voc = a * np.log(Iph / Is)
+
+    # Eqs 5-8
+    w = np.real(lambertw(np.e * Iph / Is))
+    # vmp = (1 + Rs/Rsh) * a * (w - 1) - Rs * Iph * (1 - 1/w)  # not needed
+    with np.errstate(divide='ignore', invalid='ignore'):  # zero Iph -> zero w
+        imp = Iph * (1 - 1/w) - a * (w - 1) / Rsh
+
+    vmp = a * (w - 1) - Rs * imp
+
+    vmp = np.where(Iph > 0, vmp, 0)
+    voc = np.where(Iph > 0, voc, 0)
+    imp = np.where(Iph > 0, imp, 0)
+    isc = np.where(Iph > 0, isc, 0)
+
+    out = {
+        'p_mp': imp * vmp,
+        'i_mp': imp,
+        'v_mp': vmp,
+        'i_sc': isc,
+        'v_oc': voc,
+    }
+
+    # if pandas in, ensure pandas out
+    pandas_inputs = [
+        x for x in [Iph, Is, Rsh, Rs, a] if isinstance(x, pd.Series)]
+    if pandas_inputs:
+        out = pd.DataFrame(out, index=pandas_inputs[0].index)
+
+    return out

tests/ivtools/sdm/test_desoto.py

@@ -92,3 +92,27 @@ def test_fit_desoto_sandia(cec_params_cansol_cs5p_220p):
     assert_allclose(result['dEgdT'], -0.0002677)
     assert_allclose(result['EgRef'], 1.3112547292120638)
     assert_allclose(result['cells_in_series'], specs['cells_in_series'])
+
+
+def test_fit_desoto_batzelis():
+    params = {'i_sc': 15.98, 'v_oc': 50.26, 'i_mp': 15.27, 'v_mp': 42.57,
+              'alpha_sc': 0.007351, 'beta_voc': -0.120624}
+    expected = {  # calculated with the function itself
+        'alpha_sc': 0.007351,
+        'a_ref': 1.7257632483733483,
+        'I_L_ref': 15.985408866796396,
+        'I_o_ref': 3.594308384705643e-12,
+        'R_sh_ref': 389.4379947026243,
+        'R_s': 0.13181590981241956,
+    }
+    out = sdm.fit_desoto_batzelis(**params)
+    for k in expected:
+        assert out[k] == pytest.approx(expected[k])
+
+    # ensure the STC values are reproduced
+    iv = pvsystem.singlediode(out['I_L_ref'], out['I_o_ref'], out['R_s'],
+                              out['R_sh_ref'], out['a_ref'])
+    assert iv['i_sc'] == pytest.approx(params['i_sc'])
+    assert iv['i_mp'] == pytest.approx(params['i_mp'], rel=3e-3)
+    assert iv['v_oc'] == pytest.approx(params['v_oc'], rel=3e-4)
+    assert iv['v_mp'] == pytest.approx(params['v_mp'], rel=4e-3)