Dynamic Model Reduction for Large-Scale Power Systems Using Wide-Area Measurements

To perform faster than the real-time dynamic simulation of large-scale power systems, it is necessary to reduce the simulated system size by using equivalents for surrounding areas of the study area, and existing dynamic model reduction approach could provide the needed structure of the reduced area. However, further parameter optimization is required to achieve the desired accuracy. In this paper, a particle swarm optimization (PSO) based approach is used to solve the above problem. Parameters for the individual dynamic elements in the reduced system are calibrated repeatedly until the wide-area measurements of the reduced model and the original model are very similar to each other with satisfactory accuracy. Results indicate that after optimization, the dynamic response of the reduced model matches better with that of the original one than using existing methods. Under both the generator-trip event and the bus-fault event, the reduced model has a higher frequency match and less power mismatch.


I. INTRODUCTION
The dynamic response is of great significance for the analysis of a large-scale power system. However, it is computationally demanding to conduct such type of analysis using an original model, where each element is modeled in fine detail [1]. In fact, only a certain part of the power system is taken as the main focus of the study, defined as the study area. The rest part of the system is called the external area, whose structure and internal parameters are not the main concern. To speed up the simulation process, the external area should be replaced by a simplified one. On the other hand, the reduced model should hold similar dynamic characteristics with the original one when disturbances occur in the study area.
In recent years, many schemes and algorithms have aimed to reach the above goal. Toward this end, some have proposed several state-of-the-art approaches using the modalbased method. In [2], the principal component analysis is employed to linearize the original model, which serves as the The associate editor coordinating the review of this manuscript and approving it for publication was Lei Chen . first step to reduce the model using modal-based methods. Later in [3], a novel model reduction algorithm based on an extension of balanced truncation is proposed, where the active and the reactive power of the boundary is represented by a nonlinear equation of the bus voltage and the phase angle. In [4], a validation of the Krylov subspace is attempted on reducing the linearized model of the external area, by simulating a fault in the study area. Similar contributions based on the linear transfer-function [1], [5], the autoregressive model [6], etc. have shown great performance in projecting the relationship between the study area and the external area. However, the major bottleneck of the above method is the computation complexity. Also, the accuracy of the simulation results, under a different type of disturbance to the training case, needs further enhancement as the reduction error may substantially increase [7] if the training case fails to cover all types of disturbances.
Coherency-based methods provide a promising way to overcome the above shortfalls, and this type of method is the most widely-used recently. To reduce a system using this type of method, there are three steps to be implemented: the VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ identification of the coherency group, the system equivalence, and the dynamic reduction of the system. This work includes the aggregation of the network [8] and aggregations of the generators [9], [10]. At present, many types of commercial software, such as the DYNERD [11], rely on the coherencybased equivalent approach. According to [12], the coherencybased method is capable of providing a homogeneous reduced model that has similar dynamic response characteristics with the original model. Although this type of reduction approach could provide the needed structure of the reduced area, there is still space for enhancement as aggregation algorithms in this paper do not include higher-order machine models and control systems. As a result, some critical parameters for the generator may not be accurate enough. In this situation, these parameters are necessarily calibrated to ensure that the reduced system can keep similar eigenvalues and eigenvectors with the original one.
Recent developments in large-scale power systems have spawned interest in the use of the wide-area monitoring system (WAMS) utilizing phasor measurement units (PMUs). The applications of WAMS includes 1) power system monitoring, such as the post-event monitoring [13] and the state estimation [14]; 2) power system protection, such as the enhanced back-up relays [15]- [17], the adaptive out-of-step protection [18], etc.; and 3) power system control, such as the transient stability prediction [19], intelligent control strategies based on support vector machines [20], etc. The extensive application of the WAMS has been a powerful key to providing a technical approach for improving the security level of large-scale power systems. Also, the increasing deployment of PMUs is developed for better full-system observability with a large number of available data sources, thereby being the basis for the wide-area measurement-based dynamic model reduction methods.
The contribution of this paper is to develop a dynamic model reduction method by repeatedly tuning the parameters of the reduced model until the dynamic response of the reduced system matches the original wide-area measurements obtained in large-scale power systems. Conventionally, numerical methods, such as the nonlinear least squares (NLS), have proved to be powerful tools for solving problems like parameter identification [21], frequency estimation [22], stability analysis [23], etc. However, as the system is highly nonlinear, it is quite challenging to find a group of optimized parameters using numerical methods. Among several types of intelligent algorithms, particle swarm optimization (PSO) is widely applied in the power system for its excellent performance solving problems like parameter identification [24], nonconvex economic dispatch [25], etc. In comparison with the genetic algorithm (GA) that has already been successfully applied in the power system [26], PSO is more efficient because it has a smaller number of internal parameters and a shorter time to converge. This makes the PSO-based method more suitable for real-time applications.
The remainder of this paper is organized as follows. Section II introduces how to obtain a reduced structure for the external area, which is based on the previous contribution of [9]. Section III discusses how the PSO works to tune the parameters of the external area. Section IV presents the application of the proposed method to the Northeast Power Coordinating Council (NPCC) on which the merits of the proposed method are evaluated. Finally, the contributions are concluded in Section V.

II. PARAMETER TUNING OF THE EXTERNAL AREA
A. THE MODEL OF SECOND-ORDER GENERATORS Literature [9] presents an aggregation algorithm for obtaining reduced order power networks, which has already been implemented in the EPRI dynamic equivalence program DYNRED. In this method, the aggregated generators are using the second-order model, as presented below: where U , I , and E are the phasors for the terminal voltage, terminal current, and the excitation electrical motive force (EMF), respectively. R a and X d are the stator resistance and the d-axis transient reactance. H J is the inertia. D is the damping factor. T m and T e are the mechanical torque and the electromagnetic torque, respectively. ω and δ are the rotor angular frequency and the power transfer angle, respectively. According to the above, if the reduced model includes a total number of i second-order generators, there should be a total of 2i parameters to be tuned.

B. THE THEORY OF PSO
The core philosophy of the PSO is to find the best particle in the solution space. Each particle is assumed to be traveling in this space, and its location and velocity are renewed in each generation according to the desired fitness function. The following variables are defined as below: where d is the dimension of the solution space. X i represents the location of the particle i. Similarly, V i represents the velocity of particle i given by: The fitness of the particle i is given by (4). For a problem of minimization, equation (4) will have better fitness when the objective function holds a low value.
Let P i be the location of particle i when it has the best fitness, and P g be the best location of all particles in the population, i.e., the global best fitness, which are given by: Similar to other types of optimization algorithms, in each generation, the location of each particle is determined according to equation (7), which is given by: where g represents the order of generation. To renew the velocity and the location of particles, equation (8) and (9) are employed in each generation: where c 1 and c 2 are two acceleration factors, which perform as the step sizes to renew the location and velocity of each particle. r 1 and r 2 are two random numbers that uniformly distributed between 0 and 1. ω is the dynamic inertia weight. j represents the j-th dimension among 1 to d.

III. METHODOLOGY
After the application of the method in [9] to the original model, the tie lines between the study area and the external area can be left unchanged in the reduced model. For this reason, given the same disturbance that occurred in the study area, the similarity of the measured frequency on each tie line in the reduced model to that of the original model can be used as an indicator to represent the model accuracy.
As shown in Fig. 1, assume that there are a total number of T tie lines between the two areas. If the total number of disturbance scenarios that should be considered is S, we assume that the measured frequency on these tie lines are f r. 11

A. DETERMINING THE FITNESS FUNCTION
Under the same type of scenarios s, the similarity of the frequency of a certain tie line t is calculated by (10): (10) where f r.st is the mean of f r.st . F s.t is the fitness, which represents the similarity of the frequency measurement on the t-th tie line under s-th scenario. If this index is closer to 1, it indicates that the time domain response of the reduced model is in higher accordance with that of the original model, in which case both the variation characteristics (oscillation mode) and the magnitude error of the reduced model are quite close to those of the original model.
To take the frequency measurements on all tie lines into consideration, the average fitness in scenario s is given by: Thus, the fitness function of the PSO that to be optimized is given by: As indicated, the optimizing equation (12) is to find the F with a minimum value, which is equal to finding a group of particles that make the frequencies of tie lines in the reduced model matches best with the original model.

B. PARAMETER SETTING
PSO uses c 1 and c 2 to control the step sizes of the renewing of the location and the velocity of each particle. To have better efficiency, the higher step size is preferred. However, if the step size is set too high, the tuning process will be difficult to reach a convergence result. According to [27], these two parameters are set to 2.0.
The setting of ω is to obtain a balanced capability between finding a global minimum and finding a local minimum. The higher ω is, the better global convergence capability the PSO will have. Vice versa, the PSO will hold a higher local convergence capability. Therefore, at the beginning of the iteration, a higher ω is preferred; and ω should decrease with the growth of the generation. Thus, a dynamic ω is set as given by: where ω max and ω min are the maximum and the minimum of ω, which are conventionally set to 0.9 and 0.2. G is the generation size, and g represents the number of the current generation. σ is the attenuation index to determine the decline rate of ω, which is set to 1 in this study. VOLUME 8, 2020 To prevent the particles from running out of the searching space, the velocity of each particle should have some limitations, which is given by: Thus, the velocity and location of each particle are renewed according to equations (8), (9), and (14).

C. DYNAMIC MODEL REDUCTION
For illustrative purposes, the process to reduce the model is given below, as shown in Fig. 2. The details are as follows.
1) Obtaining the structure of the reduced model according to [9]. Thus, the external area should include several numbers of generators, and each of them is both modeled in secondorder. The total number of parameters to be calibrated is twice the number of generators.
2) Parameter initialization. Assume that the total population size is N , which represents the number of particles that participate in each generation. G is the total generation of iteration. The setting of S depends on the number of oscillation modes of the model to be reduced, which should be analyzed specifically.
3) In the beginning, s, n, and g are set to 1. In each generation, PSO generates N groups of parameters. For each group of parameters, S groups of simulation cases are performed using these parameters in the reduced model. Meanwhile, the same simulations are conducted in the original model. After simulations are finished, equation (10) is used to calculate the similarity between the frequency measurements of the tie lines of the reduced model and that of the original model. Then, equation (11) is used to determine the average fitness in scenario s. Finally, the average fitness of all scenarios is summed up, and F is determined as the fitness function to be optimized according to equation (12). 4) When 3) has been repeated N times using the above N groups of parameters that the PSO has generated, the proposed method will pick the group of parameters that produce the smallest F among them. Then, it will be checked if g = G.
If not, n and s will be initialized to 1 again.
Step 3) will be repeated until the simulation process ends.
5) The PSO will output the parameters obtained, and the reduction of the model will be ended.

IV. CASE STUDY A. REDUCTION OF THE NPCC MODEL
Based on the PSS/E software, the NPCC system is used as the original model, and it has 6 zones and a total number of 140 buses. In this study, Zone-1 is set as the study area that should be kept unchanged, while the structure of the other 5 zones is reduced using the widely-used inertial and slow coherency aggregation algorithm on the DYNRED software [9]. Grid conditions, such as the scale of the system and the range of changes in the grid structure, are given in Tab. 1. The structure of the NPCC system is shown in Fig. 3.
As indicated, the total number of generators in Zone 2-6 is reduced from 39 to 8. Among these zones, the only generator left in Zone-3 is not aggregated. As a result, the total number of generators to be tuned is 7. Each parameter for these generators is given by Tab. 2.  To have a higher accuracy of the reduced model under both the bus-fault event and the generator-trip event, two disturbance scenarios are included to tune the inertias and the damping constants of the above generators. They are a bus-fault event on bus 7 and a generator-trip event on bus 26. Therefore, S is set to 2. The total generation is set to G = 100, and the population size is set to N = 250. T is set to 2 because the frequencies for both of the two tie lines should be taken into consideration. In each generation, the particle with   the best fitness is recorded, and the variation of its fitness is shown in Fig. 4.
As indicated, the best fitness starts at about 0.379, ending at about 0.105 after the 32-nd generation. This demonstrates that the error between the reduced model and the original model, under both the generator-trip event and the bus-fault event, is growing lower with the optimization of these parameters. At last, the parameters are output as below in Tab. 3.

B. SIMULATION AND COMPARISON
To assess the accuracy of the reduced model in presenting the frequency response of the original model under generationtrip contingencies, the cases in Tab. 4 are conducted.
In each contingency, the corresponding generator trips at t=1 s, and the simulation ends at t=5 s. In each case, comparison results are provided among the response of the original model, the response of the reduced model using the initial parameters, and that using the optimized parameters. The results are shown in Fig. 5, Fig. 6, and Fig. 7, respectively.
For the first contingency, a generator connected at bus-23 trips at t=1 s, and the loss of active power is 276.65 MW. As a result, the frequency dips very soon after the occurrence of this event. The frequency variation of the two tie lines is quite similar as they are geographically near. As a result, they show a similar transient characteristic.
By using the parameters obtained from PSO, both the frequencies of the two tie lines match quite well with the curve from the original model. The oscillation mode of the reduced model is also in high accordance with the original model, showing that the reduced model from the PSO has better accuracy than the model obtained from DYNRED software. For instance, by using the initial parameters from DYNRED, the frequency of tie-line 1 matches well with the VOLUME 8, 2020   original one before t=2 s. However, an obvious difference can be seen by observing the second crest at about t= 3.65 s and the third trough at t=4.36 s, respectively. In contrast, by using the parameters outputted from the PSO, there is no obvious frequency mismatch. Also, the nadir frequency error also declines significantly. For contingencies 2 and 3, similar conclusions can be reached.
In comparison with generator-trip contingencies, busfault transients usually include larger magnitudes of system dynamics. To investigate if the proposed method applies to this type of fault scenarios, bus-faults are incepted when t=1 s at bus 22, 23, 36, respectively. The results are shown in Fig. 8, Fig. 9, and Fig. 10, respectively.
As indicated, for faults that occur at bus 22 and 23, the model will have higher accuracy under the bus-fault condition if the parameters of the PSO are used. By this means, the frequencies for both tie line 1 and tie-line 2 matches quite well with that of the original model. For the frequency of tieline 2 in Fig. 10, both the results come from the PSO and that from the DYNRED are not satisfactory. Therefore, further merits are needed to quantify the performance of the proposed method.

C. ERROR ANALYSIS
To quantify the performance of the proposed method, several indicators should be defined. To determine whether the   dynamic response of the reduced model is similar enough to the original one, the main focus should be placed on the frequency match and the magnitude error, which represents the similarity of the oscillation mode and that of the power mismatch. Therefore, the following indicators are defined: Frequency match. Conventionally, the Pearson correlation is a measure of the linear correlation between two variables. It has a value between 1 and −1, where 1 is a total positive linear correlation, 0 is no linear correlation, and −1 is a total negative linear correlation. According to this, equation (15) and (16) is introduced especially to quantify the frequencymatch between the reduced model and the original model. In this equation, the amplitude error has no impact on the calculation result. (15) where ρ st is the frequency-match of the tie line t under event s. N samp is the number of sampling. µ and σ are the mean and the variance of the corresponding measurement.
where ρ s is the total frequency-match of the reduced model under the disturbance scenario of s. 97868 VOLUME 8, 2020 Nadir based frequency response error. Frequency nadirs measure the minimum post contingency frequency in a presented frequency waveform, which is an essential indicator in the power system. This indicator measures the active power change divided by the change in frequency to the nadir point. The results are shown in Fig. 11.
According to the above, the power mismatch is greatly reduced if the reduced model is from the PSO. When a generator-trip event occurs, the power mismatch is no more than 40 MW. When a bus-fault event occurs, the power mismatch is no more than 20 MW. However, if we use the parameters come from DYNRED, the power mismatch can be as high as 58.50-105.3 MW in a generator-trip event, and 10.40-124.80 MW in a bus-fault event, respectively. As a result, the proposed method shows satisfactory performance in reducing the dynamic model of large-scale power systems. Using the same number of generators and the same type of structure, the proposed method reaches higher accuracy in comparison with the existing method.

D. REDUCTION OF THE TEXAS MODEL
In the above section, the parameters of c 1 and c 2 of the PSO are heuristics-based according to the experience of [27]. For this reason, the performance of the PSO in reducing other power systems with different scales should be further studied. In this section, a system with a large size is used to assess the performance of the proposed method. This system is built on the footprint of the electric reliability council of Texas. It has four voltage levels (500/230/161/115 kV) and a total generation capacity of 98 GW, a portion of which is committed to supply a load of 67 GW and 19 GVar. The built of the base case of this model is discussed in [28], and the dynamic information can be found in [29].
In Fig. 12, the zones of the Texas system are shown in different colors, and they are the far west (Zone-1), the north  (Zone-2), the west (Zone-3), the south (Zone-4), the north CE (Zone-5), the south CE (Zone-6), the coast (Zone-7), and the east (Zone-8), respectively. Among them, Zone-7 generates the largest amount of active power and also consumes the largest amount of active power, and is set as the study area. Four disturbance scenarios, including a bus-fault at bus-7045, a generator-trip at bus-7045, a bus-fault at bus-7166, and a generator-trip at bus-7166 are included to tune the inertias and the damping constants of the generators in the external area. In this scenario, S is set to 4. The total generation is set to G = 100, and the population size is set to N = 500. T is set to 3 because the frequency responses of three tielines are taken into account, and the buses of these tie lines in the study area are bus-7018, bus-7274, and bus-7389, respectively. For illustrative purposes, the changes in the parameters for the generators before and after optimization are shown in Fig. 13.
To assess the accuracy of the reduced model under different disturbance scenarios, a bus-fault contingency and a VOLUME 8, 2020  generator-trip contingency on bus-7166 are conducted. The frequency waveform of the three tie lies is recorded, as shown in Fig. 14, and Fig. 15. As seen, by using the parameters obtained from PSO, both the frequencies of the three tie lines match quite well with the curve of the original model. The oscillation mode of the reduced model is also in high  accordance with the original model, showing that the reduced model from the PSO has better accuracy than the model obtained from DYNRED software. Also, the nadir-point mismatch is significantly reduced.
To have a total view of the performance of the proposed method, in Zone-7, a total of 70 bus-fault contingencies and 70 generator-trip contingencies are conducted. Then, the frequency waveforms of the three tie-lines in the reduced model and those of the original model are compared. The results are given in Fig. 16-Fig. 19. As seen, the frequency-match index of each contingency is closer to 1 after the optimization, indicating that the proposed method works well in terms of improving the similarity of the oscillation model of the reduced model. Also, the nadir-point error is significantly reduced in most cases.

V. CONCLUSION
In this paper, the PSO is introduced to optimize the dynamic reduction model of large-scale power systems. Conclusions are as follows.
1) The implementation of the proposed method needs a reduced structure of the external area, usually obtained by the currently-used inertial and slow coherency aggregation algorithm.
2) Using the boundary measurements of the study area, parameters of the external area can be tuned and optimized under certain disturbance scenarios.
3) After the implementation of the PSO, dynamic responses of the reduced model match better with the original one under both generator-trip events and bus-fault events, indicating that the reduced model is a more accurate one in comparison with the model obtained by the inertial and slow coherency aggregation algorithm. 4) Except for the fault scenarios considered in this paper, the performance of the proposed method under asymmetric fault conditions are also needed to be assessed. Such works, including the investigation for the zero-sequence parameters, the modification of the large-scale systems used in this research, the changing of the fitness function, etc. are left for future studies.