GézaPattantyús-ÁbrahámDoctoralSchoolofMechanicalEngineeringSciences,BudapestUniversityofTechnologyandEconomics2019 Submittedto Supervisor :TamásINSPERGER,DSc. Author :DávidHAJDU Robuststabilityofdelayeddynamicalsystems

(1)

Budapest University of Technology and Economics Department of Applied Mechanics

Robust stability of

delayed dynamical systems

Author : Dávid HAJDU

Supervisor : Tamás INSPERGER, DSc.

Submitted to

Géza Pattantyús-Ábrahám Doctoral School of Mechanical Engineering Sciences, Budapest University of Technology and Economics

2019

(2)

(3)

To My Family

(4)

(5)

Acknowledgement

First of all, I would like to thankTamás Inspergerfor his leadership and support during the last six years. His patience and guidance helped me in all the time of my research. Moreover, I thank Prof. Gábor Stépán, former head of department, giving me an opportunity to be part of such an inspiring team.

I also thank all of mycolleaguesat the Department of Applied Mechanics for helping me in various problems. In particular, I would like to mention Dániel Bachrathy, Zoltán Dombóvári, Tamás Molnár, and Szabolcs Berezvai, who helped me the most to overcome the emerging diffi- culties. Hereby, I also would like to mention my foreign collaborators,Jin I. GeandGábor Orosz, who kept helping me through the last years.

Last, but not the least, I would like to thank myfamilyand mywife, Alexandra,for supporting me spiritually throughout my work.

This work has been supported by the ÚNKP-16-3-I. and ÚNKP-17-3-I.

New National Excellence Program of the Ministry of Human Capacities.

The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP/2007-2013)/ERC

Advanced Grant Agreement n. 340889.

The research reported in this paper was supported by the Higher Education Excellence Program of the Ministry of Human Capacities in the frame of Biotechnology research area of Budapest

University of Technology and Economics (BME FIKP-BIO).

v

(6)

(7)

Introduction

The mathematical investigation of time-delay systems plays an important role in applied mathematics and engineering. Delays are often used to model information propagation in biological and electrical systems, growth rate in population models, economical processes, and also mechanical systems. There exist several applications in mechanics, where the inclusion of the delayed dynamics leads to a better description of the real physical system, e.g., wheel shimmy, machining vibrations, traffic models or human motion control, just to mention a few. Over the years the number of contributions on these research areas has increased significantly and the tendency predicts a strong interest both from the academic and the industrial sides.

Mathematical representations of dynamical systems are always associated with a certain level of model simplification. The simplified models intend to capture the most important features of the real system, however, these reductions often lead to an impaired representation of the real dynamics. One of the most important question is how the model behaves if some of the inputs or parameters are slightly changed.

Uncertainties in the various inputs and parameters have significant effect on the systems’

performance. Therefore, it is necessary to investigate that in what extent the systems’ behavior depends on the uncertainties. Systems, which guarantee the required performance in case of unmodeled dynamics and uncertainties, are called robust. When the stability is preserved in the presence of parameter uncertainties, then the system is said to berobustly stable. This dissertation is devoted to the robust stability analysis of a certain class of delayed dynamical systems with engineering applications. Different models with static uncertainties are investigated, i.e., the parameter uncertainties are assumed to be constant in time. The dissertation is structured as follows.

Chapter 1 introduces the fundamental properties of time-delay systems, and the basic definitions of stability and robustness. Chapter 2 gives an introduction to the application of the classical Smith predictor, which is one of the simplest type of predictor feedback controller. It is shown that static uncertainties can drastically reduce the applicability of the classical Smith predictor. Chapter 3 gives a detailed description of a human balance control model in the frontal plane, where the uncertainties originate from anthropometric data of the human body.

Chapters 4, 5, and 6 are motivated by industrial applications: stability analysis of turning and milling operations. In these chapters time delays are introduced by the regenerative effect. In case of active chatter control the dead-time of the controller results and additional delay in the feedback loop, which interferes with the regenerative delay. Uncertainties in machining operations arise because of the partially known modal parameters of the dynamical system.

Finally, Chapter 7 presents a robust controller design for connected cruise controllers in auto- mated vehicles. Uncertainties in the control parameters and reaction time delays are considered in the human-driven vehicles, and their effect on the controller design is investigated in order to guarantee safe operation.

ix

(10)

(11)

CHAPTER 1 Mathematical background

Since the appearance of differential calculus, dynamical systems have been described by ordinary differential equations (ODE) [1], where the rate of change of state depends on the present state only.

As mentioned already, there exists several other systems, where past events also play important role in the dynamics. Delay arises in many engineering applications, such as control systems [2], wheel shimmy [3, 4], machining vibrations [5, 6], traffic models [7] or human motion control [8].

These systems are described by functional differential equations (FDE), which have an extensive literature. Large number of books and journal papers are available discussing the theories and applications of time-delay systems, see the books [1, 6, 9, 10, 11, 12, 13], just to mention a few.

FDEs can be categorized into three main groups depending on the highest derivative of the delayed states. When the rate of change of state depends only on past state, the system is called retarded and is described by retarded functional differential equations (RFDE). In the engineering literature, RFDEs are also often called simply delay-differential equations (DDE). If the rate of change of state depends on its own past values as well, then the corresponding equation is called a neutral functional differential equation (NFDE). When the rate of change of state depends also on past values of higher derivatives, then the system is described by an advanced functional differential equation (AFDE). While RFDEs and NFDEs are often used to model dynamical systems, AFDEs are rarely used in engineering applications due to their inverted causality [1].

In this introductory chapter a brief overview is given on the stability of RFDEs. First, Sec. 1.1 introduces the basic theories of autonomous and time-periodic DDEs. Second, Sec. 1.2 gives the most important definitions of robustness and robust stability based on spectral properties.

1.1 Stability of retarded functional differential equations

This present dissertation investigates models, where the governing equations are RFDEs. Systems of this class are generally written as

˙

y(t) =F(t,y_t), y(t) =˙ dy(t)

dt , (1.1)

wherey(t) ∈ Rⁿ is the vector of state variablesy_i(t), i = 1, . . . , n, dot stands for differentiation with respect to timet, n ∈ Z⁺ is the dimension of the state vector,y_t ∈ Bis a function segment defined by the time shift

y_t(ϑ) :=y(t+ϑ), ϑ∈[−ξ,0], (1.2) and ξ > 0 is the length of the delay interval. Here, B = C([−ξ,0],Rⁿ) is the Banach space of continuous real functions, provided with the supremum norm [11]. In (1.1), F: (R× B) → Rⁿ is a time-dependent functional mapping y_t to the space Rⁿ. Equation (1.1) is said to be linear if F(t,y_t) = L(t)y_t + h(t), where L(t) : (R × B) → Rⁿ is linear; linear homogeneous if

1

(12)

2 CHAPTER 1 Mathematical background Rⁿ 3 h(t) ≡ 0, and linear nonhomogeneous if h(t) 6= 0. Equation (1.1) is autonomous if F(y_t, t) =G(y_t), whereG:B →Rⁿdoes not depend on time explicitly [10].

Stationary solution of (1.1) can be a fixpoint y(t) ≡ y_eq, or a periodic orbit y(t) = y_p(t) with periodT, such that y_p(t+T) = y_p(t). Assuming small perturbationx(t) = y(t)−y_p(t) around the stationary solutiony_p(t), and neglecting the higher-order nonlinear terms yield a linear homogeneous RFDE.

1.1.1 Linear autonomous RFDEs

In this subsection we give the most important properties of linear autonomous RFDEs, and adapt the notations used by [1, 10, 12, 14, 15] and others. In this case, the general form can be written as

˙

x(t) =Lx_t, (1.3) where L : B → Rⁿ is linear. The solution can be given by introducing the solution operator T(t) : (R× B)→ B, which maps the initial function segmentx₀tox_tas

x_t =T(t)x₀. (1.4)

Note, that T(t) is a strongly continuous semigroup of operators [11, 12], which satisfies the propertyT(t+s) = T(t)T(s)andT(0) =I, whereIis the identity operator. Sometimes, this solution operator is also called the operator of translation along trajectories [12].

Let A : B → B be the infinitesimal generator of the solution operator, which introduces the operator differential equation form of (1.3) as [10, 11, 16]

x˙_t=Ax_t. (1.5)

The solution of linear systems can be sought in the form of complex exponential functions as x(t) =ce^λt, wherec∈Cⁿ\{0}. Therefore, substitution of the function segments

x_t(ϑ) =ce^λ(t+ϑ)=ce^λϑe^λt =s(ϑ)e^λt, (1.6)

x˙_t(ϑ) =λx_t(ϑ), (1.7)

into (1.5) yield an infinite-dimensional eigenvalue-eigenvector problem λI−A

s=0. (1.8)

The stability is determined by the spectrum of the infinitesimal generatorA. The elements of the spectrum are called thecharacteristic exponentsand are defined by the kernel of (1.8) as

ker λI−A

6

={0}. (1.9)

The stationary solutionx(t)≡0of (1.3) is said to be asymptotically stable if all of the characteristic exponents have negative real part, i.e.,

Re(λ_i)<0, for alli, (1.10)

where indexirefers to thei^thcharacteristic exponent.

Example 1. The autonomous time-delay system with a single point delay can be written as x(t) =˙ Ax(t) +Bx(t−τ), (1.11)

(13)

1.1 Stability of retarded functional differential equations 3 where x(t) ∈ Rⁿ, A ∈ R^n×n, B ∈ R^n×n and τ > 0 is the time delay. Substitution of the trial solutionx(t) = ce^λtyields the characteristic equationD(λ) = 0with the characteristic function

D(λ) := det λI−A−Be⁻^λτ

, (1.12)

whereIis then×nidentity matrix.

Note, that (1.12) is transcendental, thus it has infinitely many zeros, which give infinitely many characteristic exponents λ. Even though the characteristic exponents cannot be calculated analytically in general, there exists several numerical methods to investigate the exponential stability, see, e.g., [1, 15, 16]. Alternatively, theD-subdivisionmethod can be used by substitutingλ =ζ+iω (ζ ∈ R, ω ∈ R) into the characteristic function D(λ). This way, the stability boundaries often can be determined analytically [1]. The solution for the characteristic equationD(iω) = 0results the so-called transition curves, where an eigenvalue (or a complex pair of eigenvalues) cross the imaginary axis. The number of unstable characteristic exponents, i.e., the number of zeros ofD(iω) located on the open right-half complex plane, can be determined usingStepan’s formulas[6].

1.1.2 Linear time-periodic RFDEs

The general form of linear time-periodic RFDEs can be written as

˙

x(t) = L(t)x_t, (1.13) whereL(t) : (R× B)→Rⁿis linear and periodic with periodT, that is,L(t+T) =L(t). There exists also an operator differential equation form of (1.13) as

˙

x_t=A(t)x_t, and A(t+T) =A(t), (1.14) where A(t) : (R× B) → B. The solution can be given by introducing the solution operator U(t, t₀) : (R×R× B)→ B,t ≥t₀(also called theevolution operatoraccording to [16]), which maps the initial function segmentx_t₀ tox_t, i.e.,

x_t =U(t, t₀)x_t₀. (1.15)

The second parameter is necessary since the initial time has to be considered. Similarly to the solution operator of nonautonomous systems, U(t, t₀) satisfies the properties U(t, t₀) = U(t, s)U(s, t₀), t ≥ s ≥ t₀, and U(t, t) = I. Moreover, due to periodicity, U(t +T, t) = U(t, t−T)also holds [16]. According to the Floquet theory, stability of the stationary solution can be investigated through the spectral properties of the monodromy operator U(T,0), which maps the solution segmentx₀to the segmentx_T one period later, i.e.,

x_T =U(T,0)x₀. (1.16)

It is also known that the trial solution is sought in the formx(t) = p(t)e^λt, wherep(t+T) = p(t) is periodic [17]. Therefore, the solution segments can be written as

x_t(ϑ) =p(t+ϑ)e^λ(t+ϑ)=p_t(ϑ)e^λ(t+ϑ)=s(ϑ)e^λt, (1.17) x_t+T(ϑ) =p(t+T +ϑ)e^λ(t+T^+ϑ) =s(ϑ)e^λte^λT, (1.18) where s(ϑ) = p_t(ϑ)e^λϑ, and e^λT =: η ∈ C is called the characteristic multiplier of the system (1.13). Substitution of (1.17-1.18) with setting t = 0 into (1.16) yields the infinite-dimensional eigenvalue-eigenvector problem in the form

ηI−U(T,0)

s=0. (1.19)

(14)

4 CHAPTER 1 Mathematical background The characteristic multipliers are determined by the kernel of (1.19), i.e.,

ker ηI−U(T,0)

6

={0}, η6= 0. (1.20)

The trivial solution is said to be asymptotically stable if all of the eigenvalues of the monodromy operator are inside the unit circle of the complex plane, i.e.,

|η_i|<1, for alli, (1.21)

where indexirefers to thei^thcharacteristic multiplier.

Example 2. The time-periodic system with a single point delay can be written in the form x(t) =˙ A(t)x(t) +B(t)x(t−τ), A(t+T) =A(t), B(t+T) =B(t). (1.22) where x(t) ∈ Rⁿ, A(t) ∈ R^n×n, B(t) ∈ R^n×n and τ > 0 is the time delay. Stability can be investigated based on the approximation of the monodromy operator, see, e.g., [1, 15].

1.2 Sensitivity and robustness

Analysis of stability of systems subjected to uncertain dynamics requires the extension of the basic mathematical definitions. In case of linear systems when no uncertainty is present, stability can be assessed using the associated point spectrum of characteristic roots. Static uncertainty spoils this favorable property and introduces the continuous pseudospectra. Throughout the dissertation, we assume that the nominal system is fully known, moreover uncertainty is bounded and static in time.

1.2.1 Pseudospectrum

In order to describe the closeness to instability one may use a measure defined in terms of perturbations or uncertainty. The computation of pseudospectra (spectral value sets) is a key issue in characterization of sensitivity and therefore became an established tool in robust analysis [12]. In more details, pseudospectrum is defined by the set of eigenvalues in the complex plane that can be reached under bounded perturbations. In order to present the basic mathematical concepts, let us consider a matrixA ∈K^n×n(K=RorC) and define its spectrum as

Λ(A) :={λ∈C: det(λI−A) = 0}. (1.23) A square matrix is said to be Hurwitz stable if all of its eigenvalues have negative real part, and said to be unstable, if there exists at least one with positive real part. Similarly, a square matrix is said to be Schur stable if all eigenvalues have magnitude less than one and unstable otherwise.

As a generalization, a matrixAis said to beCg-stable, if all of the eigenvalues are confined to the domainCg ⊂Cand unstable otherwise. The boundary of the domainCg is denoted by∂Cg.

When matrixAis subjected to uncertaintyA, i.e.,˜ A→A+ ˜A, the-pseudospectrum can be introduced as

Λ(A) :={λ ∈C: det λI−(A+ ˜A)

= 0, A˜ ∈K^n×n,kA˜k< }, (1.24) where k · k is the spectral norm [18]. The complex unstructured pseudospectrum is the most conservative estimation on sensitivity since it does not consider any topology in the matrix structure and uncertainty is assumed to be complex. This yields the definition

Λ^C(A) :={λ∈C: det λI−(A+ ˜A)

= 0, A˜ ∈C^n×n,kA˜k< }, (1.25)

(15)

1.2 Sensitivity and robustness 5 which is equivalent to

Λ^C(A) = {λ∈C:kR(λ,A)k>1/}, (1.26) where R(λ,A) := (λI−A)⁻¹ denotes the resolvent of A at λ. This provides a computable formula for the complex unstructured pseudospectrum.

Sensitivity analysis with respect to real perturbation plays a more important role in engineering applications, therefore we introduce the definition ofreal unstructured pseudospectrumas

Λ^R(A) :={λ∈C: det λI−(A+ ˜A)

= 0, A˜ ∈Rⁿ^×ⁿ,kA˜k< }. (1.27) The calculation procedure of the real pseudospectra of matrix A is more challenging, but computable formula is provided by [19]:

Λ^R(A) ={λ ∈C:υ R(λ,A)

>1/}, (1.28)

where

υ R(λ,A)

= inf

γ∈(0,1]σ₂

Re R(λ,A)

−γIm R(λ,A) γ⁻¹Im R(λ,A)

Re R(λ,A)

, (1.29)

andRe(·)andIm(·)denote real and imaginary parts, andσ₂(·)denotes the second largest singular value, see [18, 19] for more details. It is shown that function υ is unimodal inγ, therefore local minima can be easily found.

Example 3. In order to visualize the spectral properties and the different spectra with uncer- tainties, let us consider the Demmel matrixE,[E]_ij =E_ij,i, j = 1. . . , ndefined as [20]

E_ij =

(−b^j⁻ⁱ ifj ≥i,

0 otherwise, (1.30)

where b = 10andn = 5 [18]. In this case, the eigenvalues ofE areλ_1,2,3,4,5 =−1. Results are presented in Fig. 1.1, where panels (a-b) illustrate the unstructured complex and real pseudospectra for different -levels computed by definitions(1.25)and (1.28), respectively. Panels (c-d) present the eigenvalues computed for 1000 random perturbations bounded by ≤ 10⁻⁴. Note that the topology of spectra and therefore the sensitivity of system highly depends on the structure of the perturbation.

The application of spectral value sets is widely applied in mathematics, however, has not found many engineering applications yet [21]. Consider [22] as an exception, which provided a pseudospectral approach to calculate the spectral abscissa for perturbed dynamical systems. For further applications, see [12, 23, 24, 25], just to mention a few.

1.2.2 Stability radius

The definition of stability radius [26] is strongly related to the concept of pseudospectrum. It is equivalent to the smallest perturbation level at which one of the perturbed matrices become unstable. This definition gives a measure of robustness of the nominal matrix A subjected to perturbation A. Recall the definition of unstructured complex pseudospectrum (1.25), then the˜ evaluation along the (stability) boundary∂Cg yields thecomplex unstructured stability radius

r_C(A) = sup

λ∈∂Cg

kR(λ,A)k

!−1

. (1.31)

(16)

6 CHAPTER 1 Mathematical background

0 1 2

-2 -1

-2 -1 0

0 1 2

-2 -1

-2 -1 0

0 1 2

-2 -1

-2 -1 0

0 1 2

-2 -1

-2 -1 0

ǫ≤10⁻⁴ ǫ≤10⁻⁴

Re(λ) Re(λ)

Im(λ)

Im(λ) Im(λ)

Im(λ)

Complex pseudospectrum Real pseudospectrum (a)

(b)

(c)

(d)

Fig. 1.1. Complex and real pseudospectra of a Demmel matrix [20] (50.5) in [18] (kE˜k =,= 10⁻³,10⁻⁴,. . ., 10⁻⁸): (a) Complex pseudospectrum; (b) Eigenvalues of 1000 complex random perturbations with≤10⁻⁴; (c) Real pseudospectrum; (d) Eigenvalues of 1000 real random perturbations with≤10⁻⁴.

For continuous-time autonomous systems∂Cg = iω, ω ∈ R(based on Hurwitz stability) and for discrete-time (or time-periodic systems, where the stability is determined by the spectral properties of the monodormy matrix) ∂Cg = e^iϕ, ϕ ∈ [0,2π) (based on Schur stability). Similarly to real perturbation structures, thereal unstructured stability radiusis defined as

r_R(A) = sup

λ∈∂Cg

υ R(λ,A)!₋1

. (1.32)

While the definitions (1.31-1.32) are introduced for matrixA(more precisely for its resolvent), it can be easily extended to general nonlinear matrix functions [12].

The unstructured stability radii defined by formulae (1.31) and (1.32) can be conservative since the structure of matrix A is not taken into account. This feature is important in engineering applications, where canonical forms are used, as demonstrated by the following example.

Example 4. Consider the unforced single-degree-of-freedom system (spring-mass-damper system) governed by

mx(t) +¨ cx(t) +˙ kx(t) = 0, (1.33) where m is the mass, c is the damping coefficient, k is the stiffness and x(t) is the generalized

(17)

1.2 Sensitivity and robustness 7 coordinate. The first-order representation reads

x(t)˙

¨ x(t)

=

0 1

−_m^k −_m^c

| {z } A

x(t)

˙ x(t)

, (1.34)

where the system matrix A is real. The complex unstructured stability radius defined by (1.31) is associated with the complex full block perturbation of matrix A, and it therefore leads to conservative robustness bounds (i.e., the entries0and1are also perturbed in the first row).

In order to reduce the gap between the real parametric uncertainties and unstructured perturbations, the so-called structured perturbation can be introduced [23]. In many cases, the perturbation structure can be taken into account by the form

A+ ˜A =A+B∆C, (1.35)

whereA ∈Kⁿ^×ⁿ,∆∈K^p^×^qis the perturbation matrix with dimensionspandq, whileB∈Kⁿ^×^p and C ∈ K^q×n are scaling (or weight) matrices with appropriate dimensions. The complex structured stability radius of matrix A with respect to perturbations in the form B∆C can be introduced as [23]

r_C^∆(A,B,C) := sup

λ∈∂Cg

kC(λI−A)⁻¹Bk

!−1

, (1.36)

wherer^∆

C(A,B,C)gives the norm of the smallest∆∈C^p^×^q, such thatA+B∆Cis notCg-stable.

This can be extended to real structured perturbations by applying the already introduced formulae as

r^∆_R(A,B,C) := sup

λ∈∂Cg

υ C(λI−A)⁻¹B!−1

. (1.37)

Example 5. Consider an unforced single-degree-of-freedom system with structured parametric perturbationsδ₁ andδ₂ in the following form

x(t)˙

¨ x(t)

=

0 1

−_m^k −_m^c

+ 0

1

|{z}

B

δ₁ δ₂

| {z }

∆

1 0 0 1

| {z } C

! x(t)

˙ x(t)

=

0 1

−_m^k +δ₁ −_m^c +δ₂

x(t)

˙ x(t)

, (1.38) then the uncertainty affects only the second row of matrixA. In this formulation the uncertainty of parametersm,c, andk must be replaced by two other uncertain parameters: δ₁ andδ₂.

The definition of stability radius can be generalized to matrix polynomial functions and even to time-delay systems, see [27, 28]. In order to visualize structured uncertainty, let us assume additive perturbations in the system parameters denoted bym˜, c˜, ˜k. This case cannot be treated by (1.38) due to the multiple divisions bym˜.

Example 6. Let the perturbed characteristic equation of single spring-mass-damper system be given as

D(λ) + ˜D(λ) = (m+ ˜m)λ²+ (c+ ˜c)λ+ (k+ ˜k) = 0. (1.39)

(18)

8 CHAPTER 1 Mathematical background It can be shown, that(1.39)can be reformulated as

1 + mλ²+cλ+k−1

˜

mλ²+ ˜cλ+ ˜k

= 1 + mλ²+cλ+k−1

| {z } D(λ)⁻¹

λ² λ 1

| {z } w(λ)





˜ m

˜ c

˜k





| {z }

∆

= 0,

(1.40) wherew(λ)is a weight matrix and∆is the perturbation matrix (now a vector). The complex and real structured stability radii are defined as [29]

r^∆_C D(λ),w(λ)

=



sup

λ= iω ω≥0

kD(λ)⁻¹w(λ)k





−1

, (1.41)

r^∆_R D(λ),w(λ)

=



sup

λ= iω ω≥0

υ D(λ)⁻¹w(λ)





−1

, (1.42)

where r^∆_C and r^∆_R give the norm of the smallest complex or real ∆, such that the perturbed characteristic function has no zero on the imaginary axis.

1.2.3 Structured singular value

Applications of pseudospectra and stability radii are found to be effective tools, when the perturbation structures are relatively simple. However, there exist many cases when the structures are complicated and for real engineering applications tight robustness bounds cannot be computed. In- stead, in order to handle a wider class of uncertain systems, the structured singular value (µ-value) was introduced by [30], and it became one of the most effective mathematical technique in robust analysis. In many applications it is possible to represent mixed structured uncertainty problems (real and complex) in a block-diagonal form and analyze its effect on the uncertain system. The perturbation matrix∆is therefore not necessarily limited to full block matrix structures as opposed to the structured stability radii defined by (1.36) and (1.37).

Here, we define the structured singular valueµof matrixM∈C^n×mas

µ(M) :=







∆min∈X_K

¯

σ(∆) : det (I−M∆) = 0 ₋1

0, if no∆∈ XKmakesI−M∆singular,

(1.43) where∆ ∈ C^m^×ⁿ, XK is the set of allowed perturbations and σ(¯ ·) = k · krepresents the largest singular value. In other words, µ is the inverse of the largest singular value of the smallest perturbation matrix∆∈ XKfrom the set of possible structures, such thatI−M∆is singular. Or alternatively, 1 is the element ofΛ(M∆). The more general definition ofµincluding computational issues is given in Appendix A.

In order to highlight the differences between the representation of uncertainties in Sec. 1.2.2 and in Sec. 1.2.3, let us recall the single-degree-of-freedom model and assume the same parametric uncertainties in the scalar coefficients.

Example 7. Consider the forced single-degree-of-freedom system (spring-mass-damper system) governed by

(m+ ˜m)¨x(t) + (c+ ˜c) ˙x(t) + (k+ ˜k)x(t) =f(t), (1.44)

(19)

1.2 Sensitivity and robustness 9

y1

y2

y3

u1

u2

u3

(a) (b)

k

˜k c 1 m

˜ c

˜ m w=F(s)

w z=X(s)

z

u y

∆=



m˜ 0 0 0 c˜ 0 0 0 k˜





∆

N(s) 1

s 1

s

Fig. 1.2. Uncertain spring-mass-damper system: (a) Block-diagram with uncertainties; (b)N-∆uncertain intercon- nection structure.

where f(t)is the forcing. The block-diagram including the uncertainties is presented in Fig. 1.2, where w = F(s) is the input and w = X(s) is the output (F(s) and X(s) are the Laplace transforms of f(t) and x(t), respectively), moreover, the inputs and outputs of each uncertainty block (m,˜ ˜c, ˜k) are denoted by y_i and u_i. The equations can be reformulated by constructing a multiple-input and multiple-output system (MIMO) in the form





 y₁ y₂ y₃ w





=

z N(s)}| { 1

ms²+cs+k







−s² −s² −s² s²

−s −s −s s

−1 −1 −1 1











 u₁ u₂ u₃ z





, (1.45)



u₁ u₂ u₃



=





˜

m 0 0 0 ˜c 0 0 0 ˜k





| {z }

∆



y₁ y₂ y₃



. (1.46)

The transfer function between the inputz and outputwcan be given by the upper linear fractional transformation (upper LFT Fu). Let the generalized transfer function matrixN(s)be partitioned as

y w

=

M(s) N₁₂(s) N₂₁(s) N₂₂(s)

| {z } N(s)

u z

and u =∆y, (1.47)

then

w=Fu N(s),∆

z = N₂₂(s) +N₂₁(s)∆(I−M(s)∆)⁻¹N₁₂(s)

z, (1.48)

moreover

Fu N(s),∆

= 1

(m+ ˜m)s² + (c+ ˜c)s+ (k+ ˜k). (1.49) The upper LFT exists if I−M(s)∆ is invertible in (1.48). Due to the required stability of the nominal system (hence the notion of robust stability), it is sufficient to show thatI−M(iω)∆

(20)

10 CHAPTER 1 Mathematical background is invertible along the imaginary axis for all ω ≥ 0, i.e.,det(I−M(iω)∆) 6= 0for any allowed perturbation ∆ and ω ≥ 0, cf. (1.43). Notice the difference between formulations of the perturbation matrix∆ in (1.40) and (1.46). The structured singular value analysis is not limited to full block perturbation, instead, repeated scalar, full block matrix (non-parametric uncertainty), mixed perturbation structures can simultaneously be considered in the block-diagonal matrix∆, see Appendix A.

While the stability radius and structured singular value are defined in different ways, the relations between them in special cases can be easily given. Let us recall the structured complex perturbations A+B∆C, withA∈ Kⁿ^×ⁿ(stable),B ∈Kⁿ^×^p,C∈K^q^×ⁿand∆∈C^p^×^q(complex full block).

In this case the following equivalence can be given

r^∆

C(A,B,C) = sup

λ∈∂Cg

kC(λI−A)⁻¹Bk

!₋1

= sup

λ∈∂Cg

µ M(λ)!₋1

, (1.50)

whereM(λ) = C(λI−A)⁻¹B. The relation in (1.50) holds only if∆∈ C^p^×^q is a complex full block matrix (see also the upper bound in (A.4) in Appendix A).

(21)

CHAPTER 2 Smith predictor

In many control applications the presence of time delay detrimentally affects the performance of the closed loop system. The finite speed of information propagation, digital sampling and data processing are typically listed as sources of this phenomenon. Time delay is also considered to be one of the main reason of unstable behavior of control systems, which therefore has to be eliminated, or alternatively, compensated in the feedback loop.

The Smith predictor [2], introduced in 1957, is probably the best known technique to overcome the destabilizing effect of feedback delays in control systems. This chapter is devoted to stability analysis of the original Smith predictor, which has been studied deeply a number of times in the literature, see [12, 31, 32]. While many modifications of the original concept became accessible during the decades, researchers still remained interested about the initial attempt to compensate the effects of delay more adequately.

In the present study a time-domain-based analysis is performed to highlight the possible sources of low robustness of stability in the presence of model uncertainty. It is demonstrated through two case studies that the original Smith predictor is not capable to compensate the effect of time delay when the predictor model is slightly mistuned. The structure of this chapter is as follows. Section 2.1 gives a more detailed introduction to the Smith predictor using frequency and time-domain- based approaches. Section 2.2 details the mathematical solution using the D-subdivision method, while Sec. 2.2.2 and 2.2.3 present two examples using a marginally stable and an asymptotic stable system. In Sec. 2.3 the original Smith predictor is compared to the finite spectrum assignment controller. Main results are highlighted in Sec. 2.4.

2.1 The Smith predictor

It is known that state prediction is a fundamental concept for systems with feedback delay [33].

Predictive controllers intend to eliminate delay from the feedback loop by using a prediction of the actual state based on an internal model of the plant. In the ideal case, stability properties of the system subject to predictor feedback match the stability properties of the same system without delay.

A main concern about prediction-based controllers is that they require knowledge about the plant and the feedback delay. The Smith predictor [2] uses a mathematical model of the plant implemented in the controller to predict the actual and delayed states. The slightest mismatch between the internal model used for the prediction and the actual system may lead to poor performance or may even destabilize the closed-loop system in extreme cases. The conditions for practical stability (i.e. the preservation of stability for infinitesimal modeling mismatches) for the Smith predictor was given in [34] and for some special cases in [35, 36, 37]. The effect of delay mismatches was investigated in [32]. It is known that in case of a proper but not strictly proper delay-free transfer function, the sensitivity to infinitesimal parameter mismatches can be explained by the discontinuity of the

11

(22)

12 CHAPTER 2 Smith predictor

r u x

d

xd ˆxd

ˆ x

C(s) P(s)

P(s)ˆ

e⁻^{τ s} e^−ˆ^{τ s}

CSP(s)

Fig. 2.1.The block diagram of the Smith predictor.

spectrum of the associated difference equation [32, 38].

Here, we consider a system with a strictly proper delay-free transfer function and show that sensitivity to infinitesimal parameter mismatches can still occur for a marginally stable plant. Note that this sensitivity is of different nature from those analyzed in [32, 38]. The sensitivity presented here is related to the uncertainties of the plant parameters and the feedback delay. Moreover, we visualize the sensitivity of the Smith predictor to parameter mismatches by constructing stability diagrams for a marginally stable and for an asymptotically stable second-order plant. Transition between the Smith predictor and delayed state feedback is illustrated by a series of stability diagrams.

2.1.1 Frequency-domain representation

The Smith predictor is usually represented in frequency domain either by its block diagram or by its transfer function [2, 12, 31, 32]. The block diagram is shown in Fig. 2.1. As mentioned in Sec. 2.1, the point of the Smith predictor is that the feedback delay is eliminated from the control loop using a prediction of the actual and the delayed states based on an internal model of the plant.

Let us denote the transfer function of the plant by P(s), the transfer function of the plant used by the internal model byPˆ(s), the transfer function of the primary controller byC(s), the actual feedback delay byτ, the delay used by the internal model byτˆ, moreoversis the Laplace variable.

In practice, the internal model is not perfectly accurate, thereforeP(s)6= ˆP(s)andτ 6= ˆτ. The transfer function from the plant input disturbancedto the outputxcan be given as

W_dx(s) = P(s)(1 +C(s) ˆP(s)−C(s) ˆP(s)e⁻^{τ s}^ˆ )

1 +C(s) ˆP(s)−C(s) ˆP(s)e^−ˆ^{τ s}+C(s)P(s)e^{−τ s}. (2.1) If the plant and the controller are factorized as

P(s) = B₁(s)

A₁(s), Pˆ(s) = Bˆ₁(s)

Aˆ₁(s), C(s) = B₂(s)

A₂(s), (2.2)

then the transfer function reads W_dx(s) =

Aˆ₁(s)A₂(s)B₁(s) +B₁(s) ˆB₁(s)B₂(s)(1−e⁻^{τ s}^ˆ )

×

A₁(s)A₂(s) ˆA₁(s) +A₁(s)B₂(s) ˆB₁(s)(1−e^−ˆ^{τ s}) + ˆA₁(s)B₁(s)B₂(s)e^{−τ s}₋1

. (2.3) Clearly, ifAˆ₁(s) = A₁(s), then the poles of the open-loop system (which are the zeros ofA₁(s)) are the poles of the closed-loop system, too. This deduction is often used to conclude that the Smith predictor cannot stabilize an unstable plant. However, ifAˆ₁(s)6=A₁(s)then this conclusion is not valid.

(23)

2.2 Stability analysis of a second-order system 13

2.1.2 Time-domain representation

Time-domain representations are rarely discussed in the literature (see Eq. (2.45) in the book [33]

for an exception). In this subsection, the time-domain equations of the Smith predictor are presented based on its block diagram.

Time-domain representation of the Smith predictor in case of a state feedback controller can be written as

x(t) =˙ Ax(t) +Bu(t), (2.4)

˙ˆ

x(t) = ˆAˆx(t) + ˆBu(t), (2.5) u(t) = K x(t−τ)−x(tˆ −τˆ) + ˆx(t)

, (2.6)

where x ∈ Rⁿ is the vector of actual state variables,xˆ ∈ Rⁿ is the auxiliary vector of predicted state variables,A∈ Rⁿ^×ⁿ andAˆ ∈Rⁿ^×ⁿare the actual and the model state matrices,B ∈Rⁿ^×^m and Bˆ ∈ R^n×m are the actual and the model input matrices, and matrix K ∈ R^m×n contains the control gains. It is assumed now, for simplicity, that the state vectorxis fully available and there is no need for an observer.

Without loss of generality it can be assumed that the reference inputris zero (ifris not zero then the variational system around the reference inputrhas the form of Eqs. (2.4)-(2.6)). Equation (2.4) describes the actual system while Eq. (2.5) corresponds to the internal model. The corresponding control law can be given in integral form as

u(t) =K x(t−τ)−

t−ˆτ

Z

0

e^A(t^ˆ ⁻^τ^ˆ⁻^θ)Bu(θ)ˆ dθ+ Zt

0

e^A(t^ˆ ⁻^θ)Bu(θ)ˆ dθ

!

. (2.7)

Thus, the control law involves integrals of the control input over the interval[0, t]. The closed-loop system can be described by a system of retarded functional differential equations (RFDEs) with two delays (τ andτˆ) as

˙

x(t) =Ax(t) +BK x(t−τ)−x(tˆ −τˆ) + ˆx(t)

, (2.8)

˙ˆ

x(t) = ˆAx(t) + ˆˆ BK x(t−τ)−x(tˆ −τˆ) + ˆx(t)

. (2.9)

The characteristic equation is D(s) = det

sI−A−BKe⁻^{τ s} −BK(1−e⁻^{τ s}^ˆ )

−BKeˆ ^{−τ s} sI−Aˆ −BK(1ˆ −e^−ˆ^{τ s})

= 0. (2.10)

It can be seen that the closed-loop system subjected to the Smith predictor is equivalent to a delayed state feedback system with an augmented state vector of dimension double that of the open-loop system.

2.2 Stability analysis of a second-order system

Here, we consider a general stable second-order plant with control input of the form

¨

x(t) +bx(t) +˙ ax(t) =u(t), u(t) = f(x(t−τ),x(t),ˆ ˆx(t−τˆ)), (2.11) where b ≥ 0 is the damping parameter and a > 0 is called the system parameter. When the predictor is utilized in the controller, one needs to estimate the parameters in the system, that isˆb,

(24)

14 CHAPTER 2 Smith predictor ˆ

aand ˆτ are estimations ofb, aandτ, respectively. Accordingly, the system matrices of (2.4) are given by

x(t) = x(t)

˙ x(t)

, A =

0 1

−a −b

, B= 0

1

, (2.12)

the predictor model in (2.5) is given by ˆ

x(t) = x(t)ˆ

˙ˆ x(t)

, Aˆ =

0 1

−ˆa −ˆb

, Bˆ = 0

1

, (2.13)

moreover the matrix of control gains in (2.6) reads K=

−k_p −k_d

, (2.14)

wherek_p andk_d are the proportional and the derivative gains.

2.2.1 Special case: delayed state feedback

First, let us consider for simplicity that the system is undamped, i.e.,b= 0and (2.11) degrades to

¨

x(t) +ax(t) =u(t), u(t) = f(x(t−τ),x(t),ˆ x(tˆ −ˆτ)). (2.15) A crane model is presented in Fig. 2.2(a) as a example, where the damping is often very small or completely neglected [39].

If Aˆ = 0 andBˆ = 0then the predicted state is constant in time and (2.8) gives the delayed state feedback

x(t) =˙ Ax(t) +BKx(t−τ). (2.16) The associated characteristic function according to (2.10) is (omitting the internal model completely)

D(s) =s²+a+ (k_p+sk_d)e^{−τ s}. (2.17) The stability diagram of this special case can be constructed by the D-subdivision method [6]. The D-curves, which are associated with pure imaginary characteristic exponents of the forms = iω, can be given in the parametric form

ifω = 0 : k_p =−a, k_d ∈R, (2.18)

ifω 6= 0 : k_p = (ω²−a) cos(ωτ), k_d = ω²−a

ω sin(ωτ). (2.19) The number of unstable characteristic exponents in the domains separated by the D-curves can be determined using Stepan’s formulas (see Eqs. (2.17) and (2.19) in [6]). The corresponding stability diagrams for a marginally stable system (b= 0,a >0) is shown in Fig. 2.2(b) withτ = 1.

2.2.2 Smith predictor: a marginally stable system

Let us assume that the open loop system has no internal damping (b = ˆb = 0), and the Smith predictor is employed to eliminate the effect of delay from the feedback loop according to (2.4- 2.6). The stability diagrams can be determined with the same procedure as used for delayed state feedback. Based on (2.10), the characteristic function of the system is

D(s) = (s²+ ˆa)(s² +a) + (k_p+sk_d) (s²+a) + (s²+ ˆa)e^{−τ s}−(s²+a)e^−ˆ^{τ s}

. (2.20)

(25)

2.2 Stability analysis of a second-order system 15 (a) Crane model

u(t)

g

m l

x

(b) Stability diagram

kp

−20 0 20 40 −1 0

−1 1

1

1 1

−a kd

kd

−5

0 0 0

5 5 10

3

3 4

2

ω 2 ω

ω= 0 ω= 0

Stable

Unstable

Fig. 2.2. (a) Crane model as an example for a marginally stable open-loop system. (b) Stability diagrams with the number of the unstable characteristic exponents for a marginally stable open-loop system (a= 0.5,b= 0) subjected to delayed state feedback (τ= 1). Regions of asymptotic stability are indicated by gray shading.

Application of the D-subdivision method gives the D-curves in the following form. Ifω = 0, then

k_p =−a, k_d∈R. (2.21)

This result is similar to that obtained for a delayed state feedback controller. Ifω6= 0, then

k_p=

(a−ω²)(ˆa−ω²)

(a−ω²) cos(ωˆτ)−1

−(ˆa−ω²) cos (ωτ)

g₀ , (2.22)

k_d=

(a−ω²)(ˆa−ω²)

(a−ω²) sin(ωτˆ)−(ˆa−ω²) sin(ωτ)

ωg₀ , (2.23)

where

g₀ = 2a²+ â²−2(2a+ â)ω²+ 3ω⁴ + 2(a−ω²) (â−ω²)

cos(ωτ)−cos (τ −τˆ)ω

−(a−ω²) cos(ωˆτ)

!

. (2.24) The robustness of a marginally stable plant (b = 0, a > 0) in case of parameter mismatches can be investigated using parameter sweeping. A series of stability diagrams with the number of unstable characteristic exponents are presented in Fig. 2.3 for different mismatches of the system parameter and the feedback delay. These diagrams can be considered projections of the four- dimensional stability chart in the parameter space (kp, k_d,â,τ). For the ideal case, whenˆ aˆ = a andτˆ = τ, the stability boundaries are given byk_p > −a andk_d > 0, which corresponds to the stability condition of the delay-free system. In case ofâ=a, the rightmost characteristic exponents are purely imaginary (see (2.20)), consequently, the system is marginally stable (see striped gray regions in Fig. 2.3). This property is inherited from the open-loop system [40]. However, ifâ6=a then the system becomes asymptotically stable (see gray shading regions in Fig. 2.3). It can be seen that stability properties change radically ifâ > a. Even the slightest overestimation of the system parameter can destabilize the system. Slight underestimation of the system parameter, however, does not significantly affect the stability properties. The effect of the uncertainty in the time delay is smaller than that of the system parameter. In general, mismatches in the delay cause a slight reduction of the stable domains. Note that the D-curvek_p = −adoes not change for any system parameter mismatch.

(26)

16 CHAPTER 2 Smith predictor

4 4

2 2

3 2

4

3 2

4

2 2

3 2

4

3 2

4

2 2

3 2

4

3 2

4

2 2

3 2

4

3 2

4

2 2 2

3 2

4

3 2

4 2

Stable

Unstable

kp

kdkdkdkdkd

0 0

0

0 0

0 0 0 0

0 0

−1

−1 1 1 1 1 1

1

1 1 1 1

1 1 1

1 1

1

1 1 1

ˆ

a= 0.8a â= 0.9a â=a â= 1.1a â= 1.2a

ˆτ=0.8τˆτ=0.9τˆτ=τˆτ=1.1τˆτ=1.2τ

Fig. 2.3. Stability diagrams with the number of the unstable characteristic exponents for small parameter mismatches (a= 0.5,b= ˆb= 0,τ = 1). Regions of asymptotic stability and marginal stability are indicated by solid gray shading and by striped gray shading, respectively.

2

4

2 4

2

4

2 4

Stable

Unstable

kp

0 0 0

0 0

0 −0.05 −0.05 −0.05 −0.05

−0.05

−0.05 0.05 0.05 0.05 0.05

0.05

kd

ˆ

a= 0.9a aˆ= 0.99a â=a â= 1.01a â= 1.1a

Fig. 2.4. Transition of the stability boundaries (a= 0.5,b= ˆb= 0,τˆ=τ = 1). Regions of asymptotic stability and marginal stability are indicated by solid gray shading and by striped gray shading, respectively.

Figure 2.4 shows the transition of the stability diagram for small mismatches of the estimated system parameterâ, whena, τandτâre kept constant. Forâ= 0.9athere is a small loop attached to the origin associated with 4 unstable characteristic exponents. Asâ→a, the loop gets smaller and smaller and disappears ataˆ=a. Ifâis just larger thana, then the stability boundaries turns inside out, the small loop becomes stable and the domain which was stable forâ ≤abecomes unstable.

This demonstrates that the Smith predictor is sensitive to infinitesimal parameter uncertainties for marginally stable plants.

(27)

2.3 Stability analysis of a second-order system 17

2.2.3 Smith predictor: an asymptotically stable plant

In the previous section, it was demonstrated that the Smith predictor is not robust stable to infinitesimal parameter mismatches for a marginally stable plant. Dynamical systems in engineering applications often involve damping in the governing equations, which results in an asymptotically stable plant. In this section we analyze the same second-order plant with a slight damping term.

We consider the system (2.4-2.6) with A =

0 1

−a −b

and Aˆ =

0 1

−ˆa −ˆb

. (2.25)

Robustness of this system is investigated using the same concept as before. A series of stability diagrams are presented in Fig. 2.5 for different mismatches of the system parameter and the damping parameter. The damping parameter is b = 0.05, which gives the relative damping ratio of about 3.5%. As opposed to the marginally stable plant, the damped system is not sensitive to infinitesimal parameter mismatches. Still, it is sensitive to finite parameter mismatches: the stable domain shrinks significantly for finite overestimation of the system parametera.

Figure 2.6 shows the transition of the stability diagram as the estimated system parameter â changes. It can be seen that the stable domains shrinks radically to a finite region whileâchanges from1.105ato1.115a. At the critical valueâ_cr, the loop with four unstable poles turns inside out, and the stable domain starts shrinking radically. The stability boundaries associated toτˆ= τ and ˆb =bare given by the parametric curves

k_p = 1 g₁

aˆa−(a+ ˆa+b²)ω²+ω⁴

ω²−a+ (a−ˆa) cos(ωτ) +

bω(−a−ˆa+ 2ω²) bω+ (a−ˆa) sin(ωτ) (2.26) k_d = 1

ωg₁

−bω(a²−2aω²+b²ω²+ω⁴) + (a−a)ˆ bω(a+ ˆa−2ω²) cos(ωτ)+

(aˆa−(a+ ˆa+b²)ω² +ω⁴) sin(ωτ)

, (2.27)

where

g₁ = ω²−a+ (a−ˆa) cos(ωτ)2

+ bω+ (a−ˆa) sin(ωτ)2

. (2.28)

The loop turns inside out whenlim_ω_→_ω_cr₋₀ k_p = ∞andlim_ω_→_ω_cr₋₀k_d =∞. This happens when the denominatorg₁ becomes zero, i.e.

ω_cr² −a+ (a−ˆa_cr) cos(ω_crτ) = 0 (2.29) bω_cr+ (a−ˆa_cr) sin(ω_crτ) = 0, (2.30) which gives a nonlinear system of equations in the form

bω_crcos(ω_crτ) + (a−ω_cr²) sin(ω_crτ) = 0 (2.31) a+ (ω_cr² −a) cos(ω_crτ) +bω_crsin(ω_crτ) = ˆa_cr. (2.32) For the case τˆ = τ = 1, ˆb = b = 0.05and a = 0.5, numerical solution of (2.31) gives ω_cr = 0.73528, and (2.32) yieldsaˆ_cr = 0.5548. Therefore the critical overestimation isˆa_cr/a= 1.1096, as shown in Fig. 2.6. Similarly to the marginally stable plant, the stability diagram tends to that of the delayed state feedback controller (DSF, denoted by dashed line) if the system parameter is strongly overestimated.

GézaPattantyús-ÁbrahámDoctoralSchoolofMechanicalEngineeringSciences,BudapestUniversityofTechnologyandEconomics2019 Submittedto Supervisor :TamásINSPERGER,DSc. Author :DávidHAJDU Robuststabilityofdelayeddynamicalsystems

Robust stability of

delayed dynamical systems

Author : Dávid HAJDU

Supervisor : Tamás INSPERGER, DSc.

Submitted to

Géza Pattantyús-Ábrahám Doctoral School of Mechanical Engineering Sciences, Budapest University of Technology and Economics

2019

Acknowledgement

Contents

Introduction

CHAPTER 1

Mathematical background

1.1 Stability of retarded functional differential equations

1.1.1 Linear autonomous RFDEs

1.1.2 Linear time-periodic RFDEs

1.2 Sensitivity and robustness

1.2.1 Pseudospectrum

1.2.2 Stability radius

1.2.3 Structured singular value

CHAPTER 2

Smith predictor

2.1 The Smith predictor

2.1.1 Frequency-domain representation

2.1.2 Time-domain representation

2.2 Stability analysis of a second-order system

2.2.1 Special case: delayed state feedback

2.2.2 Smith predictor: a marginally stable system

2.2.3 Smith predictor: an asymptotically stable plant