 Original research
 Open Access
 Published:
Parallel computing of multicontingency optimal power flow with transient stability constraints
Protection and Control of Modern Power Systems volume 3, Article number: 20 (2018)
Abstract
To deal with the high dimensionality and computational density of the Optimal Power Flow model with Transient Stability Constraints (OTS), a credible criterion to determine transient stability is proposed based on swing curves of generator rotor and the characteristics of transient stability. With this method, the swing curves of all generator rotors will be independent one another. Therefore, when a parallel computing approach based on the MATLAB parallel toolbox is used to handle multicontingency cases, the calculation speed is improved significantly. Finally, numerical simulations on three test systems including the NE39 system, the IEEE 300bus system, and 703bus systems, show the effectiveness of the proposed method in reducing the computing time of OTS calculation.
Introduction
With the development of smart grid, the operation of power system is close to its stability limit state frequently. In the event of a contingency or a large disturbance, largescale blackouts or even the whole system crashes would happen. However, Optimal Power Flow with Transient Stability Constraints (OTS) can integrate the economic and dynamic security of the system into the same framework, which has become an effective tool to get the safety of the system.
Lots of research work has been devoted to the classical OPF [1,2,3,4]. However, OTS has become a challenge from the model description to the design of algorithms [5, 6], due to the transient stability constraints of ordinary differential equations. In [7], the ordinary differential equations are differentiated into algebraic equations, while the transient stability constraints are discretized into inequality constraints. As a result, the OTS problem is solved by the conventional optimization method, but the number of variables and equations are increased significantly, resulting in long computation time and difficult convergence. However, references [5, 8] are two classical articles dealing with OTS problems practically. The theoretical development in [5] is that dynamic equations are converted to numerically equivalent algebraic equations. And in [8], OTS gets the same size of the OPF problem by an equivalent transformation from differential equation constraints to the initial value constraints.
However, in the Jacobian matrix of OTS, a large number of dynamic sensitivity equations must be dealt with in each iteration of the optimization process [9]. The dynamic sensitivity equation which describes the dynamics of power systems is a set of differential algebraic equations (DAE) with the same size and can only be solved by numerical integral methods with massive calculation.
The work of this paper is based on [8, 9], where a credible criterion is proposed based on swing curves of generator rotor and the characteristics of transient stability. Highperformance science and engineering computing [10], represented by parallel technology, has become an important approach of increasing productivity in various industries [11,12,13,14]. By using MATLAB parallel computing toolbox, the OTS problem is solved within reasonable computing resources with the computing speed reaching to the practical or even online level [15, 16] compared with the existing methods proposed in [17,18,19].
OTS model
Foundational model
Combined with the precontingence transient stability constraints into OPF, the model can be described as:

1)
Objective function is minimum active power loss of power system.
$$ \min \kern0.5em F(x)=\sum \limits_{i\in {S}_G}{P}_{Gi}\hbox{} \sum \limits_{i\in {S}_n}{P}_{Di} $$(1)
Where P_{Gi} refers to the active power output of i^{th} power source; S_{G} is active power sources set;S_{n} is the set of all nodes;P_{Di} is the active power load of i^{th} power source.

2)
Equality constraints.
There are two parts in this section, including power flow equations and swing equations of generator.

a.
Power flow equations.
$$ \left\{\begin{array}{l}P{}_{Gi}{P}_{Di}{V}_i\sum \limits_{j=1}^n{Y}_{ij}{V}_j\cos \left({\theta}_i{\theta}_j{\alpha}_{ij}\right)=0\\ {}Q{}_{Ri}{Q}_{Di}{V}_i\sum \limits_{j=1}^n{Y}_{ij}{V}_j\sin \left({\theta}_i{\theta}_j{\alpha}_{ij}\right)=0\end{array}\right.\kern0.5em i\in {S}_n $$(2)
Where Q_{Ri} refers to the reactive power output of i^{th} power source; Q_{Di} is reactive power load at bus i; V_{i}, θ_{i} are the magnitude and phase of voltage at bus i separately; Y_{ij}, α_{ij} are magnitude and phase of transfer admittance between buses i,j.

b.
Swing equations of generator.
For simplicity, a classic generator model is used in this paper. The difference equation is as follows:
Where δ_{i} refers to the rotor angle of i^{th} generator; ω_{i} is the rotor angular speed of i^{th} generator where ω_{0} defines the synchronous; D_{i} is damping value of i^{th} generator; M_{i} stands for the moment of inertia of i^{th} generator; P_{ei} is electromagnetic power of i^{th} generator:

3)
Inequality constraints.

a.
Operating constraints:
$$ \left\{\begin{array}{l}{\underline{P}}_{Gi}\le {P}_{Gi}\le {\overline{P}}_{Gi}\kern2.5em i\in {S}_G\\ {}{\underline{Q}}_{Ri}\le {Q}_{Ri}\le {\overline{Q}}_{Ri}\kern2em i\in {S}_R\\ {}{\underline{V}}_i\le {V}_i\le {\overline{V}}_i\kern3.25em i\in {S}_n\\ {}{\underline{P}}_{ij}\le {P}_{ij}\le {\overline{P}}_{ij}\kern2.75em i\in {S}_{CL}\end{array}\right. $$(5)

a.
Where \( {\overline{P}}_{Gi},{\underline{P}}_{Gi} \), \( {\overline{Q}}_{Ri},{\underline{Q}}_{Ri} \) are upper bound and lower bound of active and reactive power sources; \( {\overline{V}}_i,{\underline{V}}_i \), \( {\overline{P}}_{ij},{\underline{P}}_{ij} \) are upper bound and lower bound of voltage magnitude at bus i and power flow at line ij; S_{R}, S_{CL} are reactive power sources set and lines set.

b.
Stability constraints:
Stability criteria are as follows:
Where \( \overline{\delta},\underline{\delta} \) upper bound and lower bound of rotor relative swing are angle; δ_{COI} is inertia center angle, which is defined as follows:
Where M_{i} stands for the moment of inertia of the i^{th} generator.
Transient stability constraints
The transient process of power system is described by differential algebraic equations. The related transient stability preventive control problem is an optimization problem with differential algebraic equations. Its general form can be described as follows:
Formula (8) is the objective function. Formula (9) is a power flow equation. And formula (10) represents static security constraints, including the voltage limit, generator output limit, etc. Formula (11, 12) represent transient stability equality constraints. And formula (13) represents the transient stability inequality constraints under each contingency.
In this paper, based on the functional transformation technique, transient stability constraints (11–13) formulas are converted into a unique inequality equation.
The traditional timedomain simulation of transient stability analysis considers that the rotor of the unit and the center of inertia relative swing angle does not exceed a given limit angle, which can be considered the system run synchronously without loss of stability in the rolling process.
The formula is as follows:
Whether the system transient stability requirement is satisfied or not can be determined by calculating whether the sum of the excess area of each generator’s swing curve is less than or equal to zero.
The swing curve of the i^{th} generator is shown in Fig. 1.
To maintain transient stability, there is:
And the formula (11–13) is transformed into formula (15). However, when formula (15) is added as an inequality constraint, the convergence is not satisfactory with low computational efficiency.
Therefore, based on the contents of Fig. 1, this paper presents a simple and effective method to determine the transient stability of the system under the given power angel limit.
In the relative swing curve of all generator rotors, there must be at least one point deviated from the highest point. If the point satisfies formula (14), the transient stability condition of the system can be considered to be satisfied. As is shown in Fig. 2.
The t_{m} point of the generator 1 represents the highest point that deviates from the center of inertia in all the generators. Once the generator 1 is constrained, the other generators, including the generator 2 and the generator 3, naturally satisfy the constraint conditions.
Assuming that the point is the t_{m} time of the i^{th} generator, the constraint expression can be obtained:
Formula (16) can be used as the transient stability inequality constraint instead of the formula (14). The model after conversion is as follo6ws:
Parallel algorithm design
In order to simplify the calculation, this paper uses the classical model of multi machine system and uses the modern interior point algorithm to solve the formula (17–20). In the model, Jacobi and Hessian arrays of formula (18–19) can be obtained directly. However, the Jacobi and Hessian array of formula (20) can be obtained as follows:
Where \( {g}_{t_m} \) is an explicit expression for \( {x}_{t_m} \), so \( \partial {g}_{t_m}/\partial {x}_{t_m} \) and \( {\partial}^2{g}_{t_m}/\partial {x_{t_m}}^2 \) can be derived directly. While \( {x}_{t_m} \) is not an explicit expression for x_{0}, \( \partial {x}_{t_m}/\partial {x}_0 \) and \( {\partial}^2{x}_{t_m}/\partial {x_0}^2 \) can be calculated as follows:
Formula (21) can be solved by any method that can effectively solve differential equations and the fourthorder RungeKutta method is used in this paper. According to the formula (21–22), for each contingency, the swing curve of all the generator rotors can be calculated separately without affecting each other.
In solving the multicontingency problem, each contingency requires repeated calculations, which is the most timeconsuming place as well. As the size of the system expands, the time in this place even occupies more than 90% of the total computing time. Therefore, aiming at the main timeconsuming part of the program, a parallel algorithm is designed in this paper to allocate the tasks according to the expected multicontingency.
A diagram for the proposed approach is shown in Fig. 3. That is, the parallel loops can be decomposed and processed by client and worker modes. Client refers to a CPU kernel that performs and assignments parallel tasks. Worker refers to the other CPU cores that run parallel code, and all workers work at the same time to achieve accelerated computing.
The parallel architecture of this paper relies on MATLAB parallel computing toolbox. The key method is running parallel forloops on workers within contingencies.
And there are three main variables which should be classified before the parallel forloops.

1)
Loop Index Variable
The loop variable represents the number of times the loop is executed. The parallel strategy in this paper is used in the expected multicontingency.

2)
Sliced Variable
The large variable array contains all contingencies can be segmented into each worker by using sliced variable in the implementation of the parallel forloops program i and therefore the data transmission between worker and client is reduced. The array contains all the contingency information is actually a contingency variable in this paper.

3)
Broadcast Variable.
Data arrays contains nodal admittance matrix, generator operating data and parameters are classified as broadcast variable which will be shared by all the workers.
Multicontingency OTS problem parallel computing process is as follows:

1)
First of all, the conventional optimal power flow calculation is carried out without considering the transient stability constraints of formula (20), only the calculation of formula (17–19) are solved. Then taking the optimal solution of the conventional optimal power flow as the initial value, the transient stability analysis is carried out to determine whether the formula (20) is satisfied. If the optimal solution is satisfied, the calculation is terminated; otherwise, the next step is taken.

2)
The initial values of the parameters under each contingency are passed to each worker, and each worker calculates the Jacobian and Hessian matrix independently to complete the iterative calculation of the interior point method.

3)
In the iteration process, it’s necessary to judge repeatedly whether the solution of formula (20) and formula (17–19) are satisfied. If both satisfied, the calculation is terminated; otherwise, the calculation process above is repeated.
Numerical results
Test systems summary
In this section, the NE39 system, the IEEE 300bus system and 703bus systems are used for simulation from the aspects including computing time, speedup and efficiency. Furthermore, the simulations are run on 32 Cores 1.8 GHz with 16 GB of RAM. Set the simulation time T = 2s, the integral step size Δt = 0.01s or 0.02s, and the swing angle limit is ±100°. A numerical summary of the three systems is listed in Table 1.
Simulation analysis of single contingency OTS
Formula (16) is used as transient stability constraint to solve the OTS model (17–20). Test results of three system are shown in Tables 2, 3 and 4.
From Table 1 to Table 4, the calculation results of contingency lines of branch 8–9 in NE39 system, branch 5–1 in IEEE300 system and branch 681–337 in C703 system including the objective function and iterations is completely consistent with the conventional OPF calculation, which shows that the swing curve can meet the condition of angel limit ±100° under the transient stability analysis without adjustment when the conventional optimal power flow result is used as the initial value, also called ‘inactive contingency’.
In the meantime, the calculation results of contingency lines of branch 21–22 in NE39 system, branch 97–96 in IEEE300 system and branch 10–442 in C703 system has a larger objective function and more iterations compared with the conventional OPF calculation, indicating that the swing angle of generators reaches the upper limit when expected contingency are calculated, which means that the contingency must be unstable without constraints, also called ‘active contingency’. In order to meet the requirements of transient stability, economy is sacrificed properly, and more iteration is added to adjust the system operation mode meanwhile.
The swing curves of three test systems are shown in Figs. 4, 5 and 6.
Simulation analysis of multicontingency OTS
Test results of three systems are shown in Tables 5, 6 and 7.
From Table 5 to Table 7, iterations of the NE39 system and the IEEE300 system changes from 17 to 55 and from 26 to 56 respectively and iterations of the C703 system changes little when considering multicontingency. In particular, with the same number of iterations, the calculation time increases approximate linearly as the number of contingency increases, reflecting the good computational characteristics of the model in this paper.
According to the results of C703 system, when considering three contingency, the results of objective function is consistent with the results of the first contingency which shows that the first contingency plays a dominant position when considering the previous three contingency. The system adjustment mainly aims to the constraints of the first contingency which can meet the constraints of other contingency as well. However, when considering 6 contingency, the 4th contingency becomes the dominant contingency. Therefore, if the leading contingency in the sets can be found, a great amount of computation of OTS problem can be reduced effectively, which can be another research goal in the future.
The swing curves of three test systems are shown in Figs. 7, 8 and 9.
Simulation analysis of parallel computing
Performance index
Speedup refers to the ratio of computation time that the same task running on a singleprocessor and parallel processor system, which is used to measure the performance of parallel programs. In the field of parallel computing, the speedup ratio is defined as
Where T_{s} and T_{p} refers to the serial execution time and parallel execution time, respectively.
Speedup can reflect the overall performance of the program, but can’t describe the role of each processor, so the concept of parallel efficiency is proposed:
q is the core number involved in the calculation, S_{p} is the speedup ratio mentioned above.
Results and analysis
Benchmark results are shown in Table 8. Speedup ratio and parallel efficiency are shown in Table 9.
For all the test systems, it’s found that the speedup ratio increase with the number of contingency and the approach tends to perform better on large systems even up to 7.28 times. It is also observed that the acceleration is not proportional to the number of contingencies. This phenomenon mainly has the following two reasons:

1)
Communication Loss. As the size of the system expands, parallel tasks may increase the number of iterations and then the consumption of communication operation time of total computation time increased in a global iteration of the program’s global variables. As a result, the parallel efficiency decreased with the increase of contingency. In fact, in the same number of contingencies, the smaller the system size, the corresponding efficiency tends to be lower as shown in Table 9.

2)
Serial Calculation. The parallel method is used to solve the coefficient matrix of the modified equation merely, so the other process of the algorithm still uses the serial solution. In general, as the number of contingencies increases, the amount of data required for serial computation is increasing as well, which is another reason the speedup ratio tends to be saturated.
Nevertheless, the speedup ratio still shows a good upward trend. It is shown that the parallel method of this paper can adapt to the calculation and analysis of large power systems, which can effectively reduce the computing time of OTS.
Furthermore, compared with the existing approach in [17,18,19], for example, the computing time in 2 contingencies in IEEE 300 system is 68.264 s. And it’s 5.94 s in this paper.
In addition, compared with the existing approach, the proposed algorithm has the following advantages:

1)
A stronger acceleration
The method proposed in this paper has better acceleration effect, especially in larger systems or more contingencies, which shows excellent calculation performance.

2)
More contingencies can be calculated
This algorithm can handle even more than 32 contingencies. Theoretically, as long as the number of cores is sufficient, more contingencies can be processed at the same time and a higher speedup ratio can be achieved.

3)
Better convergence
The simulation data also shows that the algorithm has strong convergence for any number of contingency.

4)
Easy upgrading and extending
The parallel part of the algorithm is mainly based on parallel computing technology. In hardware, without changing the algorithm structure, the increase in the number of cores and the frequency can easily achieve the performance extension.
Conclusions
In this paper, the OTS model based on the criteria according to the swing curves of generator rotor and the characteristics of transient stability analysis is proposed. Taking into account the existing computing power, a parallel method based on MATLAB toolbox is used as well. Tests results in three different systems shown that the computing time of OTS is reduced efficiently without obvious changes on the number of iterations and the optimal solution.
The algorithm proposed in this paper is simple in constraint and good in convergence. The multicore processors parallel computing adopted in this paper is basically the same as the single CPU because of its hardware structure, which can improve the computing capacity without changing the hardware structure. Therefore, the method is versatile and scalable and is not limited to a specific algorithm or architecture.
Abbreviations
 DAE:

Differential algebraic equations
 OTS:

Optimal Power Flow model with Transient Stability Constraints
References
 1.
ElHawary, M. E. (1996). Optimal power flow: Solution technologies, requirement and challenges [R]. Piscataway: IEEE Tutoral Service, IEEE Service Center.
 2.
Tong, X. J., Wu, F. F., Zhang, Y. P., et al. (2007). A semismooth Newton method for solving optimal power flow. Journal of Industrial and Management Optimization, 3(3), 553–567.
 3.
Luo, K., Lin, M. G., & Tong, X. J. (2006). Decoupled semismooth Newton algorithm for optimal power flow problems. Control and Decision, 21(5), 580–584.
 4.
Luo, K., & Tong, X. J. (2006). New model and algorithm for solving the KKT system of optimal power flow. Control Theory and Applications, 23(2), 245–250.
 5.
Deqiang, G., Thomas, R. J., & Zimmerman, R. D. (2000). Stability constrained optimal power flow. IEEE Transactions on Power Systems, 15(2), 535–540.
 6.
Momoh, J. A., Koesslwer, R. J., & Bond, M. S. (1997). Challenges to optimal power flow. IEEE Transactions on Power Systems, 12(1), 444–447.
 7.
La Scala, M., Trovato, M., & Antonelli, C. (1998). Online dynamic preventive control: An algorithm for transient security dispatch [J]. IEEE Transactions on Power Systems, 13(2), 601–610.
 8.
Chen, L., Tada, Y., Okamoto, H., et al. (2001). Optimal operation solutions of power systems with transient stability constraints. IEEE Transactions on Circuits and Systems, 48(3), 327–339.
 9.
Ming bo, L., Yan, X., & Jie, W. (2003). Calculation of available transfer capability with transient stability constraints. Proceedings of the CSEE, 23(9), 28–33.
 10.
Vecchiola, C., Pandey, S., & Buyya, R. (2009). HighPerformance Cloud Computing: A View of Scientific Applications. International Symposium on Pervasive Systems, Algorithms, and Networks (pp. 4–16). IEEE Computer Society.
 11.
Gomez, A., & Betancourt, R. (1990). Implementation of the fast decoupled load flow on a vector computer [J]. IEEE Transactions on Power Apparatus and Systems, 5(3), 977–983.
 12.
Sasaki, H., Aoki, K., & Yokoyama, R. (1987). A parallel computation algorithm for static state estimation by means of matrix inversion lemma. IEEE Power Engineering Review, PER7(8), 40–41.
 13.
Li, Y., & McCalley, J. D. (2009). Decomposed SCOPF for improving efficiency. IEEE Transactions on Power Apparatus and Systems, 24(1), 494–495.
 14.
Abur, A. (1988). A parallel scheme for the forward/backward substitutions in solving sparse linear equations. IEEE Transactions on Power Apparatus and Systems, 3(4), 1471–1478.
 15.
Sutter, H., & Larus, J. (2005). Software and the concurrency revolution. Q focus: Multiprocessors, 3(7), 54–62.
 16.
J Y, C., & X H, M. (2015). Multiobjective optimal power flow considering transient stability based on parallel NSGAII. IEEE Transactions on Power Apparatus and Systems, 30(2), 857–866.
 17.
Jiang, Q., & Geng, G. (2010). A reducedspace interior point method for transient stability constrained optimal power flow. IEEE Transactions on Power Apparatus and Systems, 25(3), 1232–1240.
 18.
Mo, N., Zou, Z. Y., Chan, K. W., et al. (2007). Transient stability constrained optimal power flow using particle swarm optimisation. Iet Generation Transmission & Distribution, 1(3), 476–483.
 19.
Geng, G., Jiang, Q., & Sun, Y. (2016). Parallel transient stabilityconstrained optimal power flow using GPU as coprocessor[J]. IEEE Transactions on Smart Grid, PP(99), 1.
Funding
This work was supported by the National Natural Science Foundation of China (Grant No.51577085).
Author information
Affiliations
Contributions
YD Yang conceived and designed the study. YD Yang and AJ Song performed the experiments and simulations. YD Yang, AJ Song and H Liu wrote the paper. YD Yang, AJ Song, H Liu, ZJ Qin, J Deng and JJ Qi reviewed and edited the manuscript. All authors read and approve the final manuscript.
Corresponding author
Ethics declarations
Competing interests
The authors declare that they have no competing interests.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
About this article
Cite this article
Yang, Y., Song, A., Liu, H. et al. Parallel computing of multicontingency optimal power flow with transient stability constraints. Prot Control Mod Power Syst 3, 20 (2018). https://doi.org/10.1186/s416010180095z
Received:
Accepted:
Published:
Keywords
 Power system
 Transient stability
 Multicontingency
 Optimal power flow
 Parallel computing