Implement iam.fedis and iam.schlick

[kandersolar]

• [x] Closes #1564
• [x] I am familiar with the contributing guidelines
• [x] Tests added
• [x] Updates entries in docs/sphinx/source/reference for API changes.
• [x] 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`).
• [x] New code is fully documented. Includes numpydoc compliant docstrings, examples, and comments where necessary.
• [x] Pull request is nearly complete and ready for detailed review.
• [x] Maintainer: Appropriate GitHub Labels (including remote-data) and Milestone are assigned to the Pull Request and linked Issue.

Yu Xie from NREL presented a new IAM model at the recent PVPMC meeting (slides) and wanted to contribute it to pvlib. This one is a little unusual in that it calculates the IAM factor as a relative quantity, taking the indices of refraction of both the PV surface and the accompanying pyranometer into account.

[kandersolar]

I suppose we could move the diffuse calculations to a separate function like was done with martin_ruiz. It's a little strange to return scalar sky & ground IAM but time series direct IAM as the current code does for a fixed tilt simulation.

[kandersolar]

Offline discussion with @cwhanse and the model's author determined that the parameter this PR calls n_ref may have some niche uses but should be set to the same value as n for standard PV modeling workflows (unlike as assumed in the reference). That is now the default behavior.

I also added pvlib.iam.schlick in this PR as suggested in #1564.

[adriesse]

Should/does fedis direct be equal to schlick? If not, why not?

[kandersolar]

Should/does fedis direct be equal to schlick? If not, why not?

Close-ish but not identical. FEDIS is using the true Fresnel equations for the direct component so it cannot match exactly. The Schlick approximation only comes in for the diffuse components: in the absence of an analytical integration of the real Fresnel equations, FEDIS's diffuse components are instead an integration of the Shlick approximation. It's the alternative to Marion's integration method -- whereas Marion's is an approximate integration of an exact integrand, FEDIS is an exact integration of an approximate integrand.

[adriesse]

Close-ish but not identical. FEDIS is using the true Fresnel equations for the direct component

I see, but we already have that in iam.physical, right?

[kandersolar]

I see, but we already have that in iam.physical, right?

Yes, with two minor differences: FEDIS does not consider extinction (i.e. it assumes K=0 in iam.physical), and FEDIS has the option to normalize the IAM profile with respect to another device with a different index of refraction. The latter is what prevents iam.physical from being a superset of FEDIS's direct IAM calculation. A separate index of refraction for the normal incidence normalization coefficient is probably only useful in unusual modeling scenarios, but nevertheless it is part of the model so I include it here.

For context, I have the sense that the FEDIS paper to some extent views its diffuse IAM equations as the primary scientific contribution with the direct component being included for completeness.

[adriesse]

FEDIS has the option to normalize the IAM profile with respect to another device with a different index of refraction.

Normally IAM is a property or characteristic of one device, so we should try to avoid confusion. If a ratio is needed for some unusual modeling scenarios, this can be calculated with any pair of existing IAM models and parameters, I would think.

[kandersolar]

The behavior as currently implemented is to default to equal refractive index so that everything is for a single device. I'm open to removing the second refractive index altogether so long as we include sufficient justification in the docstring, but this would be a substantial deviation from the reference which, for better or worse, places a lot of focus on the multiple refractive index aspect.

[adriesse]

Should/does fedis direct be equal to schlick? If not, why not?

Close-ish but not identical. FEDIS is using the true Fresnel equations for the direct component so it cannot match exactly. The Schlick approximation only comes in for the diffuse components: in the absence of an analytical integration of the real Fresnel equations, FEDIS's diffuse components are instead an integration of the Shlick approximation. It's the alternative to Marion's integration method -- whereas Marion's is an approximate integration of an exact integrand, FEDIS is an exact integration of an approximate integrand.

This is a good way to put the distinction. Based on my most recent re-read, FEDIS is an exact integration of an approximate integrand plus a scaling/correction factor w to make up for the fact that the integrating the approximate integrand (Schlick) doesn't actually produce the desired result.

[kandersolar]

Maybe of interest... comparison with analogous calculations from iam.marion_integrate (compare with Fig 3 in the reference):

Source

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

n = 1.6
aoi = np.linspace(0, 90)
surface_tilt = np.linspace(0, 90, num=20)

def physical_no_extinction(aoi):
    return pvlib.iam.physical(aoi, n=n, K=0)


fig, axes = plt.subplots(3, 1)

# %%

physical = physical_no_extinction(aoi)
fedis = pvlib.iam.fedis(aoi, 0, n=n)['direct']

axes[0].plot(aoi, fedis, label='fedis')
axes[0].plot(aoi, physical, label='physical (no extinction)')
axes[0].set_ylabel('IAM (direct)')
axes[0].set_xlabel('aoi [degrees]')
axes[0].legend()

# %%

fedis = pvlib.iam.fedis(0, surface_tilt, n=n)['sky']

polycoeffs = [2.77526e-09, 3.74953, -5.18727, 3.41186, -1.08794, 0.136060]
w = np.polynomial.polynomial.polyval(n, polycoeffs)

marion_schlick = pvlib.iam.marion_integrate(pvlib.iam.schlick,
                                            surface_tilt=surface_tilt, region='sky', num=1000)
marion_physical = pvlib.iam.marion_integrate(physical_no_extinction,
                                             surface_tilt=surface_tilt, region='sky', num=1000)

axes[1].plot(surface_tilt, fedis, label='fedis')
axes[1].plot(surface_tilt, marion_schlick * w, label='marion + schlick (with $w$)')
axes[1].plot(surface_tilt, marion_physical, label='marion + physical')
axes[1].set_ylabel('IAM (sky)')
axes[1].set_xlabel('surface_tilt [degrees]')
axes[1].legend()

# %%

fedis = pvlib.iam.fedis(0, surface_tilt, n=n)['ground']
marion_schlick = pvlib.iam.marion_integrate(pvlib.iam.schlick,
                                            surface_tilt=surface_tilt, region='ground', num=1000)
marion_physical = pvlib.iam.marion_integrate(physical_no_extinction,
                                             surface_tilt=surface_tilt, region='ground', num=1000)

axes[2].plot(surface_tilt, fedis, label='fedis')
axes[2].plot(surface_tilt, marion_schlick * w, label='marion + schlick (with $w$)')
axes[2].plot(surface_tilt, marion_physical, label='marion + physical')
axes[2].set_ylabel('IAM (ground)')
axes[2].set_xlabel('surface_tilt [degrees]')
axes[2].legend()

[cwhanse]

@kanderso-nrel that bottom panel is an eyebrow-raiser. I had not seen that result before. Although ground-reflected irradiance onto the front surface is a very minor part of plane-of-array, that figure suggests it should be further reduced by roughly a factor of 2, for many installations.

[kandersolar]

This PR now implements four separate functions (schlick/fedis, direct/diffuse) as @adriesse suggested. For interest here is a plot similar to above but comparing two indices of refraction to show the effect of FEDIS's weight coefficient w:

Source

import matplotlib.pyplot as plt
import numpy as np
from pvlib.iam import physical, fedis, fedis_diffuse, schlick, schlick_diffuse, marion_integrate

aoi = np.linspace(0, 90)
surface_tilt = np.linspace(0, 90, num=20)

def physical_no_extinction(aoi):
    return physical(aoi, n=n, K=0)


fig, axes = plt.subplots(3, 2)

for i, n in enumerate([1.32, 1.5]):

    iam_physical = physical_no_extinction(aoi)
    iam_fedis = fedis(aoi, n=n)

    axes[0, i].plot(aoi, iam_fedis, label='fedis')
    axes[0, i].plot(aoi, iam_physical, label='physical (no extinction)')
    axes[0, i].set_ylabel('IAM (direct)')
    axes[0, i].set_xlabel('aoi [degrees]')
    axes[0, i].legend()

    iam_schlick = schlick_diffuse(surface_tilt)[0]
    iam_fedis = fedis_diffuse(surface_tilt, n=n)[0]
    marion_schlick = marion_integrate(schlick, surface_tilt=surface_tilt, region='sky', num=1000)
    marion_physical = marion_integrate(physical_no_extinction, surface_tilt=surface_tilt, region='sky', num=1000)

    axes[1, i].plot(surface_tilt, iam_fedis, label='fedis_diffuse')
    axes[1, i].plot(surface_tilt, iam_schlick, label='schlick_diffuse')
    axes[1, i].plot(surface_tilt, marion_schlick, label='marion + schlick')
    axes[1, i].plot(surface_tilt, marion_physical, label='marion + physical')
    axes[1, i].set_ylabel('IAM (sky)')
    axes[1, i].set_xlabel('surface_tilt [degrees]')
    axes[1, i].legend()

    iam_schlick = schlick_diffuse(surface_tilt)[1]
    iam_fedis = fedis_diffuse(surface_tilt, n=n)[1]
    marion_schlick = marion_integrate(schlick, surface_tilt=surface_tilt, region='ground', num=1000)
    marion_physical = marion_integrate(physical_no_extinction, surface_tilt=surface_tilt, region='ground', num=1000)

    axes[2, i].plot(surface_tilt, iam_fedis, label='fedis_diffuse')
    axes[2, i].plot(surface_tilt, iam_schlick, label='schlick_diffuse')
    axes[2, i].plot(surface_tilt, marion_schlick, label='marion + schlick')
    axes[2, i].plot(surface_tilt, marion_physical, label='marion + physical')
    axes[2, i].set_ylabel('IAM (ground)')
    axes[2, i].set_xlabel('surface_tilt [degrees]')
    axes[2, i].legend()

    axes[0, i].set_title(f'n={n}')

[adriesse]

Nicely done. Let's see what other comments appear here over the coming days...

On a related note, I think it would be great if the original author could submit two examples to pvlib-python: one example illustrating an important use case for n_ref ; and one to replicate the validation comparisons that are presented in the publication.

[adriesse]

Perhaps we should invite a additional reviewer on this one?

[adriesse]

It turns out that term1 in fedis_diffuse() is exactly the same as the weighting function w in fedis(), just in simplified form. I suggest making the code the same in both functions and referring to this ratio/weight/multiplier as w0. (Or the calculation could even be put in a separate helper function.)

For clarity, I would suggest using the long version because it shows what's actually going on, and putting the short version in a comment with reference to the paper.

[kandersolar]

For clarity, I would suggest using the long version because it shows what's actually going on, and putting the short version in a comment with reference to the paper.

I have kept the equations as-is so they continue being directly tied to the reference, but the latest commit at least leaves a comment noting the equivalence for any interested readers. Hopefully that achieves the same goal?

Perhaps we should invite a additional reviewer on this one?

I wouldn't mind input from additional reviewers here, but the three of us (@adriesse, @cwhanse, myself) reaching consensus would be good enough for me.

[kandersolar]

Following offline discussion, we have decided to move forward here with only the two schlick functions. I'll wait a few days to merge this to leave opportunity to review for anyone who wants to.

[kandersolar]

Thanks again for the thoughtful and productive discussion @adriesse and @cwhanse

[adriesse]

No problem. I learned a few things in the process!