Laleh Ravanbod

Nathalie Verdière

Carine Jauberthie

Determination of set-membership identifiability sets

Keywords:

HAL is

Introduction

The concept of identifiability is important for insuring that a mathematical model is well-defined in the sense that there exists an unambiguous mapping between its parameters and the output trajectories. This concept insures that parameters estimation procedure will succeed and will give trustable results [START_REF] Verdière | Identifiability and estimation of pharmacokinetic parameters of ligands of macrophage mannose receptor[END_REF]. However, uncertainties on the knowledge of measures or parameters for example, have to be considered in the numerical procedure. The stochastic framework can take into account these uncertainties but with appropriate assumptions about noise and model error distributions. However, some uncertainties are better characterized by simple bounds. For example, tolerances on parameter values are provided by the manufacturer as lower and upper bounds corresponding to the inherent variability of production line. That is why set-membership (SM) models, in which errors and disturbances affecting a dynamical system are represented by prior bounds, represent an advantageous alternative to stochastic framework. Moreover, success of operational estimation methods developed these last years on SM models confirm their interest, for example [START_REF] Verdière | Set-membership identifiability of nonlinear models[END_REF][START_REF] Raıssi | Set-membership state and parameter estimation for systems described by nonlinear differential equations[END_REF]. Even theoretical concepts can be treated by interval analysis as the injectivity of a differentiable function. For example, in [START_REF] Lagrange | Injectivity analysis using interval analysis : Application to structural identifiability[END_REF], a numerical algorithm based on interval analysis is presented to determine domains on which a differentiable function is injective and thus, can permit to study the structural identifiability of parametric models. Their algorithms partition the study domain into two domains: an undetermined one on which the function is not proved to be injective and a domain on which the function is injective.

The concept of identifiability of bounded-error uncertain models has been introduced in ([2], [START_REF] Jauberthie | Set-membership identifiability: definitions and analysis[END_REF]). In [START_REF] Jauberthie | Set-membership identifiability: definitions and analysis[END_REF], two definitions named SM and µ-SM-identifiability have been proposed. The first one is purely conceptual whereas the second one subsumes classical identifiability and can be put in correspondence with interval based parameter estimation methods. A SM-identifiable set is a connected set of the parameter domain such that the trajectories that it generates are distinct to trajectories arising from its complementary. Thus, the search of global (µ-)SM-identifiable sets permits to determine regions in the parameter space leading to an unique behavior of the corresponding trajectories compared to the rest of parameter space. In practice, it can be important to know if the feasible parameter set is reduced to one connected set or not and in which regions of the parameter space it is located, while taking into account the uncertainties. For example, this knowledge can be useful for doing SM identification and parameter estimation (see for instance [START_REF] Walter | Estimation of parameter bounds from bounded-error data: a survey[END_REF][START_REF] Walter | Bounding Approaches to System Identification[END_REF][START_REF] Milanese | Estimation theory for nonlinear models and set membership uncertainty[END_REF][START_REF] Raıssi | Set-membership state and parameter estimation for systems described by nonlinear differential equations[END_REF]) or detection [START_REF] Jauberthie | Fault detection and identification relying on setmembership identifiability[END_REF]. The aim of this paper is to determine numerically these SM-identifiable connected sets that correspond to different output behaviors for the considered system. For doing this, a method based on differential algebra is proposed. This method consists in linking the µ-SM-identifiability/SMidentifiability definitions to the partial injectivity/restricted-partial injectivity of a real rational function. The definition of partial injectivity has been introduced in [START_REF] Lagrange | Injectivity analysis using interval analysis : Application to structural identifiability[END_REF] and characterizes perfectly the concept of µ-SM-identifiability. We have completed this definition in introducing the notion of restricted-partial injectivity for characterizing the concept of SM-identifiability. In order to determine numerically these SM-identifiable sets, we have implemented an algorithm composed of three steps ( [START_REF] Jaulin | Applied Interval Analysis, with examples in parameter and state estimation, Robust control and robotics[END_REF]): the mince step, the evaluate step and the regularize step. Indeed, the method consists in calculating the image of the real rational function obtained by the theoretical part. However, it is well-known that it is not always an easy task since interval analysis does not provide any inclusion test for the point test y ∈ f (X). In [START_REF] Jaulin | Applied Interval Analysis, with examples in parameter and state estimation, Robust control and robotics[END_REF], they propose three steps in order to obtain a guaranteed inclusion box containing the image of a continuous function. The first step, named mincing consists in building a non-minimal regular subpaving. The second one, the evaluate step, computes the image of a box by the inclusion function. The third one, the regularize step, computes a regularize subpaving containing the image of the inclusion function. However, If we keep their algorithm, the obtained boxes, supposed to correspond to SM-identifiable sets will be composed of non desired parts, that is parts of non SM-identifiable sets. That is why, an algorithm to find a regular internal subpaving is proposed in this step.

The paper is organized as follows. After having recalled the SM-identifiability definitions in the second section, a method based on differential algebra is proposed in order to determine the SMidentifiable sets. This method is based on a proposition linking the (µ-)SM-identifiability and the (restricted-)partial injectivity of a real rational function. In Section 3, we present our algorithm based on guaranteed numerical computation to find the SM-identifiable sets in the feasible parameters set of a model. In the fourth part, illustrative examples show the efficiency of our algorithm.

The concept of set-membership identifiability

In this section, we recall the notion of SM-identifiability published in [START_REF] Jauberthie | Set-membership identifiability: definitions and analysis[END_REF][START_REF] Jauberthie | Fault detection and identification relying on setmembership identifiability[END_REF] for the class of systems formalized by (2.1) :

Γ =            ẋ(t, p) = g(x(t, p), u(t), p), y(t, p) = h(x(t, p), p), x(t 0 , p) = x 0 ∈ X 0 , p ∈ P ⊂ U P , t 0 ≤ t ≤ T, (2.1)
where :

• x(t, p) ∈ R n and y(t, p) ∈ R m denote the state variables and the outputs at time t respectively.

• u(t) ∈ R r is the input vector at time t.

• the initial conditions x 0 , if any, are supposed to belong to a bounded set X 0 .

• the vector of parameters p belongs to a connected set P supposed to be included in U P where U P is an a priori known set of admissible parameters. U P is either included in R p or equal to R p . • the functions g and h are real and analytic1 on M , where M is an open set of R n such that x(t, p) ∈ M for every t ∈ [t 0 , T ] and p ∈ P . T is a finite or infinite time bound. Afterwards, P c will design the complementary of P in the considered space.

Definitions

The definitions that are proposed are given in the case of controlled systems but they can be stated similarly in the case of uncontrolled systems. In these definitions, Y Γ (P, u) denotes the set of outputs, solution of Γ with the input u, where P represents a connected set of R p . Y Γ (P, u) is also called the output behavior of Γ arising from P . P * is a connected set of R p . Definition 2.1. The model Γ given by (2.1), P * = ∅, P * ⊆ U P , is globally SM-identifiable if there exists an input u such that Y Γ (P * , u) = ∅ and Y Γ (P * , u)∩Y Γ ( P , u) = ∅, P ⊆ U P =⇒ P * ∩ P = ∅.

The previous definition expresses that a connected set P * is globally SM-identifiable if the output behavior of Γ arising from P * , i.e. the output behavior of Γ for any p ∈ P * , is distinguishable from the output behavior of Γ arising from its complementary set P * c , i.e. Y Γ (P * , u) and Y Γ (P * c , u) do not share any identical output trajectory.

Let us now consider a nonempty bounded connected set Π of R p and d a classical metric on R p [START_REF] Bourbaki | Elements of Mathematics[END_REF], [START_REF] Munkres | Topology a first courses[END_REF]. On the metric space (Π, d), let µ be a continuous map from Π to Π. µ is a contraction if there is a nonnegative number k < 1 such that for all π 1 , π 2 in Π, d(µ(π 1 ), µ(π 2 )) < kd(π 1 , π 2 ) [START_REF] Munkres | Topology a first courses[END_REF]. Let us also define the diameter of Π by the least upper bound of {d(π 1 , π 2 ), π 1 , π 2 ∈ Π}.

In the following definition, the set P * is supposed to be a bounded connected set and µ is a contraction from P * to P * . Definition 2.2. The nonempty bounded connected set P * ⊆ U P is globally µ-SM-identifiable if, for all contractions µ from P * to P * , µ(P * ) is globally SM-identifiable.

Under the conditions of definition 2.2, we may equivalently say that the model Γ given by (2.1) is globally µ-SM-identifiable with respect to P * . Definition 2.2 differs from definition 2.1 in the sense that the set P * may be reduced as small as desired by the contraction µ while still retaining the property of SM-identifiability. This is true by Banach fixed-point theorem, which implies that the diameter of µ(P * ) tends to zero [START_REF] Munkres | Topology a first courses[END_REF]. In this case, µ-SM-identifiability meets classical identifiability and, interestingly, it means that classical identifiability holds for any p ∈ P * . In [START_REF] Jauberthie | Set-membership identifiability: definitions and analysis[END_REF], it has been proven that the µ-SM-identifiability is equivalent to the classical identifiability (see [START_REF] Ljung | On global identifiability for arbitrary model parametrizations[END_REF] for a good survey of this definition), that the µ-SM-identifiability implies the global SM-identifiability but the reciprocal is not true.

In Figure 1, the sets A i , i = 1, . . . , 5 are SM-identifiable and our aim is to determine

F={F 0 = A 0 , F 1 = A 1 , F 2 = A 3 }.
This list constitutes a list of SM-identifiable sets such that the trajectories generated by these connected sets are distinct to their complementary.

Partial injectivity

The definition of partial injectivity of a function has been introduced in [START_REF] Lagrange | On sufficient conditions of injectivity, development of a numerical test via interval analysis[END_REF]. This notion perfectly characterizes µ-SM-identifiability. A second definition named restricted-partial injectivity completes the first one and characterizes the global SM-identifiability.

Definition 2.3. Consider a function f : A → B and any set

A 1 ⊆ A. The function f is said to be a partial injection of A 1 over A, noted (A 1 , A)-injective, if ∀a 1 ∈ A 1 , ∀a ∈ A, a 1 = a ⇒ f (a 1 ) = f (a). A0 A1 A4 A5 A2 a1 A3 a0 FIGURE 1. Example of SM-identifiable sets f is said to be A-injective if it is (A, A)-injective.
In [START_REF] Lagrange | Injectivity analysis using interval analysis : Application to structural identifiability[END_REF], an algorithm based on interval analysis for testing the injectivity of a given differentiable function is presented and a solver called ITVIA (Injectivity Test Via Interval Analysis) implemented in C++ is mentioned. From a given function, the solver partitions a given box into two domains: a domain on which the function is partially injective and an indeterminate domain on which the function may or not be injective. When the latter is empty, the function is injective over the initial box. In order to characterize the global SM-identifiability, the notion of restricted-partial injectivity is introduced. The algorithm proposed in [START_REF] Lagrange | On sufficient conditions of injectivity, development of a numerical test via interval analysis[END_REF], can be easily adapted for testing this new definition. In this definition, A c 1 designs the complementary of A 1 in A. Definition 2.4. Consider a function f : A → B and any set A 1 ⊆ A. The function f is said to be a restricted-partial injection of A 1 over A c 1 or a (A 1 , A)-restricted injection if :

∀a 1 ∈ A 1 , ∀a ∈ A c 1 , f (a 1 ) = f (a).
In section 3, the algorithm proposed in [START_REF] Lagrange | Injectivity analysis using interval analysis : Application to structural identifiability[END_REF] is taken again and adapted to determine the non µ-SM-identifiable parameters domains.

Method for analyzing the (µ-)SM-identifiability of a connected set

In [START_REF] Jauberthie | Set-membership identifiability: definitions and analysis[END_REF][START_REF] Verdière | Set-membership identifiability of nonlinear models[END_REF], we have proposed two methods for testing the global (µ-)SM-identifiability definition: the first one is based on the Taylor Series approach, the second one on differential algebra. The latter presents many advantages. The first one is that it reduces the SM-identifiability study of the model to the study of a real rational function whose range is the parameter domain of the model. The second one is that this approach leads to a numerical method given a first estimation of the parameters of the model ( [START_REF] Verdière | Identifiability and estimation of pharmacokinetic parameters of ligands of macrophage mannose receptor[END_REF][START_REF] Verdière | Set-membership identifiability of nonlinear models[END_REF]). In [START_REF] Jauberthie | Set-membership identifiability: definitions and analysis[END_REF][START_REF] Verdière | Set-membership identifiability of nonlinear models[END_REF], equivalent conditions between the definition of µ-SM-identifiability and the injectivity of a real rational function has been proved and is recalled below. In this paper, we present an equivalent condition between the definition of SM-identiability and the restricted-partial injectivity of the same real rational function.

The differential algebra [START_REF] Kolchin | Differential algebra and algebraic groups[END_REF] permits to obtain differential polynomials linking outputs, inputs and parameters which can be expressed as [START_REF] Verdière | Identifiability and estimation of pharmacokinetic parameters of ligands of macrophage mannose receptor[END_REF]:

w i (y, u, p) = m i 0 (y, u) + ni k=1 θ i k (p)m i k (y, u), i = 1, . . . , m, (2.2)
where

(θ i k (p)) 1≤k≤ni are rational in p, θ i u = θ i v (u = v), (m i k (y, u)) 0≤k≤ni
are differential polynomials with respect to y, u and m i 0 (y, u) ≡ 0.

{θ i k (p) 1≤k≤ni } is called the exhaustive summary of w i .
The size of the system is the number of outputs. For the time being, we suppose that i = 1, that is there is one output and n 1 = n, w 1 = w, m 1 k (y, u) = m k (y, u). The case of fewer outputs will be seen at the end of this section.

Consider t + 0 the right limit of t 0 2 and l the higher order derivative of y in (2.2). Hereafter, ∆w(y, u) will design the functional determinant formed from the {m k (y, u)} 1≤k≤n . The following theorem permits to obtain necessary and sufficient conditions for having global SM-identifiability or µ-SM-identifiability.

Theorem 2.5. Assume that the functional determinant ∆w(y, u) is not identically equal to zero. Consider P * a connected subset of U P . If the function φ : p = (p 1 , . . . , p p ) → (θ 1 (p), . . . , θ n (p), y(t + 0 , p), . . . , y (l-1) (t + 0 , p)) is (P * , U P )restricted injective (resp. (P * , U P )-injective) then P * is globally SM-identifiable (resp. µ-SM-identifiable).

In the two cases, if the coefficient of y (l) in (2.2) is not equal to 0 at t 0 , then the reciprocal is valid 3 .

Proof -The proof is an easily extension of the one proposed in [START_REF] Jauberthie | Set-membership identifiability: definitions and analysis[END_REF]. Remark -In the case of m outputs, the procedure is the following. For each of the m obtained differential polynomials w i (y, u, p), the functional determinant is evaluated. If it is not identically equal to zero, the associated exhaustive summary is added to the image of the function φ whose the (restricted-)partial injectivity has to be studied.

Testing the SM-identifiability

The steps for proving the SM-identifiability's definitions are summed up below.

1. Finding the differential polynomials w i in using, for example, the package DifferentialAlgebra of Maple. 2. Evaluate the functional determinants and construct the function φ.

Verify the (restricted-)partial injectivity of the function φ.

The proposed method consists in reducing the study of SM-identifiability's definitions to the study of the (restricted)-partial injectivity of the function φ whose range corresponds to the admissible parameter set of the model. The domains on which the function is restricted-partial injective but not partial injective will correspond to global SM-identifiable sets and domains on which the function φ is partially injective will correspond to µ-SM-identifiable sets.

In the following section, an algorithm to determine the SM-identifiable sets in the parameter space is proposed. It is based on interval analysis tools whose main advantage is that it guarantees numerical solutions provided as sets.

2 t + 0 is considered to ensure the existence of derivatives. 3 When initial conditions are not considered, the function φ becomes φ : p = (p 1 , . . . , pp) → (θ 1 (p), . . . , θn(p)) and the reciprocal of the theorem is not yet valid.

Determination of the (µ-)SM-identifiable sets

The proposed algorithms find the connected sets in the range of a continuous differentiable function such that, for each of them, their image is distinct to the one obtained from their complementary in the parameter space. Afterwards, they will be named SM-identifiable sets to be put in correspondence with the SM-identifiable sets defined from the model Γ. It takes back the algorithm of Lagrange et al. ([8] [9]) called ITVIA (for Injectivity Test Via Interval Analysis) which partitions the domain of a differentiable function into two domains: an undetermined domain not proved partially injective and a partial injective domain. The latter corresponds to a list of µ-SM-identifiable sets. To complete the partition, we have to find the SM-identifiable sets in the undetermined domain.

Denote f : R n → R m (f ∈ C 1 ) a differentiable vector function defined over a given n-dimensional box [x] ∈ IR n , where IR is the set of all intervals and [f ] its inclusion function [START_REF] Jaulin | Applied Interval Analysis, with examples in parameter and state estimation, Robust control and robotics[END_REF]. Regarding to [START_REF] Chabert | Contractor programming[END_REF], a subpaving of [x] is a set of non-overlapping boxes included in [x]. A subpaving can be considered either as a collection (list) of boxes K = [x] (1) , [x] (2) , ... or as a union (1) [x] (2) .... Hence, a subpaving can either be viewed as discrete subset of IR n or as a convex subset of R n . Subpavings permit to approximate convex sets with arbitrary precision. A regular subpaving is a subpaving generated by successive bisections and it can be easily represented by a binary tree.

K = [x]
After obtaining a paving of µ-SM-identifiable sets by the algorithm ITVIA, we determine its complementary and, to complete the work, we proceed as in Image evaluation [START_REF] Jaulin | Applied Interval Analysis, with examples in parameter and state estimation, Robust control and robotics[END_REF]. Firstly, we consider a subpaving with boxes of width smaller than (mince step). Secondly, we find SM-identifiable sets as a list of non-overlapped boxes (evaluate step). Finally, we transform this list to a regular subpaving (regularize step). We begin by the case m = 1 and then we generalize our algorithm to the multidimensional case.

Case of unidimensional functions

Let {P} a collection (list) of non overlapped boxes included in [x], we introduce the following notations. N corresponds to the length of the list {P}. max ). Then, i Max and i Min are the minimal corresponding indexes. Finally, is the maximal diameter of the boxes considered in the list {P}. Mathematically, these definitions can be rewritten under the form:

• length({P}) = N , • max [f ],{P} = {max([f ]({P} i )), 1 ≤ i ≤ N } ∈ R N , • min [f ],{P} = {min([f ]({P} i )), 1 ≤ i ≤ N } ∈ R N , • Max [f ],{P} + = max(max [f ],{P} ) ∈ R, • Max [f ],{P} -= max(min [f ],{P} ) ∈ R, • Min [f ],{P} -= min(min [f ],{P} ) ∈ R, • Min [f ],{P} + = min(max [f ],{P} ) ∈ R,
• = w({P}) i=1,...,N = max 1≤i≤n (w({P} i )) where w represents paving width. 4 Example 1:

In order to illustrate these notations, the following function f defined by: f (x) = sin(10(x -0.1) [START_REF] Jaulin | Applied Interval Analysis, with examples in parameter and state estimation, Robust control and robotics[END_REF] for the definition of connected box), we add it to {L} and we continue. If there exists a box [x] temp not connected with one of the box of {L}, one gets the box containing the second maximum of the function f . The second maximum is named disconnected afterwards, an example is shown at with its complementary image. Since {L}

+ [f ]
,{P},a is supposed to be connected, it constitutes a SMidentifiable set. Similarly, one gets the following proposition. Proposition 3.5.

In (3.3), for a = SMin [f ],{P} -, if the list {L} - [f ]
,{P},a is a connected set, then it is a SM-identifiable set. We denote it by {L} Algorithm 2 is used to find the connected sets verifying Proposition 3.4. In this algorithm, two lists {L 1 } and {L 2 } verify {L 1 } = {L 2 }, if they are composed of the same elements even in different order. Moreover, the notation [x a ] designates the component {P} i of the list {P} for which min(

[f ]({P} i )) = a. Algorithm 2 {L} + SM I,[f ],{P} = ConnectedSet + ({P} , max [f ],{P} , min [f ],{P} ) Input:{P}, max [f ],{P} , min [f ],{P} Output: {L} + SM I,[f ],{P} Initialization: a = SMax [f ],{P} + , {L 1 } = ∅, {L 2 } = {L} + [f ],{P},a while {L 2 } = {L 1 } do {L 1 } = {L 2 } a = SMax - [f ],{L1} (in using Algorithm 1) {L 2 } = {L} + [f ],{L1},a if a is a disconnected maximum then {L 2 } = {L 2 } -[x a ] end if end while {L} + SM I,[f ],{P} = {L 2 }
Another algorithm can be deduced in substituting + and max byand min in order to find {L} -SM I,[f ],{P} .

3.1.3. Final algorithms. Algorithms 3 and 4 are the final algorithms which determine the SMidentifiable sets in the definition domain of a function and take again the previous Algorithms 1 and 2. {U} SM I,[f ],{P} corresponds to the set of the two lists {L} As it can be seen, these algorithms are different only in the order of treating the partitions {U + } and {U -}. In addition, line 1, consisting in the subtraction of two lists, prevents the apparition of the same boxes in {U + } and {U -}.

+ SM I,[f ],{P} and {L} - SM I,[f ],{P} . Algorithm 3 {U} SM I,[f ],{P} =SMISet1([f ], {P}) Input: [f ], {P} Output: {U} SM I,[f ],{P} Initialization: {U + } = ∅, {U -} = ∅, {L 1 } = {P} {U + } = {L} + SM I,[f ],{L1} (in using Algorithm 2) 1 : {L 2 } = {L 1 } -{U + } {U -} = {L} - SM I,[f ],{L2} (

The case of multidimensional functions

In this section, we consider the case m > 1. We suppose that an interval vector function

[f ] = [ [f 1 ], ..., [f m ] ]
f i ] and {L 1 } ) {L 2 } = {L 1 } -{U + } i {U -} i = {L} - SM I,[fi],{L2} (in using Algorithm 2 in substituting [f ] and {P} by [f i ] and {L 2 } ) {L 1 } = {L 2 } -{U -} i end for {U} SM I,[f ],{P} = {{U + } , {U -}}
However, as in the case of m = 1, other algorithms can be considered, for example the following algorithm:

Input: [f ], {P} Output: {U} SM I,[f ],{P} Initialization: {U + } = ∅, {U -} = ∅, {L 2 } = {P} for i=1:m do {U -} i = {L} - SM I,[fi],{L2} {L 1 } = {L 2 } -{U -} i {U + } i = {L} + SM I,[fi],{L1} {L 2 } = {L 1 } -{U + } i end for {U} SM I,[f ],{P} = {{U + } , {U -}}
Changing the order of the interval functions, [f i ], i = 1, ..., m and also the order of computing partitions {U -} i and {U + } i , altogether, there can be (2m)! algorithms like Algorithms 5 and 6.

Regularization

In [START_REF] Jaulin | Applied Interval Analysis, with examples in parameter and state estimation, Robust control and robotics[END_REF], in order to evaluate the image of a function at a domain set, three main sequential steps, mince, evaluate and regularize are introduced. In Sections 3.3 and 3.4, the two first steps, the mince and the evaluate ones have been explained. In this section, the last step, i.e. the regularization is realized. Figure 4 illustrates an example of the regularization step based on the algorithm presented in [START_REF] Jaulin | Applied Interval Analysis, with examples in parameter and state estimation, Robust control and robotics[END_REF]. This algorithm can also be applied to a list composing of overlapped or non overlapped boxes. However, as it can be seen in Figure 4, the initial list, here consisting of overlapped boxes in Figure 4 (a), is inside the regular subpaving of the Figure 4 (b), or, in other words, the subpaving is an external regular subpaving. If we apply this algorithm to a list of SM-identifiable sets, some parts of the initial domain which are not SM-identifiable will be included in the obtained regular subpaving. To cope with this problem, the Algorithm 7 is proposed to find a regular subpaving inside a list. However, the list must only consist of non overlapped boxes which is fortunately the case of a list of SM-identifiable sets found by the method presented in the previous sections. In this algorithm, the function named volume and applied to an interval box

[x] = [[x 1 ] [x 2 ] ...[x n ]] T ∈ IR n is defined as: volume([x]) = (x 1 -x 1 ) × (x 2 -x 2 ) × ... × (x n -x n ) with [x i ] = [x i , x i ] for i = 1, . . . , n.
To explain how this algorithm works, we consider the simple case n = 2. In this case, since [x] is a rectangle, the function volume consists in evaluating its surface. Beginning by the convex In that case, the function φ is defined by φ :

R 2 → R 2 , φ(p) = p 2 1 -p 2 2
3. Algorithm 5 is used with a regular paving such that its resolution is = 0.05. The results are illustrated in Figure 6 after regularization. Light gray, dark gray and black color boxes correspond to SM-identifiable sets but which are not µ SM-identifiable. The state gray color box represents the µ-SM-identifiable sets. With Algorithm 6, we do not obtain such SM-identifiable sets. 3. Choosing = 0.01 and Algorithm 5, the results illustrated in Figure 7 are obtained.

Conclusion

In this paper, we propose a guaranteed method in order to determine, in the parameter space of a model, the (µ-)SM-identifiable sets, each of them corresponding to different output behaviors from their complementary. For doing this, differential algebra tools have been used for obtaining differential polynomials. From them, a real function depending only on the parameters has been constructed. We show that the domains on which this function is partially injective correspond to µ-SM-identifiable sets, otherwise if the function verifies only the restricted-partial injectivity, the corresponding domain determines SM-identifiable sets. An algorithm based on interval analysis tools has been proposed to determine the SM-identifiable sets in the function's range, the latter corresponding to parameter space. Finally, Two examples have been provided in order to study the SM-identifibility of the model.

  [f ],{P} (resp. min [f ],{P} ) denotes the sequence of the maximal value (resp. the minimal value) of [f ] on each box of {P} and the two following, Max [f ],{P} + and Min [f ],{P} + (resp. Max [f ],{P} -and Min [f ],{P} -), the maximal and minimal values of the sequence max [f ],{P} (resp. min [f ],{P}

Figure 2 .FIGURE 3 .

 23 FIGURE 3. Connected and non connected sets.

-

  SM I,[f ],{P} .

Figure 3

 3 Figure 3 providing by Equation (3.1) illustrates the importance of the connected property verification. In the first Figure, one gets a = SMax [f ],{P} + = max([f ]({P} 7 )), {L} + [f ],{P},a = {{P} 2 , {P} 3 , {P} 4 , {P} 5 } and it can be easily seen that {L} + [f ],{P},a is the list {L} + SM I,[f ],{P} since it is connected. Notice that SMax [f ],{P} -is the lower bound of a sub-list. However, in the second Figure, {L} + [f ],{P},a = {{P} 2 , {P} 4 , {P} 5 } is no more connected and can

  in using Algorithm 2 and in substituting + and max byand min){U} SM I,[f ],{P} = {{U + } , {U -}} Algorithm 4 {U} SM I,[f ],{P} =SMISet2([f ], {P}) Input: [f ], {P} Output: {U} SM I,[f ],{P} Initialization: {U + } = ∅, {U -} = ∅, {L 2 } = {P} {U -} = {L} -SM I,[f ],{L2}(in using Algorithm 2 in substituting + and max byand min) 1 :{L 1 } = {L 2 } -{U -} {U + } = {L} + SM I,[f ],{L1} (in using Algorithm 2) {U} SM I,[f ],{P} = {{U + } , {U -}}

  and a subpaving {P} of [x] are available. The main idea of Algorithm 5 consists in applying Algorithms 3 or 4 consecutively on each component of the vector [f ]. For example, Algorithm 3 is applied with [f 1 ] and the list {P}. Suppose that it gives the list {U} SM I,[f1],{P} of SM-identifiable sets. Algorithm 3 is then reapplied with [f 2 ] and {P} -{U } SM I,[f1],{P} . The process goes on until the last component of the function [f ]. Based on Algorithms 3 and 4, Agorithm 5 has been implemented. Algorithm 5 {U} SM I,[f ],{P} =SMISetN1([f ], {P}) Input: [f ], {P} Output: {U} SM I,[f ],{P} Initialization: {U + } = ∅, {U -} = ∅, {L 1 } = {P} for i=1:m do {U + } i = {L} + SM I,[fi],{L1} (in using Algorithm 2 in substituting [f ] and {P} by [

2 FIGURE 4 .

 24 FIGURE 4. (a): list of overlapped boxes , (b): external regular subpaving

FIGURE 5 .

 5 FIGURE 5. (a): list of non overlapped boxes, (b): internal regular subpaving (c): external regular subpaving

2 FIGURE 6 .

 26 FIGURE 6. Example 2: State gray: µ-SM-identifiable set, light Grays and Black: SM-identifiable sets.

Example 3 : 3 ) 2 .

 332 Consider the following example:   ẋ1 (t) = x 1 (t) + x 2 (t) + u(t), x 1 (0) = p 1 (1 + sin(p 1 ) -p 2 sin(p 1 )) + p 2 cos(p 1 ), ẋ2 (t) = -x 2 1 (t) -((1 -p 2 )p 1 cos(p 1 ) -p 2 sin(p 1 ) -2p 2 )x 1 x 2 , x 2 (0) = 0, y = x 1 .(4.2) 1. By setting c 1 = cos(p 1 ) and c 2 = sin(p 1 ), the Rosenfeld-Groebner algorithm, implemented in the package DifferentialAlgebra of Maple gives the following differential polynomial:w(y, u, p) = -u -ẏ + ÿ + y 2 + ((1 -p 2 )p 1 cos(p 1 ) -p 2 sin(p 1 ) -2p 2 )( ẏy -uy -y 2 ). (4.Clearly the associated functional determinant is not identically equal to zero. Suppose too that u(0) = 0, then the following function φ : R 2 → R 2 can be considered:φ(p) = (1 -p 2 )p 1 cos(p 1 ) -p 2 sin(p 1 ) -2p 2 p 1 (1 + sin(p 1 ) -p 2 sin(p 1 )) + p 2 cos(p 1 ) .

2 FIGURE 7 .

 27 FIGURE 7. Example 3: State gray: µ-SM-identifiable sets, Black: SM-identifiable sets.

  For example, consider that M = max([f ]({P})). In this algorithm, we find Max [f ],{P} + , that is the maximum value of the list max([f ]({P})) and we initialize a list named {L} by the corresponding box of Max [f ],{P} + in {P}. Then, we find the next maximum value of max([f ]({P})). If its corresponding box [x] temp in {P} is connected to one of the component of {L}, that is it intersects one of the box of {L} (see

		f(x)=sin(10(x-0.1) 2 )/x, dashed: max([f](.)), solid: min([f](.))
	2.5	Max +
	2		Max -
	1.5
						SMax +
	1	SMax -
	0.5
	0
	0.2 SMin SMax [f ],{P} 0.4 0.6 0.8 1 -1.5 -1 -0.5 Min + + (see Figure 2).	x	1.2 Min -	1.4	1.6 SMin -	1.8	2
	While Max [f ],{P}	+ and Max [f ],{P}	-can be directly evaluated from max([f ]({P})), we pro-
	pose Algorithm 1 for evaluating SMax [f ],{P}	+ or SMax [f ],{P}	-. There, SMax [f ],{P}	+ and
	SMax [f ],{P}

2 )/x, where [x] ∈ [0.1, 2], (3.1) is represented at Figure 2. {P} consists of equal boxes of width = 0.05 in [0.1, 2]. +

FIGURE 2. Illustration of the notations, to simplify, the subscript {[f ], {P}} is omitted.

3.1.1. Second maximum and second minimum. The aim of the two following sections is to isolate the domain around the maximum (resp. minimum) of f which constitutes a SM-identifiable set as the interval [0.2; 0.6] at Figure

2

. For doing this, one has to determine first the second maximum of f (resp. minimum), in particular the lower and the upper values of the box containing this second maximum (resp. minimum). The lower value will be denoted SMax [f ],{P} -and the upper value, -are computed in replacing M ∈ R N by max([f ]({P})) and min([f ]({P})), respectively.

This assumption is important at section

2.3 for the use of differential algebra.

Acknowledgment

This work was supported by the French National Research Agency (ANR) in the framework of the project ANR-11-INSE-006 (MAGIC-SPS).

This work was supported by the French National Research Agency (ANR) in the framework of the project ANR-11-INSE-006 (MAGIC-SPS).

Relation between partitioning and second maximum, minimum. For the final step, we construct the SM-identifiable set defined around the maximum or minimum value of f . Indeed, when we search the second maximum or minimum, we obtain the list of boxes {L} whose image, by the function f , is not connected with its complementary. However, we are not insure to obtain a connected set. For being convinced, see for example the second figure of

The notion of connected lists formed of intervals is introduced before giving a sufficient condition for having a connected SM-identifiable set in [x]. Definition 3.3. A list {L} formed of intervals is said connected if the union of its components is an interval of R. Proposition 3.4.

,{P},a is a connected set, then it is a SM-identifiable set. We denote it by {L}

Proof -The list {L} + [f ],{P},a corresponds, by construction, to the list of boxes not having a common intersection with its complementary. Concretely, the image {L}

union of the paving {P}, we construct a list by consecutive bisections as explained in that follows. At each iteration, we get the first element of the last list, i.e. [x], and we evaluate its surface, i.e. v x . We compute also the surface of this element which is covered by the paving, i.e. v t . If v t exists, and if its difference with v x is very small comparing to v x , then [x] is accepted and is added to the response list {P R }. If v x exists but the difference is large and [x] is wide compared to the precision , then we bisect [x] and we add the results to the end of the last list. The new list is produced. Recall that the elements of {P} do not intersect and consequently {P} i ∩ [x] do not have intersection with {P} j ∩ [x] for all i, j from 1 to length({P}), i = j. This property enables us to compare sum of the surfaces {P} i ∩ [x], i.e. v t , with the surface of [x], i.e. v x .

Figure 5 shows the difference between internal and external regular subpavings found for a list. As it can be seen, the list consists of non overlapped boxes (Figure 8 (a)), the internal regular subpaving is inside the list (Figure 8

Examples


Hélène Frankowska
email: frankowska@math.jussieu.fr

Daniela Tonon
email: dtonon@imperial.ac.uk

† Hélène

The Goh necessary optimality conditions for the Mayer problem with control constraints

The well known Goh second order necessary optimality conditions in optimal control theory concern singular optimal controls taking values in the interior of a set of controls U. In this paper we investigate these conditions for the Mayer problem when U is a convex polytope or a closed subset of class C 2 for an integrable optimal control ū(•) that may take values in the boundary of U. This is indeed a frequent situation in optimal control and for this reason the understanding of this issue is crucial for the theory of second order optimality conditions. Applying the Goh transformation we derive necessary conditions on tangent subspace to U at ū(t) for almost all t's.

In the presence of an endpoint constraint, if the Mayer problem is calm, then similar second order necessary optimality conditions are satisfied whenever the maximum principle is abnormal. If it is normal, then analogous results hold true on some smaller subspaces.

Introduction

Let U ⊂ R m , K ⊂ R n be closed, x 0 ∈ R n , f : [0, 1] × R n × R m → R n
and ϕ : R n → R be sufficiently regular. Consider the Mayer optimal control problem min ϕ(x(1)) :

x(•) ∈ S [0,1] (x 0 ), x(1) ∈ K , (1)
where S [0,1] (x 0 ) is the set of all absolutely continuous solutions of the control system

x (t) = f (t, x(t), u(t)), u(t) ∈ U, x(0) = x 0 . ( 2
)
Define the Hamiltonian H :

[0, 1] × R n × R n × R m → R : H (t, x, p, u) = p, f (t, x, u) .
We first consider the Mayer problem without endpoint constraints:

min ϕ(x(1)) : x(•) ∈ S [0,1] (x 0 ) . (3)
Let ( x(•), ū(•)) be a strong local minimizer of [START_REF] Schättler | Geometric optimal control[END_REF]. By the maximum principle, for an absolutely continuous p : [0, 1] → R n satisfying -p (t) = H x (t, x(t), p(t), ū(t)) for a.e. t ∈ [0, 1] (

and -p(1) = ∇ϕ( x( 1)) [START_REF] Krener | The high order maximal principle[END_REF] the following maximum principle holds true for a.e. t : p(t), f (t, x(t), ū(t)) = max u∈U p(t), f (t, x(t), u) . [START_REF] Krener | The high order maximal principle and its application to singular extremals[END_REF] Thus for a.e. t ∈ [0, 1] such that the optimal control ū(t) lies in the interior of U, H u [t] = 0 and the Legendre-Clebsch condition H uu [t] ≤ 0 is satisfied, where [t] := (t, x(t), p(t), ū(t)). If ū(•) is singular, in the sense that H uu [•] = 0, the Legendre-Clebsch condition doesn't bring any additional information about optimal controls. Different, but related, definitions of singular controls and discussions on their importance can be found for instance in [START_REF] Bonnard | Singular trajectories and their role in control theory[END_REF][START_REF] Goh | Optimal singular rocket and aircraft trajectories[END_REF][START_REF] Schättler | Geometric optimal control[END_REF], and the references therein.

Assume that H uu [t] = 0 a.e. in [t 1 ,t 2 ]. According to [START_REF] Goh | Necessary conditions for singular extremals involving multiple control variables[END_REF], if U is an open convex polytope, then the matrix

H xu [t] f u [t] is symmetric and R(t) := 1 2 f u [t] * H xx [t] f u [t] - 1 2 d dt H xu [t] f u [t] -H xu [t] f x [t] f u [t] + 1 2 H xu [t] d dt f u [t]
is negative definite for a.e. t ∈ [t 1 ,t 2 ]. The first property, known in the literature as the Goh condition, in the case of affine control systems reduces to a condition on the Lie brackets of the flux functions, as we recall in Remark 3.2. The second property can be seen as a generalized Legendre-Clebsch condition.

To deal with closed control constraints Goh also suggested to use variations of the optimal control which are equal to zero when ū(t) ∈ ∂U. However, this can't be applied to recover second-order optimality conditions when ū(t) ∈ ∂U almost everywhere.

We present here a different approach, based on non-zero variations of the optimal control even when ū(•) is a boundary control.

For every u ∈ U let us consider the maximal subspace P u contained in the tangent cone to U at u.

Below, Theorems 4.1 and 5.2 provide pointwise second order necessary optimality conditions respectively in the case when U is a convex polytope or a subset of class C 2 for a singular optimal control ū(•) being merely integrable. Considering H xu [t] f u [t] and R(t) as bilinear operators on the subspace P ū(t) , in both cases, we still have symmetry of

H xu [t] f u [t].
In the former case we also get the negative definiteness of R(t) on P ū(t) .

In Corollary 4.2 we prove that if U is a closed convex polytope in R m and f is affine with respect to u:

f (x, u) = f 0 (x) + g(x)u, (7)
with g(•) ∈ C 2 (R n ; M(n × m) and if g k (x) denotes the k-th column of the matrix g(x), then the matrix valued mapping Θ : [0, 1] → M(m × m), given by

θ k,l (t) := p(t), [g k , g l ]( x(t)) ∀ k, l ∈ {1, . . . , m},
enjoys the following property: for a.e. t ∈ [0, 1], the subspace Θ(t) P ū(t) (the image of P ū(t) by the linear operator Θ(t)) is orthogonal to the subspace P ū(t) . In particular we recover the well known Goh condition in the geometric optimal control: Θ(t) = 0 for a.e. t such that ū(t) ∈ Int U.

For the Mayer problem with the endpoint constraint K, we assume calmness at x(•). This allows to reduce the problem to the one without endpoint constraints but involving a nondifferentiable cost function and so the classical first and second order necessary conditions are not satisfied. Still we are able to derive second order conditions through a variational approach similar to the one without end point constraints.

Since the very first work by Goh, [START_REF] Goh | Necessary conditions for singular extremals involving multiple control variables[END_REF], the Goh condition was investigated by many researchers in the geometric optimal control theory following fundamental works by Krener, cf. [START_REF] Krener | The high order maximal principle[END_REF] and [START_REF] Krener | The high order maximal principle and its application to singular extremals[END_REF], on the Goh condition (and its higher order extensions) when the optimal control is piecewise C ∞ . His approach was taken over by several authors and the Goh condition became mostly known in the form involving Lie brackets, as a second-order necessary optimality condition for the abnormal maximum principle. See for instance [START_REF] Agrachev | Control theory from the geometric viewpoint[END_REF][START_REF] Agrachev | The principle of second order optimality for time-optimal problems[END_REF] and [3, pp.314-319].

Several authors made an extensive use of the Goh transformation to recover necessary and sufficient second order optimality conditions in an integral and a pointwise form, see for instance [START_REF] Aronna | Quadratic order conditions for bang-singular extremals[END_REF], [START_REF] Aronna | Singular solutions in optimal control: second order conditions and a shooting algorithm[END_REF] and the bibliographies contained therein. However, in these latter works many additional structural restrictions on the control system, control sets and also on the optimal control are imposed, as for instance in Assumptions 1, 2 from [START_REF] Aronna | Quadratic order conditions for bang-singular extremals[END_REF] concerning the structure optimal controls .

The outline of this paper is as follows. In Section 2 we introduce some notations and provide some preliminary results. Section 3 is devoted to the second order necessary optimality conditions for an arbitrary set U and any integrable optimal control at times when it belongs to the interior of U. Section 4 deals with the case when U is a convex polytope and Section 5 with the case when U is of class C 2 . The last section adresses the second order necessary optimality conditions for the Mayer problem involving endpoint constraints.

Most of results presented here are given without proofs, because they are quite long and technical and will be published elsewhere.

Preliminaries and notations

The norm in R n is denoted by • and the inner product by •, • . Let B := {x ∈ R n : x < 1} be the open unit ball in R n . For a nonempty subset K of R n , ∂ K denotes its boundary, Int K its interior and K c its complement. The distance function

d K : R n → R + is defined by d K (x) := inf y∈K x -y for all x ∈ R n and the oriented distance b K : R n → R from K = R n by b K (x) := d K (x)-d K c (x) for all x ∈ R n . We set b K (•) = 0 if K = R n .
Finally, by M(k × r) we denote the space of matrices having k lines and r columns. Definition 2.1 Let C : R m → R m be a linear mapping and F be a subspace of R m . We say that C is symmetric on F if Cu, v = u,Cv for all u, v ∈ F and that C is negative definite on F if Cu, u ≤ 0 for all u ∈ F.

Let U be a nonempty closed subset of R m . For u 0 ∈ U the adjacent tangent cone to U at u 0 is defined by

T U (u 0 ) := u ∈ R n : lim h→0+ d U (u 0 + hu) h = 0 .
The second-order adjacent subset to U at (u 0 , u) ∈ R m × R m is the set defined by

T (2) U (u 0 , u) := v ∈ R n : lim h→0+ d U (u 0 + hu + h 2 v) h 2 = 0 .
For any u ∈ U we denote by P u the largest subspace contained in T U (u).

Let f : [0, 1] × R n × R m → R n , f (t, x, u) be the deriva- tive of f (t, •, •) at (x, u), f x (t,
x, u) and f u (t, x, u) be the partial derivatives of f with respect to x and u; f (t, x, u) denotes the Hessian of f (t, •, •) at (x, u), while f xx , f xu and f uu denote the second-order partial derivative of f with respect to x and u. Fix a trajectory-control pair ( x(•), ū(•)) of ( 2) and set

f [t] = f (t, x(t), ū(t)). Let f u [t], f x [t] and f [t]
be defined in a similar way. Throughout the whole paper we impose the following assumptions : For some δ > 0

(a) f (t, •, •) ∈ C 2 on x(t) + δ B × ū(t) + δ B , ∀ t; (b) f u [•] , f [•] , f [•] ∈ L 1 ([0, 1]; R + ) and ∃ a(•) ∈ L 1 ([0, 1]; R + ) such that the mappings f (t, •, •) and f (t, •, •) are a(t)-Lipschitz on ( x(t) + δ B) × ( ū(t) + δ B) for a.e. t ∈ [0, 1]; (c) f (t, x,U) is closed for all t ∈ [0, 1], x ∈ x(t) + δ B; (d) ϕ(•) ∈ C 2 on x(1) + δ B. Recall that ( x(•), ū(•)) is a strong local minimizer for (3) (resp. (1)) if there exists η > 0 such that x(•) minimizes ϕ over all x(•) ∈ S [0,1] (x 0 ) (resp. over all x(•) ∈ S [0,1] (x 0 ) with x(1) ∈ K) satisfying x(•) -x(•) ∞ ≤ η.
Remark 2.2 Results contained in Sections 3-5 below are stated only for strong local minimizers. However they are also valid for the so called weak local minimizers. In contrast, Section 6 concerns only strong local minimizers. Not to switch between these different notions, in this paper we investigate only strong local minimizers.

Theorem 2.3 (Maximum Principle, [START_REF] Pontryagin | The mathematical theory of optimal processes[END_REF]) Let ( x(•), ū(•)) be a strong local minimizer of the Mayer problem [START_REF] Schättler | Geometric optimal control[END_REF]. Then there exists an absolutely continuous p : [0, 1] → R n satisfying (4), ( 5) and [START_REF] Krener | The high order maximal principle and its application to singular extremals[END_REF]. Furthermore, for a.e. t ∈ [0, 1] and for any u ∈ P ū(t) we have

H u [t]u = 0. Define A := {t ∈ [0, 1] : ū(t) ∈ ∂U, H u [t] = 0} and U A = {u(•) ∈ L ∞ ([0, 1]; R m ) : ∃ c > 0, h 0 > 0 s.t. ∀ h ∈ [0, h 0 ], ∀ t ∈ A, d U ( ū(t) + hu(t)) ≤ ch 2 }.
Note that for any u(•) ∈ U A , u(t) ∈ T U ( ū(t)) a.e. in A.

(•), v(•) ∈ L ∞ ([0, 1]; R m ) such that u(•) ∈ U A and v(t) ∈ T (2)
U ( ū(t), u(t)) for a.e. t ∈ A, and for p(•) as in the maximum principle either

1 0 H u [t]u(t)dt < 0 or H u [t]u(t) = 0 a.e. in [0, 1] and 0 ≥ 1 0 H u [t]v(t) + 1 2 y(t) * H xx [t]y(t) +u(t) * H xu [t]y(t) + 1 2 u(t) * H uu [t]u(t))dt -1 2 y(1) * ϕ ( x(1))y(1) , (8)
where y(•) solves the linear system

y (t) = f x [t]y(t) + f u [t]u(t) y(0) = 0. ( 9
)
The proof follows the same lines as the one in [START_REF] Hoehener | Variational approach to second-order optimality conditions for control problems with pure state constraints[END_REF].

Second order conditions in the interior of control constraints

Consider a strong local minimizer ( x(•), ū(•)) of the Mayer problem (3) and p(•) as in the maximum principle of Theorem 2.3. Let the matrix valued mapping R(t) be defined as in the introduction.

Theorem 3.1 Assume that for some 0 ≤ t 1 < t 2 ≤ 1, ū(•) is integrable on [t 1 ,t 2 ] and H uu [t] = 0 a.e. in [t 1 ,t 2 ]; H xu [•], f u [•] are Lipschitz on [t 1 ,t 2 ]; H xx [•], f x [•] are essentially bounded on [t 1 ,t 2 ]. (10)
Then, for a.e. t ∈ [t 1 ,t 2 ] such that ū(t) ∈ Int U, the matrix

H xu [t] f u [t] is symmetric and R(t) ≤ 0.
The statement of the above theorem is proved for the Lebesgue points t of ( ū(•), R(•)) such that ū(t) ∈ Int U. The proof is long and technical and will be published elsewhere. It involves the Goh transformation and some advanced properties of Lebesgue points. Remark 3.2 If f is affine with respect to controls, the symmetry of the matrix H xu [t] f u [t] implies that the adjoint state is orthogonal to the Lie brackets of the flux. Indeed let f be as in [START_REF] Agrachev | Control theory from the geometric viewpoint[END_REF], where g(•) ∈ C 2 (R n ; M(n×m)) and U is a closed nonempty subset of R m and let g k (x) denote the k-th column of g(x). Then for all t such that H xu [t] f u [t] is symmetric and for any k, l ∈ {1, . . . , m},

0 = (H xu [t] f u [t] -f u [t] * H xu [t] * ) k,l = n ∑ j=1 p j (t) n ∑ s=1 ∂ g j,k ∂ x s [t]g s,l [t] - ∂ g j,l ∂ x s [t]g s,k [t] ,
where g j,k is the j-th element of g k . Hence for any k, l ∈ {1, . . . , m},

p(t), ∂ g k ∂ x [t]g l [t] - ∂ g l ∂ x [t]g k [t] = p(t), [g k , g l ][t] = 0.

Second order conditions on convex polytopes

A nonempty subset U ⊂ R m is called a convex polytope in R m if U is an intersection of a finite family of affine half-spaces of R m . Then the adjacent tangent cone T U (u 0 ) is equal to the tangent cone of convex analysis to U for any u 0 ∈ U. The proof follows the same lines as the proof of Theorem 3.1, but a change of base is involved and the control variations take values in the tangent to U at ū(t).

The following corollary generalizes the classical Goh condition to boundary controls. Corollary 4.2 Under all the assumptions of Theorem 4.1 suppose that f is as in [START_REF] Agrachev | Control theory from the geometric viewpoint[END_REF], where g(•) ∈ C 2 (R n ; M(n × m)) and U is a closed convex polytope in R m and let g k (x) denote the k-th column of g(x). Consider the matrix valued mapping Θ : [0, 1] → M(m × m), whose elements are

θ k,l (t) := p(t), [g k , g l ]( x(t)) ∀ k, l ∈ {1, . . . ,
[t] f u [t] -(H xu [t] f u [t]) * )b = 0, we get 0 = m ∑ k=1 m ∑ l=1 a k (H xu [t] f u [t] -f u [t] * H xu [t] * ) k,l b l . Consequently p(t), ∑ m k=1 ∑ m l=1 a k [g k , g l ]( x(t))b l = 0 and therefore 0 = ∑ m k=1 ∑ m l=1 a k p(t), [g k , g l ]( x(t)) b l = a, Θ(t)b , as claimed.

Second order conditions for smooth control constraints

In this section we consider a proper closed set U ⊂ R m of class C 2 , see for instance [START_REF] Delfour | Shapes and geometries[END_REF]. For a strong local minimizer ( x(•), ū(•)) of the Mayer problem (3) define the following matrix valued mapping C(•) on [0, 1]

C(t) = b U ( ū(t)) if ū(t) ∈ ∂U 0 otherwise Remark 5.1 If U is of class C 2 , then a) ∀ u 0 ∈ ∂U there exists a neighborhood W (u 0 ) of u 0 such that the oriented distance b U (•) ∈ C 2 (W (u 0 )), see for instance Theorem 4.3 in [13]; b) ∀ u 0 ∈ ∂U, T U (u 0 ) = {v ∈ R m : ∇b U (u 0 ), v ≤ 0} and the subspace P u 0 := {v ∈ R m : ∇b U (u 0 ), v = 0} = ∂ T U (u 0 )
has the dimension equal to m -1.

c) ∀ u 0 ∈ ∂U, ∇b U (u 0 ) is the unit outward normal to U at u 0 and, applying the inverse mapping theorem to b U , it is not difficult to show that for any u

∈ ∂ T U (u 0 ) we have -1 2 b U (u 0 )u, u ∇b U (u 0 ) ∈ T (2)
U (u 0 , u).

For a trajectory-control pair ( x(•), ū(•)) and for every t ∈ [0, 1] such that ∇b U ( ū(t)) does exist, define the m-1 dimensional subspace Assume that for some

S(t) := {u ∈ R m : ∇b u ( ū(t)), u = 0}.
0 ≤ t 1 < t 2 ≤ 1, ū(•) is integrable on [t 1 ,t 2 ].
Then for a.e. t ∈ [t 1 ,t 2 ] such that ū(t) ∈ ∂U and for some δ t > 0

u * (H uu [s] -H u [t] b U ( ū(s))) u = 0
for all u ∈ S(s) and for a.e. s ∈ [tδ t ,t

+ δ t ] satisfying || ū(t) -ū(s)|| ≤ δ t ; H xu [•], f u [•] are Lipschitz on [t -δ t ,t + δ t ]; H xx [•], f x [•] are essentially bounded on [t -δ t ,t + δ t ],
we have H xu [t] f u [t] is symmetric on the subspace P ū(t) .

Mayer's problem with the endpoint constraints

Let K be a proper closed nonempty subset of R n of class C 2 . Then b K (•) ∈ C 2 on a neighborhood of ∂ K, cf. [START_REF] Delfour | Shapes and geometries[END_REF]. We study here the Mayer optimal control problem (1) involving the end point constraints.

Let ( x(•), ū(•)) be a strong local minimizer of the above problem and P ū(t) be defined as in Section 2. If x(1) belongs to Int K, then the previously developed second order variational analysis does apply as though there were no endpoint constraints. This is due to the fact that small perturbations of ū(•) lead to feasible trajectories. For this reason we have to investigate only the case x(1) ∈ ∂ K. We deduce, for instance from [START_REF] Vinter | Optimal control[END_REF], the following maximum principle. Theorem 6.1 Consider a strong local minimizer ( x(•), ū(•)) of (1) with x(1) ∈ ∂ K. Then there exists an absolutely continuous function p : [0, 1] → R n and λ ∈ {0, 1}, µ ≥ 0, satisfying λ + µ > 0, (4), ( 6) and the transversality condition

-p(1) = λ ∇ϕ( x(1)) + µ∇b K ( x(1)) .
The above maximum principle may be abnormal, that is λ = 0. In this case p(1) may be taken equal to -∇b K ( x(1)). Lemma 6.2 Consider a strong local minimizer ( x(•), ū(•)) of (1) with x(1) ∈ ∂ K for which an abnormal maximum principle of Theorem 6.1 holds true for some µ > 0. Then for every measurable bounded selection u(t) ∈ P ū(t) a.e. in [0, 1], the solution y(•) of ( 9)

satisfies ∇b K ( x(1)), y(1) = 0 and p(t), f u [t]u(t) = 0 a.e. in [0, 1].
If the maximum principle is abnormal, then the cost is not involved into the first order optimality conditions. To apply the variational approach developed in the previous sections, we use the exact penalization which allows to deal with a normal maximum principle.

Definition 6.3 The Mayer problem (1) is called calm at x(•) if there exist ε > 0 and M > 0 such that ∀ ν ∈ (0, ε) and x(•) ∈ S [0,1] (x 0 ) satisfying d K (x(1)) ≤ ν and ||x(•) -x(•)|| ∞ ≤ ε we have ϕ(x(1)) ≥ ϕ( x(1)) -Mν.
If the Mayer problem (1) is calm at x(•), then by Proposition 6.4.3 from [START_REF] Clarke | Optimization and nonsmooth analysis[END_REF], there exists k > 0 such that ( x(•), ū(•)) is a strong local minimizer of the following minimization problem:

min ϕ(x(1)) + kd K (x(1)) : x(•) ∈ S [0,1] (x 0 ) . (11) Since ϕ(•) + kd K (•) is no longer C 2 ,
we need the following (non smooth) maximum principle that can be deduced, for instance, from [START_REF] Vinter | Optimal control[END_REF]. Theorem 6.4 Consider a strong local minimizer ( x(•), ū(•)) of ( 11) with x(1) ∈ ∂ K. Then there exists an absolutely continuous function p : [0, 1] → R n satisfying (4), [START_REF] Krener | The high order maximal principle and its application to singular extremals[END_REF] and the transversality condition

-p(1) ∈ ∇ϕ( x(1)) + k[0, ∇b K ( x(1))], (12)
where [0, ∇b K ( x( 1))] = {µ∇b K ( x(1)) : µ ∈ [0, 1]}.

Denote by q 1 (•) the solution of (4) for p(1) = -∇ϕ( x(1)) and by q 2 (•) the solution of (4) for p(1) = -∇ϕ( x(1)) -k∇b K ( x(1). Note that q i (1) are extremal points in the transversality condition [START_REF] Hoehener | Variational approach to second-order optimality conditions for control problems with pure state constraints[END_REF]. The mappings q i (•) do not satisfy, in general, the maximum principle [START_REF] Krener | The high order maximal principle and its application to singular extremals[END_REF] but are used below to express the second order necessary optimality conditions. Define for i = 1, 2

ϕ 1 (•) := ϕ(•), ϕ 2 (•) := ϕ(•) + kb K (•), H i xx [•] := H xx (•, x(•), q i (•), ū(•)). The mappings H i xu [•]
and H i uu [•] are defined similarly. Theorem 6.5 Consider a strong local minimizer ( x(•), ū(•)) of ( 11) with x(1) ∈ ∂ K and p(•) as in Theorem 6.4. Then for every u(•) ∈ U A such that H u [•]u(•) = 0 and the corresponding solution y(•) of ( 9) verifies ∇b K ( x(1)), y(1) = 0, and for any v(

•) ∈ L ∞ ([0, 1]; R m ) satisfying v(t) ∈ T (2)
U ( ū(t), u(t)) for a.e. t ∈ A, there exists i ∈ {1, 2} such that

0 ≥ 1 0 H i u [t]v(t) + 1 2 y(t) * H i xx [t]y(t) + u(t) * H i xu [t]y(t) + 1 2 u(t) * H i uu [t]u(t))dt - 1 2 y(1) * (ϕ i ) ( x(1))y(1). ( 13
)
Consider the solution q(•) of the linear system (4) such that q(1) = -∇b K ( x(1)) and define the subspaces T t ⊂ P ū(t)

T t := {u ∈ P ū(t) : q(t), f u [t]u = 0}.
By Lemma 6.2, if the abnormal maximum principle of Theorem 6.1 holds true, then T t = P ū(t) a.e. in [0, 1] and q(•) satisfies the maximum principle. Observe next that for every measurable essentially bounded selection u(t) ∈ T t and the corresponding solution y(•) of ( 9) we have ∇b K ( x(1)), y(1) = 0. Therefore, by Theorem 6.5, (13) holds true for every u(•) ∈ U A such that u(t) ∈ T t a.e. in [0, 1]. This and similar, but lengthly proofs as those of results from Sections 3,4,5 allow to deduce from (13) the following three theorems.

Theorem 6.6 Consider a strong local minimizer

( x(•), ū(•)) of (1) such that x(1) ∈ ∂ K. Assume that (1) is calm at x(•), that for some 0 ≤ t 1 < t 2 ≤ 1, ū(•) is integrable on [t 1 ,t 2 ] and for j = 1, 2 H j uu [t] = 0 a.e. in [t 1 ,t 2 ]; H j xu [•], f u [•] are Lipschitz on [t 1 ,t 2 ]; H j xx [•], f x [•] are essentially bounded on [t 1 ,t 2 ]. (14)
Then for a.e. t ∈ [t 1 ,t 2 ] with ū(t) ∈ IntU there exists

i ∈ {1, 2} such that H i xu [t] f u [t] is symmetric on T t and R i (t) := 1 2 f u [t] * H i xx [t] f u [t] -1 2 d dt H i xu [t] f u [t] -H i xu [t] f x [t] f u [t] + 1 2 H i xu [t] d dt f u [t] is negative definite on T t .
Theorem 6.7 Let U be a convex polytope and ( x(•), ū(•)) be a strong local minimizer of (1) such that x(1) ∈ ∂ K. Assume that (1) is calm at x(•) and for some 0 Corollary 6.8 Under all the assumptions of Theorem 6.7 suppose that f is as in [START_REF] Agrachev | Control theory from the geometric viewpoint[END_REF], where g(•) ∈ C 2 (R n ; M(n × m)) and U is a closed convex polytope in R m and let g k (x) denote the k-th column of g(x). For i = 1, 2 define the matrix valued mappings Θ Remark 6.10 Similar results can be obtained when K is an intersection of a finite family of closed sets K j of class C 2 . Then the assumptions and statements will involve Hamiltonians H j for each active index j.

≤ t 1 < t 2 ≤ 1, ū(•) is integrable on [t 1 ,
i : [0, 1] → M(m × m), where θ i k,l (t) := q i (t), [g k , g l ]( x(t)) ∀ k, l ∈ {1, . . . ,

Conclusions and future works

For a nonlinear Bolza optimal control problem we proposed pointwise second order necessary optimality conditions extending those of Goh to the case of a closed control constraint given by a polytope or by a C 2 subset of R m and for integrable optimal controls that may take values on the boundary of control constraints. Furthermore, we have shown that on the interior of the control constraints Goh conditions hold true for integrable optimal controls without imposing any structural properties on the control set. This was done by investigating first the free endpoint problem and then, the calmness assumption allowed us to reduce the Mayer problem with endpoint constraints to the one without them, but involving a nondifferentiable cost function. For this equivalent problem we obtained similar pointwise conditions using the developed variational approach for the free endpoint problems. Our future work will concern a larger class of end-point constraints and more general control sets.

Theorem 2 . 4 (

 24 Integral second order conditions) Consider a strong local minimizer ( x(•), ū(•)) of the Mayer problem (3).Then for all u

Theorem 4 . 1

 41 Let U ⊂ R m be a convex polytope in R m . Consider a strong local minimizer ( x(•), ū(•)) of the Mayer problem (3) and p(•) as in the maximum principle of Theorem 2.

Theorem 5 . 2

 52 Let U be of class C 2 . Consider a strong local minimizer ( x(•), ū(•)) of the Mayer problem (3) and p(•) as in the maximum principle of Theorem 2.3.

  3. Assume that ū(•) is integrable on [t 1 ,t 2 ] for some 0 ≤ t 1 < t 2 ≤ 1 and (10) holds true. Then H xu [t] f u [t] is symmetric on the subspace P ū(t) and R(t) is negative definite on P ū(t) for a.e. t ∈ [t 1 ,t 2 ].

  m}. Then for a.e. t ∈ [t 1 ,t 2 ], the subspace Θ(t) P ū(t) is orthogonal to P ū(t) . In particular, Θ(t) is equal to zero for a.e. t ∈ [t 1 ,t 2 ] such that ū(t) ∈ Int U. Proof. By Theorem 4.1 for a.e. t ∈ [t 1 ,t 2 ], H xu [t] f u [t] is symmetric on P ū(t) . Fix such t. Since for all a, b ∈ P ū(t) , a

* (H xu

  t 2 ] and (14) holds true for j = 1, 2. Then for a.e. t ∈ [t 1 ,t 2 ] there exists i ∈ {1, 2} such that H i xu [t] f u [t] is symmetric on T t and R i (t) is negative definite on T t .

  m}.Then for a.e. t ∈ [t 1 ,t 2 ] there exists i ∈ {1, 2} such thatΘ i (t) (T t ) is orthogonal to T t . Let U be of class C 2 and ( x(•), ū(•)) be a strong local minimizer of the Mayer problem (1) such that x(1) ∈ ∂ K. Assume that (1) is calm at x(•) and for some 0 ≤ t 1 < t 2 ≤ 1, ū(•) is integrable on [t 1 ,t 2 ].Then for a.e. t ∈ (t 1 ,t 2 ) such that ū(t) ∈ ∂U and for some δ t > 0 and j = 1, 2 the following relations are satisfied:u * H j uu [s] -H j u [s] b U ( ū(s)) u = 0 ∀ u ∈T s and a.e. s ∈ (tδ t ,t + δ t ) with || ū(t)ū(s)|| ≤ δ t ; H j xu [•], f u [•] are Lipschitz on (tδ t ,t + δ t ); H j xx [•], f x [•] are essentially bounded on (tδ t ,t + δ t ) we can find i ∈ {1, 2} for which H i xu [t] f u [t] is symmetric on T t .

	Consider subspaces S(s) defined as in Section 5.
	Theorem 6.9

* This work was co-funded by the European Union under the 7th Framework Programme "FP7-PEOPLE-2010-ITN", grant agreement number 264735-SADCO.


Keywords: Antenna characterization, Cauchy's method, complex natural resonance, matrix pencil method, poles, residues, singularity expansion method

In this paper, the accuracy of the Singularity Expansion Method (SEM) used for antenna characterization is investigated. A well-known limitation of the SEM is that pole extraction is very sensitive to noise. A comparison between two main methods of pole extraction is presented. The Matrix Pencil (MP) method and the Cauchy's method are used to extract poles from the radiated fields of a dipole antenna and two bowtie antennas. Results are presented for simulated fields and the robustness to a white Gaussian noise is also analyzed. We show that the MP method allows working with lower SNR than Cauchy's method and is more accurate for field reconstruction.

I. INTRODUCTION

OR many years, the Singularity Expansion Method (SEM) [START_REF] Baum | On the singularity expansion method for the solution of electromagnetic interaction problems[END_REF], introduced by C. E. Baum in 1971, has been widely used in the radar domain. The SEM represents a solution of an electromagnetic problem in terms of singularities (poles) in the complex frequency plane. Since singularities are independent of the direction of the incoming wave, the SEM has been widely studied for target identification [START_REF] Baum | The singularity expansion method and its application to target identification[END_REF] [START_REF] Chauveau | Selection of contributing natural poles for the characterization of perfectly conducting targets in resonance region[END_REF]. The information on poles can give some indications on the general shape and the constitution of the illuminated target. Moreover, the SEM has been used in the antenna domain such as in [START_REF] Tesche | On the analysis of scattering and antenna problems using the singularity expansion technique[END_REF] where the SEM formalism has been applied for modeling the time response of the current of a thin wire antenna. In [START_REF] Barnes | On the singularity expansion method as applied to the EMP analysis and simulation of the cylindrical dipole antenna[END_REF], Barnes analyses the dipole antenna response to an electromagnetic pulse using the SEM. More recently, this method has been applied in both time and frequency domains to model antenna Manuscript received December 10, 2012. This work was financially supported in part by the Direction Générale de l"Armement (DGA).

François Sarrazin and Ala Sharaiha are with the Institute of Electronics and Telecommunications of Rennes, University of Rennes 1, Rennes, France (corresponding author to provide phone: +33-223-233-374; e-mail: francois.sarrazin@univ-rennes1.fr, ala.sharaiha@univ-rennes1.fr).

Philippe Pouliguen is with the Direction Générale de l"Armement (DGA), MRIS, Paris, France, (e-mail: philippe.pouliguen@dga.defense.gouv.fr).

Patrick Potier and Janic Chauveau are with the Direction Générale de l"Armement (DGA), MI, Bruz, France, (e-mail: patrick.potier@dga.defense.gouv.fr, janic.chauveau@dga.defense.gouv.fr).

effective length, in order to fully describe the antenna radiation patterns, directivity and gain using only a few sets of parameters (poles and residues) [START_REF] Licul | Unified frequency and time-domain antenna modeling and characterization[END_REF] [START_REF] Marchais | Compact characterisation of ultrawideband antenna responses from frequency measurements[END_REF]. In [START_REF] Licul | Unified frequency and time-domain antenna modeling and characterization[END_REF] [START_REF] Marchais | Compact characterisation of ultrawideband antenna responses from frequency measurements[END_REF], the authors extract poles from the whole antenna response. But, according to the theory of C.E. Baum [START_REF] Baum | On the singularity expansion method for the solution of electromagnetic interaction problems[END_REF], physics phenomena can be observed only in the late time antenna response. In this paper, we focus on the physical approach of the SEM and only the late time response is considered.

The well-known limitation of the SEM is that pole extraction is very sensitive to noise. In transient domain, there are two main methods to extract poles from the antenna electromagnetic response. The first one is the well-known Prony"s method introduced by Baron de Prony in 1795 [START_REF] Prony | Essai expérimental et analytique sur les lois de la dilatabilité de fluides élastiques et sur celles de la force expansive de la vapeur d"alcool[END_REF]. This method has been modified to be used on noisy data using a Least Square (LS) approach in 1950 [START_REF] Householder | On Prony's method of fitting exponential decay curves and multiple-hit survival curves[END_REF] and a Total Least Square (TLS) approach in 1987 [START_REF] Rahman | Total least square approach for frequency estimation using linear prediction[END_REF]. In 1990, Hua and Sarkar suggested the Matrix Pencil (MP) method [START_REF] Hua | Matrix pencil method for estimating parameters of exponentially damped/undamped sinusoids in noise[END_REF] also based on a TLS approach. They have compared Prony and MP methods and have shown that the MP method is more robust to noise than Prony"s algorithm [START_REF] Sarkar | Using the Matrix Pencil method to estimate the parameters of a sum of complex exponentials[END_REF]. This has been verified on noisy antenna responses [START_REF] Sarrazin | Comparison between matrix pencil and Prony methods applied on noisy antenna responses[END_REF]. Moreover, the MP method is computationally more efficient [START_REF] Hua | Generalized pencil-of-function method for extracting poles of an EM system from its transient response[END_REF]. More details on MP are given in appendix I. In the frequency domain, the main way to extract poles of a transfer function is the Cauchy"s method developed in 1821 by Cauchy [START_REF] Cauchy | Sur la formule de Lagrange relative à l"interpolation[END_REF]. A TLS approach has also been used to improve its robustness to noise [START_REF] Adve | The effect of noise in the data on the Cauchy method[END_REF]. More details on Cauchy"s method are given in appendix II. MP and Cauchy methods are the two main efficient methods of pole extraction, in transient and frequency domains, respectively, but there are few works dealing with their comparison [START_REF] Sarrazin | Analysis of two methods of poles extraction for antenna characterization[END_REF].

The objective of this paper is to determine which method is the most appropriate to extract poles for antenna characterization by using either the Total Least Square Matrix Pencil (TLS MP) method in the time domain or the Total Least Square Cauchy (TLS Cauchy) method in the frequency domain. In section II, the SEM is presented. Next, in section III, the two methods are applied to noiseless fields radiated by three different antennas: a narrow band dipole antenna and two Ultra Wide Band (UWB) bowtie antennas with different flare angles. To study the robustness of these methods in the presence of noise, two different kinds of noise are added to the simulated fields and results are compared in section IV. The SEM was developed to describe the global behavior of an object"s response excited by an electromagnetic wave. In the time domain, this response is composed of two successive parts. The first one is called the early time response and is mainly due to the excitation impulse. The duration E T of this early time response depends on the pulse duration P T and the greatest dimension D of the antenna as

Accuracy of Singularity Expansion Method in

P E T c D T   /
where c is the speed of light. The second part, called the late time response, occurs after the early time response and is only due to the radiation of the induced current propagated on the antenna after its illumination. The SEM allows modeling the late time response of an object as a decaying exponential sum as:

, ) (

    M n n n s s R s H (2)
where s is the Laplace complex variable. Each pole is defined as:

n n n j s     ( 3
)
where n  is the negative damping coefficient of the n th pole and n  is the resonant pulsation of the n th pole.

III. POLE EXTRACTION FROM NOISELESS RADIATED ANTENNA RESPONSES

This section presents a comparison between the two methods applied to fields radiated by three antennas: a dipole and two bowtie antennas shown in Fig. 1. Their lengths are L = 33.75 mm. The diameter of the dipole is D = 1.12 mm, so the ratio L/D = 30. Its gap length is 1.12 mm and the impedance of its lumped port is 73 Ω. Bowtie antennas have two different flare angles: 20° and 45° and the impedance of their lumped port is 200 Ω. These antennas are simulated using CST Microwave Studio [START_REF]The CST website[END_REF] with transient solver. A power source is used and the excitation signal is a Gaussian pulse. A far field probe is used to measure the electric field in the boresight direction at a distance R = 2 m.

A. Dipole Antenna

The field radiated by the dipole antenna is shown in Fig. algorithm is used on the frequency radiated field obtained by a Fast Fourier Transform (FFT). Results of both extractions are shown in Fig. 3 and Fig. 4 for poles and residues, respectively. Three pairs of poles and residues are extracted by the two methods and their agreement is very good (relative difference less than 1 %). Resonant frequencies of the poles are represented using black lines on its input impedance in Fig. 5. We can see that these poles correspond almost to resonances of the dipole antenna at 2 /  , 2 / 3 and 2 / 5 , where  is the wavelength in free space. So all extracted poles have physical meanings. Late time radiated fields are reconstructed using these two different sets of poles and residues, extracted with MP and Cauchy methods, and (1). Results are presented in Fig. 6. The three curves are close; the Normalized Mean Square Error (NMSE) is less than 3 % for both methods. In Fig. 6, only the late time response of the dipole antenna is presented whereas in Fig. 2 the entire radiated field is shown. The early time duration is 0.2 ns. In Fig. 3, the marker size depends on the weight of each pole. The weight is computed as the ratio between residue and damping coefficient of the pole and is normalized by the maximum weight. This representation is a good way to estimate the contribution of each pole [START_REF] Chauveau | Selection of contributing natural poles for the characterization of perfectly conducting targets in resonance region[END_REF]. In Table I, NMSE of the dipole responses reconstructed using only some of the physical poles extracted are presented. It shows that the pole 2 is the most important contribution to the late time dipole response, then the pole 1 and finally the pole 3. MP and Cauchy methods allow extracting the same physical poles from the noiseless field radiated by a dipole antenna. Moreover, it is possible to reconstruct the radiated field with a very small error. It means that a 6-poles set is enough to accurately model the late time response of a dipole antenna.

B. Bowtie Antenna with 20° flare angle° Fig. 7 shows the field radiated by the bowtie antenna with flare angle equal to 20°. This antenna has a wider band than the dipole antenna. In fact, larger flare angles correspond to wider bandwidths. Since the antenna is less resonant, it is more difficult to define which poles have physical meanings when extracted by the two methods. In order to select poles extracted in the time domain, we use the Window Moving Technique (WMT) also known as Time-Frequency analysis [START_REF] Van Blaricum | An analysis of existing Prony"s method techniques[END_REF] [START_REF] Manteghi | Short time Matrix Pencil for chipless RFID detection applications[END_REF]. The idea of the WMT is to move a time window with a given duration through the entire signal by small time steps. The minimum time step is equal to the sample period. MP method is applied on each time window. The assumption is that, depending on the window, the position of the mathematical poles will change from window to window, whereas the physical poles will remain essentially unchanged.

To improve pole extraction in the frequency domain, the frequency radiated field is split into two responses (from 1 to 7 GHz and from 7 to 23 GHz). In this paper, subbands have been selected empirically after several tests. However, a systematic approach relating the resonant behavior of the antenna to the needed number and size of the subbands is under investigation. The Cauchy"s algorithm is applied separately on these two subbands. Each subband contains only a few poles so it is easier to extract them accurately. Residues are then computed again using the extracted poles in each subband and the entire response in order to accurately model the complete response. Poles extracted with the two methods are shown in Fig. 8. Three pairs of poles are extracted by the two methods and they are in very good agreement. A maximum relative error of 1 % is found for the resonant frequencies of these poles whereas a relative error varying from 1 to 3 % is obtained for the damping coefficients.

Resonant frequencies of the poles extracted with the MP algorithm are represented with a black line on the antenna"s input impedance in Fig. 9. We can see that these frequencies mainly correspond to the resonances of the antenna input impedance, i.e. when its imaginary part is close to zero. Associated residues are presented in Fig. 10. Late time radiated field is reconstructed using poles and residues extracted with MP and Cauchy methods and results are shown in Fig. 11. The three curves are overlapped. NMSE with the two reconstructed fields are less than 4 %. Therefore, it is possible to model the late time response of this wideband antenna using only a set of 6 poles and residues. NMSE of the reconstructed fields using only some poles are presented in the Table II. It shows that the pole 2 is still the dominant one followed by the pole 3. Note that even if the pole 1 has the strongest damping coefficient, it provides the smallest contribution to the radiated field.

C. Bowtie Antenna with 45° flare angle

The field radiated by the bowtie antenna with flare angle equal to 45° is shown in Fig. 12. This UWB bowtie antenna has a much wider band than the previous one. Therefore, we need to use the WMT in the transient domain and split the frequency response into three subbands in order to extract the physical poles presented in Fig. 13. The number of extracted poles is unchanged but damping coefficients are much larger in absolute value than for the previous bowtie antenna, especially for the third pole (around -35.10 9 Neper/s). The resonant frequencies, extracted with the MP method, are represented with a black line on the impedance of the antenna in Fig. 14. Associated residues are presented in Fig. 15 and they are in very good agreement. Late time responses reconstructed using the whole poles are presented in Fig. 16. The NMSE is less than 3% for both methods. NMSE of the reconstructed fields using only a few poles are presented in Table III. They confirm the weight of each pole presented in Fig. 13, i.e. the pole 2 is the dominant one, followed by the poles 3 and 1.

IV. POLE EXTRACTION FROM NOISY RADIATED ANTENNA RESPONSES

Behaviors of MP and Cauchy methods in presence of noise are now compared. We consider two different approaches as presented in the flowcharts in Fig. 17. The first one considers that the original data is in the transient domain and a noise is added to the transient simulated field to obtain a noisy response with a desired Signal to Noise Ratio (SNR). Then, poles are extracted directly with MP and using a FFT. In the second one, it is considered that the original data is in the frequency domain, so the noise is added to the frequency simulated field. The Cauchy"s method is then directly applied and an IFFT is performed to use MP. Two different kinds of noise are considered: a White Gaussian Noise (WGN) and a mixed noise composed of WGN, impulse noise and single carrier noise.

Due to the noise and the overestimation of the number of poles to be extracted, some poles can be very different of those extracted from the noiseless case. These poles are called "mathematical" because they do not have a physical sense.

Since the aim of this study is to define which method is the best to extract physical poles of an antenna, we only kept poles close to those extracted from the noiseless field. We consider a pole well extracted when its resonant frequency is around 5 % of the original one, i.e. poles extracted from noiseless data, and damping coefficient around 30 %. Using these poles extracted from noisy data, the radiated field is reconstructed and the NMSE compared to the simulated field is computed. These operations (noise addition, pole extraction, field reconstruction and NMSE computation) are repeated for each SNR value from -10 to 60 dB with a 5 dB step. For each SNR value, these operations are repeated 100 times to limit the random effect of the noise. An average NMSE is then computed using the 100 NMSE values. Results for the WGN case for the first and second approach are given in Fig. 18 and Fig. 19, respectively. The two approaches provide very similar results. In fact, adding the WGN in the frequency or the transient domain does not change the results. Therefore, the FFT and the IFFT seem to not disturb pole extraction with the two methods. We can notice that, for a given SNR, the NMSE is generally more important for the 45° bowtie than the 20° bowtie and the dipole antennas. Therefore, we may conclude that for wideband antennas, it is difficult to extract poles with a good accuracy whatever the method used. For the dipole antenna case, the difference between the two methods is around 5 dB to obtain the same NMSE. However, for the two others antennas, the difference is much higher. In fact, for the 20° bowtie, one needs a 45 dB SNR to obtain a 10 % NMSE with the Cauchy"s method, unlike 20 dB sufficient with the MP method. For the 45° bowtie, there is a 20 dB difference. So the MP method allows dealing with signals with SNR 5 to 25 dB lower than using Cauchy"s method.

The same analysis is done using the mixed noise. Since the two approaches give very similar results, only those obtained from approach 1 are presented in Fig. 20. As for the previous analysis, it is easier for both methods to extract poles from the dipole antenna response than for the bowtie antennas. Otherwise, even if the results are close between MP and Cauchy methods for the dipole response, the MP algorithm allows obtaining a lower NMSE than the Cauchy"s method for the two bowties. As an example, the dipole response, noiseless and in presence of mixed noise for SNR = 10 dB, and the reconstructed field are presented in Fig. 21. The NMSE of the reconstructed field is 10 %.

From these analyses, we can conclude that the MP method is less sensitive to noise than the Cauchy"s one, especially when applied to wideband or UWB antenna responses. Nevertheless, the SNR needed to obtain poles with accuracy is quite high. During measurements, it will be necessary to increase the SNR as much as possible using especially a time gating to filter unwanted echoes [START_REF] Marchais | Compact characterisation of ultrawideband antenna responses from frequency measurements[END_REF]. V. CONCLUSION

The well-known limitation of the SEM is that pole extraction is high sensitive to noise. Therefore, one has to extract poles very carefully in order to obtain physical ones. In this article, two of the best numerical methods of pole extraction, MP and Cauchy, have been applied to fields radiated by three different antennas, a dipole and two bowtie antennas with different flare angles. In the noiseless case, both methods allow extracting physical poles with a good accuracy. It follows that the late time response of these antennas can be reconstructed by using only a 6-poles set. We also compared the robustness of MP and Cauchy methods in the presence of two different kinds of noises. For the simple case of the dipole antenna, results are close but the MP method is more accurate than Cauchy"s method when SNR becomes low. For the two bowtie antennas, the difference between the two methods is more significant. Indeed, the MP method allows dealing with signal with a SNR 20-25 dB lower than that needed for Cauchy"s method.

APPENDIX I MATRIX PENCIL METHOD

The data samples k y are defined as:

                   1 1 1 2 1 1 0 K L K L K L L y y y y y y y y y Y        (5)
L is called the Pencil parameter and is very important to filter noise. It is usually chosen between K/3 and K/2 [START_REF] Hua | Matrix pencil method for estimating parameters of exponentially damped/undamped sinusoids in noise[END_REF]. In fact, the variance of the extracted poles is the lowest for these values.    

      H V U Y ' ' ' '   where:     , ' * ) ( 2 1 M L K M u u u U    (6)   , 0 0 0 0 0 0 0 0 0 0 0 0 ' * 2 1 M M M                   (7)
. ' * ) 1 ( 2 1 M L M v v v V    (8) From   ' V , it is possible to define two submatrices   ' 1 V and   ' 2 V as:   , ' ' * ) 1 ( 1 1 M L V l V         (9)
and

  , ' ' * ) 1 ( 1 2 M L L l V V          (10)
where i l is the i th line of   ' V . Using the TLS approach [21], poles are obtained from the eigenvalues of       H H V V ' '

APPENDIX II CAUCHY"S METHOD

The main idea of this method is to approximate the transfer function ) (s H of an antenna into a ratio of two polynomials P and Q as:

A QR decomposition of ( 13) is made such as:

0 0 22 12 11              b a R R R . ( 14
)
So these two following equations are obtained: [START_REF] Cauchy | Sur la formule de Lagrange relative à l"interpolation[END_REF] . 0 12 11   b R a R [START_REF] Adve | The effect of noise in the data on the Cauchy method[END_REF] A SVD of 22 R is done such as:

0 22  b R
      , 22 H V U R   (17)
Using a TLS approach [START_REF] Van Huffel | Analysis of the total least square problem and its use in parameter estimation[END_REF], coefficients b are obtained from

  1   Q V b
. Coefficients a are then computed using [START_REF] Adve | The effect of noise in the data on the Cauchy method[END_REF].

Poles are now found by computing squares of the Q polynomial.

      Q n n Q n Q n n s s b s b 1 0
) ( [START_REF]The CST website[END_REF] where n s are the poles and s the Laplace variable. The transfer function can now be written as: [START_REF] Manteghi | Short time Matrix Pencil for chipless RFID detection applications[END_REF]

  2. The MP algorithm is directly applied on the late time transient radiated field of the dipole antenna whereas the Cauchy"s a. Dipole b. 20° Bowtie c. 45° Bowtie Fig. 1. Geometry of the three considered antennas.

Fig. 2 .

 2 Fig. 2. Simulated transient electric far field of the dipole antenna in the boresight direction.

Fig. 3 .

 3 Fig. 3. Poles extracted from transient simulated electric far field of the dipole antenna with MP and Cauchy methods in a complex plane.

Fig. 4 .

 4 Fig. 4. Residues extracted from transient simulated electric far field of the dipole antenna with MP and Cauchy methods in a complex plane.

Fig. 7 .

 7 Fig. 7. Simulated transient electric far field of the bowtie antenna in the boresight direction.

Fig. 6 .

 6 Fig. 6. Radiated fields of the dipole antenna, simulated and reconstructed using Cauchy poles and Matrix Pencil poles.

Fig. 8 .

 8 Fig. 8. Poles extracted from simulated transient electric far field of the bowtie 20° antenna with MP and Cauchy methods in a complex plane.

Fig. 5 .

 5 Fig. 5. Complex impedance of the dipole antenna. Black lines represent resonant frequencies of extracted poles.

Fig. 11 .

 11 Fig. 11. Radiated fields of the bowtie antenna, simulated and reconstructed using Cauchy poles and Matrix Pencil poles.

Fig. 10 .

 10 Fig. 10. Residues extracted from simulated transient electric far field of the bowtie antenna with MP and Cauchy methods in a complex plane.

Fig. 9 .

 9 Fig. 9. Complex impedance of the bowtie antenna. Black lines represent resonant frequencies of extracted poles.

Fig. 12 .

 12 Fig. 12. Simulated transient electric far field of the 45° bowtie antenna in the boresight direction.

Fig. 15 .

 15 Fig. 15. Residues extracted from simulated transient electric far field of the bowtie antenna with MP and Cauchy in a complex plane.

Fig. 16 .Fig. 17

 1617 Fig. 16. Radiated fields of the bowtie antenna, simulated and reconstructed using Cauchy poles and Matrix Pencil poles.

  Fig. 17. The flowcharts of the two approaches (1: transient, 2: frequency).

Fig. 13 .Fig. 14 .

 1314 Fig. 13. Poles extracted from simulated transient electric far field of the bowtie 45° antenna with MP and Cauchy methods in a complex plane.

  a. Matrix Pencil method b. Cauchy method Fig. 19. NMSE of the reconstructed fields versus SNR for approach 2 and for white Gaussian noise. a. Matrix Pencil method b. Cauchy method Fig. 18. NMSE of the reconstructed fields versus SNR for approach 1 and for white Gaussian noise.

T

  is the sampling period. A data matrix   Y is constructed from the data samples.



  c. Matrix Pencil method d. Cauchy method Fig. 20. NMSE of the reconstructed fields versus SNR for approach 1 and for mixed noise.

Fig. 21 .

 21 Fig. 21. Field radiated by the dipole antenna, noiseless, in presence of mixed noise for an SNR of 10 dB and reconstructed with NMSE of 10 %.

1 .

 1 It is now possible to compute residues from (4).

  Time and Frequency Domains to Characterize Antennas in Presence of Noise

	II. SEM THEORY
	François Sarrazin, Janic Chauveau, Philippe Pouliguen, Patrick Potier and Ala Sharaiha, Senior
	Member, IEEE
	F

TABLE I NMSE

 I OF THE RECONSTRUCTED DIPOLE RESPONSE

	Poles	1-2-3	1-2	1-3	2-3	2
	NMSE (%)	2.7	6.2	92	19.5	24.5

TABLE II

 II

	NMSE OF THE RECONSTRUCTED BOWTIE 20° RESPONSE
	Poles	1-2-3	1-2	1-3	2-3	2
	NMSE (%)	4	26	80	12	33
			TABLE III
	NMSE OF THE RECONSTRUCTED BOWTIE 45° RESPONSE
	Poles	1-2-3	1-2	1-3	2-3	2
	NMSE (%)	2.5	54	98	4.5

  Then, a Singular Value Decomposition (SVD) is applied to this matrix as       H Y contains exactly M nonzero eigenvalues corresponding to the M poles of the system. However, in the noisy case, the other eigenvalues are not exactly equal to zero. So it is necessary to filter these eigenvalues. The smallest ones, minor to a threshold ε, are set to zero. A new matrix  

	Y		U		V	where H defines
	the Hermitian transpose,   U and   V are unitary matrices,
	composed of the eigenvectors of    H Y Y	and     Y Y H	,
	respectively, and    is a diagonal matrix containing the
	singular values of   Y . In the noiseless case, the matrix  
						' Y can be written
	as

François Sarrazin was born in Poitiers, France, in 1987. He received the M.S. degree in electronics and electrical engineering from Polytech"Nantes (Ecole polytechnique de l"université de Nantes), Nantes, France, in 2010. Since 2010, he has been working toward the Ph.D. degree at the Institute of Electronics and Telecommunications of Rennes (IETR), University of Rennes 1, Rennes, France. His current research activities are in the field of antenna characterization using the Singularity Expansion Method (SEM) applied on Radar Cross Section (RCS) measurement. Janic Chauveau was born in Saint-Nazaire, France, in 1981. He received the M.S. degree in electronics and electrical engineering from Polytech"Nantes (Ecole polytechnique de l"université de Nantes), Nantes, France, in 2004. He received the Ph.D. degree from the University of Nantes, France, in 2007. He is currently working in the field of radar signatures at the Direction Générale de l"Armement (DGA), DGA-Maîtrise de l"Information, Bruz, France.

His research interests include electromagnetic scattering and inverse scattering problems.

Philippe Pouliguen received the M.S. degree in signal processing and telecommunications, the Doctoral degree in electronic and the "Habilitation à Diriger des Recherches" degree from the University of Rennes 1, France, in 1986, 1990 and 2000. In 1990, he joined the Direction Générale de l"Armement (DGA) at the Centre d"Electronique de l"Armement (CELAR), now DGA Information Superiority (DGA/IS), in Bruz, France, where he was a "DGA expert" in electromagnetic radiation and radar signatures analysis. He was also in charge of the EMC (Expertise and ElectoMagnetism Computation) laboratory of CELAR. Since December 2009, he is the head of "acoustic and radio-electric waves" scientific domain at the office for advanced research and innovation of the strategy directorate, DGA, Bagneux, France. His research interests include electromagnetic scattering and diffraction, Radar Cross Section (RCS) measurement and modeling, asymptotic high frequency methods, radar signal processing and analysis, antenna scattering problems and Electronic Band Gap Materials. In these fields, he has published more than 30 articles in refereed journals and more than 100 conference papers.

Patrick Potier received the Ph.D. degree in structure and property of the material from the University of Rennes, Rennes, France, in 1984. From September 1984 to September 1987, he was a Research Engineer with Thomson CSF, Paris, France. Since September 1987, he has been an Engineer with the Center "Maîtrise de l"information" (Information superiority) de la "Direction Générale de l"armement" (General Armaments Directorate, the French procurement agency), Bruz, France. He treats the aspects of antennas and radiation and ensures the follow-up of various studies and theses mainly on the topics broad band reflect or transmit arrays. Currently, he is a Professor at the University of Rennes 1, Rennes, France and the Co-Head of the "Complex Radiation Systems" Team in the "Antennas and Microwave Devices" Department. His main research activities include small antennas, broadband and UWB antennas, miniaturization, printed spiral and helical antennas, antennas for mobile communications, etc. Conducted and involved in more than 20 development projects for private companies and participates in the European Network of Excellence ACE (Antenna Center of Excellence) in the small antenna WP. He is the author and coauthor of 48 international papers, 150 conference presentations and holds twelve European patents.

Ala

Prof. Sharaiha is a senior member of the IEEE and is a reviewer for the IEEE APS, IEEE AWPL, the IET Letters and the IET Microwave Antennas Propagation. He was the conference Chairman of the 11th International Canadian Conference ANTEM (Antenna Technology and Applied Electro-Magnetics), held at Saint-Malo in France, 2005.