# Probabilistic Harmonic Calculation in Distribution  Networks with Electric Vehicle Charging Stations

Authors: Jianxue Wang; Yanlin Cui; Minghui Zhu
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

Integrating EV charging station into power grid will bring impacts on power system, among which the most significant one is the harmonic pollution on distribution networks. Due to the uncertainty of the EV charging process, the harmonic currents brought by EV charging stations have a random nature. This paper proposed a mathematical simulation method for studying the working status of charging stations, which considers influencing factors including random leaving factor, electricity price, and waiting time. Based on the proposed simulation method, the probability distribution of the harmonic currents of EV charging stations is obtained and used in the calculation of the probability harmonic power flow. Then the impacts of EVs and EV charging stations on distribution networks can be analyzed. In the case study, the proposed simulation and analysis method is implemented on the IEEE-34 distribution network. The influences of EV arrival rates, the penetration rate, and the accessing location of EV charging station are also investigated. Results show that this research has good potential in guiding the planning and construction of charging station.


---

## Body

## 1. Introduction

Car industry has been developing for decades at a cracking pace. In order to pursue better performance and higher utility of vehicles, many valuable researches have been done to address the issues like adaptive control, mathematical modeling, and fault diagnosis [1–3]. In recent years, the pressure of fossil fuel shortage and other environmental problems leads to high-speed development of electric vehicles (EVs) due to their emission-free characteristic and effective utilization of renewable energy [4, 5]. The promising trend of EVs attracts plenty of researchers to concentrate on the functioning and the influence of electrical vehicles. Unlike conventional vehicles, the power supply of EVs mainly comes from local power distribution system through EV charging stations [6, 7] and that will lead to impacts on distribution network. In order to study the impacts, the EV charging station should be modeled properly and many factors need to be taken into consideration.The impacts of EV charging station integration can be divided into two categories basically. The first category contains the influences induced by uncoordinated charging of large amounts of EVs such as voltage deviations, thermal overloads, increase in losses, and decrease in voltage stability [6, 8–10]. Researches have investigated aspects including the different EV penetrations on distribution network investment [7], the long-term operating capacity reserve [11], and transformer life [12]. Coordinated charging strategies are proposed [12–15] to promote the distribution gird’s ability to hold more EVs without reinforcements, improve voltage profile, and reduce power losses.The other category is the impacts brought by high rating nonlinear switching devices in EV chargers [16] and the typical one in this category is harmonic pollution. Harmonic pollution will shorten the life span of transformers, reduce the utilization of transmission lines, and interfere with the signals of measuring and protection equipment. Therefore, studying the harmonic pollution brought by EV charging stations is of great importance. Paper [17] establishes a stochastic aggregate harmonic load model to study the harmonic characteristics of EV charger. Paper [18] presents a harmonic simulation method to study the harmonic impact of EVs. Paper [19] develops an approach to decrease the harmonics in a high EV penetration grid, and the harmonic load flow is used to model the distorted grid and identify the worst bus based on the typical daily load curve.In most of the calculation, the worst situation of harmonic of charging station was used. However, the probability of the worst situation is rare in daily life. In most of these studies the harmonic currents are generated by individual EV charger. The random characteristic of charging behavior was not considered in the calculation of charging station harmonic and was hardly consistent with the actual situation in our daily life. Besides, there are limited studies focused on the system level to evaluate the harmonic distortion of charging stations.In this paper, this problem is investigated from single equipment to the entire system. First, a simplified simulation model of EV charger is used to analyze the harmonic characteristics. Second, based on the queuing theory and charger model, a mathematical simulation method to get the probability distribution of charging station state is proposed, considering the influences of random factors, electricity price, and waiting time. The statistics of the simulation results are performed to obtain the probability distribution of total harmonics. Third, through the probability harmonic power flow, the probability characteristics of probability distribution of charging stations’ influences on system harmonics are analyzed. At last, the key factors, including the arrival rate, penetration rate, and the access location of EV charging station, were analyzed based on a modified IEEE-34 system.


## 2. The Working State Model of EV Charging Station

### 2.1. Model the EV Charger

In order to analyze the harmonic characteristic of EV charging station, the basic working model of a single EV charger needs to be built first. In this paper, we consider a kind of currently the most widely used EV charger. This typical EV charger consists of a three-phase uncontrolled rectifier and a high frequency DC/DC converter which are coupled with a high frequency transformer [4]. The output voltage and current of the charger are controlled by feedback control equipment in the high frequency DC/DC converter. Since elaborately modeling the EV charger is too complicated, it consumes much time for the simulation.To reduce the complexity of the EV charger model, a simplification proposed in [20] is adopted here, where the high frequency DC/DC converter and battery load are replaced by an equivalent resistance. The simplified equivalent circuit is shown in Figure 1.Figure 1
The structure diagram of equivalent simplified mode for EV charger [20].Based on Ohm’s law, the equivalent resistance can be calculated as(1)
R
(
t
)
=
U
I
I
I
(
t
)
=
U
I
2
P
I
(
t
)
=
η
U
I
2
P
O
(
t
)
,
where R is the equivalent resistance of high frequency DC/DC converter, P
I is the input power of the high frequency DC/DC converter, P
O is the output power of EV charger, η is the efficiency of high frequency DC/DC converter, U
I is the output voltage of three-phase uncontrolled rectifier, and I
I is the output current of three-phase uncontrolled rectifier.According to (1), the equivalent resistance should be a function of the charging time to reflect the different output power in different charging stages. The empirical formula of P
O during the charging process was obtained through fitting the actual-measured data [20], which can be put as
(2)
P
O
(
t
)
=
{
0.79
P
O
max
⁡
t
0.048
,
0
<
t
<
t
′
,
P
O
max
⁡
e
-
0.021
(
t
-
t
′
)
,
t
′
<
t
<
t
max
⁡
,
where P
O is the output power of EV charger, t is charging time, t
′ is the time spot when charging mode turns from constant current charging to constant voltage charging, t
max
⁡ is the maximum charging time, and P
O
max
⁡ is the rated output power.In [20], η was treated as a constant; however, according to the simulation result in MATLAB/Simulink environment, the efficiency of the high frequency DC/DC converter η varies with the output power of the EV charger P
O. Their relationship is presented in Figure 2.Figure 2
The relationship between efficiency and output power.Using the curve-fitting technique, the expression can be obtained as(3)
η
(
P
O
)
=
-
0.011
+
193
P
O
1
+
1.785
P
O
+
0.012
P
O
2
.Through (1)–(3), we can get the value of equivalent resistance at any time spot during charging.To illustrate the effectiveness of the simplified model in harmonic analysis, a comparison between the original and the simplified model (considering the varying efficiency during charging) has been made through simulating in MATLAB/Simulink environment. Simulation results show that the maximum error of the harmonic between the original and the simplified model is merely 0.45%, while the equivalent simplified model costs 40% less time than the original one. Therefore, the simplified model appears to be more effective and suitable, especially when repetitive calculations are required.


### 2.2. Build the Working State Model for an EV Charging Station

Analysis of the harmonic impact of EV charging stations requires knowing the harmonic currents injection of the charging stations to distribution network. With the EV charger simulation model mentioned in Section2.1, the harmonic currents that one charger brings to the grid can be obtained in simulation software. The harmonic currents of the entire charging station can be calculated by adding up the harmonic currents generated by every charger in work. Therefore, it is necessary to build the working state model for the EV charging station to acquire the number and the status of working chargers.To build the working state model, the charging process in the charging station, which depends on EVs’ charging behaviors, needs to be studied. Charging behaviors of EVs have much to do with the power conditions of their batteries. Normally, the power condition of a battery is measured by its SOC, which is defined as the available capacity in percentage of full charge capacity [21]. The battery power needs to be adequate to ensure the normal operation of electric vehicles. With the battery power consuming, the SOC will keep dropping and may get down to an overdischarging state [21]. When the battery power drops below a threshold value, the EV cannot work normally and thus needs to get charged at a nearby charging station. The charging process of an EV (as shown in Figure 3) can be modeled as a queuing process based on queuing theory (introduced in [22]), which can be summarized as “first come, first served.” The EV that needs charging adds to the tail of the queue when it enters the charging station. If there is no available EV charger, the EV will be waiting for an opening; if there is an available EV charger, it will get charged immediately. The EV will leave when the SOC reaches 100% or when there are some factors to make it leave early.Figure 3
The diagram of EVs’ charging process.In [23, 24], the queuing theory is used to describe the charging process and it is proved feasible. In [24], a typical queuing model, M/M/c, is used to portray the customer streams of EV charging stations. The typical queuing theory assumes that the inter-arrival and service time are generally subject to exponential distribution. As for EV charging, it is feasible to assume that the interarrival time of two consecutive EVs is exponentially distributed but it will lead to large error to consider that the charging time also obeys exponential distribution. The charging time is a random variable influenced mainly by the initial SOC, the random leaving factor, the electricity price, and the waiting time, so it does not seem to follow any already known probability distributions. Therefore, the typical queuing model may not describe the charging behavior effectively.In order to make the model more realistic, this paper adopts a numerical simulating method to study the EV charging process. First, the states of an EV charger should be modeled. The state of charger is determined by the SOC of the battery connected to it, where the range of SOC is from 0 to 100 in percentage. The battery charging process can be treated as a Markov process since the SOC of the next time spot is only related to the current SOC and has nothing to do with the SOC before, which can be put as follows:(4)
P
{
SOC
(
t
m
-
1
)
⊂
A
∣
SOC
(
t
1
)
=
SOC
1
,
SOC
(
t
2
)
=
SOC
2
,
…
,
SOC
(
t
m
)
=
SOC
m
}
=
P
{
SOC
(
t
m
+
1
)
⊂
A
∣
SOC
(
t
m
)
=
SOC
m
}
,
where t
m is the current time and SOC
(
t
m
) is the current state of charge.In order to make it easier for study, we discretize the domain of SOC’s possible values. Figure4 shows the transition of the states during charging. The SOC of the battery goes towards a higher level over the charging time if the EV does not quit at the middle of charging. If it does quit, the SOC of the battery goes to 100% directly. It is worth noting that 100% of SOC here means that the EV leaves and the charger is free. So there are two conditions that a state is treated as a 100% state. One is that the EV gets fully charged and leaves and the other is that the EV leaves in the middle of charging.Figure 4
The diagram of charging state transition.The diagram can be described in(5)
P
m
,
m
+
Δ
SOC
=
p
(
soc
)
=
1
-
q
(
soc
)
P
m
,
m
+
Δ
SOC
=
SOC
m
,
m
=
0
P
m
,
100
=
q
(
soc
)
,
(
0
<
m
<
100
)
.P
i
,
j means the transition probability from state i to j, Δ
SOC is the incremental power in a simulation time unit, and q
(
soc
) is the probability of quitting before the EV is fully charged, which will vary due to the current value of SOC.To start the simulation, the initial SOC should be given. Setting it to zero or some fixed low value will be inappropriate because the remaining SOC of last time is random. In fact, according to [25], the battery’s initial SOC obeys two-parameter β distribution, whose probability density function is
(6)
f
(
x
;
α
,
β
)
=
x
α
-
1
(
1
-
x
)
β
-
1
∫
0
1
u
α
-
1
(
1
-
u
)
β
-
1
d
u
=
Γ
(
α
+
β
)
Γ
(
α
)
Γ
(
β
)
x
α
-
1
(
1
-
x
)
β
-
1
=
1
B
(
α
,
β
)
x
α
-
1
(
1
-
x
)
β
-
1
.Quitting probabilityq
(
soc
) needs to be pinned down to complete the model for EV chargers. Generally, it contains 3 components, including random leaving factor, waiting time q
time, and electricity price q
price.The typical threshold valueSOC
′ is set to 80%. Leaving earlier before SOC reaches SOC
′ will do some harm to the battery, so people usually do not leave earlier no matter how long they will wait or how much the electricity costs. When SOC is equal to or greater than the threshold SOC
′, waiting time factor and electricity cost factor begin to work. The whole process can be illustrated by
(7)
q
(
soc
)
=
{
q
′
,
soc
<
SOC
′
,
q
′
+
q
time
+
q
price
,
soc
≥
SOC
′
.Furthermore,q
time is determined not only by the waiting time but also by the difference between the current SOC and the threshold SOC
′. With the increasing of SOC, both the waiting time and the difference to SOC
′ will increase, which leads to a higher waiting time-related quitting probability. In (8), q
time can be thought to be in a linear relationship with the product of the charging time t and the difference between the battery’s SOC and SOC
′ in this study, where time coefficient k
′ is a constant:
(8)
q
time
=
{
k
′
t
(
soc
-
SO
C
′
)
,
soc
≥
SOC
′
,
0
,
soc
<
SOC
′
.q
price is also decided by two components. One is the difference between the real-time electricity price (p
now) and the expected electricity price (P
E), and the other is the difference between SOC and SOC
′. The expression is shown in (9), where price coefficient k
′
′ is a constant:
(9)
q
price
=
{
k
′
′
(
p
now
-
p
E
)
(
soc
-
SOC
)
,
soc
≥
SO
C
′
,
p
now
>
P
E
,
0
,
soc
<
SO
C
′
.Finally, a correction should be made to improve the model of charging process. It is noted that an EV driver may be scared off by a long queue or get impatient waiting in line and then decide to leave before getting the car charged, so the effective arrival rate is a little bit less than the arrival rate. The correctional charging process is demonstrated in Figure5.Figure 5
The correctional diagram of EVs’ charging process.Considering the waiting time, the arrival rate is fixed in (10), where λ is the real arrival rate and λ
′ is the fixed arrival rate. The directly leaving rate (q
L) is the product of the leaving factor (k
wait) and the average waiting time (T
wait) which can be acquired in the simulation of queuing process:
(10)
λ
′
=
λ
-
q
L
=
λ
-
k
wait
T
wait
.Using the model given above to simulate the activity of each EV charger, the state of the whole charging station can be simulated by the sequential Monte Carlo method. This makes the preparation for the calculation of harmonic distribution.


### 2.3. Calculate the Harmonic Distribution of Charging Station

The harmonic currents of the whole charging station can be calculated by adding up the harmonic currents of all EV chargers. However, the states of chargers are uncertain so the harmonics of the charging station cannot be calculated directly. Using the working state model built in Section2.2, the states of chargers in an EV charging station can be determined at any time spot; then the simulation model of the entire charging station can be built in simulators to calculate the total harmonics. After simulation, the probability distribution of the charging station can be obtained through the statistical analysis of simulation results.The simulation model of the charging station is built by the simplified models of EV chargers mentioned in Section2.1. In each step of the EV charging station’s mathematical simulation, all the chargers’ output power is determined; thus, the equivalent resistance can be calculated through (1). After that, the value of equivalent resistance can be set in MATLAB/Simulink and then the harmonics of the charging station can be simulated for statistical analysis to acquire the probability distribution. Figure 6 illustrates the steps in a flowchart.Figure 6
Calculating the EV charging station’s harmonic distribution.Details of the steps in Figure6 are given as follows.(1)
Initialize the valid data, for example, the simulation step time (T) and the convergence criterion. Then start the mathematical simulation process.
(2)
Update the battery’s SOC in timeT through
(11)
Δ
SOC
(
t
)
=
∫
t
t
+
T
P
O
(
t
)
d
t
∫
0
t
all
P
O
(
t
)
d
t
×
100
%
,
where T is the time of simulation step, t is the battery’s charging time, and t
all is the total time of a battery’s SOC changing from 0 to 100 in percent.
(3)
Calculate the number of EVs which quit charging in timeT, considering the random leaving factor, the price of electricity, and the waiting time.
(4)
Calculate the number of EVs which arrive in timeT using the correctional arrival rate to increase the calculation accuracy. Then update the number of EV cars in charging station.
(5)
Give the initial SOC of each EV which arrives in timeT with (1). Then find the corresponding charging time (t) and output power (P
O).
(6)
Arrange the EV cars in the queue for charging according to the “first come, first served” principle.
(7)
Calculate the output power by (2) and (3), and then the value of equivalent resistance can be obtained by (1).
(8)
Set the resistance parameters of each charger in the simulation model of EV charging station. Then start the simulation and record results of harmonic and total power in an array.
(9)
Keep simulation until the probability distribution of charging state meets the convergence criterion that Euclidean distance of probability distributions obtained using different simulation time must be within the range of allowable error, which can be put as(12)
d
i
j
=
∑
i
=
1
M
(
x
i
k
-
x
j
k
)
2
<
ϵ
,
where x
i
k, x
j
k stands for the difference of the kth component of random variable x in the ith andjth distributions obtained in different simulation time, M stands for the number of variables, and ε stands for the allowable error.
(10)
Obtain the probability distribution of the charging state and the harmonics.


## 2.1. Model the EV Charger

In order to analyze the harmonic characteristic of EV charging station, the basic working model of a single EV charger needs to be built first. In this paper, we consider a kind of currently the most widely used EV charger. This typical EV charger consists of a three-phase uncontrolled rectifier and a high frequency DC/DC converter which are coupled with a high frequency transformer [4]. The output voltage and current of the charger are controlled by feedback control equipment in the high frequency DC/DC converter. Since elaborately modeling the EV charger is too complicated, it consumes much time for the simulation.To reduce the complexity of the EV charger model, a simplification proposed in [20] is adopted here, where the high frequency DC/DC converter and battery load are replaced by an equivalent resistance. The simplified equivalent circuit is shown in Figure 1.Figure 1
The structure diagram of equivalent simplified mode for EV charger [20].Based on Ohm’s law, the equivalent resistance can be calculated as(1)
R
(
t
)
=
U
I
I
I
(
t
)
=
U
I
2
P
I
(
t
)
=
η
U
I
2
P
O
(
t
)
,
where R is the equivalent resistance of high frequency DC/DC converter, P
I is the input power of the high frequency DC/DC converter, P
O is the output power of EV charger, η is the efficiency of high frequency DC/DC converter, U
I is the output voltage of three-phase uncontrolled rectifier, and I
I is the output current of three-phase uncontrolled rectifier.According to (1), the equivalent resistance should be a function of the charging time to reflect the different output power in different charging stages. The empirical formula of P
O during the charging process was obtained through fitting the actual-measured data [20], which can be put as
(2)
P
O
(
t
)
=
{
0.79
P
O
max
⁡
t
0.048
,
0
<
t
<
t
′
,
P
O
max
⁡
e
-
0.021
(
t
-
t
′
)
,
t
′
<
t
<
t
max
⁡
,
where P
O is the output power of EV charger, t is charging time, t
′ is the time spot when charging mode turns from constant current charging to constant voltage charging, t
max
⁡ is the maximum charging time, and P
O
max
⁡ is the rated output power.In [20], η was treated as a constant; however, according to the simulation result in MATLAB/Simulink environment, the efficiency of the high frequency DC/DC converter η varies with the output power of the EV charger P
O. Their relationship is presented in Figure 2.Figure 2
The relationship between efficiency and output power.Using the curve-fitting technique, the expression can be obtained as(3)
η
(
P
O
)
=
-
0.011
+
193
P
O
1
+
1.785
P
O
+
0.012
P
O
2
.Through (1)–(3), we can get the value of equivalent resistance at any time spot during charging.To illustrate the effectiveness of the simplified model in harmonic analysis, a comparison between the original and the simplified model (considering the varying efficiency during charging) has been made through simulating in MATLAB/Simulink environment. Simulation results show that the maximum error of the harmonic between the original and the simplified model is merely 0.45%, while the equivalent simplified model costs 40% less time than the original one. Therefore, the simplified model appears to be more effective and suitable, especially when repetitive calculations are required.


## 2.2. Build the Working State Model for an EV Charging Station

Analysis of the harmonic impact of EV charging stations requires knowing the harmonic currents injection of the charging stations to distribution network. With the EV charger simulation model mentioned in Section2.1, the harmonic currents that one charger brings to the grid can be obtained in simulation software. The harmonic currents of the entire charging station can be calculated by adding up the harmonic currents generated by every charger in work. Therefore, it is necessary to build the working state model for the EV charging station to acquire the number and the status of working chargers.To build the working state model, the charging process in the charging station, which depends on EVs’ charging behaviors, needs to be studied. Charging behaviors of EVs have much to do with the power conditions of their batteries. Normally, the power condition of a battery is measured by its SOC, which is defined as the available capacity in percentage of full charge capacity [21]. The battery power needs to be adequate to ensure the normal operation of electric vehicles. With the battery power consuming, the SOC will keep dropping and may get down to an overdischarging state [21]. When the battery power drops below a threshold value, the EV cannot work normally and thus needs to get charged at a nearby charging station. The charging process of an EV (as shown in Figure 3) can be modeled as a queuing process based on queuing theory (introduced in [22]), which can be summarized as “first come, first served.” The EV that needs charging adds to the tail of the queue when it enters the charging station. If there is no available EV charger, the EV will be waiting for an opening; if there is an available EV charger, it will get charged immediately. The EV will leave when the SOC reaches 100% or when there are some factors to make it leave early.Figure 3
The diagram of EVs’ charging process.In [23, 24], the queuing theory is used to describe the charging process and it is proved feasible. In [24], a typical queuing model, M/M/c, is used to portray the customer streams of EV charging stations. The typical queuing theory assumes that the inter-arrival and service time are generally subject to exponential distribution. As for EV charging, it is feasible to assume that the interarrival time of two consecutive EVs is exponentially distributed but it will lead to large error to consider that the charging time also obeys exponential distribution. The charging time is a random variable influenced mainly by the initial SOC, the random leaving factor, the electricity price, and the waiting time, so it does not seem to follow any already known probability distributions. Therefore, the typical queuing model may not describe the charging behavior effectively.In order to make the model more realistic, this paper adopts a numerical simulating method to study the EV charging process. First, the states of an EV charger should be modeled. The state of charger is determined by the SOC of the battery connected to it, where the range of SOC is from 0 to 100 in percentage. The battery charging process can be treated as a Markov process since the SOC of the next time spot is only related to the current SOC and has nothing to do with the SOC before, which can be put as follows:(4)
P
{
SOC
(
t
m
-
1
)
⊂
A
∣
SOC
(
t
1
)
=
SOC
1
,
SOC
(
t
2
)
=
SOC
2
,
…
,
SOC
(
t
m
)
=
SOC
m
}
=
P
{
SOC
(
t
m
+
1
)
⊂
A
∣
SOC
(
t
m
)
=
SOC
m
}
,
where t
m is the current time and SOC
(
t
m
) is the current state of charge.In order to make it easier for study, we discretize the domain of SOC’s possible values. Figure4 shows the transition of the states during charging. The SOC of the battery goes towards a higher level over the charging time if the EV does not quit at the middle of charging. If it does quit, the SOC of the battery goes to 100% directly. It is worth noting that 100% of SOC here means that the EV leaves and the charger is free. So there are two conditions that a state is treated as a 100% state. One is that the EV gets fully charged and leaves and the other is that the EV leaves in the middle of charging.Figure 4
The diagram of charging state transition.The diagram can be described in(5)
P
m
,
m
+
Δ
SOC
=
p
(
soc
)
=
1
-
q
(
soc
)
P
m
,
m
+
Δ
SOC
=
SOC
m
,
m
=
0
P
m
,
100
=
q
(
soc
)
,
(
0
<
m
<
100
)
.P
i
,
j means the transition probability from state i to j, Δ
SOC is the incremental power in a simulation time unit, and q
(
soc
) is the probability of quitting before the EV is fully charged, which will vary due to the current value of SOC.To start the simulation, the initial SOC should be given. Setting it to zero or some fixed low value will be inappropriate because the remaining SOC of last time is random. In fact, according to [25], the battery’s initial SOC obeys two-parameter β distribution, whose probability density function is
(6)
f
(
x
;
α
,
β
)
=
x
α
-
1
(
1
-
x
)
β
-
1
∫
0
1
u
α
-
1
(
1
-
u
)
β
-
1
d
u
=
Γ
(
α
+
β
)
Γ
(
α
)
Γ
(
β
)
x
α
-
1
(
1
-
x
)
β
-
1
=
1
B
(
α
,
β
)
x
α
-
1
(
1
-
x
)
β
-
1
.Quitting probabilityq
(
soc
) needs to be pinned down to complete the model for EV chargers. Generally, it contains 3 components, including random leaving factor, waiting time q
time, and electricity price q
price.The typical threshold valueSOC
′ is set to 80%. Leaving earlier before SOC reaches SOC
′ will do some harm to the battery, so people usually do not leave earlier no matter how long they will wait or how much the electricity costs. When SOC is equal to or greater than the threshold SOC
′, waiting time factor and electricity cost factor begin to work. The whole process can be illustrated by
(7)
q
(
soc
)
=
{
q
′
,
soc
<
SOC
′
,
q
′
+
q
time
+
q
price
,
soc
≥
SOC
′
.Furthermore,q
time is determined not only by the waiting time but also by the difference between the current SOC and the threshold SOC
′. With the increasing of SOC, both the waiting time and the difference to SOC
′ will increase, which leads to a higher waiting time-related quitting probability. In (8), q
time can be thought to be in a linear relationship with the product of the charging time t and the difference between the battery’s SOC and SOC
′ in this study, where time coefficient k
′ is a constant:
(8)
q
time
=
{
k
′
t
(
soc
-
SO
C
′
)
,
soc
≥
SOC
′
,
0
,
soc
<
SOC
′
.q
price is also decided by two components. One is the difference between the real-time electricity price (p
now) and the expected electricity price (P
E), and the other is the difference between SOC and SOC
′. The expression is shown in (9), where price coefficient k
′
′ is a constant:
(9)
q
price
=
{
k
′
′
(
p
now
-
p
E
)
(
soc
-
SOC
)
,
soc
≥
SO
C
′
,
p
now
>
P
E
,
0
,
soc
<
SO
C
′
.Finally, a correction should be made to improve the model of charging process. It is noted that an EV driver may be scared off by a long queue or get impatient waiting in line and then decide to leave before getting the car charged, so the effective arrival rate is a little bit less than the arrival rate. The correctional charging process is demonstrated in Figure5.Figure 5
The correctional diagram of EVs’ charging process.Considering the waiting time, the arrival rate is fixed in (10), where λ is the real arrival rate and λ
′ is the fixed arrival rate. The directly leaving rate (q
L) is the product of the leaving factor (k
wait) and the average waiting time (T
wait) which can be acquired in the simulation of queuing process:
(10)
λ
′
=
λ
-
q
L
=
λ
-
k
wait
T
wait
.Using the model given above to simulate the activity of each EV charger, the state of the whole charging station can be simulated by the sequential Monte Carlo method. This makes the preparation for the calculation of harmonic distribution.


## 2.3. Calculate the Harmonic Distribution of Charging Station

The harmonic currents of the whole charging station can be calculated by adding up the harmonic currents of all EV chargers. However, the states of chargers are uncertain so the harmonics of the charging station cannot be calculated directly. Using the working state model built in Section2.2, the states of chargers in an EV charging station can be determined at any time spot; then the simulation model of the entire charging station can be built in simulators to calculate the total harmonics. After simulation, the probability distribution of the charging station can be obtained through the statistical analysis of simulation results.The simulation model of the charging station is built by the simplified models of EV chargers mentioned in Section2.1. In each step of the EV charging station’s mathematical simulation, all the chargers’ output power is determined; thus, the equivalent resistance can be calculated through (1). After that, the value of equivalent resistance can be set in MATLAB/Simulink and then the harmonics of the charging station can be simulated for statistical analysis to acquire the probability distribution. Figure 6 illustrates the steps in a flowchart.Figure 6
Calculating the EV charging station’s harmonic distribution.Details of the steps in Figure6 are given as follows.(1)
Initialize the valid data, for example, the simulation step time (T) and the convergence criterion. Then start the mathematical simulation process.
(2)
Update the battery’s SOC in timeT through
(11)
Δ
SOC
(
t
)
=
∫
t
t
+
T
P
O
(
t
)
d
t
∫
0
t
all
P
O
(
t
)
d
t
×
100
%
,
where T is the time of simulation step, t is the battery’s charging time, and t
all is the total time of a battery’s SOC changing from 0 to 100 in percent.
(3)
Calculate the number of EVs which quit charging in timeT, considering the random leaving factor, the price of electricity, and the waiting time.
(4)
Calculate the number of EVs which arrive in timeT using the correctional arrival rate to increase the calculation accuracy. Then update the number of EV cars in charging station.
(5)
Give the initial SOC of each EV which arrives in timeT with (1). Then find the corresponding charging time (t) and output power (P
O).
(6)
Arrange the EV cars in the queue for charging according to the “first come, first served” principle.
(7)
Calculate the output power by (2) and (3), and then the value of equivalent resistance can be obtained by (1).
(8)
Set the resistance parameters of each charger in the simulation model of EV charging station. Then start the simulation and record results of harmonic and total power in an array.
(9)
Keep simulation until the probability distribution of charging state meets the convergence criterion that Euclidean distance of probability distributions obtained using different simulation time must be within the range of allowable error, which can be put as(12)
d
i
j
=
∑
i
=
1
M
(
x
i
k
-
x
j
k
)
2
<
ϵ
,
where x
i
k, x
j
k stands for the difference of the kth component of random variable x in the ith andjth distributions obtained in different simulation time, M stands for the number of variables, and ε stands for the allowable error.
(10)
Obtain the probability distribution of the charging state and the harmonics.


## 3. Probabilistic Harmonic Power Flow of Distribution System Integrating EVs

The influence of integrating EVs into the distribution network has a probabilistic nature due to the random and uncertain characteristic of EVs. Based on the probability distribution of total harmonic of charging stations, the probability characteristics of charging stations’ influences on system harmonics can be analyzed through the probabilistic harmonic power flow [10, 26]. From the simulation results of harmonics, we find that there is no typical curve that can describe the probability characteristics of charging stations’ harmonic. Thus in this research, the Monte Carlo sampling is used to calculate the probabilistic harmonic power flow.Because the harmonic power brought by EV charging station is much less than the fundamental power in distribution network, the fundamental power flow and the harmonic power flow can be decoupled and calculated, respectively. Newton-Raphson method was used to calculate the fundamental power flow.The harmonic models of load, transmission line, and transformer adopt the typical harmonic models in [27]. Based on the harmonic model of electrical components and the result of fundamental power flow, n time harmonic admittance matrix (Y
n) can be obtained according to the harmonic analysis theory. Harmonic sources are treated as an ideal current source, so the n time injection current matrix (I
n) can be obtained. The n time harmonic voltage matrix (U
n) can be obtained by solving
(13)
Y
n
U
n
=
I
n
.After solving harmonic power flow at each time, the effective voltage (U
i) and total harmonic distortion of voltage (THD
v) can be calculated through
(14)
U
i
=
∑
n
=
1
H
|
U
1
n
|
2
,
THD
v
=
[
(
∑
n
=
2
H
|
U
i
n
|
2
)
0.5
|
U
i
1
|
]
×
100
%
.After the calculation of probabilistic harmonic power flow based on Monte Carlo, the probability distribution of effective voltage andTHD
v can be acquired. Then the probability characteristics of charging stations’ influences on system harmonics under different conditions can be analyzed.


## 4. Case Study

The case study is implemented on a modified IEEE-34 test system [28], which is shown in Figure 7.Figure 7
The diagram of the modified IEEE-34 test system.An EV charging station is connected to the original IEEE-34 test system at Node 23 through a transformer. The penetration power of the EV charging station is set to 15% in this case. Detailed information of the charging station is listed in Table1.Table 1
The detailed information of the charging station.
Voltage of network/kV
The number of transformers
The ratio of transformation
The number of chargers
The total charging power/kW
The number of parking spaces
10
1
10 kV/400 V
35
315
50The charging time of the EV battery from out-of-power to fully charged is 270 minutes. The parameters of the EV charger are listed in Table2.Table 2
Parameters of the EV charger.
Voltage/V
Rated output power/kW
Output voltage/V
Output current/A
L
f
(mH)
C
f
(uF)
400
9
300
30
1.5
2115The arrival rate is set to 0.1 during the day and 0.2 at night since fewer people use cars at night and the electricity price is lower. According to the load curve, the load level during the night is lower than other times. The random leaving factor, the waiting time, and electricity price are considered during the daytime. The parameters of the working state model are listed in Table3.Table 3
Parameters of the working state model.
Parameters
α
β
q
k/10−4
k
′
′
k
wait/10−4
Value
2.7
5.5
0.004
1.25
0.2
1


### 4.1. The Efficiency Analysis of Working State Model

In order to prove the effectiveness of the model proposed in this paper, the probability distribution of charging station’s state is calculated using two different models. The first is the conventional queuing theory model, in which the average service time of charger is 203 minutes. For in the typical model, the random leaving factor, the waiting time, and the price of electricity are not reflected. So in order to prove the effectiveness of the method proposed in this paper, these factors are ignored. The probability distribution of EV charging station’s state obtained using two models is shown in Figure8.Figure 8
The comparison of probability distribution of EV charger in work using typical and proposed model.From Figure8, it is obvious that the probability distribution obtained by the typical model is similar to that by the proposed model.The parameters of probability distribution of station state are listed in Table4.Table 4
Parameters of the probability distribution of station state.
Model
Expected number of chargers in work
The number of chargers with highest possibility
The probability of all chargers in work
The average length of the queue
Typical
20.3
20
0.002
20.3
Proposed model
20.3
20
0.003
20.4Among the parameters listed in Table4, the expected number of working chargers and the number of chargers with the highest probability are the same, while the other two parameters are similar in two models. The simulation step time (T) is set to be 3 minutes in this case. If it is set smaller, the value of those similar parameters will get closer. Therefore, the method proposed in this paper can be used to simulate the charging process and calculate the probability distribution of charging station’s harmonic.


### 4.2. Result of Harmonic Probability Power Flow

Based on the working state model of charging station, the distribution of charging station’s total power and harmonics can be acquired by doing statistics on simulation results. The probability distributions of charging station’s total power during daytime and night are shown in Figures9 and 10, respectively.Figure 9
The probability distribution of charging station’s power during daytime.Figure 10
The probability distribution of charging station’s power during night.From Figures9 and 10, it can be noted that the EV charging station is more heavily utilized during night than during daytime.The distribution of harmonics of charging station can also be acquired. Based on the distribution, the expectation of the current amplitude, the harmonic content, and the expectation of the current angle during day and night are calculated, which are listed in Table5. It is worth noting that only 6
k
±
1th (k
=
1,2
,
…) harmonics are listed because they are the major harmonics of EV charging station.Table 5
The detailed information of charging station’s harmonics.
Time
Parameters
Order
5
7
11
13
17
19
Daytime
The expectation of current amplitude/10−2 p.u.
3.77
1.61
0.82
0.44
0.41
0.25
The harmonic content/%
34.2
14.6
7.5
4.0
3.7
2.3
Night
The expectation of current amplitude/10−2 p.u.
6.45
2.13
1.47
0.74
0.66
0.45
The harmonic content/%
31.5
10.4
7.2
3.6
3.3
2.2The expectation of current amplitude decreases with the increase of the harmonic order. It is easy to see that the low-order harmonic current accounts for larger proportion than the high-order harmonic. The harmonics with orders higher than 19 are not listed because the amplitudes are too small; thus, the influences are neglectable.Based on the result of the probabilistic harmonic power flow, the probability distributions of root mean square (RMS) voltages of each node are obtained. The overranging nodes are listed in Table6.Table 6
The detailed information of voltage RMS over the range.
Node
8
9
10
11
12
13
22
23
24
25
26
The probability of voltage RMS over the range/10−2
0.01
0.27
4.88
5.50
6.73
7.70
0.03
0.03
2.07
0.27
13.69According to Table6, there are 11 nodes whose RMS voltages are overrange. Among these nodes, Node 26 has the highest overranging probability. The probability distribution of voltage RMS of Node 26 is shown in Figure 11.Figure 11
The probability distribution of voltage in node 26.It can be seen that the distribution has a bimodal characteristic, which is due to the peak-valley load during 24 hours.Based on the result of the probabilistic harmonic power flow, the expectedTH
D
v of each node is calculated, which is shown in Figure 12.Figure 12
The expected voltage distortion of each node.From Figure12, it can be observed that Node 23, the EV charging station node, and its neighbor node, Node 24, have very high expected TH
D
v. Generally, it is obvious that the expected TH
D
v of the node near Node 23 is high while the expected TH
D
v of the node far from Node 23 is low.The probability distribution ofTHD
v of Node 23 is shown in Figure 13.Figure 13
The probability distribution of voltage distortion in node 23.Table7 lists the nodes whose TH
D
v is over range. There are 11 nodes whose probability of THD
v is over the range; 7 of them have an overranging probability that is higher than 49%. The result shows that the integration of EV charging station significantly affects the safety and power quality of the distribution network. So filter devices need to be installed to diminish the influence.Table 7
The detailed information ofTHD
v over the range.
Node
5
6
7
8
9
10
11
22
23
24
25
The probability ofTHD
v over the range/%
0.4
29.2
49.2
53.6
49.2
29.6
3.8
68.1
76.1
76.1
49.2


### 4.3. Influence Factors of Charging Station

In this part of research, the influences of the arrival rate, the penetration rate, and the accessing location of EV charging station on the probabilistic harmonic power flow are analyzed as follows.(1) The Arrival Rate. First, the influence of arrival rate is studied. Calculations have been made when arrival rate is set to different values, 0.1 and 0.2. The results of the probabilistic harmonic power flow during daytime are listed in Table 8.Table 8
The detailed information ofTHD
v over the range in different arrival rate.
Arrival rate
Node
6
7
8
9
10
22
23
24
25
0.1
The probability ofTHD
v over the range/%
0.9
5.5
9.3
3.6
0.9
25.8
53.9
53.9
3.7
Arrival rate
Node
5
6
7
8
9
10
11
22
23
24
25
0.2
The probability ofTHD
v over the range/%
0.8
43.8
93.8
96.6
89.9
45.6
3.4
100
100
100
89.9With the increase of the arrival rate, the harmonic distortion of the network becomes worse: two new nodes, Node 5 and Node 11, becomeTHD
v overranging nodes. The influenced area becomes larger due to the increase of the arrival rate and so does the THD
v overranging probability. In fact, the overranging probability even reaches 100% at Nodes 22, 23, and 24. Therefore, it is easy to conclude that the harmonic pollution is more significant when the arrival rate increases.(2) The Penetration Rate. The harmonic impact on the distribution network is studied during 24 hours when the penetration rate increases to 30%. The results of probabilistic harmonic power flow are listed in Table 9.Table 9
The detailed information ofTHD
v over the range when penetration rate is 30%.
Node
4
5
6
7
8
9
10
11
The probability ofTHD
v over the range
0.048
0.454
0.700
0.889
0.889
0.889
0.747
0.668
Node
12
13
21
22
23
24
25
The probability ofTHD
v over the range
0.239
0.187
0.048
0.982
0.994
0.994
0.889Compared to the results in Table7, the number of THD
v overranging nodes has increased by 4, the influence area is larger than the original result, and the probability of THD
v over the range changes from 0.761 to 0.994 at Node 23. It can be seen that the impact on the power quality of the distribution network will be more significant when the penetration rate increases.(3) Accessing Location. Situations where the EV charging station accesses Nodes 10, 23, and 34 are compared. The results are listed in Table 10.Table 10
The detailed information ofTHD
v over the range in different accessing location.
Accessing location
Node
7
8
9
10
11
12
13
22
23
24
25
26
10
The probability ofTHD
v over the range/%
1.4
1.4
3.2
8
5.1
1.4
0.6
1.4
1.6
1.6
3.2
1.4
Accessing location
Node
6
7
8
9
10
22
23
24
25
23
The probability ofTHD
v over the range/%
0.9
5.5
9.3
3.6
0.9
25.8
53.9
53.9
3.7
Accessing location
Node
18
19
20
34
34
The probability ofTHD
v over the range/%
1.4
1.4
1.4
23.7From the results, two conclusions which are useful in EV charging station planning can be drawn. The influenced area when the EV charging station accesses to Node 10 is larger than that of Node 23. It can be interpreted in a way that Node 10 is the main feeder of distribution network so accessing it will result in a larger influenced area compared with accessing a branch feeder. When accessing location is Node 34, the influenced area and the degree of influence are smaller than that when accessing the other two nodes. That is because Node 34 is close to Node 15, which is connected with reactive power compensator that has an inhibitory effect on harmonics. Therefore, it can be concluded that it is better to set the charging station close to reactive power compensation devices or away from the main feeder.


## 4.1. The Efficiency Analysis of Working State Model

In order to prove the effectiveness of the model proposed in this paper, the probability distribution of charging station’s state is calculated using two different models. The first is the conventional queuing theory model, in which the average service time of charger is 203 minutes. For in the typical model, the random leaving factor, the waiting time, and the price of electricity are not reflected. So in order to prove the effectiveness of the method proposed in this paper, these factors are ignored. The probability distribution of EV charging station’s state obtained using two models is shown in Figure8.Figure 8
The comparison of probability distribution of EV charger in work using typical and proposed model.From Figure8, it is obvious that the probability distribution obtained by the typical model is similar to that by the proposed model.The parameters of probability distribution of station state are listed in Table4.Table 4
Parameters of the probability distribution of station state.
Model
Expected number of chargers in work
The number of chargers with highest possibility
The probability of all chargers in work
The average length of the queue
Typical
20.3
20
0.002
20.3
Proposed model
20.3
20
0.003
20.4Among the parameters listed in Table4, the expected number of working chargers and the number of chargers with the highest probability are the same, while the other two parameters are similar in two models. The simulation step time (T) is set to be 3 minutes in this case. If it is set smaller, the value of those similar parameters will get closer. Therefore, the method proposed in this paper can be used to simulate the charging process and calculate the probability distribution of charging station’s harmonic.


## 4.2. Result of Harmonic Probability Power Flow

Based on the working state model of charging station, the distribution of charging station’s total power and harmonics can be acquired by doing statistics on simulation results. The probability distributions of charging station’s total power during daytime and night are shown in Figures9 and 10, respectively.Figure 9
The probability distribution of charging station’s power during daytime.Figure 10
The probability distribution of charging station’s power during night.From Figures9 and 10, it can be noted that the EV charging station is more heavily utilized during night than during daytime.The distribution of harmonics of charging station can also be acquired. Based on the distribution, the expectation of the current amplitude, the harmonic content, and the expectation of the current angle during day and night are calculated, which are listed in Table5. It is worth noting that only 6
k
±
1th (k
=
1,2
,
…) harmonics are listed because they are the major harmonics of EV charging station.Table 5
The detailed information of charging station’s harmonics.
Time
Parameters
Order
5
7
11
13
17
19
Daytime
The expectation of current amplitude/10−2 p.u.
3.77
1.61
0.82
0.44
0.41
0.25
The harmonic content/%
34.2
14.6
7.5
4.0
3.7
2.3
Night
The expectation of current amplitude/10−2 p.u.
6.45
2.13
1.47
0.74
0.66
0.45
The harmonic content/%
31.5
10.4
7.2
3.6
3.3
2.2The expectation of current amplitude decreases with the increase of the harmonic order. It is easy to see that the low-order harmonic current accounts for larger proportion than the high-order harmonic. The harmonics with orders higher than 19 are not listed because the amplitudes are too small; thus, the influences are neglectable.Based on the result of the probabilistic harmonic power flow, the probability distributions of root mean square (RMS) voltages of each node are obtained. The overranging nodes are listed in Table6.Table 6
The detailed information of voltage RMS over the range.
Node
8
9
10
11
12
13
22
23
24
25
26
The probability of voltage RMS over the range/10−2
0.01
0.27
4.88
5.50
6.73
7.70
0.03
0.03
2.07
0.27
13.69According to Table6, there are 11 nodes whose RMS voltages are overrange. Among these nodes, Node 26 has the highest overranging probability. The probability distribution of voltage RMS of Node 26 is shown in Figure 11.Figure 11
The probability distribution of voltage in node 26.It can be seen that the distribution has a bimodal characteristic, which is due to the peak-valley load during 24 hours.Based on the result of the probabilistic harmonic power flow, the expectedTH
D
v of each node is calculated, which is shown in Figure 12.Figure 12
The expected voltage distortion of each node.From Figure12, it can be observed that Node 23, the EV charging station node, and its neighbor node, Node 24, have very high expected TH
D
v. Generally, it is obvious that the expected TH
D
v of the node near Node 23 is high while the expected TH
D
v of the node far from Node 23 is low.The probability distribution ofTHD
v of Node 23 is shown in Figure 13.Figure 13
The probability distribution of voltage distortion in node 23.Table7 lists the nodes whose TH
D
v is over range. There are 11 nodes whose probability of THD
v is over the range; 7 of them have an overranging probability that is higher than 49%. The result shows that the integration of EV charging station significantly affects the safety and power quality of the distribution network. So filter devices need to be installed to diminish the influence.Table 7
The detailed information ofTHD
v over the range.
Node
5
6
7
8
9
10
11
22
23
24
25
The probability ofTHD
v over the range/%
0.4
29.2
49.2
53.6
49.2
29.6
3.8
68.1
76.1
76.1
49.2


## 4.3. Influence Factors of Charging Station

In this part of research, the influences of the arrival rate, the penetration rate, and the accessing location of EV charging station on the probabilistic harmonic power flow are analyzed as follows.(1) The Arrival Rate. First, the influence of arrival rate is studied. Calculations have been made when arrival rate is set to different values, 0.1 and 0.2. The results of the probabilistic harmonic power flow during daytime are listed in Table 8.Table 8
The detailed information ofTHD
v over the range in different arrival rate.
Arrival rate
Node
6
7
8
9
10
22
23
24
25
0.1
The probability ofTHD
v over the range/%
0.9
5.5
9.3
3.6
0.9
25.8
53.9
53.9
3.7
Arrival rate
Node
5
6
7
8
9
10
11
22
23
24
25
0.2
The probability ofTHD
v over the range/%
0.8
43.8
93.8
96.6
89.9
45.6
3.4
100
100
100
89.9With the increase of the arrival rate, the harmonic distortion of the network becomes worse: two new nodes, Node 5 and Node 11, becomeTHD
v overranging nodes. The influenced area becomes larger due to the increase of the arrival rate and so does the THD
v overranging probability. In fact, the overranging probability even reaches 100% at Nodes 22, 23, and 24. Therefore, it is easy to conclude that the harmonic pollution is more significant when the arrival rate increases.(2) The Penetration Rate. The harmonic impact on the distribution network is studied during 24 hours when the penetration rate increases to 30%. The results of probabilistic harmonic power flow are listed in Table 9.Table 9
The detailed information ofTHD
v over the range when penetration rate is 30%.
Node
4
5
6
7
8
9
10
11
The probability ofTHD
v over the range
0.048
0.454
0.700
0.889
0.889
0.889
0.747
0.668
Node
12
13
21
22
23
24
25
The probability ofTHD
v over the range
0.239
0.187
0.048
0.982
0.994
0.994
0.889Compared to the results in Table7, the number of THD
v overranging nodes has increased by 4, the influence area is larger than the original result, and the probability of THD
v over the range changes from 0.761 to 0.994 at Node 23. It can be seen that the impact on the power quality of the distribution network will be more significant when the penetration rate increases.(3) Accessing Location. Situations where the EV charging station accesses Nodes 10, 23, and 34 are compared. The results are listed in Table 10.Table 10
The detailed information ofTHD
v over the range in different accessing location.
Accessing location
Node
7
8
9
10
11
12
13
22
23
24
25
26
10
The probability ofTHD
v over the range/%
1.4
1.4
3.2
8
5.1
1.4
0.6
1.4
1.6
1.6
3.2
1.4
Accessing location
Node
6
7
8
9
10
22
23
24
25
23
The probability ofTHD
v over the range/%
0.9
5.5
9.3
3.6
0.9
25.8
53.9
53.9
3.7
Accessing location
Node
18
19
20
34
34
The probability ofTHD
v over the range/%
1.4
1.4
1.4
23.7From the results, two conclusions which are useful in EV charging station planning can be drawn. The influenced area when the EV charging station accesses to Node 10 is larger than that of Node 23. It can be interpreted in a way that Node 10 is the main feeder of distribution network so accessing it will result in a larger influenced area compared with accessing a branch feeder. When accessing location is Node 34, the influenced area and the degree of influence are smaller than that when accessing the other two nodes. That is because Node 34 is close to Node 15, which is connected with reactive power compensator that has an inhibitory effect on harmonics. Therefore, it can be concluded that it is better to set the charging station close to reactive power compensation devices or away from the main feeder.


## 5. Conclusions

A new method for studying the working status of EV charging stations is proposed in this paper, while considering the random leaving factor, the electricity price, and the waiting time. This method is proved to be effective through the comparison with the conventional queuing theory model. In this paper, this method is used to analyze the characteristics of harmonics from the charging station. The probabilistic harmonic power flow is used to analyze the impact of charging station’s harmonics on the distribution network. The proposed model can be directly applied in some other fields that involve charging station’s state simulation. In addition, other factors can also be considered in the simulation process. Since the method proposed in this paper is a simulation process, further study will focus on the mathematical theory of the charging behavior.



---

*Source: 167565-2014-08-14.xml*