Mihai Pascaru
email: m_pascaru@yahoo.com

Keywords: glocalization, psycho-social impact, position groups, population fears, traditional and modern mining

Romanian glocalization: case study on the Roşia

Montană gold corporation mining project.

Introduction

The Romanian communist system was a quasi-closed one, as it is the case inside all communist regimes. Until the fall of the communist regime, in 1989, the connections to the globalization process had been limited to their commercial dimension. Everything was exported, including the basic food of the Romanian people; and import had a highly selective character. After 1989, in a more or less scheduled manner, a massive campaign Sitio web: http://territoriosposibles.fahce.unlp.edu.ar/ La Plata, 17 al 20 de octubre de 2012. 1

Equipo TAG Territorios Posibles UNLP-CONICET y Universidades asociadas del GDRI Groupe de Recherche CNRS INTI International Network of Territorial Intelligence Instituto IdhICS, Facultad de Humanidades y Ciencias de la Educación, Universidad Nacional de La Plata -CONICET meant to connect Romania to the globalization processes started. However, paradoxically, this process brought with it a significant amount of damage to the economy, including in the area of mineral resources exploitation. Essentially, within only 10 years, the most of the mines in Romania were closed. The reason was unprofitability. Moreover, after another 10 years, mining was almost forgotten or occasionally brought to public attention through the controversies surrounding possible foreign investments in the field. The Rosia Montana Gold Corporation (RMGC) mining project was and still is among the topics of long and heated controversy. This project is promoted by Rosia Montana Gold Corporation, a joint venture between the Romanian state company Minvest S.A. Deva (19,3%), Gabriel

Resources Limited from Canada (80%) and minority shareholders (0,7%). (Pascaru, 2007: 40). The project should have been implemented starting with the year 2000. The exploitation itself was scheduled to begin in 2005 and -based on the existent mineral reserves -it was meant to go on for a period of 17 years. Yet it has not started to this day.

At one point, the delay in the implementation of the project led to the intervention of Romanian president, Traian Basescu. He blamed this delay on the politicians' cowardice:

"If this project has been buried from 1997 until now, you have to know that one of the serious causes for this is the cowardice of politicians who wanted to spare themselves the headache. If the project would have been launched in 1997, the increase of the gold price would have found us in the middle of an exploitation process." (Mediafax, 29 august 2011).

To us, the RMGC project is an example of a vector of globalization that penetrated the open society that characterized Romanian after the fall of the communist regime.

Natale Ammaturo noted that, in the last decades, a large quantity of studies and research related to globalization were published and that it is difficult to organize the enormous amounts of literature dedicated to this topic. This literature passes on from

McLuhan's global village to discussions related to the globalization of the village (Ammaturo, 2004: 5) The relationship between global society and local societies is extremely complex.

Regarding this subject Schuerkens wrote: "For several decades, we have found the creations of new local social forms as an answer to the process of globalization… Thus glocalization permits the opening up of a new theoretical understanding of social transformation and change." (Schuerkens, 2003: 207).

When global and local overlap, there are three levels we can identify. These represent the levels of glocalization as well: 1) local-national level, 2) local-regional level and 3) local-community level (see Fig. 1).

Romanian Glocalization. Researchs at Rosia Montana

The vector of globalization and the position groups

Our first research on Rosia Montana was carried out in Corna and Bunta, two villages that would be taken off the map by the works in the project. Here the impact of the RMGC project, as a vector of globalization, was felt through the immergence of position groups at a small community level. The position groups in the village Corna, for instance, were built around a more or less explicit agreement or disagreement towards the RMGC project and its effects. Consequently the following factors have been considered: 1) agreement/disagreement to the idea of building a starter dam on the site of the Corna village; 2) the assessment of the projects' influence upon the area. Position groups were characterized according to the community matrix indices: 1) inter-knowledge index, 2) communication index, 3) mutual action index. These indices have values between 0 (total lack of knowledge, communication and action) and 1 (maximum amount of knowledge, communication and action). (Pascaru, 2003).

In 2003, in the village of Corna alone, 9,4% of the respondents stated that they agree, without regret, to the RMGC project's initiative to build in the Corna village a starter dam;

28,1% stated that they agree, but feel sorry about it and 29,7% do not agree. For the ones who agreed without regret the community matrix index was over the average (0,50), especially due to a high inter-knowledge index (0,82) and a high communication index (0,57). The ones who stated they did not agree also had a higher community matrix index due to a high mutual action index (0,19).

When asked "Do you think the project will have a positive influence upon this region?" 31,3% of the respondents answered YES and 48,8% answered NO. The community matrix index for those who answered YES was above the average (0,46), especially due to the inter-knowledge index (0,79), while the index value for the inhabitants who answered NO Sitio web: http://territoriosposibles.fahce.unlp.edu.ar/ La Plata, 17 al 20 de octubre de 2012. ISBN 978-950-34-0932-9

Equipo TAG Territorios Posibles UNLP-CONICET y Universidades asociadas del GDRI Groupe de Recherche CNRS INTI International Network of Territorial Intelligence Instituto IdhICS, Facultad de Humanidades y Ciencias de la Educación, Universidad Nacional de La Plata -CONICET was very close to the community average (0,73 -inter-knowledge index, 0,43communication index, 0,13 -mutual action index, 0,43 -matrix index).

In a provisional overview, we estimated that the following groups adopted a favourable position towards the RMGC project: 1) those who stated they agreed and had no regrets about the disappearance of the village and whose community matrix index was above average, especially due to their involvement in mutual actions; they represent a group whose resources could be harnessed to support the communitarian development for the project; 2) those who believed the RMGC project would have a positive influence upon the area, with a community matrix index above average largely due to the high value of the inter-knowledge index. The groups who adopted a more or less clear position towards rejecting the RMGC project were: 1) those who did not agree to the disappearance of the village, with an index value above average; 2) those who believed that the RMGC project will not have a positive influence in the area; the index value of the community matrix and its dimensions equalled the community average for this group, which can be qualified as a "mass group". The practical meaning of these classifications, as we have acknowledged at the time this research was made, could be that, in the case of a collective debate and decision, due to their position within the community, the groups who appeared to support the RMGC project would prevail upon the groups that rejected it. (Pascaru, 2007: 79-82).

Local alternatives to globalization

In 2011 we have coordinated a series of research in Rosia Montana and Abrud (a neighbouring town). The results of the research showed that the number of people who had heard about modern mining was approximately equal to the number of people who had not. Still, little over half of the interviewees (54,4%) gave affirmative answers. The majority of those who stated that they had heard about modern mining answered that it was different than traditional mining. The equivalence between the two types of mining could have been the result of the representations regarding the impact and/or the technology used in the process. 3,2% more interviewees from Abrud than those of Rosia Montana believed that there are differences between modern and traditional mining.

Sitio web: http://territoriosposibles.fahce.unlp.edu.ar/ La Plata, 17 al 20 de octubre de 2012. Equipo TAG Territorios Posibles UNLP-CONICET y Universidades asociadas del GDRI Groupe de Recherche CNRS INTI International Network of Territorial Intelligence Instituto IdhICS, Facultad de Humanidades y Ciencias de la Educación, Universidad Nacional de La Plata -CONICET (Rachieru et al., 2011: 10-11). Most answers to the question about advantages and disadvantages of traditional mining vs. modern mining contained references to the number of jobs created because of mining in the past, but also to the contribution of traditional mining to an increase of the GDP. In contrast, modern mining was only represented as a job source, where the bosses make money and where higher qualifications and responsibilities are required. In traditional mining, low work efficiency was considered a disadvantage for the economy but an advantage for the environment. However, according to the inhabitants, modern mining would have economical advantages offered by the new but this would draw higher risks of pollution. From a social point of view, advantages in the past were present in the comfort brought by a safe pension and high wages. Modern mining did not seem to promise anything like this. (Rachieru et al., 2011: 11-12). Overall, as an indirect reaction to globalization, supporting the advantages of traditional mining is in fact a manifestation of glocalization. The globalization, located in Rosia Montana, generates nostalgia.

Conclusions

The natural question at the end of this study is: What is the specific character of localising a powerful vector of globalization such as the Rosia Montana Gold Corporation project? By answering this question we will also draw the lines for the profile of Romanian glocalization. At a local-national level, we have seen that the cowardice of the politicians, as the Romanian president called it, can also be taken into account when defining Romanian glocalization. At a local-regional level, we have a paradox of rejecting globalization, even more so when its impact is less dramatic. To it, we can add the fast structuring of the position groups into organizations (associations pro and against the RMGC project) with their intense activity to delay the start of the exploitation process throughout the years. At present, there are data indicating the fact that a pro-project decision or a permanent rejection of the project would be more beneficial than an endlessly prolonged uncertainty, including for the inhabitants of Rosia Montana. In terms of Sitio web: http://territoriosposibles.fahce.unlp.edu.ar/ La Plata, 17 al 20 de octubre de 2012. Equipo TAG Territorios Posibles UNLP-CONICET y Universidades asociadas del GDRI Groupe de Recherche CNRS INTI International Network of Territorial Intelligence Instituto IdhICS, Facultad de Humanidades y Ciencias de la Educación, Universidad Nacional de La Plata -CONICET local identity preservation -a question brought forth by the position groups, but also by interventions outside the community -the RMGC project was constantly adjusted, firstly in order to preserve the historical centre of Rosia Montana, and then to value and preserve the rich archaeological patrimony. In terms of environmental protection, initially the RMGC planned rehabilitation projects at the end of the exploitation. These plans were changed for projects meant to save Rosia Montana from the pollution left behind by the traditional mining.

We will conclude this study with a hypothesis suggested by the case study on RMGC: the depth of glocalization is greater as we move from a macro (local-national) to a micro level (local-communal). We believe that this hypothesis can offer a generous opening to further research in the field of glocalization. BIBLIOGRAPHY Ammaturo, N. ( 2004

Fig. 1 .LOCAL

 1 Fig. 1. Levels of glocalization

  ), La dimensione della solidarità nella società globale, Milano, FrancoAngeli. Pascaru, M. (2003), Matricea comunitara. Cunoastere, comunicare si actiune comuna in satul contemporan, Cluj-Napoca, Presa Universitara Clujeana. Pascaru, M. (2007), Habitatul risipit de globalizare. Impactul psiho-social si rezidential al proiectului Rosia Montana Gold Corporation, Cluj-Napoca, Editura Argonaut. Rachieru, A.D., Petcu, D., Gherhes, V., Obrad, G., Suciu, S. (2011), Raportul cu privire la impactul evolutiei economice asupra calitatii vietii in localitatile Zlatna, Baia de Aries, Abrud si Rosia Montana, Alba Iulia, Centrul de Cercetare si Dezvoltare Socio-Economica "Muntii Apuseni". Schuerkens, U. (2003), Social Transformation Between Global forces and Local Life-Worlds: Introduction, Current Sociology, May/July 2003, Vol. 51, pp. 195-208, SAGE Publications. Van Der Bly C. E., M. (2005), Globalization: A Triumph of Ambiguity, in Current Sociology, November 2005, Vol. 53(6), pp. 875-893, SAGE Publications. *** "BASESCU: PROIECTUL DE LA ROSIA MONTANA A STAT INGROPAT DIN CAUZA Sitio web: http://territoriosposibles.fahce.unlp.edu.ar/ La Plata, 17 al 20 de octubre de 2012. ISBN 978-950-34-0932-9 Equipo TAG Territorios Posibles UNLP-CONICET y Universidades asociadas del GDRI Groupe de Recherche CNRS INTI International Network of Territorial Intelligence Instituto IdhICS, Facultad de Humanidades y Ciencias de la Educación, Universidad Nacional de La Plata -CONICET /WWW.MEDIAFAX.RO/SOCIAL/ BASESCU-PROIECTUL-DE-LA-ROSIA-MONTANA-A-STAT-INGROPAT-DIN-CAUZA-LASITATII-POLITICIENILOR-8668622. Sitio web: http://territoriosposibles.fahce.unlp.edu.ar/ La Plata, 17 al 20 de octubre de 2012. ISBN 978-950-34-0932-9 8

  Equipo TAG Territorios Posibles UNLP-CONICET y Universidades asociadas del GDRI Groupe de Recherche CNRS INTI International Network of Territorial Intelligence Instituto IdhICS, Facultad de Humanidades y Ciencias de la Educación, Universidad Nacional de La Plata -CONICET

	. While economists successfully defined globalization as open
	economy, Van Der Bly believed that, "sociologists might switch the debate and consider
	the idea of globalization as the idea of 'an open society', and discuss which empirical
	parameters indicate the level of 'openness' of a society." (Van Der Bly, 2005: 884).
	Sitio web: http://territoriosposibles.fahce.unlp.edu.ar/
	La Plata, 17 al 20 de octubre de 2012. ISBN 978-950-34-0932-9


Didier Georges
email: didier.georges@grenoble-inp.fr

Gildas Besançon
email: gildas.besancon@grenoble-inp.fr

Jean-François Dulhoste

A Decentralized Optimal LQ State Observer Based On An Augmented Lagrangian Approach

Keywords: Decentralized optimal LQ state observer, price decomposition-coordination, augmented Lagrangian

This paper is devoted to the design of a decentralized optimal batch LQ state observer for state estimation of large-scale interconnected systems, well suited for implementation on a sensor network. The here-proposed approach relies on both the use of an augmented Lagrangian formulation and a price-decomposition-coordination algorithm. The state estimation of an openchannel hydraulic system illustrates the effectiveness of this approach and is used to provide a comparison with alternative methods.

INTRODUCTION

The design of decentralized or distributed state estimation or control algorithms is an important research topic which has attracted a great interest for more than 30 years [START_REF] Looze | On decentralized estimation and control with application to freeway ramp metering[END_REF][START_REF] Sanders | Specific structures for large scale state estimation algorithms having information exchange[END_REF][START_REF] Siljak | Decentralized Control of Complex Systems[END_REF][START_REF] Khan | Distributing the Kalman filter for large-scale systems[END_REF][START_REF] Stankovic | Consensus based overlapping decentralized estimator[END_REF][START_REF] Farina | Moving-horizon partition-based state estimation of large scale systems[END_REF][START_REF] Hess | Distributed Kalman filter algorithms for self-localization of mobile devices[END_REF][START_REF] Maestre | Distributed receding horizon Kalman filter[END_REF][START_REF] Mansouri | Decentralized observers design for descriptor systems[END_REF][START_REF] Ishizaki | Hierarchical decentralized observer design for linearly coupled network systems[END_REF][START_REF] Martensson | A scalable modularized synthesis method for distributed Kalman filters[END_REF]. On the one hand, estimation of large-scale complex systems is of the greatest importance for monitoring applications in various engineering fields, such as power grids, road traffic networks, environmental systems (open-channel hydraulic systems, water supply networks, weather, ...). On the other hand, sensor networks [START_REF] Akyildiz | Wireless sensor networks: A survey[END_REF] are now recognized for being well suited for measurement, monitoring, tracking of distributed physical phenomena, such as environmental phenomena (weather, seismic events, wildfires, air / soil / river pollution [START_REF] Ghanem | Sensor Grids For Air Pollution Monitoring[END_REF], sound or vibration monitoring ...) or complex industrial systems. Sensor networks are defined as a collection of embedded sensors with communication capabilities. Large scale deployment of such ad hoc networks mainly relies on the availability of cheap embedded sensors. For that reason, it is necessary to consider algorithmic architectures capable of performing low complexity computations distributed on each sensor node. Here the attention is paid to sensor network-based monitoring applications for large-scale systems, in which local measurements are made and a local state estimation with low computational complexity is performed by each sensor agent while some computation results are provided to the other agents connected to it. Furthermore a systematic way to manage decentralized state estimation through the application of dual optimization theory and the deterministic interpretation of the continuous-time Kalman filter as described in [START_REF] Bornard | Observability and observers[END_REF] is provided. Unlike the approaches in [START_REF] Looze | On decentralized estimation and control with application to freeway ramp metering[END_REF][START_REF] Mansouri | Decentralized observers design for descriptor systems[END_REF], the proposed approach provides an exact decentralized solution to the centralized Kalman filter and does not require strong structural properties of the system. The here-proposed approach is related to some previous works in [START_REF] Maestre | Distributed receding horizon Kalman filter[END_REF], where dual decomposition is used. However the authors in [START_REF] Maestre | Distributed receding horizon Kalman filter[END_REF] only consider the case of coupling of the subsystems via the measurement outputs. The here-proposed approach also differs from approaches based on consensus (see [START_REF] Stankovic | Consensus based overlapping decentralized estimator[END_REF][START_REF] Khan | Distributing the Kalman filter for large-scale systems[END_REF]), since the local systems are never some overlapping subsystems. It can be compared to the approach proposed in [START_REF] Farina | Moving-horizon partition-based state estimation of large scale systems[END_REF] since decentralization leads to the solution of distributed small-size moving horizon estimators without overlapping. However to the best of our knowledge, such an augmented Lagrangian formulation is used for the first time. It will be also shown that the proposed observer algorithm can be used to derive a moving horizon state observer. The paper is now organized as follows: In section II, the state estimation problem is stated. The decentralized design, based on both the use of an augmented Lagrangian formulation and a relaxation technique, together with some backgrounds on dual optimization, is then presented in section III. Section IV will be devoted to the derivation of a moving horizon observer based on the decentralized algorithm. Section V proposes an environmental state estimation application that illustrates the effectiveness of the approach. In this section, the proposed approach will be also compared to alternative approaches. Finally, the paper ends with some conclusions.

PROBLEM STATEMENT

We consider the problem of designing a decentralized state observer for large-scale interconnected linear systems defined by:

ẋi (t) = A i x i (t) + B i u i (t) + D i N j=1 F ij x j (t), (1)
y i (t) = C i x i (t), i = 1, ..., N, (2)
where the full state vector x ∈ R n and the full input vector u ∈ R m are partitioned as

x = (x T 1 , x T 2 , ..., x T N ) T and u = (u T 1 , u T 2 , ..., u T N )
T respectively, where x i ∈ R ni and u i ∈ R mi , and each subsystem admits a measurement output y i ∈ R pi . The overall system is assumed to be observable. System (1)-( 2) may be equivalently defined as the following algebraic-differential system:

ẋi (t) = A i x i (t) + B i u i (t) + D i s i (t), (3)
s i (t) - N j=1 F ij x j (t) = 0, (4)
y i (t) = C i x i (t), i = 1, ..., N, (5)
where the s i 's represent some slack variables. We seek to derive an optimal Linear Quadratic (LQ) state observer with a finite horizon by solving the following optimal output tracking problem, based on the knowledge of

y i (t), u i (t), ∀t ∈ I = [0, T ]: min vi(t),ŝi(t),t∈I 1 2 N i=1 T 0 { ŷi (t) -y i (t) 2 R -1 i + v i (t) 2 Q -1 i }dt + 1 2 N i=1 xi (0) -x i 2 M -1 i (6) s.t. ẋi (t) = A i xi (t) + B i u i (t)i + E i v i (t) + D i ŝi (t), ŷi (t) = C i xi (t), ŝi (t) = N j=1 F ij xj (t),
where Q i , R i and M i , i = 1, ..., N are some symmetric positive-definite matrices of adequate dimensions, and E i 's are also matrices of appropriate dimensions. x i is the a priori most likely value of x i (0). This problem is an equivalent formulation to the continuous Kalman filter design problem which seeks the minimum variance state estimate of the overall state x of the system (see [START_REF] Bornard | Observability and observers[END_REF][START_REF] Kwakernaak | Linear optimal control systems[END_REF]), where Q i is interpreted as the covariance matrix of a zero mean gaussian noise affecting the state equation of x i through matrix E i , and R i is the covariance of a zero mean gaussian noise affecting the measurement output y i . M i is the covariance matrix of random initial state x i (0).

STATE OBSERVER DESIGN

The optimal output tracking problem ( 6) is closely related to the classical optimal trajectory tracking problem, except that the cost penalizes the initial states rather than the terminal states. For that reason, it appears to be very convenient to reverse time by considering the change (t → T -t) and interval I = [T, 0], and thus turning the optimal output tracking problem into: min ṽi(t),si(t),i∈I

1 2 N i=1 T 0 { C i xi (t) -ỹi (t) 2 R -1 i + ṽi (t) 2 Q -1 i }dt + 1 2 N i=1 xi (T ) -x i 2 M -1 i s.t. ẋi (t) = -A i xi (t) -B i ũi (t) -E i ṽi (t) -D i si (t), si (t) = N j=1 F ij xj (t), (7)
where xi (t

) = xi (T -t),ũ i (t) = u i (T -t), ṽi (t) = v i (T -t) and si (t) = ŝi (T -t).
Let us now consider an augmented Lagrangian formulation of the problem as follows (in the spirit of [START_REF] Cohen | Decomposition coordination methods in large scale optimization problems. The nondifferentiable case and the use of augmented Lagrangians[END_REF][START_REF] Georges | Distributed model predicted control via decomposition-coordination techniques and the use of an augmented lagrangian[END_REF] for instance):

L c (ṽ, s, μ) = 1 2 N i=1 T 0 { C i xi (t) -ỹi (t) 2 R -1 i + ṽi (t) 2 Q -1 i }dt + 1 2 N i=1 xi (T ) -x i 2 M -1 i + N i=1 T 0 < μi (t) + c 2 [s i (t) - N j=1 F ij xj (t)], si (t) - N j=1 F ij xj (t) > dt (8) s.t. ẋi (t ) = -A i xi (t) -B i ũi (t) -E i ṽi (t) -D i si (t),
where c > 0 is the coefficient of the augmented Lagrangian, which should be chosen large enough to ensure the existence of a saddle-point. ṽ, s, μ denote the vectors of all the ṽi 's, si 's and μi 's, respectively. < ., . > represents the usual scalar product. µ is the vector of Lagrange multipliers associated to the interconnection constraints. The interest of considering an augmented Lagrangian rather than the ordinary one may be found in the fact the problem is not strongly convex with respect to the slack variables s i . The augmented Lagrangian formulation may be viewed as a regularization technique which strongly convexifies the problem and therefore ensures the existence of a saddle-point of the constrained problem for non-strongly convex problems. As a consequence, the relaxation algorithms, such as Uzawa algorithm [START_REF] Arrow | Studies in Linear and Nonlinear Programming[END_REF], will converge by using a constant step gradient method, rather than a gradient method with "small steps" of type σ 1 , which tend to zero as the number of dual iterations tends to infinity (see [START_REF] Cohen | Decomposition coordination methods in large scale optimization problems. The nondifferentiable case and the use of augmented Lagrangians[END_REF]), and therefore, cannot be easily handled in practice. On the basis of the augmented Lagrangian defined by ( 9), the problem consists in finding a saddle-point (ṽ, s, μ ) of L c , defined by min ṽ,s

max μ L c (ṽ, s, μ) = max μ min ṽ,s L c (ṽ, s, μ) (9)
where w(μ) = min ṽ,s L c (ṽ, s, μ) is the so-called dual function, which has to be maximized.

Some background on dual optimization and augmented Lagrangians

In this paper, we consider iterative computation of a saddle-point solution of (9) by using a relaxation algorithm of the family of "Uzawa algorithms". In order to illustrate such an algorithm, let us consider the following constrained "separable" optimization problem:

min u J(u) = min ui N i=1 J i (u i ) (10)
s.t. θ(u) = N i=1 θ i (u i ) = 0. ( 11
)
and the ordinary Lagrangian

L(u, µ) = J(u)+ < µ, θ(u) >= N i=1 {J i (u i ) + µθ i (u i )} 1 steps ρ k with ρ k ≥ 0, +∞ i=0 ρ k = +∞ and +∞ i=0 (ρ k ) 2 < +∞.
where J is a convex functional defined on a Hilbert space U = U 1 ×...×U N and the θ i 's are some vectorial mapping from

U i to C. θ(u) = 0 means that θ(u) ∈ -C = -{0},
where C is the closed convex cone in C reduced to zero. In this case, the Uzawa algorithm becomes the so-called "price-decomposition-coordination algorithm" (PDCA) defined by:

(1) Initialization (k = 0): Choose u 0 i and µ 0 i . (2) At iteration k: Compute N independent subproblems:

u k+1 i = arg min ui J i (u i )+ < µ k , θ i (u i ) >, i = 1, ..., N .
(3) Update µ:

µ k+1 = µ k + ρ k N i=1 θ i (u k+1 i ). ( 4
) if µ k+1 -µ k < c , then stop, else k + 1 → k and go to step (2).
At step (2), the problem splits into N subproblems (decomposition procedure), which may be solved in parallel with less complexity. This interesting feature was well exploited for computation of large-scale optimization problems, for which existence of a saddle-point is ensured (otherwise convergence is not obtained).

Step (3) is the coordination procedure which consists in adjusting the Lagrange multipliers µ in order to satisfy the constraint θ(u) = 0. µ k is sent to each subproblem at step (2) in order to compute the dual functional of the problem evaluated at

µ k : w(µ k ) = min u L(u, µ k ).
µ k is used to find a consensus in the sense of constraint satisfaction. Unfortunately, introducing an augmented Lagrangian with the penalty term c 2 θ(u) 2 would induce the loss of the separability property. Of course, it would be interesting to combine the advantage of augmented Lagrangians (existence of a saddle-point) with separability (leading to parallel or decentralized computations with lower complexity). In Cohen et al. [START_REF] Cohen | Decomposition coordination methods in large scale optimization problems. The nondifferentiable case and the use of augmented Lagrangians[END_REF], a "norm-linearization technique" has been proposed to meet these two requirements. For that purpose, step (2) of the Uzawa algorithm is simply replaced by

u k+1 = arg min u J(u) + 1 2 u -u k 2 + < µ k + cθ(u k ), θ (u k ).u >, (12)
where 0 < < 1/cτ 2 , with τ being the Lipschitz constant of θ(u), and θ (u) is the Jacobian matrix of θ, while, in the same time, ρ k can be chosen constant, ρ k = ρ, with 0 < ρ < 2c.

Application to state observer design

In our case, in order to circumvent the loss of separability 2 and according to [START_REF] Akyildiz | Wireless sensor networks: A survey[END_REF], we consider the modified Lagrangian Lc (ṽ, s, μ) = 1 2

induced by c 2 si (t) - N j=1 F ij xj (t)
N i=1 T 0 { C i xi (t) -ỹi (t) 2 R -1 i + ṽi (t) 2 Q -1 i }dt + 1 2 N i=1 xi (T ) -x i 2 M -1 i + 1 2 N i=1 T 0 xi (t) -xk i (t) 2 dt + 1 2 N i=1 T 0 si (t) -sk i (t) 2 dt + N i=1 T 0 < μi (t) + c[s k i (t) - N j=1 F ij xk j (t)], si (t) - N j=1 F ij xj (t) > dt (13)
s.t. ẋi (t) = -A i xi (t) -B i ũi (t) -E i ṽi (t) -D i si (t)
where 0 < < 1/cτ 2 , with τ being the Lipschitz constant of the vector of constraints si (t) -

N j=1 F ij xj (t), i = 1, ..., N .
Lc is obviously the sum of N independant optimal LQ tracking problems, defined at each dual iteration k, and whose solutions are given by the following Lemma.

Lemma 3.1 The optimal solution to each LQ tracking problem i, i = 1, ..., N min ṽi,si,t∈I

1 2 T 0 { C i xi (t) -ỹi (t) 2 R -1 i + ṽi (t) 2 Q -1 i }dt + 1 2 xi (T ) -x i 2 M -1 i + 1 2 T 0 xi (t) -xk i (t) 2 dt + 1 2 T 0 si (t) -sk i (t) 2 dt + T 0 < μk i (t) + c[s k i (t) - N j=1 F ij xk j (t)], si (t) - N j=1 F ij xj (t) > dt (14)
s.t. ẋi (t) = -A i xi (t) -B i ũi (t) -E i ṽi (t) -D i si (t)
is given, for all t ∈ I , by

ṽi (t) = Q i E T i ( Pi (t)x i (t) + gi (t)), si (t) = s k i (t) + (D T i ( Pi (t)x i (t) + gi (t)) -μk i (t) -c[s k i (t) - N j=1 F ij xk j (t)]), (15)
where Pi (t), gi (t) and xi (t) are solutions of the canonical equations 18)

-Ṗi (t) = -A T i Pi (t) -Pi (t)A i -Pi (t)E i Q i E T i Pi (t) -Pi (t)D i D T i Pi (t) + C T i R -1 i C i + 1 I d , (16)
Pi (T ) = M -1 i , -ġi (t) = -A T i gi (t) -Pi (t)g i (t) -C T i R -1 i y i (t) -Pi (t)D i (s k i (t) -(μ k i (t) +c[s k i (t) - N j=1 F ij xk j (t)])) - N j=1 F T ji (μ k j (t) + c[s k j (t) - N i=1 F ji xk i (t)]) -Pi (t)B i ũi (t) -Pi (t)E i Q i E T i gi (t) - 1 xk i (t), (17) gi
(T ) = -M -1 i x i , ẋi (t) = -A i xi (t) -B i ũi (t) -E i ṽi (t) -D i si (t), (
xi (0) = -Pi (0) -1 gi (0). ( 19
)
Proof of Lemma 3.1. This optimal LQ tracking problem is well documented in the literature. Equations ( 16)-( 17)-( 18)-( 19) are easily obtained in closed-loop form both by application of the Pontryagin maximum principle and by considering that the adjoint state pi (t) has a separable form in both time and space pi (t) = Pi (t)x i (t) + gi (t), where Pi (t) is a symmetric matrix and gi (t) is a vector of the same dimensions as those of state xi (t) (see e.g. [START_REF] Anderson | Optimal control. Linear Quadratic Methods[END_REF]). Furthermore [START_REF] Anderson | Optimal control. Linear Quadratic Methods[END_REF] reflects the fact that each xi (0) has to minimize Lc (ṽ, s, μk ).

Theorem 3.2

The optimal trajectory solution xi (t), t ∈ I, to each LQ tracking problem ( 14), i = 1, ..., N , expressed in positive time, at each dual iteration k, is given by

Ẇi (t) = W i (t)A T i + A i W i (t) + E i Q i E T i + D i D T i -W i (t)(C T i R -1 i C i + 1 I d )W i (t), ( 20
)
W i (0) = M i . ẋi (t) = (A i -W i (t)C T i R -1 i C i )x i (t) + W i (t)C T i R -1 i y i (t) + B i u i (t) + D i s i (t) + 1 W i (t)(x k i (t) -xi (t)), + W i (t)( N j=1 F T ji (µ k j (t) + c[s k j (t) - N i=1 F ji xk i (t)])), xi (0) = x i , (21)
with

s i (t) = s k i (t) -(µ k i (t) + c[s k i (t) - N j=1 F ij xk j (t)]). ( 22
)
Proof of Theorem 3.2. Coming back to positive time and using the notations xi (t 16)- [START_REF] Alamir | Nonlinear Moving Horizon Observers: Theory & Real-Time Implementation[END_REF] given by Lemma 3.1 becomes

) = xi (T -t), P i (t) = Pi (T - t), g i (t) = gi (T -t), u i (t) = ũi (T -t), v i (t) = ṽi (T -t), s i (t) = si (T -t) and µ i (t) = μi (T -t), differential system (
Ṗi (t) = -A T i P i (t) -P i (t)A i -P i (t)E i Q i E T i P i (t) -P i (t)D i D T i P i (t) + C T i R -1 i C i + 1 I d , (23)
P i (0) = M -1 i , ġi (t) = -A T i g i (t) -P i (t)D i D T i g i (t) -C T i R -1 i y i (t) -P i (t)D i (s k i (t) -(µ k i (t) +c[s k i (t) - N j=1 F ij xk j (t)])) -P i (t)B i u i (t) - N j=1 F T ji (µ k j (t) + c[s k j (t) - N i=1 F ji xk i (t)]) -P i (t)E i Q i E T i g i (t) - 1 xk i (t), ( 24
)
g i (0) = -M -1 i x i .
If we set W i (t) = P -1 i (t), differential Riccati equation (23) becomes differential Riccati equation [START_REF] Kwakernaak | Linear optimal control systems[END_REF]. We recognize in equation [START_REF] Kwakernaak | Linear optimal control systems[END_REF] the differential Riccati equation of the optimal LQ state observer, solution of each LQ tracking problem [START_REF] Cohen | Decomposition coordination methods in large scale optimization problems. The nondifferentiable case and the use of augmented Lagrangians[END_REF], i = 1, ...N . Furthermore, the state observer equation is derived by considering estimated state xi (t) = -W i (t)g i (t), ∀t ∈ I = [0, T ], by using ( 20) and (24) (since according to [START_REF] Anderson | Optimal control. Linear Quadratic Methods[END_REF], xi (T ) = -W i (T )g i (T )), which leads to subsystem i observer ( 21), i = 1, ..., N . xi (0) = x i , since g i (0) = -M -1 i x i (that is equivalent to gi (T ) = 0 given by Lemma 3.1). Finally setting xi (t) = xi (t), and xk i (t) = xk i (t), ∀t ∈ I leads to (22).

A decentralized LQ state observer (DSO)

Finally, the decentralized optimal LQ state observer may be described as follows:

(1) Initialization (k = 0): Choose x0 i (t), s 0 i (t) and µ 0 i (t), ∀t ∈ [0, T ], i = 1, ..., N .

(2) At iteration k: For each subsystem i, i = 1, ..., N :

Integrate on the interval [0, T ], each optimal LQ state observer i, i = 1, ..., N to get estimated trajectory xk+1 i (t), by using the solution W i (t) of Riccati equation [START_REF] Kwakernaak | Linear optimal control systems[END_REF]

(computed off-line if M i is constant): ẋk+1 i (t) = (A i -W i (t)C T i R -1 i C i )x k+1 i (t) + B i u i (t) + W i (t)C T i R -1 i y i (t) + D i s k+1 i (t) + W i (t)( N j=1 F T ji (µ k j (t) + c[s k j (t) - N i=1 F ji xk i (t)])) + 1 W i (t)(x k i (t) -xk+1 i (t)), xk+1 i (0) = x i . ( 25
)
where

s k+1 i (t) = s k i (t) -(µ k i (t) + c[s k i (t) - N j=1 F ij xk j (t)]). (26)
(3) Coordination: Update the µ i 's:

µ k+1 i (t) = µ k i (t) + ρ(s k+1 i (t) - N j=1 F ij xk+1 j (t)), ( 27
)
∀t ∈ [0, T ], i = 1, ..., N .

(4) if

N i=1 µ k+1 i -µ k i 2 < 2 c /ρ 2 then stop, else k+1 → k
and go to step 2).

The following theorem can be stated:

Theorem 3.3 • Convergence: if 0 < < 1/cτ
where τ is the Lipschitz constant of the interconnection constraints and 0 < ρ < 2c, then the sequence (x k i , s k i , µ k i ), i = 1, ..., N generated by algorithm (DSO) is bounded and the related cluster point is a saddle-point of L c .

• Qualification of the solution: The obtained cluster point is solution to the decentralized optimal LQ state estimation problem [START_REF] Farina | Moving-horizon partition-based state estimation of large scale systems[END_REF], which is equivalent to a centralized finite-horizon Kalman filter for system (1)-( 2).

Proof of Theorem 3.3. The first statement is proved by using Theorem 15 in Cohen et al. [START_REF] Cohen | Decomposition coordination methods in large scale optimization problems. The nondifferentiable case and the use of augmented Lagrangians[END_REF], where the assumptions on convexity of both the cost functionals and the constraints, and the so-called Constraint Qualification Condition are satisfied in our case. The second statement is proved by verifying that the cluster point satis-fies the necessary conditions for optimality of problem (6) (due to the lack of place, this is left to the reader).

Some Remarks:

• This algorithm is well suited for decentralized state estimation through communication networks, such as sensor networks, since each of the N interconnected agents has to compute its own state estimation on the basis of information (x k j (t), µ k j (t)), ∀t ∈ [0, T ], provided by the others agents j connected to it and then has to send its computation results (x k+1 i (t), µ k+1 i (t)), ∀t ∈ [0, T ] at each iteration k, only to these connected agents j ∈ C i , where C i is the set of indices of agents connected to i. For n and N large, the computational complexity can be dramatically reduced from a n-dimensional LQ estimation problem to N n idimensional LQ estimation problem solved in parallel with additional 23), even if pairs (C i , A i ) are not observable. This assumption is often needed in other approaches (see [START_REF] Farina | Moving-horizon partition-based state estimation of large scale systems[END_REF] for instance). • The main parameter to be tuned is the regularization coefficient c. It is well known (see [START_REF] Cohen | Decomposition coordination methods in large scale optimization problems. The nondifferentiable case and the use of augmented Lagrangians[END_REF] for example) that c must be chosen not too small to ensure existence of a saddle-point in the case of non-strongly convex problems, but not too large to avoid bad numerical conditioning in solving step (2) of the pricedecomposition-coordination algorithm. • The practical implementation of this continuous-time decentralized observer may be easily performed by using a piece-wise constant discretization of the inter-

(n i + n si )-dimensional exchange data (x i (t), s i (t), µ i (t)), ∀t ∈ [0, T ], where n [1,9]si is vector s i (t) dimension. • Step (2) of
connection constraints s i (t) - N j=1
F ij xj (t) and the related Lagrange multipliers µ i (t), induced by a time discretization of the interval [0, T ], while the differential equations (including the Riccati equations) may be solved by using standard integration schemes, together with an approximation of the cost function integral by an appropriate quadrature scheme.

A MOVING HORIZON STATE OBSERVER BASED ON ALGORITHM (DSO)

To get on-line state estimation rather than batch estimation, algorithm (DSO), which provides the exact solution to a finite horizon Kalman filter, can be used as the core of a moving horizon observer approach, by con-

sidering t t-T instead of T 0
. The receding horizon observer uses only measurements on the moving horizon [t -T, t] and discards past measurements prior to the current moving horizon (see [START_REF] Alamir | Nonlinear Moving Horizon Observers: Theory & Real-Time Implementation[END_REF] for instance). Here a sampled moving horizon observer with observation horizon T = N s τ s , that updates the estimated state at time samples t l = lτ s , where τ s denotes the measurement sampling period, is defined as follows:

xi (t l ) = S i (y i t l t l-N ), i = 1, ..., N, l > N s (28)
where S i (y i t l t l-N ) denotes the estimated state obtained at time t l from iterative algorithm (DSO)2 , based on the output measurements y i (τ ), τ ∈ [t l-N , t l ], denoted as y i t l t l-N , for each subsystem i, i = 1, ..., N . At time sample t l , each observer (25) is initialized with estimated state xi (t l ) computed at time sample t l-1 (i.e. on moving interval [t l-1-N , t l-1 ]).

Convergence of the moving horizon state estimation

In this section, we will consider two particular versions of the tuning parameters T and M i , i = 1, ..., N leading to simple convergence results:

• Version 1: Deadbeat property: Following [START_REF] Kwon | A receding horizon Kalman FIR filter for linear continuous-time systems[END_REF], moving horizon observer (28) exhibits a deadbeat property provided that T > 0 and M -1 i = 0, i = 1, ..., N (assumption of infinite covariance of the initial state). The implementation of this version requires a new filter form for (25) which can be derived by using change of coordinates ν i (t) = W -1 i xi (t) (see eqn. [START_REF] Mansouri | Decentralized observers design for descriptor systems[END_REF] in [START_REF] Kwon | A receding horizon Kalman FIR filter for linear continuous-time systems[END_REF] for more details). T may be chosen arbitrarily small since system (1) is time-invariant and the initial estimated state of each observer (25) may be anything at each time sample t l .

• Version 2: Asymptotic convergence property:

Choosing T = τ s (N s = 1), x i = xi (t l-1 ), where xi (t l-1 ) denotes the local state estimate computed at previous time sample t l-1 , and

M i = W i (t l ) t l-1 , i = 1, ..., N
, where W i (t l ) t l-1 denotes the solution at time sample t l of Riccati equation ( 20) computed at previous time sample t l-1 , leads to an asymptotic convergence property. In that case, Riccati equations [START_REF] Kwakernaak | Linear optimal control systems[END_REF] have to be computed on-line, if the estimation horizon T is unknown. However the computational complexity of steps ( 2) and (3) of algorithm (DSO) is reduced to a minimum since the horizon is limited to one time sample. Since algorithm (DSO) computes the exact solution to a Kalman filter at each time sample, as l → +∞, the moving horizon observer solution tends to the solution of the infinite horizon Kalman filter, which is known for providing state estimation error dynamics stability under the assumption that the overall system is observable.

AN ILLUSTRATIVE EXAMPLE

The main objective is to build a state observer in order to estimate the water velocity distribution along an irrigation canal from distributed level measurements. In practice, using a state observer for such systems is very interesting, since it is very expensive to use velocity or flow rate sensors in addition to level sensors. It is assumed that a N + 1-sensor network has been distributed along the canal stretch [0, L], which consists in N -1 sensor-nodes equipped with only a water level sensor, while water flow velocity (or flow rate) are also available at the upstream and downstream nodes via direct measurements or hydraulic laws. The sensor-nodes have some limited computational capacities and are able to communicate with a SCADA master station using for instance an ethernet link, according to Fig. 1. A spatial discretization of Saint Venant equations governing open-channel water flow dynamics, is performed thanks to a finite-difference scheme with spatial sample dx, followed by a linearization of the discretized system around an equilibrium point of Saint Venant equations, which lead to the following observable linear system:

ḣi (t) = 1 dx (v i (t) -v i+1 (t)), i = 0, ...N -1, vi (t) = 1 dx (α i-1 h i-1 (t) + β i-1 v i-1 (t) -α i h i (t) -β i v i (t)) +δ i h i (t) + γ i v i (t), i = 1, ..., N, y i (t) = h i (t), i = 1, ..., N -1 ( 29
)
where h i (t), and v i (t) stand for the deviation values from the equilibrium of the water depth and velocities, expressed at each spatial coordinate x i . More details on this derivation are available in [START_REF] Georges | State observer design for open-channel hydraulic systems based on a sensor network[END_REF]. A natural decomposition of this system under form (3)-( 5) can be obtained for each sensor i, i = 1, ..., N , as follows:

ḣi (t) = 1 dx (v i (t) -s 1 i (t)), vi (t) = 1 dx (α i-1 s 2 i (t) + β i-1 s 3 i (t) -α i h i (t) -β i v i (t)) +δ i h i (t) + γ i v i (t), y i (t) = h i (t), i = 1, ..., N -1, s i (t) = (s 1 i (t), s 2 i (t), s 3 i (t)) T = (v i+1 (t), h i-1 (t), v i-1 (t)) T ḣ0 (t) = 1 dx (v 0 (t) -s 2 1 (t)), v0 (t) = u 1 (t), y 0 (t) = (h 0 (t), v 0 (t)) T , s 0 (t) = s 1 0 (t) = v 1 (t), ḣN (t) = u 2 (t), vN (t) = 1 dx (α N -1 s 1 N (t) + β N -1 s 2 N (t) -α N h N (t) -β N v N (t)) + δ N h N (t) + γ N v N (t), y N (t) = (h N (t), v N (t)) T , s N (t) = (s 1 N (t), s 2 N (t)) T = (h N -1 (t), v N -1 (t)) T
The interconnection variables correspond to states of the nearest neighbours of i only. Consequently, the data exchanges are limited to the nearest neighbors of each sensor. Table 1 summarizes the parameters of the test system, which corresponds to a 10 km-long canal. The decentralized observer parameters used in the simulations are given by table 2, where I n is the (n × n) identity matrix and in addition τ is the Lipschitz constant of the constraints, ρ = 1.99c and = 0.99 cτ 2 . In this cas study, version 2 of the (DSO)-based decentralized moving-horizon observer (DMHO1), (where N s = 1), and a version with a larger moving horizon N s = 4 (DMHO2) were implemented with sampling time τ s = 480s on an estimation horizon of M = 20 time samples. The simulations were performed on a desktop with Intel Core i5-2320 @ 3.00GHz and 6 Go RAM under SCILAB software 3 . Fig. 2 and 3 provide a comparison between the simulated linearized system and the two versions of the observer, with a Gaussian noise of zero mean and a standard deviation of 10 -3 added to the measured outputs. Water velocities are effectively estimated, after a 10% increase of the upstream water velocity. For both versions of the observer, the average number of dual iterations at each sample time is equal to 57.

Comparisons with some alternative methods

Using the same case study with the same covariance matrices as defined in table 2, the here-proposed approaches (DMHO1 and DMHO2) are successively compared with the Centralized Infinite-Horizon Kalman Filter (CIHKF), the Centralized Information Filter (CIF) described in [START_REF] Khan | Distributing the Kalman filter for large-scale systems[END_REF] (page 4922). O(n 3 i ) where n i is the state dimension of subsystem i, i = 1, ..., N .

Communications are also limited to connected neighbors of each subsystem. Table 3 provides some comparisons of the root-mean-square error of estimation (RM SE:

1 M M i=1
x(k) -x(k) 2 ), the average time required to compute the overall state estimation (denoted as T node ) by each sensor node (subsystem). In the case of CIF, the time needed to compute the Riccati covariance matrices over the estimation horizon (denoted as T c ) by the central location is also provided.

Table 3 Comparisons (with N = 20, τ = 480s, M = 20, and the parameters in Tables 1 and2). All the approaches exhibited similar RMSE, with a better estimation provided by CIHKF. CIF is the far less effective approach in terms of computation time needed at the central location. For both DMHO1 and DMHO2, real-time implementation is possible since each node can compute the current estimate in less than τ = 480s (the worst average computation time per node is given by 541/M = 27.05s).

CONCLUSIONS

A decentralized optimal LQ state observer design suitable for a sensor-network-based state estimation of largescale interconnected systems has been proposed. This observer design relies on both the use of an augmented Lagrangian formulation and a price-decompositioncoordination algorithm. Under mild assumptions, the algorithm converges to the solution of the centralized optimal LQ observer (finite-horizon Kalman filter). This algorithm is the core of a decentralized moving horizon state estimator. A 20 sensor-network based hydraulic canal estimation problem illustrated the effectiveness of the proposed approach, which was also compared to existing approaches. A very similar version may be derived in the framework of discrete-time linear interconnected systems. The extension to the case of nonlinearly-interconnected linear systems is also possible, with very similar communication and computation requirements. Further researches will be devoted to the analysis of the impact of communication noises or losses on the observer performance.

Fig. 1 .

 1 Fig. 1. Sensor Network Structure.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Comparison between the true water velocities and the estimated ones -with DMHO1.

  ApproachRM SE T node (in sec) Tc (in sec)

  algorithm (DSO) is equivalent to the computation of N small-scale finite horizon Kalman filters. • Coordination step (3) can be decentralized through each agent i, without needing a central coordinator agent. • Existence of a bounded and positive definite matrix P

i (t) over the interval [0, T ] is ensured by the regularization term 1 I d in the Riccati equation (

  The centralized infinitehorizon Kalman filter provides the optimal LQ state estimation in the sense of the minimal variance of the state estimation error in the presence of both state and output noises defined by the covariance matrices in table 2 or, equivalently, the solution of LQ problem (6) when T → +∞. The Centralized Information Filter is a partially distributed implementation of the finite-horizon Kalman filter: At each sensor level, the so-called local observation variables are computed and sent to a central location which computes updates of both the global state estimate and the Riccati equation covariance matrix at each time sample. When the number of sensor nodes N is large, the number of communications and the O(n 3 ) computations at the central location (where n is the state number) are some serious limitations of the approach, together with large latency and single point of failure disadvantages. In comparison, approach (DMHO) does not require a central location, and all the O(n 3 ) computations are reduced to

	N
	i=1

At each time sample t l , algorithm (DSO) can be initialized at k = 0 with any initial trajectory. However experiments show the use of the optimal solution computed at t l-1 can significantly reduce the number of dual iterations, i.e., increase the convergence speed of algorithm (DSO).

SCILAB is an open source software for numerical computation, trademark of Scilab Entreprises SAS 2013.

This work was partly supported by France-Venezuela ECOS program V11M02 / PI-2010000303. This paper was not presented at any IFAC meeting.


Thierry Miquel
email: thierry.miquel@enac.fr

A Flatness Based Speed Control Approach For Merge Behind Operations

Over the last few years, several concepts dealing the delegation to the flight crew of some tasks currently performed by the air traffic controllers have emerged. Among these new ideas, relative guidance has appeared to be capable to contribute to the enhancement of air traffic capacity though it raises difficult technical challenges. Indeed, this kind of maneuver appears difficult to perform manually, and may induce an excessive increase in flight crew workload, thus requiring new on-board automated functions. Some linear and nonlinear techniques have already been studied to design a feedback loop which performs automatically merging maneuvers and maintains station keeping behind a designated aircraft. The main contribution of this paper consists in the design of a feedback control loop using flatness control: this technique takes advantage of feedback linearization through a reference error between the ghost position of the leading aircraft and the current position of the trailing aircraft to provide a commanded speed to the autothrottle speed control system of the trailing aircraft. The robustness of the proposed design with respect to uncertainties in the dynamics of the leading aircraft and in the dynamics of the autothrottle speed control system of the trailing aircraft is achieved by the periodic update of the reference error. The reference error generation is detailed and simulation results are presented.

I. Introduction

his paper addresses a specific aspect of Air Traffic Services (ATS), namely the achievement of an orderly and expeditious flow of air traffic. Terminal control area (TMA) air traffic controllers have to merge two or more streams into a single stream by means of radar vectoring and speed instructions. The task of establishing properly spaced landing sequences is quite demanding for air traffic controllers, especially under heavy traffic conditions. An example of typical flight path for arriving aircraft at an airport is depicted hereafter: T Runway Meter fix 1 Meter fix 2

Fig. 1: Example of sequencing and merging operations for arriving aircraft at an airport

On the ground, an automation tool named Arrival Manager (AMAN) usually helps air traffic controllers to build a sequence of aircraft in order to safely and expeditiously land them [1] . Unfortunately, the airborne counterpart of the Arrival Manager, which could help the flight crew to merge its aircraft towards a meter fix according to a sequencing constraint, and ensure the station keeping after the meter fix, is not yet available. This yet to be designed onboard function will ensure the merging and station keeping tasks behind a leading aircraft designated by the air traffic controller. This new function will use the information delivered by some surveillance and communication systems and should be included in future Flight Management System (FMS) capabilities. This will enable FMS to generate a trajectory with a prescribed delay relatively to another aircraft over a specified meter fixes.

Such a new capability onboard aircraft needs on-line information on the aircraft environment, and more specifically the knowledge of the leading aircraft position and velocity. In this respect, the Automatic Dependent Surveillance-Broadcast (ADS-B) is a potential key enabler to support these surveillance requirements [2] . Aircraft equipped with ADS-B capabilities broadcast their position, velocity and identification periodically (e.g. every second). Any neighboring aircraft capable of receiving those data will therefore be able to track the surrounding traffic.

Airborne Separation Assistance System (ASAS) uses ADS-B data to provide improved airborne surveillance in support of new operational procedures that increasingly delegate tasks from air traffic controllers to pilots [4] . The goal is to achieve a higher integration between the ground and the air. By definition proposed by ICAO [4] , ASAS is an airborne system that allows the flight crew to maintain separation between their aircraft and one or more other aircraft, and provides information concerning the surrounding traffic. ASAS on-board an aircraft interacts with the existing systems (particularly the FMS) or even with the installed equipment (particularly the displays); this may also require new equipments depending on the architecture. Not only the system, but also its utilization was defined by ICAO: an ASAS application is a set of operational procedures for air traffic controllers and flight crews, which makes use of the capabilities of airborne separation assistance systems to meet a clearly-defined operational objective.

In June 2001, the FAA/EUROCONTROL published a document [3] within which 4 ASAS Application Categories were defined:

• Airborne Traffic Situational Awareness applications (ATSA)

• Airborne Spacing applications (ASPA)

• Airborne Separation applications (ASEP)

• Airborne Self-separation applications Clearly, automatic merging and station keeping is an ASAS Airborne Spacing application (ASPA) for which controllers will issue new spacing instructions like Remain, Merge behind or Radar Vector Then Merge [3] . Flight crews will be asked through those new air traffic control instructions to establish and maintain a given time or distance relative to another aircraft and will use ASAS to achieve and maintain the required spacing. On the ground, controllers will monitor and provide separation: they will intervene if necessary to ensure that the actual spacing will not decrease below the separation standard; under that aspect, an ASPA instruction is similar to the instruction clearing an aircraft to fly at a specified flight level. This type of application is clearly in the scope of the European programme SESAR [5] and the US programme NextGen. Those new technologies and procedures are part of a new paradigm that has the potential to completely transform the roles of air traffic controllers and pilots [3] .

In terms of operational benefits, this kind of application will relieve air traffic controllers of providing time consuming radar vectoring instructions to the trailing aircraft once the flight crew has accepted the relative guidance clearance. Thus, the expected benefit of such new capabilities is an increase of air traffic controller availability, which could result in increased air traffic efficiency and / or capacity. Enhancement airborne traffic situational awareness of the flight crew with associated safety benefits are also expected.

Preliminaries studies have mainly investigated the station keeping phase without taking into consideration the merging phase. This field is addressed for UAVs or military aircraft by means of linear and nonlinear techniques [6,7] . However, research for civil aircraft where safety and passenger comfort are crucial issues is still in its initial stage: in some early works, station keeping is either performed manually [8] or by means of a proportional, integral and derivative (PID) control law [9,10] .

However, operations where the merging constraint is relaxed at the meter fix have to be tackled as part of the envisioned ASPA applications and may yield in more efficient operations than applying station keeping during the merging phase. The purpose of this paper is to investigate the design of a new autopilot mode dedicated to the achievement of a specified delay / distance between aircraft at a specified meter fix. This paper considers a constant time delay as the separation objective at the meter fix; indeed it is more operational to express separation objective in terms of delay rather than distance in terminal maneuvering [11] . Eurocontrol Experimental Centre has investigated a specific method named Point Merge method for merging arrival flows in terminal area [14] . This method relies on the introduction of area navigation (RNAV, P-RNAV) and defines new route structures to revisit the merging of arrival flows. In the following, we consider that the merging and station keeping operations are achieved without change of the standard arrival routes.

American Institute of Aeronautics and Astronautics

The paper is organized as follows: in the preliminaries, reference frame and dynamics of the error to be compensated are introduced. This leads to the very simple linear state space representation. The subsequent sections present the design of the flatness based speed control which sets the commanded speed applied to the autothrottle speed control system of the trailing aircraft. Illustrative examples are then presented in order to show the effectiveness of the proposed approach.

II. Flatness Based Speed Control Design

A. Differentially flat system

Flatness was originally introduced 1995 [12] . Roughly speaking, a flat system is a square input/output system (i.e. a system having the same number of inputs and outputs) for which there exists an output vector such that all states and inputs can be expressed in terms of this output vector and its derivatives.

More precisely, a nonlinear system: ẋ= f x , u  (1) is differentially flat if one can find an output z of the form:

z= x , u , ... , u  
(2) where u (s) denotes the s-th order derivative of u with respect to time, and such that:

{ x= 0  z , ... , z s   u= 1 z , ... , z s  , z s1  (3)
Output z is called flat output. In addition, system (1) is said Lie-Bäcklund equivalent to the following system (called trivial system), where vector v is the new input, also called pseudo control:

v= z s1 (4)
Imposing a given arbitrary trajectory to z yields a trajectory for all the system variables x and u, without integrating any differential equation. Remark that the time derivatives involved in the above formulas do not imply to take derivatives of noisy signals since it involves precomputed open-loop time functions.

Fig. 2 shows the principle of a feedback loop for flatness based control: Σ denotes the system to be controlled and described by the dynamics ẋ= f  x , u , Σ -1 its inverse from which the open loop control u olc = 1  z r ,... , z r s  , z r s 1   is computed, and K a feedback gain which enables the control u to retrieve from errors between the measured state vector x and the reference trajectory zr.

Σ -1 Σ zr(t),zr (1 ) (t),… u(t) x(t) + - - + K uolc(t)

B. Reference frame and error dynamics

The reference frame is affixed to the trailing aircraft. From Fig. 3, the following notations are adopted:

• Subscript L stands for ghost position of the leading aircraft, that is the position of the leading aircraft delayed by the required delay which should be achieved at the meter fix P. This delay is set be the air traffic controller; • Subscript stands F for the current position of trailing aircraft (follower);

• dL stands for the distance between the ghost position of the leading aircraft and the meter fix P, whereas dF stands for the current distance between the trailing aircraft and the meter fix.

American Institute of Aeronautics and Astronautics

3

C.

V F dF d L V L χ L χ F P L F 0 X Y Meter fix

Fig. 3: Reference frame and notations

The dynamics of the error e between the ghost position of the leading aircraft and the current position of the trailing aircraft is simply given by: ė=V L -V F (5) In the following, the difference VL -VF will be used as a pseudo control, denoted u, to control the dynamics of the error e:

u=V L -V F ⇒ ė=u (6) The linear model which represents the dynamics of the error e between the ghost position of the leading aircraft and the current position of the trailing is then very simple: it is basically an integrator.

D. Feedback control using flat outputs

It is clear that the error e between the ghost position of the leading aircraft and the current position of the trailing aircraft is a flat output of system (6).

It is assumed in the following that an autopilot function dealing with airspeed control is available on-board the trailing aircraft. The purpose of the relative guidance control system is then to set the commanded airspeed such that the trailing aircraft reaches the meter fix at the same time than the ghost position of the leading aircraft; then station keeping has to be maintained behind the leading aircraft.

To achieve global stability, a Lyapunov based design is used: assume that there exists a scalar function L of the state e, with continuous first order derivatives such that [13] :

• L e is positive definite • L e∞ as ∥e∥ ∞ •
L e is negative definite The following quadratic candidate Lyapunov function will be used, where er stands for the reference error between the ghost position of the leading aircraft and the current position of the trailing aircraft:

L e= 1 2  e-e r  2 (7)
The pseudo control u has to be designed such that the reference error er tracked by system (6). Taking time derivative and using (6) yields:

L e= e-e r   ė-ėr = e-e r  uėr  (8) L e will be negative definite by choosing the pseudo control u as follows, where kp is a gain to be set by the designer:

u= ėr -k p e-e r ⇒ L e=-k p  e-e r  2 (9)
The commanded airspeed applied to the autopilot function dealing with airspeed control is then given by:

u≝V L -V F ⇒ V F =V L -u=V L -ėr  k p e-e r  (10)
The error e(t) between the positions of ghost aircraft and the trailing aircraft results from the integration of ( 5) American Institute of Aeronautics and Astronautics but can be directly obtained without any integration thanks to ADS-B data.

When compared with Fig. 2, the term V Lėr represents the open loop control uolc(t), whereas kp is the feedback gain K which enables the control to retrieve from errors between the measured error e and the reference error er.

The problem of the relative guidance is then viewed as a problem where explicit trajectory generation of reference er is required. Since the behavior of flat systems is determined by their flat outputs, the use of flatness to design relative guidance controller is quite natural. Indeed, trajectory in the flat output space can be planned and then mapped to appropriate inputs.

E. Reference error computation

As far as the flight plan of the leading aircraft is not known by the trailing aircraft, the duration T needed by the ghost position of the leading aircraft to reach the meter fix is simply estimated by the ratio between distance dL to the meter fix and the airspeed VL of the leading aircraft ghost position:

T ≈ d L 0 V L (11)
Referring to Fig. 3, and assuming a constant airspeed VL for the leading aircraft, the reference error er(t) between the positions of ghost aircraft and the trailing aircraft is chosen as follows, where l(t) denotes the reference distance of the trailing aircraft to the meter fix:

{ e r t≝ d F 0-l t -d L 0-V L t  ėr ≝V L -V r where { l 0 ≝0 l t ≝V t e0 =d F 0-d L 0  (12)
Then, the control law (10) reduces to the following expression, where dF(t) denotes the actual distance between the trailing aircraft and the merging point and VF the commanded airspeed applied to the autothrottle speed control system of the trailing aircraft:

{ e t≝d F t - d L 0-V L t  e r t ≝ d F 0-l t-d L 0-V L t  ėr ≝V L -V r ⇒ V F =V L -ėr k p e-e r =V r k p d F t-d F 0l t  (13)
As far as it is desired that trailing aircraft arrives at the meter fix at the same time than the ghost position of the leading aircraft, and with the same speed, the reference er between the ghost position and the trailing aircraft position has to be computed such that the following constraints are satisfied:

{ e r T =0 ėr T =V L -V r T ≝0 (14)
In the following, we will replace time t by the dimensionless parameter τ which is defined by the ratio between actual time t and the estimated duration T of the merging maneuver (see ( 11)):

0≤≝ t T ≤1 (15)
Two options to set the reference length l(τ) of the trailing aircraft to the meter fix are presented hereafter.

• The first option only consider the constraints presented in ( 14):

l =T  a 0  a 2 b atanb-1atanb  (16)
Where parameters a0, a2 and b are free parameters. Time derivation of the preceding equation leads to the expression of the trailing aircraft's reference airspeed:

V r = a 0  a 2 b-1 2 1 (17)
In order to satisfy the constraints presented in ( 14), parameters a0, a2 and b shall satisfy the following relationships:

American Institute of Aeronautics and Astronautics { e r 1≝0 ⇔l 1=d F 0⇔ a 0  atan b b

a 2 = d F 0 T ėr 1≝V L -V r 1=0⇔ a 0  a 2 =V L (18)
Assuming that parameter b is already set, parameters a0 and a2 can easily be computed from the preceding equations.

• The second option consider the constraints presented in ( 14) and also sets the initial value of the reference speed at the current speed of the trailing aircraft. We consider the following expression for reference length l(τ), where parameters a0, a1, a2 and b are free parameters:

l =T  a 0  a 1 b atan b  a 2 b atanb-1atanb  (19)
Time derivation of the preceding equation leads to the expression of the trailing aircraft's reference airspeed:

V r = a 0  a 1 b 2 1  a 2 b-1 1 (20)
In order to set the initial value of the reference speed at the current speed of the trailing aircraft and to satisfy the constraints presented in ( 14), parameters a0, a1, a2 and b shall satisfy the following relationships:

{ V r 0 ≝V 0 ⇔a 0 a 1  a 2 b1 =V 0 e r 1≝0 ⇔ l 1=d F 0 ⇔a 0  atanb b a 1 a 2 = d F 0 T ėr 1≝V L -V r 1=0 ⇔ a 0  a 1 b1 a 2 =V L (21)
Assuming that parameter b is already set, parameters a0, a1 and a2 can easily be computed from the preceding equations.

III. Illustrative examples

A. Scenarios

In this section, two scenarios are designed in order to illustrate the properties of the control laws previously designed. The meter fix point P is situated at 0 NM.

For both scenarios, the ghost position of the leading aircraft starts at dL(0) = -25 NM with initial airspeed VL(0) = 220 kts. It is supposed to broadcast its data every second (ADS-B basic assumption). The initial position of the trailing aircraft is dF(0) = -30 NM with an initial airspeed VF(0) = 210 kts.

In the first scenario, the ghost airspeed of the leading aircraft remains constant, whereas in the second scenario the ghost airspeed of the leading aircraft decreases towards 120 kts with a deceleration of 0.01g.

While the ghost position of the leading aircraft has not reached the meter fix point P, the speed control feedback loop of the trailing aircraft is in the merge behind mode: the reference errors e r t  and ėr t  are computed according to the two options previously presented, and the commanded airspeed VF which feeds the autopilot is given by (10):

V F =V L -u=V Lėr  k p e-e r  (22) Once the ghost position of the leading aircraft has passed the meter fix point P, the speed control feedback loop switches to the remain behind mode where the reference errors e r t  and ėr t  are set to zero. In such mode, the commanded airspeed VF is simply:

V F =V L  k p e (23) In both modes the gain kp has been set to 50. In addition the robustness of the proposed design with respect to uncertainties in the dynamics of the leading aircraft and the dynamics of the speed control autopilot is achieved by the periodic update of the reference error: every 30 sec for both options.

The dynamics of the speed control autopilot of the trailing aircraft is modeled by a second order linear system with a damping ratio of 0.7 and a natural frequency w0 of 0.5 rad/sec. Thus the relationship between then actual speed V of the trailing aircraft and the commanded airspeed VF is given by: American Institute of Aeronautics and Astronautics

Fig. 2 :

 2 Fig. 2: Feedback loop for flatness based control

Fig. 4 :Fig. 5 :

 45 Fig. 4: Speeds of trailing aircraft (commanded in red and actual in green) and ghost (blue) Fig. 5 presents the positions of ghost and trailing aircraft for all scenarios:

It is worth noticing that the acceleration of the trailing aircraft has been limited to 0.05g.

Positions xL of the ghost and xF of the trailing aircraft are linked to speed VL and V through an integration:

The two flatness based speed control design to compute the commanded speed VF are compared against a design where the commanded speed VF is simply computed as proportional to the actual error between the ghost position of the leading aircraft and the current position of the trailing aircraft, that is as in ( 23) even during the merge behind phase.

B. Simulation results

All the scenarios are summarized and labeled in the following [15] . Fig. 4 presents commanded true airspeed and actual airspeed of the trailing aircraft as well as the airspeed of the ghost for all scenarios whereas Fig. 5 presents the positions of ghost and trailing aircraft for all scenarios. Speeds are expressed in kts, positions in NM and time in seconds.

As expected, the scenarios involving a flatness based speed control enable the trailing aircraft to reach the meter fix at the same time than the ghost position of the leading aircraft, whereas the scenario involving a proportional based speed control leads the trailing aircraft to catch the ghost position before the meter fix. In terms of airspeed, the scenarios involving a flatness based speed control leads to lower speed deviation and speed command when compared to the scenario involving a proportional based speed control. Thus the greater complexity of the flatness based speed control due the fact that time explicitly appears in the control law and to the need of periodically update the computation of the reference error are compensated by higher performances for merge achievement.

American Institute of Aeronautics and Astronautics

IV. Conclusion

In this paper, the design of a new autopilot mode dedicated to the achievement of a specified delay between commercial aircraft at a specified meter fix has been considered.

This envisioned new capability on-board commercial aircraft takes advantage of the Automatic Dependent Surveillance-Broadcast (ADS-B) and provides new perspectives to potentially increase air traffic control efficiency. It could be the airborne counterpart of the ground based arrival manager.

The proposed approach is based on flatness control: this technique takes advantage of feedback linearization through a reference error between the ghost position of the leading aircraft and the current position of the trailing aircraft. It provides a commanded speed to the speed control autopilot of the trailing aircraft. The robustness of the proposed design with respect to uncertainties in the dynamics of the leading aircraft and in the dynamics of the speed control autopilot of the trailing aircraft is achieved by the periodic update of the reference error. When compared to a proportional based control, the greater complexity of the flatness based control due the fact that time explicitly appears in the control law and to the need of periodically update the computation of the reference error are compensated by higher performances for merge operations.

This paper has considered a 1 dimensional movement. The extension to the 3 dimensional space can easily be done by considering the curvilinear abscissa flown by aircraft.

Simulation results illustrates the efficiency of the proposed design when compared to a proportional based speed control depending only of the actual error. Nevertheless, the way to accurately get the time at which the meter fix is reach by the leading aircraft should be tackled. Furthermore, additional studies in terms of operational scenarios and taking into account the vertical motion of the aircraft as well as the effect of wind are needed in order to refine and validate the proposed design.

V. References


Maxime Galloux
email: maxime.galloux@gmail.comgondran

Alexandre Gondran

Cyril Allignol
email: allignol@recherche.enac.fr

Nicolas Barnier
email: barnier@recherche.enac.fr

Conflict resolution by speed regulation

Keywords: Air Traffic Management, Speed Regulation, Conflict Avoidance

teaching and research institutions in France or abroad, or from public or private research centers.

I. INTRODUCTION

For more than forty years, global air traffic has never ceased to increase. The current traffic control systems are reaching their structural limits, so that the traffic growth might reduce the safety level of the airspace. Thus, new methods and concepts are to be set up in order to adapt to future traffic, as advocated by research programmes such as SESAR 1 [START_REF]SESAR Concept of Operations[END_REF] in Europe or NextGen in the US. Therefore, in this paper, we propose a model for speed regulation in order to avoid a maximum number of conflicts. This work is built upon the conflict model presented in [START_REF] Barnier | 4D-Trajectory Deconfliction Through Departure Time Adjustment. 8 th USA/Europe Air Traffic Management Research and Development Seminar ATM 2009[END_REF] and uses small speed regulations like in [START_REF] Rey | Technical Report on the Minimization of Potential Air Conflicts using Speed Control[END_REF] to avoid conflicts.

A standard day of traffic within the european airspace is made of approximately 30,000 flights. Each flight follows a flight plan described as a sequence of waypoints that the aircraft will fly over. The aim of this project is to avoid the air space conflicts in advance by the use of slight modifications of flight speeds, or equivalently by constraining fly over times on waypoints (RTA 2 ), in order to decrease the controllers' workload. Our model assumes that a conflict between two flights that follow intersecting trajectories can only happen on a common waypoint of their routescatching-up flights are handled by considering multiple waypoints. After a first potential conflict detection step, resulting in a combinatorial optimization problem that generalizes the one presented in [START_REF] Barnier | 4D-Trajectory Deconfliction Through Departure Time Adjustment. 8 th USA/Europe Air Traffic Management Research and Development Seminar ATM 2009[END_REF], the RTAs are computed by a centralized system so as to avoid as many conflicts as possible and 1 Single European Sky ATM Research 2 Requested Time of Arrival to optimize operating costs. Then they are transmitted to aircraft during the flight to modify their flight plans. To take uncertainties into account, the problem would be iteratively solved over a rolling horizon (typically 20-30 min) as presented in [START_REF] Rey | Technical Report on the Minimization of Potential Air Conflicts using Speed Control[END_REF].

This combinatorial optimization problem could be solved by state-of-the-art techniques such as Constraint Programming as in [START_REF] Barnier | 4D-Trajectory Deconfliction Through Departure Time Adjustment. 8 th USA/Europe Air Traffic Management Research and Development Seminar ATM 2009[END_REF] or Local Search as in [START_REF] Allignol | Optimized Flight Level Allocation at the Continental Scale[END_REF].

II. SEPARATION CONSTRAINTS

As in [START_REF] Rey | Technical Report on the Minimization of Potential Air Conflicts using Speed Control[END_REF], this model is designed to take into account the functionalities of future Flight Management Systems (FMS) in the SESAR context, which will be able to dynamically accommodate several RTAs on the waypoints of their trajectories with an accuracy of a few seconds. Its output is therefore a set of RTAs for each flight involved in a potential conflict within the time window of the resolution, trying to minimize the number of actual conflicts and their durations.

In this context, the trajectory of a flight i is represented by a sequence of 3D points and associated times corresponding to the waypoints of its route indicated by the flight plan:

ω k i , θ k i , k ∈ [1, n i ]
where the θ k i i,k are the decision variables and n i the number of waypoints of flight i. Furthermore, we note:

• σ k i , the curvilinear abscissa (or oriented length along the trajectory) of each waypoint k of flight i;

• v k i (for all k ∈ [1, n i -1]
), the speed, considered constant, between waypoints ω k i and ω k+1 i .

DISCRETIZATION OF TRAJECTORIES

Following the conflict model presented in [START_REF] Barnier | 4D-Trajectory Deconfliction Through Departure Time Adjustment. 8 th USA/Europe Air Traffic Management Research and Development Seminar ATM 2009[END_REF], potential conflicts between two flights i and j are detected by pairwise checks on the points of a discretization of each trajectory, with a grain fine enough to ensure that even the shortest potential conflicts will be taken into account, as described in [START_REF] Barnier | Trajectory deconfliction with constraint programming[END_REF] for example. We note:

(p k ′ i , t k ′ i ), k ′ ∈ [1, m i ]
the discretization of the trajectory of flight i consisting of a sequence of m i 3D points p k ′ i and associated times t k ′ i , as illustrated in 1.

Similarly to the waypoints, we note s k ′ i the curvilinear abscissa of each point p k ′ i . These abscissae can be easily

ω l+1 j ω k-1 i ω k+1 i ω k i = ω l j ω l-1 j flight i p k ′ i p l ′ j p l ′ +1 j 5 N M flight j
Figure 1: Projection in the horizontal plane of the trajectories of flights i and j, which are in potential conflict at their common waypoint

ω k i = ω l j . Point p k ′ i is in potential conflict with points p l ′ j and p l ′ +1 j .
computed by adding the distances between the previous points of the trajectory:

s k ′ i = k ′ -1 l ′ =1 dist(p l ′ i , p l ′ +1 i )
providing that the turning points of the trajectory, i.e. the waypoints ω k i , are included in the discretization points p k ′ i . A similar relation would be true for the abscissae σ k i of the waypoints (as they indeed are the turning points), but only during cruise, as the waypoints are too distant from each other to approximate the vertical profile of the descent or climb phase precisely enough.

CONFLICT DETECTION

Two points p k ′ i and p l ′ j of flight i and j are in potential conflict if their horizontal and vertical distances both violates the separation norm (usually 5 NM horizontally and 1000 ft vertically), as depicted in figure 1 for the horizontal plane. For two such points, there would be an actual conflict between flights i and j if they are located at these points at the same time. To avoid a conflict, it is therefore necessary that:

∀k ′ ∈ [1, m i ] , ∀l ′ ∈ [1, m j ] , dist H (p k ′ i , p l ′ j ) < 5 N M ∧ dist V (p k ′ i , p l ′ j ) < 1000 f t =⇒ t k ′ i = t l ′ j (1)
where dist H is the distance in the horizontal plane and dist V in the vertical plane.

For intersecting trajectories, we assume, as in [START_REF] Rey | Technical Report on the Minimization of Potential Air Conflicts using Speed Control[END_REF], that successive potentially conflicting points are all located in the vicinity of the same waypoint. As the speed is considered constant between two consecutive waypoints, it is possible to translate these conflict inequations 1 to their closest waypoint, such that the conflict constraints can be expressed as a temporal separation at this waypoint.

The resolution will therefore consist in regulating the speed of these flights (through the issuing of RTAs at the waypoints only) so as to avoid conflicts, i.e. ensure that the corresponding inequation θ k i = θ l j holds for all the neighbouring pairs of potentially conflicting points p k ′ i and p l ′ j . However, catch-up conflicts along the same route portion (or if the angle of intersecting trajectories is very low) cannot be considered local to a single waypoint. Thus, all potentially conflicting points located between waypoints ω k-1 i and ω k+1 i are reported to a potential conflict associated with ω k i . Several such conflicts will then be defined at each successive waypoints as long as the flights follow the same route.

We can now define a generalized (intersecting or catchingup) potential conflict between two flights i and j as the set of conflicting pairs of trajectory points around a single waypoint:

C kl ij = {(k ′ , l ′ ) ∈ [1, m i ] × [1, m j ] , s.t. p k ′ i and p l ′ j are in potential conflict near ω k i = ω l j , and σ k-1 i < s k ′ i < σ k+1 i , σ l-1 j < s l ′ j < σ l+1 j
To compute the resulting contraint between θ k i and θ l j , we first need to express t k ′ i as a function of θ k i and the speed of the aircraft, which is v k-1 i before waypoint ω k i and v k i afterwards. Therefore:

∀k ′ ∈ [1, m i ] , ∃k ∈ [1, n i ] , t k ′ i =      θ k i + s k ′ i -σ k i v k i si s k ′ i ≥ σ k i θ k i + s k ′ i -σ k i v k-1 i si s k ′ i < σ k i (2)
Inequation 1 and equation 2 can then be combined, with four different cases depending on the locations of points p k ′ i and p l ′ j with respect to waypoint ω k i = ω l j . Note that if flights i and j have several distinct (non continuous) potential conflicts, there will be as many (non empty) sets C kl ij of conflicting points:

∀(k, l) ∈ [1, n i ] × [1, n j ] s.t. C kl ij = ∅, ∀(k ′ , l ′ ) ∈ C kl ij , t k ′ i = t l ′ j ⇔ θ k i -θ l j =                    s l ′ j -σ l j v l j - s k ′ i -σ k i v k i if s l ′ j ≥ σ l j and s k ′ i ≥ σ k i s l ′ j -σ l j v l j - s k ′ i -σ k i v k-1 i if s l ′ j ≥ σ l j and s k ′ i < σ k i s l ′ j -σ l j v l-1 j - s k ′ i -σ k i v k i if s l ′ j < σ l j and s k ′ i ≥ σ k i s l ′ j -σ l j v l-1 j - s k ′ i -σ k i v k-1 i if s l ′ j < σ l j and s k ′ i < σ k i
(3) However, these expressions depend on the (unknown) variable speeds of aircraft, whereas the conflict model presented in [START_REF] Barnier | 4D-Trajectory Deconfliction Through Departure Time Adjustment. 8 th USA/Europe Air Traffic Management Research and Development Seminar ATM 2009[END_REF] uses static bounds for the time difference between flights i and j at waypoint ω k i = ω l j . In the next section, we explain how the speed variations are tightly bounded in our operational context, which allows us to approximate the values of equation 3 by small intervals.

III. SMALL SPEED ADJUSTMENT AND CONFLICT APPROXIMATION

Following [START_REF] Rey | Technical Report on the Minimization of Potential Air Conflicts using Speed Control[END_REF], we consider the same kind of small speed variation as in the ERASMUS project [START_REF] Villiers | Automatisation du contrôle de la circulation aérienne "ERASMUS" -Une voie conviviale pour franchir le "mur de la capacité[END_REF]. Two ratio paramaters α ≤ 1 and α ≥ 1 are introduced to bound the possible speed adjustment, such that if we note v 0 i the reference speed of flight i, all other speed variables are restricted to take values in a small interval α v 0 i , α v 0 i :

∀i ∈ [1, n], ∀k ∈ [1, n i -1] , v k i ∈ α v 0 i , α v 0 i ( 4
)
where n is the number of flights in the instance. These ratio paramaters are typically chosen in the range α = 0.97 and α = 1.06 to limit the cost of regulation on the fuel consumption of aircraft and to have the smallest impact possible on standard Air Traffic Control (ATC) practices.

As the speed variations considered are small and bounded, we can bound the values of θ k iθ l j described in equation 3 by constants r k ′ l ′ ij (lower bound) and r k ′ l ′ ij (upper bound) defined as follows:

∀(k, l) ∈ [1, n i ] × [1, n j ] s.t. C kl ij = ∅, ∀(k ′ , l ′ ) ∈ C kl ij , r k ′ l ′ ij =                    s l ′ j -σ l j α v 0 j - s k ′ i -σ k i α v 0 i if s l ′ j ≥ σ l j and s k ′ i ≥ σ k i s l ′ j -σ l j α v 0 j + σ k i -s k ′ i α v 0 i if s l ′ j ≥ σ l j and s k ′ i < σ k i - σ l j -s l ′ j α v 0 j - s k ′ i -σ k i α v 0 i if s l ′ j < σ l j and s k ′ i ≥ σ k i - σ l j -s l ′ j α v 0 j + σ k i -s k ′ i α v 0 i if s l ′ j < σ l j and s k ′ i < σ k i r k ′ l ′ ij =                    s l ′ j -σ l j α v 0 j - s k ′ i -σ k i α v 0 i if s l ′ j ≥ σ l j and s k ′ i ≥ σ k i s l ′ j -σ l j α v 0 j + σ k i -s k ′ i α v 0 i if s l ′ j ≥ σ l j and s k ′ i < σ k i - σ l j -s l ′ j α v 0 j - s k ′ i -σ k i α v 0 i if s l ′ j < σ l j and s k ′ i ≥ σ k i - σ l j -s l ′ j α v 0 j + σ k i -s k ′ i α v 0 i if s l ′ j < σ l j and s k ′ i < σ k i
(5) For one conflict set C kl ij , bounds of a forbidden interval for the difference of times θ k iθ l j at the waypoint can then be computed by cumulating the inequations over the pairs of conflicting points (k ′ , l ′ ) ∈ C kl ij :

r kl ij = min (k ′ ,l ′ )∈C kl ij r k ′ l ′ ij r kl ij = max (k ′ ,l ′ )∈C kl ij r k ′ l ′ ij
Moreover, the bounds on the speed of each aircraft further constrain the θ k i variables for two consecutive waypoints

ω k i and ω k+1 i : ∀k ∈ [1, n i -1] , θ k+1 i -θ k i ∈ d k i α v 0 i , d k i α v 0 i (6)
with

d k i = σ k+1 i -σ k i
, the distance between ω k i and ω k+1 i . In the following, these bounds on the travel time of flight i between two consecutive waypoints are respectively noted

T k i = d k i αv 0 i and T k i = d k i αv 0 i .

IV. MODEL

The conflict detection processing of the previous sections can now be sum up with a more standard (and concise) combinatorial decision problem formulation, with variables θ k i representing the RTA of flight i on waypoint ω k i , for a set of n potentially conflicting flights over the time window considered:

Find: ∀i ∈ [1, n] , (θ k i ) k∈[1,ni] s.t.: ∀(i, j) ∈ [1, n] 2 , i < j, ∀(k, l) ∈ [1, n i ] × [1, n j ] s.t. C kl ij = ∅, θ k i -θ l j / ∈ r kl ij , r kl ij ∀i, ∀k ∈ [1, m i -1] , θ k+1 i -θ k i ∈ T k i , T k i
The first set of constraints are derived from the conflicts between two potentially conflicting flights i and j, r kl ij and r kl ij being the bounds of the forbidden values for the difference of fly times over the conflicting waypoint ω k i = ω l j . The second one characterizes the speed limitation for each aircraft as stated in equation 6.

OPTIMIZATION

Among the admissible solutions of this decision problem, the ones that minimize airlines costs should be preferred. To optimize their operating costs, airlines generally tune the Cost Index (CI) parameter (used by the FMS to optimize the flight parameters along its trajectory) for each flight, which represents the relative importance of the cost of fuel with respect to the cost of flight time, i.e.:

CI = cost time cost fuel
where cost time is in $ per time unit, and cost fuel in $ per mass unit. For a given flight i, we can therefore consider that the airline cost is the sum of the extra cost of time and of the extra cost of fuel multiplied by CI i (the Cost Index of flight i), compared to the reference trajectory at the reference speed:

cost i = CI i × cost fuel i + T 0 i -θ ni i -θ 1 i where T 0 i = ni-1 k=1 d k v 0 i
is the total flight time of flight i at its reference speed v 0 i . If we assume that v 0 i is the optimal speed of the flight and that a discrepancy from v 0 i leads to a proportional increase of fuel consumption during the time flown at this speed, then:

cost fuel i = ni-1 k=1 v k i -v 0 i v 0 i θ k+1 i -θ k i = ni-1 k=1 d k i v 0 i -θ k+1 i -θ k i
So the cost for one flight can be expressed as:

cost i = CI i ni-1 k=1 d k i v 0 i -θ k+1 i -θ k i + T 0 i -θ ni i -θ 1 i
and the total cost of one solution simply is the sum of the costs of all flights:

cost = n i=1 cost i (7)
Note that other parameters could be taken into account, like the number of speed changes, or could use a more realistic function to estimate the effect of the speed discrepancy from the reference speed, provided we could gather enough data from airlines and manufacturers.

INFEASIBILITY

If some conflict constraints cannot be satisfied, it may still be of interest to relax these constraints and search for solutions with as few violated constraints as possible. Furthermore, when a conflict constraint is violated, the longer the conflict lasts, the more hazardous is the operational situation, so it is interesting as well to search for solutions that minimize the duration of all actual conflicts, as proposed in [START_REF] Rey | Technical Report on the Minimization of Potential Air Conflicts using Speed Control[END_REF].

As the conflict constraints of our model enforce that θ k iθ l j does not belong to r kl ij , r kl ij , we have to compute the "distance" of the time difference from the center of the interval M kl ij = 1 2 r kl ij + r kl ij :

dist kl ij = θ k i -θ l j -M kl ij
The longer the conflict lasts, the closer θ k i -θ l j is to M kl ij , and the smaller is dist kl ij . So to obtain a measure of the duration of the conflict, the distance is subtracted from the half of the length of the interval r kl ij , r kl ij , i.e. from L kl ij = r kl ijr kl ij . Moreover, the cost should be 0 outside the conflict interval:

max(0, 1 2 L kl ij -(θ k i -θ l j ) -M kl ij )
The total cost over all remaining conflicts is the sum of all the actual conflicts durations, for all pairs of flights in potential conflict (i.e. ∀i < j s.t. ∃(k, l), C kl ij = ∅):

cost conflict = i<j k,l s.t. C kl ij =∅ max(0, 1 2 L kl ij -|(θ k i -θ l j ) -M kl ij |)
This cost could be combined with the cost defined in equation 7 with suitable weighting to balance the cost of remaining conflicts with the operational one, so as to take both criteria into account for overconstrained instances.

V. CONCLUSIONS AND PERSPECTIVES

To overcome the traffic growth predicted by current ATM research programs in Europe, we propose a novel deconfliction model based on small speed regulations, following the mixed integer program presented in [START_REF] Rey | Technical Report on the Minimization of Potential Air Conflicts using Speed Control[END_REF] with a conflict model that generalizes the work presented in [START_REF] Barnier | 4D-Trajectory Deconfliction Through Departure Time Adjustment. 8 th USA/Europe Air Traffic Management Research and Development Seminar ATM 2009[END_REF]. The associated centralized solver would output new RTAs to dynamically adjust the flight plans during the flight, taking operational costs for airlines and for ATC into account. The resolution would be iteratively performed over a rolling horizon (20-30 min) to handle the uncertainties inherent to trajectory prediction.

The described model is currently being implemented using Constraint Programming and Local Search as optimization techniques. Simulations will be carried out with Europe-wide traffic data.