Biobjective Optimization-Based Frequency Regulation of Power Grids with High-Participated Renewable Energy and Energy Storage Systems

Large-scale renewable energy sources connected to the grid bring new problems and challenges to the automatic generation control (AGC) of the power system. In order to improve the dynamic response performance of AGC, a biobjective of complementary control (BOCC) with high-participation of energy storage resources (ESRs) is established, with the minimization of total power deviation and the minimization of regulation mileage payment. To address this problem, the strength Pareto evolutionary algorithm is employed to quickly acquire a high-quality Pareto front for BOCC. Based on the entropy weight method (EWM), grey target decision-making theory is designed to choose a compromise dispatch scheme that takes both of the operating economy and power quality into account. At last, an extended two-area load frequency control (LFC)model with seven AGC units is taken to verify the effectiveness and the performance of the proposed method.


Introduction
Automatic generation control (AGC) is one of the important tools to maintain the contact line exchange power and realtime network frequency of the power system within the schedule when the grid experiences load disturbances [1]. Traditional AGC units mainly include thermal and hydro units, which are hard to fast track the dynamic power input commands due to their low regulation performance [2]. With the development of renewable energy, a large number of wind power and photovoltaic (PV) units are connected to the grid. On the one hand, as the outputs of wind power and PV units are regulated by power electronics equipment, they can quickly respond to the dynamic power input regulation commands. On the other hand, since the large-scale wind power and PV units are subject to climate conditions, their generation outputs have large random fluctuations, which aggravates the pressure on frequency regulation of power system. In addition, more and more new energy storage resources (ESRs) are joining the grid, such as chemical battery energy storage, electric vehicles, and grid-scale battery storage.
Currently, there is a worldwide effort to balance the intermittency of renewable energy in the grid with highcapacity batteries, in which these ESRs have considerably faster regulation performance than conventional power generators. In China, a vast 200 MW, 800 MWh Vanadium Redox Flow Battery in the Dalian High-Tech zone is currently being readied for operation, which is the largest chemical energy storage plant in the world [3]. In California, the largest battery energy storage project in the world, Gateway Energy Storage, is in progress [4]. e addition of energy storage can improve the economic operation of the system. It can help reduce the pressure on frequency regulation of power system caused by PV units. Wind power and PV units generate electricity by storing excess energy in large-capacity battery sets, which are fed back into the grid when the batteries are not generating electricity to relieve the pressure on frequency regulation of the power system. Besides, the traditional frequency regulation compensation scheme does not effectively stimulate fast response resources to be engaged in AGC because it fails to provide fair compensation in terms of actual regulation performance [5]. e Federal Energy Regulatory Commission (FERC) published Order No. 755 in 2011 to establish a more reasonably priced compensation framework [6]. is command document showed that the compensation payment in the performance-based frequency regulation market is not only determined by the regulation capacity, but also by the actual regulation volume and unit performance [7]. However, some independent system operators (ISOs) implemented 755 negatively and did not establish fair and dynamic markets, because they adopted existing rules that protect existing resources such as liquefied natural gas (LNG) units from competition that is created by ESRs with lower regulation mileage payment and faster regulation performance. And then in February 2018, the FERC issued a landmark order for ESRs, Order No. 841, setting standards for ESR participation in frequency regulation [8].
e regulation units such as ESRs with higher regulation performance will be preferred for ISOs when receiving the same power regulation command. And it will be paid more compared to the regulation unit with lower regulation performance.
So far, few research studies have addressed the complementary control between wind, PV, ESRs, and other frequency regulation resources. In [9], a new biobjective optimization model of real-time AGC dispatch (BOAD) was constructed, with the minimization of the total power deviation and the regulation mileage payment. However, it did not consider ESRs to participate in AGC. Hence, this paper presents a new biobjective of complementary control (BOCC) for AGC with the high participation of ESRs.
To solve the BOCC with high nonlinearity, the metaheuristic-based multiobjective optimization algorithms show stronger global searching ability than the traditional mathematical optimization methods [10][11][12][13][14][15]. Widely used multiobjective algorithms include the improved strength Pareto evolutionary algorithm (SPEA2) [16] and a nondominated sorting genetic algorithm II (NSGA-II) [17]. Compared with NSGA-II, the fitness function of SPEA2 is not calculated directly from the value of the objective function, which is indirectly linked to the value of the objective function but to the Pareto strength. SPEA2 ensures that the converged solution set is a set of noninferior solutions and considers the density of individuals in the solution space to maintain the diversity of the population and ensure the uniform distribution of the solution set. erefore, this paper adopts SPEA2 to solve BOCC.
In addition, since ISOs can allocate only one AGC dispatch signal to each unit at each control interval, the ISOs should select the best suitable solution from the multiple solutions obtained on the Pareto front.
is requires the application of a decision-making method to help the ISOs select a best compromise dispatch scheme from the set of obtained Pareto solutions. Grey target decision-making is one of the methods to solve multiobjective optimization by using the grey systematic theory to select a best compromise solution. In [18], a transformation operator of rewarding the good and punishing the bad was proposed to classify the index of effect sample matrix into three categories. Furthermore, some subjective weighting methods such as the analytic hierarchy process method (AHP) [19] and Delphi method [20] are usually used to determine the weights of indexes by the experts' evaluation matrix, but they easily lead to deviations in indicator weights, while the objective weighting method such as the entropy weight method (EWM) is based on inherent information of decisionmaking matric to determine the weights of indexes. And then grey target decision-making theory with EWM is designed to choose a dispatch plan that takes into account both of the operating economy and power quality, which can provide an objective decision without a subjective judgment. e remaining of this work is organized as follows: Section 2 presents the mathematical model of BOCC. Section 3 gives the detailed implementation of SPEA2 and grey target decision-making with EWM for BOCC. e simulation results and discussions are given in Section 4. Finally, Section 5 concludes the paper.

Mathematical Model of BOCC
2.1. AGC Dispatch Framework. AGC mainly consists of two operations. Figure 1 displays an extended two-area load frequency control (LFC) model. ΔP T is the junction line power exchange deviation; Δf is the real-time frequency deviation; ΔP out is the actual regulated power output; and ΔP D is the power disturbance. e goal of the first operation is to approximate the real-time power disturbance by a PI controller with the inputs of the real-time frequency deviation Δf and the junction line power exchange deviation ΔP T . In the second operation, the ISO assigns the total generation command ΔP C to all the AGC units. Note that this paper focuses on the specific allocation process of the second operation, which is addressed by the proposed SPEA2 and grey target decision-making with EWM. e specific framework of the AGC dispatch model has been given in reference [9], and the ESRs are added to participate in the AGC scheduling process.

Constraints.
In BOCC, it should consider various constraints, including the power balance constraint, dynamic response process with generation ramp constraint (GRC), regulation capacity constraint, and energy transfer constraint [9], as follows: (1) Power balance constraint: at the k th control interval, the total power regulation command output by the controller should be equal to the sum of the power regulation input signals received by all AGC units, as follows: where ΔP in i (k) denotes the input power command received by the i th unit at k th control interval and ΔP c (k) denotes the output of the PI controller.
(2) GRC: depending on the response time delay, AGC units can be sorted into different types of units [21]. Like wind power and PV units, the dynamic response model of ESRs has no a generation ramp constraint (GRC), as shown in Table 1. And the dynamic response function is shown in Figure 2. Without considering the GRC and power limiter, the actual regulated power output is relating to a Laplacian inverse transfer function, as follows: where G i (s) is the energy transfer function of the i th AGC unit; ΔT is the delay time constant of the i th unit; and M in i (k) is the regulation mileage input of the i th AGC unit at k th control interval.
If the GRC and power limiter are considered, then the output of the AGC unit can be calculated, as follows: where ΔP min i and ΔP max i are the minimum and maximum regulation capacities of the i th unit, respectively; R min i and R max i are the minimum and maximum power regulation variations of the i th unit, respectively; and ΔP rate i is the maximum ramp rate of the i th unit.

Objective Function.
Since the proposed BOCC is intended to minimize the total power deviation between the regulation command and the actual power regulation output and minimize the regulation mileage payment, the objective functions can be written as follows: where R i represents the regulation mileage payment for the i th AGC unit, as follows: where λ is the price of regulation mileage; S p i is the performance score; ΔP out i (k) is the actual regulation power output of i th unit at the k th time control interval; and M out i (k) is the regulation mileage output at the kth control interval.

Principle of SPEA2
3.1.1. Optimization Operations. In general, SPEA2 mainly contains four operations, as follows: (1) Initialization: generating a random set of solutions in the solution space X j � ΔP in i (i � 1, 2, . . . , n, j � 1, 2, . . . , N) X j forms the initialized population P 0 , and the size of the population is N. And then creating the empty archive P 0 ′ � ∅, and the size of the archive is N. In addition, at the kth control interval, the lower and upper boundary of all optimization variables are set equally to the lower and upper adjustment capacity of all units, respectively.
(2) Fitness calculation: based on the equation power balance constraint in (1) can be regarded as a known quantity during optimization. By considering the constraint of the d th unit output, if the boundary conditions are violated, the fitness function of j th individual should be given with a large penalty value. Hence, the fitness function can be designed as follows: where ΔP min d and ΔP max d are the minimum and maximum regulation capacities of the d th unit, respectively.
(3) Sorted by actual fitness function and copy: the actual fitness function can be calculated by equations (11)- (15). en, the results can be sorted by adaptation value. When the number of noninferior solutions is less than the number of archive size N, copy the first N − |P t+1 | individuals X j with F(X j ) ≥ 1 from the resulting ordered list to P t+1 [11], and | · | corresponds to the number of elements of a set. (4) Archive truncation procedure: when the number of noninferior solutions exceeds the number of archive Response time delay Governor-turbine Power limiter size N, some of the solutions should be removed from P t+1 until |P t+1 | � N.

Parameters of SPEA2. (1) Strength Value S(X i ).
Denote the number of solutions that j th individual dominates: where X j denotes the i th individual of the nondominated set A; P t denotes the t th generation population; P t denotes the t th generation archived population; | · | corresponds to the number of elements of a set; + represents the sum of two sets and the symbol; and ≻ indicates that the latter is dominated by the former.
(2) Raw Fitness R(X j ). Denote the quality of other solutions that are better than the i th solution. Its lower value means better solution. In particular, when R(X j ) is 0, it means that there is no solution strictly better than this solution in the internal and external populations, and it is a noninferior solution. It is calculated as follows: (3) Density D(x i ). A number less than 1 indicates the density of the solution around the individual, and the larger the value is, the denser the solution set is. When the number of noninferior solutions exceeds the number of archive size N, some of the solutions should be rejected in the dense area. However, the individuals in the sparser regions are retained to evolve in the next generation, so that the diversity of solutions can be maintained: where δ k j means the distance between X j individuals and the k th nearest individuals; N denotes the population size; and N denotes the archive size.
(1) Actual fitness function F(x i ): it can be calculated as follows:

Design of Grey Target Decision-Making with EWM.
In this paper, grey target decision-making based on EWM is used to filter out the compromise solutions of the Pareto solution set. A bullseye is selected in grey target decisionmaking region formed by the solution set. e distance between each solution and the bullseye is taken as an important basis for grey target decision-making, and then the solutions are ranked according to the distance. Each solution of the solution set is considered as a separate decision solution. e weights and distances to the bullseye of each solution obtained by EWM do not rely on the evaluation and preference of experts, so that the decision is credible.

Design of the Effect Sample Matrix.
e Pareto solution set X based on the SPEA2 algorithm is a matrix of n rows and m columns, which is also the set of unit outputs. Here, the absolute values of each solution in X can be taken as one of the decision indicators, as follows: To consider the reduction of total power deviation and regulation mileage payment, two objective function values F 1 and F 2 are used as one of the evaluation indicators.
In this paper, we consider the limiting the variation of the output of each unit by adding an index D, representing the Euclidean distance of each solution of X9 to the origin, as follows: erefore, there are m + 3 evaluation indicators, which are m unit outputs, two objective function values, and Euclidean distance squared D. erefore, the matrix of effect samples containing n decision options and m + 3 decision objectives is expressed as follows:

Design of the Bullseye Vector.
e operator Z j based on the principles of rewarding the best and punishing the worst is calculated as follows: Since the indicators are positive and the smaller they are, the better the program is, so this paper selects the cost-type indicator formula, the decision-making matrix V is calculated as follows: en, the decision matrix can be obtained

Design of Bullseye.
In this paper, EWM is used to calculate the weights of the decision indicators, which is an objective method that does not depend on expert judgments or decision makers' preferences and helps to filter an objective solution from numerous solutions of the Pareto solution set and use this objective solution as the best compromise dispatch scheme. e weight y ij and entropy value E j are calculated based on the index values of each program, respectively, as follows: Mathematical Problems in Engineering According to the bullseye vector v 0 � v 0 1 , v 0 2 , . . . , v 0 m+3 [13], the bullseye distances of each program can be expressed as follows: e principle of the screening program is that the closer the indicator is to the bullseye, the better the solution.

Selection of Decision Options.
When three consecutive AGC time intervals are taken, the set of three Pareto solution sets obtained by the SPEA2 algorithm is allowed to compare with the change in the decision solution, and the results are shown in Figure 3. It can be concluded from Figure 3 that the optimal rational point obtained by applying this decision is not necessarily in the geometric center of the set but is the result of combining the output of each unit as well as the biobjective value.

Calculation Flow.
e whole calculation flow of solving BOCC by SPEA2 and grey target decision-making with EWM is provided in Table 2, where the termination condition of SPEA2 is set to the maximum number of iterations.

Case Studies
In this paper, the SPEA2 and grey target decision-making based BOCC are evaluated via the static test and dynamic test on the extended two-area LFC model by expanding the original single equivalent generating unit in area A into seven AGC units. e control time cycle of AGC is equal to 4 s, and the price of regulation mileage is equal to 2 $/MW. e main parameters of transfer function parameters of the units are given in Table 3. And the main parameters of regulation of the units are given in Table 4. e population size and the maximum number of iteration steps of all the algorithms are set to 50 and 50, respectively. And the archive size of SPEA2 is set to 50.

Static Test Experiments.
In the static test, two load disturbances are applied to the dynamic model, with ΔP D � 50 MW and ΔP D � − 50 MW. Here, 10 algorithms including the multiobjective immune algorithm with nondominated neighbor-based selection (NNIA) [22], NSGA-II, strength Pareto evolutionary algorithm based on reference direction for multiobjective (SPEAR) [23], SPEA2, a decision variable clustering-based evolutionary algorithm for large-scale many-objective (LMEA) [24], a cellular genetic algorithm for multiobjective optimization (MOCELL) [25], a proposal for multiple-objective particle swarm optimization (MOPSO) [26], an evolutionary many-objective optimization algorithm using reference point-based nondominated (NSGA-III) [27], many-objective evolutionary optimization based on reference points (PREA) [28], and region-based selection in evolutionary multiobjective optimization (PESAII) [29] are selected for comparison. Figure 4 provides the comparison of the local Pareto front of the ten algorithms. Figure 4(a) shows that the solutions obtained by LMEA and PESAII deviate obviously from the ideal Pareto front, in which the Pareto front obtained by SPEA2 tends to be better than the fronts of NSGA-II and NSGA-III. Figure 4(b) shows that the solutions obtained by LMEA, SPEAR, and PESAII deviate obviously from the ideal Pareto front, while the Pareto front obtained by SPEA2 tends to be better than the fronts of NSGA-II, MOCELL, and NSGA-III. Meanwhile, the Pareto solutions of SPEA2 are evenly and widely distributed on the local Pareto front. Figure 4(c) shows that Pareto fronts obtained by 10 algorithms with 10 runs are converted into an approximate ideal Pareto surface. To better compare the convergence and diversity of the algorithms, the performance indexes of four algorithms are evaluated in Table 5, including inverted generational distance (IGD) [30], generational distance (GD) [31], pure diversity (PD) [32], hyper volume (HV) [33], DM (diversity metric) [34], spread [35], and spacing [36], where IGD [37] and HV denote the accuracy and diversity of algorithm; GD considers the accuracy of the algorithm; PD and spread denote the diversity of the algorithm; spacing denotes the evenness of the solutions obtained by algorithm; GD, IGD, spread, and spacing are the negative indexes, while DM, HV, and PD are the positive indexes [38][39][40].
It is evident that the SPEA2 outperforms the other three algorithms in the extend two-area LFC model at the two disturbances: (1) e SPEA2 has the smallest GD value among the average data of the algorithm metrics, representing that it has a good merit of convergence performance. (2) It has the largest DM and HV average, which demonstrates that SPEA2 has the largest objective space covered by the Pareto front. Mathematical Problems in Engineering AGC dispatch is less than 4 s, which means that SPEA2 can converge to the Pareto front much faster than other three algorithms.

4.2.1.
Step Load Disturbance. To better verify the superiority of SPEA2 and grey target decision-making, the following dynamic simulation tests are performed in real-time optimization of BOCC. is paper introduces a comparison using the adjustable capacity method, called the proportional (PROP) method [41]. PROP is an engineering calculation method, which distribute the dispatching signal of the unit by using the proportion of the adjustable capacity of the unit. erefore, the output of the i th AGC unit at k th time interval is calculated as follows: en, two load disturbances (see Figures 5 and 6) are applied to the extended two-area LFC model, which are ΔP D � 70 MW and ΔP D � − 50 MW. Figure 5(a) reveals that the proposed SPEA2 with grey target decision-making and EWM can get a more minor error between the total input power and total output power than PROP via the regulation coordination between ESRs and other units, as shown in Figure 5(b). Hence, it slightly reduces the peak of frequency deviation due to it can match the total regulation command closely. Moreover, it results a lower regulation mileage payment than PROP, as illustrated in Figure 7(a). Figure 6(a) shows that compared to the PROP, the proposed method obtains a less power deviation, reduces the overshoot of the total power command, and makes the total       power output curve markedly closer to the total command curve. However, it results a higher regulation mileage payment (see Figure 7(a)), because it assigns larger commands to the AGC units with fast performance scores.

Continuous
Step Load Disturbance. In the following dynamic simulation test, two continuous step load disturbances with 120 control intervals a step and 300 control intervals a step (see Figures 8(b) and 9(c)) are applied to the extended two-area LFC mode.
At the short-time interval step, Figure 8(a) illustrates that the proposed method can also gain a slightly smaller power deviation. After the 460 control interval (see Figure 8(a)), the mileage payment of PROP quickly exceeds that of SPEA2. After reaching the lower limit capacity, the outputs of the PV unit and wind power remain almost static (see Figure 8(b)). However, for PROP, those outputs continue to change dramatically during the simulation.
At the long-time interval step, Figure 9(a) shows that the proposed method also gets a smaller power deviation and a smaller frequency deviation. Moreover, it makes a lower regulation mileage payment (see Figure 7(a)) because PROP distributes more input power command to the PV unit and wind power, while the disturbance borne by wind power and PV unit is shared by ESR units after SPEA2 optimization. It can be concluded that at the short-time interval step, the proposed method can help significantly reduce the mileage payment for IOSs. On the contrary, at the long-time interval step, it can make the total output command well close to the total input command.
To further discuss the relationship between the regulation mileage payment and regulated output power of new energy units and ESRs, Figure 7(a) displays the payment changes under the four disturbance cases. In case of ideal single-step disturbance, due to the stochastic character of the calculation results of the multiobjective algorithm, SPEA2 has uncertainty in reducing the regulation mileage payment, but in the continuous step disturbance case, the algorithm plays a significant role in reducing the regulation mileage payment. In case of successive step disturbance, the two curves before and after the optimization have an intersection point (see Figure 7(a)) due to the changes of the total command assignment scheme (see Figures 7(b) and 7(c)). In the starting phase of the disturbance, because the proposed method assigns mainly units with fast performance scores and high regulation price to participate in the AGC dispatch, it can obtain a slightly higher regulation mileage payment than PROP. However, in the middle stage of disturbance, the PROP still involves wind power and PV unit in regulation frequency, and the regulation capacity of ESR units is not fully utilized. However, in the dispatch with SPEA2 participation, the disturbance borne by wind power and PV units is shared by ESR units, and the superior regulation performance of ESR units is still maximized, so the dispatching payment with SPEA2 participation will be lower than the payment before optimization. In conclusion, ESRs play an important role in AGC dispatch, and the proposed method can make full use of this advantage.
To further test the performance of the SPEA2 algorithm, this example compares the online optimization results of the four cases separately, as shown in Table 6, where accuracy is used to measure the approximate degree of the actual regulation output and the regulation command curve during the simulation time. It can be found two main points from Table 5 compared with PROP, as follows: (1) e proposed method can significantly improve the power quality for power grid by reducing the power deviation and improving each performance index of the system dynamic response, especially in the case of continuous short-time interval step disturbances. e power deviation, average |Δf|, and |ACE| reduce 66.8%, 8.0%, and 5.2%, respectively, in the case of continuous and long-time interval step load disturbance. e power deviation, average |∆f|, and |     continuous long-time interval step disturbances. Particularly, in the two continuous disturbance test, the regulation mileage payments reduce by 22.5% and 9.6%.

Conclusions
In summary, the presented work includes the following three contributions: (1) e proposed BOCC can effectively coordinate various frequency regulation resources and ESRs for AGC, which can simultaneously reduce the total power deviation and the regulation mileage payment.
(2) e presented SPEA2 can converge to a high-quality Pareto front for BOCC in the relatively short time, while grey target decision-making with EWM can effectively select a compromise dispatch scheme that makes full use of the advantage of ESRs in frequency regulation to consider both of the operating economy and power quality.
(3) e simulation results of the extended two-area LFC model verify that the combination of SPEA2 and grey target decision-making with EWM can effectively solve BOCC. It can effectively improve the response performance by reducing the |ACE|, average |Δf|, and total power deviation and simultaneously rise the operating economy for IOSs by reducing the total regulation mileage payment.
To further improve the operating economy and power quality, the future studies will focus on the higher-participated renewable energy and smaller control periodic. Furthermore, the method based on model predictive control or deep learning will be the key technology to solve BOCC.
Data Availability e data that support the findings of this study are available on request from the corresponding author. e data are not publicly available due to privacy or ethical restrictions.