Rachid Malti
email: rachid.malti@ims-bordeaux.fr

Milan R Rapaić
email: rapaja@uns.ac.rs

Vukan Turkulov
email: vukan_turkulov@uns.ac.rs

A unified framework for exponential stability analysis of irrational transfer functions in the parametric space

Keywords: Stability, Distributed parameter system, Irrational transfer function, Time-delay system, Interval arithmetics

This paper presents a unified framework for exponential stability analysis of linear stationary systems with irrational transfer functions in the space of an arbitrary number of unknown parameters. Systems described by irrational transfer functions may be of infinite dimension, typically having an infinite number of poles and/or zeros, rendering their stability analysis more challenging as compared to their finite-dimensional counterparts. The analysis covers a wide class of distributed parameter systems, time delayed systems, or even fractional systems. First, it is proven that, under mild hypotheses, new poles may appear to the right of a vertical axis of abscissa (imaginary axis, when = 0) through a continuous variation of parameters only if existing poles to the left of cross the vertical axis. Hence, by determining parametric values for which the crossing occurs, known as stability crossing sets (SCS), the entire parametric space is separated into regions within which the number of right-half poles (including multiplicities) is invariant. Based on the aforementioned result, a constraint satisfaction problem is formulated and a robust estimation algorithm, from interval arithmetics that uses contraction and bisection, is used to solve it. Applications are provided for determining the SCS of (i) a controlled parabolic 1D partial di↵erential equation, namely the heat equation, in finite and semi-infinite media, (ii) time-delay rational systems with distributed and retarded type delays, (iii) fractional systems, providing stability results even for incommensurate di↵erentiation orders.

Introduction

Many engineering systems exhibit dynamical behaviors that can be captured by partial di↵erential equations (PDE), or delayed ordinary and/or PDE. These distributed parameter systems (DPS) yield irrational transfer functions that may involve some functions of the Laplace variable like, exp, sinh, cosh, etc. A wide variety of transfer functions of DPS, solutions of PDE, is exhibited in [START_REF] Du↵y | Transform Methods for Solving Partial Di↵erential Equations[END_REF][START_REF] Curtain | Transfer functions of distributed parameter systems: A tutorial[END_REF]. Such transfer functions are usually infinite dimensional, with an infinite number of poles and/or zeros.

Context

Stability of irrational systems is a rich field of study in Automatic Control. It would be pretentious to claim an exhaustive bibliographical review covering all its facets. For time domain approaches, interested readers are referred to some recent papers [START_REF] Mondié | Lyapunov stability tests for linear time-delay systems[END_REF][START_REF] Prieur | Feedback stabilization of a 1-d linear reaction-di↵usion equation with delay boundary control[END_REF][START_REF] Katz | Constructive method for finitedimensional observer-based control of 1-d parabolic pdes[END_REF][START_REF] Aghayan | Robust stability of uncertain fractional order systems of neutral type with distributed delays and control input saturation[END_REF][START_REF] Katz | Delayed finite-dimensional observerbased control of 1-d parabolic pdes[END_REF], and the cited references therein. Several methods are developed for analyzing stability of some specific DPS such as clamped-free damped string [START_REF] Lhachemi | Input-to-state stability of a clamped-free damped string in the presence of distributed and boundary disturbances[END_REF], telegrapher's equation [START_REF] Sano | Stability analysis of the telegrapher's equations with dynamic boundary condition[END_REF], heat equation [START_REF] Li | The stability and stabilization of heat equation in non-cylindrical domain[END_REF], wave equation [START_REF] Gao | Stability of solution for uncertain wave equation[END_REF]. As compared to the mentioned references, this paper analysis stability of more general irrational systems on the basis of a frequency domain approach. Hence, a state of the art of frequency domain approaches is presented below.

Nyquist stability criterion is an important tool that has been continuously used for stability analysis of irrational transfer functions of DPS in [START_REF] Chait | A nyquist stability criterion for distributed parameter systems[END_REF][START_REF] Logemann | Circle criteria, small-gain conditions and internal stability for infinite-dimensional systems[END_REF], time-delay FoS in [START_REF] Zhang | Nyquist-based stability analysis of non-commensurate fractional-order delay systems[END_REF], and incommensurate FoS in [START_REF] Ivanova | Stability and resonance conditions of second-order fractional systems[END_REF]. Some other wellknown finite-dimensional results, such as circle criterion and small gain theorem, have been generalized to a large class of DPS in [START_REF] Logemann | Circle criteria, small-gain conditions and internal stability for infinite-dimensional systems[END_REF].

Another approach is based on root continuity when parameters vary continuously. It has been widely investigated for stability analysis of time-delay-systems (TDS) in [START_REF] Gryazina | The D-decomposition theory[END_REF][START_REF] Neimark | D-partition and robust stability[END_REF][START_REF] Lee | On the tau-decomposition method of stability analysis for retarded dynamical systems[END_REF]El'sgol'ts and Norkin, 1973), and is based on the splitting of the parametric space into multiple regions, with the number of unstable poles being invariant inside each region. Related methods applicable to TDS, have been considered in [START_REF] Walton | Direct method for tds stability analysis[END_REF]; [START_REF] Marshall | Timedelay systems Stability and Performance Criteria with Applications[END_REF], on the basis of the following parametrization of the characteristic function:

F (s, ⌧ ) = A(s) + C(s)e ⌧ s
(1)

where A(s) and C(s) 6 = 0 are polynomials of the Laplace variable s. Conditions are determined on the delay that destabilizes the system. Multiple generalizations of the the method were proposed:

• by the authors themselves, including multiple commensurate delays in [START_REF] Walton | Direct method for tds stability analysis[END_REF])

F (s, ⌧ ) = A(s) + N X n=1 C n (s)e n⌧ s (2)
• by [START_REF] Bonnet | Stabilization of some fractional delay systems of neutral type[END_REF] to fractional TDS of neutral type with a fixed commensurate di↵erentiation order ↵ > 0 where A(s) and C(s) in (1) (or C n (s) in ( 2)) are extended to pseudo-polynomials (polynomials of s ↵ -variable).

• by [START_REF] Abusaksaka | Bibo stability of some classes of delay systems and fractional systems[END_REF] to systems involving terms like e ⌧ s ↵ (sometimes referred to as fractional delays) with arbitrary, but fixed ↵, A(s) and C(s) being polynomials of s ↵ :

F (s, ⌧ ) = A(s) + C(s)e ⌧ s ↵ (3)
The Walton-Marshall method is restrictive in multiple point of views: (i) only commensurate delays are considered: all delays are multiples of a constant, (ii) the extension to FoS considers all di↵erentiation orders to be commensurate (iii) in the latter extension, the authors consider systems with e s ↵ ⌧ , where ↵ is the fixed commensurate order of A(s) and C(s) pseudo-polynomials. The method developed in this paper allows handling systems with much more complex structures, and much more complex parametrizations. The structure may include systems described by transfer functions which are arbitrary irrational functions of the Laplace variable s, with terms like e p s , cosh( p s), or sinh( p s), s ↵ and many others. The parametrization may be diverse as stability can be investigated for arbitrary combinations of parameters including (but not limited to) transfer function coe cients, incommensurate delays, incommensurate di↵erentiation orders.

A novel method, based on the application of Rouché's theorem, has recently been developed for stability analysis of (i) a wide class of systems including a class of DPS in (Turkulov et al., 2023a), (ii) TDS in (Turkulov et al., 2023b) and (iii) fractional order systems (FoS) in [START_REF] Rapaić | On stability regions of fractional systems in the space of perturbed orders[END_REF]. Given a parametrized transfer function and an arbitrary parametric point, the proposed method is capable of identifying the surrounding points/regions in the parametric space for which the number of unstable poles remains invariant. Iterating the method, allows reaching the boundary of the stability domain, however with an important time complexity when the whole stability/instability region is explored. Instead, the method developed in this paper searches straightforwardly for the domain boundaries known as the Stability Crossing Sets (SCS), i.e. a set of surfaces for which there is at least a pole crossing the imaginary axis. This is done in a unified framework of linear irrational systems that includes classes of DPS, TDS, and FoS. Such approaches have been successfully demonstrated in the particular case of TDS of retarded type with two and three independent delays [START_REF] Hale | Global geometry of the stable regions for two delay di↵erential equations[END_REF][START_REF] Gu | On stability crossing curves for general systems with two delays[END_REF]Morȃrescu et al., 2007;[START_REF] Sipahi | Extraction of 3D stability switching hypersurfaces of a time delay system with multiple fixed delays[END_REF][START_REF] Gu | Stability crossing set for systems with three delays[END_REF], providing insightful graphical representation of stability equivalence regions. To the best of authors' knowledge, such approaches have not been extended in the general framework of irrational transfer functions including DSP and FoS.

Various control problems of rational systems with uncertain parameters have been treated in the frequency domain in (Jaulin et al., 2001, section 7.4), based on interval arithmetics, such as robust stability, stability analysis, stability margins, and some basic control design. Di Loreto et al. ( 2007) extend some of these control problems to TDS. Regarding (robust) stability analysis, the authors of the latter paper formulate the problem of searching the whole stability/instability regions in the parametric space. In doing so, the problem turns out to be ill-posed in the stability region, as it requires testing the stability of all (hyper)-boxes of elementary sizes in a brute-force way as pointed out in [START_REF] Malti | Stability analysis of incommensurate elementary fractional systems using interval arithmetics[END_REF]. Our paper uses similar interval-arithmetics-techniques, based on contraction and bissection, to search however for the stability limits (SCS), instead of searching for the whole stability and instability regions, which ensures a well-posedness of the problem.

Contributions of the paper

A unified framework is proposed in this paper for exponential stability analysis of linear stationary systems in the frequency domain explicitly presented in a transfer function form. The proposed method is based on computing the SCS and applies to a wide range of irrational transfer functions [START_REF] Du↵y | Transform Methods for Solving Partial Di↵erential Equations[END_REF][START_REF] Curtain | Transfer functions of distributed parameter systems: A tutorial[END_REF].

The contributions of the paper can be summarized as follows:

• Clear hypotheses are set up under which root continuity is proven, for a wide class of transfer functions parameterized in a generic way.

• The problem of finding all characteristic-functionroots on the imaginary axis (and more generally on the axis of abscissa 2 R) is formulated as a constraint satisfaction problem (CSP).

• A robust algorithm from interval arithmetics, based on contraction and bisection [START_REF] Jaulin | Applied interval analysis[END_REF][START_REF] Jaulin | Actions of the hyperoctahedral group to compute minimal contractors[END_REF], which avoids all sources of errors due to floating point arithmetics and due to discretization, is used to solve this problem. As a consequence, a gap is bridged between algorithms from numerical analysis based on interval arithmetics and the stability analysis in automatic control.

• A comprehensive illustration of the method is provided on (i) a DPS, namely controlled heat equations in finite and semi-infinite media, (ii) a retarded and a distributed TDS, providing stability results even for incommensurate delays (not multiples of a constant), (iii) a FoS providing stability results even for incommensurate di↵erentiation orders which stability analysis is known to be more challenging than the commensurate ones.

To the best of authors' knowledge, no algorithm of similar universality has ever been reported previously in the literature for finding the SCS for such a wide class of linear irrational systems.

Notation

The set of non-negative and non-positive real numbers are respectively denoted by R + and R . C is the field of complex numbers, C + and C + are respectively the sets of all complex numbers with real part bigger than , and bigger or equal to . Consequently, C + 0 and C + 0 denote respectively the open and the closed right-half complex plane. Re(.) and Im(.) are respectively real and imaginary parts of a complex number, L {.} is the Laplace transform and ⇥ the cartesian product. The closure of a set is defined as its smallest closed superset.

Stability. Consider several subalgebras of transfer functions, as discussed in [START_REF] Curtain | An Introduction to Infinite-Dimensional Linear Systems Theory[END_REF][START_REF] Partington | H1 and bibo stabilization of delay systems of neutral type[END_REF]. A transfer function G is said to be stable in the H 1 sense if it is holomorphic and bounded in C + 0 . Next, consider the class A as the class of all transfer functions that are holomorphic in Re(s) > and continuous on the boundary Re(s) = , with finite and unique limit when s ! ± j1. For > 0, A is referred to as the subalgebra of -exponentially stable transfer functions. On the other hand, a transfer function belonging to A for some negative , but not for any 0, is referred to as exponentially unstable. Note that for 2 > 1 > 0,

A 2 ⇢ A 1 ⇢ A 0 ⇢ H 1 .
Characteristic function and stability crossing sets. When dealing with linear time-invariant (LTI) rational systems, the concept of a characteristic polynomial is well defined as transfer function denominator. It allows assessing system stability by studying root locii of the characteristic polynomial which are the only singularities of the transfer function.

When dealing with LTI, irrational systems denoted G(s, ✓), where s 2 C is the Laplace variable, and ✓ = (✓ 1 , • • • , ✓ n ) 2 T ⇢ R n is the parametric vector, the concept of characteristic function (instead of polynomial) f (s, ✓) requires some clarifications. In the present work, the characteristic function is understood as any function f (s, ✓) such that for every ✓ 2 T 1. f has no finite poles of any multiplicity, 2. f has branch points and/or essential singularities everywhere G has branch points and/or essential singularities1 , 3. all finite zeros of f match poles of G in both location and multiplicity.

The selection of the characteristic function is not unique. Several di↵erent characteristic functions may be constructed from a single transfer function, according to the points above. For example, a characteristic function of a given transfer function G(s, ✓) could be derived as

f (s, ✓) = Q(s, ✓) G(s, ✓) , (4)
with Q being an analytic function in C having zeros, with matching multiplicities, at the same locations as the zeros of G. Without loss of generality, transfer functions with poles-zeros cancellations are not handled in this paper and should be treated separately. Additionally, let ⌦ f, denote the stability crossing set (SCS) of f , i.e. a set of surfaces for which there is at least a zero of f crossing the axis Re(s) =

⌦ f, = {✓ 2 T | (9! 2 R) f ( + j!, ✓) = 0} .
(5)

Let NU f, (✓) denote the number of zeros of a characteristic function f (s, ✓) with real part bigger or equal to , where each zero is counted as many times as its multiplicity. Although multiple characteristic functions may be defined for an irrational transfer function, the stability crossing sets are well-defined and independent of the selection of f .

Problem formulation and hypotheses

Consider an LTI system with transfer function G(s, ✓) and characteristic equation f (s, ✓), ✓ 2 T ⇢ R n (with T being the parametric space). Choose 2 R such that G(s, ✓) has finite and unique limit for s ! ± j1, all singularities of f have real parts smaller or equal to , with ⌅ denoting the set of imaginary parts of all singularities of f on the axis Re(s) = . The problem under consideration is partitioning the parametric space into non-overlapping sub-regions i) f (s, ✓) has no zero on the axis Re(s) = , for all ✓ 2 T k , and

T k (k 2 1, . . . , K) such that for all k Re(s) = Im(s) Re(s) M R ✏ 1 ✏ ✏ 2 ✏ 3 + j! 2 + j! 3 + j! 1
ii) NU f, (✓ 1 ) = NU f, (✓ 2 ) for all ✓ 1 , ✓ 2 2 T k .
Formally, the following hypotheses are postulated:

(H1) For every ✓ 2 T , s 7 ! f (s, ✓) is analytic in C + 2 , except possibly in a finite set of points on the boundary, + j! for

! 2 ⌅ = {! 1 , . . . ! N }. For every (s, ✓) 2 ⇣ C + \ ( + j⌅) ⌘ ⇥ T , (s, ✓) 7 ! f (s, ✓) is con- tinuous.
(H2) For all + j! k , ! k 2 ⌅, there exists ✏ k > 0 such that f (s, ✓) 6 = 0 for all s satisfying Re(s) and Hypothesis (H1) serves two purposes. Firstly, it ensures that the transfer function is parameterized in a continuous fashion with respect to s and ✓, except possibly at a finite number of points on Re(s) = . Those singular points, illustrated by red dots in Fig. 1, are gathered within the set + j⌅. Secondly, (H1) determines the type of guaranteed stability one is able to seek for. If (H1) is satisfied 2 A function is analytic in a closed set if and only if it is analytic in every point of that set. For sets having nonempty boundary (like C + ) this actually means that the function is also analytic in some (possibly small) neighborhood of that boundary. In other words, a function analytic in a closed set is also analytic in some open superset that encloses that set. for some < 0, then it is possible to investigate for exponential stability, by localizing zeros of the characteristic function. If it is satisfied for some > 0, but not for any smaller , then the system under consideration is exponentially unstable, regardless of the position of the zeros of f . However, if (H1) is satisfied for 0, but not for any smaller , then neither exponential stability nor exponential instability can be claimed solely based on the position of roots of the characteristic equation. In this case, stability should be investigated by additional considerations, e.g. H 1 -stability can be claimed by establishing that the transfer function is bounded on C + 0 [START_REF] Partington | H1 and bibo stabilization of delay systems of neutral type[END_REF]).

|s j!|  ✏ k ,
The introduced hypotheses also limit the extent of the parametric space T . Hypothesis (H2) is a safeguard from cases in which zeros of the characteristic function in C + either originate from or terminate at the singularities located on the imaginary axis. A simple example would be a system with characteristic function f (s, K) = K + s 3 2 , having a zero with positive real part originating from the branch point at the origin when K < 0. This however is not permitted by (H2) which hence restricts the study of f (s, K) to K > 0. Hypothesis (H3) guards against cases in which zeros emerge from infinity in C + . It, for example, prevents negative values of delays (non-causal leads), and also safeguards against cases in which system structure abruptly changes for certain parametric values. An interesting practical example is the case of loop transfer function with relative degree zero, such as e.g.

G(s, K) = K 1 s s+1 , which closed-loop characteristic polyno- mial is f (s, K) = (s + 1) K(s 1)
. In this particular case, a single branch of the root locus originates at s = 1 for K = 0. By increasing K the root slides along the negative real axis, and reaches asymptotically 1 when K ! 1 , as illustrated in Fig. 2 . The value K = 1 is singular, since the order of the system changes, and for this particular value of K the characteristic polynomial has no zero. However when K > 1, a root pops up from +1, and decreases asymptotically to s = 1 as K grows. Consequently, the characteristic function of this system does not fulfill hypothesis (H3) when K 2 R. However, it is worth mentioning that (H3) allows investigating separately the stability for 0  K < 1. It is also important to stress that the hypothesis (H3) excludes a certain class of TDS of neutral type, as in [START_REF] Cooke | Discrete delay, distributed delay and stability switches[END_REF][START_REF] Baretta | Geometric stability switch criteria in delay di↵erential systems with delay-dependent parameters[END_REF][START_REF] Bonnet | Stabilization of some fractional delay systems of neutral type[END_REF][START_REF] Nguyen | A problematic issue in the walton-marshall method for some neutral delay systems[END_REF].

1.5. Some examples of system classes satisfying the hypotheses All in all, hypotheses (H1)-(H3) are the only conditions a transfer function must fulfill for the results of this paper to hold. These hypotheses are satisfied by wide classes of systems.

Some examples of such classes include (i) finite dimensional LTI systems, i.e. systems with rational transfer functions (ii) TDS of retarded type with no constraint on the delays, including those considered by [START_REF] Walton | Direct method for tds stability analysis[END_REF] (see the example in section 5.2), (iii) fractional incommensurate systems with no exponential term (see the example in section 6), (iv) fractional TDS of retarded type, similar to the ones considered in [START_REF] Bonnet | Stabilization of some fractional delay systems of neutral type[END_REF] however by releasing commensurability constraints on the time delays and/or on the fractional derivatives, (v) systems with fractional delays like in [START_REF] Abusaksaka | Bibo stability of some classes of delay systems and fractional systems[END_REF] but again without commensurability constraints (see the example in section 4.1), (vi) systems that are described by PDE [START_REF] Du↵y | Transform Methods for Solving Partial Di↵erential Equations[END_REF][START_REF] Curtain | Transfer functions of distributed parameter systems: A tutorial[END_REF], (see the example in subsection 4.1 and 4.2), (vii) systems that are solutions of fractional PDE [START_REF] Cvetićanin | Generalized timefractional telegrapher's equation in transmission line modeling[END_REF]. Transfer functions described in (vi) and (vii) typically involve hyperbolic functions of powers of rational and quasi-rational functions of the Laplace variable.

Paper layout

Within the present work and under hypothesis (H1)-(H3), we start by proving that the only way for the number of zeros of a characteristic function to change in C + , when the parameters vary continuously, is by crossing the vertical axis Re(s) = . This theoretical result, related to root continuity, is studied in section 2.

Then, in section 3 the problem of finding the SCS is formulated as a constraint satisfaction problem based on interval arithmetics, and solved using a robust algorithm based on contraction and bissection. Applications to DPS, TDS, and FoS are presented in sections 4, 5, and 6 respectively, before concluding.

Unlike many previously published results, mainly focused on time-delay systems (TDS) [START_REF] Hale | Global geometry of the stable regions for two delay di↵erential equations[END_REF][START_REF] Gu | On stability crossing curves for general systems with two delays[END_REF][START_REF] Sipahi | Extraction of 3D stability switching hypersurfaces of a time delay system with multiple fixed delays[END_REF][START_REF] Gu | Stability crossing set for systems with three delays[END_REF], a unified framework is proposed in the present paper that applies to a wide range of irrational systems, regardless of the number of investigated parameters. The only practical limitation is related to the exponential complexity of the proposed algorithm.

Root Continuity

The following claims will ensure continuity of all righthalf zeros of the characteristic function. Start by recalling the root continuity result from (Dieudonné, 1960, Theorem 9.17.4).

Theorem 1. Let D be an open set in C, T a metric space, f a continuous, complex valued function in C ⇥ T such that for each ✓ 2 T s 7 ! f (s, ✓) is analytic in D. Let D 1 be an open subset of D whose closure D 1 is compact and contained in D, and let ✓ 0 be such that no zero of f (s, ✓ 0 ) is on the boundary of D 1 . Then, there exists a neigborhood W of ✓ 0 such that for every ✓ 2 W 1. there are no zeros of f (s, ✓) on the boundary of D 1 , 2. the sum of the orders of zeros of f (s, ✓) belonging to D 1 is independent of ✓.

Dieudonné's Theorem 1 can be used to justify stability investigation by evaluating the SCS. To do so, it is necessary to prove analyticity of the characteristic function in a compact (closed and bounded) set of possible locations of its zeros in the right half complex plane. In this case, Theorem 1 guarantees that the only manner for the zeros to appear in the right-half complex plane, is by crossing the imaginary axis. Related conditions have been established for a wide class of TDS in the literature [START_REF] Cooke | Discrete delay, distributed delay and stability switches[END_REF][START_REF] Gu | On stability crossing curves for general systems with two delays[END_REF][START_REF] Bonnet | Stabilization of some fractional delay systems of neutral type[END_REF]. However, in a more general class of irrational transfer functions, Theorem 1 cannot be applied straightforwardly, because it is not trivial to find a set D 1 ⇢ C + which is simultaneously relevant for stability analysis and posses all the required properties.

The following claim extends Theorem 1 in an original way to support stability analysis of a more general class of systems described by irrational transfer functions which contain arbitrary irrational functions of the Laplace variable s, like e p s , cosh( p s), or sinh( p s). The claim ensures that zeros of the characteristic function cannot originate from infinity at the right of the axis Re(s) = , nor from the singular points lying on the stability boundary. Hence, the identification of the the SCS is su cient to reason about stability of such systems.

Theorem 2. Under hypotheses (H1)-(H3), given any connected subset S of T , if f ( + j!, ✓) 6 = 0 for all ! 2 R \ ⌅ and all ✓ 2 S, then for any ✓ 1 , ✓ 2 2 S the number of zeros of f in C + , counting multiplicities, is the same: (with ' 2 [ ⇡ 2 , ⇡ 2 ]). Let D 1 be the interior of D 1 . Due to hypotheses (H2), (H3), and the condition f ( + j!, ✓) 6 = 0 assumed in the statement of the Theorem, there can be no zeros of f on the boundary of D 1 for all ✓ 2 S. Therefore, D 1 satisfies all the assumptions of Theorem 1. Consequently, for all ✓ 0 2 S there exists a neighborhood W (✓ 0 ) ⇢ S such that the number of zeros of f inside D 1 , including multiplicities, is independent of ✓ 2 W (✓ 0 ). Notice further that this implies that the number of zeros (including multiplicities) of f is constant within C + , since due to (H2) there can be no right-half singularities within the "small" semicircles of D 1 , and due to (H3) there can be no roots in C + with modulus greater than M < R. Consequently, for every ✓ 0 2 S there exists a neighborhood W (✓ 0 ) such that the number of zeros in C + is independent of ✓ 2 W (✓ 0 ). In the remainder, it is proven by contradiction that the number of zeros of f with real part bigger or equal to , including multiplicities, does not change within the entire set S. Assume that the number of roots, counting multiplicities, is not the same for every ✓ 2 S. Then, there exists at least two parametric points ✓ 1 , ✓ 2 2 S such that the number of zeros at ✓ 1 is di↵erent from the one at ✓ 2 . Since S is connected by assumption, there exists a simple path P ⇢ S connecting ✓ 1 and ✓ 2 , and at least one point ✓ 0 2 P such that the number of zeros, counting multiplicities, is not constant in any neighborhood of ✓ 0 , which contradicts the existence of W (✓ 0 ) and completes the proof. ⇤ Definition 1. Given a parametric space T , and an LTI system with characteristic function f (s, ✓) with stability crossing set ⌦ f, , defined in (5), any two points ✓ 1 and ✓ 2 from T are said to belong to the same parametric region if there exists a continuous path within T connecting ✓ 1 with ✓ 2 which does not intersect ⌦ f, .

NU f (✓ 1 ) = NU f (✓ 2 ).
Notice that "belong to the same parametric region" is an equivalence relation. Therefore, the stability crossing set splits the parametric space T into multiple disjoint regions T k , k = 1, . . . , K  1. An immediate consequence of Theorem 2 is that given any parametric region T k of T , and any pair of points ✓ 1 , ✓ 2 2 T k , the number of zeros in C + , counting multiplicities, is the same for f (s, ✓ 1 ) and f (s, ✓ 2 ). Thus, under hypotheses (H1)-(H3), zeros of the characteristic function cannot suddenly appear or change multiplicity in C + . The only way for the number of zeros to change is by crossing the axis Re(s) = when system parameters vary continuously. Hence, the determination of the SCS, ⌦ f, defined in (5), allows characterizing all the regions in the parametric space inside which the characteristic function has the same number of zeros.

The required theoretical results for root continuity being established, the problem of finding the SCS is formulated and solved in the next section.

Robust estimation algorithm of the SCS

Since the parameters gathered in ✓ vary continuously, the problem of determining stability regions can be resumed to the problem of finding the SCS, defined in (5), for which the zeros of f (s, ✓) cross the imaginary axis. This problem is formulated as a Constraint Satisfaction Problem (CSP):

CSP : ( f ( + j!, ✓) = 0 (!, ✓) 2 {R + \ ⌅} ⇥ T . (6)
The searching space of ! is restricted to {R + \ ⌅}, because all real-valued impulse-response systems have complex conjugate poles. The solution set of all the feasible parameters of the CSP (6) is rewritten as

S = n (!, ✓) 2 R + \ ⌅ ⇥ T f ( + j!, ✓) = 0 o . (7)
A guaranteed and robust solution of this CSP can be obtained, in the searching domain, using interval arithmetics introduced by [START_REF] Moore | Interval analysis[END_REF]. The reader is also referred to [START_REF] Jaulin | Applied interval analysis[END_REF] for an introduction to interval arithmetics. The characterization of the whole solution set S in ( 7) can be formulated as a set inversion problem

S = f 1 (0) \ R + \ ⌅ ⇥ T (8)
and solved by guaranteed methods based on contraction and bisection. ⇢ [⇣]. Among the di↵erent types of contractors, [START_REF] Jaulin | Applied interval analysis[END_REF], the one used in the paper is based on the centered inclusion functions (Jaulin et al., 2001, section 2.4.3, p.33). Such a contractor is more adapted, as compared to the natural contractor (Jaulin et al., 2001, section 2.4.2, p.29), because the SCS are smooth and continuous functions with respect to the searching parameters. Hence, the mean value theorem, on which the centered inclusion function is based, applies. It is shown in [START_REF] Jaulin | Applied interval analysis[END_REF] that for this contractor the width of the inclusion function tends to the minimal inclusion function when the width of the searching interval tends to 0, which is not the case for the natural inclusion function. Recently, Jaulin (2023) developed an enhanced contractor based on the centered form, which is asymptotically minimal and hence provides even faster convergence.

Set Inversion Via Interval Analysis (SIVIA)

The SIVIA algorithm [START_REF] Jaulin | Set inversion via interval analysis for nonlinear bounded-error estimation[END_REF] has been proposed to solve constraint propagation problems using bisection. It has been used in di↵erent contexts such as: state and parameter estimation of non linear systems [START_REF] Raissi | Set membership state and parameter estimation for systems described by nonlinear di↵erential equations[END_REF], robust estimation of frequency domain models [START_REF] Khemane | Robust estimation of fractional models in the frequency domain using set membership methods[END_REF], parameter estimation in a Algorithm 1 SIVIA algorithm with only an outer enclosure end if 14: end procedure glucose model [START_REF] Herrero | Robust set-membership parameter estima-tion of the glucose minimal model[END_REF], or even to compute invariant sets of closed-loop control systems [START_REF] Romig | Using interval analysis to compute the invariant set of a nonlinear closed-loop control system[END_REF].

1: procedure SIVIA(in: [t], [⇣], ⌘, S ; out: S) 2: contractor (in: [⇣] ; out [⇣]) 3: if [t]([⇣]) = [0] then return; 4: end if 5: if w([⇣])  ⌘ then
Applying the SIVIA algorithm to (6) allows obtaining an outer enclosure S of the solution set S, if it exists3 as defined in ( 7), such that S ✓ S. SIVIA is a recursive algorithm based on partitioning of the parameter set into three regions: feasible, indeterminate and unfeasible. However, only indeterminate and unfeasible regions 3 can be obtained for the CSP (6). Hence the presentation of the SIVIA algorithm is restrained to these two cases. It uses an inclusion test [t] : {R + \ ⌅} ⇥ T ! N which is a function allowing to prove if a box [⇣] is unfeasible or undetermined. If unfeasible, the box is simply ignored. If undetermined, it is bisected and tested again unless its width w([⇣]) is less than a precision parameter ⌘ tuned by the user and which ensures that the algorithm terminates after a finite number of iterations. The outer enclosure S is then computed as a union of all undetermined boxes as indicated in line 6 of Algorithm 1. Line 2 calls the contractor at each execution of SIVIA.

In this paper, the SIVIA algorithm guarantees that all the outer enclosures of the SCS are obtained in the desired searching box which should be chosen carefully together with the precision factor ⌘ due to the exponential complexity of the algorithm. When no prior knowledge is available, the SIVIA algorithm may be initialized with a large searching box and a big precision factor ⌘ in such a way to obtain a rough estimation of the SCS. Then, the searching box can be narrowed and the precision factor enhanced to get more precise outer enclosures.

The SIVIA algorithm and the contractor are implemented using the Codac library in Python (Rohou et al., x 0 = 0.15 0 ϕ(t, 0)

L = 0.20m or L = ∞ x Figure 3: Thin rod of length L.
2022), in di↵erent contexts: a controlled heat equation in section 4, rational TDSs in section 5, and FoS is section 6.

Application to stability analysis of a distributed parameter system

Consider a one-dimensional heat di↵usion in a thin rod4 of length 0 < L  1. The rod is thermally isolated, except at its boundary cross-sections as in Fig. 3. At the initial end x = 0, the rod is subject to an adjustable thermal flux, '(t, 0) = @✓ @x (t, 0), where is the thermal conductivity. Two cases are treated: either the rod is of finite length or infinite. When the rod length is finite, the opposite end is kept at ambient temperature ✓(t, L) = 0. The corresponding model describing spatio-temporal distribution of the temperature along the rod is given by the well-known heat equation,

@ 2 ✓(t, x) @x 2 = 1 @✓(t, x) @t , (9)
where is the thermal di↵usivity of the medium. The temperature of the rod is measured at a crosssection x = x 0 = 0.15m, with 0  x 0  L. This temperature is then controlled remotely (see Fig. 4) using a proportional controller with a gain K. Additionally, time required to transmit data from the controller to the actuator, and from the sensor to the controller is denoted by ⌧ 2 , so that the total loop delay is ⌧ 0. Hence, the e↵ective transfer function of the proportional delayed controller is

G k (s) = Ke ⌧ s . ( 10
)
The objective of this example is to find the whole set of stabilizing controllers (10) in the K versus ⌧ plane in both cases: finite and semi-infinite rods. A similar problem is treated in (Morȃrescu and Niculescu, 2007), with the same delayed proportional controller G k (s), however applied to a rational system.

Semi-infinite spatial domain (L = 1)

In the semi-infinite domain, the transfer function between the input flux (s, 0) = L {'(t, 0)} and the temperature ⇥(s, x) = L {✓(t, x)} at a cross-section x is (Curtain

G(s) = ⇥(s, x) (s, 0) = e x p s p s , (11)
with the closed-loop transfer function G c`( s) =

Ke s⌧ e x p s f (s,K,⌧ )
and the closed-loop characteristic function

f (s, K, ⌧ ) = r s + Ke s⌧ e x p s . (12)
Due to the presence of p s a branch-cut is necessary. It is chosen along the negative real axis including the branch point 0 and 1. Hence, the characteristic function f (s, K, ⌧ ) is holomorphic in the complement of the branch-cut line of the complex plane and the arguments of s are restrained to

| arg(s)| < ⇡ (13)
In the special case when ⌧ = 0, all the roots of f (s, K, 0) can be determined analytically and their locus plotted versus K. The solutions of f (s, K, 0) = 0 are given by all the determinations, n = 0, ±1, ±2, . . ., of Lambert's W n function, see e.g. [START_REF] Corless | On the lambertw function[END_REF], s = x 2 W 2 n Kx , provided ( 13) is fulfilled. The root locus of f (s, K, 0) is plotted versus gain K, in Fig. 5, for the principal determination of the Lambert function W 0 (as the upper imaginary part) and for the determination W 1 (as the lower imaginary part). The other determinations of W n , n = 1, ±2, ±3, . . . are not represented as they are well beyond the scale of Fig. 5 towards 1. They may yield other poles crossing the imaginary axis for much greater values of K than the ones represented in Fig. 5, which shows that two poles pop up from the branch-cut when K ⇡ 3 ⇥ 10 3 . They cross the imaginary axis towards instability when K ⇡ 56 ⇥ 10 3 . The objective of the paper is to determine for what parametric values (here K and in a more general case K and ⌧ ) the poles cross the imaginary axis. No analytical solution is known for ⌧ 6 = 0 (and obviously K 6 = 0).

Finite spatial domain (L < 1)

In this case, the following transfer function is obtained for the heat di↵usion process under consideration [START_REF] Curtain | Transfer functions of distributed parameter systems: A tutorial[END_REF] 12) for di↵erent values of K, illustrating that poles pop up from the plane-cut. When the poles cross the imaginary axis (red squares), the system becomes unstable.

G(s) = ⇥(s, x) (0, s) = sinh((L x) p s ) p s cosh(L p s ) , (14)
Numbers on the curves correspond to values of K. The plane-cut is represented in blue along R .

with the closed-loop characteristic function

f (s, K, ⌧ ) = r s cosh ✓ L r s ◆ + Ke s⌧ sinh ✓ (L x) r s ◆ . (15)
As in the previous example, the proposed algorithm is applied with = 0, due to the presence of a branch point at the origin. \ {0}, the thermal system under consideration is slow. Additionally, one can see from Fig. 5 that the root locus of ( 12) for ⌧ = 0 crosses the imaginary axis for an ! in the interval (0.04, 0.05). Consequently, the first poles crossing the imaginary axis towards instability, in both cases (12) for any ⌧ and (15), are expected to be within the frequency interval [10 3 , 0.1] rad/sec. In case the obtained outer enclosure of [!] touches either the lower bound 10 3 or the upper bound 0.1 and/or in the case the resulting SCS appear with discontinuities, the searching interval [!] is enlarged and the SIVIA algorithm run again. The tolerance is set to an arbitrarily small value as the width of each element w(

The obtained outer enclosures of the stability crossing sets are plotted in the K versus ⌧ plane and the crossing frequencies along the imaginary axis in Fig. 6, in red for the semi-infinite rod and in blue for the finite rod. Since both transfer functions ( 11) and ( 14) are bounded on the imaginary axis, the systems under consideration are H 1stable in the parametric regions for which there are no roots in C + 0 (the lower left parts of Fig. 6, left). Hence, the red line (similarly the blue) in Fig. 6, left, splits the plane into two regions: a guaranteed stability region (the lower left containing K = 0 and ⌧ = 0) and a guaranteed instability region as the upper right one. In between, the outer enclosures of the SCS are plotted. Their width can be reduced by reducing the precision factor ⌘ in ( 16) at the price of a higher computational burden. In the case of the semi-infinite rod, for ⌧ = 0, the results of Figs.6 (red squares) are further confirmed by the intersection of the root-locus with the imaginary axis, in Fig. 5.

To the best of authors' knowledge, in both cases finite with ⌧ 6 = 0 and semi-infinite aluminium rods, there is no analytical method allowing to compute the roots of the characteristic function f in ( 12) or ( 15). The only alternative method for evaluating exponential stability of the control loop is based on a graphical representation that uses Cauchy's principle argument.

Application to stability analysis of time delay systems

Distributed delay system

Consider in this section a distributed delay system taken from (Turkulov et al., 2023b, example 9) Its characteristic function, computed by straightforward integration of the Laplace transform of (18), 16). The obtained outer enclosures of the SCS are plotted in the K versus ⌧ plane in Fig. 7, together with the outer enclosures of [!] along the axes = 0 and = 0.02 which confirms and extends the results of (Turkulov et al., 2023b, example 8) obtained using Rouché's theorem. The exterior of the blue curve corresponds to exponentially stable systems of class A 0.02 .

modeled by ẋ(t) = Z 0 ⌧ e K↵ x(t + ↵)d↵. (18)
f (s, K, ⌧ ) = s 2 + sK + 1 e ⌧ (s+K) , (19)

Time delay system of retarded-type

Consider a TDS system of retarded type, taken from (Turkulov et al., 2023b, example 10), with a characteristic function given by

f (s, ⌧ 1 , ⌧ 2 ) = s 2 + 2se ⌧1s + e ⌧2s . ( 20
)
Its stability is investigated with respect to ⌧ 1 and ⌧ 2 by setting the initial searching box to

[⇣] = ([!], [⌧ 1 ], [⌧ 2 ]) = ([0.1, 4], [0, 2], [0, 4]) .
The obtained outer enclosures of the SCS are plotted in the ⌧ 1 versus ⌧ 2 plane in Fig. 8, together with the outer enclosures of [!] along the imaginary axis.

The number of unstable poles, computed at any point inside the di↵erent regions, is indicated by a number. Fig. 8 confirms and extends again the results of (Turkulov et al., 2023b, example 10) obtained using Rouché's theorem. The interior of the blue contour indicates the region of the parametric space corresponding to exponentially stable systems of class A 0.02 .

Application to stability analysis of a fractional system

Consider in this section a fractional system, taken from (Rapaić and Malti, 2019, example 3), with a characteristic function given by f (s, ↵ 1 , ↵ 2 ) = s ↵2 + 2s ↵1 + 1 ( 21) and (↵ 1 , ↵ 2 ) 2 {R + \ {0}} 2 . The proposed algorithm is applied for determining the SCS of f (s, ↵ 1 , ↵ 2 ) with = 0. The characteristic function (21) complies with the hypotheses (H1)-(H3). The branch point at s = 0 is excluded from the searching box (⌅ = {0}), initialized at: ([0.1, 6], [0.01, 15], [0.01, 20]) with a tolerance set arbitrarily to a small value as in ( 16). The obtained outer enclosures S of the SCS, is plotted in Fig. 9, together with the outer enclosures of [!] along the axes = 0 and = 0.02, which confirms the results of (Rapaić and Malti, 2019, example 3). Hence, all fractional systems, incommensurate and commensurate, that have di↵erentiation orders in the lower-left region of Fig. 9, have no unstable poles. Considering additionally that G(s) = f 1 (s, ↵ 1 , ↵ 2 ) is bounded on C + 0 , one may also conclude that the parametric region with NU f,0 = 0 is the region of H 1 -stability.

[⇣] = ([!], [↵ 1 ], [↵ 2 ]) =

Conclusions

A unified framework has been presented in this paper for exponential stability analysis of irrational transfer func- tions in the frequency domain. First, it has been proven that the only way for the number of zeros of a characteristic function to change in the right of a vertical axis of abscissa , when its parameters vary continuously, is by crossing that vertical axis. Based on this theoretical result, the problem of finding the set of parameters for which poles cross the vertical axis has been formulated as a constraint satisfaction problem. This original formulation is universal as it applies to a wide class of irrational transfer functions. It has been solved using the robust SIVIA algorithm, based on interval arithmetics, that uses contraction and bissection. It has been successfully applied for determining the SCS of (i) a controlled heat partial di↵erential equation, in finite and semi-infinite media, (ii) time-delay systems with distributed and retarded type delays, (iii) fractional systems, providing stability results even for incommensurate di↵erentiation orders. The proposed algorithm may be used for any number of transfer function parameters. The only limitation is related to the timecomplexity of the SIVIA algorithm which is known to be exponential in terms of the number of parameters. Hence, any prior knowledge, used to reduce the size of the initial search box, increases procedure e ciency. The examples presented in the paper have been intentionally limited to two parameters for graphical representation purposes.

Figure 1 :

 1 Figure 1: Illustration of the notation involved in the hypotheses (H1)-(H3) and in the proof of Theorem 2.

  and all ✓ 2 T .(H3) There exists M > 0 such that f (s, ✓) 6 = 0 for all ✓ 2 T and all s such that Re(s) and |s | M .

Figure 2 :

 2 Figure 2: Root locus of closed-loop system with loop transfer function G(s, K) = K 1 s 1+s .

Proof 1 .

 1 According to (H1), there exists an open set D containing C + \ ( + j⌅) such that f is continuous on D ⇥ T . Let M satisfy (H3), and select an arbitrary R max{M, |! 1 |, . . . , |! N | }, with ! k 2 ⌅. Consider ✏ k , introduced in (H2), and define 0 < ✏  min k ✏ k . Introduce set D 1 as the closed set bounded by the axis Re(s) = excluding +j[! k ✏, ! k +✏], "small" semi-circles +j! +✏ j' and the "big" semi-circle + Re j'

  Let [x] = [x, x] design an interval which is a closed, bounded, and connected set of real numbers. The CSP (6) can be solved by a contractor C, which is an operator which permits to reduce the domain [⇣] = ([!], [✓]) without any bisection. Hence, contracting the box [⇣] means replacing it by a smaller box [⇣] ⇤ such that the solution set S remains unchanged, i.e. S ⇢ [⇣] ⇤

Figure 4 :

 4 Figure4: A feedback control loop, where G(s) is the transfer function between the heat flux and the temperature at x 0 , for either a semiinfinite or a finite rod.

Figure 5 :

 5 Figure5: Root locus of f (s, K, 0) in (12) for di↵erent values of K, illustrating that poles pop up from the plane-cut. When the poles cross the imaginary axis (red squares), the system becomes unstable. Numbers on the curves correspond to values of K. The plane-cut is represented in blue along R .

4. 3 .

 3 Robust estimation algorithm of the SCS Both characteristic functions (12) and (15) comply with the hypotheses (H1)-(H3). They have a branch point at s = 0, which is excluded from the searching box (⌅ = {0}), initialized at [⇣] = ([!], [K], [⌧ ]) = [10 3 , 0.1], [0, 7 ⇥ 10 4 ], [0, 100] . Although theoretically the searching interval of ! is R +

  [⇣]) divided by 2 9 ⌘ = w([⇣])/2 9 (16) ⌘ = (0.193 • 10 3 ,136.7, 0.195)

Figure 6 :

 6 Figure 6: Left: Outer enclosures of the SCS in the finite and semiinfinite aluminium rods in the K versus ⌧ plane. The number of unstable poles is indicated inside each region (by the corresponding color). Right: Outer enclosure of the crossing frequencies, along the positive imaginary axis. The red square (right and left plots) corresponds to the SCS of f (s, K, 0) in (12) as predicted in the root locus of Fig.5.

Figure 7 :

 7 Figure 7: Left: Outer enclosures of the SCS of f (s, K, ⌧ ) in (19). When = 0.02, only the external part is sketched delimiting systems guaranteed to belong to A 0.02 ; the number of unstable poles is indicated inside each red region. Right: Outer enclosure of the crossing frequencies of the axis = 0.02 in blue and of the imaginary axis in red.

  complies with the hypotheses (H1)-(H3), with no singularities in C + . Hence the boundary of instability region is investigated by choosing = 0 and = 0.02 . The searching box is initialized at [⇣] = ([!], [K], [⌧ ]) = ([0.01, 2], [0, 0.4], [0, 20]), and the tolerance factor set as in (

Figure 8 :

 8 Figure 8: Left: Outer enclosures of the the SCS of f (s, ⌧ 1 , ⌧ 2 ) in (20). When = 0.02, only the internal part is sketched delimiting systems guaranteed to belong to A 0.02 ; the number of unstable poles is indicated in each red region. Right: Outer enclosure of the crossing frequencies along the axis Re(s) = 0.02 in blue and along the imaginary axis in red.

Figure 9 :

 9 Figure 9: Left: Outer enclosures of the SCS of f (s, ↵ 1 , ↵ 2 ) in (21). The number of unstable poles is indicated inside each region. Right: Outer enclosure of the crossing frequencies, along the positive imaginary axis.

The distinction between di↵erent types of singularities is made on the basis of the shape of the Laurent series, as studied in complex analysis (see e.g.[START_REF] Lang | Graduate texts in Mathematics -Complex Analysis[END_REF]).

The original SIVIA algorithm allows computing additionally an inner enclosures: S ✓ S ✓ S. However, the CSP formulated in (6), can only yield an outer enclosure S.

For numerical application, the rod is considered of aluminium type with a thermal di↵usivity = 98.8 ⇥ 10 6 m 2 /s and a conductivity = 237 Wm 1 K 1 (see e.g.[START_REF] Baehr | Heat and Mass Transfer[END_REF]).