A New Hybrid Approach Using the Simultaneous Perturbation Stochastic Approximation Method for the Optimal Allocation of Electrical Energy Storage Systems

: This paper deals with the optimal allocation (siting and sizing) of distributed electrical energy storage systems in unbalanced electrical distribution systems. This problem is formulated as a mixed, non-linear, constrained minimization problem, in which the objective function involves economic factors and constraints address the technical limitations of both network and distributed resources. The problem is cumbersome from the computational point of view due to the presence of both constraints of an intertemporal nature and a great number of state variables. In order to guarantee reasonable accuracy-although limiting the computational efforts-a new approach is proposed in this paper: it is based on a Simultaneous Perturbation Stochastic Approximation (SPSA) method and on an innovative inner algorithm, which allows it to quickly carry out the daily scheduling (charging/discharging) of the electrical energy storage systems. The proposed method is applied to a medium voltage (Institute of Electrical and Electronics Engineers) IEEE unbalanced test network, to demonstrate the effectiveness of the procedure in terms of computational effort while preserving the accuracy of the solution. The obtained results are also compared with the results of a Genetic Algorithm and of an exhaustive procedure.


Introduction
Electrical Energy Storage Systems (EESSs) have been recognized as a viable solution for implementing the smart grid paradigm, providing features in load levelling, integrating renewable and intermittent sources, improving power quality (PQ) and reliability, reducing energy import during peak demand periods, and so on [1].
In particular, EESSs can be exploited in distribution systems to pursue several objectives that range from implementing demand response to minimizing electrical energy costs. System objectives can also be pursued; indeed, voltage profiles as well as other PQ aspects (such as unbalances) can be improved by controlling EESSs smart interfacing converters. For instance [2] proposed the compensation of unbalances by adequately using storage systems.
In this paper, a new approach based on the Simultaneous Perturbation Stochastic Approximation (SPSA) method [21,22] is proposed for optimally allocating (i.e., sizing and siting) EESSs in unbalanced distribution systems.
SPSA is widely used to solve optimization problems in several frameworks [23]. It operates by simultaneously perturbing all of the unknown problems (both continuous and integer), stochastically approximating differentiation at each iteration. Although SPSA is a stochastic optimization solution method, simultaneous population-free perturbations make the computation very effective.
In addition, the complexity effort introduced by adding constraints of an intertemporal nature (due to the EESSs daily charging/discharging stages) is overcome in this paper by applying an inner simplified approach. This approach is based on a busbar system where EESSs are exploited to reduce the energy purchase during hours with greater prices and to increase the energy during hours with smaller prices. Furthermore, the simplified strategy attempts to realize a load leveling: when the EESS is charging or discharging, the load curve is possibly flattened.
Eventually, the main features of the proposed approach are: • A distributed storage is considered to catch the potential advantages brought by EESSs in an unbalanced distribution system. • The procedure accounts for many economic and technical aspects of the EESSs allocation.

•
The implementation of the solving algorithm based on the SPSA method allows to considerably shorten the computational time while providing good-quality solutions.

•
The inner simplified approach allows it to quickly carry out the daily scheduling of the EESSs, further shortening the computational time.

•
The comparison of the obtained results with the results of a Genetic Algorithm (GA) and of an exhaustive procedure gives evidence of the accuracy and of the computational effort reduction.
This paper is organized as follows. Section 2 shows the formulation of the planning problem. Section 3 refers to the new approach used for solving the optimization problem. Section 4 briefly recalls the GA applied for comparison. Section 5 shows the results of numerical simulations and the comparisons with a GA solution and an exhaustive solution. Final considerations are reported in Section 6.

Problem Formulation
The EESSs can be sized and located in an unbalanced distribution system, aiming at optimizing the benefits for the whole system. In particular, several effects of EESSs on the distribution system can be taken into account and can be optimized.
In particular, economic and technical effects have to be considered. In fact, EESSs can modify the pattern of energy imported from the upstream grid in view of an investment for the EESSs installation. Moreover, the allocation (siting and sizing) of EESSs can significantly affect the PQ levels at all buses of the distribution systems. Indeed, an adequate allocation of EESSs can improve the voltage profile at all the system nodes and can control the currents flowing through the system lines. The impact of EESSs on the unbalanced factors is another aspect that is worth being investigated. In this respect, the EESSs can be effectively designed at the planning stage in order to obtain the best performance and to support the distribution system operator (DSO) in keeping unbalance factors under the maximum Standard allowable values.
The planning of EESSs for an unbalanced distribution system can be formulated as a mixed integer, non-linear, constrained, optimization problem, in which a proper objective function is minimized and a large number of equality/inequality constraints are met.
In this paper, the objective function f obj to be minimized depends on the cost of the energy acquired from the upstream grid over the planning period, and on the EESSs costs: where PV C E is the present value of the cost of the energy that the distribution system acquires from the upstream grid in the planning period, and PV(C ES ) is the present value of the EESSs cost. The vector of problem variables X includes the size, the allocation bus, and the allocation phase of each EESS.
The following relationship can be obtained by expanding the cost items in Equation (1): The objective function in Equation (2) accounts also for the battery replacement by means of the term r s y . When, during year y, there is a battery installation (only at y = 1) or replacement, r s y is not zero. If the battery lifetime exceeds the remaining years of the planning period, the value of r s y is less than 1 to account for the residual value of the battery at the end of the planning period. To make a reasonable economic analysis, a trend of the batteries installation costs is included in Equation (2), with an installation cost EC ES y that varies with the considered year y. We assume, moreover, not to replace the interfacing static converters during the planning period.
The decision variables of the optimization problem are the sizing (power and energy) and siting of the single-phase and three-phase storage systems. In particular, the power sizing of the storage system is assumed to be discrete and multiple of an elementary size, whereas the siting is clearly a discrete variable (the grid buses). Therefore, if the numbers of battery units connected at the phases of a generic bus are all greater than zero, then a three-phase converter will connect the battery to the distribution system. On the contrary, if one of these numbers is zero, single-phase EESSs will be installed. The energy size of an EESS; i.e., E size s in Equation (2)-is the available energy capacity. The objective function in Equation (2) has to be minimized subject to meeting a set of equality and inequality constraints, which refer to the technical limitations of both the network and its distributed resources. In the following, for the sake of simplicity, we refer only to three phases of EESS.
First of all, for each storage system and during each day of operation, the energy charged must be equal to the energy discharged: The efficiency γ s,h in Equation (3) depends also on the operation of the EESS; indeed, the values in charging and in discharging steps may be different. The power P ES s,h,d,y is positive or negative, according to the discharge or charge steps.
Due to constraints on the expected life of the batteries, each battery can be charged or discharged in assigned hours. Since the objective is to minimize the cost of electricity, the hours of charging and discharging depend on the structure of the electricity prices (tariffs).
The value of P ES s,h,d,y for each battery cannot exceed admissible ranges. In particular, during the charging hours, the following constraints apply: During the discharging hours, the following constraints apply: The size of the EESS AC/DC interfacing converter imposes constraints on the active and reactive powers that the EESS can absorb/inject. In particular, for each storage system and during each hour of operation of the energy storage, the following constraints have to be verified: Note that the apparent power S ES s coincides with the power size when the converter operates at unitary power factor.
Furthermore, for each EESS, the nominal discharging time (i.e., the ratio between the energy size E size s and the power size P size s ) has to be constrained into a range defined by the specific technology of the storage device.
With reference to the network, the three-phase load flow equations [24] have to be included. In particular, the following equations apply at each three-phase node: where P p b,h,d,y and Q p b,h,d,y are the net injected active and reactive powers. Extending Equations (7) in order to include single-phase and two-phase nodes is trivial.
With reference to the slack bus, that is set at the bus of interconnection to the upstream network (i.e., bus #1), the magnitude and the argument of phase voltages are specified: Moreover, at the interconnection bus, the apparent power flowing through the interfacing transformer is constrained by its rating S trans f ormer ; then, the following constraint has to be considered: Meeting PQ requirements at all of the buses leads to the following constraints: Eventually, the lines of the system have a specified ampacity that cannot be exceeded; then, the line phase currents are limited as well: We note that the planning problem is based on the choice of a planning period, the length (i.e., the number of years) of which depends on the expected lifetimes of the equipment. Moreover, it is evident that considering all of the days of all of the years in the planning period may push the computational effort beyond reasonable time scales; therefore, it is rational to consider a reduced set

Solving Procedure
The planning problem formulated in Section 2 is a cumbersome optimization problem that may require tremendous computational effort to be solved, particularly in unbalanced distribution systems. Therefore, developing new methods that can limit the computational effort in the solution of the optimization problem while saving the accuracy of the results, is quite essential.
Based on these considerations, a hybrid approach is proposed in this paper to quickly and accurately solve the optimization problem presented in Section 2. This approach consists in a four-step iterative procedure that includes an inner routine ( Figure 1): 1st Step: The Simultaneous Perturbation Stochastic Approximation (SPSA) method selects specific buses at which EESSs of specified power ratings are sited. 2nd Step: In the inner routine, an algorithm based on a simplified approach provides, for each day of each year, the hourly active and reactive power profiles of the EESSs and their energy sizes. 3rd Step: The network constraints (Equations (7)-(12)) are verified through load flow analyses. In this way, only feasible conditions are taken into consideration, whereas conditions that do not satisfy the network constraints are discarded. If the network constraints are never verified, the algorithm returns to step 1, otherwise it calculates (Equation (2)). 4rd Step: The convergence is checked by monitoring the number of iterations and/or the decrease of the objective function (Equation (2)). If the convergence is not achieved, the algorithm returns to step 1.
Note that, in the inner routine, a simplified algorithm minimizes the objective function (Equation (2)), satisfying the EESSs constraints. It is based on a busbar system in which EESSs are exploited to reduce the energy purchase during hours with greater prices, and to increase the energy purchase during hours with smaller prices. Furthermore, the implemented strategy attempts to realize a load leveling: when the EESS is charging or discharging, the load curve is possibly flattened.
Note also that, solving the EESSs siting and sizing planning problem by means of non-simplified procedures may require so much time that it is difficult to obtain reasonable results. Instead, the simplified step-procedure proposed in this paper allows it to dramatically shorten the computational time while saving the accuracy of the results. In fact, the algorithm based on the SPSA provides good solutions with a considerably shortened computational time; moreover, the inner simplified approach allows it in turn to quickly carry out the daily scheduling of the EESS, further shortening the computational efforts.
In Sections 3.1 and 3.2, details about the aforesaid algorithms will be provided. Moreover, Section 4 briefly recalls the micro Genetic Algorithms (µGA) applied in the numerical applications of Section 5 to compare the proposed procedure. time while saving the accuracy of the results. In fact, the algorithm based on the SPSA provides good solutions with a considerably shortened computational time; moreover, the inner simplified approach allows it in turn to quickly carry out the daily scheduling of the EESS, further shortening the computational efforts.
In Sections 3.1 and 3.2, details about the aforesaid algorithms will be provided. Moreover, Section 4 briefly recalls the micro Genetic Algorithms (μGA) applied in the numerical applications of Section 5 to compare the proposed procedure.

The Simultaneous Perturbation Stochastic Approximation Method
The SPSA method was firstly proposed by J. C. Spall [21] to solve optimization problems and, since then, it has widely been applied to several planning problems [23].
The SPSA method is based on the simultaneous perturbation of all unknown variables x i (i = 1, . . . , q) and on the differentiation approximation. This is performed in a stochastic way for each iteration [22].
Let us initially refer to a generic unconstrained optimization problem with integer variables: being the vector X composed of q integer variables: Let Energies 2018, 11, x 7 of 19

The Simultaneous Perturbation Stochastic Approximation Method
The SPSA method was firstly proposed by J. C. Spall [21] to solve optimization problems and, since then, it has widely been applied to several planning problems [23].
The SPSA method is based on the simultaneous perturbation of all unknown variables = 1, … , and on the differentiation approximation. This is performed in a stochastic way for each iteration [22]. Let us initially refer to a generic unconstrained optimization problem with integer variables: being the vector composed of q integer variables: Let be the operator that rounds a real number to the nearest integer toward positive infinity.
The elements of the vector are updated at each iteration of the SPSA algorithm; for instance, considering the kth iteration, we have: where: . be the operator that rounds a real number to the nearest integer toward positive infinity. The elements of the vector X are updated at each iteration of the SPSA algorithm; for instance, considering the kth iteration, we have: where: The constant a k is included to accelerate the convergence, and it depends on the expected maximum number of iterations, the expected step size, and the starting values of the optimization variables; ∆ k is the random perturbation vector (with dimension q), the components of which are independently generated following a Bernoulli (±1) distribution [25]; the constant c k is a positive number. More details about the values of the constants and on the generation of the perturbation vector are in [25].
If the optimization problem is constrained by a set of inequality constraints: min f (X) (17) subject to the SPSA algorithm can be applied as well, but the update of the vector X accounts for the violated constraints [26]. Note that the SPSA method does not directly handle equality constraints, whereas a set of equality constraints (Equation (7) with the positions of Equation (8)) appears in the optimization problem to be solved in the outer routine.
Fortunately, these equations can be separately considered [27] (3rd Step of the proposed procedure). In fact, once the daily scheduling of EESSs is known after performing the inner routine, Equation (7) (including, of course, the conditions stated by (8)) can be separately solved in order to determine the phase voltages (magnitudes and arguments) at all buses. The latter allows, then, to verify whether inequality constraints among Equations (9-12) are met or not, discarding unfeasible siting and sizing solutions. In the case of an unfeasible solution, the update of the vector X accounts for the violated constraints, as previously mentioned.

Inner Algorithm: the EESSs Daily Scheduling
The optimal daily scheduling of EESSs (i.e., charging/discharging cycles of EESSs) along with the energy size of the EESSs are obtained by applying a simplified but effective and fast algorithm, based on a busbar-system representation of the network (Figure 2). The losses of the distribution system are not taken into consideration, allowing it in this way to avoid time-consuming algorithms (e.g., three-phase load flow) and to adopt fast operations (we show once again that the three-phase load flow equations (and the losses) are taken into account in Step 3).
The power rating of the equivalent EESS in Figure 2 is given by where the size P size s of the EESS installed at bus s is fixed in the outer routine by SPSA for s ∈ Ω S and the equivalent load in Figure 2 incorporates all the loads of the original system under study.
system are not taken into consideration, allowing it in this way to avoid time-consuming algorithms (e.g., three-phase load flow) and to adopt fast operations (we show once again that the three-phase load flow equations (and the losses) are taken into account in Step 3).
The power rating of the equivalent EESS in Figure 2 is given by ∑ ∈ where the size of the EESS installed at bus s is fixed in the outer routine by SPSA for ∈ Ω and the equivalent load in Figure 2 incorporates all the loads of the original system under study. With reference to the system in Figure 2, the inner routine solves an optimization problem in which the objective function to be minimized accounts for the total costs of Equation (2), given by the sum of the costs of the energy bought from the upstream network (which depends on the daily scheduling of the equivalent EESSs, once the equivalent load is assigned) over the planning period and the EESSs costs. A set of EESSs constraints are also considered, as shown hereinafter. We assume that a tariff scheme is assigned as input data of the optimization problem. Each day is divided into tariff periods, each characterized by an energy cost. The tariff periods depend on the season of the year. Typically, this tariff scheme includes on-peak hours (which are characterized by the greatest energy prices), off-peak hours (characterized by the smallest energy prices) and part-peak hours (i.e., the hours during which energy prices are comprised between on-peak and off-peak prices).
The decision variables of the optimization problem (output of the inner routine) are the values of the active powers in the charging/discharging stages of the equivalent EESS for all of the hours of the typical days in years; such values are linked to the energy size of the equivalent EESS, as it is shown in the following. In order to reduce the computational efforts, some typical days (e.g., working day, Saturday, Holidays) are considered in each season of the year.
In particular, assuming the nominal discharging time of the equivalent EESS to be a discrete variable, the optimization problem of the inner routine is solved by applying an exhaustive procedure; i.e., considering all of the possible integer values of the nominal discharging time (from one hour, up to a maximum value ), and selecting the optimal value to be the one associated to the smallest value of the objective function Equation (2) involving the total costs.
The maximum nominal discharging time can be determined as the ratio of the greatest value of maximum energy that can be charged (during low-price hours) and discharged (during highprice hours) over a day without violating any constraints, to the (the integer rounded value is taken). Its value is constrained into a range defined by the EESS technology During the exhaustive procedure, once is fixed, for each typical day of each year of the planning period, the charging/discharging cycle of the equivalent EESS is determined: (i) by reducing the energy consumption during hours with greater prices, and by increasing the energy purchase during hours with smaller prices; and (ii) by implementing a strategy aiming at load levelling. In particular, when the EESS is charging or discharging, the equivalent load curve is possibly flattened.
More in detail, the most inexpensive daily operation of the equivalent EESS is obtained by moving as much energy as possible from the on-/part-peak hours to the off-peak hours, charging batteries during the off-peak hours and discharging them during the on-/part-peak hours. The peak power during the off-peak hours is straightforwardly adapted. Moreover, the inner algorithm tries to perform a load levelling, when it is possible, during the time intervals of each tariff period. With such characteristics, even if we operate in a busbar system (in which losses are neglected), we are With reference to the system in Figure 2, the inner routine solves an optimization problem in which the objective function to be minimized accounts for the total costs of Equation (2), given by the sum of the costs of the energy bought from the upstream network (which depends on the daily scheduling of the equivalent EESSs, once the equivalent load is assigned) over the planning period and the EESSs costs. A set of EESSs constraints are also considered, as shown hereinafter. We assume that a tariff scheme is assigned as input data of the optimization problem. Each day is divided into tariff periods, each characterized by an energy cost. The tariff periods depend on the season of the year. Typically, this tariff scheme includes on-peak hours (which are characterized by the greatest energy prices), off-peak hours (characterized by the smallest energy prices) and part-peak hours (i.e., the hours during which energy prices are comprised between on-peak and off-peak prices).
The decision variables of the optimization problem (output of the inner routine) are the values of the active powers in the charging/discharging stages of the equivalent EESS for all of the hours of the typical days in years; such values are linked to the energy size of the equivalent EESS, as it is shown in the following. In order to reduce the computational efforts, some typical days (e.g., working day, Saturday, Holidays) are considered in each season of the year.
In particular, assuming the nominal discharging time of the equivalent EESS to be a discrete variable, the optimization problem of the inner routine is solved by applying an exhaustive procedure; i.e., considering all of the possible integer values of the nominal discharging time t disch (from one hour, up to a maximum value t max disch ), and selecting the optimal value to be the one associated to the smallest value of the objective function Equation (2) involving the total costs.
The maximum nominal discharging time can be determined as the ratio of the greatest value of maximum energy that can be charged (during low-price hours) and discharged (during high-price hours) over a day without violating any constraints, to the P EESS (the integer rounded value is taken). Its value is constrained into a range defined by the EESS technology During the exhaustive procedure, once t disch is fixed, for each typical day of each year of the planning period, the charging/discharging cycle of the equivalent EESS is determined: (i) by reducing the energy consumption during hours with greater prices, and by increasing the energy purchase during hours with smaller prices; and (ii) by implementing a strategy aiming at load levelling. In particular, when the EESS is charging or discharging, the equivalent load curve is possibly flattened.
More in detail, the most inexpensive daily operation of the equivalent EESS is obtained by moving as much energy as possible from the on-/part-peak hours to the off-peak hours, charging batteries during the off-peak hours and discharging them during the on-/part-peak hours. The peak power during the off-peak hours is straightforwardly adapted. Moreover, the inner algorithm tries to perform a load levelling, when it is possible, during the time intervals of each tariff period. With such characteristics, even if we operate in a busbar system (in which losses are neglected), we are confident that levelling the load leads to a reduction of the mean required power and, consequently, we also expect active power losses to be reduced.
As previously evidenced, the optimization problem solved in the inner routine includes a set of constraints that have to be verified when the exhaustive procedure is applied; the EESSs constraints of Section 2 (expressed for the equivalent EESS) and the following ones: i capacity of the grid: the "updated" daily curves (provided by the sum of the load powers and of the EESS power) do not have to exceed the peak power of the loads, ii exportation is not allowed: when the EESS is discharging, power cannot flow toward the main grid, These must be verified. Two examples are reported hereinafter in order to better clarify the strategy implemented to optimally charge and discharge the equivalent EESS, with reference to the three tariff periods (i.e., off-peak, part-peak, and on-peak hours) in which the smallest price occurs during the off-peak hours.
In the next figures, P on peak,d,y and P o f f peak,d,y are the constant values of power during the on peak hours (h ∈ Ω disch,d,y ) and the off-peak hours (h ∈ Ω ch,d,y ), respectively. P on peak,d,y and P o f f peak,d,y are determined by considering the contribution of the energy charged (h ∈ Ω ch,d,y ) and discharged (h ∈ Ω disch,d,y ) by the equivalent EESS. In Figure 4, P part peak,d,y is the constant value of power during the part-peak hours.
Different cases can occur, on the basis of the rating of the equivalent EESS (i.e., P EESS ), on the basis of the value of the nominal discharging time, and, of course, on the basis of the assigned load curve. Two of these cases are discussed hereinafter, considering a specified typical day.
In the first case we assume t disch equal to 6 h; as shown in Figure 3, the equivalent EESS discharges during the peak hours in such a way that a levelled power, i.e., P on peak,d,y , is obtained. Similarly, the equivalent EESS is charged during the off-peak hours in such a way that a levelled power, i.e., P o f f peak,d,y is obtained.
In the second case we assume the power rating and the nominal discharging time of the equivalent EESS to be significantly greater than in the previous case. As shown in Figure 4, the energy stored by the equivalent EESS is greater than the energy required by the loads during the peak hours; therefore, part of this energy can be also used in the part-peak hours. As a consequence, P on peak,d,y is zero, and the EESS discharges during some hours of the part-peak period. This allows it not to draw energy from the main grid when prices are high. In particular, note that in the case represented in Figure 4, constraints on the minimum and maximum values of the "updated" load curve (the one that include the storage) are binding. In fact, during the charging phase in the off-peak hours, the "updated" load curve has a maximum equal to the peak value of the original load curve due to the condition (i) reported above. Moreover, during the discharging phase in the on-peak hours, the levelled load is zero due to the condition (ii) reported above.
Once the optimization problem of the inner routine is solved and the charging/discharging cycle of the equivalent EESS is obtained for all of the days, we assume that all of the EESSs that constitute the equivalent EESS in Figure 2 operate with the same charging/discharging starting times, and with active powers of each of them scaled of the quantity P size s /P EESS (s ∈ Ω S ). Regarding the reactive power, in order to support the network operation, the converters of the storage systems are controlled to provide the maximum reactive power compatible with the constraint Equation (6) while assuring not to inject reactive power into the upstream network.

Micro Genetic Algorithms
The μGA used in [28] is applied in this paper. The μGA, which is used to reduce the processing time required by simple GAs, explores the possibility of working with small populations. The μGA evolves with populations of only five individuals, and it uses the selection and the crossover. The mutation is not applied, since diversity is guaranteed by periodically refreshing the population. The replacement of the population (except the replacement of the best individual) is also performed [28].
In particular, the selection is based on the roulette wheel method. The amount of diversity of the population, after the application of genetic operators, is measured by counting the total number of genes that are unlike the genes of the best individual. When the diversity of the population is smaller than a selected threshold, four individual populations are deleted and replaced by new randomly-

Micro Genetic Algorithms
The μGA used in [28] is applied in this paper. The μGA, which is used to reduce the processing time required by simple GAs, explores the possibility of working with small populations. The μGA evolves with populations of only five individuals, and it uses the selection and the crossover. The mutation is not applied, since diversity is guaranteed by periodically refreshing the population. The replacement of the population (except the replacement of the best individual) is also performed [28].
In particular, the selection is based on the roulette wheel method. The amount of diversity of the population, after the application of genetic operators, is measured by counting the total number of genes that are unlike the genes of the best individual. When the diversity of the population is smaller than a selected threshold, four individual populations are deleted and replaced by new randomly-

Micro Genetic Algorithms
The µGA used in [28] is applied in this paper. The µGA, which is used to reduce the processing time required by simple GAs, explores the possibility of working with small populations. The µGA evolves with populations of only five individuals, and it uses the selection and the crossover. The mutation is not applied, since diversity is guaranteed by periodically refreshing the population. The replacement of the population (except the replacement of the best individual) is also performed [28].
In particular, the selection is based on the roulette wheel method. The amount of diversity of the population, after the application of genetic operators, is measured by counting the total number of genes that are unlike the genes of the best individual. When the diversity of the population is smaller than a selected threshold, four individual populations are deleted and replaced by new randomly-generated individuals. The next step is performed with this new population, made up of • a control on the fitness improvement provided by the next solution; or • a maximum number of generated individuals.
In the first case, the algorithm s processing ends when the best fitness value remains stable for a certain number of generations; a maximum number of generations is previewed after the algorithm ends. In particular, in order to avoid premature convergences, this verification is performed after the diversity check.
The second condition stops the algorithm at an assigned finite number of generated individuals (200 in this paper).

Case Study
The planning problem related to the sizing and siting of EESSs is solved for the IEEE unbalanced 34-bus test system shown in Figure 5. This system contains a mixture of single-and three-phase lines and loads. Note that the system lines 808-810, 816-818, 818-820, 820-822, 824-826, 854-856, 858-864, and 862-838 are single-phase and the remaining lines are three-phase. The complete network data and parameters are in [29].
The EESSs can be allocated at single-phase and three-phase buses. Time-of-Use pricing of energy bought from the utility is assumed, with the tariffs reported in Table 1 [30]. These tariffs are chosen since they are characterized by a significant spread between on-peak and off-peak prices, which nowadays is a mandatory condition to profitably use storage systems.
Three typical days (i.e., working day, Saturday, Holiday) and four seasons are considered. As a result, 24 typical days are assumed [31]. A rate of increase equal to 2% was assumed in order to account for an increase of load over the planning period. The peak powers are set at 70% of the nominal power of the load in [29].
Both the effective rate of change and the discount rate are assumed equal to 3%. Two operation modes of energy storage systems are considered: • Mode 1: The energy storage systems are used only in the summer months: during the off-peak hours they can charge, and during the rest of the day the battery can discharge. The storage systems also exchange the reactive power subject to the constraints of Equation (6). • Mode 2: For each day of the year, the energy storage systems can charge during the off-peak hours, and they can discharge during the remaining hours. Both active and reactive powers can be exchanged subject to the constraints of Equation (6).  a control on the fitness improvement provided by the next solution; or  a maximum number of generated individuals.
In the first case, the algorithmʹs processing ends when the best fitness value remains stable for a certain number of generations; a maximum number of generations is previewed after the algorithm ends. In particular, in order to avoid premature convergences, this verification is performed after the diversity check.
The second condition stops the algorithm at an assigned finite number of generated individuals (200 in this paper).

Case Study
The planning problem related to the sizing and siting of EESSs is solved for the IEEE unbalanced 34-bus test system shown in Figure 5. This system contains a mixture of single-and three-phase lines and loads. Note that the system lines 808-810, 816-818, 818-820, 820-822, 824-826, 854-856, 858-864, and 862-838 are single-phase and the remaining lines are three-phase. The complete network data and parameters are in [29].
The EESSs can be allocated at single-phase and three-phase buses. Time-of-Use pricing of energy bought from the utility is assumed, with the tariffs reported in Table 1 [30]. These tariffs are chosen since they are characterized by a significant spread between onpeak and off-peak prices, which nowadays is a mandatory condition to profitably use storage systems.
Three typical days (i.e., working day, Saturday, Holiday) and four seasons are considered. As a result, 24 typical days are assumed [31]. A rate of increase equal to 2% was assumed in order to account for an increase of load over the planning period. The peak powers are set at 70% of the nominal power of the load in [29].
Both the effective rate of change and the discount rate are assumed equal to 3%. Two operation modes of energy storage systems are considered:  Mode 1: The energy storage systems are used only in the summer months: during the off-peak hours they can charge, and during the rest of the day the battery can discharge. The storage systems also exchange the reactive power subject to the constraints of Equation (6).  Mode 2: For each day of the year, the energy storage systems can charge during the off-peak hours, and they can discharge during the remaining hours. Both active and reactive powers can be exchanged subject to the constraints of Equation (6).   With reference to the constraints, the maximum line currents are fixed at the ratings reported in [29], and the maximum value of the unbalance factor at each bus is set at 3%. The minimum and the maximum values of the voltage at each phase of each bus are set at 90% and 110% of the nominal value, respectively.

Analysis of Several Technologies
Several technologies can be considered for the batteries used in the EESSs. They have different characteristics in terms of:

•
Energy and power installation costs • Electrochemical properties (energy density, power density) • Costs evolution • Performances.
From the analysis of the data reported in [32], it is evident that NaS (Sodium-Sulfur) batteries are associated with the smallest costs, whatever their operation mode. However, the EESSs based on these batteries are available starting from 800 kVA/4.8 MWh size [33]; therefore, they are not compatible with the distributed EESSs concept analyzed in this paper for the network under study.
Na-NiCl 2 batteries can be considered as the candidate technology for EESSs due to their modularity and to economic considerations [32].
With respect to EESSs, the following assumptions stand: • The unit storage system available at any phase of each bus is assumed to come in discrete sizes of 50 kVA; • The standard value of the power/energy ratio for Na-NiCl 2 batteries is generally around 1/3 [34]; • The efficiencies in charging and discharging modes are set at 0.9 [34]; • The installation cost at year 1 is assumed to be 400 $/kWh [32]; the expected evolution of battery costs in the next years is also provided in [32]. Based on this analysis, in this application, the replacement cost was evaluated for each year of the planning period according to [32].

Results
Tables 2 and 3 report the results obtained for Mode 1 and Mode 2 operations, whereas Table 4 reports the values of the objective function for the optimal configurations of Tables 2 and 3 (in per unit (p.u.) value of the objective function without any EESSs).
We point out from the analysis of the results reported in Tables 2-4 that: • the optimal value of the power/energy ratio is always 1/6; this value is adequate for Na-NiCl 2 batteries.

•
The Mode 1 allows it to obtain smaller total costs. These results confirm that there is convenience in using the storage systems only in the most adequate conditions, due to the non-negligible cycling costs.  Table 3. Optimal location and size of energy storage systems for Mode 2.

Storage systems Location and Size
Three-phase storage systems 150 kVA/900 kWh at bus #806 150 kVA/900 kWh at bus #844 300 kVA/1800 kWh at bus #858 Single-phase storage systems 50 kVA/300 kWh at bus #818 A further analysis aims at evaluating the performance of the optimal configuration of Table 2 when the operating Mode 2 is implemented, and at evaluating the performance of the optimal configuration of Table 3 when the operating Mode 1 is implemented. The results are reported in Table 5. The EESSs have a positive effect on voltage amplitudes and on the unbalance factor as well. Of course, the voltage amplitude (as well as the unbalance factor) varies with the node, the day, and the hour. The minimum, the mean, and the maximum values of the voltage amplitudes during one year are shown in Figure 6 for all of the buses and for the three phases. The minimum, the mean, and the maximum values of the unbalance factor during one year are shown in Figure 7 for all of the three-phase buses.   The planning problem of the optimal sizing and siting of EESSs is also solved by means of the μGA. In this case, the outer routine applies the μGA, in which the EESSs sites and sizes (powers) are discrete variables. The optimal location and size for the operation of Mode 1 are reported in Table 6, and the objective function of the optimal configuration (in per unit value to the objective function without any EESS) is 0.8971 p.u. The result is very close to the one obtained by applying the SPSA, in spite of a much longer computational time. In particular, the proposed procedure (based on the application of the SPSA algorithm) has required a computational time of about 60 hours, which is 30 times lower than that requested by the μGA. These times must be handled with care because our software was not optimized for computational speed and the simulations were carried out with Matlab (MathWorks, Natick, MA, USA) programs by a 3.6 GHz-16 GB RAM PC Xeon processor E3-1280v2 (Intel corporation, Santa Clara, CA, USA). Nowadays, massive computation is becoming accessible with the new machines and configurations (parallel distributed processing and environment). Table 6. Optimal location and size of energy storage systems for mode 1 by applying the μGA. Eventually, an exhaustive search is performed for a single EESS NaS installation (800 kVA/4.8 MWh) in order to check the benefits of the distributed EESS concept. The exhaustive search results in an optimal value of the objective function equal to 0.8991 p.u., and the optimal location is at the node #832. Therefore, despite the smaller cost of the NaS technology, the distributed concept allows to obtain better performances.

Conclusions
In this paper, an optimization problem for the optimal siting and sizing of battery systems in an unbalanced electrical distribution systems is formulated and solved with a simplified hybrid algorithm. The algorithm significantly reduces the computational efforts involved in the solution of The planning problem of the optimal sizing and siting of EESSs is also solved by means of the µGA. In this case, the outer routine applies the µGA, in which the EESSs sites and sizes (powers) are discrete variables. The optimal location and size for the operation of Mode 1 are reported in Table 6, and the objective function of the optimal configuration (in per unit value to the objective function without any EESS) is 0.8971 p.u. The result is very close to the one obtained by applying the SPSA, in spite of a much longer computational time. In particular, the proposed procedure (based on the application of the SPSA algorithm) has required a computational time of about 60 hours, which is 30 times lower than that requested by the µGA. These times must be handled with care because our software was not optimized for computational speed and the simulations were carried out with Matlab (MathWorks, Natick, MA, USA) programs by a 3.6 GHz-16 GB RAM PC Xeon processor E3-1280v2 (Intel corporation, Santa Clara, CA, USA). Nowadays, massive computation is becoming accessible with the new machines and configurations (parallel distributed processing and environment). Table 6. Optimal location and size of energy storage systems for mode 1 by applying the µGA. Eventually, an exhaustive search is performed for a single EESS NaS installation (800 kVA/4.8 MWh) in order to check the benefits of the distributed EESS concept. The exhaustive search results in an optimal value of the objective function equal to 0.8991 p.u., and the optimal location is at the node #832. Therefore, despite the smaller cost of the NaS technology, the distributed concept allows to obtain better performances.

Conclusions
In this paper, an optimization problem for the optimal siting and sizing of battery systems in an unbalanced electrical distribution systems is formulated and solved with a simplified hybrid algorithm. The algorithm significantly reduces the computational efforts involved in the solution of the high dimensional mixed integer and nonlinear optimization problem by applying (i) the Simultaneous Perturbation Stochastic Approximation method and (ii) an inner simplified procedure that allows it to quickly carry out the daily scheduling (charging/discharging) of the EESSs. The numerical application to an unbalanced IEEE test network clearly shows the feasibility and effectiveness of the proposed solution in terms of accuracy and reduced computational time. The results are compared with the ones obtained by applying a micro Genetic Algorithm.
The proposed approach gives insight into the technical and economic advantages of using distributed electrical energy storage systems in an unbalanced distribution network.
In the present paper, the size of EESSs was assumed as a discrete variable; however, it can be considered as a continuous variable as well. Since this point is worth further investigation, a future work will address the comparison of the discrete and continuous sizing of EESSs.
Future research on the subject will extend the proposed analytical approach to take into account the different sources of uncertainties involved in the EESSs allocation problem.