A Methodology for Determining Permissible Operating Region of Power Systems According to Conditions of Static Stability Limit

Abstract

For power systems with long-distance ultra-high-voltage (UHV) transmission lines, power transmission limits are often determined by static stability limits. Therefore, the assessment of stability and finding solutions to improve the stability reserve are essential for the operation of the system. This article presents an analytical approach to construct limit characteristics according to static stability conditions on a power plane. Based on the approach proposed, a program is developed and tested on a system with long-distance UHV transmission lines, showing a good performance.

1. Introduction

The rapid development of economies, science, and technology in the world has led to an increasing demand for energy consuming in general and electricity in particular. Power systems are growing and expanding to link power systems of regions or to connect large power plants to consumers using ultra-high-voltage (UHV) or extra-high-voltage (EHV) transmission lines. Large interconnected power systems with long UHV or EHV lines have different features to small ones. The presence of a large amount of reactive power generated by a high voltage line during its operation influences the transmission capacity of the line and the stability of the system.

The definition and classification of power system stability are given in detail in [1]. Power system stability is broadly classified into three categories [1]: rotor angle stability [2,3,4], frequency stability [5,6] and voltage stability [7]. Voltage stability is defined as the ability of a system to maintain steady voltages at all buses in the system after being subjected to a disturbance from a given initial operating condition [1,7]. In voltage stability analysis, it is useful to further classify voltage stability into two subcategories, i.e., large disturbance and small disturbance voltage stability. The former refers to the ability of the system to maintain stable voltages followed by large disturbances such as system faults, loss of generation, or circuit contingencies. On the contrary, the latter refers to the capability of the system to continue steady voltages when subjected to small perturbations such as incremental changes in system load. Generally, the analysis for small disturbances is done as a steady-state stability analysis. In the present paper, static voltage stability is considered. In static stability analysis, the system is assumed to be operated in an equilibrium state. Static analysis evaluates the feasibility of the operating point to provide the system operators with a permissible region where the system can operate normally.

It is of great importance to operators of the system to identify how far the power system is from voltage collapse in order to carry out appropriate remedial actions to avoid unexpected results. The distance towards voltage instability can be obtained by using voltage stability indices. A detailed review of the voltage stability indices can be found in [8]. Voltage stability indices could generally be divided into two groups: Jacobian matrix and system variables based voltage stability indices [9]. Jacobian matrix based voltage stability indices could be used to predict the voltage collapse point and determine the voltage stability margin. The main drawback of such indices is that they require high computational time; hence, they are not suitable for online voltage stability studies. On the other hand, system variables based voltage stability indices are obtained by using the elements of the admittance matrix of the system and some system variables, e.g., bus voltages, power flows, etc., so they require less computational time compared to the first group. Thanks to this advantage, the second group is appropriate for online voltage stability assessment.

So far, many techniques for assessing voltage stability have been developed in the literature. Among them, P-V and Q-V curves are still widely used both in the research and in industry [7,10,11,12,13,14,15]; with these curves, the maximum active and reactive power at a system load bus can be determined when the voltage at that bus reaches the voltage collapse boundary. In particular, a P-V curve gives the relationship between the real power load and bus voltage. It can intuitively provide the operators of the system voltage stability criterion as the margin between the voltage collapse point and the current operating point. On the contrary, a Q-V curve shows the change of bus voltages with respect to reactive power injection or absorption. It gives the operators information on the maximum reactive power at a bus when reaching the minimum voltage limit. In general, both procedures for building P-V and Q-V curves are time-consuming since a large number of power flows need to be executed using conventional methods and models [16]. Owing to this drawback, they cannot to be used for online analysis. Moreover, P-V and Q-V curves could be used to analyze voltage stability only for certain increasing modes and not for the whole view of the analysis. Commonly, those procedures focus on individual buses by stressing the bus considered independently; hence, they are not able to fully reflect the real stability condition of the system.

For determining the visual characteristics of the voltage stability limit and the operating region for the system, a P-Q curve [7,17,18,19,20] is very useful. On a power plane, the operating region indicates both the possible active and reactive power in which the load bus of the system can operate normally. Once the boundary is determined, the distance to the voltage collape of the system can be directly assessed. In [17], the curve is assumed to follow a circle. In [7,18,19,20], the P-Q curve is defined by a parabolic equation. Those assumptions about the shape of the curve are not realistic and may lead to the over-estimation of the reactive power limit.

In this paper, we develop an analytical method to construct a P-Q curve point-by-point, without using any assumption about the shape of the curve on a power plane, and provide a visual operating region for voltage stability studies. The technique can be applied to large-scale power systems to assess voltage stability in real-time. In addition, the P-Q curves we build are realistic since there is no assumption about the change of load used.

The remainder of this paper is organized as follows: In Section 2, a new technique for determining the permissible operating regions of power systems according to static stability limit on the power plane is represented. In Section 3, the testing of the proposed methodology on a system with long-distance UHV or EHV transmission lines is described and the results are discussed. Finally, Section 4 concludes the paper.

2. Algorithm for Determining Permissible Operating Region of Power Systems According to Conditions of Static Stability Limit on Power Plane

2.1. Determining Points on Stability Boundary

For assessing the stability of a power system, stability margins are often used, i.e., K_P = [(P_lim − P₀)/P₀] × 100% and K_Q = [(Q_lim − Q₀)/Q₀] × 100%, where K_P and K_Q are the stability reserve ratios, P_lim and Q_lim are the values of limitation, and P₀ and Q₀ are the values at the operating points for active and reactive power, respectively (see Figure 1). The determination of P_lim and Q_lim for buses is based on completely conventional conditions; for example, P_j of bus j is considered to ascend to the limit, while all remaining parameters are constant [13]. Such determination is different from the actual variation of parameters in operation. Suppose that the characteristics of stability limit of bus j is assumed to be determined in the power plane, as in Figure 1, in which if Q_j is constant, the reserve capacity according to P_j is proportional to distance a, or if cos ϕ_i is constant, it is proportional to distance b, while the parameters may vary with distance c in reality. Thus, for the evaluation of the stability reserve, it is needed to determine the reserve according to the shortest distance d from the operating point to the stability boundary. The challenge is to calculate and draw the characteristics curve under static stable limit conditions to determine the operating region on the power plane, taking into account nonlinearities.

In this section, the pragmatic criterion dQ/dV presented in [21] is used to construct the permissible operating region on a power plane under the conditions of a static stability limit. The operating states of a power system can be assessed based on dQ/dV, i.e., if dQ/dV is negative, the system is stable; otherwise, the system is unstable. Based on this criterion, we develop a mathematical formulation to determine points on the boundary (i.e., determine a curve on the power plane) between stable and unstable states of the system, and then we use this boundary to separate the operating regions into two regions, i.e., stable and unstable regions on the power plane, and assess its stability.

Suppose that a power system has n buses in which there are k generation buses (among k generation buses, Bus 1 is the slack bus of the system). The system is represented by an equivalent diagram, and then an adopted Gaussian elimination method [22] is used to simplify the diagram to a simple one that includes k generation buses and the load bus considered, as shown in Figure 2. It is worth noting that, thanks to this technique, the proposed methodology in this paper can be effectively applied for large-scale power systems.

In Figure 2,

• E_i is the electromotive force of generators connected at the ith generation bus (i = 1 ÷ k);
  ${\dot{Z}}_i=R_i$ + j$X_i$ is the equivalent impedance of the branch connecting generation bus i and the considered load at bus t;
• P_id and Q_id are the real and reactive power transferred from generation bus i to the considered load at bus t and before equivalent impedance ${\dot{Z}}_i$, respectively;
• P_ic and Q_ic are the real and reactive power transferred from generation bus i to the considered load at bus t and after equivalent impedance ${\dot{Z}}_i$, respectively;
• P_t, Q_t, and S_t are the real, reactive, and apparent power of the considered load at bus t;
• ${\dot{Y}}_0$ is the equivalent admittance-to-ground at the considered load bus t, ${\dot{Y}}_0=1/{\dot{Z}}_0$, where ${\dot{Z}}_0=R_0$ + j$X_0$ is the corresponding equivalent impedance;
• P₀, Q₀, and S₀ are the real, reactive, and apparent power transferred through ${\dot{Y}}_0$.

The mathematical formulation for determining points on a stability boundary is as follows.

Based on the power balance conditions at the load bus, we have:

$$\{\begin{array}{l}{Q}_t+Q_0=Q_{1c}+{\displaystyle \sum _{i=2}^k{Q}_{ic}}\\ P_t+P_0=P_{1c}+{\displaystyle \sum _{i=2}^k{P}_{ic}}\end{array}$$

(1)

so:

$$\{\begin{array}{l}{Q}_t=(Q_{1c}-Q_0)+{\displaystyle \sum _{i=2}^k{Q}_{ic}}\\ P_{1c}=(P_t+P_0)-{\displaystyle \sum _{i=2}^k{P}_{ic}}\end{array}$$

(2)

where:

$$\begin{array}{l}{P}_0=G\cdot V_t^2;Q_0=B\cdot V_t^2\\ G=\frac{R_0}{R_0^2+X_0^2};B=-\frac{X_0}{R_0^2+X_0^2}\end{array}$$

(3)

The power transferred from generation bus i to the considered load bus can be calculated using electromotive force ${\dot{E}}_i$ and voltage at the load bus ${\dot{V}}_i$. The equations for calculating the active and reactive power transferred to the end of branch i are as follows

$$\{\begin{array}{}\mathrm{(4a)}& P_{ic}=-\frac{V_t^2}{Z_i}\sin {\alpha }_i+\frac{E_i{V}_t}{Z_i}\sin ({\delta }_i+{\alpha }_i)\mathrm{(4b)}& Q_{ic}=-\frac{V_t^2}{Z_i}\cos {\alpha }_i+\frac{E_i{V}_t}{Z_i}\cos ({\delta }_i+{\alpha }_i)\end{array}$$

where, E_i and V_i are magnitudes of ${\dot{E}}_i$ and ${\dot{V}}_i$, respectively; δ_i is the angle difference between ${\dot{E}}_i$ and ${\dot{V}}_i$; and α_i is the phase angle of ${\dot{Z}}_i$.

From Equations (4a) and (4b), Q_ic is deduced as follows:

$$Q_{ic}=\sqrt{\frac{E_i^2{V}_t^2}{Z_i^2}-{\left(P_{ic}+\frac{V_t^2}{Z_i}\sin {\alpha }_i\right)}^2}-\frac{V_t^2}{Z_i}\cos {\alpha }_i$$

(5)

where:

$$Z_i=\sqrt{R_i^2+X_i^2};\hspace{1em}{\alpha }_i=\frac{\pi }{2}-\mathrm{arctg}\frac{X_i}{R_i}$$

From Equations (2), (3), and (5), the reactive power at the load bus is computed as follows:

$$\begin{array}{cc}\hfill Q_t& =\sqrt{\frac{E_1^2{V}_t^2}{Z_1^2}-{\left((P_t+G\cdot V_t^2)-{\displaystyle \sum _2^k{P}_{ic}}+\frac{\sin {\alpha }_1}{Z_1}{V}_t^2\right)}^2}+\hfill \\ & +{\displaystyle \sum _2^k\sqrt{\frac{E_i^2{V}_t^2}{Z_i^2}-{\left(P_{ic}+\frac{\sin {\alpha }_i}{Z_i}{V}_t^2\right)}^2}}-\left(\frac{\cos {\alpha }_1}{Z_1}+B+{\displaystyle \sum _2^k\frac{\cos {\alpha }_i}{Z_i}}\right)V_t^2\hfill \end{array}$$

(6)

From the formula to calculate the active power at two endings of branch i, P_ic, corresponding to the parameters of the operation modes, can be deduced as:

$$P_{id}=\frac{E_i^2}{Z_i}\sin {\alpha }_i+\frac{E_i{V}_t}{Z_i}\sin ({\delta }_i-{\alpha }_i)$$

so:

$$\sin ({\delta }_i-{\alpha }_i)=(P_{id}-\frac{E_i^2}{Z_i}\sin {\alpha }_i)\frac{Z_i}{E_i{V}_t}$$

(7)

We perform some trigonometric transformations as follows:

$$\begin{array}{ll}& \sin ({\delta }_i+{\alpha }_i)=\sin \left[({\delta }_i-{\alpha }_i)+2\cdot {\alpha }_i\right]\\ \Rightarrow & \sin ({\delta }_i+{\alpha }_i)=\sin ({\delta }_i-{\alpha }_i)\cdot \cos 2{\alpha }_i+\sin 2{\alpha }_i\cdot \cos ({\delta }_i-{\alpha }_i)\\ \Rightarrow & \sin ({\delta }_i+{\alpha }_i)=\sin ({\delta }_i-{\alpha }_i)\cdot \cos 2{\alpha }_i+\sin 2{\alpha }_i\cdot \sqrt{1-{\sin }^2({\delta }_i-{\alpha }_i)}\end{array}$$

(8)

From Equations (4a), (7), and (8), we can determine P_ic as:

$$P_{ic}=-\frac{V_t^2}{Z_i}\sin {\alpha }_i+(P_{id}-\frac{E_i^2}{Z_i}\sin {\alpha }_i)\cos 2{\alpha }_i+\sin 2{\alpha }_i\sqrt{\frac{E_i^2{V}_t^2}{Z_i^2}-{\left(P_{id}-\frac{E_i^2}{Z_i}\sin {\alpha }_i\right)}^2}$$

(9)

As can be seen from Equations (4a) and (6), the reactive power at the load bus Q_t is a function of the voltage at load bus V_t:

$$Q_t=f(V_t)$$

If there is a fluctuation of voltage at the load bus, it will lead to an imbalance of power. The amount of imbalance of reactive power can be identified as follows:

$$\Delta Q=Q_{\sum G}-Q_{\sum T}=f(V_t)-Q_t$$

where Q_ΣG is the total reactive power transferred to the load bus from all generators in the system; Q_ΣT = Q_t + Q₀.

Taking the derivative of both sides of the above equation with respect to V_t, we obtain:

$$\frac{\mathrm{d}\Delta Q}{\mathrm{d}{V}_t}=\frac{\mathrm{d}f(V_t)}{\mathrm{d}{V}_t}-\frac{\mathrm{d}{Q}_t}{\mathrm{d}{V}_t}$$

Suppose that Q_t varies very little with the variation of V_t, so:

$$\frac{\mathrm{d}{Q}_t}{\mathrm{d}{V}_t}=0\Rightarrow \frac{\mathrm{d}\Delta Q}{\mathrm{d}{V}_t}=\frac{\mathrm{d}f(V_t)}{\mathrm{d}{V}_t}$$

According to pragmatic criterion dQ/dV in [14], the system reaches its static stability limit when:

$$\frac{\mathrm{d}\Delta Q}{\mathrm{d}{V}_t}=\frac{\mathrm{d}f(V_t)}{\mathrm{d}{V}_t}=0$$

(10)

The process for determining points on a stability boundary is implemented step by step as follows.

Step 1:

Suppose:

$$Q_{i\min }\le Q_{id}\le Q_{i\max },\forall i=\overline{2,k}$$

(11)

in which, Q_imin and Q_imax are the regulating limits for generating the reactive power of source i (limits of the regulator of an excitation system) and Q_id is the reactive power at the first point of branch i. The nonlinear equations in Equation (10) are solved for voltage at the load bus in limited state V_t,lim. Q_t is calculated by putting V_t,lim and a given value P_t into Equation (6). In order to conclude whether the point M (Q_t, P_t) is located on the static stability limit of the system or not, it is needed to verify Hypothesis (11).

The reactive power at the first point of branch i can be calculated by using the following formula:

$$Q_{id}=\frac{E_i^2}{Z_i}\cos {\alpha }_i-\sqrt{{\left(\frac{E_i{V}_{t,\lim }}{Z_i}\right)}^2-{\left(P_{id}-\frac{E_i^2}{Z_i}\sin {\alpha }_i\right)}^2}$$

(12)

If all the generation buses satisfy Hypothesis (11), point M (Q_t, P_t) will be located on the static stability limit of the system. In case there are m generators within their regulating limits and (k-m) ones exceeding their regulating limits,

$$\begin{array}{l}{Q}_{i\min }\le Q_{id}\le Q_{i\max }\hspace{1em}\hspace{1em}\hspace{1em}\forall i=\overline{2,m}\\ Q_{jd}\le Q_{j\min }\mathrm{or}{Q}_{jd}\ge Q_{j\max }\hspace{1em}\forall j=\overline{m+1,k}\end{array}$$

Step 2:

For source j exceeding its regulating limit:

• Q_jd ≤ Q_jmin set Q_jd = Q_jmin
• Q_jd ≥ Q_jmax set Q_jd = Q_jmax

In this case, the reactive power transferred to load bus can be approximately determined:

$$\{\begin{array}{l}{Q}_{jc}=Q_{jd}-\frac{P_{jd}^2+Q_{jd}^2}{U_{j,rated}^2}{X}_j=\mathrm{const}\\ P_{jc}=P_{jd}-\frac{P_{jd}^2+Q_{jd}^2}{U_{j,rated}^2}{R}_j=\mathrm{const}\end{array}$$

(13)

where U_j,rated denotes the rated voltage at bus j.

The formula to calculate the reactive power at load bus (6) can be changed as follows:

$$\begin{array}{cc}\hfill f(V_t)& =\sqrt{\frac{E_1^2{V}_t^2}{Z_1^2}-{\left(P_t+G{V}_t^2+\frac{\sin {\alpha }_1}{Z_1}{V}_t^2-{\displaystyle \sum _2^m{P}_{ic}}-{\displaystyle \sum _{m+1}^k{P}_{jc}}\right)}^2}\hfill \\ & +{\displaystyle \sum _2^m\sqrt{\frac{E_i^2{V}_t^2}{Z_i^2}-{\left(P_{ic}+\frac{\sin {\alpha }_i}{Z_i}{V}_t^2\right)}^2}}+{\displaystyle \sum _{k+1}^n{Q}_{jc}}\hfill \\ & -\left(B+\frac{\cos {\alpha }_1}{Z_1}+{\displaystyle \sum _2^k\frac{\cos {\alpha }_i}{Z_i}}\right)U_t^2\hfill \end{array}$$

(14)

where, P_ic is calculated in Equation (9) and P_jc and Q_jc are obtained from Equation (13).

Solving nonlinear Equation (10) we can obtain V_t,lim and f(V_t) (see (14)); continue to check m adjustable sources if there is any source violated. If there is no source violated, M(Q_t, P_t) is located on the stability limit; otherwise, we recalculate Step 2.

2.2. Development of Algorithm for Drawing Characteristics of Stability Limit

From the calculation for determining points located on the stability boundary in Section 2.1, we can build an algorithm for drawing the characteristics of the stability limit, as in Figure 3.

3. Application to Building a Program for Determining Permissible Operating Region of Power Systems on Power Plane

Based on the algorithm developed in Section 2, we build a procedure and then a computer program to determine the operating region of power systems using Delphi programming language. For the illustration, we run it for six-bus power system, as shown in Figure 4, in which Power System 2 (PS2) is connected to Power System 1 (PS1) with an unlimited capacity via a long-distance 500 kV UHV transmission line, in the middle of the line at Bus 6, to which a power plant (G2) is connected. In addition, there are two regional electricity networks with low capacity that are represented as two loads connected to Bus 5 and Bus 6. There is also a load connected to Bus 3.

PS1 has an unlimited capacity so the voltage at Bus 1 (U1) can be considered constant, and Bus 1 is considered the slack bus in the system. PS2 is represented by an electromotive force E_eq and an equivalent impedence Z_eq. Therefore, the system considered can be represented by an equivalent diagram, as shown in Figure 5. From this equivalent diagram, the algorithm in Section 2 is applied to construct the characteristics of the stability limit on the power plane. The program is run for the test system with its initial interface, as in Figure 6.

For example, when clicking on Bus 3, the resulting P-Q curve is obtained, as shown in Figure 7 in which the circular point represents the operating point. In this case, the operating point belongs to the stable region (the region covered by two axes and the P-Q curve), indicating that the system is stable. If the load at Bus 3 is increased, exceeding P-Q curve (stability boundary), a warning message will appear to notify the user that the system is unstable, as in Figure 8.

Using the program, the permissible operating region at Buses 3, 5, and 6 can be investigated. Moreover, factors that affect the operating region at a load bus can be examined such as the power change at other load buses in the system, the change in power generated (e.g., P_2f), the number of working units (change in Q_2max) of the power plant connected to Bus 6, the change in power generated (P_3f), and the limit of adjusting the reactive power Q_3max of PS2. In order to identify the critical bus and factors affecting the stability of the system, the operating parameters are assumed to be set in the program, as in

Table 1. Operating parameters for the six-bus test system.

Bus 3		Bus 5		Bus 6		Q2max (Mvar)	P2f (MW)	Q3max (Mvar)	P2f (MW)
P (MW)	Q (Mvar)	P (MW)	Q (Mvar)	P (MW)	Q (Mvar)
600	360	200	120	400	240	700	360	700	200

.

The operating regions at Buses 5, 6, and 3 are plotted in Figure 9, Figure 10 and Figure 11, respectively. By comparing the distances from the operating points to the P-Q curves for all load buses in the system, the most dangerous bus can be identified. In this system, Bus 3 is found to be the most dangerous one that needs to be carefully supervised during the operation of the system.

In particular, while surveying the operating region at Bus 3, as in Figure 11, we decrease Q_3max from 1150 Mvar to 1000 Mvar; the characteristics of stability limit will be changed, as shown in Figure 12. In other cases, we adjust P_2f from 360 MW to 240 MW, for which the characteristics of stability limit are shown in Figure 13, or we adjust the load at Bus 6 from P₆ + jQ₆ = 400 + j240 MVA to P₆ + jQ₆ = 500 + j360 MVA, so that the curve will be obtained as in Figure 14. Thus, by adjusting the operating parameters and observing the change in the P-Q curve, we can find effective solutions to improve the stability of the system.

Extensive testing on the above six-bus system indicates the good performance of the proposed approach in comparison with the results obtained by other commercial software that includes functions for power system stability studies. For example, Figure 15 shows the operating region at Bus 3 and the determination of values for calculating the stability reserve ratios, while

Table 2. Comparison of the results.

Paramaters	Comparison of Results
	Proposed Approach	PSS/E
P0 (MW)	600	600
Q0 (Mvar)	360	360
S0 (MVA)	699.7	699.7
Plim (MW)	800	833
Qlim (Mvar)	800	821
Pi (MW)	772	791
Qi (Mvar)	448	475
Si (MVA)	892.6	923
KP (%)	27.5	31.8
KQ (%)	33.4	38.8
KS (%)	122.3	128

gives the results obtained by the proposed method compared with the results from the Power System Simulator for Engineering (PSS/E) software. It should be noted that the curve we build and the stability reserve ratios that are computed are realistic because we construct the curve point-by-point without using any assumption about the shape of the curve or assumption about the change of load. As illustrated in Figure 15, when surveying the operating region at a certain load bus (e.g., Bus 3), from the operating point of the load and the characteristics curve under the static stable limit built, the stability reserve ratios can be computed as follows. At the operating point, P₀ and Q₀ are known; then S₀ is calculated as $S_0=\sqrt{P_0^2+Q_0^2}$. To determine the limitation of Q (Q_lim in Figure 15), P is kept constant (P = P₀), while Q is increased until it reaches the curve. Analogously, P_lim is determined by increasing P until it reaches, the curve whilst Q is held constant (Q = Q₀). Maintaining constant cosϕ and increasing S (i.e., increasing both P and Q along the line that crosses the origin point and the operating point in Figure 15) until it reaches the curve, P_i and Q_i are found out and S_i is calculated as $S_i=\sqrt{P_i^2+Q_i^2}$. Eventually, stability reserve ratios K_P, K_Q, and K_S are computed:

K_P = [(P_lim − P₀)/P₀] × 100%; K_Q = [(Q_lim − Q₀)/Q₀] × 100%; K_S = [(S_i − S₀)/S₀] × 100%

Using the above formulas, results are obtained as shown in Table 2.

The proposed method is implemented on an Intel Core i5 CPU 2.3 GHz/4.00 GB RAM PC (Computer Dell Vostro 14-5459, Core i5 6200U, Manufacturer: Dell Inc., Chengdu, China), taking 3 × 10⁻³ s, 5 × 10⁻³ s, and 0.11 s for the above six-bus system, IEEE (The Institute of Electrical and Electronics Engineers) nine-bus system, an IEEE 118-bus system [23], respectively. Hence, the proposed approach can be run and can quickly provide results, showing that it can be used for real-time analysis, thanks to the use of an analytical approach in developing the algorithm in this paper.

4. Conclusions

In this paper, we present a methodology for constructing the characteristics of the stability limit of the system according to the static stability conditions on a power plane. Based on the proposed algorithm, a procedure and a program are built and real-time analysis can be performed; the program can be run and quickly provide results (e.g., it needs only 0.11 s for an IEEE 118-bus system with an Intel Core i5 CPU 2.3 GHz/4.00 GB RAM PC) thanks to our use of an analytical approach to develop the algorithm.

With the above-mentioned attractive feature, the proposed approach and the program developed based on it can be applied to a practical power system to monitor its stability reserve in a real-time environment by using online data acquisition systems to provide the necessary information about the system and its operating parameters.

When surveying the operating region of the load buses in the system, by observing the distances from the operating points to the stability boundaries, the operator can easily identify the most critical buses in the system, which should then be carefully monitored during the operation. For example, for the six-bus power system considered in Section 3, Bus 3 is found to the most critical one (Figure 11 for Bus 3 compared to Figure 9 for Bus 5 and Figure 10 for Bus 6).

In addition, by observing the change in the characteristics of the stability limit (P-Q curve) when adjusting operating parameters such as the consumption capacity at the load bus, the power output from power plants, the number of operating units in the system, etc., the most efficient solution for raising the stability reserve of the system can be indicated.

One more important point that should be noted about the proposed methodology is that the curve we build, as presented in Section 2, is realistic since we construct it point-by-point without using any assumption about the shape of the curve or any assumption about the change of load either.