# Real curve analysis and stability of time-delay systems

Yacine Bouzidi, Adrien Poteaux

► To cite this version:

Yacine Bouzidi, Adrien Poteaux. Real curve analysis and stability of time-delay systems. 7th IFAC Symposium on System Structure and Control SSSC 2019, Daniela Danciu, Sep 2019, Sinaia (Romania), Romania. pp.94-98, 10.1016/j.ifacol.2019.11.033. hal-04262761

HAL Id: hal-04262761

https://inria.hal.science/hal-04262761v1

Submitted on 27 Oct 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Real curve analysis and stability of time-delay systems

Yacine Bouzidi* Adrien Poteaux**

* Inria Lille - Nord Europe, Parc Scientifique de la Haute Borne, 40 Avenue Halley, Bat A - Park Plaza, 59650 Villeneuve d'Ascq, France (e-mail: yacine.bouzidi@inria.fr).

** Univ. Lille, CNRS, INRIA, UMR 9189 - CRISTAL - Centre de Recherche en Informatique Signal et Automatique de Lille, F-59000 Lille, France (e-mail: adrien.poteaux@univ-lille.fr).

**Abstract:** In this work, we investigate the asymptotic stability of LTI differential commensurate time-delay systems whose dynamics are defined in the frequency domain by quasi-polynomials of the form $f(s, \tau) = \sum_{j=0}^{m} p_j(s) e^{-j\tau s}$. We propose a new approach for studying the stability with respect to the delay: we determine the asymptotic behavior of the roots of the quasi-polynomial near to the critical pairs of $f(s, \tau)$ by analysing the intersection of an associated algebraic space curve with sufficiently refined 3D real boxes. Compared to the existing methods, our method does not require any Puiseux series computation.

**Keywords:** Time delay systems, curve analysis, stability, certified tests, symbolic-numeric methods.

## 1. INTRODUCTION

In this paper, we investigate the asymptotic stability of LTI differential commensurate time-delay systems whose dynamics are defined by the following state-space representation

$$ \dot{x}(t) = \sum_{k=0}^{m} A_k x(t - k \tau) \quad (1) $$

where $\tau \in R_+ := \{\tau \in R \mid \tau \ge 0\}$, $m \in \mathbb{N}$ and $A_0, \dots, A_m$ are constant matrices with entries in a field $\mathbb{K}$ (e.g. $\mathbb{Q}$ or $\mathbb{R}$). The characteristic function of (1) is defined by a quasi-polynomial of the form

$$ f(s, \tau) = \sum_{j=0}^{m} p_j(s) e^{-j\tau s}, \quad (2) $$

where the $p_j$'s are polynomials in the complex variable $s$ with coefficients in $\mathbb{K}$.

Such a system is said to be asymptotically stable if all the zeros of the quasi-polynomial $f(s, \tau)$ have negative real parts, i.e., $f(s, \tau) \ne 0$ for all $s$ in the closed right half-plane $\bar{\mathbb{C}}_+ := \{s \in \mathbb{C} \mid \Re(s) \ge 0\}$. Considering the delay $\tau$ as a parameter, the stability property will obviously depend on the value of $\tau$ and a classical question is then to study the asymptotic stability with respect to the delay value $\tau$.

A classical approach for studying the asymptotic stability depending on the delay is based on the analysis of the so-called critical pairs of $f(s, \tau)$, that is the pairs $(\omega \in \mathbb{R}, \tau \in \mathbb{R}_+)$ such that $f(i\omega, \tau) = 0$; see for instance Gu et al. [2003], Li et al. [2015], Marshall et al. [1992], Niculescu [2001], Olgac and Sipahi [2002] and the references therein. If such critical pairs exist, the stability is derived from the way these critical imaginary roots behave under small variations of the time-delay $\tau$. As a result, this approach

yields a decomposition of the $\tau$-axis into a set of disjoint intervals in which the considered system is either stable or unstable.

Algorithmically, the above approach can be decomposed in two independent problems. First, the computation of the critical pairs of $f(s, \tau)$, and second, for each critical pair $(s_0, \tau_0)$, the study of the asymptotic behavior of the root $s_0$ with respect to a small variation of the delay $\tau_0$.

There exist numerous methods for computing the critical pairs of $f(s, \tau)$ (see for instance Bouzidi et al. [2016], Niculescu [2001], Olgac and Sipahi [2002]). In Bouzidi et al. [2016] for example, based on an algebraic formulation for the critical pairs of $f(s, \tau)$, the problem is reduced to solving a system of algebraic equations in two variables and classical symbolic methods are used to obtain certified numerical approximations for the critical pairs. Such approximations consist in a pair of intervals isolating the two coordinates $s$ and $\tau$ of the critical pairs.

For the study of the asymptotic behavior of critical imaginary roots, several methods have been developed for the case of simple imaginary roots (see, e.g., Gu et al. [2003], Li et al. [2015], Marshall et al. [1992], Niculescu [2001], Olgac and Sipahi [2002] and the references therein). For multiple imaginary roots (MIR), recent methods based on Puiseux series was developed in Bouzidi et al. [2016], Li et al. [2015, 2013]. These methods require in particular the symbolic computation of the Puiseux series of $f(s, \tau)$ at all the critical pairs.

**Contributions and outline of the paper.** In this paper, we assume that the critical pairs of $f(s, \tau)$ are given as certified numerical approximations (using for instance the method in Bouzidi et al. [2016]) and we focus on the analysis of the asymptotic behavior of $s$
at these critical pairs. We propose a new approach that avoids any Puiseux series computation at the critical pairs. Instead, the asymptotic behavior is derived through the analysis around the critical pairs of the real branches of a space algebraic curve obtained by local approximation of the quasi-polynomial $f(s, t)$. This analysis is done while avoiding as much as possible symbolic computations, using interval arithmetic instead. More precisely, given a critical pair $(\omega_0, \tau_0)$, our approach proceeds in the following way:

(1) We compute an upper bound $m$ for the multiplicity of $(\omega_0, \tau_0)$ by evaluating the successive derivatives of $f$ up to the point that we can certify a non-zero evaluation. See Section 3.1.

(2) We consider $f(s+\omega_0, \tau+\tau_0)$ as an element of $\mathbb{C}[[s, \tau]]$ by replacing the exponential by their local Taylor series. Note that its coefficients are given with error bounds. To make computations possible, we truncate these series modulo some powers of $s$ and $\tau$, taking at least $m+1$ as a truncation bound. This provides a polynomial $p(s, \tau)$ that is an approximation of the studied quasi-polynomial (note that approximations are coming from both numerical errors and the truncation process). See Section 3.2.

(3) Finally, we study the local behaviour of the real curve defined by $\{p_r(s_r, s_i, \tau) = p_i(s_r, s_i, \tau) = 0\}$ where $p_r$ and $p_i$ are respectively the real and the imaginary part of $p(s_r + i s_i, \tau)$. More precisely, we isolate the curve in a 3D box that contains only the germ at $(0, 0, 0)$ of the real curve; then, one can conclude about the stability of the system by looking at the intersections of real branches with this box.

**Organisation of the paper.** For a sake of clarity, we briefly recall in Section 2 the method presented in Bouzidi et al. [2016] for the computation of the critical pairs of $f(s, t)$. Then, in Section 3, we show how a quasi-polynomial can be locally approximated at a given point by a polynomial while preserving the local topology of the latter. Using this approximated curve, we present in Section 4 a strategy to compute the asymptotic behavior of this curve at a given point by looking at the intersection of the latter with some specific 3D boxes. This enables us to study the stability of the system.

## 2. CRITICAL PAIRS OF A QUASI-POLYNOMIAL

As mentioned in the introduction, to study the asymptotic stability of a system of the form (1), we first need to compute the critical pairs of the corresponding quasi-polynomial $f(s, \tau)$. These critical pairs are the intersection points of the curve defined by the quasi-polynomial, $f(s, \tau) = 0$, with the imaginary axis in $s$, that is

$$ \{(\omega \in \mathbb{R}, \tau \in \mathbb{R}_+) | f(i\omega, \tau) = 0\} $$

We now briefly recall the method given in Bouzidi et al. [2016] for the certified computation of these critical pairs.

The main idea is to use the so-called Rekasius transformation, that allows to recover the infinite number of zeros of $f(i\omega, \tau)$ from the real zeros of a zero-dimensional algebraic bivariate system. More precisely, we replace the terms $e^{-\tau i\omega}$ in the quasipolynomial $f(i\omega, \tau)$ by the rational fraction $\frac{1-\gamma i\omega}{1+\gamma i\omega}$ with $\gamma \in \mathbb{R}$. Cleaning the denominators, one obtain a complex polynomial of the form $\mathcal{R}(\omega, \gamma) +$

$i\mathcal{I}(\omega, \gamma)$. Computing the critical pairs then amounts to identifying the real zeros of the following bivariate polynomial system:

$$ \begin{cases} \mathcal{R}(\omega, \gamma) = 0, \\ \mathcal{I}(\omega, \gamma) = 0, \end{cases} \qquad (3) $$

and from each solution $(\omega, \gamma)$ of (3) the critical delays can then be obtained by:

$$ \tau_k = \frac{2}{\omega} (\arctan(\omega \gamma) + k \pi), \quad k \in \mathbb{Z}. \qquad (4) $$

The method proposed in Bouzidi et al. [2016] for computing the real solutions of the above algebraic system consists first in computing a univariate representation of the solutions (5), which is roughly speaking, a one-to-one mapping between the solutions of the system (3) and the root of a univariate polynomial $f(t)$ (see Rouillier [1999] for details), then, numerical approximations for the real solutions are obtained from this representation.

$$ \begin{cases} f(t) = 0 \\ \omega = g_{\omega}/g_{0} \\ \gamma = g_{\gamma}/g_{0} \end{cases} \qquad (5) $$

From the above representation, one can compute intervals with rational endpoints that isolate the coordinates $\omega$ and $\gamma$ of the solutions. Moreover, by substituting the obtained intervals in (4), one also obtain isolating intervals for the corresponding delays $\tau$. Another advantage of the previous representation is that it allows to control the precision of the numerical approximation of the solutions. Indeed, given an isolating interval for a root of the polynomial $f(t)$, one can refine this interval up to a desired precision while it still contains the root which results to refined intervals for $\omega, \gamma$ and $\tau$. As we will see in the sequel, this feature is very useful in the description and the design of our algorithm.

## 3. LOCAL APPROXIMATION

### 3.1 An upper bound on the multiplicities of critical pairs

An important information that we need for the description of our algorithm is the multiplicities of the critical pairs as roots of the quasi-polynomial $f(s, \tau)$. Computing explicitly such multiplicities is however a challenging problem which is yet not solved (for detail, see Boussaada and Niculescu [2016] and references therein). Therefore, instead of exact multiplicities, we only compute upper bounds on the latter, which will be sufficient for our purpose. Let first recall the definition of the multiplicity of a critical pair of a quasi-polynomial $f(s, \tau)$.

**Definition 1.** A root $(\omega_0, \tau_0)$ of the quasi-polynomial $f(i\omega, \tau)$ is said to be of multiplicity $k$ if and only if $f(\omega_0, \tau_0) = \frac{\partial f}{\partial s}(\omega_0, \tau_0) = \cdots = \frac{\partial^{k-1} f}{\partial s^{k-1}}(\omega_0, \tau_0) = 0$ and $\frac{\partial^k f}{\partial s^k}(\omega_0, \tau_0) \neq 0$

According to the previous definition, if the evaluation at the point $(\omega_0, \tau_0)$ can be achieved in an exact way, a straightforward approach for computing the multiplicity of a critical pair $(\omega_0, \tau_0)$ is to evaluate it on the successive derivatives of $f(s, \tau)$ until a non-zero value is obtained.

The key point of this paper is to avoid as much as possible to run costly symbolic computations. In particular, we are
not going to compute the multiplicities of the critical pairs, but only a bound for each of them. To do so, we will just evaluate the successive derivatives $\frac{\partial f^k}{\partial k^s}$ on the given 2D box (using interval arithmetic - see Moore [1979] for details), until we can conclude that the m-th derivative is non zero in this box. In practice, this strategy will often provide the correct multiplicity, but over-bound might happen, so that we cannot assume to know the correct multiplicity.

In particular, if we get a bound $m = 1$ for a critical pair, then this multiplicity is correct and we can conclude that the system is not stable.

## 3.2 Local approximation

Given a critical pair $(i\omega_0, \tau_0)$, the quasi-polynomial $f(s, \tau)$ can be considered in $\mathbb{K}[(s - i\omega_0, \tau - \tau_0)]$. Thus, it can be approximated by a polynomial via a truncation of the powers of $(s - i\omega_0)$ and $(\tau - \tau_0)$. In practice, we only have to use the Taylor series of the exponential to convert our quasi-polynomial in series. More precisely, it is well known that $e^x$ can be approximated by

$$ \sum_{k=0}^{n-1} \frac{x^k}{k!} $$

with an error bounded by $\frac{x^n}{n!} e^x$. We thus can approximate $f(s + i\omega_0, \tau + \tau_0)$ as

$$ \sum_{j=0}^{m} p_j(s + i\omega_0)e^{-j i \omega_0 \tau_0} e^{-j i \omega_0 \tau} e^{-j \tau_0 s} e^{-j s \tau} $$

by replacing each last three exponentials by their truncated Taylor series (with $x$ being respectively $-j i \omega_0 \tau$, $-j \tau_0 s$ and $-j s \tau$). This provides us a polynomial $p(s, \tau)$ for which there are two different approximations:

* the truncation we just made (for which we have a good majoration of the error), that can be improved by increasing the truncation order,

* the approximations we have for $\tau_0$ and $\omega_0$ ($p$ having coefficients defined by interval), that we can improve on demand from the previous computations.

Note that we will always take a truncation bound $n$ greater than the bound $m$ we have for the multiplicity of the critical pair.

# 4. ASYMPTOTIC BEHAVIOR

As shown in Section 2, one can compute a set of disjoint 2D boxes in $\mathbb{R}^2$ that isolate all the critical pairs of $f(s, \tau)$. Moreover, thanks to the symbolic univariate representation of these pairs (5), the corresponding boxes can be refined up to an arbitrary precision. Also, for a given critical pair $(\omega_0, \tau_0)$ of $f(s, \tau)$, Section 3 provides a local approximation of the quasi-polynomial $f(s, \tau)$ by a two variables polynomial $p(s, \tau)$ that preserves its local topology.

Given such a polynomial $p(s, \tau)$, as well as a critical pair of $p$ (a zero of the form $(i\omega_0, \tau_0)$), the goal of this section is to show how to analyze the asymptotic behavior of this critical pair. As a consequence, provided that this polynomial results from an approximation of a quasi-polynomial, this will naturally provide the asymptotic behavior of the critical pair as a zero of the quasi-polynomial.

## 4.1 Asymptotic behavior through curve intersection in $\mathbb{R}^3$

Let us consider a polynomial $p(s, \tau)$ where $s \in \mathbb{C}$ and $\tau \in \mathbb{R}$. In the following, we set $s = x + iy$, where $x, y \in \mathbb{R}$. For simplicity, we also denote to the variable $\tau$ by $z$. Substituting $s = x + iy$ as well as $\tau = z$ in the polynomial $p$ yields $p(x + i y, z) = g(x, y, z) + i h(x, y, z)$ where $g, h$ are polynomials in three variables with real coefficients. With this notation, we define a critical pair of $p$ as a point $(0, y_0, z_0) \in \mathbb{R}^3$ such that $p(i y_0, z_0) = 0$, or equivalently, as the intersection point of the axis $x = 0$ with curve $C$ defined by

$$ C := \begin{cases} g(x, y, z) = 0 \\ h(x, y, z) = 0 \end{cases} $$

and whose coordinate in $y$ (resp. $z$) is $y_0$ (resp. $z_0$). In the latter case, and to avoid confusion with critical points of the curve, we refer to it as a zero-crossing point. Moreover, for simplicity, we assume in the sequel that the curve $C$ does not admit a whole component lying in a plan orthogonal to the $x$-axis. This translates into the fact that for any fixed value $x_0$ of $x$, the algebraic system $\{g(x_0, y, z) = h(x_0, y, z) = 0\}$ admits a finite number of complex zeros.

Given a zero-crossing point $(0, y_0, z_0)$, the objective in the following is to analyse the asymptotic behavior of the coordinate $x = 0$ under a small variation of the coordinate $z_0$. We suppose that this zero-crossing point $(0, y_0, z_0)$ is given by a 3D box in $\mathbb{R}^3$ that isolates it from other zero-crossing points of the curve. Such a box can be obtained for example by extending the 2D box obtained from the representation (3) on the $x$-axis.

First, we need to introduce the notion of $x$-critical point of $C$. A point $(x, y, z)$ of $C$ is called $x$-critical if it is a critical point of the following projection map

$$ \Pi_x : \begin{array}{c} C \\ (x, y, z) \end{array} \to \mathbb{R} \\ \mapsto x $$

Algebraically, a point of the curve $C$ is an $x$-critical point if and only if it is a zero of the polynomial $\frac{\partial g}{\partial y} \frac{\partial h}{\partial z} - \frac{\partial g}{\partial z} \frac{\partial h}{\partial y}$. In the following, for a sake of briefness, we denote to this polynomial simply by $d(x, y, z)$.

As a consequence, we can identify two different kinds of zero-crossing points of $C$.

*Definition 2.* Let $(0, y_0, z_0)$ be a zero-crossing point of $C$. This point is called $x$-critical if and only if $d(0, y_0, z_0) = 0$. Otherwise, it is called regular.

Given a zero-crossing point $(0, y_0, z_0)$ of $C$, if the latter is a regular zero-crossing points, then the analysis of the asymptotic behavior near to it is straightforward. Indeed, at the neighborhood of these kind of points, the curve $C$ consists in one regular branch that passes through the point $(0, y_0, z_0)$. As a consequence, the variation of $x$ can be deduced from the sign of $d(0, y_0, z_0)$. As $(0, y_0, z_0)$ is not an $x$-critical point, one can compute a box that contains the latter (in fact, it is simply a 2D box in $y$ and $z$), and refine it until the interval resulting from the evaluation $d(B)$ does not contain zero.

For $x$-critical zero-crossing points, analyzing the asymptotic behavior of $C$ is more involved. In broad words, our
approach consists in computing a 3-D box B that isolate this point from other x-critical points of the curve C. Then, this box B is refined until all the intersections of C with the boundary of B are located on the sides that are orthogonal to the x-axis. This guarantees that all the branches that intersect the boundary of the refined box pass through the considered zero-crossing point. The behavior is then easily deduced from these intersections.

**Definition 3.** A box B containing a point p := (0, y₀, z₀) of C is called isolating for p if and only if all the real branches of C in B pass through the point p.

**Theorem 1.** Let B be a box that contains one x-critical zero-crossing point (0, y₀, z₀) of C and such that the two following conditions hold

* B does not contain any other x-critical point of C

* The only intersection points of C with B are located on the sides of B that are orthogonal to the x-axis.

Then, B is an isolating box for (0, y₀, z₀).

**Proof.** This is a direct consequence of the fact that B contains a single zero-crossing point of C and no other x-critical points. Indeed, since all the branches of C intersect B on the sides that are orthogonal to the x-axis, if one of these branches does not pass through (0, y₀, z₀), this implies that either it intersects the axis x = 0 at y₁ ≠ y₀, or it passes through an x-critical point, which contradicts our assumption.

**Remark 2.** Figures 1 and 2 below show the projection of two non-isolating 3D boxes on the (x - z)-plan. The first situation corresponds to a box B that contains an x-critical zero-crossing point (lying on the branch depicted in blue) as well as another zero-crossing point (lying on the branch depicted in red). The second situation corresponds to a box B that contains two different x-critical points of C among which one is a zero-crossing point.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Two zero-crossing points</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2. A zero-crossing point and an x-critical point</figcaption></figure>

According to the previous theorem, given a box that contains an x-critical zero-crossing point (0, y₀, z₀) and no other x-critical points of C, it is sufficient to refine the latter until all the intersections with C lie on the sides that are orthogonal to the x-axis, in order to get an isolating box for (0, y₀, z₀). This is always possible since we have assumed that there is no branch of C lying in a plan orthogonal to the x-axis.

## 4.2 Computational aspects

Let (0, y₀, z₀) be an x-critical zero-crossing point of C and suppose given a 3D box B that isolate the point (0, y₀, z₀) from other zero-crossing points of C.

A first computational step consists in refining the box B until (0, y₀, z₀) is the only x-critical point of C that lies in the refined box. To do so, a straightforward approach consists in computing disjoint boxes that isolate the real zeros of the following system

$$S_x := \begin{cases} g(x, y, z) = 0 \\ h(x, y, z) = 0 \\ d(x, y, z) = 0 \end{cases}$$

As for the system in (3), this can be done by computing a univariate representation of the solutions followed by a numerical approximation. This provides a set of disjoint 3D boxes that isolate the real zeros of $S_x$. Thenceforth, it suffices to refine the box B until it overlaps with only one box of this set to get a refined box that contains only the x-critical point (0, y₀, z₀). Let denote to this new box by $B'$.

The second computational step consists in refining the box $B'$ until all the intersections with C lie only on the desired sides. Let $B^*$ be a box resulting from the refinement of $B'$ up to a certain precision. We denote by $y_{min}, y_{max}$ (resp. $z_{min}, z_{max}$), the rational endpoints of its interval on y (resp. z). Then, the box $B^*$ admits intersections only on the sides that are orthogonal to the x-axis if and only if the following algebraic systems do not have real zeros.

$$\begin{cases} g(x, y_{min}, z) = 0 \\ h(x, y_{min}, z) = 0 \end{cases} \quad \begin{cases} g(x, y_{max}, z) = 0 \\ h(x, y_{max}, z) = 0 \end{cases}$$

$$\begin{cases} g(x, y, z_{min}) = 0 \\ h(x, y, z_{min}) = 0 \end{cases} \quad \begin{cases} g(x, y, z_{max}) = 0 \\ h(x, y, z_{max}) = 0 \end{cases}$$

Consequently, to compute the desired box, our algorithm consists in refining successively the initial box and checking the existence of real solutions for the above algebraic systems. The algorithm stops when no real solution is found.

## 5. CONCLUSION

We addressed in this paper some aspects related to the study of the stability of differential time-delay systems with commensurate delays. More precisely, following the approach in common use, based on the computation of the critical pairs, we proposed a new strategy for analyzing the asymptotic behavior at these points, that gets
rid of Puiseux series computations. Instead, our strategy
proceeds by locally approximating at the critical pairs, the
quasi-polynomial by a polynomial, and then by analyzing
the asymptotic behavior at the zeros of this polynomial
seen as specific points of algebraic curve in $\mathbb{R}^3$. The im-
plementation of this approach is the subject of an ongoing
work. In future work, an interesting problem is to inves-
tigate more deeply the impact of numerical errors, due to
the truncation as well as the local approximation at the
critical pairs, on the certification of our result.

REFERENCES

Boussaada, I. and Niculescu, S.I. (2016). Tracking the algebraic multiplicity of crossing imaginary roots for generic quasipolynomials: A vandermonde-based approach. *IEEE Transactions on Automatic Control*, 61(6), 1601–1606.

Bouzidi, Y., Poteaux, A., and Quadrat, A. (2016). Computer algebra methods for the stability analysis of differential systems with commensurate time-delays. *IFAC-PapersOnLine*, 49(10), 194–199.

Gu, K., Kharitonov, V.L., and Chen, J. (2003). *Stability of Time-Delay Systems*. Birkhäuser.

Li, X.G., Niculescu, S.I., and Çela, A. (2015). *Analytic Curve Frequency-Sweeping Stability Tests for Systems with Commensurate Delays*. Springer.

Li, X.G., Niculescu, S.I., Çela, A., Wang, H.H., and Cai, T.Y. (2013). On computing puiseux series for multiple imaginary characteristic roots of lti systems with commensurate delays. *IEEE Transactions on Automatic Control*, 58(5), 1338–1343.

Marshall, J.E., Gorecki, H., Korytowski, A., and Wal-
ton, K. (1992). *Time-Delay. Systems: Stability and
Performance Criteria with Applications*. London: Ellis
Horwood.

Moore, R.E. (1979). *Methods and applications of interval analysis*, volume 2. Siam.

Niculescu, S.I. (2001). *Delay Effects on Stability: A Robust Control Approach*. Number 269 in Lecture Notes in Control and Information Sciences. Springer.

Olgac, N. and Sipahi, R. (2002). An exact method for the stability analysis of time-delayed linear time-invariant (LTI) systems. *IEEE Transactions on Automatic Control*, 47(5), 793–796.

Rouillier, F. (1999). Solving zero-dimensional systems through the rational univariate representation. *J. of Applicable Algebra in Engineering, Communication and Computing*, 9(5), 433–461.