Alejandro Martínez-González
email: alejandro.martinez@l2s.centralesupelec.fr

César Fernando

Méndez Barrios

Silviu-Iulian Niculescu
email: silviu.niculescu@l2s.centralesupelec.fr

César-Fernando Méndez-Barrios

Some remarks on the regular splitting of quasi-polynomials with two delays. Characterization of double roots in degenerate cases

Keywords: Delay systems, Characteristic Roots, Stability, Weierstrass Polynomial, Asymptotic Behavior

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

The stability analysis of linear time-invariant (LTI) systems with time-delay have been extensively studied in the past decades and there exists an abundant literature covering the subject (see, for instance, [START_REF] Niculescu | Delay Effects on Stability: A Robust Control Approach[END_REF]; [START_REF] Gu | Stability of time-delay systems[END_REF]; [START_REF] Michiels | Stability, Control, and Computation for Time-Delay Systems. An Eigenvalue-Based Approach[END_REF]; [START_REF] Li | Analytic Curve Perspective for Time-Delay Systems[END_REF] and the references therein).

As discussed in Chen et al. (2010a,b), even in the the simplest case of a single constant delay, the stability tests are not easy to perform. Such a difficulty arises from the fact that that the delay systems are infinitedimensional and the corresponding characteristic function is, in fact, a quasi-polynomial that always have an infinite number of characteristic roots (see, for instance, [START_REF] Michiels | Stability, Control, and Computation for Time-Delay Systems. An Eigenvalue-Based Approach[END_REF] and the references therein). In the retarded case, by using an appropriate continuity argument1 , the problem of (exponential) stability can be reduced to (i) the detection and (ii) the analysis of the behavior of the characteristic roots located on the imaginary axis2 in the case when such roots exist3 .

In this context, one the problem particularly treated in the literature was to understand how the delay parameters affect the behavior of the critical characteristic roots, and to explicitly compute the stability domains in the delayparameter space. In the case of a single or commensurate delays, such domains reduce to a finite number of delay intervals. For more insights on the existing methods, one may refer to [START_REF] Niculescu | Delay Effects on Stability: A Robust Control Approach[END_REF]; [START_REF] Olgac | An exact method for the stability analysis of time-delayed linear time-invariant (LTI) systems[END_REF]; [START_REF] Gu | Stability of time-delay systems[END_REF]; [START_REF] Michiels | Stability, Control, and Computation for Time-Delay Systems. An Eigenvalue-Based Approach[END_REF]; [START_REF] Li | Analytic Curve Perspective for Time-Delay Systems[END_REF]. To the best of the authors' knowledge, the complete characterization of the stability domains in the incommensurate delays case is still an open problem.

Next, it is well known that the roots of polynomials are continuous functions of the coefficients as long as the leading coefficient does not vanish, see, e.g., [START_REF] Knopp | Theory of functions. Parts I and II[END_REF]. Furthermore, in the case of simple roots, these functions are also differentiable. However, in the case of multiple roots, such a property does not necessarily hold and Puiseux series may be used to perform the analysis.

In the commensurate delays case, these conclusions are also valid for quasi-polynomials of retarded type, where the characteristic roots are seen as functions of one variable -the delay parameter (see, for instance, [START_REF] Michiels | Stability, Control, and Computation for Time-Delay Systems. An Eigenvalue-Based Approach[END_REF]; [START_REF] Li | Analytic Curve Perspective for Time-Delay Systems[END_REF]). By using an operatorbased approach, such an idea was exploited by [START_REF] Chen | An eigenvalue perturbation approach to stability analysis, part II: when will zeros of time-delay systems cross imaginary axis[END_REF], where the authors characterized the asymptotic behavior of multiple critical characteristic roots in semisimple 4 and some not semi-simple 5 cases.

Based on the remarks above, the asymptotic behavior of multiple critical characteristic roots with respect to several (incommensurate) delays needs a deeper understanding (degenerate cases classification and characterization).

As mentioned before, one of the natural and standard approaches to handle such a problem is to study of the properties of the corresponding Puiseux series solutions. Such an idea, exploited in the commensurate delays case (for some insights, see [START_REF] Li | On computing Puiseux series for multiple imaginary characteristic roots of lti systems with commensurate delays[END_REF]) cannot be extended straightforwardly to the incommensurate delays. In fact, as in the multivariate polynomial case, such a method has several drawbacks and, in some situations, such Puiseux series solutions do not always exist (see Subsection 3.1 below as well as various examples proposed by [START_REF] Alonso | Local parametrization of space curves at singular points[END_REF]; [START_REF] Aroca | Puiseux parametric equations of analytic sets[END_REF]). In order to avoid such situations, the problem needs to be well-posed (see, e.g., [START_REF] Monforte | Formal Laurent series in several variables[END_REF]).

The migration of double characteristic roots depending on two parameters is studied in [START_REF] Gu | Migration of double imaginary characteristic roots under small deviation of two delay parameters[END_REF]; [START_REF] Irofti | Some insights into the migration of double imaginary roots under small deviation of two parameters[END_REF], where, without invoking Puiseux series, under appropriate assumptions 6 , the authors proposed a more conventional geometric approach to classify the asymptotic behavior of the critical characteristic roots in the parameter space. Such a method allows covering some particular asymptotic behavior in the delay-parameter space when the parameters change in the neighborhood of the critical point. Next, in the multiple delays case, [START_REF] Li | An iterative frequency-sweeping approach for stability analysis of linear systems with multiple delays[END_REF] proposed an iterative frequency-sweeping approach, by leaving one delay parameter free, while the others are fixed. Thus, the critical characteristic roots depend on only one parameter, and the analysis is performed by using the method developed by the same authors in [START_REF] Li | Analytic Curve Perspective for Time-Delay Systems[END_REF]. Finally, in the case of two parameters, some analysis was proposed by [START_REF] Maurer | Puiseux expansion for space curves[END_REF]; [START_REF] Lipman | Quasi-ordinary singularities of surfaces in C 3[END_REF], by using an appropriate parametrization of a given surface. This idea was exploited by [START_REF] Martínez-González | Some Remarks on the Asymptotic Behavior for Quasipolynomials with Two Delays[END_REF] for the computation of roots of quasi-polynomials with two incommensurate delays without any discussion on multiple critical characteristic roots.

Taking into consideration the previous discussion, the main contribution of this paper is threefold:

• Give conditions for the existence of multi-parameter Puiseux series and give a classification of its solutions; • Relax the non-degeneracy assumption proposed by [START_REF] Irofti | Some insights into the migration of double imaginary roots under small deviation of two parameters[END_REF]; • Propose an appropriate algorithm for studying the corresponding asymptotic behavior and related its splitting properties.

More precisely, under some particular conditions that guarantee the existence of well-defined multi-parameter solutions, we propose to relax the assumption for the existence of a parametrization around the multiple roots as well as its splitting properties. To complete the presentation, an explicit computation of the associated Puiseux series solution is proposed. To the best of the authors' knowledge, all these contributions represent a novelty in the open literature.

The remaining paper is organized as follows: Section 2 includes some preliminary results. In Section 3 some moti-6 "least degenerate" double critical characteristic roots vating example is outlined and the problem formulation is stated. Section 4 is devoted to the main results. In Section 5, several numerical examples are presented. Finally, some concluding remarks end the paper.

Throughout the paper, the following notations will be adopted: for z ∈ C, arg (z) ∈ [0, 2π), (z) ( (z)) denote the argument, real (imaginary) part of z, respectively. Next, R + denotes the set of positive real values, C[x] the ring of polynomials and C{x} the ring of convergent power series.

PRELIMINARY RESULTS

In the sequel, a dynamical LTI system of retarded type including two delays τ 1 and τ 2 is considered. For the sake of brevity and simplicity, its characteristic function f is given by the following quasi-polynomial:

f (s, τ 1 , τ 2 ) := p 0 (s) + p 1 (s)e -τ1s + p 2 (s)e -τ2s
, (1) where p k are polynomials given as,

p 0 (s) = s n + n-1 =0 a 0 s , p k (s) = n-1 =0 a k s , k ∈ {1, 2}.

Local Representation Around Multiple Roots

It is possible to reduce the analytic properties of f (z, x) to some algebraic ones. To this purpose, consider the following result [START_REF] Mailybaev | On the Weierstrass preparation theorem[END_REF]:

Theorem 1. (Weierstrass Preparation Theorem). Suppose that f (z, x) is an analytic function vanishing at the singular point z 0 ∈ C, x 0 ∈ C n , where z = z 0 is an m-multiple root of the equation f (z, x) = 0, i.e., f (z 0 , x 0 ) = ∂f ∂z = • • • = ∂ m-1 f ∂z m-1 = 0, ∂ m f ∂z m = 0, ( 2
) where derivatives are evaluated at (z 0 , x 0 ).

Then, there exists a neighborhood U 0 ⊂ C n+1 of the point (z 0 , x 0 ) ∈ C n+1 in which the function f (z, x) can be expressed as

f (z, x) = W (z, x) b (z, x) , (3) where W (z, x) is given by (z -z 0 ) m + w m-1 (x) (z -z 0 ) m-1 + • • • + w 0 (x) , (4)
and w 0 (x),. . . ,w m-1 (x), b (z, x) are analytic functions uniquely defined by the function f (z, x), and

w i (x 0 ) = 0, b (z 0 , x 0 ) = 0.
In [START_REF] Martínez-González | Asymptotic Analysis of Multiple Characteristics Roots for Quasi-polynomials of Retarded-Type[END_REF] the authors proposed an explicit method for computing the associated Weierstrass polynomial of (1). Remark 1. It can be seen from Theorem 1 that, since b (z, x) is an holomorphic non vanishing function at (0, 0).

Then, there must exist some neighborhood U ⊂ C n+1 at which b(z, x) preserves the same property. Hence, based on this observation, we can ensure that the roots behavior of a given quasi-polynomial f in the neighborhood U will be completely described by the roots behavior of W (x, x), see, for instance, [START_REF] Hörmander | An introduction to complex analysis in several variables[END_REF].

Splitting Properties

In the case of one parameter, the quasi-polynomial f defines an appropriate plane curve given by f = 0 (see [START_REF] Wall | Singular points of plane curves[END_REF]; [START_REF] Walker | Algebraic Curves[END_REF]), and the root locus of a multiple root of f is characterized by its branches which explicitly describe a solution curve C ∈ C 2 . This characterization is composed by a finite union of r-branches s j (τ 1/mj ) which can be parametrized as (s, τ ) = (ϕ(τ ), τ mj ). In the form of Puiseux series, the root locus gives rise to the following branching:

s jσ (τ ) = c jσ τ m j +o (|τ | mj ) , j = 0, . . . , r-1, σ = 1, . . . , m j (5) where each branch has multiplicity m j such that m = m 1 + m 2 + • • • + m r .
Definition 2.1. The trivial solution s * = 0 has a Complete Regular Splitting (CRS) property at τ * = 0 if c jσ = 0, ∀j. Next, if some of the coefficients c jσ for which m j = 1 may be equal to zero, then the trivial solution has a Regular Splitting (RS) property at τ * = 0. Finally, in the remaining cases of the coefficients c jσ , the trivial solution has a Non-Regular Splitting (NRS) property at τ * = 0.

In Martínez-González et al. (2019a), the authors proposed a methodology based on the first partial derivatives to analyzed the splitting behavior, see Figure 1. It is worth mentioning that these results cover only the commensurate delays case and the extension to incommensurate delays is not straightforward.

-0.6 -0.4 -0.2 0 0.2 0.4 ℜ(s) 0.8 0.9 1 1.1 1.2 1.3 ℑ(s)
Fig. 1. Illustration of Completely Regular Splitting (CRS) property in the case of a triple critical characteristic root s * = j

Quasi-ordinary Roots

In order to extend the above approach to the multiparameter case, the so-called Abhyankar-Jung Theorem must be taken into consideration (see, for instance, [START_REF] Abhyankar | On ramification of algebraic functions[END_REF]). To this end, the notion of quasiordinary polynomials must be introduced. Definition 2.2. Let W (z, x 1 , x 2 ) be a Weierstrass polynomial (3). Then W is a quasi-ordinary polynomial with respect to z, if the discriminant ∆ is of the form

∆(x 1 , x 2 ) = x p1 1 x p2 2 V (x 1 , x 2 ), x 1 , x 2 ∈ N, (6)
where V is an analytic function such that V (0, 0) = 0.

The discriminant ∆ of a given polynomial is a function that depends on the coefficients and it can be defined up to some constant factor as the product of the square of its roots.

In general, it can be expressed by means of the m-solutions z i = φ(x) as follows

∆(x) = i<j (φ i -φ j ) .
The discriminant of a quadratic Weierstrass polynomial is given by ∆(x) = w 2 1 (x) -4w 0 (x) (see, for instance, [START_REF] Wall | Singular points of plane curves[END_REF]). Now, the existence of parametric equation that satisfies the equation W = 0 for a quasi-ordinary polynomials, i.e., the existence of Puiseux series solutions in the multiparameter case is given by the following result (see, for instance [START_REF] Zurro | The abhyankar-jung theorem revisited[END_REF]; [START_REF] Lipman | Quasi-ordinary singularities of surfaces in C 3[END_REF]). Theorem 2. (Abhyankar-Jung Theorem). Let f ∈ C[z] be a quasi-ordinary polynomial in z with analytic coefficients in (x 1 , x 2 ). Then, there exist a natural number r such that the roots of f are given by convergent fractional power series that belong to C{x

1/r 1 , x 1/r 2 }.

Regularity Condition

Consider now the quasi-polynomials f (s, τ 1 , τ 2 ) described in (1), with a multiple root s = 0 at (τ * 1 , τ * 2 ). Furthermore, we restrict the analysis to the case when f satisfies the following assumption: Assumption 2.1. (Regularity Condition). Let s = 0 be a m-multiple of the quasi-polynomial f at (τ * 1 , τ * 2 ). Then, the following condition ∂f ∂τ i (0,0) = 0 (7)

holds for at least one τ i with i = 1, 2. Remark 2. As a consequence of the Implicit Function Theorem, if the Regularity Condition (7) holds for both τ 1 and τ 2 , then, the characteristic equation f = 0 defines the parameters (τ 1 , τ 2 ) as continuous functions of s, in some sufficiently small neighborhood of s = 0. This implies that the regularity condition presented in (7), relax the conditions imposed by the so-called non-degeneracy condition proposed by [START_REF] Irofti | Some insights into the migration of double imaginary roots under small deviation of two parameters[END_REF]. Lemma 3. If the following condition:

∂f ∂τ i s0 = 0, i = 1, 2.
holds, then the equation f (s, τ 1 , τ 2 ) = 0 is satisfied by the unique solutions τ 1 and τ 2 defined by continuous functions τ i (s 0 ) with s 0 ∈ V * δ (iω) excluding the double root s * = iω.

PROBLEM FORMULATION

The present work focuses on the analysis of the asymptotic behavior of multiple imaginary roots under multiparameter perturbations. More precisely, we consider as variable parameters the pair (τ 1 , τ 2 ).

In this vein,for the quasi-polynomial f (s, τ 1 , τ 2 ), we will focus in the following problems:

(1) derive the splitting behavior of a double root at iω under small variations of the delays (τ 1 , τ 2 ); (2) characterize conditions on the delay-parameter space (τ 1 , τ 2 ) guaranteeing the existence of Puiseux series s(τ 1 , τ 2 ) of multiple roots, together with its corresponding region of convergence, of the form

s(τ ) = c 1 τ 1/d 1 + c 2 τ 1/d 1 + c 3 τ 1/d 1 τ 1/d 1 + o(τ 1/d 1 τ 1/d n )
, where d ∈ {1, 2};

(3) characterize the existing links between the assumptions made in previous works devoted to the same problem and the condition for quasi-ordinary singularities.

Insights on Multivariate Polynomial Perturbations

In the sequel, we present some difficulties that arise when considering multi-parameter functions. In order to illustrate such arguments, we will present two examples that motivate the proposed approach.

Example 4. Consider the following polynomial:

P (z, ) = z 2 + 3 1 z + 2 2 1 + 2 2 2 , where = ( 1 , 2 ) ∈ R 2 .
Here, 1 and 2 are considered as perturbation parameters. It is clear to see, that if 1 = 2 = 0, z = 0 is a root of multiplicity two. We can see that P is a Weierstrass polynomial (3) such that the regularity condition ( 7) is not satisfied. In this case, the solutions z 1,2 ( ) are not continuous at := ( 1 , 2 ) = (0, 0). Furthermore, z 1,2 ( ) does not have a unique representation as a power series which is convergent in some punctured neighborhood of the origin. In order to illustrate this assertion on P , let us consider the case when

| 1 | < | 2 |.
In this region the solutions admit the following representation:

z 1,2 ( ) = - 1 2 (3 1 ± i4 2 ) + 1 16 1 ±i 1 2 ± i 64 1 2 3 + ± i 2048
z k ( ) = -2 k-1 1 +(-1) k 4 2 2 1 +4 2 1 3 + 32 2 1 5 + O 2 1 5
.

The previous example shows the difficulties inherent in the study of the migration of the solution for small perturbation of the parameters.

Example 5. Consider now the perturbed polynomial Q given by 8) It is easy to see that one of the roots is given by z 1 ( 1 ) = -1 , and the polynomial (z + 1 ) -1 Q(z, ) posses a double root at z 2,3 = 0 if ( ) = (0, 0). The two remaining roots can be found using the binomial expansion, with expansion given as follows

Q(z, ) = z 3 +( 1 + 1 2 )z 2 +( 1 + 2 + 2 1 2 )z+( 2 1 + 1 2 ). (
z 2,3 = ±i 1/2 2 + O 1/2 1 1/2 2
These two solutions posses the form given in (2) in the problem formulation, with c 1 = c 2 = 0 and c 3 = i. A particular case of fractional power series solution, known as quasi-ordinary singularities, can be expressed as

z 1,2 ( ) = u/2 1 v/2 2 ϕ( 1/2 1 , 1/2
2 ) such that ϕ(0, 0). This situation is guaranteed by the structure of the corresponding discriminant of Q.

MAIN RESULTS

As mention in the Introduction, our approach is based on the properties of the associated Weierstrass polynomial (4). This allows us using the algebraic properties for the root behavior analysis. Since any critical solution (s * , τ * ) can always be translated to the origin by appropriate shifts

s → s -s * , τ 1 → τ 1 -τ * 1 , τ 1 → τ 1 -τ *
1 hereinafter we will assume that (s * , τ * 1 , τ * 2 ) = (0, 0). In addition, we will consider that m ∈ N with m ≥ 2 is the algebraic multiplicity of f at (0, 0), that is,

f (0, 0) = ∂f ∂s (0,0) = • • • = ∂ m-1 f ∂z m-1 (0,0) = 0 and ∂ m f ∂z m (0,0) = 0.

Puiseux Series Solutions for Quasi-Polynomials

By means of a recursive procedure, in the commensurate delays case, Martínez-González et al. (2019a) propose a method to compute the splitting behavior of multiple root for quasi-polynomials under the variation of the delay parameter. By defining an appropriate solution surface around a multiple root, such a procedure can be extended to systems with two delay parameters. In other words, the space curve C defined by the set {(s, τ 1 , τ 2 ) ∈ C × R 2 : f = 0} can be parametrized by fractional power series in (τ 1 , τ 2 ) called Puiseux series, as in the case of one parameter delay discussed in Section 2.2. Proposition 6. Let the regularity condition (7) holds for i = 1. Then, the quasi-polynomial f admits Puiseux series solutions in τ 1 of the form

s(τ ) = c 1 (τ 2 )τ 1/m 1 + o τ 1/m
where the coefficients c k (p 2 ) can be expressed as a power fractional series in τ 2 .

The above results give some explicit representation of m-solutions which determine the solution surface. These solutions are in the form of Puiseux series, which give some insights on the splitting for a fixed τ 2 . In previous works, by using iterated Newton diagram procedure, the leading terms of the quasi-polynomial f (s, τ ) are given in an explicit manner. Now, under the regularity condition (7) the following property guarantees the existence of Puiseux series solutions in the general case. Proposition 7. Suppose that the Assumption 2.1 is satisfied for τ i with i = 1, 2. Then, the leading terms are given as

s j (τ ) = c 1 τ 1/m 1 + c 2 τ 1/β 2 + o τ 1/m , (9)
where m, β ∈ N such that β ≤ m. Proposition 8. Assume that f satisfies (7) for τ i with i ∈ {1, 2}. Let n ∈ N, such that the following partial derivatives are satisfied:

∂f ∂τ j = • • • = ∂ n-1 f ∂τ n-1 j = 0, ∂ n f ∂τ n j = 0, (10)
for n > 1, j ∈ {1, 2} and j = i. Then, the m-solutions of f are characterized as,

s k (τ 1 , τ 2 ) = τ 1/m i τ 1/ni j ϕ k (τ 1/m i , τ 1/ni j ), k = 1, 2, . . . , m,
for some n i ≤ m and ϕ k (0, 0) = 0.

Double Root of Quasi-Polynomials

The above construction does not give any details on the characterization of the root behavior in the parameter space (τ 1 , τ 2 ). Thus, to describe such a behavior, it will be shown the explicit relationship between the local form (quadratic polynomial in s) and the branches of the double root. Hence, by assuming that m = 2 and according to the previous results, the Weierstrass polynomial W (s, τ ) can be expressed by: W (s, τ ) = s 2 + w 1 (τ 1 , τ 2 )s + w 0 (τ 1 , τ 2 ). ( 11) It is possible to characterize the the root behavior by means of the Discriminant D s given by: D s (τ 1 , τ 2 ) := w 2 (τ 1 , τ 2 ) -4w 0 (τ 1 , τ 2 ). ( 12) Clearly, in the parameter space (τ 1 , τ 2 ), the condition D s = 0 guarantees the existence of a double root. Proposition 9. Let s * = iω be a double root of the quasipolynomial f (s, τ ). Then, there exists a change of variable such that the roots have the form z 2 = ϕ(τ 1 , τ 2 ), where ϕ can be expressed in the form convergent power series of order 1 in τ 1 and τ 2 .

Complete Regular Splitting

Now, to obtain information on the migration of a double root we study the properties of the branch given by the Puiseux series, but, without its explicit computation. Notice that the existence of such a series is guaranteed by Proposition 7. Hence, the migration of the double root in all cases is summarized by the following result. Proposition 10. Let W be the Weierstrass polynomial of f for the critical point (0, 0), such that D s = 0. Assume that the regularity condition ( 7) is satisfied for τ 1 and τ 2 . Then, in a neighborhood of (0, 0) the two solutions of f possess the Completely Regular Splitting Property (CRS) with respect to τ 1 and τ 2 , that is, both solutions can be expressed as:

s 1,2 (τ ) = ±c 1 τ 1/2 1 ± c 2 τ 1/2 2 + O (τ ) ,
where c 1 = 0 and c 2 = 0.

NUMERICAL EXAMPLES

In this section, we consider several numerical examples encountered in the control literature, that will allow us to illustrate the effectiveness of the proposed results. Example 11. Consider the following quasi-polynomial f (s, τ 1 , τ 2 ) = p 0 (s) + p 1 (s)e -τ1s + e -τ2s (13) where p 0 (s) := s 2 -2s + 2, p 1 (s) := 2 cos(1)s -2 (cos(1) + sin(1)) . Simple computations show that for (τ 1 , τ 2 ) = (1, 2), f has a critical root at s * = i with multiplicity two. Additionally, the first partial derivatives are given by ∂f ∂τ 1 (i,1,2) ≈ 2.91 + i0.584, ∂f ∂τ 2 (i,1,2) ≈ -0.91 + i0.416.

By applying Proposition 6, we conclude that the solutions of the Weierstrass polynomial can be expanded as a Puiseux series.

Example 12. Consider now the quasi-polynomial: f (s, τ ) = p 0 (s) + p 1 (s)e -sτ1 + p 2 (s)e -sτ2 where p 0 (s)=s 5 +s 4 + 4 + π 2 s 3 +2s 2 + 2 + π 2 s+2, (14a)

p 1 (s)=1, p 2 (s) = 2s 4 +4s 2 + 2. ( 14b
)
For (τ 1 , τ 2 ) = (π, 1), f has a double root at s = i. Let f be the quasi-polynomial derived by shifting from (i, π, 1) to the origin. Then the first non-zero partial derivatives of the quasi-polynomial at (0, 0, 0) are given by

∂ f ∂τ 1 (0,0) = i, ∂ n f ∂τ n 2 (0,0) = 0, ∀n ∈ N, (15)
∂ 2 f ∂s∂τ 1 (0,0) = 1 -iπ, ∂ n+1 f ∂s∂τ n 2 (0,0) = 0, ∀n ∈ N. ( 16
)
It is clear that the associated Weierstrass polynomial W has the following structure: W (s, τ ) := s 2 + w 1 (τ ) + w 0 (τ ). Moreover, since the regularity condition (7) holds, we only need to compute w 0 (τ ). After simple computations, one gets: w 0 (τ ) = -2i (8 + π 2 ) + i(8 -3π) + 16e -i τ 1 + o (τ ) . Thus, its solutions are given by:

s(τ ) = i ± √ 2i 3/2 √ ( 8+π
2 )+i(8-3π)+16e -i (τ 1 -π) 1/2 + O (|τ -τ * |) . Proposition 10 guarantees that the solution (i, π, 1) has the CRS property. This behavior is illustrated in Figure 2.

ℜ(s)

s0,1(τ ), τ1 < π, τ2 < 1 s0,0(τ ), τ1 < π, τ2 < 1 2. Root locus of quasi-polynomial f (s, τ ) ( 14) around (i, π, 1).

Example 13. Consider now the special case of a model of population dynamics discussed by [START_REF] Irofti | Some insights into the migration of double imaginary roots under small deviation of two parameters[END_REF]:

f (s, τ 1 , τ 2 ) = s -a -b 1 -e -sτ1 s -c 1 -e -sτ2 s . ( 17
)
Taking the same values as the ones proposed by [START_REF] Irofti | Some insights into the migration of double imaginary roots under small deviation of two parameters[END_REF], that is, a = -0.214104, b = -0.996801 and c = 0.5, we know that for (τ * 1 , τ * 2 ) ≈ (3.84003026849, 10.44866732901), f possesses a double root at s * = i. In this case, we have the following partial derivatives: ∂f ∂τ 1 (i,τ * ) = 0, ∂f ∂τ 2 (i,τ * ) = 0.

Thus, the regularity condition (7) holds for τ 1 and τ 2 , and consequently, by using Proposition 6, its solutions can be expanded as Puiseux series with fractional power of 1 2 . Furthermore, by Proposition 7, its solutions behave as: s(τ ) = i ± (0.3885 -0.3307i)(τ 2 -τ * 2 ) 1/2 + O (|τ -τ * |) .

CONCLUDING REMARKS

The paper addresses the asymptotic behavior of multiple critical roots for quasi-polynomials of retarded-type with two incommensurate delays. The proposed approach allows identifying the structure of the solutions by only computing some partial derivatives. Thus, such results reveal if the solutions behave as a power series or as a Pusieux series. Besides, the adopted approach allows to relax some conditions imposed on the quasi-polynomial in some results previously reported in the literature, and as a consequence, our method allows capturing the behavior of any quasi-polynomial of retarded type with algebraic multiplicity two.

  instead of the previous region, we consider the region | 2 | < | 1 |. Then, for k ∈ {1, 2} the solutions admit the following representation

a direct consequence of the Rouché's lemma, see, e.g.,[START_REF] Michiels | Stability, Control, and Computation for Time-Delay Systems. An Eigenvalue-Based Approach[END_REF]

Such roots are simply called critical characteristic roots.

In fact, in the case when there are no characteristic roots on the imaginary axis, the stability/instability of the dynamical system free of delays is preserved for any delay value, i.e. the delay-independent stability/instability property, see, for instance,[START_REF] Niculescu | Delay Effects on Stability: A Robust Control Approach[END_REF];[START_REF] Gu | Stability of time-delay systems[END_REF];[START_REF] Michiels | Stability, Control, and Computation for Time-Delay Systems. An Eigenvalue-Based Approach[END_REF] for further arguments.

ACKNOWLEDGEMENTS

The research of the first author is financially supported by CONACyT, Mexico and partially supported by Laboratory of Signals and Systems (University Paris-Saclay, CNRS, CentraleSupélec), France.


Cédric Escudero
email: cedric.escudero@grenoble-inp.fr

Paolo Massioni
email: paolo.massioni@insa-lyon.fr

Gérard Scorletti
email: gerard.scorletti@ec-lyon.fr

Eric Zamaï
email: eric.zamai@insa-lyon.fr

Security of Control Systems: Prevention of Aging Attacks by means of Convex Robust Simulation Forecasts

Keywords: Aging attacks, secure control theory, Lyapunov methods, sum of squares, linear matrix inequalities

Aging attacks are a class of attacks targeting closed-loop control systems connected to a network. Such attacks consist in slightly modifying the control signals to increase the wear and tear of the physical system while maintaining the delivered service. In this paper, we exploit the technique of robust simulation of linear systems for forecasting and preventing these attacks by means of convex optimization. Besides, this technique allows the computation of the time at which the normal aging is not guaranteed anymore for a given control input set. The paper end highlights an example application upon an electrical machine.

INTRODUCTION

Modern control systems consist of digital control components, such as controllers and networks, for monitoring and controlling a physical system. With their digitization, they are now suffering from vulnerabilities like any computerbased system. As a result, modern control systems have been the target of attackers [START_REF] Franck | From ICS attacks' analysis to the S.A.F.E. approach: Implementation of filters based on behavioral models and critical state distance for ICS cybersecurity[END_REF]. So, they aim at taking control over the physical system in order to damage it or degrade the control system's performances which might lead to damages. Over the past decade, the physical system integrity has became a major concern for the security of control systems [START_REF] Escudero | Processaware model based IDSs for industrial control systems cybersecurity: Approaches, limits and further research[END_REF]. In particular, deception attacks, a type of attacks, have gained attention. They consist in altering the behavior of the physical system by compromising one or multiple control components (e.g. controller, network) (Teixeira et al., 2012).

In the secure control literature, prevention and detection of deception attacks have been investigated from the standpoint of the control system's performances [START_REF] Giraldo | A survey of physics-based attack detection in cyber-physical systems[END_REF]. In other words, they are interested in attacks that alter the control system's objective which is related to the services delivered by the physical system. [START_REF] Mo | Detecting integrity attacks on scada systems[END_REF] investigate the detection and the feasibility of replay attacks in which an attacker manipulates the control signals while replaying a sequence of recorded sensor measurements. [START_REF] Bai | Security in stochastic control systems: Fundamental limitations and performance bounds[END_REF] quantify the maximum degradation of the control system's performances an attacker can induce by hijacking and replacing the control signal in the presence of an anomaly detector. With a similar attacker model, Teixeira et al. (2012) exploit the unobservable states in order to mislead the anomaly detector and the controller about the physical system's states. [START_REF] Murguia | Reachable sets of hidden cps sensor attacks: Analysis and synthesis tools[END_REF] analyze the consequences of stealthy attacks over the control system's performances by the means of positively invariant sets. Similarly, Hadizadeh Kafash et al. (2018) propose to limit the positively invariant set to satisfy properties of the physical system's services.

In this work context, we focus on a class of attacks that we call "aging attacks". Their goal is mainly to reduce the physical system's lifetime. They consist in slightly modifying the control signals reaching the actuators to maximize the wear and tear of the machinery, subsequently reducing its lifetime, causing increased corrective maintenance and potential dangers to the attacked party. At the same time, the attacker will make sure that the control system's objective is not degraded for the stealthiness purpose. As aging attack examples, the temperature increase in an electrical machine deteriorates earlier its insulation system which leads further to a failure. Also, the power consumption increase in an autonomous system leads to a fast discharge of the energy storage which ceases its operations.

Aging attacks act on dynamical systems, which makes their forecast and prevention amenable to be dealt with using the methods of systems and control theory. Previous works have been conducted in [START_REF] Escudero | Prevention of aging attacks: Malicious nature of the control signal[END_REF] to limit the available actuator output range by the means of positively invariant sets. As a continuation of these works, we propose in this paper to use robust convex simulation which is less conservative and allows to consider the time dimension. The methodology is based on the ability of forecasting the possible reachable states of the closed-loop plant under an external attack, which allows tuning the actuators range in order to prevent them. The proposed forecasting method is based on Pseudo-Lyapunov function for determining ellipsoidal bounding sets that capture the state trajectory of a dynamical system over a given time interval.

The rest of the paper is organized as follows. Section 2 presents the problem, the attacker model and the main theoretical background. Section 3 gives the main theorem which is sufficient conditions for computing ellipsoidal bounding set over a time interval. Section 4 presents the algorithms for limiting the available actuator output range such that aging attacks cannot be launched. Section 5 highlights our main contributions for preventing attacks on an example application of a brushed DC motor. Finally, concluding remarks and directions are given in Section 6.

PRELIMINARIES

Notation

Let R be the set of real numbers, and R n×m be the set of real n × m matrices. Given a vector v ∈ R n , v denotes its transpose, and [v] i is the i th element of v. Given a matrix A, A indicates its transpose, diag(a 1 , a 2 , ..., a n ) denotes the diagonal matrix with the diagonal elements a 1 , a 2 , ..., a n , I n is the identity matrix of size n, 0 n,m is the zero matrix of size n, m. The notation A 0 (A 0) indicates the positive (negative) semidefiniteness -i.e. all the eigenvalues of the symmetric matrix A are positive (negative) or equal to zero-, whereas A 0 (A ≺ 0) indicates the positive (negative) definitenessi.e. eigenvalues are strictly positive (negative).-We also define E v (A, v) as the ellipsoid of dimension n with matrix

A ∈ R n×n , A = A 0 and centered in v ∈ R n -i.e. E v (A, v) = {v ∈ R n | (v -v) A(v -v) 1}.-In addition, let H v (c, µ) be an hyperplane defined as H v (c, µ) = {v ∈ R n | c v = µ}. Finally, let R m [t]
be the set of polynomials of degree up to m in the variable t, and R n×n m [t] be the set of symmetric matrix-valued polynomials of degree up to m in the variable t of size n × n, with m ∈ N.

Physical system: In this paper, we consider dynamical systems modelled as affine time-invariant continuous-time systems as follows:

ẋ(t) = Ax(t) + Bu(t) + a (1) where x(t) ∈ R n is the state vector, u(t) ∈ R m is the input vector, A ∈ R n×n , B ∈ R n×m , and a ∈ R n .
The model in (1) above can be simplified by defining the extended state vector x(t) = [x (t), 1] , which allows reformulating it as:

ẋ(t) = Ãx(t) + Bu(t) (2) with à = A a 0 1,n 0 , B = B 0 1,m . (3)
With this notation, an ellipsoid n+1) .

E x (Q, x) = {x ∈ R n | (x - x) Q(x -x) 1} is equivalently defined by the set {x = [x , 1] | x ∈ R n , x Qx 1}, where Q = Q -Qx -x Q x Qx ∈ R (n+1)×(
(4)

Remark 1. Without loss of generality, we consider that in the model in (2) it is possible to distinguish the states related to the service delivered by the physical system x g , in other words the states involved into the control system's objective, from the states x s related to the stress of the physical system.

Anomaly detector : Let an output feedback controller satisfy a control system's objective for the physical system's services in (2) such that (ū(t), xg (t)) is a nominal trajectory for the system in (2). In addition, it is reasonable to assume the presence of sensors for monitoring the service such that the controller can comply with the control system's objective. This set of measurement is captured by y g (t) ∈ R l with y g (t) = Cx g (t).

In addition, we consider an anomaly detector that monitors deviations from a nominal trajectory at each time instant from the component-wise absolute value:

r g (t) = |ŷ g (t) -y g (t)| (5)
where ŷg (t) is the estimated measurement vector, and r g (t) ∈ R l is the residue vector with each of its component is evaluated to a threshold τ R l for determining the presence of an anomaly as in (6):

r g (t) τ. (6)

Attacker model and problem formulation

Regarding the attack space formulated in Teixeira et al. (2012), we define the attacker model for aging attacks as follows. The attacker is able to inject additional control data in the actuator channel by compromising either the controller or the actuator channel itself which is captured by having u(t) = ū(t)+u a (t). However, the attacker cannot eavesdrop on the sensor and actuator data. In other words, the attack policy assumes no online data. In addition, the attack is computed a priori corresponding to an openloop type of attack policy. The attacker has access to the physical system model in (2) and the anomaly detector model in ( 5), (6).

Therefore, the attacker's objective consists in searching for a control signal u a (t) that maximizes the stress states x s (t) while guaranteeing the service residue r g (t) remains below the anomaly detector's threshold τ (stealthiness).

In complement, the defender's objective is the prevention of aging attacks formulated as searching for the upper bound of u(t) that maximizes the service states x g (t) while guaranteeing the stress states remains below a certain theshold corresponding to a normal aging behavior.

Clearly, adding additional sensors for monitoring the stress states might be considered, but will inherently increase the attack surfaces allowing an attacker taking advantage of them. That is why we propose to prevent the aging attacks by restricting the control input u(t) received by the physical system.

The rest of the section provides the main theoretical background for understanding the results of this paper.

Sum of squares

In this paper, we are going to rely on a technique called Sum Of Squares (SOS), which allows casting several classes of polynomial problems into convex optimization problems. We recall here briefly the basic notions that are necessary for understanding the paper, the interested reader can then consult the references that are provided. Definition 1. (SOS problems, [START_REF] Parrilo | Semidefinite programming relaxations for semialgebraic problems[END_REF]). Let p(t) ∈ R 2d [t]; we call Sum Of Squares problem (SOS) the problem of finding whether there exists a finite number l of polynomials

π i (t) ∈ R d [t] such that p(t) = l i=1 π i (t) 2 . ( 7
)
If such an expression exists, then p(t) is a sum of squares (SOS), which implies that p(t) 0 for all t.

The class above can be extended to the class of Matrix SOS problems (MSOS). Definition 2. (MSOS problems, Chesi ( 2010)). Let P(t) ∈ R n×n 2d [t]; we call a Matrix Sum Of Squares problem (MSOS) the problem of finding whether there exists a finite number l of matrices of polynomials Π i (t) such that

P(t) = l i=1 Π i (t) Π i (t). (8)
If such a decomposition exists, then P (t) is a matrix sum of squares (MSOS), which implies that P (t) 0 for all t.

The definitions above are restricted to univariate polynomials, as this is the case we consider in this paper. SOS and MSOS problems are convex problems. A derived class of problems is that of feasibility problems under SOS or MSOS constraints. Definition 3. (SOS constraints feasibility, [START_REF] Parrilo | Semidefinite programming relaxations for semialgebraic problems[END_REF]).

Let

p i (t, θ) ∈ R 2d [t]
, ∀i ∈ {1, 2, ..., q}, where θ ∈ R ρ is a vector of parameters or unknowns, with p i (t, θ) affine with respect to the entries of θ. A feasibility problem under SOS constraints consists in finding, if it exists, a value of θ = θ * for which p i (t, θ * ) is SOS, for i = 1, ..., q.

(9) If such a θ * exists, then the problem is feasible; otherwise it is unfeasible. Definition 4. (MSOS constraints feasibility, Chesi (2010)).

Let P i (t, θ) ∈ R n×n 2d [t], ∀i ∈ {1, 2, ..., q}
, where θ ∈ R ρ is a vector of parameters or unknowns, and the matrices P i (t, θ) are affine with respect to the entries of θ. A feasibility problem under SOS constraints consists in finding, if it exists, a value of θ = θ * for which P i (t, θ * ) is MSOS, for i = 1, ..., q.

(10) If such a θ * exists, then the problem is feasible; otherwise it is unfeasible.

Feasibility problems under SOS and MSOS problems are convex optimization problems, and they can be reformulated as linear matrix inequality (LMI) feasibility problems; this can either be done explicitly, or relying on automated procedures, like the one available in the Yalmip toolbox [START_REF] Löfberg | Pre-and post-processing sum-ofsquares programs in practice[END_REF] under Matlab. Minimising any single affine function of the decision variables under SOS and/or MSOS constraints is also a convex optimization problem.

The generalised S-procedure

The S-procedure allows restricting some classes of inequalities to a certain given subset of the variables that are concerned [START_REF] Boyd | Linear matrix inequalities in system and control theory[END_REF]. In this paper, we rely on a general expression which can be specialised according to the cases. Lemma 1. (Generalised S-procedure). Let F (x), G(x) be (symmetric matrix) functions of the (vector) variable x. Let g(x) be a scalar function of x. The following implications hold.

F (x) 0 for G(x) 0 ⇐ F (x) -G(x)λ 0 for λ 0, ∀x (11) F (x) 0 for g(x) 0 ⇐ F (x) -g(x)Λ 0 for Λ 0, ∀x (12)
The terms Λ and λ are called multipliers, which can be chosen at one's convenience, i.e. they are decision variables subject to the positivity constraints above. When used in the context of polynomial problems as in this case, the lemma above is a direct consequence of a lemma known as Positivstellensatz [START_REF] Chesi | LMI techniques for optimization over polynomials in control: a survey[END_REF] or p-satz, of which several versions exist in the literature. The multipliers can in this case have a polynomial dependence on x, which allows satisfying the positivity constraints by means of either SOS or MSOS constraints.

Robust simulation and Pseudo-Lyapunov functions

By robust simulation [START_REF] Kantner | Robust simulation and nonlinear performance[END_REF] we mean simulation of a dynamical system for a whole set of initial conditions, under a certain number of constraints. Namely, the problem of robust simulation consists in finding a bounding set for the state x(t) at a final time t f , given a set of possible initial values for it at an initial time t 0 ; this under the hypotheses of a given dynamical equation and in the presence of additional constraints of different kind. In this work, we consider ellipsoids as bounding sets.

In the context of this work, we adapt the ideas in Ben-Talha et al. ( 2017); [START_REF] Tobenkin | Invariant funnels around trajectories using sumof-squares programming[END_REF], that are based on the search of a time-dependent pseudo-Lyapunov function. The idea is to define a positive definite function

V (x(t), t) = x(t) Q(t)x(t) > 0 ∀x(t) = 0, (13)
such that V (x(t), t) 0, where Q(t) is a time-varying symmetric matrix. Then the following lemma holds, under the hypothesis of u(t) = 0. Lemma 2. Consider the system in (2), with initial conditions at t 0 satisfying x(t 0 ) Qt0 x(t 0 ) 1. If there exists a matrix function Q

(t) = Q(t) ∈ R n×n , with Q(t 0 ) = Qt0 such that Q(t) + Ã Q(t) + Q(t) Ã 0 for t ∈ [t 0 , t f ], ( 14
) then x(t) Q(t)x(t) 1 for t ∈ [t 0 , t f ].
The lemma above can be interpreted as follows: x(t 0 ) is inside the level curve of value 1 of the function V (x(t), t).

The inequality in ( 14) implies that for all trajectories of the state according to (2), the function can only decrease; so those trajectories are bounded to stay inside the same level curve at all times t ∈ [t 0 , t f ].

MAIN THEOREM

The main theoretical result of this paper consists in a theorem providing sufficient conditions for computing the ellipsoidal bounding set over a time interval [t 0 , t f ] for a system as in (1). This theorem will be subsequently employed in a set of algorithms, whose goal is the prediction and prevention of potential aging attacks.

The main theorem is obtained by generalising Lemma 2 in order to account for the presence of a bounded input signal. In fact, we will assume

u(t) ∈ E u (R, ū), i.e. (u -ū) R(u -ū) 1 ( 15
)
with R ∈ R m×m a symmetric matrix, for which we assume R 0.

In the presence of input, the time derivative of the pseudo-Lyapunov function V (x(t), t) defined in ( 13) is not ( 14), but it becomes:

V (x(t), t) = x(t) u(t) Q(t) + Ã Q(t) + Q(t) Ã Q(t) B B Q(t) 0 m,m
x(t) u(t) .

(16) Notice that if V (x(t), t) 0, then V (x(t), t) defines timevariant ellipsoidal level sets from which the state trajectory will never go out over the given time interval; so if a set of initial conditions is within a level curve, the trajectories stemming from them will never exit the level curve of the same value. The time-variant Q(t) allows the ellipsoidal level curves to move and grow or shrink over time with the state trajectory, whether they are converging or not.

In order to make the search tractable, we are going to limit it to the set of matrix-valued polynomials, i.e. we will set Q(t) ∈ R (n+1)×(n+1) 2d

[t] for an arbitrarily chosen integer d 1. It is clear that specifying Q(t) in this way adds some conservatism, which can be progressively reduced by increasing the degree of the matrix-valued polynomial. Q(t) is then its first time-derivative, with

Q ∈ R (n+1)×(n+1) 2d-1 [t].
As the attacker wants to remain stealthy regarding the anomaly detector, we may assume that the control signals will be chosen carefully in order to not degrade the service delivered by the physical system. For this reason, we can assume that a part of the state can be constrained to belong to a given specific set, the set that will not cause any suspects of attack to arise. This ellipsoid will be defined as E x (Ξ, ξ), where Ξ in general will be rank-deficient (it will basically only constrain x g ), as it will constrain only one part of the state; E x (Ξ, ξ) can even coincide with R n×n (formally by picking Ξ = 0). As a consequence of this, we formulate the following assumption: Assumption 1. There exists an input signal u(t) ∈ E u (R, ū) such that x(t) ∈ E x (Ξ, ξ) for any t ∈ [t 0 , t f ] and for initial conditions x(t 0 ) ∈ E x (Q(t 0 ), x(t 0 )) .

If the assumption above is not satisfied, it means that no aging attack can be performed without being detected, which put us on the safe side.

Finally, just before formulating our main result, we define three additional terms which will be useful later on. First, it will be necessary at some steps to extract the term 1 from x, for which we define Γ = [0 1,n 1] (17) such that Γx = 1. Subsequently, we can also define γ(t) = (t f -t)(t -t 0 ), (18) which allows expressing the constraint of t ∈ [t 0 , t f ] equivalently as simply γ(t) 0. Last, consider

Ĩ = I n 0 n,1 0 1,n 0 , Ξ = Ξ 0 n,1 0 1,n 0, , ξ = ξ 0 . (19)
We can now formulate the main theorem.

Theorem 1. Consider an affine time-invariant continuoustime system according to (1) or equivalently (2), with matrices defined according to (3).

If for a given d

1 there exist Q(t) ∈ R (n+1)×(n+1) 2d [t], Z(t), X(t) ∈ R (n+1)×(n+1) 2d-2 [t] α(t), β(t) ∈ R 2d-2 [t]
, and a scalar ε > 0 for which the following constraints are satisfied:

-

M (t) -α(t)S -β(t)T -γ(t)Z(t) is MSOS, (20) Q(t) -γ(t)X(t) -ε Ĩ is MSOS (21) Z(t), X(t) are MSOS (22) α(t), β(t) are SOS (23) with M (t) = Q(t) + Ã Q(t) + Q(t) Ã Q(t) B B Q(t) 0 m,m S = Γ Γ -Γ ū RūΓ Γ ū R RūΓ -R T = Γ Γ -Ξ + Ξ ξΓ + Γ ξ Ξ -Γ ξ Ξ ξΓ 0 n+1,m 0 m,n+1 0 m,m
and Γ, γ(t), Ĩ defined as in ( 17), ( 18), ( 19), then

x(t 0 ) Q(t 0 )x(t 0 ) 1 ⇒ x(t) Q(t)x(t) 1 ∀t ∈ [t 0 , t f ], (24) under the constraints ∀t ∈ [t 0 , t f ]u(t) ∈ E u (R, ū), x(t) ∈ E x (Ξ, ξ).
Proof. Consider first (21); left and right multiply by x(t), and consider γ(t)X(t) as an S-procedure term by positive multiplier X(t) as in ( 22); this implies (Lemma 1): Subsequently,consider (20); left and right multiply by [x(t) , u(t) ] , and consider α(t)S, β(t)T and γ(t)Z(t) as S-procedure terms by positive multipliers α(t), β(t) and Z(t) as in ( 22) and ( 23); this implies (Lemma 1):

x(t) Q(t)x(t) ε||x(t)|| 2 (25) when γ(t) 0, i.e. V (x(t), t) = x(t) Q(t)x(t) > 0 when x(t) = 0, t ∈ [t 0 , t f ].
[x(t) , u(t) ] M (t) [x(t) , u(t) ] = V (x(t), t) 0 (26) when [x(t) , u(t) ] S [x(t) , u(t) ] 0 ⇔ u(t) ∈ E u (R, ū). ( 27
) [x(t) , u(t) ] T [x(t) , u(t) ] 0 ⇔ x(t) ∈ E x (Ξ, ξ), (28
) and when γ(t) 0, i.e. when t ∈ [t 0 , t f ]. This means that the value of V (x(t), t) can only increase under the stated constraints, i.e. V (x(t 0 ), t 0 ) 1 ⇒ V (x(t), t) 1 in the considered interval, which is the theorem statement. Remark 2. Although SOS problems are exact for the univariate case, Theorem 1 is conservative, i.e. it provides sufficient but not necessary conditions, for two main reasons: 1) the pseudo-Lyapunov function is constrained to be of a specific form, and 2) the use of the generalized S-procedure for obtaining (20). In both cases the conservatism can be reduced by increasing d, which provides more degrees of freedom to the function, as well as it reduces the conservatism of the S-procedure according to the applicable formulation of the p-satz [START_REF] Lasserre | An introduction to polynomial and semi-algebraic optimization[END_REF]. Notice also that the polynomial degrees of Q(t), Z(t), X(t), α(t), and β(t) have all been chosen in order to have all terms in (20) of the same degree; in fact, the overall cost of solving the MSOS problem depends on the overall degree in t of the expression in ( 20), so it is computationally inefficient to have terms of the sum of higher degree with respect to the others. In other words, the conservatism is minimum for a given computational cost when all terms have the same degree.

ALGORITHMS FOR ATTACK ANALYSIS AND PREVENTION

In this section, we are proposing algorithms involving Theorem 1 for preventing aging attacks over a given time interval. In other words, we want to restrain the input set

E u (R, ū) into the set E u ( R, ū) (E u ( R, ū) ⊆ E u (R, ū)) such that a dangerous set D x is avoided over a time interval [t 0 , t f ].
Firstly, we are interested in finding the minimal bouding ellipsoid for the system defined in (1) over the time interval [t 0 , t f ]. In fact, there exists many bounding ellipsoids under Theorem 1; however, we want to find the tighest ellipsoid set that encapsulates the state trajectory over the time interval. This can be obtained by maximizing the trace of Q(t f ) under Theorem 1 with initial conditions (Q(t 0 ),x(t 0 )) as stated in P 1 . In order to give more degrees of freedom to the search of such ellipsoids and without loss of generality, let t 0 be time 0 and split the final time t f in N ∈ N time steps such that N = t f ∆ -1 t with ∆ t a small time step, which simplifies the computations. Hence, Algorithm 1 is proposed for finding the minimal bounding ellipsoid for system in (2) with a control input set E u (R, ū)) over the time interval [0, ∆ t ] by using Theorem 1.

P 1 : maximize Qt,Z(t),X(t),α(t),β(t) trace( Q∆t )
subject to (20), ( 21), ( 22), (23)

Notice that the numerical procedure will have a limited numerical precision, for this reason a normalisation step (step 3.) has been introduced; typically the ellipsoids matrices might not fit into the form of ( 4), with a nonmatching lower-right entry that corresponds to defining a set {v ∈ R n | (v-v) A nn (v-v) ρ} with ρ = 1. The function Normalize( Ãnn ) transforms a non-normalized (nn) el-

Algorithm 1 : [Q ∆t , x∆t ]=Alg1(Q(t 0 ),x(t 0 ),R,ū,Ξ,ξ, I 1 ) Input: Q(t 0 ),x(t 0 ),R,ū, Ξ, ξ, I 1 :={A, B, ∆ t , d} Initial condition: Q 0 = Q(t 0 ), x0 = x(t 0 ) 1. Set Q0 = Q 0 -Q 0 x0 -x 0 Q 0 x 0 Q 0 x0 , Q(t) = 2d i=0 Qi t i 2. Solve P 1 3. [Q(∆ t ), x(∆ t )]=Normalize( Q(∆ t )) Output: Q(∆ t ), x(∆ t ) lipsoid defined as E nn v (A nn , v) = {v ∈ R n | (v-v) A nn (v- v) ρ} into a normalized one E v (A, v) (ρ = 1).
Secondly, we propose Algorithm 2, whose goal is to determine if and at which time t ∩ (with t ∩ ∈ {0, ∆ t , ..., N ∆ t }) the bounding ellipsoid hits a dangerous set D x . We consider the dangerous set as a set where the lifetime for the given system is reduced.

Consider the dangerous set as the union of i number of halfspaces H i x (c i , µ i ) defined by their boundary hyperplane in (29):

D x = {x ∈ R n : p i=1 c i x(t) ≥ µ i } (29)
where c i ∈ R n is an affine combination of states in x, µ i ∈ R is the bound. Let us formulate the overlapping distance from an ellipsoid E v (A, v) to an hyperplane [START_REF] Kurzhanskiy | Ellipsoidal Toolbox[END_REF] as follows:

H i v (c i , µ i ) denoted dist(E v (A, v), H i v (c i , µ i )) detailed in
dist(E v (A, v), H i v (c i , µ i )) = (c i A -1 c i ) 1/2 -|µ i -c i v| (c i c i ) 1/2
(30) where dist(E v (A, v), H v (c, µ)) 0 if the ellipsoid hits the hyperplane -i.e. the ellipsoid intersects the hyperplane-; otherwise it is negative. Hence, for determining if an ellipsoid E v (A, v) hits a dangerous set D x in (29), it is sufficient to verify that the distance from the ellipsoid to each hyperplane is negative, leading to the boolean function Hit(A,v,H i v ).

Algorithm 2 : [hit, t ∩ ]=Alg2(Q(t 0 ),x(t 0 ),R,ū,Ξ,ξ, I 1 ,I 2 ) Input: Q(t 0 ),x(t 0 ),R,ū,Ξ,ξ, I 1 ,I 2 := {H i x , t f } Initial condition: Q 0 = Q(t 0 ), x0 = x(t 0 ) 1. Set N = t f ∆ -1 t 2. For j = 1 to N 2.1. [Q ∆t , x∆t ]=Alg1(Q 0 ,x 0 ,R,ū,Ξ,ξ, I 1 ,I 2 ) 2.2. If Hit(Q ∆t ,x ∆t ,H i x ) = T rue Set hit = T rue, t ∩ = j∆ t , Goto Output 2.3. Set Q 0 = Q ∆t , x0 = x∆t 3. Set hit = F alse, t ∩ = ∅ Output: hit, t ∩
Finally, Algorithm 3 aims at restraining the input set to the set E u ( R, ū) such that the given dangerous set is guaranteed to be unreachable by the state trajectory for system in (1) over the time interval [0,t f ]. The algorithm implements a bisection method on the scaling term δ for R for solving the quasi-convex problem. Note that higher the trace of R is, tighter the constraints on u(t) are.

Algorithm 3 converges, with respect to a tolerance tol, to the minimal restrained input set E u ( R, ū) with R = δR 0 for a sufficiently large R 0 such that the dangerous set is not reachable by the state trajectory for E u (R 0 , ū), and a δ ∈ [0, 1] sufficiently small such that the state trajectory can reach the dangerous set for E u (δR 0 , ū).

Algorithm 3 :[diag, δ wrk ]= Alg3(R 0 ,δ,tol,Q(t 0 ),x(t 0 ),ū,Ξ,ξ, I 1 ,I 2 ) Input: R 0 ,δ,tol,Q(t 0 ),x(t 0 ),ū,Ξ,ξ, I 1 ,I 2 Initial condition: diag = N otF ound, Q 0 = Q(t 0 ), x0 = x(t 0 ) 1. Set δ = 1 2. [hit, t ∩ ]=Alg2(Q(t 0 ),x(t 0 ),δR 0 ,ū,Ξ,ξ,I 1 ,I 2 ) 3. If hit = F alse 3.1. Set diag = N everHit, δ = δ, Goto Output 4. [hit, t ∩ ]=Alg2(Q(t 0 ),x(t 0 ),δR 0 ,ū,Ξ,ξ,I 1 ,I 2 ) 5. If hit = T rue 5.1 Set diag = AlwaysHit, δ = ∅, Goto Output 6. Set δ = δ 7. While (δ -δ > tol) 7.1. δ it = (δ + δ)/2 7.2. R = δ it R 0 7.3. [hit, t ∩ ]=Alg2(Q(t 0 ),x(t 0 ),R,ū,Ξ,ξ,I 1 ,I 2 ) 7.4. If hit = T rue, Set δ = δ it 7.4. Else, Set δ = δ it , diag = F ound, δ = δ it Output: diag, δ

APPLICATION

In this section, we propose to apply the proposed algorithms for a brushed DC motor by first analyzing the potential aging an attacker can induce from the power supply limits, and lastly preventing them.

Description of the system

Consider a brushed DC motor from [START_REF] Salah | Parameters identification of a permanent magnet DC motor[END_REF] with a constant load torque T l = 0.5 given in (31) controlled by an output feedback controller issuing a voltage u(t) through a network to its armature for achieving a control objective on its angular velocity ω(t) (service); consider also a single aging factor, which is the winding temperature increase θ(t). Regarding the lifetime, we assume a normal aging behavior given for θ(t) ≤ θ n . For the sake of simplicity, we consider θ(t) = k DC i 2 (t) with k DC = 0.34 the DC gain of the thermal model which leads to -i n i(t) i n with i n = 21.6 A for guaranteeing the lifetime. Hence, the dangerous set D x is defined for H 1

x (c 1 , µ 1 ) and H 2 x (c 2 , µ 2 ) with c 1 = [0, 1] , µ 1 = 21.6, c 2 = [0, -1] , µ 2 = -21.6. Recall that the control signal reaching the brushed DC motor is expressed as u(t) = ū(t) + u a (t) with ū(t) the nominal control input and u a (t) the attack signal. The state-space matrices of the system are: In both the coming example applications, we have considered Q(t) MSOS of degree 4 (d = 2), Z(t), X(t) a matrix (MSOS of degree 0), and α(t) SOS of degree 4, and β(t) an unknown scalar (SOS of degree 0).

A = -0.

Attack analysis: forecast simulation

In this subsection, we analyze the potential aging an attacker can induce to the brushed DC motor from the power supply limits. The control input set can be defined as R = 1 u 2 bnd which is maximal (R max ) -i.e. in the sense of its size-for u bnd = 110 V and ū = 0 from the power supply limits. From this maximal control input set, we want to verify if the physical integrity of the brushed DC motor -i.e. aging-can be transgressed over the time interval [t 0 ,t f ]. Instead of stopping the simulation once a bounding ellipsoid at time t = j∆ t in Algorithm 2 hits with the dangerous set, we stop it once the bounding ellipsoids converge, in this case for t f = 1s. We apply Algorithm 2 with the inputs: R = 1 110 2 , ū = 0, Q(t 0 ) = diag(10 5 , 10 7 ), x(t 0 ) = [0, 0] , Ξ = 0, ξ = [0, 0] , ∆ t = 0.1 ms, t f = 1 s. The result is graphically shown in Fig. 1 for each bounding ellipsoids at each 10 ms. Fig. 1. Bounding ellipsoids at each 10 ms for t f = 1 s from initial conditions E x (diag(10 5 , 10 7 ), [0, 0] ) (blue star) and E u ( 1 110 2 , 0); final bounding ellipsoid (filled ellipsoid in red) ; boundaries of the dangerous set (red lines).

The bounding ellipsoid hits for the first time the dangerous set at t ∩ = 34.9 ms. From this analysis, we can note that the physical system can reach the dangerous set over the given time interval [0,1s]. As a result, an attacker is capable to induce an anomalous aging to the brushed DC motor by injecting an attack signal u a (t).

Prevention of attacks

In this subsection, we want to restrain the control input set such that the dangerous set can never be hit. We apply Algorithm 3 on the previous example for finding the maximal restrained control input set such that none of the bounding ellipsoids hits the dangerous set over the given time interval. The inputs are: R 0 = 1 30 2 , δ = 0.01, tol = 0.01, Q(t 0 ) = diag(10 5 , 10 7 ), x(t 0 ) = [0, 0] , ū = 0, Ξ = 0, ξ = [0, 0] , ∆ t = 0.1 ms, t f = 1 s. The restrained control input set E u ( R, ū) is given from the output δ = 0.6674 of Algorithm 3 and computed as R = δR 0 = 1 36.72 2 . The resulting bounding ellipsoids over the time interval [0,1 s] for the restrained control input set E u ( R, ū) computed in Algorithm 3 are shown for each 10ms in Fig. 2. Fig. 2. Bounding ellipsoids at each 10 ms for t f = 1 s from initial conditions E x (diag(10 5 , 10 7 ), [0, 0] ) (blue star) and E u ( 1 36.72 2 , 0); final bounding ellipsoid (filled ellipsoid in red) ; boundaries of the dangerous set (red lines).

As we can see in Fig. 2, the dangerous set is guaranteed unreachable for u(t) ∈ E u ( 1 36.72 2 , 0). Therefore, aging attacks are prevented by restricting the control input u(t) to the restrained input set.

In Fig. 3, it is represented the bounding ellipsoids evolution with respect to the time t ∈ [0, 1 s] for the maximal control input set E u ( 1 110 2 , 0) computed in the previous subsection (transparency) and the restrained control input set E u ( 1 36.72 2 , 0) (non transparent). Secondly, we want to prevent stealthy aging attacks that might occur while the system trajectory is constant and given for (ū,x) with ū = 20 and x = [10.07, 0.5] . The attacker aims at remaining stealthy regarding the anomaly detector in (5) tuned at τ = 20 by ensuring the states related to the service lie into E x (Ξ, ξ) with Ξ = diag( 1 τ 2 , 0) and ξ = x. Applying Algorithm 3 with R0 = 1 1 2 , δ = 0.01, tol = 0.1, and Q(t 0 ) = diag(10 5 , 10 7 ), x(t 0 ) = x, it yields the restrained input set E u ( 1 3.1 2 , 20). The resulting bounding ellipsoids at each 20 ms are given in Fig. 4. As we can observe in Fig. 4, the bounding ellipsoids never hit the dangerous set for the computed restrained input set. This ensures for the constant nominal trajectory to prevent aging attacks.

Finally, we compute the evolution of the time at which the normal aging is not guaranteed anymore with respect to the evolution of the restrained input set. In other words, we compute the hit time t ∩ in function of δu bnd with δ ∈ [0, 1] from Algorithm 2 with the following inputs: R = 1 (δu bnd ) 2 , ū = 0, Q(t 0 ) = diag(10 5 , 10 7 ), x(t 0 ) = [0, 0] , Ξ = 0, ξ = [0, 0] , ∆ t = 0.1 ms, t f = 1 s. The resulting curve is shown in Fig. 5. It represents the secured zone in which We can observe that the hit time increases with the decrease of the size of the control input set. Hence, the lower the upper bound of u(t), the longer is the period of time in which the normal aging is guaranteed. In addition, from a certain value of δ the dangerous set is not reachable for any time t ≥ 0 for the state trajectory of the brushed DC motor. The corresponding R is the one obtained previously ( R = 1 36.72 2 ) with the results shown in Fig. 2.

CONCLUSION

In this paper, we have shown how it is possible to use convex robust simulation for forecast and prevention of in function of the hit time t ∩ with t f = 1 s, initial conditions E x (diag(10 5 , 10 7 ), [0, 0] ) and E u (R, 0). aging attacks. This paper has introduced this new idea and it has presented some preliminary results on a model of brushed DC motor on a very simple attack benchmark, but we can envisage several extensions as topics of future research. First of all, we have considered as the only means of attack prevention the time-invariant restriction of the control input set, but in principle this can be made timevarying, allowing for less restrictive set. Further research will also investigate more deeply the practical application of the methodology to industrial plants.

Fig. 3 .

 3 Fig. 3. Bounding ellipsoids evolution from initial conditions E x (diag(10 5 , 10 7 ), [0, 0] ) with respect to the time t for t f = 1 s at each 40 ms for E u ( 1 110 2 , 0) (transparent) and for E u ( 1 36.72 2 , 0) (non-transparent)

Fig. 5 .

 5 Fig. 5. Secured zone: control input set with R = 1 (δu bnd ) 2


Christian A Wilhelmsen
email: nils.wilhelmsen@mines-paristech.fr

Florent Di Meglio
email: florent.dimeglio@mines-paristech.fr

An Observer for the Electrically Heated Vertical Rijke Tube with Nonlinear Heat Release

Keywords: Thermoacoustics, Rijke tube, Observers, Nonlinearities, Distributed parameter systems

This paper proposes an observer for the electrically heated vertical Rijke tube, using a model capable of inciting thermoacoustic instabilities. The model consists of linear, distributed acoustics coupled with nonlinear, lumped heat release. Using a boundary pressure measurement taken at the bottom end of the tube, the observer is designed by copying the system model and reconstructing both the bottom and top boundary conditions, the former exactly and the latter with an error that is shown to be converge to zero exponentially. It is proven that the observer produces state estimates that converge to their correct values asymptotically. Furthermore, it is shown the state estimation errors stay bounded when there is a modelling error in the boundary acoustic impedance. The proposed observer is simulated and compared to two alternative observers.

INTRODUCTION

Thermoacoustic instabilities can be a severe issue in combustion chambers, and have been observed in rocket engines, gas turbines, industrial furnaces and afterburners [START_REF] Correa | Power generation and aeropropulsion turbines: From combustion science combustion technology[END_REF]; [START_REF] Crocco | Research on combustion instability in liquid propellant rockets[END_REF]; [START_REF] Macquisten | Combustion oscillations in a twinstream afterburner[END_REF]). They occur due to a positive feedback loop between the heat release from the combustion process and the acoustics in the chamber when these are sufficiently in phase [START_REF] Putnam | Burner oscillations of the gauze-tone type[END_REF]). Oscillations can grow into a limit cycle with amplitude large enough to damage the combustion system. However, it is difficult in practice to forecast when or if thermoacoustic instabilities occur for a given combustor [START_REF] Lieuwen | Combustion instabilities in gas turbine engines: operational experience, fundamental mechanisms, and modeling[END_REF]), and hence much research effort has been devoted to understanding the phenomenon better, and also developing techniques for monitoring and mitigating the instabilities.

One common experimental setup for studying thermoacoustic instabilities, being one of the simplest arrangements capable of exhibiting the phenomenon, is the electrically heated Rijke tube. It has the advantage of being simple to model due to the absence of complicated combustion dynamics. The Rijke tube (with a flame) was first introduced in Rijke (1859), and has since been the subject of numerous studies to gain understanding of and develop methods to mitigate thermoacoustic instabilies, see [START_REF] Raun | A review of rijke tubes, rijke burners and related devices[END_REF] for a review.

Much of the previous work on the Rijke tube has consisted in characterizing its stability limits, with [START_REF] Carrier | The mechanics of the rijke tube[END_REF] pioneering linear stability analysis of the system and later [START_REF] Bayly | Onset and equilibration of oscillations in general rijke devices[END_REF] taking into account nonlinear features. In addition to experimental studies of the stability limits of the Rijke tube, studies on active control strategies applied to attenuate the thermoacoustic oscillations in the Rijke tube have been performed. A control law consisting of measuring the pressure signal upstream of the heater and subsequently sending this signal phase-shifted and amplified to a loudspeaker downstreams has been applied to a Rijke tube in [START_REF] Heckl | Active control of the noise from a rijke tube[END_REF], being one of the first studies investigating active control of the Rijke tube. More recently, a more sophisticated full-state feedback boundary control law designed via infinite-dimensional backstepping on a linearized PDE-ODE model of the electrically heated Rijke tube has been derived in de [START_REF] De Andrade | Backstepping stabilization of a linearized ode-pde rijke tube model[END_REF]. To pair with this full-state feedback control law, a corresponding boundary observer for the linearized PDE-ODE model is derived in de Andrade et al. (2018a).

A heat release model which, together with an acoustic model, captures the most important behaviour of the electrically heated Rijke tube is King's law [START_REF] King | Xii. on the convection of heat from small cylinders in a stream of fluid: Determination of the convection constants of small platinum wires with applications to hot-wire anemometry[END_REF]). For the observer design in de Andrade et al. (2018a), King's law is linearized, which makes the mathematical analysis tractable -however this linear ODE model does not reflect the full nonlinear dynamics one typically obtains in practice. A nonlinear heat release model is needed to model the saturated response one sees for large amplitudes and resultant limit cycle behaviour [START_REF] Agostino | Nonlinear identification of thermoacoustic instabilities with limit cycles in a rijke tube[END_REF]). To maintain this behaviour in the design, we propose a differ- ent observer taking into account the nonlinear features of King's law. This is the main contribution of the paper.

Our approach is as follows. We design an observer by injecting measured signals at two boundaries of an infinite dimensional observer. The zero equilibrium of the resulting error system is shown to be Globally Asymptotically Stable (GAS). Further, we show that the presence of uncertainty in the boundary acoustic impedance results in bounded errors.

The paper is organized in the following way. In Section 2 the governing equations are introduced and the problem we are solving is defined. Next, in Section 3 the design and convergence analysis of the observer is covered. Results from simulations are presented in Section 4 before some concluding remarks are offered in Section 5.

PROBLEM STATEMENT

Consider a vertical Rijke tube of length L and constant circular cross-sectional area A with an electrical heater placed at a fixed vertical position z 0 ∈ (0, L 2 ) 1 . The phenomenon we are studying evolves mainly in the vertical direction inside the Rijke tube, so we use one-dimensional (1D) first principles to model the dynamics of the air column inside the tube (de Andrade et al. ( 2018b)). To establish an acoustic model, define pressure perturbation P := P -P , velocity perturbation V := V -V , and heat fluctuation Q := Q -Q, with P total pressure, P constant atmospheric pressure, V total velocity, V constant base flow velocity, Q heat release and Q its constant mean value. Also, the mean air density is denoted ρ. We write, as in [START_REF] Epperlein | Thermoacoustics and the rijke tube: Experiments, identification, and modeling[END_REF], the linearized mass, momentum and energy conservation laws of 1D gas dynamics to obtain the following linear acoustic model, 1 A necessary condition for thermoacoustic instabilities to occur is for the heater to be in the lower half of the tube.

∂ t V (z, t) + 1 ρ ∂ z P (z, t) = 0 (1a) ∂ t P (z, t) + γ P ∂ z V (z, t) = γ A δ(z -z 0 ) Q(t) (1b) P (L, t) = Z V (L, t) + U (t) (1c) P (0, t) = -Z V (0, t) (1d) P (x, 0) = P0 (x) (1e) V (x, 0) = V0 (x),
(1f) defined over the spatial domain z ∈ [0, L] and for time t ≥ 0, initialized from P0 , V0 ∈ L 2 (0, L). Also, γ is the adiabatic constant, γ := γ -1, Z is an acoustic impedance at the tube boundaries and U is a boundary gauge pressure actuation signal.

To describe the heat release rate Q from the electrical heater, King's law [START_REF] King | Xii. on the convection of heat from small cylinders in a stream of fluid: Determination of the convection constants of small platinum wires with applications to hot-wire anemometry[END_REF]) is used, which together with dynamics described in [START_REF] Lighthill | The response of laminar skin friction and heat transfer to fluctuations in the stream velocity[END_REF] reads

Q(t) = - 1 τ Q(t) + 1 τ l w (κ + κ v |V (z 0 , t)|)(T w -T g ) (2a) Q(0) = Q 0 , (2b
) and is initialized from Q 0 ∈ R, with τ a time-constant due to the boundary layer formed around the electrical heater, l w the length of the heating wire, T w is the wire temperature, T g is the time-averaged gas temperature around z 0 , κ is the thermal conductivity of air and κ v is an empirically determined constant. The model (2) captures the two most important dynamic features required of the heat source for thermoacoustic instabilities to occur, namely a time lag between changes in the local velocity and the heat release rate, and also a positive dependency between the time derivative of the heat release rate and the local air speed. Unlike previous attempts at designing observers for the Rijke tube, as in de Andrade et al. (2018a), we do not linearize King's law in (2) to capture the full nonlinearity of its dynamics.

Hence, given the Rijke tube thermoacoustic model (1)-(2), the objective of this paper is to design an observer producing convergent state estimates P , V and Q of respectively the pressure P , velocity V and heat release Q, assuming the boundary pressure measurement Y (t) := P (0, t)

(3) is available, only. As noted in [START_REF] Epperlein | Thermoacoustics and the rijke tube: Experiments, identification, and modeling[END_REF], in practice the nodes in the Rijke tube are situated slightly outside the tube ends. Hence by placing the pressure sensor right by the tube opening or slightly inside, the node is avoided and one obtains information that can be used for state estimation. Next, we proceed with the observer design and assess its convergence properties.

OBSERVER DESIGN

Model in Riemann coordinates

For the observer design and analysis, we rewrite the linearized acoustics (1) in Riemann invariant coordinates and fold the spatial domain around z 0 to move the heat release to the system boundary. To facilitate this, we introduce the invertible linear spatial coordinate transforms

z i : x → z, z 1 (x) := z 0 (1 -x) (4a) z 2 (x) := z 0 + x(L -z 0 ), ( 4b
) with x ∈ [0, 1] and i ∈ {1, 2}. The subscript i denotes which part of the Rijke tube x is mapped to, with z 1 mapping x to points below the electrical heater and z 2 mapping x to points above the electrical heater. Next, define the Riemann coordinates

u i (x, t) := P (z i (x), t) + k V (z i (x), t) (5a) v i (x, t) := P (z i (x), t) -k V (z i (x), t). ( 5b
)
where

k := γ P ρ. ( 6
)
This allows us to rewrite the linearized acoustics (1a)-( 1b) over (x, t)

∈ (0, 1) × [0, ∞) as ∂ t u 1 (x, t) = λ 1 ∂ x u 1 (x, t), ∂ t u 2 (x, t) = -λ 2 ∂ x u 2 (x, t), (7a) ∂ t v 1 (x, t) = -λ 1 ∂ x v 1 (x, t), ∂ t v 2 (x, t) = λ 2 ∂ x v 2 (x, t), (7b)
where

λ 1 := λ z 0 , λ 2 := λ L -z 0 , λ := γ P ρ , (8)
with λ being the speed of sound inside the tube. Remark 1. Since z 0 ≤ L 2 , we have that λ 1 ≥ λ 2 . This fact is useful later on in the observer design.

Next, the acoustics boundary conditions (1c)-( 1d) are rewritten as

u 2 (0, t) = u 1 (0, t) + µ X(t), v 1 (0, t) = v 2 (0, t) + µ X(t), (9a) u 1 (1, t) = dv 1 (1, t), v 2 (1, t) = du 2 (1, t) + W (t), (9b)
where we have denoted X := Q, the boundary reflection coefficient d, the ODE boundary coefficient µ and input signal W , are respectively defined

d := Z -k Z + k , µ := kγ γA P , W (t) := (1 -d)U (t). ( 10
)
Denoting X := Q, the heat release model ( 2) is also rewritten as

Ẋ(t) = -aX(t) + b 1 |b 2 + b 3 [u 1 (0, t) -v 2 (0, t)]| + b 4 (11) with a := 1 τ , b 1 := l w (T w -T g )κ v τ , b 2 := V , (12a) b 3 := 1 2k , b 4 := l w (T w -T g )κ τ . ( 12b
)
The gauge pressure measurement (3) can in the Riemann invariant coordinates (5) be written as

Y (•) = 1 2 (u 1 (1, •) + v 1 (1, •)).
Applying the boundary condition for u 1 we see by defining the boundary measurement signal y(t) := v 1 (1, t) (13) the gauge pressure measurement can be reconstructed as Y = 1+d 2 y. These dynamics are schematically depicted in Figure 2. Notice that several feedback loops make the dynamics potentially unstable.

Proposed observer

With the measurement signal y defined in ( 13), we propose the following observer

Ẋ(t) = -a X(t) + b 1 |b 2 + b 3 [û 1 (0, t) -v2 (0, t)]| + b 4 (14a) ∂ t û1 (x, t) = λ 1 ∂ x û1 (x, t) (14b) ∂ t û2 (x, t) = -λ 2 ∂ x û2 (x, t) (14c) ∂ t v1 (x, t) = -λ 1 ∂ x v1 (x, t) (14d) ∂ t v2 (x, t) = λ 2 ∂ x v2 (x, t) (14e) û2 (0, t) = û1 (0, t) + µ( X(t) -X) (14f) v1 (0, t) = v2 (0, t) + µ( X(t) -X) (14g) û1 (1, t) = dy(t) (14h) v2 (1, t) = dy(t + λ -1 1 -λ -1 2 ) + W (t) + d û1 (0, t -λ -1
2 ) -v2 (0, t -λ -1 2 ) . (14i) We state now the main result of this paper, which pertains to the convergence properties of ( 14), before explaining the rationale behind the observer followed by a formal proof of the result. Theorem 3.1. Consider system (7)-( 12), and the state observer (14) using the measurement (13). Assume they have initial conditions (u i,0 , v i,0 , X 0 ) ∈ L 2 (0, 1) × L 2 (0, 1) × R and (û i,0 , vi,0 , X0 ) ∈ L 2 (0, 1) × L 2 (0, 1) × R, respectively. Then, the zero equilibrium of the dynamics of the estimation errors ũi := u i -ûi , ṽi := v i -vi , X := X -X is GAS.

This observer consists of a copy of the Rijke tube dynamics in Riemann coordinates ( 7)-( 12), with the exception of (14h)-(14i). While (14h) consists of injecting the measured output directly, Equation (14i) deserves more explanation. It is based on the following considerations. First, notice that substituting Equation (9a) into the general solution of u 2 in terms of the boundary at

x = 0, u 2 (1, t) rewrites u 2 (1, t) = u 2 (0, t -λ -1 2 ) = u 1 (0, t -λ -1 2 ) + µ X(t -λ -1
2 ) (15) Besides, using the expression for v 1 (0, t) from (9a) together with the measurement y 15) and ( 16) yields

(t) = v 1 (1, t) yields µ X(t) = y(t + λ -1 1 ) -v 2 (0, t) (16) Combining (9b), (
v 2 (1, t) =d u 1 (0, t -λ -1 2 ) -v 2 (0, t -λ -1 2 )
Fig. 3. Schematic of the error system signals. The only internal feedback is ṽ2 (1, •) feeding into itself.

+y(t + λ -1 1 -λ -1 2 ) + W (t) (17)
The boundary condition (14i) follows by considering estimates of u 1 (0, t-λ -1

2 ) and v 2 (0, t-λ -1 2 ) in lieu of the true values. As we show in the next Section, the resulting error converges asymptotically to zero. Remark 2. Notice that the observer is causal, in particular the signal y(t + λ -1 1 -λ -1 2 ) is available at time t due to Remark 1.

Convergence analysis

The most critical boundary conditions in the error system for stability are ũ1 (1, •), ṽ2 (1, •), so we derive their expressions first. It is trivial to see that ũ1 (1, •) = 0, while subtracting (14i) from (17

) yields ṽ2 (1, t) =d ũ1 (0, t -λ -1 2 ) -ṽ2 (0, t -λ -1 2 ) . (18)
The other terms in the observer ( 14) are copies of the corresponding terms in the original system ( 7)-( 12) and hence their corresponding error dynamics are easily computed. Therefore, the state estimation error in ũi , ṽi , X satisfies the dynamics

∂ t ũ1 (x, t) = λ 1 ∂ x ũ1 (x, t) (19a) ∂ t ũ2 (x, t) = -λ 2 ∂ x ũ2 (x, t) (19b) ∂ t ṽ1 (x, t) = -λ 1 ∂ x ṽ1 (x, t) (19c) ∂ t ṽ2 (x, t) = λ 2 ∂ x ṽ2 (x, t) (19d) ũ1 (1, t) = 0 (19e) ũ2 (0, t) = ũ1 (0, t) + µ X(t) (19f) ṽ1 (0, t) = ṽ2 (0, t) + µ X(t) (19g) ṽ2 (1, t) = d ũ1 (0, t -λ -1 2 ) -ṽ2 (0, t -λ -1 2 ) (19h) Ẋ(t) = -a X(t) + b 1 |b 2 + b 3 [u 1 (0, t) -v 2 (0, t)]| -b 1 |b 2 + b 3 [û 1 (0, t) -v2 (0, t)]|.
(19i) A schematic view of the error system is shown in Figure 3, illustrating the cascade structure of its dynamics which ensure the convergence of its states to zero. We are now ready to prove Theorem 3.1.

Proof. [Proof of Theorem 3.1] From (19a), (19e) we see that ũ1 (0, t -λ -1 2 ) = 0 for time t ≥ λ -1 1 + λ -1 2 . Hence the boundary condition (19h) simplifies after this time to ṽ2 (1, t) = -dṽ 2 (0, t -λ -1 2 ) (20) allowing us to conclude that ṽ2 (0, t) = -dṽ 2 (0, t -2λ -1 2 ).

(21) Since |d| < 1, we can conclude that ṽ2 (0) → 0 exponentially as t → ∞. Next, we can bound (19i) by the following inequality:

Ẋ(t) ≤ -a X(t) + b 1 |b 3 [ũ 1 (0, t) -ṽ2 (0, t)]|. ( 22
) Define g(t) := b 1 |b 3 ι(t)|, ι(t) := ũ1 (0, t) -ṽ2 (0, t) (23
) where we know g tends to zero as t → ∞. Also, defining the linear system in Ξ, Ξ(t) = -a Ξ(t) + g(t), (24) with initial condition Ξ(0) = X(0), then Ξ is Input-to-State Stable (ISS) with respect to g. More precisely (see [START_REF] Khalil | Nonlinear systems[END_REF]), we can establish the bound (where 0 ≤ t 0 ≤ t),

| Ξ(t -t 0 )| ≤ e -a(t-t0) | Ξ(t 0 )| + 1 a sup t0≤τ ≤t |g(τ )|. ( 25
)
Since Ẋ ≤ Ξ and X(0) = Ξ(0), we can establish X(t) ≤ Ξ(t).

(26) As g is exponentially vanishing as t → ∞, we see the right hand side of (25) goes to zero and hence X → 0 asymptotically as t → ∞, which proves the Theorem.

Robustness to modelling error in boundary condition

In practice the acoustic impedance Z appearing in Equations (1c)-( 1d) is difficult to estimate correctly, implying the boundary coefficient d appearing in (9b) for the Riemann coordinate formulation is prone to being incorrectly modelled. Denote the estimate of d as d, and define d := d -d as the modelling error.

When d is used in place of d in the observer ( 14), the expressions for ũ1 (1), ṽ2 (1), given by (19e), (19h), change to ũ1 (1, t) = dy(t) ( 27)

ṽ2 (1, t) = d(ũ 1 (0, t -λ -1 2 ) -ṽ2 (0, t -λ -1 2 ) (28) + d µX(t -λ -1 2 ) + u 1 (0, t -λ -1
2 ) with the rest of the error system (19) being unaffected. We state now a Proposition on the sensitivity of the state estimate error to errors in modelling this boundary condition. Proposition 3.2. Assume that the states of the original system are bounded, and the estimate d of d satisfies | d| < 1.

(29) Then all the error signals ũi , ṽi , X are bounded.

Proof. In the following we denote the Laplace transform of the time-domain signal f as f , i.e. f (s) = L(f (•)). With the signal h defined as h(t) := u 1 (0, t) + µX(t) (30) we find that ǔ1 (0, s)

v2 (0, s) = H(s) y(s) ȟ(s) (31)
where

H(s) :=   de -λ -1 1 s 0 de -(2λ -1 2 +λ -1 1 )s 1+ de -2λ -1 2 s e -2λ -1 2 s 1+ de -2λ -1 2 s   ( 32
)
is a transfer matrix. This allows us to reconstruct the signal ι defined in (23) in terms of signals y, h as

ι(t) = 1 -1 T L -1 H(s) y(s) ȟ(s) (33)
Since | d| < 1, Equation ( 32) implies that H is stable (see [START_REF] Niculescu | Delay effects on stability: a robust control approach[END_REF] for a more extended treatment of transfer functions for systems with time delays), and since y, h ∈ L ∞ , one has ι, g ∈ L ∞ , with g defined by ( 23).

Further, since ũ1 (1, •) = dy(•), ũ1 (1, •) ∈ L ∞ and hence all of ũ1 is bounded. Besides, one has

v2 (1, s) = de -λ -1 2 s 1 + de -λ -1 2 s dy(s)e -λ -1 1 s + ȟ(s) (34)
which again using that | d| < 1 implies that ṽ2 (1, •) ∈ L ∞ and therefore all of ṽ2 is bounded. Equation ( 22) further implies the following bound on X,

X(t) ≤ e -at | X(0)| + 1 a ||g|| ∞ . (35) Finally, (19f)-(19g) express ũ2 (0, •), ṽ1 (0, •) as the sum of bounded signals, therefore ũ2 (x, •), ṽ1 (x, •) ∈ L ∞ , ∀x ∈ [0, 1].

SIMULATIONS

In this section the observer ( 14) is tested in simulations implemented in MATLAB. The proposed observer has certain advantages, such as being guaranteed stable for all system parameters, global asymptotic convergence and being robust to modelling errors in the boundary condition d.

However, it also has certain drawbacks, such as having no tuning parameters to calibrate and possible poor measurement noise filtering properties as the measurement is injected directly.

For comparison we test, where possible, the observer from this paper with two alternative observers:

(1) A "trivial" observer identical to the observer ( 14 A total of three different simulations scenarios are presented. In the first two simulations, the Rijke tube model ( 1)-( 2) is implemented with parameters that incite thermaocoustic instabilities. No actuation signal is imposed on the top boundary so we assume U = 0. The physical parameters used are summarized in Table 1. For the observers in the first simulation, it is assumed all system parameters are known exactly, but in the second simulation an incorrect value of Z is used in the observers. Lastly, results from a toy simulation, that implements the system (7), ( 9), ( 11), together with the observer ( 14) and observers mentioned above, is presented. This demonstrates that observer (14) in this paper works for certain system parameters the alternative observers do not work for, due to them having no proven stability guarantee.

All results presented are from simulations that integrate the equations for a total of 1 second, with the observers being turned on at time t = 0.25 seconds. The PDEs are integrated using a first-order finite difference upwind scheme, while the ODE is integrated with a fourth-order Runge-Kutta scheme. The spatial discretization used is dx = 1.00 × 10 -2 m, while the temporal discretization is dt = 8.25 × 10 -8 s. All observers are initialized to û1,0 (x) = sin(2x) û2,0 (x) = cos(3x) v1,0 (x) = e x v2,0 (x) = sin(5x) -2 cos(x) X0 = 500.

Thermoacoustic simulations with no modelling error

In this simulation it is assumed the observers have perfect knowledge of the system parameters in Table 1. First, to showcase the performance of the proposed observer ( 14), in Figures 456we see respectively the estimation errors for the pressure P := P -P and velocity Ṽ := V -V , evaluated at positions z = L, z = z 0 , and the heat release Q := Q -Q, all associated with the estimates P , V , and Q produced by this observer. As predicted by Theorem 3.1, these approach the origin asymptotically as t → ∞.

Next, the transient performance of the various observers in estimating the heat release rate Q is compared. Even though the equations were integrated up to time t = 1 second, for clarity a time-interval of length 0.2 seconds around the time the observers are turned on is shown. In Figure 7 the heat release rate estimation error associated with the different observers is displayed for times between t = 0.24 and t = 0.44 seconds. This interval contains the time instant when observers are turned on, so the initial behaviour can be seen from this plot. The "trivial" observer seems to converge fastest, with the observer using backstepping gains converging second fastest and the observer ( 14) from this paper converging more slowly, but asymptotically, as is shown in Figure 6. Interestingly, in this example the observer ( 14) had overall the worst performance. This can be attributed to the fact that the error -dṽ 2 (0, t -λ -1 2 ) is introduced in the boundary for v2 (1) to ensure stability, with the trade-off being that the performance is degraded initially, before improving as time increases. Next, the performance of the observers when there is a modelling error in the boundary coefficient is compared. Remark 3. In this example, the convergence properties of the "trivial" observer and the nonlinear observer using the backstepping gains is similar, showing that the observer gains have little extra effect on the convergence properties.

W m -1 K -1 κv Empirical constant in King's law 1.5 W s 1 2 m -3 2 K -1 ρ Average air density 1.2 kg m -3 P Atmospheric pressure 10 5 N m -2 V Base flow velocity 0.35 m s -1 γ Adiabatic constant 1.4 - Table 1. Physical

Thermoacoustic simulations with modelling error in boundary

Here the observers are tested against modelling error in the reflection coefficient d. Therefore, for this simulation it is assumed the observers have correct knowledge of all the physical parameters in Table 1 except for Z. Instead the observers use an estimate Ẑ = 50, which gives an estimate d ≈ -0.78.

Overall, the effect of the incorrect knowledge of Z is that the observers have an offset in the state estimate. The magnitude of this offset for the various observers in estimating the heat release rate can be seen in Figure 8. From this plot it can be seen that the observer ( 14) has the smallest offset from the correct value, followed by the observer using backstepping gains, and lastly the "trivial" observer has the largest offset. It is interesting to note that despite having incorrect knowledge of the boundary condition, the observer ( 14) performs better during the transient for the simulation shown here compared to the simulation in Section 4.1. This is due to the smaller of d the error -2 (0, t -λ 2 )) at the boundary for v2 •) shrink faster in the transient.

The simulations suggest the design of the observer ( 14) makes it more robust to modelling errors in the boundary reflection coefficient d compared to the other two observers. For the observer using backstepping gains, the gains helped make the offset smaller than the offset from the "trivial" observer, but not to the extent that the observer ( 14) managed. In Proposition 3.2 it was proven that the estimation error produced by the observer (14) when using an incorrect boundary condition would be boundedfor the other two observers this has of yet not been proven.

Another performance aspect that is not necessarily guaranteed for the other two observers is that they converge for all possible system parameters. Next a toy simulation is offered where the observer ( 14) converges to the correct estimate but the other observers become unstable.

Toy simulations showing stability for alternative system parameters

Due to the way the boundary condition (14i) is constructed for the observer (14), it is possible to prove that the state estimates produced by this observer converge to their true values. For the trivial observer, due to it having a potentially unstable feedback loop of previous ODE states influencing the dynamics of the current ODE states via the nonlinearity present in King's law, there exists certain parameters for which the observer can become unstable, As can be seen, the estimate Xtrivial shoots up to values much larger than the range X is oscillating within. The observer (14), however, remains unscathed and produces an estimate X that tracks the state X fairly well, asymptotically.

CONCLUSIONS

An observer taking into account both the nonlinearity of the heat release model and distributed states of the tube acoustics has been derived for the electrically heated vertical Rijke tube. Instead of applying a standard observer design procedure the observer was designed simply Fig. 9. ODE state estimation error X produced by the observer ( 14) and Xtrivial produced by the trivial observer, for a toy example. by copying the system dynamics and reconstructing some of its boundary conditions from known signals.

It was proven that this observer converges to the origin asymptotically and that it is robust to modelling errors in the boundary condition. Simulations were offered to demonstrate the observer in three examples. Compared to an observer that only injects the measurement at a single boundary, the observer from this paper does not suffer the same problem of blowing up for certain system parameters.

As further work, simulating the observer against a more realistic flow model, such as the nonlinear gas dynamics model based on mass, momentum and energy balance given in Epperlein et al. would give insight about observer performance from a practical point. observer design for the Rijke tube using a flame rather than electrical heater model could be a viable research direction.

Fig. 1 .

 1 Fig. 1. Electrically heated vertical Rijke tube. The magenta arrow represents base flow induced by the electrical heater.

Fig. 2 .

 2 Fig. 2. Schematic of Rijke tube system in Riemann coordinates. The green arrows are couplings between the ODE state and PDE states, the red arrows are boundary couplings between the PDE states and the magenta arrows represent input and output signals.

  ), except the boundary condition v2 (1,•) = dû 2 (1, •) + W (•) is used in place of (14i).(2) A nonlinear observer consisting of the "trivial" observer together with observer gains from deAndrade et al. (2018a), multiplied by the output estimation error (y(•) -v1 (1, •)) and added to the observer dynamics (14a)-(14e).

Fig. 4a .Fig. 4b .

 4a4b Fig. 4a. Pressure estimation error P (L, •) at tube end opposite from where sensing takes place, for observer (14).

Fig. 5a .

 5a Fig. 5a. Velocity estimation error Ṽ (L, •) at tube end opposite from where sensing takes place, for observer (14).

Fig. 5b .

 5b Fig. 5b. Velocity estimation error Ṽ (z 0 , •) at heater location for observer (14).

Fig. 6 .

 6 Fig. 6. Heat release rate estimation error Q for observer (14).

Table 2 .

 2 Parameters for toy simulation and the observer in deAndrade et al. (2018a) was only proven to be stable by assuming that the heat release model is linearized.As a demonstration of what could happen in practice, the system (7), (9), (11) with parameters specified in Table2is implemented. Figure9shows a plot of the resultant state X, shown in blue, versus the estimate Xtrivial produced by the trivial observer, shown in yellow, and the estimate X produced by the observer (14), shown in orange. The estimate produced by the nonlinear observer using backstepping gains from de Andrade et al. (2018a) can not be shown because it becomes unstable very quickly.

	1000
	800
	600
	400
	200
	0
	-200
	0.24	0.26	0.28	0.3	0.32	0.34	0.36	0.38	0.4	0.42	0.44
	Fig. 8. Heat release rate estimation error Q produced
	by observer (14), Qtrivial produced by the "trivial"
	observer and Qbackstepping produced by the observer
	from de Andrade et al. (2018a) but using a nonlinear
	ODE, with all observers having incorrect knowledge
	of the acoustic boundary Z. The plot shows the curves
	for times t ∈ [0.24, 0.44].
				Parameter	Value
				a		500
				b 1		474750000
				b 2		3
				b 3		0.001
				b 4		8349270
				d		-0.9711
				µ		0.2925
				λ 1		975.9
				λ 2		325.3

ACKNOWLEDGEMENTS

The authors would like to thank Laurent Praly for his insightful comments on the observer design.

The work of N.C.A. Wilhelmsen has received funding through the European Union's Horizon 2020 research and innovation programme under the Marie Sk lodowska-Curie grant agreement No. 766264.