Stability, Control and Application of Time-Delay Systems

Stability, Control and Application of Time-Delay Systems gives a systematic description of these systems. It includes adequate designs of integrated modeling and control and frequency characterizations. Common themes revolve around creating certain synergies of modeling, analysis, control, computing and applications of time delay systems that achieve robust stability while retaining desired performance quality. The book provides innovative insights into the state-of-the-art of time-delay systems in both theory and practical aspects. It has been edited with an emphasis on presenting constructive theoretical and practical methodological approaches and techniques.

• Unifies existing and emerging concepts concerning time delay dynamical systems
• Provides a series of the latest results in large-delay analysis and multi-agent and thermal systems with delays
• Gives in each chapter numerical and simulation results in order to reflect the engineering practice

• Language: English
• Publisher: Elsevier Science
• Release date: Jun 27, 2019
• ISBN: 9780128149294

Book preview

China

Chapter 1

On the numerical determination of stability regions in the delay space via dominant pole estimation

Libor Pekař    Faculty of Applied Informatics, Department of Automation and Control Engineering, Tomas Bata University in Zlín, Zlín, Czechia

Abstract

This contribution intends to provide the reader with the presentation of a numerical gridding iterative algorithm to determine all stability regions within the prescribed area of the delay space. In every single grid node, the iterative estimation of the rightmost pole is computed based on the polynomial approximation of the characteristic quasipolynomial, by utilizing the knowledge of the dominant pole estimation in the nearest grid node. The polynomial approximation is made via the Taylor series-based expansion in the vicinity of the closest dominant poles estimation, and by using the bilinear transformation followed with prewarping for a discrete-time approximation. Exponential terms are subjected to a quadratic extrapolation method to get commensurate delays. Two-step Newton’s iteration method with averaging is used to detect imaginary axis crossings. Neutral delay case is concisely discussed as well. Two numerical examples demonstrate the accuracy and efficiency of the algorithm. Possible future directions of this research and algorithm modifications are proposed and discussed in brief as well.

Keywords

Bilinear transformation; Delay-dependent stability; Exponential stability; Newton’s iteration method; Quadratic extrapolation; Taylor series expansion

Chapter outline

1Introduction

2Preliminaries

2.1LTI-TDS model

2.2Stability and spectral properties

3Numerical gridding DDS algorithm

3.1Framework of the algorithm

3.2Algorithm steps in detail

3.3Continuous-time approximation

3.4Discrete-time approximation

3.5Approximation of neutral quasipolynomials

4Examples

4.1Retarded system

4.2Neutral system

5Conclusions

Acknowledgments

References

Acknowledgments

We would like to express our gratitude to the European Regional Development Fund and to the Ministry of Education, Youth and Sports that financially supported this work under the National Sustainability Program project No. LO1303 (MSMT-7778/2014).

1 Introduction

It is well known that linear time-invariant time-delay systems (LTI-TDSs) are typical representatives of infinite-dimensional systems [1]. It means that they own infinitely many system poles (or, characteristic values), the loci of which have many interesting and tricky features [2, 3]. This fact inherently implies that analytic and control tasks related to LTI-TDSs and their spectra are generally nontrivial. Therefore, some rigorous solutions suffer from rich or hardly implementable mathematical formulations.

Characteristic root loci of LTI-TDSs are closely related to system stability [or bounded-input bounded-output stable even if there exist purely imaginary poles, or the rightmost spectrum asymptotically reach the imaginary axis [2, 5, 6]. Hence, the detection of imaginary axis crossing can be used as a tool for exponential stability analysis.

The systematic searching of stabilizing regions in the delay space is called the delay-dependent stability (DDS) analysis. A couple of research results have dealt with the study of DDS. A family of indirect DDS methods utilizes Lyapunov-Krasovskii approaches and linear matrix inequalities; however, these approaches usually suffer from mathematical complexity, provide conservative results, and they are difficult to be implemented in practice—see, for instance, a discussion in Ref. [7]. Contrariwise, direct (or, frequency-domain) approaches adopt the above-introduced idea of the determination of purely imaginary poles that constitute the stability margin. Several techniques within this framework have been investigated so far. Namely, the cluster treatment of characteristic roots (CTCR) paradigm [8, 9], the direct method [10], the Puiseux series expansion technique [11], the matrix pencil method [12], the Kronecker multiplication method [13], the use of the argument principle [14], etc. Despite numerous DDS ideas and principles, there is still a lack of relatively simple, easily programmable, and practically well-implementable DDS methods that give a sufficiently accurate estimation of stabilizing regions in the delays space.

Hence, this chapter is aimed at the presentation of a gridding-based direct computational method for the determination of all multiple delays for LTI-TDSs ensuring the system (exponential) stability. The core of the proposed method lies in the iterative estimation of the rightmost (dominant) pole or a conjugate pair in every single-grid node of the discretized delay space. The delay-free dominant pole can initially be simply and exactly computed as a polynomial root. Since system pole positions are continuous with respect to delay values almost everywhere, the approximated already known dominant pole locus in the nearest grid node to the current one serves as the initial estimation for the iterative computation. The iterative estimations are computed via a polynomial approximation of the so-called characteristic quasipolynomial. Two ideas how to find the approximating polynomial are suggested. As first, the Taylor series-based expansion of a particular order in the neighborhood of the nearest dominant poles estimation is performed [15, 16]. As second, the Tustin (bilinear) transformation followed by the prewarping correction is utilized to get the discrete form of the approximating characteristic polynomial [17, 18]. Here, the delay (exponential) terms are subjected to a quadratic extrapolation method to get commensurate delays that are closely related to the discrete-time shifting operator. Once an imaginary axis crossing is detected, the so-called root tendency values are used to get a more accurate critical delay and the corresponding critical frequency estimations, by means of averaged (multistep) Newton’s zero-point extrapolation principle. Finally, the set of found critical delays is eventually augmented via quadratic regression of the found points.

This chapter, in fact, summarizes, integrates, and extends the two above-referred approaches [16, 17]. Both the ideas are compared, and an extension to a very delicate family of so-called neutral systems is suggested—strong stability, the essential spectrum, and the approximation of the associated characteristic exponential polynomial are discussed as well. An example solving the stabilization of a skater on the controlled swaying bow, modeled as an unstable retarded TDS, is given to the reader. The obtained numerical results are compared to those received by the CTCR paradigm. Another academic example suggests a possibility how to modify the algorithm in the case of neutral delays. The Quasi-Polynomial mapping Rootfinder (QPmR) software tool by Vyhlídal and Zítek [19] is used to get the reference values in the examples.

The rest of the chapter is organized as follows: in Section 2, the definition of an LTI-TDS along with the introduction of its exponential stability and selected spectral properties are given. Then, the reader is acquainted with an overview of the proposed DDS algorithm including its continuous- and discrete-time versions in Section 3. Moreover, the approximation of neutral quasipolynomials and its associated difference equation is suggested. Section 4 includes two numerical examples. The chapter is then concluded with a discussion of possible algorithm modifications and its summary.

Notation

Re(s) and Im(sis the n-dimensional Euclidean space. F(s | p) stands for a function of s with a parameter set p. The superscript T means the vector or matrix transpose, and ⌊⋅⌋ represents the floor function. Symbol ∥⋅∥s means the supreme norm.

2 Preliminaries

In this preliminary section, the class of LTI-TDSs for the purpose of this contribution is introduced first; then, exponential stability is formulated and some important properties of the spectrum of LTI-TDS characteristic values are presented.

2.1 LTI-TDS model

Definition 1

Let an LTI-TDS be formulated by state and output functional differential equations as

   (1)

are delays, and A0, Ai, B0, Bi, C, Ci, Hi express real-valued matrices of compatible dimensions.

Note that, in Definition 1, only lumped delays are considered for the simplicity. Whenever ∃i such that Hi≠0, a system is called neutral (NTDS); otherwise, it is retarded (RTDS).

The characteristic function of system (1) can be constructed using the Laplace transform as

   (2)

di(s. The so-called associated characteristic exponential polynomial , where dn, i .

2.2 Stability and spectral properties

Definition 2

Gu et al. [4], Michiels and Niculescu [3], Zhang and Sun [7]

System (1) is said to be exponentially stable if there exist a > 0, μ , ∀t ≥ 0 for all φ, where φ(θ) = x(θ) with θ ∈ [−L, 0] stands for the initial condition continuous-time function.

In order to characterize exponential stability by means of the system spectrum Σ := {s : Δ(s) = 0}, let us define the spectral abscissa first. Then, the following statement holds.

Proposition 1

Michiels and Niculescu [3]

An RTDS is exponentially stable if and only if α(⋅) < 0. An NTDS is exponentially stable if and only if there exists ε > 0 such that α(⋅) < −ε.

(or a pair) appears exactly on the imaginary axis, one can determine the switching frequency and the corresponding switching delay identifies a point of the stability border in the delay space.

In contrast to delay-free LTI systems, the (exponential) stability analysis and/or synthesis is, however, even more complicated for the delayed case. This is mainly due to the spectral properties of NTDSs. Let γ be the supreme of real parts of the so-called essential spectrum Σess := {s : dn(s) = 0}. Moreover, let us define

. Now, an overview of some LTI-TDS spectral properties follows.

Proposition 2

Hale and Lunel [1], Michiels and Niculescu [3], Vanbiervliet et al. [20]

The following statements hold:

(1)There are only finitely many characteristic roots sk ∈ Σ in the half-plane Re(s) > β , and these poles are isolated, for an RTDS. On the contrary, subsets of system and essential poles (sk, ess ∈ Σess) constitute vertical strips at high frequencies that asymptotically converge to each other, for an NTDS, and there may be located infinitely many system poles in the half-plane Re(s) > γ.

(2)For an RTDS, system poles behave continuously and smoothly with respect to system parameters and τ, represents the vector of all system delays. Whereas, the value of γ is not continuous with respect to delays for an NTDS. Moreover, function α(τ) may be nonsmooth and/or non-Lipschitz at some points.

It implies from : in such a case, a system with α < 0 cannot be fully assumed as the stable one. In addition, the shape of the rightmost infinite vertical strip of system poles has to be taken into account as well. A possibility of abrupt changes in γ gives rise to the notion of strong stability.

Definition 3

Michiels and Niculescu [3], Michiels and Vyhlídal [21]

System (1) is said to be strongly stable .

Proposition 3

System (1) is strongly stable if and only if

   (3)

It is worth noting that Michiels and Vyhlídal [21] proposed an effective estimation of the safe upper bound c , which is, inter alia, a continuous function of τ.

3 Numerical gridding DDS algorithm

Two versions of the proposed DDS determination procedure are presented in this section by means of a concise formulation via a summarizing step-by-step framework algorithm. Both the versions, continuous- and discrete-time one, have already been published in Refs. [16, 17], respectively, in detail. Herein, a unifying and clear concept extending the algorithm steps is given to the reader, followed by an idea how to cope with NTDSs.

3.1 Framework of the algorithm

The framework algorithm can be formulated in Algorithm 1.

Algorithm 1

Numerical gridding DDS algorithm

(1)The delay space is equidistantly discretized with a selected discretization step.

(2)Dominant pole(s) for the delay-free case is (are) exactly computed.

(3)Steps 4–6 are performed for every single node of the discretized delay space (one goes successively in parallel to particular delay axes).

(4)The dominant pole estimation in the closest node is set.

(5)In the neighborhood of this estimation, an updated dominant pole estimation is computed (for the particular node). If the imaginary axis is not crossed, the loop is finished (see Step 3); otherwise, go to Step 6.

(6)Approximations of switching delays and switching frequencies are computed by means of known (current and preceding) dominant poles and corresponding delays’ estimations. These values are computed successively, that is, delay by delay in the delay space.

(7)The algorithm output is represented by the sets of switching delays and switching frequencies after an additional interpolation of points obtained in Step 6.

3.2 Algorithm steps in detail

Now, particular steps of Algorithm 1 are explained in more detail excluding Step 5 which is introduced separately in Sections 3.3 and 3.4.

Regarding Step 1, let the mesh grid in the delay space (for the given Δ(s, τ)) be

with a discretization step Δτ⋅, j = τ⋅, j+1 − τ⋅, j) are initialized as the empty ones.

Step 2 can simply be solved via standard software tools as the following polynomial root-searching problem:

   (4)

Let this delay-free dominant pole obtained from Eq. ((the value is saved). The nested loops according to Step 3 are performed such that the outer loop goes through delays while the inner ones go through grid nodes for the particular delay indexed by jl, l = 1, 2, …, nτcan be found at the position

   (5)

. Hence, the corresponding preceding dominant pole estimation according to Step 4 reads

   (6)

where the subscript denotes the grid position.

Regarding Step 6, it is worth a priori noting that a pair of purely imaginary characteristic roots crosses the imaginary axis only if the speed of this pair in the real axis with respect to delays is nonzero. This speed can be quantified by the value of the root tendency function defined as

   (7)

.

The value of RT can be used to compute the switching delays by applying Newton’s zero point extrapolation principle as follows. Consider a fixed pole sk, the locus of which depends on τ, . Then, the elements of τ⁰ can be approximated as

   (8)

When the crossing is detected in Step 5 of Algorithm 1, we have preceding (sp) and current (sc) dominant pole values (and corresponding grid nodes in the delay space, τp, τc, is computed repeatedly via Eq. (8) with iteratively updated estimations of dominant poles (see Step 5).

Note that in Algorithm 1, formula (8) is performed step by step for every single nonzero delay element of vectors τp and τc indexed by l = μ, μ − 1, …, 2, 1. Thus, the value of τc =τ⁰ is successively updated followed by the computation of the updated estimation of sc in the vicinity of sp: see Sections 3.3–3.5 for more detail about Step 5 of the algorithm. The eventual delay vector τ⁰ is computed in the vicinity of sc according to Step 5. A nonzero RT are made.

Remark 1

Due to Proposition 2, ∇s(τ) , not its differentiability. The zero-point τ⁰ in points τ¹ and τ² as

   (9)

Then, Step 6 of for l = μ, μ − 1, …, 2, 1 according to Eq. (9).

In contrast to the use of Newton’s method, both the values τp, τc (and also the corresponding sp, sc) have to be necessarily updated when going through the set of particular delay elements, l = μ, μ − 1, …, 2, 1. Namely, if jl≠0 for a particular l, the following setting is made

Otherwise, whenever jl≠0, the use of Eq. ((and the value of l is decremented).

Then, the updated estimations of the corresponding dominant poles sp, sc are computed in the neighborhood of the current value of sp (note that its initial value is given by the output of Step 5 of Algorithm 1).

The algorithm supposes that sc is sufficiently close to sp while grid stepping. This assumption holds almost everywhere; however, a discontinuity or abrupt changes in the value of α(τ) can rarely appear due to Proposition 2. Thus, if there is a pair of nondominant poles with the real part very close to α(τ), or multiple dominant poles appear, or the function α(τ) is close to its minimum, numerical experiments imply that the estimation of sc ought to be reset—one can simply determine the rightmost root (rather than the closest one) as the updated dominant pole, or use the QPmR [19] in such a case.

Regarding NTDSs, as mentioned earlier, if condition (3) in Proposition 3 does not hold, the system is not strongly stable, and thus, there is no reason to analyze exponential stability since Eq. (3) does not depend on the particular value of τ. On the contrary, the safe upper bound c (defined below Proposition 1) is a function of τ. It, moreover, holds that there is only a finite number of poles right from c and γ, respectively. Only the approximated spectrum lying right from these values can be considered when making a decision whether the system is exponentially stable.

A more detailed view of Step 5 follows, that is, the task how to estimate the value of sc in the vicinity of sp by means of a characteristic quasipolynomial approximation is presented and its possible solutions are suggested.

3.3 Continuous-time approximation

This technique, in general, iteratively approximates Δ(s, τ) by a polynomial

   (10)

, where

which gives the following linear matrix equation to determine polynomial coefficients:

are unknown coefficients and C, R , which is inspired by the Rekasius transformation applied in the original CTCR algorithm [8].

Step 5 starts with the setting s0 = sp, followed by the computation of Eq. (10) with τ being is determined where sk , which yields the updated estimation of the dominant system pole. The reset s0 = s, where i is obtained as the output of Step 5 of Algorithm 1.

3.4 Discrete-time approximation

Another polynomization technique is inspired by digital filter design principles. The approximating polynomial has the form

   (11)

where z expresses the variable from the z-transform, ts stands for the sampling period, and z0 is the image of s0 in the z-plane given by the relation

   (12)

that preserves stability.

in every iteration step based on the estimated value of the dominant pole s0 corresponding to the nearest grid node. All roots zk are consequently computed within the iteration step. Then, transformation (12) is applied to every single root zk, and the updated estimation of sstems from the very last iteration step.

is done in two separate steps:

(1)Variables s in Δ(s, τ) expressing derivatives (i.e., s-powers) are subjected to the following Tustin (bilinear) transformation with prewarping:

(13)

where q means the shifting operator that corresponds to z−1 and ω0 is the frequency of undamped oscillations. In contrast to the simple bilinear transformation, the use of prewarping preserves frequencies of the original system after its digital approximation; that is, ωshould coincide. It was proved by examples that the use of Eq. (13) gives better results when estimating the dominant poles pair [17]. Note that by adopting a second-order finite-dimensional model, one can write ω0 = |s0|.

(2)Exponential terms expressing delays are subjected to

(14)

for a general case; thus, it is necessary to apply an integer approximation of the exponent (power). The following theorem can be used for the quadratic extrapolation.

Theorem 1

Pekař et al. [17] and Pekař and Prokop [16]

as

   (15)

, r = k − f ∈ (0, 1).

The remaining and a rather ambiguous task is how to set ts properly. Note that from the point of view of the z-transform (see Eq. 12) and due to the computer representation of real numbers, it is usually recommended for periodic systems that 0.2 ≤ ω0ts ≤ 0.5; however, the value of ts should be sufficiently small due to the discretization of derivatives and delta model theory [23]. In fact, the Tustin formula represents the trapezoidal (feedback-feedforward) rule of the derivative approximation, from the operator point of view. A suitable trade-off between these two opposing requirements can effectively be found via numerical experiments applied to a particular problem.

3.5 Approximation of neutral quasipolynomials

As mentioned earlier, it is desirable to know the approximate loci of the rightmost subset of Σess for a strongly stable NTDS that constitutes a vertical strip in the complex. The data can serve as the estimation of γ. This task can be formulated in terms of the searching of an approximating polynomial (or, an exponential polynomial with commensurate delays), for which the rightmost roots subset can easily be found. This endeavor is motivated by the fact that a couple of analytic results on the spectral distribution for NTDSs with commensurate delays have recently been published, see, for example, Refs. [of dn(sof Σess which cannot express the exact true essential spectrum.

has already been presented in Ref. [of an NTDS, which can be used to update the value of s0 in Step 5 of Algorithm 1.

is proposed.

3.5.1 Associated characteristic exponential polynomial approximation

is an integer multiple of τ0. In addition, τ, for some real k0 and the dominant pole s0, ess from the essential spectrum. A suitable value of k0 can be found by numerical experiments for a particular dn(s, τ).

As a consequence, the following optimization problem is to be solved:

   (16)

is determined, dn(s, τ) is repeatedly (iteratively) subjected to Eq. (, according to Eq. (12), to get commensurate delays. The eventual approximating exponential polynomial (in each iteration step) can be then expressed as

into nonnegative zand use Eq. (12) (with τs = τ0) to get the essential spectrum approximation in the s-plane. The estimation of s0, ess is then done in the same way as presented in the last paragraph of Section 3.3, that is, as the closest current dominant pole estimation to the preceding