Energy Management Method for Fast-Charging Stations with the Energy Storage System to Alleviate the Voltage Problem of the Observation Node

Abstract

Large-scale fast charging of electric vehicles (EVs) probably causes voltage deviation problems in the distribution network. Installing energy storage systems (ESSs) in the fast-charging stations (FCSs) and formulating appropriate active power plans for ESSs is an effective way to reduce the local voltage deviation problem. Some deterministic centralized strategies used for ESSs at FCSs are proposed to solve the voltage deviation problem mentioned above. However, the randomness of the EV load is very large, which can probably reduce the effects of deterministic centralized strategies. A fast and reliable centralized strategy considering the randomness of the EV load for ESSs is a key requirement. Therefore, we propose in this paper a day-ahead scheduling strategy with the aim of maximizing the probability of the nodal voltage change being smaller than a preset limit at the observation node. In the proposed strategy, the uncertainty of EV load is taken into account and the probability of the voltage change of an observation node is quantified by a proposed analytic assessment model (AMM). Furthermore, a voltage change optimization model (VCOM) based on a novel control parameter β is proposed, where β can be used as a constraint to suppress the nodal voltage change at the observation node. Finally, the IEEE 33-bus test system is used to verify the effectiveness of the proposed day-ahead ESS strategy.

1. Introduction

The adoption of EVs is regarded as having very large potential in terms of reducing carbon emissions [1,2]. The take-up rate of EVs has increased rapidly in recent years, the number of EVs registered worldwide by 2018 having surpassed 5 million [3]. The high penetration level of EVs can bring on a sharp increase in power load, which may cause adverse impacts on the distribution grid [4], especially fast-charging EVs [5]. During peak load periods, the distribution grid may easily suffer security risks such as transformer overload, voltage problems, or cable overload [6].

The voltage problems caused by increasing EV charging load are widely seen as concerning. Previous studies have confirmed that high penetrations of EVs can reduce the voltage stability margin and may lead to the voltage instability of the power grid [7]. A similar conclusion was obtained by using the IEEE 3-bus test system [8]. Since there is very little voltage instability in the distribution networks operating under normal conditions, more research focus is required on the impacts of EV charging on voltage change. The randomness and impulsive characteristics of EV recharging loads mean that the injection of power at the charging node can vary greatly, which may increase the nodal voltage change at the local node or neighboring nodes. The findings in [9] show that voltage violation may occur at 10–20% EV penetrations, and in [10] the voltage deviation reaches 3.46% at 20% EV penetrations. Furthermore, large-scale EV fast charging also can cause voltage flicker in the distribution network [11].

Many studies show that installing energy storage systems (ESS) in fast-charging stations with an appropriate energy management strategy can alleviate voltage problems. A compensation control strategy for the ESS is proposed to deal with energy imbalance in continuous operation [12]. An advanced ESS control strategy based on fuzzy control is used to reduce the local voltage fluctuation [13]. The authors of Ref. [14] studied the ESS used in the FCS, which demonstrates that the voltage profile of the transformer in the FCS can be improved by the ESS. A local control algorithm for the FCS with flywheel energy storage systems was proposed in [15], which can also guarantee the charging demand power of the EV during the period of charging control. However, the randomness of the EV load is not considered in these studies. Furthermore, the actual adjustable capacity of the ESS in an FCS is very limited so that it is generally less effective despite using the above local voltage regulation methods to solve the voltage problem in the distribution network.

To address these inadequacies, a centralized day-ahead scheduling strategy for ESSs in FCSs is proposed in this paper to reduce the adverse impact of FCSs on the voltage changes at observation nodes in the radial distribution network. This paper assumes that a nodal voltage regulation system has been established between the distribution network operator and the FCS operator to reduce the adverse impact caused by FCSs, which includes a regional-level optimization decision-making system (RODS) and station-level energy management systems (SEMSs), as shown in Figure 1. RODS is responsible for formulating the day-ahead active power plan for the ESS. SEMS i is responsible for executing the day-ahead active power plan i.

The rest of the paper is arranged as follows. Section 2 presents the problem formulation. Section 3 expounds on the AMM and the control characteristics of the control parameters. The VCOM is illuminated in Section 4. Section 5 presents the case study and discussion. Finally, Section 6 draws conclusions from the research.

2. Problem Formulation

2.1. FCS Model

The power supply system usually consists of the grid, the photovoltaic system (PV), and the ESS, as shown in Figure 2. Since the traditional load mainly made up of lighting and air conditioning load is too low in comparison with the EV load, it is reasonable to ignore the traditional load, in which case the FCS model can be expressed as,

$$\{\begin{array}{l}{P}_{j\kappa }^{\mathrm{G}}=P_{j\kappa }^{\mathrm{ch}}-P_{j\kappa }^{\mathrm{S}}\\ P_{j\kappa }^{\mathrm{ch}}=P_{j\kappa }^{\mathrm{EV}}-P_{j\kappa }^{\mathrm{PV}}\\ P_{j\kappa }^{\mathrm{G}}\le S_j^{\mathrm{Trans}}\end{array}$$

(1)

where $P_{j\kappa }^{\mathrm{G}}$ stands for the nodal grid load at the node j at κth time interval, $P_{j\kappa }^{\mathrm{PV}}$ is the output power of the photovoltaic (PV) source at κth time interval, $P_{j\kappa }^{\mathrm{S}}$ stands for the charge-discharge active power of the ESS at κth time interval, where negative sign stands for the discharge and positive sign stands for the charge, $P_{j\kappa }^{\mathrm{E}\mathrm{V}}$ is the EV’s charging load, and $S_j^{\mathrm{T}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}$ is the rated capacity of the service transformer at the FCS.

In many studies, the nodal load data at different time intervals are regarded as normal distribution variables, and the normal probability density function is usually used to depict the uncertainties of the load of each node [16,17,18]. We divide a day into 24 equal time intervals, every time interval is 1 h. Let ∆$P_{j\kappa }^{\mathrm{E}\mathrm{V}}$ stand for the EV’s charging load change in κth time interval, which can be calculated as follows,

$$\Delta P_{j\kappa }^{\mathrm{EV}}=P_{j\kappa }^{\mathrm{EV}}-\mathrm{mean}(P_{j\kappa 1}^{\mathrm{EV}},P_{j\kappa 2}^{\mathrm{EV}},\dots ,P_{j\kappa T}^{\mathrm{EV}})$$

(2)

Histograms of the EV’s charging load change at 4th, 5th, and 7th time intervals can be seen in Figure 3. Although EV charging load changes at these time intervals do not obey normal distribution, they concentrate within 20 kW. Since nodal load changes within 20 kW can hardly affect the nodal voltage significantly under normal circumstances, we use the constant power model to depict the EV’s charging load at the time interval when the EV’s charging load change is not a normal distribution variable. Let κ_j stand for a set of time intervals when the EV’s charging load does not obey the normal distribution at the FCS j. Then $P_{j\kappa }^{\mathrm{E}\mathrm{V}}$ can be formulated as,

$$\{\begin{array}{l}{P}_{j\kappa }^{\mathrm{EV}}~N\left({\mu }_{j\kappa },{\sigma }_{j\kappa }^2\right), \kappa \in {\mathsf{\kappa }}_j\\ P_{j\kappa }^{\mathrm{EV}}=\mathrm{mean}(P_{j\kappa 1}^{\mathrm{EV}},P_{j\kappa 2}^{\mathrm{EV}},\dots ,P_{j\kappa T}^{\mathrm{EV}}), \kappa \notin {\mathsf{\kappa }}_j\end{array}$$

(3)

where ${\mu }_{j\kappa }$ and ${\sigma }_{j\kappa }$ stand for the expectation and standard change of $P_{j\kappa }^{\mathrm{E}\mathrm{V}}$ respectively. It is reasonable to adopt the constant power model to express the output active power of the PV at the κth time interval, since the fact that the change range of the PV power output is usually much smaller than that of the EV’s charging load. Hence, when κ ∈ κ_j, $P_{j\kappa }^{\mathrm{c}\mathrm{h}}$ is also a normal distribution variable, N(${\mu }_{j\kappa }$−$P_{j\kappa }^{\mathrm{P}\mathrm{V}}$,${\sigma }_{j\kappa }^2$). Since the linear transformation of a normal distribution still follows a normal distribution, when κ ∈ κ_j, the FCS model can be written as,

$$\{\begin{array}{l}{P}_{j\kappa }^{\mathrm{G}}={\mu }_{jk}-P_{j\kappa }^{\mathrm{P}\mathrm{V}}+\Delta P_{j\kappa }^{\mathrm{ch}}-P_{j\kappa }^{\mathrm{S}}\\ \Delta P_{j\kappa }^{\mathrm{ch}}~N\left(0,{\sigma }_{j\kappa }^2\right)\\ P_{j\kappa }^{\mathrm{G}}\le S_j^{\mathrm{Trans}}\\ \kappa \in {\mathbf{\kappa }}_j\end{array}$$

(4)

where ∆$P_{j\kappa }^{\mathrm{c}\mathrm{h}}$ stands for the Gaussian error component of $P_{j\kappa }^{\mathrm{c}\mathrm{h}}$.

2.2. Impact on the Radial Distribution Network

Numerous studies show that large-scale EV fast charging easily leads to nodal voltage change in the radial distribution network, which can deteriorate the power quality. Therefore, it is meaningful to reduce the nodal voltage change caused by FCSs. In this paper, we use the evaluating function $F_{o\kappa }$(.) to quantify the adverse impact of the FCSs on voltage change at observation node o at κth time interval. $F_{o\kappa }$(.) can be expressed as follows,

$$\{\begin{array}{l}{F}_{o\kappa }\left({\delta }_{o\kappa },{\beta }_{o\kappa },{\alpha }_{o\kappa }\right)=\frac{1}{\mathsf{\Gamma }\left(\alpha \right)}{\displaystyle {\int }_0^{{\beta }_{o\kappa }{\delta }_{{}_{o\kappa }}^2}{t}^{{\alpha }_{o\kappa }-1}{e}^{-t}dt}\\ \mathsf{\Gamma }\left({\alpha }_{o\kappa }\right)={\displaystyle {\int }_0^{+\infty }{t}^{{\alpha }_{o\kappa }-1}{e}^{-t}dt, {\alpha }_{o\kappa }>0}\end{array}$$

(5)

where ${\delta }_{o\kappa }^{o\kappa }$ stands for the upper limit of voltage change at observation node o at κth time interval, ${\beta }_{o\kappa }$ and ${\alpha }_{o\kappa }$ are control parameter and shape parameter, respectively, and Γ(.) is a Gamma function. The detailed descriptions of $F_{o\kappa }$(.) and the characteristic of the control parameter ${\beta }_{o\kappa }$ are introduced in Section 3.2 and Section 3.3, respectively.

2.3. Optimization Model

We can build up a one-objective optimization model to describe the problem which is focused on in this paper,

$$\begin{array}{l}\max f\left(\mathbf{x}\right)\\ \mathrm{s}.\mathrm{t}.\\ c_j\le u_j\left(\mathbf{x}\right)\le d_j, j=1,\dots ,q\\ a_{\kappa }\le x_{\kappa }\le b_{\kappa }, \kappa =1,\dots ,m\end{array}$$

(6)

where f(x) stands for the objective function calculated by the elements of the optimization vector x, x = [$x_1$, …, $x_m$], u_j(x) stands for the jth inequality constraint. It is obvious that many intelligent algorithms can be used to solve Equation (6), such as particle swarm optimization (PSO) and the genetic algorithm (GA). Compared with other intelligent algorithms, the traditional GA is more robust and can effectively search the complex spaces, but it is easier to fall into a local optimum [19,20,21]. Therefore, an improved real-coded genetic algorithm (RCGA-rdn) is used to solve Equation (6), which can reduce the probability of falling into local optimum and improve the computational efficiency [22]. The constrained optimization model needs to be converted into an unconstrained optimization model before using the RCGA-rdn. In this paper we use the penalty function method to achieve this conversion, then Equation (6) can be rewritten as,

$$\{\begin{array}{l}\max \left\{P_1\left(\mathbf{x}\right)\cdot f_1\left(\mathbf{x}\right),\dots ,P_{pp}\left(\mathbf{x}\right)\cdot f_{pp}\left(\mathbf{x}\right)\right\},\\ P_g\left(\mathbf{x}\right)=1-\frac{1}{q}\sqrt{{\displaystyle \sum _{j=1}^q{\left(\frac{\Delta u_{jg}\left(\mathbf{x}\right)}{\Delta u_j^{\max }\left(\mathbf{x}\right)}\right)}^2}}\\ \Delta u_{jg}\left(\mathbf{x}\right)=\{\begin{array}{ll}\left|u_{jg}\left(\mathbf{x}\right)-d_j\right|& ,\mathrm{if} u_{jg}\left(\mathbf{x}\right)>d_j \\ \left|u_{jg}\left(\mathbf{x}\right)-c_j\right|& ,\mathrm{if} u_{jg}\left(\mathbf{x}\right)<c_j\\ 0& ,\mathrm{else}\end{array}\\ \Delta u_j^{\max }\left(\mathbf{x}\right)=\max \left\{\Delta u_{j1}\left(\mathbf{x}\right),\dots ,\Delta u_{jpp}\left(\mathbf{x}\right)\right\}\\ j=1,\dots q\\ g=1,\dots pp\end{array}$$

(7)

where P_g(.) stands for the penalty function of population g, and pp is the initial population number. For a detailed description of the optimization model, please refer to Section 4.

3. Analytical Assessment Model of the Nodal Voltage Change

3.1. The Radial Distribution Network Model

Figure 4 shows a typical radial distribution network, there is only one source node that can be regarded as a voltage source and we use current sources to depict the load nodes. Let U_s stand for the phase voltage at the source node s, U_d stand for the phase voltage at the load node d, and U_o stand for the phase voltage at the observation node o. Based on the circuit superposition theorem, the typical radial distribution network can be decomposed into three different type sub-circuits as shown in Figure 4. The voltage at the observation node o (U_o) can be expressed as,

$$U_o=U_{\mathrm{s}}-{\displaystyle \sum _{d\in \mathbf{D}}\frac{S_d^{*}}{U_d^{*}}{Z}_{od}}$$

(8)

where $S_d^{*}$ stands for the conjugate complex draw power at node d, $U_d^{*}$ is the conjugate phase voltage at node d, Z_od is the shared impedance between the node d and observation node o from the source node s, and D is a set of the load nodes.

Let the complex draw power at the node d change from S_d to S_d + ∆S_d and the nodal voltage change from U_d to U_d + ∆U_d, then the voltage change at the observation (∆U_o) can be written as,

$$\begin{array}{ll}\Delta U_o& =-{\displaystyle \sum _{d\in \mathbf{D}}\left(\frac{S_d^{*}+\Delta S_d^{*}}{U_d^{*}+\Delta U_d^{*}}-\frac{S_d^{*}}{U_d^{*}}\right)Z_{od}}\\ & ={\displaystyle \sum _{d\in \mathbf{D}}\frac{S_d^{*}\Delta U_d^{*}-\Delta S_d^{*}{U}_d^{*}}{U_d^{*}\left(U_d^{*}+\Delta U_d^{*}\right)}{Z}_{od}}\end{array}$$

(9)

Since ∆$U_d^{*}$ << $U_d^{*}$($U_d^{*}$ + ∆$U_d^{*}$), ∆$U_d^{*}$/$U_d^{*}$($U_d^{*}$ + ∆$U_d^{*}$)→0, Equation (9) can be simplified as,

$$\Delta U_o=-{\displaystyle \sum _{d\in \mathbf{D}}\frac{\Delta S_d^{*}}{U_d^{*}+\Delta U_d^{*}}{Z}_{od}}$$

(10)

Let Z_od = R_od+jX_od, ∆S_d = ∆$P_d^{\mathrm{G}}$+j∆$Q_d^{\mathrm{G}}$ and ∆U_d = ∆$U_d^{\mathrm{r}}$ + j∆$U_d^{\mathrm{i}}$. By listing the real and imaginary parts of ∆U_o separately, we can get,

$$\{\begin{array}{l}\Delta U_o^{\mathrm{r}}=-{\displaystyle \sum _{d\in \mathrm{D}}\left[\begin{array}{l}\frac{\left(\Delta P_d^{\mathrm{G}}{R}_{od}+\Delta Q_d^{\mathrm{G}}{X}_{od}\right)\left(U_d^{\mathrm{r}}+\Delta U_d^{\mathrm{r}}\right)}{{\left(U_d^{\mathrm{r}}+\Delta U_d^{\mathrm{r}}\right)}^2+{\left(U_d^{\mathrm{i}}+\Delta U_d^{\mathrm{i}}\right)}^2}\\ -\frac{\left(\Delta P_d^{\mathrm{G}}{X}_{od}-\Delta Q_d^{\mathrm{G}}{R}_{od}\right)\left(U_d^{\mathrm{i}}+\Delta U_d^{\mathrm{i}}\right)}{{\left(U_d^{\mathrm{r}}+\Delta U_d^{\mathrm{r}}\right)}^2+{\left(U_d^{\mathrm{i}}+\Delta U_d^{\mathrm{i}}\right)}^2}\end{array}\right]}\\ \Delta U_o^{\mathrm{i}}=-{\displaystyle \sum _{d\in \mathrm{D}}\left[\begin{array}{l}\frac{\left(\Delta P_d^{\mathrm{G}}{X}_{od}-\Delta Q_d^{\mathrm{G}}{R}_{od}\right)\left(U_d^{\mathrm{r}}+\Delta U_d^{\mathrm{r}}\right)}{{\left(U_d^{\mathrm{r}}+\Delta U_d^{\mathrm{r}}\right)}^2+{\left(U_d^{\mathrm{i}}+\Delta U_d^{\mathrm{i}}\right)}^2}\\ +\frac{\left(\Delta P_d^{\mathrm{G}}{R}_{od}+\Delta Q_d^{\mathrm{G}}{X}_{od}\right)\left(U_d^{\mathrm{i}}+\Delta U_d^{\mathrm{i}}\right)}{{\left(U_d^{\mathrm{r}}+\Delta U_d^{\mathrm{r}}\right)}^2+{\left(U_d^{\mathrm{i}}+\Delta U_d^{\mathrm{i}}\right)}^2}\end{array}\right]}\end{array}$$

(11)

where ∆$U_o^{\mathrm{r}}$ and ∆$U_o^{\mathrm{i}}$ stand for the real and imaginary parts of ∆U_o. Since the phase angle of the phase voltage is usually small, it is reasonable to ignore the ∆$U_d^{\mathrm{r}}$ and ∆$U_d^{\mathrm{i}}$. Let U_d = |U_d| ∠ θ_d, then Equation (11) can be rewritten as,

$$\{\begin{array}{l}\Delta U_o^{\mathrm{r}}\approx {\displaystyle \sum _{d\in \mathbf{D}}\left[A_{od}\Delta P_d^{\mathrm{G}}+B_{od}\Delta Q_d^{\mathrm{G}}\right]}\\ \Delta U_o^{\mathrm{i}}\approx {\displaystyle \sum _{d\in \mathbf{D}}\left[C_{od}\Delta P_d^{\mathrm{G}}+D_{od}\Delta Q_d^{\mathrm{G}}\right]}\\ A_{od}=-\frac{R_{od}\cos {\theta }_d-X_{od}\sin {\theta }_d}{\left|U_d\right|},B_{od}=-\frac{X_{od}\cos {\theta }_d+R_{od}\sin {\theta }_d}{\left|U_d\right|}\\ C_{od}=-\frac{X_{od}\cos {\theta }_d+R_{od}\sin {\theta }_d}{\left|U_d\right|},D_{od}=\frac{R_{od}\cos {\theta }_d-X_{od}\sin {\theta }_d}{\left|U_d\right|}\end{array}$$

(12)

Since U_d and Z_od are known variables, A, B, C, and D are constants. Therefore, ∆$U_o^{\mathrm{r}}$ and ∆$U_o^{\mathrm{i}}$ can be approximated as a linear combination of ∆$P_d^{\mathrm{G}}$ and ∆$Q_d^{\mathrm{G}}$.

3.2. The AMM

It is obvious that the probability of the voltage change (∆U_o) smaller than the limiting value (δ_o) can be used to measure the impact of the FCSs on the voltage change at the observation node o. The larger P(|∆U_o| < δ_o) is, the smaller the impact of the FCSs on the voltage change will be. In this section, we derive the AMM which is used to evaluate the impact of the FCSs on the voltage change at any node in the LV radial distribution network. In this paper, we select the average load at the node without the FCS as the initial load at this node. Let the complex draw power change at the node without the FCS be zero, then we can obtain the nodal voltage change of the observation node, which is caused by charging load changes of FCSs. The reactive power change (∆Q) can be ignored at the node with the FCS since the power factor of the node with the FCS is usually close to 1 so that the reactive power change is usually very small. Hence, for an n-node LV distribution network system, Equation (12) can be rewritten as,

$$\{\begin{array}{l}\Delta U_o^{\mathrm{r}}={\displaystyle \sum _{i=1}^n{A}_{oi}\Delta P_i^{\mathrm{G}}},A_{oi}=-\frac{R_{oi}\cos {\theta }_i-X_{oi}\sin {\theta }_i}{\left|U_i\right|}\\ \Delta U_o^{\mathrm{i}}={\displaystyle \sum _{i=1}^n{C}_{oi}\Delta P_i^{\mathrm{G}}},C_{oi}=-\frac{X_{oi}\cos {\theta }_i+R_{oi}\sin {\theta }_i}{\left|U_i\right|}\end{array}$$

(13)

where ∆$P_i^{\mathrm{G}}$ is 0 at the node without the FCS and is ∆$P_{j\kappa }^{\mathrm{G}}$ at the node with the FCS. Rewrite Equation (13) into a vector form as follows,

$$\left[\begin{array}{l}\Delta U_o^{\mathrm{r}}\\ \Delta U_o^{\mathrm{i}}\end{array}\right]=\left[\begin{array}{l}\begin{array}{ccc}{A}_{o1}& \cdots & A_{on}\end{array}\\ \begin{array}{ccc}{C}_{o1}& \cdots & C_{on}\end{array}\end{array}\right]\left[\begin{array}{c}\Delta P_1^{\mathrm{G}}\\ ⋮\\ \Delta P_n^{\mathrm{G}}\end{array}\right]$$

(14)

Let the ESS at the FCS have two output modes as follows:

$$\mathrm{Mode}1:P_{j\kappa }^{\mathrm{S}}=\left(1-k_{j\kappa }\right)P_{j\kappa }^{\mathrm{ch}},\kappa \in {\mathsf{\kappa }}_j^{}$$

(15)

$$\mathrm{Mode}2:P_{j\kappa }^{\mathrm{S}}=C_{j\kappa }^{\mathrm{S}},\kappa \notin {\mathsf{\kappa }}_j^{}$$

(16)

where the output active power of the ESS is a normal distribution variable when κ ∈ κ_j and is a constant ($C_{j\kappa }^{\mathrm{S}}$) when κ ∉ κ_j. By plugging Equation (15) into Equation (1), $P_{j\kappa }^{\mathrm{G}}$ is,

$$P_{j\kappa }^{\mathrm{G}}=k_{j\kappa }{P}_{j\kappa }^{\mathrm{ch}}, \kappa \in {\mathsf{\kappa }}_j^{}$$

(17)

Since $P_{j\kappa }^{\mathrm{c}\mathrm{h}}$ is a normal variable when κ ∈ κ_j, $P_{j\kappa }^{\mathrm{G}}$ is still a normal variable at this time. Let ${\mu }_{j\kappa }^{\mathrm{G}}$ stand for the expectation of $P_{j\kappa }^{\mathrm{G}}$ when κ ∈ κ_j. Select ${\mu }_{j\kappa }^{\mathrm{G}}$ as the initial nodal load for the FCS j, then ∆$P_{j\kappa }^{\mathrm{G}}$ is a normal variable, ∆$P_{j\kappa }^{\mathrm{G}}$~N(0, ${\gamma }_{j\kappa }^2$). By plugging Equation (16) into Equation (1), we can obtain the expression of $P_{j\kappa }^{\mathrm{G}}$ when κ ∉ κ_j,

$$P_{j\kappa }^{\mathrm{G}}=P_{j\kappa }^{\mathrm{EV}}-P_{j\kappa }^{\mathrm{PV}}-C_{j\kappa }^{\mathrm{S}},\kappa \notin {\mathsf{\kappa }}_j^{}$$

(18)

Since $P_{j\kappa }^{\mathrm{P}\mathrm{V}}$, $P_{j\kappa }^{\mathrm{E}\mathrm{V}}$, and $C_{j\kappa }$ are all constant when κ ∉ κ_j, ∆$P_{j\kappa }^{\mathrm{G}}$ is zero at this time. Then Equation (14) can be rewritten as,

$$\{\begin{array}{l}\left[\begin{array}{l}\Delta U_{o\kappa }^{\mathrm{r}}\\ \Delta U_{o\kappa }^{\mathrm{i}}\end{array}\right]=\left[\begin{array}{l}\begin{array}{ccc}{A}_{o1\kappa }& \cdots & A_{on\kappa }\end{array}\\ \begin{array}{ccc}{C}_{o1\kappa }& \cdots & C_{on\kappa }\end{array}\end{array}\right]\left[\begin{array}{c}\Delta P_{1\kappa }^{\mathrm{G}}\\ ⋮\\ \Delta P_{n\kappa }^{\mathrm{G}}\end{array}\right]\\ \Delta P_{i\kappa }^{\mathrm{G}}~N\left(0,{\gamma }_{j\kappa }^2\right),\kappa \in {\mathsf{\kappa }}_j^{}\&\&i\in \mathbf{E}\\ \Delta P_{i\kappa }^{\mathrm{G}}=0,\kappa \notin {\mathsf{\kappa }}_j^{}\left|\right|i\notin \mathbf{E}\end{array}$$

(19)

where E is a set of nodes with the FCS. According to Equation (19), ∆$U_{o\kappa }^{\mathrm{r}}$ and ∆$U_{o\kappa }^{\mathrm{i}}$ are linear combinations of ∆$P_{j\kappa }^{\mathrm{G}}$. Therefore, ∆$U_{o\kappa }^{\mathrm{r}}$ and ∆$U_{o\kappa }^{\mathrm{i}}$ can be regarded as normal distribution variables. We assume that the load at each node is independent of each other, then ∆$U_{o\kappa }^{\mathrm{r}}$ and ∆$U_{o\kappa }^{\mathrm{i}}$ can be expressed as,

$$\Delta U_{o\kappa }^{\mathrm{r}}~N\left(0,{\displaystyle \sum _{j=1}^n{A}_{oj\kappa }^2{\gamma }_{j\kappa }^2}\right),\mathrm{and} \Delta U_{o\kappa }^{\mathrm{i}}~N\left(0,{\displaystyle \sum _{j=1}^n{C}_{oj\kappa }^2{\gamma }_{j\kappa }^2}\right)$$

(20)

Since Equation (20) has the same form for each time interval, the subscript κ is omitted for convenience of writing in the following derivation. Let ∆U_o = (∆$U_o^{\mathrm{r}}$,∆$U_o^{\mathrm{i}}$) ^T, ∆U_o~N(0,C_o), where C_o is a covariance matrix of ∆U_o for observation node o. C_o can be written as,

$${\mathbf{C}}_o=\left[\begin{array}{cc}{\displaystyle \sum _{j=1}^n{A}_{oj}^2}{\gamma }_j^2& {\displaystyle \sum _{j=1}^n{A}_{oj}^{}{C}_{oj}^{}}{\gamma }_j^2\\ {\displaystyle \sum _{j=1}^n{A}_{oj}^{}{C}_{oj}^{}}{\gamma }_j^2& {\displaystyle \sum _{j=1}^n{C}_{oj}^{{}^2}}{\gamma }_j^2\end{array}\right]$$

(21)

C_o can be diagonalization by eigenvalue decomposition as:

$${\mathsf{\Lambda }}_o{=\mathrm{W}}_o{}^{\mathrm{T}}{\mathrm{C}}_o{\mathrm{W}}_o=\left[\begin{array}{cc}{\lambda }_{1o}& \\ & {\lambda }_{2o}\end{array}\right]$$

(22)

where W_o is an eigenmatrix of C_o, ${\lambda }_{1o}$ and ${\lambda }_{2o}$ are related eigenvalues and Λ_o is the diagonally similar matrix of C_o. Let V_o = $W_o^{\mathrm{T}}$∆U_o, where V_o = (V_1o,V_2o)^T. Since Cov(V_1o,V_2o) is equal to 0, V_1o and V_2o are independent mutual variables. $V_o^{\mathrm{T}}$V_o can be expanded as,

$$\begin{array}{ll}{\mathbf{V}}_o^{\mathrm{T}}{\mathbf{V}}_o& ={\mathsf{\Delta }\mathbf{U}}_o^{\mathrm{T}}{\mathbf{W}}_o{\mathbf{W}}_o^{\mathrm{T}}{\mathsf{\Delta }\mathbf{U}}_o\\ & ={\left(\Delta U_o^{\mathrm{r}}\right)}^2+{\left(\Delta U_o^{\mathrm{i}}\right)}^2\end{array}$$

(23)

Since ∆$U_o^2$ = (∆$U_o^{\mathrm{r}}$)² + (∆$U_o^{\mathrm{r}}$)², ∆$U_o^2$=$V_o^{\mathrm{T}}$V_o. Therefore, ∆$U_o^2$ is the sum of two independent weighed chi-square random variables. Then ∆$U_o^2$ can approximately obey Gamma distribution (Γ(α_o, β_o)) and the parameters of Γ (α_o, β_o) can be calculated as [22],

$${\alpha }_o=\frac{{\left({\lambda }_1+{\lambda }_2\right)}^2}{2\left({\lambda }_1^2+{\lambda }_2^2\right)}, \mathrm{and} {\beta }_o=\frac{{\lambda }_1+{\lambda }_2}{2\left({\lambda }_1^2+{\lambda }_2^2\right)}$$

(24)

Then, P(|∆U_o| < δ_o) can be calculated as,

$$F_o\left({\delta }_o,{\beta }_o,{\alpha }_o\right)=\mathrm{P}\left(\left|\Delta U_o\right|<{\delta }_o\right)={\displaystyle {\int }_0^{{\delta }_o^2}\frac{{\beta }_o^{{\alpha }_o}}{\mathsf{\Gamma }\left({\alpha }_o\right)}{x}^{{\alpha }_o-1}{e}^{-{\beta }_{o}x}dx}$$

(25)

The smaller the $F_o$(.) is, the greater the impact of the FCSs on the voltage change at the observation node o will be.

3.3. The Control Characteristic of β_o

When δ_o is a known constant, Fo(.) is decided by parameters α_o and β_o. We can increase Fo(.) by controlling the two parameters. In this section, we discuss the control characteristic of β_o, which can give theoretical support for the VCOM. According to Equation (24), α_o and β_o are functions of λ_1o and λ_2o. We next note that λ_1o and λ_2o can be expanded as,

$$\{\begin{array}{l}{\lambda }_{1o}=\frac{1}{2}\left\{{\displaystyle \sum _{j=1}^n{K}_{oj}}{\gamma }_j^2+\sqrt{{\left\{{\displaystyle \sum _{j=1}^n{K}_{oj}}{\gamma }_j^2\right\}}^2-4{\displaystyle \sum _{j=1}^{n-1}{Z}_{oj}}{\gamma }_{j\kappa }^2{\gamma }_{j+1}^2}\right\}\\ {\lambda }_{2o}=\frac{1}{2}\left\{{\displaystyle \sum _{j=1}^n{K}_{oj}}{\gamma }_j^2-\sqrt{{\left\{{\displaystyle \sum _{j=1}^n{K}_{oj}}{\gamma }_j^2\right\}}^2-4{\displaystyle \sum _{j=1}^{n-1}{Z}_{oj}}{\gamma }_{j\kappa }^2{\gamma }_{j+1}^2}\right\}\\ K_{oj}=A_{oj}^2+C_{oj}^2\\ Z_{oj}=A_{oj}{C}_{oj+1}-A_{oj+1}{C}_{oj}\end{array}$$

(26)

Further expand Z_oj and we can obtain the expansion of Z_oj as follows,

$$\begin{array}{l}{Z}_{oj}=A_{oj}{C}_{oj+1}-A_{oj+1}{C}_{oj}\\ =\frac{\left(R_{oj}{R}_{oj+1}-X_{oj}{X}_{oj+1}\right)\sin {\theta }_{jj+1}+\left(R_{oj}{X}_{oj+1}-R_{oj+1}{X}_{oj}\right)\cos {\theta }_{jj+1}}{\left|V_j\right|\left|V_{j+1}\right|}\end{array}$$

(27)

where ${\theta }_{jj+1}$ = ${\theta }_{j+1}$−θ. Since the phase difference of the voltage between the neighboring nodes in the LV radial distribution network is usually very small, it is reasonable to assume that ${\theta }_j$ ≈ ${\theta }_{j+1}$, then Equation (22) can be simplified as,

$$Z_{oj}=\frac{R_{oj}{X}_{oj+1}-R_{oj+1}{X}_{oj}}{\left|V_j\right|\left|V_{j+1}\right|}$$

(28)

The impedance angle difference between shared impedances between adjacent nodes is usually small, therefore,

$$\begin{array}{l}\tan {\phi }_{oj}-\tan {\phi }_{oj+1}=\frac{X_{oj}}{R_{oj}}-\frac{X_{oj+1}}{R_{oj+1}}\to 0\\ \Rightarrow Z_j\to 0\\ \Rightarrow {\left\{{\displaystyle \sum _{j=1}^n{K}_j}{\gamma }_j^2\right\}}^2\gg 4{\displaystyle \sum _{j=1}^{n-1}{Z}_j}{\gamma }_{j\kappa }^2{\gamma }_{j+1}^2\\ \Rightarrow {\left\{{\displaystyle \sum _{j=1}^n{K}_j}{\gamma }_j^2\right\}}^2-4{\displaystyle \sum _{j=1}^{n-1}{Z}_j}{\gamma }_{j\kappa }^2{\gamma }_{j+1}^2\to {\left\{{\displaystyle \sum _{j=1}^n{K}_j}{\gamma }_j^2\right\}}^2\end{array}$$

(29)

where φ_oj is the impedance angle of the shared impedance $Z_{oj}$. Therefore, Equation (26) can be simplified as,

$$\{\begin{array}{l}{\lambda }_1\approx {\displaystyle \sum _{j=1}^n\left[A_{oj}^2+C_{oj}^2\right]}{\gamma }_j^2\\ {\lambda }_2\approx 0\end{array}$$

(30)

By plugging Equation (30) into Equation (24), we can obtain,

$${\alpha }_o\approx 0.5, \mathrm{and} {\beta }_o=\frac{1}{2{\displaystyle \sum _{j=1}^n\left[A_{oj}^2+C_{oj}^2\right]}{\gamma }_j^2}$$

(31)

From Equation (31), α_o is approximately equal to 0.5, which is not affected by the change of the nodal load change. Let α_o = 0.5 and t = β_ox, then $F_o$(.) can be rewritten as,

$$F_o\left({\delta }_o,{\beta }_o\right)={\displaystyle {\int }_0^{{\beta }_o{\delta }_o^2}\frac{t^{-0.5}{e}^{-t}}{\mathsf{\Gamma }\left(0.5\right)}dt}$$

(32)

Since the integrand function of Equation (32) is not less than 0 and β_o > 0, $F_o$(.) increases as the increase of β_o until $F_o$(.) reaches 1. In other words, the smaller the value of β_o is, the smaller the nodal voltage change at the observation will be.

4. The VCOM

We restore the subscript κ in the rest of this paper. The objective of the day-ahead scheduling strategy for the ESS is to reduce the adverse impact of the FCSs on nodal voltage change at observation nodes. According to the analysis in Section 3, we can use the AMM to evaluate this impact. Therefore, the optimization objective function can be written as,

$$F_o\left({\delta }_o,{\beta }_o\right)={\displaystyle {\int }_0^{{\beta }_o{\delta }_o^2}\frac{t^{-0.5}{e}^{-t}}{\mathsf{\Gamma }\left(0.5\right)}dt}$$

(33)

where O is a set of observation nodes. The inequality constraint of $P_{j\kappa }^{\mathrm{S}}$ can be written as,

$$\mathrm{P}\left(P_{j\kappa }^{\mathrm{Sd}}\le P_{j\kappa }^{\mathrm{S}}\le P_j^{\mathrm{S}+}\right)\ge {\xi }_j^{\mathrm{Ps}},\kappa \in {\mathsf{\kappa }}_j$$

(34)

$$P_j^{\mathrm{S}-}\le P_{j\kappa }^{\mathrm{S}}\le P_j^{\mathrm{S}+},\kappa \in {\mathsf{\kappa }}_j$$

(35)

where ${\xi }_j^{Ps}$ is a preset constant for restraining $P_{j\kappa }^{\mathrm{S}}$, $P_j^{\mathrm{S}+}$, and $P_{j\kappa }^{\mathrm{S}\mathrm{d}}$ stand for the rated upper and dynamic lower limits of the charge–discharge active power for the ESS at the FCS j, respectively. $P_{j\kappa }^{\mathrm{S}-}$ can be expressed as,

$$P_{j\kappa }^{\mathrm{Sd}}=\{\begin{array}{l}{P}_{j\kappa }^{\mathrm{ch}}-S^{\mathrm{Trans}},P_{j\kappa }^{\mathrm{ch}}-S^{\mathrm{Trans}}>P_j^{\mathrm{S}-}\\ P_j^{\mathrm{S}-},P_{j\kappa }^{\mathrm{ch}}-S^{\mathrm{Trans}}\le P_j^{\mathrm{S}-}\end{array}$$

(36)

where $P_j^{\mathrm{S}-}$ is the rated lower limit of the charge-discharge active power for the ESS at the FCS j. Let k_jκ = ${\gamma }_{j\kappa }$/δ_jκ when κ ∈ κ_j, we can rewrite Equation (15) by plugging k_jκ = ${\gamma }_{j\kappa }$/δ_jκ in Equation (15),

$$P_{j\kappa }^{\mathrm{S}}=\left(1-\frac{{\gamma }_{j\kappa }^{}}{{\sigma }_{j\kappa }^{}}\right)P_{j\kappa }^{\mathrm{ch}},\kappa \in \mathsf{\kappa }$$

(37)

Next, according to the Bayesian probability formula, the analytical expression for the constraint (34) can be expanded as,

$$\{\begin{array}{l}\mathrm{P}\left(P_{j\kappa }^{\mathrm{Sd}}\le P_{j\kappa }^{\mathrm{S}}\le P_j^{\mathrm{S}+}\right)=\{\begin{array}{l}\left[\begin{array}{l}\Phi \left(a_{j\kappa }\right)-\left(1-\Phi \left(b_{j\kappa }\right)\right)\left(1-\Phi \left(c_{j\kappa }\right)\right)\\ -\Phi \left(c_{j\kappa }\right)\Phi \left(d_{j\kappa }\right)\ge {\xi }_j^{\mathrm{Ps}}\end{array}\right],\kappa \in {\mathsf{\kappa }}_j\\ \Phi \left(c_{j\kappa }\right)+\Phi \left(e_{j\kappa }\right)-\Phi \left(c_{j\kappa }\right)\Phi \left(e_{j\kappa }\right)\ge {\xi }_j^{\mathrm{Ps}},\kappa \notin {\mathsf{\kappa }}_j\end{array}\\ a_{j\kappa }=\frac{{\sigma }_{j\kappa }^{}{P}_j^{\mathrm{S}+}}{{\sigma }_{j\kappa }^{}\left({\sigma }_{j\kappa }^{}-{\gamma }_{j\kappa }^{}\right)}-\frac{{\mu }_{j\kappa }^{}-P_{j\kappa }^{\mathrm{PV}}}{{\sigma }_{j\kappa }^{}},b_{j\kappa }=\frac{S_j^{\mathrm{Trans}}}{k_{\kappa }{\sigma }_{j\kappa }^{}}-\frac{{\mu }_{j\kappa }^{}-P_{j\kappa }^{\mathrm{PV}}}{{\sigma }_{j\kappa }^{}}\\ c_{j\kappa }=\frac{P_{j\kappa }^{\mathrm{PV}}+P_j^{\mathrm{S}-}+S^{\mathrm{Trans}}-{\mu }_{j\kappa }^{}}{{\sigma }_{j\kappa }^{}},d_{j\kappa }=\frac{{\sigma }_{j\kappa }^{}{P}_j^{\mathrm{S}-}}{{\sigma }_{j\kappa }^{}\left({\sigma }_{j\kappa }^{}-{\gamma }_{j\kappa }^{}\right)}-\frac{{\mu }_{j\kappa }^{}-P_{j\kappa }^{\mathrm{PV}}}{{\sigma }_{j\kappa }^{}}\\ e_{j\kappa }=\frac{P_{j\kappa }^{\mathrm{PV}}+P_j^{\mathrm{S}-}+S^{\mathrm{Trans}}-{\mu }_{j\kappa }^{}}{{\sigma }_{j\kappa }^{}}\end{array}$$

(38)

where Φ(.) is the standard normal probability distribution function. Let ${\mathrm{SOC}}_{j0}$ stand for the initial state of charge (SOC) for the ESS at the FCS j at the beginning of the day, then ${\mathrm{SOC}}_{j\kappa }$ can be written as,

$${\mathrm{SOC}}_{j\kappa }={\mathrm{SOC}}_{j0}-\frac{1}{E_j^{\mathrm{S}}}{\displaystyle \sum _{i=1}^{\kappa }{P}_{ji}^{\mathrm{S}}}$$

(39)

where $E_j^{\mathrm{S}}$ stands for the rated capacity of the ESS at the FCS j. The constraint of the SOC can be written as,

$$\mathrm{P}\left({\mathrm{SOC}}_j^{-}\le {\mathrm{SOC}}_{j\kappa }\le {\mathrm{SOC}}_j^{+}\right)\ge {\xi }_j^{\mathrm{SOC}}$$

(40)

where ${\xi }_j^{\mathrm{S}\mathrm{O}\mathrm{C}}$ is a preset constant for restraining ${\mathrm{S}\mathrm{O}\mathrm{C}}_{j\kappa }$. ${\mathrm{S}\mathrm{O}\mathrm{C}}_j^{+}$ and ${\mathrm{S}\mathrm{O}\mathrm{C}}_j^{-}$ stand for the rated upper and lower limits of the SOC for the ESS at the FCS j. The analytical expression of the constraint (40) is,

$$\{\begin{array}{l}\mathrm{P}\left(\mathrm{SO}{\mathrm{C}}_j^{-}\le \mathrm{SO}{\mathrm{C}}_{j\kappa }\le \mathrm{SO}{\mathrm{C}}_j^{+}\right)=\Phi \left(w_{j\kappa }\right)-\Phi \left(r_{j\kappa }\right)\ge {\xi }_j^{\mathrm{SOC}}\\ w_{j\kappa }=\frac{E_j^{\mathrm{S}}\left(\mathrm{SO}{\mathrm{C}}_{j0}-\mathrm{SO}{\mathrm{C}}_j^{-}\right)}{{\sigma }_Q^{}}-\frac{{\displaystyle \sum _{i\notin {\mathsf{\kappa }}_j\&\&i\le \kappa }{P}_{ji}^{\mathrm{S}}}+{\mu }_Q}{{\sigma }_Q^{}}\\ r_{j\kappa }=\frac{E_j^{\mathrm{S}}\left(\mathrm{SO}{\mathrm{C}}_{j0}-\mathrm{SO}{\mathrm{C}}_j^{+}\right)}{{\sigma }_Q^{}}-\frac{{\displaystyle \sum _{i\notin \left|{\mathsf{\kappa }}_j\right|\&\&i\le \kappa }{P}_{ji}^{\mathrm{S}}}+{\mu }_Q}{{\sigma }_Q^{}}\\ {\mu }_Q={\displaystyle \sum _{i=1}^{\left|{\mathsf{\kappa }}_j\cap {\mathbf{T}}_j^{\kappa }\right|}\left(1-\frac{{\gamma }_{ji}^{}}{{\sigma }_{ji}^{}}\right)\left({\mu }_{j\kappa }^{}-P_{j\kappa }^{\mathrm{PV}}\right)}\\ {\sigma }_Q^2={\displaystyle \sum _{i=1}^{\left|{\mathsf{\kappa }}_j\cap {\mathbf{T}}_j^{\kappa }\right|}{\left(1-\frac{{\gamma }_{ji}^{}}{{\sigma }_{ji}^{}}\right)}^2{\sigma }_{ji}^2}+2\left({\displaystyle \sum _{i=1}^{\left|{\mathsf{\kappa }}_j\cap {\mathbf{T}}_j^{\kappa }\right|-1}{\displaystyle \sum _{l=i+1}^{\left|{\mathsf{\kappa }}_j\cap {\mathbf{T}}_j^{\kappa }\right|}\left(1-\frac{{\gamma }_{ji}^{}}{{\sigma }_{ji}^{}}\right)\left(1-\frac{{\gamma }_{jl}^{}}{{\sigma }_{jl}^{}}\right)\mathrm{Cov}\left(P_{ji}^{\mathrm{ch}},P_{jl}^{\mathrm{ch}}\right)}}\right)\end{array}$$

(41)

where $T_j^{\kappa }$ is a set of time intervals, $T_j^{\kappa }$ = [1, 2, …, κ]. Then the basic constraints which are used to ensure the normal operation of the ESS at the FCS j can be listed as follows:

$$\{\begin{array}{l} \left(36\right), \left(39\right)\\ P_j^{\mathrm{S}-}\le P_{j\kappa }^{\mathrm{S}}\le P_j^{\mathrm{S}+},\kappa \in {\mathsf{\kappa }}_j\end{array}$$

(42)

Using the analysis according to Section 3.3, we can control the voltage change by adjusting ${\beta }_o$, which means β_oκ can be designed as a constraint which is used to restrain the voltage change at the observation node in the optimization framework as follows:

$${\beta }_{o\kappa }\ge {\beta }_{o\kappa }^{\min }$$

(43)

where ${\beta }_{o\kappa }^{\min }$ stands for the lower limit for ${\beta }_{o\kappa }$. According to Equation (32), F_oκ(.) increases with the increase of ${\beta }_{o\kappa }$ until it reaches 1 when δ_o is a known constant, we can obtain ${\beta }_{o\kappa }^{\min }$ by the Algorithm 1 as follows:

Algorithm 1: Calculate${\mathit{\beta }}_{\mathit{o}\mathit{\kappa }}^{\mathbf{min}}$

1: INPUT: This algorithm knows the acceptable probability of the nodal voltage change at the observation node ($F_{o\kappa }^{\min }$) at time interval κ, the initial nodal load at each node at time interval κ, the adjustment step size of ${\beta }_{o\kappa }^{\min }$ (∆), basic information of the distribution network, which can be used to calculate power flow by using traditional Newton power flow method.

2: OUTPUT: ${\beta }_{o\kappa }^{\min }$ = [${\beta }_{1\kappa }^{\min }$, ${\beta }_{2\kappa }^{\min }$,…, ${\beta }_{n\kappa }^{\min }$]

3: PROCEDURE:

4: Obtain the nodal voltage at each node by using traditional Newton power flow method.

4: for o = 1 to n do

5: Obtain ${\beta }_{o\kappa }$ according to Equation (31)

6: Obtain Foκ(.) according to Equation (32)

7: for I = 1 to N do

8: if $F_{o\kappa }$(.) ≥ $F_{o\kappa }^{\min }$

9: ${\beta }_{o\kappa }^{\min }$ ←${\beta }_{o\kappa }$—∆

10: Update Foκ(.) by plugging ${\beta }_{o\kappa }^{\min }$ into Equation (32)

11: else

12 ${\beta }_{o\kappa }^{\min }$←${\beta }_{o\kappa }$

13: end if

14: end for

15: end for

Then the VCOM can be written as,

$$\{\begin{array}{l}\max {\displaystyle \sum _{\kappa =1}^m{\displaystyle \sum _{o\in \mathrm{O}}^{}{F}_{ok}\left({\delta }_{o\kappa },{\beta }_{o\kappa }\right)}}\\ \mathrm{s}.\mathrm{t}.\\ \left(38\right), \left(41\right)~\left(43\right)\\ 0\le {\gamma }_{j\kappa }^{}\le \sqrt{2}{\sigma }_{j\kappa }^{},\kappa \in {\mathsf{\kappa }}_j\\ P_j^{\mathrm{S}-}\le P_{j\kappa }^{\mathrm{S}}\le P_j^{\mathrm{S}+},\kappa \in {\mathsf{\kappa }}_j\\ j=1,\dots ,n\\ \kappa =1,\dots ,m\end{array}$$

(44)

By converting Equation (44) to the unconstrained optimization model as shown in Equation (7), RCGA-rdn can be used to solve the VCOM. The RCGA-rdn stops iterating when the number of the genetic offspring reaches the specified value and we select the best elite solution as the optimal solution.

The ratio plan can be converted into the active power plan of the ESS according to Equation (45):

$$\{\begin{array}{l}{P}_{j\kappa }^{\mathrm{S}}=P_j^{\mathrm{S}-}+v_j\left(P_j^{\mathrm{S}+}-P_j^{\mathrm{S}-}\right),\kappa \notin {\mathsf{\kappa }}_j\\ {\gamma }_{j\kappa }^{}=v_j\sqrt{2}{\sigma }_{j\kappa }^{},\kappa \in {\mathsf{\kappa }}_j\end{array}$$

(45)

where v_j is a random number in the range [0, 1].

5. Case Study and Discussion

In this section, we verify the effectiveness of the AMM and the proposed day-ahead ESS strategy. The computer used for the test adopts an Intel^® Core™ i7-6700 Processor and 8 GB of memory, in which the CPU works at the nominal frequency of 3.4 GHz. The base voltage of the system is 12.66 kV, and the base capacity is 100 MW. The reference voltage is set to 0.95.

Two FCSs are connected to nodes 10 and 25 respectively, as shown in Figure 5. The basic configurations of FCSs A and B are the same as shown in

Table 1. The basic configuration.

Parameter	Description	Unit	Value
STrans	Rated capacity of the service transformer.	kVA	1000
SPV	Rated capacity of the PV.	kWp	50
PS+	Rated upper boundary of the charge-discharge active power of the ESS.	kW	250
PS−	Rated lower boundary of the charge-discharge active power of the ESS.	kW	−250
SOC0	The initial SOC for the ESS at the beginning of the day.	%	50
SOC+	Rated upper boundary of the SOC of the ESS.	%	30
SOC−	Rated lower boundary of the SOC of the ESS.	%	80
ES	Rated capacity of the ESS.	kWh	1000

. The probability distributions of EV load changes of the FCSs A and B at different time intervals in a day are listed in

Table 2. The probability distribution of EV load change at the FCS A in different time intervals in a day.

κ	The EV Load at κth Time Interval/kW
1~4	N(789,1922)	N(532,2182)	N(167,1212)	E(40)
5~8	E(40)	E(40)	N(223,402)	N(108,742)
9~12	N(107,952)	N(190,1322)	N(294,1592)	N(406,1532)
13~16	N(547,1442)	N(572,1692)	N(575,1852)	N(352,1982)
17~20	N(207,1632)	N(306,1452)	N(312,1632)	N(309,1802)
21~24	N(328,1802)	N(350,1432)	N(494,1542)	N(690,1802)

and

Table 3. The probability distribution of EV load change at the FCS B in different time intervals in a day.

κ	The EV Load at κth Time Interval/kW
1~4	N(789,1922)	N(532,2182)	N(167,1212)	E(40)
5~8	E(40)	E(40)	E(40)	N(108,742)
9~12	N(107,952)	N(190,1322)	N(294,1592)	N(406,1532)
13~16	N(547,1442)	N(572,1692)	N(575,1852)	N(352,1982)
17~20	N(207,1632)	N(306,1452)	N(312,1632)	N(309,1802)
21~24	N(328,1802)	N(350,1432)	E(40)	N(690,1802)

respectively. Here, sets of charging power test data come from sets of 60-day historical charging power data from two different FCSs. Typical daily EV charging curves of FCSs A and B are shown in Figure 6.

Each FCS is equipped with a PV source and the typical daily PV curve is shown in Figure 7.

5.1. The Performance of the AMM

The cumulative probability distribution (CPD) of the actual voltage change is obtained by statistical analysis based on a Monte Carlo probability power flow method (MCPFL) [23], which is used as a comparison to the AMM. The sampling times of MCPFL reach 10,000. Figure 8 shows the CPD curves of the voltage change at node 23 at the time interval (17:00, 17:59). It can be seen that the CPD curve obtained by the AMM is very close to the actual CPD curve, which can verify the effectiveness of the AMM. A similar result can also be seen in Figure 9 where node 12 is selected as the observation node. Finally, the correlation coefficient is used to measure the closeness between the actual CPD curve and the CPD curve obtained by the AMM. The corresponding correlation coefficients are listed in

Table 4. The correlation coefficients.

Observation Node	The Correlation Coefficient
2~9	0.99	0.99	0.99	0.99	0.99	0.99	0.99	0.99
10~17	0.99	0.99	0.99	0.99	0.99	0.99	0.99	0.99
18~25	0.99	0.99	0.99	0.99	0.99	0.99	0.99	0.99
26~33	0.99	0.99	0.99	0.99	0.99	0.99	0.99	0.99

. The correlation coefficients all reach 0.99, which also verifies the effect of the AMM.

Table 5. Computational overheads.

Methods	Average Time, s	CPU, %
MCPFL	1.527	2
AMM	20.313	11

lists the computational overheads of assessment on voltage change at an observation node. It can be seen that the calculation time of the AMM is very short and its CPU occupancy rate is also small.

5.2. The Effect of the Proposed Day-Ahead ESS Strategy

Next, two different cases are studied to verify the proposed day-ahead ESS strategy.

(i)

Case 1: the traditional PV-ESS complementarity strategy is used for the ESS.

(ii)

Case 2: the strategy proposed in this paper is used for the ESS.

For case 1, a compensation control strategy for the ESS mentioned in [12] is used. For case 2, the ESS executes the proposed day-ahead ESS strategy. The ESS automatically enters the charging state both in cases 1 and 2 when the SOC of the ESS is lower than 30%. We select terminal node 18 as an observation node. Besides, the voltage change constraint of the observation node is set to 0.01 and the acceptable probability of the nodal voltage change of the observation node at each time interval is set to 0.95.

Figure 10 shows P(|∆U_o| < δ_o) at each node under different cases in a day. It is obvious that the day-ahead ESS strategy can increase P(|∆U_o| < δ_o) in comparison with the traditional PV-ESS complementarity strategy. Voltages under case 2 are closer to the reference voltage than other cases as shown in Figure 11a, which is more intuitively displayed in Figure 11b. Simulation results above demonstrate the effect of the proposed day-ahead ESS strategy. It can be seen from Figure 12a that the ESS charges during periods of low EV charging load and provides active power compensation by discharge during periods of large EV charging load to reduce the voltage change between actual and reference voltages. In addition, we can find that due to the constraint of the SOC of the ESS, the ESS cannot always compensate for the active power. For example, from 6.00 p.m. to 9.00 p.m., since the SOC of the ESS is close to its lower boundary, the ESS can hardly compensate the active power, which results in the failure of the day-ahead ESS strategy. A similar result can also be obtained according to Figure 13. Therefore, to give full play to the effect of the day-ahead ESS strategy, it is inseparable from reasonable energy storage system planning.

6. Conclusions

With the increasing penetration of EVs, the adverse impact of EV charging load on the nodal voltage in the distribution network can probably be aggravated. To alleviate this impact, this paper proposes a centralized day-ahead ESS strategy with full consideration of the randomness of EV charging load. This strategy can be deployed in the RODS and its output can be calculated based on the VCOM. Finally, simulation results verify the effectiveness of the proposed day-ahead ESS strategy. The main contributions of this article are as follows:

• A voltage change optimization model (VCOM) considering the randomness of the EV load is constructed to alleviate the voltage change problem caused by EV fast charging and it can be easily solved by traditional intelligent algorithms, such as GA.
• An analytical assessment model (AAM) of the nodal voltage change with shorter computational time and higher reliability is proposed.