# A Review of Some Subtleties of Practical Relevance for Time-Delay Systems of Neutral Type

Authors: Keqin Gu
Journal: ISRN Applied Mathematics
 (2012)

---

## Abstract

This paper reviews some subtleties in time-delay systems of neutral type that are believed to be of particular relevance in practice. Both traditional formulation and the coupled differential-difference equation formulation are used. The discontinuity of the spectrum as a function of delays is discussed. Conditions to guarantee stability under small parameter variations are given. A number of subjects that have been discussed in the literature, often using different methods, are reviewed to illustrate some fundamental concepts. These include systems with small delays, the sensitivity of Smith predictor to small delay mismatch, and the discrete implementation of distributed-delay feedback control. The framework prsented in this paper makes it possible to provide simpler formulation and strengthen, generalize, or provide alternative interpretation of the existing results.


---

## Body

## 1. Introduction

Time-delay systems of neutral type may be used to model a system without feedback control, such as a lossless transmission line [1]. It may also be an appropriate model for systems under feedback control, such as discrete implementation of distributed-delay feedback control [2]. Compared with systems of retarded type, analyzing systems of neutral type involves a number of rather subtle points. A thorough understanding of these subtleties may be crucial to understanding some rather surprising phenomena of practical importance.One of such phenomena is the drastic change of stability under arbitrarily small delay deviation from the nominal value. This phenomenon has been documented for decades in the control systems [3] and is known under various circumstances as practical stability [4–8], w-stability [9, 10], and robust stability under small delay [11] in the control systems circle and known as strong stability [12–14] in the more mathematical circle. Some simpler problems, such as the practical stability problem of Smith predictor under small delay mismatch [5–7], have been understood well using other methods, such as Nyquest stability criterion. However, some other problems, such as the mechanism of instability of discrete approximation of distributed-delay feedback control [10, 15–18], require more systematic approach to be understood completely.The purpose of this paper is to summarize the basic theoretical background of the stability of time-delay systems of neutral type that the author believes to be of particular relevance, and to illustrate how these theoretical results can be applied in some practical settings. In addition to the traditional formulation, the coupled differential-difference equation formulation will be used extensively.It is not the purpose of this paper to conduct a comprehensive survey of general time-delay systems. For that purpose, the readers are referred to the papers by Niculescu et al. [19], Kharitonov [20], Gu and Niculescu [21], Richard [17], Normey-Rico and Camacho [22], and Sipahi et al. [23]. Many books are devoted to different aspects of time-delay systems. Bellman and Cooke’s book [24] provides a very readable introduction with a thorough treatment of the distribution of characteristic roots. Hale and Verduyn Lunel’s book [12] provides a comprehensive mathematical coverage of time-delay systems. The book by Diekmann et al. [25] is devoted to retarded systems only, but with more thorough treatment. The book by Kolmanovskii and Myshkis [26] gives a comprehensive coverage that may be more accessible to readers with limited mathematical background. This book also contains numerous practical examples in different fields of science and engineering. The book by Malek-Zavarei and Jamshidi [27] is also a useful reference for some topics of interests that are accessible to readers with limited mathematical background. The short book by  Ėl’sgol’ts [28] covers many results that may be difficult to find in other books. A more comprehensive coverage can be found in the book by Ėl’sgol’ts and Norkin [29]. For treatments with more modern approaches and concentration on stability problems, the readers are referred to the books by Gu et al. [30] and Niculescu [31]. The books by Górecki et al. [32] and by  Kolmanovskiĭ and Shaĭkhet [33] also cover optimal control problems, among other topics. Stépán’s book [34] is devoted mainly to frequency domain approach, and the book contains a rich collection of stability charts. The book by Michiels and Niculescu [35] also uses frequency domain approach, but it also covers control problems. Mao’s book [36] covers stochastic systems. Stochastic systems are also covered in [26, 37], but the coverage is less thorough in terms of theoretical development. The book by Wu et al. [38] uses mainly simple Lyapunov-Krasovskii functional approach. While such approaches are rather conservative in general, the book also contains some more updated results on time-varying delays and linear matrix inequality manipulations. The book by Normey-Rico and Camacho [39] is devoted to process control. Erneux’s short book [40] is more concerned with a dynamical systems point of view. Kuang’s book [41] covers the applications to population dynamics along with basic theories. A comprehensive coverage of numerical methods may be found in Bellen and Zennaro [42]. A number of books are devoted to systems with only input or output delays. Zhong’s book [43] is devoted to robust control. The book by Zhang and Xie [44] is devoted to optimal control and parameter estimation. Krstic’s book [45] extends many results to nonlinear and distributed parameter systems. There are also a few books with much narrower scopes. For example, the book by Silva et al. [46] is devoted to obtaining stabilizing parameter regions of PID control of systems with an input delay. The book by Wang et al. [47] is devoted to finite spectrum assignment. Finally, a number of books, such as Curtain and Zwart [48] and Bensoussan et al. [49], also treat time-delay systems in the framework of infinite-dimensional systems.This paper is organized as follows. Section2 covers the traditional formulation of functional-differential equations and differential-difference equations. Section 3 covers the more recent coupled differential-difference equation formulation. Section 4 provides various stability definitions. Section 5 reviews the distributions of characteristic roots and their relations to stability. Section 6 covers the stability issues of difference equations. Especially, the discontinuity of its spectrum is observed, and conditions of stability under small delay deviations are presented. Section 7 discusses the continuity issue in the stability analysis of the complete time-delay system of neutral type and summarizes the principles and guidelines for practical stability analysis. Section 8 considers the retention of stability under small delays. Section 9 covers the sensitivity of the stability of the Smith  predictor subjected to a small delay deviation. Section 10 presents two classes of time-delay systems that appear in the form of neutral type but behave as retarded type systems. Section 11 gives a comprehensive coverage of the limitations of discrete implementation of distributed-delay feedback control. The effect of dependence of delay errors are examined. Section 12 covers the marginal case when an infinite number of left half plane characteristic roots approach the imaginary axis. Section 13 concludes the paper.


## 2. Functional-Differential Equations

Traditionally, a linear time-delay system of neutral type is described by the following functional-differential equation:(2.1)ddt[D(t)xt]=L(t)xt,
where D and L are ℝn-valued linear operators for each given t; xt is defined as
(2.2)xt(θ)=x(t+θ),θ∈[-r,0],
and r is the maximum delay. In other words, the notation xt represents a shift of the time function x by the amount t, and a restriction to the interval [-r,0]. In general, the linear operators are in the form of
(2.3)D(t)ϕ=ϕ(0)-∫-r0‍dθ[μ(t,θ)]ϕ(θ),(2.4)L(t)ϕ=∫-r0‍dθ[η(t,θ)]ϕ(θ),
where the subscript θ is used to indicate that θ is the integration variable, and μ and η are of bounded variation with respect to θ for each given t. For the problem to be well posed, the integral on the right hand side of (2.3) should be uniformly nonatomic at 0; that is, for any given ε>0, there exists a δ such that the total variation of μ(t,θ) as a function of θ within [-δ,0] is less than ε for any t. The system is reduced to the retarded type if μ=0.It is sufficient in this paper to consider the following special case:(2.5)∫-r0‍dθ[μ(t,θ)]ϕ(θ)=∑k=1K‍Dk(t)ϕ(-rk)+∫-r0‍D(t,θ)ϕ(θ)dθ,(2.6)∫-r0‍dθ[η(t,θ)]ϕ(θ)=∑k=0K‍Ak(t)ϕ(-rk)+∫-r0‍A(t,θ)ϕ(θ)dθ,
where
(2.7)0=r0<r1<⋯<rK=r.
Notice, the summation in (2.5) starts from 1 rather than 0 so that the problem is well posed (i.e., μ is uniformly nonatomic). We will use (2.1) only for notational simplicity. In most cases, we discuss the special case of time-invariant systems with discrete delays; in which case, (2.1) may be written as the following differential-difference equation:
(2.8)x˙(t)-∑k=1K‍Dkx˙(t-rk)=A0x(t)+∑k=1K‍Akx(t-rk),
where Ak∈ℝn×n, Dk∈ℝn×n. For a given initial time t0 the initial condition for (2.1) is given in the form of
(2.9)xt0=ϕ,
where
(2.10)ϕ∈C([-r,0],ℝn).
The initial condition (2.9) may be expressed more explicitly as
(2.11)x(t0+θ)=ϕ(θ),-r≤θ≤0.The basic theory of such systems can be found in the book by Hale and Verduyn Lunel [12] and the references therein. For example, the existence and uniqueness of solutions may be found in [50].


## 3. Coupled Differential-Difference Equations

Time-delay systems of neutral type were initially motivated by some physical systems described by partial differential equations of hyperbolic type with time and space as the independent variables. When one is only interested in certain discrete points on space, the equation can often be reduced to the form of coupled differential-difference equations. A well-known example is the lossless transmission line given by Brayton in [1]. However, similar models date back to as early as 1940s. See, for example, [51] where a steam system is modeled in this form. As will be seen later in this paper, some well-known feedback control methods, such as Smith predictor and discrete implementation of distributed-delay feedback control, also result in coupled differential-difference equations. In other words, time-delay systems of neutral type are often more naturally described by coupled differential-difference equations.A more general description of such systems is coupled differential-functional equations [52]
(3.1)x˙(t)=f(t,x(t),yt),(3.2)y(t)=g(t,x(t),yt),
where the subscript t in yt indicates a shift and restriction of y similar to xt in (2.2). In some literature, the delayed x, xt, is also included in the model [53]. However, it is always possible to transform such a model to the form given in (3.1) and (3.2) by introducing additional variables as discussed in [54]. Let
(3.3)𝒞a={(ψ,ϕ)∣ψ∈ℝn,ϕ∈C([-r,0],ℝm),ϕ(0)=g(t,ψ,ϕ)}.
If f and g satisfy certain continuity conditions, and the initial conditions
(3.4)x(t0)=ψ,y(t0+θ)=ϕ(θ),-r≤θ≤0
satisfy
(3.5)(ψ,ϕ)∈𝒞a,
then
(3.6)(x(t),yt)∈𝒞a,∀t≥t0.In this paper, we are mainly interested in linear time-invariant systems with discrete delays. Such a system can always be described by coupled differential-difference equations of the following form:(3.7)x˙(t)=Ax(t)+∑j=1K‍Bjyj(t-rj),yk(t)=Ckx(t)+∑j=1K‍Dkjyj(t-rj),k=1,2,…K,
where x(t)∈ℝn; yk∈ℝmk; A, Bj, Ck, and Dkj are real matrices of appropriate dimensions. As shown in [54], any linear time-invariant system with multiple discrete delays can be written in the above standard form. In some topics discussed in this paper, such as the discrete implementation of distributed-delay feedback control to be covered in Section 11, it is also important to consider the case when the delays are linear combinations of independent parameters. A transformation may be carried out so that the independent parameters will appear as independent delays. This process is briefly mentioned in the discussion after Corollary  6.1 in Chapter 9 of [12] and described in more detail in [55].The easiest way of arriving at the description (3.7) from a system block diagram is through a process known as “pulling out delays.” This process will be briefly described in Section 8, and it parallels the process of “pulling out uncertainties” that is described in detail by Doyle et al. in [56]. If properly modeled, the “pulling out delay” process should result in a state-space description given by (3.7) with a smaller state space than that of (2.8) or the more general description given by (3.1) and (3.2). Indeed, the initial conditions, which describe the initial state, may be specified as
(3.8)x(t0)=ψ,yk(t0+θ)=ϕk(θ),-rk≤θ<0.
Notice that, yk(σ+θ) for -r≤θ<rkis not needed for the case of rk<r. Therefore, instead of 𝒞a, the state space may be further restricted to
(3.9)𝒞b={(ψ,ϕ1,ϕ2,…,ϕK)|ψ∈ℝn,ϕk∈C([-rk,0],ℝmk),ϕk(0)=Ckψ+∑j=1K‍Dkjϕj(-rj),k=1,2,…,K}.For a time functionz, let z(rk)t denote the shifting of z by t and a restriction to the interval [-rk,0],
(3.10)z(rk)t(θ)=z(t+θ),-rk≤θ≤0.
Then the state at t is
(3.11)(x(t),y1(r1)t,y2(r2)t,…,yK(rK)t).
If the initial conditions (3.8) satisfy
(3.12)(ψ,ϕ1,ϕ2,…,ϕK)∈𝒞b,
then it is not difficult to see that
(3.13)(x(t),y1(r1)t,y2(r2)t,…,yK(rK)t)∈𝒞b,∀t≥t0.Notice that, there is only one delayrk associated with each yk. This “one-channel-one-delay” formulation permits one to obtain a rather simple general solution in terms of fundamental solutions [54] as compared to that for traditional formulation given, for example, by Henry [57]. It is also important to notice that mk are typically much smaller than n in many practical systems, which means significant reduction of computational effort in the stability analysis using Lyapunov-Krasovskii functional approach [58–60]. These facts make it desirable to choose the coupled differential-difference equation model even for systems of retarded type, in which case the matrices Dkj satisfy
(3.14)det(I-DE(s))=1,
where D and E(s) are given in (5.4) in the next section. This is one of the examples in which a system may appear as of neutral type by the form of its description, but it actually behaves as one of retarded type, as to be discussed in Section 10.It is interesting to observe that the model has been known as the “Roesser’s model” and has been studied earlier using frequency domain approaches [61, 62]. The process of “pulling out delays” was also used by Meinsma et al. [11] in the context of studying the stability of systems with small delays.In some cases, some simple substitutions allow one to transform (3.7) to the form of (2.8). Therefore, the coupled differential-difference equations are often considered as an alternative description of (2.8). In general, however, one needs to take derivative of (3.2) in order to write the whole system, described by (3.1) and (3.2) (when the system is linear), in the standard form (2.1). However, in order to make the resulting system equivalent to the original system, it is necessary to constrain the state space to some subspace [63, 64], which causes substantial complication in the analysis. Most early studies concentrate on the description (2.1). An exception is [65] where direct analysis was carried out. In recent years, there has been a substantial interest in direct analysis of coupled differential-difference equations. See, for example, [53, 66–70].


## 4. Stability

It is convenient to usez to represent the state, 𝒞 to represent the state space, and ∥·∥ to represent the norm of the state for all four descriptions given above. Specifically, let |·| refer to the 2-norm of column vectors, then for the system described by (2.1) or (2.8),
(4.1)𝒞=C([-r,0],ℝn),z(t)=xt,∥z(t)∥=max-r≤θ≤0|x(t+θ)|.
For the coupled differential-functional equations described by (3.1) and (3.2),
(4.2)𝒞=𝒞a,z(t)=(x(t),yt),∥z(t)∥=max-r≤θ≤0{|x(t)|,|y(t+θ)|}.
In the case of coupled differential-difference equations (3.7), we have
(4.3)𝒞=𝒞b,z(t)=(x(t),y1(r1)t,y2(r2)t,…,yK(rK)t),∥z(t)∥=max-rk≤θk≤01≤k≤K{|x(t)|,|y1(t+θ1)|,|y2(t+θ2)|,…,|yK(t+θK)|}.
We will also refer to a system described in any of the four descriptions as a time-delay system. Then we may give the following definition of stability.Definition 4.1.
For a time-delay system, the trivial solutionz(t)=0 is said to be stable if for any t0∈ℝ and any ε>0, there exists a δ=δ(t0,ε)>0 such that ∥z(t0)∥<δ implies ∥z(t)∥<ε for all t≥t0. It is said to be asymptotically stable if it is stable, and for any t0∈ℝ, there exists a δa=δa(t0) such that ∥z(t0)∥<δa implies limt→∞z(t)=0. It is said to be uniformly stable if it is stable and δ(t0,ε) can be chosen independently of t0. It is uniformly asymptotically stable if it is uniformly stable, and there exists a δa>0 such that for any η>0, there exists a T=T(δa,η), such that ∥z(t0)∥<δa implies ∥z(t)∥<η for t≥t0+T and t0∈ℝ. It is exponentially stable if there exist an M>0 and an α>0 such that
(4.4)∥z(t)∥≤M∥z(t0)∥e-α(t-t0).It is noted that exponential stability as defined above is also known asuniform exponential stability in some literature (see, e.g., [71]). For a linear system, it is well known that uniform asymptotic stability is equivalent to exponential stability [13, 52]. On the other hand, even for linear time-invariant systems, it is possible for a system to be asymptotically stable, but not exponentially stable, if the system is of neutral type [71–74].For linear time-invariant systems, the stability is closely related to the characteristic roots, which will be reviewed in the next section.


## 5. Characteristic Roots

The characteristic equation for the system described by (2.8) is
(5.1)Δ(s)=det(sI-A0-∑k=1K‍e-rks(sDk+Ak))=0.
Similarly, the characteristic equation for the system described by (3.7) is
(5.2)Δ(s)=det(sIn-A-BE(s)-CIm-DE(s))=0,
where
(5.3)m=∑k=1K‍mk,B=(B1B2⋯BK),C=(C1C2⋮CK),(5.4)D=(D11D12⋯D1KD21D22⋯D2K⋮⋮⋱⋮DK1DK2⋯DKK),E(s)=diag(e-r1sIm1e-r2sIm2⋯e-rKsImK).
The solutions to the characteristic equation will be known as the characteristic roots. Let
(5.5)σ-=sup{Re(s)∣Δ(s)=0}.
It is known that σ- is finite for time-delay systems of neutral type. The system is exponentially stable if σ-<0. Indeed, the system trajectories in this case can be bounded by [57, 64]
(5.6)∥z(t)∥≤M∥z(0)∥eσt
for any σ>σ-.Ifσ->0, then there exists at least one trajectory that grows exponentially, and the system is obviously unstable. This case is sometimes known as “exponentially unstable” in order to distinguish it from the case of polynomial growth associated with some cases of σ-=0. The general case for σ-=0 is rather complicated and will be discussed in Section 12 later on.In general, a time-delay system has an infinite number of characteristic roots. However, asΔ(s) is an entire function, there can only be a finite number of characteristic roots within any bounded domain [75]. These characteristic roots form root chains that are rather easy to describe [24]. For the characteristic equation (5.1) or (5.2), there are two types of root chains.The first type isretarded chains. In a retarded chain, the characteristic roots fall in the region
(5.7)|Re(s+μlogs)|<c
for sufficiently large |s|. Different retarded chains have different values of μ>0 and c>0. As a consequence, there may only be a finite number of roots on the right of the vertical line Re(s)=α in the complex plane for any given α.The other type isneutral chains. In a neutral chain, the characteristic roots are bounded by two vertical lines
(5.8)c1≤Re(s)≤c2,
for sufficiently large roots. The positions of such vertical lines are determined by the difference equation associated with the system. Specifically, in the case of (2.8), the associated difference equation is
(5.9)x(t)=∑k=1K‍Dkx(t-rk),
with the corresponding characteristic equation
(5.10)Δ0(s)=det(I-∑k=1K‍e-rksDk)=0.
The difference equation associated with (3.7) is
(5.11)yk(t)=∑j=1K‍Dkjyj(t-rj),k=1,2,…,K,
with the characteristic equation
(5.12)Δ0(s)=det(I-DE(s))=0.
If s0 is a solution of (5.10) or (5.12), then there is a series of corresponding characteristic roots of (5.1) or (5.2), sk, k=1,2,…, |sk|→∞, that approach the vertical line [57]
(5.13)Re(sk)→Re(s0).
Therefore, we can state the following.Theorem 5.1.
A necessary condition for the time-delay system to be exponentially stable is(5.14)σ-0=sup{Re(s)∣Δ0(s)=0}<0
for the associated difference equation. If σ-0>0, then the time-delay system is exponentially unstable.


## 6. Stability of Difference Equations

From the discussions in the previous section, it is important to understand the stability problem of difference equations (5.10) and (5.12). To emphasize the fact that the solution of such equations are defined on t∈ℝ rather than on discrete time, they are known as difference equations of continuous time.Let(6.1)Z={Re(s)∣Δ0(s)=0}.
Obviously,
(6.2)σ-0=maxZ-.
Similar to differential-difference equations, a difference equation is exponentially stable if σ-0<0. If σ-0>0, then the difference equation is exponentially unstable [57].An important concept in understanding the stability problem of difference equations is rational independence. The real numbersr1,r2,…,rK are said to be rationally independent if any rational (or equivalently, integer) combination may vanish only if all the coefficients vanish; that is,
(6.3)∑k=1K‍αkrk=0,αkrational(orequivalently,integers),
may be satisfied only if αk=0 for all 1≤k≤K. If r1,r2,…,rK are not rationally independent, then they arerationally dependent, in which case we may always find a group of rationally independent numbers τ1,τ2,…,τN, N<K, such that r1,r2,…,rK may be expressed as rational (or integer) combinations of τ1,τ2,…,τN. In this case, τ1,τ2,…,τN are known as a rationally independent basis. The number N is independent of the choice of the rationally independent basis. When N=1, then r1,r2,…,rK may all be expressed as integer multiples of a single number, and they are known as commensurate (or commensurable in some literature). If N≥2, then they areincommensurate (or  uncommensurable in some literature). For any given r1*,r2*,…,rK* and an arbitrarily small ε>0, we may always find rationally independent r1,r2,…,rK that satisfy
(6.4)|rk-rk*|<ε,1≤k≤K.Henry [57] obtained the set Z- in the general case, that is, when the delays are integer combinations of rationally independent basis. He also observed that σ-0 is in general a discontinuous function of the delays r1,r2,…,rK. Such discontinuity was also observed and illustrated by examples by Melvin in [76].In order to fully comprehend the issue, let’s concentrate on the difference equation (5.9) for the moment. Hale [77] formulated the conditions for stability under arbitrarily small delay variations. This formulation leads to the following theorem.Theorem 6.1.
The following four statements are equivalent.(i)
Equation (5.9) is exponentially stable for a fixed set of rationally independent delays r1>0,r2>0,…,rK>0.(ii)
For a given nominal delaysr10>0,r20>0,…,rK0>0, and a small ε>0, (5.9) is exponentially stable for all delays r1,r2,…,rK that satisfy
(6.5)|rk-rk0|<ε.(iii)
Equation (5.9) is exponentially stable for arbitrary positive delays r1>0,r2>0,…,rK>0.(iv)
The matricesD1,D2,…,DK satisfy
(6.6)sup0≤θk≤2πk=1,2,…,Kρ(∑k=1K‍eiθkDk)<1,
where ρ(·) is the spectrum radius of the matrix concerned.
If on the other hand,(6.7)sup0≤θk≤2πk=1,2,…,Kρ(∑k=1K‍eiθkDk)>1,
then (5.9) with any fixed rationally independent delays r1>0,r2>0,…,rK>0 is exponentially unstable.The first part of theorem, that is, the equivalence of the four statements, may be found in [12, Theorem  6.1 of Chapter 9]. The last part may be found in [13] with a new proof about the equivalence between (i) and (iv). In practice, there are always errors in estimating or setting delays. If the errors of different delays vary independently, then the above theorem applies. The equivalence of statements (ii) and (iii) is very disquieting; as far as the robust stability is concerned, there is no difference between the case where the delays vary within an arbitrarily small range (often known as practical stability or local strong stability) or the case where the delays are allowed to assume any positive values (delay independent stability, or stability independent of delays). This discontinuity is indeed at the root of many surprising phenomena in many systems with delays.Checking condition (6.6) is not easy in general. A practically computable condition is given by Carvalho in [78] in the form of a linear matrix inequality, which was motivated by Lyapunov functional formulation on the H2 norm. As indicated by Boyd et al. [79], efficient numerical methods based on interior point algorithm are available to solve such linear matrix inequalities. The following condition, which is equivalent to one in [78], is from [80].Proposition 6.2.
The condition (6.6) is satisfied if there exist symmetric positive definite matrices Sk, k=1,2,…,K such that
(6.8)DT∑k=1K‍SkD-diag(S1S2⋯SK)<0,
where
(6.9)D=(D1D2⋯DK).In (6.8), “<0” is used to indicate that the matrix on the left hand side is symmetric negative definite. Similarly, “>” will be used to denote positive definiteness. If the difference equation (5.9) is a scalar equation, then Dk, k=1,2,…,K are scalars, and (6.6) is reduced to
(6.10)∑k=1K‍|Dk|<1.If we apply the robust stability condition (6.6) to the system described by (5.11), we may conclude the following.Corollary 6.3.
The system described by (5.11) is exponentially stable for all rk,  |rk-rk0|<ε, k=1,2,…K if and only if
(6.11)ρ0=Δsupδk∈ℂ,|δk|=11≤k≤Kρ(DE(δ))<1,
where
(6.12)E(δ)=diag(δ1Im1δ2Im2⋯δKImK).Of course, the above is still valid if the delays are allowed to assume any positive values in view of Theorem6.1. If we relax the constraint |δk|=1 in (6.11) to |δk|≤1, we obtain
(6.13)supδk∈ℂ,|δk|≤11≤k≤Kρ(DE(δ))≥ρ0.
Those who are familiar with the structured singular value problem (see [56, 81–83]) may recognize that the left hand side of (6.13) is equal to the structured singular value of the matrix D under the structure described by the matrix E in (6.12) (multiple scalar blocks),
(6.14)μ(D)=1min{ρ[E(δ)]∣δk∈ℂ,k=1,2,…,K,det(I-DE(δ))=0}.
Therefore,
(6.15)ρ0≤μ(D),
and (6.11) can be guaranteed by
(6.16)μ(D)<1.It is well known that calculation of structured singular value is not easy. A sufficient condition in the form of linear matrix inequality is given below [82].Corollary 6.4.
The condition (6.16), and therefore (6.11), is satisfied if there exist Sk∈ℝmk×mk,
(6.17)Sk=SkT>0,k=1,2,…,K,
such that
(6.18)S-DTSD>0,
where
(6.19)S=diag(S1S2⋯Sk).The above condition may also be derived directly from (6.8) [54]. On the other hand, it is interesting to point out that (6.6) may also be guaranteed by (6.16) with
(6.20)D=(D1D2⋯DK-1DKIn0⋯000In⋯00⋮⋮⋱⋮⋮00⋯In0)
and the uncertainty structure
(6.21)E(δ)=diag(δ1Inδ2In⋯δKIn).The linear matrix inequality form of stability conditions (6.18) and (6.8) are very useful in formulating Lyapunov-Krasovskii functional stability conditions of linear time-delay systems of neutral type [54, 80]. It is also interesting that another condition given by Fridman [84] is also equivalent to (6.18) as shown in the appendix of [58]. For the stability problem of nonlinear difference equation of continuous time using Lyapunov method, an interesting study is given by Pepe [85]. It is noticed that the criteria is not as tight as (6.18) if applied to linear systems. See also the recent book by Shaikhet [86] for a more comprehensive study.


## 7. Continuity Issue and Practical Stability

Continuous dependence of characteristic roots on the system parameters is the basis of many important techniques used in the stability analysis. Examples of these techniques include root-locus [87] and D-decomposition [88] (also known as D-partition or D-subdivision). Indeed, it is well known that if the leading coefficient of a polynomial does not vanish, the roots of the polynomial depend on the coefficients continuously. When the root concerned is simple, then it is an analytical function of coefficients. Even around a multiple root, a Puiseux series around the nominal value is possible (see [89] and Part II, Chapter 5 of [90]). An example of discontinuity due to the vanishing leading coefficient is
(7.1)εs2+s+1=0.
The two roots for nonzero ε are
(7.2)s1=-1+1-4ε2ε=-1+ε+o(ε),s2=-1-1-4ε2ε=-1+1ε+o(1).
While s1 is a continuous function of ε, s2 is obviously discontinuous at ε=0. However, such points are rather easy to discover. An example of discontinuity caused by vanishing leading coefficient in time-delay systems is given in [7]. The situation for time-delay systems also has some similarity with this example; at the critical parameter values, some characteristic roots may have discontinuity, while other roots change continuously with the parameters.For time-delay systems of retarded type, although new roots may suddenly appear with infinite magnitudes near some parameter values, such new roots always appear at far left of the complex plane (i.e., with-∞ real parts) and do not affect the stability analysis. Therefore, continuous dependence of characteristic roots on the system parameters is widely used in stability analysis. One of such techniques is again D-decomposition that identifies the parameter values that correspond to the presence of imaginary roots. These values divide the parameter space into regions with fixed number of right half plane roots, from which the stable parameter regions may be easily obtained [29, 91]. Especially, delays may also be used as parameters, which is important as the case of zero delay is reduced to a polynomial equation, which is rather easy to analyze. For the single delay case, see [92, 93], and the polynomial coefficient part of [94]. For multiple commensurate delay case, see [95]. For multiple delay case, see [96–99]. Another interesting method that is matrix based is formulated by Chen et al. [100]. Indeed, the book by Stépán [34] contains many practical examples of such analysis. The more recent book by Michiels and Niculescu [35] also contains some more  recent such analysis.Attempts have also been made in recent years to extendD-decomposition method to time-delay systems of neutral type. See, for example, [101, 102]. Obviously, such analysis includes the systems of retarded type as a special case and is more general. However, as σ- defined in (5.5) may be discontinuous with respect to the delays for such systems, the theoretical basis needs to be carefully examined. Indeed, Cooke and van den Driessche [94] studied the characteristic equation of the form
(7.3)P(s)+Q(s)e-rs=0,
with P(s) a higher order polynomial than the polynomial Q(s). This obviously represents a time-delay system of retarded type. However, in an attempt to extend the result to the more general case of analytical functions, they inadvertently included the possibility of systems of neutral and even advanced type, rendering the result invalid as was shown by Boese [103] and acknowledged by Cooke [104].For time-delay systems of neutral type, as the parameters change, there is a possibility of sudden appearance of an infinite number of characteristic roots with positive real parts without going through the imaginary axis. This can be traced back to the discontinuity of the spectrum of the associated difference equation. Asσ-0 defined in (5.14) for a difference equation is in general discontinuous with respect to delays, Theorem 5.1 and the discussions preceding this theorem indicate that σ- defined in (5.5) for the complete system is also a discontinuous function of delays in general. This renders the stability analysis rather complicated. However, if we restrict ourselves to robust stability under arbitrarily small deviations of parameters, the situation is much better. We will use the term in [4–7] and call such robust stability as practical stability, as defined below.Definition 7.1.
Consider a system with parameterα. The system is said to be practically stable at α=α0 if there exists an ε>0 such that for any permissible α with ∥α-α0∥<ε, the system remains exponentially stable.The parameter above should be interpreted as a vector. Therefore, the case of multiple parameters is also covered. Obviously, the main interest in discussing the practical stability problem is for the parameter vector to include some delays. Thepermissible deviation should be specifically defined according to the specific problem. For example, at zero delay, we should define a permissible delay as positive in order to avoid creating an inherently unstable time-delay system of advanced type. In some cases, deviation of different components of parameter vector α may be constrained to satisfy certain linear relations. As shown in [55], systems with delays that are integer combinations of independent parameters may always be transformed to ones with the independent parameters as delays. This process will also be illustrated in the problem of discrete implementation of distributed-delay feedback control to be covered in Section 11. Therefore, we will only consider the case of independent parameter deviation in this section.LetΔ0(s)=0 be the characteristic equation of the difference equation that depends on the delays rk,k=1,2,…,K. Then σ-0 defined in (5.14) is in general a discontinuous function of rk and may be denoted as σ-0(r1,r2,…,rK). Define
(7.4)σ-0p(r1,r2,…,rK)=limε→0+max∥r^k-rk∥<ε1≤k≤Kσ-0(r^1,r^2,…,r^K).
Then, similar to the case discussed in [57], σ-0p(r1,r2,…,rK) is a continuous function of the delays rk. Obviously, practical stability requires
(7.5)σ-0p(r1,r2,…,rK)<0,
which is not the case for conventional stability. Therefore, practical stability analysis is much easier in terms of continuity. As all the parameters are subject to errors, requiring practical stability is also essential in practice.We may also handle the spectrum of the overall system in a similar manner and define(7.6)σ-p(r1,r2,…,rK)=limε→0+max∥r^k-rk∥<ε1≤k≤Kσ-(r^1,r^2,…,r^K).
However, it is typically more convenient to consider the spectrum related to the difference equation and the remaining characteristic roots separately. Indeed, if the delays r1,r2,…,rK are subject to independent errors, then a small deviation may render them to be rationally independent. In this case, the conditions for (7.5) can be written in a form that only depends on the coefficient matrices, and independent of the delays, as shown in Section 6. If the condition (7.5) is satisfied, then all the characteristic roots that satisfy
(7.7)Re(s)>-ε
for some small ε>0 are continuous functions of system parameters. Therefore, a continuity argument can be used for both the quantity σ-0p and those characteristic roots that satisfy (7.7).Based on the above discussion, a common practice in usingD-decomposition method for stability analysis of time-delay systems of neutral type is to guarantee the satisfaction of (7.5), and keep track of the number of right half plane characteristic roots as the parameter changes. For example, to repair the results of [94], Boese [103] suggested restricting to neutral type systems that satisfy certain additional inequality ((2.4) in [103]). It can be easily seen that this inequality guarantees the stability of the associated difference equation. Other examples of enforcing this condition are Lemma  12 in [105], (2.3) in [106], Assumption  3 in [107], and Assumption  3 in [108].


## 8. Small Delays

An immediate application of the theory developed so far is the stability of a stable delay-free system when it is subjected to small delays. Traditionally, a system with the possibility of instability under small delays is classified as “not well-posed” [3]. For a period, there had been substantial interest on this issue. This also includes infinite-dimensional systems described by partial differential equations [109, 110]. Some studies on infinite-dimensional systems [111] have very similar formulations and conclusions with finite-dimensional systems [11]. The discussions in this paper will be restricted to finite-dimensional systems.Let a linear system ofnth order be exponentially stable when it does not contain any delay. Let there be K components of the system that may be subjected to small delays r1,r2,…,rK. Let the kth delay component be mk-dimensional. We want to study if there exists a small ε>0, such that the system remains exponentially stable for all rk∈(0,ε).This problem can be written in a standard form if we use the process of “pulling out delays.” The process begins with removing all the delay elements. The output of each delay becomes an input to the remaining part of the system. Similarly, the input to each delay becomes an output to the remaining part of the system. The part of system with all the delays removed may be written in the standard state-space form as(8.1)x˙(t)=Ax(t)+∑j=1K‍Bjuj(t),yk(t)=Ckx(t)+∑j=1K‍Dkjuj(t),k=1,2,…,K,
where x(t)∈ℝn is the state variable, the input uk(t)∈ℝmk was the output of kth delay element, and the output yk(t)∈ℝmk was the input to the kth delay element. Reconnect the delay elements in their original locations. This results in
(8.2)uk(t)=yk(t-rk),k=1,2,…,K.
Equations (8.1) and (8.2) now completely describe the original system. A substitution of (8.1) by (8.2) yields the final description of the system
(8.3)x˙(t)=Ax(t)+∑j=1K‍Bjyj(t-rj),(8.4)yk(t)=Ckx(t)+∑j=1K‍Dkjyj(t-rj),k=1,2,…K,
which is in the standard form of coupled differential-difference equations.Whenrk=0 for all k, there are n characteristic roots for the system. These characteristic roots are continuous functions of delays and will remain on the left half of the complex plane when ε is small. However, as delays increase from zero, an infinite number of new characteristic roots appear. The locations of these roots depend on the associated difference equation. Therefore, the only condition that needs to be checked to guarantee exponential stability is the stability of the associated difference equation (8.4). The necessary and sufficient condition for this stability is the satisfaction of (6.11). A sufficient condition is the satisfaction of (6.16) with the uncertainty structure defined by (6.12).When all the delays are single-input single-output, Meinsma et al. [11] showed that the condition (6.12) is necessary and sufficient in the sense of input-output stability. It can be shown that this is also necessary and sufficient in the sense of exponential stability. To do so, we first need the following lemma.Lemma 8.1.
LetG∈ℂn×n be partitioned to
(8.5)G=(HP),
where H is a column vector. Then
(8.6)supδ∈ℂ,|δ|=1ρ(δHP)=supδ∈ℂ,|δ|≤1ρ(δHP).Proof.
Let(8.7)Δδ(s)=det[sI-(δHP)].
Then Δδ(s) is affine with δ for each fixed s. Let the quantity on the left hand side of (8.6) be ρ0. For a small ε>0, form a contour
(8.8)Γ={(ρ0+ε)eiθ∣θ∈[0,2π]}.
As s assumes a fixed point on this contour, the set
(8.9)Cs={Δδ(s)∣δ∈ℂ,|δ|=1}
forms a circle. This circle cannot enclose the origin. (Otherwise, for large real α, since Δδ(αs) is dominated by (αs)n, the origin must be outside of the curve {Δδ(αs)∣δ∈ℂ,|δ|=1}. By continuity, there must exists a α>1 and |δ|=1 to satisfy Δδ(αs)=0. But this implies
(8.10)supδ∈ℂ,|δ|=1ρ(δHP)≥α(ρ0+ε)>ρ0,
a contradiction.) The set
(8.11)Ds={Δδ(s)∣δ∈ℂ,|δ|≤1}
consists of all the points on and inside of Cs. As s goes around Γ, the Argument principle [75] and
(8.12)supδ∈ℂ,|δ|=1ρ(δHP)<ρ0+ε
imply that each point on Cs goes around the origin exactly n times. But this also means that each point on Ds also goes around the origin n times. Therefore, for any |δ|≤1, all n roots of Δδ(s) are inside Γ; that is,
(8.13)supδ∈ℂ,|δ|≤1ρ(δHP)<ρ0+ε.
As ε is arbitrary, we can conclude that
(8.14)supδ∈ℂ,|δ|≤1ρ(δHP)≤ρ0.
On the other hand, it is obvious that
(8.15)supδ∈ℂ,|δ|≤1ρ(δHP)≥supδ∈ℂ,|δ|=1ρ(δHP)=ρ0.
This proves (8.6).Apply the above lemma forG=D, and δ=δk with δj, j≠k fixed, and repeat the process for each k; we arrive at the following proposition.Proposition 8.2.
For the system (8.3) and (8.4) with all yk(t)∈ℝ, let
(8.16)D=(D11D12⋯D1KD21D22⋯D2K⋮⋮⋱⋮DK1DK2⋯DKK),(8.17)E(δ)=diag(δ1δ2⋯δK).
Then,
(8.18)supδk∈ℂ,|δk|=11≤k≤Kρ(DE(δ))=μ(D),
where the uncertainty structure is defined by E(δ).The conclusion is obvious from the above proposition and is stated below.Corollary 8.3.
If the system described by (8.3) and (8.4) is exponentially stable for rk=0, k=1,2,…,K, and yk(t)∈ℝ, k=1,2,…,K. Then, there exists an ε>0 such that the system is exponentially stable for all rk∈(0,ε) if and only if μ(D)<1, where the uncertainty structure is defined in (8.17).Another interesting result regarding small delay is on the stabilization of difference equations. Theoretically, such a stabilization is possible only if a derivative is used [13, 112]. Logemann and Townley [113] showed that if the open-loop difference equation is unstable, then the stabilized system can be destabilized by an arbitrarily small input delay. The instability problem caused by delay mismatch in Smith predictor to be discussed in the next section is of similar nature.


## 9. Delay Sensitivity of Smith Predictor

Smith predictor, proposed in [114], is a well-known control method for processes with a delay. This section is devoted to illustrating how the theory developed so far can be used to solve the practical stability problem of Smith predictor delay mismatch. Although this problem has been solved completely in the literature, the method presented here is much simpler. The structure of Smith predictor is shown in Figure 1, where Gp(s)e-Ts is the plant to be controlled, Gm(s) is an estimated model of the plant without delay, Gc(s) is a control designed based on the estimated model, and τ is the estimated delay.Figure 1
Smith predictor.In the ideal case, the model should be exactly equal to the plant,(9.1)Gm(s)=Gp(s),τ=T,
in which case, if the system is single-input single-output, the closed-loop transfer function can be calculated as
(9.2)Gcl(s)=Gc(s)Gp(s)e-Ts1+Gc(s)Gp(s).
The above equation shows why the Smith predictor is so attractive. One can design the controller based on the plant model without delay. The resulting closed-loop transfer function under the ideal case has exactly the same dynamics as if the control Gc(s) is applied to a system without delay Gp(s), and the delay may be considered as taken out of the feedback loop. Numerous extensions and analyses have been made on Smith predictor over the years; see [5, 115–117] and the references therein.In the single-input single-output case, it was shown by Palmor [5] that the closed-loop system may become unstable under arbitrarily small deviation of the estimated delay τ from the plant delay T. A condition for this not to happen, which Palmor called practical stability, is derived using Nyquist stability criterion. In the following, it will be shown that this phenomenon is closely related to the discontinuity of the spectrum of the associated difference equation.Let the state-space description ofe-TsGp(s) be
(9.3)x˙p(t)=Apxp(t)+Bpu(t-T),(9.4)y(t)=Cpxp(t)+Dpu(t-T),
that of (1-e-τs)Gm(s) be
(9.5)x˙m(t)=Amxm(t)+Bm[u(t)-u(t-τ)],(9.6)v(t)=Cmxm(t)+Dm[u(t)-u(t-τ)],
and that of Gc(s) be
(9.7)x˙c(t)=Acxc(t)+Bcw(t),(9.8)u(t)=Ccxc(t)+Dcw(t).
Then, the closed-loop system can be described by the above equations with the additional constraint
(9.9)w(t)=r(t)-v(t)-y(t).
This is in the form of coupled differential-difference equations. The associated difference equations are (9.4), (9.6), (9.8), and (9.9), where r(t), xp(t), xm(t), and xc(t) are the inputs to the difference equations. Eliminating the variables y(t), v(t), and w(t) in these equations yields
(9.10)u(t)=(I+DcDm)-1[DcDmu(t-τ)-DcDpu(t-T)]+(I+DcDm)-1[Ccxc(t)-DcCmxm(t)-DcCpxp(t)+Dcr(t)].
Applying the condition (6.6) to the above, we may conclude the following.Theorem 9.1.
If the system described by (9.3) to (9.9) is exponentially stable for τ=T, then there exists an ε>0 such that the system remains exponentially stable for any τ∈(T-ε,T+ε); that is, the system is practically stable under small delay deviation, if and only if
(9.11)sup0≤θ≤2πρ[(I+DcDm)-1Dc(eiθDm+Dp)]<1.In the above, we have used the fact that multiplyinge-iθ2 does not change the spectrum radius of a matrix, and
(9.12){ei(θ1-θ2)∣0≤θk≤2π,k=1,2}={eiθ∣0≤θ≤2π}.
As discussed in Section 6, the quantity on the left hand side of (9.11) may be bounded by a structured singular value. Therefore, the result can be related to the known results given in [6, 7]. We may apply (6.8) to (9.11) to obtain a sufficient condition in the form of linear matrix inequality.IfDp=Dm, then
(9.13)ρ[(I+DcDm)-1Dc(eiθDm+Dp)]=|eiθ+1|ρ[(I+DcDm)-1DcDm].
In this case, (9.11) reduces to
(9.14)ρ[(I+DcDm)-1DcDm]<12
since
(9.15){|eiθ+1|∣0≤θ≤2π}=[0,2].
The condition (9.14) is given in Proposition  12 of [8]. The sufficient part may be found in [7].IfGp(s) is single-input single-output, then all the matrices in (9.11) are scalars, and the condition (9.11) reduces to
(9.16)|Dc1+DcDm|(|Dm|+|Dp|)<1.
If Dm=Dp, then the above reduces to Proposition  4 of [8]. The sufficient part may be found in [5].


## 10. Neutral Systems That Behave as Retarded Systems

It is interesting to note that some time-delay systems expressed in the form of neutral type actually behave like one of retarded type. The first type of such systems have a difference equation that has empty spectrum. Consider the system described by (2.8). If the associated difference equation (5.9) has empty spectrum, that is,
(10.1)Δ0(s)=1,∀s∈ℂ,
then, as shown by Henry in [118], the solutions to (5.9) vanishes for t≥rn for any initial condition. As shown by Stoorvogel et al. in [119], we may make a variable transformation
(10.2)z(t)=x(t)-∑k=1K‍Dkx(t-rk).
The condition (10.1) means that the transformation (10.2) is a valid one as x(t) may be expressed by a linear combination of z(t+θ), -nr≤θ≤0. Indeed, a Laplace transform of (10.2) yields
(10.3)Z(s)=H(s)X(s),
where
(10.4)H(s)=I-∑k=1K‍e-rksDk.
This equation can be solved to yield
(10.5)X(s)=[H(s)]-1Z(s).
From
(10.6)Δ0(s)=det[H(s)]=1,
it can be seen that [H(s)]-1 is a polynomial matrix of e-rks with order not higher than n-1. Therefore, x(t) can be expressed as a linear combination of z(t) and its delayed version up to (n-1)r. Equation (2.8) can thus be transformed to a differential-difference equation of retarded type.Similarly, for the coupled differential-difference equations (3.7), if the difference equation (5.11) satisfies (10.1), y(t) can be expressed as a linear combination of Cx and its delayed version, and the system can be expressed as a differential-difference equation of retarded type in x(t).The second type of systems that behave like those of retarded type  involve  neutral distributed delay. Consider the system(10.7)ddt[x(t)-∫-r0‍D(θ)x(t+θ)dθ]=∑k=0K‍Akx(t-rk).
If D(θ) is of bounded variation, then an integration by parts yields
(10.8)ddt∫-r0‍D(θ)x(t+θ)dθ=∫-r0‍D(θ)∂∂θ[x(t+θ)]dθ=D(0)x(t)-D(-r)x(t-r)-∫-r0‍d[D(θ)]x(t+θ).
Using the above, (10.7) is transformed to one of retarded type.More fundamental understanding of the nature of systems of neutral or retarded type can be gained from an abstract formulation. Let the state of the system (10.7) at time t be z(t). Then solutions can be represented by a strongly continuous operator T(t),
(10.9)z(t)=T(t)z(0).
Similarly, the solutions to the associated difference equation (excluding the neutral distributed delay) can be represented as
(10.10)z(t)=TD(t)z(0).
As shown by Henry [57], T(t) can be viewed as a compact perturbation of TD for sufficiently large t. The operator TD(t), if the spectrum is not empty, is not compact. From this point of view, we may classify a time-delay system as of retarded type whenever the solution operator T(t) is compact for sufficiently large t. The difference is also obvious from the spectrum of T(t); while the spectrum of the infinitesimal generator of T(t) contains only eigenvalues, it is not true for the spectrum of T(t) itself. Indeed, for systems of neutral type, T(t) also contains continuous spectrum [57]. Hale and Verduyn Lunel [14] discussed the issue from the point of view of the essential spectrum of T(1) that cannot be changed by compact perturbations. The stability is determined by eigenvalues (which are continuous functions of parameters) of T(1) if and only if the essential spectrum of T(1) is less than 1.On the other hand, it is sometimes desirable to express a systems of retarded type in the form of neutral type. For example, in analyzing additional dynamics due to model transformation [120, 121], a more clear understanding can be obtained by writing it in the form of a differential equation coupled with a functional equation (distributed delay); see [122] and Section  5.3.3 of [30]. The distributed-delay feedback control to be discussed next is another such example. Nonetheless, it is important to keep in mind the fundamental difference of these two types of systems.


## 11. Discrete Implementation of Distributed-Delay Feedback Control

### 11.1. Basic Formulation

Distributed-delay feedback control is an important method to control systems with delays. When the nominal system is unstable, the Smith predictor using the architecture shown in Figure1 involves unstable pole-zero cancellation, and an alternative implementation using distributed-delay feedback control may be used to avoid this problem [115]. Another common method to control systems with input or output delays is finite spectrum assignment [2, 47, 123, 124]. The system with single input-delay and distributed-delay feedback control may be expressed as
(11.1)x˙(t)=Ax(t)+Bu(t-r),(11.2)u(t)=Fx(t)+∫-r0‍G(θ)u(t+θ)dθ.
Exact implementation of control (11.2) is very difficult, if not impossible, due to the integration term. Some literature suggests using discrete delays to approximate the distributed delay in (11.2) [2]. The basic idea is to divide the interval [-r,0] to K small intervals, [θk,θk-1], k=1,2,…,K, with
(11.3)0=θ0>θ1>θ2>⋯>θK=-r,
and to approximate the integration in each small interval [θk,θk-1]. We may choose a point in each interval
(11.4)-rk∈[θk,θk-1].
Common choices are
(11.5)-rk=θk,θk-1,orθk+θk-12.
When all the intervals are sufficiently small, it seems reasonable to approximate any u(t+θ) with θ∈[θk,θk-1] by u(t-rk). Therefore, (11.2) may be approximated by
(11.6)u(t)=Fx(t)+∑k=1K‍Gku(t-rk),
where
(11.7)Gk=∫θkθk-1‍G(θ)dθ.
Less accurate method of obtaining Gk may be used. For example, it is common to use
(11.8)Gk=G(-rk)(θk-1-θk),
leading to a rectangular rule of numerical integration. As was discussed by Zhong [18], (11.7) compares favorably with (11.8) in terms of accuracy, although the difference is not fundamental. We will call (11.6) quasi-rectangular implementation, and any choice of Gk and rk along this line will be known as a rectangular-like implementation. If G(θ)u(t+θ) is Riemann integrable, then for any rectangular-like implementation, the expression on the right hand side of (11.6) converges to that of (11.2) as hmax→0, where
(11.9)hmax=max1≤k≤K(θk-1-θk).Alternatively, we may use[u(t+θk)+u(t+θk-1)]/2 to approximate u(t+θ), θ∈[θk,θk-1]. As a result, (11.6) should be replaced by
(11.10)u(t)=Fx(t)+∑k=0K‍Hku(t-rk),
where
(11.11)H0=12∫θ10‍G(θ)dθ,Hk=12∫θk+1θk-1‍G(θ)dθ,1≤k≤K-1,HK=12∫θKθK-1‍G(θ)dθ,(11.12)rk=-θk,k=0,1,…,K.
Similarly, a less accurate choice is to use
(11.13)H0=-12(θ0-θ1)G(0),Hk=12(θk-1-θk+1)G(θk),1≤k≤K-1,HK=12(θK-1-θK)G(θK),
which becomes the trapezoidal rule of numerical integration. Again, there is no fundamental difference between the more accurate expression (11.11) and the slightly less accurate expression (11.13). However, in this case, the feedback rule (11.10) is not well posed as the nonzero coefficient H0 at r0=0 violates the uniform nonatomic requirement discussed in Section 2. This problem can be avoided if the size of the interval [θ1,0] is sufficiently small, in which case H0 is small, and I-H0 is invertible. We may solve for u(t) in (11.10) to obtain a well-posed implementation rule,
(11.14)u(t)=Cx(t)+∑k=1K‍Dku(t-rk),
where
(11.15)C=(I-H0)-1F,Dk=(I-H0)-1Hk,k=1,2,…,K,rk=-θk.
We will call (11.14) quasi-trapezoidal implementation. In practice, there may also be small deviation of rk. Any of these implementations will be known as a trapezoidal-like implementation. Obviously, if G(θ)u(t+θ) is Riemann integrable, then for any trapezoidal-like implementation, the expression on the right hand side of (11.10) or (11.14) converges to that of (11.2) as hmax→0.For a uniform gridding(11.16)θk=-krK,
it was pointed out in [15, 16] that instability was observed even in numerical simulation when control rules similar to (11.6) or (11.14) are used in place of (11.2) even when hmax becomes very small. The problem was declared by Richard in the survey paper [17] as one of the open problems in the control of time-delay systems. Significant insight has been gained since then [10, 18, 125, 126], and alternative implementation strategies have been proposed [10, 126–128]. This section will only discuss the stability property of such systems and will not discuss improvement of implementation. The material presented here is more general and many results are often more precise.


### 11.2. Fundamental Limitation for General Case

Consider the internal dynamics of the distributed-delay feedback (11.2),
(11.17)u(t)=∫-r0‍G(θ)u(t+θ)dθ.
It is exponentially stable if and only if all the roots of its characteristic equation
(11.18)Δd(s)=det(I-∫-r0‍eθsG(θ)dθ)=0
are on the strict left half plane. A fundamental limitation of discrete implementation is due to this fact.Theorem 11.1.
Let the matrix functionG(θ) be Riemann integrable in [-r,0]. If the internal dynamics of (11.17) is exponentially unstable, then the feedback control system consisting of (11.1) and a feedback control using either a rectangular-like implementation (11.6) or a trapezoidal-like implementation (11.14) is exponentially unstable for a sufficiently small hmax (defined in (11.9)).Proof.
If (11.17) is exponentially unstable, then there exists a s0, Re(s0)>0, such that (11.18) is satisfied for s=s0. The function Δd(s) is an entire function that is not a constant. Therefore, the order of root s0 is finite, and there is an a>0 such that Δd(s)≠0 for any s in the region 0<|s-s0|≤a. We may obviously make
(11.19)a<Re(s0).
Let
(11.20)Γ={s∈ℂ∣|s-s0|=a}.
Since Γ is compact, we may define
(11.21)b=min{|Δd(s)|∣s∈Γ}.
Obviously
(11.22)b>0.
In the compact set Γ, as either rectangular-like implementation or trapezoidal-like implementation can approximate the integration in (11.18) to arbitrary accuracy with sufficiently small hmax, the related quantity Δd(s) may also be approximated to arbitrary accuracy. To be definite, consider the quasi-trapezoidal implementation. There exists a h^, such that for any hmax≤h^,
(11.23)|Δim(s)-Δd(s)|<b2,fors∈Γ,
where
(11.24)Δim(s)=det(I-∑k=1K‍Dkerks).
However, according to Rouché’s theorem [75], the inequality (11.23) implies that Δim(s) has the same number of roots as Δd(s) within the region enclosed by Γ, which is completely on the right half complex plane due to (11.19). Therefore, Δim(s) has at least one right half plane root. As
(11.25)Δim(s)=0
is the characteristic equation of the difference equation
(11.26)u(t)=∑k=1K‍Dku(t-rk),
this difference equation is exponentially unstable. However, according to Theorem 5.1, this also means that the complete system consisting of (11.1) and (11.14) is exponentially unstable.The special case of the above theorem with uniform gridding and finite spectrum assignment can be found in [125]. The above theorem gives a fundamental limitation of any discrete implementation of a distributed-delay feedback control system. We may understand the situation by reversing the role of the plant dynamics (11.1) and the controller dynamics (11.2) or its discrete implementation, such as (11.14), and considering x as the input to the “plant” (11.2) or its  discretized version (11.14). In the case of (11.2), as its solution operator is compact, it is possible for an appropriate (11.1) to make the whole system stable. On the other hand, as the solution operator for (11.14) is not compact, no linear differential equation in the form of (11.1) may stabilize it as its unstable essential spectrum cannot be changed by a compact operator.


### 11.3. Commensurate Delays

The most favorable implementation is using commensurate delays. Imagine if we have a multichanneled device that delays the inputs in all channels byh. Like any device, there is bound to be errors. However, the delays for all channels will be identical. If such a device is available, then we can imagine implementing the feedback law as shown in Figure 2. In this case, when error is considered, the discrete delays become
(11.27)rk=kh,(11.28)h=rK+εh,|εh|<ε.
The proportional relation (11.27) is guaranteed by the structure of the system and does not contain any error. Although the delays still contain errors, they are not independent.Figure 2
Implementation to enforce proportional delays.Introduce the variable(11.29)y^1(t)=(y1(t)y2(t)⋮yK(t))=(u(t)u(t-h)u(t-2h)⋮u(t-(K-1)h)),y^2(t)=yK(t).
Then, (11.1) and (11.14) become
(11.30)x˙(t)=Ax(t)+B^(y^1(t-h)y^2(t-h^)),(11.31)(y^1(t)y^2(t))=C^x(t)+D^(y^1(t-h)y^2(t-h^)),
where
(11.32)h^=r-(K-1)h=rK-(K-1)εh,(11.33)B^=(00⋯0B),C^=(C0⋮0),D^=(D1D2⋯DK-1DK0I0⋯0000I⋯000⋮⋮⋱⋮⋮⋮00⋯I0000⋯I00).
Let
(11.34)D~=(D1D2⋯DK-1DKI0⋯000I⋯00⋮⋮⋱⋮⋮00⋯I0).
Then
(11.35)D^=(D~0En0),
where
(11.36)En=(00⋯0I0).
Although the delays h and h^ are not exactly independent, considering them as independent does not increase the conservatism as the last column block of D^ vanishes. Indeed, under this assumption, the condition for exponential stability of the difference equation (11.31) for all permissible small delay errors becomes
(11.37)maxδk∈ℂ|δk|=1ρ(D^E(δ1,δ2))=maxδ1∈ℂ|δ1|=1ρ(D^E(δ1,0))=ρ(D~)<1,
or(11.38)max{|λ||det(λKI-∑k=1K‍λK-kDk)=0}<1.
It is noticed that the above condition is independent of the delay error εh.To get a sense of how well this implementation works, consider the special case of the example discussed in [10, 15, 16], which is restated below.Example 11.2.
The system with scalar variables(11.39)x˙(t)=x(t)+u(t-r)
is stabilized using the following distributed-delay feedback control using finite spectrum assignment method:
(11.40)u(t)=-2[erx(t)+∫-r0‍e-θu(t+θ)dθ].This system has a single characteristic root at-1, and therefore the system is exponentially stable. To evaluate the stability of the feedback control internal dynamics
(11.41)u(t)=-2∫-r0‍e-θu(t+θ)dθ,
consider the characteristic equation
(11.42)1+2∫-r0‍e(s-1)θdθ=0,
or
(11.43)1+2(1-e-r(s-1))s-1=0.
It is easily seen that the system is stable for sufficiently small r>0. The smallest delay r for the above equation to have imaginary solutions can be obtained by letting s=jω in the above equation and solving the resulting equation. This r must be the positive solution of the equation
(11.44)2ercos(r4e2r-1)-1=0
such that
(11.45)s=i4e2r-1
is indeed a solution of (11.43). This can be obtained numerically as
(11.46)r=rmax=0.958399,
and the corresponding imaginary characteristic roots are
(11.47)s=±i5.118262.Consider now the discrete implementation of the control law (11.40) using the quasi-trapezoidal method (11.14). The coefficients can be easily calculated as
(11.48)C=-2e(K-1)r/K,Dk=-e(k-2)r/K(e2r/K-1),1≤k≤K-1,DK=-e(K-2)r/K(er/K-1).
Equation (11.37) or (11.38) may be used to obtain the maximum delay so that the difference equation (11.31) is exponentially stable. The values for various K are listed in  Table 1 as rmd.Table 1
K
2
3
5
10
15
100
1000
r
m
d
1.3863
1.0920
1.0000
0.9684
0.9628
0.9585
0.9584
r
m
a
0.921
0.920
0.909
0.926
0.935It can be seen from  Table1 that rmd indeed approaches rmax given in (11.46) as K increases.For the overall system, as the stability is not very sensitive to the retarded delayh^, we may analyze the system assuming h^=h. Theoretically the stability of such a system may be checked using such frequency domain methods as [95]. However, for a higher order system or with a large K, such analysis is rather unrealistic except some very special cases. On the other hand, a Lyapunov-Krasovskii functional method described in [52] may be used for this purpose. For the system in Example 11.2, using N=3 (according to the convention used in [52]), the maximum delays to retain exponential stability for various K up to 15 were obtained. The results are also listed in Table 1 as rma. It can be seen that as K increases, the maximum delay allowed for stability of the overall system approaches that for the difference equation. This is expected as for a larger K, the discrete delay closely approximates the distributed delay, and the only limiting factor is due to the stability of the difference equation.It should be pointed out that the theory covers only the case for largeK. No conclusion can be drawn regarding small K. Indeed, from Table 1, the cases for K=2 and 3 have larger maximum delays for stability than the case for K=5. This is not uncommon for other methods of approximation as well. For example, based on the theory presented in [129], Zhou et al. [130] used a sequence of distributed-delay operators that asymptotically approach the one obtained by finite spectrum assignment to stabilize marginally stable open-loop systems. A numerical example indicates that lower order approximations have larger stability margins.The case for commensurate delays is rather easy to implement in numerical simulation. One should use (11.30) and (11.31) instead of (11.1) and (11.14) in order to make sure all the components of y^ have the same delay. The simulation conducted in [15, 16] still slightly exceeds the limit shown in Table 1. Therefore, even if the simulation was done on (11.30) and (11.31), instability can still be expected, albeit it may take a rather long time to show up in step response.


### 11.4. Fundamental Limitation for Rationally Independent Delays

If the delays are subject to independent deviations, then no matter how small the deviations are, it is always possible to make the delays rationally independent. Therefore, in such a case, for practical stability, it is essential to consider the possibility of rationally independent delays. In this case, the stability of the system with discrete implementation is subject to a more strict fundamental limitation. Let(11.49)D(s)=∫-r0‍eθsG(θ)dθ.
Then, we may state the following result.Theorem 11.3.
If(11.50)supω=ℝρ[D(iω)]>1,
then the feedback control system consisting of (11.1) and a feedback control using either a rectangular-like implementation (11.6) or a trapezoidal-like implementation (11.14) with r1,r2,…,rK rationally independent is exponentially unstable for a sufficiently small hmax (defined in (11.9)).Proof.
From (11.50), there is a ω0 such that
(11.51)ρ[D(iω0)]>1.
As G is real, we may assume ω0>0 without loss of generality. Similar to the proof of Theorem 11.1,
(11.52)Dim(jω0)=∑k=1K‍Dkeiω0rk
can be made arbitrarily close to D(jω0) for a sufficiently small hmax. For such discrete implementation,
(11.53)ρ[∑k=1K‍Dkeiω0rk]>0.
According to the last part of Theorem 6.1, the difference equation is exponentially unstable. As a result, the whole system is exponentially unstable according to Theorem 5.1.For the special case of single-input single-output system with distributed delay determined by finite spectrum assignment, [126] reported that (11.50) may be replaced by D(0), and if the inequality is reversed, the discrete implementation is exponentially stable for sufficiently small hmax.The condition opposite to (11.50) is more strict than the stability of the internal dynamics of the distributed-delay feedback controller, as shown below.Proposition 11.4.
The roots of (11.18) are all on the strict left half plane; that is, (11.2) is exponentially stable, if
(11.54)supω=ℝρ[D(iω)]<1.Proof.
We will prove this by contradiction. Suppose(11.55)sup{Re(s)∣Δd(s)=0}≥0.
Consider the equation
(11.56)det(I-1zD(s))=0.
This defines solutions s as functions of z∈ℝ by implicit function theorem. For z=1, (11.55) indicates that there exists an s with Re(s)≥0 to satisfy (11.56). Consider this branch of solution as a continuous function of z, s=f(z). As z→+∞, it is obvious that Re(s)→-∞. By continuity, there must exists a z1≥1 such that s1=f(z1) satisfies Re(s1)=0. For this s1, since z1≥1 is an eigenvalue of D(s1),
(11.57)supω=ℝρ[D(s1)]≥1.
Therefore, (11.54) cannot be true.


### 11.5. Independent Delays

If each delay is individually implemented, then we can write(11.58)rk=krK+εk,
where εk are small errors. (Strictly speaking, if the system is multi-input multioutput, then we still require different components of yk(t) to have the same delay. If this cannot be enforced, then further complication may arise. However, the analysis may still be carried out using similar idea.) In order to formulate the problem to the standard form of (3.7), introduce the state variables
(11.59)y1(t)=u(t),yk(t)=u(t-rk-1),k=2,3,…,K,yK+1(t)=yK(t)=u(t-rK-1).
Then the system described by (11.1) and (11.14) may be written as
(11.60)x˙(t)=Ax(t)+ByK+1(t-hK+1),y1(t)=Cx(t)+∑k=1K‍Dkyk(t-hk),yk+1(t)=yk(t-hk),k=1,2,…,K-1,yK+1(t)=yK-1(t-hK-1),
where
(11.61)h1=r1=rK+ε1,hk=rk-rk-1=rK+εk-εk-1,k=2,3,…,K,hK+1=r-rk-1=rK-εK-1.
Instead of independent εk, k=1,2,…,K, we may alternatively consider
(11.62)δ1=ε1,δk=εk-εk-1,k=2,3,…,K,δK+1=-εK-1.
With this notation, the system may be written as
(11.63)x˙(t)=Ax(t)+B^(y^1(t-h1)y^2(t-h2)⋮y^K+1(t-hK+1)),(y^1(t)y^2(t)⋮y^K+1(t))=C^x(t)+D^(y^1(t-h1)y^2(t-h2)⋮y^K+1(t-hK+1)),
where the matrices B^, C^, and D^ are defined in (11.33). Again, as the last column block of the D^ vanishes, no conservatism is introduced by assuming εK+1 as independent, and Corollary 6.3 may be used to evaluate the stability of the difference equation described by D^ with partition structure mk=n, k=1,2,…,K+1. Actually, as the last column block of D^ vanishes, the dynamics of the difference equation is completely determined by the first K equations, and it is sufficient to apply Corollary 6.3 to the matrix D~ defined in (11.34) with the partition structure mk=n, k=1,2,…,K.For the overall system, a Lyapunov-Krasovskii functional method based on the description (11.60) has been formulated in [54]. Its numerical implementation is given in [58]. The formulation in [58] is directly applicable to this case.For the system given in Example11.2, since it is a scalar system, we may use (6.10) instead, which becomes
(11.64)∑k=1K-1‍e(k-2)r/K(e2r/K-1)+e(K-2)r/K(er/K-1)<1.
The above can be reduced to
(11.65)2er<1+2er/K.
The above condition shows that the exponential stability condition for the difference equation for all delays given in (11.58) in this example becomes more stringent as K increases. The maximum r that satisfies (11.65), that we denote as rmd, is listed in Table 2 for various K.Table 2
K
2
3
5
10
15
100
1000
r
m
d
0.6238
0.5256
0.4689
0.4346
0.4244
0.4081
0.4057The maximum delay for the overall system to be stable,rma, as estimated by the Lyapunov-Krasovskii functional method presented in [58] has also been calculated. It was found that the rma is identical to rmd for K up to 10. This is not surprising; only the spectrum associated with the difference equation is sensitive to small delay deviations. As shown in the commensurate delay case, the other characteristic roots would have allowed much larger rma.AsK→∞,  er/K→1, and  the condition (11.65) reduces to
(11.66)r<ln(32)≈0.4054.
As the characteristic root at -1 is continuous with respect to the delay, we can conclude the following about the discrete implementation about the system in Example 11.2.Corollary 11.5.
The system described by (11.39), (11.14), and (11.48) is practically stable for sufficiently large K if (11.66) is satisfied. On the other hand, if r>ln(3/2), then for any given ε>0, there exist a sufficiently large K and a set of εk that satisfy |εk|<ε, and the corresponding rk that satisfy (11.58) such that the system is exponentially unstable.To evaluate the system using the conditions given by Theorem11.3, we may calculate
(11.67)ρ[D(jω)]=|-2∫-r0‍e-θeiωθdθ|≤2∫-r0‍|e-θeiωθ|dθ=2∫-r0‍e-θdθ=2(er-1).
The upper bound above is actually tight as it is equal to ρ[D(jω)]ω=0. Therefore, Theorem 11.3 indicates that, for discrete feedback control with sufficiently small hmax, the complete system that is practically unstable is
(11.68)supω∈ℝρ[D(jω)]=2(er-1)>1,
or
(11.69)r>ln(32).
The analysis above indicates that this condition is tight at least for this system.It is interesting to note that the bound above in this particular scalar example was identical to the bound  obtained by Mirkin [10] based on w-stability.


### 11.6. An Intermediate Case

We may consider situations somewhere in between the two situations considered above. For example, if we want to implement the case forK=5, but we only have devices to implement 3-channeled simultaneous delays, we may consider using one device to implement the first three delays, and another to implement the remaining two delays. In this way, the closed-loop system becomes
(11.70)x˙(t)=Ax(t)+B^(y^1(t-h1)y^2(t-h2)y^3(t-h3)),(y^1(t)y^2(t)y^3(t))=C^x(t)+D^(y^1(t-h1)y^2(t-h2)y^3(t-h3)),
where y^1(t)∈ℝ3n, y^2(t)∈ℝ2n, y^3(t)=ℝn, and hk=r/5+εk, although ε3 is not independent, which is not important as it does not affect the difference equation. For the system in Example 11.2, using Corollary 6.3 yields the maximum delay for practical stability to be at 0.6563. The overall system using the method described in [58] again gives the same stability limit for the overall system, which is not surprising if we compare with the commensurate delay case.


## 11.1. Basic Formulation

Distributed-delay feedback control is an important method to control systems with delays. When the nominal system is unstable, the Smith predictor using the architecture shown in Figure1 involves unstable pole-zero cancellation, and an alternative implementation using distributed-delay feedback control may be used to avoid this problem [115]. Another common method to control systems with input or output delays is finite spectrum assignment [2, 47, 123, 124]. The system with single input-delay and distributed-delay feedback control may be expressed as
(11.1)x˙(t)=Ax(t)+Bu(t-r),(11.2)u(t)=Fx(t)+∫-r0‍G(θ)u(t+θ)dθ.
Exact implementation of control (11.2) is very difficult, if not impossible, due to the integration term. Some literature suggests using discrete delays to approximate the distributed delay in (11.2) [2]. The basic idea is to divide the interval [-r,0] to K small intervals, [θk,θk-1], k=1,2,…,K, with
(11.3)0=θ0>θ1>θ2>⋯>θK=-r,
and to approximate the integration in each small interval [θk,θk-1]. We may choose a point in each interval
(11.4)-rk∈[θk,θk-1].
Common choices are
(11.5)-rk=θk,θk-1,orθk+θk-12.
When all the intervals are sufficiently small, it seems reasonable to approximate any u(t+θ) with θ∈[θk,θk-1] by u(t-rk). Therefore, (11.2) may be approximated by
(11.6)u(t)=Fx(t)+∑k=1K‍Gku(t-rk),
where
(11.7)Gk=∫θkθk-1‍G(θ)dθ.
Less accurate method of obtaining Gk may be used. For example, it is common to use
(11.8)Gk=G(-rk)(θk-1-θk),
leading to a rectangular rule of numerical integration. As was discussed by Zhong [18], (11.7) compares favorably with (11.8) in terms of accuracy, although the difference is not fundamental. We will call (11.6) quasi-rectangular implementation, and any choice of Gk and rk along this line will be known as a rectangular-like implementation. If G(θ)u(t+θ) is Riemann integrable, then for any rectangular-like implementation, the expression on the right hand side of (11.6) converges to that of (11.2) as hmax→0, where
(11.9)hmax=max1≤k≤K(θk-1-θk).Alternatively, we may use[u(t+θk)+u(t+θk-1)]/2 to approximate u(t+θ), θ∈[θk,θk-1]. As a result, (11.6) should be replaced by
(11.10)u(t)=Fx(t)+∑k=0K‍Hku(t-rk),
where
(11.11)H0=12∫θ10‍G(θ)dθ,Hk=12∫θk+1θk-1‍G(θ)dθ,1≤k≤K-1,HK=12∫θKθK-1‍G(θ)dθ,(11.12)rk=-θk,k=0,1,…,K.
Similarly, a less accurate choice is to use
(11.13)H0=-12(θ0-θ1)G(0),Hk=12(θk-1-θk+1)G(θk),1≤k≤K-1,HK=12(θK-1-θK)G(θK),
which becomes the trapezoidal rule of numerical integration. Again, there is no fundamental difference between the more accurate expression (11.11) and the slightly less accurate expression (11.13). However, in this case, the feedback rule (11.10) is not well posed as the nonzero coefficient H0 at r0=0 violates the uniform nonatomic requirement discussed in Section 2. This problem can be avoided if the size of the interval [θ1,0] is sufficiently small, in which case H0 is small, and I-H0 is invertible. We may solve for u(t) in (11.10) to obtain a well-posed implementation rule,
(11.14)u(t)=Cx(t)+∑k=1K‍Dku(t-rk),
where
(11.15)C=(I-H0)-1F,Dk=(I-H0)-1Hk,k=1,2,…,K,rk=-θk.
We will call (11.14) quasi-trapezoidal implementation. In practice, there may also be small deviation of rk. Any of these implementations will be known as a trapezoidal-like implementation. Obviously, if G(θ)u(t+θ) is Riemann integrable, then for any trapezoidal-like implementation, the expression on the right hand side of (11.10) or (11.14) converges to that of (11.2) as hmax→0.For a uniform gridding(11.16)θk=-krK,
it was pointed out in [15, 16] that instability was observed even in numerical simulation when control rules similar to (11.6) or (11.14) are used in place of (11.2) even when hmax becomes very small. The problem was declared by Richard in the survey paper [17] as one of the open problems in the control of time-delay systems. Significant insight has been gained since then [10, 18, 125, 126], and alternative implementation strategies have been proposed [10, 126–128]. This section will only discuss the stability property of such systems and will not discuss improvement of implementation. The material presented here is more general and many results are often more precise.


## 11.2. Fundamental Limitation for General Case

Consider the internal dynamics of the distributed-delay feedback (11.2),
(11.17)u(t)=∫-r0‍G(θ)u(t+θ)dθ.
It is exponentially stable if and only if all the roots of its characteristic equation
(11.18)Δd(s)=det(I-∫-r0‍eθsG(θ)dθ)=0
are on the strict left half plane. A fundamental limitation of discrete implementation is due to this fact.Theorem 11.1.
Let the matrix functionG(θ) be Riemann integrable in [-r,0]. If the internal dynamics of (11.17) is exponentially unstable, then the feedback control system consisting of (11.1) and a feedback control using either a rectangular-like implementation (11.6) or a trapezoidal-like implementation (11.14) is exponentially unstable for a sufficiently small hmax (defined in (11.9)).Proof.
If (11.17) is exponentially unstable, then there exists a s0, Re(s0)>0, such that (11.18) is satisfied for s=s0. The function Δd(s) is an entire function that is not a constant. Therefore, the order of root s0 is finite, and there is an a>0 such that Δd(s)≠0 for any s in the region 0<|s-s0|≤a. We may obviously make
(11.19)a<Re(s0).
Let
(11.20)Γ={s∈ℂ∣|s-s0|=a}.
Since Γ is compact, we may define
(11.21)b=min{|Δd(s)|∣s∈Γ}.
Obviously
(11.22)b>0.
In the compact set Γ, as either rectangular-like implementation or trapezoidal-like implementation can approximate the integration in (11.18) to arbitrary accuracy with sufficiently small hmax, the related quantity Δd(s) may also be approximated to arbitrary accuracy. To be definite, consider the quasi-trapezoidal implementation. There exists a h^, such that for any hmax≤h^,
(11.23)|Δim(s)-Δd(s)|<b2,fors∈Γ,
where
(11.24)Δim(s)=det(I-∑k=1K‍Dkerks).
However, according to Rouché’s theorem [75], the inequality (11.23) implies that Δim(s) has the same number of roots as Δd(s) within the region enclosed by Γ, which is completely on the right half complex plane due to (11.19). Therefore, Δim(s) has at least one right half plane root. As
(11.25)Δim(s)=0
is the characteristic equation of the difference equation
(11.26)u(t)=∑k=1K‍Dku(t-rk),
this difference equation is exponentially unstable. However, according to Theorem 5.1, this also means that the complete system consisting of (11.1) and (11.14) is exponentially unstable.The special case of the above theorem with uniform gridding and finite spectrum assignment can be found in [125]. The above theorem gives a fundamental limitation of any discrete implementation of a distributed-delay feedback control system. We may understand the situation by reversing the role of the plant dynamics (11.1) and the controller dynamics (11.2) or its discrete implementation, such as (11.14), and considering x as the input to the “plant” (11.2) or its  discretized version (11.14). In the case of (11.2), as its solution operator is compact, it is possible for an appropriate (11.1) to make the whole system stable. On the other hand, as the solution operator for (11.14) is not compact, no linear differential equation in the form of (11.1) may stabilize it as its unstable essential spectrum cannot be changed by a compact operator.


## 11.3. Commensurate Delays

The most favorable implementation is using commensurate delays. Imagine if we have a multichanneled device that delays the inputs in all channels byh. Like any device, there is bound to be errors. However, the delays for all channels will be identical. If such a device is available, then we can imagine implementing the feedback law as shown in Figure 2. In this case, when error is considered, the discrete delays become
(11.27)rk=kh,(11.28)h=rK+εh,|εh|<ε.
The proportional relation (11.27) is guaranteed by the structure of the system and does not contain any error. Although the delays still contain errors, they are not independent.Figure 2
Implementation to enforce proportional delays.Introduce the variable(11.29)y^1(t)=(y1(t)y2(t)⋮yK(t))=(u(t)u(t-h)u(t-2h)⋮u(t-(K-1)h)),y^2(t)=yK(t).
Then, (11.1) and (11.14) become
(11.30)x˙(t)=Ax(t)+B^(y^1(t-h)y^2(t-h^)),(11.31)(y^1(t)y^2(t))=C^x(t)+D^(y^1(t-h)y^2(t-h^)),
where
(11.32)h^=r-(K-1)h=rK-(K-1)εh,(11.33)B^=(00⋯0B),C^=(C0⋮0),D^=(D1D2⋯DK-1DK0I0⋯0000I⋯000⋮⋮⋱⋮⋮⋮00⋯I0000⋯I00).
Let
(11.34)D~=(D1D2⋯DK-1DKI0⋯000I⋯00⋮⋮⋱⋮⋮00⋯I0).
Then
(11.35)D^=(D~0En0),
where
(11.36)En=(00⋯0I0).
Although the delays h and h^ are not exactly independent, considering them as independent does not increase the conservatism as the last column block of D^ vanishes. Indeed, under this assumption, the condition for exponential stability of the difference equation (11.31) for all permissible small delay errors becomes
(11.37)maxδk∈ℂ|δk|=1ρ(D^E(δ1,δ2))=maxδ1∈ℂ|δ1|=1ρ(D^E(δ1,0))=ρ(D~)<1,
or(11.38)max{|λ||det(λKI-∑k=1K‍λK-kDk)=0}<1.
It is noticed that the above condition is independent of the delay error εh.To get a sense of how well this implementation works, consider the special case of the example discussed in [10, 15, 16], which is restated below.Example 11.2.
The system with scalar variables(11.39)x˙(t)=x(t)+u(t-r)
is stabilized using the following distributed-delay feedback control using finite spectrum assignment method:
(11.40)u(t)=-2[erx(t)+∫-r0‍e-θu(t+θ)dθ].This system has a single characteristic root at-1, and therefore the system is exponentially stable. To evaluate the stability of the feedback control internal dynamics
(11.41)u(t)=-2∫-r0‍e-θu(t+θ)dθ,
consider the characteristic equation
(11.42)1+2∫-r0‍e(s-1)θdθ=0,
or
(11.43)1+2(1-e-r(s-1))s-1=0.
It is easily seen that the system is stable for sufficiently small r>0. The smallest delay r for the above equation to have imaginary solutions can be obtained by letting s=jω in the above equation and solving the resulting equation. This r must be the positive solution of the equation
(11.44)2ercos(r4e2r-1)-1=0
such that
(11.45)s=i4e2r-1
is indeed a solution of (11.43). This can be obtained numerically as
(11.46)r=rmax=0.958399,
and the corresponding imaginary characteristic roots are
(11.47)s=±i5.118262.Consider now the discrete implementation of the control law (11.40) using the quasi-trapezoidal method (11.14). The coefficients can be easily calculated as
(11.48)C=-2e(K-1)r/K,Dk=-e(k-2)r/K(e2r/K-1),1≤k≤K-1,DK=-e(K-2)r/K(er/K-1).
Equation (11.37) or (11.38) may be used to obtain the maximum delay so that the difference equation (11.31) is exponentially stable. The values for various K are listed in  Table 1 as rmd.Table 1
K
2
3
5
10
15
100
1000
r
m
d
1.3863
1.0920
1.0000
0.9684
0.9628
0.9585
0.9584
r
m
a
0.921
0.920
0.909
0.926
0.935It can be seen from  Table1 that rmd indeed approaches rmax given in (11.46) as K increases.For the overall system, as the stability is not very sensitive to the retarded delayh^, we may analyze the system assuming h^=h. Theoretically the stability of such a system may be checked using such frequency domain methods as [95]. However, for a higher order system or with a large K, such analysis is rather unrealistic except some very special cases. On the other hand, a Lyapunov-Krasovskii functional method described in [52] may be used for this purpose. For the system in Example 11.2, using N=3 (according to the convention used in [52]), the maximum delays to retain exponential stability for various K up to 15 were obtained. The results are also listed in Table 1 as rma. It can be seen that as K increases, the maximum delay allowed for stability of the overall system approaches that for the difference equation. This is expected as for a larger K, the discrete delay closely approximates the distributed delay, and the only limiting factor is due to the stability of the difference equation.It should be pointed out that the theory covers only the case for largeK. No conclusion can be drawn regarding small K. Indeed, from Table 1, the cases for K=2 and 3 have larger maximum delays for stability than the case for K=5. This is not uncommon for other methods of approximation as well. For example, based on the theory presented in [129], Zhou et al. [130] used a sequence of distributed-delay operators that asymptotically approach the one obtained by finite spectrum assignment to stabilize marginally stable open-loop systems. A numerical example indicates that lower order approximations have larger stability margins.The case for commensurate delays is rather easy to implement in numerical simulation. One should use (11.30) and (11.31) instead of (11.1) and (11.14) in order to make sure all the components of y^ have the same delay. The simulation conducted in [15, 16] still slightly exceeds the limit shown in Table 1. Therefore, even if the simulation was done on (11.30) and (11.31), instability can still be expected, albeit it may take a rather long time to show up in step response.


## 11.4. Fundamental Limitation for Rationally Independent Delays

If the delays are subject to independent deviations, then no matter how small the deviations are, it is always possible to make the delays rationally independent. Therefore, in such a case, for practical stability, it is essential to consider the possibility of rationally independent delays. In this case, the stability of the system with discrete implementation is subject to a more strict fundamental limitation. Let(11.49)D(s)=∫-r0‍eθsG(θ)dθ.
Then, we may state the following result.Theorem 11.3.
If(11.50)supω=ℝρ[D(iω)]>1,
then the feedback control system consisting of (11.1) and a feedback control using either a rectangular-like implementation (11.6) or a trapezoidal-like implementation (11.14) with r1,r2,…,rK rationally independent is exponentially unstable for a sufficiently small hmax (defined in (11.9)).Proof.
From (11.50), there is a ω0 such that
(11.51)ρ[D(iω0)]>1.
As G is real, we may assume ω0>0 without loss of generality. Similar to the proof of Theorem 11.1,
(11.52)Dim(jω0)=∑k=1K‍Dkeiω0rk
can be made arbitrarily close to D(jω0) for a sufficiently small hmax. For such discrete implementation,
(11.53)ρ[∑k=1K‍Dkeiω0rk]>0.
According to the last part of Theorem 6.1, the difference equation is exponentially unstable. As a result, the whole system is exponentially unstable according to Theorem 5.1.For the special case of single-input single-output system with distributed delay determined by finite spectrum assignment, [126] reported that (11.50) may be replaced by D(0), and if the inequality is reversed, the discrete implementation is exponentially stable for sufficiently small hmax.The condition opposite to (11.50) is more strict than the stability of the internal dynamics of the distributed-delay feedback controller, as shown below.Proposition 11.4.
The roots of (11.18) are all on the strict left half plane; that is, (11.2) is exponentially stable, if
(11.54)supω=ℝρ[D(iω)]<1.Proof.
We will prove this by contradiction. Suppose(11.55)sup{Re(s)∣Δd(s)=0}≥0.
Consider the equation
(11.56)det(I-1zD(s))=0.
This defines solutions s as functions of z∈ℝ by implicit function theorem. For z=1, (11.55) indicates that there exists an s with Re(s)≥0 to satisfy (11.56). Consider this branch of solution as a continuous function of z, s=f(z). As z→+∞, it is obvious that Re(s)→-∞. By continuity, there must exists a z1≥1 such that s1=f(z1) satisfies Re(s1)=0. For this s1, since z1≥1 is an eigenvalue of D(s1),
(11.57)supω=ℝρ[D(s1)]≥1.
Therefore, (11.54) cannot be true.


## 11.5. Independent Delays

If each delay is individually implemented, then we can write(11.58)rk=krK+εk,
where εk are small errors. (Strictly speaking, if the system is multi-input multioutput, then we still require different components of yk(t) to have the same delay. If this cannot be enforced, then further complication may arise. However, the analysis may still be carried out using similar idea.) In order to formulate the problem to the standard form of (3.7), introduce the state variables
(11.59)y1(t)=u(t),yk(t)=u(t-rk-1),k=2,3,…,K,yK+1(t)=yK(t)=u(t-rK-1).
Then the system described by (11.1) and (11.14) may be written as
(11.60)x˙(t)=Ax(t)+ByK+1(t-hK+1),y1(t)=Cx(t)+∑k=1K‍Dkyk(t-hk),yk+1(t)=yk(t-hk),k=1,2,…,K-1,yK+1(t)=yK-1(t-hK-1),
where
(11.61)h1=r1=rK+ε1,hk=rk-rk-1=rK+εk-εk-1,k=2,3,…,K,hK+1=r-rk-1=rK-εK-1.
Instead of independent εk, k=1,2,…,K, we may alternatively consider
(11.62)δ1=ε1,δk=εk-εk-1,k=2,3,…,K,δK+1=-εK-1.
With this notation, the system may be written as
(11.63)x˙(t)=Ax(t)+B^(y^1(t-h1)y^2(t-h2)⋮y^K+1(t-hK+1)),(y^1(t)y^2(t)⋮y^K+1(t))=C^x(t)+D^(y^1(t-h1)y^2(t-h2)⋮y^K+1(t-hK+1)),
where the matrices B^, C^, and D^ are defined in (11.33). Again, as the last column block of the D^ vanishes, no conservatism is introduced by assuming εK+1 as independent, and Corollary 6.3 may be used to evaluate the stability of the difference equation described by D^ with partition structure mk=n, k=1,2,…,K+1. Actually, as the last column block of D^ vanishes, the dynamics of the difference equation is completely determined by the first K equations, and it is sufficient to apply Corollary 6.3 to the matrix D~ defined in (11.34) with the partition structure mk=n, k=1,2,…,K.For the overall system, a Lyapunov-Krasovskii functional method based on the description (11.60) has been formulated in [54]. Its numerical implementation is given in [58]. The formulation in [58] is directly applicable to this case.For the system given in Example11.2, since it is a scalar system, we may use (6.10) instead, which becomes
(11.64)∑k=1K-1‍e(k-2)r/K(e2r/K-1)+e(K-2)r/K(er/K-1)<1.
The above can be reduced to
(11.65)2er<1+2er/K.
The above condition shows that the exponential stability condition for the difference equation for all delays given in (11.58) in this example becomes more stringent as K increases. The maximum r that satisfies (11.65), that we denote as rmd, is listed in Table 2 for various K.Table 2
K
2
3
5
10
15
100
1000
r
m
d
0.6238
0.5256
0.4689
0.4346
0.4244
0.4081
0.4057The maximum delay for the overall system to be stable,rma, as estimated by the Lyapunov-Krasovskii functional method presented in [58] has also been calculated. It was found that the rma is identical to rmd for K up to 10. This is not surprising; only the spectrum associated with the difference equation is sensitive to small delay deviations. As shown in the commensurate delay case, the other characteristic roots would have allowed much larger rma.AsK→∞,  er/K→1, and  the condition (11.65) reduces to
(11.66)r<ln(32)≈0.4054.
As the characteristic root at -1 is continuous with respect to the delay, we can conclude the following about the discrete implementation about the system in Example 11.2.Corollary 11.5.
The system described by (11.39), (11.14), and (11.48) is practically stable for sufficiently large K if (11.66) is satisfied. On the other hand, if r>ln(3/2), then for any given ε>0, there exist a sufficiently large K and a set of εk that satisfy |εk|<ε, and the corresponding rk that satisfy (11.58) such that the system is exponentially unstable.To evaluate the system using the conditions given by Theorem11.3, we may calculate
(11.67)ρ[D(jω)]=|-2∫-r0‍e-θeiωθdθ|≤2∫-r0‍|e-θeiωθ|dθ=2∫-r0‍e-θdθ=2(er-1).
The upper bound above is actually tight as it is equal to ρ[D(jω)]ω=0. Therefore, Theorem 11.3 indicates that, for discrete feedback control with sufficiently small hmax, the complete system that is practically unstable is
(11.68)supω∈ℝρ[D(jω)]=2(er-1)>1,
or
(11.69)r>ln(32).
The analysis above indicates that this condition is tight at least for this system.It is interesting to note that the bound above in this particular scalar example was identical to the bound  obtained by Mirkin [10] based on w-stability.


## 11.6. An Intermediate Case

We may consider situations somewhere in between the two situations considered above. For example, if we want to implement the case forK=5, but we only have devices to implement 3-channeled simultaneous delays, we may consider using one device to implement the first three delays, and another to implement the remaining two delays. In this way, the closed-loop system becomes
(11.70)x˙(t)=Ax(t)+B^(y^1(t-h1)y^2(t-h2)y^3(t-h3)),(y^1(t)y^2(t)y^3(t))=C^x(t)+D^(y^1(t-h1)y^2(t-h2)y^3(t-h3)),
where y^1(t)∈ℝ3n, y^2(t)∈ℝ2n, y^3(t)=ℝn, and hk=r/5+εk, although ε3 is not independent, which is not important as it does not affect the difference equation. For the system in Example 11.2, using Corollary 6.3 yields the maximum delay for practical stability to be at 0.6563. The overall system using the method described in [58] again gives the same stability limit for the overall system, which is not surprising if we compare with the commensurate delay case.


## 12. Marginal Case

Let the characteristic equation of a time-delay system be(12.1)Δ(s)=0.
Let
(12.2)σ-=sup{Re(s)∣Δ(s)=0}.
Recall from Section 5 that the system is exponentially stable if σ-<0, and it is exponentially unstable if σ->0. The situation for
(12.3)σ-=0
is more complicated. As such cases are encountered in some practical control methods [131–133], it is of interest to provide an overview of existing results.A few cases parallel to the systems without delay are rather obvious. If there are a finite number of characteristic roots that are simple and on the imaginary axis, and the remaining roots satisfy(12.4)Re(s)<-ε
for some ε>0, then the system is marginally stable, and any solution with a bounded initial condition remains bounded. A multiple root on the imaginary axis may cause instability.Complication arises when all the characteristic rootss satisfy
(12.5)Re(s)<0,
but there are a series of characteristic roots sk that asymptotically approach the imaginary axis
(12.6)Re(sk)→0.
This is possible only if
(12.7)σ-0=0,
where
(12.8)σ-0=sup{Re(s)∣Δ0(s)=0},
add Δ0(s)=0 is the characteristic equation of the associated difference equation. For a single delay system, it was shown by Hahn [134] that such a system is always asymptotically stable but not exponentially stable. For example, the following simple example given by Datko [71]
(12.9)x˙(t)-x˙(t-1)+x(t)=0
can be shown as asymptotically stable. Its characteristic equation is
(12.10)g(s)=s(1-e-s)+1=0.
However, the following system given in [74]
(12.11)d2dt2[x(t)-2x(t-1)+x(t-2)]+2ddt[x(t)-x(t-1)]+x(t)=0
is unstable. Its characteristic equation is
(12.12)[g(s)]2=[s(1-e-s)+1]2=0,
which obviously has identical distribution of characteristic roots as that for the system (12.9) except that all the roots are double roots instead. The reason is that there is a possibility of an extra t for the solutions corresponding to each characteristic root of (12.11) as compared to the solutions for (12.9).Snow [135] seems to be the first attempt to construct an unstable system that satisfies (12.5) and (12.6). Although it contributed some major ideas for the other works, it contains a mistake. A correct construction of such system with some generality is given in Brumley [73]. The constructed system has commensurate delays. The main idea is to estimate the rate of convergence for sk to the imaginary axis, and the order of multiple characteristic roots needed to drive the system unstable. Gromova [72] also provides a method of constructing such systems and also commented on the case for incommensurate delays.


## 13. Conclusions

The stability analysis of time-delay systems of neutral type requires understanding of some subtle points, especially the discontinuity of the spectrum of the associated difference equation. Especially, care must be taken when continuity of characteristic roots is used. A review of some major points are provided that integrates the coupled differential-difference equation formulation. Some pratical problems are discussed based on the theories. These include the small delay problem, the sensitivity of Smith predictor, and the discrete implementation of distributed-delay feedback control. Some derivations are simplified, some results are strengthened or extended to more general case, and new perspectives are gained.



---

*Source: 725783-2012-11-27.xml*