# An iterative frequency-sweeping approach for stability analysis of linear systems with multiple delays

An iterative frequency-sweeping approach for stability analysis of linear systems with multiple... Abstract In this article, we study the stability of linear systems with multiple (incommensurate) delays, by extending a recently proposed frequency-sweeping approach. First, we consider the case where only one delay parameter is free while the others are fixed. The complete stability w.r.t. the free delay parameter can be systematically investigated by proving an appropriate invariance property. Next, we propose an iterative frequency-sweeping approach to study the stability under any given multiple delays. Moreover, we may effectively analyse the asymptotic behaviour of the critical imaginary roots (if any) w.r.t. each delay parameter, which provides a possibility for stabilizing the system through adjusting the delay parameters. The approach is simple (graphical test) and can be applied systematically to the stability analysis of linear systems including multiple delays. A deeper discussion on its implementation is also proposed. Finally, various numerical examples complete the presentation. 1. Introduction In this article, we consider the following linear time-delay system with multiple delay parameters   $$\label{MTDSsytem} \dot x(t) = Ax(t) + \sum\limits_{\ell = 1}^L {{B_\ell }} x(t - {\tau _\ell }),$$ (1.1) where $$A$$ and $$B_{\ell}$$ are constant matrices with compatible dimensions; $$\tau_{\ell} \ge 0$$, $$\ell = 1, \ldots ,L$$, are independent delay parameters ($$L$$ denotes the number of delay parameters). The delay combination may be expressed by a vector $$\overrightarrow \tau = ({\tau _1}, \ldots, {\tau _L})$$. The stability of time-delay system (1.1) has been largely investigated in the literature, see e.g., a survey paper Sipahi et al. (2011), or the books Michiels & Niculescu (2014) and Niculescu (2001), and the references therein. In the literature, in the framework of the so-called parameter-based approach, the stability problems for time-delay systems may be roughly divided in two categories: the $$\tau$$-decomposition problem (Lee & Hsu, 1969) and the $$D$$-decomposition problem (Neimark, 1949). For the former the delays are treated as free parameters, while for the latter some system/controller parameters are considered free (the delays are all fixed). For further discussions on such topics, we refer to Michiels & Niculescu (2014) and the references therein. The problem considered in this article belongs to the so-called $$\tau$$-decomposition problem. We start by recalling some of the existing results for the systems with single delay parameter. If $$L=1$$, the system (1.1) reduces to the form   $$\label{SingleDelaysytem} \dot x(t) = Ax(t) + Bx(t - \tau ),$$ (1.2) where $$B$$ is a constant matrix and $$\tau$$ is the delay parameter. The characteristic function for the system (1.2) is $$\det (\lambda I - A - B{e^{- \tau \lambda }})$$, where $$I$$ is the identity matrix of appropriate dimensions. If the delay parameters of system (1.1) are commensurate, the system can be expressed in the form   $$\label{CommensurateDelaysytem}\dot x(t) = Ax(t) + \sum\limits_{\ell = 1}^L {{B_\ell }} x(t - \ell \tau ),$$ (1.3) with the single delay parameter $$\tau$$. The characteristic function for system (1.3) writes as $$\det (\lambda I - A - \sum\limits_{\ell = 1}^L {{B_\ell }} {e^{ - \ell \tau \lambda }})$$. For the time-delay systems (1.2) and (1.3), the corresponding characteristic functions involve only one delay parameter (this is a major distinction with the multiple-delay system (1.1)). In such a case, the complete stability problem1 was systematically solved by the frequency-sweeping approach recently proposed in Li et al. (2015). The core result lies in that the invariance property concerning the asymptotic behaviour of the critical imaginary roots (CIRs) w.r.t. the infinitely many positive critical delays (CDs) was proved. It is worth mentioning that one of the main difficulties in analysing the complete stability lies in the characterization and the corresponding classification of the case with multiple characteristic roots on the imaginary axis. In the case where the delays are not commensurate, further discussions on characterizing multiple characteristic roots on the imaginary axis as well as related properties can be found in Boussaada & Niculescu (2016a,b). However, the stability problem for the multiple-delay system (1.1) is much more complicated and may have some unexpected dynamic behaviour (for instance, the delay interference phenomenon, see Michiels & Niculescu, 2007). Roughly speaking, to the best of the authors’ knowledge, the existing studies for the multiple-delay system (1.1) can be categorized into two classes: (i) One class of studies focus on finding the stability crossing set (SCS) in the delay parameter space.2 For a class of time-delay systems without cross-terms in the characteristic functions, the SCSs for the case $$L=2$$ and for the case $$L=3$$ are reported in Gu et al. (2005) and Gu & Naghnaeian (2011), respectively. Moreover, the geometric structure of the corresponding two- and three-dimensional SCSs are studied in detail. An extension to the two-delay case with one cross-term in the characteristic function was reported later in Naghnaeian & Gu (2013). By using some different arguments (based on the properties of the Rekasius transformation), another important series of results have been reported in, e.g., Sipahi & Delice (2009), Sipahi & Delice (2011) and Sipahi & Olgac (2006) on obtaining the SCSs (also called potential stability switching curves in these references). The methods are now applicable to general linear time-delay systems. In our opinion, the main advantage of the above class of studies is that the corresponding SCS in the case $$L=2$$ or $$3$$ may be visualized by a two- or three-dimensional figure, from which we may intuitively analyse the stability. If $$L > 3$$, a visualization of the SCS is difficult. In some cases (coupling small with large delays), some visualization was also proposed in the case of four delays describing immune dynamics in leukaemia (see e.g., Niculescu et al., 2010). (ii) Another class of methods concentrate on counting the number of unstable roots under a given delay vector $$\overrightarrow \tau$$. The idea lies in applying the argument principle to the corresponding characteristic function. An advantage of this class of studies is that the number of the delays $$L$$ is allowed to be any large. A representative conclusion can be found in Stépán (1979, 1989), known as the Stépán’s formula. Other interesting results in the same spirit can also be found in, e.g., Hassard (1997), Hu & Liu (2007) and Vyhlídal & Zítek (2009). However, when the system has CIRs, it is difficult to adopt the above two classes of methodologies to analyse the asymptotic behaviour of the CIRs. The difficulty is mainly related to the treatment of the case with multiple and/or degenerate CIRs. The asymptotic behaviour analysis is of practical and theoretical importance, especially when $$L$$ is large. When CIRs appear for a given $$\overrightarrow \tau$$, we usually want to know if there is a delay vector near $$\overrightarrow \tau$$ such that the system is asymptotically stable. If so, we will further consider how to find such a delay vector. From the above discussions, we see that there is still much room for the stability analysis of linear systems with multiple (incommensurate) delays, which motivates the work of this article. A straightforward idea is to extend the mathematical results for asymptotic behaviour analysis of the single-delay problem (proposed in Li et al., 2015) to the multiple-delay problem. However, to the best of the authors’ knowledge, such an extension is rather difficult as the asymptotic behaviour analysis for multiple-parameter problems has remained open in mathematics.3 In this article, we will propose an ‘indirect’ yet effective approach, called the iterative frequency-sweeping approach. The core of this approach is an important property to be proved in this article: If the multiple-delay system (1.1) has only one free delay parameter (the other ones are fixed), the effect of a CIR’s asymptotic behaviour on the stability w.r.t. its infinitely many positive CDs is invariant. This invariance property is generalized from the one confirmed in Li et al. (2015) for the systems with commensurate delays. To the best of the authors’ knowledge, such a result was not proposed in the open literature and it gives insights in understanding the way multiple delays may affect the stability property. As mentioned in the single-delay case, one of the main difficulties of the complete characterization of the stability is related to the asymptotic behaviour analysis of the multiple and/or degenerate CIRs. The invariance property is essential for overcoming this difficulty. With the invariance property, we may study the complete stability w.r.t. any delay parameter $${\tau _\ell }$$ through some simple and quite effective frequency-sweeping test. Furthermore, for any given delay combination $$\overrightarrow {{\tau ^ \#}}$$, we may accurately compute the number of unstable roots by using $$L$$ times frequency-sweeping tests in an appropriate iterative manner. Moreover, if the multiple-delay system (1.1) has CIRs, we may easily analyse the asymptotic behaviour of the CIRs w.r.t. each delay parameter from the frequency-sweeping curves (FSCs). As a consequence, we know if a stabilizing delay combination exists near $$\overrightarrow {{\tau ^ \#}}$$ and how to find it (if it exists). Finally, we propose an algorithm for a class of multiple-delay systems such that the stability may be automatically analysed by a single computer program. In our opinion, such an algorithm is easy to understand, to implement, and to apply, as illustrated by some examples taken from the open literature. The rest part of article is organized as follows. Some preliminaries are reviewed in Section 2. In Section 3, the complete stability of a linear time-delay system including multiple delays with single free delay is analysed. In Section 4, an approach is presented for computing the number of unstable roots for a linear system with any combination of multiple delays. Numerical examples are given in Section 5. In Section 6, an algorithm for automatic stability analysis for a class of multiple-delay systems is proposed. Finally, in Section 7, some concluding remarks end the article. Notations: In the sequel, $$\mathbb{R}$$ ($$\mathbb{R}_+$$) denotes the set of (positive) real numbers and $$\mathbb{C}$$ is the set of complex numbers; $$\mathbb{C}_-$$ and $$\mathbb{C}_+$$ denote respectively the left half-plane and right half-plane in $$\mathbb{C}$$; $$\mathbb{C}_0$$ is the imaginary axis and $$\partial \mathbb{D}$$ is the unit circle in $$\mathbb{C}$$; $$\mathbb{Z}$$, $$\mathbb{N}$$, and $$\mathbb{N}_+$$ are the sets of integers, non-negative integers and positive integers, respectively. $$\varepsilon$$ is a sufficiently small positive real number. Next, $$I$$ is the identity matrix and $$\overrightarrow 0 = (0, \ldots ,0)$$ of appropriate dimensions. For $$\gamma \in \mathbb{R}$$, $$\lceil \gamma \rceil$$ denotes the smallest integer greater than or equal to $$\gamma$$. Finally, $$\det (\cdot)$$ denotes the determinant of its argument. 2. Preliminaries The characteristic function of time-system (1.1) is   $f(\lambda ,\overrightarrow \tau ) = \det \left(\lambda I - A - \sum\limits_{\ell = 1}^L {{B_\ell }{e^{ - {\tau _\ell }\lambda }}} \right)\!,$ which is a quasipolynomial. For a non-zero delay vector $$\overrightarrow \tau$$, the characteristic equation $$f(\lambda ,\overrightarrow \tau ) = 0$$ has an infinite number of roots, i.e., the time-system (1.1) has an infinite number of characteristic roots. Time-delay system (1.1) is asymptotically stable if and only if all the characteristic roots lie in the open left half-plane $$\mathbb{C}_-$$, see e.g., Hale & Verduyn Lunel (1993) and Michiels & Niculescu (2014) for more comprehensive explanation. Throughout this article, we denote by $$NU(\overrightarrow \tau)$$ the number of characteristic roots in the open right half-plane $$\mathbb{C}_+$$ for system (1.1) in the presence of a delay vector $$\overrightarrow \tau$$. Clearly, the time-delay system (1.1) is asymptotically stable if and only if $$NU(\overrightarrow \tau)=0$$ and the system has no CIRs. If at $$\overrightarrow {{\tau ^*}} = (\tau _1^*, \ldots ,\tau _L^*)$$, $$f(j\omega^* , \overrightarrow {\tau^*}) = 0$$ ($${\omega ^*} \in \mathbb{R}_+$$, $$j$$ is the imaginary unit, and such a pair $$(\,j\omega^* , \overrightarrow {\tau^*})$$ is called a critical pair), then for all $$\overrightarrow \tau = \overrightarrow {{\tau ^*}} + ({k_1}, \ldots ,{k_L})\frac{{2\pi }}{{{\omega ^*}}}$$ ($${k_\ell } \in \mathbb{Z}$$ such that $$\tau _\ell ^* + {k_\ell }\frac{{2\pi }}{{{\omega ^*}}} \ge 0$$, $$\ell = 1, \ldots ,L$$), $$f(\,j\omega^*, \overrightarrow \tau) = 0$$. That is to say, a CIR $$j\omega^*$$ repeats infinitely many times as each $$\tau_{\ell}$$ increases, with a periodicity $$\frac{{2\pi }}{{\omega^*}}$$. Remark 2.1 It is a common assumption that $$\lambda=0$$ is not a characteristic root, otherwise the system can not be asymptotically stable for any $$\overrightarrow \tau$$. Owing to the conjugate symmetry of the spectrum, it suffices to consider only the CIRs with nonnegative imaginary parts. As mentioned in Section 1, the asymptotic behaviour analysis w.r.t. multiple free delay parameters is rather complicated (corresponding to an open mathematical problem). We will propose an iterative procedure: Analysing the asymptotic behaviour w.r.t. one free delay each time. 3. Complete stability w.r.t. one delay parameter In this section, we study the case where $$L-1$$ delay parameters among $$\tau_1$$, $$\ldots$$, $$\tau_L$$ are fixed and the remaining one is ‘free’. The objective is to study the complete stability of time-delay system (1) w.r.t. the remaining free delay parameter. For simplicity, we adopt the following notation $$\delta (\ell ) = ({\delta _1}(\ell ), \ldots ,{\delta _L}(\ell ))$$ with   ${\delta _i}(\ell ) = \left\{ \begin{array}{l} 0, {\kern 1pt} {\kern 1pt} {\rm{if}}{\kern 1pt} {\kern 1pt} i \ne \ell , \\ 1,{\kern 1pt} {\kern 1pt} {\rm{if}}{\kern 1pt} {\kern 1pt} i = \ell . \\ \end{array} \right.$ Then, $$\overrightarrow \tau = \sum\limits_{\ell = 1}^L {{\tau _\ell }\delta (\ell )}$$. Suppose $$\tau_{\chi}$$ ($$\chi \in \{ 1, \ldots ,L\}$$) is the free parameter and the other $$L-1$$ fixed delay parameters are $$\tau _\ell^ \#$$ ($$\ell \ne \chi$$). In this scenario, $$\overrightarrow \tau$$ can be expressed as $$\overrightarrow \tau = {\tau _\chi }\delta (\chi) + {F_\chi}$$, with $${F_\chi } = \sum\limits_{\ell \ne \chi} {{\tau _\ell^ \#}\delta (\ell)}$$. For instance, in the case $$L=4$$, suppose $${\tau _1} = \tau _1^ \#$$, $${\tau _3} = \tau _3^ \#$$, and $${\tau _4} = \tau _4^ \#$$ are fixed, while $$\tau_2$$ is the free delay parameter. We may express $$\overrightarrow \tau = (\tau _1^ \#, {\tau _2}, \tau _3^ \#, \tau _4^ \#)$$ as $$\overrightarrow \tau = {\tau _2}\delta (2) + {F_2}$$ with $${F_2} = (\tau _1^ \#, 0, \tau _3^ \#, \tau _4^ \#)$$ and $$\delta (2) = (0,1,0,0)$$. We study the stability of time-delay system (1.1) w.r.t. $$\tau_{\chi}$$ along the whole non-negative $$\tau_{\chi}$$-axis, i.e., the complete stability problem w.r.t. the delay parameter $$\tau_{\chi}$$. More precisely, we will keep track of the number of unstable roots, denoted by $$N{U_{{F_\chi}}}({\tau _\chi})$$, for $${\tau _\chi} \in [0, + \infty )$$. In this context, the characteristic function $$f(\lambda ,\overrightarrow \tau )$$ can be rewritten as   \begin{eqnarray}\label{ftauFsecondform} f(\lambda ,{\tau _\chi},{F_\chi}) = \sum\limits_{i = 0}^{{q_\chi}} {{{\widetilde a}_{\chi i}}(\lambda ){e^{ - i{\tau _\chi}\lambda }}}, \end{eqnarray} (3.1) where $$\widetilde a_{\chi i}$$ ($$i =0, \ldots ,{q_\chi}$$) are polynomials in $$\lambda$$ and $${e^{ - \tau _\ell ^\# \lambda }}$$ ($$\ell \in \{ 1, \ldots ,L\}$$, $$\ell \ne \chi$$) with real coefficients. Furthermore, the expression (3.1) can be viewed as a polynomial of $${e^{ - {\tau _\chi}\lambda }}$$ where $${\widetilde a_{\chi 0}}(\lambda ), \ldots ,{\widetilde a_{\chi {q_\chi}}}(\lambda )$$ are treated as the coefficient functions. We introduce the following two-variable polynomial expression of $$f(\lambda ,{\tau _\chi},{F_\chi})$$:   \begin{eqnarray}\label{ptauF} p(\lambda ,{z_\chi},{F_\chi}) = \sum\limits_{i = 0}^{{q_\chi}} {{{\widetilde a}_{\chi i}}(\lambda )z_{^\chi}^i}, \end{eqnarray} (3.2) where   ${z_\chi} = {e^{ - {\tau _\chi}\lambda }}.$ Remark 3.1 The coefficient functions $${\widetilde a_{\chi i}}(\lambda)$$ ($$i=1, \ldots, q_{\chi}$$) of (3.1) as well as (3.2) are independent of $$\tau_{\chi}$$. The detection of the CIRs and the CDs for $$f(\lambda ,{\tau _\chi},{F_\chi}) = 0$$ amounts in detecting the critical pairs $$(\lambda, z_\chi)$$ ($$\lambda \in \mathbb{C}_0$$ and $$z_\chi \in \partial \mathbb{D}$$) for $$p(\lambda ,{z_\chi},{F_\chi})=0$$. Without any loss of generality, suppose that there exist $$u$$ critical pairs for $$p(\lambda, z_\chi, F_\chi)=0$$ denoted by $$({\lambda _{\chi,0}} = j{\omega _{\chi,0}},{z_{\chi,0}})$$, $$({\lambda _{\chi,1}} = j{\omega _{\chi,1}},{z_{\chi,1}})$$, $$\ldots$$, $$({\lambda _{\chi,u-1}} = j{\omega _{\chi,u-1}},{z_{\chi,u-1}})$$ with $$0 < {\omega _{\chi,0}} \le {\omega _{\chi,1}} \le \cdots \le {\omega _{\chi, u - 1}}$$. Once all the critical pairs $$({\lambda _{\chi,\alpha} },{z_{\chi,\alpha} }),\alpha = 0, \ldots ,u - 1$$, are found, all the critical pairs $$(\lambda,\tau_\chi)$$ for $$f(\lambda ,{\tau _\chi},{F_\chi}) = 0$$ can be obtained: For each CIR $$\lambda_ {\chi, \alpha}$$, the corresponding (infinitely many) CDs are given by $${\tau _{\chi, \alpha, k}} \buildrel {\it{\Delta}} \over = {\tau _{\chi, \alpha, 0}} + \frac{{2k\pi }}{{{\omega _{\chi, \alpha} }}}, k \in \mathbb{N}, {\tau _{\chi, \alpha, 0}} \buildrel {\it{\Delta}} \over = \min \{ \tau \ge 0:{e^{ - \tau {\lambda _{\chi, \alpha} }}} = {z_{\chi, \alpha} }\}$$. The pairs $$({\lambda _{\chi, \alpha} },{\tau _{\chi, \alpha ,k}}),k \in \mathbb{N}$$, define a set of critical pairs associated with $$({\lambda _{\chi, \alpha} },{z_{\chi, \alpha} })$$. Remark 3.2 According to the root continuity argument, if the time-delay system (1.1) has no CIRs for all $${\tau _\chi } \ge 0$$, $$NU({\tau _\chi}\delta (\chi) + {F_\chi}) = NU({F_\chi})$$ for all $$\tau_{\chi}>0$$. 3.1. Asymptotic behaviour of a critical pair An essential step for the stability analysis is to address the asymptotic behaviour of a critical pair $$({\lambda _{\chi, \alpha }}, {\tau _{\chi, \alpha, k}})$$. Since the CIR may be multiple and split as $$\tau_{\chi}$$ increases near a positive CD, we now introduce a general notation to describe the asymptotic behaviour for a critical pair $$({\lambda _{\chi,\alpha }} = j{\omega _{\chi,\alpha }},{\tau _{\chi,\alpha ,k}})$$ (provided $$F_{\chi}$$ is given)   $\begin{array}{l} {\it{\Delta}} NU_{{\lambda _{\chi,\alpha }}}^{{F_\chi}}({\tau _{\chi,\alpha ,k}}) \buildrel {\it{\Delta}} \over = NU(({\tau _{\chi,\alpha ,k}} + \varepsilon )\delta (\chi) + {F_\chi }) - NU(({\tau _{\chi,\alpha ,k}} - \varepsilon )\delta (\chi) + {F_\chi}), \end{array}$ which stands for the change in $$NU(\overrightarrow \tau)$$ caused by the variation of the CIR $${\lambda _{\chi ,\alpha }}$$ as $$\tau_{\chi}$$ increases from $${\tau _{\chi,\alpha ,k}} - \varepsilon$$ to $${\tau _{\chi,\alpha ,k}} + \varepsilon$$, i.e., the asymptotic behaviour of the CIR $${\lambda _{\chi,\alpha }}$$ at a positive CD $${\tau _{\chi,\alpha ,k}}$$. The value of $${\it{\Delta}} NU_{{\lambda _{\chi,\alpha }}}^{{F_\chi }}({\tau _{\chi,\alpha ,k}})$$ can be calculated by invoking the associated Puiseux series. One may refer to Chapter 4 of Li et al. (2015) for a general algorithm for invoking the Puiseux series. At this stage, it is interesting to see that the work of invoking the Puiseux series may be replaced by a graphical (frequency-sweeping) test, which will be presented later in this article. In order to analyse the complete stability, we need to specifically address the situation as $$\tau_{\chi}$$ increases from $$0$$ to $$+\varepsilon$$. In other words, we need to know $$NU({F_\chi} + \varepsilon \delta (\chi))$$. Similarly to Theorem 5.1 in Li et al. (2015), we have the following: Theorem 3.1 If the time-delay system (1.1) has no CIRs when $$\tau_{\chi}=0$$, $$NU({F_\chi} + \varepsilon \delta (\chi))= NU({F_\chi})$$. Otherwise, $$NU({F_\chi} + \varepsilon \delta (\chi)) - NU({F_\chi})$$ equals to the number of the values in $$\mathbb{C}_+$$ of the Puiseux series for all the corresponding CIRs when $$\tau_{\chi}=0$$ with $${\it{\Delta}} \tau_{\chi} = + \varepsilon$$. Proof. If $$F_\chi$$ is zero vector, the system has finitely many characteristic roots when $$\tau_{\chi}=0$$. As $$\tau_{\chi}$$ increases to a sufficiently small positive number $$+\varepsilon$$, infinitely many new characteristic roots appear at the far left of the complex plane while the original (finitely many) characteristic roots change continuously w.r.t. $$\tau_{\chi}$$. If $$F_\chi$$ is a non-zero vector, all the (infinitely many) characteristic roots change continuously as $$\tau_{\chi}$$ increases from $$0$$ to $$+\varepsilon$$. Thus, if the system has no CIRs when $$\tau_{\chi}=0$$, it is trivial to know that $$NU({F_\chi} + \varepsilon \delta (\chi)) - NU({F_\chi}) = 0$$. Otherwise, the value of $$NU({F_\chi} + \varepsilon \delta (\chi)) - NU({F_\chi})$$ is determined by the asymptotic behaviour of the corresponding CIRs, which can be explicitly analysed through the Puiseux series (as stated in Theorem 3.1). □ 3.2. Frequency-sweeping curves (FSCs) Under a given $$F_{\chi}$$, the FSCs for time-delay system (1.1) are obtained by the following procedure. Frequency-sweeping curves (FSCs): Sweep $$\omega \ge 0$$ and for each $$\lambda = j\omega$$ we have $$q_{\chi}$$ solutions of $$z_{\chi}$$ such that $$p(\lambda ,{z_\chi},{F_\chi}) = 0$$ (denoted by $${z_{\chi,1}}(\,j\omega ), \ldots ,{z_{\chi,{q_\chi}}}(\,j\omega )$$). In this way, we obtain $$q_{\chi}$$ FSCs $${{\it{\Gamma}} _{\chi,i}}(\omega ):$$$$\left| {{z_{\chi,i}}(\,j\omega )} \right|$$ vs. $$\omega$$, $$i = 1, \ldots ,{q_\chi}$$. For simplicity, we denote by $$\Im _1$$ the line parallel to the abscissa axis with ordinate equal to 1. If $$({\lambda _{\chi,\alpha }} = j{\omega _{\chi,\alpha }},{\tau _{\chi,\alpha ,k}})$$ are a set of critical pairs, some FSCs intersect $${\Im _1}$$ at $$\omega = {\omega _{\chi, \alpha}}$$. It is straightforward to see that all the CIRs and the CDs can be detected by the FSCs (more precisely, the intersection of the FSCs and the line $$\Im _1$$). We now address the asymptotic behaviour of the FSCs (such an idea was proposed and largely discussed in Li et al., 2015). For a set of critical pairs $$({\lambda _{\chi,\alpha} =j \omega_{\chi,\alpha}},{\tau _{\chi, \alpha ,k}})$$ (as usually assumed, $$\lambda _{\chi,\alpha} \ne 0$$), there must exist some FSCs such that $${z_{\chi,i}}(\,j{\omega _{\chi,\alpha }}) = {z_{\chi ,\alpha }} = {e^{ - {\tau _{\chi,\alpha ,0}}{\lambda _{\chi,\alpha }}}}$$ intersecting $${\Im _1}$$ when $$\omega = {\omega _{\chi ,\alpha }}$$. Among these FSCs, we denote the number of those above $$\Im _1$$ when $$\omega = {\omega _{\chi,\alpha }} + \varepsilon$$ ($$\omega = {\omega _{\chi,\alpha }} - \varepsilon$$) by $$NF_{{z_{\chi,\alpha }}}^{{F_\chi}}({\omega _{\chi,\alpha }} + \varepsilon )$$ ($$NF_{{z_{\chi,\alpha }}}^{{F_\chi}}({\omega _{\chi ,\alpha }} - \varepsilon )$$). We now introduce the notation $${\it{\Delta}} NF_{{z_{\chi,\alpha }}}^{{F_\chi}}({\omega _{\chi,\alpha }})$$ as defined below   ${\it{\Delta}} NF_{{z_{\chi,\alpha }}}^{{F_\chi}}({\omega _{\chi,\alpha }}) \buildrel {\it{\Delta}} \over = NF_{{z_{\chi,\alpha }}}^{{F_\chi}}({\omega _{\chi,\alpha }} + \varepsilon ) - NF_{{z_{\chi,\alpha }}}^{{F_\chi}}({\omega _{\chi,\alpha }} - \varepsilon ).$ Remark 3.3 It is a useful property that $${\it{\Delta}} NF_{{z_{\chi,\alpha }}}^{{F_\chi}}({\omega _{\chi ,\alpha }})$$ is invariant w.r.t. different CDs. Moreover, the value of $${\it{\Delta}} NF_{{z_{\chi,\alpha }}}^{{F_\chi}}({\omega _{\chi ,\alpha }})$$ may be easily obtained graphically (by observing how the corresponding FSCs intersect $$\Im _1$$): it equals to the number change of the corresponding FSCs above $${\Im _1}$$ when $$\omega = {\omega _{\chi,\alpha }} + \varepsilon$$ and $$\omega = {\omega _{\chi,\alpha }} - \varepsilon$$. For instance, for a set of critical pairs $$({\lambda _{\chi,\alpha} =j \omega_{\chi,\alpha}},{\tau _{\chi, \alpha ,k}})$$, suppose that there exists one and only one FSC such that $${z_{\chi,i}}(\,j{\omega _{\chi,\alpha }}) = {z_{\chi ,\alpha }} = {e^{ - {\tau _{\chi,\alpha ,0}}{\lambda _{\chi,\alpha }}}}$$ and that as $$\omega$$ increases near $${\omega _{\chi ,\alpha }}$$ the FSC intersects $${\Im _1}$$ from above to below. Then, $$NF_{{z_{\chi,\alpha }}}^{{F_\chi}}({\omega _{\chi,\alpha }} + \varepsilon ) = 0$$, $$NF_{{z_{\chi,\alpha }}}^{{F_\chi}}({\omega _{\chi,\alpha }} - \varepsilon )=1$$ and hence $${\it{\Delta}} NF_{{z_{\chi,\alpha }}}^{{F_\chi}}({\omega _{\chi,\alpha }}) \buildrel {\it{\Delta}} \over = NF_{{z_{\chi,\alpha }}}^{{F_\chi}}({\omega _{\chi,\alpha }} + \varepsilon ) - NF_{{z_{\chi,\alpha }}}^{{F_\chi}}({\omega _{\chi,\alpha }} - \varepsilon ) = 0 - 1 =-1$$. This is the case to be seen in Example 5.1. 3.3. Invariance property We are now in a position to present the result concerning the invariance property of the asymptotic behaviour w.r.t. the free delay parameter $$\tau_\chi$$ for time-delay system (1.1). Theorem 3.2 For a CIR $${\lambda _{\chi,\alpha }}$$ of time-delay system (1.1) with given $${F_\chi}$$, $${\it{\Delta}} NU_{{\lambda _{\chi,\alpha }}}^{{F_\chi }}({\tau _{\chi,\alpha ,k}})$$ is a constant $${\it{\Delta}} NF_{{z_{\chi ,\alpha }}}^{{F_\chi}}({\omega _{\chi,\alpha }})$$ for all $${\tau _{\chi,\alpha, k}} > 0$$. The proof is given in the Appendix. 3.4. Ultimate stability With the above invariance property, we now address the ultimate stability issue (i.e., the system stability as $${\tau _\chi} \to \infty$$). Property 3.1 For all $$i= 1, \dots, q_\chi$$, it follows that $$\left| {{z_{\chi, i}}(\,j\omega )} \right| \to \infty$$ as $$\omega \to \infty$$. Proof. As the time-system (1.1) is of retarded type and $$\left| {{e^\iota }} \right| = 1$$ for any purely imaginary number $$\iota$$, we have   $\mathop {\lim }\limits_{\omega \to \infty } \frac{{{{\widetilde a}_{\chi i}}(\,j\omega )}}{{{{\widetilde a}_{\chi 0}}(\,j\omega )}} = 0,i = 1, \ldots ,{q_\chi }.$ The proof is now complete according to (3.2).□ A critical frequency $${\omega _{\chi,\alpha }}$$ is called a crossing (touching) frequency for an FSC $${\it{\Gamma}} _{\chi,i} (\omega)$$, if $${\it{\Gamma}} _{\chi,i} (\omega)$$ crosses (touches without crossing) $${\Im _1}$$ as $$\omega$$ increases near $$\omega = \omega_{\chi,\alpha}$$. With the notions above, we have the following results. Theorem 3.3 $${{\it{\Gamma}} _{\chi ,i}}(\omega ),i = 1, \ldots ,{q_\chi }$$, have a crossing frequency, there exists a $$\tau_{\chi}^*$$ such that the time-delay system (1.1) is unstable for all $$\tau_{\chi} > \tau_{\chi}^ *$$ and $$\mathop {\lim }\limits_{{\tau _\chi } \to \infty } NU({\tau _\chi }\delta (\chi ) + {F_\chi }) = \infty$$. Theorem 3.4 A time-delay system (1.1) with given $${F_\chi}$$ must belong to the following three types: Type 1: The system has a crossing frequency and $$\mathop {\lim }\limits_{{\tau _\chi} \to \infty } NU({\tau _\chi}\delta (\chi) + {F_\chi}) = \infty$$; Type 2: The system has neither crossing frequencies nor touching frequencies and $$NU({\tau _\chi}\delta (\chi) + {F_\chi}) = NU({F_\chi})$$ for all $$\tau _ \chi >0$$; Type 3: The system has touching frequencies but no crossing frequencies and $$NU({\tau _\chi}\delta (\chi) + {F_\chi})$$ is a constant for all $$\tau_\chi$$ other than the CDs. Following the line of Theorems 9.1 and 9.2 in Li et al. (2015), we may prove Theorems 3.3 and 3.4 in light of Property 3.1. 3.5. Explicit computation of number of unstable roots Assume that $$NU({F_\chi})$$ is known (this value can be obtained by the procedure to be given in the next section). We now show that $$NU({\tau _\chi}\delta (\chi) + {F_\chi})$$ can be expressed as an explicit function of $$\tau_{\chi}$$. For each CIR $${\lambda _{\chi,\alpha }}$$, we may choose any positive CD $${\tau _{\chi,\alpha ,k}}$$ to compute $${\it{\Delta}} NU_{{\lambda _{\chi ,\alpha }}}^{{F_\chi}}({\tau _{\chi,\alpha ,k}})$$ (the value is denoted by $${U_{{\lambda _{\chi,\alpha }}}}$$), through invoking the Puiseux series. Alternatively, we may directly have that $${U_{{\lambda _{\chi,\alpha }}}} = {\it{\Delta}} NF_{{z_{\chi,\alpha }}}^{{F_\chi}}({\omega _{\chi,\alpha }})$$ from the FSCs, according to Theorem 3.2. In the light of the invariance property, the explicit expression of $$NU({\tau _\chi}\delta (\chi) + {F_\chi})$$ is as follows: Theorem 3.5 For time-delay system (1.1) and for any $$\tau_{\chi}$$ which is not a CD, it follows that   $NU({\tau _\chi}\delta (\chi) + {F_\chi}) = NU(\varepsilon \delta (\chi) + {F_\chi}) + \sum\limits_{\alpha = 0}^{u - 1} {N{U_{\chi ,\alpha }}({\tau _\chi})} ,$   ${\rm where}\ N{U_{\chi,\alpha }}({\tau _\chi}) = \left\{ {\begin{array}{*{20}{l}} {0,\tau < {\tau _{\chi,\alpha ,0}},}\\ {2{U_{{\lambda _{\chi,\alpha }}}}\left\lceil {\frac{{{\tau _\chi} - {\tau _{\chi,\alpha ,0}}}}{{2\pi /{\omega _{\chi,\alpha }}}}} \right\rceil ,\tau > {\tau _{\chi,\alpha ,0}},} \end{array}} \right.{\rm{if}}\,{\tau _{\chi,\alpha ,0}} \ne 0,$   $\left\{ {\begin{array}{*{20}{l}} {0,\tau < {\tau _{\chi,\alpha ,1}},}\\ {2{U_{{\lambda _{\chi,\alpha }}}}\left\lceil {\frac{{{\tau _\chi} - {\tau _{\chi,\alpha ,1}}}}{{2\pi /{\omega _{\chi,\alpha }}}}} \right\rceil ,\tau > {\tau _{\chi,\alpha ,1}},} \end{array}} \right.{\rm{if}}\,{\tau _{\chi,\alpha ,0}} = 0.$ We can now systematically solve the complete stability w.r.t. the free delay parameter $$\tau_\chi$$. The system is asymptotically stable if and only if $$\tau_\chi$$ lies in the interval(s) with $$NU({\tau _\chi}\delta (\chi) + {F_\chi}) =0$$ excluding the CDs. The ultimate stability is known according to Theorems 3.3 and 3.4. 4. Stability analysis for any given delay vector We now study the stability of time-delay system (1.1) for any given delay vector $$\overrightarrow {{\tau ^ \#}} = (\tau _1^ \#, \ldots ,\tau _L^ \#)$$. Based on the results of Section 3, an iterative frequency-sweeping approach is given below. An Iterative Frequency-Sweeping Approach Step 0: Set $$\chi=1$$ and $${F_1} = (0, \ldots ,0)$$. Step 1: Compute $$NU({F_\chi } + \varepsilon \delta (\chi ))$$ by Theorem 3.1. Step 2: Generate the FSCs corresponding to $$p(\lambda ,{\tau _\chi},{F_\chi}) = 0$$, denoted by $${{\it{\Gamma}} _{\chi ,i}}(\omega ),i = 1, \ldots ,{q_\chi }$$. Following the results in Section 3, we can know $$NU({F_\chi } + \tau _\chi ^ \# \delta (\chi ))$$ and detect the CIRs (if any!). Step 3: If $$\chi < L$$, let $$\chi = \chi + 1$$ and $${F_\chi } = {F_{\chi - 1}} + \tau _{\chi - 1}^ \# \delta (\chi - 1)$$. Return to Step 1. (After running Steps 1–3 for the last time (i.e., when $$\chi = L$$), we know the value of $$NU(\overrightarrow {{\tau ^ \#}} )$$ and whether the system (1.1) has CIRs when $$\overrightarrow \tau = \overrightarrow {{\tau ^ \#}}$$). Step 4: If the system has no CIRs when $$\overrightarrow \tau = \overrightarrow {{\tau ^ \#}}$$, skip to Step 5. Otherwise, generate the FSCs corresponding to $$p(\lambda ,{\tau _\ell},{F_\ell ^ \#}) = 0$$ with $${F_\ell ^ \#} = \sum\limits_{\kappa \ne \ell} {{\tau _\kappa^ \#}\delta (\kappa)}$$, denoted by $${\it{\Gamma}} _{\ell,i}^ \# (\omega ),i = 1', \ldots ,{q_\ell}$$, $$\ell = 1, \ldots ,L$$.4 (We may analyse the asymptotic behaviour of the CIRs when $$\overrightarrow \tau = \overrightarrow {{\tau ^ \#}}$$ w.r.t. each delay parameter $$\tau_{\ell}$$, $$\ell=1,\ldots,L$$, from the FSCs $${\it{\Gamma}} _{\ell,i}^ \# (\omega ),i = 1, \ldots ,{q_\ell}$$. To be more precise, we can determine if increasing or decreasing a delay parameter $$\tau_{\ell}$$ sufficiently near $$\tau _\ell ^ \#$$ may stabilize the system, according to Theorem 3.2. As a result, we know if there exists a $$\overrightarrow \tau$$ in a small neighbourhood of $$\overrightarrow {{\tau ^ \#}}$$, at which the system is asymptotically stable.) Step 5: The procedure ends. Proposition 4.1 For time-delay system (1.1) and for any given delay vector $$\overrightarrow {{\tau ^ \#}} = (\tau _1^ \#, \ldots ,\tau _L^ \#)$$, the iterative frequency-sweeping approach allows to calculate $$NU(\overrightarrow \tau )$$ with $$\overrightarrow \tau = \overrightarrow {{\tau ^ \#}}$$. Proof. It is worth to mention that the invariance property guarantees that the stability characterization is complete w.r.t. any free delay parameter. After the first iteration (Steps 1–3 for the first time), we know the stability property for $$\tau = (\tau _1^\# ,0, \ldots ,0)$$. Next, after the second iteration (Steps 1–3 for the second time), we know the stability property for $$\tau = (\tau _1^\# ,\tau _2^\# ,0, \ldots ,0)$$. In this iterative manner, we are able to examine the stability properties for the delay vectors $$(\tau _1^\# ,0, \ldots ,0)$$, $$(\tau _1^\# ,\tau _2^\# ,0, \ldots ,0)$$, $$(\tau _1^\# ,\tau _2^\# ,\tau _3^\# ,0, \ldots ,0)$$, $$\ldots$$ The iterative frequency-sweeping approach requires totally $$L$$ iterations (one per independent delay parameter). Finally, as expected, after the last iteration (i.e., Steps 1–3 for the $$L$$-th time), we can determine the stability property for $$\tau = \overrightarrow {{\tau ^ \#}} = (\tau _1^ \#, \ldots ,\tau _L^ \#)$$.□ Remark 4.1 The above asymptotic behaviour analysis (by means of Step 4) is simple and effective. The effect of each delay parameter’s asymptotic behaviour can be explicitly known. It should be pointed out that the asymptotic behaviour analysis w.r.t. all the $$L$$ delay parameters (simultaneously) is very difficult (theoretically, it corresponds to some Puiseux series for $${\it{\Delta}} \lambda$$ w.r.t. $${\it{\Delta}} {\tau _1}, \ldots ,{\it{\Delta}} {\tau _L}$$), especially when $$L$$ is large. However, as mentioned in Section 1, to the best of the authors’ knowledge, no general tool has been reported for the multiple-parameter asymptotic behaviour analysis so far. 5. Numerical examples We now give some illustrative examples concerning the properties as well as the results presented in this article. Example 5.1 (Invariance property) Consider a time-delay system involving two delay parameters $$\tau_1$$ and $$\tau_2$$ with the characteristic function   $f(\lambda ,\overrightarrow \tau ) = {e^{ - ({\tau _1} + {\tau _2})\lambda }} - ({\lambda ^6} - {\lambda ^4} + {\lambda ^2}){e^{ - {\tau _2}\lambda }} - ({\lambda ^{10}} - {\lambda ^8} + {\lambda ^6}){e^{ - {\tau _1}\lambda }} + {\lambda ^{12}}.$ Suppose $$\tau_2$$ is fixed as $$\tau_2^ \#= 2 \pi$$ and $$\tau_1$$ is the free delay parameter. We here illustrate the invariance property. The corresponding characteristic function can be expressed as $$f(\lambda ,{\tau _1},{F_1}) = {\widetilde a_{10}}(\lambda ) + {\widetilde a_{11}}(\lambda ){e^{ - {\tau _1}\lambda }}$$ with $${\widetilde a_{10}}(\lambda ) = {\lambda ^{12}} - ({\lambda ^6} - {\lambda ^4} + {\lambda ^2}){e^{ - {\tau _2^ \#}\lambda }}$$ and $${\widetilde a_{11}}(\lambda )={e^{ - {\tau _2^ \#}\lambda }} - ({\lambda ^{10}} - {\lambda ^8} + {\lambda ^6})$$. At $$\tau_1 = (2k + 1)\pi$$, $$\lambda=j$$ is a CIR and, in particular, $$\lambda=j$$ is a triple CIR at $$\tau_1 = \pi$$ (it is simple at all $$\tau_1 = (2k + 1)\pi, k \in \mathbb{N}_+$$). According to Theorem 3.2, $${\it{\Delta}} NU_j^{{F_1}}((2k + 1)\pi ) = {\it{\Delta}} NF_{ - 1}^{{F_1}}(1) = - 1$$ for all $$k \in \mathbb{N}$$, where $${\it{\Delta}} NF_{ - 1}^{{F_1}}(1)$$ can be obtained from the FSC $${{\it{\Gamma}} _{1,1}}$$ shown in Fig. 1. To verify the result, we give the Puiseux series for the critical pair $$(\,j,\pi)$$:   ${\it{\Delta}} \lambda = (0.2978 + 0.1019j){({\it{\Delta}} {\tau _1})^{\frac{1}{3}}} + o({({\it{\Delta}} {\tau _1})^{\frac{1}{3}}})$ and the Taylor series for the critical pairs $$(\,j,(2k + 1)\pi )$$ (due to the limited space, we only give the results for $$k=1$$ and $$k=2$$):   ${{\it{\Delta}} \lambda = - 0.1592j{\it{\Delta}} {\tau _1} + 0.0253j{{({\it{\Delta}} {\tau _1})}^2} + ( - 0.0113 + 0.0132j){{({\it{\Delta}} {\tau _1})}^3} + o({{({\it{\Delta}} {\tau _1})}^3}),k = 1,}$   $\begin{array}{*{20}{l}} {{\it{\Delta}} \lambda = - 0.0796j{\it{\Delta}} {\tau _1} + 0.0063j{{({\it{\Delta}} {\tau _1})}^2} + ( - 0.0007 + 0.0006j){{({\it{\Delta}} {\tau _1})}^3} + o({{({\it{\Delta}} {\tau _1})}^3}),k = 2.} \end{array}$ It is worth mentioning that our result is consistent with the above series analysis. □ Fig. 1. View largeDownload slide FSC for Example 5.1. Fig. 1. View largeDownload slide FSC for Example 5.1. Example 5.2 (Complete stability) Consider a time-delay system with the characteristic function   $f(\lambda ,\overrightarrow \tau ) = {\lambda ^4} + 2{\lambda ^2} + 3{e^{ - {\tau _1}\lambda }} - 3{e^{ - {\tau _2}\lambda }} + {e^{ - ({\tau _1} + {\tau _2})\lambda }}.$ We study the complete stability w.r.t. $$\tau_2$$ when $$\tau_1$$ is fixed as $$\tau_1^ \#=0.5$$. First, when $$\overrightarrow \tau = \overrightarrow 0$$, the system has four CIRs (both $$+j$$ and $$-j$$ are double CIRs). Owing to the conjugate symmetry of the spectrum, it suffices to consider the asymptotic behaviour of the CIR $$+j$$. Applying the method in Chapter 4 of Li et al. (2015), we have that the asymptotic behaviour of the CIR $$+j$$ w.r.t. $$\tau_1$$ at $$\overrightarrow \tau = \overrightarrow 0$$ corresponds to the Puiseux series $${\it{\Delta}} \lambda = {( - j{\it{\Delta}} {\tau _1})^{\frac{1}{2}}} + o({({\it{\Delta}} {\tau _1})^{\frac{1}{2}}})$$. Thus, in the light of Theorem 3.1, $$NU(( + \varepsilon ,0)) = + 2$$. Next, through the FSC $${{\it{\Gamma}} _{1,1}}$$ with $${F_1} = (0,0)$$ (Fig. 2a), we have that $$NU((0.5,0)) = + 2$$ (two sets of critical pairs are obtained from Fig. 2(a): $$({\lambda _{1,0}} = j,{\tau _{1,0,k}} = 2k\pi )$$ and $$({\lambda _{1,1}} = 1.9566j,{\tau _{1,1,k}} = 1.6056 + \frac{{2k\pi }}{{1.9566}})$$). Finally, we generate the FSC $${{\it{\Gamma}} _{2,1}}$$ with $${F_2} = (0.5,0)$$ (Fig. 2(b)) and we have that the system has only one set of critical pairs: $$({\lambda _{2,0}} = j,{\tau _{2,0,k}} = {\rm{0}}{\rm{.9443}} + 2k\pi )$$. According to Theorem 3.2, each time $$\tau_2$$ increases near $${\tau _{2,0,k}}$$, $${\it{\Delta}} NU_{{\lambda _{2,0}}}^{{F_2}}({\tau _{2,0,k}}) = 0$$. Thus, $$NU((0.5,{\tau _2})) = + 2$$ for all $$\tau_2 \ge 0$$ except at $$\tau_2 = {\tau _{2,0,k}}$$. Finally, it is worth mentioning that the above analysis is consistent with the SCS, shown in Fig. 3. It is interesting to see that, in the three regions partitioned by the SCS in Fig. 3, the system has the same number of unstable roots.□ Fig. 2. View largeDownload slide FSCs for Example 5.2. (a) $${{\it{\Gamma}} _{1,1}}$$. (b) $${{\it{\Gamma}} _{2,1}}$$. Fig. 2. View largeDownload slide FSCs for Example 5.2. (a) $${{\it{\Gamma}} _{1,1}}$$. (b) $${{\it{\Gamma}} _{2,1}}$$. Fig. 3. View largeDownload slide Stability crossing set for Example 5.2. Fig. 3. View largeDownload slide Stability crossing set for Example 5.2. Example 5.3 (Stability analysis for a given delay vector) Consider the time-delay system in VII.B of Gu & Naghnaeian (2011) including three delays, with the characteristic function   $f(\lambda ,\overrightarrow \tau ) = {\lambda ^3} + 3\lambda + 7 + ({\lambda ^2} + 3\lambda + 1){e^{ - {\tau _1}\lambda }} + (4\lambda + 3){e^{ - {\tau _2}\lambda }} + ({\lambda ^2} + \lambda + 0.1){e^{ - {\tau _3}\lambda }}.$ In the sequel, our objective is to analyse the stability for a given delay vector $$\overrightarrow {{\tau ^ \#}} = (0.3041,1.7,0.0504)$$. As seen from Fig. 14 in Gu & Naghnaeian (2011), this $$\overrightarrow {{\tau ^ \#}}$$ corresponds to an intersecting point of the SCS. We now apply the procedure proposed in Section 4. The first iteration (for $$\chi = 1$$, after Step 0) $$\Rightarrow$$ Step 1: $$NU(( + \varepsilon ,0,0)) = 0$$ according to Theorem 3.1 (when $$\overrightarrow \tau = \overrightarrow 0$$, the characteristic roots are $$- 0.4457{\rm{ }} \pm {\rm{ }}3.1326j$$ and $$-1.1087$$). Step 2: We generate the FSC $${{\it{\Gamma}} _{1,1}}$$ (Fig. 4(a)), from which we obtain two sets of critical pairs: $$({\lambda _{1,0}} = 2.2819j,{\tau _{1,0,k}} = 1.9052 + \frac{{2k\pi }}{{2.2819}})$$ and $$({\lambda _{1,1}} = 3.5160j,{\tau _{1,1,k}} = 0.2756 + \frac{{2k\pi }}{{3.5160}})$$. We have that $$NU(({\rm{0}}{\rm{.3041}},0,0)) = +2$$ and the system has no CIRs at $$\overrightarrow \tau = (0.3041,0,0)$$. Step 3: Let $$\chi=2$$ and go to Step 1. The second iteration (for $$\chi = 2$$) $$\Rightarrow$$ Step 1: $$NU(({\rm{0}}{\rm{.3041}}, + \varepsilon ,0)) = +2$$ according to Theorem 3.1. Step 2: We generate the FSC $${{\it{\Gamma}} _{2,1}}$$ (Fig. 4(b)), from which we obtain two sets of critical pairs: $$({\lambda _{2,0}} = 1.9073j,{\tau _{2,0,k}} = 1.7471 + \frac{{2k\pi }}{{1.9073}})$$ and $$({\lambda _{2,1}} = 3.5133j,{\tau _{2,1,k}} = 1.7577 + \frac{{2k\pi }}{{3.5133}})$$. We have that $$NU(({\rm{0}}.{\rm{3041}},1.7,0)) = +2$$ and the system has no CIRs at $$\overrightarrow \tau = (0.3041,1.7,0)$$. Step 3: Let $$\chi=3$$ and go to Step 1. The third iteration (for $$\chi = 3$$) $$\Rightarrow$$ Step 1: $$NU(({\rm{0}}.{\rm{3041}},1.7, + \varepsilon )) = +2$$ according to Theorem 3.1. Step 2: We generate the FSC $${{\it{\Gamma}} _{3,1}}$$ (Fig. 5(a)), from which we obtain two sets of critical pairs: $$({\lambda _{3,0}} = 1.9594j,{\tau _{3,0,k}} = 0.0504 + \frac{{2k\pi }}{{1.9594}})$$ and $$({\lambda {3,1}} = 3.6118j,{\tau _{3,1,k}} = 0.0504 + \frac{{2k\pi }}{{3.6118}})$$. In light of $${{\it{\Gamma}} _{3,1}}$$, each time $$\tau_3$$ increases from $${\tau _{3,0,k}} - \varepsilon$$ to $${\tau _{3,0,k}} + \varepsilon$$ ($${\tau _{3,1,k}} - \varepsilon$$ to $${\tau _{3,1,k}} + \varepsilon$$), $${\lambda _{3,0}}$$ ($${\lambda _{3,1}}$$) rosses $$\mathbb{C}_0$$ from $$\mathbb{C}_+$$ to $$\mathbb{C}_-$$ (from $$\mathbb{C}_-$$ to $$\mathbb{C}_+$$). The complete stability w.r.t. $$\tau_3$$ can be studied. We may accurately calculate $$NU((\tau _1^\#,\tau _2^ \#,{\tau _3}))$$ as a function of $$\tau _3$$, as shown in Fig. 5(b). In particular, near $$\tau_3^ \#=0.0504$$, the variation of $$\tau_3$$ does not change $$NU((\tau _1^ \#, \tau _2^ \#, {\tau _3}))$$ since it causes opposite crossing directions (w.r.t. the imaginary axis $$\mathbb{C}_0$$) for $${\lambda _{3,0}}$$ and $${\lambda _{3,1}}$$. From the above iterations, we know that at $$\overrightarrow {{\tau ^ \#}}$$ the system has two CIRs ($$1.9594j$$ and $$3.6118j$$), without characteristic roots in $$\mathbb{C}_+$$. We next study the asymptotic behaviour of the CIRs w.r.t. $$\tau_1$$ and $$\tau_2$$, respectively, in order to see if there is a stabilizing point near $$\overrightarrow {{\tau ^ \#}}$$. Step 4 $$\Rightarrow$$ The asymptotic behaviour w.r.t. $$\tau_1$$ can be studied from the FSC $${\it{\Gamma}} _{1,1}^ \#$$ (Fig. 6(a)), i.e., $${\it{\Gamma}} _{1,1}$$ when $${F_1} = (0,1.7,0.0504)$$ as defined in the procedure. From Fig. 6(a), we have that as $$\tau_1$$ increases from $$0.3041 - \varepsilon$$ to $$0.3041 + \varepsilon$$ both the CIRs ($$1.9594j$$ and $$3.6118j$$) cross $$\mathbb{C}_0$$ from $$\mathbb{C}_-$$ to $$\mathbb{C}_+$$. Therefore, decreasing $$\tau_1$$ appropriately near $$\overrightarrow {{\tau ^ \#}}$$ may stabilize the system. Next, the asymptotic behaviour w.r.t. $$\tau_2$$ can be studied from the FSC $${\it{\Gamma}} _{2,1}^ \#$$ (Fig. 6(b)), i.e., $${\it{\Gamma}} _{2,1}$$ when $${F_2} = (0.3041,0,0.0504)$$. Similar to the asymptotic behaviour w.r.t. $$\tau_3$$, the variation of $$\tau_2$$ near $$\tau _2^ \#$$ dose not change $$NU((\tau _1^ \#, {\tau _2}, \tau _3^ \#))$$. Thus, there exists a stabilizing $$\overrightarrow \tau$$ sufficiently near $$\overrightarrow {{\tau ^ \#}}$$ with $${\tau _1} > \tau _1^ \#$$. □ Fig. 4. View largeDownload slide FSCs (a) $${{\it{\Gamma}} _{1,1}}$$ and (b) $${{\it{\Gamma}} _{2,1}}$$ for Example 5.3. Fig. 4. View largeDownload slide FSCs (a) $${{\it{\Gamma}} _{1,1}}$$ and (b) $${{\it{\Gamma}} _{2,1}}$$ for Example 5.3. Fig. 5. View largeDownload slide (a) $${{\it{\Gamma}} _{3,1}}$$ and (b) $$NU((\tau _1^ \#, \tau _2^ \#, {\tau _3}))$$ vs. $$\tau_3$$ for Example 5.3. Fig. 5. View largeDownload slide (a) $${{\it{\Gamma}} _{3,1}}$$ and (b) $$NU((\tau _1^ \#, \tau _2^ \#, {\tau _3}))$$ vs. $$\tau_3$$ for Example 5.3. Fig. 6. View largeDownload slide (a) $${\it{\Gamma}} _{1,1}^ \#$$ and (b) $${\it{\Gamma}} _{2,1}^ \#$$ for Example 5.3. Fig. 6. View largeDownload slide (a) $${\it{\Gamma}} _{1,1}^ \#$$ and (b) $${\it{\Gamma}} _{2,1}^ \#$$ for Example 5.3. Finally, at the end of this section, we give some additional discussions and comments. As mentioned and as seen above, the frequency-sweeping test proposed in this article is iterative. In this context, one natural question may arise: How about a direct approach (through the complete stability analysis along an appropriate ray) in the corresponding delay parameter space? For instance, consider the stability for Example 5.2 when $$\tau _1^\# = 0.5,\tau _2^\# = 1.5$$. We may try to study the complete stability w.r.t the ray $$\tau \cdot (1,3)$$, $$\tau \in [0,\infty )$$. Here, $$\tau$$ is the only free parameter. Obviously, the point ($$\tau _1^\# = 0.5,\tau _2^\# = 1.5$$) corresponds to the case $$\tau \cdot (1,3)$$ with $$\tau = 0.5$$. The characteristic function is $${\lambda ^4} + 2{\lambda ^2} + 3z - 3{z^3} + {z^4}$$ ($$z = {e^{ - \tau \lambda }}$$). There are totally four FSCs, as shown in Fig. 7(a). Fig. 7. View largeDownload slide FSCs for complete stability analysis along a ray. (a) Ray $$\tau \cdot (1,3)$$. (b) Ray $$\tau \cdot (5,14)$$. Fig. 7. View largeDownload slide FSCs for complete stability analysis along a ray. (a) Ray $$\tau \cdot (1,3)$$. (b) Ray $$\tau \cdot (5,14)$$. If we consider the point ($$\tau _1^\# = 0.5,\tau _2^\# = 1.4$$) similarly to the previous case study, we shall address the ray $$\tau \cdot (5, 14)$$, $$\tau \in [0,\infty )$$. The corresponding characteristic function is $${\lambda ^4} + 2{\lambda ^2} + 3{z^5} - 3{z^{14}} + {z^{19}}$$ ($$z = {e^{ - \tau \lambda }}$$). As shown in Fig. 7(b), there are totally 19 FSCs! Apparently, the analysis appears to be quite involved with some (significant) increasing complexity. As seen above, if there exists a common factor among all the delay parameters, the complete stability analysis along a ray is possible in theory. However, this is not necessarily simple or practical. For instance, if we try to address Example 5.3, the number of the FSCs is 17000 (the largest common factor among the three delays is $$0.0001$$). The above analysis is consistent with the delay interference phenomenon (Michiels & Niculescu, 2007) and we are able to link it with the number of FSCs. More precisely, a very small perturbation on the delay ratio may result in an increase in the number of FSCs. This will tend to bring more CIRs and hence an increase in $$NU(\tau)$$ for a sufficiently large $$\tau$$, according to the invariance property. 6. Algorithm implementation The effectiveness and advantage of the iterative frequency-sweeping approach are illustrated and discussed in the previous sections. In this section, we will try to implement all the steps by a single program. In this way, the stability can be automatically treated. In this section, we consider the following characteristic function involving $$L$$ independent delay parameters ($${\tau _1}, \ldots ,{\tau _L}$$)   $$\label{CFMTDSwithoutCrossterms}f(\lambda ,{\tau _1}, \ldots ,{\tau _L}) = {p_0}(\lambda ) + \sum\limits_{\ell = 1}^L {{p_\ell }(\lambda ){e^{ - {\tau _\ell }\lambda }}} ,$$ (6.1) where $${p_0}(\lambda ), \ldots, {p_L}(\lambda )$$ are real-coefficient polynomials of $$\lambda$$ such that   $\deg ({p_0}(\lambda )) > \max \{ \deg ({p_1}(\lambda )), \ldots ,\deg ({p_L}(\lambda ))\} .$ The characteristic function (6.1) corresponds to a class of multiple-delay systems of retarded type without cross terms in the characteristic functions. As earlier mentioned, the stability of such multiple-delay systems have been largely studied in the literature (see e.g., Gu et al., 2005 for the case $$L=2$$ and Gu & Naghnaeian, 2011 for the case $$L=3$$. Based on the approach proposed in this article, we now present an algorithm to automatically determine $$NU(\overrightarrow {{\tau ^\# }} )$$ for a given $$\overrightarrow {{\tau ^\# }} = (\tau _1^\# , \ldots ,\tau _L^\# )$$. Algorithm for automatic implementation of iterative frequency-sweeping approach Step 0: Choose a step-length $${{\it{\Delta}} _\omega }$$. Set $$\chi=1$$ and $${F_1} = (0, \ldots ,0)$$. Step 1: Compute $$NU({F_\chi } + \varepsilon \delta (\chi ))$$ by Theorem 3.1. (Step 1 usually may be automatically performed by computer, through solving the polynomial equation $$f(\lambda ,\overrightarrow 0 )$$. Only when the spectrum contains CIRs, we need to additionally invoke the Puiseux series.) Step 2: Sweep $$\omega$$ with the step-length $${{\it{\Delta}} _\omega }$$ and solve $$z_\chi$$ for the equation $$p(\lambda ,{z _\chi},{F_\chi}) = 0$$. Detect the signs of $$\left| {{z_\chi }} \right| - 1$$. If a sign change is found at two adjacent $$\omega$$, say $${\omega '}$$ and $${\omega ''}$$, a CIR is detected as $$\frac{{(\omega ' + \omega '')}}{2}j$$. For any critical pair ($$\lambda, \tau$$) detected in this way, we have that $${\it{\Delta}} N{U_\lambda }(\tau ) = + 1$$ ($$-1$$) if the corresponding sign change is from negative to positive (from positive to negative). Then, we know $$NU({F_\chi } + \tau _\chi ^ \# \delta (\chi ))$$ and detect the CIRs (if any). (The case $${\it{\Delta}} N{U_\lambda }(\tau ) = 0$$ may be later examined from the FSCs. Note that such a case usually does not affect the value of $$NU(\overrightarrow {{\tau ^\# }} )$$). Step 4: If $$\chi < L$$, let $$\chi = \chi + 1$$ and $${F_\chi } = {F_{\chi - 1}} + \tau _{\chi - 1}^ \# \delta (\chi - 1)$$. Return to Step 1. Step 5: We obtain the value of $$NU(\overrightarrow {{\tau ^\# }} )$$ and plot the FSCs $${{\it{\Gamma}} _{\chi ,1}}$$, $$\chi=1, \ldots, L$$. The FSCs may be used to determine if there is the case with $${\it{\Delta}} N{U_\lambda }(\tau ) = 0$$ in Step 2. Also, the FSCs may verify the results of Step 2. The algorithm stops. If there are CIRs at $$\overrightarrow {{\tau ^\# }}$$, we may further generate the FSCs $${\it{\Gamma}} _{\ell ,1}^\#$$ as in Step 4 of the iterative frequency-sweeping approach, to analyse the asymptotic behaviour. Remark 6.1 By using the above algorithm, the calculation error for CIRs is kept within $$\pm \frac{{{\it{\Delta}} _\omega }}{2}j$$. Example 6.1 Consider the three-delay system in Example 5.3. We now compute $$NU(\overrightarrow \tau )$$ at some points, by using a MATLAB program (based on the algorithm proposed in this section). The step-length is chosen as $${{\it{\Delta}} _\omega }=0.01$$. In Table 1, we list the results and the computation time (on a Laptop with an Intel Core 2.50 GHz CPU with 8 G RAM). For these points, the analysis is automatically performed by computer. Table 1 $$NU((\tau _1^\# ,\tau _2^\# ,\tau _3^\# ))$$ calculation for some representative points $$(\tau _1^\# = 0.01$$ and $$\tau _2^\# = 1.7)$$   $$NU((\tau _1^\# ,\tau _2^\# ,\tau _3^\# ))$$  Computation time (s)  $$\tau _3^\# = 0.1$$  2  0.009605  $$\tau _3^\# = 0.3$$  0  0.009712  $$\tau _3^\# = 0.6$$  2  0.009579  $$\tau _3^\# = 0.9$$  0  0.009573  $$\tau _3^\# = 1.5$$  2  0.009616  $$\tau _3^\# = 2.5$$  4  0.009555  $$\tau _3^\# = 3.0$$  2  0.009604  $$\tau _3^\# = 3.7$$  0  0.009571  $$\tau _3^\# = 3.77$$  2  0.009547  $$\tau _3^\# = 3.8$$  4  0.009589    $$NU((\tau _1^\# ,\tau _2^\# ,\tau _3^\# ))$$  Computation time (s)  $$\tau _3^\# = 0.1$$  2  0.009605  $$\tau _3^\# = 0.3$$  0  0.009712  $$\tau _3^\# = 0.6$$  2  0.009579  $$\tau _3^\# = 0.9$$  0  0.009573  $$\tau _3^\# = 1.5$$  2  0.009616  $$\tau _3^\# = 2.5$$  4  0.009555  $$\tau _3^\# = 3.0$$  2  0.009604  $$\tau _3^\# = 3.7$$  0  0.009571  $$\tau _3^\# = 3.77$$  2  0.009547  $$\tau _3^\# = 3.8$$  4  0.009589  Table 1 $$NU((\tau _1^\# ,\tau _2^\# ,\tau _3^\# ))$$ calculation for some representative points $$(\tau _1^\# = 0.01$$ and $$\tau _2^\# = 1.7)$$   $$NU((\tau _1^\# ,\tau _2^\# ,\tau _3^\# ))$$  Computation time (s)  $$\tau _3^\# = 0.1$$  2  0.009605  $$\tau _3^\# = 0.3$$  0  0.009712  $$\tau _3^\# = 0.6$$  2  0.009579  $$\tau _3^\# = 0.9$$  0  0.009573  $$\tau _3^\# = 1.5$$  2  0.009616  $$\tau _3^\# = 2.5$$  4  0.009555  $$\tau _3^\# = 3.0$$  2  0.009604  $$\tau _3^\# = 3.7$$  0  0.009571  $$\tau _3^\# = 3.77$$  2  0.009547  $$\tau _3^\# = 3.8$$  4  0.009589    $$NU((\tau _1^\# ,\tau _2^\# ,\tau _3^\# ))$$  Computation time (s)  $$\tau _3^\# = 0.1$$  2  0.009605  $$\tau _3^\# = 0.3$$  0  0.009712  $$\tau _3^\# = 0.6$$  2  0.009579  $$\tau _3^\# = 0.9$$  0  0.009573  $$\tau _3^\# = 1.5$$  2  0.009616  $$\tau _3^\# = 2.5$$  4  0.009555  $$\tau _3^\# = 3.0$$  2  0.009604  $$\tau _3^\# = 3.7$$  0  0.009571  $$\tau _3^\# = 3.77$$  2  0.009547  $$\tau _3^\# = 3.8$$  4  0.009589  The results listed in Table 1 may be verified by Fig. 14 in Gu & Naghnaeian (2011) (see also Fig. 8(a) in this article). With some slight modifications of the program mentioned above, we may obtain the SCS for the three-delay system. Here, we generate the SCS for the case $$\tau_2 = 1.7$$, see Fig. 8(a). Figure 8(a) obtained here is same as Fig. 14 of Gu & Naghnaeian (2011), which is generated by a different approach using some geometric arguments. The computation time to generate the data for Fig. 8(a) by our algorithm is 2.566361 s (on the same Laptop). It is worth to mention that the $$NU(\overrightarrow \tau )$$ distribution is directly obtained by our program (without extra calculation), since the asymptotic behaviour analysis for CIRs is covered by the approach in this article. Moreover, we can further determine the 3-D SCS in the $$(\tau_1, \tau_2, \tau_3)$$-space. For a clear illustration, we here give the SCS for a domain $$({\tau _1},{\tau _2},{\tau _3}) \in [1.5, 2.5] \times [0, 3] \times [0, 3]$$. The 3-D SCS is shown in Fig. 8(b). The computation time to obtain the data for Fig. 8(b) is 35.552933 s (on the same Laptop).□ Fig. 8. View largeDownload slide Stability crossing set for Example 6.1. (a) Cross section with $$\tau_2 = 1.7$$. (b) In $$(\tau_1, \tau_2, \tau_3)$$-space. Fig. 8. View largeDownload slide Stability crossing set for Example 6.1. (a) Cross section with $$\tau_2 = 1.7$$. (b) In $$(\tau_1, \tau_2, \tau_3)$$-space. Remark 6.2 As seen in Gu & Naghnaeian (2011), it is not trivial to determine the SCS for the multiple-delay system (6.1). In this section, we propose a different approach for this task. As illustrated, the asymptotic behaviour analysis for CIRs is included by our algorithm and hence the distribution of $$NU(\tau)$$ can be directly examined. In our opinion, this is one of the advantages of the approach proposed in our article. We think that this observation will open some new perspectives for further research. For instance, we may study how to obtain the SCSs for more general multiple-delay systems (e.g., when cross terms $${e^{ - ({\tau _1} + {\tau _2})\lambda }}$$, $${e^{ - ({\tau _1} + {\tau _3})\lambda }}$$, $${e^{ - ({\tau _2} + {\tau _3})\lambda }}$$, $${e^{ - ({\tau _1} + {\tau _2} + {\tau _3})\lambda }}$$, $$\ldots$$, are added in (6.1)) in the future. But this scope is out of our purposes in this article. 7. Conclusion We analyse the stability for linear systems with multiple (incommensurate) delay parameters. As the asymptotic behaviour analysis of the CIRs w.r.t multiple delay parameters corresponds to an open problem, we propose an indirect yet effective methodology called the iterative frequency-sweeping approach. We first study the complete stability in the case where only one delay parameter is free while the others are fixed. The invariance property (regarding the asymptotic behaviour of the CIRs) in this case is confirmed by extending the frequency-sweeping framework recently proposed for studying linear systems with single delay parameter. As a result, the complete stability problem can be easily studied by employing a frequency-sweeping test. Based on the above results, we next present an iterative frequency-sweeping approach to analyse the stability for any given combination of multiple delays. Using this approach, we may accurately compute the number of unstable characteristic roots. Furthermore, if the system has CIRs, we may analyse the asymptotic behaviour of the CIRs w.r.t each delay parameter. Consequently, we may determine if there exists a stabilizing combination of multiple delays sufficiently close to the given one, at which the system is asymptotically stable. Finally, we develop an algorithm with which the stability for a class of multiple-delay systems may be easily implemented. This work opens some new perspectives for further research. Funding National Natural Science Foundation of China (61473065); and Fundamental Research Funds for the Central Universities of China (N160402001), in part. Footnotes 1 Roughly speaking, the complete stability problem refers to finding exhaustively the stability interval(s) for the delay parameter $$\tau$$ along the whole $$\tau$$-axis $$\tau \in [0,\infty )$$. 2 A system has CIRs iff the delay parameters lie in the SCS. 3 Some results for specific two-free-parameter problems can be found in, e.g., Section 6.2.1 of Arnold et al. (2012), Beringer & Richard-Jung (2003) and Soto & Vicente (2011). For more general problems, the analysis will become much more complicated and no general result has been reported so far. 4 We do not need to specifically generate the FSCs $${\it{\Gamma}} _{L, i}^ \# (\omega )$$ since they are exactly the FSCs obtained by Step 3 when $$\chi = L$$. References Arnold V. I., Gusein-Zade S. M. & Varchenko A. N. ( 2012) Singularities of Differentiable Maps.  vol. 2: Monodromy and Asymptotics of Integrals. New York: Springer. Beringer F. & Richard-Jung F. ( 2003) Multi-variate polynomials and Newton-Puiseux expansions. Symbolic and Numerical Scientific Computation . ( Winkler F. & Langer U. eds). Berlin: Springer, pp. 240– 254. Google Scholar CrossRef Search ADS   Boussaada I. & Niculescu S. -I. ( 2016a) Characterizing the codimension of zero singularities for time-delay systems: a link with Vandermonde and Birkhoff incidence matrices. Acta Appl. Math. , 145, 47– 88. Google Scholar CrossRef Search ADS   Boussaada I. & Niculescu S. -I. ( 2016b) Tracking the algebraic multiplicity of crossing imaginary roots for generic quasipolynomials: A Vandermonde-based approach. IEEE Trans. Autom. Control , 61, 1601– 1606. Google Scholar CrossRef Search ADS   Cooke K. L. & van den Driessche P. ( 1986) On zeros of some transcendental equations. Funkcialaj Ekvacioj , 29, 77– 90. Delice I. I. & Sipahi R. ( 2012) Delay-independent stability test for systems with multiple time-delays. IEEE Trans. Autom. Control,  57, 963– 972. Google Scholar CrossRef Search ADS   Gu K. & Naghnaeian M. ( 2011) Stability crossing set for systems with three delays. IEEE Trans. Autom. Control,  56, 11– 26. Google Scholar CrossRef Search ADS   Gu K., Niculescu S. -I. & Chen J. ( 2005) On stability crossing curves for general systems with two delays. J. Math. Anal. Appl. , 311, 231– 253. Google Scholar CrossRef Search ADS   Hale J. K. & Verduyn Lunel S. M. ( 1993) Introduction to Functional Differential Equations.  New York: Springer. Google Scholar CrossRef Search ADS   Hassard B. D. ( 1997) Counting roots of the characteristic equation for linear delay-differential systems. J. Differ. Equ. , 136, 222– 235. Google Scholar CrossRef Search ADS   Hu G. -D. & Liu M. ( 2007) Stability criteria of linear neutral systems with multiple delays. IEEE Trans. Autom. Control , 52, 720– 724. Google Scholar CrossRef Search ADS   Lee M. S. & Hsu C. S. ( 1969) On the $$\tau$$-decomposition method of stability analysis for retarded dynamical systems. SIAM J. Control , 7, 242– 259. Google Scholar CrossRef Search ADS   Li X. -G., Niculescu S. -I. & Çela A. ( 2015) Analytic Curve Frequency-Sweeping Stability Tests for Systems with Commensurate Delays.  London: Springer. Google Scholar CrossRef Search ADS   Li X. -G., Niculescu S. -I., Çela A, Zhang L. & Li X. ( 2017) A frequency-sweeping framework for stability analysis of time-delay systems. IEEE Trans. Autom. Control , 62, 3701– 3716. Google Scholar CrossRef Search ADS   Michiels W. & Niculescu S. -I. ( 2007) Characterization of delay-independent stability and delay interference phenomena. SIAM J. Control Optim. , 45, 2138– 2155. Google Scholar CrossRef Search ADS   Michiels W. & Niculescu S. -I. ( 2014) Stability, Control, and Computation for Time-Delay Systems: An Eigenvalue-Based Approach.  Philadelphia: SIAM. Google Scholar CrossRef Search ADS   Naghnaeian M. & Gu K. ( 2013) Stability crossing set for systems with two scalar-delay channels. Automatica , 49, 2098– 2106. Google Scholar CrossRef Search ADS   Neimark J. ( 1949) D-subdivisions and spaces of quasi-polynomials. Prikl. Mat. Meh. , 13, 349– 380. Niculescu S. -I. ( 2001) Delay Effects on Stability: A Robust Control Approach.  Heidelberg: Springer. Niculescu S. -I., Kim P., Gu K., Lee P. & Levy D. ( 2010) Stability crossing boundaries of delay systems modeling immune dynamics in leukemia. Discrete Contin. Dyn. Syst.-Ser. B (DCDS-B) , 15, 129– 156. Sipahi R. & Delice I. I. ( 2009) Extraction of 3D stability switching hypersurfaces fo a time delay system with multiple fixed delays. Automatica , 45, 1449– 1454. Google Scholar CrossRef Search ADS   Sipahi R. & Delice I. I. ( 2011) Advanced clustering with frequency sweeping methodology for the stability analysis of multiple time-delay systems. IEEE Trans. Autom. Control , 56, 467– 472. Google Scholar CrossRef Search ADS   Sipahi R., Niculescu S. -I., Abdallah C. T., Michiels W. & Gu K. ( 2011) Stability and stabilization of systems with time delay. IEEE Control Syst. Mag. , 31, 38– 65. Google Scholar CrossRef Search ADS   Sipahi R. & Olgac N. ( 2006) A unique methodology for the stability robustness of multiple time delay systems. Syst. Control Lett. , 55, 819– 825. Google Scholar CrossRef Search ADS   Soto M. J. & Vicente J. L. ( 2011) The Newton procedure for several variables. Linear Algebra Appl. , 435, 255– 269. Google Scholar CrossRef Search ADS   Stépán G. ( 1979) On the stability of linear differential equations with delay. Coll. Math. Soc. J. Bolyai , 30, 971– 984. Stépán G. ( 1989) Retarded Dynamical Systems: Stability and Characteristic Functions.  Harlow: Longman Scientific & Technical. Vyhlídal T. & Zítek P. ( 2009) Modification of Mikhaylov criterion for neutral time-delay systems. IEEE Trans. Autom. Control , 54, 2430– 2435. Google Scholar CrossRef Search ADS   Appendix: Proof of Theorem 3.2 Consider the following characteristic function   \begin{eqnarray} \label{GeneralQuasipolynomial}f(\lambda ,\tau ) = {a_0}(\lambda ) + {a_1}(\lambda ){e^{ - \tau \lambda }} + \cdots + {a_q}(\lambda ){e^{ - q\tau \lambda }}, \end{eqnarray} (A.1) where the coefficient functions $${a_0}(\lambda ), \ldots ,{a_q}(\lambda )$$are only required to be analytic in $$\mathbb{C}_0 \backslash \{ 0\}$$ (have in mind that usually we preclude the trivial case where $$\lambda = 0$$ is a characteristic root). The above characteristic function (A.1) is called a general quasipolynomial, corresponding to a broad class of time-delay systems. Note that the characteristic function (A.1) reduces to the widely-studied quasipolynomials, corresponding to the retarded-type and the neutral-type time-delay systems, if the coefficient functions $${a_0}(\lambda ), \ldots ,{a_q}(\lambda )$$ are restricted to be polynomials of $$\lambda$$. Recently, the invariance property of the CIRs for the general quasipolynomial (A.1) was confirmed in Li et al. (2017). Clearly, the characteristic function $$f(\lambda ,{\tau _\chi },{F_\chi })$$ (3.1) is in the form of (A.1) as   ${\widetilde a_{\chi 0}}(\lambda ) + {\widetilde a_{\chi 1}}(\lambda ){e^{ - {\tau _\chi }\lambda }} + \cdots + {\widetilde a_{\chi {q_\chi }}}(\lambda ){e^{ - {q_\chi }{\tau _\chi }\lambda }},$ where the coefficient functions $${\widetilde a_{\chi 0}}(\lambda ), \ldots ,{\widetilde a_{\chi {q_\chi }}}(\lambda )$$ are polynomials in $$\lambda$$ and $${e^{ - \tau _\ell ^\# \lambda }}$$ ($$\ell \in \{ 1, \ldots ,L\}$$, $$\ell \ne \chi$$). We may now prove Theorem 3.2 as the characteristic function $$f(\lambda ,{\tau _\chi },{F_\chi })$$ falls in the class of general quasipolynomial (A.1), since the coefficient functions $${\widetilde a_{\chi 0}}(\lambda ), \ldots ,{\widetilde a_{\chi {q_\chi }}}(\lambda )$$ are analytic in $$\mathbb{C}$$. © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png IMA Journal of Mathematical Control and Information Oxford University Press

# An iterative frequency-sweeping approach for stability analysis of linear systems with multiple delays

, Volume Advance Article – Nov 22, 2017
20 pages

/lp/ou_press/an-iterative-frequency-sweeping-approach-for-stability-analysis-of-6i6iB2PROf
Publisher
Oxford University Press
ISSN
0265-0754
eISSN
1471-6887
D.O.I.
10.1093/imamci/dnx050
Publisher site
See Article on Publisher Site

### Abstract

Abstract In this article, we study the stability of linear systems with multiple (incommensurate) delays, by extending a recently proposed frequency-sweeping approach. First, we consider the case where only one delay parameter is free while the others are fixed. The complete stability w.r.t. the free delay parameter can be systematically investigated by proving an appropriate invariance property. Next, we propose an iterative frequency-sweeping approach to study the stability under any given multiple delays. Moreover, we may effectively analyse the asymptotic behaviour of the critical imaginary roots (if any) w.r.t. each delay parameter, which provides a possibility for stabilizing the system through adjusting the delay parameters. The approach is simple (graphical test) and can be applied systematically to the stability analysis of linear systems including multiple delays. A deeper discussion on its implementation is also proposed. Finally, various numerical examples complete the presentation. 1. Introduction In this article, we consider the following linear time-delay system with multiple delay parameters   $$\label{MTDSsytem} \dot x(t) = Ax(t) + \sum\limits_{\ell = 1}^L {{B_\ell }} x(t - {\tau _\ell }),$$ (1.1) where $$A$$ and $$B_{\ell}$$ are constant matrices with compatible dimensions; $$\tau_{\ell} \ge 0$$, $$\ell = 1, \ldots ,L$$, are independent delay parameters ($$L$$ denotes the number of delay parameters). The delay combination may be expressed by a vector $$\overrightarrow \tau = ({\tau _1}, \ldots, {\tau _L})$$. The stability of time-delay system (1.1) has been largely investigated in the literature, see e.g., a survey paper Sipahi et al. (2011), or the books Michiels & Niculescu (2014) and Niculescu (2001), and the references therein. In the literature, in the framework of the so-called parameter-based approach, the stability problems for time-delay systems may be roughly divided in two categories: the $$\tau$$-decomposition problem (Lee & Hsu, 1969) and the $$D$$-decomposition problem (Neimark, 1949). For the former the delays are treated as free parameters, while for the latter some system/controller parameters are considered free (the delays are all fixed). For further discussions on such topics, we refer to Michiels & Niculescu (2014) and the references therein. The problem considered in this article belongs to the so-called $$\tau$$-decomposition problem. We start by recalling some of the existing results for the systems with single delay parameter. If $$L=1$$, the system (1.1) reduces to the form   $$\label{SingleDelaysytem} \dot x(t) = Ax(t) + Bx(t - \tau ),$$ (1.2) where $$B$$ is a constant matrix and $$\tau$$ is the delay parameter. The characteristic function for the system (1.2) is $$\det (\lambda I - A - B{e^{- \tau \lambda }})$$, where $$I$$ is the identity matrix of appropriate dimensions. If the delay parameters of system (1.1) are commensurate, the system can be expressed in the form   $$\label{CommensurateDelaysytem}\dot x(t) = Ax(t) + \sum\limits_{\ell = 1}^L {{B_\ell }} x(t - \ell \tau ),$$ (1.3) with the single delay parameter $$\tau$$. The characteristic function for system (1.3) writes as $$\det (\lambda I - A - \sum\limits_{\ell = 1}^L {{B_\ell }} {e^{ - \ell \tau \lambda }})$$. For the time-delay systems (1.2) and (1.3), the corresponding characteristic functions involve only one delay parameter (this is a major distinction with the multiple-delay system (1.1)). In such a case, the complete stability problem1 was systematically solved by the frequency-sweeping approach recently proposed in Li et al. (2015). The core result lies in that the invariance property concerning the asymptotic behaviour of the critical imaginary roots (CIRs) w.r.t. the infinitely many positive critical delays (CDs) was proved. It is worth mentioning that one of the main difficulties in analysing the complete stability lies in the characterization and the corresponding classification of the case with multiple characteristic roots on the imaginary axis. In the case where the delays are not commensurate, further discussions on characterizing multiple characteristic roots on the imaginary axis as well as related properties can be found in Boussaada & Niculescu (2016a,b). However, the stability problem for the multiple-delay system (1.1) is much more complicated and may have some unexpected dynamic behaviour (for instance, the delay interference phenomenon, see Michiels & Niculescu, 2007). Roughly speaking, to the best of the authors’ knowledge, the existing studies for the multiple-delay system (1.1) can be categorized into two classes: (i) One class of studies focus on finding the stability crossing set (SCS) in the delay parameter space.2 For a class of time-delay systems without cross-terms in the characteristic functions, the SCSs for the case $$L=2$$ and for the case $$L=3$$ are reported in Gu et al. (2005) and Gu & Naghnaeian (2011), respectively. Moreover, the geometric structure of the corresponding two- and three-dimensional SCSs are studied in detail. An extension to the two-delay case with one cross-term in the characteristic function was reported later in Naghnaeian & Gu (2013). By using some different arguments (based on the properties of the Rekasius transformation), another important series of results have been reported in, e.g., Sipahi & Delice (2009), Sipahi & Delice (2011) and Sipahi & Olgac (2006) on obtaining the SCSs (also called potential stability switching curves in these references). The methods are now applicable to general linear time-delay systems. In our opinion, the main advantage of the above class of studies is that the corresponding SCS in the case $$L=2$$ or $$3$$ may be visualized by a two- or three-dimensional figure, from which we may intuitively analyse the stability. If $$L > 3$$, a visualization of the SCS is difficult. In some cases (coupling small with large delays), some visualization was also proposed in the case of four delays describing immune dynamics in leukaemia (see e.g., Niculescu et al., 2010). (ii) Another class of methods concentrate on counting the number of unstable roots under a given delay vector $$\overrightarrow \tau$$. The idea lies in applying the argument principle to the corresponding characteristic function. An advantage of this class of studies is that the number of the delays $$L$$ is allowed to be any large. A representative conclusion can be found in Stépán (1979, 1989), known as the Stépán’s formula. Other interesting results in the same spirit can also be found in, e.g., Hassard (1997), Hu & Liu (2007) and Vyhlídal & Zítek (2009). However, when the system has CIRs, it is difficult to adopt the above two classes of methodologies to analyse the asymptotic behaviour of the CIRs. The difficulty is mainly related to the treatment of the case with multiple and/or degenerate CIRs. The asymptotic behaviour analysis is of practical and theoretical importance, especially when $$L$$ is large. When CIRs appear for a given $$\overrightarrow \tau$$, we usually want to know if there is a delay vector near $$\overrightarrow \tau$$ such that the system is asymptotically stable. If so, we will further consider how to find such a delay vector. From the above discussions, we see that there is still much room for the stability analysis of linear systems with multiple (incommensurate) delays, which motivates the work of this article. A straightforward idea is to extend the mathematical results for asymptotic behaviour analysis of the single-delay problem (proposed in Li et al., 2015) to the multiple-delay problem. However, to the best of the authors’ knowledge, such an extension is rather difficult as the asymptotic behaviour analysis for multiple-parameter problems has remained open in mathematics.3 In this article, we will propose an ‘indirect’ yet effective approach, called the iterative frequency-sweeping approach. The core of this approach is an important property to be proved in this article: If the multiple-delay system (1.1) has only one free delay parameter (the other ones are fixed), the effect of a CIR’s asymptotic behaviour on the stability w.r.t. its infinitely many positive CDs is invariant. This invariance property is generalized from the one confirmed in Li et al. (2015) for the systems with commensurate delays. To the best of the authors’ knowledge, such a result was not proposed in the open literature and it gives insights in understanding the way multiple delays may affect the stability property. As mentioned in the single-delay case, one of the main difficulties of the complete characterization of the stability is related to the asymptotic behaviour analysis of the multiple and/or degenerate CIRs. The invariance property is essential for overcoming this difficulty. With the invariance property, we may study the complete stability w.r.t. any delay parameter $${\tau _\ell }$$ through some simple and quite effective frequency-sweeping test. Furthermore, for any given delay combination $$\overrightarrow {{\tau ^ \#}}$$, we may accurately compute the number of unstable roots by using $$L$$ times frequency-sweeping tests in an appropriate iterative manner. Moreover, if the multiple-delay system (1.1) has CIRs, we may easily analyse the asymptotic behaviour of the CIRs w.r.t. each delay parameter from the frequency-sweeping curves (FSCs). As a consequence, we know if a stabilizing delay combination exists near $$\overrightarrow {{\tau ^ \#}}$$ and how to find it (if it exists). Finally, we propose an algorithm for a class of multiple-delay systems such that the stability may be automatically analysed by a single computer program. In our opinion, such an algorithm is easy to understand, to implement, and to apply, as illustrated by some examples taken from the open literature. The rest part of article is organized as follows. Some preliminaries are reviewed in Section 2. In Section 3, the complete stability of a linear time-delay system including multiple delays with single free delay is analysed. In Section 4, an approach is presented for computing the number of unstable roots for a linear system with any combination of multiple delays. Numerical examples are given in Section 5. In Section 6, an algorithm for automatic stability analysis for a class of multiple-delay systems is proposed. Finally, in Section 7, some concluding remarks end the article. Notations: In the sequel, $$\mathbb{R}$$ ($$\mathbb{R}_+$$) denotes the set of (positive) real numbers and $$\mathbb{C}$$ is the set of complex numbers; $$\mathbb{C}_-$$ and $$\mathbb{C}_+$$ denote respectively the left half-plane and right half-plane in $$\mathbb{C}$$; $$\mathbb{C}_0$$ is the imaginary axis and $$\partial \mathbb{D}$$ is the unit circle in $$\mathbb{C}$$; $$\mathbb{Z}$$, $$\mathbb{N}$$, and $$\mathbb{N}_+$$ are the sets of integers, non-negative integers and positive integers, respectively. $$\varepsilon$$ is a sufficiently small positive real number. Next, $$I$$ is the identity matrix and $$\overrightarrow 0 = (0, \ldots ,0)$$ of appropriate dimensions. For $$\gamma \in \mathbb{R}$$, $$\lceil \gamma \rceil$$ denotes the smallest integer greater than or equal to $$\gamma$$. Finally, $$\det (\cdot)$$ denotes the determinant of its argument. 2. Preliminaries The characteristic function of time-system (1.1) is   $f(\lambda ,\overrightarrow \tau ) = \det \left(\lambda I - A - \sum\limits_{\ell = 1}^L {{B_\ell }{e^{ - {\tau _\ell }\lambda }}} \right)\!,$ which is a quasipolynomial. For a non-zero delay vector $$\overrightarrow \tau$$, the characteristic equation $$f(\lambda ,\overrightarrow \tau ) = 0$$ has an infinite number of roots, i.e., the time-system (1.1) has an infinite number of characteristic roots. Time-delay system (1.1) is asymptotically stable if and only if all the characteristic roots lie in the open left half-plane $$\mathbb{C}_-$$, see e.g., Hale & Verduyn Lunel (1993) and Michiels & Niculescu (2014) for more comprehensive explanation. Throughout this article, we denote by $$NU(\overrightarrow \tau)$$ the number of characteristic roots in the open right half-plane $$\mathbb{C}_+$$ for system (1.1) in the presence of a delay vector $$\overrightarrow \tau$$. Clearly, the time-delay system (1.1) is asymptotically stable if and only if $$NU(\overrightarrow \tau)=0$$ and the system has no CIRs. If at $$\overrightarrow {{\tau ^*}} = (\tau _1^*, \ldots ,\tau _L^*)$$, $$f(j\omega^* , \overrightarrow {\tau^*}) = 0$$ ($${\omega ^*} \in \mathbb{R}_+$$, $$j$$ is the imaginary unit, and such a pair $$(\,j\omega^* , \overrightarrow {\tau^*})$$ is called a critical pair), then for all $$\overrightarrow \tau = \overrightarrow {{\tau ^*}} + ({k_1}, \ldots ,{k_L})\frac{{2\pi }}{{{\omega ^*}}}$$ ($${k_\ell } \in \mathbb{Z}$$ such that $$\tau _\ell ^* + {k_\ell }\frac{{2\pi }}{{{\omega ^*}}} \ge 0$$, $$\ell = 1, \ldots ,L$$), $$f(\,j\omega^*, \overrightarrow \tau) = 0$$. That is to say, a CIR $$j\omega^*$$ repeats infinitely many times as each $$\tau_{\ell}$$ increases, with a periodicity $$\frac{{2\pi }}{{\omega^*}}$$. Remark 2.1 It is a common assumption that $$\lambda=0$$ is not a characteristic root, otherwise the system can not be asymptotically stable for any $$\overrightarrow \tau$$. Owing to the conjugate symmetry of the spectrum, it suffices to consider only the CIRs with nonnegative imaginary parts. As mentioned in Section 1, the asymptotic behaviour analysis w.r.t. multiple free delay parameters is rather complicated (corresponding to an open mathematical problem). We will propose an iterative procedure: Analysing the asymptotic behaviour w.r.t. one free delay each time. 3. Complete stability w.r.t. one delay parameter In this section, we study the case where $$L-1$$ delay parameters among $$\tau_1$$, $$\ldots$$, $$\tau_L$$ are fixed and the remaining one is ‘free’. The objective is to study the complete stability of time-delay system (1) w.r.t. the remaining free delay parameter. For simplicity, we adopt the following notation $$\delta (\ell ) = ({\delta _1}(\ell ), \ldots ,{\delta _L}(\ell ))$$ with   ${\delta _i}(\ell ) = \left\{ \begin{array}{l} 0, {\kern 1pt} {\kern 1pt} {\rm{if}}{\kern 1pt} {\kern 1pt} i \ne \ell , \\ 1,{\kern 1pt} {\kern 1pt} {\rm{if}}{\kern 1pt} {\kern 1pt} i = \ell . \\ \end{array} \right.$ Then, $$\overrightarrow \tau = \sum\limits_{\ell = 1}^L {{\tau _\ell }\delta (\ell )}$$. Suppose $$\tau_{\chi}$$ ($$\chi \in \{ 1, \ldots ,L\}$$) is the free parameter and the other $$L-1$$ fixed delay parameters are $$\tau _\ell^ \#$$ ($$\ell \ne \chi$$). In this scenario, $$\overrightarrow \tau$$ can be expressed as $$\overrightarrow \tau = {\tau _\chi }\delta (\chi) + {F_\chi}$$, with $${F_\chi } = \sum\limits_{\ell \ne \chi} {{\tau _\ell^ \#}\delta (\ell)}$$. For instance, in the case $$L=4$$, suppose $${\tau _1} = \tau _1^ \#$$, $${\tau _3} = \tau _3^ \#$$, and $${\tau _4} = \tau _4^ \#$$ are fixed, while $$\tau_2$$ is the free delay parameter. We may express $$\overrightarrow \tau = (\tau _1^ \#, {\tau _2}, \tau _3^ \#, \tau _4^ \#)$$ as $$\overrightarrow \tau = {\tau _2}\delta (2) + {F_2}$$ with $${F_2} = (\tau _1^ \#, 0, \tau _3^ \#, \tau _4^ \#)$$ and $$\delta (2) = (0,1,0,0)$$. We study the stability of time-delay system (1.1) w.r.t. $$\tau_{\chi}$$ along the whole non-negative $$\tau_{\chi}$$-axis, i.e., the complete stability problem w.r.t. the delay parameter $$\tau_{\chi}$$. More precisely, we will keep track of the number of unstable roots, denoted by $$N{U_{{F_\chi}}}({\tau _\chi})$$, for $${\tau _\chi} \in [0, + \infty )$$. In this context, the characteristic function $$f(\lambda ,\overrightarrow \tau )$$ can be rewritten as   \begin{eqnarray}\label{ftauFsecondform} f(\lambda ,{\tau _\chi},{F_\chi}) = \sum\limits_{i = 0}^{{q_\chi}} {{{\widetilde a}_{\chi i}}(\lambda ){e^{ - i{\tau _\chi}\lambda }}}, \end{eqnarray} (3.1) where $$\widetilde a_{\chi i}$$ ($$i =0, \ldots ,{q_\chi}$$) are polynomials in $$\lambda$$ and $${e^{ - \tau _\ell ^\# \lambda }}$$ ($$\ell \in \{ 1, \ldots ,L\}$$, $$\ell \ne \chi$$) with real coefficients. Furthermore, the expression (3.1) can be viewed as a polynomial of $${e^{ - {\tau _\chi}\lambda }}$$ where $${\widetilde a_{\chi 0}}(\lambda ), \ldots ,{\widetilde a_{\chi {q_\chi}}}(\lambda )$$ are treated as the coefficient functions. We introduce the following two-variable polynomial expression of $$f(\lambda ,{\tau _\chi},{F_\chi})$$:   \begin{eqnarray}\label{ptauF} p(\lambda ,{z_\chi},{F_\chi}) = \sum\limits_{i = 0}^{{q_\chi}} {{{\widetilde a}_{\chi i}}(\lambda )z_{^\chi}^i}, \end{eqnarray} (3.2) where   ${z_\chi} = {e^{ - {\tau _\chi}\lambda }}.$ Remark 3.1 The coefficient functions $${\widetilde a_{\chi i}}(\lambda)$$ ($$i=1, \ldots, q_{\chi}$$) of (3.1) as well as (3.2) are independent of $$\tau_{\chi}$$. The detection of the CIRs and the CDs for $$f(\lambda ,{\tau _\chi},{F_\chi}) = 0$$ amounts in detecting the critical pairs $$(\lambda, z_\chi)$$ ($$\lambda \in \mathbb{C}_0$$ and $$z_\chi \in \partial \mathbb{D}$$) for $$p(\lambda ,{z_\chi},{F_\chi})=0$$. Without any loss of generality, suppose that there exist $$u$$ critical pairs for $$p(\lambda, z_\chi, F_\chi)=0$$ denoted by $$({\lambda _{\chi,0}} = j{\omega _{\chi,0}},{z_{\chi,0}})$$, $$({\lambda _{\chi,1}} = j{\omega _{\chi,1}},{z_{\chi,1}})$$, $$\ldots$$, $$({\lambda _{\chi,u-1}} = j{\omega _{\chi,u-1}},{z_{\chi,u-1}})$$ with $$0 < {\omega _{\chi,0}} \le {\omega _{\chi,1}} \le \cdots \le {\omega _{\chi, u - 1}}$$. Once all the critical pairs $$({\lambda _{\chi,\alpha} },{z_{\chi,\alpha} }),\alpha = 0, \ldots ,u - 1$$, are found, all the critical pairs $$(\lambda,\tau_\chi)$$ for $$f(\lambda ,{\tau _\chi},{F_\chi}) = 0$$ can be obtained: For each CIR $$\lambda_ {\chi, \alpha}$$, the corresponding (infinitely many) CDs are given by $${\tau _{\chi, \alpha, k}} \buildrel {\it{\Delta}} \over = {\tau _{\chi, \alpha, 0}} + \frac{{2k\pi }}{{{\omega _{\chi, \alpha} }}}, k \in \mathbb{N}, {\tau _{\chi, \alpha, 0}} \buildrel {\it{\Delta}} \over = \min \{ \tau \ge 0:{e^{ - \tau {\lambda _{\chi, \alpha} }}} = {z_{\chi, \alpha} }\}$$. The pairs $$({\lambda _{\chi, \alpha} },{\tau _{\chi, \alpha ,k}}),k \in \mathbb{N}$$, define a set of critical pairs associated with $$({\lambda _{\chi, \alpha} },{z_{\chi, \alpha} })$$. Remark 3.2 According to the root continuity argument, if the time-delay system (1.1) has no CIRs for all $${\tau _\chi } \ge 0$$, $$NU({\tau _\chi}\delta (\chi) + {F_\chi}) = NU({F_\chi})$$ for all $$\tau_{\chi}>0$$. 3.1. Asymptotic behaviour of a critical pair An essential step for the stability analysis is to address the asymptotic behaviour of a critical pair $$({\lambda _{\chi, \alpha }}, {\tau _{\chi, \alpha, k}})$$. Since the CIR may be multiple and split as $$\tau_{\chi}$$ increases near a positive CD, we now introduce a general notation to describe the asymptotic behaviour for a critical pair $$({\lambda _{\chi,\alpha }} = j{\omega _{\chi,\alpha }},{\tau _{\chi,\alpha ,k}})$$ (provided $$F_{\chi}$$ is given)   $\begin{array}{l} {\it{\Delta}} NU_{{\lambda _{\chi,\alpha }}}^{{F_\chi}}({\tau _{\chi,\alpha ,k}}) \buildrel {\it{\Delta}} \over = NU(({\tau _{\chi,\alpha ,k}} + \varepsilon )\delta (\chi) + {F_\chi }) - NU(({\tau _{\chi,\alpha ,k}} - \varepsilon )\delta (\chi) + {F_\chi}), \end{array}$ which stands for the change in $$NU(\overrightarrow \tau)$$ caused by the variation of the CIR $${\lambda _{\chi ,\alpha }}$$ as $$\tau_{\chi}$$ increases from $${\tau _{\chi,\alpha ,k}} - \varepsilon$$ to $${\tau _{\chi,\alpha ,k}} + \varepsilon$$, i.e., the asymptotic behaviour of the CIR $${\lambda _{\chi,\alpha }}$$ at a positive CD $${\tau _{\chi,\alpha ,k}}$$. The value of $${\it{\Delta}} NU_{{\lambda _{\chi,\alpha }}}^{{F_\chi }}({\tau _{\chi,\alpha ,k}})$$ can be calculated by invoking the associated Puiseux series. One may refer to Chapter 4 of Li et al. (2015) for a general algorithm for invoking the Puiseux series. At this stage, it is interesting to see that the work of invoking the Puiseux series may be replaced by a graphical (frequency-sweeping) test, which will be presented later in this article. In order to analyse the complete stability, we need to specifically address the situation as $$\tau_{\chi}$$ increases from $$0$$ to $$+\varepsilon$$. In other words, we need to know $$NU({F_\chi} + \varepsilon \delta (\chi))$$. Similarly to Theorem 5.1 in Li et al. (2015), we have the following: Theorem 3.1 If the time-delay system (1.1) has no CIRs when $$\tau_{\chi}=0$$, $$NU({F_\chi} + \varepsilon \delta (\chi))= NU({F_\chi})$$. Otherwise, $$NU({F_\chi} + \varepsilon \delta (\chi)) - NU({F_\chi})$$ equals to the number of the values in $$\mathbb{C}_+$$ of the Puiseux series for all the corresponding CIRs when $$\tau_{\chi}=0$$ with $${\it{\Delta}} \tau_{\chi} = + \varepsilon$$. Proof. If $$F_\chi$$ is zero vector, the system has finitely many characteristic roots when $$\tau_{\chi}=0$$. As $$\tau_{\chi}$$ increases to a sufficiently small positive number $$+\varepsilon$$, infinitely many new characteristic roots appear at the far left of the complex plane while the original (finitely many) characteristic roots change continuously w.r.t. $$\tau_{\chi}$$. If $$F_\chi$$ is a non-zero vector, all the (infinitely many) characteristic roots change continuously as $$\tau_{\chi}$$ increases from $$0$$ to $$+\varepsilon$$. Thus, if the system has no CIRs when $$\tau_{\chi}=0$$, it is trivial to know that $$NU({F_\chi} + \varepsilon \delta (\chi)) - NU({F_\chi}) = 0$$. Otherwise, the value of $$NU({F_\chi} + \varepsilon \delta (\chi)) - NU({F_\chi})$$ is determined by the asymptotic behaviour of the corresponding CIRs, which can be explicitly analysed through the Puiseux series (as stated in Theorem 3.1). □ 3.2. Frequency-sweeping curves (FSCs) Under a given $$F_{\chi}$$, the FSCs for time-delay system (1.1) are obtained by the following procedure. Frequency-sweeping curves (FSCs): Sweep $$\omega \ge 0$$ and for each $$\lambda = j\omega$$ we have $$q_{\chi}$$ solutions of $$z_{\chi}$$ such that $$p(\lambda ,{z_\chi},{F_\chi}) = 0$$ (denoted by $${z_{\chi,1}}(\,j\omega ), \ldots ,{z_{\chi,{q_\chi}}}(\,j\omega )$$). In this way, we obtain $$q_{\chi}$$ FSCs $${{\it{\Gamma}} _{\chi,i}}(\omega ):$$$$\left| {{z_{\chi,i}}(\,j\omega )} \right|$$ vs. $$\omega$$, $$i = 1, \ldots ,{q_\chi}$$. For simplicity, we denote by $$\Im _1$$ the line parallel to the abscissa axis with ordinate equal to 1. If $$({\lambda _{\chi,\alpha }} = j{\omega _{\chi,\alpha }},{\tau _{\chi,\alpha ,k}})$$ are a set of critical pairs, some FSCs intersect $${\Im _1}$$ at $$\omega = {\omega _{\chi, \alpha}}$$. It is straightforward to see that all the CIRs and the CDs can be detected by the FSCs (more precisely, the intersection of the FSCs and the line $$\Im _1$$). We now address the asymptotic behaviour of the FSCs (such an idea was proposed and largely discussed in Li et al., 2015). For a set of critical pairs $$({\lambda _{\chi,\alpha} =j \omega_{\chi,\alpha}},{\tau _{\chi, \alpha ,k}})$$ (as usually assumed, $$\lambda _{\chi,\alpha} \ne 0$$), there must exist some FSCs such that $${z_{\chi,i}}(\,j{\omega _{\chi,\alpha }}) = {z_{\chi ,\alpha }} = {e^{ - {\tau _{\chi,\alpha ,0}}{\lambda _{\chi,\alpha }}}}$$ intersecting $${\Im _1}$$ when $$\omega = {\omega _{\chi ,\alpha }}$$. Among these FSCs, we denote the number of those above $$\Im _1$$ when $$\omega = {\omega _{\chi,\alpha }} + \varepsilon$$ ($$\omega = {\omega _{\chi,\alpha }} - \varepsilon$$) by $$NF_{{z_{\chi,\alpha }}}^{{F_\chi}}({\omega _{\chi,\alpha }} + \varepsilon )$$ ($$NF_{{z_{\chi,\alpha }}}^{{F_\chi}}({\omega _{\chi ,\alpha }} - \varepsilon )$$). We now introduce the notation $${\it{\Delta}} NF_{{z_{\chi,\alpha }}}^{{F_\chi}}({\omega _{\chi,\alpha }})$$ as defined below   ${\it{\Delta}} NF_{{z_{\chi,\alpha }}}^{{F_\chi}}({\omega _{\chi,\alpha }}) \buildrel {\it{\Delta}} \over = NF_{{z_{\chi,\alpha }}}^{{F_\chi}}({\omega _{\chi,\alpha }} + \varepsilon ) - NF_{{z_{\chi,\alpha }}}^{{F_\chi}}({\omega _{\chi,\alpha }} - \varepsilon ).$ Remark 3.3 It is a useful property that $${\it{\Delta}} NF_{{z_{\chi,\alpha }}}^{{F_\chi}}({\omega _{\chi ,\alpha }})$$ is invariant w.r.t. different CDs. Moreover, the value of $${\it{\Delta}} NF_{{z_{\chi,\alpha }}}^{{F_\chi}}({\omega _{\chi ,\alpha }})$$ may be easily obtained graphically (by observing how the corresponding FSCs intersect $$\Im _1$$): it equals to the number change of the corresponding FSCs above $${\Im _1}$$ when $$\omega = {\omega _{\chi,\alpha }} + \varepsilon$$ and $$\omega = {\omega _{\chi,\alpha }} - \varepsilon$$. For instance, for a set of critical pairs $$({\lambda _{\chi,\alpha} =j \omega_{\chi,\alpha}},{\tau _{\chi, \alpha ,k}})$$, suppose that there exists one and only one FSC such that $${z_{\chi,i}}(\,j{\omega _{\chi,\alpha }}) = {z_{\chi ,\alpha }} = {e^{ - {\tau _{\chi,\alpha ,0}}{\lambda _{\chi,\alpha }}}}$$ and that as $$\omega$$ increases near $${\omega _{\chi ,\alpha }}$$ the FSC intersects $${\Im _1}$$ from above to below. Then, $$NF_{{z_{\chi,\alpha }}}^{{F_\chi}}({\omega _{\chi,\alpha }} + \varepsilon ) = 0$$, $$NF_{{z_{\chi,\alpha }}}^{{F_\chi}}({\omega _{\chi,\alpha }} - \varepsilon )=1$$ and hence $${\it{\Delta}} NF_{{z_{\chi,\alpha }}}^{{F_\chi}}({\omega _{\chi,\alpha }}) \buildrel {\it{\Delta}} \over = NF_{{z_{\chi,\alpha }}}^{{F_\chi}}({\omega _{\chi,\alpha }} + \varepsilon ) - NF_{{z_{\chi,\alpha }}}^{{F_\chi}}({\omega _{\chi,\alpha }} - \varepsilon ) = 0 - 1 =-1$$. This is the case to be seen in Example 5.1. 3.3. Invariance property We are now in a position to present the result concerning the invariance property of the asymptotic behaviour w.r.t. the free delay parameter $$\tau_\chi$$ for time-delay system (1.1). Theorem 3.2 For a CIR $${\lambda _{\chi,\alpha }}$$ of time-delay system (1.1) with given $${F_\chi}$$, $${\it{\Delta}} NU_{{\lambda _{\chi,\alpha }}}^{{F_\chi }}({\tau _{\chi,\alpha ,k}})$$ is a constant $${\it{\Delta}} NF_{{z_{\chi ,\alpha }}}^{{F_\chi}}({\omega _{\chi,\alpha }})$$ for all $${\tau _{\chi,\alpha, k}} > 0$$. The proof is given in the Appendix. 3.4. Ultimate stability With the above invariance property, we now address the ultimate stability issue (i.e., the system stability as $${\tau _\chi} \to \infty$$). Property 3.1 For all $$i= 1, \dots, q_\chi$$, it follows that $$\left| {{z_{\chi, i}}(\,j\omega )} \right| \to \infty$$ as $$\omega \to \infty$$. Proof. As the time-system (1.1) is of retarded type and $$\left| {{e^\iota }} \right| = 1$$ for any purely imaginary number $$\iota$$, we have   $\mathop {\lim }\limits_{\omega \to \infty } \frac{{{{\widetilde a}_{\chi i}}(\,j\omega )}}{{{{\widetilde a}_{\chi 0}}(\,j\omega )}} = 0,i = 1, \ldots ,{q_\chi }.$ The proof is now complete according to (3.2).□ A critical frequency $${\omega _{\chi,\alpha }}$$ is called a crossing (touching) frequency for an FSC $${\it{\Gamma}} _{\chi,i} (\omega)$$, if $${\it{\Gamma}} _{\chi,i} (\omega)$$ crosses (touches without crossing) $${\Im _1}$$ as $$\omega$$ increases near $$\omega = \omega_{\chi,\alpha}$$. With the notions above, we have the following results. Theorem 3.3 $${{\it{\Gamma}} _{\chi ,i}}(\omega ),i = 1, \ldots ,{q_\chi }$$, have a crossing frequency, there exists a $$\tau_{\chi}^*$$ such that the time-delay system (1.1) is unstable for all $$\tau_{\chi} > \tau_{\chi}^ *$$ and $$\mathop {\lim }\limits_{{\tau _\chi } \to \infty } NU({\tau _\chi }\delta (\chi ) + {F_\chi }) = \infty$$. Theorem 3.4 A time-delay system (1.1) with given $${F_\chi}$$ must belong to the following three types: Type 1: The system has a crossing frequency and $$\mathop {\lim }\limits_{{\tau _\chi} \to \infty } NU({\tau _\chi}\delta (\chi) + {F_\chi}) = \infty$$; Type 2: The system has neither crossing frequencies nor touching frequencies and $$NU({\tau _\chi}\delta (\chi) + {F_\chi}) = NU({F_\chi})$$ for all $$\tau _ \chi >0$$; Type 3: The system has touching frequencies but no crossing frequencies and $$NU({\tau _\chi}\delta (\chi) + {F_\chi})$$ is a constant for all $$\tau_\chi$$ other than the CDs. Following the line of Theorems 9.1 and 9.2 in Li et al. (2015), we may prove Theorems 3.3 and 3.4 in light of Property 3.1. 3.5. Explicit computation of number of unstable roots Assume that $$NU({F_\chi})$$ is known (this value can be obtained by the procedure to be given in the next section). We now show that $$NU({\tau _\chi}\delta (\chi) + {F_\chi})$$ can be expressed as an explicit function of $$\tau_{\chi}$$. For each CIR $${\lambda _{\chi,\alpha }}$$, we may choose any positive CD $${\tau _{\chi,\alpha ,k}}$$ to compute $${\it{\Delta}} NU_{{\lambda _{\chi ,\alpha }}}^{{F_\chi}}({\tau _{\chi,\alpha ,k}})$$ (the value is denoted by $${U_{{\lambda _{\chi,\alpha }}}}$$), through invoking the Puiseux series. Alternatively, we may directly have that $${U_{{\lambda _{\chi,\alpha }}}} = {\it{\Delta}} NF_{{z_{\chi,\alpha }}}^{{F_\chi}}({\omega _{\chi,\alpha }})$$ from the FSCs, according to Theorem 3.2. In the light of the invariance property, the explicit expression of $$NU({\tau _\chi}\delta (\chi) + {F_\chi})$$ is as follows: Theorem 3.5 For time-delay system (1.1) and for any $$\tau_{\chi}$$ which is not a CD, it follows that   $NU({\tau _\chi}\delta (\chi) + {F_\chi}) = NU(\varepsilon \delta (\chi) + {F_\chi}) + \sum\limits_{\alpha = 0}^{u - 1} {N{U_{\chi ,\alpha }}({\tau _\chi})} ,$   ${\rm where}\ N{U_{\chi,\alpha }}({\tau _\chi}) = \left\{ {\begin{array}{*{20}{l}} {0,\tau < {\tau _{\chi,\alpha ,0}},}\\ {2{U_{{\lambda _{\chi,\alpha }}}}\left\lceil {\frac{{{\tau _\chi} - {\tau _{\chi,\alpha ,0}}}}{{2\pi /{\omega _{\chi,\alpha }}}}} \right\rceil ,\tau > {\tau _{\chi,\alpha ,0}},} \end{array}} \right.{\rm{if}}\,{\tau _{\chi,\alpha ,0}} \ne 0,$   $\left\{ {\begin{array}{*{20}{l}} {0,\tau < {\tau _{\chi,\alpha ,1}},}\\ {2{U_{{\lambda _{\chi,\alpha }}}}\left\lceil {\frac{{{\tau _\chi} - {\tau _{\chi,\alpha ,1}}}}{{2\pi /{\omega _{\chi,\alpha }}}}} \right\rceil ,\tau > {\tau _{\chi,\alpha ,1}},} \end{array}} \right.{\rm{if}}\,{\tau _{\chi,\alpha ,0}} = 0.$ We can now systematically solve the complete stability w.r.t. the free delay parameter $$\tau_\chi$$. The system is asymptotically stable if and only if $$\tau_\chi$$ lies in the interval(s) with $$NU({\tau _\chi}\delta (\chi) + {F_\chi}) =0$$ excluding the CDs. The ultimate stability is known according to Theorems 3.3 and 3.4. 4. Stability analysis for any given delay vector We now study the stability of time-delay system (1.1) for any given delay vector $$\overrightarrow {{\tau ^ \#}} = (\tau _1^ \#, \ldots ,\tau _L^ \#)$$. Based on the results of Section 3, an iterative frequency-sweeping approach is given below. An Iterative Frequency-Sweeping Approach Step 0: Set $$\chi=1$$ and $${F_1} = (0, \ldots ,0)$$. Step 1: Compute $$NU({F_\chi } + \varepsilon \delta (\chi ))$$ by Theorem 3.1. Step 2: Generate the FSCs corresponding to $$p(\lambda ,{\tau _\chi},{F_\chi}) = 0$$, denoted by $${{\it{\Gamma}} _{\chi ,i}}(\omega ),i = 1, \ldots ,{q_\chi }$$. Following the results in Section 3, we can know $$NU({F_\chi } + \tau _\chi ^ \# \delta (\chi ))$$ and detect the CIRs (if any!). Step 3: If $$\chi < L$$, let $$\chi = \chi + 1$$ and $${F_\chi } = {F_{\chi - 1}} + \tau _{\chi - 1}^ \# \delta (\chi - 1)$$. Return to Step 1. (After running Steps 1–3 for the last time (i.e., when $$\chi = L$$), we know the value of $$NU(\overrightarrow {{\tau ^ \#}} )$$ and whether the system (1.1) has CIRs when $$\overrightarrow \tau = \overrightarrow {{\tau ^ \#}}$$). Step 4: If the system has no CIRs when $$\overrightarrow \tau = \overrightarrow {{\tau ^ \#}}$$, skip to Step 5. Otherwise, generate the FSCs corresponding to $$p(\lambda ,{\tau _\ell},{F_\ell ^ \#}) = 0$$ with $${F_\ell ^ \#} = \sum\limits_{\kappa \ne \ell} {{\tau _\kappa^ \#}\delta (\kappa)}$$, denoted by $${\it{\Gamma}} _{\ell,i}^ \# (\omega ),i = 1', \ldots ,{q_\ell}$$, $$\ell = 1, \ldots ,L$$.4 (We may analyse the asymptotic behaviour of the CIRs when $$\overrightarrow \tau = \overrightarrow {{\tau ^ \#}}$$ w.r.t. each delay parameter $$\tau_{\ell}$$, $$\ell=1,\ldots,L$$, from the FSCs $${\it{\Gamma}} _{\ell,i}^ \# (\omega ),i = 1, \ldots ,{q_\ell}$$. To be more precise, we can determine if increasing or decreasing a delay parameter $$\tau_{\ell}$$ sufficiently near $$\tau _\ell ^ \#$$ may stabilize the system, according to Theorem 3.2. As a result, we know if there exists a $$\overrightarrow \tau$$ in a small neighbourhood of $$\overrightarrow {{\tau ^ \#}}$$, at which the system is asymptotically stable.) Step 5: The procedure ends. Proposition 4.1 For time-delay system (1.1) and for any given delay vector $$\overrightarrow {{\tau ^ \#}} = (\tau _1^ \#, \ldots ,\tau _L^ \#)$$, the iterative frequency-sweeping approach allows to calculate $$NU(\overrightarrow \tau )$$ with $$\overrightarrow \tau = \overrightarrow {{\tau ^ \#}}$$. Proof. It is worth to mention that the invariance property guarantees that the stability characterization is complete w.r.t. any free delay parameter. After the first iteration (Steps 1–3 for the first time), we know the stability property for $$\tau = (\tau _1^\# ,0, \ldots ,0)$$. Next, after the second iteration (Steps 1–3 for the second time), we know the stability property for $$\tau = (\tau _1^\# ,\tau _2^\# ,0, \ldots ,0)$$. In this iterative manner, we are able to examine the stability properties for the delay vectors $$(\tau _1^\# ,0, \ldots ,0)$$, $$(\tau _1^\# ,\tau _2^\# ,0, \ldots ,0)$$, $$(\tau _1^\# ,\tau _2^\# ,\tau _3^\# ,0, \ldots ,0)$$, $$\ldots$$ The iterative frequency-sweeping approach requires totally $$L$$ iterations (one per independent delay parameter). Finally, as expected, after the last iteration (i.e., Steps 1–3 for the $$L$$-th time), we can determine the stability property for $$\tau = \overrightarrow {{\tau ^ \#}} = (\tau _1^ \#, \ldots ,\tau _L^ \#)$$.□ Remark 4.1 The above asymptotic behaviour analysis (by means of Step 4) is simple and effective. The effect of each delay parameter’s asymptotic behaviour can be explicitly known. It should be pointed out that the asymptotic behaviour analysis w.r.t. all the $$L$$ delay parameters (simultaneously) is very difficult (theoretically, it corresponds to some Puiseux series for $${\it{\Delta}} \lambda$$ w.r.t. $${\it{\Delta}} {\tau _1}, \ldots ,{\it{\Delta}} {\tau _L}$$), especially when $$L$$ is large. However, as mentioned in Section 1, to the best of the authors’ knowledge, no general tool has been reported for the multiple-parameter asymptotic behaviour analysis so far. 5. Numerical examples We now give some illustrative examples concerning the properties as well as the results presented in this article. Example 5.1 (Invariance property) Consider a time-delay system involving two delay parameters $$\tau_1$$ and $$\tau_2$$ with the characteristic function   $f(\lambda ,\overrightarrow \tau ) = {e^{ - ({\tau _1} + {\tau _2})\lambda }} - ({\lambda ^6} - {\lambda ^4} + {\lambda ^2}){e^{ - {\tau _2}\lambda }} - ({\lambda ^{10}} - {\lambda ^8} + {\lambda ^6}){e^{ - {\tau _1}\lambda }} + {\lambda ^{12}}.$ Suppose $$\tau_2$$ is fixed as $$\tau_2^ \#= 2 \pi$$ and $$\tau_1$$ is the free delay parameter. We here illustrate the invariance property. The corresponding characteristic function can be expressed as $$f(\lambda ,{\tau _1},{F_1}) = {\widetilde a_{10}}(\lambda ) + {\widetilde a_{11}}(\lambda ){e^{ - {\tau _1}\lambda }}$$ with $${\widetilde a_{10}}(\lambda ) = {\lambda ^{12}} - ({\lambda ^6} - {\lambda ^4} + {\lambda ^2}){e^{ - {\tau _2^ \#}\lambda }}$$ and $${\widetilde a_{11}}(\lambda )={e^{ - {\tau _2^ \#}\lambda }} - ({\lambda ^{10}} - {\lambda ^8} + {\lambda ^6})$$. At $$\tau_1 = (2k + 1)\pi$$, $$\lambda=j$$ is a CIR and, in particular, $$\lambda=j$$ is a triple CIR at $$\tau_1 = \pi$$ (it is simple at all $$\tau_1 = (2k + 1)\pi, k \in \mathbb{N}_+$$). According to Theorem 3.2, $${\it{\Delta}} NU_j^{{F_1}}((2k + 1)\pi ) = {\it{\Delta}} NF_{ - 1}^{{F_1}}(1) = - 1$$ for all $$k \in \mathbb{N}$$, where $${\it{\Delta}} NF_{ - 1}^{{F_1}}(1)$$ can be obtained from the FSC $${{\it{\Gamma}} _{1,1}}$$ shown in Fig. 1. To verify the result, we give the Puiseux series for the critical pair $$(\,j,\pi)$$:   ${\it{\Delta}} \lambda = (0.2978 + 0.1019j){({\it{\Delta}} {\tau _1})^{\frac{1}{3}}} + o({({\it{\Delta}} {\tau _1})^{\frac{1}{3}}})$ and the Taylor series for the critical pairs $$(\,j,(2k + 1)\pi )$$ (due to the limited space, we only give the results for $$k=1$$ and $$k=2$$):   ${{\it{\Delta}} \lambda = - 0.1592j{\it{\Delta}} {\tau _1} + 0.0253j{{({\it{\Delta}} {\tau _1})}^2} + ( - 0.0113 + 0.0132j){{({\it{\Delta}} {\tau _1})}^3} + o({{({\it{\Delta}} {\tau _1})}^3}),k = 1,}$   $\begin{array}{*{20}{l}} {{\it{\Delta}} \lambda = - 0.0796j{\it{\Delta}} {\tau _1} + 0.0063j{{({\it{\Delta}} {\tau _1})}^2} + ( - 0.0007 + 0.0006j){{({\it{\Delta}} {\tau _1})}^3} + o({{({\it{\Delta}} {\tau _1})}^3}),k = 2.} \end{array}$ It is worth mentioning that our result is consistent with the above series analysis. □ Fig. 1. View largeDownload slide FSC for Example 5.1. Fig. 1. View largeDownload slide FSC for Example 5.1. Example 5.2 (Complete stability) Consider a time-delay system with the characteristic function   $f(\lambda ,\overrightarrow \tau ) = {\lambda ^4} + 2{\lambda ^2} + 3{e^{ - {\tau _1}\lambda }} - 3{e^{ - {\tau _2}\lambda }} + {e^{ - ({\tau _1} + {\tau _2})\lambda }}.$ We study the complete stability w.r.t. $$\tau_2$$ when $$\tau_1$$ is fixed as $$\tau_1^ \#=0.5$$. First, when $$\overrightarrow \tau = \overrightarrow 0$$, the system has four CIRs (both $$+j$$ and $$-j$$ are double CIRs). Owing to the conjugate symmetry of the spectrum, it suffices to consider the asymptotic behaviour of the CIR $$+j$$. Applying the method in Chapter 4 of Li et al. (2015), we have that the asymptotic behaviour of the CIR $$+j$$ w.r.t. $$\tau_1$$ at $$\overrightarrow \tau = \overrightarrow 0$$ corresponds to the Puiseux series $${\it{\Delta}} \lambda = {( - j{\it{\Delta}} {\tau _1})^{\frac{1}{2}}} + o({({\it{\Delta}} {\tau _1})^{\frac{1}{2}}})$$. Thus, in the light of Theorem 3.1, $$NU(( + \varepsilon ,0)) = + 2$$. Next, through the FSC $${{\it{\Gamma}} _{1,1}}$$ with $${F_1} = (0,0)$$ (Fig. 2a), we have that $$NU((0.5,0)) = + 2$$ (two sets of critical pairs are obtained from Fig. 2(a): $$({\lambda _{1,0}} = j,{\tau _{1,0,k}} = 2k\pi )$$ and $$({\lambda _{1,1}} = 1.9566j,{\tau _{1,1,k}} = 1.6056 + \frac{{2k\pi }}{{1.9566}})$$). Finally, we generate the FSC $${{\it{\Gamma}} _{2,1}}$$ with $${F_2} = (0.5,0)$$ (Fig. 2(b)) and we have that the system has only one set of critical pairs: $$({\lambda _{2,0}} = j,{\tau _{2,0,k}} = {\rm{0}}{\rm{.9443}} + 2k\pi )$$. According to Theorem 3.2, each time $$\tau_2$$ increases near $${\tau _{2,0,k}}$$, $${\it{\Delta}} NU_{{\lambda _{2,0}}}^{{F_2}}({\tau _{2,0,k}}) = 0$$. Thus, $$NU((0.5,{\tau _2})) = + 2$$ for all $$\tau_2 \ge 0$$ except at $$\tau_2 = {\tau _{2,0,k}}$$. Finally, it is worth mentioning that the above analysis is consistent with the SCS, shown in Fig. 3. It is interesting to see that, in the three regions partitioned by the SCS in Fig. 3, the system has the same number of unstable roots.□ Fig. 2. View largeDownload slide FSCs for Example 5.2. (a) $${{\it{\Gamma}} _{1,1}}$$. (b) $${{\it{\Gamma}} _{2,1}}$$. Fig. 2. View largeDownload slide FSCs for Example 5.2. (a) $${{\it{\Gamma}} _{1,1}}$$. (b) $${{\it{\Gamma}} _{2,1}}$$. Fig. 3. View largeDownload slide Stability crossing set for Example 5.2. Fig. 3. View largeDownload slide Stability crossing set for Example 5.2. Example 5.3 (Stability analysis for a given delay vector) Consider the time-delay system in VII.B of Gu & Naghnaeian (2011) including three delays, with the characteristic function   $f(\lambda ,\overrightarrow \tau ) = {\lambda ^3} + 3\lambda + 7 + ({\lambda ^2} + 3\lambda + 1){e^{ - {\tau _1}\lambda }} + (4\lambda + 3){e^{ - {\tau _2}\lambda }} + ({\lambda ^2} + \lambda + 0.1){e^{ - {\tau _3}\lambda }}.$ In the sequel, our objective is to analyse the stability for a given delay vector $$\overrightarrow {{\tau ^ \#}} = (0.3041,1.7,0.0504)$$. As seen from Fig. 14 in Gu & Naghnaeian (2011), this $$\overrightarrow {{\tau ^ \#}}$$ corresponds to an intersecting point of the SCS. We now apply the procedure proposed in Section 4. The first iteration (for $$\chi = 1$$, after Step 0) $$\Rightarrow$$ Step 1: $$NU(( + \varepsilon ,0,0)) = 0$$ according to Theorem 3.1 (when $$\overrightarrow \tau = \overrightarrow 0$$, the characteristic roots are $$- 0.4457{\rm{ }} \pm {\rm{ }}3.1326j$$ and $$-1.1087$$). Step 2: We generate the FSC $${{\it{\Gamma}} _{1,1}}$$ (Fig. 4(a)), from which we obtain two sets of critical pairs: $$({\lambda _{1,0}} = 2.2819j,{\tau _{1,0,k}} = 1.9052 + \frac{{2k\pi }}{{2.2819}})$$ and $$({\lambda _{1,1}} = 3.5160j,{\tau _{1,1,k}} = 0.2756 + \frac{{2k\pi }}{{3.5160}})$$. We have that $$NU(({\rm{0}}{\rm{.3041}},0,0)) = +2$$ and the system has no CIRs at $$\overrightarrow \tau = (0.3041,0,0)$$. Step 3: Let $$\chi=2$$ and go to Step 1. The second iteration (for $$\chi = 2$$) $$\Rightarrow$$ Step 1: $$NU(({\rm{0}}{\rm{.3041}}, + \varepsilon ,0)) = +2$$ according to Theorem 3.1. Step 2: We generate the FSC $${{\it{\Gamma}} _{2,1}}$$ (Fig. 4(b)), from which we obtain two sets of critical pairs: $$({\lambda _{2,0}} = 1.9073j,{\tau _{2,0,k}} = 1.7471 + \frac{{2k\pi }}{{1.9073}})$$ and $$({\lambda _{2,1}} = 3.5133j,{\tau _{2,1,k}} = 1.7577 + \frac{{2k\pi }}{{3.5133}})$$. We have that $$NU(({\rm{0}}.{\rm{3041}},1.7,0)) = +2$$ and the system has no CIRs at $$\overrightarrow \tau = (0.3041,1.7,0)$$. Step 3: Let $$\chi=3$$ and go to Step 1. The third iteration (for $$\chi = 3$$) $$\Rightarrow$$ Step 1: $$NU(({\rm{0}}.{\rm{3041}},1.7, + \varepsilon )) = +2$$ according to Theorem 3.1. Step 2: We generate the FSC $${{\it{\Gamma}} _{3,1}}$$ (Fig. 5(a)), from which we obtain two sets of critical pairs: $$({\lambda _{3,0}} = 1.9594j,{\tau _{3,0,k}} = 0.0504 + \frac{{2k\pi }}{{1.9594}})$$ and $$({\lambda {3,1}} = 3.6118j,{\tau _{3,1,k}} = 0.0504 + \frac{{2k\pi }}{{3.6118}})$$. In light of $${{\it{\Gamma}} _{3,1}}$$, each time $$\tau_3$$ increases from $${\tau _{3,0,k}} - \varepsilon$$ to $${\tau _{3,0,k}} + \varepsilon$$ ($${\tau _{3,1,k}} - \varepsilon$$ to $${\tau _{3,1,k}} + \varepsilon$$), $${\lambda _{3,0}}$$ ($${\lambda _{3,1}}$$) rosses $$\mathbb{C}_0$$ from $$\mathbb{C}_+$$ to $$\mathbb{C}_-$$ (from $$\mathbb{C}_-$$ to $$\mathbb{C}_+$$). The complete stability w.r.t. $$\tau_3$$ can be studied. We may accurately calculate $$NU((\tau _1^\#,\tau _2^ \#,{\tau _3}))$$ as a function of $$\tau _3$$, as shown in Fig. 5(b). In particular, near $$\tau_3^ \#=0.0504$$, the variation of $$\tau_3$$ does not change $$NU((\tau _1^ \#, \tau _2^ \#, {\tau _3}))$$ since it causes opposite crossing directions (w.r.t. the imaginary axis $$\mathbb{C}_0$$) for $${\lambda _{3,0}}$$ and $${\lambda _{3,1}}$$. From the above iterations, we know that at $$\overrightarrow {{\tau ^ \#}}$$ the system has two CIRs ($$1.9594j$$ and $$3.6118j$$), without characteristic roots in $$\mathbb{C}_+$$. We next study the asymptotic behaviour of the CIRs w.r.t. $$\tau_1$$ and $$\tau_2$$, respectively, in order to see if there is a stabilizing point near $$\overrightarrow {{\tau ^ \#}}$$. Step 4 $$\Rightarrow$$ The asymptotic behaviour w.r.t. $$\tau_1$$ can be studied from the FSC $${\it{\Gamma}} _{1,1}^ \#$$ (Fig. 6(a)), i.e., $${\it{\Gamma}} _{1,1}$$ when $${F_1} = (0,1.7,0.0504)$$ as defined in the procedure. From Fig. 6(a), we have that as $$\tau_1$$ increases from $$0.3041 - \varepsilon$$ to $$0.3041 + \varepsilon$$ both the CIRs ($$1.9594j$$ and $$3.6118j$$) cross $$\mathbb{C}_0$$ from $$\mathbb{C}_-$$ to $$\mathbb{C}_+$$. Therefore, decreasing $$\tau_1$$ appropriately near $$\overrightarrow {{\tau ^ \#}}$$ may stabilize the system. Next, the asymptotic behaviour w.r.t. $$\tau_2$$ can be studied from the FSC $${\it{\Gamma}} _{2,1}^ \#$$ (Fig. 6(b)), i.e., $${\it{\Gamma}} _{2,1}$$ when $${F_2} = (0.3041,0,0.0504)$$. Similar to the asymptotic behaviour w.r.t. $$\tau_3$$, the variation of $$\tau_2$$ near $$\tau _2^ \#$$ dose not change $$NU((\tau _1^ \#, {\tau _2}, \tau _3^ \#))$$. Thus, there exists a stabilizing $$\overrightarrow \tau$$ sufficiently near $$\overrightarrow {{\tau ^ \#}}$$ with $${\tau _1} > \tau _1^ \#$$. □ Fig. 4. View largeDownload slide FSCs (a) $${{\it{\Gamma}} _{1,1}}$$ and (b) $${{\it{\Gamma}} _{2,1}}$$ for Example 5.3. Fig. 4. View largeDownload slide FSCs (a) $${{\it{\Gamma}} _{1,1}}$$ and (b) $${{\it{\Gamma}} _{2,1}}$$ for Example 5.3. Fig. 5. View largeDownload slide (a) $${{\it{\Gamma}} _{3,1}}$$ and (b) $$NU((\tau _1^ \#, \tau _2^ \#, {\tau _3}))$$ vs. $$\tau_3$$ for Example 5.3. Fig. 5. View largeDownload slide (a) $${{\it{\Gamma}} _{3,1}}$$ and (b) $$NU((\tau _1^ \#, \tau _2^ \#, {\tau _3}))$$ vs. $$\tau_3$$ for Example 5.3. Fig. 6. View largeDownload slide (a) $${\it{\Gamma}} _{1,1}^ \#$$ and (b) $${\it{\Gamma}} _{2,1}^ \#$$ for Example 5.3. Fig. 6. View largeDownload slide (a) $${\it{\Gamma}} _{1,1}^ \#$$ and (b) $${\it{\Gamma}} _{2,1}^ \#$$ for Example 5.3. Finally, at the end of this section, we give some additional discussions and comments. As mentioned and as seen above, the frequency-sweeping test proposed in this article is iterative. In this context, one natural question may arise: How about a direct approach (through the complete stability analysis along an appropriate ray) in the corresponding delay parameter space? For instance, consider the stability for Example 5.2 when $$\tau _1^\# = 0.5,\tau _2^\# = 1.5$$. We may try to study the complete stability w.r.t the ray $$\tau \cdot (1,3)$$, $$\tau \in [0,\infty )$$. Here, $$\tau$$ is the only free parameter. Obviously, the point ($$\tau _1^\# = 0.5,\tau _2^\# = 1.5$$) corresponds to the case $$\tau \cdot (1,3)$$ with $$\tau = 0.5$$. The characteristic function is $${\lambda ^4} + 2{\lambda ^2} + 3z - 3{z^3} + {z^4}$$ ($$z = {e^{ - \tau \lambda }}$$). There are totally four FSCs, as shown in Fig. 7(a). Fig. 7. View largeDownload slide FSCs for complete stability analysis along a ray. (a) Ray $$\tau \cdot (1,3)$$. (b) Ray $$\tau \cdot (5,14)$$. Fig. 7. View largeDownload slide FSCs for complete stability analysis along a ray. (a) Ray $$\tau \cdot (1,3)$$. (b) Ray $$\tau \cdot (5,14)$$. If we consider the point ($$\tau _1^\# = 0.5,\tau _2^\# = 1.4$$) similarly to the previous case study, we shall address the ray $$\tau \cdot (5, 14)$$, $$\tau \in [0,\infty )$$. The corresponding characteristic function is $${\lambda ^4} + 2{\lambda ^2} + 3{z^5} - 3{z^{14}} + {z^{19}}$$ ($$z = {e^{ - \tau \lambda }}$$). As shown in Fig. 7(b), there are totally 19 FSCs! Apparently, the analysis appears to be quite involved with some (significant) increasing complexity. As seen above, if there exists a common factor among all the delay parameters, the complete stability analysis along a ray is possible in theory. However, this is not necessarily simple or practical. For instance, if we try to address Example 5.3, the number of the FSCs is 17000 (the largest common factor among the three delays is $$0.0001$$). The above analysis is consistent with the delay interference phenomenon (Michiels & Niculescu, 2007) and we are able to link it with the number of FSCs. More precisely, a very small perturbation on the delay ratio may result in an increase in the number of FSCs. This will tend to bring more CIRs and hence an increase in $$NU(\tau)$$ for a sufficiently large $$\tau$$, according to the invariance property. 6. Algorithm implementation The effectiveness and advantage of the iterative frequency-sweeping approach are illustrated and discussed in the previous sections. In this section, we will try to implement all the steps by a single program. In this way, the stability can be automatically treated. In this section, we consider the following characteristic function involving $$L$$ independent delay parameters ($${\tau _1}, \ldots ,{\tau _L}$$)   $$\label{CFMTDSwithoutCrossterms}f(\lambda ,{\tau _1}, \ldots ,{\tau _L}) = {p_0}(\lambda ) + \sum\limits_{\ell = 1}^L {{p_\ell }(\lambda ){e^{ - {\tau _\ell }\lambda }}} ,$$ (6.1) where $${p_0}(\lambda ), \ldots, {p_L}(\lambda )$$ are real-coefficient polynomials of $$\lambda$$ such that   $\deg ({p_0}(\lambda )) > \max \{ \deg ({p_1}(\lambda )), \ldots ,\deg ({p_L}(\lambda ))\} .$ The characteristic function (6.1) corresponds to a class of multiple-delay systems of retarded type without cross terms in the characteristic functions. As earlier mentioned, the stability of such multiple-delay systems have been largely studied in the literature (see e.g., Gu et al., 2005 for the case $$L=2$$ and Gu & Naghnaeian, 2011 for the case $$L=3$$. Based on the approach proposed in this article, we now present an algorithm to automatically determine $$NU(\overrightarrow {{\tau ^\# }} )$$ for a given $$\overrightarrow {{\tau ^\# }} = (\tau _1^\# , \ldots ,\tau _L^\# )$$. Algorithm for automatic implementation of iterative frequency-sweeping approach Step 0: Choose a step-length $${{\it{\Delta}} _\omega }$$. Set $$\chi=1$$ and $${F_1} = (0, \ldots ,0)$$. Step 1: Compute $$NU({F_\chi } + \varepsilon \delta (\chi ))$$ by Theorem 3.1. (Step 1 usually may be automatically performed by computer, through solving the polynomial equation $$f(\lambda ,\overrightarrow 0 )$$. Only when the spectrum contains CIRs, we need to additionally invoke the Puiseux series.) Step 2: Sweep $$\omega$$ with the step-length $${{\it{\Delta}} _\omega }$$ and solve $$z_\chi$$ for the equation $$p(\lambda ,{z _\chi},{F_\chi}) = 0$$. Detect the signs of $$\left| {{z_\chi }} \right| - 1$$. If a sign change is found at two adjacent $$\omega$$, say $${\omega '}$$ and $${\omega ''}$$, a CIR is detected as $$\frac{{(\omega ' + \omega '')}}{2}j$$. For any critical pair ($$\lambda, \tau$$) detected in this way, we have that $${\it{\Delta}} N{U_\lambda }(\tau ) = + 1$$ ($$-1$$) if the corresponding sign change is from negative to positive (from positive to negative). Then, we know $$NU({F_\chi } + \tau _\chi ^ \# \delta (\chi ))$$ and detect the CIRs (if any). (The case $${\it{\Delta}} N{U_\lambda }(\tau ) = 0$$ may be later examined from the FSCs. Note that such a case usually does not affect the value of $$NU(\overrightarrow {{\tau ^\# }} )$$). Step 4: If $$\chi < L$$, let $$\chi = \chi + 1$$ and $${F_\chi } = {F_{\chi - 1}} + \tau _{\chi - 1}^ \# \delta (\chi - 1)$$. Return to Step 1. Step 5: We obtain the value of $$NU(\overrightarrow {{\tau ^\# }} )$$ and plot the FSCs $${{\it{\Gamma}} _{\chi ,1}}$$, $$\chi=1, \ldots, L$$. The FSCs may be used to determine if there is the case with $${\it{\Delta}} N{U_\lambda }(\tau ) = 0$$ in Step 2. Also, the FSCs may verify the results of Step 2. The algorithm stops. If there are CIRs at $$\overrightarrow {{\tau ^\# }}$$, we may further generate the FSCs $${\it{\Gamma}} _{\ell ,1}^\#$$ as in Step 4 of the iterative frequency-sweeping approach, to analyse the asymptotic behaviour. Remark 6.1 By using the above algorithm, the calculation error for CIRs is kept within $$\pm \frac{{{\it{\Delta}} _\omega }}{2}j$$. Example 6.1 Consider the three-delay system in Example 5.3. We now compute $$NU(\overrightarrow \tau )$$ at some points, by using a MATLAB program (based on the algorithm proposed in this section). The step-length is chosen as $${{\it{\Delta}} _\omega }=0.01$$. In Table 1, we list the results and the computation time (on a Laptop with an Intel Core 2.50 GHz CPU with 8 G RAM). For these points, the analysis is automatically performed by computer. Table 1 $$NU((\tau _1^\# ,\tau _2^\# ,\tau _3^\# ))$$ calculation for some representative points $$(\tau _1^\# = 0.01$$ and $$\tau _2^\# = 1.7)$$   $$NU((\tau _1^\# ,\tau _2^\# ,\tau _3^\# ))$$  Computation time (s)  $$\tau _3^\# = 0.1$$  2  0.009605  $$\tau _3^\# = 0.3$$  0  0.009712  $$\tau _3^\# = 0.6$$  2  0.009579  $$\tau _3^\# = 0.9$$  0  0.009573  $$\tau _3^\# = 1.5$$  2  0.009616  $$\tau _3^\# = 2.5$$  4  0.009555  $$\tau _3^\# = 3.0$$  2  0.009604  $$\tau _3^\# = 3.7$$  0  0.009571  $$\tau _3^\# = 3.77$$  2  0.009547  $$\tau _3^\# = 3.8$$  4  0.009589    $$NU((\tau _1^\# ,\tau _2^\# ,\tau _3^\# ))$$  Computation time (s)  $$\tau _3^\# = 0.1$$  2  0.009605  $$\tau _3^\# = 0.3$$  0  0.009712  $$\tau _3^\# = 0.6$$  2  0.009579  $$\tau _3^\# = 0.9$$  0  0.009573  $$\tau _3^\# = 1.5$$  2  0.009616  $$\tau _3^\# = 2.5$$  4  0.009555  $$\tau _3^\# = 3.0$$  2  0.009604  $$\tau _3^\# = 3.7$$  0  0.009571  $$\tau _3^\# = 3.77$$  2  0.009547  $$\tau _3^\# = 3.8$$  4  0.009589  Table 1 $$NU((\tau _1^\# ,\tau _2^\# ,\tau _3^\# ))$$ calculation for some representative points $$(\tau _1^\# = 0.01$$ and $$\tau _2^\# = 1.7)$$   $$NU((\tau _1^\# ,\tau _2^\# ,\tau _3^\# ))$$  Computation time (s)  $$\tau _3^\# = 0.1$$  2  0.009605  $$\tau _3^\# = 0.3$$  0  0.009712  $$\tau _3^\# = 0.6$$  2  0.009579  $$\tau _3^\# = 0.9$$  0  0.009573  $$\tau _3^\# = 1.5$$  2  0.009616  $$\tau _3^\# = 2.5$$  4  0.009555  $$\tau _3^\# = 3.0$$  2  0.009604  $$\tau _3^\# = 3.7$$  0  0.009571  $$\tau _3^\# = 3.77$$  2  0.009547  $$\tau _3^\# = 3.8$$  4  0.009589    $$NU((\tau _1^\# ,\tau _2^\# ,\tau _3^\# ))$$  Computation time (s)  $$\tau _3^\# = 0.1$$  2  0.009605  $$\tau _3^\# = 0.3$$  0  0.009712  $$\tau _3^\# = 0.6$$  2  0.009579  $$\tau _3^\# = 0.9$$  0  0.009573  $$\tau _3^\# = 1.5$$  2  0.009616  $$\tau _3^\# = 2.5$$  4  0.009555  $$\tau _3^\# = 3.0$$  2  0.009604  $$\tau _3^\# = 3.7$$  0  0.009571  $$\tau _3^\# = 3.77$$  2  0.009547  $$\tau _3^\# = 3.8$$  4  0.009589  The results listed in Table 1 may be verified by Fig. 14 in Gu & Naghnaeian (2011) (see also Fig. 8(a) in this article). With some slight modifications of the program mentioned above, we may obtain the SCS for the three-delay system. Here, we generate the SCS for the case $$\tau_2 = 1.7$$, see Fig. 8(a). Figure 8(a) obtained here is same as Fig. 14 of Gu & Naghnaeian (2011), which is generated by a different approach using some geometric arguments. The computation time to generate the data for Fig. 8(a) by our algorithm is 2.566361 s (on the same Laptop). It is worth to mention that the $$NU(\overrightarrow \tau )$$ distribution is directly obtained by our program (without extra calculation), since the asymptotic behaviour analysis for CIRs is covered by the approach in this article. Moreover, we can further determine the 3-D SCS in the $$(\tau_1, \tau_2, \tau_3)$$-space. For a clear illustration, we here give the SCS for a domain $$({\tau _1},{\tau _2},{\tau _3}) \in [1.5, 2.5] \times [0, 3] \times [0, 3]$$. The 3-D SCS is shown in Fig. 8(b). The computation time to obtain the data for Fig. 8(b) is 35.552933 s (on the same Laptop).□ Fig. 8. View largeDownload slide Stability crossing set for Example 6.1. (a) Cross section with $$\tau_2 = 1.7$$. (b) In $$(\tau_1, \tau_2, \tau_3)$$-space. Fig. 8. View largeDownload slide Stability crossing set for Example 6.1. (a) Cross section with $$\tau_2 = 1.7$$. (b) In $$(\tau_1, \tau_2, \tau_3)$$-space. Remark 6.2 As seen in Gu & Naghnaeian (2011), it is not trivial to determine the SCS for the multiple-delay system (6.1). In this section, we propose a different approach for this task. As illustrated, the asymptotic behaviour analysis for CIRs is included by our algorithm and hence the distribution of $$NU(\tau)$$ can be directly examined. In our opinion, this is one of the advantages of the approach proposed in our article. We think that this observation will open some new perspectives for further research. For instance, we may study how to obtain the SCSs for more general multiple-delay systems (e.g., when cross terms $${e^{ - ({\tau _1} + {\tau _2})\lambda }}$$, $${e^{ - ({\tau _1} + {\tau _3})\lambda }}$$, $${e^{ - ({\tau _2} + {\tau _3})\lambda }}$$, $${e^{ - ({\tau _1} + {\tau _2} + {\tau _3})\lambda }}$$, $$\ldots$$, are added in (6.1)) in the future. But this scope is out of our purposes in this article. 7. Conclusion We analyse the stability for linear systems with multiple (incommensurate) delay parameters. As the asymptotic behaviour analysis of the CIRs w.r.t multiple delay parameters corresponds to an open problem, we propose an indirect yet effective methodology called the iterative frequency-sweeping approach. We first study the complete stability in the case where only one delay parameter is free while the others are fixed. The invariance property (regarding the asymptotic behaviour of the CIRs) in this case is confirmed by extending the frequency-sweeping framework recently proposed for studying linear systems with single delay parameter. As a result, the complete stability problem can be easily studied by employing a frequency-sweeping test. Based on the above results, we next present an iterative frequency-sweeping approach to analyse the stability for any given combination of multiple delays. Using this approach, we may accurately compute the number of unstable characteristic roots. Furthermore, if the system has CIRs, we may analyse the asymptotic behaviour of the CIRs w.r.t each delay parameter. Consequently, we may determine if there exists a stabilizing combination of multiple delays sufficiently close to the given one, at which the system is asymptotically stable. Finally, we develop an algorithm with which the stability for a class of multiple-delay systems may be easily implemented. This work opens some new perspectives for further research. Funding National Natural Science Foundation of China (61473065); and Fundamental Research Funds for the Central Universities of China (N160402001), in part. Footnotes 1 Roughly speaking, the complete stability problem refers to finding exhaustively the stability interval(s) for the delay parameter $$\tau$$ along the whole $$\tau$$-axis $$\tau \in [0,\infty )$$. 2 A system has CIRs iff the delay parameters lie in the SCS. 3 Some results for specific two-free-parameter problems can be found in, e.g., Section 6.2.1 of Arnold et al. (2012), Beringer & Richard-Jung (2003) and Soto & Vicente (2011). For more general problems, the analysis will become much more complicated and no general result has been reported so far. 4 We do not need to specifically generate the FSCs $${\it{\Gamma}} _{L, i}^ \# (\omega )$$ since they are exactly the FSCs obtained by Step 3 when $$\chi = L$$. References Arnold V. I., Gusein-Zade S. M. & Varchenko A. N. ( 2012) Singularities of Differentiable Maps.  vol. 2: Monodromy and Asymptotics of Integrals. New York: Springer. Beringer F. & Richard-Jung F. ( 2003) Multi-variate polynomials and Newton-Puiseux expansions. Symbolic and Numerical Scientific Computation . ( Winkler F. & Langer U. eds). Berlin: Springer, pp. 240– 254. Google Scholar CrossRef Search ADS   Boussaada I. & Niculescu S. -I. ( 2016a) Characterizing the codimension of zero singularities for time-delay systems: a link with Vandermonde and Birkhoff incidence matrices. Acta Appl. Math. , 145, 47– 88. Google Scholar CrossRef Search ADS   Boussaada I. & Niculescu S. -I. ( 2016b) Tracking the algebraic multiplicity of crossing imaginary roots for generic quasipolynomials: A Vandermonde-based approach. IEEE Trans. Autom. Control , 61, 1601– 1606. Google Scholar CrossRef Search ADS   Cooke K. L. & van den Driessche P. ( 1986) On zeros of some transcendental equations. Funkcialaj Ekvacioj , 29, 77– 90. Delice I. I. & Sipahi R. ( 2012) Delay-independent stability test for systems with multiple time-delays. IEEE Trans. Autom. Control,  57, 963– 972. Google Scholar CrossRef Search ADS   Gu K. & Naghnaeian M. ( 2011) Stability crossing set for systems with three delays. IEEE Trans. Autom. Control,  56, 11– 26. Google Scholar CrossRef Search ADS   Gu K., Niculescu S. -I. & Chen J. ( 2005) On stability crossing curves for general systems with two delays. J. Math. Anal. Appl. , 311, 231– 253. Google Scholar CrossRef Search ADS   Hale J. K. & Verduyn Lunel S. M. ( 1993) Introduction to Functional Differential Equations.  New York: Springer. Google Scholar CrossRef Search ADS   Hassard B. D. ( 1997) Counting roots of the characteristic equation for linear delay-differential systems. J. Differ. Equ. , 136, 222– 235. Google Scholar CrossRef Search ADS   Hu G. -D. & Liu M. ( 2007) Stability criteria of linear neutral systems with multiple delays. IEEE Trans. Autom. Control , 52, 720– 724. Google Scholar CrossRef Search ADS   Lee M. S. & Hsu C. S. ( 1969) On the $$\tau$$-decomposition method of stability analysis for retarded dynamical systems. SIAM J. Control , 7, 242– 259. Google Scholar CrossRef Search ADS   Li X. -G., Niculescu S. -I. & Çela A. ( 2015) Analytic Curve Frequency-Sweeping Stability Tests for Systems with Commensurate Delays.  London: Springer. Google Scholar CrossRef Search ADS   Li X. -G., Niculescu S. -I., Çela A, Zhang L. & Li X. ( 2017) A frequency-sweeping framework for stability analysis of time-delay systems. IEEE Trans. Autom. Control , 62, 3701– 3716. Google Scholar CrossRef Search ADS   Michiels W. & Niculescu S. -I. ( 2007) Characterization of delay-independent stability and delay interference phenomena. SIAM J. Control Optim. , 45, 2138– 2155. Google Scholar CrossRef Search ADS   Michiels W. & Niculescu S. -I. ( 2014) Stability, Control, and Computation for Time-Delay Systems: An Eigenvalue-Based Approach.  Philadelphia: SIAM. Google Scholar CrossRef Search ADS   Naghnaeian M. & Gu K. ( 2013) Stability crossing set for systems with two scalar-delay channels. Automatica , 49, 2098– 2106. Google Scholar CrossRef Search ADS   Neimark J. ( 1949) D-subdivisions and spaces of quasi-polynomials. Prikl. Mat. Meh. , 13, 349– 380. Niculescu S. -I. ( 2001) Delay Effects on Stability: A Robust Control Approach.  Heidelberg: Springer. Niculescu S. -I., Kim P., Gu K., Lee P. & Levy D. ( 2010) Stability crossing boundaries of delay systems modeling immune dynamics in leukemia. Discrete Contin. Dyn. Syst.-Ser. B (DCDS-B) , 15, 129– 156. Sipahi R. & Delice I. I. ( 2009) Extraction of 3D stability switching hypersurfaces fo a time delay system with multiple fixed delays. Automatica , 45, 1449– 1454. Google Scholar CrossRef Search ADS   Sipahi R. & Delice I. I. ( 2011) Advanced clustering with frequency sweeping methodology for the stability analysis of multiple time-delay systems. IEEE Trans. Autom. Control , 56, 467– 472. Google Scholar CrossRef Search ADS   Sipahi R., Niculescu S. -I., Abdallah C. T., Michiels W. & Gu K. ( 2011) Stability and stabilization of systems with time delay. IEEE Control Syst. Mag. , 31, 38– 65. Google Scholar CrossRef Search ADS   Sipahi R. & Olgac N. ( 2006) A unique methodology for the stability robustness of multiple time delay systems. Syst. Control Lett. , 55, 819– 825. Google Scholar CrossRef Search ADS   Soto M. J. & Vicente J. L. ( 2011) The Newton procedure for several variables. Linear Algebra Appl. , 435, 255– 269. Google Scholar CrossRef Search ADS   Stépán G. ( 1979) On the stability of linear differential equations with delay. Coll. Math. Soc. J. Bolyai , 30, 971– 984. Stépán G. ( 1989) Retarded Dynamical Systems: Stability and Characteristic Functions.  Harlow: Longman Scientific & Technical. Vyhlídal T. & Zítek P. ( 2009) Modification of Mikhaylov criterion for neutral time-delay systems. IEEE Trans. Autom. Control , 54, 2430– 2435. Google Scholar CrossRef Search ADS   Appendix: Proof of Theorem 3.2 Consider the following characteristic function   \begin{eqnarray} \label{GeneralQuasipolynomial}f(\lambda ,\tau ) = {a_0}(\lambda ) + {a_1}(\lambda ){e^{ - \tau \lambda }} + \cdots + {a_q}(\lambda ){e^{ - q\tau \lambda }}, \end{eqnarray} (A.1) where the coefficient functions $${a_0}(\lambda ), \ldots ,{a_q}(\lambda )$$are only required to be analytic in $$\mathbb{C}_0 \backslash \{ 0\}$$ (have in mind that usually we preclude the trivial case where $$\lambda = 0$$ is a characteristic root). The above characteristic function (A.1) is called a general quasipolynomial, corresponding to a broad class of time-delay systems. Note that the characteristic function (A.1) reduces to the widely-studied quasipolynomials, corresponding to the retarded-type and the neutral-type time-delay systems, if the coefficient functions $${a_0}(\lambda ), \ldots ,{a_q}(\lambda )$$ are restricted to be polynomials of $$\lambda$$. Recently, the invariance property of the CIRs for the general quasipolynomial (A.1) was confirmed in Li et al. (2017). Clearly, the characteristic function $$f(\lambda ,{\tau _\chi },{F_\chi })$$ (3.1) is in the form of (A.1) as   ${\widetilde a_{\chi 0}}(\lambda ) + {\widetilde a_{\chi 1}}(\lambda ){e^{ - {\tau _\chi }\lambda }} + \cdots + {\widetilde a_{\chi {q_\chi }}}(\lambda ){e^{ - {q_\chi }{\tau _\chi }\lambda }},$ where the coefficient functions $${\widetilde a_{\chi 0}}(\lambda ), \ldots ,{\widetilde a_{\chi {q_\chi }}}(\lambda )$$ are polynomials in $$\lambda$$ and $${e^{ - \tau _\ell ^\# \lambda }}$$ ($$\ell \in \{ 1, \ldots ,L\}$$, $$\ell \ne \chi$$). We may now prove Theorem 3.2 as the characteristic function $$f(\lambda ,{\tau _\chi },{F_\chi })$$ falls in the class of general quasipolynomial (A.1), since the coefficient functions $${\widetilde a_{\chi 0}}(\lambda ), \ldots ,{\widetilde a_{\chi {q_\chi }}}(\lambda )$$ are analytic in $$\mathbb{C}$$. © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved.

### Journal

IMA Journal of Mathematical Control and InformationOxford University Press

Published: Nov 22, 2017

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations