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

Authors: Tingyi He; Shengnan Li; Shuijun Wu; Chuangzhi Li; Biao Xu
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

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.


---

## Body

## 1. Introduction

Automatic generation control (AGC) is one of the important tools to maintain the contact line exchange power and real-time 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 high-capacity 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]. The 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]. The Federal Energy Regulatory Commission (FERC) published Order No. 755 in 2011 to establish a more reasonably priced compensation framework [6]. This 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]. The 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–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. Therefore, 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. This 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 decision-making 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.The remaining of this work is organized as follows: Section2 presents the mathematical model of BOCC. Section 3 gives the detailed implementation of SPEA2 and grey target decision-making with EWM for BOCC. The simulation results and discussions are given in Section 4. Finally, Section 5 concludes the paper.


## 2. Mathematical Model of BOCC

### 2.1. AGC Dispatch Framework

AGC mainly consists of two operations. Figure1 displays an extended two-area load frequency control (LFC) model. ΔPT is the junction line power exchange deviation; Δf is the real-time frequency deviation; ΔPout is the actual regulated power output; and ΔPD is the power disturbance. The 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 ΔPT. In the second operation, the ISO assigns the total generation command ΔPC 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. The 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.Figure 1
Framework of AGC on the extended two-area LFC model.


### 2.2. 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 thekth 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:(1)∑i=1nΔPiink−ΔPck=0,whereΔPiink denotes the input power command received by the ith unit at kth control interval and ΔPck 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:(2)ΔPioutt=L−1Giss1+Tdis∑k=1Ne−ΔTk−1sMiink,(3)ΔPioutk=ΔPioutt=k⋅ΔT,(4)Miink=ΔPiink−Piink−1,where Gis is the energy transfer function of the ith AGC unit; ΔT is the delay time constant of the ith unit; and Miink is the regulation mileage input of the ith AGC unit at kth control interval.Figure 2
Dynamic response models: (a) conventional units; (b) renewable units and ESRs.
(a)(b)Table 1
Types of transfer functions for various units.
TypeTransfer functionGsNon-reheat steam unit1/1+T1sReheat steam unit1+T2s/1+T3s1+T4s1+T5sHydro unit1−T6s1+T7s/1+0.5T6s1+T8sWT and PV1/1+T9sESRs1/1+T10sIf the GRC and power limiter are considered, then the output of the AGC unit can be calculated, as follows:(5)ΔPioutk=ΔPioutk−1+Rimin,ifΔPioutk<Rimin,ΔPioutk,ifRimin≤ΔPioutk≤Rimax,ΔPioutk−1+Rimax,ifΔPioutk>Rimin,(6)Rimin=0,ifMck≥0,max−ΔPirate⋅ΔT,ΔPimin−ΔPioutk−1,ifMck<0,(7)Rimax=minΔPirate⋅ΔT,ΔPimax−ΔPioutk−1,ifMck≥0,0,ifMck<0,where ΔPimin and ΔPimax are the minimum and maximum regulation capacities of the ith unit, respectively; Rimin and Rimax are the minimum and maximum power regulation variations of the ith unit, respectively; and ΔPirate is the maximum ramp rate of the ith unit.


### 2.3. 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:(8)minf1=∑k=1NΔPck−∑i=1nΔPioutk+1,minf2=∑i=1nRi,where Ri represents the regulation mileage payment for the ith AGC unit, as follows:(9)Ri=∑k=1NγSiPMioutk,(10)Mioutk=ΔPioutk−ΔPioutk−1,where λ is the price of regulation mileage; Sip is the performance score; ΔPioutk is the actual regulation power output of ith unit at the kth time control interval; and Mioutk is the regulation mileage output at the kth control interval.


## 2.1. AGC Dispatch Framework

AGC mainly consists of two operations. Figure1 displays an extended two-area load frequency control (LFC) model. ΔPT is the junction line power exchange deviation; Δf is the real-time frequency deviation; ΔPout is the actual regulated power output; and ΔPD is the power disturbance. The 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 ΔPT. In the second operation, the ISO assigns the total generation command ΔPC 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. The 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.Figure 1
Framework of AGC on the extended two-area LFC model.


## 2.2. 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 thekth 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:(1)∑i=1nΔPiink−ΔPck=0,whereΔPiink denotes the input power command received by the ith unit at kth control interval and ΔPck 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:(2)ΔPioutt=L−1Giss1+Tdis∑k=1Ne−ΔTk−1sMiink,(3)ΔPioutk=ΔPioutt=k⋅ΔT,(4)Miink=ΔPiink−Piink−1,where Gis is the energy transfer function of the ith AGC unit; ΔT is the delay time constant of the ith unit; and Miink is the regulation mileage input of the ith AGC unit at kth control interval.Figure 2
Dynamic response models: (a) conventional units; (b) renewable units and ESRs.
(a)(b)Table 1
Types of transfer functions for various units.
TypeTransfer functionGsNon-reheat steam unit1/1+T1sReheat steam unit1+T2s/1+T3s1+T4s1+T5sHydro unit1−T6s1+T7s/1+0.5T6s1+T8sWT and PV1/1+T9sESRs1/1+T10sIf the GRC and power limiter are considered, then the output of the AGC unit can be calculated, as follows:(5)ΔPioutk=ΔPioutk−1+Rimin,ifΔPioutk<Rimin,ΔPioutk,ifRimin≤ΔPioutk≤Rimax,ΔPioutk−1+Rimax,ifΔPioutk>Rimin,(6)Rimin=0,ifMck≥0,max−ΔPirate⋅ΔT,ΔPimin−ΔPioutk−1,ifMck<0,(7)Rimax=minΔPirate⋅ΔT,ΔPimax−ΔPioutk−1,ifMck≥0,0,ifMck<0,where ΔPimin and ΔPimax are the minimum and maximum regulation capacities of the ith unit, respectively; Rimin and Rimax are the minimum and maximum power regulation variations of the ith unit, respectively; and ΔPirate is the maximum ramp rate of the ith unit.


## 2.3. 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:(8)minf1=∑k=1NΔPck−∑i=1nΔPioutk+1,minf2=∑i=1nRi,where Ri represents the regulation mileage payment for the ith AGC unit, as follows:(9)Ri=∑k=1NγSiPMioutk,(10)Mioutk=ΔPioutk−ΔPioutk−1,where λ is the price of regulation mileage; Sip is the performance score; ΔPioutk is the actual regulation power output of ith unit at the kth time control interval; and Mioutk is the regulation mileage output at the kth control interval.


## 3. Design of SPEA2 and Grey Target Decision-Making for BOCC

### 3.1. 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 spaceXj=ΔPiini=1,2,…,n,j=1,2,…,NXj forms the initialized population P0, and the size of the population is N. And then creating the empty archive P0′=∅, 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), the variable xd=ΔPck−∑i=1i≠dnΔPiink can be regarded as a known quantity during optimization. By considering the constraint of the dth unit output, if the boundary conditions are violated, the fitness function of jth individual should be given with a large penalty value. Hence, the fitness function can be designed as follows:(11)F2Xj=F1Xj=xd−ΔPdmaxxd−ΔPdmin+1⋅108,ifxd>ΔPdmaxorxd<ΔPdmin,F1Xj=f1Xj,F2Xj=f2Xj,else,whereΔPdmin and ΔPdmax are the minimum and maximum regulation capacities of the dth unit, respectively.(3)
Sorted by actual fitness function and copy: the actual fitness function can be calculated by equations (11)–(15). Then, 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 Xj with FXj≥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 sizeN¯, some of the solutions should be removed from P¯t+1 until P¯t+1=N¯.


#### 3.1.2. Parameters of SPEA2

(1) Strength ValueSXi. Denote the number of solutions that jth individual dominates:(12)SXj=Xi|Xi∈Pt+P¯t∧Xj≻Xi,where Xj denotes the ith individual of the nondominated set A; Pt denotes the tth generation population; P¯t denotes the tth 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 FitnessRXj. Denote the quality of other solutions that are better than the ith solution. Its lower value means better solution. In particular, when RXj 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:(13)RXj=∑Xi∈Pt+P¯t,Xj≻XiSXi.(3) DensityDxi. 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:(14)DXj=1δjk+2,(15)k=N+N¯,where δjk means the distance between Xj individuals and the kth nearest individuals; N denotes the population size; and N¯ denotes the archive size.(1) Actual fitness functionFxi: it can be calculated as follows:(16)Fxi=Rxi+Dxi.


### 3.2. 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 decision-making region formed by the solution set. The 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. The 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.


#### 3.2.1. Design of the Effect Sample Matrix

The Pareto solution setX 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:(17)X′i,j=Xi,j,i=1,2,..,n,j=1,2,…,N.To consider the reduction of total power deviation and regulation mileage payment, two objective function valuesF1 and F2 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 indexD, representing the Euclidean distance of each solution of X′ to the origin, as follows:(18)Di=∑j=1mX′i,j2.Therefore, there arem + 3 evaluation indicators, which are m unit outputs, two objective function values, and Euclidean distance squared D. Therefore, the matrix of effect samples containing n decision options and m + 3 decision objectives is expressed as follows:(19)X″=X′F1F2D.


#### 3.2.2. Design of the Bullseye Vector

The operatorZj based on the principles of rewarding the best and punishing the worst is calculated as follows:(20)Zj=1n∑i=1nX″i,j,j=1,2,…,m+3.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 matrixV is calculated as follows:(21)vij=zij−xijmaxmax1≤i≤nxij−zj,zj−min1≤i≤nxij.Then, the decision matrix can be obtainedV=vijn×m+3. vj0=maxvji|1≤i≤n,j=1,2,….,m+3。Therefore, the selected bullseye vector is v0=v10,v20,…,vm+30.


#### 3.2.3. 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. The weightyij and entropy value Ej are calculated based on the index values of each program, respectively, as follows:(22)yij=xij∑i=1nxij,xij≥0,(23)Ej=−1lnn∑i=1nyijlnyij,Ej>0,(24)ωj=1−Ej∑j=1m+31−Ej.According to the bullseye vectorv0=v10,v20,…,vm+30 [13], the bullseye distances of each program can be expressed as follows:(25)di=vi−v0=∑j=1m+3ωjvij−vj02.The principle of the screening program is that the closer the indicator is to the bullseye, the better the solution.


#### 3.2.4. 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 Figure3. 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.Figure 3
The selection of decision options.


### 3.3. Calculation Flow

The whole calculation flow of solving BOCC by SPEA2 and grey target decision-making with EWM is provided in Table2, where the termination condition of SPEA2 is set to the maximum number of iterations.Table 2
The execution procedure of SPEA2 for BOCC.
Input: The total power regulation command for the current control interval.1: Initialize the parameters of SPEA2;2: Initialize the parent population and archive of SPEA2;3: Calculate the initial dispatch solution by equations (1)–(7);4: Compute the objective function for each individual by equations (8)–(10);5: Evaluate the fitness function for each individual by equation (11);6: Merge the generation population and the generation archived population into one.7: Calculate the parameters of the new set by equations (12)–(16);8: Compare the number of noninferior solutions and archive size, and then perform the truncation or reproduction to fill the archive set.9: Generate the offspring population by applying recombination and mutation operators to the set10: If the termination condition of SPEA2 is satisfied, then go to the next step; otherwise, return to Step 4;11: Find the best compromise solution from multiple solutions in Pareto solution set by grey target decision-making theory in equation (17)–(25).Output: The best power regulation inputs for all AGC units.


## 3.1. 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 spaceXj=ΔPiini=1,2,…,n,j=1,2,…,NXj forms the initialized population P0, and the size of the population is N. And then creating the empty archive P0′=∅, 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), the variable xd=ΔPck−∑i=1i≠dnΔPiink can be regarded as a known quantity during optimization. By considering the constraint of the dth unit output, if the boundary conditions are violated, the fitness function of jth individual should be given with a large penalty value. Hence, the fitness function can be designed as follows:(11)F2Xj=F1Xj=xd−ΔPdmaxxd−ΔPdmin+1⋅108,ifxd>ΔPdmaxorxd<ΔPdmin,F1Xj=f1Xj,F2Xj=f2Xj,else,whereΔPdmin and ΔPdmax are the minimum and maximum regulation capacities of the dth unit, respectively.(3)
Sorted by actual fitness function and copy: the actual fitness function can be calculated by equations (11)–(15). Then, 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 Xj with FXj≥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 sizeN¯, some of the solutions should be removed from P¯t+1 until P¯t+1=N¯.


### 3.1.2. Parameters of SPEA2

(1) Strength ValueSXi. Denote the number of solutions that jth individual dominates:(12)SXj=Xi|Xi∈Pt+P¯t∧Xj≻Xi,where Xj denotes the ith individual of the nondominated set A; Pt denotes the tth generation population; P¯t denotes the tth 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 FitnessRXj. Denote the quality of other solutions that are better than the ith solution. Its lower value means better solution. In particular, when RXj 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:(13)RXj=∑Xi∈Pt+P¯t,Xj≻XiSXi.(3) DensityDxi. 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:(14)DXj=1δjk+2,(15)k=N+N¯,where δjk means the distance between Xj individuals and the kth nearest individuals; N denotes the population size; and N¯ denotes the archive size.(1) Actual fitness functionFxi: it can be calculated as follows:(16)Fxi=Rxi+Dxi.


## 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 spaceXj=ΔPiini=1,2,…,n,j=1,2,…,NXj forms the initialized population P0, and the size of the population is N. And then creating the empty archive P0′=∅, 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), the variable xd=ΔPck−∑i=1i≠dnΔPiink can be regarded as a known quantity during optimization. By considering the constraint of the dth unit output, if the boundary conditions are violated, the fitness function of jth individual should be given with a large penalty value. Hence, the fitness function can be designed as follows:(11)F2Xj=F1Xj=xd−ΔPdmaxxd−ΔPdmin+1⋅108,ifxd>ΔPdmaxorxd<ΔPdmin,F1Xj=f1Xj,F2Xj=f2Xj,else,whereΔPdmin and ΔPdmax are the minimum and maximum regulation capacities of the dth unit, respectively.(3)
Sorted by actual fitness function and copy: the actual fitness function can be calculated by equations (11)–(15). Then, 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 Xj with FXj≥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 sizeN¯, some of the solutions should be removed from P¯t+1 until P¯t+1=N¯.


## 3.1.2. Parameters of SPEA2

(1) Strength ValueSXi. Denote the number of solutions that jth individual dominates:(12)SXj=Xi|Xi∈Pt+P¯t∧Xj≻Xi,where Xj denotes the ith individual of the nondominated set A; Pt denotes the tth generation population; P¯t denotes the tth 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 FitnessRXj. Denote the quality of other solutions that are better than the ith solution. Its lower value means better solution. In particular, when RXj 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:(13)RXj=∑Xi∈Pt+P¯t,Xj≻XiSXi.(3) DensityDxi. 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:(14)DXj=1δjk+2,(15)k=N+N¯,where δjk means the distance between Xj individuals and the kth nearest individuals; N denotes the population size; and N¯ denotes the archive size.(1) Actual fitness functionFxi: it can be calculated as follows:(16)Fxi=Rxi+Dxi.


## 3.2. 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 decision-making region formed by the solution set. The 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. The 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.


### 3.2.1. Design of the Effect Sample Matrix

The Pareto solution setX 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:(17)X′i,j=Xi,j,i=1,2,..,n,j=1,2,…,N.To consider the reduction of total power deviation and regulation mileage payment, two objective function valuesF1 and F2 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 indexD, representing the Euclidean distance of each solution of X′ to the origin, as follows:(18)Di=∑j=1mX′i,j2.Therefore, there arem + 3 evaluation indicators, which are m unit outputs, two objective function values, and Euclidean distance squared D. Therefore, the matrix of effect samples containing n decision options and m + 3 decision objectives is expressed as follows:(19)X″=X′F1F2D.


### 3.2.2. Design of the Bullseye Vector

The operatorZj based on the principles of rewarding the best and punishing the worst is calculated as follows:(20)Zj=1n∑i=1nX″i,j,j=1,2,…,m+3.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 matrixV is calculated as follows:(21)vij=zij−xijmaxmax1≤i≤nxij−zj,zj−min1≤i≤nxij.Then, the decision matrix can be obtainedV=vijn×m+3. vj0=maxvji|1≤i≤n,j=1,2,….,m+3。Therefore, the selected bullseye vector is v0=v10,v20,…,vm+30.


### 3.2.3. 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. The weightyij and entropy value Ej are calculated based on the index values of each program, respectively, as follows:(22)yij=xij∑i=1nxij,xij≥0,(23)Ej=−1lnn∑i=1nyijlnyij,Ej>0,(24)ωj=1−Ej∑j=1m+31−Ej.According to the bullseye vectorv0=v10,v20,…,vm+30 [13], the bullseye distances of each program can be expressed as follows:(25)di=vi−v0=∑j=1m+3ωjvij−vj02.The principle of the screening program is that the closer the indicator is to the bullseye, the better the solution.


### 3.2.4. 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 Figure3. 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.Figure 3
The selection of decision options.


## 3.2.1. Design of the Effect Sample Matrix

The Pareto solution setX 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:(17)X′i,j=Xi,j,i=1,2,..,n,j=1,2,…,N.To consider the reduction of total power deviation and regulation mileage payment, two objective function valuesF1 and F2 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 indexD, representing the Euclidean distance of each solution of X′ to the origin, as follows:(18)Di=∑j=1mX′i,j2.Therefore, there arem + 3 evaluation indicators, which are m unit outputs, two objective function values, and Euclidean distance squared D. Therefore, the matrix of effect samples containing n decision options and m + 3 decision objectives is expressed as follows:(19)X″=X′F1F2D.


## 3.2.2. Design of the Bullseye Vector

The operatorZj based on the principles of rewarding the best and punishing the worst is calculated as follows:(20)Zj=1n∑i=1nX″i,j,j=1,2,…,m+3.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 matrixV is calculated as follows:(21)vij=zij−xijmaxmax1≤i≤nxij−zj,zj−min1≤i≤nxij.Then, the decision matrix can be obtainedV=vijn×m+3. vj0=maxvji|1≤i≤n,j=1,2,….,m+3。Therefore, the selected bullseye vector is v0=v10,v20,…,vm+30.


## 3.2.3. 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. The weightyij and entropy value Ej are calculated based on the index values of each program, respectively, as follows:(22)yij=xij∑i=1nxij,xij≥0,(23)Ej=−1lnn∑i=1nyijlnyij,Ej>0,(24)ωj=1−Ej∑j=1m+31−Ej.According to the bullseye vectorv0=v10,v20,…,vm+30 [13], the bullseye distances of each program can be expressed as follows:(25)di=vi−v0=∑j=1m+3ωjvij−vj02.The principle of the screening program is that the closer the indicator is to the bullseye, the better the solution.


## 3.2.4. 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 Figure3. 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.Figure 3
The selection of decision options.


## 3.3. Calculation Flow

The whole calculation flow of solving BOCC by SPEA2 and grey target decision-making with EWM is provided in Table2, where the termination condition of SPEA2 is set to the maximum number of iterations.Table 2
The execution procedure of SPEA2 for BOCC.
Input: The total power regulation command for the current control interval.1: Initialize the parameters of SPEA2;2: Initialize the parent population and archive of SPEA2;3: Calculate the initial dispatch solution by equations (1)–(7);4: Compute the objective function for each individual by equations (8)–(10);5: Evaluate the fitness function for each individual by equation (11);6: Merge the generation population and the generation archived population into one.7: Calculate the parameters of the new set by equations (12)–(16);8: Compare the number of noninferior solutions and archive size, and then perform the truncation or reproduction to fill the archive set.9: Generate the offspring population by applying recombination and mutation operators to the set10: If the termination condition of SPEA2 is satisfied, then go to the next step; otherwise, return to Step 4;11: Find the best compromise solution from multiple solutions in Pareto solution set by grey target decision-making theory in equation (17)–(25).Output: The best power regulation inputs for all AGC units.


## 4. 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. The control time cycle of AGC is equal to 4 s, and the price of regulation mileage is equal to 2 $/MW. The main parameters of transfer function parameters of the units are given in Table3. And the main parameters of regulation of the units are given in Table 4. The 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.Table 3
Transfer function parameters of AGC units.
Generation typeParameters (s)Coal-firedT2 = 5, T3 = 0.08, T4 = 10, T5 = 0.3LNGT2 = 2, T3 = 0.05, T4 = 5, T5 = 0.2HydroT6 = 1, T7 = 5, T8 = 0.513WT, PVT1 = 0.01ESRsT1 = 2Table 4
Main parameters of AGC units in area A of the two-area LFC model.
Unit no.TypeTd (s)ΔPrate (MW/min)ΔPmax (MW)ΔPmin (MW)G1Coal-fired603050−50G2LNG201830−30G3Hydro515020−10G4WT1—15−5G5PV1—10−10G6ESRs11—8−12G7ESRs21—7−13


### 4.1. Static Test Experiments

In the static test, two load disturbances are applied to the dynamic model, withΔPD=50MW and ΔPD=−50MW. 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.Figure4 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
Comparison of the Pareto front surface obtained by 10 algorithms on the extended two-area LFC model: (a)ΔPD=70MW, (b)ΔPD=70MW, and (c) select the approximate ideal Pareto surface.
(a)(b)(c)Figure4(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–40].Table 5
Comparison of performance metrics of algorithms.
ΔPDFunctionIGDGDPDHVDMSpreadSpacingT (s)50 WNNIAAve8.501.013.36E + 050.5720.7260.6636.047.33E − 02Std2.960.297.18E + 040.0040.0590.0851.061.34E − 03NSGAIIAve11.230.942.99E + 050.5710.6840.6964.966.50E − 02Std5.160.265.21E + 040.0080.0470.0640.692.58E − 03SPAR2Ave11.080.642.85E + 050.5730.7520.4753.626.49E − 02Std5.330.177.60E + 040.0070.0760.0900.791.92E − 03SPEARAve12.200.982.47E + 050.5620.6340.8269.016.63E − 02Std4.170.343.47E + 040.0060.0570.0762.062.98E − 03−50 WNNIAAve12.540.662.18E + 050.4850.6890.7355.447.79E − 02Std11.870.152.47E + 040.0030.0650.0701.235.13E − 03NSGAIIAve14.680.372.15E + 050.4850.6630.7594.736.48E − 02Std11.370.143.38E + 040.0070.0790.0561.052.87E − 03SPAR2Ave22.350.451.74E + 050.4790.6920.5232.966.47E − 02Std25.490.113.52E + 040.0090.1510.1430.872.46E − 03SPEARAve13.120.612.15E + 050.4870.6120.8797.766.70E − 02Std6.040.154.79E + 040.0040.0620.0672.302.84E − 03It is evident that the SPEA2 outperforms the other three algorithms in the extend two-area LFC model at the two disturbances:(1)
The 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.(3)
It has the smallest spread and spacing average, i.e., which indicates that the Pareto solutions obtained by SPEA2 are evenly and widely distributed on the Pareto front.(4)
It also has the smallest running average time and satisfies the requirement that the time interval of 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. Dynamic Test

#### 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. This 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. Therefore, the output of the ith AGC unit at kth time interval is calculated as follows:(26)ΔPioutk=ΔPck⋅ΔPimax∑i=1nΔPimax,ifΔPck≥0,ΔPck⋅ΔPimin∑i=1nΔPimin,ifΔPck<0.Then, two load disturbances (see Figures5 and 6) are applied to the extended two-area LFC model, which are ΔPD=70MW and ΔPD=−50MW.Figure 5
Real-time optimization results obtained on the extended two-area LFC model whenΔPD=70MW, (a) overall power deviation, (b) regulation power output obtained by SPEA2 with grey target decision-making, and (c) frequency deviation.
(a)(b)(c)Figure 6
Real-time optimization results obtained on the extended two-area LFC model whenΔPD=−50MW, (a) overall power deviation, (b) regulation power output obtained by SPEA2 with grey target decision-making, and (c) frequency deviation.
(a)(b)(c)Figure5(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 7
(a) The payment under four disturbance cases. Comparison of regulated output power of new energy units and ESRs before and after optimization: (b) under a continuous and short-time interval step load disturbance and (c) under a continuous and long-time interval step load disturbance.
(a)(b)(c)Figure6(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.


#### 4.2.2. 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 Figures8(b) and 9(c)) are applied to the extended two-area LFC mode.Figure 8
Online optimization results obtained by SPEA2 under a continuous and short-time interval step load disturbance in the extended two-area LFC model, (a) overall power deviation, (b) regulation power output obtained by SPEA2 with grey target decision-making, and (c) frequency deviation.
(a)(b)(c)Figure 9
Online optimization results obtained by SPEA2 under a continuous and long-time interval step load disturbance in the extended two-area LFC model, (a) overall power deviation, (b) regulation power output obtained by SPEA2 with grey target decision-making, and (c) frequency deviation.
(a)(b)(c)At the short-time interval step, Figure8(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, Figure9(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, Figure7(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 Table6, 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)
The 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. The 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. The power deviation, average |∆f|, and |ACE| reduce 45.9%, 5.9%, and 6.4%, respectively, in the case of continuous and short-time interval step load disturbance.(2)
The proposed method can obviously reduce the regulation mileage payment and thus rise the economy operating for ISOs, especially in the case of continuous long-time interval step disturbances. Particularly, in the two continuous disturbance test, the regulation mileage payments reduce by 22.5% and 9.6%.Table 6
Result comparison of online optimization under different disturbances.
ΔPDMethod|ACE| (MW)Δf (Hz)CPS1 (%)Deviation (MW)Accuracy (%)Payment ($)AvgMaxAvgMaxAvgMin70 MWPROP0.8911.174.78E − 045.73E − 03199.99199.82550.0681.34165.58SPEA20.8811.064.52E − 045.67E − 03199.99199.82383.3681.87130.09−50 MWPROP0.648.093.46E − 044.15E − 03199.99199.90426.1481.1997.87SPEA20.616.933.04E − 043.58E − 03199.99199.9353.3682.87126.95Disturb 1 (see Figure8(a))PROP4.7211.332.38E − 035.81E − 03199.96199.812854.3874.23816.35SPEA24.429.952.24E − 035.12E − 03199.96199.851544.5782.73632.56Disturb 2 (see Figure9(a))PROP1.928.098.71E − 044.15E − 03199.99199.901240.6970.31311.12SPEA21.826.958.01E − 043.59E − 03199.99199.93411.5678.75281.25


## 4.1. Static Test Experiments

In the static test, two load disturbances are applied to the dynamic model, withΔPD=50MW and ΔPD=−50MW. 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.Figure4 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
Comparison of the Pareto front surface obtained by 10 algorithms on the extended two-area LFC model: (a)ΔPD=70MW, (b)ΔPD=70MW, and (c) select the approximate ideal Pareto surface.
(a)(b)(c)Figure4(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–40].Table 5
Comparison of performance metrics of algorithms.
ΔPDFunctionIGDGDPDHVDMSpreadSpacingT (s)50 WNNIAAve8.501.013.36E + 050.5720.7260.6636.047.33E − 02Std2.960.297.18E + 040.0040.0590.0851.061.34E − 03NSGAIIAve11.230.942.99E + 050.5710.6840.6964.966.50E − 02Std5.160.265.21E + 040.0080.0470.0640.692.58E − 03SPAR2Ave11.080.642.85E + 050.5730.7520.4753.626.49E − 02Std5.330.177.60E + 040.0070.0760.0900.791.92E − 03SPEARAve12.200.982.47E + 050.5620.6340.8269.016.63E − 02Std4.170.343.47E + 040.0060.0570.0762.062.98E − 03−50 WNNIAAve12.540.662.18E + 050.4850.6890.7355.447.79E − 02Std11.870.152.47E + 040.0030.0650.0701.235.13E − 03NSGAIIAve14.680.372.15E + 050.4850.6630.7594.736.48E − 02Std11.370.143.38E + 040.0070.0790.0561.052.87E − 03SPAR2Ave22.350.451.74E + 050.4790.6920.5232.966.47E − 02Std25.490.113.52E + 040.0090.1510.1430.872.46E − 03SPEARAve13.120.612.15E + 050.4870.6120.8797.766.70E − 02Std6.040.154.79E + 040.0040.0620.0672.302.84E − 03It is evident that the SPEA2 outperforms the other three algorithms in the extend two-area LFC model at the two disturbances:(1)
The 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.(3)
It has the smallest spread and spacing average, i.e., which indicates that the Pareto solutions obtained by SPEA2 are evenly and widely distributed on the Pareto front.(4)
It also has the smallest running average time and satisfies the requirement that the time interval of 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. Dynamic Test

### 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. This 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. Therefore, the output of the ith AGC unit at kth time interval is calculated as follows:(26)ΔPioutk=ΔPck⋅ΔPimax∑i=1nΔPimax,ifΔPck≥0,ΔPck⋅ΔPimin∑i=1nΔPimin,ifΔPck<0.Then, two load disturbances (see Figures5 and 6) are applied to the extended two-area LFC model, which are ΔPD=70MW and ΔPD=−50MW.Figure 5
Real-time optimization results obtained on the extended two-area LFC model whenΔPD=70MW, (a) overall power deviation, (b) regulation power output obtained by SPEA2 with grey target decision-making, and (c) frequency deviation.
(a)(b)(c)Figure 6
Real-time optimization results obtained on the extended two-area LFC model whenΔPD=−50MW, (a) overall power deviation, (b) regulation power output obtained by SPEA2 with grey target decision-making, and (c) frequency deviation.
(a)(b)(c)Figure5(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 7
(a) The payment under four disturbance cases. Comparison of regulated output power of new energy units and ESRs before and after optimization: (b) under a continuous and short-time interval step load disturbance and (c) under a continuous and long-time interval step load disturbance.
(a)(b)(c)Figure6(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.


### 4.2.2. 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 Figures8(b) and 9(c)) are applied to the extended two-area LFC mode.Figure 8
Online optimization results obtained by SPEA2 under a continuous and short-time interval step load disturbance in the extended two-area LFC model, (a) overall power deviation, (b) regulation power output obtained by SPEA2 with grey target decision-making, and (c) frequency deviation.
(a)(b)(c)Figure 9
Online optimization results obtained by SPEA2 under a continuous and long-time interval step load disturbance in the extended two-area LFC model, (a) overall power deviation, (b) regulation power output obtained by SPEA2 with grey target decision-making, and (c) frequency deviation.
(a)(b)(c)At the short-time interval step, Figure8(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, Figure9(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, Figure7(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 Table6, 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)
The 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. The 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. The power deviation, average |∆f|, and |ACE| reduce 45.9%, 5.9%, and 6.4%, respectively, in the case of continuous and short-time interval step load disturbance.(2)
The proposed method can obviously reduce the regulation mileage payment and thus rise the economy operating for ISOs, especially in the case of continuous long-time interval step disturbances. Particularly, in the two continuous disturbance test, the regulation mileage payments reduce by 22.5% and 9.6%.Table 6
Result comparison of online optimization under different disturbances.
ΔPDMethod|ACE| (MW)Δf (Hz)CPS1 (%)Deviation (MW)Accuracy (%)Payment ($)AvgMaxAvgMaxAvgMin70 MWPROP0.8911.174.78E − 045.73E − 03199.99199.82550.0681.34165.58SPEA20.8811.064.52E − 045.67E − 03199.99199.82383.3681.87130.09−50 MWPROP0.648.093.46E − 044.15E − 03199.99199.90426.1481.1997.87SPEA20.616.933.04E − 043.58E − 03199.99199.9353.3682.87126.95Disturb 1 (see Figure8(a))PROP4.7211.332.38E − 035.81E − 03199.96199.812854.3874.23816.35SPEA24.429.952.24E − 035.12E − 03199.96199.851544.5782.73632.56Disturb 2 (see Figure9(a))PROP1.928.098.71E − 044.15E − 03199.99199.901240.6970.31311.12SPEA21.826.958.01E − 043.59E − 03199.99199.93411.5678.75281.25


## 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. This 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. Therefore, the output of the ith AGC unit at kth time interval is calculated as follows:(26)ΔPioutk=ΔPck⋅ΔPimax∑i=1nΔPimax,ifΔPck≥0,ΔPck⋅ΔPimin∑i=1nΔPimin,ifΔPck<0.Then, two load disturbances (see Figures5 and 6) are applied to the extended two-area LFC model, which are ΔPD=70MW and ΔPD=−50MW.Figure 5
Real-time optimization results obtained on the extended two-area LFC model whenΔPD=70MW, (a) overall power deviation, (b) regulation power output obtained by SPEA2 with grey target decision-making, and (c) frequency deviation.
(a)(b)(c)Figure 6
Real-time optimization results obtained on the extended two-area LFC model whenΔPD=−50MW, (a) overall power deviation, (b) regulation power output obtained by SPEA2 with grey target decision-making, and (c) frequency deviation.
(a)(b)(c)Figure5(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 7
(a) The payment under four disturbance cases. Comparison of regulated output power of new energy units and ESRs before and after optimization: (b) under a continuous and short-time interval step load disturbance and (c) under a continuous and long-time interval step load disturbance.
(a)(b)(c)Figure6(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.


## 4.2.2. 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 Figures8(b) and 9(c)) are applied to the extended two-area LFC mode.Figure 8
Online optimization results obtained by SPEA2 under a continuous and short-time interval step load disturbance in the extended two-area LFC model, (a) overall power deviation, (b) regulation power output obtained by SPEA2 with grey target decision-making, and (c) frequency deviation.
(a)(b)(c)Figure 9
Online optimization results obtained by SPEA2 under a continuous and long-time interval step load disturbance in the extended two-area LFC model, (a) overall power deviation, (b) regulation power output obtained by SPEA2 with grey target decision-making, and (c) frequency deviation.
(a)(b)(c)At the short-time interval step, Figure8(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, Figure9(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, Figure7(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 Table6, 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)
The 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. The 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. The power deviation, average |∆f|, and |ACE| reduce 45.9%, 5.9%, and 6.4%, respectively, in the case of continuous and short-time interval step load disturbance.(2)
The proposed method can obviously reduce the regulation mileage payment and thus rise the economy operating for ISOs, especially in the case of continuous long-time interval step disturbances. Particularly, in the two continuous disturbance test, the regulation mileage payments reduce by 22.5% and 9.6%.Table 6
Result comparison of online optimization under different disturbances.
ΔPDMethod|ACE| (MW)Δf (Hz)CPS1 (%)Deviation (MW)Accuracy (%)Payment ($)AvgMaxAvgMaxAvgMin70 MWPROP0.8911.174.78E − 045.73E − 03199.99199.82550.0681.34165.58SPEA20.8811.064.52E − 045.67E − 03199.99199.82383.3681.87130.09−50 MWPROP0.648.093.46E − 044.15E − 03199.99199.90426.1481.1997.87SPEA20.616.933.04E − 043.58E − 03199.99199.9353.3682.87126.95Disturb 1 (see Figure8(a))PROP4.7211.332.38E − 035.81E − 03199.96199.812854.3874.23816.35SPEA24.429.952.24E − 035.12E − 03199.96199.851544.5782.73632.56Disturb 2 (see Figure9(a))PROP1.928.098.71E − 044.15E − 03199.99199.901240.6970.31311.12SPEA21.826.958.01E − 043.59E − 03199.99199.93411.5678.75281.25


## 5. Conclusions

In summary, the presented work includes the following three contributions:(1)
The 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)
The 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)
The 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 theACE, 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.



---

*Source: 5526492-2021-03-01.xml*