Probabilistic Modeling of Partial Shading of Photovoltaic Arrays Yaw-JuenWang and

The probabilistic modeling technique is used in this paper to ascertain the performance of energy conversion of a photovoltaic (PV) array under the influence of partial shading. Three cases of increasing complexity are studied in this paper. The first case is to find the probability density function (pdf) of the maximum power output (MPO) of a single PV module exposed to a random level of solar irradiance. Given that the probability distribution of solar irradiance is known a priori, the pdf of the MPO of the PVmodule is derived analytically. A Monte Carlo simulation is then conducted to validate the analysis. This is followed in the second case by studying the MPO of an array composed of two series PVmodules, each exposed to a random and independent level of irradiance. The third case further involves the effects of bypass diodes which are commonly installed to reduce partial shading losses.


Introduction
Partial shading is a commonly encountered problem in applications of photovoltaic (PV) energy generation.Partial shading is avoided in the design stage by considering all the possible shadows caused by nearby objects.However, unexpected shadows covering a part of PV modules often happen.Examples are shadows caused by bird dung, unmelted snow, fallen tree leaves, newly growing tree branches, and so forth.Therefore, partial shading is not easily avoided even for very professional PV projects.
Methods for analysis of partial shading are well documented in literature.Alonso-García et al. [1] involved solving equations of PV networks using the Newton-Raphson method.On the other hand, Wang and Hsu [2] tackled this problem by developing a piecewise linear circuit model that allowed almost all circuit simulation software to model a partially shaded PV system.
As partial shading exhibits some degree of uncertainty, efforts have been made to extend traditional analysis into cases where the pattern, number, and shading percentage of shaded PV modules may vary at random.The works of Gautam and Kaushika [3][4][5][6] considered several shading patterns in PV arrays of different connection configurations.A randomly generated shading pattern was tested in Wang and Hsu's study [7].Ramaprabha and Mathur [8] applied 15 random patterns to arrays of 10 different sizes.Although the studies [7,8] have introduced the idea of randomization, further works still need to be undertaken, particularly for deeper theoretical development of a fundamental probabilistic model of partial shading.
This paper is intended to ascertain how the maximum power output (MPO) of a PV array is distributed when its modules are subject to randomly varying levels of solar irradiance.Three cases of increasing complexity are studied to explain the procedure of analytically deriving the probability density function (pdf) of the MPO.The first case, which is a univariate analysis, is to find the pdf of the MPO of a PV module exposed to randomly varying irradiance.Once the univariate problem can be solved, it is extended to a more complicated bivariate case where a PV array composed of two series modules, each exposed to independent and randomly varying irradiance, is analyzed.In the third case, the modules are assumed to be equipped with bypass diodes, which are commonly used to reduce partial shading losses.

Modeling of a Single PV Module
where   is the ambient temperature, NOCT refers to the nominal operating cell temperature, and  is the solar irradiance in kW/m 2 .The current-voltage (-) relation of the module can be written as an implicit function as where   and   are the current and voltage of the equivalent diode.The diode voltage is given by The diode current is given by Shockley's equation: where  0 is the reverse saturation current of the equivalent diode,  the Boltzmann constant (= 1.381 × 10 −23 J/K),  the electron charge (= 1.602 × 10 −19 C),   the ideality factor of the module, and   the module temperature in Kelvin which depends on the irradiance  and can empirically be found by (1).It is noted in (4) that  0 is a function of   , and   is a function of , so both  0 and   can be expressed as functions of  in (4).Therefore, the implicit function   becomes a function of three variables   ,   , and  and is equal to zero.
For a given irradiance level  and a known module voltage   , the module current   can be obtained by solving (2) using a numerical method such as the Newton-Raphson method.The electric power outputted by the module is given by  (  , ) =   ⋅   (  , ) . (5) Obviously, the power  is a function of   and  and can be traced in -  plane.A curve of  for a known irradiance level as a function of   is called the - relation of a module.The voltage at which a PV module can output a maximum power at a given level of irradiance can be found by equating the derivative of the - relation to zero: Practically, the maximum power point of a - curve is tracked by a maximum power point tracker (MPPT).Several algorithms are available for realizing the function of MPPT such as the perturbation and observation method and the incremental conductance method.Let the solution to (6) be  *  .The corresponding maximum power  is We note that  *  is indeed a function of , so  is simply a function of .It is denoted by ℎ() in (7).function (  , ) in - plane (solid curves) for different levels of irradiance.The locus of maximum power in dashed curve is also traced in the plane.It is noted that the locus passes through the peak of each solid curve standing for (  , ) at a different value of .The maximum power  is a monotonically increasing function of , denoted by ℎ(), which is shown in Figure 3 and can be compared with the dashed curve depicted in Figure 2. In Figure 3, the scattered circles are calculated by solving (6) for  *  and subsequently substituting it into (7).The curve in Figure 3 is obtained by a curve-fitting algorithm which is to be discussed in the next section.

Probabilistic Modeling.
Shading bears some degree of uncertainty.When a module is shadowed, the solar irradiance it receives reduces to a level that depends on the percentage of shading.Shading is modeled by a randomly varying intensity of solar irradiance in this study.

Analytical Modeling.
We assume that  is a uniformly distributed random variable (RV) between   and   .
The assumption of a uniform RV has some practical and theoretical basis.Since  must be bounded, a uniform RV satisfies this feature.The probability density function (pdf) of a uniform RV is the easiest one to manipulate, largely simplifying the procedure of mathematical derivation.The work is then to find the pdf of . Figure 3 shows that  and  have a one-to-one correspondence.Hence the probability that  is between  and + must be equal to the probability that  is between  and  + .Mathematically, this relation can be written as [10] where   () and   () are the pdfs of  and , respectively.Since  is uniform between   and   , its pdf is written as It follows from ( 8) that the pdf of  can be given by where ℎ −1 () refers to the inverse function of ℎ, namely,  expressed in terms of .The function ℎ() in ( 7) is indeed very complicated and is impossible to be written in an analytical closed-form, let alone its inverse function.In this paper, we make a detour to avoid this tricky problem by finding an analytical function () that closely approximates ℎ(): The curve shown in Figure 3 is obtained by (11) with the parameters  = 113.7749, = −1.1734,and  = 1.5690, and the scattered dots are obtained from (7).It can be seen that () agrees very well with ℎ().From (11), the inverse function ℎ −1 () can be approximately written as Substituting ( 12) into (10) allows the pdf   () to be derived as where   = ℎ(  ) and   = ℎ(  ) are the lower and upper bounds of .
International Journal of Photoenergy

Monte Carlo Simulation.
A more straightforward way to obtain the pdf of the module MPO is to resort to the Monte Carlo (MC) simulation method.The MC method involves generating a uniformly distributed solar irradiance   .For each   , the corresponding maximum power   is calculated numerically by solving (2), (6), and (7).That is, a series of random numbers  1 ,  2 ,  3 , . . .,   is generated, and a series of random samples  1 ,  2 ,  3 , . . .,   is recorded.The data set { 1 ,  2 ,  3 , . . .,   } is employed to make an empirical probability function (epf) that is the relative frequency histogram corresponding to the data set.The acquired epf thus approximates the theoretical pdf when the number of sampling data  is large enough.For a rough profile of the pdf, several thousands of samples are enough.However, for a stable and accurate result, several millions of simulation runs may be needed.

Numerical Examples.
A simple case of a SM-50 PV module receiving a uniformly distributed irradiance between   = 0.1 kW/m 2 and   = 0.9 kW/m 2 is studied as an example.Figure 4 shows the pdf of the maximum power  obtained from (13).The cumulative distribution function (CDF)   () is also calculated and depicted in Figure 4 to ensure that the total area under the pdf equals unit.Further verification has been carried out by performing the MC simulation.The results of the MC simulations of 10 4 and 10 7 trials are also shown by the scattered dots in Figures 4 and 5.The pdf found by the 10 4 -trial MC is close enough to the analytical model with perceptible fluctuations around the analytical curve.The 10 7 -trial MC shows nearly perfect agreement with the analytical curve.To more precisely compare the pdfs obtained by the analytical and the MC methods, an amplified graph has been traced in Figure 6 in which the curve first goes down a little bit and then rises up steeply.The scattered dots standing for the histogram found by the MC simulation closely stick the analytical curve, showing very nice agreement.

Array Composed of Two Series Modules
Partial shading can be modeled by two modules connected in series that receive different levels of irradiance  and  as shown by the equivalent circuit of Figure 7 where   and   refer to the array current and voltage, respectively.

Analysis of Equivalent Circuit.
The array current   and module voltages  1 and  2 shown in Figure 7 satisfy the following three equations: where the implicit function   has been defined in (2).The parameters of the SM50 module listed in    of the maximum power is also shown by a dashed curve.We notice that the array power   is a function of , , and   .When  and  are constants,   (  ) is a function with only one peak.The value of   that gives a maximum   at given values of  and  is denoted by  *  and is given by solving the equation Practically, the value of  *  of a PV system is continuously regulated by an MPPT.It is noted that  *  is a function of  and .The array maximum output power  can be expressed as where  is apparently a function of  and  and is depicted as a three-dimensional graph as shown in Figure 9. Figure 10 shows the corresponding contours in - plane for  varying from 10 to 80 W. It is noted that the contours are symmetric to the 45-degree line, that is, the line - = 0 in - plane.This symmetry property is important when making a double integral over an area in - plane.

Analytical Modeling.
In the probabilistic analysis, the irradiances ( and ) striking on the two PV modules are considered as two uniformly distributed RVs having their respective pdfs: The joint pdf of  and  can then be written as The array MPO  is a function of two RVs  and , so  is also a RV.The probability that  is included between  and  + Δ must be equal to the probability that  and  are in the region  enclosed by two curves ℎ  (, ) =  and ℎ  (, ) =  + Δ as sketched by the hatched area in Figure 11.We can write It follows that the pdf of  can be determined by The pdf expressed by (20) has been evaluated numerically.
To calculate the double integral in (20) for a given array power , the lower and upper boundaries ℎ  (, ) =  and ℎ  (, ) =  + Δ must be computed first.The function ℎ  (, ) =  cannot be expressed analytically.Only numerical data are available.When evaluating the double integral, linear interpolation is done to find the correct value of  for a given  so that the area  can be sliced into many small strips for numerical integration.The integration is made easier and faster by making full use of the symmetry.Since the isopower curves are symmetric to the 45-degree line, only the integration to right of the 45-degree line needs to be evaluated, which enhances the computational speed by a factor of two and also enhances the accuracy since the integration to the left of the 45-degree line is more difficult because of steeper slopes of the isopower curves.

Monte Carlo Simulation.
MC simulation of the twomodule PV system has been performed to verify the calculation result of the analytical model.In the MC simulation trials, uniformly distributed random numbers have been generated according to (17) to simulate the irradiances  and  and then calculate the corresponding array maximum power .A number of 10 6 simulation trials have been carried out, which allows the epf of  to be obtained.The pdfs of  calculated using the analytical model (curve) and MC simulation (circle dots) are depicted in Figure 12.It can be seen that good agreement has been achieved.Also shown in the figure are the CDFs   () calculated using the analytical and the MC methods.It is interesting to see how the CDF curve can be used to evaluate the array performance under the assumed random shading conditions.The CDF gives the probability that  is equal to or less than a certain value.For example, the CDF curve in Figure 12 shows that the array has a probability of 0.8 to supply a maximum power equal to or less than 54 W. A more useful application of the CDF is to find the thpercentile power.For example, the 80th-percentile power of the array is 54 W. Comparison of the th-percentile powers of two arrays allows their performance to be evaluated in a probabilistic and quantitative sense.

Array Composed of Two Series Modules with Bypass Diodes
A more general case that takes into account the effects of bypass diodes is studied in this paragraph.The circuit shown in Figure 7 is extended by connecting a bypass diode to each module.The bypass diodes are modeled by a series branch containing an ideal diode, a resistance   , and a voltage source   standing for the diode barrier voltage as shown in Figure 13 in which  1 and  2 refer to the module currents.13 is a bit more complicated than that shown in Figure 7 because of insertion of bypass diodes.The following equations can be written as per KVL and KCL:

Analysis of Equivalent Circuit. The circuit in Figure
where the implicit function   is defined in (2).Equations ( 21) and ( 22) describe the nonlinearity of the PV modules that receive irradiance levels  and .Equation (24) implies that the module current minus the bypass diode current equals the array current.In ( 24) and (25), the resistance   of the bypass diodes has been modeled as a voltage-controlled resistance allowing for rectification.The value of   is controlled by the voltage across the bypass diode   by which means that   has a low value when forward biased and a high value when reversely biased.The voltages across the two bypass diodes are − 1 and − 2 , respectively.The forward voltage drop   of the bypass diodes is assumed to be 0.7 V.The - characteristics of the array when  stays equal to 0.8 kW/m 2 and  varies from 0.1 to 0.9 kW/m 2 are shown in Figure 14.These curves are obtained by solving (21)-( 25).The locus of maximum power has also been traced.The effects of bypass diodes cause the - curves to possess two peaks.The locus of maximum power that goes down through the peaks suddenly shifts to the left because the peaks on the right-hand side are no longer the global maxima.The method used in ( 6) and ( 15) that takes the derivative of a - curve does not apply to the current case because it may mistakenly find a peak that is not the global MPO.
To find the global maximum power, an algorithm that searches for the curve peak from both the left-and righthand sides has been developed.Comparison of the left-and right-hand side peaks allows the global maximum to be correctly determined.The voltage at which the array outputs a maximum power is denoted by  †  that can be determined by the aforementioned algorithm.It is noted that  †  is a function of  and .The maximum power  delivered by the array is also a function of  and  and can be given by  = ℎ  (, ) . ( The function ℎ  reflects the relation that connects , , and . For given values of  and , the array maximum power  can be found by solving (21)-( 27). Figure 15 shows the 3D graph of function ℎ  (, ).The curved surface depicted in Figure 9 slopes downward from the ridge line and directly reaches the bottom while the 3D surface in Figure 15 16.Unlike the contours in Figure 11 that look like a family of hyperbolas, the curves in Figure 16 bend nearly vertically toward the horizontal and vertical axis when an irradiance level is about twice the other irradiance level.When a contour curve begins to bend, it implies that a bypass diode starts to conduct.The pdf   () can be determined by Equation ( 29) looks very similar to (20), but the double integral over the region Φ is indeed much trickier because the contour curves turn nearly perpendicular to their original directions when a bypass diode starts to conduct, which requires that the integration step size be reduced after the turning point to achieve satisfactory accuracy.The evaluation of (29) also takes advantage of the symmetry of the contours about the 45-degree line to save the calculation time.

Monte Carlo Simulation.
The MC simulation method has been used to find the pdf and CDF of the array maximum power  when the two randomly varying irradiances  and  are uniformly distributed between   and   .One hundred thousand random pairs of (, ) have been generated, and the corresponding array maximum power  has been calculated.The epf has been found to approximate the pdf   ().Figure 18 depicts the function   () obtained analytically by (29) (curve) and by the MC simulation (scattered dots).It can be seen that there is a sudden drop, that is, a discontinuous point, near  = 38 W. Good agreement between the results of MC simulation and the analytical model is observed.The CDF   () found by both methods are also sketched in Figure 18.The applications of the CDF curve shown in Figure 18 are similar to those explained in Section 3.2 for Figure 12.For example, the CDF in Figure 18 reveals that the array supplies a power  not exceeding 30 W with a probability of 0.25.Equivalently, this implies that the array has a probability of 0.75 to provide more than 30 W of power.The pdf curve has also several applications, one of which is to find the expected value of , which has the physical meaning of long-run power output or average power output.We can compare International Journal of Photoenergy the expected values obtained from pdf curves in Figure 12 (without bypass diode) and in Figure 18 (with bypass diodes) and show the advantages of installing bypass diodes.The pdf in Figure 12 gives the expected value of  of 36.77W, against the expected value obtained from the pdf in Figure 18 of 39.58 W. We can conclude that for the same condition of random shading, on average or in the long run, the array equipped with bypass diodes supplies 2.81 W more power than the array without bypass diodes.

Conclusion
Partial shading is modeled in this study by a randomly varying level of solar irradiance.A method that employs an analytical function to approximate the relationship between the irradiance and the MPO and hence allows an analytical closed-form of the pdf of the MPO to be derived has been presented.MC simulation has been conducted to verify the result of the proposed method.Good agreement has been achieved between the proposed method and MC simulation.The results of the MC simulation of 10 4 and 10 7 trials have also been compared.As expected, the simulation with a higher number of simulation trials has obtained a better agreement with the analytical model.
The study has been extended to a PV array composed of two series modules, each exposed to a random level of irradiance.The MPO of the array becomes a function of two random variables.A bivariate analysis has then been conducted to find the pdf of the array MPO using an analytical method that involves double integral over a region enclosed by two adjacent isopower contours.MC simulation runs of the two-module array have also been carried out to verify the result of the bivariate analysis and good agreement has been obtained.
The bivariate study on the two-module array has gone further to the case where bypass diodes are considered.An analytical method has again successfully been developed to find the pdf of the array MPO.Results of MC simulation and the analytical model have been compared and the analysis has been validated.Function of the maximum power of a module ℎ  (, ): Function of the maximum power of an array without bypass diode ℎ  (, ): Function of the maximum power of an array with bypass diodes ℎ −1 (): Inverse function of ℎ() , :

Principal Symbols and Abbreviations
Solar irradiance received by modules : M a x i m u mpo w e ro u t p u to fam od ul e : M a x i m u mp o w e ro u t p u to fa na r r a y CDF: Cumulative distribution function epf: Empirical probability function MC: Monte Carlo MPO: Maximum power output pdf: Probability density function RV: Random variable.

Figure 4 :Figure 5 :
Figure 4: pdf and CDF of the maximum power obtained by the analytical method (curves) and pdf obtained by the Monte Carlo method of 10 4 simulation trials (scattered dots).

Figure 6 :
Figure 6: Amplified graph for the comparison of the pdfs obtained by the analytical method (curve) and the Monte Carlo method (scattered dots) of 10 7 simulation trials.

Figure 7 :Figure 8 :
Figure 7: Equivalent circuit of an array composed of two series PV modules.

Figure 9 :𝛼 (kW/m 2 )Figure 10 :
Figure 9: Three-dimensional view of the array maximum power  as a function of  and .

Figure 11 :Figure 12 :
Figure11: Domain  between two isopower curves over which a double integration is to be carried out.

Figure 13 :
Figure 13: Equivalent circuit of a PV array composed of two series modules, each protected by a bypass diode.

Figure 14 :
Figure 14: - curves of the array for  = 0.8 kW/m 2 and different values of .The - curves have two peaks because of the effects of bypass diodes.

Figure 15 :Figure 16 :
Figure 15: Three-dimensional view of the array maximum power  as a function of  and .

4. 2 .Figure 17 :
Figure17: Region Φ between two isopower curves over which a double integral is to be evaluated.

Figure 18 :
Figure 18: pdf and CDF of the array maximum power obtained from the analytical model (curves) and from the MC simulations (scattered dots).
[2].Analysis of Equivalent Circuit.The equivalent circuit of PV module is shown in Figure1where   is the light-generating current,   and   are the module voltage and current,   and   are the parallel and series resistances, and   is the diode voltage.The analysis of the equivalent circuit has been studied in the work of Wang and Hsu[2]in detail.Here only a brief review is given.
[9].1.Maximum Power from a PV Module.The module temperature can be estimated by an empirical formula[9]

Table 1 :
Specification of SM50 PV module at standard test conditions (STC).

Table 2 :
Parameters of the solar cells used in SM50 PV module at STC.
Maximum Power.The electrical specifications of SM50 PV module are taken in this study to illustrate the maximum power characteristic.The SM50 module is composed of 36 identical PV cells connected in series.Table1lists the module's electrical specifications at the standard test conditions (STC).The values of the parameters used in this study are summarized in Table2.Figure2depicts the

Table 2
are again used for numerical applications.When  is fixed at 0.8 kW/m 2 , curves of the array output power   as a function of the array voltage   for different values of  are obtained by solving (14) and are depicted in Figure8in which the locus (): pdf of the solar irradiance   (): pdf of the maximum power of a module   (): pdfofthemaximumpowerofanarray ℎ():