Bishop88 functions sometimes return incorrect results or errors when using d2mutau and NsVbi parameters

[kandersolar]

Describe the bug
The bishop88 functions can sometimes fail, or return incorrect results, when using the d2mutau and NsVbi parameters.

To Reproduce

Set-up:

import pvlib
import numpy as np
import matplotlib.pyplot as plt

sdm_params = {
    'alpha_sc': 0.001, 'gamma_ref': 1.409, 'mu_gamma': 0.001, 'I_o_ref': 3.96e-11,
    'R_sh_ref': 4900, 'R_sh_0': 8800, 'R_sh_exp': 5.5,
    'R_s': 6.76, 'cells_in_series': 264, 'EgRef': 1.121, 'I_L_ref': 2.278883384981941
}
d2mutau = 1.5
NsVbi = 237.6

G = np.linspace(0, 1200)
T = 25
sde_params = pvlib.pvsystem.calcparams_pvsyst(G, T, **sdm_params)

method='brentq' errors with ValueError: f(a) and f(b) must have different signs.

Click to show traceback

Traceback (most recent call last):

  Cell In[16], line 17
    mpp = pvlib.singlediode.bishop88_mpp(*sde_params, d2mutau=d2mutau, NsVbi=NsVbi, method='brentq')

  File ~\software\miniconda3\envs\dev\lib\site-packages\pvlib\singlediode.py:577 in bishop88_mpp
    vd = vec_fun(voc_est, *args)

  File ~\software\miniconda3\envs\dev\lib\site-packages\numpy\lib\function_base.py:2372 in __call__
    return self._call_as_normal(*args, **kwargs)

  File ~\software\miniconda3\envs\dev\lib\site-packages\numpy\lib\function_base.py:2365 in _call_as_normal
    return self._vectorize_call(func=func, args=vargs)

  File ~\software\miniconda3\envs\dev\lib\site-packages\numpy\lib\function_base.py:2450 in _vectorize_call
    ufunc, otypes = self._get_ufunc_and_otypes(func=func, args=args)

  File ~\software\miniconda3\envs\dev\lib\site-packages\numpy\lib\function_base.py:2410 in _get_ufunc_and_otypes
    outputs = func(*inputs)

  File ~\software\miniconda3\envs\dev\lib\site-packages\pvlib\singlediode.py:572 in <lambda>
    vbr_exp: brentq(fmpp, 0.0, voc,

  File ~\software\miniconda3\envs\dev\lib\site-packages\scipy\optimize\_zeros_py.py:809 in brentq
    r = _zeros._brentq(f, a, b, xtol, rtol, maxiter, args, full_output, disp)

ValueError: f(a) and f(b) must have different signs

method='newton' gives incorrect results for some conditions:

mpp_newton = pvlib.singlediode.bishop88_mpp(*sde_params, d2mutau=d2mutau, NsVbi=NsVbi, method='newton')
plt.plot(G, mpp_newton[2])
plt.ylabel('$P_{mp}$ [W]')
plt.xlabel('Irradiance [W/m2]')

Versions:

• pvlib.__version__: 0.10.4
• scipy.__version__: 1.12.0

Additional context
I think the problem has something to do with the following:

1. The SDE with the thin film recombination current term has a singularity at NsVbi = Vd, and our code only returns valid results for Vd < NsVbi (ref Adding New Single Diode Model For CdTe #163 (comment))
2. estimate_voc can sometimes return a value on the other side of the singularity (Vd > NsVbi), where our bishop88-style calculations are not valid.
3. bishop88_mpp uses estimate_voc to calculate either an upper bound for the bracketed optimization (with method='brentq') or the initial guess (with method='newton'), meaning we are sometimes providing invalid starting points to the root finder.

Here is an example of the estimated Voc being greater than NsVbi:

import pvlib

sdm_params = {
    'alpha_sc': 0.001, 'gamma_ref': 1.409, 'mu_gamma': 0.001, 'I_o_ref': 3.96e-11,
    'R_sh_ref': 4900, 'R_sh_0': 8800, 'R_sh_exp': 5.5,
    'R_s': 6.76, 'cells_in_series': 264, 'EgRef': 1.121, 'I_L_ref': 2.278883384981941
}
d2mutau = 1.5
NsVbi = 237.6

G = 1100
T = 25

sde_params = pvlib.pvsystem.calcparams_pvsyst(G, T, **sdm_params)
voc_est = pvlib.singlediode.estimate_voc(sde_params[0], sde_params[1], sde_params[-1])
print(voc_est)  # 237.6945087345138, greater than NsVbi=237.6

cc @leliadeville

[smmeredith]

Thank you for bringing this issue up. This is something I've noticed recently as well. I have tested a fair amount of CdTe modules with bishop88 and have seen "discontinuities" like the one you pointed out above - in the form of extremely high or low MPPs, or just getting nan values returned at those points.

A minor question regarding the example you gave - was that module supposed to CdTe? If so, shouldn't the energy bandgap be 1.5eV?

[kandersolar]

was that module supposed to CdTe? If so, shouldn't the energy bandgap be 1.5eV?

Yes, good point. I guess my brain and fingers were on autopilot when putting together an example set of parameters. Interestingly, updating the parameters to use 1.5 eV doesn't seem to change the simulation results at all... am I overlooking something?

Also @smmeredith, if you can, it would be great if you would assemble a collection of parameter sets that you have observed to produce incorrect simulations. Whenever someone fixes this bug, the fix could be evaluated on that collection.

[smmeredith]

updating the parameters to use 1.5 eV doesn't seem to change the simulation results at all

I think this is because only saturation current uses it, which has no irradiance dependence.

it would be great if you would assemble a collection of parameter sets that you have observed to produce incorrect simulations

No problem, just let me know how I can share it.

So I've done a little more investigating on the 25 CdTe modules that I've been testing with bishop88, and I think you're right about the problem being related to the Voc estimate. For each module, I:

1. Calculated the estimated Voc using singlediode.estimate_voc
2. Ran the module through the singlediode.bishop_mpp model (using the newton method) at a range of temperatures from -25 to +74 C in steps of 1 C, and a range of irradiances from 20 to 1100 W/m2 in steps of 20 W/m2.
3. Calculated the efficiency based on the output p_mp and identified problematic outputs where the efficiency was either nan, less than zero, or greater than 30%.
4. Calculated NsVbi - Voc

For every module that I looked at, whenever a problematic output occurred, NsVbi - Voc was negative. This is not to say that ALL negative NsVbi - Voc caused a problematic output - in fact there were many outputs where this expression was negative and the output seemed reasonable.

[smmeredith]

Attaching a JSON of parameters that were used in my analysis.

I removed seven parameter sets that came from a third-party lab, because my single-diode model does not match PVsyst as well for those modules - and that is a separate issue that needs to be resolved by me.

sd_params.json

[cwhanse]

Modifying bishop88_mpp here from

        x0, args, method_kwargs = \
            _prepare_newton_inputs(voc_est, args, method_kwargs)

to

        xp = np.where(voc_est < NsVbi, voc_est, 0.99*NsVbi)
        x0, args, method_kwargs = \
            _prepare_newton_inputs(xp, args, method_kwargs)

fixes all np.isnan(p_mp) | (p_mp < 0) in the modeling cases described by @smmeredith (irradiance 20W - 1100W by 20W, cell temp. -25C to +74C, parameters from the posted sd_params.json file). I don't have the module areas so can't tell if the high efficiency values are also fixed, but I suspect yes.

We want to start the Newton iteration a little less than NsVbi, where the derivative dP/dV should be negative (and large in magnitude), and then I think Newton's should be reliable (as reliable as Newton's can be). Illustration of the power vs. diode voltage curve in the neighborhood of the asymptote at Vd=NsVbi, we're trying to find the max power point (red square)

.

[kandersolar]

@cwhanse how about with method='brentq'?

[cwhanse]

@cwhanse how about with method='brentq'?

Same change to the right end point also fixes the failures using method='brentq'.

Edit: ...and prevents a number of warnings that are currently emitted.