Abstract
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 PV module 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 PV modules, 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.
1. 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–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.
2. Modeling of a Single PV Module
2.1. Analysis of Equivalent Circuit
The equivalent circuit of PV module is shown in Figure 1 where 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.
2.1.1. Maximum Power from a PV Module
The module temperature can be estimated by an empirical formula [9]where is the ambient temperature, NOCT refers to the nominal operating cell temperature, and is the solar irradiance in kW/m2. The current-voltage (-) relation of the module can be written as an implicit function aswhere and are the current and voltage of the equivalent diode. The diode voltage is given byThe diode current is given by Shockley’s equation:where is the reverse saturation current of the equivalent diode, the Boltzmann constant ( J/K), the electron charge ( 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 is a function of , and is a function of , so both 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 byObviously, 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 isWe note that is indeed a function of , so is simply a function of . It is denoted by in (7).
2.1.2. Locus of the 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. Table 1 lists the module’s electrical specifications at the standard test conditions (STC). The values of the parameters used in this study are summarized in Table 2. Figure 2 depicts the 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.
2.2. 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.
2.2.1. 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 bywhere 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 , , and , and the scattered dots are obtained from (7). It can be seen that agrees very well with . From (11), the inverse function can be approximately written as Substituting (12) into (10) allows the pdf to be derived aswhere and are the lower and upper bounds of .
2.2.2. 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 is generated, and a series of random samples is recorded. The data set 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.
2.3. Numerical Examples
A simple case of a SM-50 PV module receiving a uniformly distributed irradiance between kW/m2 and kW/m2 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 104 and 107 trials are also shown by the scattered dots in Figures 4 and 5. The pdf found by the 104-trial MC is close enough to the analytical model with perceptible fluctuations around the analytical curve. The 107-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.
3. 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.
3.1. Analysis of Equivalent Circuit
The array current and module voltages and 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 Table 2 are again used for numerical applications. When α is fixed at 0.8 kW/m2, 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 Figure 8 in which the locus 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 aswhere 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.
3.2. Probabilistic Analysis
3.2.1. 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 asThe 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 writeIt follows that the pdf of can be determined byThe 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.
3.2.2. Monte Carlo Simulation
MC simulation of the two-module 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 106 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 th-percentile 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.
4. 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 and refer to the module currents.
4.1. Analysis of Equivalent Circuit
The circuit in Figure 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: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 bywhich means that has a low value when forward biased and a high value when reversely biased. The voltages across the two bypass diodes are and , 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/m2 and β varies from 0.1 to 0.9 kW/m2 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 right-hand 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 byThe 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 looks like a paper glider with lateral wings extending outward toward both sides of the ridge line. It merits noting that the projections of the ridge lines in both Figures 9 and 15 are indeed the 45-degree lines in Figures 10, 11, 16, and 17. The isopower contours are depicted in Figure 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.
4.2. Probabilistic Analysis
4.2.1. Analytical Model
The marginal pdfs of α and β remain unchanged as shown by (17). The joint pdf also remains unchanged as (18). The probability that the array maximum power is included between and must be equal to the probability that the coordinate appears within the region Φ that is enclosed by the curves and as sketched by the hatched area in Figure 17. We can writeThe pdf can be determined byEquation (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.
4.2.2. 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 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 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.77 W, 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.
5. 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 104 and 107 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.
Principal Symbols and Abbreviations
| : | Array voltage |
| : | Array current |
| : | Module voltage |
| : | Module current |
| : | Solar irradiance |
| : | Implicit function of , , and |
| : | Module power |
| : | Voltage at maximum power |
| : | pdf of the solar irradiance |
| : | pdf of the maximum power of a module |
| : | pdf of the maximum power of an array |
| : | 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 |
| : | Inverse function of |
| , : | Solar irradiance received by modules |
| : | Maximum power output of a module |
| : | Maximum power output of an array |
| CDF: | Cumulative distribution function |
| epf: | Empirical probability function |
| MC: | Monte Carlo |
| MPO: | Maximum power output |
| pdf: | Probability density function |
| RV: | Random variable. |
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
This work was supported by research grants from the Ministry of Science and Technology of Taiwan under the Contract no. MOST 103-2221-E-224-028.