Irradiance modelling for individual cells of shaded solar photovoltaic arrays

Developments in Photovoltaic (PV) design software have progressed to modelling the string or even the module as the smallest system unit but current methods lack computational eﬃciency to fully consider cell mismatch eﬀects due to partial shading. This paper presents a more eﬃcient shading loss algorithm which generates an irradiance map of the array for each time step for individual cells or cell portions. Irradiance losses are calculated from both near and far obstructions which might cause shading of both beam and diﬀuse irradiance in a three-dimensional reference ﬁeld. The irradiance map output from this model could be used to calculate the performance of each solar cell individually as part of an overarching energy yield model. A validation demonstrates the calculation of shading losses due to a chimney with less than one percent error when compared with measured values.


Introduction
Developers of solar photovoltaic (PV) systems would benefit from more accurate prediction of energy yield and internal rate of return (IRR). A reduction in prediction uncertainty to within given confidence limits would help secure lower cost financing, contributing to the drive for grid parity in the solar industry. There is a lack of consensus in the industry regarding how much separation should be left between PV arrays and near shading objects such as chimneys and dormer windows. This manifests in anecdotal reports of poorly designed systems with modules heavily shaded by obstructions, where closer attention to design would have made significant improvements to energy yield. Mismatch effects of shading are not normally considered in system modelling because the computation time would be too great.
Shading can be the most detrimental factor on performance for a domestic system. The impact of shading on performance varies depending on the electrical series and parallel arrangement of cells within a module and modules within an installed array. Whilst many approaches to shading analysis have been proposed, computational efficiency is not reported despite being of high importance when incorporating shading algorithms into an overall energy yield model. The lack of consideration of the non-linear impacts of shading on smaller systems for example means that the shading loss is significantly underestimated, especially from supposedly small obstacles such as antennas or chimneys. As an example, the system shown in Fig. 1 illustrates the case where the installer may have attested a shade loss factor close to unity under UK microgeneration guidelines (Microgeneration Certification Scheme, 2013), i.e. negligible, but the performance of the system is severely compromised due to the non-linear cell mismatch effects. An effective shading sub-model therefore needs to give feedback to inform decisions of array layout in the proximity of obstructions but must not rely on high power computing.
The algorithms to calculate shading losses within an overall PV system energy yield model can be divided into two main sub-models: (a) The shaded irradiance sub-model -which calculates irradiance incident on the cells, using spatial location data for shading objects. (b) The array electrical sub-model -which calculates current & voltage for each string, taking mismatch into consideration using cell irradiance calculated in the shaded irradiance sub-model (Bishop, 1988;Quaschning and Hanitsch, 1996;Overstraeten and Mertens, 1986;Liu et al., 2011).
This paper is concerned primarily with the shaded irradiance calculation, the output of which can be interfaced with any electrical mismatch model. Shaded irradiance models fall into two main categories, those which view Models in category (a) commonly use rendering to generate a three dimensional view of a building or district, with shading used to indicate zones of varying irradiation. (Mardaljevic and Rylatt, 2003;Compagnon, 2004;Levinson et al., 2009;Lukač et al., 2013). A key challenge of this approach is the computation time to model irradiance for each surface segment and for each hourly sun position. A logical optimisation is to bin sun positions into zones of similar irradiation, for example from 4000 hourly sun positions above the horizon into 250 bins (Mardaljevic and Rylatt, 2003 simulations at neighbourhood level using building coordinates from GIS or LiDAR coordinates. Category (b) models use the concept of the sunpath diagram to analyse whether the surface or portion of it is shaded for a given time-step (Quaschning and Hanitsch, 1998;Robinson and Stone, 2004;Cheung and Chung, 2007;Cellura et al., 2012;Drif et al., 2008;Crocker and Sullivan, 2013).
PV design tools such as PV*Sol and PVsyst model shading in PV arrays with strings or modules as the smallest unit (Mermoud and Lejeune, 2010;Valentin, 2013). Modelling the energy yield cell-by-cell is considered too slow with currently available methods (MacAlpine et al., 2012), simulation time has been identified as a key factor for users, after complexity, cost and CAD integration (Horvat and Dubois, 2012).
The key challenge of the shaded irradiance sub-model is therefore to efficiently generate a matrix of irradiance values for each portion of the PV array for each time-step. The array portion could be a module, sub-module, cell or even sub-cell. In this work the array portions are defined as individual cells. Generation of this irradiance map is potentially a complex task with significant computation time, which could be problematic to incorporate into an energy yield model. The main challenge of this work is to identify a method which achieves the required accuracy whilst operating a minimal total number of calculations.
This paper presents a model which is computationally light enough to operate on a typical laptop computer but calculates irradiance losses due to shading with an accuracy of 1% or less. The model combines elements of the specifically developed sky-patch (Tregenza, 1987) concept of irradiance distributions and the polygon point containment method (Quaschning and Hanitsch, 1998;Mortenson, 1990) for shape analysis to create a new efficient shadow calculation algorithm.

Overview of the model
This shading sub-model consists of two main phases of operation: The geometry preparation and then the time series simulation.

Phase 1, Horizon analysis and sky-dome preparation
The horizon as 'seen' by each cell in the array remains fixed for the duration of the simulation, therefore it can be digitised as a dome of sky-patches and stored in memory. The only exception to this rule is for non-mature trees, which can be accommodated by repeating the simulation for each year of tree growth. Tree growth rate is not included in the methods presented.
Each sky-patch has a Boolean value to indicate whether the cell is shaded or unshaded when the sun position is within it. This preparation is done only once before the time series simulation, so that computation in each loop of the time series is minimised. Two main types of sky-dome configurations were considered. The first type has patches equally spaced in elevation and azimuth but with varying size, as shown in Fig. 2, the second type is described in Fig. 9, Section 5.
The flowchart in Fig. 3 presents an overview of the preparation phase of the model. In the first step the horizon line surveyed for a single datum point on the array must be translated to a horizon line for every cell in the array. An array of Boolean sky-patches is then created for each cell. Each sky-patch is given a 'shaded' Boolean value if it is below the horizon line for that cell or the horizon line passes through it. The diffuse irradiance shading loss factor can then be calculated for each cell by counting the number of shaded sky-patches and calculating the cosine losses.

Phase 2, Time series simulation
For each time-step in the simulation, the irradiance calculation is a simple Boolean check if the cell is shaded for that sun position in the predefined sky-dome. The beam and diffuse loss factors are applied to the beam and diffuse irradiance respectively.
The output from the model is an array of irradiance values for each cell and for each time-step, this array can then be input into a cell performance and interaction model as part of a detailed array energy yield simulation.
3. Detailed description of the model 3.1. Preparation for simulation: surveying horizon coordinates PV system designers are familiar with recording horizon data as (azimuth, elevation) coordinates onto a two dimensional sunpath diagram as found in the user interface of various energy yield tools (PV GIS, 2013) and as shown in Fig. 4. A similar approach was taken for the user interface of this model.
To model shading on individual cells requires the surveyed shading objects and the array to be described in a three dimensional space. For simplicity of calculation, the coordinates of each point on the surveyed shading objects are initially described by a spherical vector [a H , e H , d], where a H , e H and d represent the azimuth angle, elevation angle and distance from a datum point on the array to a given point on the horizon as shown in Fig. 5.
The datum position of the array for the model is defined as the centre of the middle cell of the lowest row of cells as shown in Fig. 5 so the horizon must be surveyed from this point or translated to it.
Use of the spherical coordinate system means that input data has some compatibility with the simple 2D horizon plots on a sun-path diagram consisting of an azimuth angle and elevation for each point as shown in Fig. 4. Each 2D point is converted to this 3D frame by adding a distance to each point. Cartesian coordinates from 3D CAD packages could also be imported by simply converting to spherical coordinates.

Translating the horizon coordinates from the datum cell to all other cells in the PV array
Once the user has input the horizon data surveyed from the datum position of the array, the datum horizon line can be translated to a horizon line as 'seen' by every cell in the array using the methodology below. For increased resolution, the same approach could be applied for multiple points within cells, or for increased speed, a sub module with one bypass diode could be considered the smallest unit in the array. For this paper the solar cell is defined as the smallest physical and electrical unit in the array.
All geometric formulae in the model use right handed Cartesian or spherical coordinate systems as described in ISO standard 80000-2 Pt 16 (ISO, 2009) , where: the x-axis is positive southwards; the y-axis is positive eastwards and the z-axis is positive upwards, as shown in the reference frame in Fig. 5. The spherical coordinate system in the user interface follows the same convention as other PV software (PVGIS, Meteonorm, PVsyst, PV*Sol) but differs from ISO80000-2 so the horizon coordinates must first be converted as follows: Azimuth aðISOÞ ¼ ÀAzimuth aðGUIÞ ð 2Þ Each horizon point is then converted from polar to Cartesian co-ordinates using the following formulae: The centre of each cell in the array is identified with coordinates [i, j] as shown in Fig. 6. The spacing between cells centres is defined as For square, octagonal, or round solar cells: Fig. 4. Drawing of the horizon as a polygon, viewed from sun towards array.
For the purpose of this description, the additional spacing between modules and the module frame is omitted for clarity. However in practice for larger arrays these would be included.
The Cartesian vectors describing the position of each cell relative to its nearest neighbour for an array plane with azimuth of zero and tilt angle e a are calculated as: Eq. (5) assumes crystalline PV Cells, which are usually square or non-regular octagonal in shape.
The position of any point on the horizon relative to any cell in the array can then be calculated from the datum horizon using the following translation formulae: Each horizon point is then converted from Cartesian to polar co-ordinates using the following formulae: The calculations described in Eqs.
(1)-(7) ultimately calculate a small translation of the coordinates for a given horizon point from the datum cell to any other cell in the array as shown in Fig. 7, which shows the elevation angle of the horizon from the datum cell e H1 translated to the elevation angle of the same horizon point from another cell e H2 .

Digitising the horizon line to an array of sky-patches
The sky is described as a dome of sky-patches, through elevation of 0-90°and azimuth 0-360° (Fig. 9). Note that for computational simplicity, the sky-patch's azimuth is from 0°to 360°with south as 180° (Fig. 8). The user interface for this model uses the À180°to 180°ISO solar convention with south as 0°in common with PV design packages, conversion is simply.
Two distinct sky-dome configurations were compared for use in the model. The first approach has an arrangement of patches which is symmetrical in both axes, as shown in Fig. 2. This approach allows any sky-patch resolution to be used according to the computing power available and desired accuracy. The number of sky-patches n is defined as For example, a sky-dome with 1-degree 2 resolution would have 324,000 sky patches. A sky-dome with 10-degree resolution is shown in Fig. 2. This sky-dome has the disadvantage of uneven weighting between irradiation in lower and upper portions of sky, so it was rejected as being unsuitable.
The sky dome definition used in this paper has patches of equal size and equally spaced in elevation but differing azimuthal spacing for each row of patches, similar to the BRE/CIE sky-dome (Fig. 9). The advantages of this approach are that the patches are equally sized, so have equal weighting, and it has potential to use the Perez anisotropic diffuse irradiance model (Robinson and Stone, 2004; Perez et al., 1993) which uses the BRE/CIE sky-dome arrangement. The limitation of the pure BRE/CIE skydome is that it would not account for the times when the sun position aligns with a gap between circular patches. This source of bias error is resolved by using an amended version of the BRE/CIE sky-dome, with quadrilateral instead of round patches but with the same centre points as shown in Fig. 10. Each patch has the same planar angle in azimuth and elevation axes and therefore the same solid angle. The advantage of quadrilateral patch shape over for example hexagonal, is that each patch is described by only four angles, the minimum and maximum azimuth and elevation, allowing for simple and fast computation within the time series loop.
The coordinates of sky-patches are calculated using the values in Table 1, The BRE/CIE sky-dome has a total of 151 patches, arranged in horizontal bands, with the band centres spaced with a 12°elevation angle.
Sky-domes with 4, 3, 2 and 1 degree band centre angular spacings (Table 2) are also compared in the validation.
Since the horizon does not change during the hourly simulation, the sky-dome for each cell needs to be set up once only. For each sky-patch, the model must check whether: (a) The sky-patch is in front of the array. (b) The cell has an unobstructed line of sight to the sky-patch.
If both these are true, then the sky-patch can 'have beam irradiance' -stored as a Boolean value in an array which is then called for each time step. The array does not need to contain actual irradiance values since these are called later during time series simulation. This approach minimises the calculation required for each step during the time series simulation, since for a given hour only a single memory location is accessed to verify if beam irradiance is present for the sky-patch in question. Fig. 11a shows a view of a house shaded by a tree which will be surveyed and converted to a sky-dome based shading map.
Digitisation from a horizon of vectors to an array of Boolean sky-patches is achieved using the point in polygon containment test (PPCT). For the purpose of the PPCT the sky-dome is projected onto a two dimensioned view of the dome, in the same way the Earth is plotted on maps in the Mercator projection. Fig. 11b shows sky-patch A below the horizon line and sky-patch B above the horizon line. This two dimensional Fig. 9. Drawing showing the BRE/CIE sky-dome of equally sized skypatches (Tregenza, 1987). Fig. 10. Amended BRE/CIE sky-dome with quadrilateral patches. Table 1 Spacing of patches in the BRE/CIE sky-dome, which has 12 degree angular spacing between bands.

Band index
Elevation angle of band centre (°)

Number of patches in band
Azimuth angle between patches in band (°) 1 6 32 11.25 2 1 8 3 0 1 2 3 30 28 12.85714 4 4 2 2 4 1 5 5 5 4 1 8 2 0 6 6 6 1 2 3 0 7 7 8 6 6 0 8 90 1 360 view of the sky-dome can be considered to have an extreme left and right side. The PPCT requires a horizontal test line to be created from each point to be examined (A&B) to a virtual point on the extreme left of the sky-dome with the same elevation angle (in unshaded space) as shown in Fig. 11b. This test line is checked against each portion of the horizon line to identify if they intersect. The number of intersections is counted for each sky-patch. If there is an odd number of intersections, the point is below/within the horizon (as with Point 'A'). If the number of intersections is zero or even, the point is above/outside the horizon (as with Point 'B') Quaschning and Hanitsch, 1998;Mortenson, 1990. For this model to work for all cases, the extreme left point on the horizon must have zero elevation. If the origin of the horizon line is not at zero elevation then an additional point with zero elevation must be added, to the extreme left of the view. A vertical test line could also be used, but would require the same total number of computations to check against each section of the polygon.
Whilst simpler methods could be used to check if points are above or below the horizon, the PPCT has been shown to work reliably for multiple polygons with any shape and any number of points. Taking the example of a horizon with overhanging sections as shown in Fig. 12, the model must recognise that patch A would have beam irradiance despite being below part of the horizon.
Each sky-patch is also assigned a Boolean value for 'in front of array'. Storing this information separately from the 'Unshaded' Boolean value allows the flexibility to extend the model for bifacial modules. To identify if a patch is behind the array it is necessary to define the formula for the arc across the sky-dome where an infinite PV array would intersect with it (Fig. 13).
The elevation angle of a point on an arc where an infinitely large array would intersect the skydome e B , is defined as where a P is the azimuth angle of the sky patch and e A is the tilt angle of the PV array. In this scenario, the sky patch contributes irradiance only if: 4. Time series simulation

Diffuse irradiance loss calculation
For fixed array systems both the horizon and array are static. If an isotropic sky is assumed, a simple loss factor for diffuse irradiance due to shading can be calculated once for each cell based on the number of shaded sky-patches. This is then applied to the diffuse irradiance for each time-step.
where L D,I,J = the diffuse loss factor for cell with position ij in the array and P SP st ¼ Unshaded is the number of unshaded sky-patches as visible from cell I, J. P SP st is the total number of sky-patches.
The model described here assumes an isotropic sky, but the BRE/CIE sky-dome definition is the same as that used in the Perez anisotropic diffuse irradiance model (Perez et al., 1993), so these 2 models could be combined to form an anisotropic shading model by incorporating the weighting factors to the sum in Eq. (12).

Beam irradiance loss calculation
For each time step, the model must: (a) Get beam and diffuse irradiance from stored values or in plane irradiance model.  (d) Check whether the sky-patch is in front of the array (using previously stored Boolean values for each sky-patch). (e) Check if the sky-patch is unshaded (using previously stored Boolean values for each sky-patch).
The sky-patch indexes s & t are calculated simply as: where R X = sky-dome resolution in the X (azimuth) axis R Y = sky-dome resolution in the Y (elevation) axis.

Validation
The model was validated against a PV array of 7 modules with array azimuth 0°(south) and tilt angle 35°. Two scenarios were used: (1) Chimney to the side of the array (Fig. 14a).
(2) Chimney to the front of the array (Fig. 14b).
The resulting representation of the horizon in skypatches is shown in Fig. 15. The chimney appears asymmetric due to the increasing width of each patch with increasing elevation.   The shaded/unshaded property of each cell was recorded in a spreadsheet, as shown in Fig. 16. Each square represents a quarter cell in this image to increase the validation resolution.
This data was then compared with the irradiance map for that time as generated by the model.
There are two possible ways to calculate the error of the model (1) The difference in the number of shaded cells between the model and the rooftop measurements.
(2) The number of cells which are incorrectly modelled as shaded plus the number of cells which were incorrectly modelled as unshaded.
The latter approach is important because even if the absolute number of cells shaded is correct, if they are in incorrect positions this may impact on the modelling of electrical mismatch losses between cells. Graphs Figs. 17a and 17b show how the errors vary with different shadow shapes but were at or below 1% for both scenarios for all times of day.
The choice of patch resolution is necessarily a compromise between accuracy and computation time. Fig. 18 indicates that below 2 degrees, smaller patch sizes make only marginal improvements in accuracy but cause exponential increases in computation time. The unit of computing time for a one degree patch resolution was 2.9 ls to run 8760 hourly simulations on a typical laptop with an Intel i5 processor running at 2.7 GHz.

Conclusions
Shading was identified as a critical area for improvement of the accuracy of PV energy yield modelling. Energy yield models for system designers need to generate fast results for designers so they can quickly optimise designs to maximise financial yields.
The proposed model uses established methods of the sky-dome of patches and the polygon point containment method in a new configuration to model shading more accurately and more efficiently. It generates fast results by reducing the computation in each time step to an absolute minimum yet generates a map of irradiance on each cell for each time step.
The Irradiance map output of the model is designed for ease of connection to cell performance sub models as part of an overall energy yield model.
The model allows for any number of shading objects of any shape since the datum horizon line can contain any number of points. These could easily be generated from CAD drawings for near objects and from GIS data for far objects.
The modelling of diffuse irradiance here assumes an isotropic sky, but use of the sky-dome concept enables anisotropic diffuse models to be integrated in further work, in particular the BRE/CIE sky-dome is the same definition as used in the Perez anisotropic model. This model was designed for single tilted arrays, where cells are arranged in a uniform spacing in both directions, in that the module layout is symmetrical in both directions. Capability to model asymmetric arrays could be added by storing an array of actual cell position coordinates, rather than assuming regular spacing.
Identification of shaded cells was achieved with an error of less than 1%. Further improvements are possible with high accuracy surveying methods and above average computing power.

Acknowledgements
This work has been supported by a joint UK-India initiative in solar energy through a joint project 'Stability and