Modify initial Voc for bishop88 functions

pvlib/singlediode.py

@@ -304,6 +304,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
@@ -317,7 +320,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)
@@ -431,6 +434,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"
@@ -449,10 +455,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)
@@ -561,6 +567,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]
@@ -574,12 +583,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,39 @@ 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):
+    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)
+    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()