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).

[markcampanelli]

I wrote up the derivation, which led me to realize an even simpler expression for the integral:

import numpy
from scipy.integrate import quad

def integrand(phi_rad, theta_tilt_rad, iam_fun):
    return iam_fun(numpy.degrees(numpy.arccos(numpy.sin(phi_rad) * numpy.sin(theta_tilt_rad)))) * numpy.sin(phi_rad) / 2

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

Marion_Horizon_Integral_Alternative.pdf

@cwhanse I am happy to collaborate on a publication by adding the view-factor correction to the ground diffuse integrals.

@cbcrespo I struggled to follow Marion's original derivation. Maybe this helps in understanding? Also, would you benefit by publishing a paper on this? (See also comment above for Cliff.)

[ramaroesilva]

@markcampanelli, it's meritful how you got your hands dirty and made the effort to understand what was happening under the hood. while it was easy for me to understand the practical challenge of "marion + tilt timeseries" and how interpolation could easily, efficiently, and accurately handle this, when more complicated geometry and some calculus matters it becomes out of my reach 😅

personally, and I have yet to discuss this with @AdamRJensen about it, I think that a brief technical note could be a good diffusion mechanism for a potential new, interpolation-based IAM timeseries generation (which @cbcrespo is now testing).

however, while that is a high-level modification (in the sense that it doesn't change the base code) that only helps for trackers, yours is a low-level one which benefits all cases and - from what you say - even corrects a potential gap in marion's ground calculations. so I think that in a scenario where both cases could be publishable, I would keep them separate.