Sizing and Sitting of Battery Energy Storage Systems in Distribution Networks with Transient Stability Consideration

: In this study, the capacity and location of battery energy storage systems (BESSs) in a distribution network were evaluated to increase the stability and reliability of power systems by applying the proposed transient stability indicators. The search capability of particle swarm optimization (PSO) combined with Pareto optimality in MATLAB R2019b was used to solve multi-objective problems in steady-state conditions, including active power loss, voltage deviation, and total operating cost. The BESS capacity and location were then derived using the Manhattan distance method. Subsequently, the BESS control model was constructed using the commercial software DIgSILENT PowerFactory 2021 to evaluate generator power smoothing and short-term voltage stability index during the transient response and determine the ﬁnal capacity and location of BESS in distribution networks. Finally, the accuracy of the proposed method was veriﬁed by considering the daily load and solar power generation curves using the IEEE 33-bus radial power distribution network.


Introduction
With the rise in environmental awareness in recent years and the growing problem of global energy supply, governments around the world are actively developing renewable energy sources. The capacity of renewable energy installations worldwide is increasing rapidly, with solar and wind power being the mainstream. According to the Renewables 2020 Global Energy Status Report, the capacity of renewable energy installations grew by more than 200 GW in 2019. Currently, renewable energy has become an integral part of power generation. The International Renewable Energy Agency (IRENA) [1] analyzed that 260 GW of new renewable energy installations was added globally in 2020, showing a significant increase from the 176 GW of new renewable energy installations in 2019. The increasing number of installations is because renewable energy can be used to improve peak electricity demand and reduce long-distance transmission line losses and transmission line congestions. However, grid integration of numerous renewable energy sources may bring unforeseen problems to the frequency, voltage, and stability of power systems. Some countries have drawn up relevant regulations to improve the impact of grid-connected renewable energy sources; such regulations were compiled in [2]. Moreover, the Institute of Electrical and Electronics Engineers (IEEE) proposed specifications for the grid-connected features, safety considerations, and island operations of distributed energy resources (DERs) at primary or secondary distribution voltage levels in IEEE Std. 1547 [3]. Additionally, Germany introduced VDE-AR-N 4105 as the new standard for low voltage grid-connected features [4]. The application guide is used for the connection of DERs with low voltage levels and capacities equal to or less than 100 kVA. Moreover, China used GB-T 19964 as output, the active power was injected into the system through proportional-integral (PI) and proportional-integral-derivative (PID) control. Thereafter, the IEEE-33 bus test system and BESS control model were constructed for the scenario study. The final results were obtained by generator power smoothing, and the short-term voltage stability index for the BESS capacity and location was derived from the Manhattan distance method. Subsequently, the daily load and solar power generation curves were considered to verify the simulation results. The schematic diagram of the proposed MATLAB and DIgSILENT simulations is shown in Figure 1. The main contributions of this work are as follows.
1. The sizing and sitting of a BESS in distribution networks with connected renewable energy sources were provided by considering the active power loss, voltage deviation, and total operating cost as the objective functions to minimize; 2. Generator power smoothing and short-term voltage stability index were considered using the transient response as the selection strategy to determine the final capacity and location of the BESS.
Finally, the rest of the paper is organized as follows: In Section 2, we introduce the network component models by DIgSILENT; in Section 3, we present the objective functions and constraints; in Section 4, we describe the proposed algorithm and solution steps; in Section 5, we illustrate the analytical results and comparisons; and in Section 6, we present our conclusions.

Network Component Modeling by DIgSILENT PowerFactory
PowerFactory 15.2 is a software developed by DIgSILENT GmbH, independent software and consulting company based in Gomaringen, Germany, that is suitable for power generation, transmission, distribution, and industrial grid analyses. It combines its system modeling capabilities with advanced algorithm concepts. In this study, the DIgSILENT Simulation Language (DSL) model in PowerFactory 15.2 was applied to construct the BESS model, which included battery and smart inverter models. The schematic diagram is shown in Figure 2.  The main contributions of this work are as follows.

1.
The sizing and sitting of a BESS in distribution networks with connected renewable energy sources were provided by considering the active power loss, voltage deviation, and total operating cost as the objective functions to minimize; 2.
Generator power smoothing and short-term voltage stability index were considered using the transient response as the selection strategy to determine the final capacity and location of the BESS.
Finally, the rest of the paper is organized as follows: In Section 2, we introduce the network component models by DIgSILENT; in Section 3, we present the objective functions and constraints; in Section 4, we describe the proposed algorithm and solution steps; in Section 5, we illustrate the analytical results and comparisons; and in Section 6, we present our conclusions.

Network Component Modeling by DIgSILENT PowerFactory
PowerFactory 15.2 is a software developed by DIgSILENT GmbH, independent software and consulting company based in Gomaringen, Germany, that is suitable for power generation, transmission, distribution, and industrial grid analyses. It combines its system modeling capabilities with advanced algorithm concepts. In this study, the DIgSILENT Simulation Language (DSL) model in PowerFactory 15.2 was applied to construct the BESS model, which included battery and smart inverter models. The schematic diagram is shown in Figure 2.
output, the active power was injected into the system through proportional-integral (PI) and proportional-integral-derivative (PID) control. Thereafter, the IEEE-33 bus test system and BESS control model were constructed for the scenario study. The final results were obtained by generator power smoothing, and the short-term voltage stability index for the BESS capacity and location was derived from the Manhattan distance method. Subsequently, the daily load and solar power generation curves were considered to verify the simulation results. The schematic diagram of the proposed MATLAB and DIgSILENT simulations is shown in Figure 1.  The main contributions of this work are as follows.
1. The sizing and sitting of a BESS in distribution networks with connected renewable energy sources were provided by considering the active power loss, voltage deviation, and total operating cost as the objective functions to minimize; 2. Generator power smoothing and short-term voltage stability index were considered using the transient response as the selection strategy to determine the final capacity and location of the BESS.
Finally, the rest of the paper is organized as follows: In Section 2, we introduce the network component models by DIgSILENT; in Section 3, we present the objective functions and constraints; in Section 4, we describe the proposed algorithm and solution steps; in Section 5, we illustrate the analytical results and comparisons; and in Section 6, we present our conclusions.

Network Component Modeling by DIgSILENT PowerFactory
PowerFactory 15.2 is a software developed by DIgSILENT GmbH, independent software and consulting company based in Gomaringen, Germany, that is suitable for power generation, transmission, distribution, and industrial grid analyses. It combines its system modeling capabilities with advanced algorithm concepts. In this study, the DIgSILENT Simulation Language (DSL) model in PowerFactory 15.2 was applied to construct the BESS model, which included battery and smart inverter models. The schematic diagram is shown in Figure 2.

User Defined
Step 2: Model Construction

Selector
Step 1: Element Selection

Battery Model
Lead-acid batteries are often used in DERs because of their advanced and low-cost technology with the longest development history. The electrode is composed mainly of lead Mathematics 2022, 10, 3420 4 of 25 (Pb) and oxygen (O), and the sulfuric acid solution is used as the battery electrolyte. The BESS is mainly charged and discharged through the oxidation-reduction reactions at the positive and negative poles of the batteries. In the charge state of a lead-acid battery, the main component in the positive plate is lead dioxide, and the negative plate is Pb. While in the discharge state, the main component in the positive and negative plates is lead sulfate, as shown in Figure 3.

Battery Model
Lead-acid batteries are often used in DERs because of their advanced and low-cost technology with the longest development history. The electrode is composed mainly of lead (Pb) and oxygen (O), and the sulfuric acid solution is used as the battery electrolyte. The BESS is mainly charged and discharged through the oxidation-reduction reactions at the positive and negative poles of the batteries. In the charge state of a lead-acid battery, the main component in the positive plate is lead dioxide, and the negative plate is Pb. While in the discharge state, the main component in the positive and negative plates is lead sulfate, as shown in Figure 3.  As shown in Figure 2, the DSL is constructed by selecting the network elements and structuring composite and common models. In this study, the composite model of the battery is created using the current signal module (StaImea), common model module (ElmDSL), and network element module (ElmDcu), and its connections are shown in Figure 4. At the same time, the common model of the battery contains the control strategy model and the initial condition definition for the input and output parameters, and its main control framework is illustrated in Figure 5.   As shown in Figure 2, the DSL is constructed by selecting the network elements and structuring composite and common models. In this study, the composite model of the battery is created using the current signal module (StaImea), common model module (ElmDSL), and network element module (ElmDcu), and its connections are shown in Figure 4. At the same time, the common model of the battery contains the control strategy model and the initial condition definition for the input and output parameters, and its main control framework is illustrated in Figure 5.

Smart Inverter Model
In recent years, BESSs have been integrated rapidly and heavily into power grids with renewable energy resources. BESSs have a significant impact on system voltage and frequency. In response to this situation, research is being conducted on smart inverters. Smart inverter outputs are generated by pulse-width modulation (PWM) based on a synchronized rotary axis composed of two control loops. In this study, a voltage-source converter was used, assuming that the DC side is an ideal voltage source. The switch is controlled by the current inner-loop method, and the outer loop is controlled by the active and reactive power. The d-axis signal controls the active power, and the q-axis signal controls the reactive power by synchronizing the rotary axis coordinates. Furthermore, the architecture of the voltage source inverter controller is shown in Figure 6. In an inner-loop control, the three-phase currents and voltages at the output of the converter are captured. The abc components are then converted into αβ0 components using the Clarke transformation, as shown in Equations (1) and (2).

Smart Inverter Model
In recent years, BESSs have been integrated rapidly and heavily into power grids with renewable energy resources. BESSs have a significant impact on system voltage and frequency. In response to this situation, research is being conducted on smart inverters. Smart inverter outputs are generated by pulse-width modulation (PWM) based on a synchronized rotary axis composed of two control loops. In this study, a voltage-source converter was used, assuming that the DC side is an ideal voltage source. The switch is controlled by the current inner-loop method, and the outer loop is controlled by the active and reactive power. The d-axis signal controls the active power, and the q-axis signal controls the reactive power by synchronizing the rotary axis coordinates. Furthermore, the architecture of the voltage source inverter controller is shown in Figure 6.

Smart Inverter Model
In recent years, BESSs have been integrated rapidly and heavily into power grids with renewable energy resources. BESSs have a significant impact on system voltage and frequency. In response to this situation, research is being conducted on smart inverters. Smart inverter outputs are generated by pulse-width modulation (PWM) based on a synchronized rotary axis composed of two control loops. In this study, a voltage-source converter was used, assuming that the DC side is an ideal voltage source. The switch is controlled by the current inner-loop method, and the outer loop is controlled by the active and reactive power. The d-axis signal controls the active power, and the q-axis signal controls the reactive power by synchronizing the rotary axis coordinates. Furthermore, the architecture of the voltage source inverter controller is shown in Figure 6. In an inner-loop control, the three-phase currents and voltages at the output of the converter are captured. The abc components are then converted into αβ0 components using the Clarke transformation, as shown in Equations (1) and (2). In an inner-loop control, the three-phase currents and voltages at the output of the converter are captured. The abc components are then converted into αβ0 components using the Clarke transformation, as shown in Equations (1) and (2). (1) Mathematics 2022, 10, 3420 Subsequently, the αβ0 components are transformed into dq0 components using the Park transformation, as shown in Equations (3) and (4).
In the case of a fixed frequency, ε is calculated using Equation (5), where ω 0 is the AC frequency, and ε 0 is a constant and is the initial phase angle between the αβ0 and dq0 components.
Moreover, the active and reactive powers in the three-phase coordinate system are given by Equation (6).

BESS Model
The BESS consists of a battery, a smart inverter, a battery management system, and an energy management system. BESS can increase the system voltage, stabilize the system frequency, and improve the transient stability, thus enhancing the reliability of power systems. In this study, the battery was operated in the discharge mode when active power was injected into the power system. The BESS composite model is shown in Figure 7. In the Pf controller, the input is the active power or system frequency, as shown in Figure 8a, whereas, in the QV controller, the input is the reactive power or bus voltage, as shown in Figure 8b. The measurement signal is first calculated by filtering the oscilloscope outputs such that the result is more precise when calculated with the reference value. Subsequently, the error signals are compensated using PID and PI controllers. Finally, the output results are presented as id_ref_in and iq_ref_in, respectively.

Problem Formulations
BESS has become an essential component to support energy transitions in many countries, and Taiwan is no exception. BESS is expected to improve power loss and increase the system bus voltage. However, the implementation of a BESS may be expensive. Therefore, it is crucial to consider system responses and cost factors for BESS installation. Accordingly, the objective functions and system constraints of the BESS used in this study are described as follows.

Control Variables
In order to satisfy the objective function of the planned control variables, the feasible solution interval of the objective function and the restriction condition of the inequality are expressed as min f n (x, u), n ∈ {1, 2, · · · N} (7) Here, f n (x, u) is the nth objective function; g i (x, u) is the ith constraint; x is the control variable; and u is the dependent variable. This study focused on the design of the capacity and location of a BESS to satisfy the designed objective functions. Therefore, the control variable x and the dependence variable u based on MOPSO can be expressed as Equations (10) and (11), respectively.
x = BESS loc,1 , · · · , BESS loc,n , BESS cap,1 , · · · , BESS cap,n T (10) u = V 1 , · · · , V i , V 1 , · · · , V j , P B,1 , · · · , P B,n , P loss,1 , · · · P loss,n T (11) BESS loc,n is the nth location of the BESS; BESS cap,n is the nth capacity of the BESS; V i and V j are the voltages at bus-i and bus-j, respectively; P B,n and P loss,n are the nth active power output and the nth active power loss by the BESS, respectively.

Objective Functions
In this study, the objective functions are (1) active power loss, (2) voltage deviation, and (3) total operating cost.

Active Power Loss
Active power loss is an aspect of the capacity and location optimization problem that cannot be ignored; therefore, the first objective function is defined using Equation (12).
where N b is the number of buses in the power system; n is the node; g ij is the admittance on branch-ij; V i and V j are the voltages at bus-i and bus-j, respectively; and θ ij is the voltage phase displacement on branch-ij.

Voltage Deviation
Voltage deviation is an important indicator for evaluating safety and power quality [22]; hence, it is considered the second objective function. To evaluate the voltage deviation of bus-i, Equation (13) was applied.
where N b is the number of buses in the power system; V rated is the nominal voltage at a bus, which is defined as 1.00 p.u.; and V i is the voltage at bus-i in the power system.

Total Operating Cost
As important as the minimization of active power loss and voltage deviation, the minimization of the total operating cost (TOC) also becomes the objective function, and the cost evaluation was performed using Equation (14) [23].
where k b is the power generation factor of the BESS and k i is the loss cost factor. TOC includes BESS generation cost, as shown in Equation (15), and loss cost, as shown in Equation (16), where P B is the active power output provided by the BESS; a 1 , a 2 , and a 3 are quadratic cost coefficients; and P loss is the active power loss obtained by power flow calculations. Furthermore, K 1 is the annual load cost; K 2 is the price of energy loss; L f is the loss factor evaluated using Equation (17), where k and l f are the loss factor constants.

Constraints
In the following subsection, we explained the implementation of voltage, capacity, and power flow constraints to comply with the power-system operation standards after connecting the BESS to the grid.

Voltage Constraints
The voltages at all buses shall remain within the allowable range shown in Equation (16).
where V min is the lower boundary of the voltage, defined as 0.95 p.u.; V max is the higher boundary of the voltage, defined as 1.05 p.u.; and V i is the voltage magnitude at bus-i.

BESS Technical Constraints
The installation capacity of each BESS depends on its type and operating conditions. The Taiwan Power (Taipower) company classifies BESS capacity according to the voltage level, and the installation capacity shall remain between 100 and 10,000 kW. The operation of the BESS should also be considered in terms of energy constraints. The permissible levels of the BESS are shown in Equations (19) and (20).
P B is the active power output of the BESS; P B,min and P B,max are its boundary limits; and SoC is the battery charging status with SoC min and SoC max as the boundary limits.

Branch Power Flow Constraints
The power provided by the BESS grid connection may increase the power flow in transmission lines. In order to avoid line overloading, the power flowing through the transmission lines must be within the allowable thermal limit shown in Equation (21).
where S line is the power flow through the transmission line and S max,line is the related line thermal capacity.

Proposed Algorithm
Currently, optimization algorithms are widely used in various fields, such as engineering, applied science, and mathematics. Typically, an optimization algorithm is eminent for its systematic search features with a random search mechanism. An optimization algorithm can avoid convergence into regional best solutions and effectively improve the search performance for global best solutions (Gbest). In recent years, many optimization algorithms have been developed, such as the genetic algorithm (GA) [24], particle swarm algorithm (PSO) [25], and the ant colony algorithm [26].

Pareto Optimality
Most practical engineering problems have multiple design objectives. Therefore, it is essential to consider multi-objective conditions simultaneously to determine the best solution for each combination of objective designs. When satisfying all design objectives without simplifying any objective function, some feasible solutions in the solution set are superior to others. If no solutions exist superior to a particular solution in all objective functions, then the solution is called a non-dominated solution. The set of all combinations of non-dominated solutions is called the set of Pareto-optimal solutions. In terms of minimization objectives, domination can be expressed as in Equation (22).
where x for all objectives is better than x because x is the dominant solution and x is the dominated solution. Thus, a non-dominated solution is obtained when no other solution is found to be better than the dominant solution. Then, all non-dominated solutions in the solution set space can be expressed as in Equation (23), where x i is the ith non-dominated solution.

Particle Swarm Optimization (PSO)
In 1995, Kennedy and Eberhart proposed a new type of algorithm inspired by the flocks of birds in nature called the particle swarm algorithm (PSO). In the PSO algorithm, the possible solution positions in the solved space are treated as bird movement spaces. The possibility for better solutions during the solution process is gradually improved by transmitting messages from individuals to each other.
In the PSO algorithm, an individual flies at a certain speed in the search space, and the direction and distance of the flight are controlled by particle velocity. According to the flight experience of the individual and group, the flight speed is dynamically adjusted, and the individual gradually moves to a better area of the problem space. The PSO algorithm is initialized for a group of random particles. By following the present best particle, the search position is continuously iterated to determine the best or most efficient solution. In each iteration, each particle is updated after determining the personal best solution (the present best solution of the particle) and Gbest (the present best solution found by the entire population). The velocity increment during the flight process is closely related to the flight experience of the individual and swarm and is limited by the maximum flight speed, as shown in Figure 9. The MOPSO algorithm, adopted from [27], applies the search capabilities of the PSO algorithm combined with Pareto optimality. A flowchart of the MOPSO algorithm is shown in Figure 10, and the steps of the proposed algorithm are as follows.
Step (1) Initialize parameters and population: set the number of particle swarms, external archive, number of iterations, and other parameters; Step (2) Generate randomly: generate the position and velocity of each particle randomly; Step (3) Select fitness: the fitness of each particle is calculated according to the objective function, and the non-dominated solution is obtained from the particle swarms and stored in the external archive; Step (4) Select Gbest: start the iteration and select a non-dominant solution as Gbest from the external archive; Step (5) Update: use the velocity and position of the particle update formula to obtain the new position and velocity of the particles; Step (6) Evaluate the fitness of the particles: evaluate the fitness of each particle based The MOPSO algorithm, adopted from [27], applies the search capabilities of the PSO algorithm combined with Pareto optimality. A flowchart of the MOPSO algorithm is shown in Figure 10, and the steps of the proposed algorithm are as follows.
Step (1) Initialize parameters and population: set the number of particle swarms, external archive, number of iterations, and other parameters; Step (2) Generate randomly: generate the position and velocity of each particle randomly; Step (3) Select fitness: the fitness of each particle is calculated according to the objective function, and the non-dominated solution is obtained from the particle swarms and stored in the external archive; Step (4) Select Gbest: start the iteration and select a non-dominant solution as Gbest from the external archive; Step (5) Update: use the velocity and position of the particle update formula to obtain the new position and velocity of the particles; Step (6) Evaluate the fitness of the particles: evaluate the fitness of each particle based on the objective function; Step (7) Select non-dominated solutions: find new non-dominated solutions from the population; Step (8) Check the external archive: determine whether the external archive is full or not; Step (9)  Step 4. Store the optimal solution in an external archive and select Gbest Step 5. Update position and velocity Step 6. Evaluate the fitness of particles Step 7. Evaluate nondominated solutions Step 9.

Compare nondominated solutions
Terminal condition?

End
Step 3. Find the optimal solution based on the objective functions Step 1. Set initial parameters Step 2. Randomly generate particle position and velocity Step 10.

Minimum Manhattan Distance Method
In order to find the ideal result, Chiu et al. [28] proposed the minimum Manhattan distance (MMD) to select the final solution among a set of Pareto fronts called multi-criteria decision-making (MCDM). The minimum distance from the normalized ideal vector is called the MMD.
After MOPSO evaluation, many non-dominated solutions in a Pareto front can be obtained, and the minimum value of each objective function is found in the Pareto front. The ideal vector can be calculated by normalizing the minimum value of each objective function, as shown in Equation (24).
The minimum value of the objective function is then calculated using Equation (25).

Minimum Manhattan Distance Method
In order to find the ideal result, Chiu et al. [28] proposed the minimum Manhattan distance (MMD) to select the final solution among a set of Pareto fronts called multi-criteria decision-making (MCDM). The minimum distance from the normalized ideal vector is called the MMD.
After MOPSO evaluation, many non-dominated solutions in a Pareto front can be obtained, and the minimum value of each objective function is found in the Pareto front.
The ideal vector y opt can be calculated by normalizing the minimum value of each objective function, as shown in Equation (24).
The minimum value of the objective function is then calculated using Equation (25).
where ρ N is the minimum value of the Nth objective function in the Pareto front PS A for the set of Pareto optimal solution sets and x i ∈ PS A is the ith solution in the set of all feasible solutions. Furthermore, L N is calculated as the deviation between the maximum and minimum values of each objective function, as shown in Equation (26).
By using the deviation mentioned above, the vector y(x i ) that normalizes all the solutions is represented using Equation (27).
Thus, the minimum Manhattan distance is the shortest distance from the ideal Pareto front among all the solutions, expressed using Equation (28).

Final Solution Selection Index
In the Pareto optimal solution based on the results of the Manhattan distance method, the final capacity and location of the BESS are determined by considering th.

Generator Power Smoothinge Generator Power Smoothing and Short-Term Voltage Stability Index
Solar power can provide additional external generation resources to a power grid and reduce the pressure on power units to supply electricity during the peak period. However, its uncontrollable and intermittent characteristics are problematic. For example, solar power disappears in the evening, and the reduced power generation during this period must be compensated by conventional generators. However, the slow response time of conventional generators may worsen the scheduling difficulty. Therefore, to ensure the safety and stability of the power grid, BESS is particularly important as an ancillary service. In this study, we observed the changes in a conventional generator output power by incorporating BESSs with different capacities and in various locations when the hourly load is increased or decreased, using Equation (29).
where P d,t is the power deviation of the conventional generators at hour t and P s,t and P e,t are the generator power when the load fluctuates and under normal load conditions, respectively, at hour t.

Short-Term Voltage Stability Index
The root-mean-square voltage-dip severity index (RVSI), implemented from [29], was used to evaluate the transient voltage performance for a detailed evaluation of the voltage characteristics after an accidental event. This study introduces an accidental event as a three-phase short-circuit fault in a power system because this event is the most severe case of a short-circuit fault. Short-circuit faults cause a sudden drop in the system voltage, and the fault current flowing through the conductor causes the temperature to rise, causing overheating and equipment damage. Therefore, it is essential to calculate the RSVI using Equations (30) and (31).
T s is the fault start time; T e is the time taken for voltage stabilization after the fault is cleared; V i,t is the voltage magnitude of bus-i at time t; V i,t 0 is the voltage magnitude of bus-i at time t 0 ; t 0 is the initial time; and N b is the number of buses.

Proposed Solution Steps
In this study, the MOPSO algorithm was applied to consider the daily load and solar power generation curves in a region in Taiwan. First, the generator output is used as the first indicator. The final BESS capacity and location are then obtained by simulating the short-term voltage stability index. A flowchart of the solution process is shown in Figure 11, and the steps of the proposed solution are as follows.
Step (1) Build MATLAB benchmark platform: The benchmark-related parameters, such as line impedance, bus voltage magnitude and phase angle, and active and reactive power of the load bus, are built in MATLAB; Step (2) Read parameters and mathematical models: The daily load curve and solar power generation curve parameters are incorporated into the benchmark platform. Active power loss, voltage deviation, and cost functions are constructed; Step (3) Set constraint condition: Set bus voltage, BESS technical, and branch power-flow constraints; Step (4) Initialize: Initialize the state variables and set the particle swarm-related parameters, such as the number of particle swarms, number of iterations, and acceleration factors; Step (5) Run the MOPSO algorithm: After the first iteration, the result of each particle is compared with the result of the previous iteration. If the present result is better, then the personal best solution is updated. Moreover, as the global best solution, if the present Gbest is better than the previous iteration result, then Gbest is updated; Step (6) Determine the convergence condition: If the convergence condition is satisfied, the output Gbest is obtained; otherwise, repeat Step 5 until the simulation converges; Step (7) Run Manhattan distance method: After Gbest is obtained using the MOPSO algorithm, the capacity and location of the BESS are determined using the Manhattan distance method; Step (8) Build the DIgSILENT benchmark platform: The benchmark and BESS components, such as generators, loads, transformers, inverters, and batteries, are built in DIgSILENT; Step (9) Set models and parameters: Daily load and solar power generation curves are inserted into smart inverter and battery control models; Step (10) Select the best neighboring solutions: Select five groups of BESS capacities and locations from the solutions obtained in Step 7; Step (11) First transient screening: Set the output from Step 10 into the DIgSILENT. First, generator power smoothing is evaluated when the load changes, then three sets with better BESS capacities and positions are selected; Step (12) Second transient screening: Set the output from Step 11 into the DIgSILENT. In the event of a three-phase short-circuit fault at each bus in the power system, the final BESS capacity and location are selected by considering the short-term voltage stability index.

Read parameters and models
Step 3. Set constraints Step 4. Initialization Step 5.

Run MOPSO algorithm
Step 7. Run Manhattan distance method

MATLAB DIgSILENT
No Start Transient analysis Step 11. First transient screening Step 8. Build DIgSILENT platform Step 9. Set model and parameters

End
Step 10. Select best neighboring solutions

Yes
Step 12. Second transient screening Step 6. Terminal condition? Figure 11. Flowchart of the proposed solution.

Sample System
In this study, the IEEE 33-bus system was used as a benchmark to verify the effectiveness of the proposed method. The test system consisted of 33 buses, 32 branches, and 32 loads. Furthermore, the rated voltage of the power system was 12.66 kV. The shortcircuit capacity was 100 MVA, and the total real power and reactive power loads were 3.715 MW and 2.3 MVAR, respectively [28]. The single-line diagram of the network system is shown in Figure 12. In addition, a solar power generation system was installed in the test system. Therefore, the solar power generation and load models on each bus were considered with the relevant data provided by Taipower. Finally, the simulation was performed in MATLAB R2019b and DIgSILENT PowerFactory 2021 using a computer with an Intel ® Core(TM)i7-4770 CPU @3.940 GHz, RAM 16 GB, and Windows 10 Pro-64-bit operating system.

Sample System
In this study, the IEEE 33-bus system was used as a benchmark to verify the effectiveness of the proposed method. The test system consisted of 33 buses, 32 branches, and 32 loads. Furthermore, the rated voltage of the power system was 12.66 kV. The short-circuit capacity was 100 MVA, and the total real power and reactive power loads were 3.715 MW and 2.3 MVAR, respectively [28]. The single-line diagram of the network system is shown in Figure 12. In addition, a solar power generation system was installed in the test system. Therefore, the solar power generation and load models on each bus were considered with the relevant data provided by Taipower. Finally, the simulation was performed in MATLAB R2019b and DIgSILENT PowerFactory 2021 using a computer with an Intel ® Core(TM)i7-4770 CPU @3.940 GHz, RAM 16 GB, and Windows 10 Pro-64-bit operating system.

Solar Power Generation Model
The hourly solar power generation output can be obtained using Equation (32).
where , is the solar power generation at hour t, is the rated solar power defined as 1 MW, and is the coefficient of power generation at hour t [30]. Because solar power exists from 6 a.m. to 7 p.m. and there is no sunshine for the rest of the day, the hourly generation coefficient is calculated, as shown in Table 1. The daily solar power generation output, based on the hourly coefficient, is then shown in Figure 13.

Solar Power Generation Model
The hourly solar power generation output can be obtained using Equation (32).
where P PV,t is the solar power generation at hour t, P PV is the rated solar power defined as 1 MW, and µ t is the coefficient of power generation at hour t [30]. Because solar power exists from 6 a.m. to 7 p.m. and there is no sunshine for the rest of the day, the hourly generation coefficient is calculated, as shown in Table 1. The daily solar power generation output, based on the hourly coefficient, is then shown in Figure 13.  Figure 12. Single-line diagram of IEEE 33-bus system.

Solar Power Generation Model
The hourly solar power generation output can be obtained using Equation (32).
where , is the solar power generation at hour t, is the rated solar power defined as 1 MW, and is the coefficient of power generation at hour t [30]. Because solar power exists from 6 a.m. to 7 p.m. and there is no sunshine for the rest of the day, the hourly generation coefficient is calculated, as shown in Table 1. The daily solar power generation output, based on the hourly coefficient, is then shown in Figure 13.

Load Model
The hourly load can be obtained using Equation (33).
where P n,t is the load demand of the nth bus at time t, P n is the rated load demand, and δ t is the load coefficient at time t. As shown in Table 2, the load coefficient gradually increased after 9 a.m. because the demand load continued to increase and gradually decreased at 11 p.m. as the demand load diminished. The daily load curve is shown in Figure 14.

Load Model
The hourly load can be obtained using Equation (33).
where , is the load demand of the nth bus at time t, is the rated load demand, and is the load coefficient at time t. As shown in Table 2, the load coefficient gradually increased after 9 a.m. because the demand load continued to increase and gradually decreased at 11 p.m. as the demand load diminished. The daily load curve is shown in Figure 14.

Setting Parameters
The objective functions used in this study are the active power loss, voltage deviation, and total operating cost. They are minimized by applying the search capability of PSO combined with Pareto optimality. The parameters of the MOPSO used in this study are listed in Table 3, and the related parameters in the multi-objective function are listed in Table 4. The simulations are carried out for one-BESS connection (case 1) and two-BESS connection (case 2).

Setting Parameters
The objective functions used in this study are the active power loss, voltage deviation, and total operating cost. They are minimized by applying the search capability of PSO combined with Pareto optimality. The parameters of the MOPSO used in this study are listed in Table 3, and the related parameters in the multi-objective function are listed in Table 4. The simulations are carried out for one-BESS connection (case 1) and two-BESS connection (case 2).

One BESS Case Study
The final capacity and location of the BESS were derived based on the hourly load and solar power generation results. The hourly power generation and RVSI values are presented in Table 5. During the period 1-7 a.m., the change in power generation was relatively small. Although the power generation is better between 1 and 6 a.m. than between 6 and 7 a.m., the difference is insignificant. RVSI was considered because failure is an abnormal phenomenon. However, RVSI values were better for the rest of the time, and the power generation changed more. Moreover, since the period 6-7 a.m. has a better RVSI than the period 1-6 a.m., the final capacity and location of the BESS were obtained as 409 kW on bus 33 between 6 and 7 a.m.  By considering the period 6-7 a.m., the five best neighboring solutions were obtained using the Manhattan distance method and shown in Table 6. If the result was chosen considering only the shortest distance, which is 0.7886, then the 417 kW BESS on bus 31 was selected. However, by using the strategy proposed in this paper, it was found that this MMD result may not be the best. Therefore, the five best neighboring solutions are necessary to evaluate the generator output and RVSI. Thus, the result with better generator output and RVSI was chosen, which offers the optimum BESS capacity as 409 kW on bus 33. Assume that a three-phase short-circuit fault occurs in a bus. Observing bus 33, where solar power is generated, if a 417 kW BESS on bus 31 is installed, the voltage drops to 0.6791 p.u., as shown in Figure 15a. Otherwise, if a 409 kW BESS on bus 33 is installed, the voltage becomes 0.7101 p.u., as shown in Figure 15b. In addition, in terms of the RVSI indicator, the 409 kW BESS on bus 33 performs better than the 417 kW BESS on bus 31. By considering the period 6-7 a.m., the five best neighboring solutions were obtained using the Manhattan distance method and shown in Table 6. If the result was chosen considering only the shortest distance, which is 0.7886, then the 417 kW BESS on bus 31 was selected. However, by using the strategy proposed in this paper, it was found that this MMD result may not be the best. Therefore, the five best neighboring solutions are necessary to evaluate the generator output and RVSI. Thus, the result with better generator output and RVSI was chosen, which offers the optimum BESS capacity as 409 kW on bus 33. Assume that a three-phase short-circuit fault occurs in a bus. Observing bus 33, where solar power is generated, if a 417 kW BESS on bus 31 is installed, the voltage drops to 0.6791 p.u., as shown in Figure 15a. Otherwise, if a 409 kW BESS on bus 33 is installed, the voltage becomes 0.7101 p.u., as shown in Figure 15b. In addition, in terms of the RVSI indicator, the 409 kW BESS on bus 33 performs better than the 417 kW BESS on bus 31. Based on the results of the MOPSO algorithm, the proposed method combined with the Manhattan distance method chooses a 409 kW BESS on bus 33. A comparison of its results with those obtained using the pure PSO algorithm, which is 111 kW BESS on bus 28, is shown in Table 7. The MOPSO algorithm outperforms the PSO algorithm in each Based on the results of the MOPSO algorithm, the proposed method combined with the Manhattan distance method chooses a 409 kW BESS on bus 33. A comparison of its results with those obtained using the pure PSO algorithm, which is 111 kW BESS on bus 28, is shown in Table 7. The MOPSO algorithm outperforms the PSO algorithm in each index, highlighting the effectiveness of the MOPSO algorithm. Furthermore, the comparisons of the proposed strategy with the MMD results in terms of the generator output and RVSI are shown in Figure 16a,b, respectively. For the generator output, the result from the proposed method is slightly higher than the MMD solution at 6 p.m. The generator output is 0.000016 MW more than the MMD solution, which is insignificant. At 12 p.m., the results of the proposed method allow the unit to produce 0.001976 MW less than the MMD solution. Whereas, for the rest of the time, it is either equal to or lower than the MMD solution. Obviously, the proposed method can determine a solution similar to or even better than MMD in the generator output index. However, for the RVSI, solutions equal to the MMD are found at 2, 8, and 9 p.m. For the rest of the time, a better solution than MMD can be found. For example, at 3 a.m., the proposed method results in an RVSI that allows a total average voltage solution below the MMD of 2.463 p.u. The results demonstrate the effectiveness of the proposed strategy.  Furthermore, the comparisons of the proposed strategy with the MMD results in terms of the generator output and RVSI are shown in Figure 16a,b, respectively. For the generator output, the result from the proposed method is slightly higher than the MMD solution at 6 p.m. The generator output is 0.000016 MW more than the MMD solution, which is insignificant. At 12 p.m., the results of the proposed method allow the unit to produce 0.001976 MW less than the MMD solution. Whereas, for the rest of the time, it is either equal to or lower than the MMD solution. Obviously, the proposed method can determine a solution similar to or even better than MMD in the generator output index. However, for the RVSI, solutions equal to the MMD are found at 2, 8, and 9 p.m. For the rest of the time, a better solution than MMD can be found. For example, at 3 a.m., the proposed method results in an RVSI that allows a total average voltage solution below the MMD of 2.463 p.u. The results demonstrate the effectiveness of the proposed strategy.

Two BESSs Case Study
When connecting two BESSs, the final hourly BESS capacity and location are derived based on the hourly load and solar power generation results. The hourly power generation

Two BESSs Case Study
When connecting two BESSs, the final hourly BESS capacity and location are derived based on the hourly load and solar power generation results. The hourly power generation and RVSI values are listed in Table 8. Between 0 and 8 a.m., the power generation change is relatively small. After considering the RVSI, it was found that the best RVSI is between 2 and 3 a.m., even though the power generation between 3 and 4 a.m. is better. Therefore, the final solution, 453 kW BESS on bus 33 and 237 kW BESS on bus 14, was found between 2 and 3 a.m.  Because the final solution is between 2 and 3 a.m., the five best neighboring solutions were obtained using the Manhattan distance method. The results are shown in Table 9. For the case of connecting two BESSs, if the MMD solution is chosen considering the shortest distance, which is 0.6876, then the selected combination of the two BESSs is 463 kW BESS on bus 33 and 127 kW BESS on bus 15. However, after evaluating the generator output and RVSI, it is evident that this MMD solution may not be the best because even though it results in higher generator output, it does not give the best RVSI. Then, a combination of 453 kW BESS on bus 33 and 237 kW BESS on bus 14 was selected as the final solution because the result is still within the acceptable limits, even though it is slightly more expensive. Assume that a three-phase short-circuit fault occurs in a bus. By observing bus 33, it can be shown that using a combination of 463 kW BESS on bus 33 and 127 kW BESS on bus 15 adjusts the voltage to 0.7148 p.u., as shown in Figure 17a. However, when a combination of 453 kW BESS on bus 33 and 237 kW BESS on bus 14 is installed, the voltage is adjusted to 0.7887 p.u., as shown in Figure 17b. It is proven that when the same fault occurs, installing one more BESS can reduce the voltage amplitude during the fault.
Assume that a three-phase short-circuit fault occurs in a bus. By observing bus 33, it can be shown that using a combination of 463 kW BESS on bus 33 and 127 kW BESS on bus 15 adjusts the voltage to 0.7148 p.u., as shown in Figure 17a. However, when a combination of 453 kW BESS on bus 33 and 237 kW BESS on bus 14 is installed, the voltage is adjusted to 0.7887 p.u., as shown in Figure 17b. It is proven that when the same fault occurs, installing one more BESS can reduce the voltage amplitude during the fault.  Based on the results of the MOPSO algorithm, the proposed method combined with the Manhattan distance method chooses a combination of 453 kW BESS on bus 33 and 237 kW BESS on bus 14. This result was then compared with the results of the pure PSO algorithm, which chooses a combination of 100 kW BESS on bus 32 and 119 kW BESS on bus 16. Subsequently, the comparison of the results is shown in Table 10. The results of the MOPSO algorithm show that regardless of whether one BESS or two BESSs are connected, the search results are better than those of the PSO algorithm. Finally, Figure 18a,b show a comparison between the proposed strategy and MMD solution in terms of generator output and RVSI, respectively. For the generator output, the results of the proposed method are equal to or lower than the solution of the MMD at all times. For example, at 1 a.m., 12 p.m., and 2 p.m., the results of the proposed method are similar to those of the MMD solution, and at 1 p.m., it is possible to make the unit output 0.0046 MW less than the MMD solution. The proposed method can determine a solution similar to or even better than MMD in the generator output index. However, for RVSI, a solution equal to the MMD was found at 1 a.m. For the rest of the time, a better solution than MMD can be found. For example, at 3 a.m., the proposed method results in an RVSI that allows a total average voltage solution below the MMD of 2.5693 p.u. The results show that implementing two BESSs is better than implementing one BESS. However, the proposed strategy still outperforms the MMD solutions in terms of both generator output and RVSI, regardless of whether the power system is connected with one BESS or two BESSs. solution than MMD can be found. For example, at 3 a.m., the proposed method results in an RVSI that allows a total average voltage solution below the MMD of 2.5693 p.u. The results show that implementing two BESSs is better than implementing one BESS. However, the proposed strategy still outperforms the MMD solutions in terms of both generator output and RVSI, regardless of whether the power system is connected with one BESS or two BESSs.

Conclusions
In this study, the capacity and location of BESSs in a distribution network were evaluated using the proposed evaluation indices during the transient response. In order to evaluate the multi-objective problem, the MOPSO algorithm was applied to minimize active power loss, voltage deviation, and total operating cost. Moreover, DIgSILENT was used to build the BESS model and the control strategy. Furthermore, to investigate the performance of the proposed strategy, the proposed method was compared with the MMD and PSO algorithm results. In case studies with one BESS and two BESSs, the results demonstrate the superior performance of the proposed strategy. In addition, in terms of generator output power smoothing and RVSI, the simulation results show that the proposed strategy outperforms the MMD results with better search capability in terms of both the transient response indices. Subsequently, the proposed strategy was verified using the IEEE 33-bus system, demonstrating the excellent performance of the proposed strategy.

Conclusions
In this study, the capacity and location of BESSs in a distribution network were evaluated using the proposed evaluation indices during the transient response. In order to evaluate the multi-objective problem, the MOPSO algorithm was applied to minimize active power loss, voltage deviation, and total operating cost. Moreover, DIgSILENT was used to build the BESS model and the control strategy. Furthermore, to investigate the performance of the proposed strategy, the proposed method was compared with the MMD and PSO algorithm results. In case studies with one BESS and two BESSs, the results demonstrate the superior performance of the proposed strategy. In addition, in terms of generator output power smoothing and RVSI, the simulation results show that the proposed strategy outperforms the MMD results with better search capability in terms of both the transient response indices. Subsequently, the proposed strategy was verified using the IEEE 33-bus system, demonstrating the excellent performance of the proposed strategy.