Development of a Static Equivalent Model for Korean Power Systems Using Power Transfer Distribution Factor-Based k-Means++ Algorithm

This paper presents a static network equivalent model for Korean power systems. The proposed equivalent model preserves the overall transmission network characteristics focusing on power flows among areas in Korean power systems. For developing the model, a power transfer distribution factor (PTDF)-based k-means++ algorithm was used to cluster the bus groups in which similar PTDF characteristics were identified. For the reduction process, the bus groups were replaced by a single bus with a generator or load, and an equivalent transmission line was determined to maintain power flows in the original system model. Appropriate voltage levels were selected, and compensation for real power line losses was made for the correct representation. A Korean power system with more than 1600 buses was reduced to a 38-bus system with 13 generators, 25 loads, and 74 transmission lines. The effectiveness of the developed equivalent model was evaluated by performing power flow simulations and comparisons of various characteristics of the original and reduced systems. The simulation comparisons show that the developed equivalent model maintains inter-area power flows as close as possible to the original Korean power systems.


Introduction
Confidentiality issues of critical energy infrastructure information limit public access to realistic transmission network models without the dedicated permission requirements. Moreover, a free educational version of a power system simulator provides power system solutions only up to a limited number of buses. However, in Korea, no research efforts have been made to develop a reduced network model, and no simplified equivalent model exists to approximate the overall characteristics of a large practical power system. Therefore, without a network model to study, the potential benefits of power system education and research are impossible.
The power system equivalencing technique represents a large system using a simplified system that approximates the behavior of the original system [1]. Traditionally, power system equivalents were used to speed up power system analysis by reducing the system network size [2]. This network-reduction technique is divided into static and dynamic equivalents [3,4]. The static equivalent methods reduce a system model for power flow studies, whereas the dynamic equivalents reduce the computational demands for transient stability analysis.
Classical Ward and REI equivalents are well-known techniques for realizing static network equivalents [1,5]. These techniques partition the electric power system into internal, boundary, and external system buses. The power system model size is decreased by replacing the external system with small equivalents, whereas the area of interest (the internal system) is unchanged. These techniques focus on maintaining the characteristics of a specific portion of the system; thus, the reduced model might be unsuitable for large-scale network analysis where power flow interactions among areas are required.
Recently, a new network-reduction technique was proposed for system-planning analysis. The method, extensively discussed by Xu and Overbye [6], uses a power transfer distribution factor (PTDF) matrix based on the DC power flow model to identify and aggregate the bus groups. A reduced equivalent system correctly approximates the PTDF characteristics of the original system so that the reduced system maintains the line power flows between areas as close as possible to the original system. Thus, the structural properties of the overall system can be retained in the reduced system [7][8][9].
The correct approximation achieved using the reduced system depends on how well the buses are grouped in the original system. For a large power system, identifying the bus groups is complicated because of the huge dimensions of the PTDF matrix. Shi and Tylavsky [10] evaluated various bus-clustering methods based on data-processing algorithms, with k-means++ methods showing the best performance. The k-means++ method overcomes the limitations of the k-means method by advanced initial bus selection and the Euclidean distance of the PTDF matrix.
For the reduction in large power systems, all the characteristics of the original system cannot be preserved; thus, the reduction processes generally focus on specific elements of the power system. In [11], the authors noted that the south/southeastern part of Brazil's power system can be reduced to a seven-bus system with five generators. The equivalent system correctly represents the electromechanical oscillations and has been used in several small-signal stability studies. Reduction in southern and eastern Australia's power system to a 14-generator system was achieved by Gibbard [12] by focusing on modal control problems. The reduced model was used to validate the control effects of a power system stabilizer in multimachine power systems. In [13], the simplified two-area interconnected Hokkaido-Honshu power system of eastern Japan was used in order to investigate the effectiveness of fuzzy-PID controller for improved frequency stability. In this manner, with clear objectives for use, the reduced equivalent models can work effectively.
In this study, a static equivalent model preserving the overall system topology of Korean power systems was developed using the PTDF-based k-means++ bus-clustering method. The original system of more than 1600 buses was reduced to a 38-bus system with 13 generators, 25 loads, and 74 transmission lines. The overall system topology was maintained by developing the PTDF matrix and using the k-means++ algorithm for bus-clustering. The bus groups were replaced by a single bus with a generator or load depending on the net power injection in the group. The equivalent impedance of the transmission lines was determined to preserve the power flows among the identified areas as close as possible to the original system. For an advanced representation of the original system, an appropriate voltage magnitude of a bus was selected, and real power losses, ignored in DC power flow models, were compensated.
The next three sections of this paper are organized as follows. Section 2 describes the method for static model reduction using the PTDF matrix and k-means++ algorithm and introduces the detailed procedures for determining transmission line impedances between the areas and line loss compensation. Section 3 presents the developed equivalent model for Korean power systems and performs simulation comparisons with original systems to evaluate the performance of the reduced system. Section 4 presents the conclusions.

PTDF Matrix
PTDF is a linear approximation of the active power flow sensitivity in a specific transmission line regarding the shift of real power [14]. In other words, when 1 MW of real power is injected into a specific bus, a change in the power flow in each line can be approximated using the PTDF calculation. A PTDF matrix represents the difference in the real power flow over any line because of the power injection at all system buses.
The PTDF matrix is fully based on the DC power model. The corresponding power flows and power injections are linearly related to a susceptance matrix and bus voltage angles, respectively [15,16]. The PTDF matrix, shown in Figure 1, can be expressed as follows: where B branch and B bus are the branch and bus susceptance matrices with dimensions of L × N and N × N, respectively. L is the number of transmission lines, N is the number of buses, except for a slack bus, x is the transmission line reactance vector, and C is the incidence matrix with dimensions of L × N.
Energies 2020, 13, x FOR PEER REVIEW 3 of 12 calculation. A PTDF matrix represents the difference in the real power flow over any line because of the power injection at all system buses. The PTDF matrix is fully based on the DC power model. The corresponding power flows and power injections are linearly related to a susceptance matrix and bus voltage angles, respectively [15,16]. The PTDF matrix, shown in Figure 1, can be expressed as follows: where ℎ and are the branch and bus susceptance matrices with dimensions of L × N and N × N, respectively. L is the number of transmission lines, N is the number of buses, except for a slack bus, is the transmission line reactance vector, and C is the incidence matrix with dimensions of L × N. Each column ℎ ( = 1,2, … , ) in this matrix has a corresponding PTDF value, and each column corresponds to a bus and each row corresponds to a line in the system.

Bus-Clustering Using the Power Transfer Distribution Factor-Based k-means++ Algorithm
Using the formulated PTDF matrix, the bus groups showing similar PTDF characteristics can be identified by comparing the elements in the column of the matrix. The Euclidean distance measures the similarity, and the entire system can be partitioned into different areas or clusters, as given in Equation (2): where ℎ is the ℎ column of the PTDF matrix, ϕ is the component in the ℎ row and ℎ column of the PTDF matrix, and L is the number of transmission lines. The k-means++ algorithm is an improved k-means algorithm used to form the bus clusters according to similarity measures from the Euclidean distances of each ℎ in the PTDF matrix. In the k-means algorithm, the initial centroids in the clusters were randomly selected, and the Euclidean distance between ℎ and the cluster's centroid was calculated. The cluster members were identified using the closest Euclidean distance, and the centroid of each cluster was updated by averaging all vectors within the cluster. These steps were repeated until the desired number of clusters was determined. However, for a large power system network, the k-means approach has limitations of nonconvergence and the appearance of many empty clusters because the k-means algorithm randomly selects the initial centroids without knowing the appropriate number of clusters [17]. Problems associated with the k-means algorithm can be overcome using a k-means++ algorithm in the iteration process, improving the convergence speed and resulting in better cluster formations [17,18].
For the bus-clustering with the k-means++ algorithm, the first centroid was randomly selected Each column h i (i = 1, 2, . . . , N) in this matrix has a corresponding PTDF value, and each column corresponds to a bus and each row corresponds to a line in the system.

Bus-Clustering Using the Power Transfer Distribution Factor-Based k-Means++ Algorithm
Using the formulated PTDF matrix, the bus groups showing similar PTDF characteristics can be identified by comparing the elements in the column of the matrix. The Euclidean distance measures the similarity, and the entire system can be partitioned into different areas or clusters, as given in Equation (2): where h i is the ith column of the PTDF matrix, φ ij is the component in the ith row and jth column of the PTDF matrix, and L is the number of transmission lines. The k-means++ algorithm is an improved k-means algorithm used to form the bus clusters according to similarity measures from the Euclidean distances of each h i in the PTDF matrix. In the k-means algorithm, the initial centroids in the clusters were randomly selected, and the Euclidean distance between h i and the cluster's centroid was calculated. The cluster members were identified using the closest Euclidean distance, and the centroid of each cluster was updated by averaging all vectors within the cluster. These steps were repeated until the desired number of clusters was determined. However, for a large power system network, the k-means approach has limitations of nonconvergence and the appearance of many empty clusters because the k-means algorithm randomly selects the initial centroids without knowing the appropriate number of clusters [17]. Problems associated with the k-means algorithm can be overcome using a k-means++ algorithm in the iteration process, improving the convergence speed and resulting in better cluster formations [17,18].
For the bus-clustering with the k-means++ algorithm, the first centroid was randomly selected from all PTDF vectors, followed by calculation of its Euclidean distance (d i1 , d i2 , . . . , d ik , . . .) from the centroid to all existing noncentroids for each vector h i . For the remaining noncentroid vector, the one with the largest D(h i ) = max{d i1 , d i2 , . . . , d ik , . . .} was selected and designated as the new centroid. This process was continued until the desired number of centroids was obtained. Once the iteration was completed, the same method as in the k-means algorithm was used to identify clustering information. Figure 2 illustrates the clustering process using the k-means++ algorithm for a set of data in a two-dimensional Euclidean space. By using the k-means++ algorithm, a converged solution was obtained, and the number of clusters could be controlled [17].
Energies 2020, 13, x FOR PEER REVIEW 4 of 12 process was continued until the desired number of centroids was obtained. Once the iteration was completed, the same method as in the k-means algorithm was used to identify clustering information. Figure 2 illustrates the clustering process using the k-means++ algorithm for a set of data in a twodimensional Euclidean space. By using the k-means++ algorithm, a converged solution was obtained, and the number of clusters could be controlled [17].

Determining the Equivalent Transmission Line Impedance between Clusters
During the reduction process, a single bus replaced each cluster and the aggregated generator and load was attached to the bus. The equivalent impedance of the line connecting the identified clusters was determined. The reduced system preserves the inter-area power flows. For the original system, the inter-area power flows could be calculated by summing all line power flows over the transmission lines connecting the corresponding areas as follows: where − is an × 1 vector containing the inter-area power flows, is the number of transmission lines in the reduced system, Π is an × matrix summing the line flows, Φ is the PTDF matrix, and is the bus-injection vector. After aggregating the bus groups, the reduced network's power flow is different from that of the original network. The formula for the reduced model is given by Equation (4), which is obtained from the definition of PTDF with a reduced transmission line reactance vector and reduced busbranch incidence matrix : where subscript R denotes the reduced system. The principal constraint associated with network reduction is to find a set of transmission line

Determining the Equivalent Transmission Line Impedance between Clusters
During the reduction process, a single bus replaced each cluster and the aggregated generator and load was attached to the bus. The equivalent impedance of the line connecting the identified clusters was determined. The reduced system preserves the inter-area power flows. For the original system, the inter-area power flows could be calculated by summing all line power flows over the transmission lines connecting the corresponding areas as follows: where P inter−area f low is an L R × 1 vector containing the inter-area power flows, L R is the number of transmission lines in the reduced system, Π f low is an L R × L matrix summing the line flows, Φ is the PTDF matrix, and P inj is the bus-injection vector.
After aggregating the bus groups, the reduced network's power flow is different from that of the original network. The formula for the reduced model is given by Equation (4), which is obtained from the definition of PTDF with a reduced transmission line reactance vector x R and reduced bus-branch incidence matrix C R : where subscript R denotes the reduced system. The principal constraint associated with network reduction is to find a set of transmission line reactance minimizing the power flow difference between the original and reduced systems. The reactance was identified using the optimization method proposed by Shi and Tyvalsky [19] and can be expressed as follows: where x i R is the ith component of the transmission reactance vector and · represents the Euclidean norm. In the reduced model, the injections at each bus are the sum of the injections at all buses within each cluster in the original system. For example, if the amount of real power generation in the cluster is greater than the real power load, the extra real power injection is represented by a generator. For the opposite case, a load represents the extra real power consumption in the cluster. Equation (6) expresses the reduced system power injection.
where Π g is the sum of the bus injections with dimensions of an N R × N matrix, and each inter-area power flow is the summation of all contributions from different bus injections.

Selecting an Appropriate Bus Voltage Magnitude
During the reduction process, when a single bus replaced the bus groups, an appropriate bus voltage magnitude was determined, maintaining the overall voltage characteristics of the original system. For developing the equivalent method, the average value of voltage magnitudes within a cluster was selected, as given by Equation (7). For a generator bus, the set voltage was adjusted to match the average value, whereas for a load bus, a shunt capacitor was provided, supporting the required voltage level.
where V i average is the average voltage of the ith area, V k is the bus voltage magnitude at bus k, and n is the number of buses in the ith area.

Real Power Loss Compensation
The DC power flow model used for developing the PTDF matrix neglects real power losses. Although real power losses on a transmission line are relatively small, their cumulative sum could cause a significant loss of simulation accuracy. In the equivalent model, real power losses in the original system were compensated by adding resistance to the transmission line connecting two areas. Resistances were added to minimize changes in the transmission real power flows of the equivalent model. The appropriate resistance was determined using Equation (8): where R i,j compensated is the resistance of the transmission line between areas i and j, I i,j R is the magnitude of the current over the transmission line between areas i and j in the equivalent model, P k loss denotes the transmission line real power losses in the original system, and k is a transmission line connecting areas i and j.

Summary of the Procedures
Developing an efficient network-reduction model involves five steps, namely: (1) bus-clustering using a PTDF-based k-means++ algorithm, (2) aggregating the bus groups and converting them into a single equivalent bus, (3) determining the equivalent line impedance, (4) adjusting the bus voltage magnitude, and (5) compensating for real power losses on the transmission line. This study proposes the last two steps for better approximation using the developed equivalent model. Figure 3 illustrates these steps.

Developing an Equivalent Model for Korean Power Systems
The reduction approach was implemented using MATLAB and applied to a practical Korean power system. Table 1 presents detailed information about the original system used for the development.

Developing an Equivalent Model for Korean Power Systems
The reduction approach was implemented using MATLAB and applied to a practical Korean power system. Table 1 presents detailed information about the original system used for the development. The number of buses for the reduction process can be selected by considering the appropriate representation of the overall Korean power system topology. The limitation of the number of buses was also considered for the free use of a commercial power system simulator. Easy access to the equivalent model was expected for undergraduate and graduate students. Thirty-eight clusters were selected using the Euclidean distance of the PTDF matrix given in Equation (2) and the k-means++ algorithm. The largest cluster includes 280 buses, and the smallest one has only 11 buses. Most clusters have between 20 and 50 buses. These bus groups in the original network were replaced by single buses in the reduced system. Figure 4 presents the corresponding buses and their geographic locations.  The power injections from all buses within an identified cluster were summed up, and a generator or load, depending on positive or negative injections, was connected to the corresponding bus. The generation and load were scaled down to one-tenth for simplicity. Figure 4 shows that the generator buses are marked as G, whereas the unmarked buses are the corresponding load buses. The buses were grouped around cities with large loads and near power generation complexes.
When there was at least one transmission line between clusters, a corresponding line was created in the reduced model. In total, 74 transmission lines were identified. The equivalent impedance of a line can be calculated by solving the optimization problem in Equation (5). The PTDF matrix for the reduced system has dimensions of 74 × 37, except for a slack bus. The reduced model comprises 13 generator buses, 25 load buses, and 74 transmission lines. Figure 5 shows a one-line diagram of the reduced Korean power system, and Table 2 presents the information for the reduced system. Table 2. Reduced Korean power system information. The power injections from all buses within an identified cluster were summed up, and a generator or load, depending on positive or negative injections, was connected to the corresponding bus. The generation and load were scaled down to one-tenth for simplicity. Figure 4 shows that the generator buses are marked as G, whereas the unmarked buses are the corresponding load buses. The buses were grouped around cities with large loads and near power generation complexes.
When there was at least one transmission line between clusters, a corresponding line was created in the reduced model. In total, 74 transmission lines were identified. The equivalent impedance of a line can be calculated by solving the optimization problem in Equation (5). The PTDF matrix for the reduced system has dimensions of 74 × 37, except for a slack bus. The reduced model comprises 13 generator buses, 25 load buses, and 74 transmission lines. Figure 5 shows a one-line diagram of the reduced Korean power system, and Table 2 presents the information for the reduced system. The performance of the equivalent model was evaluated by comparing the corresponding power flows through the original and reduced models. Figure 6 presents a simulation comparison. The power flows in the original model are a summation of power flows between the corresponding areas, with an average difference of 0.4 MW, and a maximum value of 22.5 MW. Figure 6 shows that the equivalent model correlates well with the original model. Therefore, the characteristics of the overall power flow among these areas can be maintained using the developed equivalent model. The performance of the equivalent model with a changed operating point was validated by injecting an additional 100 MW of real power into bus 3 in the reduced system. In Figure 7   The performance of the equivalent model was evaluated by comparing the corresponding power flows through the original and reduced models. Figure 6 presents a simulation comparison. The power flows in the original model are a summation of power flows between the corresponding areas, with an average difference of 0.4 MW, and a maximum value of 22.5 MW. Figure 6 shows that the equivalent model correlates well with the original model. Therefore, the characteristics of the overall power flow among these areas can be maintained using the developed equivalent model. The performance of the equivalent model was evaluated by comparing the corresponding power flows through the original and reduced models. Figure 6 presents a simulation comparison. The power flows in the original model are a summation of power flows between the corresponding areas, with an average difference of 0.4 MW, and a maximum value of 22.5 MW. Figure 6 shows that the equivalent model correlates well with the original model. Therefore, the characteristics of the overall power flow among these areas can be maintained using the developed equivalent model. The previous evaluation for considering the changed operating point was extended to all the system buses. With the additional 100 MW injection or withdrawal at a bus, the differences of the real power flow variations between the original and the reduced systems were compared. The average errors with respect to the 100 MW change at the corresponding bus are presented in Figure 8. As shown in the Figure 8, the maximum value of the average error with the injection is 2.7 MW at bus 32 and the minimum value is 0.11 MW at bus 6. With the withdrawal, the maximum error is 2.7 MW at bus 32 and the minimum one is 0.02 MW at bus 33. The average errors are not large enough to produce a significant difference between two systems. Thus, even with the changed operating points from both injection and withdrawal, the power flow characteristics of the original system can be preserved with the reduced system. The average voltage of each group in the original system was compared with the corresponding bus voltage of the reduced system. For a generator bus in the reduced system, the voltage magnitude was set to the average cluster voltage magnitude in the original system. For a load bus, an additional shunt capacitor was provided, matching the average voltage level. Ten shunt capacitors were added The previous evaluation for considering the changed operating point was extended to all the system buses. With the additional 100 MW injection or withdrawal at a bus, the differences of the real power flow variations between the original and the reduced systems were compared. The average errors with respect to the 100 MW change at the corresponding bus are presented in Figure 8. As shown in the Figure 8, the maximum value of the average error with the injection is 2.7 MW at bus 32 and the minimum value is 0.11 MW at bus 6. With the withdrawal, the maximum error is 2.7 MW at bus 32 and the minimum one is 0.02 MW at bus 33. The average errors are not large enough to produce a significant difference between two systems. Thus, even with the changed operating points from both injection and withdrawal, the power flow characteristics of the original system can be preserved with the reduced system. The previous evaluation for considering the changed operating point was extended to all the system buses. With the additional 100 MW injection or withdrawal at a bus, the differences of the real power flow variations between the original and the reduced systems were compared. The average errors with respect to the 100 MW change at the corresponding bus are presented in Figure 8. As shown in the Figure 8, the maximum value of the average error with the injection is 2.7 MW at bus 32 and the minimum value is 0.11 MW at bus 6. With the withdrawal, the maximum error is 2.7 MW at bus 32 and the minimum one is 0.02 MW at bus 33. The average errors are not large enough to produce a significant difference between two systems. Thus, even with the changed operating points from both injection and withdrawal, the power flow characteristics of the original system can be preserved with the reduced system. The average voltage of each group in the original system was compared with the corresponding bus voltage of the reduced system. For a generator bus in the reduced system, the voltage magnitude was set to the average cluster voltage magnitude in the original system. For a load bus, an additional shunt capacitor was provided, matching the average voltage level. Ten shunt capacitors were added to the reduced system. Figure 9 presents a comparison of the bus voltage magnitudes in the original The average voltage of each group in the original system was compared with the corresponding bus voltage of the reduced system. For a generator bus in the reduced system, the voltage magnitude was set to the average cluster voltage magnitude in the original system. For a load bus, an additional shunt capacitor was provided, matching the average voltage level. Ten shunt capacitors were added to the reduced system. Figure 9 presents a comparison of the bus voltage magnitudes in the original and reduced systems. The maximum difference in the voltage magnitude was less than 0.005 p.u. For comparison, the voltage magnitudes before and after adding the shunts were analyzed. The results show that the reduced system without the shunt compensation shows a greater deviation from the original system voltage characteristics, whereas the reduced system with the compensation exhibits a negligible difference.
Energies 2020, 13, x FOR PEER REVIEW 10 of 12 show that the reduced system without the shunt compensation shows a greater deviation from the original system voltage characteristics, whereas the reduced system with the compensation exhibits a negligible difference.    Table 3 shows the comparison of real power errors between the original and reduced models. In Table 3, all real power flows and losses over transmission lines for the entire system are summed. The average error in inter-area power flows is less than 1%. The total real power losses have an error of 2.6 MW, which is a small amount compared to the total line power flows.  show that the reduced system without the shunt compensation shows a greater deviation from the original system voltage characteristics, whereas the reduced system with the compensation exhibits a negligible difference.    Table 3 shows the comparison of real power errors between the original and reduced models. In Table 3, all real power flows and losses over transmission lines for the entire system are summed. The average error in inter-area power flows is less than 1%. The total real power losses have an error of 2.6 MW, which is a small amount compared to the total line power flows.  Table 3 shows the comparison of real power errors between the original and reduced models. In Table 3, all real power flows and losses over transmission lines for the entire system are summed.
The average error in inter-area power flows is less than 1%. The total real power losses have an error of 2.6 MW, which is a small amount compared to the total line power flows.

Conclusions
This study developed a static equivalent model for Korean power systems using a PTDF-based k-means++ algorithm. The reduced system represents the power flows between areas as close to the original system as possible. The equivalent model was created following five steps, namely: (1) identification of the clusters using Euclidean distance with PTDF matrix and k-means++ algorithm, (2) aggregating the bus groups and converting them into a single generator or load bus, (3) determining the equivalent line impedance, (4) setting the voltage magnitude to the average value of the cluster, and (5) compensating for real power line losses. A large practical Korean power system with more than 1600 buses and over 85 GW of generation capacity was reduced to a 38-bus system with 13 generators, 25 loads, and 74 transmission lines. The performance of the reduced model was evaluated regarding its ability to preserve the characteristics of the original system by comparing the power flows, voltage magnitudes, and real power losses of the reduced and original models. The average power flow error was 0.4 MW, and even with the changed operating point from additional 100 MW injection or withdrawal at all buses, the average error was not large enough. Therefore, the developed equivalent model maintains the overall power flow characteristics of the original Korean power system. The developed equivalent Korean power system can be used for research and education purposes without confidentiality issues of critical energy infrastructure information. Undergraduate and graduate students will have more learning opportunities and experience of the specific characteristics of Korean power systems. As a continuation of this work, a dynamic reduced model focusing on preserving inter-area oscillations will be developed based on the static equivalent model.