A Novel Frequency-Domain Approach for the Exact Range of Imaginary Spectra and the Stability Analysis of LTI Systems With Two Delays

This paper presents a novel frequency-domain approach to reveal the exact range of the imaginary spectra and the stability of linear time-invariant systems with two delays. First, an exact relation, i.e., the Rekasius substitution, is used to replace the exponential term caused by the delays in order to transform the transcendental characteristic equation to a quasi-polynomial. Second, this quasi-polynomial is uniquely tackled by our proposed Dixon resultant and discriminant theory, leading to the elimination of delay-related elements and the revelation of the exact range of the frequency spectra of the original system of interest. Then, by sweeping the frequency over this obtained range, the stability switching curves are declared exhaustively. Last, we deploy the cluster treatment of characteristic roots (CTCR) paradigm to reveal the exact and complete stability map. The proposed methodologies are tested and verified by a numerical method called Quasi-Polynomial mapping-based Root finder (QPmR) over an example case.


I. INTRODUCTION
In this paper, the frequency spectra and the stability of a general class of linear-time invariant systems with two delays are analyzed from a new perspective: where x ∈ R n is the state vector; A, B 1 , B 2 are constant, known matrices in R n×n ; and τ 1 , τ 2 are rationally independent delays. In this paper, the boldface capital notation denotes vector and matrix quantities. Apparently, the characteristic equation of (1) is where s is the Laplace variable, and τ = (τ 1 , τ 2 ) ∈ R 2 is the delay vector that contains rationally independent delays.
The associate editor coordinating the review of this manuscript and approving it for publication was Zhiguang Feng .
The investigation of the stability of (2) in the domain of the delays is a basic yet challenging problem among researchers in the control community [1]- [5]. The notorious delay terms in (1) bring about the exponential terms rendering (2) transcendental and thus extremely difficult to be tackled. One of the most crucial and interesting problems is to reveal the stability feature of the time-delay systems in the domain of the delays. For this, a number of works have been proposed and verified successfully [6]- [9]. Based on the obtained stability maps, various control concepts have been developed with the aim to increase the delay robustness of systems against large delays, for example, the sign inverting [10], [11] and delay scheduling controls [12], [13]. These large delays can be several orders of magnitude larger than the sampling period for a particular system [14], [15], in comparison with the conventionally treated small delays that are a few sampling periods long [16].
However, among the prementioned successful attempts to tackle the stability analysis of time-delay systems, none proposed a direct algorithm, to the best of the authors' knowledge, to reveal the exact range of the frequency spectra of the general time-delay system in (1). Several other attempts have been indeed made for the stability investigation of time-delay systems with more than two delays and the bounds of the frequency spectra have been revealed as a by-product [17]- [19]. However, we wish to clarify here that the obtained frequency bounds in the cited works are by no means for the original time-delay system being dealt with, but rather for the two-dimensional cross section of the system of interest in the domain of two arbitrarily selected delays. In other words, the obtained frequency bounds are, in fact, a subset of the original frequency spectra. It is noteworthy that the revelation of the exact frequency range of (1) lays the foundation for its stability analysis, which serves as the main contribution of the paper. A series of techniques are proposed to determine the stability map including the Rekasius substitution [20], Dixon resultant [21] and discriminant [22], frequency sweeping technique [23], and the cluster treatment of characteristic roots CTCR) paradigm [8], [24], [25]. Specifically, the Rekasius substitution is a handy tool to convert the transcendental equation (2) into a quasi-polynomial, which, in turn, facilitates the deployment of numerous well-known linear-system techniques. The Dixon resultant and discriminant are then used successively to eliminate delay-related elements in the characteristic eqution to obtain the exact range of the frequency spectra of the original system (1). By sweeping such frequency range, the complete set of stability switching curves are determined. With this knowledge, the CTCR paradigm is utilized to declare the exact stability map. Finally, the proposed technique is verified over a published numerical case using a numerical tool called Quasi-Polynomial mapping-based Root-finder (QPmR) [26]. This numerical method is based on mapping the quasi-polynomial and on utilizing asymptotic properties of the chains of the roots, which are determined as the intersection points of the zero-level curves of the real and imaginary parts of the quasi-polynomial in a certain region of the complex plane.
The rest of the paper is structured as follows: Section II reviews some crucial definitions of the CTCR paradigm and the Rekasius substitution [20]. The main novel results on determining the range of the frequency spectra and the complete stability map for the time-delay system (1) are presented in Section III. Section IV shows the capability and strength of the proposed approach over a published example case in the literature using the QPmR method.

II. PRELIMINARIES
A general class of time-delay system (1) is declared asymptotically stable if and only if all its infinitely many characteristic roots reside in the left half of the complex plane. As the continuity of these characteristic roots with respect to delays has been proved previously [27], the domain of the delays is partitioned into stable and unstable regions if one restricts the discussion within the delay-dependent stability problem rather than the delay-independent one [28]. The boundaries separating the stable and unstable regions are the stability switching curves consisting of both the kernel and offspring curves (KOH) [29], which are, in fact, the only loci where the characteristic equation (2) exhibits purely imaginary roots. The complete knowledge of the KOH is essential for the declaration of the exact stability map of (1) in the domain of the delays. Actually, the KOH is a fundamental concept for the CTCR paradigm.

A. REVIEW OF THE CTCR PARADIGM
Before introducing the formal definitions of the KOH, we first present a key definition of the complete set of the frequency spectra, denoted as , of (2) as follows, where i is the imaginary unit, and the τ , ω notation implies that a purely imaginary root, ω ∈ R + of (2) exists for a specific combination of delays, τ ∈ R 2+ . Note that special cases in which ω = 0 or zero delays are not discussed within this paper. With this, the following definitions arise. Definition 1 (Kernel Hypersurfaces ℘ 0 ): The hypercurves in the R 2+ domain that consist of all the points τ ∈ R 2+ that cause a purely imaginary root of (2), i.e., s = ±ωi and satisfy the constraint 0 < τ k ω < 2π, k = 1, 2 are called the kernel hypercurves. The points on these hypercurves contain the smallest delay values that incur the given pair of imaginary roots at the frequency ω.
Definition 2 (Offspring Hypersurfaces ℘): The hypercurves obtained from the kernel according to the following point-wise nonlinear transformation: where p 1 , p 2 are not zero simultaneously, are called the offspring hypercurves. They are defined due to the periodicity of the purely imaginary roots with respect to the time delays.
It should be noted that the stability switching curves, i.e., KOH are actually ℘ 0 ∪ ℘. Another key concept follows these two definitions as below.
Definition 3 (Root Tendency (RT)): At any point τ ∈ ℘ 0 ∪ ℘, an infinitesimal increase in any of the delays τ j (j = 1, 2) creates either a leftward or rightward transition of the purely imaginary root across the imaginary axis in the complex plain. RT indicates the direction of the crossing as only one of the delays, τ j , increases by , 0 < 1, while the other remains fixed: With these key definitions, we will then present a basic yet crucial mathematical manipulation, i.e., Rekasius substitution, to initially convert the transcendental equation (2) into a quasi-polynomial.

B. REKASIUS SUBSTITUTION
As a preparatory step for the main contribution in the next section, we present an exact mathematical relation called Rekasius substitution [20], to replace the exponential term with a rational polynomial as follows: where T k ∈ R is called the agent parameter. The above equation is an exact relation where s = ±ωi, ω ∈ R + with the following connection between τ k and T k :

C. CONVERSION OF THE CHARACTERISTIC EQUATION
The substitution of (6) into (2) results in a rational polynomial , which can then be converted into a polynomial in terms of s with parameters in T as follows, It is notable that only the imaginary spectra of (8) are identical to those of (2) while the remaining complex roots of (8) are completely different from those of (2). With these preliminary results, we then present the main contribution of the paper in the next section.

III. MAIN RESULT
We first define a similar frequency set as in (3) for the converted characteristic equation (8) as: Next, we present a theorem to show the equivalence of these two frequency sets.
Theorem 1 (Identical Imaginary Spectra): Between the characteristic equation (2) for the original time-delay system (1) and the converted characteristic equation (8), the imaginary spectra are identical. In other words,ˆ = Proof: From the properties of the Rekasius substitution, (6) is an exact equation for purely imaginary roots, s = ±ωi. Therefore, the purely imaginary roots obtained from the original characteristic equation (2) should be the same as those calculated from the converted equation (8) using the Rekasius substitution. Notice that is the set consisting of all the purely imaginary roots of (2) andˆ is the set composed of the exhaustive imaginary roots of (8), one can claim that these two sets are identical.
We will later show a series of unique steps to reveal the setˆ , then according to Theorem 1, we can reveal the complete imaginary spectra of the original time-delay system, i.e., . For now, we revisit the converted characteristic equation (8). It is trivial to observe that the new form of the equivalent characteristic equation (8) is, in fact, a complex polynomial in terms of three variables, i.e., s = ±ωi, T 1 , and T 2 . If a solution ω ∈ R + exists for (8), it should make both the real and imaginary parts of this complex equation vanish at the same time, where the coefficients m j s and n j s are polynomials in T 1 with real coefficients in ω. In order for (10) and (11) to have a common T 2 solution, we deploy the Dixon resultant and discriminant approach to obtain the necessary and sufficient condition.

A. DIXON RESULTANT AND DISCRIMINANT
The Dixon resultant method is widely used to obtain the necessary and sufficient condition for a set of polynomial equations to share nontrivial common solutions [21], [30] in addition to peer methodologies, such as Sylveter [31], Macaulay [32], and Sparse [33] resultant. Dixon resultant is selected over the others due to its efficiency in computation. Next, we present fundamental formulations of the Dixon resultant theory for the complete logic flow. To simplify the notation, we denote the polynomials (10) and (11) with coefficients parameterized in ω and T 1 , We then formulate the Dixon resultant based on these two polynomials to get the necessary and sufficient condition for a common T 2 solution. Consider the following polynomial in T 2 where p k (α) , k = 1, 2 stands for uniformly replacing T 2 with a dummy variable α. The polynomial in (14) is called the Dixon polynomial. It is obvious that δ is symmetric with respect to α and T 2 , for the expression remains unchanged when exchanging α and T 2 , i.e., δ (T 2 , α) = δ (α, T 2 ). Besides, δ is of degree d max in α where d max = max deg (p 1 , T 2 ) , deg (p 2 , T 2 ) and the notation deg (p k , T 2 ) , k = 1, 2 stands for the highest degree of T 2 in polynomial p k . Due to the fact that each common zero T 2 of p 1 and p 2 will make δ (T 2 , α) zero regardless of α values, the coefficient of each power product of α in δ (T 2 , α) must be identically zero at this common T 2 root. This property leads to d max equations in T 2 based on the VOLUME 8, 2020 coefficients of α k (k = 0, 1, . . . , d max − 1). The coefficient matrix F (ω, T 1 ) ∈ R d max ×d max of these d max equations is called the Dixon matrix, For a nontrivial solution of (15) to exist, the Dixon matrix The above equation yields the necessary and sufficient condition on the coefficients of p 1 and p 2 for the existence of a common T 2 solution of (10) and (11). We formalize the conclusion by the following theorem, proof of which has just been discussed. Theorem 2 ( [34]): The necessary and sufficient condition for two polynomials p 1 (T 2 ) and p 2 (T 2 ) to share a nontrivial common zero is that the corresponding Dixon matrix The notation R T 2 represents the Dixon resultant obtained by eliminating T 2 . Notice that R T 2 is a function of ω and T 1 . If the range of ω values is known, one can scan it within this range and calculate the corresponding T 1 ∈ R values from R T 2 . With the obtained pairs (ω, T 1 ) ∈ R 2 , the remaining variable T 2 ∈ R can be calculated from either (10) or (11). Then the delay values which form the KOH could be reconstructed from (7) with feasible triplets (ω, T 1 , T 2 ) ∈ R 3 . The complete and exact knowledge of the frequency range of (1) is the foundation for these steps, which serves as the main contribution of the paper. To achieve this, we propose the novel methodologies in the next subsection.

B. DETERMINATION OF THE EXACT RANGE OF IMAGINARY SPECTRA
Due to the continuity and the differentiability of ω with respect to agent parameters [19], we start the process by checking the following extremum condition We will show later that the above condition leads to another equation in ω and T 1 . Notice that both (17) and this new equation have to be satisfied simultaneously, and we end up having two equations with two unknowns. Taking the differential of (17) leads to which can be converted to Considering that the general term ∂R T 1 ∂ω is not identically zero and the condition in (18)

holds, the following equation is valid
We now end up having two equations (17) and (21) with two unknowns ω and T 1 . Next, we bring up the following definition by applying the Dixon resultant procedure again on (17) and (21) to eliminate T 1 . Definition 4 ( [22]) (Discriminant of R T 2 (ω, T 1 )): The discriminant of R T 2 (ω, T 1 ) with respect to T 1 , denoted as D T 1 (ω), is the resultant formulated to eliminate T 1 by using two polynomial equations R T 2 (ω, T 1 ) = 0 and is a function of only one variable, i.e., ω. Solving its real roots yields candidates for the lower and upper positive bounds ofˆ for the case where the imaginary spectra are strictly positive, i.e., the lower bound is nonzero. Notice that special cases with zero imaginary spectrum is excluded from this paper. We will present the main contribution as the following theorem.

Theorem 3: The exact lower and upper bounds ofˆ for the time-delay system (1) are the minimum and maximum positive real roots of the discriminant D T 1 (ω) that correspond to
Proof: The imaginary spectraˆ can be found exhaustively by calculating the purely imaginary roots of the system's characteristic equation (2), which exhibits real and imaginary parts being (10) and (11), respectively. Notice that both these two equations are in terms of three variables, i.e., ω, T 1 , and T 2 with the last two being agent parameters. The objective to find the exact lower and upper bounds of is reduced to the detection of the minimum and maximum positive real values of ω corresponding to the common solution triplets (ω, T 1 , T 2 ) ∈ R 3 for (10) and (11). A common T 2 solution of these two equations enforces R T 2 in (17) being zero. This equation and the extremum condition of ω in (18) lead to the zero discriminant, the solution of which leads to the minimum and maximum positive real values of ω.

C. A NEW METHODOLOGY FOR THE GENERATION OF THE STABILITY MAP
As pointed out previously, the discriminant D T 1 (ω) is a function of only one variable, i.e., ω, and therefore its positive real roots can be calculated precisely. Denote the sequence of all the positive roots of D T 1 asω 1 <ω 2 < . . . <ω q−1 <ω q , where q ∈ Z + is a positive integer. For the calculation of the upper bound ofˆ , the following procedure is presented with the initial condition k = q, that is, the largest positive real root of D T 1 (ω) is used: (1) Plugω k into R T 2 ω k , T 1 = 0 and solve for T 1 ∈ R. If the real solution set is empty, reduce the counter k by 1 and repeat step (1). Otherwise, denote the T 1 ∈ R solution sequence as T 11 < T 12 < . . . < T 1(r−1) < T 1r , where r ∈ Z + . 36598 VOLUME 8, 2020 (2) For j = 1, 2, . . . , r − 1, r, plug the composition ω k , T 1j ∈ R 2 into (10) or (11) to solve for T 2 ∈ R. If no such T 2 ∈ R root exists, increase the counter j by 1 and restart step (2). Otherwise, declareω k as the upper bound ofˆ and denote it as ω before exit. If j reaches r and still no T 2 ∈ R exists for (10) or (11), reduce k by 1 and return to step (1).
(3) If k reaches 1 and still no T 2 ∈ R exists for (10), then the system is delay-independent stable or unstable, that is, the stability of the time-delay system is not affected by time delays.
For the calculation of the lower bound, denoted as ω, follow the same procedure from step (1) to step (3) except beginning with k = 1 and increasing it by 1 in steps (1) and (2) until the first (T 1 , T 2 ) ∈ R 2 are obtained. Label the currentω i as ω. We can then claim that the lower and upper bounds ofˆ are ω and ω, respectively. With this crucial knowledge, the following procedure is performed to reveal numerically the exact and complete stability map of the time-delay system (1): (a) Scan ω within the range ω, ω using a sufficiently small step size and solve R T 2 (ω, T 1 ) = 0 for T 1 ∈ R.
(b) For each (ω, T 1 ) ∈ R 2 pair obtained in (a), solve (10) or (11) for T 2 ∈ R. If such T 2 ∈ R root exists, go to the next step. Otherwise, return to step (a) with the next ω value to be swept.
The above steps (a) to (c) produce the exact set of ℘ 0 . This set, together with the generated ℘, constitute the complete KOH, i.e., the stability switching curves, in the domain of (τ 1 , τ 2 ) ∈ R 2+ . Then the stable and unstable regions in the domain of the delays can be determined effectively by the CTCR paradigm. We will then test the proposed approach over an example case in the next section.

IV. CASE STUDY
In order to demonstrate the effectiveness and strength of the proposed method, we take an example case from a previous published work [10] The corresponding characteristic equation is f (s, τ ) = s 2 − 0.4s + 6.48 + (−0.192s − 1.0788) e −τ 1 s + (0.384s + 0.1296) e −τ 2 s + 0.039024e −(τ 1 +τ 2 )s . (23) Note that for this particular example (22), we have rank (B 1 ) = rank (B 2 ) = 1, and rank (B 1 + B 2 ) = 2. We then suppress the steps involving the Dixon resultant as in (17), and discriminant as in Definition 4 for the sake of space constraints. Due to Theorem 3, this system's exact range of the imaginary spectra, i.e.,ˆ , is obtained as ω, ω = [2.1690, 2.9165], which is in accordance with the  numerical findings in the cited paper. We then follow the steps (a) to (c) in the last section to get the ℘ 0 and then the complete set of ℘ as per (4). It is noticeable that the union of ℘ 0 and ℘ is the KOH. Then the CTCR paradigm is implemented to reveal the exact stability map of the system as shown in Fig. 1. The obtained stability map is identical to that in [10]. For the purpose of further verification, a numerical algorithm called QPmR is deployed to calculate the dominant (rightmost) roots of (1). The real parts of the dominant roots are plotted with respect to the delay values τ 1 and τ 2 as in Fig. 2. The stable regions that fall below the Re [s dom ] = 0 plane precisely match those in Fig. 1, which further verifies the proposed methodology.

V. CONCLUSION
In this paper, a novel frequency-domain approach is proposed for the determination of the exact frequency range and the stability feature of the time-delay system (1). Different from previous works in the literature, this is the first attempt, to the best of the authors' knowledge, to extract the exact range of the imaginary spectra of the system of interest. By sweeping VOLUME 8, 2020 the obtained range, the complete set of the stability switching curves could be effectively determined, which facilitates the deployment of the CTCR paradigm for the stability analysis. For future work, we would like to extend the proposed methodology for systems with more delays and a higher level of complexity in dynamics, including but not limited to multi-agent systems [35] and switched systems [36], [37]. He was a Postdoctoral Fellow with the College of Automation Engineering, Nanjing University of Aeronautics and Astronautics, China. He is currently an Associate professor with the College of Electronics and Information Engineering, Suzhou University of Science and Technology, China, and a Visiting Scholar with the School of Electrical Engineering, Korea University, Seoul, South Korea. His research interests include observer design, model-based fault detection, and fault-tolerant control. VOLUME 8, 2020