Fast approach for transient stability constrained optimal power flow based on dynamic reduction method

: The main challenge to solve the transient stability constrained optimal power ﬂ ow (TSC-OPF) problem is its huge dimension due to numerous discretised transient variables and constraints. This problem becomes more serious when large power systems are considered. This study presents a fast approach to realise a global TSC-OPF based on dynamic reduction method, which decomposes the power system into several coherent areas and represents the original system by a reduced equivalent system. In this approach, the single transient stability constraint is obtained by simulating the reduced system instead of the full system. The new approach reduces the simulation execution time and thus increases the ef ﬁ ciency of the TSC-OPF. Two case studies indicate that the proposed approach can remarkably reduce the CPU time of the TSC-OPF procedure, compared with the TSC-OPF based on the full-system simulation. The new approach is very practical in solving the TSC-OPF problem in large power systems where numerous machines are coherent.


Introduction
Optimal power flow (OPF) is a well-established tool in power system planning and operation. But it is not suitable for modern power systems operating under market rules, where system operators tend to operate electric power systems close to their security limits due to economic constraints and intensified transactions [1]. Under such operation situations, security-constrained OPFs are more appreciated. One of them is the transient stability constrained optimal power flow (TSC-OPF), which assures the power system operating stably under contingency conditions. Mathematically, the TSC-OPF is a non-linear programming problem with variables and constraints in time domain. To solve TSC-OPF problems, the time-varying variables and constraints are discretised and are included into the OPF formulation. By solving this optimisation problem, an optimal generation dispatch can be determined [2][3][4][5][6][7][8][9][10][11]. However, the corresponding programming model has huge dimension and requires heavy computation. When large power systems are considered, a substantial increase of the programming dimension will make the computation time and the memory consumption unacceptable. Some authors propose to reduce the system's discretised dimension by special discretisation methods [2,3]. Others eliminate unnecessary constraints and variables by using reduced admittance matrix of power systems [4,5]. Owing to their efficiency and robustness, modern optimisation techniques are also proposed to solve the TSC-OPF problems [6][7][8][9][10][11]. Recently, some inspiring approaches based on single-machine equivalent (SIME) model have also been proposed [12][13][14][15][16]. The SIME model aggregates the multi-machine system angle trajectories to a one-machine infinite bus (OMIB) equivalent angle trajectory [17]. Thus, the stability constraints in the TSC-OPFs are reduced to a single trajectory of the OMIB. The required simulation time is also reduced by setting the simulation time at the time to instability.
Based on the SIME method, a new approach is proposed to realise the TSC-OPF using an independent dynamics simulation [18,19]. In this approach, the dynamic time-varying constraints are taken out of the optimisation formulation and are implemented by an independent simulation algorithm. Based on the simulation results, a single transient stability constraint is derived and is then included into the conventional OPF optimisation model. Using this approach, the dimension of optimisation model of the TSC-OPF is largely reduced, while the non-linear characteristic of power system is taken into account. Although this approach succeeds in reducing the dimension of optimisation model, the total computation time of the TSC-OPF is still very long for large power systems due to the time-consuming dynamics simulation.
It is well known that, in large power systems, some generators have similar swing patterns when the system is suffering a disturbance. In order to reduce computation memory and simulation time, the concept of coherence of generators is established and the dynamic reduction method is proposed [20]. Aggregated or reduced system is commonly used by operators and engineers for power system transient stability studies. Successful implementation of the dynamic reduction method can contribute towards significant computation efficiency and even real-time security assessment. A review of dynamic reduction method is summarised in [21].
The principle of the SIME relies on the observation that the loss of synchronism of a multi-machine system originates from the machine separation into two groups: (i) the group of machines that loose synchronism (critical machines, CM) and (ii) all other machines (non-critical machines, NM). The SIME uses the physical parameters and the time-varying data of the CM and the NM groups to generate an equivalent machine for each group. The two groups of machines are further reduced to an equivalent OMIB model, from which the equal critical area can infer the transient stability bound of the original multi-machine system. Broadly speaking, the SIME method is also an aggregation method, which aggregates the machines of the CM and NM groups separately into two equivalent machines. But this aggregation procedure follows a dynamics simulation procedure and therefore does not have any influence on the simulation efficiency.
In the CM or NM groups, it is possible that some machines have similar swing patterns and can be aggregated into an equivalent machine. When the dynamic reduction method is used with the SIME method, a fast transient stability assessment will be obtained. Consequently, the efficiency of the TSC-OPF will be improved.
In this paper, a fast approach is proposed to realise the TSC-OPF by a dynamics simulation based on the dynamic reduction method. In this approach, the power system is decomposed into several coherent areas and is represented by a smaller equivalent system. Using the reduced equivalent system, a transient stability constraint is derived quickly and is included into the optimisation algorithm. Comparing with the TSC-OPF procedure based on full-system simulation, the new approach will reduce more the computation time of the TSC-OPF and is very practical for large power systems, which generally have wide coherent areas.
The rest of the paper is structured as follows: Section 2 introduces a TSC-OPF approach using an independent dynamics simulation. Section 3 describes the proposed approach. The description of the dynamic reduction method to assess the transient stability is presented. Section 4 discusses the results obtained from the application of the proposed approach to the IEEE 10-machine and IEEE 50-machine benchmark systems. The comparisons with the TSC-OPF based on the full-system simulation are given, demonstrating the benefits of the proposed approach. Finally, the contributions of this paper are summarised and highlighted in Section 5.
2 TSC-OPF with a single stability constraint

Formulation of TSC-OPF
Assuming that there is only one contingency, the general formulation of the TSC-OPF, is expressed as follows 3. Transient stability constraints where f is the cost function of the generator active powers P g .
x(t) are the state variables of power system such as machine angles, velocities, and so on. y(t) are the algebraic variables such as bus voltage magnitudes and angles. G s (x 0 , y 0 , P g ) and H s (x 0 , y 0 , P g ) are the steady-state equality and inequality constraints, such as power balance equations and system physical and operation limits. x 0 and y 0 are the steady-state values (or initial conditions) of x(t), y(t). Note that (1) and (2)

TSC-OPF formulation with an independent dynamics simulation
The TSC-OPF formulated by (1)-(4) is a non-linear optimisation problem with variables and constraints in time domain. If a totally explicit global approach is used, the above TSC-OPF problem is very difficult to solve due to its huge dimension, which comes mainly from the discretisation of constraints (3) and (4). Based on the SIME model, in [13], it is proposed to reduce the multi-machine stability constraints (4) to one single constraint at time step t u as where δ UT (t u ), δ CT (t u ) are the OMIB angular deviation at t u on the unstable and the critical trajectory, considering the original fault clearing time t f and the computed CCT, respectively. Δδ tu is a desired threshold. In (5), δ CT (t u )i sa fixed scalar value. δ UT (t u ) is a function of the CM and NM angles at t u and is expressed as where M i and δ i are the inertia and the angle of machine i, respectively. The inertia coefficients of the CM and NM groups are An important dimension reduction in the TSC-OPF optimisation model is achieved when (4) is replaced by (5), since now it is only required to stabilise the single rotor trajectory δ UT (t u ), instead of stabilising each machine rotor trajectory δ i (t u ) in the original TSC-OPF formulation.
Moreover, the sensitivity analysis shows that the simulation time (or control time) t ctrl can be applied at any selected point during period [t 0 , t u ] [ 15,16]. A practical approach is achieved in [16] by setting t ctrl to t 0 . The dynamic constraints (3) are completely eliminated from the optimisation model. The size of the resulting optimisation problem is reduced to one very similar to that of the conventional OPF.
To include the non-linearity of power system into the transient stability assessment, in [19], the authors propose to set t ctrl to t u and use an independent dynamics simulation to realise the dynamic constraints (3). In this way, the dynamic constraints are also eliminated from the optimisation model. This approach has the same dimension of optimisation model as in [16] and is described in brief as follows.
When the single transient stability (5) is used, the TSC-OPF formulation consists of (1)-(3) and (5), where the constraints (2), (3) and (5) are coupled. The steady-state constraints (2) establish the relationship between the steady-state variables x 0 , y 0 and P g . The dynamic constraints (3) link the initial values x 0 , y 0 and the final values x(t u ), y(t u ), whereas the variable δ UT (t u ) in (5) is deduced from x(t u ), y(t u ) by using the OMIB model (6). Their relationships are shown in Fig. 1.
From Fig. 1, it can be seen that the variable δ UT (t u ) explicitly depends on the control variable P g . During optimisation procedure, once the control variable P g changes, the stability variable δ UT (t u ) changes too, which can be realised by simulating the dynamic constraint (3) from [t 0 , t u ] in an independent way. On the other hand, since the procedure of calculating δ UT (t u ) includes automatic implementation of the dynamic constraints (3), the dynamic constraints (3) can be eliminated from the TSC-OPF formulation as follows Min fP g s.t.
G s x 0 , y 0 , P g = 0 In the above formulation, δ UT (t u ) is the result of the independent dynamics simulation procedure shown in Fig. 1. δ CT (t u )i safixed scalar value. Comparing with a conventional OPF formulation, the above formulation adds only one single constraint. The formulation (7) can be solved by a conventional OPF algorithm with small modification. The advantages and implementation of the formulation (7) are given in [19].
3 Transient stability constraint with dynamic reduction method

Principle of dynamic reduction method
The procedure to establish the relationship between x(t u ), y(t u ) and x 0 , y 0 is in fact the dynamics simulation of the power system subjected to a disturbance. This simulation requires large amount of computation efforts in terms of time and memory when a large power system is considered. Using the dynamic reduction method, a large power system can be represented by a smaller equivalent system while maintaining the essential dynamic behaviour of the original system. The dynamic reduction process can be divided into two main steps [20]: 1. Determination of coherent generator groups 2. Aggregation of generators and reduction of the network The most important part of the dynamic reduction is to identify the groups of coherent generators. Generators are said to be 'coherent' when they have similar swing patterns following a disturbance. The coherency condition generally simplifies to Generators i and j that satisfy (8) are coherent and can be aggregated into a single equivalent generator. There are several methods to identify the group of coherent generators. They include electrical distance method, time-domain approach and frequency-domain approach, which are described in [21]. In our paper, the slow coherency method [22] is chosen, due to its advantages such as almost independence of disturbance. This characteristic is very convenient for the TSC-OPF, which requires the simplest equivalent system under different fault conditions to improve the efficiency.
Slow coherency is an application of the singular perturbation or two-time-scale method in power systems [23]. The method assumes that the state variables of an nth-order system are divided into r slow states, and (n − r) fast states, in which the r slowest states represent r groups with the slow coherency. To find the slow coherent areas, the method requires the calculation of the slow eigenvector matrix of the electromechanical model of the power system. The r (number of areas) most linearly independent rows of the eigenvector matrix will become the reference generators. A grouping algorithm is then applied to group non-reference generators (n − r) to the reference generators (r). The grouping algorithm is summarised as follows [23]: 1. Choose the number of groups r based on the slowest modes 2. Compute the eigenvector matrix U for the r smallest eigenvalues 3. Apply Gaussian elimination with complete pivoting on U to obtain r reference machines. Each group will have only one reference machine. 4. Order the first r rows of U (called U 1 ) according to the order found in step 3 and compute L as 5. Use the L matrix to assign other machines to the coherent areas according to the largest entry in each row of L.
After the groups of coherent generators are identified, the power network can be reduced by aggregating the coherent machines into a single machine. There are many aggregation techniques, classified by the machine model. For the classical model, there are principally two methods: inertial aggregation and slow coherency aggregation. For a detailed model, the synchronic model equivalence and the structure preservation techniques can be applied [24]. In this paper, the inertial aggregation technique is used since it requires less computation and represents very well the dynamics of the coherent generators in short time term. The aggregation procedure includes two steps: 1. Nodal aggregation: An equivalent bus replaces the terminal buses of the coherent generators such as (10) where V eq is the voltage of the equivalent bus. V i , S i are the voltage and MVA of bus i connected with coherent generator i. α i is the complex turn's ratio between the original voltage and the equivalent voltage. 2. Machine aggregation: Using the inertia techniques, the parameters of the equivalent machine are calculated by where S eq , H eq and X′ deq are the MVA, inertia and transient reactance of the equivalent generator. S i , H i and X′ di are the MVA, inertia and transient reactance of coherent generator i.

SIME model with dynamic reduction method
As shown in (6), the angle δ UT (t u ) of the OMIB is obtained from the angles of two defined groups: CM and NM, which are calculated from the corresponding individual machine's angle δ i (t u ) if the full-system simulation is used. When the dynamic reduction method is introduced, the angle δ UT (t u ) is calculated from the machine angles of the reduced system. The procedure to integrate the dynamic reduction method in the TSC-OPF is summarised as follows: 1. The proposed TSC-OPF is initialised using a conventional OPF algorithm consisting of minimising (1) subjected to the steady-state constraints (2). 2. For a given unstable contingency, run time-domain simulation to identify the CM and NM groups in the system. Construct the OMIB model to find t u , CCT and δ CT (t u ). t u is determined by the criteria described in [17].
3. Identify the groups of coherent generators in the system by using the slow coherency method. The identification of coherent groups of machines is based on the coherency map matrix [25] where v i is the ith row of U whose columns are the eigenvectors corresponding to the lowest frequency mode. A tolerance C tol is selected, so that machines i and j are said to be coherent if 4. Aggregate the coherent machines in the NM or CM groups. Each coherent group is aggregated into one equivalent machine by using (10) and (11). 5. Construct the dynamic equations of the reduced systeṁ where x red (t), y red (t) are the state and algebraic variables of the reduced systems. P g_red is the active powers of the machines in the reduced system. 6. Construct the new OMIB model using the reduced equivalent system.
where  [26]) of the OMIB angle trajectories of the reduced system as where t u is the time to instability of the SIME model of the full system. δ UT_red (t) and δ UT_f (t) are the OMIB angle trajectories using the reduced and the full system, respectively. The maximum AAE are defined as that produced by the reduced system which has the same t u and CCT as the full system in tolerance error (one time step).
8. Construct the stability constraint based on the reduced system simulation as δ UT_red (t u ) is calculated as in the same procedure as shown in Fig. 1, except that the variables P g red , x 0 red , y 0 red of the reduced system are used. 9. Form the TSC-OPF formulation based on the reduced-system simulation as 10. Solve the TSC-OPF formulation (18) using a conventional OPF algorithm. Note that, in (18), the optimisation variables P g , x 0 and y 0 are the variables of the full system. The outputs of the optimisation algorithm are the same as the TSC-OPF based on the full-system simulation.   Fig. 2 depicts the total procedure to integrate the dynamic reduction method into the TSC-OPF. This procedure can be divided into two main steps: the 'off-line' construction of the reduced system and the 'on-line' simulation of the reduced system. The 'off-line' steps are the above steps before the optimisation process starts in Fig. 2 and are used to construct a proper reduced equivalent system for the full system. The 'on-line' steps are the simulation of the reduced system during the optimisation process and are used to offer the stability constraint (17) for the optimisation algorithm. During the optimisation process, the reduced system is simulated iteratively to estimate the derivate ∇δ tu_red (P g ) in the perturbed Karush-Kuhn-Tucker (KKT) optimality conditions of (18) [19].
The efficiency of this method depends on the ratio between the number of machines in the full system and that in the reduced system. The more the machines swing coherently in the original system, the less there is machines in the reduced system, and the more efficient this method will be. Obviously, the dynamic reduction method is not worthwhile for small systems with just a few coherent machines. However, for large power systems with wide coherent areas, this method may become very attractive.

Case studies
In this section, two well-known power systems, the IEEE 39-bus 10-machine system and the IEEE 145-bus 50-machine system, are used to validate the proposed method. All the simulations are carried out using the MATLAB 7.14.0 platform in a PC with Intel Core™ i7 3930 K processor, 3.2 GHz, 8 GB RAM. The steady-state OPF for the two systems are obtained using MATPOWER [27]. The loads of the power systems are considered as constant impedances. For simplification, the machine classical model is used. The power system toolbox (PST) [28] is employed to find the groups of coherent machines.

IEEE 39-bus 10-machine system
The diagram of the test system is shown in Fig. 2. The data of this system and the machine parameters come from [29]. The generator ratings and quadratic generation cost function are taken from [30]. The steady-state voltage magnitude limits in p.u. for generators and load buses are 0.95 < V g < 1.09 and 0.95 < V l < 1.07, respectively. The base case operating point is presented in Table 1.
By choosing C tol = 0.95, the 10-machine system is partitioned into six coherent areas and is shown by the dotted lines in Fig. 3. The detailed grouping information is shown in Table 2. This partitioning is the same as in [31].
Case 1: The first contingency scenario is defined as a three-phase short-circuit fault, occurring at bus 29 and cleared after 100 ms by tripping line 28-29, indicated as fault 28-29*. Fig. 4 shows the machine swing patterns for this contingency, which is transient unstable. The CM group includes one machine M 9 , while the remaining machines are non-critical machines. The SIME analysis confirms that the system is transient unstable with t u = 420 ms and CCT = 70 ms.
From the full-system simulation, it is observed that M 1 and M 8 of area I are no longer coherent due to this contingency.
But the machines in areas III, IV, V are less disturbed and have nearly the same swing patterns. Thus, these machines are considered as a coherent group. When the dynamic reduction method is applied for areas III, IV, V, the full system is reduced to a small system with five machines, of which four retained machines (M 1 , M 8 −M 10 ) and one equivalent machine for the remaining six machines (M 2 −M 7 ). Fig. 5 shows the comparison of the swing angles of the critical machine and the equivalent machine, in the full system and in the reduced system separately. The simulation results show that the reduced system represents very well the original system. Table 3 shows the transient stability parameters of the OMIB model using the full and reduced systems. Fig. 6 shows the OMIB angle trajectories using the two systems separately. It can be seen that the reduced system represents well the transient stability of the full system, with AAE = 0.01. The reduction of simulation execution time by the reduced system is 46.4%, which is nearly equal to the proportion of machines in the two systems.
Using the transient stability parameters in Table 3, the TSC-OPF is carried out by simulating the full and the reduced systems separately. The results are shown in Tables 4 and 5.
From Tables 4 and 5, it is observed that the two approaches produce very similar results. The SIME analysis confirms that the two operating points both satisfy the transient stability conditions and have the same CCT = 100 ms. The reduction of TSC-OPF converging time using the reduced system simulation is 45.94%, which is nearly the same as the simulation time reduction. Case 2: A three-phase short-circuit fault at bus 7 is cleared after 350 ms by tripping line 6-7, indicated as fault 6-7*. The full-system simulation shows that the system is transient unstable with M 1 −M 9 loosing synchronism, as shown in Fig. 7a. Now, the CM group includes nine machines M 1 −M 9 , while M 10 is the non-critical machine. The dynamic reduction method is applied to the CM groups, which can be represented by two equivalent machines (M 1 , M 8 , M 9 ), (M 4 −M 7 ) and two retained machines (M 2 , M 3 ). Thus, the full system will become a reduced system with five machines, of which three retained machines (M 2 , M 3 , M 10 ) and two equivalent machines.   . Thus, the full system will become a reduced system with four machines, of which the retained machine M 10 and three equivalent machines. Case 6: A three-phase short-circuit fault at bus 10 is cleared after 250 ms by tripping line 10-11, indicated as fault 10-11*. The full-system simulation shows the system is transient unstable with M 1 -M 9 loosing synchronism, as shown in Fig. 7e. Now, the CM group includes nine machines M 1 -M 9 , while the NM group consists of one machine M 10 . The dynamic reduction method is applied for the CM groups, which is represented by two equivalent machines (M 1 , M 8 , M 9 ), (M 4 −M 7 ) and two retained machines M 2 , M 3 . Thus, the full system will become a reduced system with five machines, of which three retained machine M 2 , M 3 , M 10 and two equivalent machines. It is noted that, in the last three cases, although the nine machines M 1 -M 9 are all CMs, they are grouped in different ways. It is shown that the fault may change the coherency relationship between machines, especially the machines in the faulted area. Thus, before the dynamic reduction method is applied, a full-system simulation is required to verify the coherency relationship between machines. Table 6 shows the machine grouping information, the parameters of the OMIB model, the results and the execution time of the TSC-OPF procedure, using the full system and the reduced system separately. The machines in parenthesis in the third and fourth columns are the coherent machines which are reduced to an equivalent machine. In all cases, the TSC-OPF with dynamic reduction method gives similar results as with the full-system simulation and is more efficient.

IEEE 145-bus 50-machine system
In this section, we will demonstrate the efficiency of the proposed approach on the IEEE 145-bus, 50-machine benchmark system [32] provided by the PST software. The classical model is used for all the machines. The base-case operating conditions and the dynamic data of this system can be found in [28]. The limits of the bus voltage magnitude, 0.9 < V g < 1.1 are used for generator buses and 0.8 < V l < 1.2 for the remaining buses. The objective of TSC-OPF is to minimise the total cost associated with generation rescheduling as [6] min 0.1 ×      where P x,o is the initial active power of machine x. P x is the OPF results to be solved. N is the total number of machines. The software PST is first used to identify the coherent areas. By choosing C tol = 0.95, the system is divided into seven areas, as shown in Table 7. Note that area VII is a group of large, less disturbed generators although they are not strictly coherent.

Case 1:
The first contingency scenario is defined as a three-phase short-circuit fault, occurring at bus 7 and cleared after 130 ms by tripping line 6-7, indicated as fault 6-7*. Fig. 8 shows the machine swing patterns under this contingency, where M 2 and M 6 in area II loose synchronism.
Although M 2 and M 6 now are the critical machines, the full-system simulation shows that the remaining machines in area II still have same swing patterns as for the machines in areas I and III. Thus, these machines can be reduced to an equivalent machine. The coherent relationship between the machines in areas IV, V, VI and VII is not disturbed by the fault. The 50 machines now can be divided into five coherent groups and two retained generators as shown in Table 9. When the dynamic reduction is applied, the full system becomes a reduced system where each coherent group is simplified to a single machine. Fig. 9 shows the comparison of the swing patterns of the two retained machines 2, 6 in the full and the reduced systems. The good matching between the solid and the dot lines shows that the reduced system represents well the full system. The power response of the OMIB model in Fig. 10, using the full and the reduced systems separately, also confirms the validity of the reduced system. Table 8 gives the transient stability parameters of the OMIB for the full and the reduced systems separately, showing that the reduced system with seven machines represents well the full system with 50 machines. But the reduced system improves largely the simulation efficiency. Under the same conditions, the execution time of the reduced system is 5.6% that of the full system.
Using the transient stability parameters in Table 8, the TSC-OPF is carried out by including the full and the reduced systems simulation separately. The results in Table 10 show that the two approaches give nearly the same value. The swing angles of the 50 machines after the TSC-OPF procedure are plotted in Fig. 11, which shows that the system becomes stable. The SIME analysis shows that the optimal operating point satisfies the stability conditions with CCT = 130 ms.  Table 9. Case 3: A three-phase short-circuit fault at bus 101 is cleared after 280 ms by tripping line 69-101, indicated as fault 69-101*. The full-system simulation shows that the system is      Table 9.  Table 9. Table 10 shows the parameters of the OMIB model, the results and the execution time of the TSC-OPF procedure, using the full system and the reduced system separately. In all cases, the TSC-OPF with the reduced-system simulation gives as good results as with the full-system simulation but is much more efficient.
Comparing with the 10-machine system, the cases of the 50-machine system show that the proposed approach is more efficient for large systems. It is clear that the efficiency of the proposed approach depends on the number of coherent machines. If more machines are coherent, the smaller the reduced system will be, the more efficient the proposed method.

Conclusion
This paper presents an efficient TSC-OPF procedure by using an independent simulation based on the dynamic reduction method. In this new approach, based on the principle of slow coherency, the dynamics of the original system is represented by a reduced equivalent system. By using the SIME method based on the reduced-system simulation, a fast transient stability assessment is derived and is included into the TSC-OPF procedure. The new approach reduces the dynamic constraints in the TSC-OPF problem and thus reduces the computation burden. Two power systems test cases have shown the efficiency of the method. Since the new approach depends on the slow coherency areas in power systems, it presents more advantages for large power systems where the coherent machines are numerous and the slow coherency areas are large.