Convergence aspects of the matrix Lambert W function method for solving systems of delay differential equations (DDEs) are considered. Recent research results show that convergence problems can occur with certain DDEs when using the well-established Q-iteration approach. A complementary, and recently proposed, W-iteration approach is shown to converge even on systems where the Q-iteration fails. Furthermore, the role played by the branch numbers k = -∞ .. -2, -1, 0, 1, 2 , .. ∞ of the matrix Lambert W function, Wk, in terms of initializing the iterative solutions, is also discussed and elucidated. Several second order examples, known to have convergence problems with Q-iteration, are readily solved by W-iteration. Examples of third and fourth order DDEs show that the W-iteration method is also effective on higher-order systems. 1. Introduction This paper considers the matrix Lambert W function approach, developed during 2003-2006, to the solution of systems of delay differential equations (DDEs) [Asl & Ulsoy 2003, Yi & Ulsoy 2006]. The method solves linear systems of DDEs with a single known delay. It can be used to assess stability, obtain free and forced solutions, determine observability and controllability, and for controller/observer design via eigenvalue assignment [Yi et al 2010a]. The method, explained in more detail in the next section, relies on a numerical search for an unknown matrix Q that satisfies certain conditions. In this paper this original numerical iteration is referred to as the Q-iteration approach. The method has been successfully applied to a variety of time delay systems as previously documented in numerous publications [Yi et al 2007, 2008, 2010a]. Recent studies have, however, demonstrated convergence problems with the Q-iteration approach on some specific problems [Wei et al 2014, Cepeda-Gomez & Michiels 2015]. Researchers have also recently proposed a complementary approach, which replaces an iteration to find the unknown matrix Q with iteration 1 Professor Emeritus, Department of Mechanical Engineering, ulsoy@umich.edu 2 Academic Affiliate, Department of Mathematics, ritagtk@umich.edu DRAFT - 11/5/18 2 of the matrix Lambert W function itself [Ivanoviene & Rimas 2015]. This complementary approach will be denoted in this paper as the W-iteration approach. Furthermore, there are recent research results that demonstrate that it is not necessary to use all the branches of the multi-valued matrix Lambert W function [Cepeda-Gomez & Michiels 2015]. The complex multi-valued matrix Lambert W function, Wk, has an infinite number of branches denoted by the branch numbers k = -∞ .. -2, -1, 0, 1, 2 , .. ∞. The work by Cepeda-Gomes and Michiels shows that if the initial value of Q is chosen appropriately, then it is sufficient to only use the branches k = 0 and -1. The purpose of this paper is to better understand these convergence issues and their potential resolution, as well as to discuss and elucidate the role that the branches of the multi-valued matrix Lambert W function play in these iterative solutions. This paper demonstrates, for the first time, that the well-documented convergence problems in the literature can be resolved by using the complementary W-iteration rather than the original Q-iteration method. Furthermore, the role of multiple branches, rather than only the k = 0 and -1 branches, in initializing these iterative solutions is elucidated based upon the multi-valued nature of the Lambert W function. 2. The Lambert W Function and Related Approaches to Solving DDEs 2.1 The Lambert W function Introduced in the 1700's by Lambert and Euler, the Lambert W function is defined to be any function, W(H), that satisfies [Corless et al 1996]: W(H)eW(H)=H (1) In general the argument H as well as the function Wk , which has an infinite number of branches denoted by the branch numbers k = -∞ .. -2, -1, 0, 1, 2 , .. ∞, are complex valued. The W function is similar to the complex logarithm function, P(H) = log H, and its inverse, H = eP(H), which has an infinite number of branches, k, and has branch cuts which are horizontal lines in the complex P – plane, defined by P = (2k+1)πi, where i = -∞ .. -2, -1, 0, 1, 2 , .. ∞. The curves which define the branch cuts for the Lambert W function are a subset of the so-called Quadratrix of Hippias, as discussed in more detail in [Corless et al 1996]. Commercial software packages, such as MATLAB, Maple, and Mathematica, have embedded commands for numerically computing the Lambert W function [Corless et al 1996]. Plots of the complex Lambert W function, showing its repetitive nature (similar to the complex logarithm function) are given in Fig. 1 for branches k = -2, -1, 0, 1, 2. Extension of the Lambert W function to a matrix Lambert W function, denoted as W(H) has been presented in several studies based on the use of the DRAFT - 11/5/18 3 Jordan canonical form of the complex matrix argument, H and is summarized here in Appendix A [Asl & Ulsoy 2003, Jarlebring & Damm 2007, Yi et al 2010a]. Figure 1. Three Dimensional Plot of the (a) Real Part and (b) Imaginary Part of the scalar Lambert W Function, Wk(x+yi) for Branches k = -2, -1, 0, 1, 2 Figure 2. Re(Wk) vs x = Re(Hk) for branches k = -2, -1, 0, 1, 2 100500xreal(lambertw(-2, x + y i)) -50-100-100-80-60-40-200y20406080-1.5-1-0.500.511.522.5 100 (a) Re(Wk) 100500ximag(lambertw(2, x + y i)) -50-100-100-80-60-40-200y20406080-202468101214 100 (b) Im(Wk) -100-80-60-40-20020406080100-3-2-101234 x = Re(Hk) Re(Wk) k=0 k=1 k=2 k=-1 k=-2 DRAFT - 11/5/18 4 Figure 2 is a simplified 2-dimensional plot, based on Fig. 1 (a), that shows the real part of the complex scalar Lambert W function, Re(Wk), versus the real part of its complex scalar argument, x = Re(Hk), for the branches k = -2, -1, 0, 1 and 2. Because of the multi-valued nature of this function a single value of the function (e.g., Re(Wk) ≈ 2) can be obtained from various branches with different arguments (e.g., W0(10), W1(35), W-1(-40), W2(65), W-2(-70), etc.). A range of values for Re(Wk), such as 1.5 to 3, can be obtained by either: (1) using a single branch (e.g., W0) with a range of arguments (e.g., H0 = 5 to 40), or (2) using multiple branches (e.g., k = 0, 1, 2) with a single argument (e.g., H ≈ 40). This multi-valued nature of the Lambert W function is relevant to the initialization of the numerical solution methods (i.e, the Q-iteration and W-iteration) discussed in the following sections. 2.2 Solution of DDEs and the Q iteration approach [Yi et al 2010a] The Lambert W function in Eq. (1) can be used to solve the scalar DDE: !x(t)=ax(t)+adx(t−h) (2) where h is a single known delay. Assuming the exponential solution form: x(t)=estC (3) and substituting Eq. (3) into Eq. (2) one obtains the characteristic equation: s−a−adesh=0 (4) Multiply both sides of Eq. (4) by heshe−ah to obtain: h(s−a)eh(s−a)=hade−ha (5) Note that Eq. (5) is of the same form as Eq. (1) if one defines W=h(s−a) (6) and the argument H=hade−ha . Consequently, the solution to the characteristic equation (i.e., the eigenvalues) are found by solving Eq. (6) for s for each branch k = -∞ .. -2, -1, 0, 1, 2 , .. ∞. sk=(1/h)Wk(hade−ha)+a (7) It has been proven in [Shinozaki & Mori 2006] that the infinite eigenspectrum of Eq. (2) is given by Eq. (7) and that the principal (i.e., k=0) branch gives the rightmost eigenvalue in the s-plane and, thus, determines system stability. For a system of DDEs, analogous to Eq. (2), represented as: DRAFT - 11/5/18 5 !x(t)=Ax(t)+Adx(t−h) (8) following the same approach, and corresponding to Eq. (5), one obtains: h(S−A)ehSe−hA=hAde−hA (9) If S and A commute, then ehSe−hA=eh(S-A) and we can write the solution, analogous to Eq. (7) in terms of the matrix Lambert W function, as: S=(1/h)W(hAde−hA)+A (10) However, in general S and A do not commute, and Eq. (10) is not a valid solution to Eq. (9). To address this issue, one can introduce an unknown matrix Q to satisfy: h(S−A)eh(S−A)=hAdQ (11) Thus the solution matrix, S, comparing Eq. (1) to Eq. (11), can be obtained as: S=(1/h)W(hAdQ)+A (12) Substituting Eq. (12) into Eq. (9), one obtains the following condition that the unknown matrix Q must satisfy: W(hAdQ)eW(hAdQ)+hA=hAd (13) Consequently, the system of nonlinear equations in Eq. (13) can be initialized as Q = Qini and then solved iteratively for Q. In this paper this solution method is referred to as the Q-iteration approach. Since Eq. (13) is a nonlinear system of equations with multiple solutions, the choice of initial condition, Qini, will determine the specific solution Q that is obtained. In many problems reported in the literature Qini=e−hA is used, and this corresponds to the value that Q would take if S and A did commute (which they do not). If the same argument (i.e., same value of Qini) is used with different branches then one can obtain different values of Wk corresponding to different eigenvalues. Alternatively one can obtain those different eigenvalues by using only one branch, but varying the argument. This follows from the multi-valued nature of the Lambert W function and from our previous discussion of the simplified plot in Fig. 2. For example, using only the fixed argument Re(Hk) = 20, and all the branches (i.e., k = -2, -1, 0, 1, 2) one obtains values for Re(Wk) in the range of 0.8 to 2.5. Alternatively, using only branch k = 0 one can obtain the same range of values for Re(W0) by using different arguments Re(H0) in the range 5 to 20. Note that when Ad is rank deficient, which is often the case in natural and engineered time delay systems, the matrix Q obtained as the solution of Eq. (13) will not be unique. This can lead to numerical convergence issues in some particular DRAFT - 11/5/18 6 cases, as documented in the literature [Wei et al 2014, Cepeda-Gomez & Michiels 2015]. Finally, as discussed in the previous subsection, the matrix Lambert W function is a multi-valued complex function, with complex argument H, and an infinite number of branches denoted by the branch numbers k = -∞ .. -2, -1, 0, 1, 2 , .. ∞ [Corless et al 1996, Asl & Ulsoy 2003]. Corresponding to each branch, k, of the Lambert W function, Wk, there is a solution Qk of Eq. (13), and an argument Hk = hAdQk. The infinite eigenspectrum of the system of DDEs in Eq. (8) is obtained from the eigenvalues of the solution matrices: Sk=(1/h)Wk(hAdQk)+A (14) For the scalar case of Eq. (2), the rightmost eigenvalues can be determined by the principal branch, i.e., k = 0, of the scalar Lambert W function [Shinozaki & Mori 2006]. Such a proof can readily be extended to the matrix case when A and Ad in Eq. (8) commute [Jarlebring & Damm, 2007]. No such proof is currently available for the general case of matrix DDEs. However, in all the examples considered in the literature, it has been observed that the rightmost eigenvalue is always obtained when using only the first k = -m, -(m-1), .., 0, .., (m-1), m branches, where m is the nullity of Ad [Yi et al 2010a, Duan et al 2012]. 2.3 Complementary W-iteration approach [Ivanoviene & Rimas 2015] The essence of the approach proposed in [Ivanoviene & Rimas 2015], which will be denoted here as the W-iteration approach, is that Eq. (13) can be simplified and written as: DkeDk+hA=hAd (15) and iteratively solved for an unknown matrix Dk, rather than an unknown matrix Qk , for each branch k. The iteration is initialized, using the matrix Lambert W function for each branch k, as Dini = Wk(hAde-hA). As noted previously, this choice of initial condition matrix corresponds to the value that D = W would take if S and A did commute (which they do not). The iteration for Dk converges to the matrix Wk in Eq. (14). Consequently this complement to the original Q-iteration method is referred to here as the W-iteration method and has the following advantages [Ivanoviene & Rimas 2015]: 1. The initial matrix Dini is always chosen in the same way as Dini = Wk(hAde-hA). So the issues of how to initialize Q and whether it will converge are avoided. 2. The problem of non-uniqueness of Q in the original Q-iteration method, when the matrix Ad is rank deficient, is also avoided. DRAFT - 11/5/18 7 3. Numerical Examples Several numerical examples are considered in this section to show that: (1) the W-iteration method works well on problems that were previously shown to have convergence problems with the Q-iteration method, and (2) the W-iteration method works well, not only on second order, but also higher order systems of DDEs. The basic MATLAB code used to solve these numerical examples is given in Appendix B. 3.1 Example from [Cepeda-Gomez & Michiels 2015] Consider the system of DDEs, as in Eq. (8), with a delay h = 5 and A=01−5−1⎡⎣⎢⎤⎦⎥;Ad=00−3−0.6⎡⎣⎢⎤⎦⎥ (16) For this system the Q-iteration method was shown to not converge, with initialization Qini=e−hA [Cepeda-Gomes & Michiels 2015]. However, if the solution matrix S is known a priori, then an appropriate initialization of Q can be found for which the Q-iteration method will converge. Unfortunately, if the solution S is known a priori then, of course, there is no need for the Q-iteration at all! When the W-iteration method, with initialization Dini = Wk(hAde-hA) is used, the solution is readily obtained. The dominant (i.e., rightmost) eigenvalues of the system are the complex conjugate pair located at λ1,2=0.0377±1.7911i. As noted previously, in all the examples considered in the literature it has been observed that the rightmost eigenvalue is always obtained by using only the first k = -m, -(m-1), .., 0, .., (m-1), m branches, where m is the nullity of Ad [Yi et al 2010a, Duan et al 2012]. Since rank(Ad) = 1, i.e., m = nullity of Ad is 1, then we need only to consider the branches k = -1, 0, and 1 to ensure that these dominant (i.e., rightmost) eigenvalues are found. The MATLAB code in the Appendix quickly converges and yields the results summarized in Table 1 below. k Initial Condition for Wk Solution Matrix, Sk Eig(Sk) -1 000.66+1.28i−2.20−4.23i⎡⎣⎢⎤⎦⎥ 01−4.962+0.141i0.017−0.979i⎡⎣⎢⎤⎦⎥ 0.038 + 1.791i -0.020 - 2.771i 0 00−0.110.37⎡⎣⎢⎤⎦⎥ 15.1403.512−72.433−16.3953⎡⎣⎢⎢⎤⎦⎥⎥ -0.628 + 2.403i +1 000.66−1.28i−2.20+4.23i⎡⎣⎢⎤⎦⎥ 01−4.962−0.141i0.017+0.979i⎡⎣⎢⎤⎦⎥ 0.038 - 1.791i -0.020 + 2.771i Table 1. W-iteration solution for example from [Cepeda-Gomes & Michiels 2015] DRAFT - 11/5/18 8 It is evident from these results that the convergence problems observed in [Cepeda-Gomez & Michiels 2015] when using the Q-iteration method have been eliminated by using the W-iteration method. Furthermore, the observation that only k = -m, .. 0. .., m branches of Wk are needed to obtain the dominant (i.e., rightmost) eigenvalues is once again confirmed, since m =1 and the dominant eigenvalues are 0.038+1.791i. Note also that the initial Dini = Wk(hAde-hA) matrix is complex, as is the resulting solution matrix, Sk, for branch k = i and their conjugates are found for k = -i. The eigenvalues of the system of DDEs are either real or occur in complex conjugate pairs, thus, both branches i and –i must be considered to find those complex conjugate pairs. Thus, the eigenvalues obtained in sequence from right to left in the s-plane are 0.038 + 1.791i, -0.020 + 2.771i, and -0.628 + 2.403i. 3.2 Example from [Wei et al 2014] Consider, with h = 1, the system of DDEs defined by: A=01−2.52.5⎡⎣⎢⎤⎦⎥;Ad=002.50⎡⎣⎢⎤⎦⎥ (17) This system corresponds to a region in the parameter plane where the Q-iteration method does not converge [Wei et al 2014]. However, when the W-iteration method is used it converges quickly to the results summarized in Table 2. Again, for this example m = 1, so we show results here only for the branches k = -1, 0, 1. k Initial Condition for Wk Solution Matrix, Sk Eig(Sk) -1 00−0.252+2.804i0.432−4.802i⎡⎣⎢⎤⎦⎥ 0.710011.266+0.823i−15.156−1.159i⎡⎣⎢⎢⎤⎦⎥⎥ 0.710 -15.156 – 1.159i 0 00−0.9111.560⎡⎣⎢⎤⎦⎥ 0100.710⎡⎣⎢⎤⎦⎥ 0 0.710 +1 00−0.252−2.804i0.432+4.802i⎡⎣⎢⎤⎦⎥ 0.710011.266−0.823i−15.156+1.159i⎡⎣⎢⎢⎤⎦⎥⎥ 0.710 -15.156 + 1.159i Table 2. W-iteration solution for example from [Wei et al 2014] The dominant, or rightmost, eigenvalue is 0.710 and is found from all the three branches considered. As in the previous example, the initial condition matrix and the solution matrix are complex, and those for branch k = 1 are the conjugates of the ones for branch k = -1. Thus, the eigenvalues obtained in sequence from right to left in the s-plane are 0.710, 0, and -15.156 + 1.159i. DRAFT - 11/5/18 9 3.3 Higher order examples from [Yi et al 2010, Ulsoy 2015] In addition to the two second-order system examples in Sections 3.1 and 3.2, this section includes examples of third and fourth order systems of DDEs [Yi et al 2010a, Ulsoy 2015]. First, consider, with h = 0.06, the system of DDEs defined by: A=−27−0.009769.5999−40.2750−40.6578018.06084.1480⎡⎣⎢⎢⎢⎤⎦⎥⎥⎥;Ad=0002100000⎡⎣⎢⎢⎢⎤⎦⎥⎥⎥ (18) This third order system corresponds to a closed loop control system for a diesel engine as described in detail in [Yi et al 2010b]. The rank of Ad is one, so m = 2. The eigenvalues of the solution matrix Sk, for each branch k = -2, -1, 0, 1, 2 are shown in Table 3. Branch, k Eigenvalues of Sk -2 -21.56 + 23.71i, -145.5 – 208.33i -1 -21.56 + 23.71i, -114.4 – 90.52i 0 -10.0, -21.56 + 23.71i 1 -21.56 + 23.71i, -114.4 + 90.52i 2 -21.56 + 23.71i, -145.5 + 208.33i Table 3 W-iteration solution for example from [Yi et al 2010b] The eigenvalues obtained, ordered from right to left in the s-plane, are -10, -21.56 + 23.71i , -114.4 + 90.52i, and -145.5 + 208.33i. The rightmost eigenvalue, -10, is found for the principal branch, k = 0. Although not shown here, as in the previous examples, the initial condition matrix and the solution matrix are complex, and those for branch k = 1, 2 are the conjugates of the ones for branch k = -1, -2 respectively. Next, consider, with h = 0.01, the system of DDEs defined by: A=00100001−500−117.8−67.11−10.1−0.1⎡⎣⎢⎢⎢⎢⎤⎦⎥⎥⎥⎥;Ad=000000000010000000⎡⎣⎢⎢⎢⎢⎤⎦⎥⎥⎥⎥ (19) This fourth order system corresponds to closed loop active control of a two-degree-of-freedom mechanical vibration problem as described in detail in [Ulsoy 2015]. The rank of Ad is one, so m = 3 and the eigenvalues of the solution matrix Sk, for the branches k = -3, -2, -1, 0, 1, 2, 3 are considered as shown in Table 4. The rightmost eigenvalue is at -1.75+0.34i, and is found for all the branches considered. Although not shown here, as in the previous examples, the initial condition matrix and the solution matrix are complex, and those for branches k = 1, 2, 3 are the conjugates of the ones for branches k = -1, -2, -3 respectively. DRAFT - 11/5/18 10 Branch, k Eigenvalues of Sk -3 -1.75 + 0.34i, -285 - 1718i, -1728 – 1.426i -2 -1.75 + 0.34i, -2.69 – 0.70i, -239.3 - 1088.4i -1 -1.75 + 0.34i, -2.69 – 0.70i, -153.6 - 463.5i 0 -1.75 + 0.34i, -2.69 + 7.00i 1 -1.75 + 0.34i, -2.69 + 0.70i, -153.6 + 463.5i 2 -1.75 + 0.34i, -2.69 + 0.70i, -239.3 + 1088.4i 3 -1.75 + 0.34i, -285 + 1718i, -1728 + 1.426i Table 4 W-iteration solution for example from [Ulsoy 2015] In addition to the examples presented in this section, the W-iteration method was also evaluated on several other examples from [Wei et al 2014, Cepeda-Gomes & Michiels 2015, Ivanoviene & Rimas 2015] and found to readily converge, and to always yield the dominant eigenvalues for branches k = -m, .., -1, 0, 1, .., m. Due to space limitations, those additional results are not included here. 4. Summary, Conclusions and Future Work This paper considered the solution of systems of DDEs using the matrix Lambert W function. The original Q-iteration and the complementary W-iteration methods were considered and applied to problems where the Q-iteration was known to have convergence problems. On those non-converging example problems, and on others, the W-iteration method was shown to converge to the expected results. Consequently, the solution of systems of DDEs using the matrix W function is shown to be a reliable engineering tool when W-iteration is employed. Furthermore, by considering the multi-valued nature of the complex Lambert W function, it was illustrated that the necessary iterations can be initialized by using only the principal branch or by using multiple branches of the function. However, using only the principal branch requires that the desired solution be know a priori, thus, making that approach impractical. Using multiple branches allows for multiple initializations of the iterative solution method, which then converge to different solutions in the infinite eigenspectrum. Based on the examples presented one can conclude that W-iteration works when Q-iteration does not, and should be preferred. One can use either multiple branches, or different arguments for the principal branch, to initialize the iteration. However, the second option is not practical, since it requires that the solution be known a priori. Furthermore, the long-standing observation that at most m (where m = nullity of Ad) branches are needed to ensure that one finds the rightmost eigenvalues is confirmed in all the examples considered in this paper. While these are important and encouraging results, which demonstrate that the W-iteration approach to solving systems of DDEs can be a powerful and reliable DRAFT - 11/5/18 11 engineering tool, some basic research questions remain. For example, the important and useful observation that only m branches are needed to ensure that one finds the rightmost eigenvalues has not been proven for general A and Ad. Another key question is whether the selection of different branches to initialize the W-iteration can ensure that all solutions will indeed be found. In short, there is a need to place this practical engineering method on a more sound theoretical foundation. 5.