Impact of participation ratios on the stability delay margins computed by direct method for multiple-area load frequency control systems with demand response

ABSTRACT This article studies the effect of dynamic demand response (DR) control on stability delay margins of load frequency control (LFC) systems including communication time delays. DR control is a significant tool to control the responsive loads and increase the reliability of LFC system. The DR control effort on the frequency regulation is provided to each control area of LFC system, called as LFC-DR system. Although the DR control provides some benefits to power grid, communication networks equipped in LFC systems cause communication time delays that degrade dynamic stability of the LFC systems resulting in exponential terms in the characteristic equation of LFC-DR system. This study utilizes an exact method to eliminate the exponential terms without any approximation and transform it into a regular polynomial. The method is utilized to identify stability delay margins for various proportional–integral gains and participation ratios of the secondary and DR control loops for the LFC-DR system. The delay margin values obtained are confirmed by time-domain simulations and a root finder algorithm based on quasi-polynomial mapping. Results indicate that the DR control significantly increases stability delay margins and improves the frequency response of the system as compared with conventional frequency regulation methods.


Introduction
This article investigates the effect of demand response (DR) control on stability margins of load frequency control (LFC) systems including communication time delays. The objective of LFC systems is to regulate the frequency around the nominal value and to maintain scheduled power exchanges of the tie-line connecting control areas when any mismatches between the generation and load demand occurs. In conventional LFC systems, the frequency regulation is generally achieved by adjusting power outputs of conventional thermal or hydro power plants [1]. It is expected that renewable energy (RE) sources including wind power and photovoltaic (PV) systems will have a significant share of power generation in the smart power grid prospect [2,3]. Because of this penetration, the frequency regulation is becoming a challenging task as conventional LFC systems get more complex in terms of frequency regulation. In addition, highly variable generation of RE sources is inadequate to regulate the system frequency.
To expand the integration of RE sources consider- ing the technical and economic feasibilities into power systems, a mapping approach based on the matching between the demand and the supply is reported [4]. This study determines the optimal RE sources capacities with/without energy storage system (ESS) and the location of RE sources using the multi-objective optimization technique on the matching between the demand and the supply. Furthermore, adaptive model predictive control-based linear time-varying Kalman filter [5] and a novel optimized controller-based multiverse optimization algorithm [6] are developed to increase the robustness and stability margin of LFC control integrated with RE sources.
Energy storage devices such as electric vehicles (EVs) [7,8] and responsive loads for dynamic demand control [9][10][11] are becoming promising tools for the frequency control and power grids stability because of the shortcomings of RE sources including high costs, low efficiency, and intermittent nature of their power generations. Owing to the high cost of ESSs, real-time smart active participation of controllable loads known as DR has become an essential tool to balance between peak load demand and generation. Vehicle to grid services and thermostatically controlled loads such as ventilation, heaters, air-conditioners are a few examples of such controllable loads. To utilize a faster and more reliable method for maintaining a balance between generation and demand sides, this concept was first introduced in 1980 by [12]. The change in regular consumption patterns of electricity usage by demand-side resources is referred as DR. The DR programs can be divided into two different types: time-based programs (non-dispatchable) and incentive-based programs (dispatchable) depending upon the factors like economy and network resiliency. DR offers diverse services in this era of smart grids. It can help in incentivizing the customers financially, as well as to benefit the utility companies [13]. The impacts of intermittency of RE sources can also be neutralize using DR [14,15]. To mitigate the frequency and/or voltage fluctuations, ancillary services can be provided using DR [16][17][18]. It can be employed for various other purposes such as planning of transmission expansion [19] and improvement in utilization of transformer [20].
Owing to its fast response, flexibility and economic efficiency, the DR control is a useful compensation for conventional frequency regulation techniques in power system. Therefore, there exists several studies devoted to examine the effect of DR on the frequency regulation for the conventional LFC and automatic generation control schemes. In Reference [21], a DR control loop having communication time delay was first introduced to the traditional single-area LFC model. Results presented in that study clearly illustrated that LFC systems enhanced by a DR control loop (LFC-DR) have a relatively better dynamic performance than that of conventional LFC. Systems. In Reference [22], the DR control loop including time delay was implemented into each control area of a two-area thermal LFC system and the compliant control action between LFC and DR loops was shown to be sufficient for ensuring minimum frequency deviation. The tie-line power is taken as an additional input signal of DR control loop, and genetic algorithm was used to determine optimal controller gains for a quick frequency stabilization in different control areas [17]. To provide robustness to the load disturbances, uncertainties in system parameters and multiple delays in the secondary and DR control loops were considered and a robust Proportional-Integral-Derivative (PID)-type controller for a multi-area LFC-DR system in a deregulated multiarea power system was proposed in Reference [23]. In addition, a single-area LFC system was modified by adding both DR and virtual inertia control loops with associated communication time delays to improve frequency dynamics and the impact of various parameters of DR and virtual control loops such as time delays, their power-sharing factors and frequency dead band was comprehensively analyzed in Reference [24]. An intelligent DR scheme was presented in Reference [14] to determine the control area where the disturbances occurred and to apply the DR exactly to that control area. In addition, a fuzzy-PI-based supervisory controller was proposed as a coordinator between secondary frequency control and DR avoiding large frequency undershoots/overshoots instigated by time delays. A thermostatic load control strategy employing ventilation, heating, and air-conditioning units was proposed in Reference [16] for primary as well as secondary frequency regulation. It was shown that a relatively stable frequency reserve could be provided by considering daily demand profile of thermostatic loads. To decrease frequency detection error and communication delay, a hybrid control approach was developed as a combination of distributed and centralized control methods to control the flexible loads in Reference [25]. To increase the frequency robustness in presence of uncertainties and load deviation, active disturbance rejection control for a single-area LFC system enhanced by DR control loop was proposed in Reference [26]. Moreover, an adaptive delay compensator was designed in the same study to reduce the effect of time delays on the frequency stability. Recently, authors in Reference [27] extended their earlier work reported in Reference [24] to develop a mathematical model for stability and sensitivity analysis of the frequency response of the system considering important parameters associated with DR and virtual control loops. Results presented in that work revealed that the stability as well as the performance of the closed-loop system is sensitive to variations in DR, supplementary, and virtual inertia controls. In practice, conventional controllers like PI, PD, PID controllers having their gains fixed are adopted for the LFC-DR model. However, considering dynamic performance and robustness of LFC systems, the intelligent control techniques associated with linear quadratic regulator (LQR) approach were discussed in Reference [28]. To reduce the frequency deviation caused by time delays, LQR controller with linear matrix inequality was proposed as a coordinator between the DR and secondary control loop in a single-area thermal power system integrated with wind power generation. In Reference [29], the combination of the developed model predictive control method and an improved particle swarm optimization method was proposed to optimize a residential energy consumption system assuming time-variable retail pricing. The studies reported in References [30] and [31] proposed the decentralized active DR (DADR) and stochastic DADR system models in the form of a stochastic control algorithm. In the proposed models, load reduction was initiated depending on the frequency fluctuations measured at the point of the DADR installation.
Communication time delays have become a matter of great concern in the dynamic performance of traditional LFC systems as they can decrease damping performance of the control system and even they can cause instability if they go beyond the stability delay marg in References [32][33][34]. These delays can occur in LFC systems if an open communication network is employed to receive the measured data from central controller to power plant and vice versa. With the increasing integration of RE sources, EVs and DR control, such delays have even become much more significant [21][22][23][24][35][36][37][38]. Even though there exist various studies on the delay-dependent stability and the computation of stability margin in the conventional LFC systems, studies focusing on the effect of both time delays along with the utilization of DR control on the frequency regulation are very limited. For example, studies reported in References [21][22][23][24]27] recognized the importance of time delays observed in the DR control loop on the frequency regulation. However, time delays in the secondary control loop were neglected. It is well known in the literature that communication delay in the secondary loop is larger than one in the DR control loop, which in turn significantly affects the frequency stability [32][33][34]. Moreover, those studies have utilized an approximation technique for the exponential type transfer function representing the time delay in DR control loop. Such an approximation does not reflect the true characteristic of time delays and their impact on the stability. Besides that, it also increases the system dimension depending on the order used in the Padé approximation. More importantly, the exact computation of stability delay margins of LFC-DR systems and analytical studies on the impact of DR control loop on stability delay margins were not presented. Similarly, in Reference [23], stability delay margins were obtained by trial-and-error simulation method rather than using an exact method. In References [25,26], various compensation schemes were proposed to decrease the frequency deviation in the presence of time delays in the DR control loop. In those studies, authors did not present any qualitative/quantitative analysis to determine the impact of the DR controls with time delays on the delay-dependent stability of LFC-DR systems.
To overcome shortcomings of existing studies on the time-delayed LFC-DR systems, this article aims to calculate stability delay margins of LFC-DR systems having time delay on the secondary control loop. Various stability delay margin calculation methods are available in the existing literature that can be classified as timedomain indirect methods and frequency-domain direct method. Identification of the complex roots of the quasi-polynomial characteristic equation on the imaginary axis can be easily done by using these frequencydomain approaches. This group of methods includes the following: (i) direct method based on elimination of exponential terms [39], (ii) Rekasius substitution [40,41] and (iii) frequency sweeping test [36]. Amongst these approaches, the direct method proposed by Walton and Marshall [39] was efficiently implemented to determine stability delay margins of time-delayed twoarea LFC systems and microgrid not including DR control loop [32,42] and single-area LFC system with EVs aggregator denoted by LFC-EVs system [38]. Moreover, the Rekasius substitution was applied to delay margin computation of microgrid LFC system [43], two-area LFC system [44], two-area LFC system having DR control loop [45] and two-area LFC-EVs system [37]. The frequency sweeping test was applied to the delay margin computation of a single-area LFC system with EVs aggregator [36]. The existing studies clearly show that the frequency-domain direct methods accurately compute stability margins of electrical power systems with communication time delays. However, the main drawback of these approaches is that they are applicable to constant delay only. Pekař and Gao in Reference [46] presented a thorough literature review on the approaches used for stability margin assessment of continuous-time linear timeinvariant systems having constant time delays. When mentioned, these methods are compared with each other; some observations have been remarked. Both the direct method and Rekasius substitution aim to eliminate the delay-dependent transcendental (exponential) terms in the quasi-characteristic polynomials using different approaches. The direct method employs a recursive procedure without using any approximation and obtains a regular polynomial without exponential terms whose positive real roots exactly match to the complex roots on the jω-axis of the original quasi-polynomial [32,39]. On the other hand, Rekasius substitution is an exact transformation for the roots lying on the imaginary axis. Therefore, during the substitution, these purely complex roots of the characteristic polynomial with delay-dependent exponential terms are preserved and the system characteristic equation is converted to an ordinary single-variable equation not including any exponential terms. Routh stability criterion then could be easily used to compute purely imaginary roots of the single-variable regular polynomial [40,41]. However, with the construction of the Routh table, high order power system characteristic equation exposes the cumbersome symbolic process. From the computational point of view, this increases the mathematical complexity in terms of computational burden and error when Rekasius substitution is implemented for multiarea LFC-DR systems. The comparisons between Rekasius substitution and direct method are explained in details in Section 3. In addition, a detailed comparison of these methods as applied to the two-area LFC system were presented in Reference [44]. Finally, the frequency sweeping test consisting of a combination of the binary iteration algorithm and frequency sweeping also computes exact delay margin results. However, the selection of the frequency range for the sweeping test requires undesired computational effort [36].
The time-domain indirect approaches rely on linear matrix inequalities and Lyapunov stability theory. These methods can only determine the sufficient conditions for the system stability and there exist various studies focusing on the reduction of its conservativeness [33,34]. Numerous inequalities were proposed in recent years such as Wirtinger inequality, Jensen inequality [47] free-matrix-based inequality [48] and Bessel-Legendre inequality [49]. Time-domain indirect methods were used to calculate stability margins in multi-area LFC systems without DR control loops [33,34]. These approaches can be used for both timevarying and constant delay cases. Although there exist tremendous efforts to reduce the conservativeness of this approach, it is well known in the literature that frequency-domain direct methods give more accurate and less conservative stability delay margins than timedomain indirect methods [32][33][34].
Existing literature has investigated the robust frequency performance of LFC-DR system and the effect of time delays of only DR control loop on frequency stability. Moreover, many studies address the delay-dependent stability analysis of LFC-DR system using the trial-and-error simulation method and Padé approximation rather than an exact analytical method. The simulation study reported in Reference [21] indicates that the largest delay value on DR control loop for larger power system having high inertia and slower response is 500 msec. Whereas, larger delays on secondary control loop depending on communication structure for such power systems are observed. Therefore, not considering the delays of DR control loop, this study examines the impact of DR control on stability delay margins of the LFC system via an analytically elegant method, the elimination of exponential terms [32,38,42].
The main contents of this article are as following: First, for the selected power-sharing factors between the generators and DR control loop, delay margins are computed for a wide-range of PI controller gain values to evaluate the impact of the controller gains. Second, delay margins are determined for various powersharing scenarios between the conventional generator and DR control loop to assess how the participation of controllable responsive loads affect stability delay margins. Third, some case studies for the efficiency of the proposed method and relationship between participation ratio and delay margin are presented. The implementation of an exact analytical method to calculate delay margin in two-area LFC-DR system and a comprehensive analysis of the impact of the DR control loop on the stability delay margin and the frequency regulation are the major contributions of this study. The best of authors' knowledge, an investigation on the impact of DR control on the stability delay margins of LFC systems has not been reported in the literature. Finally, a root finder algorithm based on quasipolynomial mapping known as QPmR algorithm [50] along with time-domain simulations [51] are used to validate the correctness of stability margin results. The comparison of stability delay margins of LFC-DR system with those of LFC system not including DR control loop clearly illustrates that delay margins significantly increase as the participation of the DR control loop into the frequency regulation increases. More importantly, simulation studies prove that the inclusion of the DR control loop reduces undesired oscillations on the frequency response and stabilizes the LFC system including time delays.

LFC-DR system model
Two-area LFC system model with a DR control loop in each control area is presented in Figure  1. It can be observed that the conventional twoarea LFC system is illustrated by solid lines while the DR loop in shown by dashed lines. In this figure, f i , X gi , P mi , P gi , P DR,i , and P Li , (i = 1, 2) denotes the deviation in the frequency, position of valve, mechanical power output, power output of generator, DR control loop power output and load disturbance in each control area, respectively. Moreover, M i , D i , R i , T gi , T ci , T ri , F Pi , β i , ACE i , and T 12 (i = 1, 2) represent the inertia constant, load damping constant, speed regulation constant, time constants of governor, reheat and turbine, fraction of total turbine power, frequency bias factor, area control error, and the tie-line coefficient of each control area, respectively. Please note that a proportional-integral (PI) controller G ci (s) = K Pi + K Ii s is used as LFC controller and DR controller, where K Pi and K Ii are the proportional and integral controller gains, respectively.
With the integration of the DR control into the twoarea LFC system, the required control effort called is shared between DR and secondary loops in each control area as follows [21]: where α 0i and α 1i represent participation ratio of the secondary and DR control loops with α 0i + α 1i = 1. Finally, the measurement and data transfer time delays in τ are lumped. Multiple constant or time-varying delays are generally observed in multi-area LFC systems. As explained in Reference [32], in an open communication network, delay can arise during: (1) transmission of area control error (ACE) signals from the control center to the individual generation units and (2) from a telemetry delay when remote terminal units send the telemetry signals to the control center. Assuming that the control center waits to receive the telemetered values, the analysis for each delay case is identical. Therefore, all delays are generally aggregated into a single constant or time-varying delay from the control center. In the study, time delays in DR control loop are neglected since they are much smaller than those observed in the secondary control loop [22]. It is assumed that the delays in the secondary control loops of both control areas are considered to be equal, and they are constant delays denoted by an exponential function of e −sτ 1 = e −sτ 2 = e −sτ as depicted in Figure 1. With these assumptions, the stability delay margins of LFC-DR system are obtained for various PI controller parameters of DR control and secondary control loops and the participation factors.
The characteristic polynomial of the two-area system in Figure 1 needs to be determined to examine the delay-dependent stability of LFC-DR system and particularly to evaluate the DR control loop effect on delay margins. This polynomial is transcendental type and is given as where P(s), Q(s) and R(s) polynomials of s whose coefficients depend on system parameter and are determined as P(s) = p 13 s 13 + p 12 s 12 + p 11 s 11 + p 10 s 10 + p 9 s 9 + p 8 s 8 + p 7 s 7 + p 6 s 6 + p 5 s 5 + p 4 s 4 + p 3 s 3 Q(s) = q 10 s 10 + q 9 s 9 + q 8 s 8 + q 7 s 7 + q 6 s 6 + q 5 s 5 + q 4 s 4 + q 3 s 3 + q 2 s 2 R(s) = r 7 s 7 + r 6 s 6 + r 5 s 5 + r 4 s 4 + r 3 s 3 + r 2 s 2 (3) The coefficients of polynomials in Equation (3) are not represented here due to insufficient space.

Stability margin identification: direct method
The main purpose of stability analysis of time-delayed dynamical system is to determine that if the system is delay-dependent stable or delay-independent stable. If the system is delay-independent stable, this implies that the system will remain stable for all finite delays. For a system to be delay-dependent stable, the system will be stable for τ < τ * and unstable for τ > τ * where τ * represents the stability delay margin for selected system parameters.
For the LFC-DR system whose characteristic equation is given in Equation (2) to be stable, all roots of this equation have to be located in the stable left half of the s-plane. Observe that the characteristic polynomial of the LFC-DR system in Equation (2) is a quasi-polynomial because of exponential terms e −τ s and e −2τ s . Therefore, the polynomial has infinite many roots and calculating the roots is a challenging task. However, for the delay-dependent stability analysis, all roots are not required and it is essential to find time delay for which the quasi-polynomial has roots located on jω-axis. If the original characteristic polynomial of Equation (2) has a solution of s = jω c , then ( − s, τ ) = 0 being complex conjugate will have the same solution: (−s, τ ) = P(−s) + Q(−s)e τ s + R(−s)e 2τ s = 0 (4) The exponential terms in Equations (2) and (4) have to be removed so that a new characteristic polynomial having no transcendentality can be found. The new characteristic equation can be defined using an iterative procedure [32,39] given below: (1) From Equations (5) and (6), it can be observed that the root s = jω c of Equations (2) and (4) is also a root of the following new characteristic polynomial: (1) It can be noticed that the new characteristic equations in Equation (7) have only a single e −τ s or e τ s . This shows that the previous commensuracy degree of 2 is now reduced to 1. Moreover, after the elimination of e −2τ s term in Equation (2), the degrees of P (1) (s) and Q (1) (s) polynomials have become 26 and 23, respectively. The procedure can be simply repeated to remove exponential terms, e −τ s and e τ s in Equation (7) and an augmented characteristic equation having no transcendental terms can be obtained as follows: (2) (s, τ ) = P (2) (s) = 0 ( 9 ) where It is to be observed that for some τ , the imaginary root of Equation (2) is also a root of Equation (9) because the imaginary roots of the original characteristic polynomial of Equation (2) are well-preserved by the elimination procedure. Following polynomial in ω c 2 is obtained by substituting s = jω c into Equations (9) and (10): This new characteristic polynomial of Equation (11) has a degree of 52. The communication delay value corresponding to the stability margin for which the roots of Equation (2) cross the jω-axis is calculated by [32,39]: Moreover, an expression is required to find the root tendency (RT) of the roots to assess the direction of roots transition at s = jω c while observing an increment ( τ ) in τ values. The root s = jω c crosses the jωaxis either to the stable left half plane for RT = −1 or to the unstable right half plane for RT = +1. It is shown in Reference [32] that the RT of a root of (s, τ ) = 0 is similar to that of the corresponding root of (2) (s, τ ) = 0 when the following condition is satisfied for s = jω c crossing root: The RT is then determined by: where The minimum of the communication delays τ * min found by Equation (12) with RT = +1 is the stability margin of the LFC-DR system as per the definition of stability delay margin. It will be useful to compare the direct method based on elimination of the exponential terms with other frequency-domain Rekasius substitution method [40] applied to the delay margin computation of LFC-DR system with constant communication delays [45]. Equations (5)-(11) clearly elaborate the elimination procedure of the exponential terms presented in the characteristic equation (4) by employing direct method, whereas Rekasius substitution method uses an exact transformation described as e −sτ = (1 − Ts)/(1 + Ts) (T ∈ is called pseudo-delay) to eliminate the exponential terms. Using this substitution, the transcendental characteristic equation (4) is converted into a polynomial without transcendentality similar to one given in Equation (11) such that its purely imaginary roots determined by Routh stability criterion coincide with the purely imaginary roots of the characteristic equation exactly. The corresponding delay margin is determined by [40]: 2 ω c (Tan −1 (ω c T) ± rπ), r = 0, 1, 2, . . . (16) It should be noticed that both methods use different substitutions that are exact and ultimately obtain an augmented characteristic polynomial without transcendentality. However, in the proposed method, the real roots (if there exist any) of this new polynomial coincide with the imaginary roots ω c of the characteristic equation exactly. Whereas, Rekasius substitution method presented in Reference [40] requires the introduction of a pseudo-delay T and an additional computational burden in form of Routh Hurwitz stability criterion to determine the pseudo-delay T (The imaginary roots of the characteristic equation ω c ). The construction of the Routh's array requires the cumbersome symbolic process for high-order power system characteristic equations. The study reported in Reference [45] constructs a Routh's array with columns of 7 and rows of 16 depending on rational polynomials of T for the two-area LFC-DR system. Thus, significantly increasing the mathematical complexity of the study in terms of computational burden and errors. Moreover, the proposed method in this article is a practical and handy criterion to determine whether the system is delay-dependent or independent stable.

Results
This section gives results for stability margins of the LFC-DR system and verification studies of theoretical stability delay margins by QPmR algorithm. System parameters in each area are as follows [35]: T gi = 0.2, T ci = 0.3, T ri = 12,

Step-by-step analysis
The procedure of the delay margin computation is based on the following steps: Step 1: Two-area LFC-DRs system parameters are selected and state-space equation model of the time-delayed system is found.
Step 2: K P and K I controller gains set is selected.
Step 3: For the given PI controller gains, the characteristic polynomial of the two-area LFC-DRs system is determined using Equation (2).
Step 4: The augmented characteristic polynomial of the two-area LFC-DRs system is obtained using Equation (11) and its real positive roots; {ω c } = {ω c1 , ω c2 , . . . , ω cq } are computed. The resulting polynomial of the two-area LFC-DRs system have a degree of n.2 p = Equation (13)2 2 = 52 (p = 2 is the degree of commensuracy and n = 13 is order of the characteristic polynomial).
Step 5: The corresponding root tendencies for all real positive roots are determined using Equation (14).
Step 6: By employing Equation (12), the corresponding stability delay margins for real positive roots found in Step (4) having positive root tendency RT = +1 are computed.
Step 7: The stability margin of the system is chosen as the smallest value τ * min among those computed in Step (6).
Step 8: The accuracy of theoretically computed stability margins is then validated by QPmR algorithm together with time-domain simulations.

Theoretical results of two-area LFC-DRs system
The stability delay margin values using steps 1-8 for various PI controller gains and different participation ratios of DRs are computed and presented in Tables  1-3, respectively. For all cases, stability delay margins decrease as K I increases for the fixed K P , indicating a less stable LFC-DR system, whereas stability delay margins increase for almost all values of K I with an increase in K P values resulting in a more stable LFC-DR system. It should be mentioned here that stability delay margins are not determined for values of K P and K I when the delay-free LFC-DR system (τ = 0) is unstable. The corresponding positions in are labelled as ( * ) in tables. More importantly, the comparison of delay margins in Tables 2 and 3 with ones in Table 1 clearly reveals the fact that stability delay margins significantly increase when a DR control loop with a corresponding participation ratio is included in LFC system for frequency regulation. For example, one can see from Table 1 that for K P = 0.5, K I = 0.3, stability margin is determined as τ * = 2.3049 s when the participation of the secondary control loop is 60% for area 1 and 70% for area 2 (40% participation of the DR control loop for area 1 and 30% participation of the DR control loop for area 2, α 01 = 0.6, α 11 = 0.4 and α 01 = 0.7, α 11 = 0.3). When the participation of the secondary control loop is reduced to 60% for area 1 and area 2 (40% participation of the DR control loop for area 1 and area 2, α 01 = α 02 = 0.6, α 11 = α 12 = 0.4), the stability delay margin increases from τ * = 2.3049 s to τ * = 2.6176 s, which illustrates 13.57% increase. Furthermore, the control effort of the secondary control loop is reduced to 60% for area 1 and 50% for area 2 (40% participation of the DR control loop for area 1 and 50% for area 2) that is α 01 = 0.6, α 11 = 0.4 and α 01 = 0.5, α 11 = 0.5. The stability margin for this scenario rises to τ * = 3.0820 s, which corresponds to 33.71% increase with respect to τ * = 2.3049 s. Finally, theoretical delay margin results computed by proposed method are compared with Rekasius substitution method. It should be noticed that for the all set of PI controller values, delay margin results presented in Table 2 exactly match with those given in Table 3 of our previous study of Reference [45], where stability delay margins have been computed using Rekasius substitution method of Reference [40].

Verification of the results by QPmR algorithm and simulation studies
The accuracy of the theoretical delay margin (τ * ) and root crossing (±jω c ) is verified by QPmR algorithm and using time-domain simulations. Note from Table  2 that the stability delay margin for K P = 0.5, K I = 0.3 is computed as τ * = 2.6176 s. Figure 2 shows the position of dominant roots of the characteristic equation and frequency response of LFC-DR system around the stability delay margin that is for τ = 2.5176 s, τ * = 2.6176 s and τ = 2.7176 s. Figure 2(a) shows the dominant roots and the frequency response for τ = 2.5176 s, which is less than the delay margin, τ * = 2.6176 s It is obvious that all roots of characteristic polynomial in Equation (1) are present in the stable left half of the s − plane and oscillations in the frequency are decaying. Therefore, the LFC-DR system is stable. For τ * = 2.6176 s, a pair of complex roots is now on the imaginary axis as depicted in Figure 2(b), and the frequency response of the LFC-DR system has sustained oscillations, indicating the marginal stability. It should be observed that QPmR algorithm gives the same purely imaginary roots of s = ±jω c = ±0.3811 rad/s as obtained by the proposed method. Finally, a pair of complex roots at τ = 2.7176 s crosses the imaginary axis towards the unstable right half plane, which causes unstable growing oscillations in the frequency response as shown in Figure 2(c).

DR control loop effect on the frequency response
The effect of the participation of the DR control loop on the frequency performance of the LFC-DR system is also investigated for three different cases of time delays when PI controller gains are chosen as K P = 0.5, K I = 0.3 and the participation of the DR control loop for area 1 is selected as α 01 = 0.6, α 11 = 0.4. In the first case, the time delay is chosen as τ = 2.0 s and at first, the participation ratios of the secondary control loop and the DR control loop for area 2 is selected as α 02 = 0.7, α 12 = 0.3. Please note from Table 1 that this case is stable case since the selected delay is less than the stability delay margin, (τ = 2.0 s < τ * = 2.3049 s). The frequency response is shown in Figure 3. It can be seen that the frequency response has undesired oscillations even though it is stable. Next, the participation ratios of the secondary control loop and the DR control loop is selected as α 02 = 0.6, α 12 = 0.4 and α 02 = 0.5, α 12 = 0.5. Figure 3 also compares the frequency response with that of LFC system for α 02 = 0.7, α 12 = 0.3. It is evident from Figure 3 that undesired oscillations quickly damp out when 50% or 70% of the control effort is provided by the DR control loop. This case clearly illustrates that the inclusion of DR control loop significantly improves the frequency response performance of the LFC system. In the second case, time delay is selected as τ = 2.3049 s for which the participation ratios of the secondary control loop and the DR control loop for area 2 is selected as α 02 = 0.7, α 12 = 0.3, and the system is marginally stable as presented in Table 1. Figure 4 compares the frequency responses for the same participation ratios of the first case. It is clear that with the inclusion of the DR control loop, the marginally stable LFC system becomes stable. In the third case, time delay is chosen as τ = 2.4 s for which the participation ratios of the secondary control loop and the DR control loop for area 2 is selected as α 02 = 0.7, α 12 = 0.3, and the system is unstable since this delay is larger than the stability delay margin of τ * = 2.3049 s. Figure 5 clearly shows that the inclusion of DR control loop with a participation ratio of 50% or 70% makes the unstable LFC system stable. The second and third cases clearly illustrates the stabilizing effect of DR control loop for LFC system in the presence of communication time delays.
Finally, the impact of the participation ratio of the DR control on the stability delay margin is analyzed for three different PI controller gains, namely (K P = 0.5, K I = 0.3), (K P = 0.5, K I = 0.5) and (K P = 0.5, K I = 0.7). Figure 6 shows that the stability delay margin remarkably increases for all controller gain values when participation ratio of the DR control loop is

Conclusions
This article has comprehensively investigated the effect of DR control loop on the stability margin of LFC systems having time delays. A frequency-domain direct method has been implemented to compute stability margins in the two-area LFC-DR system for a widerange of PI controller gains and different participation ratios of the secondary and DR control loops. Moreover, theoretical delay margins have been confirmed by the QPmR algorithm and time-domain simulations as well. Following comments could be made from the results: (i) The presented results reveal that the selection of PI controller gains plays an important role in computing communication delay margins and thereby ensuring the stability of the two-area LFC-DR system. One can select the appropriate controller gains using Tables 1-3 for the allowable communication delay value of the two-area LFC-DR systems. (ii) A small contribution of DR share on the frequency control effort significantly increases the stability delay margins for all PI controller gains and thus the DR control loop expands the stability margin of LFC system. (iii) The higher DR shares significantly increase the stabilizing effect on LFC systems. Simulation studies indicate that the increase in DR participation factor reduces the undesired oscillations on the frequency responses and quickly minimizes the frequency deviations of LFC system with time delays. (iv) The comparisons between the direct method and Rekasius substitution method indicate that the direct method for commensurate or single-delay cases does not require a symbolic computational burden and the proposed method could be easily applied to high order power systems experiencing communication delays.
Considering the uncertainties in communication time delays and multi-area LFC-DR system parameters, the impact of DR control loop on controller parameter design and robust controller performance of the LFC system should be investigated in the future studies. The following future works for the LFC systems including DR control have been put in perspective: (i) Identification of stability regions based on stability boundary locus method in the controller parameters space and (ii) Computation of robust stability regions in the controller parameter space using the complex Kharitonov's theorem.