- Email: [email protected]
Real Curve Analysis and Stability of Time-delay Systems
7th IFAC Symposium on System Structure and Control 7th IFAC Symposium on System Structure and Control Sinaia, Romania, September 9-11,Structure 2019 7th IFAC Symposium on System System and Control Control 7th IFAC Symposium on and Available online at www.sciencedirect.com Sinaia, Romania, September 9-11,Structure 2019 7th IFAC Symposium on System and Control Sinaia, Romania, September 9-11, 2019 Sinaia, Romania, September 9-11,Structure 2019 Sinaia, Romania, September 9-11, 2019
ScienceDirect
IFAC PapersOnLine 52-17 (2019) 94–98
Real Real Real Real
Curve Analysis and Stability Curve Analysis and Stability Curve Analysis and Stability Time-delay Systems Curve Analysis and Stability Time-delay Systems Time-delay Systems Time-delay Systems ∗∗ Yacine Bouzidi ∗∗ Adrien Poteaux ∗∗
of of of of
Yacine Bouzidi ∗∗ Adrien Poteaux ∗∗ Yacine Bouzidi Bouzidi ∗ Adrien Adrien Poteaux Poteaux ∗∗ Yacine Yacine Bouzidi Adrien Poteaux ∗∗ ∗ ∗ Inria Lille - Nord Europe, Parc Scientifique de la Haute Borne, 40 Inria Lille Nord Europe, Parc Scientifique de la Haute Borne, 40 ∗ ∗ Inria Lille -- Nord Europe, Parc Scientifique de Borne, 40 Avenue Halley, Bat A -- Park Plaza, 59650 Villeneuve d’Ascq, France Nord Europe, Parc Scientifique de la la Haute Haute Borne, 40 ∗ Inria Lille Avenue Halley, Bat A Park Plaza, 59650 Villeneuve d’Ascq, France Inria Lille Nord Europe, Parc Scientifique de la Haute Borne, 40 Avenue Halley, Bat A Park Plaza, 59650 Villeneuve d’Ascq, France (e-mail: [email protected]). Avenue Halley, Bat A - Park Plaza, 59650 Villeneuve d’Ascq, France (e-mail: [email protected]). ∗∗ Avenue Halley, Bat A - INRIA, Park Plaza, 59650 Villeneuve d’Ascq, France (e-mail: [email protected]). UMR 9189 CRIStAL Centre de ∗∗ Univ. Lille, CNRS, (e-mail: [email protected]). INRIA, UMR 9189 - CRIStAL - Centre de ∗∗ Univ. Lille, CNRS, (e-mail: [email protected]). ∗∗ Univ. Lille, CNRS, INRIA, UMR 9189 -- CRIStAL CRIStAL Centre de Recherche en Informatique Signal et de F-59000 CNRS, INRIA, UMR 9189 -- Centre de ∗∗ Univ. Lille, Recherche en Informatique Signal et Automatique Automatique de Lille, Lille, F-59000 Univ.Lille, Lille, CNRS,(e-mail: INRIA, UMR 9189 - CRIStAL - Centre de Recherche en Informatique Signal et Automatique de Lille, F-59000 France [email protected]). Recherche en Informatique Signal et Automatique de Lille, F-59000 Lille, France (e-mail: [email protected]). Recherche en Informatique Signal et Automatique de Lille, F-59000 Lille, France (e-mail: [email protected]). Lille, France (e-mail: [email protected]). Lille, France (e-mail: [email protected]). Abstract: In this work, we investigate investigate the the asymptotic asymptotic stability stability of of LTI LTI differential differential commensurate commensurate Abstract: In this work, we Abstract: In this work, we investigate the asymptotic stability of LTI differential commensurate time-delay systems whose dynamics are defined in the frequency domain by quasi-polynomials of Abstract: systems In this work, we investigate the asymptotic stability ofdomain LTI differential commensurate time-delay whose dynamics defined in the frequency by quasi-polynomials of m we −j τ sare Abstract: Inτthis work, investigate the asymptotic stability ofdomain LTI differential commensurate time-delay systems whose dynamics are defined in the frequency domain by quasi-polynomials of m the form f (s, ) = p (s) e . We propose a new approach for studying the stability with −j τ s j time-delay systems whose dynamics are defined in the frequency by quasi-polynomials of j=0 the form f (s, τ) = (s) e−j . Wedefined propose new approachdomain for studying the stability with m τs sare m j=0 p time-delay systems whose in aaathe frequency bythe quasi-polynomials of τ the formtoff (s, (s, = pjjjdynamics (s) ee−j . We propose new approach for studying the stability stability with with respect the we determine asymptotic behavior of the roots of quasi-polynomial the form ττdelay: )) = p (s) propose new approach for studying the j=0 m −j τ the s . We j=0 respect to the delay: we determine the asymptotic behavior of the roots of the quasi-polynomial the form f (s, τ ) = p (s) e . We propose a new approach for studying the stability with j j=0 respect to the the delay:pairs we determine determine theanalysing asymptotic behavior of the the roots of the the quasi-polynomial quasi-polynomial near to critical of the intersection of an algebraic respect to delay: we the asymptotic of of near to the the critical of ff (s, (s, ττ )) by by thebehavior intersection of roots an associated associated algebraic space space respect to the delay:pairs werefined determine theanalysing asymptotic behavior ofthe the roots ofmethods, the quasi-polynomial near to the critical pairs of ff (s, ττ )) by analysing the intersection of an associated algebraic space curve with sufficiently 3D real boxes. Compared to existing our near to the critical pairs of (s, by analysing the intersection of an associated algebraic space curvetowith sufficiently refined 3D realanalysing boxes. Compared to theofexisting methods, our method method near the critical pairs of f (s, τ ) by the intersection an associated algebraic space curve with sufficiently refined 3D real boxes. Compared to the existing methods, our method does not require any Puiseux series computation. curve with sufficiently refinedseries 3D real boxes. Compared to the existing methods, our method does not require any Puiseux computation. curve with sufficiently refined 3D real boxes. Compared to the existing methods, our method does not not require require any any Puiseux Puiseux series series computation. computation. does Copyright 2019. The Authors. Published Elsevier Ltd. All rights reserved. does not © require any Puiseux seriesby computation. Keywords: Time delay systems, curve Keywords: Time delay systems, curve analysis, analysis, stability, stability, certified certified tests, tests, symbolic-numeric symbolic-numeric Keywords: Time delay systems, curve analysis, stability, certified tests, methods. Keywords: Time delay systems, curve analysis, stability, certified tests, symbolic-numeric symbolic-numeric methods. Keywords: Time delay systems, curve analysis, stability, certified tests, symbolic-numeric methods. methods. methods. 1. INTRODUCTION INTRODUCTION yields 1. yields aa decomposition decomposition of of the the ττ -axis -axis into into aa set set of of disjoint disjoint 1. yields aa decomposition decomposition of the the ττ -axis -axis intoisaaeither set of ofstable disjoint intervals in which the considered system or 1. INTRODUCTION INTRODUCTION yields of into set disjoint intervals in which the considered system is either stable or 1.investigate INTRODUCTION yields a decomposition of the τ -axis intois aeither set ofstable disjoint intervals in which the considered system or unstable. In this paper, we the asymptotic stability of intervals in which the considered system is either stable or In this paper, we investigate the asymptotic stability of unstable. intervals in which the considered system is either stable or unstable. In this paper, we investigate the asymptotic stability of LTI differential commensurate time-delay systems whose unstable. In this paper, we investigate the asymptotic stability of Algorithmically, the above approach can be decomposed LTIthis differential commensurate time-delay systems whose unstable. In paper, we investigate the asymptotic stability of the above approach can be decomposed LTI commensurate time-delay systems dynamics are defined defined by the the following following state-space repre- Algorithmically, LTI differential differential commensurate time-delay systems whose whose Algorithmically, the above can be in First, the of dynamics are by state-space repreAlgorithmically, the problems. above approach approach be decomposed decomposed LTI differential commensurate time-delay systems whose in two two independent independent problems. First, can the computation computation of dynamics are defined by the following state-space representation dynamics are defined by the following state-space repreAlgorithmically, the above approach can be decomposed in two independent problems. First, the computation of the critical pairs of f (s, τ ), and second, for each critical sentation m in two independent the for computation of dynamics the critical pairs of fproblems. (s, τ ), andFirst, second, each critical mthe following state-space representation sentation are defined by in two independent problems. First, the computation of the critical pairs of f (s, τ ), and second, for each critical pair (s , τ ), the study of the asymptotic behavior of the m 0 0 the critical pairs of f (s, τ ), and second, for each critical x(t) ˙ = A x(t − k τ ) (1) m sentation k pair (s , τ ), the study of the asymptotic behavior of the 0 0 x(t) ˙ = (1) the critical pairs of f (s, τ ), and second, for each critical m Ak x(t − k τ ) pair (s , τ ), the study of the asymptotic behavior of the root ss000 with to small variation the τthe 0 .. , τ00 ), respect the study the asymptotic x(t) ˙˙ = (1) x(t) = k=0 A Akk x(t x(t − −k k ττ )) (1) pair root (s with to aaof small variation of ofbehavior the delay delayof 0 pair (s , τ0 ), respect the study the asymptotic behavior ofτττthe x(t) ˙{τ ∈=Rk=0 A≥ −mk ∈ τ )N := {0, 1, . (1) root ss0000 with respect to aaofsmall variation of the delay k x(t 0 .. root with respect to small variation of the delay k=0 0 := | τ 0}, . .} where τ ∈ R k=0 + + There exist numerous methods for computing the critical root s with respect to a small variation of the delay τ | τ ≥ 0}, m ∈ N+ := {0, 1, . . .} There 0exist numerous methods for computing the critical where τ ∈ R+ := {τ ∈ Rk=0 0. :=are {τ ∈ ∈ R || ττ ≥ ≥matrices 0}, m m∈ ∈ with N+ :=entries {0, 1, 1, in . .. .} .} where τ, .∈ ∈. .R R, A and A constant aa There + 0τ m exist numerous methods for computing the critical := {τ R 0}, N := {0, . where pairs of f (s, τ ) (see for instance Bouzidi et al. [2016], + + There exist numerous methods for computing the critical and A , . . . , A are constant matrices with entries in 0τ ∈ R+ m pairs of f (s,numerous τ ) (see for instance Bouzidi et the al. critical [2016], :=are {τ ∈ RThe | τ ≥matrices 0}, m ∈ with N+function :=entries {0, 1,ofin . . .} where There exist methods for[2002]). computing and A A are constant matrices with entries in field K K Q, pairs of ff[2001], (s, for instance Bouzidi et al. and ,, .. .. .. ,, A constant aa Niculescu Olgac and Sipahi Bouzidi et m= pairs of (s, ττ )) (see (see for instance BouzidiIn et al. [2016], [2016], field A K000(e.g., (e.g., Km =are Q, R). R). The characteristic characteristic function ofin(1) (1) Niculescu [2001], Olgac and Sipahi [2002]). In Bouzidi et al. al. and A , . . . , A constant matrices with entries a m pairs of f (s, τ ) (see for instance Bouzidi et al. [2016], field K (e.g., K = Q, R). The characteristic function of (1) is defined by a quasi-polynomial of the form Niculescu [2001], Olgac and Sipahi [2002]). In Bouzidi et for al. field K (e.g., K = Q, R). The characteristic function of (1) [2016] for example, based on an algebraic formulation Niculescu [2001], Olgac and Sipahi [2002]). In Bouzidi et al. is defined by K a quasi-polynomial of the formfunction of (1) [2016] for example, based on an algebraic formulation for field K (e.g., = Q, R). The characteristic Niculescu [2001], Olgac and Sipahi [2002]). In Bouzidi et al. m is defined by a quasi-polynomial of the form [2016] for example, example, based on an algebraic formulation for is defined by a quasi-polynomial of the form the critical pairs of f (s, τ ), the problem is reduced to [2016] for based on an algebraic formulation for m the critical pairs ofbased f (s, τon ), an thealgebraic problemformulation is reduced for to −j τform s is defined by a quasi-polynomial of ethe m pj (s) [2016] for example, , (2) ff (s, τ ) = m −j τ s the critical pairs of f (s, τ ), the problem is reduced to solving a system of algebraic equations in two variables the critical pairs of f (s, τ ), the problem is reduced to p (s) e , (2) (s, τ ) = solving a system of algebraic equations in two variables j −j τ s m −j τ s the critical pairs of f (s, τ ), the problem is reduced to pjj (s) (s) ee−j τ s ,, (2) and (s, ττ )) = = j=0 p solving a system of algebraic equations in two variables classical symbolic methods are used to obtain certi(2) ff (s, solving a system of algebraic equations in two variables j=0 and classical symbolic methods are used to obtain certi, (2) solving f (s, τ ) = j=0 pj (s) e a system of algebraic equations in two variables and classical symbolic methods are used to obtain certified numerical approximations for the critical pairs. Such j=0 and classical symbolic methods are used to obtain certiwhere in numerical approximations for the critical pairs. certiSuch where the the p pjj ’s ’s are are polynomials polynomials in the the complex complex variable variable ss fied j=0 and classical symbolic methods are used to obtain fied numerical numerical approximations for the critical pairs. Such Such approximations consist in aa pair of intervals isolating the approximations for the critical pairs. where the polynomials in with in K. wherecoefficients the p pjj ’s ’s are are polynomials in the the complex complex variable variable ss fied approximations consist in pair of intervals isolating the with coefficients in K. numerical approximations for the critical pairs. Such wherecoefficients the pj ’s are in the complex variable s fied approximations consist in aathe pair of intervals isolating the two coordinates s and τ of critical pairs. approximations consist in pair of intervals isolating with coefficients in polynomials K. with in K. two coordinates consist s and τinofathe critical pairs. isolating the Such a system is said to be asymptotically stable if all the approximations pair of intervals the with coefficients in K. two coordinates s and τ of the critical pairs. coordinates sthe and τ of the critical pairs. Such a system is said to be asymptotically stable if all the two For the study of asymptotic behavior of critical imagSuch a system is said to be asymptotically stable if all the zeros of the quasi-polynomial f (s, τ ) have negative real two coordinates s and τ of the critical pairs. Such aofsystem is said to be asymptotically if all real the For the study of the asymptotic behavior of critical imagzeros the quasi-polynomial f (s, τ ) havestable negative Such system to be asymptotically stable if all real the inary For the study of asymptotic behavior of imagroots, have been the For studyseveral of the themethods asymptotic behavior of critical criticalfor imagzeros aof of the quasi-polynomial (s, ττclosed haveright negative real parts, i.e., ff (s, ττis)) said = 0 for all ss in half-plane zeros the quasi-polynomial ffthe (s, )) have negative inarythe roots, several methods have been developed developed for the parts, i.e., (s, = 0 for all in the closed right half-plane For the study of the asymptotic behavior of critical imagzeros of the quasi-polynomial f (s, τ ) have negative real inary roots, several methods have been developed for the case of simple imaginary roots (see, e.g., Gu et al. [2003], inaryofroots, several methods the parts, i.e., (s, = 0 0 for for all ssConsidering in the the closed closedthe right half-plane {s ∈ C ≥ 0}. delay ττ as aa case C parts, i.e., ff (s, ττ|| ))(s) = all in right half-plane simple imaginary rootshave (see,been e.g.,developed Gu et al. for [2003], + := inary roots, several methods have been developed for the C := {s ∈ C (s) ≥ 0}. Considering the delay as + parts, i.e., f (s, τ ) = 0 for all s in the closed right half-plane case of simple imaginary roots (see, e.g., Gu et al. [2003], Li et al. [2015], Marshall et al. [1992], Niculescu [2001], case of simple imaginary roots (see, e.g., Gu et al. [2003], := {s ∈ C | (s) ≥ 0}. Considering the delay τ as a C+ parameter, the stability property will obviously depend on Li et al. [2015], Marshall et al. [1992], Niculescu [2001], := {s ∈ C | (s) ≥ 0}. Considering the delay τ as a C + case of simple imaginary roots (see, e.g., Gu et al. [2003], parameter, the stability property will obviously depend on Li et et al. al. [2015], Marshall et al. al. [1992], Niculescu [2001], Olgac and Sipahi [2002] the references therein). For := {s ∈the |stability (s) ≥ property 0}. Considering the delay τ ason a Li C+ value et Niculescu [2001], parameter, will depend the ττCand a question is study Olgac and[2015], SipahiMarshall [2002] and and the[1992], references therein). For parameter, stability property will obviously obviously on Li et al. [2015], Marshall et al. [1992], Niculescu [2001], the value of ofthe and a classical classical question is then then to todepend study the the Olgac and Sipahi [2002] and the references therein). For multiple imaginary roots (MIR), recent methods based on parameter, the stability property will obviously depend on Olgac and Sipahi [2002] and the references therein). For the value of τ and a classical question is then to study the asymptotic stability with respect to the delay value τ . multiple imaginary roots (MIR), recent methods based on the value of τ and a classical question is then to study the Olgac and Sipahi [2002] and the references therein). For asymptotic stability with respect to the delaytovalue τ . the Puiseux multiple imaginary roots (MIR), recent methods based series was developed in Bouzidi et al. [2016], Li the value of τ and a classical question is then study multiple imaginary roots (MIR), recent methods based on on asymptotic stability stability with with respect respect to to the the delay delay value value ττ .. Puiseux series was developed in Bouzidi et al. [2016], Li asymptotic multiple imaginary roots (MIR), recent methods based on A classical approach for studying the asymptotic stability Puiseux series was developed in Bouzidi et al. [2016], Li et al. [2015, 2013]. These methods require in particular the asymptotic stability with respect to the delay value τ . Puiseux series was developed in Bouzidi et al. [2016], Li A classical approach for studying the asymptotic stability et al. [2015, 2013]. These methods require in particular the Puiseux series wasThese developed in Bouzidi etofparticular al. [2016], Li A classical approach for studying the asymptotic stability depending the delay based on the analysis of the soet al. al. [2015, [2015, 2013]. These methods require in particular the symbolic computation of the Puiseux series f (s, τ ) at all A classical on approach for is studying the asymptotic stability et 2013]. methods require in the depending on the delay is based on the analysis of the sosymbolic computation of methods the Puiseux seriesinofparticular f (s, τ ) atthe all A classical approach for the asymptotic stability et al. [2015, 2013]. These require depending on the delay delay isstudying based on the analysis of the socalled critical pairs of f (s, τ ), that is the pairs (ω ∈ R, τ ∈ symbolic computation of the Puiseux series of f (s, τ ) at all the critical pairs. depending on the is based on the analysis of the socomputation of the Puiseux series of f (s, τ ) at all called critical of f (s, τ ), that theanalysis pairs (ωof ∈ the R, τso∈ symbolic the critical pairs. depending on pairs the delay is=based onis the computation of the Puiseux series of f (s, τ ) at all called critical pairs of τff)(s, (s, ), that isfor theinstance pairs (ω (ωGu ∈ R, R,etττal. ∈ symbolic R such that ff (i ω, 0 ;; see the critical pairs. +) called critical pairs of ττ ), that is the pairs ∈ ∈ the critical pairs. R ) such that (i ω, τ ) = 0 see for instance Gu et al. + Contributions and called critical pairs of f (s, τ ), that is the pairs (ω ∈ R, τ ∈ the critical pairs. R ) such that f (i ω, τ ) = 0 ; see for instance Gu et al. [2003], Li et al. [2015], Marshall et al. [1992], Niculescu + ) such that f (i ω, τ ) = 0 ; see for instance Gu et al. Contributions and outline outline of of the the paper. paper. In In this this R+ [2003], Li that et al.f [2015], Marshall etforal.instance [1992], Niculescu Contributions and outline of the paper. In paper, we assume that the critical pairs of f (s, ττ )) this are R ) such (i ω, τ ) = 0 ; see Gu et al. Contributions and outline of the paper. this + [2003], Li et al. [2015], Marshall et al. [1992], Niculescu [2001], Olgac and Sipahi [2002] and the references therein. paper, we assumeand thatoutline the critical pairspaper. of f (s,In are [2003], Li et al. [2015], Marshall etthe al.references [1992], Niculescu Contributions of the In this [2001], Olgac and Sipahi [2002] and therein. paper, we assume that the critical pairs of f (s, τ ) given as certified numerical approximations (using for [2003], Li et al. [2015], Marshall et al. [1992], Niculescu paper, aswecertified assume numerical that the critical pairs of f(using (s, τ ) are are [2001], Olgac andpairs Sipahi [2002] and the references references therein. If such critical exist, the stability is derived from given approximations for [2001], Olgac and Sipahi [2002] and the therein. paper, we assume that the critical pairs of f (s, τ ) are If such critical pairs exist, the stability is derived from given as certified numerical approximations (using for instance the method in Bouzidi et al. [2016]) and we [2001], Olgac and Sipahi [2002] and the references therein. given as certified numerical approximations (using for If such such critical pairs imaginary exist, the the roots stability is derived derived from instance the way these critical behave under small the method in Bouzidi et al. [2016]) and we If critical pairs exist, stability is from given as the certified numerical approximations (using for the way critical these critical imaginary roots behave under small instance method in Bouzidi et al. [2016]) and we focus on the analysis of the asymptotic behavior of If such pairs exist, the stability is derived from instance the method in Bouzidi et al. [2016]) and wess the way these critical imaginary roots behave under small variations of the time-delay τ . As a result, this approach focus on the analysis of the asymptotic behavior of the way these critical imaginary roots behave under small instance the method in Bouzidi et al. [2016]) and wes variations of the time-delay τ . As a result, this approach focus on the analysis of the asymptotic behavior of the way these critical imaginary roots behave under small focus on the analysis of the asymptotic behavior of variations of the time-delay τ . As a result, this approach variations of the time-delay τ . As a result, this approach focus on the analysis of the asymptotic behavior of ss variations of the time-delay τ . As a result, this approach 2405-8963 Copyright © 2019. The Authors. Published by Elsevier Ltd. All rights reserved.
Copyright © 2019 IFAC 172 Copyright 2019 IFAC 172 Control. Peer review© under responsibility of International Federation of Automatic Copyright © 172 Copyright © 2019 2019 IFAC IFAC 172 10.1016/j.ifacol.2019.11.033 Copyright © 2019 IFAC 172
2019 IFAC SSSC Sinaia, Romania, September 9-11, 2019
Yacine Bouzidi et al. / IFAC PapersOnLine 52-17 (2019) 94–98
at these critical pairs. We propose a new approach that avoids any Puiseux series computation at the critical pairs. Instead, the asymptotic behavior is derived through the analysis around the critical pairs of the real branches of a space algebraic curve obtaind by local approximation of the quasi-polynomial f (s, t). This analysis is done while avoiding as much as possible symbolic computations, using interval arithmetic instead. More precisely, given a critical pair (ω0 , τ0 ), our approach proceeds in the following way: (1) We compute an upper bound m for the multiplicity of (ω0 , τ0 ) by evaluating the successive derivatives of f up to the point that we can certify a non-zero evaluation. See Section 3.1. (2) We consider f (s + ω0 , τ + τ0 ) as en element of C[[s, τ ]] by replacing the exponential by their local Taylor series. Note that its coefficients are given with error bounds. To make computations possible, we truncate these series modulo some powers of s and τ , taking at least m + 1 as a truncation bound. This provides a polynomial p(s, τ ) that is an approximation of the studied quasi-polynomial (note that approximations are coming from both numerical errors and the truncation process). See Section 3.2. (3) Finally, we study the local behaviour of the real curve defined by {pr (sr , si , τ ) = pi (sr , si , τ ) = 0} where pr and pi are respectively the real and the imaginary part of p(sr + i si , τ ). More precisely, we isolate the curve in a 3D box that contains only the germ at (0, 0, 0) of the real curve; then, one can conclude about the stability of the system by looking at the intersections of real branches with this box. Organisation of the paper. For a sake of clarity, we briefly recall in Section 2 the method presented in Bouzidi et al. [2016] for the computation of the critical pairs of f (s, t). Then, in Section 3, we show how a quasipolynomial can be locally approximated at a given point by a polynomial while preserving the local topology of the latter. Using this approximated curve, we present in Section 4 a strategy to compute the asymptotic behavior of this curve at a given point by looking at the intersection of the latter with some specific 3D boxes. This enables us to study the stability of the system. 2. CRITICAL PAIRS OF A QUASI-POLYNOMIAL As mentioned in the introduction, to study the asymptotic stability of a system of the form (1), we first need to compute the critical pairs of the corresponding quasipolynomial f (s, τ ). These critical pairs are the intersection points of the curve defined by the quasi-polynomial, f (s, τ ) = 0, with the imaginary axis in s, that is {(ω ∈ R, τ ∈ R+ ) | f (i ω, τ ) = 0} We now briefly recall the method given in Bouzidi et al. [2016] for the certified computation of these critical pairs. The main idea is to use the so-called Rekasius transformation, that allows to recover the infinite number of zeros of f (i ω, τ ) from the real zeros of a zero-dimensional algebraic bivariate system. More precisely, we replace the terms e−τ i ω in the quasipolynomial f (i ω, τ ) by the rational iω fraction 1−γ 1+γ i ω with γ ∈ R. Cleaning the denominators, one obtain a complex polynomial of the form R(ω, γ) + 173
95
i I(ω, γ). Computing the critical pairs then amounts to identifying the real zeros of the following bivariate polynomial system: R(ω, γ) = 0, (3) I(ω, γ) = 0, and from each solution (ω, γ) of (3) the critical delays can then be obtained by: 2 (4) τk = (arctan(ω γ) + k π), k ∈ Z. ω The method proposed in Bouzidi et al. [2016] for computing the real solutions of the above algebraic system consists first in computing a univariate representation of the solutions (5), which is roughly speaking, a one-to-one mapping between the solutions of the system (3) and the root of a univariate polynomial f (t) (see Rouillier [1999] for details), then, numerical approximations for the real solutions are obtained from this representation.
f (t) = 0 ω = gω /g0 γ = gγ /g0
(5)
From the above representation, one can compute intervals with rational endpoints that isolate the coordinates ω and γ of the solutions. Moreover, by substituting the obtained intervals in (4), one also obtain isolating intervals for the corresponding delays τ . Another advantage of the previous representation is that it allows to control the precision of the numerical approximation of the solutions. Indeed, given an isolating interval for a root of the polynomial f (t), one can refine this interval up to a desired precision while it still contains the root which results to refined intervals for ω, γ and τ . As we will see in the sequel, this feature is very useful in the description and the design of our algorithm. 3. LOCAL APPROXIMATION 3.1 An upper bound on the multiplicities of critical pairs An important information that we need for the description of our algorithm is the multiplicities of the critical pairs as roots of the quasi-polynomial f (s, τ ). Computing explicitly such multiplicities is however a challenging problem which is yet not solved (for detail, see Boussaada and Niculescu [2016] and references therein). Therefore, instead of exact multiplicities, we only compute upper bounds on the latter, which will be sufficient for our purpose. Let first recall the definition of the multiplicity of a critical pair of a quasi-polynomial f (s, τ ). Definition 1. A root (ω0 , τ0 ) of the quasi-polynomial f (i ω, τ ) is said to be of multiplicity k if and only if ∂ k−1 f f (ω0 , τ0 ) = ∂f ∂s (ω0 , τ0 ) = · · · = ∂sk−1 (ω0 , τ0 ) = 0 and ∂k f (ω0 , τ0 ) ∂sk
= 0
According to the previous definition, if the evaluation at the point (ω0 , τ0 ) can be achieved in an exact way, a straightforward approach for computing the multiplicity of a critical pair (ω0 , τ0 ) is to evaluate it on the successive derivatives of f (s, τ ) until a non-zero value is obtained. The key point of this paper is to avoid as much as possible to run costly symbolic computations. In particular, we are
2019 IFAC SSSC 96 Sinaia, Romania, September 9-11, 2019
Yacine Bouzidi et al. / IFAC PapersOnLine 52-17 (2019) 94–98
not going to compute the multiplicities of the critical pairs, but only a bound for each of them. To do so, we will just k evaluate the successive derivatives ∂f on the given 2D box ∂k s (using interval arithmetic - see Moore [1979] for details), until we can conclude that the m-th derivative is non zero in this box. In practice, this strategy will often provide the correct multiplicity, but over-bound might happen, so that we cannot assume to know the correct multiplicity. In particular, if we get a bound m = 1 for a critical pair, then this multiplicity is correct and we can conclude that the system is not stable. 3.2 Local approximation Given a critical pair (i ω0 , τ0 ), the quasi-polynomial f (s, τ ) can be considered in K[[s − i ω0 , τ − τ0 ]]. Thus, it can be approximated by a polynomial via a truncation of the powers of (s − i ω0 ) and (τ − τ0 ). In practice, we only have to use the Taylor series of the exponential to convert our quasi-polynomial in series. More precisely, it is well known that ex can be approximated by n−1 xk k!
k=0 n by xn! ex .
with an error bounded We thus can approximate f (s + i ω0 , τ + τ0 ) as m pj (s + i ω0 )e−j i ω0 τ0 e−j i ω0 τ e−j τ0 s e−j s τ j=0
by replacing each last three exponentials by their truncated Taylor series (with x being respectively −j i ω0 τ , −j τ0 s and −j s τ ). This provides us a polynomial p(s, τ ) for which there are two different approximations: • the truncation we just made (for which we have a good majoration of the error), that can be improved by increasing the truncation order, • the approximations we have for τ0 and ω0 (p having coefficients defined by interval), that we can improve on demand from the previous computations. Note that we will always take a truncation bound n greater than the bound m we have for the multiplicity of the critical pair. 4. ASYMPTOTIC BEHAVIOR As shown in the previous section, one can compute a set of disjoint 2D boxes in R2 that isolate all the critical pairs of f (s, τ ). Moreover, thanks to the symbolic univariate representation of these pairs (5), the corresponding boxes can be refined up to an arbitrary precision. Also, for a given critical pair (ω0 , τ0 ) of f (s, τ ), Section 3 provides a local approximation of the quasi-polynomial f (s, τ ) by a two variables polynomial p(s, τ ) that preserves its local topology. Given such a polynomial p(s, τ ), as well as a critical pair of p (a zero of the form (i ω0 , τ0 )), the goal of this section is to show how to analyze the asymptotic behavior of this critical pair. As a consequence, provided that this polynomial results from an appximation of a quasi-polynomial, this will naturally provide the asymptotic behavior of the critical pair as a zero of the quasi-polynomial. 174
4.1 Asymptotic behavior through curve intersection in R3 Let us consider a polynomial p(s, τ ) where s ∈ C and τ ∈ R. In the following, we set s = x + i y, where x, y ∈ R. For convenience, we also denote to the variable τ by z. Substituting s = x+i y as well as τ = z in the polynomial p yields p(x + i y, z) = g(x, y, z) + i h(x, y, z) where g, h are polynomials in three real variables with real coefficients. With this notation, we define a critical pair of p as a point (0, y0 , z0 ) ∈ R3 such that p(i y0 , z0 ) = 0, or equivalently, as the intersection point of the axis x = 0 with curve C defined by g(x, y, z) = 0 C := h(x, y, z) = 0
and whose coordinate in y (resp. z) is y0 (resp. z0 ). In the latter case, and to avoid confusion with critical points of the curve, we refer to it as a zero-crossing point. Moreover, for simplicity, we assume in the sequel that the curve C does not admit a whole component lying in a plan orthogonal to the x-axis. This translates into the fact that for any fixed value x0 of x, the algebraic system {g(x0 , y, z) = h(x0 , y, z) = 0} admits a finite number of complex zeros. Given a zero-crossing point (0, y0 , z0 ), the objective in the following is to analyse the asymptotic behavior of the coordinate x = 0 under a small variation of the coordinate z0 . We suppose that this zero-crossing point (0, y0 , z0 ) is given by a 3D box in R3 that isolates it from other zero-crossing points of the curve. Such a box can be obtained for example by extending the 2D box obtained from the representation (3) on the x-axis. First, we need to introduce the notion of x-critical point of C. A point (x, y, z) of C is called x-critical if it is a critical point of the following projection map Πx :
C →R (x, y, z) → x
Algebraically, a point of the curve C is an x-critical point ∂g ∂h ∂g ∂h 1 if and only if it is a zero of the polynomial ∂y ∂z − ∂z ∂y . In the following, for a sake of briefness, we denote to this polynomial simply by d(x, y, z). As a consequence, we can identify two different kinds of zero-crossing points of C. Definition 2. Let (0, y0 , z0 ) be a zero-crossing point of C. This point is called x-critical if and only if d(0, y0 , z0 ) = 0. Otherwise, it is called regular. Given a zero-crossing point (0, y0 , z0 ) of C, if the latter is a regular zero-crossing points, then the analysis of the asymptotic behavior near to it is straightforward. Indeed, at the neighborhood of these kind of points, the curve C consists in one regular branch that passes through the point (0, y0 , z0 ). As a consequence, the variation of x can be deduced from the sign of d(0, y0 , z0 ). As (0, y0 , z0 ) is not an x-critical point, one can compute a box that contains the latter (in fact, it is simply a 2D box in y and z), and 1
It should be stressed that in our setting and according to this algebraic characterization, a singular point of C is also considered as an x-critical.
2019 IFAC SSSC Sinaia, Romania, September 9-11, 2019
Yacine Bouzidi et al. / IFAC PapersOnLine 52-17 (2019) 94–98
97
refine it until the interval resulting from the evaluation d(B) does not contain zero. For x-critical zero-crossing points, analyzing the asymptotic behavior of C is more involved. In broad words, our approach consists in computing a 3-D box B that isolate this point from other x-critical points of the curve C. Then, this box B is refined until all the intersections of C with the boundary of B are located on the sides that are orthogonal to the x-axis. This guarantees that all the branches that intersect the boundary of the refined box pass through the considered zero-crossing point. The behavior is then easily deduced from these intersections. Definition 3. A box B containing a point p := (0, y0 , z0 ) of C is called isolating for p if and only if all the real branches of C in B pass through the point p. Theorem 1. Let B be a box that contains one x-critical zero-crossing point (0, y0 , z0 ) of C and such that the two following conditions hold • B does not contain any other x-critical point of C • The only intersection points of C with B are located on the sides of B that are orthogonal to the x-axis. Then, B is an isolating box for (0, y0 , z0 ).
Fig. 2. A zero-crossing point and an x-critical point According to the previous theorem, given a box that contains an x-critical zero-crossing point (0, y0 , z0 ) and no other x-critical points of C, it is sufficient to refine the latter until all the intersections with C lie on the sides that are orthogonal to the x-axis, in order to get an isolating box for (0, y0 , z0 ). This is always possible since we have assumed that there is no branch of C lying in a plan orthogonal to the x-axis. 4.2 Computational aspects
Proof. This is a direct consequence of the fact that B contains a single zero-crossing point of C and no other xcritical points. Indeed, since all the branches of C intersect B on the sides that are orthogonal to the x-axis, if one of these branches does not pass through (0, y0 , z0 ), this implies that either it intersects the axis x = 0 at y1 = y0 , or it passes through an x-critical point, which contradicts our assumption. Remark 2. Figures 1 and 2 below show the projection of two non-isolating 3D boxes on the (x − z)-plan. The first situation corresponds to a box B that contains an x-critical zero-crossing point (lying on the branch depicted in blue) as well as another zero-crossing point (lying on the branch depicted in red). The second situation corresponds to a box B that contains two different x-critical points of C among which one is a zero-crossing point.
Let (0, y0 , z0 ) be an x-critical zero-crossing point of C and suppose given a 3D box B that isolate the point (0, y0 , z0 ) from other zero-crossing points of C.
A first computational step consists in refining the box B until (0, y0 , z0 ) is the only x-critical point of C that lies in the refined box. To do so, a straighforward approach consists in computing disjoint boxes that isolate the real zeros of the following system
Sx :=
g(x, y, z) = 0 h(x, y, z) = 0 d(x, y, z) = 0
As for the system in (3), this can be done by computing a univariate representation of the solutions followed by a numerical approximation. This provides a set of disjoint 3D boxes that isolate the real zeros of Sx . Thenceforth, it suffices to refine the box B until it overlaps with only one box of this set to get a refined box that contains only the x-critical point (0, y0 , z0 ). Let denote to this new box by B . The second computational step consists in refining the box B until all the intersections with C lie only on the desired sides. Let B be a box resulting from the refinement of B up to a certain precision. We denote by ymin , ymax (resp. zmin , zmax ), the rational endpoints of its interval on y (resp. z). Then, the box B admits intersections only on the sides that are orthogonal to the x-axis if and only if the following algebraic systems do not have real zeros.
Fig. 1. Two zero-crossing points 175
g(x, ymin , z) = 0 h(x, ymin , z) = 0
g(x, ymax , z) = 0 h(x, ymax , z) = 0
g(x, y, zmin ) = 0 h(x, y, zmin ) = 0
g(x, y, zmax ) = 0 h(x, y, zmax ) = 0
2019 IFAC SSSC 98 Sinaia, Romania, September 9-11, 2019
Yacine Bouzidi et al. / IFAC PapersOnLine 52-17 (2019) 94–98
Consequently, to compute the desired box, our algorithm consists in refining successively the initial box and checking the existence of real solutions for the above algebraic systems. This can be done for example by computing a rational univariate representation (5) for each of these systems and then checking the existence of real roots of univariate polynomials using classical dedicated methods. The algorithm stops when no real solution is found. 5. CONCLUSION We addressed in this paper some aspects related to the study of the stability of differential time-delay systems with commensurate delays. More precisely, following the approach in common use, based on the computation of the critical pairs, we proposed a new strategy for analyzing the asymptotic behavior at these points, that gets rid of Puiseux series computations. Instead, our strategy proceeds by locally approximating at the critical pairs, the quasi-polynomial by a polynomial, and then by analyzing the asymptotic behavior at the zeros of this polynomial seen as specific points of algebraic curve in R3 . The implementation of this approach together with a complete complexity study will be the subject of an ongoing work. In particular, attention will be put on the comparison from the complexity point of view between our approach and the Puiseux-based counterpart. In future work, an interesting problem is to investigate more deeply the impact of numerical errors, due to the truncation as well as the local approximation at the critical pairs, on the certification of our result. REFERENCES Boussaada, I. and Niculescu, S.I. (2016). Tracking the algebraic multiplicity of crossing imaginary roots for generic quasipolynomials: A vandermonde-based approach. IEEE Transactions on Automatic Control, 61(6), 1601–1606. Bouzidi, Y., Poteaux, A., and Quadrat, A. (2016). Computer algebra methods for the stability analysis of differential systems with commensurate time-delays. IFACPapersOnLine, 49(10), 194–199. Gu, K., Kharitonov, V.L., and Chen, J. (2003). Stability of Time-Delay Systems. Birkh¨ auser. Li, X.G., Niculescu, S.I., and C ¸ ela, A. (2015). Analytic Curve Frequency-Sweeping Stability Tests for Systems with Commensurate Delays. Springer. Li, X.G., Niculescu, S.I., C ¸ ela, A., Wang, H.H., and Cai, T.Y. (2013). On computing puiseux series for multiple imaginary characteristic roots of lti systems with commensurate delays. IEEE Transactions on Automatic Control, 58(5), 1338–1343. Marshall, J.E., Gorecki, H., Korytowski, A., and Walton, K. (1992). Time-Delay. Systems: Stability and Performance Criteria with Applications. London: Ellis Horwood. Moore, R.E. (1979). Methods and applications of interval analysis, volume 2. Siam. Niculescu, S.I. (2001). Delay Effects on Stability: A Robust Control Approach. Number 269 in Lecture Notes in Control and Information Sciences. Springer. Olgac, N. and Sipahi, R. (2002). An exact method for the stability analysis of time-delayed linear time-invariant 176
(LTI) systems. IEEE Transactions on Automatic Control, 47(5), 793–796. Rouillier, F. (1999). Solving zero-dimensional systems through the rational univariate representation. J. of Applicable Algebra in Engineering, Communication and Computing, 9(5), 433–461.
ScienceDirect
IFAC PapersOnLine 52-17 (2019) 94–98
Real Real Real Real
Curve Analysis and Stability Curve Analysis and Stability Curve Analysis and Stability Time-delay Systems Curve Analysis and Stability Time-delay Systems Time-delay Systems Time-delay Systems ∗∗ Yacine Bouzidi ∗∗ Adrien Poteaux ∗∗
of of of of
Yacine Bouzidi ∗∗ Adrien Poteaux ∗∗ Yacine Bouzidi Bouzidi ∗ Adrien Adrien Poteaux Poteaux ∗∗ Yacine Yacine Bouzidi Adrien Poteaux ∗∗ ∗ ∗ Inria Lille - Nord Europe, Parc Scientifique de la Haute Borne, 40 Inria Lille Nord Europe, Parc Scientifique de la Haute Borne, 40 ∗ ∗ Inria Lille -- Nord Europe, Parc Scientifique de Borne, 40 Avenue Halley, Bat A -- Park Plaza, 59650 Villeneuve d’Ascq, France Nord Europe, Parc Scientifique de la la Haute Haute Borne, 40 ∗ Inria Lille Avenue Halley, Bat A Park Plaza, 59650 Villeneuve d’Ascq, France Inria Lille Nord Europe, Parc Scientifique de la Haute Borne, 40 Avenue Halley, Bat A Park Plaza, 59650 Villeneuve d’Ascq, France (e-mail: [email protected]). Avenue Halley, Bat A - Park Plaza, 59650 Villeneuve d’Ascq, France (e-mail: [email protected]). ∗∗ Avenue Halley, Bat A - INRIA, Park Plaza, 59650 Villeneuve d’Ascq, France (e-mail: [email protected]). UMR 9189 CRIStAL Centre de ∗∗ Univ. Lille, CNRS, (e-mail: [email protected]). INRIA, UMR 9189 - CRIStAL - Centre de ∗∗ Univ. Lille, CNRS, (e-mail: [email protected]). ∗∗ Univ. Lille, CNRS, INRIA, UMR 9189 -- CRIStAL CRIStAL Centre de Recherche en Informatique Signal et de F-59000 CNRS, INRIA, UMR 9189 -- Centre de ∗∗ Univ. Lille, Recherche en Informatique Signal et Automatique Automatique de Lille, Lille, F-59000 Univ.Lille, Lille, CNRS,(e-mail: INRIA, UMR 9189 - CRIStAL - Centre de Recherche en Informatique Signal et Automatique de Lille, F-59000 France [email protected]). Recherche en Informatique Signal et Automatique de Lille, F-59000 Lille, France (e-mail: [email protected]). Recherche en Informatique Signal et Automatique de Lille, F-59000 Lille, France (e-mail: [email protected]). Lille, France (e-mail: [email protected]). Lille, France (e-mail: [email protected]). Abstract: In this work, we investigate investigate the the asymptotic asymptotic stability stability of of LTI LTI differential differential commensurate commensurate Abstract: In this work, we Abstract: In this work, we investigate the asymptotic stability of LTI differential commensurate time-delay systems whose dynamics are defined in the frequency domain by quasi-polynomials of Abstract: systems In this work, we investigate the asymptotic stability ofdomain LTI differential commensurate time-delay whose dynamics defined in the frequency by quasi-polynomials of m we −j τ sare Abstract: Inτthis work, investigate the asymptotic stability ofdomain LTI differential commensurate time-delay systems whose dynamics are defined in the frequency domain by quasi-polynomials of m the form f (s, ) = p (s) e . We propose a new approach for studying the stability with −j τ s j time-delay systems whose dynamics are defined in the frequency by quasi-polynomials of j=0 the form f (s, τ) = (s) e−j . Wedefined propose new approachdomain for studying the stability with m τs sare m j=0 p time-delay systems whose in aaathe frequency bythe quasi-polynomials of τ the formtoff (s, (s, = pjjjdynamics (s) ee−j . We propose new approach for studying the stability stability with with respect the we determine asymptotic behavior of the roots of quasi-polynomial the form ττdelay: )) = p (s) propose new approach for studying the j=0 m −j τ the s . We j=0 respect to the delay: we determine the asymptotic behavior of the roots of the quasi-polynomial the form f (s, τ ) = p (s) e . We propose a new approach for studying the stability with j j=0 respect to the the delay:pairs we determine determine theanalysing asymptotic behavior of the the roots of the the quasi-polynomial quasi-polynomial near to critical of the intersection of an algebraic respect to delay: we the asymptotic of of near to the the critical of ff (s, (s, ττ )) by by thebehavior intersection of roots an associated associated algebraic space space respect to the delay:pairs werefined determine theanalysing asymptotic behavior ofthe the roots ofmethods, the quasi-polynomial near to the critical pairs of ff (s, ττ )) by analysing the intersection of an associated algebraic space curve with sufficiently 3D real boxes. Compared to existing our near to the critical pairs of (s, by analysing the intersection of an associated algebraic space curvetowith sufficiently refined 3D realanalysing boxes. Compared to theofexisting methods, our method method near the critical pairs of f (s, τ ) by the intersection an associated algebraic space curve with sufficiently refined 3D real boxes. Compared to the existing methods, our method does not require any Puiseux series computation. curve with sufficiently refinedseries 3D real boxes. Compared to the existing methods, our method does not require any Puiseux computation. curve with sufficiently refined 3D real boxes. Compared to the existing methods, our method does not not require require any any Puiseux Puiseux series series computation. computation. does Copyright 2019. The Authors. Published Elsevier Ltd. All rights reserved. does not © require any Puiseux seriesby computation. Keywords: Time delay systems, curve Keywords: Time delay systems, curve analysis, analysis, stability, stability, certified certified tests, tests, symbolic-numeric symbolic-numeric Keywords: Time delay systems, curve analysis, stability, certified tests, methods. Keywords: Time delay systems, curve analysis, stability, certified tests, symbolic-numeric symbolic-numeric methods. Keywords: Time delay systems, curve analysis, stability, certified tests, symbolic-numeric methods. methods. methods. 1. INTRODUCTION INTRODUCTION yields 1. yields aa decomposition decomposition of of the the ττ -axis -axis into into aa set set of of disjoint disjoint 1. yields aa decomposition decomposition of the the ττ -axis -axis intoisaaeither set of ofstable disjoint intervals in which the considered system or 1. INTRODUCTION INTRODUCTION yields of into set disjoint intervals in which the considered system is either stable or 1.investigate INTRODUCTION yields a decomposition of the τ -axis intois aeither set ofstable disjoint intervals in which the considered system or unstable. In this paper, we the asymptotic stability of intervals in which the considered system is either stable or In this paper, we investigate the asymptotic stability of unstable. intervals in which the considered system is either stable or unstable. In this paper, we investigate the asymptotic stability of LTI differential commensurate time-delay systems whose unstable. In this paper, we investigate the asymptotic stability of Algorithmically, the above approach can be decomposed LTIthis differential commensurate time-delay systems whose unstable. In paper, we investigate the asymptotic stability of the above approach can be decomposed LTI commensurate time-delay systems dynamics are defined defined by the the following following state-space repre- Algorithmically, LTI differential differential commensurate time-delay systems whose whose Algorithmically, the above can be in First, the of dynamics are by state-space repreAlgorithmically, the problems. above approach approach be decomposed decomposed LTI differential commensurate time-delay systems whose in two two independent independent problems. First, can the computation computation of dynamics are defined by the following state-space representation dynamics are defined by the following state-space repreAlgorithmically, the above approach can be decomposed in two independent problems. First, the computation of the critical pairs of f (s, τ ), and second, for each critical sentation m in two independent the for computation of dynamics the critical pairs of fproblems. (s, τ ), andFirst, second, each critical mthe following state-space representation sentation are defined by in two independent problems. First, the computation of the critical pairs of f (s, τ ), and second, for each critical pair (s , τ ), the study of the asymptotic behavior of the m 0 0 the critical pairs of f (s, τ ), and second, for each critical x(t) ˙ = A x(t − k τ ) (1) m sentation k pair (s , τ ), the study of the asymptotic behavior of the 0 0 x(t) ˙ = (1) the critical pairs of f (s, τ ), and second, for each critical m Ak x(t − k τ ) pair (s , τ ), the study of the asymptotic behavior of the root ss000 with to small variation the τthe 0 .. , τ00 ), respect the study the asymptotic x(t) ˙˙ = (1) x(t) = k=0 A Akk x(t x(t − −k k ττ )) (1) pair root (s with to aaof small variation of ofbehavior the delay delayof 0 pair (s , τ0 ), respect the study the asymptotic behavior ofτττthe x(t) ˙{τ ∈=Rk=0 A≥ −mk ∈ τ )N := {0, 1, . (1) root ss0000 with respect to aaofsmall variation of the delay k x(t 0 .. root with respect to small variation of the delay k=0 0 := | τ 0}, . .} where τ ∈ R k=0 + + There exist numerous methods for computing the critical root s with respect to a small variation of the delay τ | τ ≥ 0}, m ∈ N+ := {0, 1, . . .} There 0exist numerous methods for computing the critical where τ ∈ R+ := {τ ∈ Rk=0 0. :=are {τ ∈ ∈ R || ττ ≥ ≥matrices 0}, m m∈ ∈ with N+ :=entries {0, 1, 1, in . .. .} .} where τ, .∈ ∈. .R R, A and A constant aa There + 0τ m exist numerous methods for computing the critical := {τ R 0}, N := {0, . where pairs of f (s, τ ) (see for instance Bouzidi et al. [2016], + + There exist numerous methods for computing the critical and A , . . . , A are constant matrices with entries in 0τ ∈ R+ m pairs of f (s,numerous τ ) (see for instance Bouzidi et the al. critical [2016], :=are {τ ∈ RThe | τ ≥matrices 0}, m ∈ with N+function :=entries {0, 1,ofin . . .} where There exist methods for[2002]). computing and A A are constant matrices with entries in field K K Q, pairs of ff[2001], (s, for instance Bouzidi et al. and ,, .. .. .. ,, A constant aa Niculescu Olgac and Sipahi Bouzidi et m= pairs of (s, ττ )) (see (see for instance BouzidiIn et al. [2016], [2016], field A K000(e.g., (e.g., Km =are Q, R). R). The characteristic characteristic function ofin(1) (1) Niculescu [2001], Olgac and Sipahi [2002]). In Bouzidi et al. al. and A , . . . , A constant matrices with entries a m pairs of f (s, τ ) (see for instance Bouzidi et al. [2016], field K (e.g., K = Q, R). The characteristic function of (1) is defined by a quasi-polynomial of the form Niculescu [2001], Olgac and Sipahi [2002]). In Bouzidi et for al. field K (e.g., K = Q, R). The characteristic function of (1) [2016] for example, based on an algebraic formulation Niculescu [2001], Olgac and Sipahi [2002]). In Bouzidi et al. is defined by K a quasi-polynomial of the formfunction of (1) [2016] for example, based on an algebraic formulation for field K (e.g., = Q, R). The characteristic Niculescu [2001], Olgac and Sipahi [2002]). In Bouzidi et al. m is defined by a quasi-polynomial of the form [2016] for example, example, based on an algebraic formulation for is defined by a quasi-polynomial of the form the critical pairs of f (s, τ ), the problem is reduced to [2016] for based on an algebraic formulation for m the critical pairs ofbased f (s, τon ), an thealgebraic problemformulation is reduced for to −j τform s is defined by a quasi-polynomial of ethe m pj (s) [2016] for example, , (2) ff (s, τ ) = m −j τ s the critical pairs of f (s, τ ), the problem is reduced to solving a system of algebraic equations in two variables the critical pairs of f (s, τ ), the problem is reduced to p (s) e , (2) (s, τ ) = solving a system of algebraic equations in two variables j −j τ s m −j τ s the critical pairs of f (s, τ ), the problem is reduced to pjj (s) (s) ee−j τ s ,, (2) and (s, ττ )) = = j=0 p solving a system of algebraic equations in two variables classical symbolic methods are used to obtain certi(2) ff (s, solving a system of algebraic equations in two variables j=0 and classical symbolic methods are used to obtain certi, (2) solving f (s, τ ) = j=0 pj (s) e a system of algebraic equations in two variables and classical symbolic methods are used to obtain certified numerical approximations for the critical pairs. Such j=0 and classical symbolic methods are used to obtain certiwhere in numerical approximations for the critical pairs. certiSuch where the the p pjj ’s ’s are are polynomials polynomials in the the complex complex variable variable ss fied j=0 and classical symbolic methods are used to obtain fied numerical numerical approximations for the critical pairs. Such Such approximations consist in aa pair of intervals isolating the approximations for the critical pairs. where the polynomials in with in K. wherecoefficients the p pjj ’s ’s are are polynomials in the the complex complex variable variable ss fied approximations consist in pair of intervals isolating the with coefficients in K. numerical approximations for the critical pairs. Such wherecoefficients the pj ’s are in the complex variable s fied approximations consist in aathe pair of intervals isolating the two coordinates s and τ of critical pairs. approximations consist in pair of intervals isolating with coefficients in polynomials K. with in K. two coordinates consist s and τinofathe critical pairs. isolating the Such a system is said to be asymptotically stable if all the approximations pair of intervals the with coefficients in K. two coordinates s and τ of the critical pairs. coordinates sthe and τ of the critical pairs. Such a system is said to be asymptotically stable if all the two For the study of asymptotic behavior of critical imagSuch a system is said to be asymptotically stable if all the zeros of the quasi-polynomial f (s, τ ) have negative real two coordinates s and τ of the critical pairs. Such aofsystem is said to be asymptotically if all real the For the study of the asymptotic behavior of critical imagzeros the quasi-polynomial f (s, τ ) havestable negative Such system to be asymptotically stable if all real the inary For the study of asymptotic behavior of imagroots, have been the For studyseveral of the themethods asymptotic behavior of critical criticalfor imagzeros aof of the quasi-polynomial (s, ττclosed haveright negative real parts, i.e., ff (s, ττis)) said = 0 for all ss in half-plane zeros the quasi-polynomial ffthe (s, )) have negative inarythe roots, several methods have been developed developed for the parts, i.e., (s, = 0 for all in the closed right half-plane For the study of the asymptotic behavior of critical imagzeros of the quasi-polynomial f (s, τ ) have negative real inary roots, several methods have been developed for the case of simple imaginary roots (see, e.g., Gu et al. [2003], inaryofroots, several methods the parts, i.e., (s, = 0 0 for for all ssConsidering in the the closed closedthe right half-plane {s ∈ C ≥ 0}. delay ττ as aa case C parts, i.e., ff (s, ττ|| ))(s) = all in right half-plane simple imaginary rootshave (see,been e.g.,developed Gu et al. for [2003], + := inary roots, several methods have been developed for the C := {s ∈ C (s) ≥ 0}. Considering the delay as + parts, i.e., f (s, τ ) = 0 for all s in the closed right half-plane case of simple imaginary roots (see, e.g., Gu et al. [2003], Li et al. [2015], Marshall et al. [1992], Niculescu [2001], case of simple imaginary roots (see, e.g., Gu et al. [2003], := {s ∈ C | (s) ≥ 0}. Considering the delay τ as a C+ parameter, the stability property will obviously depend on Li et al. [2015], Marshall et al. [1992], Niculescu [2001], := {s ∈ C | (s) ≥ 0}. Considering the delay τ as a C + case of simple imaginary roots (see, e.g., Gu et al. [2003], parameter, the stability property will obviously depend on Li et et al. al. [2015], Marshall et al. al. [1992], Niculescu [2001], Olgac and Sipahi [2002] the references therein). For := {s ∈the |stability (s) ≥ property 0}. Considering the delay τ ason a Li C+ value et Niculescu [2001], parameter, will depend the ττCand a question is study Olgac and[2015], SipahiMarshall [2002] and and the[1992], references therein). For parameter, stability property will obviously obviously on Li et al. [2015], Marshall et al. [1992], Niculescu [2001], the value of ofthe and a classical classical question is then then to todepend study the the Olgac and Sipahi [2002] and the references therein). For multiple imaginary roots (MIR), recent methods based on parameter, the stability property will obviously depend on Olgac and Sipahi [2002] and the references therein). For the value of τ and a classical question is then to study the asymptotic stability with respect to the delay value τ . multiple imaginary roots (MIR), recent methods based on the value of τ and a classical question is then to study the Olgac and Sipahi [2002] and the references therein). For asymptotic stability with respect to the delaytovalue τ . the Puiseux multiple imaginary roots (MIR), recent methods based series was developed in Bouzidi et al. [2016], Li the value of τ and a classical question is then study multiple imaginary roots (MIR), recent methods based on on asymptotic stability stability with with respect respect to to the the delay delay value value ττ .. Puiseux series was developed in Bouzidi et al. [2016], Li asymptotic multiple imaginary roots (MIR), recent methods based on A classical approach for studying the asymptotic stability Puiseux series was developed in Bouzidi et al. [2016], Li et al. [2015, 2013]. These methods require in particular the asymptotic stability with respect to the delay value τ . Puiseux series was developed in Bouzidi et al. [2016], Li A classical approach for studying the asymptotic stability et al. [2015, 2013]. These methods require in particular the Puiseux series wasThese developed in Bouzidi etofparticular al. [2016], Li A classical approach for studying the asymptotic stability depending the delay based on the analysis of the soet al. al. [2015, [2015, 2013]. These methods require in particular the symbolic computation of the Puiseux series f (s, τ ) at all A classical on approach for is studying the asymptotic stability et 2013]. methods require in the depending on the delay is based on the analysis of the sosymbolic computation of methods the Puiseux seriesinofparticular f (s, τ ) atthe all A classical approach for the asymptotic stability et al. [2015, 2013]. These require depending on the delay delay isstudying based on the analysis of the socalled critical pairs of f (s, τ ), that is the pairs (ω ∈ R, τ ∈ symbolic computation of the Puiseux series of f (s, τ ) at all the critical pairs. depending on the is based on the analysis of the socomputation of the Puiseux series of f (s, τ ) at all called critical of f (s, τ ), that theanalysis pairs (ωof ∈ the R, τso∈ symbolic the critical pairs. depending on pairs the delay is=based onis the computation of the Puiseux series of f (s, τ ) at all called critical pairs of τff)(s, (s, ), that isfor theinstance pairs (ω (ωGu ∈ R, R,etττal. ∈ symbolic R such that ff (i ω, 0 ;; see the critical pairs. +) called critical pairs of ττ ), that is the pairs ∈ ∈ the critical pairs. R ) such that (i ω, τ ) = 0 see for instance Gu et al. + Contributions and called critical pairs of f (s, τ ), that is the pairs (ω ∈ R, τ ∈ the critical pairs. R ) such that f (i ω, τ ) = 0 ; see for instance Gu et al. [2003], Li et al. [2015], Marshall et al. [1992], Niculescu + ) such that f (i ω, τ ) = 0 ; see for instance Gu et al. Contributions and outline outline of of the the paper. paper. In In this this R+ [2003], Li that et al.f [2015], Marshall etforal.instance [1992], Niculescu Contributions and outline of the paper. In paper, we assume that the critical pairs of f (s, ττ )) this are R ) such (i ω, τ ) = 0 ; see Gu et al. Contributions and outline of the paper. this + [2003], Li et al. [2015], Marshall et al. [1992], Niculescu [2001], Olgac and Sipahi [2002] and the references therein. paper, we assumeand thatoutline the critical pairspaper. of f (s,In are [2003], Li et al. [2015], Marshall etthe al.references [1992], Niculescu Contributions of the In this [2001], Olgac and Sipahi [2002] and therein. paper, we assume that the critical pairs of f (s, τ ) given as certified numerical approximations (using for [2003], Li et al. [2015], Marshall et al. [1992], Niculescu paper, aswecertified assume numerical that the critical pairs of f(using (s, τ ) are are [2001], Olgac andpairs Sipahi [2002] and the references references therein. If such critical exist, the stability is derived from given approximations for [2001], Olgac and Sipahi [2002] and the therein. paper, we assume that the critical pairs of f (s, τ ) are If such critical pairs exist, the stability is derived from given as certified numerical approximations (using for instance the method in Bouzidi et al. [2016]) and we [2001], Olgac and Sipahi [2002] and the references therein. given as certified numerical approximations (using for If such such critical pairs imaginary exist, the the roots stability is derived derived from instance the way these critical behave under small the method in Bouzidi et al. [2016]) and we If critical pairs exist, stability is from given as the certified numerical approximations (using for the way critical these critical imaginary roots behave under small instance method in Bouzidi et al. [2016]) and we focus on the analysis of the asymptotic behavior of If such pairs exist, the stability is derived from instance the method in Bouzidi et al. [2016]) and wess the way these critical imaginary roots behave under small variations of the time-delay τ . As a result, this approach focus on the analysis of the asymptotic behavior of the way these critical imaginary roots behave under small instance the method in Bouzidi et al. [2016]) and wes variations of the time-delay τ . As a result, this approach focus on the analysis of the asymptotic behavior of the way these critical imaginary roots behave under small focus on the analysis of the asymptotic behavior of variations of the time-delay τ . As a result, this approach variations of the time-delay τ . As a result, this approach focus on the analysis of the asymptotic behavior of ss variations of the time-delay τ . As a result, this approach 2405-8963 Copyright © 2019. The Authors. Published by Elsevier Ltd. All rights reserved.
Copyright © 2019 IFAC 172 Copyright 2019 IFAC 172 Control. Peer review© under responsibility of International Federation of Automatic Copyright © 172 Copyright © 2019 2019 IFAC IFAC 172 10.1016/j.ifacol.2019.11.033 Copyright © 2019 IFAC 172
2019 IFAC SSSC Sinaia, Romania, September 9-11, 2019
Yacine Bouzidi et al. / IFAC PapersOnLine 52-17 (2019) 94–98
at these critical pairs. We propose a new approach that avoids any Puiseux series computation at the critical pairs. Instead, the asymptotic behavior is derived through the analysis around the critical pairs of the real branches of a space algebraic curve obtaind by local approximation of the quasi-polynomial f (s, t). This analysis is done while avoiding as much as possible symbolic computations, using interval arithmetic instead. More precisely, given a critical pair (ω0 , τ0 ), our approach proceeds in the following way: (1) We compute an upper bound m for the multiplicity of (ω0 , τ0 ) by evaluating the successive derivatives of f up to the point that we can certify a non-zero evaluation. See Section 3.1. (2) We consider f (s + ω0 , τ + τ0 ) as en element of C[[s, τ ]] by replacing the exponential by their local Taylor series. Note that its coefficients are given with error bounds. To make computations possible, we truncate these series modulo some powers of s and τ , taking at least m + 1 as a truncation bound. This provides a polynomial p(s, τ ) that is an approximation of the studied quasi-polynomial (note that approximations are coming from both numerical errors and the truncation process). See Section 3.2. (3) Finally, we study the local behaviour of the real curve defined by {pr (sr , si , τ ) = pi (sr , si , τ ) = 0} where pr and pi are respectively the real and the imaginary part of p(sr + i si , τ ). More precisely, we isolate the curve in a 3D box that contains only the germ at (0, 0, 0) of the real curve; then, one can conclude about the stability of the system by looking at the intersections of real branches with this box. Organisation of the paper. For a sake of clarity, we briefly recall in Section 2 the method presented in Bouzidi et al. [2016] for the computation of the critical pairs of f (s, t). Then, in Section 3, we show how a quasipolynomial can be locally approximated at a given point by a polynomial while preserving the local topology of the latter. Using this approximated curve, we present in Section 4 a strategy to compute the asymptotic behavior of this curve at a given point by looking at the intersection of the latter with some specific 3D boxes. This enables us to study the stability of the system. 2. CRITICAL PAIRS OF A QUASI-POLYNOMIAL As mentioned in the introduction, to study the asymptotic stability of a system of the form (1), we first need to compute the critical pairs of the corresponding quasipolynomial f (s, τ ). These critical pairs are the intersection points of the curve defined by the quasi-polynomial, f (s, τ ) = 0, with the imaginary axis in s, that is {(ω ∈ R, τ ∈ R+ ) | f (i ω, τ ) = 0} We now briefly recall the method given in Bouzidi et al. [2016] for the certified computation of these critical pairs. The main idea is to use the so-called Rekasius transformation, that allows to recover the infinite number of zeros of f (i ω, τ ) from the real zeros of a zero-dimensional algebraic bivariate system. More precisely, we replace the terms e−τ i ω in the quasipolynomial f (i ω, τ ) by the rational iω fraction 1−γ 1+γ i ω with γ ∈ R. Cleaning the denominators, one obtain a complex polynomial of the form R(ω, γ) + 173
95
i I(ω, γ). Computing the critical pairs then amounts to identifying the real zeros of the following bivariate polynomial system: R(ω, γ) = 0, (3) I(ω, γ) = 0, and from each solution (ω, γ) of (3) the critical delays can then be obtained by: 2 (4) τk = (arctan(ω γ) + k π), k ∈ Z. ω The method proposed in Bouzidi et al. [2016] for computing the real solutions of the above algebraic system consists first in computing a univariate representation of the solutions (5), which is roughly speaking, a one-to-one mapping between the solutions of the system (3) and the root of a univariate polynomial f (t) (see Rouillier [1999] for details), then, numerical approximations for the real solutions are obtained from this representation.
f (t) = 0 ω = gω /g0 γ = gγ /g0
(5)
From the above representation, one can compute intervals with rational endpoints that isolate the coordinates ω and γ of the solutions. Moreover, by substituting the obtained intervals in (4), one also obtain isolating intervals for the corresponding delays τ . Another advantage of the previous representation is that it allows to control the precision of the numerical approximation of the solutions. Indeed, given an isolating interval for a root of the polynomial f (t), one can refine this interval up to a desired precision while it still contains the root which results to refined intervals for ω, γ and τ . As we will see in the sequel, this feature is very useful in the description and the design of our algorithm. 3. LOCAL APPROXIMATION 3.1 An upper bound on the multiplicities of critical pairs An important information that we need for the description of our algorithm is the multiplicities of the critical pairs as roots of the quasi-polynomial f (s, τ ). Computing explicitly such multiplicities is however a challenging problem which is yet not solved (for detail, see Boussaada and Niculescu [2016] and references therein). Therefore, instead of exact multiplicities, we only compute upper bounds on the latter, which will be sufficient for our purpose. Let first recall the definition of the multiplicity of a critical pair of a quasi-polynomial f (s, τ ). Definition 1. A root (ω0 , τ0 ) of the quasi-polynomial f (i ω, τ ) is said to be of multiplicity k if and only if ∂ k−1 f f (ω0 , τ0 ) = ∂f ∂s (ω0 , τ0 ) = · · · = ∂sk−1 (ω0 , τ0 ) = 0 and ∂k f (ω0 , τ0 ) ∂sk
= 0
According to the previous definition, if the evaluation at the point (ω0 , τ0 ) can be achieved in an exact way, a straightforward approach for computing the multiplicity of a critical pair (ω0 , τ0 ) is to evaluate it on the successive derivatives of f (s, τ ) until a non-zero value is obtained. The key point of this paper is to avoid as much as possible to run costly symbolic computations. In particular, we are
2019 IFAC SSSC 96 Sinaia, Romania, September 9-11, 2019
Yacine Bouzidi et al. / IFAC PapersOnLine 52-17 (2019) 94–98
not going to compute the multiplicities of the critical pairs, but only a bound for each of them. To do so, we will just k evaluate the successive derivatives ∂f on the given 2D box ∂k s (using interval arithmetic - see Moore [1979] for details), until we can conclude that the m-th derivative is non zero in this box. In practice, this strategy will often provide the correct multiplicity, but over-bound might happen, so that we cannot assume to know the correct multiplicity. In particular, if we get a bound m = 1 for a critical pair, then this multiplicity is correct and we can conclude that the system is not stable. 3.2 Local approximation Given a critical pair (i ω0 , τ0 ), the quasi-polynomial f (s, τ ) can be considered in K[[s − i ω0 , τ − τ0 ]]. Thus, it can be approximated by a polynomial via a truncation of the powers of (s − i ω0 ) and (τ − τ0 ). In practice, we only have to use the Taylor series of the exponential to convert our quasi-polynomial in series. More precisely, it is well known that ex can be approximated by n−1 xk k!
k=0 n by xn! ex .
with an error bounded We thus can approximate f (s + i ω0 , τ + τ0 ) as m pj (s + i ω0 )e−j i ω0 τ0 e−j i ω0 τ e−j τ0 s e−j s τ j=0
by replacing each last three exponentials by their truncated Taylor series (with x being respectively −j i ω0 τ , −j τ0 s and −j s τ ). This provides us a polynomial p(s, τ ) for which there are two different approximations: • the truncation we just made (for which we have a good majoration of the error), that can be improved by increasing the truncation order, • the approximations we have for τ0 and ω0 (p having coefficients defined by interval), that we can improve on demand from the previous computations. Note that we will always take a truncation bound n greater than the bound m we have for the multiplicity of the critical pair. 4. ASYMPTOTIC BEHAVIOR As shown in the previous section, one can compute a set of disjoint 2D boxes in R2 that isolate all the critical pairs of f (s, τ ). Moreover, thanks to the symbolic univariate representation of these pairs (5), the corresponding boxes can be refined up to an arbitrary precision. Also, for a given critical pair (ω0 , τ0 ) of f (s, τ ), Section 3 provides a local approximation of the quasi-polynomial f (s, τ ) by a two variables polynomial p(s, τ ) that preserves its local topology. Given such a polynomial p(s, τ ), as well as a critical pair of p (a zero of the form (i ω0 , τ0 )), the goal of this section is to show how to analyze the asymptotic behavior of this critical pair. As a consequence, provided that this polynomial results from an appximation of a quasi-polynomial, this will naturally provide the asymptotic behavior of the critical pair as a zero of the quasi-polynomial. 174
4.1 Asymptotic behavior through curve intersection in R3 Let us consider a polynomial p(s, τ ) where s ∈ C and τ ∈ R. In the following, we set s = x + i y, where x, y ∈ R. For convenience, we also denote to the variable τ by z. Substituting s = x+i y as well as τ = z in the polynomial p yields p(x + i y, z) = g(x, y, z) + i h(x, y, z) where g, h are polynomials in three real variables with real coefficients. With this notation, we define a critical pair of p as a point (0, y0 , z0 ) ∈ R3 such that p(i y0 , z0 ) = 0, or equivalently, as the intersection point of the axis x = 0 with curve C defined by g(x, y, z) = 0 C := h(x, y, z) = 0
and whose coordinate in y (resp. z) is y0 (resp. z0 ). In the latter case, and to avoid confusion with critical points of the curve, we refer to it as a zero-crossing point. Moreover, for simplicity, we assume in the sequel that the curve C does not admit a whole component lying in a plan orthogonal to the x-axis. This translates into the fact that for any fixed value x0 of x, the algebraic system {g(x0 , y, z) = h(x0 , y, z) = 0} admits a finite number of complex zeros. Given a zero-crossing point (0, y0 , z0 ), the objective in the following is to analyse the asymptotic behavior of the coordinate x = 0 under a small variation of the coordinate z0 . We suppose that this zero-crossing point (0, y0 , z0 ) is given by a 3D box in R3 that isolates it from other zero-crossing points of the curve. Such a box can be obtained for example by extending the 2D box obtained from the representation (3) on the x-axis. First, we need to introduce the notion of x-critical point of C. A point (x, y, z) of C is called x-critical if it is a critical point of the following projection map Πx :
C →R (x, y, z) → x
Algebraically, a point of the curve C is an x-critical point ∂g ∂h ∂g ∂h 1 if and only if it is a zero of the polynomial ∂y ∂z − ∂z ∂y . In the following, for a sake of briefness, we denote to this polynomial simply by d(x, y, z). As a consequence, we can identify two different kinds of zero-crossing points of C. Definition 2. Let (0, y0 , z0 ) be a zero-crossing point of C. This point is called x-critical if and only if d(0, y0 , z0 ) = 0. Otherwise, it is called regular. Given a zero-crossing point (0, y0 , z0 ) of C, if the latter is a regular zero-crossing points, then the analysis of the asymptotic behavior near to it is straightforward. Indeed, at the neighborhood of these kind of points, the curve C consists in one regular branch that passes through the point (0, y0 , z0 ). As a consequence, the variation of x can be deduced from the sign of d(0, y0 , z0 ). As (0, y0 , z0 ) is not an x-critical point, one can compute a box that contains the latter (in fact, it is simply a 2D box in y and z), and 1
It should be stressed that in our setting and according to this algebraic characterization, a singular point of C is also considered as an x-critical.
2019 IFAC SSSC Sinaia, Romania, September 9-11, 2019
Yacine Bouzidi et al. / IFAC PapersOnLine 52-17 (2019) 94–98
97
refine it until the interval resulting from the evaluation d(B) does not contain zero. For x-critical zero-crossing points, analyzing the asymptotic behavior of C is more involved. In broad words, our approach consists in computing a 3-D box B that isolate this point from other x-critical points of the curve C. Then, this box B is refined until all the intersections of C with the boundary of B are located on the sides that are orthogonal to the x-axis. This guarantees that all the branches that intersect the boundary of the refined box pass through the considered zero-crossing point. The behavior is then easily deduced from these intersections. Definition 3. A box B containing a point p := (0, y0 , z0 ) of C is called isolating for p if and only if all the real branches of C in B pass through the point p. Theorem 1. Let B be a box that contains one x-critical zero-crossing point (0, y0 , z0 ) of C and such that the two following conditions hold • B does not contain any other x-critical point of C • The only intersection points of C with B are located on the sides of B that are orthogonal to the x-axis. Then, B is an isolating box for (0, y0 , z0 ).
Fig. 2. A zero-crossing point and an x-critical point According to the previous theorem, given a box that contains an x-critical zero-crossing point (0, y0 , z0 ) and no other x-critical points of C, it is sufficient to refine the latter until all the intersections with C lie on the sides that are orthogonal to the x-axis, in order to get an isolating box for (0, y0 , z0 ). This is always possible since we have assumed that there is no branch of C lying in a plan orthogonal to the x-axis. 4.2 Computational aspects
Proof. This is a direct consequence of the fact that B contains a single zero-crossing point of C and no other xcritical points. Indeed, since all the branches of C intersect B on the sides that are orthogonal to the x-axis, if one of these branches does not pass through (0, y0 , z0 ), this implies that either it intersects the axis x = 0 at y1 = y0 , or it passes through an x-critical point, which contradicts our assumption. Remark 2. Figures 1 and 2 below show the projection of two non-isolating 3D boxes on the (x − z)-plan. The first situation corresponds to a box B that contains an x-critical zero-crossing point (lying on the branch depicted in blue) as well as another zero-crossing point (lying on the branch depicted in red). The second situation corresponds to a box B that contains two different x-critical points of C among which one is a zero-crossing point.
Let (0, y0 , z0 ) be an x-critical zero-crossing point of C and suppose given a 3D box B that isolate the point (0, y0 , z0 ) from other zero-crossing points of C.
A first computational step consists in refining the box B until (0, y0 , z0 ) is the only x-critical point of C that lies in the refined box. To do so, a straighforward approach consists in computing disjoint boxes that isolate the real zeros of the following system
Sx :=
g(x, y, z) = 0 h(x, y, z) = 0 d(x, y, z) = 0
As for the system in (3), this can be done by computing a univariate representation of the solutions followed by a numerical approximation. This provides a set of disjoint 3D boxes that isolate the real zeros of Sx . Thenceforth, it suffices to refine the box B until it overlaps with only one box of this set to get a refined box that contains only the x-critical point (0, y0 , z0 ). Let denote to this new box by B . The second computational step consists in refining the box B until all the intersections with C lie only on the desired sides. Let B be a box resulting from the refinement of B up to a certain precision. We denote by ymin , ymax (resp. zmin , zmax ), the rational endpoints of its interval on y (resp. z). Then, the box B admits intersections only on the sides that are orthogonal to the x-axis if and only if the following algebraic systems do not have real zeros.
Fig. 1. Two zero-crossing points 175
g(x, ymin , z) = 0 h(x, ymin , z) = 0
g(x, ymax , z) = 0 h(x, ymax , z) = 0
g(x, y, zmin ) = 0 h(x, y, zmin ) = 0
g(x, y, zmax ) = 0 h(x, y, zmax ) = 0
2019 IFAC SSSC 98 Sinaia, Romania, September 9-11, 2019
Yacine Bouzidi et al. / IFAC PapersOnLine 52-17 (2019) 94–98
Consequently, to compute the desired box, our algorithm consists in refining successively the initial box and checking the existence of real solutions for the above algebraic systems. This can be done for example by computing a rational univariate representation (5) for each of these systems and then checking the existence of real roots of univariate polynomials using classical dedicated methods. The algorithm stops when no real solution is found. 5. CONCLUSION We addressed in this paper some aspects related to the study of the stability of differential time-delay systems with commensurate delays. More precisely, following the approach in common use, based on the computation of the critical pairs, we proposed a new strategy for analyzing the asymptotic behavior at these points, that gets rid of Puiseux series computations. Instead, our strategy proceeds by locally approximating at the critical pairs, the quasi-polynomial by a polynomial, and then by analyzing the asymptotic behavior at the zeros of this polynomial seen as specific points of algebraic curve in R3 . The implementation of this approach together with a complete complexity study will be the subject of an ongoing work. In particular, attention will be put on the comparison from the complexity point of view between our approach and the Puiseux-based counterpart. In future work, an interesting problem is to investigate more deeply the impact of numerical errors, due to the truncation as well as the local approximation at the critical pairs, on the certification of our result. REFERENCES Boussaada, I. and Niculescu, S.I. (2016). Tracking the algebraic multiplicity of crossing imaginary roots for generic quasipolynomials: A vandermonde-based approach. IEEE Transactions on Automatic Control, 61(6), 1601–1606. Bouzidi, Y., Poteaux, A., and Quadrat, A. (2016). Computer algebra methods for the stability analysis of differential systems with commensurate time-delays. IFACPapersOnLine, 49(10), 194–199. Gu, K., Kharitonov, V.L., and Chen, J. (2003). Stability of Time-Delay Systems. Birkh¨ auser. Li, X.G., Niculescu, S.I., and C ¸ ela, A. (2015). Analytic Curve Frequency-Sweeping Stability Tests for Systems with Commensurate Delays. Springer. Li, X.G., Niculescu, S.I., C ¸ ela, A., Wang, H.H., and Cai, T.Y. (2013). On computing puiseux series for multiple imaginary characteristic roots of lti systems with commensurate delays. IEEE Transactions on Automatic Control, 58(5), 1338–1343. Marshall, J.E., Gorecki, H., Korytowski, A., and Walton, K. (1992). Time-Delay. Systems: Stability and Performance Criteria with Applications. London: Ellis Horwood. Moore, R.E. (1979). Methods and applications of interval analysis, volume 2. Siam. Niculescu, S.I. (2001). Delay Effects on Stability: A Robust Control Approach. Number 269 in Lecture Notes in Control and Information Sciences. Springer. Olgac, N. and Sipahi, R. (2002). An exact method for the stability analysis of time-delayed linear time-invariant 176
(LTI) systems. IEEE Transactions on Automatic Control, 47(5), 793–796. Rouillier, F. (1999). Solving zero-dimensional systems through the rational univariate representation. J. of Applicable Algebra in Engineering, Communication and Computing, 9(5), 433–461.