Swarm-based mean-variance mapping optimization for optimal placement of energy storage with synthetic inertia control on a low inertia power grid

Utilizing energy storage systems equipped with virtual/synthetic inertia control has emerged as one of the solutions to solve low inertia issues in power systems with mas-sive penetration of renewables. However, while the aspect regarding the control of the synthetic inertia control has been widely investigated, the optimal placement of synthetic inertia control is rarely investigated. Therefore, the optimal placement of energy storage system with synthetic inertia control based on swarm-based mean-variance mapping optimization is proposed. The performance of the proposed method is tested by considering several cases, both single and multiple contingencies. By utilizing the proposed method, the optimal placement of energy storage system equipped with synthetic inertia control could be obtained, in which the resulting frequency response is better compared to randomly placed ones. The proposed method could also be applied by considering multiple contingencies in different locations. Furthermore, the resulting ﬁtness values for different independent optimization repetitions for various study cases are consistent, leading to the efﬁcacy to be applied


INTRODUCTION
Due to the environmental problem caused by power generation plants based on fossil fuels, in the last few decades, distributed generators based on renewable energy sources (RESs) have been rapidly developed and integrated into the power systems. RESs, most notably photovoltaic (PV) and wind energy gradually increase their share in the energy mix and replacing the traditional power sources, for example, coal and nuclear [1]. It can be seen from the global installation of PV and wind generation which achieves a total of 627 and 650 GW, respectively, by 2019 [2]. However, while the significant increase of RESs penetration is a good indicator of improved RESs utilization as the efforts to reduce the effect of greenhouse gas emissions, it negatively affects the frequency stability of power systems [3][4][5], due to the reduction of overall system inertia caused by the utilization of inverter-based generation. In power systems with a large number of synchronous generators (SGs), This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. © 2021 The Authors. IET Generation, Transmission & Distribution published by John Wiley & Sons Ltd on behalf of The Institution of Engineering and Technology the inertial response is considered an inherent attribute of a power system and hence, was not treated as part of ancillary service. However, due to higher integration of inertialess RESs-based generators utilizing inverter which in turn, replacing the traditional synchronous generators (SGs), the system inertial response reduces. Due to the significant penetration of RESs which result in a significant system inertia reduction, the power system operators started to consider the potential value of the inertial response, which could be provided by using synchronous condensers, wind power plants, and also energy storage systems (ESS) via synthetic/virtual inertia emulation [6].
To enable the synthetic inertia support from the RESs, shortterm energy storage is used, along with an inverter and a proper control system [7]. The concept is referred to as a virtual synchronous generator (VSG) [8] or a virtual synchronous machine [9]. The term synchronverter [10], also refers to the same approach. In this paper, the term "VSG" is used for addressing the aforementioned concepts.
Recently, various research works to augment the frequency stability of power systems by using VSG have been conducted, for example, [3,[11][12][13]. However, most of the research works on VSG focus on the improvement in its control system to achieve a better frequency response. The aspect related to the placement of VSG is rarely investigated. The frequency stability of the power grid is related, not only to the level of inertia in the grid but also to where it is located in the power system [14]. The placement of VSG is also essential in an interconnected power system since it could affect the damping, particularly the interarea mode [15]. In addition, disturbance location and the inertia placement would determine the resiliency of the power system [16].
Thus, the optimal placement of VSG is one of the important aspects of applying VSG. Since the VSG requires energy storage, the problem of optimal VSG placement is closely linked to the energy storage placement. In refs. [17,18], particle swarm optimization (PSO) and whale optimization, respectively, are implemented to optimally determine the placement and sizing of the battery energy storage system (BESS) to reduce power loss. In ref. [19], the method to determine the size and placement of the storage system to support both distribution and transmission networks is proposed. In ref. [20], genetic algorithm is implemented to obtain the location of superconducting magnetic energy storage to enhance voltage stability. In ref. [21], optimal placement of distributed ESS to improve voltage profile and minimize line loading and line losses using an artificial bee colony algorithm is proposed.
However, while the problem of the placement of the ESS in general has been widely investigated, the problem of placement of the ESS equipped with synthetic inertia control (SIC) has not been thoroughly elaborated. In ref. [16], the method for determining the optimal allocation and placement of virtual inertia based on the minimization of the H 2 norm of the power system is proposed. However, this method involves the complexity in the gradient computation [22]. Therefore, in this paper, a novel optimal placement of ESS equipped with SIC to augment the frequency stability of a power system with low inertia using swarm-based mean-variance mapping optimization (MVMOs) [23] is proposed. The advantages of the MVMO include improved performance in terms of convergence speed [24] and the minimum reach [25]. Another advantage lies in the use of normalized search space, helping to ensure precise compliance with the defined boundaries [26]. Using MVMOs, as the swarm-variant of the MVMO, the exploration of the solution space would be performed more aggressively compared to that of MVMO to obtain a better fitness value.
The structure of this paper is described as follows. Section I provides the preface to the discussed topic. Section II explains the employed VSG model. The problem formulation and the proposed optimal placement of VSGs based on MVMOs are elaborated in Section III. In Section IV, the simulation results are presented and discussed. Finally, the conclusion is given in Section V.
The contributions of this paper are: The general concept of virtual synchronous generator (i) Proposing a novel approach for determining the optimal placement of VSG for system frequency enhancement using optimization technique (i.e. the MVMOs). (ii) The candidate locations for VSG placement could be determined by also considering different types and locations of disturbance that might occur in the system. (iii) Using the proposed method, multi-contingencies could also be considered in determining optimal locations of VSG to obtain more appropriate candidate locations of VSG placement in power system planning.

ENERGY STORAGE AND VIRTUAL SYNCHRONOUS GENERATOR MODEL
The concept of a VSG, in general, is depicted in Figure 1. It is formed by energy storage with an inverter controlled by a proper SIC system. The main purpose of the synthetic inertia algorithm is to replicate the inertia property of an SG into the inverter [3].
Typically, a VSG consists of energy storage parts (i.e. the battery cell), a converter interface called power conversion systems (PCS), and the control parts. The model of the BESS utilized in this paper is depicted in Figure 2.
For frequency stability studies performed in this paper, the PCS of the BESS is as depicted in Figure 3(a). It utilizes P-Q decoupled control which enables the independent regulation of active and reactive powers according to their respective setpoint. Typically, the response of a power electronics device is much faster than the dynamics of an SG. Thus, T ES as the time constant of the PQ controller is set as 0.02 s.
The synthetic inertia controller employed in this paper utilizes the derivative technique [11], which adjusts the active power output by using a frequency deviation signal (Δf). Subject to the frequency deviation, the VSG will calculate the required inertia power as follows [11]: where J VI and D VI are the synthetic/virtual inertia and damping constant, respectively. The diagram of this controller is depicted  Figure 3(b). Finally, the power output of the VSG is given as ref. [27]: The flowchart of the mean-variance mapping optimization algorithm [23,29,34]

PROPOSED OPTIMIZATION METHOD
The MVMO is one of the methods of evolutionary optimization algorithms. In MVMO, all optimization variables work within a normalized range and are defined by a repeatedly updated archive that stores the best solutions that have been obtained. Then, a mapping function that accounts for the variance and the mean of the optimization variables is employed for mutation. MVMO projects randomly selected variables onto the related mapping function which will guide the solution following the best solution set that has been obtained [28].

Original mean-variance mapping optimization
The MVMO is one of the population-based optimization techniques [29]. MVMO has been implemented to solve various problems, for example, optimal transmission expansion planning [30], optimal reactive power dispatch [26], placement and tuning of power system stabilizer (PSS) [25] and parameter identification [31].
The flowchart of the MVMO algorithm is presented in Figure 4, with the basic steps are as described below [28,29,[32][33][34]: Step 1: Initialization and Variable Normalization The parameters to be initialized include the initial value of the shape factor, the shaping scaling factor, the number of problem dimensions and dynamic population size.
The search space utilized for all variables is in the range of [0,1]. However, the fitness evaluation is performed using the actual values. Hence, de-normalization is performed in each iteration during the optimization.
Step 2: De-normalize each variable and evaluate its fitness The variables are de-normalized to their original range and then applied to the formulated optimization problem to evaluate its fitness based on the objective function and applied constraints.
Step 3: Solution archive The archive in MVMO stores n number of best child solutions that have been obtained in descending order based on the fitness so that the first ranked solution refers to the best solution that has been found. When the archive is filled up and the fitness of the newly generated solution is better than the ones in the solution archive, an update is performed. During the iterations, the members of the archive keep changing with the improved fitness. The mean x i and variance v i are then examined for each optimization variable x i after every archive update by using (4) and (5), respectively. k is the population size. x i is initialized with the initial x i value, while v i is set as 1 [29].
Then, the shape variable s i is computed as: The shape variable s i , that characterize the mapping function, is calculated by also considering the scaling factor f s , which could regulate the shape of the mapping function and also the search process. The details on the mapping function and the effects of different mean and shape variable values are elaborated in refs. [23,29].
Step 4: Offspring generation The individual with the best fitness f best (i.e. the first ranked position in the archive) and its related optimization values x best are used as the parent of the population for the next evaluation. The parent is used for creating the next offspring.
Among k variables in the optimization problem, m variables are chosen for mutation. The m selected variables are then mutated via the mapping function. The new value of each chosen dimension is determined by using the following equations [29]: in which x * i is randomly varied variable with uniform distribution within [0,1], and h x , h 1 and h 0 are the inputs of the mapping function given by: where h x refers to mapping function, and is described as:

Swarm-based mean-variance mapping optimization
To improve the ability for global search, the swarm-based MVMOs is proposed in ref. [23]. MVMOs has been implemented to solve various problems related to power systems optimization problems, for example, optimal placement and tuning of PSS [35] and optimal reactive power dispatch [36]. Unlike the original MVMO, the swarm type (i.e. the MVMOs) allows for more aggressive exploration of the solution space by initiating the search with a set of particles where each particle has its own defined memory with respect to the correlating archive and mapping function. At first, each particle independently performs m steps to gather a set of individual solutions. After that, the particles begin to exchange information [23,33]. The flowchart of MVMOs is depicted in Figure 5.
To avoid redundancy caused by the overly close distance between particle, the normalized distance between the local best solution of each particle and the global best solution is calculated by [23,36]: where n is the number of considered optimization variables. The i-th particle is removed if the D i is less than a pre-defined threshold D Min .
After an independent evaluation, the global best solution will guide the solution search by assigning x gbest rather than x l best,i as a parent for every offspring. The rest of the steps are similar to the original MVMO. In Figure 5, i and k, respectively, are the counter for function evaluation and particle, while m and Np are the maximum numbers of independent runs and the number of particles.

Objective function
To augment the frequency stability of a power grid, SIC is implemented to minimize the maximum frequency deviation when a significant power imbalance occurs and also to maintain the rate-of-change-of-frequency (RoCoF) at acceptable values.

FIGURE 5
The flowchart of the mean-variance mapping optimizations algorithm [23] Hence, to be able to obtain the desired purposes (i.e. minimizing frequency deviation and maintaining RoCoF within limits) by adjusting the location and the SIC parameters, the optimization problem in this paper utilizes two objective functions which are formulated as follows. loc and Par are the location of VSG and its parameter, respectively.
In Equation (11), the first objective function is defined as a function of frequency nadir/zenith, while the second objective function is related to RoCoF.
In Equations (11)-(13), f nom is nominal frequency, f nadir∕zenith is frequency nadir/zenith subject to the contingencies in low inertia condition, RoCoF refer to the RoCoF (in Hz/s), w 1 is weighting factor for Obj 1 , w 2 is weighting factor for Obj 2 , and x represents the decision variables vector consisting of activated VSG locations and their parameters.
By utilizing the formulated objective function, it is expected that the obtained optimal placement locations of the VSG are the location that resulting in the minimum frequency deviation while maintaining RoCoF within limits.
The optimization constraints in this study grouped into three parts: 1. Frequency deviation and RoCoF constraints.
2. The number of ESS units to be activated in the test system.
3. BESS operational and capacity constraints.
J VI,min ≤ J VI,i ≤ J VI,max.
where f nadir,min is the minimum acceptable frequency nadir (set as 49 Hz), f zenith,max is the maximum acceptable frequency zenith (set as 51 Hz), RoCoF max is the maximum tolerable RoCoF value (set as 0.6 Hz s −1 ), N bus is the number of buses of the test system, P ES,min and P ES,max are the lower and upper

Implementation in DIgSILENT PowerFactory
The simulations in this paper are performed by using DIgSI-LENT PowerFactory. To implement the swarm-based MVMOs to obtain the optimal placement of ESS equipped with SIC, DIgSILENT Programming Language (DPL) is used. The general procedure for the implementation of MVMOs in DIgSI-LENT PowerFactory is shown in Figure 6. DPL is utilized to perform the search procedures in the MVMOs, that is, initialization, normalization-denormalization, archiving n-best solution, parent assignment and offspring generation. From the performed optimization, certain optimization values would be obtained. These obtained values are then used in the dynamic time-domain simulation and the resulting system frequency response is evaluated. The fitness values are calculated by scanning frequency nadir at each bus, along with the RoCoF. With this simulation-based optimization framework, the nonlinearities of the power system model can be taken into account in determining the optimal placement of ESS equipped with SIC.

RESULTS AND DISCUSSION
The proposed method is validated by using a variant of the PST 16-machine test system in ref. [37]. The transmission network uses voltages of 380, 220 and 110 kV. The single-line diagram of the test system is provided in Figure 7. In this study, to represent low inertia conditions, the inertia constant of SGs in the system is halved from the initial values. Each applied VSG will have a rating of 100 MW. The simulations and optimization process in this paper are conducted by using DIgSILENT PowerFactory.
There are two optimization scenarios carried out in this paper. The first scenario only considers one contingency, which consists of one case with a contingency resulting in the overfrequency (Case 1) and another case with a contingency resulting in the underfrequency (Case 2). In the second scenario (Case 3), several contingencies are considered to obtain more appropriate placement locations of VSG.

Placement results for generator outage case (Case 1)
The first scenario to validate the optimal placement of VSGs in this study considers only one contingency. The first study case is the N-1 contingency, in which a switch event is applied at t = 0 s on bus A1b to simulate the 1000 MW loss of generation. The MVMOs is executed for 10 independent repetitions. The number of installed VSGs is set as three units. The resulting optimal placement using the MVMOs is shown in Table 1 and the convergence curve is shown in Figure 8. From the results of ten optimization repetitions, three buses (i.e. A2, A3 and A6) are obtained as the candidates for the VSG location.

Optimal placement using particle swarm optimization (Case 1)
The results obtained from MVMOs are also compared with PSO to verify the result obtained from MVMOs. PSO is chosen as a comparison method due to its similarities with MVMOs, in which both methods use a set of particles (swarm). For a fair comparison, the PSO is stochastically run for ten times with a population size set as 10 particles and also 500 iterations. The result is presented in Table 2 and the convergence curve is shown in Figure 9. From the results of the    optimization, the obtained VSG location candidates using PSO are the same as the MVMOs (i.e. bus A2, A3 and A6). The same results obtained from a well-known PSO algorithm verifies the validity and efficacy of the MVMOs to solve the optimal VSG placement problem. Furthermore, comparing the result obtained from MVMOs and PSO in Figures 8 and 9, respectively, it can be seen that the fitness value obtained by using MVMOs from each FIGURE 10 The frequency response of the test system (Case 1) under: -high inertia condition (blue) -low inertia condition without virtual synchronous generators (red) -low inertia condition with optimally placed virtual synchronous generators (green) -low inertia condition with randomly placed virtual synchronous generators (purple) independent run is more consistent compared to that of PSO's. Hence, in terms of the optimization result, the obtained optimal placement from MVMOs is more robust compared to PSO.
The VSGs are then placed at the optimal locations and the frequency response subject to the loss of 1000 MW generation is analyzed. The frequency response of the test system under conditions of high inertia (blue line), low inertia without VSGs (red line), low inertia with optimally placed VSGs (green line) and low inertia with randomly placed VSGs (purple line) is depicted in Figure 10.
It can be seen from the graph that VSG utilization can augment the frequency stability of the low inertia power grid, where the frequency deviation of systems equipped with VSG is smaller with lower RoCoF, compared to those without VSG. Figure 10 also shows the comparison of the frequency responses between the optimally placed VSGs using MVMOs and the randomly placed VSGs. The results in Figure 10 shows that VSG placement at the optimal locations will result in lower frequency deviation (i.e. higher frequency nadir) and also lower RoCoF compared with VSG placement at random locations. Therefore, the system frequency response subject to the disturbance could be improved by optimally place the VSGs in the system.

Placement results for step load case (Case 2)
In the second case, the MVMOs is tested in the case of overfrequency. A switch event is applied at bus B10 to simulate the loss of the largest load (800 MW). The optimization using MVMOs is then executed for ten independent repetitions. The result is shown in Table 3 and the convergence curve is shown in Figure 11. Based on the results, different optimal locations are obtained compared to Case 1 (i.e. loss of generation). From this result, the determination of the placement location of VSGs should consider various scenarios. However, from the optimization results, bus A3 emerges as a strong candidate for the placement location of VSG since it appears as one of the candidates in both Case 1 and Case 2.

Placement results for different number of virtual synchronous generators
In the previous simulation, the number of VSGs in each optimization process is always set equal to 3. In this section, the effect of the different numbers of installed VSGs and their optimal placement locations are investigated. The optimization processes are performed considering the contingency in Case 1. The results are depicted in Table 4. These results indicate that bus A3 has always been a location candidate for VSG placement. It also appears from this result that a higher number of VSGs will provide a smaller fitness value, which means it will produce a smaller frequency deviation and a smaller RoCoF for the same type of disturbance.
Meanwhile, the comparison of the frequency response for optimized VSG locations with different VSG numbers for Case

FIGURE 12
The frequency response for optimized virtual synchronous generator locations with different virtual synchronous generator numbers 1 is depicted in Figure 12. From Figure 12, it is evident that the installation of more VSGs results in better frequency response, as shown by lower frequency deviation (i.e. higher frequency nadir).

Placement results for different weighting factor
In the previous calculations, the weighting factor w 1 and w 2 is set equal to 1 and 0.5, respectively. It means that the calculation takes more into account the frequency deviation compared to the RoCoF. To see the effects of changes in weighting factors on the results of the VSG placement, the optimization problem for Case 1 is then recalculated with variations in w 1 and w 2 considering the same contribution from frequency deviation and RoCoF (w 1 = w 2 = 0.5) and with more emphasis on RoCoF rather than frequency deviation (w 1 = 0.5,w 2 = 1). The results are shown in Table 5. From the optimization results with the variations above, the same candidate locations are obtained for Case 1 (i.e. bus A2, A3 and A6), either for calculations with the same weighting factor or when a greater concern is given to one of the objective functions. Thus, the optimal placement locations of VSG obtained by using MVMOs are less affected by the selection of the weighting factors.

Results for multi-scenario optimization
In practice, the power system is essentially time-varying, not only depending on a single operating point. Thus, in realistic power planning procedure, the robustness of the placement results need to be considered against different types of scenario variations, for example, different location of disturbance or when the seasonal loading level changes. Then, the solution from the previous single-scenario-based optimization might not be optimal or feasible anymore. In this scenario (Case 3), the case involving different locations of disturbance is considered as an example of the multi-scenario optimization applied in the optimization problem. This scenario uses a 250 MW mobile load injected at one of the pre-selected buses. The pre-selected buses are all 380 kV buses. Then, the process advances to the subsequent bus. The fitness value is examined for each contingency set. The optimizations are then performed for five independent repetitions with the number of particles = 10 and maximum iteration = 250. The placement results are shown in Table 6 and the convergence curve is depicted in Figure 13.
From these five optimization repetitions, there are three buses (A1, A2 and A3) that always appear as candidates for the VSG location. The average fitness value is 0.1021, with a standard deviation of 1.6787 × 10 -5 . The average calculation time for optimal VSG placement with 250 iterations is 1647 minutes and the average number of iterations needed to converge is 115. Based on the consistency of the obtained fitness value from different optimization repetitions, the MVMOs could be implemented as a reliable method to determine the placement locations of VSG considering multi-scenarios.
By taking into account various types of contingency locations that might occur in the power grid, it is expected that a more The convergence curve for mean-variance mapping optimizations for multi-scenario optimization reliable objective value and more credible candidate locations could be obtained to be considered as energy storage locations with SIC to augment frequency stability of the power system. The optimal locations of VSG would be particularly important when the considered contingencies are large contingencies, in which the benefit of the optimal placement of VSGs would be more significant compared to the cases of smaller contingencies.

CONCLUSION
Based on the performed simulations, the optimal VSG placement problem can be effectively solved using the swarm-based MVMO (MVMOs). The obtained optimized placement of VSG (i.e. BESS with SIC) can improve the overall frequency stability in the simulated system by better suppressing both frequency deviation and RoCoF compared to the randomly placed ones. The results of the optimal placement locations of VSGs obtained from MVMOs are generally similar to those obtained from PSO. However, the fitness value obtained from MVMOs is more consistent compared to that of PSO. Hence, the MVMOs is viable to solve the optimal placement problem of VSG.
From the results of the performed simulations, the VSG placement is important in applying VSG to augment the frequency stability of a power system with low system inertia. The results of optimal placement in two cases of a single contingency scenario show that the obtained optimal location for VSG placements is influenced by the type and location of the disturbance being taken into consideration. Hence, it is important to consider various system events to obtain more appropriate locations for VSG placement.
To be used in realistic power planning procedures, the robustness of the placement results needs to be considered against various types of scenario variations. Therefore, the proposed optimization framework is also designed and has been successfully applied to a multi-contingencies scenario, which can lead to more credible locations candidate and more reliable objective values.