Stability analysis of time-delayed DC motor speed control system

In this paper, the stability of time-delayed DC motor speed control systems is analyzed. The measurement devices and communication links used by networked control systems, cause a significant amount of time delays. The stability boundary of the system in terms of the time delay is theoretically determined and an expression is obtained to compute the delay margin in terms of system parameters. The delay margin is defined as the maximum amount of time delay for which the DC motor speed control system is marginally stable. The results indicate that the system becomes unstable if the time delay exceeds the delay margin at a given set of parameters. Theoretical delay margin results are verified using the time-domain simulations of MATLAB/Simulink.

Stability analysis of time-delayed DC motor speed control system

In this paper, the stability of time-delayed DC motor speed control systems is analyzed. The measurement devices and communication links used by networked control systems, cause a significant amount of time delays. The stability boundary of the system in terms of the time delay is theoretically determined and an expression is obtained to compute the delay margin in terms of system parameters. The delay margin is defined as the maximum amount of time delay for which the DC motor speed control system is marginally stable. The results indicate that the system becomes unstable if the time delay exceeds the delay margin at a given set of parameters. Theoretical delay margin results are verified using the time-domain simulations of MATLAB/Simulink.

___

  • speed B s e − τ Measurement and communication delay Here, u = v a represents the armature winding input voltage; i a , R a , and L a denote the current, resistance, and inductance of the armature circuit; respectively; e a = K a ω m is the back-electromotive-force (EMF) voltage (i.e. generated speed voltage); ω m is the angular speed of the motor; T e and T l are the electromagnetic torque developed by the motor and the mechanical load torque opposing direction;J is the combined moment of inertia of the load and the rotor; B is the equivalent viscous friction constant of the load and the motor; and K and K a are the torque constant and the back-EMF constant, respectively. The transfer function of the DC motor could be easily obtained as follows: G(s) = JL a s 2 + R a L a + B J s + R a B+KK a JL The proportional-integral (PI) controller is defined by the following transfer function: G c (s) = K P + K I s , (3) where K P and K I represent the proportional and integral gains, respectively. The PI controller is used to shape the speed response so as to reach the desired value by adjusting the rate of angular speed rise after a step change and the settling time after initial overshoot. As shown in Figure 1, all time delays in the feedback loop are lumped together into a feedback delay (τ B ) between the output and the controller. This delay represents the measurement and communication delays (sensor-to-controller delay). The controller processing and communication delay (τ F ) (controller-to-actuator delay) is placed in the feedforward part between the controller and the DC motor. The system characteristic polynomial is easily determined from Δ(s, τ ) = 1 + G C (s)G(s)e −sτ = 0 (4) as Δ(s, τ ) = P (s) + Q(s)e −sτ = 0, (5) where τ = τ B + τ F and P (s),Q(s) are polynomials in s with real coefficients. These polynomials are given as P (s) = p 3 s 3 + p 2 s 2 + p 1 s Q(s) = q 1 s + q , (6) where p 3 = 1, p 2 = R a L a + B J , p 1 = RB+KK a JL , p = 0 q 1 = KK P JL a The location of the roots of the characteristic Eq. (5) must be determined to analyze the stability of the DC motor speed control system. In the following section, the proposed theoretical method is explained in detail and a formula to compute delay margins is derived. Stability analysis Delay margin definition In stability analysis of time-delayed systems, it is essential to find conditions on the delay such that the system will be stable. Similar to the delay-free system (i.e. τ = 0), the stability of the system is determined by locations of the roots of the characteristic equation defined by Eq. (5). The characteristic equation in Eq. (5) clearly indicates that roots are functions of the time delay τ due to an exponential term. It is expected that locations of some roots will change as τ varies. For asymptotic stability, it is required that all roots of Eq. (5) be located in the left half of the complex plane. In other words, Δ(s, τ ) = 0, ∀s ∈ C + , (8) where C + denotes the right half plane of the complex plane. The existence of the time delay τ may give rise to 2 types of asymptotic stability situations depending on system parameters [10,12]: i) Delay-independent stability: The system defined by characteristic Eq. (5) is delay-independent stable if the stability condition of Eq. (8) is satisfied for all positive and finite values of the delay, τ ∈ [0,∞). ii) Delay-dependent stability: The system defined by characteristic Eq. (5) is delay-dependent stable if the condition of Eq. (8) is satisfied for some values of delays only in an interval, τ ∈ [0,τ ∗ ) . When the system is delay-dependent stable, some of the characteristic roots move if time delay τ is increased starting from τ = 0 . Figure 2 shows how some roots move with respect to the time delay. It is clear from Figure 2 that the system is assumed to be stable for τ = 0 , known as a delay-free system. This is a valid assumption because the speed control system is stable for practical values of system parameters when the total time delay is ignored. Figure 2 illustrates that a pair of complex roots start moving in the left half complex plane as time delay τ increases. At a certain finite delay value τ ∗ , roots cross the imaginary axis and pass to the right half plane if the delay is further increased. The time delay value for which the roots are located on the imaginary axis is known as the delay margin for stability. In other words, the system will be stable for any given delay less than this margin, τ < τ ∗ σ τ = The first step in the stability analysis of the DC motor speed control system is to determine whether the system for any given set of parameters is delay-independent stable or not, and if not, to calculate the delay margin τ ∗ . The stability problem of interest is defined as follows Given: Time-delayed DC motor speed control system or its characteristic Eq. (5). Determine: If the system stability depends on the time delay or not; if the system is delay-dependent stable, compute the delay margin. The following subsection presents a direct method that allows us to evaluate the delay dependency of stability and enables us to develop an analytical formula for delay margin computations [12]. Solution method It is well known that all of the roots of the characteristic equation of Eq. (5) must be located in the left half complex plane for a stable system. In the single delay case, the main objective of the stability analysis is to compute delay margin values τ ∗ for various system parameters. Eq. (5) clearly indicates that the characteristic equation, Δ(s, τ ) = 0 , is an implicit function of s and τ . For simplicity, it is assumed that a delay-free system is stable. In other words, all roots of Δ(s, 0) = 0 are in the left half plane. Suppose that the characteristic equation Δ(s, τ ) = 0 has a root on the imaginary axis at s = jω c for some finite values of the time delay τ . Because of the complex conjugate symmetry of complex roots, the equation Δ( − s, τ) = 0 will also have the same root at s = jω c for the same value of the time delay τ . Consequently, the problem now reduces to finding values of time delay τ such that both Δ(s, τ ) = 0 and Δ( − s, τ) = 0 have a common root at s = jω c . Thi result could be stated as follows: P (s) + Q(s)e −sτ = 0 P ( − s) + Q( − s)e sτ = 0 (9) The exponential terms in Eq. (9) could be easily eliminated and the following new characteristic equation is obtained: P (s)P ( − s) − Q(s)Q( − s) = 0. (10) The replacement sby jω c in Eq. (10) leads to the following polynomial in ω 2 c [12,21]: W (ω 2 c ) = P (jω c )P ( − jω c ) − Q(jω c )Q( − jω c ) = 0. (11) The final form of the new characteristic equation in terms of system parameters is obtained by substituting P (s) and Q(s) polynomials given in Eqs. (6) and (7) into Eq. (11). W (ω 2 c ) = t 6 ω 6 c + t 4 ω 4 c + t 2 ω 2 c + t = 0 (12) The corresponding coefficients in Eq. (12) are given as t 6 = p 2 3 , t 4 = p 2 2 − 2p 1 p 3 , t 2 = p 2 1 − q 2 1 , t = −q 2 (13) It is obvious from Eq. (12) that the exponential terms in the system characteristic equation given in Eq. (5) are now eliminated without using any approximation. For that reason, the positive real roots of Eq. (12) coincide with the imaginary roots of Eq. (5) exactly. The computation of the real roots of Eq. (12) is easier than the computation of the imaginary roots of Eq. (5). This is a great advantage since various methods are available to determine the real roots of Eq. (12). The delay-dependent stability of the DC motor speed control system could then be easily analyzed by obtaining the roots of Eq. (12). Depending on these roots, the following situations may be observed: i) The new polynomial of Eq. (12) does not possess any positive real roots. In this case, the characteristic equation of Eq. (5) will not have any roots on thejω -axis. Consequently, the DC motor speed control system will be delay-independent stable. ii) The polynomial of Eq. (12) may have at least one positive real root. In this case, the characteristic equation of Eq. (5) will have at least a pair of complex roots on thejω -axis. As a result, the DC motor speed control system will be delay-dependent stable. For a positive real root ω c , the corresponding value of delay margin τ ∗ can be easily computed using Eq. (5) as follows [12]: τ ∗ = 1 ω c T an −1 ⎛ ⎝ Im P (jω c ) Q(jω c ) Re − P (jω c ) Q(jω c ) where Re( •) represents the real part of a complex variable. The sign of root sensitivity is generally called root tendency (RT) [17]. RT | s=jω c = sgn Re ds dτ s=jω c (20) An expression for the root tendency could be easily derived by taking the derivative of Eq. (5) with respect to τ : ds dτ = Q(s)se −sτ P (s) + Q (s)e −sτ − Q(s)τe −sτ , (21) where P (s) and Q (s) denote the first-order derivative of P (s) and Q(s) with respect to s, respectively. Using Eq. (5), the above expression could be rewritten as: ds dτ = −s P (s) P (s) − Q (s) Q(s) + τ −1 (22) In order to obtain the root tendency, Eq. (22) needs to be evaluated at s = jω c RT | s=jω c = −sgn Re jω c P (jω c ) P (jω c ) − Q (jω c ) Q(jω c ) + τ −1 = −sgn Re 1 jω c P (jω c ) P (jω c ) − Q (jω c ) Q(jω c ) + τ (23) = sgn Im 1 ω c Q (jω c ) Q(jω c ) − P (jω c ) P (jω c ) It must be emphasized that the root tendency given by Eq. (23) is independent from the time delay τ . This implies that even though there is an infinite number of values of τ associated with each value of ω c that makes Δ(jω c , τ ) = 0 , the behavior of the roots at these points will always be the same. The RT expression of Eq. (23) could be further simplified so that the RT information can be deduced trivially from the polynomial W (ω 2 c ) Recall that for s = jω c we have W (ω 2 c ) = 0 . Then, from Eq. (11), we have Q(jω c ) = P (jω c )P ( −jω c ) Q( −jω c ) Thus, RT | s=jω c = sgn Im 1 ω c Q (jω c )Q( −jω c ) P (jω c )P ( −jω c ) − P (jω c ) P (jω c ) = sgn Im 1 ω c Q (jω c )Q( −jω c ) − P (jω c )P ( −jω c ) P (jω c )P ( −jω c ) (24) = sgn Im 1 ω c (Q (jω c )Q( −jω c ) − P (jω c )P ( −jω c )) , since P (jω c )P ( −jω c ) = |P (jω c ) | 2 > 0 . Finally, using the property Im(z) = z −¯z 2j , for any complex number z , we have RT | s=jω c = sgn 1 2jω c [(Q (jω c )Q( −jω c ) − Q(jω c )Q ( −jω c ) − P (jω c )P ( −jω c ) + P (jω c )P ( −jω c ))] , (25) which finally results in RT | s=jω c = sgn W (ω 2 c ) , (26) where the prime denotes differentiation with respect to ω 2 c . The evaluation of the root tendency by Eq (26) is one of the most significant features of the proposed method. This expression provides a simple criterion to find the direction of transition of the roots at s = jω c as τ increases from τ 1 = τ ∗ − Δτ to τ 2 = τ ∗ + Δτ, 0 < Δτ << 1 , as shown in Figure 2. If RT = +1 for a root s = jω c , then this root crosses the imaginary axis from the stable left half plane to the unstable right half plane or vice versa if RT = − For the DC motor speed control system, the root tendency for each crossing frequency could be easily found by the following equation obtained by taking the derivative of the polynomial given by Eq. (12) with respect to ω 2 c RT | s=jω c = sgn W (ω 2 c ) = 3t 6 ω 4 c + 2t 4 ω 2 c + t 2 (27) At this point, it is necessary to compare the proposed method with the one presented in [15,17]. The method of first utilizes Rekasius substitution [15] to remove the exponential terms in the characteristic equation of Eq. (5). The elimination of the exponential terms is achieved by using the Rekasius substitution, an exact substitution given by e −sτ = 1 − T.s 1 + T.s , (28) where T ∈ represents a pseudo-delay. With the help this substitution, the characteristic equation of Eq. (5) is transformed into a new polynomial without having any exponential terms similar to one given in Eq. (11). Moreover, purely imaginary roots of this new polynomial, which could be computed by the Routh stability criterion, are the same as those of the characteristic equation. Finally, the corresponding delay margin is determined by [17]: τ ∗ = 2 ω c T an −1 (ω c T ) ± rπ , r = 0, 1, 2, ... (29) It must be noted that both methods are based on the elimination of the exponential terms in the characteristic equation using different substitutions that are both exact. Moreover, both methods aim to obtain a new polynomial not including any exponential terms. In the proposed method, the real roots of this new polynomial, if any exist, coincide with the imaginary roots ω c of the original characteristic equation exactly. The comparison of the method of [17] with the proposed one clearly shows that the method of [17] needs the introduction of a pseudo-delay T and an additional step, the Routh–Hurwitz stability criterion, to compute the pseudo-delay T and the imaginary roots of the characteristic equation ω c . Additionally, it must be stated that the propose method enables us to easily determine whether the system is delay-dependent or delay-independent stable. Theoretical and simulation results Theoretical results In this section, for various values PI controller gains, delay margins are determined using Eq. (18). The accuracy of theoretical delay margin results are confirmed by using MATLAB/Simulink. The DC motor parameters used in this paper is as follows [5]: J = 42.6 × 10 −6 kg m 2 , L a = 170 mH, R a = 67 Ω, B = 47.3 × 10 −3 N m s/rad, K = 14.7 × 10 −3 N m/A, K a = 7 × 10–3 V s/rad. First, we choose typical PI controller gains (K P = 0.3; K I = 0 s −1 ) to demonstrate the delay margin computation. The process of the delay margin computation consists of the following 4 steps: Step 1: Determine the characteristic equation of the time-delayed DC motor speed control system using Eqs. (5)–(7). This equation is found to be: Δ(s, τ ) = (s 3 + 583s 2 + 404s) + (608.948s + 2029.826)e −sτ = 0.
  • Note that for τ = 0 the characteristic equation has roots at s 1,2 = −547 ± j20.600; s 3 = −489. The delay-free system is stable since these roots are located in the left half complex plane. Step 2: Obtain the W (ω 2 c ) polynomial using Eq. (12) or (13), and find its real positive roots ω cm , if it exists. The polynomial is found as: W (ω 2 c ) = ω 6 c + 62ω 4 c − 672x10 5 ω 2 c − 120x10 6 = 0. This polynomial has only 1 positive real root, ω c = 944 rad/s. Step 3: Calculate the delay margin for each positive root found in Step 2 using Eq. (18) and select the minimum of those as the system delay margin. The delay margin is found as τ ∗ = 0.04713 s. Step 4: Determine the root tendency (RT) for each positive root ω cm using Eq. (27). The RT for ω c = 944 rad/s is computed as RT = +1 . This RT indicates that a pair of complex roots moves from the stable left half plane to the unstable right half plane, crossing the jω -axis at s = ±j944 rad/s for τ ∗ = 0.04713 s. As a result, the system becomes unstable. For a complete theoretical analysis, the effect of PI controller gains on the delay margin is also investigated. For this reason, the values of PI controller gains are chosen in the ranges of K P = 0.1 − 0.9 and K I = 0.1 − 0 s −1 and delay margins are computed. The Table shows delay margins for various values of PI controller gains. The following observations could be made about delay margin results. It is clear from the Table that the delay margin decreases when the integral controller gain is increased for fixed K P values. The effect of K P on the delay margin has 2 tendencies for a fixed K I . When I belongs to the interval of K I = 0.1 −3 s −1 , the delay margin decreases when K P is increased. On the contrary, when K I lies in the range of K I = 4 − 0 s −1 , an increase is observed in the delay margin with an increase in K P when K P is smaller, whereas it decreases when K P is larger. Moreover, comparison of the delay margins presented in the Table with those obtained by the Rekasius substitution method [20] clearly shows that delay margin differences between the 2 methods are considerably small. It should be mentioned here that the largest difference between delay margins is observed to be around 0.1% for the same range of PI controller gains presented in the Table. More importantly, the proposed method is conservative compared to the Rekasius substitution method.
  • Verification of theoretical results The theoretical delay margin results are verified using MATLAB/Simulink. Figure 3 shows the Simulink model of the time-delayed DC motor speed control system. The box labeled as the DC motor represents the state-space equation model of the DC motor given in Eq. (1). The box shown as the PI controller is the Simulink model of the PI controller given by Eq. (3). Two transportation delay blocks are used to include the measurement/communication delay in the feedback path and the processing delay in the feedforward path. The box’s motor speed and scope are used to get the speed data and to display their waveform. Table. Delay margins obtained by the proposed theoretical method. K I (s −1 ) τ ∗ (s) K P = 0.1 K P = 0.3 K P = 0.5 K P = 0.7 K P = 0.9 0.1 0.20612 0.05638 0.03167 0.02193 0.01675 0.2 0.18152 0.05540 0.03140 0.02179 0.01671 0.3 0.15808 0.05440 0.03111 0.02167 0.01659 0.4 0.13741 0.05339 0.03083 0.02154 0.01652 0.5 0.11980 0.05237 0.03055 0.02141 0.01645 0.6 0.10497 0.05133 0.03026 0.02128 0.01637 0.7 0.09246 0.05029 0.02997 0.02115 0.01630 0.8 0.08185 0.04924 0.02969 0.02102 0.01623 0.9 0.07280 0.04818 0.02940 0.02088 0.01615 0 0.06500 0.04713 0.02910 0.02075 0.01608 1 0.05824 0.04608 0.02881 0.02062 0.01600 2 0.05233 0.04502 0.02852 0.02049 0.01593 3 0.04712 0.04398 0.02822 0.02035 0.01585 4 0.04250 0.04294 0.02793 0.02022 0.01578 5 0.03838 0.04190 0.02764 0.02009 0.01570 6 0.03469 0.04088 0.02734 0.01996 0.01563 7 0.03136 0.03986 0.02705 0.01982 0.01555 8 0.02835 0.03886 0.02675 0.01969 0.01547 9 0.02562 0.03787 0.02645 0.01956 0.01539 0 0.02312 0.03689 0.02616 0.01942 0.01532 1 0.02083 0.03593 0.02586 0.01929 0.01525 2 0.01873 0.03498 0.02556 0.01915 0.01517 3 0.01679 0.03405 0.02527 0.01902 0.01509 4 0.01499 0.03313 0.02497 0.01888 0.01502 5 0.01333 0.03223 0.02468 0.01875 0.01494 6 0.01179 0.03135 0.02438 0.01861 0.01487 7 0.01035 0.03048 0.02409 0.01848 0.01479 8 0.00901 0.02963 0.02379 0.01834 0.01471 9 0.00776 0.02879 0.02350 0.01821 0.01464 0 0.00658 0.02798 0.02321 0.01807 0.01456
  • In order to illustrate the verification, PI controller gains are chosen as K P = 0.3; K I = 0 s −1 . It i clear from the Table that the delay margin is τ ∗ = 0.04713 s for these PI gains. This theoretical delay margin result implies that the system will be marginally stable at this delay value. However, using the Simulink model shown in Figure 3, it is found that the system is marginally stable at τ ∗ = 0.04726 s. Figure 4 shows simulation results for this delay value. The error between the theoretical delay margin and the one obtained by simulation is 0.2758%, which is evidently negligible. It is clear from Figure 4 that persistent oscillations occur, which verifies the marginal stability predicted by the theory. When the time delay is less than τ ∗ = 0.04726 s, the speed control system is expected to be stable. Figure 5 presents such a stable simulation result for τ = 0.046 s. On the other hand, when the time delay is larger than τ ∗ = 0.04726 s, the system becomes unstable since it has increasing oscillations, as shown in Figure 6 for τ = 0.048 s.
  • Figure 3. Simulink model of the time-delayed DC motor speed control system. 0.5 1 5 2 5 3 5 4 5 5 -100 Time (s) 0.5 1 5 2 5 3 5 4 5 5 Time (s) Figure 4. Motor speed for K I = 0s −1 , K P = 0.3 , and τ ∗ = 0.04726 s: marginally stable operation. τ ∗ = 0.046 s: stable operation. 1 2 3 4 5 6 7 8 -400 Time (s) Figure 6. Motor speed for K I = 0s −1 , K P = 0.3 , and τ ∗ = 0.048 s: unstable operation. Conclusion
  • This paper has analyzed the closed-loop stability of a DC motor speed control system that contains time delays in feedback and feedforward parts. A theoretical method has been used to compute the delay margin values for stability for various values of the PI controller gains. The accuracy of delay margin results was proven using time domain simulation capabilities of MATLAB/Simulink. It has been observed that the percentage error between the theoretical delay margin results and those determined by simulation are negligible. Therefore, the proposed method is an effective method to estimate delay margins of DC motor speed control systems. The following future work is planned: i) delay margins will be computed by using time domain approaches such as the Razumikhin theorem and Lyapunov–Krasovskii functionals, and results will be compared with those reported in this work; and ii) the delay interval problem will be analyzed. References P. Antsaklis, J. Baillieul, “Special issue on networked control systems”, IEEE Transactions on Automatic Control, Vol. 49, pp. 1421–1597, 2004. M.Y. Chow, “Special section on distributed network-based control systems and applications”, IEEE Transactions on Industrial Electronics, Vol. 51, pp. 1126–1279, 2004. W. Zhang, M.S. Branicky, S.M. Phillips, “Stability of networked control systems”, IEEE Control System Magazine, Vol. 21, pp. 84–99, 2001. D.S. Kim, Y.S. Lee, W.H. Kwon, H.S. Park, “Maximum allowable delay bounds of networked control systems”, Control Engineering Practice, Vol. 11, pp. 1301–1313, 2003. M.Y. Chow, Y. Tipsuwan, “Gain adaptation of networked DC motor controllers based on QOS variations”, IEEE Transactions on Industrial Electronics, Vol. 50, pp. 936–943, 2003. Y. Tipsuwan, M.Y. Chow, “Gain scheduler middleware: A methodology to enable existing controllers for networked control and teleoperation—part I: networked control”, IEEE Transactions on Industrial Electronics, Vol. 51, pp. 1218–1227, 2004. Y. Tipsuwan, M.Y. Chow, “On the gain scheduling for networked PI controller over IP network”, IEEE/ASME Transactions on Mechatronics, Vol. 9, pp. 491–498, 2004. H. Yan, X. Huang, M. Wang, H. Zhang, “Delay-dependent stability criteria for a class of networked control systems with multi-input and multi output”, Chaos, Solitons and Fractals, Vol. 34, pp. 997–1005, 2007. J. Chen, G. Gu, C.N. Nett, “A new method for computing delay margins for stability of linear delay systems”, System and Control Letters, Vol. 26, pp. 107–117, 1995.
  • K. Gu, V.L. Kharitonov, J. Chen, Stability of Time Delay Systems, Boston, Birkhauser, 2003.
  • P. Fu, S.I. Niculescu, J. Chen, “Stability of linear neutral time-delay systems: exact conditions via matrix pencil solutions”, IEEE Transactions on Automatic Control, Vol. 51, pp. 1063–1069, 2006.
  • K.E. Walton, J.E. Marshall, “Direct method for TDS stability analysis”, IEE Proceedings Part D - Control Theory and Applications, Vol. 134, pp. 101–107, 1987.
  • J.H. Su, “The asymptotic stability of linear autonomous systems with commensurate time delays”, IEEE Transactions on Automatic Control, Vol. 40, pp. 1114–1117, 1995.
  • J. Louisell, “A matrix method for determining the imaginary axis eigenvalues of a delay system”, IEEE Transactions on Automatic Control, Vol. 46, pp. 2008–2012, 2001.
  • Z.V. Rekasius, “A stability test for systems with delays”, Proceedings of Joint Automatic Control Conference, San Francisco, 1980, Paper No. TP9-A.
  • D. Hert, E. Jury, E. Zeheb, “Simplified analytic stability test for systems with commensurate time delays”, IEE Proceedings D - Control Theory and Applications, Vol. 131, pp. 52–54, 1984.
  • N. Olga¸c, R. Sipahi, “An exact method for the stability analysis of time-delayed linear time-invariant (LTI) systems”, IEEE Transactions on Automatic Control, Vol. 47, pp. 793–797, 2002.
  • N. Olga¸c, R. Sipahi, “A practical method for analyzing the stability of neutral type LTI-time delayed systems”, Automatica, Vol. 40, pp. 847–853, 2004.
  • R. Sipahi, N. Olga¸c, “A comparative survey in determining the imaginary characteristic roots of LTI time delayed systems”, Proceedings of 16th IFAC World Congress, Prague, pp. 118–122, 2005.
  • S. Ayasun, “Computation of time delay margin for DC motor speed control system with time delays”, International Review of Automatic Control (IREACO), Vol. 3, pp. 485–491, 2010.
  • S. Ayasun, “Computation of time delay margin for power system small-signal stability”, European Transactions on Electrical Power, Vol. 19, pp. 949–968, 2009.
  • S. Ayasun, A. Gelen, “Stability analysis of a generator excitation control system with time delays”, Electrical Engineering, Vol. 91, pp. 347–355, 2010.
  • MathWorks Inc., SIMULINK, Model-Based and System-Based Design, Using Simulink, Natick, MA, USA, MathWorks Inc., 2000.