# Energy Losses and Voltage Stability Study in Distribution Network with Distributed Generation

Authors: Hongwei Ren; Congying Han; Tiande Guo; Wei Pei
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

With the distributed generation technology widely applied, some system problems
such as overvoltages and undervoltages are gradually remarkable, which are
caused by distributed generations like wind energy system (WES) and photovoltaic
system (PVS) because of their probabilistic output power which relied on natural conditions.
Since the impacts of WES and PVS are important in the distribution system
voltage quality, we study these in this paper using new models with the probability
density function of node voltage and the cumulative distribution function of
total losses. We apply these models to solve the IEEE33 distribution system to be
chosen in IEEE standard database. We compare our method with the Monte Carlo
simulation method in three different cases, respectively. In the three cases, these
results not only can provide the important reference information for the next stage
optimization design, system reliability, and safety analysis but also can reduce amount
of calculation.


---

## Body

## 1. Introduction

Electric power systems have been originally designed based on the unidirectional power flow. Nevertheless, in the last years the conception of distributed generation (DG) such as wind power generation and solar power generation has led to new consideration on the distribution networks (DN) [1]. As the ratio of the distributed generation in power system expanded, to study the effect of distributed generation on the system steady run is more and more important. Because wind power and solar power are with stochastic volatility, the penetration of DG may impact the operation of DN in both beneficial and detrimental ways [2–9]. The positive impacts of DG may possibly be voltage support, power loss reduction, support of ancillary services, and improved reliability, whereas, negative ones included protection coordination, dynamic stability, and islanding. Numerous researchers have dealt with the issue of size and site of DG into DNs. A group of articles optimize sizing and/or siting of DG units in order to obtain maximum benefits, such as maximum loss reduction or reliability and minimum cost [10–18]. However, the above-mentioned papers use power flow analysis either for a certain loading condition or for a few specific scenarios (e.g., seasonal loadings) based on measured data or default test cases [3–5, 9, 10, 14–21]. And some have not considered stochastic volatility of distributed generation. This paper is mainly about energy losses and voltage stability assessment in distribution network with distributed generation considering stochastic volatility.With the development of science and technology and the raised awareness of environmental protection, every country is becoming more and more interested in the renewable energy sources specifically because they are reproducible and nonpolluting. These technologies include hydro-, wind [9], solar [10], biomass [11], and tidal technology. Among these renewable energies, wind and solar technology have evolved very rapidly over the past decade and the reduction of capital costs, the improvement of reliability, and the efficiency have made the wind and solar power be able to compete with conventional power generation [12]. The renewable DG technologies like wind and solar have special characteristics due to their main source of energy. Obviously, the primary energy source of a wind turbine is wind. The wind speed is not a constant quantity during the operation of wind turbine and is highly dependent on climate condition of the area where wind turbine is installed. The solar technology is also dependent on the climate and geographical location. It is the reason that they exhibit uncertainty and variability in their output [13]. Some methods are proposed to model the impact of these uncertainties on distribution network performance; the general and powerful tool is based on Monte Carlo simulation for simulating the uncertainties, but this method needs too large amount of calculation, especially the power flow calculation. Hence, this paper proposed a tool to reduce the amount of calculation about uncertainties for distribution network.This paper is organized as follows: in Section2, preliminary theory, the fuzzy theory, Monte Carlo simulation, and the half invariant method used for the stochastic analysis are presented. The model is given in Section 3. The proposed solution algorithm and simulation result are presented in Section 4. Some general conclusions are presented in Section 5.


## 2. Preliminary Theory

The value of load in each bus and DG generation are controllable with decisions of their owners. In this section, we prepare the preliminary theory used in this paper. One part is about fuzzy mathematics theory; the other is Monte Carlo and half invariants method.


### 2.1. Fuzzy Mathematics Theory

In fuzzy mathematics theory, a membership function is defined which describes how much each element belongs to a fuzzy set.μ
A
(
x
) is a membership function that takes values in the interval [
0,1
]. For each element x
∈
A, A is a fuzzy set of universe of discourse U. In this paper, fuzzy trapezoidal membership (FTM) with a notation A
=
(
a
min
⁡
,
a
L
,
a
U
,
a
max
⁡
) is used as shown in Figure 1.Figure 1
Fuzzy trapezoidal number.In engineering problems, the question is, knowing the uncertain input variablesx
i, how to give the membership function of y. The α
-cut method [16] answers this question in the following way: for a given fuzzy set A defined on universe of discourse, that is, U, the crisp set A
α is defined as all elements of U which have membership degree to A, greater than or equal to α, as calculated in the following equation:
(1)
A
α
=
{
x
∈
U
∣
μ
A
(
x
)
⩾
α
}
.
The α
-cut of each input variable, that is, x
i
α, is calculated using (1); then the α
-cut of y, that is, Y
α, is calculated as follows:
(2)
Y
α
=
(
y
α
,
y
α
)
,
(3)
y
α
=
min
⁡
f
(
x
)
,
(4)
y
α
=
max
⁡
f
(
x
)
,
(5)
x
∈
(
x
α
,
x
α
)
.
This means for each α-cut, (3) and (4) are solved. The upper bound of y
α is obtained by (4), and the lower bound of y
α is obtained by (3). The defuzzification is a mathematical process for converting a fuzzy number into a crisp one [16]. In [19], the centroid method is used for defuzzification of fuzzy numbers. The defuzzified value of a given fuzzy quantity, that is, A
*, is calculated as follows:
(6)
A
*
=
∫
‍
μ
A
(
x
)
x
d
x
∫
‍
μ
A
(
x
)
d
x
.
Transforming fuzzy variables into random variables is commonly used in engineering. This paper adopts a conversion method defined as follows:
(7)
f
(
x
)
=
μ
A
(
x
)
∫
‍
μ
A
(
x
)
d
x
.
This method not only retained the distribution information of fuzzy variables membership functions but also met the completeness and the nonnegativity of the probability density function.


### 2.2. Monte Carlo Simulation and Half Invariants (Cumulants) Method

The main concept of Monte Carlo simulation (MCS) method is described as follows: suppose a multivariable function, namely,y, y
=
f
(
Z
), where Z
=
(
Z
1
,
…
,
Z
m
), in which Z
1 to Z
m are random variables with their own probability distribution function (PDF). In [21], the MCS acts as follows: first of all, it will generate a value, that is, Z
i
Θ, for each input variable Z
i using its own PDF and form Z
Θ
=
(
Z
1
Θ
,
…
,
Z
m
Θ
) and then calculate the value of y
Θ using y
Θ
=
f
(
Z
Θ
). This process will be repeated for a number of iterations. The trend of the output, that is, y, will determine its PDF.Some of uncertain input parameters follow from PDF, such as the value of wind which follows a Weibull PDF [20]. MCS is a powerful tool for analyzing the uncertainties which follow any PDF. But MCS calculation is too big, hence this paper presents a new method—half invariants (cumulants) method instead of Monte Carlo simulation method.Just as expectation and variance, half invariant is also a numerical characteristic of random variableX, which can be calculated by the moment of the random variable X.IfY is a linear function of random variables X
1 and X
2 with their own PDF, the problem is, knowing the PDFs of all variables X
1 and X
2, what would be the PDF of Y? The half invariants (cumulants) method applied half invariant properties and Gram-Charlier series theorem: first, calculate the half invariant of random variable Y using the half invariant of random variables X
1 and X
2: then get the PDF of random variable Y using Gram-Charlier series theorem. This can avoid complicated convolution operation and a large number of Monte Carlo simulations.Gram-Charlier Series Theorem. Suppose X is a random variable; then the PDF of X can be expressed as follows:
(8)
PDF
(
x
)
=
φ
(
x
)
+
c
1
φ
(
1
)
(
x
)
+
c
2
φ
(
2
)
(
x
)
+
⋯
,
where φ
(
x
) is the PDF of the standard normal distribution, φ
(
i
)
(
x
) is the ith derivative, and c
i is the coefficient which can be calculated by the moment of the random variable X.


## 2.1. Fuzzy Mathematics Theory

In fuzzy mathematics theory, a membership function is defined which describes how much each element belongs to a fuzzy set.μ
A
(
x
) is a membership function that takes values in the interval [
0,1
]. For each element x
∈
A, A is a fuzzy set of universe of discourse U. In this paper, fuzzy trapezoidal membership (FTM) with a notation A
=
(
a
min
⁡
,
a
L
,
a
U
,
a
max
⁡
) is used as shown in Figure 1.Figure 1
Fuzzy trapezoidal number.In engineering problems, the question is, knowing the uncertain input variablesx
i, how to give the membership function of y. The α
-cut method [16] answers this question in the following way: for a given fuzzy set A defined on universe of discourse, that is, U, the crisp set A
α is defined as all elements of U which have membership degree to A, greater than or equal to α, as calculated in the following equation:
(1)
A
α
=
{
x
∈
U
∣
μ
A
(
x
)
⩾
α
}
.
The α
-cut of each input variable, that is, x
i
α, is calculated using (1); then the α
-cut of y, that is, Y
α, is calculated as follows:
(2)
Y
α
=
(
y
α
,
y
α
)
,
(3)
y
α
=
min
⁡
f
(
x
)
,
(4)
y
α
=
max
⁡
f
(
x
)
,
(5)
x
∈
(
x
α
,
x
α
)
.
This means for each α-cut, (3) and (4) are solved. The upper bound of y
α is obtained by (4), and the lower bound of y
α is obtained by (3). The defuzzification is a mathematical process for converting a fuzzy number into a crisp one [16]. In [19], the centroid method is used for defuzzification of fuzzy numbers. The defuzzified value of a given fuzzy quantity, that is, A
*, is calculated as follows:
(6)
A
*
=
∫
‍
μ
A
(
x
)
x
d
x
∫
‍
μ
A
(
x
)
d
x
.
Transforming fuzzy variables into random variables is commonly used in engineering. This paper adopts a conversion method defined as follows:
(7)
f
(
x
)
=
μ
A
(
x
)
∫
‍
μ
A
(
x
)
d
x
.
This method not only retained the distribution information of fuzzy variables membership functions but also met the completeness and the nonnegativity of the probability density function.


## 2.2. Monte Carlo Simulation and Half Invariants (Cumulants) Method

The main concept of Monte Carlo simulation (MCS) method is described as follows: suppose a multivariable function, namely,y, y
=
f
(
Z
), where Z
=
(
Z
1
,
…
,
Z
m
), in which Z
1 to Z
m are random variables with their own probability distribution function (PDF). In [21], the MCS acts as follows: first of all, it will generate a value, that is, Z
i
Θ, for each input variable Z
i using its own PDF and form Z
Θ
=
(
Z
1
Θ
,
…
,
Z
m
Θ
) and then calculate the value of y
Θ using y
Θ
=
f
(
Z
Θ
). This process will be repeated for a number of iterations. The trend of the output, that is, y, will determine its PDF.Some of uncertain input parameters follow from PDF, such as the value of wind which follows a Weibull PDF [20]. MCS is a powerful tool for analyzing the uncertainties which follow any PDF. But MCS calculation is too big, hence this paper presents a new method—half invariants (cumulants) method instead of Monte Carlo simulation method.Just as expectation and variance, half invariant is also a numerical characteristic of random variableX, which can be calculated by the moment of the random variable X.IfY is a linear function of random variables X
1 and X
2 with their own PDF, the problem is, knowing the PDFs of all variables X
1 and X
2, what would be the PDF of Y? The half invariants (cumulants) method applied half invariant properties and Gram-Charlier series theorem: first, calculate the half invariant of random variable Y using the half invariant of random variables X
1 and X
2: then get the PDF of random variable Y using Gram-Charlier series theorem. This can avoid complicated convolution operation and a large number of Monte Carlo simulations.Gram-Charlier Series Theorem. Suppose X is a random variable; then the PDF of X can be expressed as follows:
(8)
PDF
(
x
)
=
φ
(
x
)
+
c
1
φ
(
1
)
(
x
)
+
c
2
φ
(
2
)
(
x
)
+
⋯
,
where φ
(
x
) is the PDF of the standard normal distribution, φ
(
i
)
(
x
) is the ith derivative, and c
i is the coefficient which can be calculated by the moment of the random variable X.


## 3. Half Invariants Modeling

The assumptions for modeling the two types of uncertainties, constraints, and the objective functions are described as follows.


### 3.1. Uncertainty Modeling

#### 3.1.1. Load

It is assumed that the values of load in each bus are controllable with decisions of their owners. In this paper, we assume that the distribution network operators (DNO) can just describe them with a membership function as follows:(9)
S
h
,
i
D
=
S
i
,
f
D
×
DL
F
h
×
(
ξ
min
⁡
,
ξ
L
,
ξ
U
,
ξ
max
⁡
)
,
where S
i
,
f
D is the apparent forecasted value of peak load in bus i and DL
F
h is the demand level factor at demand level h which takes values between 0 and 1. Finally, S
h
,
i
D is the fuzzy value of demand in bus i and demand level h.DG generation pattern: the amount of energy which a controllable DG unit injects into the network is uncertain, and usually it depends on the decisions of DG owner so the DNO cannot have a PDF of it if there is not much historic data about it. The output power of a controllable DG unit is modeled using a membership function as follows:(10)
P
h
,
i
dg
=
C
i
,
f
dg
×
(
ζ
min
⁡
,
ζ
L
,
ζ
U
,
ζ
max
⁡
)
,
where C
i
,
f
dg is the capacity of DG unit installed in bus i and P
h
,
i
dg is the active power of a DG unit in bus i in demand level h.Photovoltaic generation pattern: the amount of solar radiation that reaches the ground, besides on the daily and yearly apparent motion of the sun, depends on the geographical location (latitude and altitude) and on the climatic conditions (e.g., cloud cover). The generation schedule of a photovoltaic generation pattern highly depends on the irradiance in the site. The variation of irradiance, that is,E, can be modeled using a beta PDF as follows:
(11)
PDF
(
E
)
=
Γ
(
p
+
q
)
Γ
(
p
)
Γ
(
q
)
(
E
E
max
⁡
)
p
-
1
(
1
-
E
E
max
⁡
)
q
-
1
,
where E and E
max
⁡ are the actual light intensity and maximum light intensity and p and q are the shape of the beta distribution parameters.The generated power of the photovoltaic generation is determined as follows:(12)
P
h
,
i
sloar
=
E
A
γ
,
where A is the area of the solar panels and γ is the photoelectric conversion efficiency.Wind turbine generation pattern: the generation schedule of a wind turbine highly depends on the wind speed in the site. The variation of wind speed, that is,v, can be modeled using a Weibull PDF [22] and its characteristic function which relates the wind speed and the output of a wind turbine [23] as follows:
(13)
PDF
(
v
)
=
(
β
α
)
(
v
α
)
β
-
1
exp
⁡
(
-
(
v
α
)
β
)
,
where β is the shape parameters and α is the scale parameters of the Weibull PDF of wind speed in the zone under study. The generated power of the wind turbine is determined using its characteristics as follows:
(14)
P
h
,
i
wind
=
{
p
i
,
r
wind
,
v
in
c
<
v
<
v
rate
,
v
-
v
in
c
v
out
c
-
v
in
c
p
i
,
r
wind
,
v
rate
<
v
<
v
out
c
,
0
,
else
,
where P
i
,
r
wind is the rated power of wind turbine installed in bus i, P
h
,
i
wind is the generated power of wind turbine in bus i and demand level h, v
out
c is the cut out speed, v
in
c is the cut in speed, and v
rate is the rated speed of the wind turbine.


### 3.2.  Active Losses

The power flow equations must be satisfied in each demand levelh and at each bus i as follows:
(15)
P
h
,
i
net
=
V
h
,
i
∑
‍
(
Y
i
j
V
h
,
i
(
G
i
j
cos
⁡
θ
i
j
+
B
i
j
sin
θ
i
j
)
)
,
Q
h
,
i
net
=
V
h
,
i
∑
‍
(
Y
i
j
V
h
,
i
(
G
i
j
cos
⁡
θ
i
j
-
B
i
j
sin
θ
i
j
)
)
,
where P
h
,
i
net and Q
h
,
i
net are the net active and reactive power injected to the network in bus i at level h. The above equation can be written in the matrix form as W
=
f
(
X
). The power flow equations at each branch are given as follows:
(16)
P
h
,
i
j
=
V
i
V
j
(
G
i
j
cos
⁡
θ
i
j
+
B
i
j
sin
θ
i
j
)
+
t
i
j
G
i
j
V
i
2
,
Q
h
,
i
j
=
-
V
i
V
j
(
G
i
j
cos
⁡
θ
i
j
-
B
i
j
sin
θ
i
j
)
+
(
B
i
j
-
b
i
j
0
)
G
i
j
V
i
2
.
The above equation can be written in the matrix form as Z
=
g
(
X
). Make the first order Taylor expansion as follows:
(17)
W
=
f
(
X
0
)
+
J
|
X
0
Δ
X
,
Z
=
g
(
X
0
)
+
G
|
X
0
Δ
X
.
Then we get the following linear relationship:
(18)
Δ
X
=
S
0
·
Δ
W
,
Δ
Z
=
G
0
·
Δ
X
=
G
0
·
S
0
·
Δ
W
=
T
0
Δ
W
.
The total active loss of the network in each demand level is equal to the sum of all active power injected to each bus as follows:
(19)
Los
s
h
=
∑
‍
(
P
h
,
i
j
+
P
h
,
j
i
)
·
τ
h
Los
s
h
=
H
(
X
)
.
By the Taylor expansion and linearization, the total active loss of the network is equal to the sum of all active power injected to each bus, that is,
(20)
Los
s
net
=
∑
‍
Los
s
h
,
where Los
s
h
=
H
(
X
0
)
+
Δ
Los
s
h
=
Los
s
h
,
0
+
Δ
Los
s
h.


## 3.1. Uncertainty Modeling

### 3.1.1. Load

It is assumed that the values of load in each bus are controllable with decisions of their owners. In this paper, we assume that the distribution network operators (DNO) can just describe them with a membership function as follows:(9)
S
h
,
i
D
=
S
i
,
f
D
×
DL
F
h
×
(
ξ
min
⁡
,
ξ
L
,
ξ
U
,
ξ
max
⁡
)
,
where S
i
,
f
D is the apparent forecasted value of peak load in bus i and DL
F
h is the demand level factor at demand level h which takes values between 0 and 1. Finally, S
h
,
i
D is the fuzzy value of demand in bus i and demand level h.DG generation pattern: the amount of energy which a controllable DG unit injects into the network is uncertain, and usually it depends on the decisions of DG owner so the DNO cannot have a PDF of it if there is not much historic data about it. The output power of a controllable DG unit is modeled using a membership function as follows:(10)
P
h
,
i
dg
=
C
i
,
f
dg
×
(
ζ
min
⁡
,
ζ
L
,
ζ
U
,
ζ
max
⁡
)
,
where C
i
,
f
dg is the capacity of DG unit installed in bus i and P
h
,
i
dg is the active power of a DG unit in bus i in demand level h.Photovoltaic generation pattern: the amount of solar radiation that reaches the ground, besides on the daily and yearly apparent motion of the sun, depends on the geographical location (latitude and altitude) and on the climatic conditions (e.g., cloud cover). The generation schedule of a photovoltaic generation pattern highly depends on the irradiance in the site. The variation of irradiance, that is,E, can be modeled using a beta PDF as follows:
(11)
PDF
(
E
)
=
Γ
(
p
+
q
)
Γ
(
p
)
Γ
(
q
)
(
E
E
max
⁡
)
p
-
1
(
1
-
E
E
max
⁡
)
q
-
1
,
where E and E
max
⁡ are the actual light intensity and maximum light intensity and p and q are the shape of the beta distribution parameters.The generated power of the photovoltaic generation is determined as follows:(12)
P
h
,
i
sloar
=
E
A
γ
,
where A is the area of the solar panels and γ is the photoelectric conversion efficiency.Wind turbine generation pattern: the generation schedule of a wind turbine highly depends on the wind speed in the site. The variation of wind speed, that is,v, can be modeled using a Weibull PDF [22] and its characteristic function which relates the wind speed and the output of a wind turbine [23] as follows:
(13)
PDF
(
v
)
=
(
β
α
)
(
v
α
)
β
-
1
exp
⁡
(
-
(
v
α
)
β
)
,
where β is the shape parameters and α is the scale parameters of the Weibull PDF of wind speed in the zone under study. The generated power of the wind turbine is determined using its characteristics as follows:
(14)
P
h
,
i
wind
=
{
p
i
,
r
wind
,
v
in
c
<
v
<
v
rate
,
v
-
v
in
c
v
out
c
-
v
in
c
p
i
,
r
wind
,
v
rate
<
v
<
v
out
c
,
0
,
else
,
where P
i
,
r
wind is the rated power of wind turbine installed in bus i, P
h
,
i
wind is the generated power of wind turbine in bus i and demand level h, v
out
c is the cut out speed, v
in
c is the cut in speed, and v
rate is the rated speed of the wind turbine.


## 3.1.1. Load

It is assumed that the values of load in each bus are controllable with decisions of their owners. In this paper, we assume that the distribution network operators (DNO) can just describe them with a membership function as follows:(9)
S
h
,
i
D
=
S
i
,
f
D
×
DL
F
h
×
(
ξ
min
⁡
,
ξ
L
,
ξ
U
,
ξ
max
⁡
)
,
where S
i
,
f
D is the apparent forecasted value of peak load in bus i and DL
F
h is the demand level factor at demand level h which takes values between 0 and 1. Finally, S
h
,
i
D is the fuzzy value of demand in bus i and demand level h.DG generation pattern: the amount of energy which a controllable DG unit injects into the network is uncertain, and usually it depends on the decisions of DG owner so the DNO cannot have a PDF of it if there is not much historic data about it. The output power of a controllable DG unit is modeled using a membership function as follows:(10)
P
h
,
i
dg
=
C
i
,
f
dg
×
(
ζ
min
⁡
,
ζ
L
,
ζ
U
,
ζ
max
⁡
)
,
where C
i
,
f
dg is the capacity of DG unit installed in bus i and P
h
,
i
dg is the active power of a DG unit in bus i in demand level h.Photovoltaic generation pattern: the amount of solar radiation that reaches the ground, besides on the daily and yearly apparent motion of the sun, depends on the geographical location (latitude and altitude) and on the climatic conditions (e.g., cloud cover). The generation schedule of a photovoltaic generation pattern highly depends on the irradiance in the site. The variation of irradiance, that is,E, can be modeled using a beta PDF as follows:
(11)
PDF
(
E
)
=
Γ
(
p
+
q
)
Γ
(
p
)
Γ
(
q
)
(
E
E
max
⁡
)
p
-
1
(
1
-
E
E
max
⁡
)
q
-
1
,
where E and E
max
⁡ are the actual light intensity and maximum light intensity and p and q are the shape of the beta distribution parameters.The generated power of the photovoltaic generation is determined as follows:(12)
P
h
,
i
sloar
=
E
A
γ
,
where A is the area of the solar panels and γ is the photoelectric conversion efficiency.Wind turbine generation pattern: the generation schedule of a wind turbine highly depends on the wind speed in the site. The variation of wind speed, that is,v, can be modeled using a Weibull PDF [22] and its characteristic function which relates the wind speed and the output of a wind turbine [23] as follows:
(13)
PDF
(
v
)
=
(
β
α
)
(
v
α
)
β
-
1
exp
⁡
(
-
(
v
α
)
β
)
,
where β is the shape parameters and α is the scale parameters of the Weibull PDF of wind speed in the zone under study. The generated power of the wind turbine is determined using its characteristics as follows:
(14)
P
h
,
i
wind
=
{
p
i
,
r
wind
,
v
in
c
<
v
<
v
rate
,
v
-
v
in
c
v
out
c
-
v
in
c
p
i
,
r
wind
,
v
rate
<
v
<
v
out
c
,
0
,
else
,
where P
i
,
r
wind is the rated power of wind turbine installed in bus i, P
h
,
i
wind is the generated power of wind turbine in bus i and demand level h, v
out
c is the cut out speed, v
in
c is the cut in speed, and v
rate is the rated speed of the wind turbine.


## 3.2.  Active Losses

The power flow equations must be satisfied in each demand levelh and at each bus i as follows:
(15)
P
h
,
i
net
=
V
h
,
i
∑
‍
(
Y
i
j
V
h
,
i
(
G
i
j
cos
⁡
θ
i
j
+
B
i
j
sin
θ
i
j
)
)
,
Q
h
,
i
net
=
V
h
,
i
∑
‍
(
Y
i
j
V
h
,
i
(
G
i
j
cos
⁡
θ
i
j
-
B
i
j
sin
θ
i
j
)
)
,
where P
h
,
i
net and Q
h
,
i
net are the net active and reactive power injected to the network in bus i at level h. The above equation can be written in the matrix form as W
=
f
(
X
). The power flow equations at each branch are given as follows:
(16)
P
h
,
i
j
=
V
i
V
j
(
G
i
j
cos
⁡
θ
i
j
+
B
i
j
sin
θ
i
j
)
+
t
i
j
G
i
j
V
i
2
,
Q
h
,
i
j
=
-
V
i
V
j
(
G
i
j
cos
⁡
θ
i
j
-
B
i
j
sin
θ
i
j
)
+
(
B
i
j
-
b
i
j
0
)
G
i
j
V
i
2
.
The above equation can be written in the matrix form as Z
=
g
(
X
). Make the first order Taylor expansion as follows:
(17)
W
=
f
(
X
0
)
+
J
|
X
0
Δ
X
,
Z
=
g
(
X
0
)
+
G
|
X
0
Δ
X
.
Then we get the following linear relationship:
(18)
Δ
X
=
S
0
·
Δ
W
,
Δ
Z
=
G
0
·
Δ
X
=
G
0
·
S
0
·
Δ
W
=
T
0
Δ
W
.
The total active loss of the network in each demand level is equal to the sum of all active power injected to each bus as follows:
(19)
Los
s
h
=
∑
‍
(
P
h
,
i
j
+
P
h
,
j
i
)
·
τ
h
Los
s
h
=
H
(
X
)
.
By the Taylor expansion and linearization, the total active loss of the network is equal to the sum of all active power injected to each bus, that is,
(20)
Los
s
net
=
∑
‍
Los
s
h
,
where Los
s
h
=
H
(
X
0
)
+
Δ
Los
s
h
=
Los
s
h
,
0
+
Δ
Los
s
h.


## 4. Algorithm and Simulation Results

This paper supposed that the injection power at each bus is independent and made power flow equation linearization. The uncertainties are neither random variables nor fuzzy variables. Transforming fuzzy variables into random variables with the method is proposed ahead.


### 4.1. Algorithm

Step 1.
Input feeder data,h
=
1.Step 2.
Read load and DG data ath level.Step 3.
Run power flow with Newton-Raphson ath level.Step 4.
Calculate the half invariance of the injection power flow of the generator and load ath level.Step 5.
Calculate the half invariance ofΔ
W.Step 6.
Calculate the half invariance ofΔ
X according to S
h
,
0.Step 7.
Calculate the PDF of variateΔ
X according to Gram-Charlier Series theorem.Step 8.
Calculate the half invariance ofΔ
Z according to T
h
,
0.Step 9.
Calculate the half invariance ofΔ
Loss according to M
h
,
0.Step 10.
Calculate the PDF ofΔ
Los
s
h according to Gram-Charlier series theorem at h level and store the data.Step 11.
According Step7, calculate the probability of out-of-limit voltage at h level and store the data.Step 12.
h
=
h
+
1; if h
⩽
24, then turn to Step 2; else turn to Step 13.Step 13.
According to the data stored in Step10, fit the PDF of Δ
Loss of one year.Step 14.
According to the data stored in Step11, estimate the time of out-of-limit voltage of one year.Step 15.
End.


### 4.2. Simulation Results

The proposed methodology is applied to a IEEE33 distribution network which is shown in Figure2 and joined with two distributed generations (i.e., node34 and node35).Figure 2
IEEE33 distribution network.The first one is the Monte Carlo simulation method and the second is the method proposed by us in this paper. The two methods, respectively, are used in three cases.Case
1 is not joined with any distributed generation (node34 and node35 are out of work); case
2 is joined with two wind turbine generations (node34 and node35), while case
3 is joined with one wind turbine generation (node34) and one photovoltaic generation (node35). The simulation results are as shown in Figures 4 and 5. Simulation parameters are given in Table 1. It is assumed that there are 24 demand levels in each year with equal duration of τ
h
=
365 h. The variations of demand level factors are depicted in Figure 3.Table 1
Data used in this paper.
Parameters
Unit
Value
c
8.78
V
min
⁡
pu
0.95
v
in
c
m/s
3
v
rate
m/s
13
ξ
L
0.925
ξ
max
⁡
1.15
ζ
L
0.9
ζ
max
⁡
1
k
1.75
V
max
⁡
pu
1.05
v
out
c
m/s
25
ξ
min
⁡
0.850
ξ
U
1.075
ζ
min
⁡
0
ζ
U
1
τ
h
h
365Figure 3
The variations ofDL
F
h in each demand level.Figure 4
Voltage probability density function of the 33rd node.Figure 5
Cumulative distribution function of total losses.Some cumulants are given in Tables1, 2, and 3. The first order cumulant is expected of node voltage; the second order cumulant is the variance. Figure 4 is the probability density function of the voltage amplitude of node33. From this figure, we can know that, in case
1, the random fluctuations of voltage are similar to normal distribution; probability of out-of-limit voltage is almost zero. When adding distributed generation, on one hand, line voltage had improved significantly and load node voltage rises. On the other hand, because of the randomness of the wind power and photovoltaic power, the node voltage fluctuation and the probability of out-of-limit voltage significantly increased. By comparing with case
2 and case
3, we can also find that, for wind and solar hybrid power systems, due to its complementarity, its impact on system voltage fluctuation is relatively smaller compared with single wind power and the probability of out-of-limit voltage is reduced obviously. The experimental results show that the proposed fuzzy variables are effective to describe load instead of normal variable (Table 4).Table 2
Moment of bus voltage in case 1.
Node
The 1st order cumulant
The 2nd order cumulant
The 3rd order cumulant
The 4th order cumulant
5
1.0178
0.3181
×
1
0
-
6
0
0
10
0.9893
0.3482
×
1
0
-
5
0
0
18
0.9786
0.8142
×
1
0
-
5
0
0
22
0.9756
0.8796
×
1
0
-
5
0
0Table 3
Moment of bus voltage in case 2.
Node
The 1st order cumulant
The 2nd order cumulant
The 3rd order cumulant
The 4th order cumulant
5
1.0211
0.5012
×
1
0
-
5
0.1928
×
1
0
-
10
-
0.1012
×
1
0
-
16
10
1.0097
0.6032
×
1
0
-
4
0.1074
×
1
0
-
8
-
0.1404
×
1
0
-
15
18
0.9676
0.1542
×
1
0
-
5
0.4143
×
1
0
-
8
-
0.1542
×
1
0
-
11
22
0.9458
0.1696
×
1
0
-
5
0.4387
×
1
0
-
8
-
0.4542
×
1
0
-
10Table 4
Moment of bus voltage in case 3.
Node
The 1st order cumulant
The 2nd order cumulant
The 3rd order cumulant
The 4th order cumulant
5
8.0207
0.5181
×
1
0
-
5
0.4085
×
1
0
-
11
-
0.1685
×
1
0
-
17
10
7.9887
0.4582
×
1
0
-
4
0.6875
×
1
0
-
10
-
0.4735
×
1
0
-
16
18
6.9693
0.1124
×
1
0
-
3
0.8655
×
1
0
-
9
-
0.4076
×
1
0
-
13
22
6.9876
0.1685
×
1
0
-
3
0.8725
×
1
0
-
9
-
0.1628
×
1
0
-
13Figure5 is the distribution of network loss; the calculation results of two methods are almost the same, but calculated amount of method 2 is far less than method 1. In other words, the method of this paper is effective.


## 4.1. Algorithm

Step 1.
Input feeder data,h
=
1.Step 2.
Read load and DG data ath level.Step 3.
Run power flow with Newton-Raphson ath level.Step 4.
Calculate the half invariance of the injection power flow of the generator and load ath level.Step 5.
Calculate the half invariance ofΔ
W.Step 6.
Calculate the half invariance ofΔ
X according to S
h
,
0.Step 7.
Calculate the PDF of variateΔ
X according to Gram-Charlier Series theorem.Step 8.
Calculate the half invariance ofΔ
Z according to T
h
,
0.Step 9.
Calculate the half invariance ofΔ
Loss according to M
h
,
0.Step 10.
Calculate the PDF ofΔ
Los
s
h according to Gram-Charlier series theorem at h level and store the data.Step 11.
According Step7, calculate the probability of out-of-limit voltage at h level and store the data.Step 12.
h
=
h
+
1; if h
⩽
24, then turn to Step 2; else turn to Step 13.Step 13.
According to the data stored in Step10, fit the PDF of Δ
Loss of one year.Step 14.
According to the data stored in Step11, estimate the time of out-of-limit voltage of one year.Step 15.
End.


## 4.2. Simulation Results

The proposed methodology is applied to a IEEE33 distribution network which is shown in Figure2 and joined with two distributed generations (i.e., node34 and node35).Figure 2
IEEE33 distribution network.The first one is the Monte Carlo simulation method and the second is the method proposed by us in this paper. The two methods, respectively, are used in three cases.Case
1 is not joined with any distributed generation (node34 and node35 are out of work); case
2 is joined with two wind turbine generations (node34 and node35), while case
3 is joined with one wind turbine generation (node34) and one photovoltaic generation (node35). The simulation results are as shown in Figures 4 and 5. Simulation parameters are given in Table 1. It is assumed that there are 24 demand levels in each year with equal duration of τ
h
=
365 h. The variations of demand level factors are depicted in Figure 3.Table 1
Data used in this paper.
Parameters
Unit
Value
c
8.78
V
min
⁡
pu
0.95
v
in
c
m/s
3
v
rate
m/s
13
ξ
L
0.925
ξ
max
⁡
1.15
ζ
L
0.9
ζ
max
⁡
1
k
1.75
V
max
⁡
pu
1.05
v
out
c
m/s
25
ξ
min
⁡
0.850
ξ
U
1.075
ζ
min
⁡
0
ζ
U
1
τ
h
h
365Figure 3
The variations ofDL
F
h in each demand level.Figure 4
Voltage probability density function of the 33rd node.Figure 5
Cumulative distribution function of total losses.Some cumulants are given in Tables1, 2, and 3. The first order cumulant is expected of node voltage; the second order cumulant is the variance. Figure 4 is the probability density function of the voltage amplitude of node33. From this figure, we can know that, in case
1, the random fluctuations of voltage are similar to normal distribution; probability of out-of-limit voltage is almost zero. When adding distributed generation, on one hand, line voltage had improved significantly and load node voltage rises. On the other hand, because of the randomness of the wind power and photovoltaic power, the node voltage fluctuation and the probability of out-of-limit voltage significantly increased. By comparing with case
2 and case
3, we can also find that, for wind and solar hybrid power systems, due to its complementarity, its impact on system voltage fluctuation is relatively smaller compared with single wind power and the probability of out-of-limit voltage is reduced obviously. The experimental results show that the proposed fuzzy variables are effective to describe load instead of normal variable (Table 4).Table 2
Moment of bus voltage in case 1.
Node
The 1st order cumulant
The 2nd order cumulant
The 3rd order cumulant
The 4th order cumulant
5
1.0178
0.3181
×
1
0
-
6
0
0
10
0.9893
0.3482
×
1
0
-
5
0
0
18
0.9786
0.8142
×
1
0
-
5
0
0
22
0.9756
0.8796
×
1
0
-
5
0
0Table 3
Moment of bus voltage in case 2.
Node
The 1st order cumulant
The 2nd order cumulant
The 3rd order cumulant
The 4th order cumulant
5
1.0211
0.5012
×
1
0
-
5
0.1928
×
1
0
-
10
-
0.1012
×
1
0
-
16
10
1.0097
0.6032
×
1
0
-
4
0.1074
×
1
0
-
8
-
0.1404
×
1
0
-
15
18
0.9676
0.1542
×
1
0
-
5
0.4143
×
1
0
-
8
-
0.1542
×
1
0
-
11
22
0.9458
0.1696
×
1
0
-
5
0.4387
×
1
0
-
8
-
0.4542
×
1
0
-
10Table 4
Moment of bus voltage in case 3.
Node
The 1st order cumulant
The 2nd order cumulant
The 3rd order cumulant
The 4th order cumulant
5
8.0207
0.5181
×
1
0
-
5
0.4085
×
1
0
-
11
-
0.1685
×
1
0
-
17
10
7.9887
0.4582
×
1
0
-
4
0.6875
×
1
0
-
10
-
0.4735
×
1
0
-
16
18
6.9693
0.1124
×
1
0
-
3
0.8655
×
1
0
-
9
-
0.4076
×
1
0
-
13
22
6.9876
0.1685
×
1
0
-
3
0.8725
×
1
0
-
9
-
0.1628
×
1
0
-
13Figure5 is the distribution of network loss; the calculation results of two methods are almost the same, but calculated amount of method 2 is far less than method 1. In other words, the method of this paper is effective.


## 5. Conclusions

A method combining half invariants and fuzzy mathematics theory is proposed for evaluation of active losses in the distribution network and the time of out-of-limit voltage. By comparing the simulation, it can be found that the proposed method can completely replace the Monte Carlo simulation method and, moreover, reduce a large amount of calculation. The model considers probabilistic presentation of wind speed using a Weibull PDF and probabilistic description of loads using normal distribution.On the other hand, as the reduction of calculation, this method can not only conveniently be used to calculate the network loss and the time of out-of-limit voltage and to assess the effects of a distributed generation to the distribution network, but also be used as the distributed power optimization index when considering size and site.



---

*Source: 939482-2014-08-21.xml*