Fakhteh Saadatniaki

Anton Korniienko

Gerard Scorletti

Usman A Khan

On the input-output approach towards distributed estimation

Keywords: Distributed control, Consensus control and estimation, Agents and autonomous systems

In this paper, we study distributed estimation of continuous-time, linear time-invariant systems monitored by a network of agents communicating over a graph. We assume that no agent may possess enough measurements in its neighborhood to estimate the entire state vector on its own. In this context, we provide a networked Kalman-type estimator that combines prediction and innovation with information fusion among the neighboring agents and consider an approach based on designing static estimator gains. The main contribution of this paper is to analyze the estimation error using the notions of dissipativity and the input-output approach, which enable us to formulate stability and performance arguments as quasiconvex optimization problems involving linear matrix inequalities. We show that the resulting estimation error is stable and further ensures a given level of performance regarding noise rejection. Simulations illustrate the concepts described in this paper.

I. INTRODUCTION

Recently, there has been a significant amount of work in the area of distributed estimation of physical and social phenomena [START_REF] Khan | On the stability and optimality of distributed Kalman filters with finite-time data fusion[END_REF]- [START_REF] Doostmohammadian | Graph-theoretic inference in social networks[END_REF]. Some prime examples include tracking in multi-agent networks [START_REF] Oh | Tracking and coordination of multiple agents using sensor networks: System design, algorithms and experiments[END_REF], [START_REF] Safavi | An opportunistic algorithm for localization in mobile networks[END_REF], estimation or optimization in large-scale systems [START_REF] Khan | Cooperation for aggregating complex electric power networks to ensure system observability[END_REF], [START_REF] Xi | ADD-OPT: Accelerated distributed directed optimization[END_REF], or studying opinion formation or voting models in complex social networks [START_REF] Rahimian | Group decision making and social learning[END_REF]. Estimating the underlying dynamics in such systems is a challenging problem as measurements are typically distributed over a network of geographically-dispersed agents (sensors, robots, individuals). Distributed estimation, see e.g., relevant work in [START_REF] Casbeer | Distributed information filtering using consensus filters[END_REF]- [START_REF] Battistelli | Kullback-Leibler average, consensus on probability densities, and distributed state estimation with guaranteed stability[END_REF], thus enables estimating the dynamics without collecting measurements at a central location, a practically infeasible task, but by enabling local estimation and interagent information exchange. In this paper, we focus on the distributed estimation of a Continuous-Time, Linear Time-Invariant (CT-LTI) dynamical system whose measurements are distributed over a network of sparsely-connected agents.

Relevant work on decentralized estimation of CT-LTI systems can be found in [START_REF] Stankovic | Decentralized parameter estimation by consensus based stochastic approximation[END_REF]- [START_REF] Zhou | Distributed estimation for moving target based on state-consensus strategy[END_REF]. For instance, Ref. [START_REF] Stankovic | Decentralized parameter estimation by consensus based stochastic approximation[END_REF] considers decentralized state estimation of linear stochastic systems based on a combination of local Kalman filters and a dynamic consensus scheme among the agents. Ref. [START_REF] Olfati-Saber | Coupled distributed estimation and control for mobile sensor networks[END_REF] proposes a theoretical framework for coupled distributed estimation and motion control of mobile sensor networks based on distributed Kalman filtering for collaborative target † FS and UAK are with the Department of Electrical and Computer Engineering at Tufts University, Medford, MA 02155, USA, {fakhteh.saadatniaki,khan@ece.}@tufts.edu. His work is supported by an NSF Career award # CCF-1350264.

‡ AK and GS are with Laboratoire Ampère, Université de Lyon, École Centrale de Lyon, 69134 Ecully Cedex, France, {first name}.{last name}@ec-lyon.fr. Their work is supported by a grant from la Région Rhône-Alpes.

tracking. Ref. [START_REF] Nascimento | Continuous-time distributed estimation[END_REF] extends adaptive diffusion models to continuous-time setups. Moreover, a recursive algorithm is presented in [START_REF] Zhou | Distributed estimation for moving target based on state-consensus strategy[END_REF] for the distributed estimation of a moving target under switching interconnection topologies, where the stability of the algorithm is analyzed under mild observability and connectivity assumptions. Of significant relevance to this paper is Ref. [START_REF] Korniienko | Performance control for interconnection of identical systems: Application to PLL network design[END_REF], which tackles the problem of control law design for interconnected systems, ensuring global stability and a certain pre-specified performance criteria using the notion of dissipativity, focusing on certain input-output properties of dynamical systems related to the conservation, dissipation, and transport of energy loosely referred to as dissipative properties [START_REF] Brogliato | Dissipative Systems Analysis and Control: Theory and Applications[END_REF].

In this paper, we assume that the no agent possesses enough measurements in its neighborhood to be able to estimate the entire state-vector on its own. In this sense, standard estimation techniques are not applicable and we propose a consensus plus innovation estimator that recovers the observability at each agent with the help of consensus on the neighboring predictions. To keep the estimator design simple and develop an estimator that requires minimal computation and coordination, we use a static gain that is same for both of the consensus and the innovation terms. In this context, we use the machinery developed in Ref. [START_REF] Korniienko | Performance control for interconnection of identical systems: Application to PLL network design[END_REF] and cast the distributed estimation problem in a similar decentralized control framework, where based on dissipativity characterization of the dynamical systems, we show that the gain design problem can be reduced to two quasiconvex optimization problems under Linear Matrix Inequality (LMI) constraints. In particular, we apply the bounded-real lemma and matrix elimination procedure to obtain existence conditions and a design procedure for a static estimator gain, which guarantees stable error dynamics while also satisfying a given performance level in terms of an upper bound on the H ∞ -norm of the global error transfer function. Our design is based on the input-output approach, i.e., given a power spectral density (PSD) specification on the input, the H ∞ -norm is closely related to the PSD of the output, see e.g., [START_REF] Zhou | Essentials of robust control[END_REF], [START_REF] Skogestad | Multivariable Feedback Control: Analysis and Design[END_REF].

We now describe the rest of the paper. Section II formulates the problem and describes related concepts on dissipativity and the input-output approach. Sections III and IV study the estimation error analysis and present our main results on the design of a static estimator gain ensuring stability and performance specifications, respectively. Section V illustrates our results and Section VI concludes the paper.

Notation: The superscripts ' ' and ' * ' denote transpose and complex-conjugate transpose, respectively. The N × N identity matrix is denoted by I N and the n×m zero matrix is denoted by 0 n×m . The indexes are dropped if the context is clear. The diagonal aggregation of two matrices A and B is denoted by diag(A, B). With matrix G partitioned into four blocks G 11 , G 12 , G 21 , and G 22 , G K denotes the Redheffer product, [START_REF] Doyle | Review of LFTs, LMIs, and µ[END_REF], i.e.,

G K = G 11 + G 12 K (I -G 22 K) -1 G 21 ; similarly, K G = G 22 + G 21 K (I -G 11 K) -1 G 12 .
We use {α} to denote the real part of a complex number α ∈ C. For a stable LTI system represented by [START_REF] Zhou | Essentials of robust control[END_REF]. The ith eigenvalue of a real square matrix P is denoted by λ i (P ). For symmetric matrices, ' ' (' ') denotes positive semi-definiteness (positivedefiniteness). Finally, R + represents the set of nonnegative real numbers.

transfer matrix G, G ∞ sup { G(jω) | ω ∈ R } denotes its H ∞ -norm,

II. PROBLEM FORMULATION

Consider the following CT-LTI dynamics:

ẋ(t) = Ax(t), t ≥ 0, ( 1
)
where t is the continuous-time variable, x(t) ∈ R n is the state at time t, and A is the system matrix. The measurements for this system are distributed over a network of agents, communicating according to a graph G = (V, E), where V = {1, 2, . . . , N } is the set of agents and E is the set of edges (links) between the agents. The neighborhood at the ith agent is defined as N i {i} ∪ {j|(i, j) ∈ E} and each agent i has the following observation model:

y i (t) = H i x(t) + r i (t), (2)
where y i (t) ∈ R pi is the local measurement, H i is the local observation matrix, and r i (t) is the zeromean additive noise whose PSD is bounded by 1. Letting A ∈ {0, 1} N ×N to denote the corresponding adjacency matrix and

D = diag(|N 1 \ {1}|, • • • , |N N \ {N }|
) to denote the degree matrix, the normalized adjacency matrix is defined as Ā = D -1 A. Given the dynamical system described by Eq. ( 1) and local observations in Eq. ( 2), the goal is to estimate the state x(t) of the system in a distributed manner. We assume that the system is globally observable, i.e., with the system matrix, A, and the collection of all measurement matrices, H i 's. However, no agent is assumed to be observable in its neighborhood, i.e., the system matrix, A, is not necessarily observable with the measurement matrices, H j , j ∈ N i , at any agent in the network. Clearly, due to a lack of (neighborhood) observability, no agent can implement the standard Kalman filtering equations to estimate the state, x(t), even if it collects all of the measurements in its neighborhood. The distributed estimation problem is thus to estimate the state of the system in this setting. Finally, we note that there is no process noise in Eq. ( 1). Nevertheless, the state estimation problem considered here is of relevance because the initial conditions are not known, the measurements are corrupted by noise, and no agent is locally observable.

We propose a variant of the Networked Kalman-type Estimator (NKE), initially introduced in [START_REF] Khan | On connectivity, observability, and stability in distributed estimation[END_REF], which follows a prediction + innovation framework based on state and output exchange among agents with an additional consensus term from inter-agent interactions to address the distributed nature of the problem. Not only does this force agreement but it also diffuses information on the state variables for which there may be no measurement in the neighborhood:

˙ x i (t) = A x i (t) predictor + j∈Ni K ij ( x j (t) -x i (t)) state exchange: consensus + j∈Ni B ij H j y j (t) -H j x i (t) measurement exchange: innovation , (3)
where x i (t) is the local state estimate at agent i, K ij ∈ R n×n are consensus gain matrices, and B ij ∈ R n×n are innovation gain matrices between agents i and j. Since each agent estiamtes the same dynamical system, we assume the gain matrices to be constant over the consensus and innovation terms and across each agent, i.e., K ij = K |Ni\{i}| , B ij = K, where K ∈ R n×n ; this also serves as a baseline with respect to the solutions obtained with higher degrees of freedom.

A. Error Dynamics

Defining the local estimation error as e i (t)

x i (t)-x(t), it is straightforward to show that ėi(t) = ˙ xi(t) -ẋ(t) = Aei(t) + K |Ni \ {i}| j∈N i [ej(t) -ei(t)] -K   j∈N i H j Hj   ei(t) + K j∈N i H j rj(t) external input due to noise . (4)
The matrix K can be thought of as a static estimator gain, designing which is the main goal of this paper such that the estimation error is stable at each agent and meets certain performance criteria with regard to noise rejection. The block diagram corresponding to the graphical representation of the local estimation error dynamics is depicted in Fig. 1.

Fig. 1: Local error dynamics

The signals enclosed in the dashed octagon represent the estimation errors from the neighborhood N i of agent i, constructing the auxiliary variable e Ni (t) as labeled in Fig. 1. The aggregate auxiliary signal from the entire network, denoted as ē(t), reveals how the network topology affects the global estimation error dynamics and is given by

ē(t) =    eN 1 (t) . . . eN N (t)    = ( Ā ⊗ In)    e1(t) . . . eN (t)    ( Ā ⊗ In)e(t). (5)
Next we investigate the effect of agent measurements on the global estimation error dynamics:

ȳ(t) =    y N1 (t) . . . y N N (t)    = D H    e 1 (t) . . . e N (t)    , (6)
where D H is the block-diagonal matrix with j∈Ni H j H j being its ith block, and

r(t) =    r N1 (t) . . . r N N (t)    =    j∈N1 H j r j (t) . . . j∈N N H j r j (t)    , (7)
which is the external input due to measurement noise. Note that T i (s) (as labeled in Fig. 1) can be rewritten as G(s)K(I + G(s)K) -1 based on the push-through rule from [START_REF] Skogestad | Multivariable Feedback Control: Analysis and Design[END_REF]Chapter 3]; consequently, we can express T i (s) in terms of the Redheffer product P(s) K with

P(s) = P 11 (s) P 12 (s) P 21 (s) P 22 (s) = 0 G(s) I -G(s) . (8)
For the global error process, concatenating u i (t) (as labeled in Fig. 1) in a vector u(t) results in the internal signal

u(t) = r(t) + ē(t) -ȳ(t) = r(t) + [( Ā ⊗ In) -DH ]e(t). ( 9
)
The resulting global error process with r(t) as the input and e(t) as the output is represented in Fig. 2. This system can be interpreted as the interconnection F of N identical subsystems T i (s) = T (s), mathematically denoted as (I N ⊗ T (s)) F , where e(s) = (I N ⊗ T (s)) u(s). As indicated in Fig. 2, the interconnection F is static; i.e., a gain matrix, which means that the relation [in Fig. 2] holds in both time and frequency domains. We now describe the main problem addressed in this paper.

Fig. 2: Global error dynamics

Problem 1: (Distributed Estimation): Given the above formulation, find the estimator gain K such that given the local dynamics P as defined in Eq. ( 8) and the interconnection matrix F in Fig. 2, the corresponding global error dynamics is stable and satisfies the performance specification

(I N ⊗ T (s)) F ∞ < η, (10)
for a given η > 0. This objective directly affects noise rejection properties of the global error by setting an upper bound η on the noise amplification from r to e, see Fig. 2. To this aim, we will design procedures based on the dissipativity properties of the interconnection matrix F and the local system T (s). The advantage of the dissipativity framework lies in the fact that it reduces the complexity of global stability and performance analyses by studying the local [sub]systems, which is the main motivation behind inputoutput approaches in the study of networked systems, see e.g., [START_REF] Moylan | Stability criteria for large-scale systems[END_REF]- [START_REF] Korniienko | Control law design for distributed multi-agent systems[END_REF]. The analysis is carried out in Sections III and IV for stability and performance, respectively. Preliminaries and essential results related to dissipativity are described next.

B. Preliminaries on Input-output Approach

Consider a general state space

ẋ = f (x, u), u ∈ U = R m , y = h(x, u), y ∈ Y = R p , (11)
where x ∈ R n is the state and (u, y) is the input-output pair. The concept of dissipativity describes the system by relating the internally stored energy of the system to a generalized energy supply rate function s(u(t), y(t)). There are many different notions of dissipativity introduced in various publications, see e.g., [START_REF] Brogliato | Dissipative Systems Analysis and Control: Theory and Applications[END_REF]Chapter 4] for a detailed discussion. In our analysis of distributed estimation, we use the following notion of dissipativity in frequency domain from [START_REF] Korniienko | Performance control for interconnection of identical systems: Application to PLL network design[END_REF]:

Definition 1 (Dissipativity): A causal, stable LTI system T is [strictly] (X, Y, Z)-dissipative,
for almost every ω ∈ R + .

The following result from [START_REF] Korniienko | Performance control for interconnection of identical systems: Application to PLL network design[END_REF] serves as the basis for our proposed solution to the problem of distributed estimation explained in the next two sections.

Theorem 1: Given η > 0, a stable interconnection F , a local plant P(s), and the real matrices X = X ≺ 0, Y , and Z = Z of appropriate dimensions, if there exist (i) a positive-definite matrix S ∈ R N ×N such that F is {diag(S ⊗X, -η 2 I), diag(S ⊗Y, 0), diag(S ⊗Z, I)}dissipative, and (ii) a local controller K such that T (s) = P(s) K is strictly {-Z, -Y , -X}-dissipative, then the large-scale system I N ⊗T (s) F , depicted in Fig. 2 is stable and

I N ⊗ T (s) F ∞ ≤ η.
Theorem 1 can be specialized to analyze only the internal stability by slightly modifying condition (i) as below.

(i ) a positive-definite matrix S such that M = F 11 , defined in Fig. 2, is {S ⊗ X, S ⊗ Y, S ⊗ Z}-dissipative.

In this paper, we assume diagonal matrices of the form

X = xI n , Y = yI n , Z = zI n , (13)
for the choice of X, Y , and Z in Theorem 1 [and its modified version] and find x, y, and z such that conditions (i) [(i )] and (ii) are both satisfied. Matrices of more general structures will be the focus of future studies.

III. ERROR STABILITY

Recall that our goal is to design the estimator gain matrix K such that the global error dynamics are internally stable and satisfy a given level of performance. In this section, we consider the stability of the error process, which ensures that the estimation error would vanish if no measurement noise is present across the network; the performance will be discussed in the next section. We only focus on the case where the system matrix A is not Hurwitz since otherwise, the error process is stable with K = 0 n . We now present our main result on the error stability.

Theorem 2 (Stability): Given the dynamics in Eq. ( 1) with non-Hurwitz system matrix A, the global estimation error process (I N ⊗T (s)) F , and M = F 11 = Ā⊗I n -D H , there exists an estimator gain matrix K ∈ R n×n ensuring internal stability if there exist a positive-definite matrix S ∈ R N ×N , and real values χ > 0 and c such that

M (S ⊗ I n S )M χ(cM -I nN ) S(cM -I nN ), (14) χ < min{ 1 c 2 , 1 (c -1) 2 }. ( 15
)
The proof of Theorem 2 is beyond the scope of this paper and is omitted because of space limitations. We have included a sketch of the proof here.

Proof Sketch : Directly applying the modified version of Theorem 1 and choosing x = -1, y = c, and z = 1 χ -c 2 in Eq. ( 13), Eq. ( 14) is equivalent to condition (i ). Moreover, condition (ii) in Theorem 1 can be expressed as

χ (T (jω) -cI) * (T (jω) -cI) -I ≺ 0, (16)
which is equivalent to T ∞ < 1 for the auxiliary system T (s) = √ χ (T (s) -cI). Applying Bounded-Real Lemma [START_REF] Boyd | Linear Matrix Inequalities in System and Control Theory[END_REF]Chapter 10] to the state-space representation of T (s) and introducing the change of variables Q = P K, T ∞ < 1 if and only if there exists a positive-definite matrix P ∈ R n×n such that

A P + P A -(Q + Q) + χI Q -cχI Q -cχI (c 2 χ -γ 2 )I ≺ 0, (17)
for γ 2 = 1. The proof is completed by deriving the equivalence between Eq. ( 15) and the feasibility of Eq. ( 17) through the application of the Elimination Procedure for Matrix Variables [START_REF] Boyd | Linear Matrix Inequalities in System and Control Theory[END_REF]Chapter 2] on Eq. ( 17).

Theorem 2 establishes the existence conditions for a stabilizing estimator gain matrix K. In the remainder of this section, we provide an algorithm to compute K based on this theorem and the resulting observations made in the following Remarks.

Remark 1: In the case of a single-input and single-output system T , Eq. ( 16) means that the Nyquist plot of the transfer function T (s) lies inside the circle centered at c+j0 with radius 1 √ χ . Similar interpretation applies to the multiple-input and multiple-output case as well. To relax this constraint for a given center c, the radius has to be maximized, which is equivalent to minimizing χ.

The key idea of the algorithm is based on the fact that minimizing χ, as mentioned in Remark 1, would strengthen the constraint in Eq. ( 14); therefore, Eqs. ( 14) and ( 16) cannot be simultaneously relaxed. As a result, the following optimization problem can be solved to satisfy the conditions of Theorem 2, leading to the computation of the estimator gain matrix K ensuring global stability.

minimize χ,S χ subject to S 0, S = S ⊗ I n , M SM χ (cM -I nN ) S (cM -I nN ) , (18)
which is a generalized eigenvalue problem for a fixed value of c (see Remark. 2), with S ∈ R N ×N and χ ∈ R + as the decision variables, [START_REF] Boyd | Linear Matrix Inequalities in System and Control Theory[END_REF].

Remark 2:

The center c has to be fixed in order to have a tractable optimization problem in Eq. ( 18); it must be chosen in a way that Eq. ( 14) is feasible. The choice of c can be based on the spectrum of M ; in fact, it can be shown that Eq. ( 14) is feasible as long as

|c| ≤ 1 2 max i λ i (M ) . ( 19
)
Remark 3: Having formulated the optimization in Eq. ( 18) as a generalized eigenvalue problem, we cannot include Eq. ( 15) as an explicit constraint; consequently, we proceed in the following iterative manner: Denoting the values obtained from the above procedure as c * and χ * , we use these values to verify whether the constraint in Eq. ( 15) is met, in which case the static estimator gain is computed as K = (P f ) -1 Q f , where P f and Q f are solutions to the feasibility problem in Eq. ( 17) for γ = 1, c = c * , and χ = χ * . If not, we fix a new value for c from the interval obtained in Eq. ( 19) and repeat the earlier steps. Given a sufficiently large number of samples from the interval in Eq. ( 19) (gathered in the sample set C), one approach is to do a linear search on c until the conditions of Theorem 2 are met or the sample set is exhausted, the latter implying that a static estimator gain K, ensuring stability of the global error, was not found.

The above procedure to compute the static estimator gain is summarized in Algorithm 1.

Remark 4: In case of Hurwitz (stable) system matrix A, the condition in Eq. ( 15) simplifies to χ < 1 c 2 .

Algorithm 1 Compute Static Estimator Gain K given M , C = {samples from the interval in (19)}, repeat c * ← a value in C χ * ← solve the generalized eigenvalue problem in Eq. ( 18)

for c = c * if χ * < min{ 1 c * 2 , 1 (c * -1) 2 } then P f , Q f ←
solve the feasibility problem in Eq. [START_REF] Brogliato | Dissipative Systems Analysis and Control: Theory and Applications[END_REF] for

c = c * , χ = χ * K ← P f -1 Q f return K end if until C has been exhausted. return Static estimator gain K not found IV. PERFORMANCE
In this section, we present our main result on H ∞ performance characteristics of the global error process and compute estimator gain matrix K satisfying the H ∞ performance objective specified in Eq. ( 10) in addition to stability.

χ < min{ 1 c 2 , 1 (c -1) 2 }. (21)
Proof Sketch : The proof is based on Theorem 1 and follows the same procedure as the proof of Theorem 2.

Similar to the case of internal stability, we present a systematic procedure to compute an estimator gain matrix K satisfying the performance specification in Eq. [START_REF] Khan | Collaborative scalar-gain estimators for potentially unstable social dynamics with limited communication[END_REF] in addition to internal stability based on Theorem 3. An analogous observation to the one made in Remark 1, reveals that condition (ii) in Theorem 1 (equivalent to Eq.( 16)) and Eq. ( 20) cannot be simultaneously relaxed, leading to the following optimization problem:

minimize χ,S χ subject to S 0, S = S ⊗ I n , Φ      0 0 0 0 0 0 0 0 0 0 1 χ S 0 0 0 0 I      Φ χΦ        1 χ S 0 0 0 0 η 2 χ I 0 0 0 0 0 0 0 0 0 0        (22)
with Φ as defined earlier.

To formulate the optimization in Eq. ( 22) as a generalized eigenvalue problem, we first solve the stability problem. Using the resulting values for c, χ and (I N ⊗ T (s)) F ∞ (as a baseline for η), we introduce the change of variable ς = 1 χ S and parameter β = η 2 χ . The resulting optimization problem, presented in Eq. ( 23), is a generalized eigenvalue problem with ς ∈ R N ×N and χ ∈ R + as the decision variables, [START_REF] Boyd | Linear Matrix Inequalities in System and Control Theory[END_REF].

minimize χ,ς χ subject to ς 0, ς = ς ⊗ I n , Φ    0 0 0 0 0 0 0 0 0 0 ς 0 0 0 0 I    Φ χΦ     ς 0 0 0 0 βI 0 0 0 0 0 0 0 0 0 0     Φ. (23)
In the optimization problem above, β acts as a tuning parameter; if the optimization problem has no solution, β should be increased to relax the last constraint of the problem. Also note that for a given value of β, minimizing χ results in the minimization of the upper bound on noise amplification gain η according to the relation η 2 = χβ.

V. SIMULATIONS

In this section, we illustrate the concepts introduced in this paper and evaluate the effectiveness of the proposed static gain distributed estimator in terms of internal stability and H ∞ performance, experimentally. We consider a network of N = 20 agents communicating over an undirected graph based on a 3-circulant graph with a few more links generated at random. The agents' goal is to estimate the state of an unstable n = 3-dimensional CT-LTI dynamical system evolving according to the model in Eq. ( 1) with a randomly chosen initial condition x(0) using Eq. (3). We have considered a bare minimum measurement setup where 3 non-neighboring agents i 1 , i 2 , and i 3 are selected at random such that each of these three agents measures one of the system states in a mutually exclusive and collectively exhaustive manner while the rest make no measurement. Note that in this setup, observability is nonexistent from both individual and neighborhood perspectives.

In the remainder of this section, we demonstrate the stability and performance characteristics of the static gain estimator through a series of Monte Carlo simulations.

A. Internal Stability

To illustrate the stability of the estimaion error process, agents are assumed to have perfect (noise-free) observations; i.e., r i (t) = 0. Given zero initial estimate x i (0), each agent makes an estimate of the system state using Eq. ( 3) for t ∈ [0 , 10] seconds. This procedure is repeated over 100 Monte Carlo simulations, each with randomly generated sets of non-Hurwitz system matrix A, communication graph G, and observer agents {i 1 , i 2 , i 3 }. Consistent with the theory, the distributed estimator in Eq. ( 3) results in asymptotically stable error as illustrated in Fig. 3a, where the normalized 2-norm of the global error, over each trial and averaged over 100 Monte Carlo simulations, has been plotted.

B. Performance

In this subsection, the performance level η in Eq. ( 10) is set equal to 11 and the performance characteristics of the global error process are demonstrated in presence of noise signal r(t) defined in Eq. ( 7), for the case where agent observations are subject to Additive White Gaussian Noise with PSD equal to σ 2 = 0.25. This procedure is repeated over 100 Monte Carlo simulations with a fixed set of system matrix A, communication graph G, and the bare minimum observation scheme explained above. The resulting normalized 2-norm of the global error signal is plotted in Fig. 3b, averaged over 100 Monte Carlo simulations. We note that the error process is stable and the H ∞ -norm of the resulting global error process (I N ⊗ T (s)) F ∞ is 10.6803, satisfying the performance specification in Eq.( 10).

VI. CONCLUSIONS

In this paper, we describe a static estimator gain design method to solve the distributed estimation problem for CT-LTI systems monitored by a network of agents. In this context, we express the problem of distributed estimation as corresponding optimization problems and apply techniques based on the input-output approach to perform stability and performance analyses. We derive existence conditions for a static estimator gain matrix resulting in stable error dynamics, i.e., estimation error vanishing over time in the absence of measurement noise; additionally, we derive existence conditions regarding static estimator gain design, guaranteeing desirable performance characteristics for the error process in terms of noise rejection. The approach described in this paper serves as the foundation of future investigation towards extending the current work to dynamic estimator design meeting certain frequency-dependent performance objectives in addition to stability and noise rejection.

  where X = X , Y , and Z = Z are real matrices of appropriate dimensions such that X Y Y Z is full-rank, if there exists ε = 0 [> 0] such that I T (jω) * X Y Y Z I T (jω) + εI 0.

Theorem 3 (

 3 Performance): Given the global error process (I N ⊗T (s)) F , there exists an estimator gain matrix K ensuring performance objective (I N ⊗ T (s)) F ∞ < η, if there exist a positive-definite matrix S ∈ R N ×N , and real values χ > 0 and c such that Φ

Fig. 3 :

 3 Fig. 3: Normalized error 2-norm e(t) 2 N .


A Korniienko
email: anton.korniienko@ec-lyon.fr

X Bombois
email: xavier.bombois@ec-lyon.fr

H Hjalmarsson
email: hakan.hjalmarsson@ee.kth.se

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

Hierarchical Robust Analysis for Identified Systems in Network

Keywords: Identification for Control, Interconnected systems, Network, Multi-Agents

This paper considers worst-case robustness analysis of a network of locally controlled uncertain systems with uncertain parameter vectors belonging to the ellipsoid sets found by identification procedures. In order to deal with computational complexity of large-scale systems, an hierarchical robustness analysis approach is adapted to these uncertain parameter vectors thus addressing the trade-off between the computation time and the conservatism of the result.

INTRODUCTION

In this paper, the problem of worst-case robustness analysis of a network of locally controlled uncertain Linear Time Invariant (LTI) subsystems is under consideration. The uncertainty of each subsystem is an uncertain real vector that belongs to an ellipsoid: an uncertainty set in the model parameter space typically obtained after identification.

This work is motivated by recent technological advances in Microelectronics, Computer Sciences, Robotics, and related topics in the field of the Multi-Agent systems [START_REF] Cao | An overview of recent progress in the study of distributed multi-agent coordination[END_REF]. The control of these network systems is usually decentralized and in order to compute controllers achieving high performance level, the model of the subsystems needs to be known. An efficient method to build the appropriate models is system identification [START_REF] Ljung | System identification[END_REF]. However, due to the presence of the noise and since the identification experiment is limited in time, the model parameters can only be identified within some prescribed uncertainty region which is typically an ellipsoid. For these reasons, in order to ensure that the computed controllers achieve the performance not only for the nominal identified model but for the true network system, it is important to take into account these uncertainties. The evaluation of the uncertainty effects on the system stability and performance is called robustness analysis.

The large scale of today's systems raises additional challenges on identification, controller design as well as on the robustness analysis. This paper focuses on the robustness analysis in the context of large-scale network systems.

In the 80's-90's, µ-analysis [START_REF] Doyle | Analysis of feedback systems with structured uncertainties[END_REF]; [START_REF] Safonov | Stability margin of diagonaly perturbed multivariable feedback systems[END_REF] was developed to investigate the performance of LTI systems in the presence of structured uncertainties. The performance is evaluated in the frequency domain Skogestad

This work originally was supported by a grant from the Région Rhône-Alpes.

and [START_REF] Skogestad | Multivariable Feedback Control, Analysis and Design[END_REF]. This approach is based on the computation of the structured singular value µ of the frequency dependent matrices, which was proved to be NPhard [START_REF] Braatz | Computational complexity of µ calculation[END_REF]. Fortunately, upper bound on µ can be efficiently computed, see [START_REF] Fan | Robustness in the presence of mixed parametric uncertainty and unmodeled dynamics[END_REF], ensuring a certain level of performance with some conservatism. By efficient, it is understood that the computation time is bounded by a polynomial function of the problem size [START_REF] Garey | Computers and Intractability: A guide to the theory of NP-Completeness[END_REF]. An adaptation of these results to classes of the uncertainties obtained by identification can be found in [START_REF] Bombois | Robustness analysis tools for an uncertainty set obtained by prediction error identification[END_REF]; [START_REF] Scorletti | Improved efficient analysis for systems with uncertain parameters[END_REF]; [START_REF] Barenthin | Identification for control of multivariable systems: Controller validation and experiment design via LMIs[END_REF].

Nevertheless, even if the computation time of the µ upper bound is polynomial function, it can still be important in the case of uncertain large-scale systems. The purpose of this paper is to extend the results [START_REF] Bombois | Robustness analysis tools for an uncertainty set obtained by prediction error identification[END_REF]; [START_REF] Scorletti | Improved efficient analysis for systems with uncertain parameters[END_REF]; [START_REF] Barenthin | Identification for control of multivariable systems: Controller validation and experiment design via LMIs[END_REF] to the context of large-scale interconnected systems, addressing the trade-off between computation time and conservatism. To do so, we adapt the hierarchical robustness analysis approach of [START_REF] Dinh | Embedding of uncertainty propagation: Application to hierarchical performance analysis[END_REF][START_REF] Dinh | Convex hierarchical analysis for the performances of uncertain large-scale systems[END_REF]; [START_REF] Laib | Phase IQC for the hierarchical performance analysis of uncertain large scale systems[END_REF][START_REF] Laib | Hierarchical robust performance analysis of uncertain large scale systems[END_REF], initially proposed in [START_REF] Safonov | Propagation of conic model uncertainty in hierarchical systems[END_REF], to the class of uncertainties obtained from system identification. A similar subject is presented in [START_REF] Bombois | Optimal identification experiment design for the interconnection of locally controlled systems[END_REF]. The main contribution of this paper is, however, a deeper investigation of the robustness analysis aspects, allowing, in contrast to [START_REF] Bombois | Optimal identification experiment design for the interconnection of locally controlled systems[END_REF], for several types of embedding and their combinations.

The next section of the paper formulates the problem under consideration, while the third section presents the main result of the paper, the hierarchical analysis approach. The computationnal complexity of the hierarchical analysis is discussed in section four while the fifth section is dedicated to the numerical illustration example. The last section concludes the paper.

Notations We denote by R and C the sets of real complex numbers and by H M the transfer function

M 22 + M 21 H (I -M 11 H)
-1 M 12 with M ij being appropriate partitions of M and standing for the Redheffer star product. The matrix   

X 1 0 0 0 . . . 0 0 0 X N   
is denoted as diag i (X i ) with (block-)diagonal elements X i (i = 1, ..., N ). For a complex number y, we denote yy * by y 2 while σ(A) denotes the maximal singular value of a complex matrix A. The notation ( ) * + X stands for X * + X for any expression X right next to the sum sign.

PROBLEM STATEMENT Let us consider a network of

N mod Single-Input Single- Output (SISO) subsystems S i (i = 1...N mod ) operated in closed loop with a SISO decentralized controller K i (i = 1...N mod ): S i (θ i ) : y i (t) = G i (s, θ i )u i (t) + v i (t) (1)
u i (t) = K i (s)(r i (t) -y i (t)) (2) r(t) = A ȳ(t) + B ref (t)
(3) where s, in order to keep the discussion as general as possible and to consider both cases, defines the Laplace variable s in the continuous time domain or the shift variable z in the discrete time domain. The vector θ i ∈ R n θ i represents the parameter vector of the ith system. We will distinguish hereafter between a variable θ i ∈ R n θ i , its unknown true value, θ i,0 ∈ R n θ i , and its estimated value, θi ∈ R n θ i . Let us also define

θ = [θ 1 , . . . , θ N ] T ∈ R n θ , θ 0 ∈ R n θ and θ ∈ R n θ :
the stacked version of the previous parameter vectors, with n θ = i n θi . The signal u i is the input applied to the system S i and y i is the measured output. This output is made up of a contribution of the input u i and of a disturbance term v i that represents both process and measurement noises and is modeled as a stochastic random process [START_REF] Ljung | System identification[END_REF]. The different true systems are thus described by transfer functions G i (s, θ i,0 ). Moreover, the vector v ∆ = (v 1 , v 2 , ..., v N mod ) T is assumed to have mutually independent components v i .

In this paper, the interconnection form used in formation control or multi-agent systems (see e.g. [START_REF] Fax | Information flow and cooperative control of vehicle formations[END_REF]; [START_REF] Korniienko | Performance control for interconnection of identical systems: Application to PLL network design[END_REF]) is under consideration. Therefore, the signal r i is a locally available reference signal that will be computed via (3). The matrices A, called the normalized adjacency matrix, and B in (3) represent the interconnection (flow of information) physically present in the network and they can be easily obtained for any interconnection topologies. Furthermore, r, ȳ are defined in the same way as v above. A possible main global objective of the network could be the tracking performance: each output y i (t) has to approach in a specified time the reference signal: ref (t) which is generally only available (throughout r i ) at one or a few nodes of the network, which is defined by the matrix B.

As an example, let us consider the network in Fig. 1 with N mod systems connected in a chain, all of the form (1) and all with a decentralized controller K i , see (2). Using (3), it is possible to define the local tracking error signals e i = r i -y i and it can be proven that such To do so, let us define general performance specifications that cover the expressed tracking performance objective but also other additional specifications. Performance input w and output z are added with a (possible dynamic) interconnection matrix M such that r z = M ȳ w (4)

Different components of the matrix M depend on the information flow in the network, i.e. matrices A and B, as well as on the specific performance measure, as will be detailed in Section 4. In this article, we focus on the performance specifications expressed in the frequency domain, see [START_REF] Skogestad | Multivariable Feedback Control, Analysis and Design[END_REF]. For this purpose, let us further define the local, independent from the network, transfer function T i and the global transfer function of the network T w→z between local (r i → y i ) and global ( w → z) signals respectively. Based on (1), ( 2) and (4) the following expression are obtained :

T i (s, θ i ) = K i (s)G i (s, θ i ) 1 + K i (s)G i (s, θ i ) T w→z (s, θ) = diag i (T i (s, θ 1 )) M
The global performance specification will be deemed satisfactory if: ∀ω, σ (T w→z ( , θ)) < W (ω) (5) where defines jω in the continuous time domain or e jω in the discrete time domain.

It is thus necessary to design (or redesign) the local controllers in order to ensure (or improve) the network performance and respect (5) with θ = θ 0 . However, since θ 0 is unknown, it will be necessary to identify a model for each of the systems S i (θ i,0 ). We assume that there is an identification procedure leading to a consistent parameter vector estimate θi of each subsystem true parameter vector θ i,0 as well as an estimate of the corresponding covariance matrices P θi . Such an identification procedure exists in open or closed-loop for each module independently, see [START_REF] Ljung | System identification[END_REF]; [START_REF] Barenthin | Identification for control of multivariable systems: Controller validation and experiment design via LMIs[END_REF], or when the modules are connected to the network [START_REF] Bombois | Optimal identification experiment design for the interconnection of locally controlled systems[END_REF]. It implies with some probability that the true parameter vector θ i,0 belongs to some uncertainty set U i defined as :

U i = {θ i | (θ i -θi ) T P -1 θi (θ i -θi ) < χ} (
6) with a constant χ given the probability level we would like to ensure and the number of elements in θ i,0 .

We also assume that there is a design procedure allowing to compute local controllers K i (s) such that the nominal global transfer function T w→z (s, θ), with θ = θ an estimate of θ 0 , respects the frequency dependent bound (5). Such design procedures could be found in [START_REF] Scorletti | An LMI approach to decentralized H ∞ control[END_REF]; [START_REF] Korniienko | Performance control for interconnection of identical systems: Application to PLL network design[END_REF].

Of course since θ is not necessarily equal to θ 0 this will not necessarily ensure the constraint (5) for the true system. In order to ensure the performance of the true system, in this article we would like to solve the following worst-case robustness analysis problem. Since θ i,0 ∈ U i for all i, it is possible to ensure (5) with θ = θ 0 by computing the worstcase gain of T w→z ( , θ), evaluated in terms of maximum singular values, ∀θ i ∈ U i . Similarly to the robustness analysis approaches [START_REF] Doyle | Analysis of feedback systems with structured uncertainties[END_REF]; [START_REF] Safonov | Stability margin of diagonaly perturbed multivariable feedback systems[END_REF]; [START_REF] Fan | Robustness in the presence of mixed parametric uncertainty and unmodeled dynamics[END_REF], this computation will be performed frequency by frequency ∀ω ∈ Ω = {ω 1 , . . . , ω Nω }. Even though it is possible to propose a method of an appropriate frequency gridding choice, in this paper, this question is not under consideration. We therefore assume that the properties ensured ∀ω j ∈ Ω imply that they are ensured ∀ω ∈ R. Problem 2.1. Given system (1)-( 3),(4), given uncertainty sets (6), compute for each ω j ∈ Ω:

min θi∈Ui(i=1...N mod ) γ(ω j ) subject to σ (T w→z ( j , θ)) < γ(ω j ) ( 7
)
where j is defined similar to .

If the minimal solution of the previous problem respects γ(ω j ) ≤ W (ω j ) for all ω j , then the computed controllers ensure that the true system T w→z ( , θ 0 ) respects the frequency dependent bound in (5) and thus the global performance.

Problem 2.1 is close to the well-known problem of worstcase robustness analysis (or µ-analysis procedure) from the Robust Control Community [START_REF] Skogestad | Multivariable Feedback Control, Analysis and Design[END_REF]. However the uncertainty sets (6), representing ellipsoids in parameter space, are not the traditional ones considered in this field. The adaptation of traditional worst-case robust analysis methods to the case of the uncertainty set obtained from the identification can be found in [START_REF] Bombois | Robustness analysis tools for an uncertainty set obtained by prediction error identification[END_REF]; [START_REF] Scorletti | Improved efficient analysis for systems with uncertain parameters[END_REF]; [START_REF] Barenthin | Identification for control of multivariable systems: Controller validation and experiment design via LMIs[END_REF]. However direct application of these results in the case of a large-scale network system, i.e. when N mod is large, is not possible due to the high system complexity implying prohibitive computation time. As was mentioned in the introduction, the main contribution of this paper is to extend these methods to the network context i.e. to derive tractable robustness performance analysis conditions while keeping computation time reasonable.

HIERARCHICAL ANALYSIS APPROACH

Keys ideas

The main idea of the hierarchical approach [START_REF] Dinh | Embedding of uncertainty propagation: Application to hierarchical performance analysis[END_REF][START_REF] Dinh | Convex hierarchical analysis for the performances of uncertain large-scale systems[END_REF]; [START_REF] Laib | Phase IQC for the hierarchical performance analysis of uncertain large scale systems[END_REF][START_REF] Laib | Hierarchical robust performance analysis of uncertain large scale systems[END_REF] is to decompose the network into two or more hierarchical levels and to perform the robustness analysis level by level by propagating the analysis results from one level to another. For some network systems such decomposition appears naturally, as for example for the system under consideration in this paper :

(i) local hierarchical level : subsystem dynamics T i (s, θ i ) defined by ( 1) and ( 2) and (ii) global hierarchical level : the global information exchange (3) and (4). The robustness analysis at each hierarchical level allows to embed the subsystem dynamics with a possibly complex non-convex (and non-linear) dependence on the uncertainty, into a much simpler subsystem description with a convex dependence on the uncertainty. We will call it the embedding procedure in the sequel. Then in the next hierarchical level, the subsystem is replaced by this simple description and the procedure is repeated once again until reaching the last hierarchical level. The last step consists in the worstcase robustness analysis based on the propagated subsystem descriptions in order to evaluate the global network performance i.e. solve the Problem 2.1. The complexity and time computation reduction is ensured thanks to the embedding procedures and by the fact that all embeddings at each hierarchical level are independent and thus can be easily performed in parallel.

In this paper, a two level hierarchical structure (local and global) is under consideration. Before formalizing this approach separately for the local and global hierarchical levels, let us first define what we mean by subsystem dynamics and subsystem description.

Since the performance measure in this paper is expressed in the frequency domain, see (5), the subsystem dynamics are defined by the structured frequency response set T s i (ω) of the subsystem transfer function at frequency ω:

T s i (ω) = {T i ( , θ i ) | θ i ∈ U i } (8)
The subsystem description in turn is defined by an uncertainty set T i (x i (ω), y i (ω), z i (ω)) of complex numbers ∆ i (ω) ∈ C that respects a frequency dependent quadratic constraint imposed by

x i (ω) ∈ R, y i (ω) ∈ C, z i (ω) ∈ R: T i (x i (ω), y i (ω), z i (ω)) = {∆ i (ω) | ∆ i (ω) 1 * x i (ω) y i (ω) y i (ω) * z i (ω) ∆ i (ω) 1 ≤ 0 (9)
Let us introduce the following definition characterizing the frequency response of a system : Definition 3.1. (Dissipativity). An LTI system H(s) is {x(ω), y(ω), z(ω)} -dissipative at ω for some x(ω) ∈ R, y(ω) ∈ C, z(ω) ∈ R, if its frequency response H( ) respects the following quadratic constraint at ω:

H( ) 1 * x(ω) y(ω) y(ω) * z(ω) H( ) 1 ≤ 0.
If the following additional constraint is imposed on x(ω), then the corresponding quadratic constraint defines a convex set :

x(ω) ≥ 0. ( 10
)
Please note that, in the case of x(ω) > 0, by Definition 3.1 and the Schur complement [START_REF] Boyd | Linear Matrix Inequalities in Systems and Control Theory[END_REF], the following constraint is implied: y 2 (ω) ≥ x(ω)z(ω). When x(ω) = 0, no constraint is imposed on y(ω) and z(ω). In order to reduce the computational complexity, the convexity constraint (10) will be used in the sequel.

If each subsystem T i (s, θ i ) is {x i (ω), y i (ω), z i (ω)} -dissipative for some frequency dependent x i (ω), y i (ω), z i (ω) and for all θ i ∈ U i and ∀ω, we then obtain the following embedding

T s i (ω) ⊂ T i ((x i (ω), y i (ω), z i (ω)), ∀ω;
and the frequency responses of the uncertain subsystems T 1 (s, θ i ), . . . , T N mod (s, θ i ) generated by varying θ i ∈ U i , can be replaced in the global hierarchical level by the corresponding subsystem description T i ((x i (ω), y i (ω), z i (ω)).

Of course, since the set T i ((x i (ω), y i (ω), z i (ω)) is in general larger than the set T s i (ω) the result of the corresponding worst-case analysis might be conservative. In order to reduce this conservatism, it is important to choose suitable x i (ω), y i (ω), z i (ω) for each subsystem. It is also possible to compute several complementary triplets x k i (ω), y k i (ω), z k i (ω) for k = 1 . . . N d defining therefore N d dissipativity properties for each subsystem. It allows to define for each subsystem a basis of dissipativity properties (a set of subsystem descriptions) and propagate it to the global hierarchical level. Such a suitable choice in the context of the uncertainty set (6) obtained through an identification procedure is presented in the next subsection while Subsection 3.3 presents how the embeddings are combined and propagated in a global hierarchical step in order to efficiently solve Problem 2.1. It is clear that the more dissipativity characterizations are used for each subsystem, the more the conservatism is reduced. Of course, the price to pay for this is the increase of computation time. For this reason it is important to find appropriate triples x k i (ω), y k i (ω), z k i (ω) at each hierarchical step ensuring as tight embedding as possible.

Local Step

In this subsection we present how to efficiently compute different dissipativity triplets x, y, z at a given frequency ω j ∈ Ω such that an uncertain system T (s, θ i ) is {x(ω j ), y(ω j ), z(ω j )} -dissipative, ∀θ i ∈ U i with U i defined in (6).

For this purpose let us define the following factorization of the transfer function T (s, θ i ), suitable for the system identification [START_REF] Bombois | Robustness analysis tools for an uncertainty set obtained by prediction error identification[END_REF]:

T (s, θ i ) = e(s) + Z N (s)θ i 1 + Z D (s)θ i (11)
with θ ∈ R n θ and then present the following Lemma. Lemma 1. Given the uncertain LTI system T (s, θ i ) in (11), it is {x(ω j ), y(ω j ), z(ω j )} -dissipative for all θ i ∈ U i and for given ω j , x(ω j ) ∈ R, y(ω j ) ∈ C, z(ω j ) ∈ R respecting (10), if and only if (i) in the case of x(ω j ) > 0 :

-α(ω j ) λ(ω j ) λ * (ω j ) -A 1 (ω j ) -ξ(ω j )B + jX (ω j ) ≤ 0 (12) (ii) in the case of x(ω j ) = 0 : ( ) * + y * (ω j )A 2 (ω j ) + A 1 (ω j )z(ω j ) -ξ(ω j )B + jX (ω j ) ≤ 0(13) with λ(ω j ) = Z N ( j ) + y(ω j ) x(ω j ) Z D ( j ) e( j ) + y(ω j ) x(ω j ) , A 1 (ω j ) = Z * D ( j )Z D ( j ) Z * D ( j ) Z D ( j ) 1 , A 2 (ω j ) = Z * D ( j )Z N ( j ) Z * D ( j )e( j ) Z N ( j ) e( j ) , α(ω j ) = y 2 (ω j )
x 2 (ω j ) -z(ω j )

x(ω j ) , B =

P -1 θ i -P -1 θ i θi -θT i P -1 θ i θT i P -1 θ i θi -χ and some ξ(ω j ) ≥ 0 ∈ R, X (ω j ) = -X T (ω j ) ∈ R n θ ×n θ .
Proof. See the detailed proof in Korniienko et al. (2018) Please note that the sufficiency of Lemma 1 can be proved using the result of [START_REF] Dinh | Convex hierarchical analysis for the performances of uncertain large-scale systems[END_REF] (see Corollary 2.2).

As is shown in the proof of Lemma 1, the result of [START_REF] Dinh | Convex hierarchical analysis for the performances of uncertain large-scale systems[END_REF] is adapted to the case of uncertain vectors that belong to an ellipsoid which recovers sufficient and necessary conditions of {x(ω j ), y(ω j ), z(ω j )} -dissipativity. This lemma is an extension of the robustness analysis result of [START_REF] Bombois | Optimal identification experiment design for the interconnection of locally controlled systems[END_REF] and will be used to generate different types of embeddings.

We will now consider two types of embedding: the disc and the band embedding, and formulate a convex optimization problem to compute them. Please note that thanks to Lemma 1, it is possible to study other types of embedding, as for example cone embedding [START_REF] Laib | Phase IQC for the hierarchical performance analysis of uncertain large scale systems[END_REF], half planes etc. Disc Embedding Given system T (s, θ) in ( 11), its frequency response set {T ( j , θ i )

| θ i ∈ U i } is embedded in a disc set at ω j if |T ( j , θ i ) -c(ω j )| ≤ ρ(ω j ), ∀θ i ∈ U i ( 14
)
where c(ω j ) ∈ C is the center of the disc and ρ(ω j ) ∈ R is its radius, see [START_REF] Dinh | Convex hierarchical analysis for the performances of uncertain large-scale systems[END_REF]. The size measure of this embedding is the radius of the disc, and the problem of the computation of the tightest embedding can be formulated as follows, assuming appropriate gridding Ω. Problem 3.1. Given system (11) and its uncertainty set ( 6), compute for each ω j ∈ Ω: min ρ 2 (ω j ) s.t. ( 12) holds and ( 15)

x(ω j ) = 1, y(ω j ) = -c(ω j ), z(ω j ) = c 2 (ω j ) -ρ 2 (ω j ).
Proof. See Korniienko et al. (2018) for the proof. 2

Band Embedding

Given system T (s, θ i ) in ( 11), its frequency response set 16) where n(ω j ) ∈ C is the complex number which defines the vector -→ n = [Re(n), Im(n)] T giving the band orientation in complex plain (it is perpendicular to both band hyperplanes) and a 1 (ω j ), a 2 (ω j ) ∈ R are the signed distances of the two band hyperplanes to the origin multiplied by |n|, see [START_REF] Dinh | Convex hierarchical analysis for the performances of uncertain large-scale systems[END_REF]. The size measure of this embedding is the band width d(ω j ) = a 1 (ω j )-a 2 (ω j ) and the problem of computation of the tightest embedding can be formulated as follows, assuming again appropriate gridding Ω. Problem 3.2. Given system (11), its uncertainty sets (6), compute for each ω j ∈ Ω: min n(ωj ),a1(ωj ),a2(ωj )

{T ( j , θ i ) | θ i ∈ U i } is embedded in a band set at ω j if ∀θ i ∈ U i 2a 2 (ω j ) ≤ ( ) * + n * (ω j )T ( j , θ i ) ≤ 2a 1 (ω j ), (
a 1 (ω j ) -a 2 (ω j ) subject to (16)
This problem is efficiently solved as follows.

Theorem 3. (Band embedding). Problem 3.2 is solved by the following LMI optimization problem: min a1(ωj ),a2(ωj ),n(ωj )

a 1 (ω j ) -a 2 (ω j ) (17)
such that (13) holds with x 1 (ω j ) = 0, y 1 (ω j ) = n(ω j ), z 1 (ω j ) = -2a 1 (ω j ) and ( 13) holds with

x 2 (ω j ) = 0, y 2 (ω j ) = -n(ω j ), z 2 (ω j ) = 2a 2 (ω j ).

Proof. See Korniienko et al. (2018) for the proof. 2

Global Step

In this subsection, we assume that for all ω j ∈ Ω and for each subsystem T i (s, θ i ), several embeddings are found in the local step. We thus obtain N d dissipativity triplets

x k i (ω j ), y k i (ω j ), z k i (ω j ) for k = 1, .
. . , N d , for each subsystem i = 1, . . . , N mod and for all ω j ∈ Ω. The next theorem allows to compute an upper bound γ U B (ω j ) on the maximum amplification γ(ω j ) of Problem 2.1. Theorem 4. Given system (1)-( 3),(4), a frequency ω j and given 8) and ( 9))

x k i (ω j ), y k i (ω j ), z k i (ω j ) such that T s i (ω j ) ⊂ T i ((x k i (ω j ), y k i (ω j ), z k i (ω j )) for k = 1, . . . , N d , i = 1, . . . , N mod (see (
The upper bound γ U B (ω j ) on the maximum amplification γ(ω j ) of Problem 2.1 is the solution of the following LMI optimization problem:

min γ2 (ω j ),T k ω ,(k=1...N d ) γ2 (ω j ) s.t. M ( j ) I * N (γ 2 (ω j )) M ( j ) I > 0, with (18)
N (γ 2 (ω j ))

∆ =      T k ω Z k d 0 0 -I T k ω Y k d 0 0 0 * T k ω Y k d 0 0 0 T k ω X k d 0 0 γ2 (ω j )I     
with real strictly definite positive diagonal matrices for k

= 1 . . . N d , T k ω = diag i (τ k i (ω j ))
, and

X k d = diag i (x k i (ω j )), Y k d = diag i (y k i (ω j )), Z k d = diag i (z k i (ω j ))
Proof. See Korniienko et al. (2018) for the proof. 2

The detailed analysis of the computational complexity of the proposed Hierarchical approach can be found in [START_REF] Laib | Hierarchical robust performance analysis of uncertain large scale systems[END_REF]; Korniienko et al. (2018). The next section rather focus on the numerical example.

NUMERICAL EXAMPLE

Let us now consider an illustration example of an Automated Highway System (AHS): a platoon of autonomous cars following external reference signals as in [START_REF] Seiler | Disturbance in vehicle strings[END_REF]. Each car's simplified model dynamics is described by (1), with G i (s, θ) = ki s 2 (τis+1) and true parameter vector θ i,0 = [τ i , k i ] T where τ i , k i were randomly chosen around 0.105 and 0.95 respectively with uniform ±10% distribution. Each system is controlled by the same initial decentralized controller K init (s) = 2s+1 0.05s+1 taken from [START_REF] Seiler | Disturbance in vehicle strings[END_REF], see (2). There are N mod = 5 cars in the network which are allowed to exchange information according to bidirectional chain topology, see [START_REF] Seiler | Disturbance in vehicle strings[END_REF], as depicted in Fig. 1 and defined by (3).

The main objective of the network is that each car follows a ramp reference signal ref (t), available only for the first car, shifted by a constant value δ i = iδ, ∀i to avoid colisions, while keeping string instability (oscillation propagation through the network) limited [START_REF] Seiler | Disturbance in vehicle strings[END_REF] has a slope of +40 dB/dec at low frequency range, then the tracking performance is ensured, see [START_REF] Skogestad | Multivariable Feedback Control, Analysis and Design[END_REF]. Moreover, a lower gain ensures a better tracking speed and the resonance peak limitation reduces the effects of string instability [START_REF] Seiler | Disturbance in vehicle strings[END_REF].

The maximal singular value of the true system T w→z (s, θ o ) with initial controller K init is represented by orange dashdotted line in Fig. 2. In order to improve the tracking performance of the network and to reduce the oscillation effects provoked by the string instability, let us impose the frequency constraint (5) with W (ω) represented in Fig. 2 by the red dashed line. To satisfy this constraint, first an identification procedure is performed leading to a consistent parameter vector estimate θi of each subsystem true parameter vector θ i,0 as well as an estimation of the corresponding covariance matrices P θi ensuring (6). Due to the presence of a double integrator in the car transfer function model, this identification experiment has to be performed in closed loop with a stabilizing controller either independently for each module (see [START_REF] Ljung | System identification[END_REF]; [START_REF] Barenthin | Identification for control of multivariable systems: Controller validation and experiment design via LMIs[END_REF]) or in the network (see [START_REF] Bombois | Optimal identification experiment design for the interconnection of locally controlled systems[END_REF]). The results of the latter method are presented in Fig. 3 where the controllers were chosen as K i (s) = K init (s), ∀i. Different discrete-time white noise excitation signals of length N id = 1000, sampling time T s = 0.01 sec and variance 10 are added via a zero order hold to the references r i of each closedloop systems T i (s, θ i,0 ). The measured discrete signal y i is also perturbed by generated mutually independent white noise discrete signals v i with variance of 4 each modeling the measurement noise effects. A standard, predictionerror identification criterion is used, see [START_REF] Ljung | System identification[END_REF]. Notice that in this example the continuous transfer function parameters k i and τ i and the corresponding covariance matrices could be directly identified since the car transfer function model is rather simple. It ensures that the nominal global transfer function T w→z (s, θ), with θ = θ, respects the frequency dependent bound in (5).

Our problem is now to efficiently test if the constraint is satisfied by the true system by solving Problem 2.1 for properly chosen Ω. To do so, the proposed hierarchical approach is used. The results of the local step embeddings for the first system and at 0.15 Hz are presented in Fig. 4 where the borders of the minimum radius disc embedding (green full circle) and of the tightest band (red full lines) are presented. For the sake of illustration, we show the borders of the structured uncertainty set T s 1 (red dots), the estimated T 1 ( θ1 ) (blue cross) and the true T 1 (θ 1,0 ) (black round) value of the corresponding frequency responses evaluated at ω = 0.15 Hz. Notice that disk center c(ω) = T 1 (jω, θ 1,0 ) = T 1 (jω, θ1 ). The results are found by solving the LMI optimization problems (15) and ( 17). Similar results are obtained for other subsystems and other frequencies from Ω. The global step analysis results are presented in Fig. 5 for two cases : computed γ U B based on the propagation of (i) disc embedding only N d = 1 (blue rounds) and of (ii) disc and band embeddings N d = 3 (red dots). Fig. 5 also presents some Monte-Carlo samples i.e. the maximal singular value of T w→z (s, θ) for randomly chosen θ i ∈ U i . As we can see, the worst-case bounds are respected. Surprisingly even though the disc embedding set is much bigger than the intersection of disc and band sets (see Fig. 4), the overall upper bound γ U B is not improved a lot, see Table 1. It is due to the fact that, in this application, the phase uncertainty information, mostly captured by the band embedding, is much less important than the gain uncertainty information, mostly captured by the disc embedding. The corresponding computation times are also given in Table 1 for both serial and parallel computation of local embeddings. Finally, maximal singular values of the true system T w→z (s, θ o ) with the new controller are represented by the blue solid line in Fig. 2. As we can see, the global performance is ensured for the true system.

CONCLUSIONS

In this paper we proposed robustness analysis method adapted to the uncertainty sets constructed by identification in a network context. The type of network in this system is usual in the literature of multi-agent systems and the size of the network plays a crucial role in the robustness analysis complexity. In order to manage the trade-off between the computation time and the precision of the obtained result, the hierarchical robustness analysis approach was proposed and illustrated in the case of SISO subsystems. Future extension is the MIMO subsystem case with an appropriate choice of hierarchical structure (with possibly more than two hierarchical levels) in order to even better address the mentioned trade-off.

  ρ(ωj ),c(ωj ) ρ(ω j ) subject to(14) This problem is efficiently solved by the following theorem. Theorem 2. (Disc embedding). Problem 3.1 is solved by the following LMI optimization problem: min ρ 2 (ωj ),c(ωj )

Fig. 2 .

 2 Fig. 2. Maximal singular value of the true system T w→z (s, θ o ) for initial controller (orange dash-dotted line), improved controller (blue solid line) and imposed frequency constraint W (ω) (red dashed line).

Fig. 4 .

 4 Fig. 4. Local step embedding results. Borders of structured uncertainty set T s 1 (red dots), of the minimum radius disc embedding (green full circle), of the tightest band (red full lines), circle center (green dot), estimated frequency response (blue cross) and true frequency response (black round).

Fig. 5 .

 5 Fig. 5. Global step analysis results. Upper bounds computed by (i) propagation of disc embedding only (blue rounds), by (ii) propagation of disc and band embeddings (red dots), Monte-Carlo samples of maximal singular value of T w→z (s, θ), for some θ i ∈ U i .

  . It can be shown that this tracking performance specification is equivalent to the ability of each car to track the same ramp signal ref (t) ensuring that all local tracking errors e i = r i -y i go to zero in steady-state. If the maximum singular value of T w→z (s, θ o )

	As a consequence, let us define performance input w(t) =
	ref (t) and performance output z(t) = r(t) -ȳ(t). It thus
	determines the interconnection topology (4) with M =
	A B A -I B	.

Table 1 .

 1 Hierarchical analysis results

		disc only	disc + band	difference
	Dissipativity number	N d = 1	N d = 3
	γ U B @ 0.13 Hz	-11.83 dB	-12.04 dB	1.8%
	γ U B @ 0.15 Hz	-12.64 dB	-12.92 dB	2.2%
	γ U B @ 0.17 Hz	-14.28 dB	-14.44 dB	1.1%
	Overall Time	15.33 sec	19.04 sec	-24.2%
	Overall Time (Parallel)	11.85 sec	14.43 sec	-21.8%


Daniel Pizarro
email: pizarro@depeca.uah.es

Manuel Mazo
email: mazo@depeca.uah.es

Enrique Santiso
email: santiso@depeca.uah.es

Hideki Hashimoto
email: hashimoto@iis.u-tokyo.ac.jp

Mobile Robot Geometry Initialization from Single Camera

. Such proposal relies on a proper initialization of pose and structure information of the robot. The present paper proposes a solution to initialization which consists of retrieving 3D structure and motion of a rigid object from a set of point matches measured by the camera. A batch Structure from Motion (SFM) approach is proposed along a short path. By incorporating odometry information available in the robot, the ambiguity generated by a single view in the solution is solved. We propose to describe robot's motion and image detection as statistical processes in which the uncertainty is properly modelled. Using a Gaussian equivalence of the processes involved, the SFM cost function is expressed as a Maximum Likelihood optimization. The paper shows the improvements of the approach in the presence of the usual odometry drift noise, compared with those using Euclidean distance as a likelihood. The proposed method is assessed on synthetic and real data.

Introduction

The usage of external cameras to obtain the pose of a robot is inspired in the concept of 'Intelligent Space' [START_REF] Lee | Intelligent space and mobile robots[END_REF]. The aim is to provide the environment with intelligence and to reduce its necessity in the robots. As a consequence the space is composed by a set of sensors (e.g. Cameras ), a distributed intelligence, and a network layer, which allows communication with agents (e.g. robots). Under this framework, robot devices derive their localization and the high level navigation tasks from the space.

In this paper the usage of a single camera, external to the robot is proposed as a sensor device for robot localization. The ability of operating with a single sensor allows the system to work with few cameras, reducing heavily the cost and avoiding complex multiple camera configurations. There exist several approaches with the same aim that the present paper, by using artificial landmarks, either active [START_REF] Hada | Trajectory tracking control of a nonholonomic mobile robot using igps and odometry[END_REF] or passive ones [START_REF] Morioka | Human-following mobile robot in a distributed intelligent sensor network[END_REF]. In other works a geometrical model [START_REF] Sogo | Acquisition of qualitative spatial representation by visual observation[END_REF], is learnt previously to the tracking task in a supervised step. The common point of many of the proposals comes from the fact that high amount of prior knowledge is used for a posterior tracking.

In [START_REF] Pizarro | Simultaneous localization and structure reconstruction of mobile robots with external cameras[END_REF], a sequential algorithm that recovers the pose (with respect a global coordinate origin) and structure (set of 3D points from robot's structure) of a mobile robot is presented. The main novelty of this proposal is that it uses only as prior information the rigidity assumption in the geometry and the odometry readings present in the robot. The approach is based on retrieving the robot's pose by using a motion model and a sparse geometrical model which is learnt online. Such method depends heavily on a proper initialization of pose and structure before applying the online method.

In this work the initialization or self-calibration required in the online method is solved using one camera and odometry information. The initialization consists of a batch sequence of images of the robot taken during a short path which are synchronized with the odometry readings. In this paper a Structure from Motion (SFM) algorithm is proposed to solve the pose and structure by minimizing a likelihood function of the residuals in image frame. As the position estimation with odometry has a growing error the algorithm must tackle properly the uncertainty in the odometry metric reference.

In [START_REF] Chen | Self-calibration of environmental camera for mobile robot navigation[END_REF] the authors obtain the orientation and position of a camera with respect to a robot by using artificial markers and a set of predefined motions for calibration. The error in odometry position estimation and its influence in the solution is studied and the predefined motions are selected to minimize calibration error.

In the present paper we propose to describe robot's motion and image detection as statistical processes in which the uncertainty is properly modelled. Using a Gaussian equivalence of the processes involved, the SFM cost function is expressed as a Maximum Likelihood optimization. The main contribution of the paper is to show the improvements of this approach in the presence of the usual odometry drift noise, compared with those using Euclidean distance as a likelihood [START_REF] Triggs | Bundle adjustment -a modern synthesis[END_REF]. Besides, in this paper the natural structure of the robot is retrieved jointly with its pose.

We state the problem in §2. In §3 the likelihood distribution is analytically obtained from measurements and inputs. The final M.L cost function is presented in §4. Results on synthetic and real images are presented in §5. Finally, conclusions are presented in §6.

Problem Statement

As was commented before, the problem to solve is to achieve Structure from Motion by using one camera and the given odometry information which has growing error with robot's motion. To formalize the particularities of the problem, in this section it will be roughly commented on the notation and information involved.

Firstly, a set of discrete processes concerning pose, structure and measurements are presented. The number of time samples k N is considered to be finite and the proposed method runs as a batch process, so all the measurements and inputs are available from the beginning.

Robot's pose at time k is described by a vector x k . Usually for 3D motion (6 D.O.F) the vector is composed of 3 position components and 3 orientation angles x k = (r x k , r y k , r z k , α, β, γ). For wheeled robots which its motion lie on a plane, the vector pose x k is reduced to 3 components (r x k , r y k , φ). Motion model x k = f (x k-1 , u k ) obtains actual position with respect to previous time and some input u k given by odometry (i.e. angular speed and linear speed of the robot).

The geometry of the robot is composed by a sparse set of N 3D points

M = {m 1 , • • • , m M }
referred from a local coordinate origin described by robot's pose x k . The rigidity in geometry involves that the points M are static in time and no temporal subindex is necessary. Function m i x = t(x k , m i ) uses actual pose x k to express m i in the global coordinate origin.

The camera sensor follows a perfect "pin-hole" described by its 3x4 projection matrix P , which encodes intrinsic and extrinsic parameters. The camera projection model is expressed by the non-homogeneous transformation y = h(m x , P ) which converts a 3D point m x expressed in a global coordinate origin into its 2D projection y in the image plane using camera parameters P . The set of measurements which correspond to structure points is encoded in vector

Y k = {y i k i = 1, • • • , M } where y i k = h(t(x k , m i ), P ).
Each measurement in image plane is obtained by using a fiducial detector which also is able to establish correspondences between images.

The objective proposed in this paper is to retrieve the initial pose x 0 and the set M given a set of image measurements Y 1:NK and odometry input information u 1:N k .

In the next sections a statistical model of all processes is proposed. A random process is defined by its P.D.F and it is expressed in bold typography. (i.e. x k = N ( xk , Σ k ) is the Gaussian random process and x k is a realization of it).

Uncertainty Modeling

In this section, the objective is to obtain the statistical description of the processes stated in §2. By describing motion and image measurements as random processes, the uncertainty present in them is properly tackled and the odometry unbound error is directly modeled. In this paper it is shown how such representation improves the solution obtained in the reconstruction algorithm compared with a simple quadratic cost function.

The statistical models appear by including uncertainty noise variables into projection (h) and motion (f ) models. We propose to approximate all statistical processes using Gaussian density distributions, which gives a meaningful description of the solution (mean) and uncertainty (covariance). The main difficulties come from the fact that the motion and projection models are in general non-linear equations, which transform the input Gaussian process and involves solving integral equations.

A first order approximation of the non linear function is used for propagating Gaussian processes. This approach has been widely used in Extended versions of Kalman Filtering [START_REF] Davison | SLAM with a single camera[END_REF]. It is known to suffer from bias [START_REF] Van Der Merwe | The unscented particle filter[END_REF], however it offers an analytical solution for the mean and covariance which will be necessary to minimize the likelihood cost function.

In the following sections the processes involved will be tackled in detail.

Motion uncertainty

Initial solution for pose x 0 is represented by a Gaussian process x 0 = N (x 0 , Σ 0 ) with a mean vector x0 and a covariance matrix Σ 0 . Between time samples, x k is consider a Markov Process with a transition kernel p(x k |x k-1 , u k ) ruled by the motion model f . For simplicity, motion model used in this paper corresponds to a simple wheeled differential robot which moves over a ground plane. Pose parameters x k = (r k , φ k ) will be composed of position in the ground plane r k with respect a global coordinate origin O and orientation φ k :

x k = x k-1 +   (V l k + W 1 k ) cos(φ k-1 + Ω k + W 2 k ) (V l k + W 1 k ) sin(φ k-1 + Ω k + W 2 k ) Ω k + W 2 k   , (1)
Input from odometry is represented by the vector u k = (V l k , Ω k ) with linear and angular speed components. Noise from odometry readings is modeled by a Gaussian process

W k = (W 1 k , W 2 k ) which is statistically independent given x k-1 . Gaussian approximation of the distribution p(x k |x 0 ) = N ( xk , Σ k )
is recovered by propagating the last pose p(x k-1 |x 0 ) and the actual odometry u k through the transition kernel p(x k |x k-1 , u k ). The result is a Gaussian process with continuously growing variance usually presented in Dead-Reckoning [START_REF] Borenstein | Measurement and correction of systematic odometry errors in mobile robots[END_REF] processes (See Figure 1).

Measurement uncertainty

Measurement uncertainty depends on the method used to track fiducial points in the image. Measurements y i k are jointly attached to form a vector Y k which will be described as a measurement process Y k . Each separate measurement y i k is contaminated by additive Gaussian noise v i k = N (0, Σ v ). The noise components are independent between each other and are included in the noise process V k .

Likelihood distribution p(Y 1:kN |x 0 , M) represents the joint distribution of the whole set of measurements in the sequence. It is obtained by combining projection model from every measurement y i k into a global measurement equation h:

Y 1:kN = h(x 1:kN , Y 1:kN , V 1:kN ) V k = {V i k i = 1, • • • , M } h = h(t(x 1 , m 1 ), P ) + V 1 1 , • • • , h(t(x kN , m M ), P ) + V M kN ∈ Ê 2•M•kN (2)
A Gaussian approximation of p(Y 1:kN |x 0 , M) = N ( Ŷ1:kN , Σ Y ) is then obtained. It is used in this paper as the likelihood from which the solution to initial pose and structure is retrieved.

Usually, Σ Y ∈ Ê nY ×nY where n Y = 2M k N , is a full correlated covariance matrix, in which all cross-correlation terms are present. The cross correlation terms will play an important role in the reconstruction step.

Assuming independence between each x k , the resulting covariance will be approximated by a box-diagonal one, where only terms from the same frame are correlated. In the same way, assuming that at time k there exists independence between different measurements, the resulting covariance matrix will be box-diagonal where only coordinates of the same measurement can be correlated. (See (3)).

In the next section, the differences and the importance of correlated terms will be explained. In §5, some synthetic experiments are run to test the approaches selected for Gaussian propagation.

Σ Y = 0 B B B B B B B B B B B B B B @ • • • . . . . . . . . . • • • 1 C C C C C C C C C C C C C C A ≈ 0 B B B B B B B B B B B B B B @ . . . 1 C C C C C C C C C C C C C C A ≈ 0 B B B B B B B B B B B B B B @ . . . 1 C C C C C C C C C C C C C C A (3)

Maximum Likelihood Solution

Rigid structure from motion algorithms are at the moment in a mature state [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]. Thanks to the advances made in the last ten years, it can be considered a solved problem given that the correspondence between frames is correct. The usage of non-linear 'Bundle Adjustment' [START_REF] Triggs | Bundle adjustment -a modern synthesis[END_REF] with robust cost functions allows to obtain a solution for the SFM problem. 'Bundle Adjustment' makes use of reprojection errors to get robot initial pose x 0 and the set M:

min x0,M N k k M i ρ(y i k -h(t(x k , m i ), P )), (4)
where function ρ strictly depends on the reprojection error statistical properties. In this paper the cost function ρ is naturally obtained by maximizing the logarithm of the joint distribution p(Y 1:kN |x 0 , M). By using Maximum Likelihood optimization, the cost function penalizes those residuals y i kh(t(x k , m i ), P ) which has more error due to odometry drift.

Maximum Likelihood weighted Least Squares

By taking logarithm, the likelihood function proposed p(Y 1:kN |x 0 , M), and approximated by a Gaussian function, is maximized by minimizing the following quadratic cost function:

min x0,M (Y 1:kN -Ŷ1:kN ) T (Σ Y ) -1 (Y 1:kN -Ŷ1:kN ) (5)
It must be noticed that both Σ Y and Ŷ1:kN are functions of the parameters to minimize x 0 and M.

It is straightforward to see that when the covariance matrix Σ Y is the identity, the resulting cost function is directly the one presented in (4) with square of residuals as ρ.

Expression ( 5) can be optimized by using any non-linear iterative optimization method such as Gauss-Newton or Region-Trust algorithms as Levenberg Mardquardt.

Sparse Hessian

Each basic step of optimization involve the calculation of the following operations: We define as Φ n = {x 0 , M} n as the value of the parameters to minimize at iteration number n.

Given Φ n-1 the solution for the increment Φ n = Φ n-1 + ∆Φ at iteration n is obtained by solving the following linear system:

J T (Σ Y ) -1 J∆Φ = J T (Σ Y ) -1 (Y 1:kN -Ŷ1:kN (Φ n-1 )), (6)
where J represents the first derivative of Ŷ1:kN evaluated at Φ n-1 .

Matrix H = J T (Σ Y ) -1 J which is a square matrix of the size of Φ represents the approximated Hessian.

At each iteration system ( 6) is solved by inverting H.

∆Φ = (H) -1 b, (7)
If P Y is a diagonal matrix, usually H have a sparse structure. This allows a very important improvement in the computational complexity for its inversion. The situation is specially critic when M contains a large number of 3D points.

As we can see in Figure 2, the inclusion of a full covariance Σ Y or only a time correlated one, destroys the natural sparsity of the problem. The direct effect is that the systems become of higher complexity.

Initial guess

As a non-linear iterative algorithm, the proposed optimization requires an initial guess for vector Φ. The assumption is that if robot height is small compared to the distance to the camera the points from robot structure can be supposed to lie on the ground plane. Using the set of measurements from the first frame Y 0 , and knowing the position of the ground plane, a set of 3D positions from global reference frame can be obtained, namely X 0 . The set M is referred to the initial pose x 0 . By using centroid of the set X 0 the position r 0 is easily obtained.

Results

In this section some results on synthetic generated data and real images are presented.

Synthetic data

The experiment conditions are the following:

• A set of M=30 points m i k initialized as a random Gaussian jointly with robot true initial position r 0 . The geometry has a deviation around the center of rotation of σ = 200(mm) • Camera calibration and ground position from a real configuration.

• Odometry noise variance Σ W is varied from σ 2 I 2x2 σ = 0.1, • • • , 20
(m/s,rad/s). • Measurement noise Σ V with standard deviation of around 0.1 pixels.

The following algorithms will be tested:

Real Data

A real system composed of a low cost camera and a real mobile robot is tested. The error initialization is the same used in the synthetic data.

The result of an initial calibration procedure is shown in Fig. 4 (b). As can be see five coplanar points are considered as initial state vector. The path required to obtain the estimation is shown in Fig. 4 (a).

Conclusions

In this paper a method for computing robot geometry and pose from single camera is presented, with the aim of properly initializing a sequential approach. The inclusion of odometry allows metric reconstruction with single camera. Knowledge about the statistical properties of odometry readings and camera detection processes allows to propose a Maximum Likelihood approach which penalizes properly information which has more uncertainty as the last positions given by a dead-reckoning process. The equivalent approach is a weighted quadratic cost functions. The reconstruction accuracy is proved to be higher in the case of M.L versions. As can be seen in Figure 3 error is higher when no statistical information is used. Between the two approximations of the covariance matrix, the one which has more cross-correlated terms (F.D) achieves better performance when the error grows. The main drawback is that the more cross correlated terms the less sparse the Hessian becomes, which increases problem complexity. As a conclusion, for low error in odometry the Measurement Box Diagonal (M.D) approach will have the same behavior with the advantage that the sparsity of the problem remains equal to a Element Diagonal approach (E.D). Real results with coplanar points are also presented, achieving good accuracy with real odometry readings and image matching. This work was supported by the Ministry of Science and Technology under RESELAI project (reference TIN2006-14896-C02-01).

  Fig. 1.

( a )Fig. 2 .

 a2 Fig. 2. Hessian Configurations

1 .•

 1 Element Diagonal covariance estimation (E.D): By supposing that Σ Y is diagonal. No weighted optimization. 2. Frame Box Diagonal covariance matrix (F.D): See Figure 2(a). 3. Measurement Box Diagonal covariance (M.D): See Figure 2(b). Experiments: Reconstruction mean error vs. Odometry noise variance. (Figure 3(a)) • Reconstruction mean error vs. % of outliers. (Figure 3 (b))

Fig. 3 .

 3 Fig. 3. Experiment 1