Improved method for computing Marion's IAM-correction to horizon diffuse POA irradiance

[markcampanelli]

As mentioned in #1402, here is a simpler and ~2x faster (and more memory efficient?) way to compute Marion's IAM-correction to horizon diffuse POA irradiance. Copied this from a Jupyter notebook:

import timeit

from matplotlib import pyplot
import numpy
from pvlib.iam import marion_integrate, physical
from scipy.integrate import quad

def integrand_bot(phi_rad, theta_tilt_rad):
    return -numpy.sin(phi_rad) * numpy.sin(theta_tilt_rad)

def integrand_top(phi_rad, theta_tilt_rad, iam_fun):
    cos_aoi = integrand_bot(phi_rad, theta_tilt_rad)
    return iam_fun(numpy.degrees(numpy.arccos(cos_aoi))) * cos_aoi

def marion_integrate_horizon(theta_tilt_rad, iam_fun):
    two_pi = 2*numpy.pi
    # quad returns result and accuracy as a tuple, so take first element.
    return quad(integrand_top, numpy.pi, two_pi, args=(theta_tilt_rad, iam_fun))[0] / (
        quad(integrand_bot, numpy.pi, two_pi, args=theta_tilt_rad)[0]
    )

theta_tilts_deg = numpy.linspace(1, 90, 90)
theta_tilts_rad = theta_tilts_deg * numpy.pi / 180
marion_integrate_horizon_results_old = []
marion_integrate_horizon_results_new = []

for theta_tilt_deg, theta_tilt_rad in zip(theta_tilts_deg, theta_tilts_rad):
    marion_integrate_horizon_results_old.append(
        marion_integrate(physical, theta_tilt_deg, "horizon")
    )
    marion_integrate_horizon_results_new.append(
        marion_integrate_horizon(theta_tilt_rad, physical),
    )

fig, ax = pyplot.subplots()
ax.plot(
    theta_tilts_deg, marion_integrate_horizon_results_old, label="pvlib-python",
)
ax.plot(
    theta_tilts_deg, marion_integrate_horizon_results_new, "--", label="new method",
)
ax.legend(loc="lower right")
ax.set_xlabel("Module Tilt [deg]")
ax.set_ylabel("IAM")
ax.set_title("Comparison of Marion Diffuse-Horizon IAM Corrections")

timeit.timeit(
    lambda: marion_integrate(physical, theta_tilts_deg[24], "horizon"), number=10000
)
timeit.timeit(
    lambda: marion_integrate_horizon(theta_tilts_rad[24], physical), number=10000
)

The key outputs are:

• Time of existing algorithm for 10000 calculations: 6.168394709005952 s
• Time of new algorithm for 10000 calculations: 3.1387881669797935 s
• A plot of the results comparison:

The new algorithm still expects radians and doesn't yet handle the special case where POA tilt is zero.

I am working on a proper math writeup to show how the two integrands are derived.

cc @cwhanse @kandersolar

[markcampanelli]

From this comment:

(*) I think the calculation may be wrong when applied to ground-reflected irradiance, since an integrating element is defined by a solid angle, not by an equal area portion of the viewed ground. I think that a proper view factor is needed in the integrand. The sky diffuse calculation doesn't have this problem because it is assumed that the irradiance source is the hemisphere rather than the ground plane.

@cwhanse Does pvlib.irradiance.get_ground_diffuse properly account for the different view factor of the POA "looking" at the ground plane (as opposed to the surface of a sphere)? (I admit that I am not well versed in view-factor computations, but I'm wanting to learn.)

I'd consider incorporating the above horizon-integral implementation improvement in a bigger paper on a view-factor improved ground-diffuse correction, assuming your hunch is correct that it is missing in the Marion method.

[cwhanse]

get_ground_diffuse (and all ground diffuse models, I think) computes ground-reflected irradiance as

GHI x albedo x VF

where VF is the view factor from a (possibly finite) planar element to the half-plane: VF = (1 - cos(tilt)) / 2.

This formula is usually derived as the complement of the view factor from the receiving tilted surface to the sky dome.

Several assumptions are built into that VF formula:

1. the radiating half-plane is semi-infinite
2. radiation from the half-plane is Lambertian (isotropic)
3. the two surfaces share an edge (i.e., the module's lower edge sits on the ground).