Modify initial Voc for bishop88 functions

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

@@ -15,6 +15,10 @@ Enhancements
 
 Bug fixes
 ~~~~~~~~~
+* Improved reliability of :py:func:`pvlib.singlediode.bishop88_mpp`,
+  :py:func:`pvlib.singlediode.bishop88_i_from_v` and
+  :py:func:`pvlib.singlediode.bishop88_v_from_i` by improving the initial
+  guess for the newton and brentq algorithms. (:issue:`2013`, :pull:`2032`)
 * Corrected equation for Ixx0 in :py:func:`pvlib.pvsystem.sapm` (:issue:`2016`, :pull:`2019`)
 * Fixed :py:func:`pvlib.pvsystem.retrieve_sam` silently ignoring the `path` parameter
   when `name` was provided. Now an exception is raised requesting to only provide one

pvlib/singlediode.py

@@ -310,6 +310,9 @@ def fv(x, v, *a):
     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
@@ -323,7 +326,7 @@ def vd_from_brent(voc, v, iph, isat, rs, rsh, gamma, d2mutau, NsVbi,
                           **method_kwargs)
 
         vd_from_brent_vectorized = np.vectorize(vd_from_brent)
-        vd = vd_from_brent_vectorized(voc_est, voltage, *args)
+        vd = vd_from_brent_vectorized(xp, voltage, *args)
     elif method == 'newton':
         x0, (voltage, *args), method_kwargs = \
             _prepare_newton_inputs(voltage, (voltage, *args), method_kwargs)
@@ -443,6 +446,9 @@ def bishop88_v_from_i(current, photocurrent, saturation_current,
 
     # 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"
@@ -461,10 +467,10 @@ def vd_from_brent(voc, i, iph, isat, rs, rsh, gamma, d2mutau, NsVbi,
                           **method_kwargs)
 
         vd_from_brent_vectorized = np.vectorize(vd_from_brent)
-        vd = vd_from_brent_vectorized(voc_est, current, *args)
+        vd = vd_from_brent_vectorized(xp, current, *args)
     elif method == 'newton':
         x0, (current, *args), method_kwargs = \
-            _prepare_newton_inputs(voc_est, (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)
@@ -579,6 +585,9 @@ def bishop88_mpp(photocurrent, saturation_current, resistance_series,
 
     # 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]
@@ -592,12 +601,13 @@ def fmpp(x, *a):
                                   vbr_a, vbr, vbr_exp),
                             **method_kwargs)
         )
-        vd = vec_fun(voc_est, *args)
+        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(voc_est, 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)

pvlib/tests/test_singlediode.py

@@ -575,3 +575,53 @@ def test_bishop88_pdSeries_len_one(method, bishop88_arguments):
     bishop88_i_from_v(pd.Series([0]), **bishop88_arguments, method=method)
     bishop88_v_from_i(pd.Series([0]), **bishop88_arguments, method=method)
     bishop88_mpp(**bishop88_arguments, method=method)
+
+
+def _sde_check_solution(i, v, il, io, rs, rsh, a, d2mutau=0., NsVbi=np.inf):
+    vd = v + rs * i
+    return il - io*np.expm1(vd/a) - vd/rsh - il*d2mutau/(NsVbi - vd) - i
+
+
+@pytest.mark.parametrize('method', ['newton', 'brentq'])
+def test_bishop88_init_cond(method):
+    # GH 2013
+    p = {'alpha_sc': 0.0012256,
+         'gamma_ref': 1.2916241612804187,
+         'mu_gamma': 0.00047308959960937403,
+         'I_L_ref': 3.068717040806731,
+         'I_o_ref': 2.2691248021217617e-11,
+         'R_sh_ref': 7000,
+         'R_sh_0': 7000,
+         'R_s': 4.602,
+         'cells_in_series': 268,
+         'R_sh_exp': 5.5,
+         'EgRef': 1.5}
+    NsVbi = 268 * 0.9
+    d2mutau = 1.4
+    irrad = np.arange(20, 1100, 20)
+    tc = np.arange(-25, 74, 1)
+    weather = np.array(np.meshgrid(irrad, tc)).T.reshape(-1, 2)
+    # with the above parameters and weather conditions, a few combinations
+    # result in voc_est > NsVbi, which causes failure of brentq and newton
+    # when the recombination parameters NsVbi and d2mutau are used.
+    sde_params = pvsystem.calcparams_pvsyst(weather[:, 0], weather[:, 1], **p)
+    # test _mpp
+    result = bishop88_mpp(*sde_params, d2mutau=d2mutau, NsVbi=NsVbi)
+    imp, vmp, pmp = result
+    err = np.abs(_sde_check_solution(
+        imp, vmp, sde_params[0], sde_params[1], sde_params[2], sde_params[3],
+        sde_params[4], d2mutau=d2mutau, NsVbi=NsVbi))
+    bad_results = np.isnan(pmp) | (pmp < 0) | (err > 0.00001)  # 0.01mA error
+    assert not bad_results.any()
+    # test v_from_i
+    vmp2 = bishop88_v_from_i(imp, *sde_params, d2mutau=d2mutau, NsVbi=NsVbi)
+    err = np.abs(_sde_check_solution(imp, vmp2, *sde_params, d2mutau=d2mutau,
+                                     NsVbi=NsVbi))
+    bad_results = np.isnan(vmp2) | (vmp2 < 0) | (err > 0.00001)
+    assert not bad_results.any()
+    # test v_from_i
+    imp2 = bishop88_i_from_v(vmp, *sde_params, d2mutau=d2mutau, NsVbi=NsVbi)
+    err = np.abs(_sde_check_solution(imp2, vmp, *sde_params, d2mutau=d2mutau,
+                                     NsVbi=NsVbi))
+    bad_results = np.isnan(imp2) | (imp2 < 0) | (err > 0.00001)
+    assert not bad_results.any()