add tolerance to tools._golden_sect_Dataframe

[cwhanse]

• Closes Apparent numerical instability in I_mp calculation using PVsyst model #1087
• I am familiar with the contributing guidelines
• Tests added
• [ ] Updates entries to docs/sphinx/source/api.rst for API changes.
• Adds description and name entries in the appropriate "what's new" file in docs/sphinx/source/whatsnew for all changes. Includes link to the GitHub Issue with :issue:`num` or this Pull Request with :pull:`num`. Includes contributor name and/or GitHub username (link with :ghuser:`user`).
• New code is fully documented. Includes numpydoc compliant docstrings, examples, and comments where necessary.
• Pull request is nearly complete and ready for detailed review.
• Maintainer: Appropriate GitHub Labels and Milestone are assigned to the Pull Request and linked Issue.

Sticking with tools._golden_sect_DataFrame because it vectorizes single-variable optimization, which looks clumsy (if possible at all) with existing scipy tools. Absolute tolerance was hardwired at 0.01V. Setting a default absolute tolerance of 1e-8V takes care of the reported instability.

[wholmgren]

Sticking with tools._golden_sect_DataFrame because it vectorizes single-variable optimization, which looks clumsy (if possible at all) with existing scipy tools.

If I recall correctly, @mikofski contributed a vectorized newton optimizer to scipy in large part so that we could remove this function.

[cwhanse]

He did, and it's in scipy now, thanks for that reminder. I'll have to think a bit about ways to get an initial estimate of Vmp that minimizes risk of divergence. That's why the bisection algorithm was originally in PVLib for Matlab, which led to the golden section function.

It also appears that at least some of the previous "known" solutions for single diode equations using the 'lambertw' method were not converged.

This PR will take more work than I initially thought.

[mikofski]

I added 2 SciPy methods, the vectorized newton method which uses NumPy, and a cythonized api to the bounded search methods, with examples to call them with cythonize prange and nogil so they're two different things. Altho easier, the newton method would be useful here, instead I'd recommend the cythonized BrentQ

[mikofski]

cython_optimize_brentq_benchmark.ipynb

[mikofski]

Toms 748 is another bounded search algorithm but idk if it's vectorized

[cwhanse]

This is ready for review.

I recalculated the IV curves used for several tests to add precision. I believe the old expected results were done by me, in Matlab, before we improved convergence in the PV_Lib Matlab code.

[wholmgren]

before we improved convergence in the PV_Lib Matlab code.

I'm a little confused... the linked PR adds a warning to the user that iteration limit was exceeded, but you're removing it in this PR. Maybe you meant to point to a different PR? Still wondering if removing that exception is a good idea though. Perhaps a warning would be better than an exception and maybe more consistent with other solvers.

[cwhanse]

I'm a little confused

Understandable. That work in PV_LIB for Matlab wasn't done via PRs, now that I look again at the one I linked.

Theoretically, the exit criterion (err < atol) will be satisfied before interations reaches the computed iterlimit. I suppose some very pathological func could be engineered to make this fail, but I couldn't find one so I remove the exception to satisfy codecov. I'm OK putting it back, either as an exception or warning, but that line won't be hit by a test.

[wholmgren]

but that line won't be hit by a test.

You can add # pragma: no cover to prevent the coverage tools from complaining.

I'm ok moving forward without the check if others are too.

[kandersolar]

I'm ok moving forward without the check if others are too.

I'm not familiar enough with what's happening here to vote. I don't object, anyway.

[wholmgren]

Shall we merge?

I am wondering if there's a meaningful difference in execution time. If so, might be worth noting in the what's new.

[cwhanse]

For the golden section function, the old default tolerance was 1e-2. With the new default of 1e-8, convergence takes 3x longer when applied to the test function. That's probably typical of what will be seen when solving single diode equations. I don't think the increase will be visible in a ModelChain simulation (most of the execution time is in the solar position algorithm).

from pvlib.tools import _golden_sect_DataFrame
import timeit

def _obj_test_golden_sect(params, loc):
    return params[loc] * (1. - params['c'] * params[loc]**params['n'])


def test__golden_sect_DataFrame_atol(atol):
    params = {'c': 0.2, 'n': 0.3}
    expected = 89.14332727531685
    v, x = _golden_sect_DataFrame(
        params, 0., 100., _obj_test_golden_sect, atol=atol)


t8 = timeit.timeit('test__golden_sect_DataFrame_atol(1e-8)', number=10000, globals=globals())
t2 = timeit.timeit('test__golden_sect_DataFrame_atol(1e-2)', number=10000, globals=globals())

t8 = 0.064 ms, and t2= 0.025 ms, per execution.

[kandersolar]

@cwhanse I think this needs to use nanmax to correctly handle nighttime. iterlimit is ending up as nan in the docs build in #1394.

[cwhanse]

I agree. But handling nan in the bounds and where returned by func is a bit more involved. Currently, I think the optimization would run to its iteration limit.