Optimal Capacity Allocation of Energy Storage in Distribution Networks Considering Active/Reactive Coordination

Abstract

Energy storage system (ESS) has been advocated as one of the key elements for the future energy system by the fast power regulation and energy transfer capabilities. In particular, for distribution networks with high penetration of renewables, ESS plays an important role in bridging the gap between the supply and demand, maximizing the benefits of renewables and providing various types of ancillary services to cope the intermittences and fluctuations, consequently improving the resilience, reliability and flexibility. To solve the voltage fluctuations caused by the high permeability of renewables in distribution networks, an optimal capacity allocation strategy of ESS is proposed in this paper. Taking the life cycle cost, arbitrage income and the benefit of reducing network losses into consideration, a bilevel optimization model of ESS capacity allocation is established, the coordination between active/reactive power of associate power conversion system is considered, and the large scale nonlinear programming problem is solved using genetic algorithm, simulated annealing and mixed integer second-order cone programming method. The feasibility and effectiveness of the proposed algorithm have been verified.

1. Introduction

The inevitable energy revolution inflicted by the depletion of fossil fuels and global warming has attracted significant attention worldwide. In order to maintain the sustainable development of power system, the penetration of renewable distributed generation (DG) in the power system is increasing continuously. The uncertainty, variability and non-dispatchability nature of the renewables have challenged the secure and reliable operation of power systems, in particular on distribution networks (DNs) [1,2,3]. In order to mitigate the negative impacts rendered by the integration of DG without unnecessary DG power curtailment and load shedding, energy storage systems (ESSs) have been identified [3]. With proper planning and management, the ESSs also can offer diverse applications in DNs, including (i) peak shifting and load leveling [4,5], (ii) renewable power output smoothing and time shifting [6], (iii) frequency modulation [7], (iv) dynamic reactive power compensation and voltage regulation [8], (v) grid expansion and upgrade deferral [9], (vi) network congestion management [10], (vii) network reconfiguration/islanding [11,12], (viii) black start and reserve [13], etc. Among various means of energy storage technologies, battery has been selected as the energy storage unit in this study due to the technology maturity and applicability. At present, the main obstacle to the large-scale deployment of batteries is the relatively low cycle time and high cost. Therefore, the optimal planning of ESS capacity has attracted considerable attention worldwide [14,15]. With the emerging technology development, cost reduction, life extension and energy conversion efficiency improvement, the economy of battery technology will gradually become prominent [16].

The capacity allocation and operation optimization of ESS are interdependent. In recent years, research has been conducted and some valuable outcomes have been achieved. Taking the prediction error effect on ESS capacity determination into consideration, setting the total investment income model with the maximum economic benefit as the objective function, an optimal allocation scheme was developed using a dual complementary topological structure [17]. Reference [18] optimized the ESS capacity configuration from the perspective of voltage quality improvement, line losses reduction, and DG permeability promotion. In [19,20,21], capacity optimization models considering ESS cost minimization and life extension were established. In particular, the ESS capacity obtained by the optimization methods proposed in [21] was designed to meet the requirements of smoothing wind power output while maximizing the economic benefits gained from selling the dispatched power to the grid against the cost of the storage. An optimal ESS configuration model considering the benefits of ESS in energy conservation, network loss reduction and environmental protection was proposed in [22] to achieve the maximum revenue during the lifetime of the ESS. A sensitivity based algorithm was proposed for optimal sitting and sizing of the ESS to reduce the power losses on DNs. However, the daily power fluctuations of DG were not taking into account, therefore the characteristics of DG intermittency cannot be fully reflected [23]. The research mentioned above mainly concentrated on the ESS allocation methods from the perspective of the ESS cost minimization or revenue maximization, the impact of the optimal operation of the ESS, in particular the reactive power of the ESS has not been fully studied.

The power conversion system (PCS) of ESS can generate or absorb reactive power, which can be considered as an effective approach for reactive power compensation of DNs. In order to exploit the ESS’s potential capability of voltage regulations in DNs, an active/reactive power coordinated optimization method was proposed in [8] utilizing the particle swarm algorithm to realize power losses reduction and voltage fluctuation minimization in DNs. In [24], a multiobjective optimization method for battery sitting and sizing on low voltage DNs was proposed from the perspective of voltage regulation, peak power reduction and annual cost minimization. The method was taking the active and reactive power from batteries and inverters into consideration on a near-future scenario for a real residential feeder. Optimal operation methods with coordination between active and reactive power of DG and ESSs were also developed to improve the economy of the overall DNs [25,26,27]. It is worth noting that the reactive power regulations play important roles in optimal operation strategy development considering the possible ancillary services in DNs, therefore affecting the network configuration, operation and the life cycle cost (LCC), further studies are required to investigate the cooperation between the active and reactive power from the overall life time of the system.

The operation of ESS has been considered as location and time specific, and the dimension for solving the optimization problem increases rapidly with the duration of time and the number of installed locations. Therefore, the operation optimization of ESS becomes a large-scale nonlinear programming problem. At present, the general reduced gradient method [28], dynamic programming recursive method [29], genetic algorithm [30] and Pareto optimal algorithm [31] are used to solve the problem. However, the time and difficulty of calculation are closely related to the size and spatial dimension of the DNs. The convex relaxation of power flow equations can usually reduce the computational burden [32]. In [33,34], the relaxation of the power flow equation to a second-order cone has been explained theoretically and mathematically, and the accuracy of the relaxation has been proved. In [32,35] the authors explored the process of the sitting and sizing of storage units through relaxed power flow equations without considering the influence of ESS’s reactive power, and the charging/discharging model is relatively simple, which cannot fully reflect the efficiency of the PCS.

Consequently, an interactive bilevel planning strategy that determines the optimal operation and the capacity of the ESS in DNs, utilizing the genetic algorithm (GA) and simulated annealing (SA) combined with mixed integer second-order cone programming (MISOCP) is presented in this paper. With active/reactive power support from the ESSs in four quadrants, the proposed ESS planning method can both reduce the LCC and increase arbitrage revenue while providing ancillary services of voltage support in DNs.

The main contributions and salient features of this paper can be highlighted as:

(1) A bilevel ESS optimal capacity allocation model for DNs is proposed to suppress voltage fluctuations and line losses caused by photovoltaic (PV) integration, and to accommodate the increasing penetration level of renewables in the near future. From the perspective of life cycle, the proposed method considers overall cost, arbitrage income and the benefit of reducing line losses simultaneously to quantify the system economics.

(2) The active/reactive power coordination strategy through PCS can manage the voltage fluctuations without additional investment whilst exploiting the arbitrage income, which is relatively rare in recent literature.

(3) The iterative solving method utilizing GA and SA combined with MISOCP is proposed to ensure the efficiency and accuracy of the entire complicated calculation process.

The remainder of this paper is organized as follows: the overall system model including charging/discharging model, LCC model, economic benefit model and capacity fading model of the ESS, and the Distflow power flow calculation are introduced in Section 2. Section 3 presents the bilevel ESS optimization model with the active/reactive coordinated operation strategy of ESS. Case studies are carried out in Section 4, and the relevant concluding remarks are delivered in Section 5.

2. System Modeling

2.1. ESS Modeling

Due to the high energy density, long life-cycle, low self-discharging rate and high comprehensive efficiency, the Li-ion battery has been chosen as the energy storage device in this paper [36]. In order to simulate the charging and discharging characteristics of ESS accurately and quantify the cost and benefit of ESS precisely, charging/discharging model, LCC model, economic benefit model and capacity fading model of the ESS are established as follows.

2.1.1. The Charging/Discharging Model of ESS

A typical ESS mainly consists of a PCS unit, a storage unit and the connecting facility of the plant. The PCS unit is a voltage source bidirectional inverter designed to operate as either an inverter while discharging or as a rectifier while charging [26]. Due to the mismatch of the voltage levels between the battery and the grid, the PCS needs to be connected to the DNs and the battery side for voltage matching, voltage isolation and charging/discharging power control [37]. The PCS permits the ESS to generate both active and reactive power in all four quadrants as illustrated by the capability curve in Figure 1.

In Figure 1, the unit circle represents the capacity of PCS as S_ESS, and the operating point of the ESS can be any point on the unit circle. Assuming the power output of the ESS to DNs is positive, therefore the maximum discharging active power of the ESS, ${\mathrm{P}}_{\mathrm{ESS}.\mathrm{dis}}^{\max }$, has a positive value and the maximum charging active power, ${\mathrm{P}}_{\mathrm{ESS}.\mathrm{ch}}^{\max }$, is negative. ${\mathrm{Q}}_{\mathrm{ESS}}^{\max }$ is the maximum reactive power of the ESS and Q_ESS(t) is the actual operating reactive power at the time t. P_ESS.ch(t) is the actual charging active power of ESS at time t. When ESS is running in the first and fourth quadrant, P_ESS.ch(t) will be replaced by P_ESS.dis(t), which represents the actual discharging active power of ESS at time t. When the power factor of PCS is 1, it is equal to S_ESS. ESS can only be in the state of discharging or charging at the same time. So, 0-1 variables, μ_dis(t) and μ_ch(t), representing the discharging/charging sign of the ESS at time t, were introduced respectively. When the ESS is charging, μ_ch(t) is 1, μ_dis(t) is 0 and vice versa. The operating power of the ESS is required to satisfy the following constraints.

$$\sqrt{{(P_{ESS.dis}(t)+P_{ESS.ch}(t))}^2+Q_{ESS}{\left(t\right)}^2}\le S_{ESS}$$

(1)

$$\left\{\begin{array}{l}0\le P_{ESS.dis}\left(t\right)\le S_{ESS}{\mu }_{dis}\left(t\right)\hfill \\ -S_{ESS}{\mu }_{\mathrm{ch}}\left(t\right)\le P_{ESS.ch}\left(t\right)\le 0\hfill \\ {\mu }_{dis}\left(t\right)+{\mu }_{\mathrm{ch}}(t)=1\hfill \\ -S_{ESS}\le Q_{ESS}\left(t\right)\le S_{ESS}\hfill \end{array}\right.$$

(2)

In addition, the residual capacity of the battery is represented by the state-of-charge (SOC), and its value can be obtained by comparing the residual capacity value of the battery with the battery capacity value. SOC was accumulated and calculated strictly according to the magnitude of the charging and discharging power in time sequence, and its value has absolute continuity in the time sequence [38], as shown in Equation (3). The SOC value of each time t should be within the upper and lower limits of SOC, as shown in Equation (4).

$$SOC(t)=SOC(t-\mathit{1})·(1-\delta )-\frac{P_{ESS.ch}\left(t\right){\eta }_C{\eta }_{PCS}\Delta t}{E_{\mathrm{rate}}}-\frac{P_{ESS.dis}\left(t\right)\Delta t}{E_{\mathrm{rate}}{\eta }_D{\eta }_{PCS}}$$

(3)

$${\mathrm{SOC}}_{\min }\le SOC\left(t\right)\le {\mathrm{SOC}}_{\max }$$

(4)

where η_C and η_D are the charging and discharging efficiency of the battery and η_PCS represents the associate PCS efficiency. SOC_max and SOC_min are the upper and lower limits of SOC and SOC(t) is the SOC at time t. E_rate is the rated capacity of the battery and δ is the self-discharging rate of the battery.

In real operation, the ESS usually works at a fixed cycle T. In order to maintain the continuous operation of the ESS, the final SOC within each cycle T should be consistent with the initial SOC, as shown in Equation (5).

$$SOC\left(0\right)=SOC\left(T\right)$$

(5)

where SOC(0) and SOC(T) are the SOC of ESS at the initial and ending time of each operation cycle.

2.1.2. The LCC Model of ESS

According to the international standard IEC 60300-3-3, the LCC of a certain product should contain six parts: concept and definition, design and development, manufacturing, installation, operation and maintenance and disposal [39]. All of the costs considered in this study are illustrated as follows.

$$\mathit{Cost}=C_1+C_2+C_3+C_4+C_5$$

(6)

where C₁ represents the investment cost, C₂ indicates the replacement cost, C₃ is the operation and maintenance cost, C₄ is the disposal cost and C₅ is the recovery value, which is negative to represent revenue.

The hardware of the ESS are composed of three main components: the storage unit, the PCS and the connecting facilities of the plant. So, the investment cost can be presented as (7):

$$C_1= \left({\mathrm{c}}_{\mathrm{E}}{E}_{rate}+{\mathrm{c}}_{\mathrm{P}}{P}_{rate}+{\mathrm{c}}_{\mathrm{B}}{E}_{rate}\right)\frac{\mathsf{\sigma }{\left(1+\mathsf{\sigma }\right)}^{\mathrm{Y}}}{{\left(1+\mathsf{\sigma }\right)}^{\mathrm{Y}}-1}$$

(7)

where c_E ($/kW·h) is the unit capacity price of ESS, c_P ($/kW) is the PCS cost, which are expressed per unit of the ESS rated power capacity and c_B ($/kW·h) is the connecting facility of the plant cost per kW·h. E_rate and P_rate are the rated capacity and rated power of the ESS, separately. In this paper, the maximum charging/discharging power of the ESS was considered as the rated power. Then, P_rate and S_ESS were equal. Y is the project cycle (year) and σ is the discount rate (%).

Typically, the lifetime of the battery and the PCS cannot meet the needs of the entire project cycle without replacement. So, annual replacement cost is expressed as follows.

$$C_2={\mathrm{c}}_{\mathrm{E}}{E}_{rate}{\displaystyle \sum _{\epsilon =1}^k\left(\frac{{\left(1-\mathsf{\beta }\right)}^{\epsilon n}}{{\left(1+\mathsf{\sigma }\right)}^{\epsilon n}}\frac{\mathsf{\sigma }{(1+\mathsf{\sigma })}^{\mathrm{Y}}}{{(1+\mathsf{\sigma })}^{\mathrm{Y}}-1}\right)}+{\mathrm{c}}_{\mathrm{P}}{P}_{rate}\frac{{\left(1-\mathsf{\beta }\right)}^{10}}{{\left(1+\mathsf{\sigma }\right)}^{10}}\frac{\mathsf{\sigma }{(1+\mathsf{\sigma })}^{\mathrm{Y}}}{{(1+\mathsf{\sigma })}^{\mathrm{Y}}-1}$$

(8)

where k is the total number of times of battery replacement ((round up, (k = Y/n−1)), n is the lifetime of the battery and ε is the sequence of replacement. β is the average annual decline rate of the ESS investment cost. According to the current literature, the lifetime of the PCS is considered to have a fixed lifetime of 10 years in this study [40].

The operation and maintenance cost C₃ include the labor and management cost, which are related to the rated power, cf is the operation and maintenance cost per kW, then

$$C_3={\mathrm{c}}_{\mathrm{f}}{P}_{rate}$$

(9)

The ESS should be decommissioned at the end of their life. The disposal cost can be formulated as [41]:

$$C_4={\mathrm{c}}_{\mathrm{d}}{P}_{rate}{\displaystyle \sum _{\epsilon =1}^k\left(\frac{{\left(1-\mathsf{\beta }\right)}^{\epsilon n}}{{\left(1+\mathsf{\sigma }\right)}^{\epsilon n}}\frac{\mathsf{\sigma }{(1+\mathsf{\sigma })}^{\mathrm{Y}}}{{(1+\mathsf{\sigma })}^{\mathrm{Y}}-1}\right)}$$

(10)

where c_d ($/kW) is a specific disposal cost of the ESS.

In general, part of the ESS equipment can be recycled, then the recovery value can be expressed as follows:

$$C_5=-{\mathrm{c}}_{\mathrm{res}}(C_1+C_2)$$

(11)

where c_res is the recover residual rate, usually selected between 3% and 5% [42].

2.1.3. The Economic Benefit Model of ESS

In this paper, the economic benefit of the ESS mainly comes from two folds. The first part is the arbitrage revenue brought by energy time shifting between the off-peak and peak prices, and the second part is the income carried by reducing line losses.

In the context of the electricity market, the ESS can achieve arbitrage through charging at off-peak time and discharging at peak time. Dividing the entire day into 24 time slots, the arbitrage income can be calculated as:

$$B_{arb}=\mathrm{G}{\displaystyle \sum _{t=1}^{24}{\displaystyle \sum _{i=1}^{{\mathrm{N}}_{\mathrm{ESS}}}(P_{ESS.dis,i}(t)}}+P_{ESS.ch,i}(t))\times price(t)$$

(12)

where B_arb is the annual arbitrage income of ESS and G is the number of days that ESS operates in a year. P_ESS.dis,i(t) and P_ESS.ch,i(t) are the discharging and charging active power of the i-th ESS, respectively. N_ESS is the number of ESS in the DNs and the price(t) is the electricity price at time t.

The allocation of ESS in DNs can effectively reduce the line losses and obtain benefits, which are shown as follows:

$$B_{\mathrm{loss}}={\mathrm{C}}_0-C_{\mathrm{loss}}$$

(13)

$$C_{\mathrm{loss}}=\mathrm{G}{\displaystyle \sum _{t=1}^{24}{\displaystyle \sum _{i=1}^{{\mathrm{N}}_{\mathrm{ac}}}{\displaystyle \sum _{j\in {\Omega }_i}^{}\frac{1}{2}({\mathrm{r}}_{\mathrm{ij}}\times I_{ij}^2(t)\times price(t))}}}$$

(14)

where B_loss is the annual income from reducing line losses. C₀ and C_loss are the annual line losses cost of DNs before and after ESS configuration separately. Additionally, N_ac is the number of system nodes and Ω_i is the set of neighbors of node i. The r_ij and I_ij(t) are the resistance and current of branch ij respectively. It is worth noting that a factor of 1/2 is applied in Equation (14) to avoid repeated line losses consideration for each line.

2.1.4. The Capacity Fading Model of ESS

Battery lifetime is usually defined as the cycle life or calendar life corresponding to the actual capacity degradation to 80% of the rated capacity. The capacity degradation is mainly caused by the corrosion of lithium and the increase of the internal resistance of the battery, which are related to the charging/discharging power, the depth of discharging, the SOC fluctuations, the number of cycles and the working temperature [43]. As shown in [44], the residual value of the Li-ion battery can be determined by the state of health (SOH), which changes from 1 to 0 when capacity declined from rated capacity to 80%. The capacity fading rate is related to the average SOC (SOC_avg) and the SOC deviation (SOC_dev), and the higher the SOC_avg or SOC_dev value, the higher the capacity fading rate [44]. The total fading capacity is a summation of the fading capacity that occurs under the experienced operating conditions, considering temperature, the charging energy, and the SOC_avg and SOC_dev, the fading capacity can be expressed as:

$$\Gamma ={\displaystyle \sum _m^{{\tau }_{sum}}(({\mathrm{x}}_{1}SO{C}_{dev,m}}·{\mathrm{e}}^{({\mathrm{x}}_2·SO{C}_{avg,m})}+{\mathrm{x}}_3{\mathrm{e}}^{{\mathrm{x}}_{4}SO{C}_{dev,m}}){\mathrm{e}}^{(-\frac{\mathrm{Ea}}{\mathrm{R}}(\frac{1}{T_m}-\frac{1}{{\mathrm{T}}_{\mathrm{ref}}}))})A{h}_m$$

(15)

where Γ is the total fading capacity and the parameters can be chosen as: x₁ is −4.092 × 10⁻⁴, x₂ is −2.167, x₃ is −1.408 × 10⁻⁵ and x₄ is 6.130, Ea is 78.06 kmol/J and R is 8.314 J/(mol·K) according to the existing literature [44]. One sampling period is represented by m and sampling duration is τ, which is selected as 24 h in this study, and τ_sum is the total time period. Ah_m, T_m, SOC_avg,m and SOC_dev,m represent the charging energy, temperature of the ESS, average SOC and SOC deviation during time period m. T_ref is the reference temperature, which is normally set as 25 °C. The temperature of Li-ion battery cells, Tm, can be formulated as Equation (16) according to the approximate linear relationship between the charging/discharging power and temperature change for a given charging profile, R_th is the empirical constant as shown in [45].

$$T_m={\int }_0^{\tau }({\mathrm{T}}_{\mathrm{ref}}+{\mathrm{R}}_{\mathrm{th}}\cdot \left|P_{dis}\left(t\right){\mu }_{dis}+P_{ch}\left(t\right){\mu }_{ch}\right|)dt/\tau$$

(16)

The SOH of the battery can be calculated as Equation (17) and the lifetime of the battery n can be calculated by the τ_sum in Equation (15) (round, n = τ_sum/365).

$$SOH=(1-\frac{\Gamma }{0.2E_{rate}})·100\%$$

(17)

2.2. Distflow Branch Power Flow Model of DNs

Due to the large resistance and reactance ratio of the DNs and difficulties in active/reactive power decoupling, the power flow modeling methods commonly used in power transmission systems are no longer applicable to DNs. More accurate power flow analysis methods are essential in developing the optimization operation of DNs. The basic constraints of DNs are modeled in the form of branch flow [33,34]. The distflow branch flow method uses active power, reactive power, branch current amplitude and node voltage amplitude to describe the power flow of DNs. It can be mathematically expressed as:

(1) AC branch flow constraints:

$$\left\{\begin{array}{c}{\displaystyle \sum _{\mathrm{h}\in {\phi }_{\mathrm{i}}}}{P}_{ih}(t)={\displaystyle \sum _{j\in {\varphi }_i}}\left(P_{ji}(t)-{\mathrm{R}}_{\mathrm{ji}}(I_{ji}^2(t))\right)+P_i(t)\\ {\displaystyle \sum _{\mathrm{h}\in {\phi }_{\mathrm{i}}}}{Q}_{ih}(t)={\displaystyle \sum _{j\in {\varphi }_i}}\left(Q_{ji}(t)-{\mathrm{X}}_{\mathrm{ji}}(I_{ji}^2(t))\right)+Q_i(t)\end{array}\right.$$

(18)

$$V_j^2(t)=V_i^2(t)-{2(\mathrm{R}}_{\mathrm{ij}}{P}_{ij}(t)+{\mathrm{X}}_{\mathrm{ij}}{Q}_{ij}(t\left)\right)+{(\mathrm{R}}_{\mathrm{ij}}^2+{\mathrm{X}}_{\mathrm{ij}}^2\left)\right(I_{ij}^2(t\left)\right)$$

(19)

$$I_{ij}^2(t)=\frac{(P_{ij}^2(t)+Q_{ij}^2(t))}{V_i^2(t)}$$

(20)

where φ_i is the set of branch head nodes with node i as the end node and Φ_i is the set of branch end nodes with node i as the head node; R_ji and X_ji are the resistance and reactance of branch ij. P_ij(t) and Q_ij(t) are the active and reactive power flowing between node i to node j at time t. I_ij(t) and V_i(t) are the current flowing from node i to node j and the voltage of node i at time t. P_i(t) and Q_i(t) are the sum of the active power and reactive power injected on the node i at time t respectively. Their expressions are as follows:

$$\left\{\begin{array}{l}{P}_i\left(t\right)=P_{PV,i}^{}\left(t\right)+P_{ESS.\mathrm{dis},i}^{}\left(t\right)+P_{ESS.ch,i}^{}\left(t\right)-P_{LOAD,i}^{}\left(t\right)\hfill \\ Q_i\left(t\right)=Q_{PV,i}^{}\left(t\right)+Q_{ESS,i}^{}\left(t\right)-Q_{LOAD,i}^{}\left(t\right)\hfill \end{array}\right.$$

(21)

where P_PV,i(t) and Q_PV,i(t) are the active and reactive power injected by the PV on the node i at time t respectively. Q_ESS,i(t) is the reactive power injected by ESS on node i. P_LOAD,i(t) and Q_LOAD,i(t) are the active and reactive power consumed by the load on the node i at time t.

(2) Constraint of operating voltage level:

$${({\mathrm{V}}_{\mathrm{i}}^{\min })}^2\le (V_i^2(t))\le {({\mathrm{V}}_{\mathrm{i}}^{\max })}^2$$

(22)

where ${\mathrm{V}}_{\mathrm{i}}^{\max }$ and ${\mathrm{V}}_{\mathrm{i}}^{\min }$ are the upper and lower limits of the voltage of node i, respectively.

(3) Branch current constraint:

$$0\le (I_{ij}^2(t))\le {({\mathrm{I}}_{\mathrm{ij}}^{\max })}^2$$

(23)

where ${\mathrm{I}}_{\mathrm{ij}}^{\max }$ is the maximum current of branch ij.

3. The Bilevel ESS Optimization Model

The ESS planning issue is closely related to the operation pattern at different time scales. The bilevel optimization model (as shown in Figure 2) can effectively decouple the energy storage capacity planning issue and the operation optimization problem. The long term planning issue can be solved in the outer optimization, while the optimal operation problem can be resolved in the inner optimization process. The outer model mainly concentrates on the capacity of the chosen ESS and the economic operation of the system is settled by the inner optimization with the ESS capacity, rated power and initial SOC determined by the outer optimization process. The internal optimization model optimizes the operation of DNs under the given ESS parameters of the external layer, obtains the operation power of the ESS, the cost of line losses and the revenue from the ESS arbitrage, then returns to the external layer and solves the problem iteratively through the internal and external optimization. The calculation scale and difficulty of the model can be effectively reduced by the iterative approach.

3.1. The Outer Optimization Model

3.1.1. The Outer Optimization Model Establishment

As shown in Equation (24), the proposed outer optimization model takes the total LCC of the ESS minus the annual average benefit of line loss reduction and arbitrage income as the optimization object. The average annual arbitrage income of the ESS and average annual line losses cost are calculated by the internal optimization model and returned to the external layer. The optimization variables are the rated ESS capacity, rated power and the initial SOC value (SOC0) of each ESS. Since SOC0 has a certain impact on the operation of the ESS, thereby affecting the service life and arbitrage benefits of the ESS, thus it is used as the optimization variable.

$$\mathrm{F}=\min (C\mathrm{ost}-B_{arb}-B_{loss})$$

(24)

where the cost is obtained from Equation (6). B_arb and B_loss are described by Equations (12) and (13) respectively.

The outer optimization model must meet the following constraints:

$$\left\{\begin{array}{l}{\mathrm{P}}_{\mathrm{ESS}}^{\min }\le P_{rate}^{}\le {\mathrm{P}}_{\mathrm{ESS}}^{\max }\hfill \\ {\mathrm{E}}_{\mathrm{ESS}}^{\min }\le E_{rate}^{}\le {\mathrm{E}}_{\mathrm{ESS}}^{\max }\hfill \\ SO{C}_{\min }\le SOC0\le {\mathrm{SOC}}_{\max }\hfill \end{array}\right.$$

(25)

where ${\mathrm{P}}_{\mathrm{ESS}}^{\min }$ and ${\mathrm{P}}_{\mathrm{ESS}}^{\max }$ are the minimum and maximum investment power of ESS respectively. ${\mathrm{E}}_{\mathrm{ESS}}^{\min }$ and ${\mathrm{E}}_{\mathrm{ESS}}^{\max }$ are the minimum and maximum investment capacity of ESS. P_rate is the maximum active power that can be released when ESS works at the unit power factor.

3.1.2. The Outer Optimization Solution

GA is a heuristic optimization algorithm, which converges to the optimal solution successively through iterations and can solve the constrained nonlinear problems effectively, however it may fall into a local optimum [46]. SA is a general probabilistic algorithm for finding the optimal solution of a proposition in a large searching domain, which makes it as a good combination with GA to determine the global optimal solution. The probability acceptance function of SA can be calculated according to (26) [47].

$$P={\mathrm{e}}^{-\frac{(fit-fit\prime )}{{\mathsf{\lambda }}^m{\mathrm{T}}_0}}$$

(26)

where, fit is the individual fitness, which is the reciprocal of the objective function in Equation (24) without the annealing operation, and fit’ is the individual fitness with SA. λ is the annealing coefficient. T₀ is the initial temperature of annealing and m is the number of times of annealing.

The variables of outer optimization model are the rated capacity, rated power and the SOC0 of the ESS. The outer optimization process is presented in Figure 3.

3.2. The Inner Optimization Model

3.2.1. Operation Strategy of the ESS

As illustrated in Section 2.1.1, ESS can operate in any quadrant as shown in Figure 1. Through the power control of the PCS, the DNs operation can be adjusted to a certain extent. Therefore, the ESS in the DNs can further improve the stability of the system and solve the voltage issues caused by the high permeability of renewable generation in the system. In this model, the PV sends active power according to a unit factor to maximize the use of renewable energy. The PCS of the ESS has a certain reactive power capacity, which can absorb or release a certain amount reactive power for reactive compensation to support the voltage profiles [48].

As shown in Figure 4, V and V_S are the voltages of grid-connected point and DN voltage, R and X are the line resistance and reactance; P_PV and Q_PV are the active and reactive power of the PV, respectively, P_ESS and Q_ESS are the active and reactive power of the ESS and P_S and Q_S are the sum of the injected active and reactive power of PV and ESS, respectively.

Therefore:

$$\left\{\begin{array}{c}{P}_S=P_{PV}+P_{ESS}\\ Q_S=Q_{PV}+Q_{ESS}\end{array}\right.$$

(27)

The grid connected point voltage V after PV and ESS is:

$$V=V_S+\frac{P_S\mathrm{R}+Q_S\mathrm{X}}{V}$$

(28)

The voltage variation ΔV_l before and after the PV and ESS access can be calculated as follows:

$$\Delta V_l=\frac{P_S\mathrm{R}+Q_S\mathrm{X}}{V}=\frac{(P_{PV}+P_{ESS})\mathrm{R}+(Q_{PV}+Q_{ESS})\mathrm{X}}{V}$$

(29)

When the ESS is connected without reactive power, the voltage variations are caused by the active power injection of PV. At this time, P_ESS is negative, and the ESS works as a load to reduce P_S. It can be seen from Equation (29) that the smaller the P_S, the smaller ΔV_l, thereby regulating the voltage to a safe range. On the other hand, at a time that the PV output is large, the ESS is used as a load to absorb power from the network to regulate the voltage within limits, the reactive power of the ESS can be utilized to reduce the unnecessary active power purchase from the utility during the period of high electricity price, therefore improve the arbitrage income of the ESS. It can be easily seen from Equation (30), the reactive power contribution from ESS, Q_ESS, is essential for voltage regulation while reducing unnecessary ESS’s active power exchange with the grid, consequently improve the arbitrage income.

$$\Delta V_l=\frac{(P_{PV}+P_{ESS})\mathrm{R}}{V}=\frac{(P_{PV}+P_{ESS}^{\prime })\mathrm{R}+Q_{ESS}^{}\mathrm{X}}{V}$$

(30)

As a result, when the PV power output causes the voltage excursions at a high electricity price period, the PCS with ESS preferentially provide reactive power to regulate the voltage level until reaching the reactive power adjustment capacity margin, then active power support is applied.

3.2.2. Inner Optimization Model Establishment

The voltage levels can be divided into the dead zone, uneconomic zone and emergency zone, the best operating state for DNs is the dead zone, however the uneconomic zone operation does not affect the security of the system [49]. When the system voltage operates in the emergency zone, the safety and reliability of the system cannot be guaranteed. According to the operation strategy of ESS illustrated in the previous section, the objective function of the inner optimization model can be constructed as Equations (31) and (32) to ensure the best operation status of voltage, the losses minimization and the revenue maximization.

$$\mathrm{f}=\min ({\mathsf{\lambda }}_1{C}_{\mathrm{loss}}^{}-{\mathsf{\lambda }}_2{B}_{arb}^{}+{\mathsf{\lambda }}_3\Delta V)$$

(31)

$$\left\{\begin{array}{l}\Delta \mathrm{V}=\mathrm{G}{\displaystyle \sum _{\mathrm{t}=1}^{24}{\displaystyle \sum _{\mathrm{i}=1}^{{\mathrm{N}}_{\mathrm{ac}}}\left|{\mathrm{V}}_{\mathrm{i}}^2\left(\mathrm{t}\right)-1\right|}}\hfill \\ \begin{array}{cc}s.t.& ({\mathrm{V}}_{\mathrm{i}}^{}\left(\mathrm{t}\right)\ge {\mathrm{V}}^{\mathrm{thr},\max }\Vert {\mathrm{V}}_{\mathrm{i}}^{}\left(\mathrm{t}\right)\le {\mathrm{V}}^{\mathrm{thr},\min })\end{array}\hfill \end{array}\right.$$

(32)

where C_loss is the annual line losses cost of the DNs after ESS configuration, represented as Equation (12) and Barb is the annual arbitrage income of ESS, which is represented as Equation (14). ΔV is the voltage variation and λ₁, λ₂ and λ₃ are the weighting coefficients. V_i(t) is the voltage amplitude of node i at time t and V^thr,max, V^thr,min are the upper and lower limits of the optimal interval of nodal voltage amplitude, respectively. On the premise that the ESS performs loss reduction and arbitrage to the greatest extent to improve economic efficiency, the voltage is approached or maintained within the optimized interval [V^thr,min, V^thr,max] as much as possible.

The inner optimization model is essentially a multiperiod optimal power flow problem at the operating layer. So, the distflow branch flow method of DN is used in this paper. The constrains of the DNs are explained by Equations (18)–(23). Additionally, the operating constraints of ESS are illustrated in Equations (1)–(5).

3.2.3. Inner Optimization Model Solution

The continuous power regulation of ESS expands the operation optimization problem of DNs from a single period of time to a continuous time series, which brings a heavy burden to the calculation. Compared with other common algorithms, second-order cone programming features with faster solving speed and stronger optimization capability [50]. In this method, the nonconvex nonlinear problem is transformed into a second-order cone programming problem, which can be solved efficiently by convex relaxation, and then the global optimal solution of the original problem can be obtained.

However, the requirements of SOCP for the mathematical model are very strict, it mainly can be summarized as the following two points [51]:

(1) The objective function of the optimization model must be a linear function about the decision variables.

(2) The feasible domain of the optimization model must be composed of nonlinear cone constraints and linear equation or inequality constraints.

The optimization model of the DNs with ESSs contains the quadratic term and the product term of the decision variable. Therefore, it is necessary to convert the corresponding model to meet the requirements of the linear objective function and cone searching space. Since Equation (2) contains 0–1 variables, the second-order cone relaxation technique is used to the power flow equations in this optimization model and then the original model is converted to a mixed integer second-order cone programming (MISOCP) model, which can be solved efficiently by CPLEX [52].

The ESS optimization model of DNs is transformed, which is divided into three steps as follows:

(1)

Linear transformation of objective function

For Equations (14) and (32) containing the second order terms $I_{ij}^2(t)$ and $V_i^2(t)$, $I_{ij,2}^{}(t)$ and $V_{i,2}^{}(t)$ are utilized to replace them. Then the objective function can be linearized as Equations (33) and (34)

$$C_{\mathrm{loss}}=\mathrm{G}{\displaystyle \sum _{t=1}^{24}{\displaystyle \sum _{i=1}^{{\mathrm{N}}_{\mathrm{ac}}}{\displaystyle \sum _{j\in {\Omega }_i}\frac{1}{2}({\mathrm{r}}_{\mathrm{ij}}\times I_{ij,2}^{}(t)\times price(t))}}}$$

(33)

$$\left\{\begin{array}{l}\Delta \mathrm{V}=\mathrm{G}{\displaystyle \sum _{\mathrm{t}=1}^{24}{\displaystyle \sum _{\mathrm{i}=1}^{{\mathrm{N}}_{\mathrm{ac}}}\left|{\mathrm{V}}_{\mathrm{i},2}^{}\left(\mathrm{t}\right)-1\right|}}\hfill \\ \begin{array}{cc}s.t.& ({\mathrm{V}}_{\mathrm{i}}^{}\left(\mathrm{t}\right)\ge {\mathrm{V}}^{\mathrm{thr},\max }\Vert {\mathrm{V}}_{\mathrm{i}}^{}\left(\mathrm{t}\right)\le {\mathrm{V}}^{\mathrm{thr},\min })\end{array}\hfill \end{array}\right.$$

(34)

Since Equation (34) contains the absolute value term $\left|V_{i,2}^{}\left(t\right)-1\right|$, the auxiliary variable ${\mu }_i^{}\left(t\right)=\left|V_{i,2}^{}\left(t\right)-1\right|$ is introduced and the constraint is added as follows [53]:

$$\left\{\begin{array}{l}{\mu }_i^{}\left(t\right)\ge 0\hfill \\ {\mu }_i^{}(t)\ge V_{i,2}^{}(t)-{({\mathrm{V}}^{\mathrm{thr},\max })}^2\hfill \\ {\mu }_i^{}(t)\ge -V_{i,2}^{}(t)+{({\mathrm{V}}^{\mathrm{thr},\min })}^2\hfill \end{array}\right.$$

(35)

After the above transformation, the objective function becomes the linear function of the decision variable.

(2)

Cone transformation of system operation constraints

The operation constraints of the system also replace the second order term for linearization. Replace the second order terms $I_{ij}^2(t)$ and $V_i^2\left(t\right)$ in Equations (18)–(20), (22) and (23) with $I_{ij,2}^{}\left(t\right)$ and $V_{i,2}^{}\left(t\right)$, as shown in (36) to (40).

$$\left\{\begin{array}{c}{\displaystyle \sum _{\mathrm{h}\in {\phi }_{\mathrm{i}}}{P}_{ih}(t)={\displaystyle \sum _{j\in {\varphi }_i}\left(P_{ji}\left(t\right)-{\mathrm{R}}_{\mathrm{ji}}(I_{ji,2}^{}(t))\right)}}+P_i(t)\\ {\displaystyle \sum _{\mathrm{h}\in {\phi }_{\mathrm{i}}}{Q}_{ih}(t)={\displaystyle \sum _{j\in {\varphi }_i}\left(Q_{ji}\left(t\right)-{\mathrm{X}}_{\mathrm{ji}}(I_{ji,2}^{}(t))\right)}}+Q_i(t)\end{array}\right.$$

(36)

$$V_{j,2}^{}(t)=V_{i,2}^{}(t)-{2(\mathrm{R}}_{\mathrm{ij}}{P}_{ij}(t)+{\mathrm{X}}_{\mathrm{ij}}{Q}_{ij}(t\left)\right)+{(\mathrm{R}}_{\mathrm{ij}}^2+{\mathrm{X}}_{\mathrm{ij}}^2\left)\right(I_{ij,2}^{}(t\left)\right)$$

(37)

$$I_{ij,2}^{}(t)=\frac{(P_{ij}^{}{\left(t\right)}^2+Q_{ij}^{}{\left(t\right)}^2)}{V_{i,2}^{}\left(t\right)}$$

(38)

$${({\mathrm{V}}_{\mathrm{i}}^{\min })}^2\le (V_{i,2}^{}(t))\le {({\mathrm{V}}_{\mathrm{i}}^{\max })}^2$$

(39)

$$0\le (I_{ij,2}^{}(t))\le {({\mathrm{I}}_{\mathrm{ij}}^{\max })}^2$$

(40)

Further relaxing the constraint condition in (38), as shown in Equation (41), then, it is transformed into the form of standard second-order cone as shown in Equation (42):

$$I_{ij,2}^{}(t)\ge \frac{(P_{ij}{\left(t\right)}^2+Q_{ij}{\left(t\right)}^2)}{V_{i,2}^{}\left(t\right)}$$

(41)

$${‖\begin{array}{c}2P_{ij}\left(t\right)\\ 2Q_{ij}\left(t\right)\\ I_{ij,2}^{}\left(t\right)-V_{i,2}^{}\left(t\right)\end{array}‖}_2\le I_{ij,2}^{}(t)+V_{i,2}^{}(t)$$

(42)

(3)

Cone transformation of ESS operation constraints

The charging and discharging constraints of the ESS are nonlinear, which do not meet the requirements of the cone optimization model for linear constraints or the constraint form of the rotating cone, it is necessary to do other model transformation. So, Equation (1) is transformed into a rotating cone constraint:

$${(P_{ESS.dis}(t)+P_{ESS.ch}(t))}^2+Q_{ESS}{\left(t\right)}^2\le 2\frac{S_{ESS}}{\sqrt{2}}\times \frac{S_{ESS}}{\sqrt{2}}$$

(43)

Through the equivalent transformation, the variables in Equation (43) can constitute the Cartesian product form of the rotating cone, which meets the strict requirements of the searching space in the convex cone range. Equations (2)–(5) are linear constraints, which meet the requirements of cone optimization for linear constraints.

After the above steps, the cone conversion process of the optimization model is completed, the objective function is linearized and the nonlinear constraints are transformed into linear constraints, second-order cone constraints, or rotating cone constraints. According to [54,55], if the second-order cone relaxation is to be established in the power system optimization model, certain conditions must be met, including: (1) Optimization problems must have feasible solutions. (2) The objective function must be a strict increasing function of branch current. (3) The network topology can be represented by a tree-like connectivity graph. The infinite norm defining the slack deviation of the second-order cone is shown in Equation (44) to verify the accuracy of Equation (42) at the optimal solution after relaxation [56]. When the slack deviation is small enough, the second-order cone relaxation can be considered to meet the calculation accuracy requirements.

$$gap={‖I_{ij}^{}(t)-\frac{{(P_{ij}^{}(t))}^2+{(Q_{ij}^{}(t))}^2}{V_i^{}\left(t\right)}‖}_{\infty }$$

(44)

The objective function of the inner model includes the cost of network losses, the arbitrage income and the voltage deviation, and the MISOCP is no more accurate when the optimal power flow is solved by the cone relaxation under the scene of high permeability of PV access. After the nth solution of the optimization model, if the gap_n representing the slack deviation is small enough, then the solution accuracy can be considered satisfied, otherwise the cutting plane constraint can be added when the nth solution is solved. The cutting plane constraints of the multiperiod optimal power flow are shown in Equation (45) [57].

$${\displaystyle \sum _{t=1}^{24}{\displaystyle \sum _{i=1}^{{\mathrm{N}}_{\mathrm{ac}}}{\displaystyle \sum _{j\in {\Omega }_i}^{}{\mathrm{R}}_{\mathrm{ij}}{I}_{ij,2}^{}\left(t\right)}}}\le {\displaystyle \sum _{t=1}^{24}{\displaystyle \sum _{i=1}^{{\mathrm{N}}_{\mathrm{ac}}}{\displaystyle \sum _{j\in {\Omega }_i}^{}{\mathrm{R}}_{\mathrm{ij}}\frac{{(P_{ij,n}^{}(t))}^2+{(Q_{ij,n}^{}(t))}^2}{V_{i,2,n}^{}\left(t\right)}}}}$$

(45)

where P_ij,n(t) and Q_ij,n(t) are the active and reactive power transmitted between branches ij at time t, which are calculated from the nth iterations. V_i,2,n(t) is the voltage square of node i at time t calculated for the nth time and I_ij,2(t) is the current square of branch ij at time t calculated for the nth time.

The solution flow of the MISOCP used in this paper is shown in Figure 5:

4. Case Studies

4.1. Profiles

In this section, the proposed method is implemented on a medium voltage DN as shown in Figure 6 for ESS planning. The parameters of the case study network are shown in

Table A1. System data.

Node i	Node j	R(Ω)	X(Ω)	Node i	Node j	R(Ω)	X(Ω)
1	2	0.062207	0.133984	63	64	0.195590	0.421271
1	3	0.088866	0.191404	64	65	0.170937	0.368172
1	7	0.106640	0.229686	65	66	0.022633	0.048749
3	4	0.071094	0.153125	160	67	0.065168	0.140363
3	5	0.115527	0.248828	67	68	0.051511	0.110947
5	6	0.088866	0.191404	67	72	0.046828	0.100860
7	8	0.210548	0.453489	67	97	0.065168	0.140363
8	12	0.039990	0.086133	68	69	0.077018	0.165885
8	9	0.039990	0.086133	69	70	0.065168	0.140363
8	13	0.146249	0.314998	70	71	0.065168	0.140363
9	14	0.151073	0.325388	72	73	0.037463	0.080689
13	34	0.229168	0.493593	72	76	0.409082	0.881100
14	11	0.088866	0.191404	73	74	0.094792	0.204167
14	10	0.088866	0.191404	74	75	0.075833	0.163334
15	16	0.133301	0.287110	76	77	0.280567	0.604298
15	17	0.124414	0.267969	77	78	0.182271	0.392585
34	15	0.349665	0.753125	78	79	0.042656	0.091875
149	1	0.071094	0.153125	78	80	0.450259	0.969790
149	150	0.154535	0.332844	80	81	0.426561	0.918747
152	52	0.072428	0.156000	81	82	0.041471	0.089323
52	53	0.022633	0.048749	82	83	0.052135	0.112290
53	54	0.049794	0.107249	81	84	0.042656	0.091875
54	55	0.063915	0.137663	84	85	0.042656	0.091875
54	57	0.049794	0.107249	97	98	0.071093	0.153124
55	56	0.131489	0.283207	98	99	0.213283	0.459378
57	58	0.135803	0.292499	99	100	0.521350	1.122909
57	60	0.131489	0.283207	100	450	0.056874	0.122499
58	59	0.101548	0.218719	450	451	0.149851	0.322757
60	61	0.150452	0.324052	13	152	0.053320	0.114844
60	620	0.105318	0.226839	The following is the interconnection switch
620	62	0.210636	0.453678
62	63	0.255770	0.550891	60	160	0.0936	0.0936

, Appendix A. In this system, the base voltage was 10 kV, with the tight upper/lower limits of the voltage being 1.03 p.u. and 0.97 p.u., respectively. The upper/lower limits of the optimal range of nodal voltage amplitude were set as 1.015 p.u. and 0.985 p.u. The node 451 and the 150 were slack buses. The maximum capacity of each branch was 500 A and the sampling time was 1 h. The daily profiles of the loads and PV are shown in Figure 7. Time of use tariff was: 0.139$/(kW·h) from 08:00 to 21:00 (peak), 0.042$/(kW·h) from 00:00 to 07:00 and 22:00 to 24:00 (off-peak).

In this study, the Li-ion battery was selected as the ESS and the parameters of the ESS are listed in

Table 1. Technical and economic parameters of ESS.

Technical parameters	ηC/ηD	ηPCS	SOCmin	SOCmax
	0.976	0.95	0.1	0.9
Economic parameters	CE	CP	CB	Cf
	156$/kW·h	154$/kW	0$/kW·h	22$/kW·year
	Cd	Cres	σ	β
	224$/kW	5%	10%	0%

according to the current literature [40,41,42,58,59]. The project life cycle was determined as 20 years. For the weighting coefficient λ1 of the distribution system network loss cost, λ2 of arbitrage income and λ3 of voltage level deviation in the objective function, the analytic hierarchy process (AHP) was used to obtain λ1 = 0.7, λ2 = 0.15 and λ3 = 0.15. The parameters of the GA and SA are listed in

Table 2. Parameters of the genetic algorithm (GA) and the simulated annealing (SA).

Parameters of GA				Parameters of SA
Maximum generation	60	Variation probability	0.05	Annealing coefficient	0.9
Size of population	30	Crossover rate	0.7	Initial temperature	100

. When the maximum gap_n was 10⁻⁴, it could be considered to meet the requirements of slack deviation.

4.2. Simulation Results and Analysis

The preselected locations of the ESS include the critical load nodes, the end of the feeder, the DG bus bars, the primary substation, etc. So, the ESS configuration nodes in this medium voltage case study distribution system are presented in Figure 6. Three sets of simulation are carried out in this part: case 1 uses the bilevel optimization model proposed in this paper to configure the ESS of the medium voltage case study distribution system considering the active and reactive power of ESS at the same time; case 2 compare the configuration results of the same model as scenario 1 without reactive power from PCS and case 3 compare the optimization result with the authors’ previous work.

Case 1: ESS configuration considering active and reactive power.

The simulation results of the bilevel optimization model considering both active and reactive power are shown in

Table 3. Optimization results of the bilevel ESS optimization model.

ESS Location	Prate (kW)	Qrate (kvar)	Erate (kW·h)	SOC0 (%)	LCC ($/year)
ESS1	500	500	3500	70	108,414
ESS2	600	600	3000	70	104,443
ESS3	600	600	3500	70	119,036

. The voltage variations of each typical nodes before and after the ESS configuration are separately presented as (a) and (b) in Figure 8. Additionally, the output of the ESS is shown in Figure 9.

According to the results shown in Figure 8, the voltage issues could be effectively solved by coordinated operation of active/reactive power of the ESS. Before ESS operation, the voltage amplitude of each typical node exceeded the upper limit from 10:00 to 14:00, and the voltage amplitude of nodes 83 and 84 exceeded the lower limit at 18:00 to 22:00 as shown in Figure 8a. Additionally, the voltage profiles shown in Figure 8b of each typical node was well kept within the optimization range after the ESS participation. In general, the charging/discharging active power of the ESS were consistent with the power supply and demand, and the arbitrage income from discharging/charging at different electricity price were carried out within the voltage safety margin, as shown in Figure 9. Overall, the arbitrage income was $177,531 and the benefit of reducing line losses is $2567 annually.

Case 2: ESS configuration considering active power only.

The results of the bilevel optimization model of the ESS without reactive power are given in

Table 4. Optimization results of the bilevel ESS optimization model without reactive power.

ESS Location	Prate (kW)	Qrate (kvar)	Erate (kW·h)	SOC0 (%)	LCC ($/year)
ESS1	500	0	4000	60	143,855
ESS2	500	0	3500	60	100,864
ESS3	500	0	5000	70	132,998

. Comparing the results of Table 3 and Table 4, it can be seen that the rated capacity and LCC of the configuration scheme without reactive power increased. As shown in Figure 10 and Figure 11, the voltage issues could be effectively solved by active power only, and the arbitrage income from energy time shifting was achieved while maintaining the voltage within the safety margin. Compared with Figure 9 and Figure 10, it can be seen that in the period with high PV injection, ESS absorb more active power with a high electricity price to solve the voltage issues, thus reducing the economy of the overall scheme. Additionally the voltage fluctuation is greater than the result shown in Figure 8b. Due to the constraints of SOC illustrated in (4)–(5), although the voltage levels are low between 16:00 and 22:00, ESS2 cannot discharge at maximum capacity to achieve better revenues during this peak hours, and has to discharge between 22:00 and 23:00 during this off-peak hours. Therefore, the overall benefit of the scheme was reduced. The arbitrage income was $106,177 and the benefit of reducing line losses was $14,069 annually.

To better quantify the economic benefit of Case 1 and Case 2, the LCC of the ESS with/without reactive power are shown in Figure 12. ESS1, ESS2 and ESS3 represent each ESS with reactive power operation and ESS1′, ESS2′ and ESS3′ represent each ESS without considering reactive power. For ESS1 and ESS3, it can be seen that reactive power can effectively reduce the LCC, and the investment and replacement cost obviously. For ESS2, the LCC of ESS2 was slightly higher than ESS2′, but a much larger arbitrage revenue was also obtained. In general, the total cost of ESS with reactive power was about $105,676 lower than that without reactive power annually. In other words, the allocation scheme of ESS with reactive power compensation saved 41% of the average annual cost compared with the scheme without considering reactive power contribution.

Case 3: ESS configuration comparison with two-stage heuristic optimization method.

Comparing the optimization method with the authors’ previous research proposed in [60], results are shown in

Table 5. Optimization results of the two-stage heuristic optimization method.

ESS Location	Prate (kW)	Qrate (kvar)	Erate (kW·h)	SOC0 (%)	LCC ($/year)
ESS1	500	0	2290	57	146,966
ESS2	375	0	2362	10	118,994
ESS3	500	0	2633	51	168,590

. The two models were designed to solve the problem of DNs voltage fluctuations caused by the increase of DG penetration, but ESS arbitrage operation and the reactive power of PCS are not considered in [60]. It worth noting that the configuration capacity of each ESS in the two-stage heuristic optimization method was smaller than the model proposed in this paper, however the overall cost was much higher. This is mainly because the ESS of the two-stage heuristic optimization method did not carry arbitrage operation and added the electricity bill into the LCC. The electricity bill was calculated in the same way as arbitrage income and the electricity bill of each ESS was $69,479, $45,509 and $82,112 respectively. Compared with the model proposed in reference [60], the model proposed in this paper saved 65% of the annual average cost, so the economy of ESS can be greatly improved through the arbitrage operation of ESS.

Figure 13 shows the comparison of the optimization results of each ESS under different scenario. It can be seen that although the rated power and capacity of each ESS in Case 3 is the smallest of the three cases, the economic benefit of each ESS in Case 3 was the worst because there was no arbitrage income and the LCC was the highest. Under the configuration method of Case 3, each ESS only had active power to regulate the system power flow without arbitrage operation. In addition, it can be seen from Figure 7 that the output of PV was high in the period peak hour. At this time, the ESS performed load operation to absorb power and suppress PV output, resulting in a high electricity bill.

From Figure 13a, it can be seen that in Case 1, the regulating pressure of its active power was reduced due to the reactive power output of ESS1, thus the lifetime was increased compared with Case 2. Although the arbitrage income of ESS1 in Case 1 and Case 2 was similar, the economic benefit of ESS1 in case1 was the best because of its long lifetime, less rated capacity and lower LCC.

For ESS2, as shown in Figure 13b, the lifetime in Case 1 was less than that in Case 2. This is because the reactive power of ESS2 in case1 regulated the voltage so that its active power can carry out arbitrage operation with greater active power output to increase the arbitrage income, thus its SOC fluctuation was greater and its lifetime was reduced. Although the LCC of ESS1 in Case1 was slightly higher than that of Case 2, its arbitrage income was significantly higher than that of Case 2. Therefore, the economic benefit of ESS2 in Case 1 had the best results among others.

It can be seen from Figure 13c that in Case 1, the regulation pressure of its active power was reduced due to the reactive power output of ESS3, and its rated capacity was greatly reduced compared with Case 2. Similar to ESS2, because the reactive power of ESS3 in Case 1 regulated the voltage, its active power could carry out arbitrage operation with greater active power output to increase the arbitrage income, so its SOC fluctuation was greater and its lifetime was reduced. According to Figure 12 and Figure 13c, although the lifetime of ESS3 in Case1 was shorter than that of Case2, its rated capacity was less, and the increase value of replacement cost caused by the decrease of its lifetime was far less than the decrease value of installation cost. The reactive power operation of ESS3 not only increased the arbitrage income but also reduced the LCC, so the economic benefit of ESS3 in Case1 was the best compared with others.

Therefore, the model proposed in this paper was more beneficial than the two-stage heuristic optimization method in terms of economy and voltage profile improvement.

5. Conclusions

To solve the voltage fluctuations in DNs caused by the proliferation of PV, an interactive bilevel ESS planning strategy was proposed considering the factors of LCC, the arbitrage income and the benefit of reducing line losses. The outer model was solved by GA combined with SA, and the inner optimization was solved by MISOCP. The simulation results demonstrated the feasibility and effectiveness of the proposed planning method and the contributions of this study could be summarized as follows.

(1) The intermittence of the distributed renewable generation integrated in the distribution grids led to the risk of voltage fluctuations, energy storage could effectively manage the impacts without added carbon footprint. The technoeconomic analysis carried out in this study verified that the overall LCC of ESS could be significantly reduced by providing ancillary services, reducing line losses and gaining revenues from arbitrage.

(2) The coordinated operation of active and reactive power from ESS and associate PCS could effectively reduce the LCC of ESS without further investment.

(3) The accuracy and the efficiency could be achieved in solving the internal and external optimization iteratively by GA and SA, combined with MISOCP.

(4) Compared with the model in [60], the economy of this model was greatly improved after taking into account the ESS arbitrage operation and reactive power.