Support ideal PV devices (#324) #340
Conversation
markcampanelli
changed the title
|
I updated pvsystem.i_from_v() to use a shunt conductance derived from the shunt resistance input, without changing the function interface. Rsh=np.inf with Rs>0 test case should now be passing. Note that computing short-circuit current with IL=0 (from zero irradiance) also appears to be improved. I'm going to start on the analogous changes to pvsystem.v_from_i(), and hopefully get some stability tests written for parameters approaching the ideal cases. |
|
Mark,
Looks good so far. Thanks for taking this on.
|
|
Thanks @thunderfish24. I don't have time to comment on the substance at this point, but here are some more style-oriented things that I noticed with a quick glance:
flake8 is a good tool for tracking down some style issues. |
|
@wholmgren Thanks for the additional tips. |
|
Here is an example of computational instability in the existing v_from_i() as Rsh->inf, and its correction in v_from_i_alt() . For this particular parameter set, the Rsh value has to be fairly large before the voltage computation goes off the rails. import pvlib
import numpy as np
Rsh = 190.
Rs = 1.065
nNsVth = 2.89
I0 = 7.05196029e-08
IL = 10.491262
I = 0.
Voc_ideal = pvlib.pvsystem.v_from_i_alt(Rsh=np.inf, Rs=Rs, nNsVth=nNsVth, I=I, I0=I0, IL=IL)
print('Voc_ideal_V = ' + str(Voc_ideal))
print()
Rsh_pvlib_limit_vec = Rsh*np.logspace(0, 15, num=15)
for Rsh_mod in Rsh_pvlib_limit_vec:
V_pvlib = pvlib.pvsystem.v_from_i(Rsh_mod, Rs, nNsVth, 0., I0, IL)
print('Rsh_mod_Ohm = ' + str(Rsh_mod) + ', V_pvlib_V = ' + str(V_pvlib) + ': current_sum_A = ' + str(pvlib.pvsystem.current_sum_at_diode_node(V_pvlib, I, IL, I0, nNsVth, Rs, Rsh_mod)))
print()
for Rsh_mod in Rsh_pvlib_limit_vec:
V_pvlib_alt, meta_dict = pvlib.pvsystem.v_from_i_alt(Rsh=Rsh_mod, Rs=Rs, nNsVth=nNsVth, I=I, I0=I0, IL=IL, return_meta_dict=True)
print('Rsh_mod_Ohm = ' + str(Rsh_mod) + ', V_pvlib_alt_V = ' + str(V_pvlib_alt) + ': current_sum_A = ' + str(meta_dict['current_sum_at_diode_node']) + ', inf_Rsh_idx = ' + str(meta_dict['inf_Rsh_idx']))
print('inf_Rsh_idx = ' + str(meta_dict['inf_Rsh_idx']))
print()
print()
V = 0.
Isc_ideal = pvlib.pvsystem.i_from_v_alt(Rsh=np.inf, Rs=Rs, nNsVth=nNsVth, V=V, I0=I0, IL=IL)
print('Isc_ideal_A = ' + str(Isc_ideal))
print()
Rsh_pvlib_limit_vec = Rsh*np.logspace(0, 15, num=15)
for Rsh_mod in Rsh_pvlib_limit_vec:
V_pvlib = pvlib.pvsystem.v_from_i(Rsh_mod, Rs, nNsVth, Isc_ideal, I0, IL)
print('Rsh_mod_Ohm = ' + str(Rsh_mod) + ', V_pvlib_V = ' + str(V_pvlib) + ': current_sum_A = ' + str(pvlib.pvsystem.current_sum_at_diode_node(V_pvlib, Isc_ideal, IL, I0, nNsVth, Rs, Rsh_mod)))
print()
for Rsh_mod in Rsh_pvlib_limit_vec:
V_pvlib_alt, meta_dict = pvlib.pvsystem.v_from_i_alt(Rsh=Rsh_mod, Rs=Rs, nNsVth=nNsVth, I=Isc_ideal, I0=I0, IL=IL, return_meta_dict=True)
print('Rsh_mod_Ohm = ' + str(Rsh_mod) + ', V_pvlib_alt_V = ' + str(V_pvlib_alt) + ': current_sum_A = ' + str(meta_dict['current_sum_at_diode_node']) + ', inf_Rsh_idx = ' + str(meta_dict['inf_Rsh_idx']))
print('inf_Rsh_idx = ' + str(meta_dict['inf_Rsh_idx'])) |
|
@wholmgren @cwhanse I think this is in a state where a preliminary review would be helpful now. Some notes--
|
|
I don't think a grace period is needed if we bump from version 0.4.x to 0.5.0. We would want to make sure that the new functions pass all of the old tests, though. Is the 10x performance degradation important? From only looking at the code I am guessing that it's still plenty fast. For example, 10x slower is probably acceptable if it doesn't significantly impact the runtime of Is I don't think you need to worry about the pandas testing for these functions. It would help if you could handle most of the errors reported when you run It makes things easier in the long run. Some of your longer lines could be shortened by using positional arguments instead of keyword arguments. |
|
@wholmgren Thanks for the feedback. I have had to put this on hold for a bit while I finish up a paper. |
|
I have decided that I would rather forgo including benchmarks and stability tests into the code, and hopefully use this PR as documentation of sufficient testing. It looks like the new code is roughly 2x slower: Rsh = 20.
Rs = 0.1
nNsVth = 0.5
I = 3.
I0 = 6.e-7
IL = 7.
N = 10000
t_start = time.time()
for k in range(N):
pvlib.pvsystem.v_from_i(Rsh, Rs, nNsVth, I, I0, IL)
t_end = time.time()
reg_delta_t = t_end - t_start
print('pvlib.pvsystem.v_from_i: average seconds per computation ' + str(reg_delta_t/N))
t_start = time.time()
for k in range(N):
pvlib.pvsystem.sdm_v_from_i(Rsh=Rsh, Rs=Rs, nNsVth=nNsVth, I=I, I0=I0, IL=IL)
t_end = time.time()
alt_delta_t = t_end - t_start
print('pvlib.pvsystem.sdm_v_from_i: average seconds per computation ' + str(alt_delta_t/N))
print('pvlib.pvsystem.sdm_v_from_i is {} times slower'.format(alt_delta_t/reg_delta_t))gives pvlib.pvsystem.v_from_i: average seconds per computation 7.457981109619141e-05
pvlib.pvsystem.sdm_v_from_i: average seconds per computation 0.00014008252620697023
pvlib.pvsystem.sdm_v_from_i is 1.8782901719380174 times slowerand Rsh = 20.
Rs = 0.1
nNsVth = 0.5
V = 0.
I0 = 6.e-7
IL = 7.
N = 10000
t_start = time.time()
for k in range(N):
pvlib.pvsystem.i_from_v(Rsh, Rs, nNsVth, V, I0, IL)
t_end = time.time()
reg_delta_t = t_end - t_start
print('pvlib.pvsystem.i_from_v: average seconds per computation ' + str(reg_delta_t/N))
t_start = time.time()
for k in range(N):
pvlib.pvsystem.sdm_i_from_v(Rsh=Rsh, Rs=Rs, nNsVth=nNsVth, V=V, I0=I0, IL=IL)
t_end = time.time()
alt_delta_t = t_end - t_start
print('pvlib.pvsystem.sdm_i_from_v: average seconds per computation ' + str(alt_delta_t/N))
print('pvlib.pvsystem.sdm_i_from_v is {} times slower'.format(alt_delta_t/reg_delta_t))gives pvlib.pvsystem.i_from_v: average seconds per computation 5.51548957824707e-05
pvlib.pvsystem.sdm_i_from_v: average seconds per computation 0.00012906100749969482
pvlib.pvsystem.sdm_i_from_v is 2.3399737352184955 times slower |
|
@wholmgren "Is return_meta_dict helpful for anything besides debugging? It seems to me that the few cases where someone would want that data could be better handled by manually computing it or using the debugger. Or did I miss something?" The |
|
@wholmgren @cwhanse Assuming this passes CI, it is ready for review for merging into master. I am strongly in favor of having a deprecation period for the existing |
markcampanelli
changed the title
| [3.01079860e+00, 2.88414114e+00, | ||
| 3.10862447e-14], | ||
| [6.00726296e+00, 5.74622046e+00, | ||
| 0.00000000e+00]]))]) |
There was a problem hiding this comment.
Be gone nan's from zero irradiance calculations!
There was a problem hiding this comment.
I suspect the NaNs arise from this line in calcparams_desoto:
Rsh = Rsh_ref * (irrad_ref / poa_global)
We should change to use np.divide here
There was a problem hiding this comment.
This is "fixed" in the sense that Rsh is computed as +inf when poa_global is zero, and now the i_from_v and v_from_i code can properly handle this infinity.
|
Looking into these failures on only some python configurations.
Update: Looks like the same failure might happen in master ATM. |
|
@wholmgren Cliff had concerns about documenting my changes to the published algorithm. Any suggestions here? I'd like to get a thorough review after posting a description somewhere good. |
pvlib/pvsystem.py
Outdated
| w = logargW | ||
| for _ in range(0, 3): | ||
| w = w * (1 - np.log(w) + logargW) / (1 + w) | ||
| lambertwterm[idx_w] = w |
There was a problem hiding this comment.
Looking at the coverage report, we never hit lines 1937-1949, even though one of the pre-existing unit tests (copied to the fixture in this PR) is supposed to cover this branch of the code. It may be that recasting in terms of Gsh helps avoid the overflow. @cwhanse Any ideas here?
There was a problem hiding this comment.
The unit test was from #226, I think.
There was a problem hiding this comment.
Not sure why the parameters from #225 don't cause entry into that code block. They cause overflow on the np.exp() call.
def test_v_from_i_bigger():
# 1000 W/m^2 on a Canadian Solar 220M with 20 C ambient temp
# github issue 225
output = pvsystem.v_from_i(190, 1.065, 2.89, 0, 7.05196029e-08, 10.491262)
'''
There was a problem hiding this comment.
So it appears that the Gsh recast helps the existing test's computation stay finite:
In[5]: pvlib.pvsystem.v_from_i(190, 1.065, 2.89, 0, 7.05196029e-08, 10.491262)
Using LambertW with argW=[ 1.64103617e+294] and lambertwterm=[ 670.94665556]
Out[5]: 54.303958833791285
However, the original parameters from #225 serve as a good test, I think:
In[5]: pvlib.pvsystem.v_from_i(381.68, 1.065, 2.681527737715915, 0, 1.8739027472625636e-09, 5.1366949999999996)
Using LambertW with argW=[ inf] and lambertwterm=[ inf]
/Users/markcampanelli/Documents/Programming/pvlib-python/pvlib/pvsystem.py:1926: RuntimeWarning: overflow encountered in exp
np.exp((-I[idx] + IL[idx] + I0[idx]) / (Gsh[idx]*nNsVth[idx]))
Out[4]: 58.19323124611128
Note that in followup work I might be interested in checking how accurate lambertw() is for such large inputs.
I will add this parameter set as a unit test and verify the coverage
|
@cwhanse @wholmgren Ready again for review. Looking for some feedback about testing LambertW input overflow in |
| 'I': 0., | ||
| 'I0': 1.8739027472625636e-09, | ||
| 'IL': 5.1366949999999996, | ||
| 'V_expected': 58.19323124611128 |
There was a problem hiding this comment.
@cwhanse Would you mind running this parameter set though the Matlab version to "independently" verify my V_expected value?
There was a problem hiding this comment.
Awesome! Thanks.
|
Note: Some informal followup performance testing suggests that the additional computation and indexing logic causes these new computations to take slightly more than 2x time. |
|
@wholmgren Still waiting on the Landscape result to verify increased code coverage with the new test. Other than that, I think it good to go except for any further reviewers you might try to muster. |
|
Makes sense –roughly a fixed computation cost to solve the ideal (or non-ideal case) using numpy vectorized operations.
Maybe we should switch to have the non-ideal as the default branch, and the ideal as the exception? Currently the code always computes the ideal branch then checks for the non-ideal.
My guess is that the ideal branch isn’t the typical use.
|
|
So I currently use the ideal computation first as an "easy" way to get the correct shape of the output. It's probably worth investigating if a bunch of "input shape logic" is ultimately faster here. |
|
Coverage of my changes is good now. |
|
@cwhanse It looks like I may be able to use numpy.broadcast to some advantage here. |
|
@cwhanse I did not see the speedup that I hoped for (still about 2-2.5x slower), but I the code is cleaner now. I think the improved computational range and stability are worth it at this point. |
|
@wholmgren Note that the test failure before the last build appears to be unrelated to my changes. |
|
@mikofski Would you be able to do a quick review of this PR? |
|
@thunderfish24 I'll try to take a quick look next week |
|
@wholmgren Can you provide a second code review so that perhaps we can get this merged? |
|
Yes, sorry, I have a number of pvlib things I’m aiming to get to by the end
of the week.
…On Tue, Sep 26, 2017 at 7:22 PM Mark Campanelli ***@***.***> wrote:
@wholmgren <https://github.com/wholmgren> Can you provide a second code
review so that perhaps we can get this merged?
—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
<#340 (comment)>,
or mute the thread
<https://github.com/notifications/unsubscribe-auth/AELiRxRgPJ_7FOVq4NFinCTLohojpRtMks5smbF9gaJpZM4N6u62>
.
|
|
@thunderfish24 this looks great! Thank you for your patience and flexibility as we worked this out. |
|
Thanks @wholmgren. |
Based on feedback, I decided to wait on implementing my full algorithm that uses the absolute residual in the sum of currents at the diode node to determine when an explicit computation for an ideal device is a better approximation for a very near ideal device. This PR documents one really edgy case where this appears to have been a better numerical computation.
By reformulating some of the LambertW formulas, this PR does make many computations more stable, so it appears to be harder to get into a poorly behaved edge case.
In future work, I propose changing the shunt resistance parameter Rsh to a shunt conductance parameter Gsh=1./Rsh, which appears to be much better behaved numerically for ideal and near-ideal devices, i.e., Gsh is approaching/at zero. This should avoid a lot of logic bloat in considering yet another limiting case in pvsystem.i_from_v() where Rsh is so large that it is essentially infinite and the shunt term in the model must be dropped out. (Kudos to Paul Basore for recommending that I switch to Gsh in my early modeling days at NREL.)
git diff upstream/master -- pvlib | flake8 --diffdocs/sphinx/source/whatsnewfile for all changes and thedocs/sphinx/source/api.rstfor new API