singlediode error with very low effective_irradiance

[pasquierjb]

Describe the bug

Since pvlib 0.9.4 release (#1606) I get an error while running the single-diode model with some very low effective irradiance values.

To Reproduce

from pvlib import pvsystem

effective_irradiance=1.341083e-17
temp_cell=13.7 

cec_modules = pvsystem.retrieve_sam('CECMod')
cec_module = cec_modules['Trina_Solar_TSM_300DEG5C_07_II_']

mount = pvsystem.FixedMount()
array = pvsystem.Array(mount=mount,
                       module_parameters=cec_module)

system = pvsystem.PVSystem(arrays=[array])

params = system.calcparams_cec(effective_irradiance, 
                               temp_cell)

system.singlediode(*params)

    303 """
    304 Vectorized golden section search for finding maximum of a function of a
    305 single variable.
   (...)
    342 pvlib.singlediode._pwr_optfcn
    343 """
    344 if np.any(upper - lower < 0.):
--> 345     raise ValueError('upper >= lower is required')
    347 phim1 = (np.sqrt(5) - 1) / 2
    349 df = params

ValueError: upper >= lower is required

Expected behavior
This complicates the bifacial modeling procedure as run_model_from_effective_irradiance can be called with very low irradiance values estimated by pvfactors (at sunrise or sunset for instance).

Versions:

• pvlib.__version__: 0.9.4
• pandas.__version__: 1.5.3
• python: 3.10

Additional context

v_oc is negative in this case which causes the error.

from pvlib.singlediode import _lambertw_v_from_i
photocurrent = params[0]
saturation_current = params[1]
resistance_series = params[2]
resistance_shunt = params[3]
nNsVth = params[4]
v_oc = _lambertw_v_from_i(resistance_shunt, resistance_series, nNsVth, 0.,
                              saturation_current, photocurrent)

[cwhanse]

That's unexpected, thanks for reporting.

I'll note that the negative Voc results from taking the difference of two very large but nearly equal numbers. It's likely limited to the CEC model, where the shunt resistance is inversely proportional to irradiance, which would be about 1e19 at photocurrent of 1e-17 for this case.

[cwhanse]

Now this gets strange: the Voc value is positive with pvlib v0.9.3. The function involved pvlib.singlediode._lambertw_v_from_i hasn't changed for many releases. In both pvlib v0.9.3 and v0.9.4, in this calculation of Voc, the lambertw term overflows so the Voc value is computed using only python arithmetic operators and numpy.log.

I'm starting to think the error depends on python and numpy versions.

[kandersolar]

The difference between 0.9.3 and 0.9.4 here may be due to slightly different values returned by calcparams_cec. Compare the output of print(list(map(str, params))); I get slightly different saturation current values for the given example. Maybe the changed Boltzmann constant in #1617 is the cause?

[cwhanse]

+1 to #1617 as the likely culprit. I get the positive/negative Voc values with the same python and numpy versions but different pvlib versions.

[cwhanse]

To illustrate the challenge, this line computes the Voc.

Stripping out the indexing the computation is

    V = (IL + I0 - I) / Gsh - \
        I * Rs - a * lambertwterm

With pvlib v0.9.4, Io is 7.145289906185543e-12. a is not affected, since a value of the Boltzmann contant is inherent in the a_ref value from the database. (IL + I0 - I) / Gsh is 107825185636.40567, I * Rs is 0, and a * lambertwterm is 107825185636.40569

With pvlib v0.9.3, Io is 7.145288699667595e-12. (IL + I0 - I) / Gsh is 107825167429.58397, I * Rs is 0, and a * lambertwterm is 107825167429.58395

The difference defining Voc is in the least significant digit.

[cwhanse]

Increasing the iterations that solve for lambertwterm doesn't fix this issue.

[kandersolar]

This smells to me like the inevitable error from accumulated round-off.

FWIW, negative Voc can be achieved in 0.9.3 as well -- try the given example but with effective_irradiance=1.e-18. The difference is that before #1606, it led to nans and warnings instead of raising an error.

[cwhanse]

@pasquierjb I recommend intercepting the effective irradiance and setting values to 0 which are below a minimum on the order of 1e-9 W/m2. That will propagate to shunt resistance = np.inf, which changes the calculation path in pvlib.singlediode and gives Voc=0.

I'm not sure we'll be able to extend the numerical solution of the single diode equation to be accurate at very low but non-zero values of photocurrent (and/or enormous but finite values of shunt resistance.)

I note that pvlib.pvsystem.calcparams_desoto doesn't like effective_irradiance=0. but is OK with effective_irradiance=np.array([0.]). Has to do with trapping and ignoring division by zero warnings and errors.

[mikofski]

Have you tried setting method='newton' instead of 'lambertw'? https://pvlib-python.readthedocs.io/en/stable/reference/generated/pvlib.pvsystem.singlediode.html#pvlib-pvsystem-singlediode

[cwhanse]

Setting method='newton' gets a solution to this case. method isn't available as a parameter of the PVSystem.singlediode method so @pasquierjb would need to change his workflow to use it. Something for us to consider adding.

[SonGorku]

My workaround for this issue was to first filter very low effective_irradiance values (<1e-8), and then filter photocurrent and saturation_current parameters when effective_irradiance=0 and made them =0. This assures that you won't get negative v_oc values.