A Hybrid Method for Optimal Siting and Sizing of Battery Energy Storage Systems in Unbalanced Low Voltage Microgrids

Abstract

This paper deals with the problem of optimal allocation (siting and sizing) of storage resources in unbalanced three-phase low voltage microgrids. The siting and sizing problem is formulated as a mixed, non-linear, constrained optimization problem whose objective function deals with economic issues and whose constraints involve technical limitations of both network and distributed resources. Emphasis is given to the power quality issue with particular attention to unbalance reduction and voltage profile improvement. Technological issues, such as those related to the preservation of batteries’ lifetime, were also taken into account. The planning problem is solved by means of a genetic algorithm which includes an inner algorithm based on sequential quadratic programming. In order to limit the processing time while maintaining reasonable accuracy, the genetic algorithm search space is significantly reduced identifying a subset of candidate buses for the siting of the storage resources. The Inherent Structure Theory of Networks and the Loading Constraints Criterion were used to identify the candidate buses. The proposed method has been applied to a low voltage test network demonstrating the effectiveness of the procedure in terms of computational burden while also preserving the accuracy of the solution.

1. Introduction

One of the most important drivers of a sustainable electrical power grid is the optimal integration of its distributed energy resources. The optimal integration requires that the value proposition of the grid resources be attractive, that is their applications must provide multiple benefits. This should occur in the case of all the types of distributed resources and particularly when these resources are electrical energy storage systems for which high installation cost must be faced. The use of storage in the frame of distribution networks provides multiple benefits, such as energy cost reduction and power quality (PQ) enhancement [1,2,3,4]. The PQ improvement can be even more useful and evident in the case of low voltage (LV) distribution networks where one of the most severe PQ problem is the voltage unbalances caused by the large presence of unbalanced lines and single-phase loads and resources. Then, recently, particular attention has been devoted to the application of storage devices to reduce unbalances [5,6,7,8].

In the relevant literature it was shown that the use of distributed energy storage systems (DESSs) can become more attractive only if they are optimally sited and sized so that their planning stage becomes a strategic starting point for the best exploitation of their technology. An increasing interest has been also addressed to the formulation of models and algorithms finalized to the optimal sizing and siting of storage based resources with reference to different applications, [9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30]. More specifically, a method used to determine the cost of batteries that makes profitable the use of DESSs in distribution networks is discussed in [9]. Ref. [10] proposes a method for the optimal allocation and economic analysis of battery energy systems in microgrids. The method proposed in [11] allows allocating DESSs while improving network operation and minimizing costs. In [12] the life cycle cost of distribution networks in presence of batteries is minimized. In [13] two methods are discussed to optimize the allocation of batteries for frequency regulation and peak shaving in transmission and distribution networks, respectively. In [14] the allocation of batteries is performed with the aim of reducing the energy exchanged at the distribution substation while maximizing the economic benefits of the distributor. A procedure for finding the minimum storage sizes in different locations of distribution networks is proposed in [15] with the aim of reducing the curtailment from distributed generation (DG). In [16] DESS allocation is aimed at alleviating the negative impact of high penetration of photovoltaic systems in distribution networks. The optimal allocation proposed in [17] is focused on the voltage support, losses minimization, cost of energy and management of congestions. In [18] a method is proposed to reduce losses and deferring network upgrading. In [19] the storage systems are optimally sized and sited in order to provide active and reactive support in view of reactive power/voltage regulation, losses reduction and cost minimization. In [20] an approach is proposed for the optimal planning of power networks in the presence of storage devices used to control variable energy resources. In [21], the planning problem of active distribution networks is dealt with by integrating energy storage systems in the network considering both long-term investment cost and reliability of operation. A two-stage stochastic model is proposed in [22] for the optimal planning of DESSs in power systems under the uncertainties due to wind generators which integrates a tradeoff analysis between the investment and operating cost. In [23] a stochastic optimization approach is used in the framework of the optimal planning of a network in presence of electrical energy storage devices and aimed at maximizing the total expected net benefits for the system operator. A cost minimization is used in [24] to find the optimal battery size in active distribution networks. Based on a multi-temporal constrained optimization problem, in [25] the sizing and placement of storage devices are addressed by minimizing investment and operation costs. In [26,27,28,29,30,31] the planning problem has been treated with specific reference to LV grids. Particularly, these papers aim to mitigate the impact of renewable generation [26,27,28], to improve reliability [29] and to reduce costs [30,31]. In [27,28,29,30] the problem of the unbalances is also discussed within the framework of the storage planning issue.

In this paper, an unbalanced LV microgrid is considered and the allocation (siting and sizing) problem is formulated as a mixed, non-linear, constrained optimization problem whose objective function deals with economic issues and whose constraints take into account technical requirements of both network and distributed resources. This planning problem is quite complex and involves mathematical models of high dimension. As discussed in [10], in fact, the use of storage devices requires to face with problems of intertemporal nature: the operation of the storage at a time step will affect its operation at the following time steps. Moreover, computational complexity is particularly cumbersome in case of LV unbalanced grids, for which three phase representations are required for all the grid components, and when load patterns during the years have to be taken into account in the planning problem. This implies that simplified approaches or approximations in the models have to be considered in order to make the study possible even in realistic applications where the problem is further complicated by the presence of large number of buses.

Analysis of new methods that can limit the computational burden of the optimization problem is of great importance. In our opinion, a possible solution to this problem can be to reduce the computational burden by containing the solution’s search space of the planning problem into a limited number of buses of the grid, which are considered as the “best” for allocating DESSs (“allocation bus set”). Therefore, in this paper the computational complexity is overcome by solving the mixed, non-linear, constrained optimization problem of DESS allocation, through a two-step hybrid procedure: (step 1) determine the limited subset of buses to be included in the reduced search space (allocation bus set); and (step 2) find the optimal solution in this reduced region by applying a genetic algorithm (GA).

The buses of the allocation bus set (step 1) are obtained by selecting among all the network buses suitable for the DESS siting (i.e., theoretically all the network buses) a subset of buses that includes only the buses whose injection has the expected better favorable impact on the network constraints satisfaction. These buses are selected by applying the Inherent Structure Theory of Networks (ISTN) and the Loading Constraints Criterion (LCC). The ISTN is based on the reformulation of network admittance and impedance matrices with respect to their eigenvalues and eigenvectors and was first proposed by Laughton in [32,33] and successively extended to solve other important problems in the field of power systems [34,35,36,37,38,39]. The LCC was proposed in the frame of the European project ADDRESS to identify different sub-networks of distribution systems, referred to as Load Areas, by evaluating the influence of prosumers’ injections on the relevant network constraints [40,41]. These techniques have never been applied to the DESS allocation problem so far.

The GA (step 2) includes an inner algorithm that performs an optimization aimed at finding the optimal daily charge/discharge cycles of storage batteries.

The main original contribution of this paper is the proposal of a complex and detailed planning problem for DESS allocation in an unbalanced LV microgrid which drastically reduces the computational burden. The complexity is mainly related to the use of a three-phase model used to solve Power Quality problems such as voltage deviations and unbalances, and to meet technical constraints such as limits on line currents, transformer and converter capability, storage’s rating and lifetime requirements. All of these complexities are handled in an optimization framework which minimizes the overall operation and planning costs. The significant reduction of the computational burden is obtained by containing the solution’s search space of the planning problem. This allows preserving the accuracy of the model which is based on a mixed, non-linear, constrained optimization problem.

The remainder of the paper is organized as follows. Section 2 shows how the planning problem has been formulated. Section 3 refers to the algorithm used for solving the problem. Section 4 shows the results of numerical simulations. Final considerations are reported in Section 5.

2. Formulation of the Planning Problem

In what follows we refer to an unbalanced, three-phase LV network, connected to the medium voltage (MV) distribution grid through an MV/LV transformer. The LV network includes single-phase and three-phase loads and photovoltaic DG units connected through power converters. It is supposed that all the microgrid’s distributed resources are owned by a unique entity or belong to a consortium of different owners. A microgrid operator manages the grid to provide benefits to the single or clustered owners. The resources are supposed to be coordinated by a centralized control system. Also, it is assumed that the DG units can be controlled by the network operator in terms of reactive power.

In this framework, the proposed procedure is aimed at allocating both single-phase and three phase DESSs in the LV network. The DESSs considered in this paper are batteries connected to the grid through interfacing power converters that can be controlled in terms of active and reactive powers.

The planning of DESSs is formulated as a mixed integer, non-linear, constrained, minimization problem [30]:

$$\min f_{obj}(x,u),$$

(1)

s.t.

$${\psi }_{\mathsf{\nu }}(x,u)=0\mathsf{\nu }=1,\dots ,n_{eq},$$

(2)

$${\mathsf{\varphi }}_{\mathsf{\gamma }}(x,u)\le 0\mathsf{\gamma }=1,\dots ,{\mathrm{n}}_{\mathrm{ineq}},$$

(3)

where ${\psi }_{\mathsf{\nu }}$ and ${\mathsf{\varphi }}_{\mathsf{\gamma }}$ are the νth equality and γth inequality constraints, respectively, $x$ is the vector of state variables (e.g., magnitude and argument of phase voltages at all buses), and $u$ is the DESS unit vector which includes the number of single batteries installed at each phase of each bus of the system:

$$u=\left[n_{DES{S}_1,}^1{n}_{DES{S}_1}^2,\dots ,n_{DES{S}_i}^p,\dots ,n_{DES{S}_{ng}}^3\right],$$

(4)

being $n_{DES{S}_i}^p$ the number of battery units connected at bus i with phase p. If the numbers of battery units connected at the phases of the bus i are all greater than zero (i.e., $n_{DES{S}_i}^1>0,n_{DES{S}_i}^2>0,n_{DES{S}_i}^3>0$), then a three-phase converter will connect the battery to the distribution system. On the contrary, if any of $n_{DES{S}_i}^p$ is equal to zero, one or two single-phase DESSs will be installed. Therefore, both the allocation and typologies of DESSs (single-phase or three-phase) are included in (4).

The planning horizon of the optimization problem refers to a specified number of years ($ny$). Each year y is characterized by $n{d}_y$ typical days that identically apply $N_{y,d}$ times in the year. Each typical day is divided into $nt$ time intervals of duration $\Delta t$. Growth of both load demand and power productions are assumed in the planning time horizon.

2.1. Objective Function

The objective function in (1) to be minimized is the total cost sustained for the installation of DESSs and for the operation of the LV network during the planning period:

$$f_{obj}(x,u)=C_{op}(x,u)+C_{inst}(x,u),$$

(5)

where $C_{op}(x,u)$ is the cost for the purchase of energy from the upstream grid over the planning period and $C_{inst}(x,u)$ includes the cost for the installation and replacement of the DESS. Other cost items such as maintenance cost and replacement cost can be easily included within these two terms [42].

The operation cost can be formulated as follows

$$C_{op}(x,u)={\displaystyle \sum _{y=0}^{ny-1}}{(\frac{1+{\alpha }_L}{1+a})}^y{\displaystyle \sum }_{d=1}^{n{d}_y}{N}_{y,d}\left[{\displaystyle \sum }_{k=1}^{nt}({\displaystyle \sum }_{p=1}^3{P}_{1_{k_{(y,d)}}}^p)\Delta tP{r}_{k_{(y,d)}}\right].$$

(6)

This cost item depends on the price of energy and on the strategy adopted by the grid operator for the management of the connected DESSs including the way they are operated. Furthermore, this cost item refers to the case of the LV network that is not allowed to sell energy to the upstream grid. In the case the LV network can sell energy to the upstream grid, specific feed-in tariffs should be considered.

The second term of the objective function is expressed as

$$C_{inst}(x,u)={\displaystyle \sum }_{i=1}^{ng}\left[(I{C}_{DESS}+r_{i}R{C}_{DESS})E_{DES{S}_i}^{size}\right],$$

(7)

which depends on the size of the installed DESSs and on the need of battery replacement during the planning period. The latter is considered in the above equation by means of $r_i$ which assumes value 0 or 1 depending on the requirement or not to replace the battery. It has to be noted that number of charging/discharging cycles and depth of discharge are factors influencing the need or not to replace the battery, as it will be shown later.

2.2. Installation Constraints

The optimization problem involves constraints related to both the installation of the devices and the operation of both the systems and the LV distribution network. The formers are specifically related to the planning problem and deal with the number of DESSs as well as their size and are reported in

Table 1. Constraints on the Installation of distributed energy storage systems (DESSs).

Number of DESSs
$0\le n_{DES{S}_i}^p\le n_{max}i=1,\dots ,n_{g}p=1,2,3$	(a)
${{\displaystyle \sum }}_{i=1}^{ng}{\displaystyle \sum _{p=1}^3}{n}_{DES{S}_i}^p\le N_{max}$	(b)
$E_{DES{S}_i}^{size}=E_{base}^{size}({{\displaystyle \sum }}_{p=1}^3{n}_{DES{S}_i}^p)i=1,\dots ,n_g$	(c)

Equation (a) imposes a maximum number of DESSs that can be allocated at each phase of each bus whereas Equation (b) imposes a maximum value for the number of DESS elements to be connected into the grid. The total sizes of DESSs allocated at bus i is given by Equation (c) in Table 1 which is the size of the DESS when it is three-phase. In case of single phase DESS placed at phase ${\xi }_i$ of bus i, the battery size $E_{base,i}^{size,{\xi }_i}$ is given by $E_{DES{S}_i}^{size,{\xi }_i}=n_{DES{S}_i}^{{\xi }_i}{E}_{base}^{size}$.

2.3. Operation Constraints

Operation constraints refer to the DESSs, DG units and buses of the network. For ease of reading, only one phase DESS has been considered in the following constraints. In case of two single-phase DESSs, the constraints should refer to both the involved phases. In

Table 2. Constraints on the Operation of DESSs.

Net Energy Charged and Discharged
${\displaystyle \sum }_{k=1}^{nt}\Delta t{\gamma }_{i,k_{(y,d)}}^{{\xi }_i}{P}_{DES{S}_{i,k_{(y,d)}}}^{{\xi }_i}=0iϵ{\mathsf{\Omega }}_{DESS},k=1,\dots ,nt,$ where ${\gamma }_{i,k_{(y,d)}}^{{\xi }_i}=\{\begin{array}{cc}\frac{1}{{\eta }_{dch,i}^{{\xi }_i}}& P_{DES{S}_{i,k_{(y,d)}}}^{{\xi }_i}\ge 0,\\ {\eta }_{ch,i}^{{\xi }_i}& P_{DES{S}_{i,k_{(y,d)}}}^{{\xi }_i}<0.\end{array}$ y = 1, …, ny; d = 1, …, ndy.	(a)
Charging/discharging periods of the day
$\begin{array}{ccc}-P_{DES{S}_{i,max}}^{{\xi }_i,ch}\le P_{DES{S}_{i,k_{(y,d)}}}^{{\xi }_i}\le 0& iϵ{\mathsf{\Omega }}_{DESS}& k\in {\mathsf{\Omega }}_{ch,i_{(y,d)}}^{{\xi }_i},\\ 0\le P_{DES{S}_{i,k_{(y,d)}}}^{{\xi }_i}\le P_{DES{S}_{i,max}}^{{\xi }_i,dch}& iϵ{\mathsf{\Omega }}_{DESS}& k\in {\mathsf{\Omega }}_{dch,i_{(y,d)}}^{{\xi }_i}.\end{array}$ y = 1, …, ny; d = 1, …, ndy.	(b)
Energy stored
$E_{DES{S}_i}^{\mathrm{lb},{\xi }_i}\le E_{DES{S}_{i_{(y,d)}}}^{\mathrm{in},{\xi }_i}-{\displaystyle \sum }_{k=1}^{\tau }\Delta t{\gamma }_{i,k_{(y,d)}}^{{\xi }_i}{P}_{DES{S}_{i,k_{(y,d)}}}^{{\xi }_i}\le E_{DES{S}_i}^{\mathit{ub},{\xi }_i}$, $iϵ{\Omega }_{DESS},$ y = 1, …, ny; d = 1, …, ndy; $\tau =1,\dots ,nt.$	(c)
Active and reactive power
$\sqrt{{(P_{DES{S}_{i,k_{(y,d)}}}^{{\xi }_i})}^2+{(Q_{DES{S}_{i,k_{(y,d)}}}^{{\xi }_i})}^2}\le S_{DES{S}_i}^{{\xi }_i}iϵ{\Omega }_{DESS}$, k = 1, …, nt; y = 1, …, ny; d = 1, …, ndy.	(d)

, with reference to the DESS connected to the bus i ($iϵ{\mathsf{\Omega }}_{DESS}$) and phase ${\xi }_i$, Equation (a) in Table 2 imposes that the net energy charged and discharged during each typical day (d = 1, …, $n{d}_y$) has to be zero.

Equation (b) in Table 2 imposes specified time periods of the day in which each DESS can only charge (${\Omega }_{ch,i_{(y,d)}}^{{\xi }_i}$) or discharge (${\Omega }_{dch,i_{(y,d)}}^{{\xi }_i}$). This allows determining the number of daily charging/discharging cycles which have to be opportunely selected to preserve the lifetime of the battery and according to the energy tariff. Battery’s lifetime, in fact, is related to charging/discharging cycles by the relation $L_b=N_{cycles}/(365\nu )$ where $L_b$ is the battery’s lifetime expressed in years, $N_{cycles}$ is the total number of charging/discharging cycles corresponding to a specific maximum depth of discharge (DoD) (The relationship between $N_{cycles}$ and DoD is provided by the battery’s manufacturer and depends on the battery’s chemistry.); ν is the number of daily charging/discharging cycles. The operational strategy chosen for the battery defines the value of ν.

Here, it is assumed that the considered energy price refers to a Time of Use (ToU) tariff, thus the above time intervals are easily defined by assuming the charging during the off-peak price and discharging during the on-peak price. Alternatively, the charging/discharging time periods could be handled as optimization variables depending on the adopted strategy.

Maximum admissible ranges are imposed to the energy stored in each battery by means of Equation (c) in Table 2 where $E_{DES{S}_i}^{lb,{\xi }_i}$ is related to the maximum admissible DoD, (whose value is still imposed to preserve the battery lifetime) and $E_{DES{S}_i}^{ub,{\xi }_i}$ is the battery size which is given by Equation (c) in Table 1. The aging effect on the capacity of the battery can be taken into account by imposing a decrease of the term $E_{DES{S}_i}^{\mathit{ub},{\xi }_i}$ along the years. The percentage of decrease depends on the aging mechanisms of the specific battery technology [43]. Since DESSs are connected through AC/DC interfacing converters, the size of these last poses limits on the active and reactive powers of the DESSs which are imposed by means of Equation (d) in Table 2. It has to be noted that, in case of three phase DESSs, Equations in Table 2 have to be modified to consider the powers resulting from all of the three phases and the whole three phase energy.

Constraints on DG units (

Table 3. Constraints on distributed generation (DG) units.

Reactive Power Provided by the DG Unit
$\sqrt{{(P_{D{G}_{i,k_{(y,d)}}}^{{\xi }_i,sp})}^2+{(Q_{D{G}_{i,k_{(y,d)}}}^{{\xi }_i})}^2}\le S_{D{G}_i}^{{\xi }_i}$ $iϵ{\Omega }_{DG}$ k = 1, …, nt; y = 1, …, ny; d = 1, …, ndy.

) refer to limits imposed to the provided reactive power which is related to the size of the interfacing converter. Equation in Table 3 refers to a single-phase DG unit connected to phase ${\xi }_i$ of bus i ($iϵ{\mathsf{\Omega }}_{DG}$). It is trivial to derive the case of three phase DG units.

The constraints related to the whole grid buses are reported in

Table 4. Constraints on the network buses.

Load Flow Equations
$P_{i,k_{(y,d)}}^p=V_{i,k_{(y,d)}}^p{\displaystyle \sum }_{j=1}^{ng}{\displaystyle \sum }_{m=1}^3{V}_{j,k_{(y,d)}}^m\left[G_{ij}^{pm}\cos ({\vartheta }_{i,k_{(y,d)}}^p-{\vartheta }_{j,k_{(y,d)}}^m)+B_{ij}^{pm}\sin ({\vartheta }_{i,k_{(y,d)}}^p-{\vartheta }_{j,k_{(y,d)}}^m)\right]$ $Q_{i,k_{(y,d)}}^p=V_{i,k_{(y,d)}}^p{\displaystyle \sum }_{j=1}^{ng}{\displaystyle \sum }_{m=1}^3{V}_{j,k_{(y,d)}}^m\left[G_{ij}^{pm}\sin ({\vartheta }_{i,k_{(y,d)}}^p-{\vartheta }_{j,k_{(y,d)}}^m)-B_{ij}^{pm}\cos ({\vartheta }_{i,k_{(y,d)}}^p-{\vartheta }_{j,k_{(y,d)}}^m)\right],$	(a)
$P_{i,k_{(y,d)}}^p=P_{D{G}_{i,k_{(y,d)}}}^{p,sp}+P_{DES{S}_{i,k_{(y,d)}}}^p+P_{L_{i,k_{(y,d)}}}^{p,sp},$	(b)
$Q_{i,k_{(y,d)}}^p=Q_{D{G}_{i,k_{(y,d)}}}^p+Q_{DES{S}_{i,k_{(y,d)}}}^p+Q_{L_{i,k_{(y,d)}}}^{p,sp},$	(c)
i = 1, …, ng; p = 1, 2, 3; k = 1, …, nt; y = 1, …, ny; d = 1, …, ndy.
Inter-connection bus constraints
$V_{1,k_{(y,d)}}^p=V^{slack},$	(d)
${\vartheta }_{1,k_{(y,d)}}^p=\frac{2}{3}\mathsf{\pi }(1-p),$	(e)
p = 1, 2, 3; k = 1, …, nt; y = 1, …, ny; d = 1, …, ndy,
$\sqrt{{\left({\displaystyle {\displaystyle \sum }_{p=1}^3}{P}_{1,k_{(y,d)}}^p\right)}^2+{\left({\displaystyle {\displaystyle \sum }_{p=1}^3}{Q}_{1,k_{(y,d)}}^p\right)}^2}\le S_{MV},$	(f)
k = 1, …, nt; y = 1, …, ny; d = 1, …, ndy.
All network busbars
$V_{min}\le V_{i,k_{(y,d)}}^p\le V_{max},$	(g)
$k{d}_{i,k_{(y,d)}}\le k{d}_{max},$	(h)
$I_{k_{(y,d)}}^l\le I_z^l$,	(i)
i = 1, …, ng; p = 1, 2, 3; k = 1, …, nt; y = 1, …, ny; d = 1, …, ndy; l = 1, …, nl.

. Particularly, the three phase load flow equations are formulated for each phase p of each bus i which in Equation (a) in Table 4 are expressed with reference to each time interval of each typical day. Also, the active $(P_{i,k_{(y,d)}}^p)$ and reactive $(Q_{i,k_{(y,d)}}^p)$ powers at the left side of Equation (a) in Table 4 refer to injected powers. Their value must comply with Equations (b) and (c) in Table 4 where, in case of absence of load, DG or DESS at the phase p of bus i, the relative power terms at the right side assume value zero. In case of linear loads, they could also be modeled by including the specified load shunt admittance in the network’s admittance matrix.

The constraints related to the whole grid buses are reported in Table 4. Particularly, the three phase load flow equations are formulated for each phase p of each bus i which in Equation (a) in Table 4 are expressed with reference to each time interval of each typical day. Also, the active $(P_{i,k_{(y,d)}}^p)$and reactive $(Q_{i,k_{(y,d)}}^p)$ powers at the left side of Equation (a) in Table 4 refer to injected powers. Their value must comply with Equations (b) and (c) in Table 4 where, in case of absence of load, DG or DESS at the phase p of bus i, the relative power terms at the right side assume value zero. In case of linear loads, they could also be modeled by including the specified load shunt admittance in the network’s admittance matrix.

Equations (d) and (e) in Table 4 refer to the bus no.1, which is assumed as slack bus. This bus is also the point of common coupling thus a constraint is imposed here on the active and reactive powers flowing through the interfacing transformer whose values are limited by its size Equation (f) in Table 4

Finally, with reference to all the grid’s buses, admissible ranges are imposed on the phase-voltage magnitudes in Equation (g) in Table 4, the unbalance factor in Equation (h) in Table 4, and the line phase currents in Equation (i) in Table 4. The unbalance factors in Equation (h) in Table 4 and currents in Equation (i) in Table 4 are expressed as function of the phase voltages.

3. Solution Method

The problem of DESS siting and sizing shown in Section 3 is a large-scale, mixed, non-linear, constrained, optimization problem that can be solved by a GA algorithm. However, the computational efforts required can be excessive and increase significantly with the size of the distribution systems, especially when unbalances and load patterns are taken into account, as they are in this paper. To reduce the processing time while maintaining reasonable accuracy, in this paper we propose to reduce the GA algorithm search space.

For this reason, the allocation problem algorithm, whose flow chart is shown in Figure 1, is dealt with in two successive steps:

• Step 1. application of suitable techniques to identify a reduced feasible region for DESS siting (allocation bus set, ABS);
• Step 2. application of the GA to find the optimal solution in the ABS.

With reference to the reduction of the feasible region (Step 1, Figure 1), it has to be noted that DESS sizing significantly influences both the objective function and the satisfaction of constraints whereas DESS siting significantly influences the satisfaction of constraints but has only a slight influence on the objective function (basically, in terms of grid losses). Thus, the reduction of the feasible region has to be searched with reference only to the DESS siting in order to avoid affecting too negatively the problem solution by acting also on the DESS sizing. More specifically, the reduction of the feasible region is obtained by selecting, among all the network buses suitable for the DESS siting (i.e., theoretically all the network buses), a subset of buses that includes only the buses whose injection has the expected better favorable impact on the network constraints satisfaction (in Equations in Table 4.

Considering that the network operating constraints to be met include the bus voltages in Equation (g) in Table 4, the voltage unbalance factors in Equation (h) in Table 4 and the line currents in Equations in Table 4, in what follows a procedure is set up to obtain the searched ABS. Each DESS should be sited so that it will (potentially) provide the maximum contribution to the improvement of one of the considered constraints. Specifically, initially three different procedures are taken into account to individuate the candidate buses with respect to the impact of the storage system on (i) voltage profile, (ii) unbalances and (iii) line currents, respectively. Among several techniques, the Inherent Structure Theory of Networks (ISTN) [32,33,34,35,36,37,38,39,40,41] is applied for voltage profile and unbalances while the Loading Constraints Criterion (LCC) shown in [40] is applied for lines currents. Then, a criterion to match the candidate buses of all the three procedures is identified to find the ABS.

In the following subsections, the three different techniques are described and the criterion to find the ABS is illustrated firstly (Step 1 above), then, how to apply the GA to find the optimal solution in the ABS is shown (Step 2 above).

3.1. Selection of Candidate Buses for the Improvement of Voltage Profile

Let us consider initially the allocation of one single-phase DESS that will cause the variation of the phase voltages at all buses of the network except the slack bus ($\overline{\Delta \mathrm{V}})$. Having in mind that the network has been represented with a three-phase model, this variation can be evaluated by:

$$\Vert \overline{\Delta V}\Vert =\Vert {(\dot{Y})}^{-1}\overline{\Delta I}\Vert ,$$

(8)

where $\dot{Y}$ is the nodal admittance matrix, $\overline{\Delta I}$ is the vector of the variation of the injected currents and the symbol $\Vert ·\Vert$ indicates the 2-norm of the vector. In this case, the vector $\overline{\Delta I}$ has only one non-zero element, that is the one corresponding to the phase of the bus where the DESS is connected. The non-zero element is the current injected by the single-phase DESS.

Based on the ISTN applied to unbalanced systems [36], if the eigensystem of the admittance matrix is known and assuming that the eigenvalues are all distinct and non-zero and an adequate normalization of eigenvalues is applied, (8) can be exploited as:

$$\Vert \overline{\Delta V}\Vert =\Vert {\displaystyle \sum }_{\mathrm{i}=1}^{3(\mathrm{ng}-1)}\frac{1}{\dot{{\mathsf{\lambda }}_{\mathrm{i}}}}{\dot{X}}_i{\dot{\mathrm{\Gamma }}}_i^{\prime }\overline{\Delta I}\Vert ,$$

(9)

where ${\dot{\mathrm{\Gamma }}}_i^{\prime }$ (${\dot{X}}_i$) is the left (right) eigenvector associated to the ith eigenvalue of $\dot{Y}$ ($\dot{{\mathsf{\lambda }}_{\mathrm{i}}}$).

As it is well known [32,33,34], the decomposition in (9) can be simplified as:

$$\Vert \overline{\Delta V}\Vert \approx \Vert \frac{1}{\dot{{\mathsf{\lambda }}_{\mathrm{k}}}}{\dot{X}}_k{\dot{\mathrm{\Gamma }}}_k^{\prime }\overline{\Delta I}\Vert ,$$

(10)

where $\dot{{\mathsf{\lambda }}_{\mathrm{k}}}$ is the eigenvalue of minimum modulus of $\dot{\mathrm{Y}}$. When the normalization of the eigenvector is used to have unitary norm, (10) can be further simplified as:

$$\Vert \overline{\Delta V}\Vert \approx \frac{1}{{\mathsf{\lambda }}_{\mathrm{k}}}\left|{\dot{\mathrm{\Gamma }}}_k^{\prime }\overline{\Delta I}\right|.$$

(11)

The DESS should be allocated at a bus such that the maximum value of the voltage variation is obtained. Then, if we refer to (11), we can expect that the maximum value of $\Vert \overline{\Delta V}\Vert$ will be experimented if the single-phase DESS is allocated at the bus corresponding to the maximum element (in absolute value) of ${\dot{\mathrm{\Gamma }}}_k$.

Since the selection is based on a simplified equation, that is (10), and it could be necessary to adopt the procedure to allocate more than one DESS (not only single-phase, but also three-phase systems), to determine the best allocations it is advisable to select a subset of candidate buses (instead of a single bus) corresponding to the highest elements of the left eigenvector ${\dot{\mathrm{\Gamma }}}_k$. Hereinafter, this subset will be indicated as ${\mathsf{\Omega }}_1$. The dimension of the subset ${\mathsf{\Omega }}_1$ depends on the number of buses included in it. The selection of the number of candidate buses can be performed by choosing the buses corresponding to the elements i (in absolute value) of the left eigenvector ${\dot{\mathrm{\Gamma }}}_k$ that meet the condition:

$$\frac{{\mathrm{\Gamma }}_{\mathrm{k},\mathrm{MAX}}-{\mathrm{\Gamma }}_{k,i}}{{\mathrm{\Gamma }}_{\mathrm{k},\mathrm{MAX}}}100\le {\mathrm{m}}_{100},$$

(12)

where ${\mathrm{\Gamma }}_{\mathrm{k},\mathrm{MAX}}$ is the maximum element (in absolute value) of ${\dot{\mathrm{\Gamma }}}_k$ and ${\mathrm{m}}_{100}$ is a proper percentage threshold value.

3.2. Selection of Candidate Buses for the Reduction of Unbalances

As it is well-known, the unbalance factor is linked to the inverse-sequence voltage values. Therefore, to obtain the greatest reduction of the unbalance factor, the DESS has to be allocated in such a way that the inverse-sequence voltages at all buses will experience the greatest reduction.

Let us now suppose that, starting from a reference operating point, a modification of the current vector is made. The variation of the inverse-sequence voltages $\Vert {\overline{\Delta V}}^{inv}\Vert$ (at all buses except the slack bus) can be evaluated as (see Appendix A):

$$\Vert {\overline{\Delta V}}^{inv}\Vert =\Vert {\dot{Z}}^{inv-zero}{\overline{\Delta I}}^{zero}+{\dot{Z}}^{inv-dir}{\overline{\Delta I}}^{dir}+{\dot{Z}}^{inv-inv}{\overline{\Delta I}}^{inv}\Vert ,$$

(13)

where ${\overline{\Delta I}}^{zero}$, ${\overline{\Delta I}}^{dir}$ and ${\overline{\Delta I}}^{inv}$ are the vectors of the variation of the injected currents at zero, direct, and inverse sequence and ${\dot{Z}}^{inv-zero}$, ${\dot{Z}}^{inv-dir}$, ${\dot{Z}}^{inv-inv}$ are the coupling impedance sub-matrices between the inverse and zero sequences, the inverse and direct sequences and the inverse and the inverse sequences, respectively (the definitions are reported in Appendix A).

As it is well known [44], it can be stated that:

$$\begin{array}{l}\Vert {\dot{Z}}^{inv-zero}{\overline{\Delta I}}^{zero}+{\dot{Z}}^{inv-dir}{\overline{\Delta I}}^{dir}+{\dot{Z}}^{inv-inv}{\overline{\Delta I}}^{inv}\Vert \\ \le \Vert {\dot{Z}}^{inv-zero}{\overline{\Delta I}}^{zero}\Vert +\Vert {\dot{Z}}^{inv-dir}{\overline{\Delta I}}^{dir}\Vert +\Vert {\dot{Z}}^{inv-inv}{\overline{\Delta I}}^{inv}\Vert .\end{array}$$

(14)

If only inverse-sequence currents are supposed to change, the variation of the inverse-sequence voltages (13) is simplified as:

$$\Vert {\overline{\Delta V}}^{inv}\Vert =\Vert {\dot{Z}}^{inv-inv}{\overline{\Delta I}}^{inv}\Vert .$$

(15)

Assuming once again that the eigenvalues of $inv({\dot{Z}}^{inv-inv})$ are all distinct and non-zero, the relation above can be exploited as:

$$\Vert {\overline{\Delta V}}^{inv}\Vert =\Vert {\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{ng}-1}\frac{1}{{\dot{\mathsf{\lambda }}}_{\mathrm{i}}^{\mathrm{inv}-\mathrm{inv}}}{\dot{X}}_i^{inv-inv}{({\dot{\mathrm{\Gamma }}}_i^{inv-inv})}^{\prime }{\overline{\Delta I}}^{inv}\Vert ,$$

(16)

with the same assumptions and meaning of the quantities as in (9), but in this case the eigensystem refers to the matrix ${({\dot{Z}}^{inv-inv})}^{-1}$.

Let ${\dot{\mathsf{\lambda }}}_{\mathrm{k}}^{\mathrm{inv}-\mathrm{inv}}$ be the eigenvalue with minimum modulus; the 2-norm of the inverse-sequence voltage variation vector (16) can be approximated as:

$$\Vert {\overline{\Delta V}}^{inv}\Vert \approx \frac{1}{{\mathsf{\lambda }}_{\mathrm{k}}^{\mathrm{inv}-\mathrm{inv}}}|{({\dot{\mathrm{\Gamma }}}_k^{inv-inv})}^{\prime }{\overline{\Delta I}}^{inv}|.$$

(17)

In (17), it has been assumed, once again, that the 2-norm of the vector ${\dot{X}}_k^{inv-inv}$ is unitary. Similar to the case of voltage profile (Section 4.1), to obtain the maximum variation of $\Vert {\overline{\Delta V}}^{inv}\Vert$ with the same injection of the inverse-sequence current, the buses of the inverse-sequence current injections have to be searched in the set of buses corresponding to the highest values (in absolute value) of elements of ${\dot{\mathrm{\Gamma }}}_k^{inv-inv}$.

Based on (13), if direct-sequence, inverse-sequence and zero-sequence currents change simultaneously, to obtain the maximum variation of $\Vert {\overline{\Delta V}}^{inv}\Vert$, the analysis above reported should be performed taking into account (13) and, then, with respect to the spectral decomposition of all coupling sub-matrices${({\dot{Z}}^{inv-zero})}^{-1}$, ${({\dot{Z}}^{inv-dir})}^{-1}$, ${({\dot{Z}}^{inv-inv})}^{-1}$, and searching in the set of busses associated with the highest values not only of elements of ${\dot{\mathrm{\Gamma }}}_k^{inv-inv}$ but also of elements of ${\dot{\mathrm{\Gamma }}}_k^{inv-dir}$ and ${\dot{\mathrm{\Gamma }}}_k^{inv-zero}$, with obvious meaning of the symbols. In what follows, the three subsets are indicated as ${\mathsf{\Omega }}_2,{\mathsf{\Omega }}_3$ and ${\mathsf{\Omega }}_4$. Also in this case, to determine the subsets ${\mathsf{\Omega }}_2,{\mathsf{\Omega }}_3$ and ${\mathsf{\Omega }}_4$ the criterion expressed by (12) can be applied and their dimensions will depend on the choice of the percentage threshold value.

3.3. Selection of Candidate Buses for the Reduction of Line Currents

The identification of candidate buses for the reduction of line currents can be carried out following the idea of load areas first introduced in [40] as part of the ADDRESS project [41]. The buses can be identified in two steps. In the first step, either historical data of measurements or results of load flow analysis can be used to recognize the most significant critical events on the network under study in terms of current flowing in the network branches. Each critical event can be identified considering either an overload event (the line current is greater than the line ampacity) or a condition meeting the following relationship (for the generic line l):

$$\frac{{\mathrm{I}}_{\mathrm{z}}^{\mathrm{l}}-{\mathrm{I}}_{\max }^{\mathrm{l}}}{{\mathrm{I}}_{\mathrm{z}}^{\mathrm{l}}}100\le {\mathrm{n}}_{100}$$

(18)

where ${\mathrm{n}}_{100}$ is a proper percentage threshold value.

In the second step, for the most common case of radial passive networks, for each critical event, a list of all the buses downstream the component involved is identified by means of simple graph navigating techniques. Obviously, in active distribution systems with DG units that can reverse the current flows in some lines under some circumstances, we can include in the list all the buses downstream the direction of the current flow. In what follows, the subset of buses identified as above is indicated as ${\mathsf{\Omega }}_5$.

3.4. Overall Selection of Allocation Bus Set for the Distributed Energy Storage Systems Siting

Once defined the criteria for the selection of the candidate buses with respect to the impact of the DESS on voltage profile, unbalances and line currents, it is needed to establish (i) the values of the thresholds such as ${\mathrm{m}}_{100}$ and ${\mathrm{n}}_{100}$ that define the dimensions of the subsets ${\mathsf{\Omega }}_1\dots {\mathsf{\Omega }}_5$ and (ii) a criterion to match all the subset ${\mathsf{\Omega }}_1\dots {\mathsf{\Omega }}_5$ to obtain the final list of candidate buses to be included in the ABS for the DESS siting.

The choice of the value of the thresholds such as ${\mathrm{m}}_{100}$ and ${\mathrm{n}}_{100}$ and the selection of the matching criterion express a trade-off between the need of a high precision analysis (the higher are the thresholds, the larger the final set of ABS becomes) and the reduction of computational burden (the higher are the thresholds, the higher the computational efforts becomes). In our experience, threshold values lower than or equal to 10% seem to be a suitable compromise solution. With reference to the matching criterion, the final ABS ${\mathsf{\Omega }}_{\mathrm{ABS}}$ of candidate buses was obtained as:

$${\mathsf{\Omega }}_{\mathrm{ABS}}={\mathsf{\Omega }}_1{{\displaystyle \cup }}^{\text{}}{\mathsf{\Omega }}_2{{\displaystyle \cup }}^{\text{}}{\mathsf{\Omega }}_3{{\displaystyle \cup }}^{\text{}}{\mathsf{\Omega }}_4{{\displaystyle \cup }}^{\text{}}{\mathsf{\Omega }}_5.$$

(19)

3.5. Hybrid Genetic Algorithm—Sequential Quadratic Programming Algorithm

Once the ABS has been identified, to solve the optimization problem a hybrid approach was used, including a GA and an inner algorithm based on sequential quadratic programming (SQP) optimization [19]. As shown in the flow chart of Figure 1, the GA creates an initial population of individuals characterized by specific buses where siting DESSs of specified sizes, according to the constraints of Table 1.

Once generated the initial population, the inner algorithm solves an optimization problem whose optimization variables are continuous. More specifically, with reference to each typical day of the whole planning period, the solution of the inner algorithm provides the daily curve of power required from the upstream grid, of the reactive power of the DG units, of the active and reactive powers of DESSs, each corresponding to the minimization of the objective function (6) while satisfying the constraints of Table 2, Table 3 and Table 4.

For each individual in the population, the fitness function (5) is calculated by using the outputs resulting from the inner optimization. Regarding the stopping criterion of the GA, it ends when the best fitness function value remains constant over an assigned number of populations’ generations. In order to impose a limitation on the computational time, a maximum number of individuals to be investigated can also be set.

4. Numerical Application

The procedure described above was applied to a LV distribution network which is a Cigrè benchmark system whose scheme is shown in Figure 2 and whose details are reported in [45]. The grid’s nominal voltage is 400 V, it includes 41 busses and is connected to a 20 kV MV grid by means of three transformers of 500, 150 and 300-kVA.

The three transformers supply a residential, a commercial and an industrial feeder, respectively. The price of energy, which is based on an existing tariff, is reported in

Table 5. Time of Use (ToU) Tariff.

	Summer Tariff		Winter Tariff
	Period	Price ($/MWh)	Period	Price ($/MWh)
On peak	12:00–18:00	542.04	8:30–21:30	161.96
Part Peak	8:30–12:00 18:00–21:30	252.90
Off Peak	21:30–8:30	142.54	21:30–8:30	132.54

[46]. DG units are supposed to be connected to the grid. They are four single-phase 5 kW PV units connected to the buses no.15 (phase no.2), no.19 (phase no.1), no.37 (phase no.3), no.40 (phase no.3) and one three phase 20 kW PV unit located at bus no.21.

It is supposed that both single-phase and three-phase DESSs can be allocated whose base size is 4 kWh having a nominal discharging time of 5 h. DESSs are based on Li-Ion batteries technology having an efficiency of 90% in charging mode and 92% in discharging mode. The batteries are assumed having a maximum DoD of 80%. To maximize their lifetime, the batteries can have only one charging/discharging cycle per day. More specifically, they can charge during the off-peak hours indicated in Table 5 and discharge during the rest of the day.

The battery’s installation cost is supposed 600$/kWh whereas the replacement cost is assumed 250$/kWh. Two typical days were considered for the cost evaluation, one winter day and one summer day. Both effective rate of change and discount rate are assumed equal to 3%. An annual 2% growth is assumed for the load demand. As an example, Figure 3 shows the per unit value profiles of three loads of the residential, industrial and commercial feeder, respectively. The locations of the loads and their peak power are reported in

Table 6. Peak Active Power and location of loads.

Bus	Phase	Active Power [kW]	Bus	Phase	Active Power [kW]
Residential Feeder			Commercial Feeder
12	1	2.4	33	1	0.6
	2	1.9		2	0.6
	3	5.2		3	3.6
16	1	17.4	34	1	1.2
	2	14		2	1.2
	3	38		3	7.2
17	1	8.7	35	1	1.5
	2	7		2	1.5
	3	19.2		3	9.0
18	1	5.5	38	1	0.72
	2	4.4		2	0.72
	3	12.2		3	4.3
19	1	7.4	39	1	1.5
	2	6		2	1.5
	3	16.4		3	9.0
Industrial feeder			40	1	0.48
21	1	14.2		2	0.48
	2	11.3		3	2.9
	3	31.2	41	1	0.48
				2	0.48
				3	2.88

. These values were derived from [45] and properly modified to obtain critical conditions in terms of voltages and currents. Constant power factors of 0.95 (residential), 0.85 (industrial) and 0.9 (commercial) have been assumed during the simulations.

In what follows three case studies have been reported referring to:

• a reference case without DESS;
• DESSs are allocated by means of only a GA;
• DESSs are allocated by means of the proposed procedure (GA applied for ABS).

4.1. Case Study 1

This case study is useful to show the operation of the considered network without the support of distributed resources, that is DGs do not provide reactive power support and there is no DESS. In this case, the network is characterized by the presence of unbalances that sometimes exceed the imposed maximum limit (2%). As an example, Figure 4 shows the unbalance factor at each hour of a typical day of the summer season of the last year of the planning period (All of the figures reported hereinafter refer to the same typical day.). In this case it can be observed that the commercial line is characterized by quite high values of the unbalance factor which in two buses, in some hours of the day, even exceeds the limit. A similar behavior is observed also in one bus of the residential line.

Line currents and bus voltages also show in some cases critical conditions. In Figure 5 the values of the line currents are reported in p.u. with reference to hour 21.00 showing a violation of the ampacity limits in some lines of the residential feeder. Figure 6 shows the voltage magnitude of the residential feeder at hour 21.00 of the typical day. Also in this case a violation of the limit can be observed. Critical conditions have been also observed in the other two feeders (here not reported for the sake of brevity) where, however, no violations are encountered.

4.2. Case Study 2

In this case the support of the distributed resources is considered and DESSs are allocated by means of a GA applied over the set constituted by all network’s buses. A maximum value of 440 kWh was imposed to the capacity that could be allocated in the whole network and a maximum value of 12 kWh was imposed to the capacity that could be placed at the single bus. Both single-phase and three-phase DESSs were allocated whose location and capacity are reported in

Table 7. Case Study 2: DESS Allocation results.

Bus no.	Phase No.	DESS Capacity (kWh)	Bus No.	Phase No.	DESS Capacity (kWh)
Residential feeder			16	2	12
4	3	12	16	3	12
5	1	4	17	1,2,3	20
5	3	4	Commercial feeder
6	2	8	23	1,2,3	28
6	3	8	25	3	8
7	3	8	27	3	12
8	3	8	28	1	8
10	2	8	29	1,2,3	36
10	3	8	33	1,2,3	24
11	1	12	34	3	4
11	2	4	36	1,2,3	24
13	3	4	39	1	4
14	1,2,3	24	41	1	12
15	1	4	41	2	12

. Six three-phase DESSs were allocated, two of them in the residential feeder and four in the commercial feeder. Many single-phase devices were also allocated in both these feeders whereas no DESS were allocated in the industrial feeder. The overall DESS’s capacity resulted equal to 332 kWh that allowed a reduction of the total cost of more than 5% compared to case 1.

Moreover, thanks to the support of the resources, the operation of the grid improved in terms of unbalances, line currents and bus voltages. This can be observed in Figure 7, Figure 8 and Figure 9 which show, still referring to the summer typical day of the last year, the hourly unbalance factor at all buses, the line currents and the voltage magnitude (both referring to the residential feeder at hour 21.00), respectively.

Figure 7 clearly shows a general reduction of the unbalance factor in all the network buses and, particularly, in those buses of the commercial feeder where in case 1 a violation of the unbalance factor’s limits appeared. Analogous improvement can be seen with reference to the line currents that, as can be observed in Figure 8, are always within their ampacity limit and to the voltage magnitude that, as can be seen in Figure 9, is always within the limits.

4.3. Case Study 3

In this case the support of the distributed resources is considered and DESSs are allocated by means of the proposed hybrid approach. At this aim, the ISTN was applied to identify the set of candidate buses with reference to voltage deviation and unbalance factor and assuming a threshold limit of 10%. Regarding the line currents, the loading constraints criterion was adopted to identify the candidate buses where the critical events were identified by considering the overload events of line current greater than the line ampacity. The resulting ABSs are shown in

Table 8. Case Study 3: DESS allocation bus set (ABS).

Location Bus Set	Residential Feeder Bus No. (Phase No.)	Industrial Feeder Bus No. (Phase No.)	Commercial Feeder Bus No. (Phase No.)
${\mathsf{\Omega }}_1$ —Voltage deviation	10(2), 11(2), 16(2), 18(2), 19(1,2,3)	21(1,2,3)	29(1,2), 30(1,2,3), 40(1,2), 41(1,2,3)
${\mathsf{\Omega }}_2$—Unbalance factor	10, 11, 18, 19	21	29, 30, 40, 41
${\mathsf{\Omega }}_3$—Unbalance factor	9, 10, 11, 15, 16, 18, 19	21	29, 30, 37, 38, 40, 41
${\mathsf{\Omega }}_4$—Unbalance factor	9, 10, 11, 18, 19	21	29, 30, 40, 41
${\mathsf{\Omega }}_5$—Line currents	3(3), 4(3), 13(3), 14(3), 15(3), 16(3)	-	-

, with reference to each feeder. In

Table 9. Case Study 3: DESS Allocation results.

Bus No.	Phase No.	DESS Capacity (kWh)	Bus No.	Phase No.	DESS Capacity (kWh)
Residential Feeder			Industrial Feeder
9	1,2,3	36	21	1,2,3	20
10	1,2,3	36	Commercial feeder
11	1,2,3	36	29	2	8
15	1,2,3	36	29	3	8
16	1,2,3	28	30	1,2,3	32
18	1	8	37	1,2,3	36
18	2	8	40	3	8
			41	1	8

the results of the proposed approach in terms of location and capacity of the allocated DESSs are reported.

Eight three-phase DESSs were allocated, five of them in the residential feeder, two in the commercial feeder and one in the industrial feeder. Six single-phase devices were also allocated in both residential and commercial feeders whereas no single-phase DESSs were allocated in the industrial feeder. The overall DESS’s capacity resulted equal to 344 kWh that allowed a reduction of the total cost of more than 6% compared to case 1.

By comparing Table 8 and Table 9, it is interesting to note that DESSs are placed in almost all of the candidate buses. Moreover, Table 9, compared to Table 7, shows that the size of the DESSs allocated in case 3 is generally larger than that of the DESSs allocated in case 2. This allowed a significant reduction of the objective function compared to case 2 still preserving the network constraints, as shown in Figure 10, Figure 11 and Figure 12. These figures show the hourly unbalance factor at all buses, the line currents and the bus voltages, respectively.

As an example of further results obtained, Figure 13, Figure 14 and Figure 15 show the active and reactive power of the controllable distributed resources. More specifically, still referring to the summer typical day of the last year, the phase active power daily profiles of the DESS placed at bus no.37 are reported in Figure 13. The corresponding stored energy profile is show in Figure 14. Figure 15 reports the reactive power support provided by the DESS placed at bus no.37 (only phase 1 is shown) and the reactive power supplied by the DG unit located at the bus no.15 (phase 2). As clearly shown in these figures, the constraints imposed on the resources are satisfied. Furthermore, even if the figure is not reported here for the sake of conciseness, the constraints imposed by the rated powers of the interfacing converters of both DESSs and DG units are satisfied.

Finally, aimed at highlighting the advantages that can be drawn by the use of the proposed hybrid ABS and GA solution method, in Figure 16 the best values of the fitness function for all of the populations created and analyzed by the GA (case 2) and the ABS/GA (case 3) are reported versus the number of individuals investigated (particularly, the values are reported in per unit with reference to the value assumed in case 1). By observing the figure, it clearly appears that the ABS/GA solution method allows identifying better solutions in a faster way, that practically means, in a lower number of iterations (i.e., individuals investigated). In fact, the best fitness function values in case of ABS/GA are always lower than those evaluated through the GA, thus meaning that the convergence criteria are met faster with the proposed approach. As stopping criterion, in our case study, we assumed that the procedure ended when the best fitness function maintained a constant value over ten consecutive iterations. With this stopping criterion, as shown in Figure 16, GA/ABS converged after analyzing 260 individuals reaching a best fitness function value of $4.78 M. It is interesting to note that in the case of GA (case 2), after analyzing 260 individuals the best fitness function value is greater (i.e., $4.85 M). Moreover, in case 2, even with higher number of iterations, the GA provides worse solutions. In the case of Figure 16, the GA was stopped after having investigated 500 individuals, without meeting the stopping criterion and without obtaining a better solution, compared to the case 3 since the fitness function best value is still higher (i.e., $4.82 M). As a consequence of that, we can conclude that in the examined case GA/ABS allows identifying in a fast way a good solution with a drastically reduced computational time compared to the case of applying GA. We note also that the time requested to find the solution also decreases because in the case of GA/ABS the inner SPQ solution is requested to solve optimal load flow problems with a reduced number of DESSs of higher sizes, thus reducing the number of the optimization variables and so, requiring a lower computational time.

In order to give an idea of the computational time of the simulations, by using a commercial personal computer (2.93 GHz processor, 4 GB RAM), the simulation time for each individual ranged from seconds to minutes depending on the number and size of DESSs corresponding to the individual. By increasing the number of DESSs, in fact, the optimization variables (i.e., stored energy, active and reactive power at each of the 24 h of the typical days) significantly increase. So high computational time is not discouraging if we consider that we are in a planning stage and that we used a detailed and complex model of the system under study to represent both the unbalances and the behavior of DESSs with the consequence that each GA individual implies many optimization variables. Anyway, as obvious, the use of more sophisticated hardware and software arrangements will allow reducing dramatically the computation time.

As a final consideration we outline that in order to check both the proposed procedure and the usefulness of adopting DESSs to face unbalances and PQ problems, numerous further simulations were performed, not reported here for the sake of conciseness. More specifically, further investigations were made by applying the procedure imposing (i) a higher value to the maximum battery’s capacity that can be installed in the network; (ii) more stringent constraints to the PQ requirements. By analyzing the results, we observed that by allocating a higher number of DESSs further economic benefits can be obtained. In case of more stringent constraints on the PQ requirements, the proposed hybrid GA/ABS procedure still resulted more efficient than the GA in terms of computational time and accuracy of the results. As an example, when reducing the limit on unbalance factor, the set of the candidate buses provided by the ISTN still allowed a faster identification of the location buses able to satisfy the constraints and provided a better objective function value than that obtained with the GA.

5. Conclusions

In this paper a new methodology for the allocation of energy storage devices in LV distribution networks was proposed. The methodology was based on an accurate modeling of both three-phase network and components, including loads, distributed generation units and distributed energy storage devices. In particular, both single-and three-phase configurations were modeled and used to properly formulate the optimization strategy in terms of cost minimization and satisfaction of constraints on the network and on its components. Attention was paid to typical power quality problems such as unbalance factor and voltage limits violation. The complexity of the considered planning problem as well as the need of managing large numbers of storage devices in large distribution networks lead to problems related to the computational burden. Aimed at addressing this issue, the proposed approach was faced in this paper by means of a hybrid solution method based on the reduction of the algorithm search space. The proposed method is able to manage the problem complexity and to identify solutions with a high degree of accuracy while containing the computational burden so making the solution method feasible and effective. The proposed method was also applied to a test LV network typically adopted to study the impact of distributed resources under smart grid context. The results clearly showed the feasibility and effectiveness of the proposed solution in terms of accuracy and reduced computational burden. The procedure was applied to critical scenarios, where the planning problem has to face violation of some technical constraints. Based on the proposed method, our current research is focused on other complexities, such as the uncertainties in loads and generations