Yield calculations

This section describes the equations applied in SunSolve P90 to determine the energy yield of a PV plant.

1. General approach

This approach to solving yield is much simpler than the sophisticated approach taken by SunSolve Yield. For example, it solves the optics by the view-factor method rather than ray tracing, and the electronic losses with multipliers rather than SPICE solving. But while it cannot determine the P50 with high fidelity, it is it captures the major trends and relationship pertinent to solving the uncertainty, and thus the P ratios (e.g. P90/P50 or P95/P50). For, as explained in LINK, yield uncertainty is dominated by the major sources of uncertainty and their interrelation.

Specifically, SunSolve P90 computes the yield with the view-factor model for the optics, the Faiman for the thermal behaviour, and with multipliers for other losses. It is similar to the approach taken by PVsyst [ref] and SAM [ref]. The equations can be computed sufficiently quickly such that thousands of simulations can be solved in less than a minute.

2. Notes on inputs

A table of user inputs is available here. It lists whether the inputs can be attributed with uncertainty distributions.

Most of the variables in these equations are user inputs.

Many of the variables can also be assigned with an uncertainty distribution. If a distribution is not permitted, it is either because (i) the related error would be negligible in any reasonable yield simulation (e.g., the module dimensions), or (ii) the uncertainty could just as well be applied to another very similar variable.

Where an input has a distribution, its value is determined at the beginning of the simulation, year and time step, as described in LINK. That is, the randomness introduced into the variables is applied before the calculations are made with the equations that follow—with one exception. The occasion when the distribution is applied to an input after it has been used in a calculation, is in the first step (calculating fD from GHI), as will be explained.

Many variables are assigned the symbol f or k. These are modifying factors and loss factors, respectively, where we multiply results by either f or (1 – k). The default value for f is usually 1, and the default value for k is usually 0.

3. Inputs in weather file (and other time-dependent inputs)

The user loads a weather file containing 8760 timesteps, one step for each hour of a 365-day year. This file would typically be a TMY file.¹

Each timestep in the weather file contains

• solar position — i.e., the solar zenith θs and solar azimuth φs;

• the GHI and DHI; and

• the wind speed w and ambient temperature Ta.

There are also some optional time-dependent inputs that can be included in the file. If they are not included, their value will be assigned their ‘fallback’ values.

• module tilt β,

• albedo α,

• spectral correction fλ, and

• structural shading kS and transmission factor kT for bifacial systems.

4. Solar irradiance

The calculations begin by finding the diffuse fraction,

$$f_D = \frac{DHI}{GHI}$$

and then introducing the error into GHI and the diffuse fraction $f_D$ (i.e., applying the uncertainty distribution). This is the only occasion where the error is applied to a variable after it has been used.

The DNI is then given by

$$DNI = \frac{GHI \cdot (1 - f_D)}{\cos(\theta_s)}$$

Next, we find the extraterrestrial irradiance on a perpendicular plane,

$$ENI = GSC \cdot f_{ES}(t)$$

where GSC is the solar constant (assumed 1361.1 W/m2 [Gue18]) and $f_{ES}(t)$ is the earth–sun correction factor [ref]; and the extraterrestrial irradiance on the horizontal plane,

$$EHI = ENI \cdot \cos(\theta_s)$$

so that we can determine the sky clearness index,

$$K_t = \frac{GHI}{EHI}$$

and Hay’s sky clearness parameter (sometimes referred to as the anisotropy index),

$$K_b = \frac{DNI}{ENI} \cdot f_C$$

Here, we have introduced the correction factor, $f_C$. Since the Hay–Davies sky irradiance model assigns $K_b$ to be the fraction of diffuse light to be treated as the circumsolar fraction,² $f_C$ represents the error in the expected circumsolar fraction.

Finally, we can calculate the isotropic irradiance to the horizontal IHI, which is the diffuse light excluding the circumsolar fraction,

$$IHI = f_D \cdot GHI \cdot (1 - K_b)$$

and the beam irradiance to the horizontal BHI,

$$BHI = GHI - IHI$$

5. Module calculations

The module area $A_m$ is calculated from the module length $L_m$ and width $W_m$,

$$A_m = L_m \cdot W_m$$

It is also useful to determine difference in azimuth between the module orientation $\varphi_m$ and the solar azimuth $\varphi_s$,

$$\Delta\varphi = \varphi_m - \varphi_s$$

The module orientation (i.e., the direction the module faces) is defined in the same way as the solar azimuth, being the clockwise angle from due north.

6. Irradiance incident to module

We introduce a useful parameter related to incident angle and direct shading,

$$h_s = \cos(\beta) + \tan(\theta_s) \cdot \sin(\beta) \cdot \cos(\Delta\varphi)$$

and the distance from the top of the module to the shadow of the direct light,

$$L_s = P \cdot h_s$$

We then find the incident angle to the front side of the module,

$$\theta_i = \cos(\theta_s) \cdot h_s$$

If this angle is between $\pi/2$ and $\pi$ then the direct light is incident to the rear side with an angle of $\pi - \theta_i$.

We next calculate the impact of shading from neighbouring modules (assuming an infinitely large field with no lateral spacing). The shaded fraction of the beam on the front side is

$$f_{\text{shaded},F} = \begin{cases} 0 & L_s \leq 0 \\ L_s / L_m & 0 < L_s < L_m \\ 1 & L_s \geq L_m \end{cases}$$

and the shaded fraction on the rear is

$$f_{\text{shaded},R} = \begin{cases} 0 & L_s \geq 0 \\ -L_s / L_m & -L_m < L_s < 0 \\ 1 & L_s \leq -L_m \end{cases}$$

and thus, the beam irradiance that falls on the front without reflecting from the ground is³

$$\Phi_{F,B} = f_{\text{shaded},F} \cdot h_s \cdot BHI \cdot f_{\Phi,B}^F$$

and the beam irradiance that intercepts the rear without reflecting from the ground is

$$\Phi_{R,B} = f_{\text{shaded},R} \cdot h_s \cdot BHI \cdot f_{\Phi,B}^R$$

where $f_{\Phi,B}^F$ and $f_{\Phi,B}^R \geq 1$; these correction factors account for the increase in beam irradiance due to reflections from other modules or objects in the scene (but not the ground). They default to unity. Like many of the optical correction factors applied here, they are not needed in ray tracing programs.

The fraction of isotropic irradiance incident to the module (without reflecting from the ground) is

$$f_{\text{iso}} = \frac{L_m - \sqrt{L_m^2 - 2L_m P \cos(\beta) + P^2} + P}{2L_m}$$

where one uses $\beta$ to compute $f_{F,\text{iso}}$ and $(\pi - \beta)$ to compute $f_{R,\text{iso}}$. Hence, the isotropic irradiance on the front is

$$\Phi_{F,I} = IHI \cdot f_{F,\text{iso}} \cdot f_{\Phi,I}^F$$

and on rear is

$$\Phi_{R,I} = IHI \cdot f_{R,\text{iso}} \cdot f_{\Phi,I}^R$$

where $f_{\Phi,I}^F$ and $f_{\Phi,I}^R$ are correction factors to account for the increase or decrease in beam irradiance due to reflections from other modules or objects in the scene (but not the ground), and shadings from near-field objects (excluding other modules).

Following the conventional view-factor approach, we also determine the fraction of beam and isotropic irradiance that reflects from the ground onto the front and rear surface of the module. These fractions are $f_{B,G}^F$, $f_{I,G}^F$, $f_{B,G}^R$ and $f_{I,G}^R$ and they assume unity albedo, no structural shading, and no transmission through or between the modules.

We can then determine the irradiance that reflects from the ground onto the front,

$$\Phi_{F,G} = \rho \cdot (BHI \cdot f_{B,G}^F + IHI \cdot f_{I,G}^F)$$

and onto the rear,

$$\Phi_{R,G} = \rho \cdot (BHI \cdot f_{B,G}^R + IHI \cdot f_{I,G}^R)$$

Hence, the total irradiance on the front of the module would be the sum of the three components,

$$\Phi_F = \Phi_{F,B} + \Phi_{F,I} + \Phi_{F,G}$$

and an equivalent equation exists for the rear total irradiance.

7. Effective irradiance incident to the modules

We now apply three correction factors to the irradiance to account for soiling $k_\sigma$, spectral effects $f_\lambda$ and incident angle IAM,

$$\Phi_{F,\text{eff}} = f_\lambda \cdot (1 - k_\sigma) \cdot (IAM(\theta_i) \cdot \Phi_{F,B} + IAM_s \cdot (\Phi_{F,I} + \Phi_{F,G}))$$

In this formula, the same spectral correction and soiling loss is applied to all front irradiance, whereas a different IAM loss is applied to the components. Specifically, an IAM is computed for the incident angle $\theta_i$ of the beam, whereas a modifier for scattered light of $IAM_s = 0.97$ is applied for the isotropic and ground-reflected light.

An equivalent equation is applied for the rear irradiance, which contains the

$$\Phi_{R,\text{eff}} = f_\lambda \cdot (1 - k_R) \cdot (IAM(\theta_i) \cdot \Phi_{R,B} + IAM_s \cdot (\Phi_{R,I} + \Phi_{R,G}))$$

noting that the same loss factors have been applied except for the soiling loss (i.e., a separate soiling is applied for front and rear irradiance).

The $IAM(\theta_i)$ can be loaded from a PAN file and interpolated, but we believe it is more accurately fitted by modelling. In either case, it can be fitted with a cosh function without much trouble.

The total effective irradiance incident to the module then becomes:

$$\Phi_{\text{eff}} = f_E \cdot (\Phi_{F,\text{eff}} + \Phi_{R,\text{eff}})$$

Where we have included one extra irradiance multiplier $f_E$, whose real purpose is to allow the user to add uncertainty to the optical calculation with one single variable. This might arise from optimising the tracker angle for diffuse-sky conditions.

8. Module power and thermal loss

We determine a single U-value following the Faiman model

$$U = U_c + w \cdot U_v$$

then module efficiency at STC,

$$\eta_{\text{STC}} = \frac{P_{m,\text{STC}}}{A_m \cdot 1000\,\text{W/m}^2}$$

and a useful parameter that we’ll use again,

$$D = \frac{U \cdot \eta_{\text{eff}}}{\eta_{\text{STC}} - B}$$

where $B$ is the thermal coefficient of power. (Note that a value of $B = 0.003\,\text{K}^{-1}$ means that output power from a module increases by 0.3% per K.)

We next solve the module power at the operational temperature,

$$P_{m,\text{noMM}} = A_m \cdot U \cdot (1 - B) \cdot (T_a - T_{\text{STC}}) - \frac{\alpha \cdot \Phi_{\text{eff}}}{B \cdot D}$$

where $\alpha$ is the absorptance (and typically assumed equal to 0.9).

We then account for cell-to-cell mismatch by applying the user-defined mismatch factor,

$$P_m = P_{m,\text{noMM}} \cdot f_{\text{MC}}$$

Although it’s not necessary for SunSolve P90, we can also calculate the operating temperature of the module,

$$T_m = T_{\text{ref}} + \frac{1}{B} \left( 1 - \frac{P_m}{A_m \cdot \Phi_{\text{eff}} \cdot \eta_{\text{STC}}} \right)$$

and the DC module efficiency accounting for temperature,

$$\eta_m = \eta_{\text{STC}} \cdot (1 - B \cdot (T_m - T_{\text{STC}})) \cdot f_{\text{MC}}$$

9. String losses

We calculate the DC power of the string accounting for the number of modules in a string $N_m$, as well as loss factors for module-to-module mismatch $k_{\text{MM}}$, string wiring losses $k_{\text{WS}}$, and max-power tracking losses (module not held at MPP) $k_{\text{MPT}}$.

$$P_{s,\text{DC}} = N_m \cdot P_m \cdot (1 - k_{\text{WS}}) \cdot (1 - k_{\text{MM}}) \cdot (1 - k_{\text{MPT}})$$

10. Inverter losses

We then compute the power entering the inverter, accounting for the number of strings per inverter, the inverter wiring losses $k_{\text{WI}}$, the string-to-string mismatch $k_{\text{MS}}$, and inverter clipping $P_{\text{clip}}$,

$$P_{i,\text{DC}} = \min(P_{\text{clip}}, N_s \cdot P_{s,\text{DC}} \cdot (1 - k_{\text{WI}}) \cdot (1 - k_{\text{MS}}))$$

before converting the power from DC into AC, accounting for the inverter efficiency $\eta_I$,

$$P_{i,\text{AC}} = P_{i,\text{DC}} \cdot \eta_I$$

11. Field losses

We calculate the AC power from the field $P_{f,\text{AC}}$, accounting for the number of inverters $N_i$ and inverter-to-inverter mismatch $k_{\text{MI}}$,

$$P_{f,\text{AC}} = N_i \cdot P_{i,\text{AC}} \cdot (1 - k_{\text{MI}})$$

12. Annual yield

The optimal annual yield is the sum of all the hourly datapoints from the field for each year without the inclusion of annual yield losses:

$$Y'_y = \Delta t \sum_{\text{year}} P_{f,\text{AC}}$$

where $\Delta t$ is the time step (assumed to be hourly) and $Y'_y$ is in Wh.

13. Annual yield losses

The following annual yield losses are then applied: DC health $f_{\text{DCH}}$, availability $f_{\text{avail}}$, and curtailment $f_{\text{curt}}$, as well as an annual degradation rate $d$.⁴

For each progressive year $y$, we follow the typical linear degradation function

$$d_y = (y - 0.5) \cdot d$$

None of these inputs include hour-to-hour variability because they’re applied at the annual level.

$$Y_y = Y'_y \cdot (1 - k_{\text{DCH}}) \cdot f_{\text{avail}} \cdot (1 - k_{\text{curt}}) \cdot (1 - d_y)$$

14. Lifetime yield

The lifetime yield is

$$Y_L = \sum_y Y_y$$

Footnotes

1. No interpolation of the hourly data is performed (e.g., to calculate sunrise and sunset). And since the position of the sun is loaded at each time step, it is irrelevant whether the time steps represent the start of hour, middle of hour, or end of hour. ↩

2. Thus, the more light that passes through the atmosphere, the clearer the day, and the smaller the isotropic contribution to the DHI. ↩

3. This works because cosine of the incident angle is cosθi = cosθs∙hs, and the irradiance incident to the module Φ = BNI ∙ cosθi, where BNI = BHI / cosθs. Hence Φ = BHI ∙ hs. ↩

4. We do not follow a compound degradation formula, dy=i=1ydi. ↩