Minimising the Levelised Cost of Electricity for Bifacial Solar Panel Arrays using Bayesian Optimisation

Bifacial solar module technology is a quickly growing market in the photovoltaics (PV) sector. By utilising light impinging on both, front and back sides of the module, actual limitations of conventional monofacial solar modules can be overcome at almost no additional costs. Optimising large-scale bifacial solar power plants with regard to minimum levelised cost of electricity (LCOE), however, is challenging due to the vast amount of free parameters such as module inclination angle and distance, module and land costs, character of the surroundings, weather conditions and geographic position. We present a detailed illumination model for bifacial PV modules in a large PV field and calculate the annual energy yield exemplary for two locations with different climates. By applying the Bayesian optimisation algorithm we determine the global minimum of the LCOE for bifacial and monofacial PV fields at these two exemplary locations considering land costs in the model. We find that currently established design guidelines for mono- and bifacial solar farms often do not yield the minimum LCOE. Our algorithm finds solar panel configurations yielding up to 23 % lower LCOE compared to the established configuration with the module tilt angle equal to the latitude and the module distance chosen such that no mutual shading of neighboring solar panels occurs at winter solstice. Our algorithm enables the user to extract clear design guidelines for mono- and bifacial large-scale solar power plants for most regions on Earth and further accelerates the development of competitively viable photovoltaic systems.


Introduction
The record power conversion efficiency (PCE) of monofacial silicon solar cells -currently the dominant solar-cell technology -is 26.7% 1 and approaches the physical limit of around 29.4%, which was calculated by Richter et al. 2 Photovoltaic (PV) systems consisting of bifacial solar modules can generate a significantly higher annual energy yield (EY) than systems using conventional monofacial PV modules, because bifacial solar modules not only utilize light impinging onto their front, but also illumination onto their rear side. 3,4 Furthermore, advanced solar-cell concepts such as PERC, PERT, PERL (passivated emitter rear contact/totally-diffused/locally-diffused) and IBC (interdigitated back contact) can easily be manufactured as bifacial solar cells. 5 Kopecek and Libal see bifacial solar cells as the concept with the 'highest potential to increase the output power of PV systems at the lowest additional cost'. 3 Indeed, the bifacial solar cell market has been gathering pace for a couple of years and several major PV companies, such as Sanyo, 6 Yingli, 7 PVG solutions, bSolar/SolAround, 8 and Trina Solar 9 introduced bifacial modules. The tenth edition of the International Technology Roadmap for Photovoltaics (ITRPV) predicts a global market share of more than 50% for bifacial modules in 2029. 10 Large-scale bifacial PV power plants already have been realised and showed a higher energy yield than their monofacial counterparts. 11 The levelized cost of electricity (LCOE) is a very relevant economic metric of a solar power plant. 12 The performance of bifacial solar modules is heavily affected by their surroundings, because they can accept light from almost every direction. Hence, a vast amount of parameters influence the resulting LCOE, for example the module and land costs, module distance and inclination angle, albedo of the ground, geographical position and the weather conditions at the location of the solar farm. Liang et al. recently identified comprehensive simulation models for energy yield analysis as one of the key enabling factors. 4 As an example, we briefly discuss how only two free parameters -land cost and module distance -affect the resulting LCOE, which makes it challenging to identify the sweet spot yielding a minimum LCOE: If two rows of tilted solar modules are installed close to each other, many modules can be installed per area. However, at too small distances shadowing will limit the rear side irradiance and consequently the total energy yield. 13,14 In contrast, putting the rows of modules far apart from each other maximizes the irradiance at the rear side and the energy yield per module. The number of modules installed per area, however, is lower and the overall energy yield of the solar farm decreases. The module inclination angle is a third free parameter, closely connected to the two aforementioned module distance and land cost, and obviously affects shadowing of neighboring solar panel rows and hence energy yield and LCOE of a bifacial solar farm, too.
Historically, the module inclination angle was usually set to the geographical latitude of the solar farm location, and the module distance was either set to a fixed value based on experience 15 or to the minimum module distance without mutual shadowing on the day of winter solstice at 9 am 16 or noon. 17 However, it has turned out that these rule-of-thumb estimates often do not lead to a minimised LCOE. 18 One reason is that these models did not consider the cost of land. Recently Patel et al. considered land costs when optimising bifacial solar farms. 16 However, also in this study the module distance and inclination angle were preset according to above mentioned winter solstice rule. Considering the enormous market growth of bifacial solar cell technology, finding the optimum configuration yielding minimum LCOE is highly desired. With the PV system costs in $ per Watt peak (Wp), land costs in $ per area and the geographic location of the solar farm as known input variables, inversely finding the optimal geometrical configuration of a bifacial PV field is a computational challenging multi-dimensional optimisation task.
In this study, we apply a multi-parameter Bayesian optimisation in order to minimise the LCOE of large-scale bifacial solar power plants. We present a comprehensive illumination model for bifacial solar arrays and calculate the annual energy yield (EY) based on TMY3 (Typical Meteorological Year 3) data for two exemplary locations near Dallas and Seattle. We calculate optimal module inclination angles and module distances yielding minimal LCOE for various module to land cost ratios. We find that our calculated optima strongly depend on both the module to land cost ratio and the geographical location. We conclude that currently used rule-of-thumb estimates for optimal module distance and tilting angle must be reconsidered. Our method enables the user to extract clear design guidelines for mono-and bifacial large-scale solar power plants principally anywhere on Earth.

Illumination model
With the illumination model we calculate the irradiance onto a solar module, which is placed somewhere in a big PV-field. We assume this field to be so big that effects from its boundaries can be neglected. Further, we assume the modules to be homogeneous: we neglect effects from the module boundaries or module space in between the solar cells. Hence, we can treat this problem as 2-dimensional with periodic boundary conditions, as illustrated in Fig. 1. A similar approach was pursued for example by Marion et al. 19 In the current model we assume the solar modules to be completely black, which means they do not reflect any light which could reach another module.
The PV field is irradiated from direct sunlight under the Direct Normal Irradiance (DNI) 1 and the direction n S , which is determined by the solar azimuth φ S and the solar zenith θ S . The latter is connected to the solar altitude a S (the height above the ground) via a S = 90 • − θ S . Further, the PV field receives diffuse light from the sky, which is given as Direct Horizontal Irradiance (DHI). However, for calculating the total irradiance onto the module, also light reflected from the ground and shadowing by the other modules must be taken into account. Figure 1 shows the different components of light, which can reach the front of a PV module at point P m . The numbers 1.-4. correspond to the numbers in the figure -illumination on the sky is w.l.o.g. indicated for module #2 while illumination from the ground is indicated w.l.o.g. for module #5. 1 The irradiance or intensity is the radiant power a surface receives per area.
2 Table 1: The four irradiance components which constitute the illumination of a solar module in dependence of the position P m on the module as defined in Fig. 1, where s is the distance BP m . These components have to be considered for front and back sides -hence eight components in total. The numbers correspond to the numbers in Fig. 1 1. Direct sunlight hits the modules under the direction n S . It leads to the irradiance component I sky dir, f (s) = DNI · cos σ mS , where s is the distance between the lower end of the module B 2 and P m , s = B 2 P m , and σ mS is the angle between the module surface normal and the direct incident sunlight.
2. Diffuse skylight I sky diff, f (s) hits the module at P m from directions within the wedge determined by D 1 P m D 2 . Diffuse light does not only reach the module from directions within the xz-plane but from a spherical wedge, which is closely linked to the sky view factor as for example used by Calcabrini et al. 20 3. I gr.
dir, f (s) denotes direct sunlight that hits the module after it was reflected from the ground. 4. Finally, I gr.
diff, f (s) denotes diffuse skylight that hits the module after it was reflected from the ground.
All four components are summarized in Table 1. Table 2 denotes all parameters that are used as input to the model.
The total irradiance (or intensity) on front is given by and similar for the back side with a subscript b instead of f . In total, we hence consider eight illumination components on our module.
As noted above, the incident light is given as DNI and DHI. The nonuniform irradiance distribution on the module front and back surfaces has to be considered. 21,22 For the further treatment, it is therefore convenient to define unit-less geometrical distribution functions as for the components arising from direct sunlight and diffuse skylight, respectively. Here we omitted the superscripts "sky" and "gr.". The calculation of the components ι gr. dir, f (s) and ι gr. diff, f (s) requires the integration over geometrical distribution functions on the ground γ dir (x g ) and γ diff (x g ), where x g is the coordinate of the point P g on the ground. x' x z y N θ m s x g 1. 3.

2.
I. II. II. II. Figure 1: Illustrating the geometrical configuration of a (periodic) PV field and the illumination components, which reach each module on the front. The modules are labeled with #1-#5. At #1, the geometrical parameters h, , d and θ m are illustrated -d is the horizontal length of a unit cell. At #2, the two irradiance components illuminating the module from the sky at P m are indicated: 1. direct and 2. diffuse. Below #3, the I. direct and II. diffuse illumination of point P g on the ground are illustrated -here diffuse illumination origins from three angular intervals. On #5 the angular range of light reaching P m from the ground is indicated. It consists of 3. direct and 4. diffuse light being reflected from the ground. Components 1.-4. are summarized in Table 1. Here, we assume w.l.o.g. that the PV system is located on the northern hemisphere and oriented towards South.
In particular, we have where we omitted the subscripts "diff" and "dir". The coordinate x g (s, α), on which γ dir and γ diff are evaluated, is defined such that the angle between the line P g P m and the module normal n m is equal to α -the integration parameter. In Fig. 1 the fractions of the ground, which are illuminated by direct sunlight, are marked in orange. Figure 2 shows an example for illumination onto the ground: subfigure (a) illustrates the position of the solar modules #1 and #2. Subfigure (b) shows the geometrical distribution functions on the ground. γ diff is minimal below the module where the angle covered by the module is largest; and maximal at x , because here the ground sees least shadow from module #1.
Depending on the geometrical module parameters and the position of the Sun, the directly illuminated area (1) may lay completely within the unit cell as in the examples in Fig. 1 and Fig. 2, (2) it may extend from one unit cell into the next or (3) no direct light can reach the ground. The latter can occur when the module spacing d decreases or when the solar altitude a S is low. Figure 3 shows the eight geometrical distribution functions ι corresponding to the irradiance components hitting the PV module on its front and back sides. While the functions originating from the sky (a) are stronger on the front side, the components originating from the ground (b) are stronger on the back side. This can be understood by the opening angles: the opening angle towards the sky is larger on the front side, but the opening angle of the ground is larger at the back.
All the calculations presented in this work were performed with Python using numpy as numerical library for fast tensor operations.

Calculating the energy yield
We calculate the annual electrical energy yield EY by feeding the illumination model described in section 2 with irradiance data. To demonstrate the features of the model, we use TMY3 (Typical Meteorological Year 3) data for this work. TMY3 data is well suited to estimate the solar energy yield for thousands of different locations. 23 Amongst other parameters, the TMY3 data contain hourly DHI(t) and DNI(t) values. The overall EY is the sum of the energy yields and EY b with a subscript b instead of f . TMY3 data is available at the time stamps t i and ∆t is the time between two time stamps, which is typically 1 h for TMY3 data. η f and η b denote the power conversion efficiency for light impinging on the front and back sides of the solar module, respectively. For this work, we use η f = 0.20 and η b = 0.18, hence a bifaciality factor of 0.9. 5 By setting η f = η b = 1, eq. (4) delivers the annual radiant exposure on the front H f (and similarly H b ).
The ι-functions are evaluated on the positionŝ i ∈ P m , where P m = {s 1 , s 2 , . . . , s Nm } is the set of all considered positions along the module. In a conventional PV module, all cells are electrically connected in series and therefore the cell generating the lowest current limits the overall module current. To take this into account, we determineŝ i such that for all s ∈ P m . This means that the position on the module with the lowest irradiance, which is proportional to the solar cell current, determines the overall module performance. For high-end solar modules, the module performance might be higher depending on how bypass diodes are implemented. Therefore, our condition establishes a lower bound of the module performance under certain illumination conditions.
For the diffuse irradiation components on the module, the corresponding geometrical distribution functions ι diff (s) need to be calculated only once because we assume the incoming diffuse radiation to be isotropic for all timestamps. For the components arising from direct sunlight, also the geometrical distribution functions ι dir (s, t i ) are time-dependent, because they depend on the position of the Sun (θ S,i , φ S,i ), 2 which we calculate using the Python package Pysolar. 25

Results and Discussion
As an example, we discuss results for two locations with different climates:  Figure 9 shows climate diagrams for these two locations. We see that H generally increases with the module spacing. However, it is not economical to have a too large distance between the rows as we will see when considering the electricity cost in Section 4.
For Dallas, the optimal angle for monofacial modules, which only can utilize front illumination, is about 28°; it is mainly determined by direct sunlight. For back illumination, H increases significantly with the module inclination angle θ m : hardly any direct light reaches the module at the back, but contributions from diffuse sky and reflected from the ground increase with θ m . Increasing the module tilt further reduces the shaded area on the ground and therefore increases ground illumination. The optimal module tilt for a bifacial module is a compromise between the optimal tilt for the front and beneficial higher tilt angles for back contribution. Overall, the optimal module tilt for bifacial modules is significantly higher than for monofacial modules. Here it is about 36°.
Overall, the trends for Seattle are comparable to those for Dallas. However, we can identify differences: the overall radiant exposure is much lower because Seattle sees around 2170 annual Sun hours, compared to about 2850 h in Dallas. 27 Further, the optimal tilt for monofacial and bifacial modules is 32°and 44°, respectively, which is explained by the higher latitude of Seattle.    7 months (with higher elevation angles of the sun) benefits from lower tilt angle θ m values this can explain the difference of latitude to optimal tilt angles. The higher fraction of diffuse light in Seattle that also benefits the radiant exposure on the front side for small θ m might additionally increase this effect. Figure 5 shows how much the different irradiation components contribute to the annual radiant exposure for a bifacial module with d = 10 m module spacing, θ m = 34°tilt and albedo A = 30% in Dallas: about 64% of the total exposure arises from direct sunlight impinging onto the module front, 22% are due to diffuse skylight impinging onto the front but the fraction of light that reaches the front from the ground is almost negligible. However, of the 12% exposure received by the back, around 85% is reflected from the ground. Hence, the albedo only has little influence onto the energy yield of monofacial modules but is very relevant for bifacial modules. Figure 12 shows corresponding results for Seattle. Compared to Dallas, Seattle shows a nearly two per cent larger contribution by the back side. While the front side receives radiation with a ratio of nearly 3:1 of direct to diffuse light, for the back side, this ration is close to 1:1.
These results show that four factors drive the gain of bifacial modules instead of monofacial modules: the albedo of the ground, the module tilt angle, the module spacing and the overall fraction of diffuse light.
Also the mounting height h affects the bifacial gain. Increasing the mounting height monotonically raises the energy yield. Therefore it is difficult to optimise this parameter without knowing additional technical and commercial constraints. However, we find that the bifacial gain starts to saturate for a height above 0.5 m, which is in agreement with work from Kreinin et al. 17 Since a mounting height of h = 0.5 m seem realistic all simulations in our work are performed with this mounting height.

Minimising the electricity cost
In section 3 we discussed how to calculate the annual electrical energy yield and we analysed how the annual radiant exposure on the modules depends on the module spacing and tilt for two examples: Dallas and Seattle. In this section, we are going to derive a simple model for the electricity cost and perform some cost optimisations.

Levelied cost of electricity
As a measure for the electricity cost we use the levelized cost of electricity (LCOE), which is a key metric for electricity generation facilities. In the simplest case, the LCOE is given as the total cost C F spent in the facility during its lifetime T (in years) divided by the total amount of electric energy E total generated in that time, where E F is the electric energy generated by the PV field in one year. In more involved models also costs of capital and discount rates are taken into account. 3 The total cost can be split into two components, associated with the peak power C P (including modules, inverters, mounting etc.) and the land consumption C L (lease, fences, cables etc.) of the facility.
By considering a facility with a PV-field of M rows with N modules each the costs can be calculated per unit cell, The peak-power related costs per module C P,m are calculated with where c P denotes the peak-power associated costs given in [c p ] = $/kWp, which we use as input parameter. I P = 1 kW/m 2 is the peak irradiance as used for standardized PV characterization 28 and w and denote the module width and length, respectively. η f denotes the power conversion efficency on the front side of the solar cell.
The cost of land consumption per module depends on module width w and spacing d, with the land cost c L given in [c L ] = $/m 2 , which is an input parameter. 3   The annual generated electric energy of the PV field is given by with the annual yield EY according to eq. 4.
Combining eqs. (6)- (11) and simplifying leads to the expression which is independent of the field dimensions M and N and the module width w.
In this study, we assume for the overall costs of the PV system c m = 1000 $/kWp, which includes all costs over the lifetime of the solar park, such as PV module investment, balance of system cost, planning, capital cost and others. The land cost is not included in this quantity. The lifetime is assumed to be T = 25 years, a typical time span for the power warranty of solar cell modules. 24 In our optimization, we aim to minimize the LCOE as parameter of the module spacing d and the solar module tilt θ m . We perform the optimization for five land-cost scenarios c L , in which we assume to include all costs that are related to an increase of area such as lease, cables, fences etc. Table 3 gives an overview of the cost scenarios and the resulting fraction of the land costs on the total costs, (C L /C F ).

Optimisation method
As optimisation method we use Bayesian optimisation, which is well suited to find a global minimum of black box functions, which are expensive to evaluate. 29 Bayesian optimisation has been used in a wide variety of applications such as robotics, 30 hyper parameter tuning 31 or physical systems. 32,33 In principle, Bayesian optimisation consists of two components: a surrogate model that approximates the black box function and its uncertainty (based on previously evaluated data points) and an acquisition function that determines the next query point from the surrogate model. After evaluating the function for the queried data point the surrogate model is updated and the next step can be computed with the acquisition function. This cycle is repeated until a specified number of steps is reached or a convergences criteria is reached. We use the implementation from scikit-optimize with Gaussian process as surrogate model and expected improvement as acquisition function. 34

Optimization results
Traditionally, the optimal tilt and module spacing are often estimated with the winter solstice rule. 35,36 The optimal distance between two rows of modules is defined as the shortest distance for which the shadow of a row of modules does not hit the next row of modules in a specified solar time window (e.g. 9 am -3 pm) on winter solstice. As a rule of thumb the tilt is often chosen to be equivalent to the latitude of the facility location. However these rules do not consider the economic trade off between land costs and energy yield or typical weather patterns (e.g. foggy winters) that vary for different locations. Figures 6 and 7 shows the optimisation results for a field of (a) bifacial and (b) monofacial PV modules in Dallas and Seattle, respectively. Black dots mark evaluated data points, the red dot marks the found optimum and the color map shows the interpolation of the LCOE by the Gaussian process. The blue line indicates the winter solstice rule (9 am).
We see that the optimum shifts to smaller module spacing with increasing land cost. Further, also the optimal module tilt decreases in order to compensate for increased shadowing because of less module spacing. Overall, bifacial installations show large module spacing and higher tilt angles in optimal configurations compared to monofacial technology. With increasing land costs and therefore reduced optimal module spacing the cost landscape gets increasingly steep. The  5 10 module spacing (m) 5 10 module spacing (m) 5 10 module spacing (m)  sensitivity of the optimized parameters increases and using non-optimal geometrical configurations results in increasing yield loss. Seattle shows the same trends for optimal configuration in different cost scenarios. Compared to Dallas optimal tilt and spacing are higher.
Our optimisation results differ significantly from the geometric parameters obtained from the winter solstice rule. For Dallas the winter solstice rule only provides comparable optimal parameters for c L = 5$/m 2 . In Seattle, the optimal distances are shorter and the optimal module tilts are larger than expected from the winter solstice rule for all cost scenarios. This can be understood when considering the large share of diffuse light during the Seattle winter, which mitigates shading losses significantly.   5 10 module spacing (m) 5 10 module spacing (m) 5 10 module spacing (m)    scenario we see a reduction of LCOE of up to 23%. The rule-of-thumb approach shows its weakness especially in Seattle at high cost scenarios, where the cost landscape is very steep (see Figures 6 and 7).
From these results it is clear that the winter solstice rule is not able to properly reflect different economic trade-offs or different illumination conditions over the course of the year. This is especially true when setting the tilt angle to the latitude of the location. For a minimal LCOE module tilt and spacing should be optimised independently from each other. Further, typical weather patterns and the economic situation location must be taken into account.

Discussion
The results of all optimisations are summarised in Fig. 8 and Table 5. We see that the optimal LCOE increases slightly with the land cost. Further, in Seattle the LCOE difference between mono-and bifacial modules is larger as in Dallas. This is caused by the larger module tilt and diffuse light share in Seattle, which increases the fraction of illumination at the module back. As discussed above, the optimal module tilt decreases with increased land consumption cost c L . In general, we see that for a utility scale solar cell plant both, the module tilt and the distance between rows, affect the annual energy yield. Increasing the distance increases the energy yield and the costs per module while tilt can be optimized cost-neutral. The optimal distance between rows is a compromise between increasing costs with higher land use for higher distances and lower energy yield due to shading for lower distances. This is true also for monofacial modules but due to the increased relevance of light reflected from the ground it is more relevant for bifacial modules.
The optimal configuration for bifacial solar cells depends on the radiation conditions and the albedo of the facility location. With increasing latitude (and therefore lower solar elevation angles), albedo and diffuse light contribution the bifacial gain will be increased and therefore make this type of PV technology more attractive for utility scale developers.
Cost optimisations for PV installation are quickly outdated because PV module prices have been decreasing for many years and land cost is very volatile. However the optimal installation geometry only depends on the ratio of land cost related to total costs and not absolute values. Hence, at a scenario of c L = 10 $/m 2 and c P = 1500 $/kWp yields the same optimization result as c L = 5 $/m 2 and c P = 750 $/kWp.

Conclusions and Outlook
We developed a detailed model to calculate the irradiation onto both sides of a PV module, which is located in a large PV field. With this model, we could estimate the annual energy yield for monofacial and bifacial PV modules as a function of the module spacing and the module tilt. A simple model to calculate the levelized cost of electricity combined with a Bayesian optimisation algorithm allowed us to minimise the LCOE as a function of module spacing and module tilt for different land consumption costs.
Our results show that the bifacial gain and optimal geometry depend on the specific location and cost scenario. The bifacial gain can be expected to increase for locations with higher latitude and higher diffuse light share.
The usually used rule of thumb, no shadowing at winter-solstice and module tilt angle equal to the geographical latitude, leads to suboptimal module spacing and tilt combinations, because it does not account for economic trade-offs and the influence of the local climate. In contrast, optimising the parameters in Seattle can lead to a 23% reduction of LCOE for high land cost scenarios. This shows the significance of site-specific optimisation and helps users to identify the configurations yielding minimal LCOE.

Conflicts of interest
There are no conflicts to declare.

Supporting Information
A Climate diagrams  C Further results Figure 12 shows how much the different irradiation components contribute to the total energy yield for a bifacial module with d = 10 m module spacing and θ m = 42°tilt in Seattle: About 69.5% of the total energy arises from direct sunlight impinging onto the module front, 27.0% are due to diffuse skylight impinging onto the front but the fraction of light 16