Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

ENH: add methods and tests for a explicit IV curve calculation of single-diode model #409

Merged
merged 76 commits into from
Jul 10, 2018

Conversation

mikofski
Copy link
Member

@mikofski mikofski commented Jan 27, 2018

pvlib python pull request guidelines

Thank you for your contribution to pvlib python!

You may submit a pull request with your code at any stage of completion, however, before the code can be merged the following items must be addressed:

  • Closes singlediode is slow #408 and closes Single diode model helper function to wrap Lambert w and other methods #410
  • Fully tested. Added and/or modified tests to ensure correct behavior for all reasonable inputs. Tests must pass on the TravisCI and Appveyor testing services.
  • Code quality and style is sufficient. Passes git diff upstream/master -u -- "*.py" | flake8 --diff and/or landscape.io linting service.
  • New code is fully documented. Includes sphinx/numpydoc compliant docstrings and comments in the code where necessary.
  • Updates entries to docs/sphinx/source/api.rst for API changes.
  • Adds description and name entries in the appropriate docs/sphinx/source/whatsnew file for all changes.

Please don't hesitate to ask for help if you're unsure of how to accomplish any of the above. You may delete all of these instructions except for the list above.

Brief description of the problem and proposed solution (if not already fully described in the issue linked to above):

Hi all,

I implemented a fast version of the single diode using a damped newton step with pre-calculated derivatives. I think that is probably the main reason for the speedup. This is still a WIP. I still have to implement the IX, IXX, and the number of points, but I'm confident we will still see the same level of speedup even after those calculations. For the number of points, I'll use the explicit method that Bishop used in the 1988 paper, since it's way faster. Also, I've started a generic newton_step method, but I haven't substituted in where I did a lot of copy & paste (aka: WET) code.

Also, I found that the v_mp, i_mp, and p_mp values were consistently slightly off in the existing method, so in the assertion tests I had to use i_mp = i_from_v(vmp) and v_mp = v_from_i(imp) where imp and vmp are the new faster values. I think this is because the faster method finds a slightly higher pmp value. Of course you can adjust the tolerance, to get even faster convergence, if you're okay with slightly inaccurate values. The newton step iterates about 4 times. The default tolerance is the min EPS of the system, a bit ridiculous, so perhaps lowering it to 1e-5 is better ?

Also note that this method so far, does not depend on SciPy at all, so perhaps it can pass the minimum python travis ci tests?

There is a logging function and a gradient check that we can remove in the final version assuming this is approved. Enjoy!

If you run test_way_faster.py you can get a flair for the speedup:

In [1]: %run pvlib/test/test_way_faster.py
DEBUG:__main__:single diode elapsed time = 0.108363[s]
DEBUG:__main__:way faster elapsed time = 0.000411[s]
DEBUG:__main__:spr_e20_327 speedup = 263.657
DEBUG:__main__:single diode elapsed time = 0.007208[s]
DEBUG:__main__:way faster elapsed time = 0.000314[s]
DEBUG:__main__:fs_495 speedup = 22.9554
In [2]: %run pvlib/test/test_way_faster.py
DEBUG:__main__:single diode elapsed time = 0.00827[s]
DEBUG:__main__:way faster elapsed time = 0.000384[s]
DEBUG:__main__:spr_e20_327 speedup = 21.5365
DEBUG:__main__:single diode elapsed time = 0.008368[s]
DEBUG:__main__:way faster elapsed time = 0.000252[s]
DEBUG:__main__:fs_495 speedup = 33.2063
In [3]: %run pvlib/test/test_way_faster.py
DEBUG:__main__:single diode elapsed time = 0.006553[s]
DEBUG:__main__:way faster elapsed time = 0.000386[s]
DEBUG:__main__:spr_e20_327 speedup = 16.9767
DEBUG:__main__:single diode elapsed time = 0.00638[s]
DEBUG:__main__:way faster elapsed time = 0.000386[s]
DEBUG:__main__:fs_495 speedup = 16.5285
In [4]: %run pvlib/test/test_way_faster.py
DEBUG:__main__:single diode elapsed time = 0.006512[s]
DEBUG:__main__:way faster elapsed time = 0.000348[s]
DEBUG:__main__:spr_e20_327 speedup = 18.7126
DEBUG:__main__:single diode elapsed time = 0.010524[s]
DEBUG:__main__:way faster elapsed time = 0.000345[s]
DEBUG:__main__:fs_495 speedup = 30.5043
In [5]: %run pvlib/test/test_way_faster.py
DEBUG:__main__:single diode elapsed time = 0.00794[s]
DEBUG:__main__:way faster elapsed time = 0.0003[s]
DEBUG:__main__:spr_e20_327 speedup = 26.4667
DEBUG:__main__:single diode elapsed time = 0.007225[s]
DEBUG:__main__:way faster elapsed time = 0.000256[s]
DEBUG:__main__:fs_495 speedup = 28.2227

@mikofski
Copy link
Member Author

Adding the i-v curve doesn't affect the speedup at all, and the curve output visually looks very similar to the existing singlediode output, except that it extends slightly to the 2nd quadrant because Vdiode is always less than Vcell, so although at Voc they are the same, at Isc, they are different. Here's the plot with singlediode as blue dashes, sorry I didn't add a legend.

faster_way-v-singlediode

@markcampanelli
Copy link
Contributor

Thanks for trying to move the ball on this Mark. I rember Bishop’s paper from my NREL days! One of the first papers that I read on the single diode model (SDM).

FWIW, for my own SDM i-from-v and v-from-i calculations, I have resorted to using Scipy’s Newton solver using an explicit ideal diode calculation for initial condition. (I currently use the solver in Halley’s method mode, but I’m not sure that’s better overall given the extra second derivative calculation cost.) One downside to the canned solver is that it is harder to cache parts of the calculation that can be reused. I have also been able to “vectorize” the canned solvers, but I haven’t determined the overhead for this either.

I like straight up Newton solving because the logic around handling ideal (or near-ideal) devices is much more straightforward than LambertW approaches. (It also seems to translate well to the corresponding double-diode model calculations, making for a more consistent code base.) In an earlier pvlib PR of mine, the question was raised about using the residual calculation on the sum of currents at the diode node in the convergence metric as opposed to the LambertW approach. It seemed like a interval calculation of the residual on several “edgy” test cases might be in order to help settle the matter. (If we cannot calculate the residuals accurately for all practical cases, then that would be very interesting indeed.)

I will try to do a similar benchmark against my current SDM code and let you know where the chips fall.

LOGGER.debug('fs_495 speedup = %g', dt_slow / dt_fast)
assert np.isclose(pvs['i_sc'], isc)
assert np.isclose(pvs['v_oc'], voc)
# the singlediode method doesn't actually get the MPP correct
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this kind of a bigish deal?

Copy link
Contributor

@markcampanelli markcampanelli Jan 27, 2018

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I calculate a ~0.02% (absolute) difference between my I_mp_A and V_mp_V values (pvfit) and the pvlib ones for the SunPower module.

pvlib calculation:
OrderedDict([('i_sc', 5.857594583158265), ('v_oc', 57.779106188588685), ('i_mp', 5.374898118822274), ('v_mp', 47.473396494199655), ('p_mp', 255.16466951077766), ('i_x', 5.770970666689348), ('i_xx', 4.048918484984828)])

pvfit calculation:
ff=0.7539291803981164, i_sc_A=5.857594583158266, i_mp_A=5.375969556815134, p_mp_W=255.16475223820953, v_mp_V=47.46395037054039, v_oc_V=57.77910618858893

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Doesn't this change in this PR, because default method is no longer LambertW?

@markcampanelli
Copy link
Contributor

My pvfit code is only 1.5-2x faster than the pvlib singlediode function (pvfit also NOT calculating the i_x and i_xx values). Not using Halley mode didn't seem to make much difference, so I am guessing that I pay a huge penalty for how I vectorize the scipy.optimize.newton() using numpy.frompyfunc().

Overall, I think the Newton approach in this PR shows great promise. (Imagine a 10x speed up in Monte Carlo simulations!) Only issue I see with this/Bishop's implementation using diode voltage in solver is the difficulty/inability specifying specific terminal voltages at which to compute current.

@adriesse
Copy link
Member

I'm not sure whether I should comment here or in the corresponding issue, but anyway...

Would it be possible to put in some slow, brute-force functions into PVLIB that we trust will give precisely the right answers? Then all optimized methods can be compared to them for accuracy.

@markcampanelli
Copy link
Contributor

I have some ideas for a group effort to create a “gold” dataset. I think this calls for a separate issue. I can open one later this weekend if I can find some time. I won’t be upset if someone else beats me to it.

@mikofski
Copy link
Member Author

@thunderfish24 and @adriesse thanks for your comments!

@adriesse I'll address yours first:

Would it be possible to put in some slow, brute-force functions into PVLIB that we trust will give precisely the right answers? Then all optimized methods can be compared to them for accuracy.

That's the beauty of Bishop's method. It's like having your cake and eating it too because it's ...

  • an explicit, "brute force" method that works for single or double diode modes and never BLOWS UP!
  • laser fast and internally vectorized, I suppose you could easily add Numba or NumExpr for even further speedups using llvm or compiletime optimizations.

But this is only used for calculation i and v when ivcurve_pnts > 0. I am still stuck guessing what the value of i will be given any specific constraint on operating condition. Given the output, I usually interpolate to find specific values. I did try to use the canned fminbound method which is essentially a bisection or golden section method, but it was much slower, bummer.

However I did put in a very brutish estimate of voc. And I also put in a limit on the Newton step, because we know at least that operating conditions will be limited to voltages in between voc and zero. I can do a quick check on the convergence rate, to make sure it is always positive, ie: getting closer to the answer. Furthermore we can store the last 2 to 3 newton steps and fit a line or parabola to back search a la Numerical Recipe's in C, chapter 9.6 and 9.7. I've implemented this before in a TMW MATLAB FEX contribution, Newton-Raphson, there is low overhead and it does improve stability better then just damping.

@mikofski
Copy link
Member Author

@thunderfish24 it would be great to find the edgy cases and put them in test_way_faster.py to see how the newton step handles them. Then I'll have a better idea what contingencies we need to take. I will still implement back tracking and line-search from Numerical Recipes in C, chapters 9.6 and 9.7 to improve overall stability.

Also, on a different note, I was wondering, since you mentioned that you use the double-diode method, what temperature and irradiance corrections do you apply, similar to the pvlib.pvsystem.calcparams_desoto method? I think we should leave the double-diode in PVMismatch but should have these corrections if they exist.

Thanks @adriesse and @thunderfish24 for your help and insight!

@adriesse
Copy link
Member

Thank you, Mark, for the contributions!

I'm thinking of a "gold" standard for the function rather than a "golden" data set. Rather than replacing the existing i_from_v and v_from_i with what we think are better/faster algorithms, keep an old one in the library that may be slow but can be relied upon to give the right answer, just for testing.

@mikofski
Copy link
Member Author

@thunderfish24 great idea regarding a "gold" data set. I like it!

Also +1 on caching or memoizitation I just started caching using a hash table for pvmismatch

@markcampanelli
Copy link
Contributor

@mikofski Sadly, I currently have very limited access to/budget for journal article downloads. I have been toying with auxiliary equations for the reverse saturation currents as defined on page 14 here: https://www.slideshare.net/sandiaecis/30-years-of-pv-modeling. An interesting question is if it is "worth it" to use (and perhaps also fit from data) a bandgap temperature coefficient like in the DeSoto SDM. Also, my diode modeling approach doesn't use short-circuit temperature coefficients or air-mass-modifiers to get photocurrent from irradiance (at what spectrum/climate/etc., right?!?), but instead an effective irradiance ratio, see http://ieeexplore.ieee.org/document/7329928/.

I am actively trying to get some inter-comparison I-V curve measurement matrix datasets to do better model comparison/selection, esp. for the double-diode model, as my latest efforts have shown the single-diode model to be not so great for modeling over a range of effective irradiance and temperature (paper in submission). If you know of willing data providers, then I'm happy to share that out too. (My initial proposal is for higher quality x-Si cell measurements from calibration laboratories with well controlled junction temperature and the ability to measure quantum efficiency as a function of temperature for better spectral corrections that are used to calculate the effective irradiance ratio. Outdoor data would most likely use a more/less well-matched PV reference device.)

I want to make it clear that I am not opposed to open-sourcing some of my algorithms, but it would be nice to receive compensation for the significant effort required to make them presentable, sufficiently optimized, well tested, documented, etc., before unleashing them into the wild.

@adriesse
Copy link
Member

adriesse commented Jan 27, 2018

@thunderfish24 We apparently share considerable interests (second paragraph) and desires (last paragraph). Something we could explore in another forum perhaps...

@mikofski
Copy link
Member Author

@adriesse I would propose that bishop88() and est_voc() be our proposed gold-standard methods, see my comments above: they are explicit, ie: brute-force, and reliable because bishop88() uses the same formulas as singlediode. Also both are extensible to 2-diode.

@mikofski
Copy link
Member Author

@adriesse to clarify:

  • bishop88() is brute force, explicit, and at least, if not more accurate than the existing methods
  • est_voc() is a rough guess that bounds the forward bias IV curve
  • Given the two methods you can interpolate or use a bisection method, eg: scipy.optimize.fminbound(func=bishop88, x1=0, x2=est_voc(*args), args) , to get any point you want with any desired tolerance and these methods will never blow up, ever!
  • And these methods can be extended to double diode model (and already are I'm PVMismatch)

I separately implemented the Newton step because I found the canned fminbound method slow, but still faster than the existing singlediode method.

I personally believe the LambertW approach is unstable and less accurate compared to Bishop's method.

f, j = fjx(x0, *x)
newton_step = f / j
# don't let step get crazy
if np.abs(newton_step / x0) > damp:
Copy link
Member Author

@mikofski mikofski Jan 28, 2018

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Change this to (x0 + 1.0) in case x0 is zero.

@mikofski
Copy link
Member Author

mikofski commented Jan 28, 2018

@adriesse and others, I've added 4 new "slow" but reliable ways that use fminbound but are still "faster" than existing singlediode method.

To summarize:

  • The "slow" way is about 2-3x faster than the existing singlediode method.
  • It is bounded by zero diode junction voltage and v_oc so it is guaranteed to be stable always
  • It is reliable because it uses the exact same governing equation for single-diode model that the existing method does
In [1]: %run pvlib/test/test_way_faster.py
DEBUG:__main__:single diode elapsed time = 0.007702[s]
DEBUG:__main__:way faster elapsed time = 0.000358[s]
DEBUG:__main__:spr_e20_327 speedup = 21.514
DEBUG:__main__:single diode elapsed time = 0.007688[s]
DEBUG:__main__:way faster elapsed time = 0.000482[s]
DEBUG:__main__:fs_495 speedup = 15.9502
In [2]: %paste
spr_e20_327 = CECMOD.SunPower_SPR_E20_327
x = pvsystem.calcparams_desoto(
    poa_global=POA, temp_cell=TCELL,
    alpha_isc=spr_e20_327.alpha_sc, module_parameters=spr_e20_327,
    EgRef=1.121, dEgdT=-0.0002677)
il, io, rs, rsh, nnsvt = x
## -- End pasted text --
In [3]: from pvlib import way_faster
In [4]: pvsystem.singlediode(*x)
Out[4]: 
OrderedDict([('i_sc', 5.857594583158265),
             ('v_oc', 57.779106188588685),
             ('i_mp', 5.374898118822274),
             ('v_mp', 47.473396494199655),
             ('p_mp', 255.16466951077766),
             ('i_x', 5.770970666689348),
             ('i_xx', 4.048918484984828)])
In [5]: %timeit pvsystem.singlediode(*x)
100 loops, best of 3: 7.09 ms per loop
In [6]: way_faster.slower_way(*x)
Out[6]: 
OrderedDict([('i_sc', 5.857594581400468),
             ('v_oc', 57.779105959441722),
             ('i_mp', 5.3759693228572862),
             ('v_mp', 47.463952436133873),
             ('p_mp', 255.16475223821305),
             ('i_x', None),
             ('i_xx', None)])
In [7]: %timeit way_faster.slower_way(*x)
100 loops, best of 3: 2.47 ms per loop

* add slow_i_from_v(v, *x) which uses bishop88 and the canned bisecton
 method from scipy.optimize.fminbound to guarantee convergence
* ditto for slow_v_from_i(i, *x)
* ditto for slow_mppt(*x)
* add slower_way(*x) which is the equivalent of faster_way(*x) but using
 the canned bisection method from scipy.optimize.fminbound instead of
 the generic newtom step method.
* also add some better docstrings and comments

Signed-off-by: Mark Mikofski <bwana.marko@yahoo.com>
@mikofski
Copy link
Member Author

@mikofski be DRY and try using canned methods from scipy.optimize like Nelder-Mead (aka simplex) before going it alone with customer Newton step. (nb: this is a note to myself)

@mikofski
Copy link
Member Author

Another note to self, check out the custom golden section method that is already in pvlib.pvsystem vs. scipy.optimize?

@markcampanelli
Copy link
Contributor

markcampanelli commented Jan 29, 2018

@mikofski So when supplying the optional ivcurve_pnts parameter, the current singlediode() function appears to return an equally spaced array of voltages with corresponding current values. I'm guessing that changing this could break consumers. Is your intent to match this behavior with your most recent changes?

@adriesse If we can get reasonable and consistent voltage distributions, then I'm all for using the Bishop method for the "gold function", because accuracy is more important for that purpose. That said, I think it is still useful to have a committed "gold dataset" used for assertions in the test suite, because this can catch regressions if the same gold function starts producing different values. Also, I still think we should attempt some sort of an interval analysis on the gold file/dataset as part of the "blessing".

Lastly, what is pvlib's policy on using compiled code? I have always assumed the library is pure python. I've been meaning to ask the scipy community about better "vectorized" support for the solver libraries, like newton and fminbnd. In fact, I may be missing a more efficient way to do this presently.

UPDATE: According to #410, JIT compiling using numba is a thing for pvlib-python!

``d2p/dv/dvd``
"""
a = np.exp(vd / nnsvt)
b = 1.0 / rsh
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why I vote for moving to gsh := 1/rsh for easier handling of ideal devices.

x = (il, io, rs, rsh, nnsvt) # collect args
# first bound the search using voc
voc_est = est_voc(il, io, nnsvt)
vd = fminbound(lambda vd: (v - bishop88(vd, *x)[1])**2, 0.0, voc_est)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It seems like the most straightforward option for current computation at specific terminal voltages might be to directly implement a vectorized, derivative-free brent algorithm with a sensible root bracket that will "guarantee" convergence. I am investigating better ways to vectorize, say, scipy's brentq for this purpose. It would be a shame to have to rewrite it from scratch just to get performant vectorization.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I haven't reviewed all of this, but just want to suggest that simplicity should be a higher objective than speed for the reference implementation.

Copy link
Member Author

@mikofski mikofski Jan 29, 2018

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I believe what I'm suggesting is the simplest method guaranteed to give the most accurate answer with the governing equations:

  1. Use an explicit method
  2. Estimate open circuit voltage by assuming shunt resistance is infinite
  3. Bound the solution between vd = 0 and vd = voc_est
  4. Use bisection, a "slow but stable" method, guaranteed to find the desired operating condition between the bounds

Further ideological simplifications:

  • use off-the-shelf methods from a mature, familiar, stable library ie: scipy.optimize.fminbound
  • use the SDM equations without algebraic manipulation ie: no Lambert-W
  • methods (with the exception of the faster_way and newton_solver()) are short and simple.

And I'm working on simplifying both of those methods too. E.g. I think we can use the canned scipy.optimize.minimize_scalar() helper function instead of the custom newton_solver in faster_way. Also of course, these are temporary names assuming #410 is implemented.

I am very open to changes. I want to make this as simple and robust as possible. Thanks!

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry, @mikofski , I was not commenting on your code at all but rather referring to @thunderfish24 suggestion to vectorize things.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, I understand now. Thanks @adriesse I'm sorry too. @thunderfish24 May I suggest we move proposals for customized improvements to scipy.optimize subroutines to another issue?

Copy link
Contributor

@markcampanelli markcampanelli Jan 30, 2018

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think Bishop method is a sufficient and extremely simple "ground truth" algorithm at the "uncontrolled" terminal voltages that it produces. However, I am sufficiently skeptical that I would want to see a consensus of 2+ algorithms produce the same result (up to some tolerance) at any other voltage values that require an iterative solver, be it scipy.optimize's fminbnd, newton, or brentq, or lambertw with argument overflow protection logic. If this is not established by algorithm consensus, then something like interval analysis on the model residuals would be convincing. I received push back for doing something even less daring with model residual computations in an early commit to #340. At the very least, for the "slow but stable" method, we should make sure all convergence flags are true, and maybe even make assertions on the absolute residuals for the substituted solution (unless the solver has already verified that, as some solvers use other stopping criteria).

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

To illustrate brentq() approach more concretely: #412

@mikofski
Copy link
Member Author

@thunderfish24 I added a test for numerical errors from floating point arithmetic using SymPy which keeps track of numerical errors and adjusts the precision as needed, see: Numerical Evaluation - Accuracy and error handling. Although this test is slow, is shows that for this particular module there are no floating point errors associated with using np.float64 which has a precision of np.finfo('float64').eps = 2.2204460492503131e-16 on my Mac OSX i5. I doubt we can make this test part of the final PR because it adds SymPy, but I hope that it will satisfy any further concerns regarding sufficient numerical precision.


if output_is_scalar:
return np.asscalar(I)
if method.lower() == 'lambertw':
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same length concern/suggestion as in v_from_i

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

see #497

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

✔️ closes #497 moved indented/nested if: else code to private function

lambda voc, i, iph, isat, rs, rsh, gamma:
brentq(fi, 0.0, voc, args=(i, iph, isat, rs, rsh, gamma))
)
vd = vec_fun(voc_est, current, *args)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

consider the same pattern as suggested in bishop88_i_from_v.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

✔️

* update api.rst and what's new, and test_pvsystem.py everywhere where
it says mpp()
* remove logging, LOGGER, timeing, and CLI for test_singlediode_methods
* remove p from out in pvsystem.singlediode if ivcurve_npts is provided
since not part of original api
* so remove p from test_singlediode_methods.py too
* add singlediode_methods to pvlib/__init__.py
* just assign out once for all methods, outside of the if/else blocks in
pvsystem.singlediode
* since IDE complains, and future maintainer may not realize, initialize
ivcurve_* with NotImplemented instead of None or nothing at all so this
is explicit, you must calculate this or else
* need to create values for i_sc, i_x, and i_xx in bishop88 methods to
fill out
if shape is not None:
v0 = np.broadcast_to(voc_est, shape).copy()
else:
v0 = voc_est
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

maybe this could be used in both [i,v]_from_[v,i] functions? feel free to reject

def _prepare_newton_inputs(i_or_v, args, v0):
        size, shape = _get_size_and_shape((i_or_v,) + args)
        if size > 1:
            args = [np.asarray(arg) for arg in args]
        # newton uses initial guess for the output shape
        # copy v0 to a new array and broadcast it to the shape of max size
        if shape is not None:
            v0 = np.broadcast_to(v0, shape).copy()
        else:
            v0 = voc_est
    return args, v0

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

now I see this can be used in mpp as well. seems like a win. maybe a _prepare_brentq_inputs would help too, but that's a bit less boilerplatey

Copy link
Member Author

@mikofski mikofski Jun 30, 2018

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can we add this to another PR please? see #498

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

✔️ closes #498

@@ -553,43 +553,112 @@ def test_v_from_i(fixture_v_from_i):
assert_allclose(V, V_expected, atol=atol)


@requires_scipy
def test_v_from_i_brentq(fixture_v_from_i):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I believe you can remove the second test if you do this:

@requires_scipy
@pytest.mark.parametrize('method,atol', [('brentq', 1e-11), ('newton', 1e-8)])
def test_v_from_i_methods(fixture_v_from_i, method, atol):
# remove atol spec and set method=method

Copy link
Member Author

@mikofski mikofski Jun 30, 2018

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I will add a follow up pr to address this see #499

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

✔️ closes #499



@requires_scipy
def test_i_from_v_brentq(fixture_i_from_v):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same simplification as suggested for v_from_i method tests above

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

agree

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

please see #499

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

✔️ closes #499

try:
from scipy.optimize import _array_newton
except ImportError:
from pvlib import tools
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

is tools used anywhere in this module? I don't see it... if not, I think the style guide says to remove the import.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

✔️

* remove extra 'p' from out, in test_pvsystem
* make tests in test_singlediode_methods more descriptive, add docstring
and change "fast"->"newton" and "slow"->"brentq"
* use suggested language for notes in max_power_point()
* replace confusing lambdas with real function definitions that are
easier to understand in bishop88_i_from_v and v.v., also add a comment
explaining why this is necessary
* remove import of tools
* add FIXME and comment explaining that we need to remove _array_newton
at some point in the future, and the reason why we have it here
@mikofski
Copy link
Member Author

I hope I've addressed enough of your comments. Perhaps we can address any remaining ones in a few follow on PR's. I've added issues #497, #498, #499, and #500, to address them. Thanks!

@mikofski
Copy link
Member Author

Copy link
Member

@wholmgren wholmgren left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@mikofski I certainly appreciate your desire to wrap this up as soon as possible! I'm not 100% opposed to merging this as-is. My main concern is that some of these changes detract from the long-term maintainability of the code (in my opinion). I tried to propose relatively complete and quick fixes to avoid asking too much more of your time, but I'm sorry if that didn't come across right. Would you be willing to review a pull request to your branch if I try to implement some of these ideas? I think this would be better than working through more PRs to pvlib master. If you're against that, and if other people are ok with the code as is, then I will go ahead and merge.

@@ -1694,8 +1696,8 @@ def sapm_effective_irradiance(poa_direct, poa_diffuse, airmass_absolute, aoi,


def singlediode(photocurrent, saturation_current, resistance_series,
resistance_shunt, nNsVth, ivcurve_pnts=None):
r'''
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@mikofski can you confirm that the equation still renders properly? If I remember correctly, I added r to some of the doc strings at one point because it was necessary for the equations to render properly.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yes, rendered okay for me - I read somewhere that you can avoid using r""":math:`\frac{numerator}{denominator}`""" by escaping the backslashes, and it works for me:

"""
.. math::
   \\frac{numerator}{denominator}
"""

image

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for the check. That sounds familiar. Maybe I added "r" in a few places because I thought it was easier to do than escape the backslashes. One could argue that it improves readability too, but let's not fret over it.

"""
# make IDE happy
ivcurve_v = NotImplemented
ivcurve_i = NotImplemented
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I disagree with your IDE. I think it's better if an exception is thrown on failed assignment if for some reason the nested if statements fail to do the right thing.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

✔️

This function uses Brent's method by default because it is guaranteed to
converge.
"""
mpp_fun = partial(singlediode_methods.bishop88_mpp, method=method.lower())
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

are you using partial here because you're thinking ahead to this function supporting additional algorithms for determining the max power point? If so, would we have a model kwarg to go along with method? If not, why use partial at all here?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

✔️ no, it was leftovers, seemed convenient at the time, but I think it was just sloppy - I've removed it


if output_is_scalar:
return np.asscalar(V)
if method.lower() == 'lambertw':
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

As is, the increased indentation makes the code break the line length limit. But fixing that would actually make it harder to implement #497. And the diff shows all these lines that are not actually changing, which is a small but non-zero impact on long term maintainability. I think a private function is a good compromise.

resistance_shunt : numeric
shunt resitance [ohms]
nNsVth : numeric
product of thermal voltage ``Vth`` [V], diode ideality factor ``n``, and
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this line and a handful of other documentation lines are longer than 79 characters. We're not as strict about 79 characters for docstrings, particularly for equations/examples, but many of these should be fixed.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

✔️

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

✔️

"""
a = np.exp(diode_voltage / nNsVth)
b = 1.0 / resistance_shunt
i = photocurrent - saturation_current * (a - 1.0) - diode_voltage * b
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just a thought. It's fine with me if it stays as is.

if __name__ == '__main__':
expected = generate_numerical_precision()
expected.to_csv(DATA_PATH)
test_numerical_precision()
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would still like some clarification on how we envision this working with #426.

* closes pvlib#497
* remove lines setting ivcurve points to NotImplemented, so that they
raise an error if they are not set
* replace nested code in if: else condition for lambertw with new
private method "_lambertw" in singlediode_methods.py
* need to break out returns from private methods, there are so many
* move _golden_sect_DataFrame to pvlib.tools
* replace nested code in if: else condition for lambertw with new
private methods _lambertw_v_from_i and _lambertw_i_from_v in
singlediode_methods.py
* move _pwr_optfcn to singlediode_methods.py
* explain how to generate the high-precision test data by running the
file from the command line
* closes pvlib#500
* also add a comment that we are creating some temporary values to make
calculations simpler
* also use more descriptive names instead of a, b, c, use v_star, g_sh,
and g_diode
* closes pvlib#498
* adds _prepare_newton_inputs(i_or_v_tup, args, v0)
* also replace verbose comment with more descriptive one to copy v0 if
it's an array to use for initial guess
* replace partial functions v_from_i_fun, i_from_v_fun, and mpp_fun with
full names, and extra keyword argument method=method.lower() everywhere
@mikofski
Copy link
Member Author

mikofski commented Jul 10, 2018

@wholmgren I hope I have responded to your comments and suggestions sufficiently. If you or anyone else has more feedback, please let me know. Otherwise, if it is ready to merge, then I hope it will be a useful addition to pvlib. Thanks!

PS: I left a comment regarding the numerical precision test. If the test is redundant, it's okay with me to remove it. I don't have any comments regarding #426, as I haven't reviewed it, and honestly I don't really understand it. Sorry.

PPS: My changes should close #497, close #498, close #499, and close #500. I added commit messages so they should close automatically

Copy link
Member

@wholmgren wholmgren left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks great @mikofski! Thank you for this major contribution and thank you for your patience!

@@ -1694,8 +1696,8 @@ def sapm_effective_irradiance(poa_direct, poa_diffuse, airmass_absolute, aoi,


def singlediode(photocurrent, saturation_current, resistance_series,
resistance_shunt, nNsVth, ivcurve_pnts=None):
r'''
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for the check. That sounds familiar. Maybe I added "r" in a few places because I thought it was easier to do than escape the backslashes. One could argue that it improves readability too, but let's not fret over it.


if output_is_scalar:
return np.asscalar(V)
if method.lower() == 'lambertw':
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

thanks.

if __name__ == '__main__':
expected = generate_numerical_precision()
expected.to_csv(DATA_PATH)
test_numerical_precision()
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, let's keep this test. I think my confusion came out of some of the discussion in #411 and my mind locked into the wrong idea.

@wholmgren wholmgren merged commit 2498040 into pvlib:master Jul 10, 2018
@mikofski mikofski deleted the faster_way branch July 10, 2018 18:51
@mikofski mikofski changed the title ENH: WIP: add methods and tests for a explicit IV curve calculation of single-diode model ENH: add methods and tests for a explicit IV curve calculation of single-diode model Aug 3, 2018
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Single diode model helper function to wrap Lambert w and other methods singlediode is slow
5 participants