Time-Delay Estimation using the Characteristic Roots of Delay Differential Equations

Problem statement: For ordinary dynamic systems (i.e., non-delayed), various methods such as linear least-squares, gradient-weighted lea st-squares, Kalman filtering and other robust techniques have been widely used in signal processi ng, robotics, civil engineering. On the other hand, time-delay estimation of systems with unknown timedelay is still a challenging problem due to difficulty in formulation caused. Approach: The presented method makes use of the Lambert W function and analytical solutions of scalar first-o rder Delay Differential Equations (DDEs). The Lambert W function has been known to be useful in s olving delay differential equations. From the solutions in terms of the Lambert W function, the d ominant characteristic roots can be obtained and used to estimate time-delays. The function is alrea dy embedded in various software packages (e.g., MATLAB) and thus, the presented method can be readi ly used for time-delay systems. Results: The presented method and the provided examples show eas e of formulation and accuracy of time-delay estimation. Conclusion: Estimation of time-delays can be conducted in an a nalytical way. The presented method will be extended to general system s of DDEs and application to physical systems.


INTRODUCTION
Parameter estimation deals with the problem of obtaining mathematical models that represent dynamic systems based on observed data. For ordinary dynamic systems (i.e., non-delayed), various methods are available in the literature, such as linear least-squares gradient-weighted least-squares, Kalman filtering and other robust techniques (Ljung, 1999;Panich, 2010). They have been widely used in signal processing, robotics, civil engineering. On the other hand, timedelay estimation of systems with unknown time delay is a challenging problem that has attracted continuous attention during the last four decades (Ren, 2005). Even though considerable efforts have been made on parameter estimation, there are still many open problems in time-delay identification due to difficulty in formulation (Richard, 2003;Belkoura et al., 2009;Khan, 2011;El-Fallah and El-Sallam, 2011;Khan and Khan, 2010). Obviously, one can expect that the more accurate the knowledge on the delay is, the higher the achievable control performances will be (Richard, 2003). If the time-delay is known, then, many powerful control techniques (e.g., Smith predictor, Finite Spectrum Assignments (FSA)) can be applied.
The method presented here is inspired from a wellknown technique for Ordinary Differential Equations (ODEs), which is available in various system dynamics textbooks (Palm, 2009). The technique is based on analytical solutions to ODEs in terms of exponential functions. The method has not been feasible to use for systems of DDEs due to lack of analytical solutions to DDEs. The Lambert W function has been known to be useful to solve DDEs (Yi et al., 2010a;Corless et al., 1996). Through the Lambert W function-based method, it is possible to extend the estimation technique for ODEs to DDEs. This study shows the estimation of time-delays based on the Lambert W function with illustrative examples.
Unlike ODEs, two initial conditions need to be specified for DDEs: A reshape function, g (t), for −h≤ t <0 and initial point, x 0 , at time, t = 0. The quantity, h, denotes the time-delay. The characteristic equation of Eq. 1 is given by: The exponential term, e −sh , induced by the timedelay term x (t-h), makes the characteristic equation transcendental (i.e., infinite dimensional and nonlinear). Thus, it is not feasible to find roots of Eq. 2, which has an infinite number of roots. The Lambert W function is defined as (Corless et al., 1996): For detailed explanations on definition and properties of the function, refer to (Corless et al., 1996). Based on the Lambert W function, a new approach has been developed for systematic analysis and control of systems with time-delays (Yi et al., 2010a). The function allows handling of the exponential term in the characteristic equation thanks to its unique definition in Eq. 3. Then, using the Lambert W function, the characteristic equation in Eq. 2 is solved as: As seen in Eq. 4, the characteristic root, s, is obtained analytically in terms of parameters, a, a d and h. Then, the solution to Eq. 2 is given by using the Lambert W function. Then, considering the initial conditions, the solution to Eq. 1 is given by using the Lambert W function as (Asl and Ulsoy, 2003): The coefficient, I k C , is determined analytically (Yi et al., 2010a) or numerically (Asl and Ulsoy, 2003) from initial conditions, g(t) and x 0 . Note that, unlike results by other existing methods for DDEs, the solution in Eq. 5 has an analytical form expressed in terms of the parameters of the DDE in Eq. 1, i.e., a, a d and h. One can explicitly determine how the parameters are involved in the solution and, furthermore, how each parameter affects each characteristic root. That enables one to formulate estimation of time-delays in an analytic way. Also, each eigenvalue is distinguished by k, which indicates the branch of the Lambert W function as seen in Eq. 5. The function is already embedded in various software packages, such as Matlab.
Rightmost characteristic root and the lambert w function: In estimating time-delays, the Lambert W function is used to find the rightmost eigenvalues. For 1st order scalar DDEs, it has been proved that the rightmost characteristic roots are always obtained by using the principal branch, k = 0 (Shinozaki and Mori, 2006). For first-order scalar delay differential equations as in Eq. 1, one has to consider two possible cases for rightmost characteristic roots: characteristic equations of DDEs as in Eq. 2 can have one real dominant root or two complex conjugate dominant roots. Thus, when estimating time-delays using characteristic roots, it is required to decide whether it is the former or the latter. The condition for decision is derived using the Lambert W function (Yi et al., 2010b). The branches of the Lambert W function and the specific range of each branch are helpful in deriving the condition (Corless et al., 1996). As seen in Fig. 1, when the argument of the function, x, is greater than -1/e the values of W 0 (x) is real. However, if x < -1/e, the values of W 0 (x) is complex and, thus, the rightmost characteristic roots are two and complex conjugates. Using the fact, we get the conclusion regarding existence of the complex conjugates rightmost characteristic roots as Eq.
For ODEs, an estimation technique using the logarithmic decrement provides an effective way to estimate the damping ratio, ζ (Palm, 2009). The technique makes use of the form: 2 n n s j 1 = −ξω ± ω − ξ For the characteristic root, s, of 2nd order ODEs. Using the characteristic root in Eq. 7, the well-known exact solution is obtained as: From measured responses and the solution form in Eq. 8, the parameters ζ and ω n are estimated (Palm, 2009). The concept is applied to DDEs as below. Example I: Logarithmic decrement for estimation (Oscillatory): Consider the system in Eq. 1 with parameters, a = −1, a d = −1.5 and h = 1, with x 0 = 1 and g (t) = 1. As explained in Introduction, DDEs have an infinite number of characteristic roots (e.g., see Fig. 2). Thus, the solution to DDEs has the infinite series of the exponential of the characteristic roots as seen in Eq. 5. Note that the rightmost roots are always obtained by using the principal branch (k = 0) and/or k = −1 (Shinozaki and Mori, 2006). Instead, even though it is not an exact solution, using only the two rightmost (i.e., dominant) complex conjugate characteristic roots (s 0 = −0.3070 + 1.9176i and s −1 = −0.3070 -1.9176i), the free solution can be approximated as: Then using algebraic manipulation, one can get Eq. 12: With the numbers in Eq. 11, the resulting damping ratio is ζ = 0.1581. Also, the damped natural frequency: On the other hand, the response of the system in Eq. 1 can obtained by running simulation using Matlab/Simulink as shown in Fig. 3 for the same parameters. Then, the result response, x(t), is depicted in Fig. 4. As seen in the figure, the peaks are:  The logarithmic decrement δ is defined as the natural logarithm of the ratio of two successive amplitudes; that is: peak1 1 peak 2 2 peak 2 2 peak3 3 In In where, P is the period of the response, x(t). Then, the damping ratio can be obtained by using Eq. 10 and 14 as (Palm, 2009): For the values x peak1 and x peak2 in Eq. 13, δ = 1.0043 and, thus, ζ = 0.1578. The frequency ω d is also estimated from the response in Fig. 4 as ω d = 2π/P = 1.9156 where the period is P = t 2 −t 1 = 3.28. The values obtained by using the Lambert W function and measuring peaks of the response are close.
Estimation of time-delays: Assume that the timedelay, h, is unknown and the response is measured as in Eq. 13. Then, using the analysis the unknown parameter, h, can be estimated by following the steps: Step 1: Measure peak amplitudes and times (e.g., Eq. 13 Step 2: Calculate ζ and ω d (e.g., Eq. 15) Step 3: Calculate the root 2 n n 1 i ξω + ω − ξ Step 4: Solve the equation Note that in Step 4, the equation is nonlinear function with only one unknown, h. The equation can be solved nonlinear solvers (e.g., fsolve in Matlab). As shown in Fig. 5, when the results in Eq. 13 are used (see Fig. 5), the estimated time-delay is 1.0014 and close to the actual value (h = 1).

Example II: Logarithmic Decrement for Estimation (Non-Oscillatory)
Consider the system in Eq. 1 with parameters, a = −1, a d = 0.5 and h = 0.5. Then the simulated response of the system is shown in Fig. 6 The estimation method is very similar to the previous example. However, unlike that of Example 1 in Fig. 4, the response in Fig. 6 decays simply over time with any peak. Thus, instead of peak values as in Eq. 13 the value of x(t) is measured for every one second. Then, the values are: The decay rate γ is obtained as the average of the natural logarithm of the ratio of two successive amplitudes; that is Eq. 17: Using the values in Eq. 16, the value decay rate is γ= -0.3961. Because the rightmost characteristic root the system is single and real, the free solution can be approximated as Eq. 18: The Eq. 19 has one unknown and, thus, can be solved using nonlinear solvers. When using fsolve in Matlab with 1 as the initial guess for the time-delay, h, the iteration is shown in Fig. 7. The estimated timedelay is 0.4766 and close to the actual value (h = 0.5).

CONCLUSION
This study has presented a new tool for the timedelay estimation of 1st order time-delay systems. Use of the Lambert W function enables solving for the characteristic roots and identifying the rightmost ones. From the responses in the time-domain the real and imaginary parts of the rightmost characteristic roots are approximated. Then by equating them to the rightmost characteristic roots in terms of the Lambert W function, one can formulate estimation of time-delays. The provided examples illustrate ease of use of the method and effectiveness considering difficulty in estimating time-delays. The presented method will be extended to general systems of DDEs (higher than 1st order) using the matrix Lambert W function (Yi et al., 2010a) and will be applied to delay problems in network systems and fault detection of actuators by authors.