# Probabilistic Modeling of Partial Shading of Photovoltaic Arrays

Authors: Yaw-Juen Wang; Ruey-Long Sheu
Journal: International Journal of Photoenergy
 (2015)

---

## 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 knowna 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.


---

## Body

## 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 Figure1 where I
L
m is the light-generating current, V
m and I
m are the module voltage and current, R
p
m and R
s
m are the parallel and series resistances, and V
d
m 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.Figure 1
Equivalent circuit of PV module.


#### 2.1.1. Maximum Power from a PV Module

The module temperature can be estimated by an empirical formula [9](1)
T
m
=
T
a
+
NOCT
-
20
0.8
·
G
,where T
a is the ambient temperature, NOCT refers to the nominal operating cell temperature, and G is the solar irradiance in kW/m2. The current-voltage (I-V) relation of the module can be written as an implicit function as(2)
f
m
I
m
,
V
m
,
G
=
I
m
-
I
L
G
-
I
d
m
G
,
V
d
m
-
V
d
m
R
p
m
=
0
,where I
d
m and V
d
m are the current and voltage of the equivalent diode. The diode voltage is given by(3)
V
d
m
=
V
m
+
I
m
R
s
m
.The diode current is given by Shockley’s equation:(4)
I
d
m
G
,
V
d
m
=
I
0
m
T
m
·
exp
⁡
q
V
d
m
n
m
k
T
m
G
-
1
,where I
0
m is the reverse saturation current of the equivalent diode, k the Boltzmann constant (=
1.381
×
1
0
-
23 J/K), q the electron charge (=
1.602
×
1
0
-
19 C), n
m the ideality factor of the module, and T
m the module temperature in Kelvin which depends on the irradiance G and can empirically be found by (1). It is noted in (4) that I
0
m is a function of T
m, and T
m is a function of G, so both I
0
m and T
m can be expressed as functions of G in (4). Therefore, the implicit function f
m becomes a function of three variables I
m, V
m, and G and is equal to zero.For a given irradiance levelG and a known module voltage V
m, the module current I
m 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)
P
V
m
,
G
=
V
m
·
I
m
V
m
,
G
.Obviously, the power P is a function of V
m and G and can be traced in P-V
m plane. A curve of P for a known irradiance level as a function of V
m is called the P-V 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 P-V relation to zero:(6)
d
d
V
m
P
V
m
,
G
=
0
.Practically, the maximum power point of a P-V 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 V
m
∗. The corresponding maximum power λ is(7)
λ
=
P
V
m
∗
,
G
=
h
G
.We note that V
m
∗ is indeed a function of G, so λ is simply a function of G. It is denoted by h
(
G
) 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. Table1 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 P
(
V
m
,
G
) in P-V 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 P
(
V
m
,
G
) at a different value of G. The maximum power λ is a monotonically increasing function of G, denoted by h
(
G
), 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 V
m
∗ 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.Table 1
Specification of SM50 PV module at standard test conditions (STC).
Electrical characteristic
Value (unit)
Cell number
36
Rated power
50.2 (W)
Current at the MPP
2.9 (A)
Voltage at the MPP
17.3 (V)
Short-circuit current
3.3 (A)
Open-circuit voltage
21.9 (V)Table 2
Parameters of the solar cells used in SM50 PV module at STC.
Parameter
Value (unit)
Reverse saturation current at 25°C (I
0)
2 × 10−7 (A)
Ideality factor (n)
1.4
Series resistance (R
s)
0.031 (Ω)
Parallel resistance (R
p)
15 (Ω)
Light-generated current atG
0 (I
L
0)
3.3 (A)
Reference solar irradiance (G
0)
1.0 (kW/m2)Figure 2
P-V curves of the module. The locus of maximum power is traced by a dashed curve that passes through the peak of each P-V curve.Figure 3
MPO of the moduleλ against the irradiance level G, obtained by a numerical method (scattered circles) and by an approximate function (curve).


### 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 thatG is a uniformly distributed random variable (RV) between G
L and G
H. The assumption of a uniform RV has some practical and theoretical basis. Since G 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 G have a one-to-one correspondence. Hence the probability that G is between g and g
+
d
g must be equal to the probability that λ is between λ and λ
+
d
λ. Mathematically, this relation can be written as [10](8)
f
g
g
·
d
g
=
f
λ
λ
·
d
λ
,where f
g
(
g
) and f
λ
(
λ
) are the pdfs of G and λ, respectively. Since G is uniform between G
L and G
H, its pdf is written as (9)
f
g
g
=
1
G
H
-
G
L
,
G
L
≤
g
≤
G
H
0
,
elsewhere.It follows from (8) that the pdf of λ can be given by(10)
f
λ
λ
=
f
g
g
·
d
g
d
λ
=
f
g
g
·
d
d
λ
h
-
1
λ
,where h
-
1
(
λ
) refers to the inverse function of h, namely, G expressed in terms of λ. The function h
(
G
) 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 η
(
G
) that closely approximates h
(
G
):(11)
λ
=
h
G
≅
η
G
=
a
1
+
G
/
c
b
.The curve shown in Figure 3 is obtained by (11) with the parameters a
=
113.7749, b
=
-
1.1734, and c
=
1.5690, and the scattered dots are obtained from (7). It can be seen that η
(
G
) agrees very well with h
(
G
). From (11), the inverse function h
-
1
(
λ
) can be approximately written as (12)
h
-
1
λ
=
G
λ
≅
c
a
λ
-
1
1
/
b
.Substituting (12) into (10) allows the pdf f
λ
(
λ
) to be derived as(13)
f
λ
λ
=
1
G
H
-
G
L
·
a
c
a
/
λ
-
1
1
/
b
-
1
-
b
λ
2
,
λ
L
≤
λ
≤
λ
H
0
,
elsewhere,where λ
L
=
h
(
G
L
) and λ
H
=
h
(
G
H
) 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 irradianceg
i. For each g
i, the corresponding maximum power λ
i is calculated numerically by solving (2), (6), and (7). That is, a series of random numbers g
1
,
g
2
,
g
3
,
…
,
g
N is generated, and a series of random samples λ
1
,
λ
2
,
λ
3
,
…
,
λ
N is recorded. The data set {
λ
1
,
λ
2
,
λ
3
,
…
,
λ
N
} 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 N 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 betweenG
L
=
0.1 kW/m2 and G
H
=
0.9 kW/m2 is studied as an example. Figure 4 shows the pdf of the maximum power λ obtained from (13). The cumulative distribution function (CDF) F
λ
(
λ
) 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.Figure 4
pdf and CDF of the maximum power obtained by the analytical method (curves) and pdf obtained by the Monte Carlo method of 104 simulation trials (scattered dots).Figure 5
pdf and CDF of the maximum power obtained by the analytical method (curves) and pdf obtained by the Monte Carlo method of 107 simulation trials (scattered dots).Figure 6
Amplified graph for the comparison of the pdfs obtained by the analytical method (curve) and the Monte Carlo method (scattered dots) of 107 simulation trials.


## 2.1. Analysis of Equivalent Circuit

The equivalent circuit of PV module is shown in Figure1 where I
L
m is the light-generating current, V
m and I
m are the module voltage and current, R
p
m and R
s
m are the parallel and series resistances, and V
d
m 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.Figure 1
Equivalent circuit of PV module.


### 2.1.1. Maximum Power from a PV Module

The module temperature can be estimated by an empirical formula [9](1)
T
m
=
T
a
+
NOCT
-
20
0.8
·
G
,where T
a is the ambient temperature, NOCT refers to the nominal operating cell temperature, and G is the solar irradiance in kW/m2. The current-voltage (I-V) relation of the module can be written as an implicit function as(2)
f
m
I
m
,
V
m
,
G
=
I
m
-
I
L
G
-
I
d
m
G
,
V
d
m
-
V
d
m
R
p
m
=
0
,where I
d
m and V
d
m are the current and voltage of the equivalent diode. The diode voltage is given by(3)
V
d
m
=
V
m
+
I
m
R
s
m
.The diode current is given by Shockley’s equation:(4)
I
d
m
G
,
V
d
m
=
I
0
m
T
m
·
exp
⁡
q
V
d
m
n
m
k
T
m
G
-
1
,where I
0
m is the reverse saturation current of the equivalent diode, k the Boltzmann constant (=
1.381
×
1
0
-
23 J/K), q the electron charge (=
1.602
×
1
0
-
19 C), n
m the ideality factor of the module, and T
m the module temperature in Kelvin which depends on the irradiance G and can empirically be found by (1). It is noted in (4) that I
0
m is a function of T
m, and T
m is a function of G, so both I
0
m and T
m can be expressed as functions of G in (4). Therefore, the implicit function f
m becomes a function of three variables I
m, V
m, and G and is equal to zero.For a given irradiance levelG and a known module voltage V
m, the module current I
m 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)
P
V
m
,
G
=
V
m
·
I
m
V
m
,
G
.Obviously, the power P is a function of V
m and G and can be traced in P-V
m plane. A curve of P for a known irradiance level as a function of V
m is called the P-V 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 P-V relation to zero:(6)
d
d
V
m
P
V
m
,
G
=
0
.Practically, the maximum power point of a P-V 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 V
m
∗. The corresponding maximum power λ is(7)
λ
=
P
V
m
∗
,
G
=
h
G
.We note that V
m
∗ is indeed a function of G, so λ is simply a function of G. It is denoted by h
(
G
) 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. Table1 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 P
(
V
m
,
G
) in P-V 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 P
(
V
m
,
G
) at a different value of G. The maximum power λ is a monotonically increasing function of G, denoted by h
(
G
), 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 V
m
∗ 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.Table 1
Specification of SM50 PV module at standard test conditions (STC).
Electrical characteristic
Value (unit)
Cell number
36
Rated power
50.2 (W)
Current at the MPP
2.9 (A)
Voltage at the MPP
17.3 (V)
Short-circuit current
3.3 (A)
Open-circuit voltage
21.9 (V)Table 2
Parameters of the solar cells used in SM50 PV module at STC.
Parameter
Value (unit)
Reverse saturation current at 25°C (I
0)
2 × 10−7 (A)
Ideality factor (n)
1.4
Series resistance (R
s)
0.031 (Ω)
Parallel resistance (R
p)
15 (Ω)
Light-generated current atG
0 (I
L
0)
3.3 (A)
Reference solar irradiance (G
0)
1.0 (kW/m2)Figure 2
P-V curves of the module. The locus of maximum power is traced by a dashed curve that passes through the peak of each P-V curve.Figure 3
MPO of the moduleλ against the irradiance level G, obtained by a numerical method (scattered circles) and by an approximate function (curve).


## 2.1.1. Maximum Power from a PV Module

The module temperature can be estimated by an empirical formula [9](1)
T
m
=
T
a
+
NOCT
-
20
0.8
·
G
,where T
a is the ambient temperature, NOCT refers to the nominal operating cell temperature, and G is the solar irradiance in kW/m2. The current-voltage (I-V) relation of the module can be written as an implicit function as(2)
f
m
I
m
,
V
m
,
G
=
I
m
-
I
L
G
-
I
d
m
G
,
V
d
m
-
V
d
m
R
p
m
=
0
,where I
d
m and V
d
m are the current and voltage of the equivalent diode. The diode voltage is given by(3)
V
d
m
=
V
m
+
I
m
R
s
m
.The diode current is given by Shockley’s equation:(4)
I
d
m
G
,
V
d
m
=
I
0
m
T
m
·
exp
⁡
q
V
d
m
n
m
k
T
m
G
-
1
,where I
0
m is the reverse saturation current of the equivalent diode, k the Boltzmann constant (=
1.381
×
1
0
-
23 J/K), q the electron charge (=
1.602
×
1
0
-
19 C), n
m the ideality factor of the module, and T
m the module temperature in Kelvin which depends on the irradiance G and can empirically be found by (1). It is noted in (4) that I
0
m is a function of T
m, and T
m is a function of G, so both I
0
m and T
m can be expressed as functions of G in (4). Therefore, the implicit function f
m becomes a function of three variables I
m, V
m, and G and is equal to zero.For a given irradiance levelG and a known module voltage V
m, the module current I
m 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)
P
V
m
,
G
=
V
m
·
I
m
V
m
,
G
.Obviously, the power P is a function of V
m and G and can be traced in P-V
m plane. A curve of P for a known irradiance level as a function of V
m is called the P-V 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 P-V relation to zero:(6)
d
d
V
m
P
V
m
,
G
=
0
.Practically, the maximum power point of a P-V 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 V
m
∗. The corresponding maximum power λ is(7)
λ
=
P
V
m
∗
,
G
=
h
G
.We note that V
m
∗ is indeed a function of G, so λ is simply a function of G. It is denoted by h
(
G
) 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. Table1 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 P
(
V
m
,
G
) in P-V 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 P
(
V
m
,
G
) at a different value of G. The maximum power λ is a monotonically increasing function of G, denoted by h
(
G
), 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 V
m
∗ 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.Table 1
Specification of SM50 PV module at standard test conditions (STC).
Electrical characteristic
Value (unit)
Cell number
36
Rated power
50.2 (W)
Current at the MPP
2.9 (A)
Voltage at the MPP
17.3 (V)
Short-circuit current
3.3 (A)
Open-circuit voltage
21.9 (V)Table 2
Parameters of the solar cells used in SM50 PV module at STC.
Parameter
Value (unit)
Reverse saturation current at 25°C (I
0)
2 × 10−7 (A)
Ideality factor (n)
1.4
Series resistance (R
s)
0.031 (Ω)
Parallel resistance (R
p)
15 (Ω)
Light-generated current atG
0 (I
L
0)
3.3 (A)
Reference solar irradiance (G
0)
1.0 (kW/m2)Figure 2
P-V curves of the module. The locus of maximum power is traced by a dashed curve that passes through the peak of each P-V curve.Figure 3
MPO of the moduleλ against the irradiance level G, obtained by a numerical method (scattered circles) and by an approximate function (curve).


## 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 thatG is a uniformly distributed random variable (RV) between G
L and G
H. The assumption of a uniform RV has some practical and theoretical basis. Since G 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 G have a one-to-one correspondence. Hence the probability that G is between g and g
+
d
g must be equal to the probability that λ is between λ and λ
+
d
λ. Mathematically, this relation can be written as [10](8)
f
g
g
·
d
g
=
f
λ
λ
·
d
λ
,where f
g
(
g
) and f
λ
(
λ
) are the pdfs of G and λ, respectively. Since G is uniform between G
L and G
H, its pdf is written as (9)
f
g
g
=
1
G
H
-
G
L
,
G
L
≤
g
≤
G
H
0
,
elsewhere.It follows from (8) that the pdf of λ can be given by(10)
f
λ
λ
=
f
g
g
·
d
g
d
λ
=
f
g
g
·
d
d
λ
h
-
1
λ
,where h
-
1
(
λ
) refers to the inverse function of h, namely, G expressed in terms of λ. The function h
(
G
) 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 η
(
G
) that closely approximates h
(
G
):(11)
λ
=
h
G
≅
η
G
=
a
1
+
G
/
c
b
.The curve shown in Figure 3 is obtained by (11) with the parameters a
=
113.7749, b
=
-
1.1734, and c
=
1.5690, and the scattered dots are obtained from (7). It can be seen that η
(
G
) agrees very well with h
(
G
). From (11), the inverse function h
-
1
(
λ
) can be approximately written as (12)
h
-
1
λ
=
G
λ
≅
c
a
λ
-
1
1
/
b
.Substituting (12) into (10) allows the pdf f
λ
(
λ
) to be derived as(13)
f
λ
λ
=
1
G
H
-
G
L
·
a
c
a
/
λ
-
1
1
/
b
-
1
-
b
λ
2
,
λ
L
≤
λ
≤
λ
H
0
,
elsewhere,where λ
L
=
h
(
G
L
) and λ
H
=
h
(
G
H
) 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 irradianceg
i. For each g
i, the corresponding maximum power λ
i is calculated numerically by solving (2), (6), and (7). That is, a series of random numbers g
1
,
g
2
,
g
3
,
…
,
g
N is generated, and a series of random samples λ
1
,
λ
2
,
λ
3
,
…
,
λ
N is recorded. The data set {
λ
1
,
λ
2
,
λ
3
,
…
,
λ
N
} 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 N 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.2.1. Analytical Modeling

We assume thatG is a uniformly distributed random variable (RV) between G
L and G
H. The assumption of a uniform RV has some practical and theoretical basis. Since G 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 G have a one-to-one correspondence. Hence the probability that G is between g and g
+
d
g must be equal to the probability that λ is between λ and λ
+
d
λ. Mathematically, this relation can be written as [10](8)
f
g
g
·
d
g
=
f
λ
λ
·
d
λ
,where f
g
(
g
) and f
λ
(
λ
) are the pdfs of G and λ, respectively. Since G is uniform between G
L and G
H, its pdf is written as (9)
f
g
g
=
1
G
H
-
G
L
,
G
L
≤
g
≤
G
H
0
,
elsewhere.It follows from (8) that the pdf of λ can be given by(10)
f
λ
λ
=
f
g
g
·
d
g
d
λ
=
f
g
g
·
d
d
λ
h
-
1
λ
,where h
-
1
(
λ
) refers to the inverse function of h, namely, G expressed in terms of λ. The function h
(
G
) 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 η
(
G
) that closely approximates h
(
G
):(11)
λ
=
h
G
≅
η
G
=
a
1
+
G
/
c
b
.The curve shown in Figure 3 is obtained by (11) with the parameters a
=
113.7749, b
=
-
1.1734, and c
=
1.5690, and the scattered dots are obtained from (7). It can be seen that η
(
G
) agrees very well with h
(
G
). From (11), the inverse function h
-
1
(
λ
) can be approximately written as (12)
h
-
1
λ
=
G
λ
≅
c
a
λ
-
1
1
/
b
.Substituting (12) into (10) allows the pdf f
λ
(
λ
) to be derived as(13)
f
λ
λ
=
1
G
H
-
G
L
·
a
c
a
/
λ
-
1
1
/
b
-
1
-
b
λ
2
,
λ
L
≤
λ
≤
λ
H
0
,
elsewhere,where λ
L
=
h
(
G
L
) and λ
H
=
h
(
G
H
) 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 irradianceg
i. For each g
i, the corresponding maximum power λ
i is calculated numerically by solving (2), (6), and (7). That is, a series of random numbers g
1
,
g
2
,
g
3
,
…
,
g
N is generated, and a series of random samples λ
1
,
λ
2
,
λ
3
,
…
,
λ
N is recorded. The data set {
λ
1
,
λ
2
,
λ
3
,
…
,
λ
N
} 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 N 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 betweenG
L
=
0.1 kW/m2 and G
H
=
0.9 kW/m2 is studied as an example. Figure 4 shows the pdf of the maximum power λ obtained from (13). The cumulative distribution function (CDF) F
λ
(
λ
) 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.Figure 4
pdf and CDF of the maximum power obtained by the analytical method (curves) and pdf obtained by the Monte Carlo method of 104 simulation trials (scattered dots).Figure 5
pdf and CDF of the maximum power obtained by the analytical method (curves) and pdf obtained by the Monte Carlo method of 107 simulation trials (scattered dots).Figure 6
Amplified graph for the comparison of the pdfs obtained by the analytical method (curve) and the Monte Carlo method (scattered dots) of 107 simulation trials.


## 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 I
a and V
a refer to the array current and voltage, respectively.Figure 7
Equivalent circuit of an array composed of two series PV modules.


### 3.1. Analysis of Equivalent Circuit

The array currentI
a and module voltages V
m
1 and V
m
2 shown in Figure 7 satisfy the following three equations:(14)
f
m
I
a
,
V
m
1
,
α
=
0
f
m
I
a
,
V
m
2
,
β
=
0
V
m
1
+
V
m
2
=
V
a
,where the implicit function f
m 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 P
a as a function of the array voltage V
a 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 P
a is a function of α, β, and V
a. When α and β are constants, P
a
(
V
a
) is a function with only one peak. The value of V
a that gives a maximum P
a at given values of α and β is denoted by V
a
∗ and is given by solving the equation (15)
∂
∂
V
a
P
a
V
a
,
α
,
β
=
0
.Practically, the value of V
a
∗ of a PV system is continuously regulated by an MPPT. It is noted that V
a
∗ is a function of α and β. The array maximum output power ζ can be expressed as(16)
ζ
=
h
a
α
,
β
=
P
a
V
a
∗
,
α
,
β
,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.Figure 8
P-V curves of array when α is fixed at 0.8 kW/m2. The locus of maximum power passes through the peak of each P-V curve.Figure 9
Three-dimensional view of the array maximum powerζ as a function of α and β.Figure 10
Isopower contours 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:(17)
f
α
α
=
1
G
H
-
G
L
,
G
L
≤
α
≤
G
H
0
,
elsewhere
f
β
β
=
1
G
H
-
G
L
,
G
L
≤
β
≤
G
H
0
,
elsewhere.The joint pdf of α and β can then be written as(18)
f
α
β
α
,
β
=
f
α
α
f
β
β
=
1
G
H
-
G
L
2
,
G
L
≤
α
≤
G
H
,
G
L
≤
β
≤
G
H
0
,
e
l
s
e
w
h
e
r
e
.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 D enclosed by two curves h
a
(
α
,
β
)
=
ζ and h
a
(
α
,
β
)
=
ζ
+
Δ
ζ as sketched by the hatched area in Figure 11. We can write(19)
P
ζ
≤
ζ
≤
ζ
+
Δ
ζ
=
P
ζ
≤
h
a
α
,
β
≤
ζ
+
Δ
ζ
=
P
α
,
β
∈
D
=
f
ζ
ζ
Δ
ζ
=
∬
D
f
α
β
α
,
β
·
d
α
d
β
.It follows that the pdf of ζ can be determined by(20)
f
ζ
ζ
=
1
Δ
ζ
∬
D
f
α
β
α
,
β
·
d
α
d
β
.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 h
a
(
α
,
β
)
=
ζ and h
a
(
α
,
β
)
=
ζ
+
Δ
ζ must be computed first. The function h
a
(
α
,
β
)
=
ζ 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 D 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.Figure 11
DomainD between two isopower curves over which a double integration is to be carried out.


#### 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 F
ζ
(
ζ
) calculated using the analytical and the MC methods.Figure 12
pdf and CDF ofζ obtained from the analytical method (curves) and the MC simulations (scattered dots).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 nth-percentile power. For example, the 80th-percentile power of the array is 54 W. Comparison of the nth-percentile powers of two arrays allows their performance to be evaluated in a probabilistic and quantitative sense.


## 3.1. Analysis of Equivalent Circuit

The array currentI
a and module voltages V
m
1 and V
m
2 shown in Figure 7 satisfy the following three equations:(14)
f
m
I
a
,
V
m
1
,
α
=
0
f
m
I
a
,
V
m
2
,
β
=
0
V
m
1
+
V
m
2
=
V
a
,where the implicit function f
m 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 P
a as a function of the array voltage V
a 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 P
a is a function of α, β, and V
a. When α and β are constants, P
a
(
V
a
) is a function with only one peak. The value of V
a that gives a maximum P
a at given values of α and β is denoted by V
a
∗ and is given by solving the equation (15)
∂
∂
V
a
P
a
V
a
,
α
,
β
=
0
.Practically, the value of V
a
∗ of a PV system is continuously regulated by an MPPT. It is noted that V
a
∗ is a function of α and β. The array maximum output power ζ can be expressed as(16)
ζ
=
h
a
α
,
β
=
P
a
V
a
∗
,
α
,
β
,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.Figure 8
P-V curves of array when α is fixed at 0.8 kW/m2. The locus of maximum power passes through the peak of each P-V curve.Figure 9
Three-dimensional view of the array maximum powerζ as a function of α and β.Figure 10
Isopower contours 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:(17)
f
α
α
=
1
G
H
-
G
L
,
G
L
≤
α
≤
G
H
0
,
elsewhere
f
β
β
=
1
G
H
-
G
L
,
G
L
≤
β
≤
G
H
0
,
elsewhere.The joint pdf of α and β can then be written as(18)
f
α
β
α
,
β
=
f
α
α
f
β
β
=
1
G
H
-
G
L
2
,
G
L
≤
α
≤
G
H
,
G
L
≤
β
≤
G
H
0
,
e
l
s
e
w
h
e
r
e
.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 D enclosed by two curves h
a
(
α
,
β
)
=
ζ and h
a
(
α
,
β
)
=
ζ
+
Δ
ζ as sketched by the hatched area in Figure 11. We can write(19)
P
ζ
≤
ζ
≤
ζ
+
Δ
ζ
=
P
ζ
≤
h
a
α
,
β
≤
ζ
+
Δ
ζ
=
P
α
,
β
∈
D
=
f
ζ
ζ
Δ
ζ
=
∬
D
f
α
β
α
,
β
·
d
α
d
β
.It follows that the pdf of ζ can be determined by(20)
f
ζ
ζ
=
1
Δ
ζ
∬
D
f
α
β
α
,
β
·
d
α
d
β
.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 h
a
(
α
,
β
)
=
ζ and h
a
(
α
,
β
)
=
ζ
+
Δ
ζ must be computed first. The function h
a
(
α
,
β
)
=
ζ 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 D 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.Figure 11
DomainD between two isopower curves over which a double integration is to be carried out.


### 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 F
ζ
(
ζ
) calculated using the analytical and the MC methods.Figure 12
pdf and CDF ofζ obtained from the analytical method (curves) and the MC simulations (scattered dots).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 nth-percentile power. For example, the 80th-percentile power of the array is 54 W. Comparison of the nth-percentile powers of two arrays allows their performance to be evaluated in a probabilistic and quantitative sense.


## 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:(17)
f
α
α
=
1
G
H
-
G
L
,
G
L
≤
α
≤
G
H
0
,
elsewhere
f
β
β
=
1
G
H
-
G
L
,
G
L
≤
β
≤
G
H
0
,
elsewhere.The joint pdf of α and β can then be written as(18)
f
α
β
α
,
β
=
f
α
α
f
β
β
=
1
G
H
-
G
L
2
,
G
L
≤
α
≤
G
H
,
G
L
≤
β
≤
G
H
0
,
e
l
s
e
w
h
e
r
e
.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 D enclosed by two curves h
a
(
α
,
β
)
=
ζ and h
a
(
α
,
β
)
=
ζ
+
Δ
ζ as sketched by the hatched area in Figure 11. We can write(19)
P
ζ
≤
ζ
≤
ζ
+
Δ
ζ
=
P
ζ
≤
h
a
α
,
β
≤
ζ
+
Δ
ζ
=
P
α
,
β
∈
D
=
f
ζ
ζ
Δ
ζ
=
∬
D
f
α
β
α
,
β
·
d
α
d
β
.It follows that the pdf of ζ can be determined by(20)
f
ζ
ζ
=
1
Δ
ζ
∬
D
f
α
β
α
,
β
·
d
α
d
β
.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 h
a
(
α
,
β
)
=
ζ and h
a
(
α
,
β
)
=
ζ
+
Δ
ζ must be computed first. The function h
a
(
α
,
β
)
=
ζ 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 D 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.Figure 11
DomainD between two isopower curves over which a double integration is to be carried out.


## 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 F
ζ
(
ζ
) calculated using the analytical and the MC methods.Figure 12
pdf and CDF ofζ obtained from the analytical method (curves) and the MC simulations (scattered dots).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 nth-percentile power. For example, the 80th-percentile power of the array is 54 W. Comparison of the nth-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 Figure7 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 R
b, and a voltage source E
b standing for the diode barrier voltage as shown in Figure 13 in which I
m
1 and I
m
2 refer to the module currents.Figure 13
Equivalent circuit of a PV array composed of two series modules, each protected by a bypass diode.


### 4.1. Analysis of Equivalent Circuit

The circuit in Figure13 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:(21)
f
m
I
m
1
,
V
m
1
,
α
=
0
(22)
f
m
I
m
2
,
V
m
2
,
β
=
0
(23)
V
m
1
+
V
m
2
=
V
a
(24)
I
m
1
-
V
m
1
+
E
b
R
b
-
V
m
1
=
I
m
2
-
V
m
2
+
E
b
R
b
-
V
m
2
(25)
I
a
=
I
m
1
-
V
m
1
+
E
b
R
b
-
V
m
1
,where the implicit function f
m 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 R
b of the bypass diodes has been modeled as a voltage-controlled resistance allowing for rectification. The value of R
b is controlled by the voltage across the bypass diode V
b by(26)
R
b
V
b
=
0.2
Ω
,
V
b
>
0
1
0
8
Ω
,
V
b
≤
0which means that R
b has a low value when forward biased and a high value when reversely biased. The voltages across the two bypass diodes are -
V
m
1 and -
V
m
2, respectively. The forward voltage drop E
b of the bypass diodes is assumed to be 0.7 V. The P-V 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 P-V 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 P-V curve does not apply to the current case because it may mistakenly find a peak that is not the global MPO.Figure 14
P-V curves of the array for α = 0.8 kW/m2 and different values of β. The P-V curves have two peaks because of the effects of bypass diodes.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 byV
a
† that can be determined by the aforementioned algorithm. It is noted that V
a
† is a function of α and β. The maximum power ζ delivered by the array is also a function of α and β and can be given by(27)
ζ
=
h
b
α
,
β
.The function h
b 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 h
b
(
α
,
β
). 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.Figure 15
Three-dimensional view of the array maximum powerζ as a function of α and β.Figure 16
Isopower contours inα-β plane.Figure 17
Region Φ between two isopower curves over which a double integral is to be evaluated.


### 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 h
b
(
α
,
β
)
=
ζ and h
b
(
α
,
β
)
=
ζ
+
Δ
ζ as sketched by the hatched area in Figure 17. We can write(28)
P
ζ
≤
ζ
≤
ζ
+
Δ
ζ
=
P
ζ
≤
h
b
α
,
β
≤
ζ
+
Δ
ζ
=
P
α
,
β
∈
Φ
=
f
ζ
ζ
Δ
ζ
=
∬
Φ
f
α
β
α
,
β
·
d
α
d
β
.The pdf f
ζ
(
ζ
) can be determined by(29)
f
ζ
ζ
=
1
Δ
ζ
∬
Φ
f
α
β
α
,
β
·
d
α
d
β
.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.


#### 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 G
H and G
L. 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 f
ζ
(
ζ
). Figure 18 depicts the function f
ζ
(
ζ
) 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 F
ζ
(
ζ
) found by both methods are also sketched in 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).The applications of the CDF curve shown in Figure18 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.


## 4.1. Analysis of Equivalent Circuit

The circuit in Figure13 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:(21)
f
m
I
m
1
,
V
m
1
,
α
=
0
(22)
f
m
I
m
2
,
V
m
2
,
β
=
0
(23)
V
m
1
+
V
m
2
=
V
a
(24)
I
m
1
-
V
m
1
+
E
b
R
b
-
V
m
1
=
I
m
2
-
V
m
2
+
E
b
R
b
-
V
m
2
(25)
I
a
=
I
m
1
-
V
m
1
+
E
b
R
b
-
V
m
1
,where the implicit function f
m 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 R
b of the bypass diodes has been modeled as a voltage-controlled resistance allowing for rectification. The value of R
b is controlled by the voltage across the bypass diode V
b by(26)
R
b
V
b
=
0.2
Ω
,
V
b
>
0
1
0
8
Ω
,
V
b
≤
0which means that R
b has a low value when forward biased and a high value when reversely biased. The voltages across the two bypass diodes are -
V
m
1 and -
V
m
2, respectively. The forward voltage drop E
b of the bypass diodes is assumed to be 0.7 V. The P-V 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 P-V 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 P-V curve does not apply to the current case because it may mistakenly find a peak that is not the global MPO.Figure 14
P-V curves of the array for α = 0.8 kW/m2 and different values of β. The P-V curves have two peaks because of the effects of bypass diodes.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 byV
a
† that can be determined by the aforementioned algorithm. It is noted that V
a
† is a function of α and β. The maximum power ζ delivered by the array is also a function of α and β and can be given by(27)
ζ
=
h
b
α
,
β
.The function h
b 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 h
b
(
α
,
β
). 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.Figure 15
Three-dimensional view of the array maximum powerζ as a function of α and β.Figure 16
Isopower contours inα-β plane.Figure 17
Region Φ between two isopower curves over which a double integral is to be evaluated.


## 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 h
b
(
α
,
β
)
=
ζ and h
b
(
α
,
β
)
=
ζ
+
Δ
ζ as sketched by the hatched area in Figure 17. We can write(28)
P
ζ
≤
ζ
≤
ζ
+
Δ
ζ
=
P
ζ
≤
h
b
α
,
β
≤
ζ
+
Δ
ζ
=
P
α
,
β
∈
Φ
=
f
ζ
ζ
Δ
ζ
=
∬
Φ
f
α
β
α
,
β
·
d
α
d
β
.The pdf f
ζ
(
ζ
) can be determined by(29)
f
ζ
ζ
=
1
Δ
ζ
∬
Φ
f
α
β
α
,
β
·
d
α
d
β
.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.


### 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 G
H and G
L. 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 f
ζ
(
ζ
). Figure 18 depicts the function f
ζ
(
ζ
) 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 F
ζ
(
ζ
) found by both methods are also sketched in 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).The applications of the CDF curve shown in Figure18 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.


## 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 h
b
(
α
,
β
)
=
ζ and h
b
(
α
,
β
)
=
ζ
+
Δ
ζ as sketched by the hatched area in Figure 17. We can write(28)
P
ζ
≤
ζ
≤
ζ
+
Δ
ζ
=
P
ζ
≤
h
b
α
,
β
≤
ζ
+
Δ
ζ
=
P
α
,
β
∈
Φ
=
f
ζ
ζ
Δ
ζ
=
∬
Φ
f
α
β
α
,
β
·
d
α
d
β
.The pdf f
ζ
(
ζ
) can be determined by(29)
f
ζ
ζ
=
1
Δ
ζ
∬
Φ
f
α
β
α
,
β
·
d
α
d
β
.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.


## 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 G
H and G
L. 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 f
ζ
(
ζ
). Figure 18 depicts the function f
ζ
(
ζ
) 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 F
ζ
(
ζ
) found by both methods are also sketched in 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).The applications of the CDF curve shown in Figure18 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.



---

*Source: 863637-2015-12-07.xml*