Next Article in Journal
Comparison of Tender Criteria for Electric and Diesel Buses in Poland—Has the Ongoing Revolution in Urban Transport Been Overlooked?
Next Article in Special Issue
DC Admittance Model of VSCs for Stability Studies in VSC-HVDC Systems
Previous Article in Journal
Multi-Stage Bargaining of Smart Grid Energy Trading Based on Cooperative Game Theory
Previous Article in Special Issue
Overview of Various Voltage Control Technologies for Wind Turbines and AC/DC Connection Systems
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Unique Symbolic Factorization for Fast Contingency Analysis Using Full Newton–Raphson Method

1
Mathematic Computer and Engineering Department, University of Quebec at Rimouski, Rimouski, QC G5L 3A1, Canada
2
Green Teck Institute, Mohammed VI Polytechnic University, Benguerir 43150, Morocco
*
Author to whom correspondence should be addressed.
Energies 2023, 16(11), 4279; https://doi.org/10.3390/en16114279
Submission received: 3 April 2023 / Revised: 3 May 2023 / Accepted: 17 May 2023 / Published: 23 May 2023
(This article belongs to the Special Issue Advanced Electric Power System 2023)

Abstract

:
Contingency analysis plays an important role in assessing the static security of a network. Its purpose is to check whether a system can operate safely when some elements are out of service. In a real-time application, the computational time required to perform the calculation is paramount for operators to take immediate actions to prevent cascading outages. Therefore, the numerical performance of the contingency analysis is the main focus of this current research. In power flow calculation, when solving the network equations with a sparse matrix solver, most of the time is spent factorizing the Jacobian matrix. In terms of computation time, the symbolic factorization is the costliest operation in the LU (Lower-upper) factorization process. This paper proposes a novel method to perform the calculation with only one symbolic factorization using a full Newton–Raphson-based generic formulation and modular approach (GFMA). The symbolic factorization retained can be used during the iterations of any power flow contingency scenario. A computer study demonstrates that reusing the same symbolic factorization greatly reduces computation time and improves numerical performance. Power system security assessment under N-1 and N-2 contingency conditions is performed for the IEEE standard 54-bus and 108-bus to evaluate the numerical performance of the proposed method. A comparison with the conventional power flow method shows that the time required for the analysis is shortened considerably, with a minimum gain of 228%. The comparative analysis demonstrates that the proposed solution has better numerical performance for large-scale networks.

1. Introduction

Contingency analysis is regarded as an integral part of power system security analysis as it determines the integrity of the power system and verifies the post-contingency equilibrium state in terms of overloads and voltage deviations [1,2]. Its main purpose is to check whether the system can operate safely without any limit violation (i.e., equipment overloading, under- or over-voltage) after the occurrence of a contingency [3].
In N-1 contingency analysis, the main focus is to evaluate the effects of a single equipment outage on power system operating conditions once the single component (such as the loss of a generator, line, and transformer) has been retrieved from the network. Traditional N-1 contingency analysis performs numerous power-flow runs with a total number of N scenarios. After solving the power flow problem for each contingency scenario, the limits are verified for branches and nodes. The system is N-1 secure if, for all scenarios, the following constraints are fulfilled:
I I m a x U min U U max
The extension of N-1 contingency analysis considers simultaneous component outages (N-k contingencies). N-k security refers to the ability of the system to maintain secure operation after k components have failed in a system with N components [4,5]. For a large-scale network, strict N-k checking requires a large amount of power flow calculation. However, in real-time security assessment, the time window for system operators to analyze the problem and take corrective actions is quite limited. Therefore, the main challenge is to find an efficient method to significantly reduce the resolution time without compromising solution accuracy [6].
In contingency analysis, the power flow method used to perform the contingency analysis has a direct impact on the solution time. The Newton–Raphson power flow-based approach is generally preferred for its superior convergence characteristics and accuracy. The most computationally demanding part of this method is solving the linear equations with LU (Lower-upper) factorization at each iteration. In terms of computation time, the symbolic factorization is the costliest operation in the LU factorization process. This paper focuses on using the Newton–Raphson method to solve contingency analysis. The central question is whether it is possible to avoid the repetitive symbolic factorizations to effectively shorten the analysis time.
We will now highlight the strengths and weaknesses of the most important methods used to reduce the computational time of contingency analysis. This critical literature review will make evident the effectiveness of the proposed solution in simulating the contingency analysis of large-scale networks (Transmission and distribution systems).
The issue of improving computational speed in contingency analysis has been discussed in many works. In [7,8], the compensation theorem is applied to simulate the changes in the passive elements of the network without changing the Jacobian matrix structure. However, this method is well suited to applications involving linear network equations, which is not the case for analysis based fully on Newton’s method.
The main alternative to the compensation method for network matrix modification is to perform partial refactorization [9]. This method consists of updating the LDU (Lower-diagonal-upper) factors to reflect changes in some elements of the Jacobian matrix. The weakness of this method is that the refactorization effort is affected by the position of the modified matrix elements; thus, little or no savings can be obtained if one of the modified elements is near the top of the matrix.
The DC power flow method has been widely used to accelerate the computation of contingency analysis and improve its numerical performance [10]. The DC power flow simplifies the AC power flow to a linear circuit problem. This method is convergent and non-iterative but less accurate than AC power flow solutions.
A reduction in computational time can be achieved by reducing the number of contingency cases to be tried [11,12]. This approach involves ranking contingencies in descending order of severity index. The Contingencies can then be simulated, starting with the most severe and continuing until there is no overloading caused by the outage of branches. The main shortcoming of this method is that the contingency ranking by the index may introduce some errors due to the screening effect.
In [13], the graph theory was used to demonstrate that the base case symbolic factorization can be reused to solve contingency cases. The key assumption of this method is that there is no branch to be put in service during a contingency. This limitation prevents practical use in distribution systems where the introduction of new lines to reroute flows occurs in the most interesting contingencies.
Recent research in the area includes optimization-based methods [14,15,16,17,18,19,20]. These methods do not enumerate every contingency but attempt to find the most severe contingency to reduce the computational burden.
Table 1 demonstrates a summary of the reviewed classical methods for managing contingency analysis, with the major strengths and weaknesses of each one.
Additionally, the aforementioned methods have only addressed cases where the change in the sparsity pattern of the Jacobian matrix is caused by the change in network topology. However, in an automatic adjusting solution [21,22] a new symbolic factorization is not only required for a change involving network topology but also during Newton’s iterative process (due to control adjustment); hence, the symbolic factorization of the base case Jacobian matrix cannot be systematically reused in contingency cases.
In this paper, an original modeling technique has been proposed to avoid the occurrence of change in the sparsity pattern within the same power flow run (due to control adjustments) and after evaluating a new contingency scenario (due to a change in network topology). The proposed solution will overcome the limitations identified in classical methods (Table 1) and allow for the systematic and fast solving of contingency analysis. More specifically, we propose to perform the contingency analysis using only one symbolic factorization. The proposed approach is completely generic and can easily accommodate the simultaneous outages of k-branches in the network (i.e., N-k contingency analysis) with the possibility to include any remedial actions that consist of putting in new service elements.
This paper uses the power flow solver based on GFMA formulation to solve the power flow problem [23]. The algorithm of GFMA is based on the rigorous full Newton–Raphson approach and relies on the automatic adjustment technique, which provides an accurate and fast power flow solution. The high flexibility of this formulation offers the possibility of re-using the base case symbolic factorization of the Jacobian matrix to recompute the power flow of any other scenario. The basic idea behind the proposed approach is to imbue the GFMA power-flow algorithm with the concept of the dynamic parameters to preserve the same sparsity pattern of the Jacobian matrix. The dynamic parameters concept is the main contribution of this research paper.
This paper starts with a brief presentation of the GFMA method for solving power flow problems. The second part presents the modeling concept of dynamic parameters. In the third part, the dynamic parameters are introduced into the component models to perform a unique symbolic factorization in the contingency analysis. The last part presents illustrative simulation examples to evaluate saved time for relatively large transmission networks.

2. GFMA Formulation

The formulation of GFMA is recalled here to establish the basis for the following algorithms. GFMA provides a generic power flow formulation capable of handling complex component models and arbitrary network topology. In this method, each component of a system is modeled autonomously as a subsystem using a non-causal mathematical representation. The implicit representation of the system and the modular approach avoid many theoretical modeling complications, as the variables and load flow constraints can be arbitrarily defined. The GFMA formulation is generic, extensible, and can handle arbitrary component models without any known limitation.
Each component is described using the model representation in the form of Equation (1). For the convenience of the reader, matrices and vectors are indicated in bold type characters.
g k ( u k , y k , x k , λ k ) = 0
where:
  • uk are the model input variables.
  • yk are the model output variables.
  • xk are the model internal variables.
  • λk are the model variables.
Each model introduces m inputs, n outputs, and h internal state variables. These variables are expressed in real units and can be arbitrarily defined. The voltages are defined as inputs to establish the electrical connection with the connected bus and the currents as outputs to formulate Kirchhoff’s Current Law. The inputs and outputs represent the interface between the model and the system. The readers can refer to [23] for more details.
The component models are aligned along the diagonal to form the equipment Jacobian matrix. This block diagonal matrix is augmented with the linking equations to establish the connection between components. These extra-equations are external to the components model and have a linear form as follows:
i , j ( c i , k . y i , k + d j , l . u j , l ) = 0
where y i , k is the i t h output variable of the k t h model; u j , l is the j t h input variable of the l t h model; c i , k and d j , n are the coefficients associated with y i , k and u j , l , respectively.
Herein, the system is initialized using the single-iteration FP solution [24]. The estimation of the initial guess of the state vector is obtained by converting (1) to the linearized version and solving the resulting linear system of equations in the form of
A . x = b
A comprehensive example is shown in Figure 1 to explore the structure of the Jacobian matrix. As can be seen from Figure 2, the upper and lower block of the Jacobian matrix represent the network equipment and linking equations, respectively. The obtained Jacobian matrix is sparse and unsymmetrical. It is important to notice that the system buses do not contribute to the Jacobian matrix of the system with additional nonzero elements (Figure 2).

3. Dynamic Parameters Approach

A sparse matrix solver provides efficient numerical solutions to solve linear equations. The KLU solver is suitable for solving unsymmetrical and sparse linear systems [25] and is considered one of the most efficient packages designed for electrical circuits. Roughly speaking, the KLU solver performs two steps in sequence:
1.
Symbolic factorization: performs optimal permutations and pre-ordering to reduce the fill-ins (number of non-zero elements).
2.
Numerical factorization: Computes the factorization sub-matrices L and U.
The change in the Jacobian matrix structure is due to the presence of conditional statements (“if-else”) in the code corresponding to different operating conditions. These conditions are evaluated, and the appropriate constraint equation is executed (see Equation (4)).
g = f 1 i f C o n d i t i o n ( 1 ) g = f i e l s e   i f C o n d i t i o n ( i ) g = f n e l s e C o n d i t i o n ( n )
The equation can be formulated as one compact mismatch equation.
g = λ 0 f 0 + i = 1 i = n λ i f i
where:
  • fi and λi are the non-linear constraint equation and the dynamic parameter for the ith condition, respectively.
  • f0 and λ0 are the linear constraint equation and the dynamic parameter used for the ini-tialization. It is worth noting that the term λ0f0 can be dropped from Equation (5) if the model is linear.
The dynamic parameters approach consists of updating the parameters λ i from iteration to iteration to enforce the appropriate mismatch equation. Specifically, when a new condition is flagged, the parameter associated with the enforced constraint equation is set to one, while the other parameters are all set to zero. This modeling technique prevents a change in the sparsity pattern of the Jacobian matrix since the same mathematical expression holds for all operating conditions.
Equations (4) and (5) are mathematically equivalent but numerically different. The main advantage of formulation (5) resides in its high computational performance since the repetitive symbolic factorizations are systematically avoided.
The dynamic parameters approach has been successfully applied to the power flow problem and will be extended here for contingency analysis. It will be demonstrated in the next section that the use of this modeling technique finds useful applications in network problems involving a change in network topology.

4. Outage Modeling

In this section, the dynamic parameters approach presented in [23] is extended to avoid the symbolic factorizations related to the change in network topology. Herein, the constraint equation enforcing the outage of the equipment ( f c ) and its associated parameter ( λ c ) are incorporated in the original formulation as follows:
g = λ c f c + ( 1 λ c ) . ( λ 0 f 0 + i = 1 i = n λ i f i )
The first term of Equation (6) prevents change in the nonzero pattern of the Jacobian matrix caused by the change in the network topology. To simulate a new contingency scenario, the parameter λ c associated with the localized equipment outage must be set to one. For the base case scenario, the same parameter is set to zero for all elements in the network.
A flowchart depicting the updates in the dynamic parameters in the contingency algorithm is shown in Figure 3.
From a computational point of view, the constraint equation f c is considered a judicious choice if the following conditions are fulfilled:
  • Condition 1: Ideally, the constraint equation fc should not introduce additional non-zero elements in the Jacobian matrix other than those initially generated from fi.
  • Condition 2: A linear constraint equation is preferred, as the Newton–Raphson converges more robustly due to its more linear formulation.
In the following sections, the equations for each model are developed using a compact mismatch equation according to the dynamic parameter approach. For the sake of clarity, the model equations are expressed in complex form, and the partial derivatives are not detailed here, as they can be obtained using the symbolical computation [26].

4.1. Line Outage

The outage of the line is represented by a linear constraint equation forcing the entering currents to zero. For the sake of clarity, the line model of Figure 4 is presented as an example.
The vector mismatch equation of the line formulated with the dynamic parameters is given by
g = λ c . I 1 I 2 f c + ( 1 λ c ) . ( I 1 I 2 Y . V 1 V 2 f 1 )
It should be noted that the linear constraint equation f 0 used for the initialization is not required here since the line model is inherently linear.

4.2. Generator Outage

The synchronous generator (SG) is modeled as variable susceptance and conductance ( G and B ) [23]. The latter two ( G and B ) are defined in the model as internal state variables which are solved to satisfy the power flow constraints of the machine. This model is numerically more efficient and converges faster than the conventionnel internal voltage behind reactance [27,28].
The SG state variables added to the system’s variables are as follows:
X S G = [ V k R , a , V k R , b , V k R . c V k I , a , V k I . b , V k I , c i n p u t s , G , B I n t e r n a l , I k R , a , I k R , b , I k R . c I k I , a , I k I . b , I k I , c o u t p u t s ]
This model introduces eight unknown variables (six for output currents and two for internal variables) required for the solvability of the same number of equations.
For a machine with PV constraint, the mismatch equation written with dynamic parameters is given by
g = λ c . G B f c + ( 1 λ c ) . ( λ 0 . G G 0 B B 0 f 0 + λ 1 . P G P d e s V G V d e s f 1 + λ 2 . 0 Q G Q min f 2 + λ 3 . 0 Q G Q max f 3 )
where the active and reactive power is expressed as a function of the state variables
P G + j Q G = I a b c * V a b c
The sign indicates the complex-conjugate operator.
This mismatch equation can capture all possible operating conditions of the machine, including the initialization and equipment outage.
In the initialization stage, the parameter λ 0 is set to one while all other parameters are set to zero, and the non-linear constraint equations of the generator are reverted to the linearized version. The linear constraint equation f 0 is specifically introduced in the mismatch equation to estimate the initial value of the susceptance and conductance. For a machine with PV constraint, G 0 and B 0 are estimated, assuming a nominal power factor of the machine and the desired active power.
G 0 + j B 0 = c o n j ( P d e s + j Q d e s ) V L L n o m 2
The parameters λ 1 , λ 2 , λ 3 are dynamically updated during the iterative power flow process to flag all operating conditions of the generator. This allows the machine to switch from PV to PQ when its reactive power limit is violated.
The outage of the generator is represented by an equation forcing the susceptance and conductance to zero. This constraint equation ( f c ) is linear and does not introduce any additional nonzero elements in the mismatch equation other than those generated from f 0 , f 1 , f 2 , and f 3 . The generator is considered effectively out of service if the parameter is set to one.
Six equations are also introduced from the equation relating the terminal voltages and currents.
I a b c = ( A Y 012 A 1 ) . V a b c
where A is the Fortescue’s transformation matrix; V a b c and I a b c are the vector of generator voltages and current in phase domain, respectively.
Y 012 = Y 0 0 0 0 G + j B 0 0 0 Y 2
With the introduction of the dynamic parameters, the constraint Equation (13) becomes
I a b c = A Y 0 0 0 0 λ c . ( G 0 + j B 0 ) + ( 1 λ c ) . ( λ 0 . ( G 0 + j B 0 ) + ( λ 1 + λ 2 + λ 3 ) . ( G + j B ) ) 0 0 0 Y 2 A 1 V a b c

4.3. Transformer Outage

In this section, the transformer outage model is implemented in a three-phase power-flow solution method based on the GFMA formulation. Therein, an ideal transformer connected between nodes k and n is modeled by describing the voltage and current equations relating the primary and secondary windings. The node m is inserted to include the three-phase impedance transformer. This solution is generic and can be used to handle arbitrary transformer connections. The transformer model depicted in Figure 5 is used to derive the constraint equations in the phase domain for an arbitrary three-phase transformer connection.
The constraint equations for the transformer are given by
g = λ c . I m a b c I n a b c f c + ( 1 λ c ) . ( I m a b c I n a b c Y t r Y t r Y t r Y t r V m a b c V n a b c f 1 ) V m a b c D ( r a t i o ) . V k a b c I k a b c + D T ( r a t i o ) . I m a b c
where Y t r is the three-phase transformer admittance matrix; r a t i o is the turn ratio of the master transformer; D is the dependency that is expressed as a function of turn ratio; and the sign T indicates the transpose matrix operator.
The transformer outage is simulated by forcing the output currents at both sides of the transformer to zero. This is achieved by setting the parameter λ c to one.

4.4. Switching Device Outage

For an ideal switch connected between nodes k and m , the model equation of the switch device depicted in Figure 6 is given by
g = λ c . I k m f c + ( 1 λ c ) . ( V k V m f 1 )
In this equation, the switch is considered in-service if λ c is set to zero (i.e., closed position). The same equation can be used to simulate the outage of the ideal switch by setting λ c to one (i.e., opened position).

4.5. Load Shedding

An optimal load-shedding scheme is necessary under contingency conditions to prevent cascade outages and complete black-out [29]. The load is formulated using the current-mismatch equations [30,31]. The load mismatch equation is given by
g = λ c . I k m f c + ( 1 λ c ) . ( λ 0 . ( I k m Y L ( V k V m ) f 0 ) + λ 1 . ( I k m P d e s j Q d e s V k * V m * ) f 1 )
where V k and V m are the voltages at nodes k and m , respectively; I k m is the injected load current; and Y L is the estimated load admittance used for the initialization.

5. General Aspects of the Algorithm

With the introduction of the dynamic parameters in the formulation of mismatch equations, the sparsity pattern of the Jacobian matrix corresponding to the base case scenario becomes identical to that of any contingency case. The unique symbolic factorization retained in all scenarios of contingency analysis corresponds to that performed for the system initialization of the base case scenario.
In the base case scenario, all the network equipment is initially in service. However, any equipment that was initially out of service but expected to be connected as part of remedial action must be included in the base case scenario and modeled as an out-of-service element. From a modeling point of view, adding and removing an element from the system are treated similarly.
Simultaneous outages (N-2) are treated as single outages (N-1). Both can reuse the symbolic factorization of the base case scenario.
The main drawback of the proposed approach is that the nonzero elements and the size of the Jacobian matrix in the contingency cases are slightly overestimated. To illustrate this point, let us consider an outage of a transmission line that is modeled using dynamic parameters. The Jacobian elements are obtained by replacing the dynamic parameter λ c in Equation (7) with one and then deriving the partial derivatives with respect to the state variables. Expressed in real form, the Jacobian elements of the positive sequence model of the line are given by
g I 1 R , I = 1 g I 2 R , I = 1 g V 1 R , I = 0 g V 2 R , I = 0
As can be seen from Equation (18), eight nonzero elements contribute to removing a branch element from the system. These elements will exist only if the symbolic factorization of the base case scenario is reused for the contingency scenario. This is evident from the fact that all elements in the base case are initially considered in service. It is worth noting that the KLU solver does not drop numerically zero entries from its sparse matrix; therefore, all the entries that are present in the data structure of the Jacobian matrix are not dropped, even though they are numerically zero.
The contingency analysis is performed with the algorithm depicted in Figure 7.

6. Validation

The proposed solution has been tested for the IEEE 57-bus and the IEEE 118-bus test systems. These networks are used to evaluate the numerical performance of the proposed solution. The objective is to compare the computational cost of the proposed solution (GFMA) with the contingency analysis method implemented in CYME power engineering software. This conventional method consists of performing extensive simulations and the strict checking of all possible contingencies using a MANA (Modified-augmented-nodal-analysis)-based power flow solution [32,33]. MANA and GFMA are two different formulations, but both are implemented in a-b-c reference frame and solved using the full-Newton algorithm and KLU solver.
Herein, the simulations are performed for the N-1 and N-2 security criteria. The N-1 criterion considers a unique contingency that is cleared by the primary protection, in this case, the faulted element is disconnected from the system without any further impact on the system. However, the N-2 criterion considers that the first outage in the system results in a second component failure. It is worth noting that the N-k security criterion leads to C k N possible contingencies ( C k N = N ! k ! N k ! ).
To evaluate the resolution time of the contingency analysis, timers were inserted into the code sections, where the symbolic and numerical factorizations are performed. The normalized CPU time with respect to the proposed solution is presented at the end of this section, and a comparison between the numerical efficiency of these two methods serves as a closing statement.
All algorithms have been programmed and executed using the MATLAB Platform, which runs on a dedicated machine (3.4 GHz i7-2600 CPU with 16 GB of RAM). A tolerance of 0.01% on the voltage mismatch is used as the stopping criterion for each power flow run.

6.1. Case-1: IEEE 57-Bus Test System

The IEEE 57-bus system shown in Figure 8 is used to test the numerical performance of the proposed method and to assess the gain in computation time. The network contains 57 buses, 63 lines, 7 generators, 15 two-winding transformers, and 42 loads. The total number of components that can potentially fail in the system is 85 ( N = 63 + 7 + 15 ). Contingency is carried out using the GFMA formulation proposed in this paper and MANA formulation implemented in CYME software.
The strict N-l of the IEEE 57-bus system needs 85 ( C 1 85 ) outage analyses. Table 2 shows that both methods require the same number of numerical factorizations; however, the number of symbolic factorizations required in the MANA is 190, while the proposed solution performs the same calculation with only one symbolic factorization. Moreover, when the nonzero pattern is already known, the depth-first search used in Gilbert/Peierls method can be skipped. This means that numerical factorization becomes computationally less expensive when the sparsity pattern is not changed. Therefore, the cost of the numerical factorization is smaller for GFMA than MANA. In this case, the GFMA is about 2.30 faster than MANA (Table 2).
For N-2 contingency analysis, the IEEE 57-bus system needs 3570 scenarios ( C 2 85 ). As expected, the CPU time is significantly reduced since one symbolic factorization is performed to evaluate this large number of outage scenarios. In this case, the GFMA is about 2.28 faster than MANA (Table 3). Although the performance ratio depends on the system size and the number of controller devices, this case study clearly demonstrates that the time required for the analysis is shortened considerably with the use of the proposed method.
The sparsity pattern of the Jacobian matrix of the base case scenario has a dimension of 3078 X 3078 with 16,242 nonzero elements (see Figure 9). The matrix structure of the base case scenario is identical to any contingency case.

6.2. Case-2: IEEE 108-Bus Test System

This IEEE 108-bus test system was used in various studies for contingency evaluation and voltage security assessment. This system contains 19 generators, 35 synchronous condensers, 177 lines, 9 transformers, and 91 loads (Figure 10). The total number of components is 240 ( N = 19 + 35 + 177 + 9 ). The strict N-l and N-2 checking of the IEEE 108-bus system needs 240 and 3570 outage analyses, respectively.
The comparative analysis presented in Table 2 and Table 3 leads to the same findings as those observed in the IEEE 54-bus system. Once more, the proposed solution proves to be faster than MANA, with a significant reduction in resolution time.
In addition, the proposed method performs better for larger systems and its efficiency increases as the system increases in size. This can be observed from the results of Table 2 and Table 3, where the detailed CPU time is provided for both networks (IEEE-54 and IEEE-108 systems). This is explained by the fact that the number of symbolic factorizations required to cover all contingencies increases drastically with the size of the network.

7. Conclusions

The proposed modeling technique presented in this paper exploits the flexibility of the GFMA formulation to avoid the repetitive symbolic factorizations required during the iterative power flow process and after the modification of network topology. Herein, the contingency analysis is performed with only one symbolic factorization using a Newton-like method and KLU sparse matrix solver.
The proposed approach is completely generic and can be applied to simulate simultaneous outages of k-branches elements in the network and any remedial actions that consist of putting in service new elements. The proposed solution is generic and can be applied to the strict N-k checking without any known limitation.
The new method has been tested on IEEE 57-bus and IEEE 108-bus test systems, with all the details represented and simulated. The power system security assessment under N-1 and N-2 contingency conditions demonstrates that a significant computational gain has been achieved by applying the concept of dynamic parameters. The results also showed that the method performs better for larger systems and that its efficiency increases as the system increases in size.

Author Contributions

Conceptualization, H.B.; methodology, H.B.; software, H.B.; validation, H.B.; formal analysis, H.B.; investigation, H.B.; resources, H.B.; data curation, H.B.; writing—original draft preparation, H.B.; writing—review and editing, H.B. and A.C.; visualization, H.B.; supervision, A.C. and A.E.O.; project administration, A.C. and A.E.O.; funding acquisition, A.E.O. All authors have read and agreed to the published version of the manuscript.

Funding

Financial support from the University of Quebec at Rimouski under grant AEO-760700.

Data Availability Statement

Not applicable.

Acknowledgments

The financial support from the University of Quebec at Rimouski under grant AEO-760700 is gratefully acknowledged.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Kip, M.; Lei, W.; Prabha, K. Power system security assessment. IEEE Power Energy Mag. 2004, 2, 30–39. [Google Scholar]
  2. Mohammed, S.; William, T.; Yong, F. Impact of Security on Power Systems Operation. Proc. IEEE 2005, 93, 2013–2025. [Google Scholar]
  3. Bulat, H.; Franković, D.; Vlahinić, S. Enhanced Contingency Analysis—A Power System Operator Tool. Energies 2021, 14, 923. [Google Scholar] [CrossRef]
  4. Liana, C.; Xuan, L.; Zuyi, L. Preventive Mitigation Strategy for the Hidden N-k Line Contingencies in Power Systems. IEEE Trans. Reliab. 2018, 67, 1060–1070. [Google Scholar]
  5. Babalola, A.; Belkacemi, R.; Zarrabian, S. Real-time cascading failures prevention for multiple contingencies in smart grids through a multi-agent system. IEEE Trans. Smart Grid 2018, 9, 373–385. [Google Scholar] [CrossRef]
  6. Hassan, H.; Mehrdad, T.; Pedram, S. Accurate, simultaneous and Real-Time screening of N-1, N-k, and N-1-1 contingencies. Int. J. Electr. Power Energy Syst. 2022, 136, 107592. [Google Scholar] [CrossRef]
  7. William, F.T. Compensation Methods for Network Solutions by Optimally Ordered Triangular Factorization. IEEE Trans. Power Appar. Syst. 1972, PAS-91, 123–127. [Google Scholar]
  8. Ongun, A.; Brian, S.; William, F.T. Sparsity-Oriented Compensation Methods for Modified Network Solutions. IEEE Trans. Power Appar. Syst. 1983, PAS-102, 1050–1060. [Google Scholar]
  9. Sherman, M.C.; Vladimir, B. Partial Matrix Refactorization. IEEE Trans. Power Syst. 1986, 1, 193–199. [Google Scholar]
  10. John, K.; Felix, W. Analysis of linearized decoupled power flow approximations for steady-state security assessment. IEEE Trans. Circuits Syst. 1984, 31, 623–636. [Google Scholar]
  11. Mikolinnas, T.A.; Bruce, W. An Advanced Contingency Selection Algorithm. IEEE Power Eng. Rev. 1981, 2, 25–26. [Google Scholar] [CrossRef]
  12. Halpin, T.F.; Robert, F.; Fink, R. Analysis of automatic contingency selection algorithms. IEEE Power Eng. Rev. 1984, 5, 29–30. [Google Scholar]
  13. Elliott, M.C.; Jingjin, W.; Yiting, Z.; Renchang, D.; Guangyi, L. Symbolic Factorization Re-utilization for Contingency Analysis. In Proceedings of the IEEE Power & Energy Society General Meeting (PESGM), Portland, OR, USA, 5–10 August 2018; pp. 1–5. [Google Scholar]
  14. Pudi, S.; Sanjeeb, M. Power system contingency ranking using Newton Raphson load flow method. In Proceedings of the Annual IEEE India Conference (INDICON), Mumbai, India, 13–15 December 2013; pp. 1–4. [Google Scholar]
  15. Rocco, C.M.; Ramirez-Marquez, J.E.; Salazar, D.E.; Yajure, C. Assessing the Vulnerability of a Power System Through a Multiple Objective Contingency Screening Approach. IEEE Trans. Reliab. 2011, 60, 394–403. [Google Scholar] [CrossRef]
  16. Ding, T.; Li, C.; Yan, C.; Li, F.; Bie, Z. Bilevel Optimization Model for Risk Assessment and Contingency Ranking in Transmission System Reliability Evaluation. IEEE Trans. Power Syst. 2017, 32, 3803–3813. [Google Scholar] [CrossRef]
  17. Street, A.; Oliveira, F.; Arroyo, J.M. Contingency-Constrained Unit Commitment with n–K Security Criterion: A Robust Optimization Approach. IEEE Trans. Power Syst. 2011, 26, 1581–1590. [Google Scholar] [CrossRef]
  18. Wang, Q.; Watson, J.; Guan, Y. Two-stage robust optimization for N-k contingency-constrained unit commitment. IEEE Trans. Power Syst. 2013, 28, 2366–2375. [Google Scholar] [CrossRef]
  19. Levitin, G. Optimal Defense Strategy Against Intentional Attacks. IEEE Trans. Reliab. 2007, 56, 148–157. [Google Scholar] [CrossRef]
  20. Levitin, G.; Hausken, K. Redundancy vs. Protection vs. False Targets for Systems Under Attack. IEEE Trans. Reliab. 2009, 58, 58–68. [Google Scholar] [CrossRef]
  21. Stott, B. Review of load-flow calculation methods. Proc. IEEE 1974, 62, 916–929. [Google Scholar] [CrossRef]
  22. Norris, M.P.; Scott, M. Automatic Adjustment of Transformer and Phase-Shifter Taps in the Newton Power Flow. IEEE Trans. Power Appar. Syst. 1971, PAS-90, 103–108. [Google Scholar]
  23. Hakim, B.; Ahmed, C.; Abderrazak, E. A generic three-phase power flow formulation for flexible modeling and fast solving for large-scale unbalanced networks. Int. J. Electr. Power Energy Syst. 2023, 148. [Google Scholar]
  24. Wilsun, X.; Hermann, W.D.; José, R.M. Generalized three-phase power flow method for the initialization of EMTP simulations. In Proceedings of the POWERCON ’98. 1998 International Conference on Power System Technology. Proceedings (Cat. No.98EX151), Beijing, China, 18–21 August 1998; Volume 2, pp. 875–879. [Google Scholar]
  25. Timothy, A.D.; Eknathan, P.N. Algorithm 907: KLU, a direct sparse solver for circuit simulation problems. ACM Trans. Math. Softw. 2010, 37, 3. [Google Scholar]
  26. Canizares, C.A. Applications of symbolic computation to power system analysis and teaching. In Proceedings of the 2006 PES Power Systems Conference and Exposition, Atlanta, GA, USA, 29 October–1 November 2006; p. 361. [Google Scholar]
  27. Chen, T.H.; Chen, M.S.; Inoue, T.; Kotas, P.; Chebli, E.A. Three-phase cogenerator and transformer models for distribution system analysis. IEEE Trans. Power Deliv. 1991, 6, 1671–1681. [Google Scholar] [CrossRef]
  28. Tatianna, B.; Penido, D.R.R.; Araujo, L.R. A new synchronous DG model for unbalanced multiphase power flow studies. IEEE Trans. Power Syst. 2020, 35, 803–813. [Google Scholar]
  29. Bhuvanagiri, R.; Mohan, K.; Surya, V. Priority Based Optimal Load Shedding in a Power System Network under Contingency Conditions. In Proceedings of the International Conference for Advancement in Technology (ICONAT), Goa, India, 21–22 January 2022; pp. 1–5. [Google Scholar]
  30. Sereeter, B.; Vuik, K.; Witteveen, C. Newton Power Flow Methods for Unbalanced Three-Phase Distribution Networks. Energies 2017, 10, 1658. [Google Scholar] [CrossRef]
  31. Garcia, P.; Pereira, J.; Carneiro, S.; Da Costa, V.; Martins, N. Three-phase power flow calculations using the current injection method. IEEE Trans. Power Syst. 2000, 15, 508–514. [Google Scholar] [CrossRef]
  32. Ilhan, K.; Jean, M.; Ulas, K.; Gurkan, S.; Omar, S. Multiphase Load-Flow Solution for Large-Scale Distribution Systems Using MANA. IEEE Trans. Power Deliv. 2014, 29, 908–915. [Google Scholar]
  33. Jean, M. Régimes transitoires électromagnétiques: Simulation. In Réseaux Électriques et Applications; TECHNIQUES DE L’INGENIEUR: Saint-Denis, France, 2015. [Google Scholar]
Figure 1. Test system example.
Figure 1. Test system example.
Energies 16 04279 g001
Figure 2. Matrix structure.
Figure 2. Matrix structure.
Energies 16 04279 g002
Figure 3. Updates in the dynamic parameters in the contingency case power flow solution.
Figure 3. Updates in the dynamic parameters in the contingency case power flow solution.
Energies 16 04279 g003
Figure 4. Line model.
Figure 4. Line model.
Energies 16 04279 g004
Figure 5. Three-phase transformer.
Figure 5. Three-phase transformer.
Energies 16 04279 g005
Figure 6. Switch device.
Figure 6. Switch device.
Energies 16 04279 g006
Figure 7. Contingency algorithm.
Figure 7. Contingency algorithm.
Energies 16 04279 g007
Figure 8. IEEE 57-bus System.
Figure 8. IEEE 57-bus System.
Energies 16 04279 g008
Figure 9. Sparsity pattern of the Jacobian matrix of the base case scenario (IEEE 57-bus System).
Figure 9. Sparsity pattern of the Jacobian matrix of the base case scenario (IEEE 57-bus System).
Energies 16 04279 g009
Figure 10. IEEE 118-bus test system.
Figure 10. IEEE 118-bus test system.
Energies 16 04279 g010
Table 1. Summary of the classical methods of contingency analysis.
Table 1. Summary of the classical methods of contingency analysis.
MethodStrengthsWeaknessesReference
Compensation method
  • Fast linearized power flow outage analysis.
  • Requires less computer storage.
  • Include the assessment of simultaneous outages of k-branches (N-k contingencies).
  • Nonlinear networks are not supported.
  • Simulate only the outage of passive elements.
[7,8]
Partial refactorization
  • Factorization is performed only on the affected part of the matrix.
  • Numerical performance depends on the position of the modified elements.
  • Numerical performance decreases if the number of changes in the matrix is small.
[9]
DC power flow
  • Very fast.
  • Non-iterative approach.
  • Risk of divergence.
  • Less accurate.
[10]
Severity index
  • Number of scenarios to be evaluated is greatly reduced.
  • Errors due to the screening effect.
  • Do not enumerate every contingency.
[11,12]
Graph theory
  • Possibility of re-using the base case symbolic factorization.
  • Fast.
  • Unsuited to simulate remedial actions.
  • Efficient for contingencies involving only a few branch elements.
[13]
Optimization methods
  • Reduction in CPU time.
  • Remedial actions can be included in the optimization function.
  • Include the assessment of simultaneous outages of k-branches.
  • Do not enumerate every contingency.
  • High Complexity.
[14,15,16,17,18,19,20]
Table 2. CPU Timings for the N-1 contingency analysis using MANA and GFMA.
Table 2. CPU Timings for the N-1 contingency analysis using MANA and GFMA.
IEEE 57-Bus SystemIEEE 108-Bus System
Proposed Solution
(GFMA)
CYME
(MANA)
Proposed Solution
(GFMA)
CYME
(MANA)
Number of Scenarios85240
Reading and Treatment of network fileNot includedNot included
Build and Update Jacobian Matrix 1.84 s1.76 s5.57 s5.13 s
Symbolic FactorizationsSymbolic Factorizations
Number of Symbolic Factorizations11901531
Timing of Symbolic Factorizations0.010 s2.09 s0.013 s6.90 s
Numerical FactorizationsNumerical Factorizations
Number of Numerical Factorizations425 *4251177 *1177
Timing of Numerical Factorizations1.06 * s2.97 s3.01 s *8.43 s
Solve0.084 s0.082 s0.12 s0.10 s
Total CPU Time2.99 s
(100%)
6.91 s
(230%)
4.24 s
(100%)
20.56 s
(484%)
* Numerical Refactorization.
Table 3. CPU Timings for the N-2 contingency analysis using the MANA and GFMA.
Table 3. CPU Timings for the N-2 contingency analysis using the MANA and GFMA.
IEEE 57-Bus SystemIEEE 108-Bus System
Proposed Solution
(GFMA)
CYME
(MANA)
Proposed Solution
(GFMA)
CYME
(MANA)
Number of Scenarios357028,680
Reading and Treatment of network fileNot includedNot included
Build and Update Jacobian Matrix 76.28 s72.46 s687.40 s651.43 s
Symbolic FactorizationsSymbolic Factorizations
Number of Symbolic Factorizations18925171,690
Timing of Symbolic Factorizations0.011 s93.71 s0.013 s931.97 s
Numerical FactorizationsNumerical Factorizations
Number of Numerical Factorizations16,950 *16,950136,170 *136,170
Timing of Numerical Factorizations42.37 s *110.17 s *267.00 s *695.00 s
Solve3.65 s3.54 s31.86 s29.45 s
Total CPU Time122.31 s
(100%)
279.89 s
(228%)
986.27 s
(100%)
2308.00 s
(234%)
* Numerical Refactorization.
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Bennani, H.; Chebak, A.; El Ouafi, A. Unique Symbolic Factorization for Fast Contingency Analysis Using Full Newton–Raphson Method. Energies 2023, 16, 4279. https://doi.org/10.3390/en16114279

AMA Style

Bennani H, Chebak A, El Ouafi A. Unique Symbolic Factorization for Fast Contingency Analysis Using Full Newton–Raphson Method. Energies. 2023; 16(11):4279. https://doi.org/10.3390/en16114279

Chicago/Turabian Style

Bennani, Hakim, Ahmed Chebak, and Abderrazak El Ouafi. 2023. "Unique Symbolic Factorization for Fast Contingency Analysis Using Full Newton–Raphson Method" Energies 16, no. 11: 4279. https://doi.org/10.3390/en16114279

APA Style

Bennani, H., Chebak, A., & El Ouafi, A. (2023). Unique Symbolic Factorization for Fast Contingency Analysis Using Full Newton–Raphson Method. Energies, 16(11), 4279. https://doi.org/10.3390/en16114279

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop