{'fraction_non_alphanumeric': 0.042117014258892777, 'fraction_numerical': 0.01862282437712488, 'mean_word_length': 4.2234296315928965, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Some Special Cases in the Stability Analysis of Multi-Dimensional Time-Delay Systems Using The Matrix Lambert W function

This paper revisits a recently developed methodology based on the matrix Lambert W function for the stability analysis of linear time invariant, time delay systems. By studying a particular, yet common, second order system, we show that in general there is no one to one correspondence between the branches of the matrix Lambert W function and the characteristic roots of the system. Furthermore, it is shown that under mild conditions only two branches suffice to find the complete spectrum of the system, and that the principal branch can be used to find several roots, and not the dominant root only, as stated in previous works. The results are first presented analytically, and then verified by numerical experiments.


Introduction
The stability analysis and control synthesis of Linear Time Invariant-Time Delay Systems (lti-tds), represented by Delay-Differential Equations (ddes), is a wide open area of research. The difficulty of this problem arises from the the fact that the delay makes this class of systems infinite dimensional. A nice review of the recent results and challenges int his area can be found in [17]. Several avenues have been followed by different researchers to find solutions to the stability and control questions. Some works study the absolute stability regions with respect to the time delay [15,18], and lead to control strategies that use the time delay as a stabilizing tool [14]. Other researchers focus on the numerical computation of the characteristic roots of the system. These works include approaches based on the discretization of the solution operator [3,5,8] or its infinitesimal generator [4,21], and methods based on root finding of the characteristic equation [19]. From here methods have been proposed to optimize the location of the dominant roots to guarantee a certain performance [12,13]. Finally, a Krylov method for computing characteristic roots of large-scale problems has been proposed in [11].
In the past decade a framework for analyzing ddes based on the Lambert W function has been developed [1,23,27]. It expands the earlier work [20]. The main idea of the methodology is to express the solution of a dde as the sum of a series of infinitely many exponential functions. The characteristic roots of the system are found analytically in terms of the Lambert W function. While the problem remains infinite dimensional, a one-to-one correspondence between the characteristic roots of the system and the branches of this multi-valued function is assumed. The stability question is then solved by earmarking the dominant characteristic root of the system with the principal branch of the Lambert W function. Therefore, only one branch has to be considered to determine whether a solution is stable or not. Furthermore, the existence of an explicit solution expressed in terms of a power series, allows analyzing structural properties like observability and controllability for systems of ddes [28], the development of pole placement techniques for control synthesis [24][25][26]31] and other applications like the estimation of decay rates [7].
The basic foundation of his methodology, i.e., the assumption that the principal branch of the Lambert W function defines the stability of the system, is well established for first order systems [1,16]. For higher order systems, however, this result has not been extended with the same rigor, and is rather based upon observations [29,30]. In this paper, we show that the main assumption does not hold in general. By studying particular second order delay systems, with a structure that is very common in applications, we show that the full spectrum of a lti-tds of this class can be found using only two branches of the Lambert W function, i.e., there is no one-to-one correspondence between the characteristic roots and the branches of the Lambert W function. This is due to the fact that an important nonlinear equation in the approach does not have a unique solution. Furthermore, we show that the principal branch can be used to find not only the dominant root of the system, but some other roots too.
The organization of the paper is as follows. In section 2 we present the Lambert W function and its matrix version, required to work with higher order systems. Section 3 reviews the methodology to solve ddes using the Lambert W function, both in the scalar and vector cases. Section 4 presents the analysis of a second order system and contains the main results. This analysis is illustrated by numerical examples in Section 5.
In the remainder of the paper, scalar quantities are denoted by italic symbols (a, b, λ) whereas vectors and matrices are represented by bold face lowercase (x) and uppercase (A, B) letters, respectively. The notation e b is used to represent the exponential function of a scalar and exp (A) represents a matrix exponential function. In a similar way, W k (z) represents the k-th branch of the Lambert W function of a scalar number, and W k (A) is the matrix Lambert W function.

The Lambert W Function
The Lambert W function is a function W (z), C → C, defined as the solution to the equation This is a multi-valued function, that is, for a z ∈ C there are infinitely many solutions to (1). To identify these values a branch number is assigned, and we refer to W k (z) as the k-th branch of the Lambert W function of z.
The branch cuts are defined in such a way that each branch has a precisely defined range [6]. For z ∈ R, only two of the branches are real valued. The principal branch, W 0 (z) is real for z ≥ −1/e and its range is the interval [−1, ∞). The branch W −1 (z) is real for −1/e ≤ z < 0, and its range is (−∞, −1]. These branches are shown in Figure 1. A comprehensive study of the history, definition and properties of the Lambert W function is found in [6]. The matrix Lambert W function is defined now. Consider a matrix H ∈ C n×n , which has the Jordan canonical decomposition H = ZJZ −1 , with J = diag (J 1 (λ 1 ) , J 2 (λ 2 ) , . . . , J p (λ p )). Following one of the standard definitions for a function of a matrix [9], the Matrix Lambert W function for a Jordan block of size m is defined as: and the matrix Lambert W function of H is defined as: Every matrix defined by (3), for k = 0, ±1, ±2 . . ., is one particular solution to the matrix equation The above, standard definition implies that the same branch of the Lambert W function is used in each Jordan block. This is not necessary to have a solution of (4). Since W 0 (0) = 0 and W 0 (0) = ∞ for k = 0, we deviate from the standard definition to avoid the infinite value. This special case, called the hybrid branch case, is defined in [23,27]. More precisely, there W k (H) uses the actual value of k for those Jordan blocks with λ = 0 and k = 0 for those blocks in which λ = 0. This definition is used in the remainder of the paper.
In a similar way, it is assumed that e −1 is not eigenvalue of H corresponding to a Jordan block of dimension higher than 1 when the principal branch is being computed. This is required to overcome the difficulty represented by the fact that W 0 e −1 is not defined. This limitation reduces the elegance of the definition of the matrix Lambert W function, but does not affect its usage [10].
3 Solution of Delay-Differential Equations using the Lambert W Function

Scalar Case
For a scalar, homogeneous dde, the characteristic equation is The solution to (6) can be expressed in terms of the Lambert W function following simple steps [1,6,23,27]. This solution has the form where k = 0, ±1, ±2, . . . indicates the branch of the Lambert W function to be used. Each one of the infinitely many roots of (6) corresponds to one of the branches of this function. It has been proven [16] that among all solutions in (7), the one that corresponds to the principal branch, k = 0, always has the largest real part and, therefore, introduces the dominant mode to the solution of the equation. To study the stability of the solution to a one dimensional dde as (5), it is necessary and sufficient to find only the solution of (6) in (7) corresponding to k = 0.

Higher Order case
Consider now a higher order DDE described by: The following steps, introduced in [1] and extended in [23,27,29,30], aim at computing characteristic roots using the matrix Lambert W function. The proposed method is based on finding a solution of the equation where S ∈ C n×n . A matrix Q is introduced, such that is satisfied. Let us define M := τ BQ. Then, from (10) and (4), with k ∈ Z, is a solution of (10). By substituting (11) into (9) the following expression is obtained, Therefore, the steps to compute characteristic roots are given by the following algorithm [23,27].
For system (8) to be stable, all characteristic roots must have negative real parts. Calculating the solution for all the branches is not possible. To work around this difficulty, the proposers of this methodology assume that for any branch k of the matrix Lambert W function, there is a unique solution M k to (12) and a corresponding S k matrix. This assumption, based on observations from many examples, leads to a stronger conjecture: the characteristic roots with largest real part correspond to the S 0 matrix, found using the principal branch of the matrix Lambert W function in Algorithm 1. This conjecture is formally stated in [27] and it is the basis for several derivative works [7,22,[24][25][26]31].
In the following section, we show that this conjecture does not hold in general. By studying a particular, yet very common, control system, we can show that (12) may have multiple solutions and there is no one to one correspondence between the branches of the matrix Lambert W function and the characteristic roots of the system. Furthermore, we show that is in particular cases is possible to find all characteristic roots using only two branches, those corresponding to k = 0 and k = −1. This contradicts the conjecture regarding stability.

A Common Special Case in Second Order Systems
We consider a single input, second order system under time delayed PD feedback. This type of control system is ubiquitous in applications. Without loss of generality, we assume that the system matrix is in companion form. This system is represented by (8) with the following A and B matrices: From the structure of B, we can see that M k = τ BQ k , for any given Q k , has the form: Applying the definition of the matrix Lambert W function, using the hybrid branch case because M k has one eigenvalue equal to zero, we get for m 22 = 0: and from here we obtain In case m 22 = 0, m 21 = 0, a simple computation yields where we used W 0 (0) = 1. We are now ready to state the main results of the paper.
Proposition 1 Let A and B be given by (13). Let {λ,λ} be any pair of complex conjugate characteristic roots of the system defined by (13). Assume their multiplicity is one. Then for either k = 0 or k = −1 there exists a real solution of (12), such that, if this solution and corresponding value of k are selected in the first step of Algorithm 1, the characteristic roots λ andλ are found in the last step of the algorithm.

Proof.
The key idea is to perform the steps of Algorithm 1 in reverse order, in the course of which k is selected. From the pair (λ,λ) we first construct a real matrix S k , of which they are the eigenvalues, namely Subsequently, we construct M k from (17), where we make distinction between two cases.  (16) and (17), we can take M k of the form (14), where m 21 ∈ R and m 22 ∈ R are chosen such that the following equations are satisfied: Such a choice is always possible for either k = 0 or k = −1, which then fix k. Equation (18) namely implies that W k (m 22 ) must be a real number.
As mentioned earlier, the definition of the branch cuts of the Lambert W function makes this function to have two real branches: k = 0, the principal branch, and k = −1, and the union of the ranges of W 0 and W −1 for real valued arguments includes R, see Figure 1. Furthermore, in the case considered, m 22 computed from (18) is different from zero, justifying the use of (16).
Case 2: 2 (λ) = a 22 . We can freely choose k ∈ {0, −1} and take Finally, it remains to show that the constructed pair (k, M k ) is a solution of (12). Let v, respectivelyv be the eigenvector corresponding to λ, respectivelyλ. If the characteristic root λ is simple, then the pair (V, Λ) is an invariant pair of (8), where

As a consequence it satisfies
see [2]. The eigenvalue decomposition of S k in (17) takes the form It follows that Substituting the latter in (20) yields, as V is invertible, Proposition 2 Let A and B be given by (13). Let λ 1 = λ 2 be two real, simple characteristic roots of the system defined by (13). Then for either k = 0 or k = −1 there exists a real solution of (12), such that, if this solution and corresponding value of k are selected in the first step of Algorithm 1, the characteristic roots λ 1 and λ 2 are found in the last step of the algorithm.
Proof. The proof is completely analogous to the proof of Proposition 1.
The differences are that we start by defining and that the two cases to be considered are λ 1 + λ 2 = a 22 and λ 1 + λ 2 = a 22 . 2 From Propositions 1-2 the following Corollary can be derived.

Corollary 1 Let
A and B be given by (13). If all characteristic roots of (13) are simple and if the number of real characteristic roots is different from one, then all characteristic roots can be found using only two branches of the matrix Lambert W function in Algorithm 1, namely k = 0 and k = −1.
Moreover, one can restrict to the real solutions of (12).

Multiple Solutions
The following example uses a system defined by matrices: and a value τ = 5 for the time delay. We tried to calculate the characteristic roots of the system following the matrix Lambert W function approach. For this we used the LambertDDE toolbox [22]. However, the software was not able to find a solution of (12) for any value of k with the default settings. For the numerical solution, the toolbox uses, as an initial estimation of Q k , matrix exp −Aτ . In this case, that value turns out not to be in the region of attraction of a solution of (12).
In order to obtain an a posteriori guess for Q k , we reverse-engineered the solution, as in the proof of Proposition 1. First, we use the QPmR algorithm [19] to find the characteristic roots of the system in a region close to the origin of the complex plane. The roots found are shown in Figure 2.
The dominant roots of the system are λ = 0.0377±j1.7911. Corresponding to these roots, we create the following S matrix:  .
This shows that W (m 22 ) ∈ [−1, ∞), which is the range of the principal branch of the Lambert W function. There is, therefore, a matrix M for which (23) is satisfied for k = 0, and that matrix is Since B and M are singular, there are an infinite number of Q 0 matrices that satisfy M 0 = τ BQ 0 , for (24). One of such matrices is When this Q 0 is used as starting value in the LambertDDE toolbox for k = 0, the numerical solution of (12), corresponding to the dominant roots, is found at the first iteration, as expected. Furthermore, if the matrix is slightly perturbed, the method still converges to the same solution after a few iterations. Now, let us consider a non-dominant pair of roots: λ = −0.4113 ± j6.4803. Following a similar reasoning, we obtain the following matrices: For this case, we again have W (m 22 ) ∈ [−1, ∞). This implies that this pair of roots can also be found using the principal branch of the matrix Lambert W function. Indeed, if a matrix close to this one, which was created as in the previous case, is used as initial condition, the numerical routine within the LambertDDE toolbox converges to the solution using the principal branch.
In fact, we have observed that using suitable initial conditions, the 11 pairs of roots presented as blue squares in Figure 2 can be found using the principal branch of the matrix Lambert W function.
If we consider now the pair of roots λ = −0.6169 ± j14.0734, the corresponding S and W k (M k ) are In this case, W (m 22 ) ∈ (−∞, −1], which is the range of the branch indexed by k = −1. Using this branch, we obtain a Q −1 matrix, which is a solution to (12) for k = −1. Selecting an initial condition close to this matrix guarantees convergence to this solution. This procedure can be repeated for all the roost marked as black squares in Figure 2, as well as for roots further to the left of the complex plane, always using the branch corresponding to k = −1.
This example shows how the whole spectrum of a system with the structure given in (13) can be calculated using only two branches of the matrix Lambert W function and properly selected initial conditions for the solution of the nonlinear equation (12).

One Real Characteristic Root
Consider the systems with whose characteristic equation is given by From (30) it can be seen that there is only one real characteristic root, at the origin and with multiplicity one. Furthermore, this is the rightmost root, as can be seen in Figure 3. An analysis of the applicability of the Lambert W approach to this example leads the following observations.
• For this particular system B has a repeated zero eigenvalue. According to the observations in [23,27], in this case the dominant root should be found using the principal branch or the branches k = ±1.
• The "reverse-engineering" approach, behind Propositions 1-2 and based on constructing real S k and M k , cannot be used starting from the characteristic root at zero. The proof is by contradiction. Assume that there is a real solution M k , and corresponding real S k , defined by (11), having eigenvalue λ = 0. By (11) and (12) this implies that S k satisfies (9). This implies on its turn that (I, S) is an invariant pair of (8), i.e., the eigenspace corresponding to the real characteristic roots is two-dimensional.
• When using the LambertDDE toolbox no solution is found for k = 0, despite the root at λ = 0 is the dominant one. However, the root at the origin is reported, together with one of the complex roots, when branches k / ∈ {0, −1} are used. This contradicts the observations mentioned in the first point.
Whether the two phenomena are related is a topic of further research.

Concluding Remarks
The methodology for the stability analysis and control synthesis of lti-tds presented in [23,27], based on the matrix Lambert W function, assumes that there is a one to one correspondence between the characteristic roots of the system and the branches of the matrix Lambert W function. It also assumes that this correspondence is such that using the principal branch always leads to the dominant roots of the system. This paper illustrates that such a correspondence does not exist in all the cases, and that for a particular, albeit very common structure for the system, the branches corresponding to k = 0 and k = −1 can be used to find all the characteristic roots of the system. It is also illustrated that the correspondence cannot always be restored when selecting the particular initial condition exp(−Aτ ) in solving (12).
In our examples matrix B had reduced rank. An interesting path of research is to find out whether the correspondence between characteristic roots and branches of the Lambert W function, holding for scalar systems, can be extended to higher-order systems, provided additional structural conditions on the system are assumed.
The MATLAB code used to create the examples is available at https://db.tt/mSI3VwbO or can be requested via email to the authors.