*2.1. Fitness Function Assessment in the NSGA-II*

A typical genetic algorithm simulates the mechanism of the natural evolution, by applying systematically the three genetic operators of *selection* (that gives more chance of reproduction to the better individuals), *crossover* (that generates offspring solutions by mixing the genetic characteristics of parents), and *mutation* (that implements random alterations of the genetic characteristics) to evolve the population towards the global optimum. Often, *elitism* (i.e., only the better solutions are considered in the evolution process) is implemented in the formation of the new population in order to increase the effectiveness of the optimization procedure and speed up the convergence of the algorithm.

Taking account of these aspects, the mechanism implemented for comparing two solutions and identifying the better one is crucial for any evolutionary algorithm. The classification procedure used by the NSGA-II algorithm is based on the definition of two attributes associated with each individual: the *non-domination rank* and the *crowding distance*. The *non-domination rank* organizes the candidates into subsets of individuals non-dominated by any other in the same subset (fronts of non-dominance). Therefore, the first non-dominated front is formed by all the individuals of the Pareto set for which no other solutions are either equal to or better than them on all of the search objectives, and it is marked with the lowest (best) rank. The other fronts are sorted depending on the number of subsets from which they

are dominated: e.g., the second front is dominated only by the first and dominates all the remaining fronts, the third front is dominated only by the first two and so on. In order to guarantee the diversity in each front, a second attribute is introduced, the *crowding distance* (CD), based on the cardinality of the solution sets and their distance to the solution boundaries. It is calculated by summing the absolute normalized differences in the function values of two adjacent solutions; the sum is extended to each objective function (*OF*), as indicated in (1):

$$CD\_i = \sum\_{k=1}^{N\_{\text{OF}}} \frac{OF\_k^{i+1} - OF\_k^{i-1}}{OF\_k^{\text{max}} - OF\_k^{\text{min}}} \tag{1}$$

where *NOF* is the number of *OF*s and *OF<sup>i</sup> <sup>k</sup>* is the *<sup>k</sup>*th *OF* value of the *<sup>i</sup>* th generic individual. The higher the *CD* of a solution, the less crowded the corresponding area of the front and, hence, the finer its fitness value.

The application of these two attributes allows finding a better solution in a pairwise comparison, assuring the correct implementation of selection operator and elitism in the evolution of the population.

#### *2.2. Solution Coding for the BESS Location Problem*

The traditional binary coding makes the GAs particularly suitable for solving sizing and siting problems of different resources, like (for power distribution systems) capacitors, distributed generators, or BESS. However, in the case of storage devices, the effectiveness of their support to network operation depends also on their hourly state of charge (SoC). Therefore, the solution coding has to represent the BESS daily scheduling as well, so that the NSGA-II can optimize the design and normal operation of storage devices simultaneously, finding full compromise solutions among different point of views. Furthermore, as it is explained in the next paragraph, part of the BESS capacity is reserved to DSO for local distribution network support. Consequently, also this capacity share has to be included among the chromosome information of the generic solution. An example of the solution coding referred to a single BESS is reported in Figure 1. The whole chromosome of a generic individual is obtained by repeating this schema for a prefixed maximum number of BESS (*NBESS*).

**Figure 1.** Chromosome section for a storage device.

In summary, the genetic parameters optimized simultaneously by the NSGA-II are:

#### 1. **BESS position in the distribution network**

It corresponds to any MV network node, excluding the primary substation busbars because in those locations the BESS may affect only the loading of the HV/MV transformers, but it cannot provide support to local network contingencies. For real network studies, it is easy to exclude also additional MV nodes due to any kind of constraint the DSO planner wants to consider. The value of the specific gene is an integer number in the interval [1, *Nnode*], where *Nnode* is the total number of available MV nodes.

#### 2. **BESS rate**

It is defined by two genes: one for the nominal power (*Pn*) and one for the nominal duration (*dn*), so identifying by their product the nominal BESS capacity (*Cn* = *Pn* · *dn*). Nominal power and duration are

the integer number within a minimum and maximum values ([*Pn\_min*, *Pn\_max*] and [*dn\_min*, *dn\_max*]), given among the input data of the problem. Generally, *dn\_min* is the elementary time step used to represent the load/generation daily patterns (1 h), *Pn\_min* is the minimum power required by the national regulation to participate in the ancillary services market, and *Pn\_max* is the maximum power permitted for the direct installation of BESS (or any other generation unit) on the MV electric distribution system.

#### 3. **BESS daily energy schedule**

It is represented by the sequence of the State of Charge (SoC) in each interval of the day. By using a time step (Δ*t*) of 1 h, there are 24 genes for each typical profile used to represent the yearly customer behavior. Therefore, if more than one typical day is used (e.g., two semester profiles, four seasonal profiles, or multiple profiles for seasons, working days and weekends), the sequence of 24 genes is repeated accordingly. The SoC values are real numbers within the interval [0%, 100%]. However, these percentages are not related to *Cn*, but to the remaining capacity curtailed of the share for DSO (*Cown\_TSO* in Figure 2). This choice allows preserving the soundness of the chromosome representation for the offspring solutions, built with the application of the genetic crossover and mutation operators. Indeed, by so doing, the new sequences of SoC remain always within the operational limits, avoiding overlapping with the semi-bands of capacity reserved to DSO.

**Figure 2.** Decoding of the chromosome of Figure 1: BESS with *Pn* = 1 MW and dn = 4 h; *CDSO* = 800 kWh (20% of Cn) reserved for DSO; upper and lower semi-bands of 400 kWh; depending on the SoC, the virtual storage can be used as a BESS of 800 kW × 1 h or 400 kW × 1 h.

#### 4. **Share of BESS capacity for DSO**

It establishes the rate of the "virtual storage" for DSO. From this value (CDSO), expressed as a percentage of Cn, the upper and lower semi-bands of capacity are calculated (Cmin\_DSO = CDSO/2; Cmax\_DSO = Cn − CDSO/2) and always reserved to DSO in order to guarantee in any moment the availability of the virtual storage capacity in case of contingencies that could bring the distribution system to a critical emergency configuration (Figure 2). The residual capacity, Cown\_TSO = Cn − CDSO, is managed by the private owner for energy trading and for offering ancillary services to TSO. The maximum available power of the virtual storage (Pmax\_DSO) may change hour by hour because it depends on the current SoC. If it is around half of Cn, it is reasonable the potential exploitation of the highest value (Pmax\_DSO = CDSO/Δt, with Pmax\_DSO ≤ Pn), while when the SoC is coincident to the lower or upper bound (Cmin\_DSO or Cmax\_DSO), the maximum available power takes its lowest value (Pmax\_DSO = CDSO/2Δt, with Pmax\_DSO ≤ Pn). Again, for this gene, the value is a real number within the interval [0%, 100%].

It is worth noticing that the decision variables (e.g., BESS SoC, share of Cn, rated power, and energy or duration) are continuous and, therefore, also the MO optimization solution space cannot be binary, as in the simplest version of the algorithm in the literature. For this reason, it has been decided to adopt a more efficient real codification. Indeed, this choice eliminates the accuracy problems related to the discretization of the solution space proper of the binary formulation that may not result adequately precise [27]. Furthermore, the formulation of genetic algorithms with the real instead of the binary coding is well suited for dealing with the gradual trend exhibited by continuous variable functions (i.e., small variations in the variables slightly alter the function values). Its implementation has required some adaptations to the typical operators of genetic algorithms (i.e., crossover and mutation), as detailed in [8].

#### **3. Objective Functions**

The advent of flexible resources in the distribution system during the deregulated power system era introduces new challenges to be faced, among which the correct cooperation between TSO and DSO is a cornerstone. Indeed, flexibility is essential for preserving the secure operation of the power system, but it is becoming central also for development and management of upcoming smart distribution networks. Moreover, the uncoordinated exploitation of flexibility from DER for TSO requirements may cause severe technical challenges to the distribution system, not originally planned for this scope.

Diverse TSO/DSO cooperation models are possible and have been analyzed by several scientific working groups [28,29], like various services market models as well [6]. In the paper, it has been assumed the opening of the ancillary services market (ASM) to the BESS located in the distribution system for the provision of primary and secondary frequency reserve. No local ASM has been supposed on the distribution level, but the DSO may request to reserve part of the BESS capacity for facing local contingencies. Two stakeholders have been considered with their goals:


In the MO formulation, the DSO's objective function (OF) has been defined as the annual risk to violate any technical constraints (to be minimized), whereas for the BESS owners the OF is the cost/benefit ratio (to be minimized as well).

By optimizing the two OFs, it is possible also to analyze the limitation on the availability of the flexibility services for the TSO caused by the bottlenecks in distribution network. Indeed, if the location of a particular BESS is convenient for the distribution system, the DSO will reserve part of the available flexibility to resolve local contingencies, so reducing the amount of flexibility accessible to TSO.

## *3.1. DSO Objective Function—Risk of Technical Constraints Violation*

Modern planning tools for designing the upcoming smart distribution networks need to overcome the traditional and extra conservative "fit and forget" approach, built around the fulfilment of the worst-case scenario, to move towards a risk-based procedure, that can correctly correlate the effectiveness of planning choices to the probability and the seriousness of possible contingencies [2].

An explicit and detailed assessment of the annual risk to violate technical constraints in a given distribution network requires probabilistic load flow (PLF) calculations, solved for each time step of the typical days used to represent the yearly behavior of distribution network customers [30]. The load/generation behavior in each hour has been assumed normally distributed. The risk assessment procedure starts with the definition of the impedance matrix [Z]*b*, relative to the *b*th network configuration

in the N-1 security analysis (where *b* = 0 means the network in normal configuration and *b* > 0 means the network reconfigured without the *bt*<sup>h</sup> element), and the acquisition of the nodal current matrix [Inode] *<sup>f</sup>* in the *f* th hour of the typical daily profile. The results of the PLF are the nodal voltage [Vnode] *<sup>f</sup>* and the branch current [Ibranch] *<sup>f</sup>* matrixes (expressed in terms of mean value, μ, and standard deviation, σ, of a normal distribution), by which the probability (*ptcv*) to overcome the voltage regulation band or the conductor thermal limit is estimated. Only the *Nc* operating conditions with non-negligible probability to violate the technical constraints (*ptcv >* 0) are stored (Figure 3—case A), while the cases on which the extremes values of [Vnode] *<sup>f</sup>* and [Ibranch] *<sup>f</sup>* (assumed equal to <sup>μ</sup> ± <sup>3</sup>σ) do not exceed the technical limits (minimum and maximum nodal voltages, *Vlim\_min* and *Vlim\_max*, and maximum branch current, *Ilim\_max*) are disregarded (Figure 3—case B).

**Figure 3.** Identification of potential contingencies (*ptcv*) and procedure for the total risk assessment.

In order to calculate the risk of technical constraints violation (*Rbf*) when the *b*th network configuration is in force during the *f* th hour of the *dth* typical day, the *ptcv* is multiplied by the occurrence probability (*pbf*) of the relative operating conditions. Such probability can be determined by simply multiplying the forced outage rate of the *b*th network element (*FORb*) and the occurrence probability of the specific customers' operating conditions (*pfd*), because these two probabilities can be considered independent:

$$FOR\_b = \frac{MTTR\_b}{MTTF\_b + MTTR\_b} \qquad p\_{fd} = \frac{nh\_{fd}}{8760},\tag{2}$$

where:

• *MTTRb* is the "Mean Time To Repair" of the bth network element, indicated in the paper with the symbol *τ<sup>b</sup>* and assumed equal to 5 h for an overhead line and 8 h for a buried cable;


For a distribution network, it is evident that *MTTF >> MTTR* (years compared to few hours). Consequently, *MTTR* can be disregarded in the denominator of the first of Equation (2), and it is acceptable to assess the occurrence probability *pb f* with the following approximated formula:

$$p\_{bf} = \left(\tau\_b \cdot \frac{\lambda\_b}{8760}\right) \cdot \left(\frac{nh\_{fd}}{8760}\right). \tag{3}$$

When the normal configuration is examined (*b* = 0), *FORb* is assumed equal to 1 and *pbf* = *pfd*. Finally, the risk component *Rbf* is expressed in hours of violation per year:

$$R\_{bf} = \left[ p\_{bf} \cdot p\_{tcv} \cdot 8760 \left[ \frac{\text{hours}}{\text{year}} \right] \right]. \tag{4}$$

The sum of all these risk components, determined for each configuration in each interval of the typical days, represents the total risk (*RTOT*) that characterizes the distribution network examined, i.e., the number of hours per year when it is possible to overcome a technical constraint.

When the total risk is greater than the acceptable limit fixed by the DSO planner, *RA*, planning solutions have to be put in place to reduce *RTOT* below *RA* and make the distribution network robust enough for the whole planning period considered. Obviously, this system development has to be optimized by comparing the conventional network reinforcement with the exploitation of flexibility services from DERs, not only in terms of costs but also of residual risk. In the paper, only the resort to services from BESS has been considered, by exploiting when necessary the share of storage capacity (CDSO%) reserved to DSO. A simple linear programming optimal power flow (OPF), capable of finding the optimal scheduling of BESS, is used for tackling the particular contingency and nullifying the residual risk [3]. If the BESS location and/or the share of capacity for DSO are ineffectual for solving the contingency, the existing risk component will be reduced according to the flexibility resource availability.

#### *3.2. BESS Owners' Objective Function—Services for TSO*

The recent trend in opening the ancillary services market to resources located in the distribution system is creating new business opportunities for private investors. Indeed, the large diffusion of small non-dispatchable renewable generation units and the concurrent decommissioning of traditional fossil-fuelled plants are increasing the need for additional flexible and fully controllable resources, like BESS. Therefore, the introduction of new actors in the power system is expected, whose goal is to make a profit with storage devices mainly by providing ancillary services to power system operators. Obviously, the participation to the energy market is available as well, and it has been included in the following cost-benefit analysis (CBA), even if its weight is lower than the income from ancillary services provision.

In the paper, with the local ASM for the distribution system not yet available, the ancillary services provision has been assumed devoted only to TSO, even if DSO may limit this operating mode by acquiring part of the BESS capacity. Thus, the BESS owners' sources of income considered in the calculations are related to the following three applications:


These BESS applications are mutually dependent, besides limited by the DSO reserve capacity, because they may be provided simultaneously, and their value relies on the active power exchanged with the grid. By optimizing the BESS daily scheduling of the stored energy, the multi-objective procedure adopted in the paper finds the optimal share of the BESS active power between arbitrage (and the influence on the distribution network operating conditions) and TSO services. For the further partition between the two frequency ancillary services, the heuristic choice of equally dividing the residual active power was taken.

By assessing the incomes from these three applications, in terms of net present values (NVP) within the whole planning period, and the storage investment (*CBESS*), the BESS owners' point of view has been expressed as a cost-benefit ratio:

$$OF\_{BESS\\_ovuures} = \frac{\mathcal{C}\_{BESS}}{(B\_{arb} + B\_{FCR} + B\_{aFRR})} \,. \tag{5}$$

In the following, more details on the determination of cost and benefits are provided.
