Antoine Drouin
email: antoine.drouin@enac.fr

Jules G Slama
email: julesslama@yahoo.com.br

Félix A C Mora-Camino

F Mora

Hybrid Backstepping Control for Rotorcraft Guidance

The purpose of this communication is to display a non-linear control approach based on backstepping for the positioning and orientation for a four-rotor aircraft. Realistic rotorcraft flight dynamics are introduced and the effectiveness of its control channels is analyzed. Then two complementary implementations of the backstepping control approach are considered. The compatibility of these control approaches with a two-layer control structure devoted to the guidance of the rotorcraft is displayed. The resulting control laws are detailed and their expected performances are discussed. A simulation study is performed where the performances of the proposed control structure are compared with those of a classical non linear inverse control solution.

I. INTRODUCTION

N the last years a large interest has risen for the fourrotor concept since it appears to present simultaneously hovering, orientation and trajectory tracking capabilities of interest for many practical applications [START_REF] Hoffmann | The Standford Tetsbed of Autonomous Rotorcraft for Multi-Agent Control[END_REF]. The flight mechanics of this rotorcraft are highly non-linear and different control approaches [START_REF] Khalil | Nonlinear Systems[END_REF], [START_REF] Sutton | Reinforcement learning: an introduction[END_REF], have been considered with little success to achieve either only autonomous hovering and orientation or also trajectory tracking. In this paper we consider the flight dynamics of a fourrotor aircraft with fixed pitch blades. The control problem of interest is the design of flight control laws enabling autonomous positioning and orientation for this class of rotorcraft. This study investigates the solution of this problem using a backstepping control approach. Here it is required that a single continuous control law performs the whole maneuver while to achieve it by manual control it appears necessary to go through a succession of elementary maneuvers since the system is underactuated. The backstepping control approach, which has also been applied to airships flight control law design [START_REF] Singh | Nonlinear decoupled control synthesis for maneuvring aircraft[END_REF], is introduced. Two different design techniques are developed following the main guidelines of this approach. These two implementations appear of direct interest for the design of a new two-layer control structure based on backstepping control laws. The proposed control approach is compared with a nonlinear inverse control approach introduced in [START_REF] Drouin | Nonlinear Control Structures for Rotorcraft Positionning[END_REF].

II. ROTORCRAFT FLIGHT DYNAMICS

The considered system is shown in figure 1 where rotors one and three are clockwise while rotors two and four are counter clockwise. The main simplifying assumptions adopted with respect to flight dynamics in this study are a rigid cross structure, constant wind, negligible aerodynamic contributions resulting from translational speed, no ground effect as well as small air density effects and negligible response times for the rotors. It is then possible to write the rotorcraft flight equations as follows [START_REF] Etkin | Dynamics of Flight-Stability and Control[END_REF].

A. Rotorcraft Flight Equations

The rotor forces and moments are given by:

2 i i f F ω = (1-1) { 4 , 3 , 2 , 1 ∈ i } } 2 i i i f k F k M ω = = i (1-2) { 4 , 3 , 2 , 1 ∈
Where f and k are positive constants and ω i is the rotational speed of rotor i. These speeds and forces satisfy the constraints:

max 0 ω ω ≤ ≤ i i { 4 , 3 , 2 , 1 ∈ i } (2-1) { } 4 , 3 , 2 , 1 0 max max ∈ = ≤ ≤ i f F F i i i ω (2-2)
Since the inertia matrix of the rotorcraft can be considered diagonal with I xx = I yy , the roll, pitch and yaw moment equations may be written as:

xx I r q k F F l p / ) ) ( ( 2 2 4 + - = & (3-1) yy I r p k F F l q / ) ) ( ( 4 3 1 + - = & (3-2) I zz I F F F F k r / )) ( ( 3 4 1 2 - + - = & (3-3)
Where p, q and r are the roll, pitch and yaw body angular rates. Here and , where

I ) ( 2 yy zz I I k - = ) ( 4 zz xx I I k - =
xx , I yy and I zz are the inertia moments in body-axis, and l is the length of the four arms of the rotorcraft.

Let φ, θ and ψ be respectively the bank, pitch and heading angles, then the Euler equations relating the derivatives of the attitude angles to the body angular rates, are given by: ) cos )(sin ( r q tg p φ φ θ φ

+ + = & (4-1) r q φ φ θ sin cos - = & (4-2) θ φ φ ψ cos / ) cos (sin r q + = & (4-
y y d F m a - - = φ ψ φ θ ψ (5-2) ) ) cos( ) )(cos( / 1 ( z z d F m g a + - = φ θ (5-3)
where x, y and z are the centre of gravity coordinates, m is the total mass of the rotorcraft and:

4 3 2 1 F F F F F + + + = (6)
Here the drag force

( )' z y x d d d d =
is given by: ( )' where ρ is the volumetric mass of the air, S is the reference surface for the aerodynamics of the rotorcraft and C d is its dimensionless drag factor. The components of the rotorcraft airspeed are:

) ( ) ( ) (
, and

x a w x u - = & y a w y v - = & z a w z w - = & (8)
The equilibrium conditions (hovering) in an horizontal wind (w z =0) with a given heading ψ are such that:

, 0 , 0 , 0 = = = r q p 0 , 0 , 0 = = = z y x & & & (9-1) )) sin )(
φ θ = = = = (9-4)

B. Analysis of Rotorcraft Flight Dynamics

Here we are interested in controlling the four-rotor aircraft so that its centre of gravity reaches and stays hovering at a predefined position while its heading acquires and maintains a given orientation. Many potential applications require this capability to be available in UAVs' while this problem can be also considered as a first step towards the design of more efficient trajectory tracking systems. The manoeuvre under study is, when performed manually through direct radio control of the four engine thrusts (see picture 1), quite difficult to be achieved in one step. Experimentally it appears that no direct approach is feasible and that much depends on the rotorcraft attitude angles φ and θ specially when considering the control of its horizontal position error (x-x c , y-y c ). Equations (5-1) and show that to get any horizontal acceleration, it is necessary to have a non zero attitude ( 0 ≠ φ or 0 ≠ θ ), they show also that the orientation of the acceleration is dependent of the heading angle ψ. Equations (3-3) with (4-3) and show that given the attitude angles φ and θ, it is easy to master the heading angle error (ψ-ψ c ) and the vertical position error (z-z c ). From equations (3-1), (3-2) and (3-3), it appears that the effectiveness of the rotor actuators is much larger with respect to the roll and pitch axis than with respect to the yaw axis. Then we consider that attitude piloting is involved with controlling the angles θ and φ. In equations (3-1) and (3-2), the effect of the rotor forces appears as differences so, we define new attitude inputs u q and u p as:

3 1 F F u q - = (10.1) 2 4 F F u p - =
In the heading and position dynamics, the effects of rotor forces and moments appear as sums, so we define new guidance inputs u ψ and u z as:

) ( ) ( 3 1 4 2 F F F F u + - + = ψ 4 3 2 1 F F F F F u z + + + = = (10.2) ' ] [ ' ] [ 4 3 2 1 z q p u u u u u F F F F F ψ = = (11-3)
Equations (3-1), (3-2) and (3-3) are rewritten:

(12-1)

xx p I r q k u l p / ) ( 2 + = & yy q I r p k u l q / ) ( 4 + = & (12-2) zz I u k r / ψ = & (12-3)
It appears that u q and u p can be made to vary significantly with u ψ and u z remaining constant. Attitude angles φ and θ can be seen as virtual controls for the horizontal position of the rotorcraft. Here the attitude dynamics are considered to be the fast dynamics , they are at the heart of the control system. The heading and height dynamics are intermediate while the dynamics of the horizontal position coordinates are the slower. This can lead to multilevel closed-loop control structures.

III. BACKSTEPPING CONTROL

A. The Backstepping Control Approach

The backstepping technique is a rather recent non-linear control technique, which applies to cascaded systems. The main idea is to use intermediate state variables as virtual inputs to take advantage of the causality relations displayed by the cascaded state representation. The convergence of the output variables towards their target values is obtained by the construction, step by step, of an auxiliary Lyapunov function. This general idea can be developed in different ways, as it will be shown in the next sub-section. The main interest of the backstepping approach is that the stability of the controlled system as well as the convergence of the outputs towards their reference values can be guaranteed without inducing, like in the case of the nonlinear control approach, the decoupling of the outputs dynamics. Indeed, it can be considered that the decoupling of the outputs dynamics demands an additional effort from the control channels with then a higher possibility of saturation for the actuators, either in position or speed, resulting in downgraded performances. Finally, another advantage of this approach is that several matrices of parameters are introduced while constructing the control law, providing a large variety of possibilities to shape conveniently the outputs dynamics as well as the control signals.

B. Direct Implementation of the Backstepping

Consider a cascaded system whose state representation is given by:

2 1 x x = & and ) , ( 2 2 U x g x = & (13) where n R x ∈ 1 , n R x ∈ 2 are state variables and n R u ∈
is the control input and g is a smooth diffeomorphism with respect to u. The control objective here is to design a control law such that the state 1

x can be stabilized at c x 1 .

Here also, 2

x can be regarded as a virtual control input for the dynamics of 1

x while the dynamics of 2 x are controlled by the real control input U . Now, suppose that there exists a control law ) , (

1 1 2 c x x G x =
such that the dynamics of 1

x can be stabilized at c x 1 while we can find a Lyapunov function ) (

1 1 1 c x x V - , which satisfies the condition: ) ( ) , ( )' / ( ) ( 1 1 1 1 1 1 1 1 1 c c c x x W x x G x V x x V - - ≤ ∂ ∂ = - & (14) where ) ( 1 1 c x x W - is a positive definite function of x 1 . A possible choice is: ) ( ) , ( 1 1 1 1 c c x x x x G - Λ - = (15)
where Λ is a positive definite symmetric matrix. Then in this case:

) ( )' ( 2 1 ) , ( ) ( 1 1 1 1 1 1 1 1 1 c c c c x x x x x x W x x V - - = = - (16)
The whole dynamics can be expressed as:

z x x G x c + = ) , ( 1 1 1 & and w z = & (17-1)
where ) , (

1 1 2 c x x G x z - = (17-2) and 2 1 2 ) / ( ) , ( x x G u x g w ∂ ∂ - = (17-3)
Then a candidate Lyapunov function of the full system is given by:

z z x x V z x x V c c ' 2 / 1 ) ( ) , ( 1 1 1 1 1 + - = - (18) ) , ( 1 1 z x x V c
is given by: The time derivative of (19)

w z z x x G x V z x x V c c ' ) ) , ( ( ' ) / ( ) , ( 1 1 1 1 1 1 + + ∂ ∂ = - & then: w z z x V x x W z x x V c c ' )' / ( ) ( ) , ( 1 1 1 1 1 1 + ∂ ∂ + - - ≤ - & (20)
and by an adequate choice of w , such as :

z x V w Ω - ∂ ∂ - ) / ( 1 1 = (21)
where Ω is a symmetric positive definite matrix, the full system is globally asymptotically stable since it satisfies the following condition:

(22) z z x x W z x x V c c Ω - - - ≤ - ' ) ( ) , ( 1 1 1 1

&

Finally, the effective control input is given by:

( ) ) , ( ) / (( ) ( 1 1 2 1 1 2 1 c x x G x x V x g U - Ω + ∂ ∂ - = - ) (23)

C. Indirect Implementation of Backstepping

Now we consider the case where the cascaded system cannot be written easily in the form (13), but it obeys to:

) , ( U X g X = & with () X h Y = (24)
where X∈R n , U∈R m , Y∈R m , g is a smooth vector field of X and U and h is a smooth vector field of X. The system has, with respect to each independent output Y , a relative degree r i i

( , i = 1, … , m) around the state n r m i i ≤ + ∑ =1 ) 1 (
X if the output dynamics can be written as:

0 ) , ( ) ( ) 1 ( ) 1 ( 1 1 U X B X A Y Y m r m r + = ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ + + M (25)
Here we assume that , where the jacobian of B with respect to the control inputs is invertible.

1 2 1 = = = = m r r r L
In that case, two auxiliary outputs can be defined:

Y Y Y L Z c & + - = ) ( 1 and Y Z & = 2 (26)
where L is a positive definite symmetric matrix. A candidate Lyapunov function is then given by:

) ( 2 1 2 ' 2 1 ' 1 2 Z Z Z Z V + = (27)
The time derivative of V 2 is such as:

2 2 1 1 2 ' ' Z Z Z Z V + = & (28-1) or Y Y L Y Y L Y Y Y L V ' ) ( ' ) 2 ) ( ( 2 - + + - = (28-2)
Choosing a control such as :

) 2 ( Y Y L Y Y L & & & & + Λ - = + (29)
where Λ is another symmetric positive definite matrix. We have: (34) which is the case when φ and θ remains small with respect to ±π/2. We introduce now the two R 4x4 symmetric positive definite matrices L and Λ and adopt the control law (31

Y L Y Y Y Y L Y Y Y L V c c ψ θ φ = and B(X, u) is such as: u Y J u X B ) ( ) , ( = (33-1) with ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢
) with ' )] / ( [ ) ( m d g g P P X A z q p - = ψ (35)
to compute the current input vector u .

B. Horizontal backstepping control layer

The state representation of the horizontal dynamics is given by: (36-1)

x v x = & y v y = & (36-2) ) )) sin( ) sin( ) cos( ) sin( ) )((cos( / 1 ( x x d F m v - + = φ ψ φ θ ψ & (36-3) ) )) sin( ) cos( ) cos( ) sin( ) )((sin( / 1 ( y
y d F m v - - = φ ψ φ θ ψ & (36-4)
where ψ and F are defined by the inner control loop. This state representation corresponds to the one studied in the case of (13-1) and (13-2). Then following the corresponding backstepping approach, we get with V 1 chosen according to relation (23) the following reference values for φ and θ :

) / ) ( cos ) ( n arcsin((si z y y x x c u d m d m + - + = ε ψ ε ψ φ (37-1) ) cos / / ) ( sin / ) ( s arcsin((co c z y y z x x c u d m u d m φ ε ψ ε ψ θ + + + = (37-2)
where ε x and ε y are given by:

⎥ ⎦ ⎤ ⎢ ⎣ ⎡ Ω - ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ - - Λ Ω + - = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ y x y y x x I c c y x & & ) ( 2 ε ε (37-3)
where Λ and Ω are symmetric positive definite matrices. Then, the horizontal position of the rotorcraft follows the linear dynamics:

⎥ ⎦ ⎤ ⎢ ⎣ ⎡ - - Λ Ω + + ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ Ω + ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ c c y y x x I y x y x ) ( 2 0 ) ( ) 2 ( ) ( ) ( 1 1 = - + Λ + Λ + Λ - - c Y Y Y L L Y L & & & (39)
When matrices Λ and L are diagonal, these dynamics are decoupled and the poles of the decoupled dynamics are the roots of the m different characteristic polynomials: and where s is the Laplace variable. In this case, since the λ i and the μ i are positive real, we get always real negative roots. In the general case, the dynamics modes of the outputs will be characterized by the solutions of the global characteristic polynomial:

0 ) ) 2 ( det( = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ Λ + + Λ - L sI L I sI m m m (41)
Since this last relation is independent of the application it is possible to study once for all the reachable pole sets within the left half complex plane. Remark: In the case of an horizontal wind, a necessary condition for final convergence and equilibrium, independent of the choice of matrices Λ and L or even of the control approach is given by: where ψ c is the heading reference value.

V. SIMULATION RESULTS

The selected gains for the backstepping control law are displayed on Table 1 while the selected dynamics for the attitude, the altitude, the heading and the horizontal position are second order linear dynamics characterized by their respective damping coefficients and natural frequencies. These values are reported on Table 2. The produced figures display different time responses of the rotorcraft under either the backstepping control law or a reference non-linear inverse control law. The comparison is performed in two stages: first, the responses of the rotorcraft to a step in attitude (either φ or θ) for each control law, are evaluated. Then, for each control law , the responses of the rotorcraft to a step in position (x, y or z) are evaluated and compared. The evaluation of the inner attitude control loop is important since the guidance capability of the rotorcraft, an under actuated device, is directly dependent of the controllability of its attitude angles.

K θ = -2 K q =-0.23 K φ = -2 K p = -0.23 K ψ = -0.02 K r = -0.025 K z = 0.12 K z & =0.15 K x = 0.137 K x &=0.183 K y = 0.137 K = 0.183 y & TABLE II SELECTED DYNAMIC PARAMETERS θ = 0.8 ω θ = 10 rad/s ζ φ = 0.8 ω φ = 10 rad/s ζ ψ = 0.8 ω ψ = 2 rad/s ζ z = 0.8 ω z = 1.5 rad/s ζ x = 0.8 ω x = 1.5 rad/s ζ y = 0.8 ω y = 1.5 rad/s
The results (see figures 2 to 5) show that the two control laws, in both levels, present equivalent performances. However, while the non linear inverse control law produces second order linear dynamics for the attitude angles and the position and heading outputs, the backstepping control law produces clearly a non linear behaviour for these variables.

In particular (figure 4) , since the final convergences of the backstepping control law is rather slow, the non linear inverse solution can produce, for a same response time, a less input demanding solution. However, as shown in the following figures, other parameters settings may lead to responses where the backstepping approach is slightly superior. Other simulation studies should be performed in particular to show clearly the advantage of using advanced non-linear control law instead of empirical-intuitive ones. Also, the realisation of simulation studies should be of interest to explore the impact of actuator saturations on the flight domain and feasible manoeuvres of the rotorcraft.

VI. CONCLUSION In this communication the applicability of a non-linear control approach to the positioning and orientation of a rotorcraft has been treated. Since this system is highly nonlinear, naturally unstable and rather under-actuated, the design of a unique control law to perform safely the whole manoeuvre is not straightforward and a multilevel control approach must be considered. So a multilevel control structure has been introduced. It appeared that the direct application of the backstepping control approaches was not desirable and that it was more judicious to realize two different implementations of the backstepping guidelines to insure first the internal stability and then guidance of the controlled system.

REFERENCES FIGURES

Fig. 1 .

 1 Fig. 1. Four rotor aircraft

Fig. 3 .Fig. 2 Fig. 4 .Fig. 5 .

 3245 Fig.3. Step response of φ (ω φ =15 rad/s) x x c

& & & & & &

& & & & & & &

& & & & &

& & & & & = 0 (38)Since the actuator settings are determined by the inner control loop, let us have a look at the corresponding dynamics. The outputs φ, θ, ψ and z of the inner closed loop follows the dynamics given by:

and it is straightforward to show that the system is globally asymptotically stable. The corresponding control law is given by :

)

IV. APPLICATION OF BACKSTEPPING TO ROTORCRAFT FLIGHT CONTROL Analyzing relations (3-i), (4-i) and (5-i), i =1 to 3, it appears that their equations can be separated into two sets: one relative to the slower dynamics, the horizontal dynamics, and corresponding to the first case considered in the previous section and one relative to other dynamics and corresponding to the second case. The above approaches of backstepping are now applied to each of the control layer necessary to perform attitude control and guidance of the rotorcraft.

A. Control of Rotorcraft Attitude and Level

The attitude and altitude dynamics can be given by the state equations:

(cos sin ) /(cos cos

where the exact expressions of


Cyril Douthe
email: cyril.douthe@ifsttar.fr

Olivier Baverel
email: baverel@enpc.fr

Morphological and mechanical investigation of double layer reciprocal structures

Keywords: Spatial structures, structural morphology, form-finding, dynamic relaxation, linear elastic behaviour

From a technological point of view, reciprocal structures have the advantage of simplifying considerably connections in the sense that members are connected only by pairs. This constructional advantage has two main consequences: technological constrains are replaced by geometrical constrains and bending moments are increased through the non convergence of members to nodes with multiple connections. Geometrical constrains can be dealt by using form-finding methods such as dynamic relaxation algorithm. However resisting bending moments to gain stiffness is difficult without increasing the weight of the structure. For standard reticulated structures, common strategies consist in introducing curvature in the structure and/or modifying the structure into a double layer space structure. The proposed paper is thus an attempt to apply these strategies to reciprocal structures and to develop spherical domes with structural thickness. Several configurations will be investigated and compared in term of geometrical feasibility and structural performance.

INTRODUCTION

From a technological point of view, reciprocal structures have the advantage of simplifying considerably connections in the sense that members are connected only by pairs. This constructional advantage has two main consequences: technological constrains are replaced by geometrical constrains and bending moments are increased through the non convergence of members at nodes with multiple connections. Indeed it has been shown in [START_REF] Douthe | Design of nexorades or reciprocal frame systems with the dynamic relaxation method[END_REF] that the reciprocal structures are softer and generate higher stresses than their conventional counterpart. In order to increase the performances of reciprocal structures, it is proposed here to investigate two basic strategies used for reticulated systems: to introduce curvature and/or to increase the structural thickness. Various configurations of double layer systems will thus be shown in this paper. In the following section, simple configurations where structural thickness is obtained by the elements themselves (Vierendeel trusses or any other space trusses) will be examined. A simple construction made of ladder elements will be shown. Then, in the third section, a solution for building a double layer reciprocal system with a flat configuration is developed based on analytical calculations [START_REF] Sénéchal | Analytical investigations on elementary nexorades[END_REF]. In the fourth section, a double layer dome configuration is investigated. The form-finding method using the dynamic relaxation algorithm [START_REF] Douthe | Design of nexorades or reciprocal frame systems with the dynamic relaxation method[END_REF] is detailed and a comparison of the structural behaviour with a simple layer configuration is made. Finally the authors conclude on the potential of double layer reciprocal systems and illustrate their feasibility on a prototype.

DOUBLE LAYER SYSTEMS WITH TRUSS ELEMENTS

A possibility for creating double layer reciprocal systems is to use truss elements in standard single layer configuration similar to those shown in [3,[START_REF] Baverel | Configuration processing of nexorades using genetic algorithms[END_REF]. In the illustrations shown from figure 1 to 3, the truss is in fact a Vierendeel truss. A simple fan is shown in figure 1. In this figure , one can see that the higher and the lower parts of each truss are laying one the top of the other. As a consequence, the trusses do not stand vertically but slightly tilted. This tilt depends on the height of the beam and on the diameter of the elements used for connections (which is here zero as the ladder are simply laying on each other reciprocally). It is important to note that the configuration of figure 2 and 3 has only a morphological purpose and that it does not attempt to demonstrate the mechanical behaviour of such structures which would need further investigations. Nevertheless, the structure is more than 10 m span and stands it self-weight, which already prove its potential for practical application.

Moreover, it was demonstrated in [START_REF] Douthe | Design of nexorades or reciprocal frame systems with the dynamic relaxation method[END_REF] and [START_REF] Baverel | Nexorades: a family of interwoven space structures[END_REF] that when loaded, the fan have a slight rotation in torsion that generates a moment in the weak axis of the members. Further investigation of the mechanical behaviour of such structures should also consider using members made of spatial trusses instead of simple Vierendel trusses.

DOUBLE LAYER SYSTEMS WITH SLENDER MEMBERS FOR FLAT CONFIGURATIONS

In the preceding section, the method proposed to obtain a "double layer" structure was to start from a single layer reciprocal structure and to replace slender members by space trusses. In the second method presented here, the double layer reciprocal system is obtained by transformation of a standard double layer spatial structure in a similar way to the transformations studied in [START_REF] Sénéchal | Analytical investigations on elementary nexorades[END_REF][START_REF] Baverel | Configuration processing of nexorades using genetic algorithms[END_REF][START_REF] Baverel | Nexorades Based on Regular Polyhedra[END_REF]. There are several way to compute practically the final configuration of the reciprocal system after the transformation: some analytical methods have been developed by rotation of members in [START_REF] Sénéchal | Analytical investigations on elementary nexorades[END_REF] and [START_REF] Baverel | Nexorades Based on Regular Polyhedra[END_REF] and also some numerical methods using genetic algorithm [START_REF] Baverel | Configuration processing of nexorades using genetic algorithms[END_REF] or dynamic relaxation [START_REF] Douthe | Design of nexorades or reciprocal frame systems with the dynamic relaxation method[END_REF]. This section will focus on the analytical method while the next section will make use of the dynamic relaxation algorithm.

The analytical developments presented in the following concentrate on one specific configuration shown in figure 4 [START_REF] Bangash | Elements of Spatial Structures: Analysis and Design[END_REF] which consists in a two way space grid with square pyramids. When transformed into a reciprocal system, the members are rotated so that multiple connections disappear and that members are only connected by pairs. Considering the symmetries of the spatial structure, all the nodes and all the fans are identical (excepting of course fans laying on the perimeter of the spatial frame). The configuration of a typical fan after transformation is shown in figure 5. The green line represents the members that are on the upper layer or lower layer and the red members represent the diagonals element that connect both layers in figure 4. The parameters of the model are: the length and engagement length of the members of the upper (or lower) layer, those of the diagonal members and the eccentricity of the connection (or the distance between the axis of two connected members). There are thus 5 design parameters. The strategy for determining analytically the values of these parameters is the same as in [START_REF] Sénéchal | Analytical investigations on elementary nexorades[END_REF]:

-write that the direction of the eccentricity between two members is perpendicular to each member, -write that the distance between two members is equal to the eccentricity, -write that the distance between two contact points on a member is equal to the engagement length.

-write geometric compatibility between upper and lower fand (which is equivalent to adjusting the diagonal length so that upper and lower fans can effectively be connected).

To ease calculations (which are detailed in [START_REF] Sénéchal | Nexorades: geometrical analyses and constructions[END_REF]), it is supposed that the relative engagement length are equals and that the eccentricity is constant. Under those conditions, it is found that, out-of five parameters, only two are independent: the length of the upper (or lower) members and the eccentricity. Neglecting homothetic configurations, the geometry of the basic cell of the reciprocal system shown in figure 6 is thus entirely determined by the ratio of the eccentricity and the length of a member. Especially, the height of the double layer system which is a key parameter of the mechanical behaviour is given by that ratio: the higher the ratio eccentricity over length, the higher the ratio height over length.

Introducing some practical data by taking a member length of 100 cm with an eccentricity of 4 cm, it is found that the necessary engagement length is close to 3.6 cm and the height of the structure of 107 cm (see figure 6). Extending this configuration over several fans, the reciprocal system shown in figure 7 is obtained. It is higher and denser than the initial configuration of figure 4. As it is, the reciprocal structure is thus relatively difficult to build (engagement length and eccentricity being of the same order of magnitude) and also less efficient from a mechanical point of view as it uses more material for the diagonals which contribute very few to the overall stiffness. The only solution for improving the mechanical behaviour of this double layer reciprocal system is thus to give up the hypothesis of equal relative engagement length of members and diagonals, so that the height becomes again an independent design parameter. This however add one more parameter in a complex non-linear analytical problem which is still under investigation.

DOUBLE LAYER SYSTEMS WITH SLENDER MEMBERS FOR CURVED CONFIGURATIONS

Form-finding of double layer grids

The configuration investigated in this section is inspired by a scheme already illustrated by R. Buckminster Fuller for spatial structures that he had called "magnesium-framed geodesic dome" [7]. The idea is to connect two superimposed triangular grids with two different lengths (the larger one been √3 times longer than the other) as shown in figure 8. The nodes of the larger grid coincide with nodes of the smaller grid, so that a total number of twelve members converge at these nodes which form the upper layer of the spatial systems. The nodes of the smaller grid that do not coincide with nodes of the larger grid are connected directly with each others through an hexagonal mesh with nodes to which six members converge. These nodes form the lower layer of the spatial systems. The members connecting the lower hexagonal mesh to the upper triangular mesh insure the transfer of shear forces between the two layers. These three types of member (upper layer, lower layer and diagonals) are visible on figure 8. The objective here is thus to transform this flat grid in which multiple members converge at nodes into a double layer reciprocal system. In previous section, the system was extremely symmetric and the irregularities of the system at boundaries were neglected, so that an analytic solution for the transformation could be found. This time, symmetries are more complex and practical solutions for the interruption of the grid at boundary are desired. The transformation into a reciprocal system will thus be done numerically with a form-finding method based on the dynamic relaxation algorithm which had been described by the authors in [START_REF] Douthe | Design of nexorades or reciprocal frame systems with the dynamic relaxation method[END_REF]. The principle of this method relies on the introduction of a fictitious mechanical problem where members initially converging at nodes are bent to form a reciprocal system. The solution of the problem is a structure in equilibrium where the bending energy of all members has been dissipated and where members are straight. The fictitious mechanical properties are defined according to the designer needs as will be illustrated in next paragraph. Then, as in the final configuration, members are straight, there is no need for catching finely the curvature. Convergence speed is thus increased by using the least necessary number of nodes to describe the members geometry, which means only four nodes per member (one node at each connection). Convergence is also improved by the use of a routine which has been added to the method presented in [START_REF] Douthe | Design of nexorades or reciprocal frame systems with the dynamic relaxation method[END_REF] and which consists in increasing gradually the bending stiffness of members and updating their reference length as they straighten.

The geometrical parameters of the form-finding consist in a set of two configurations: a reference configuration which is stress free and used to evaluate at each time step the stresses in the members and an initial configuration from which calculations are started. As by definition the members are straight, the reference configuration reduces to a set of reference lengths of the members. The reference length for the members of the hexagonal grid is thus set to 100 cm with an engagement length of 10 cm. That of the members of the larger triangular grid is consequently of 173.2 cm while their engagement lengths are the same and taken to 10 cm. Then, as the intermediate members connecting the upper and lower grid will be used as spacers between the two layers, their length is fixed to 122.1 cm (100 cm in projection on the grids planes and 70 cm vertically), so that after the form-finding, the structural height should be close to 70 cm.

At the external nodes located on the boundary of the grid, the regularity of the reciprocal system is necessarily interrupted: the number of members converging at nodes is different from nodes in the current part of the grid. To keep a certain regularity for these external fans, some adjustments have been made by modifying locally the engagement length of a member or by adding members with suited lengths (see the locations of the modifications on figure 8). An example of each adjustment is illustrated in figure 9a and9b. Figure 9a shows an external fan to which nine members are converging instead of the expected twelve and where it is necessary to increase the engagement length of one member to close the truncated dodecagon. Figure 9b shows then an external fan in which a new member has been added. Initially five members (three members from the upper layer and two interconnecting members) converge toward this external node which would have a dodecagonal form if located in the current part of the grid. A sixth member with the suited length is thus added to close the open polygon respecting the geometry of the dodecagon. This construction principle is applied to all incomplete hexagons and dodecagons of the grid of figure 8, so that in total 18 additional members are introduced and 12 engagement lengths modified. For the initial configuration, the choice of lengths is free provided that geometrical constraints are satisfied, that boundary conditions and connectivity of nodes are respected. The initial configuration is thus chosen to guide the formfinding toward a solution that suits the designer needs. Here, the initial configuration is defined by projecting the grid of figure 8 on a sphere with a 5 m radius and by bending the members extremities to satisfy the reciprocal system connectivity with a constant eccentricity for all connexions of 4 cm. (see figure 10). To achieve this transformation of the initial grid into a reciprocal structure by bending, it is necessary to define how the members initially converging at the same node will rely on each other in the reciprocal structure. The transformation of each node is characterised by two parameters that had been called in [3] the disposition and the style of the fan. Here the dispositions of fans are chosen so that, after transformation, fans of the lower layer form downward fans while fans of the upper layer form upward fans (see zooms in figure 10). Both layer can hence be easily connected. Then, it had been seen in [START_REF] Douthe | Design of nexorades or reciprocal frame systems with the dynamic relaxation method[END_REF] that, from a mechanical point of view, alternating fans styles increases significantly the rotations of the fans and the bending moments in the members. All fans are thus defined with the same style: rightward.

Fig. 10: Initial configuration with bent members

Finally the fictitious mechanical properties of the members are set, so that their bending stiffness is much larger than their axial stiffness, so that the lengths of members are susceptible of finite variations. They can hence adapt to the various geometrical constraints and still stay straight (the out-of straightness after form-finding should be less than 10 -3 degree).

Given all these parameters, the result of the form-finding or the final configuration is shown in figure 11. The two grids (the lower hexagonal and the upper triangular) lay on two different pseudo-spheres with a radius around 6.0 m. The distance between the layers is of approximately 70 cm as expected. During the form-finding, the lengths of the members have slightly changed. The members of the larger grid (forming the upper layer) have an average length of 183.0 cm with an engagement length around 10.5 cm. Those of the hexagonal grid have a length of 103.5 cm with an engagement length of 9.5 cm. Concerning the interconnecting members, their length is also around 123 cm but their engagement lengths vary depending on the layer of the extremity: around 6.0 cm in the upper layer and 17.5 cm in the lower layer. These lengths are slightly different at the extremities of the grid where the regularity of the scheme is only partly conserved.

Fig. 11: Final configuration of the double layer reciprocal system

Considering the number of parameters entering the form-finding procedure (namely the geometry of members, the disposition and style of the fans, the initial configuration and the boundary conditions), the reciprocal double layer structure shown in figure 11 is only one among an infinity of possible structures. At this stage where only geometrical issues have been discussed, the criterion for choosing a suitable structure is constructability. Practical experiences from previous realisations of regular polyhedra [START_REF] Sénéchal | Analytical investigations on elementary nexorades[END_REF][START_REF] Baverel | Nexorades Based on Regular Polyhedra[END_REF] have proved that engagement lengths larger than 1.5 times the eccentricity between members axes are sufficient to insure feasibility and placing of connectors. As this is here the case, the reciprocal structure seems suitable. Further optimisation should thus be made relying on the mechanical behaviour which will be investigated in section 4c.

Form-finding of a single layer equivalent grid

The equivalent single layer grid is obtained by transformation into a reciprocal structure of the smaller triangular grid in figure 8 with no additional members on the boundary. The geometrical parameters of the form-finding are thus the following: an identical reference length of 100 cm with an engagement length of 20 cm (twice higher than previously because the number of members converging at a node is twice lower) and an eccentricity of 4 cm. All fans have upward disposition and rightward style. The initial configuration of the grid is such that the nodes lay on a spherical dome with a radius of 5 m. It must be here specified that this time the position of all external nodes is fixed. Like previously the fictitious mechanical properties of the members are set, so that their bending stiffness is much larger than their axial stiffness; their lengths are thus adaptable. The result of the form-finding is shown in figure 12.

In the final configuration of figure 12, the reciprocal structure has a height of 1.94 m for a 7 m span which corresponds to a spherical dome with a radius of 4.1 m. Like for the double layer structure, during the form-finding, the lengths of the members have changed and vary between 143 cm at crown to 116.5 cm close to the supports. The engagement lengths vary between 34.5 cm at crown to 21.5 cm close to the supports. This dome structure is indubitably constructible and due to its double curved shape should prove good mechanical behaviour despite a relatively large engagement length.

Mechanical behaviour of the two configurations

During the form-finding step, fictitious mechanical properties had been used to find the geometry of the previous reciprocal systems. In the geometries shown in figure 10 or 11, the members are straight but stretched (elongated) with a prestress that is not physical. Before conducting any mechanical analyses under external load, it is thus necessary to relieve the systems from inner axial stresses. This is done by setting every members lengths at rest (their reference lengths) equal to their lengths in the final configuration obtained by form finding. The removing of these stresses in the single layer or double layer structure has no consequence on the geometry of the structure (or, more precisely, the changes in the form are smaller than the out-of straightness tolerated in the form-finding).

Then, for the study of the mechanical behaviour, realistic mechanical and geometrical properties are introduced. Members are thus taken as steel circular hollow sections with an external diameter of 40 mm and a thickness of 3 mm. Boundary conditions consist in blocking the positions of nodes located on the periphery of both structures. For the single layer dome, the boundary conditions are thus identical to those used for the form-finding. For the double layer structure, it is not easy to define precisely external nodes considering the number of members touching the periphery, so that it was chosen to support the structure on the additional members shown in figure 9 whose 6 degrees of freedom have been fixed.

Concerning external loads, two cases have been investigated: one symmetric uniformly distributed load and one nonsymmetrical distributed load one half of the structure. To ease comparison between systems, these distributed load have been introduced via vertical point loads applied at each end of every member of the smaller grid (this means every member of the single layer system and only the shorter members of the double layer system). The intensity of each point load is 0.5 kN which corresponds approximately to 3.5 kN/m 2 (which is an acceptable density for a roof structure).

From a numerical point of view, non-linear analyses are conducted using the same algorithm as for the form-finding (reproducing hence the methodology used in [START_REF] Douthe | Design of nexorades or reciprocal frame systems with the dynamic relaxation method[END_REF]). Convergence is reached when the initial kinetic energy of the structure has been dissipated (practically when the current kinetic energy peak is bellow 1/100 th of the initial/maximum peak). It must here be noted, that the poor meshing of the form-finding has been slightly improved for the accuracy of results and that the density of the mesh has been refined to reach six nodes per member: one for each end connection, one for each engagement connection and then two intermediate nodes. Mechanical behaviour is assessed through: the overall stiffness of the structure (measured through the average and maximal displacement of the loaded points and the rotation of the fans), the intensity of reactions and inner stresses in the members (essentially axial forces and bending moments).

≈ 0.0 ≈ 0.1 ≈ 0.0 0.1 Single layer 0.2 0.5 ≈ 0.0 0.5
Table 1 shows the displacements under the two loadings. It is first noticed that under both loadings the displacements are low and would satisfy typical slope criterion for roofing structures: 1/200 th of the span is here equal to 35 mm. It is also observed that, under dissymmetric loading, the maximum displacements are of the same magnitude than under symmetric loading which denote very sane mechanical behaviour. Comparing both structures, it is remarked that the displacements in the single layer structure are higher than those of the double layer structure which assesses that the system has effectively a significant structural height. Moreover, it is remarked that fan rotations (which are caused by the presence of eccentric normal forces converging at the fan) have very low magnitudes, especially in the double layer system. This difference between both systems should be related to the fact that the single layer structure works mainly as a membrane structure with forces that flow within the surface, whereas the double layer structure works rather as a thick plate where forces flow within the thickness of the plate. That being said, the order of magnitude of the rotations in the single layer dome are comparable to those observed in the dome of [START_REF] Douthe | Design of nexorades or reciprocal frame systems with the dynamic relaxation method[END_REF] (taking into account the fact that here the Young modulus is twice higher, one finds 1.2 degree in [START_REF] Douthe | Design of nexorades or reciprocal frame systems with the dynamic relaxation method[END_REF] and 0.5 here) which confirms that the influence of the fan rotations on the general behaviour can be considered as a second order effect. Table 2 shows then extreme values of the efforts in both structures. As expected, it is observed that the horizontal thrust relatively to the vertical reaction is higher in the dome structure and that the inclination of the total reaction is approximately that the last member. This confirms that the single layer system carries load mainly within its surface, that its behaviour is dominated by membrane actions. These membrane actions correspond to normal forces in the members of the reciprocal systems and, as these forces are not converging toward the same points, they create high bending stresses. Their values however remains such that the construction of the structure seems feasible by using a steel grade of S355 or by slightly increasing the inertia of the members.

The distributions of the normal forces under symmetric and dissymmetric loading can be seen respectively in figure 13a and 13b (in which the left half of the structure is the loaded one). In both loading cases, the magnitude of forces are comparable. Under symmetric loading, all members are in compression except a kind of tension ring which overtakes hoop forces above supports as would be the case in a spherical shell. Under dissymmetric loading, compression dominates the loaded half of the structure while a more complex distribution of tension and compression is observed in the other half.

Concerning the bending moments, their distributions are shown in figure 14a and14b. In both loading cases, it is remarked that they develop mainly in the loaded areas and that their magnitude decreases as the member are closer to the support. However, it is noticed that, despite a lower total load under dissymmetric loading, the highest value of the bending moment is 15% higher than under symmetric loading. Once again, this is a phenomenon that can be observed in continuous spherical shell in which dissymmetric loadings often generate higher stress than symmetric ones.

The general behaviour of the double layer system is quite different as it is globally dominated by a thick shell behaviour and bending effects. The loads which act perpendicular to the surface are resisted through compression in the upper layer (which consists in a triangular grid) and tension in the lower layer with some compression as it consists in an hexagonal grid. This is illustrated in figure 15 where hot colours correspond to compression. Moreover, table 2 shows that compression forces are approximately of the same magnitude in the double layer structure as in the single layer dome but that tension forces are higher. This is related to the fact that, this time, tension forces do not correspond to hoop action but to shear forces within the shell which are logically higher at the supports. An another illustration of the fact that the reciprocal system works effectively like a double layer thick shell can be found when comparing the magnitude of forces and bending moments under symmetric and dissymmetric loading. They are almost equal, the double layer system resists in the same way to both loadings which typical from structures that carry loads by bending.

CONCLUSIONS

Different configurations of double layers reciprocal systems where investigated in this paper. The first configuration showed the potential of using spatial members to create reciprocal configurations with structural thickness. The second configuration then illustrated how planar double layer reciprocal systems can be obtained by transformation of standard spatial grids. The results of analytical developments on a two ways square grid have been presented. However the too restrictive hypotheses made for the analytical calculations of the geometry leaded to engagement length that are quite small compared to the diameter of the members so that make the construction difficult. Further development with more design parameters are thus currently under progress.

In the last chapter, two spatial configurations of reciprocal systems have been investigated. The first one is a double layer dome based on the superposition of two triangular grids, the second one a single layer dome based on a triangular grid. A form-finding method for the definition of geometries compatible with reciprocal arrangement of members has been illustrated. Then a mechanical study of the two structures has been conducted. On the one hand, the double layer reciprocal structure works effectively like a spatial structure, like a thick plate in bending with high bending stiffness and low sensitivity to dissymmetric loading. On the other hand, the single layer dome works like a thin shell dominated by membrane forces and is therefore more sensible to dissymmetric loading. However it was shown that both structures are efficient and could satisfy standard design criterion for roofing structures.

Finally to illustrate the practical feasibility of the double layer reciprocal system, a prototype has been built based on the configuration of Figure 11. The structure is shown in figure 17, it has a span of almost 3 m. The assembly was relatively easy but the chosen members were to soft to really conduct mechanical tests on this prototype.

Fig. 1 :Fig. 2 :Fig. 3 :

 123 Fig. 1: View of a fan composed with three truss elements

Fig. 4 :Fig. 5 :

 45 Fig. 4 : Two way space grid with square pyramid (square on square offset)

Fig. 6 :

 6 Fig. 6: Perspective view and plan view of a configuration with 5 fans

Fig. 7 :

 7 Fig. 7: Plane double layer reciprocal system with 81 fans

Fig. 8 :

 8 Fig. 8: Principle scheme of the double layer grid

Fig. 9 :

 9 Fig. 9: Adjustment of external fans a) engagement length b) additional member.

Fig. 12 :

 12 Fig. 12: Equivalent single layer reciprocal structure

Fig. 13 :

 13 Fig. 13: Normal forces in the single layer dome under a) symmetric b) dissymmetric loading

Table 1 :

 1 Displacements of single and double layer reciprocal systems.

	Vertical	Symmetric loading	Dissymmetric
	displ.			loading
	[in mm]	Average Maximu	Average Average
			m
	Double	7.5	11.5	3.0	6.0
	layer
	Single layer	12	30	5.5	29
	Fan	Symmetric loading	Dissymmetric
	rotation			loading
	[in deg]	Average Maximu	Average Maximu
			m		m
	Double
	layer

Table 2 :

 2 Forces in single and double layer reciprocal systems.

		N min	N max	M max	σ min	σ max
		[kN]	[kN]	[kN.	[MPa	[MPa
				m]	]	]
			Symmetric loading
	Double	-6.0	4.3	0.57	-188	189
	layer
	Single layer	-7.2	2.2	0.71	-245	228
			Dissymmetric loading
	Double	-9.5	9.0	0.59	-194	199
	layer
	Single layer	-8.3	2.6	0.82	-277	269


Minesh Poudel
email: minesh.poudel@airbus.fr

Rakia Kaffa-Jackou

Felipe França
email: felipe@cos.ufrj.br

Félix Mora-Camino

A Multi Agent Systems Approach for Aircraft Emergency Evacuation

In this communication we consider the problem of optimal control of a multi agents system evolving over a bi dimensional grid. Starting from an initial state, the aim of the problem is to find a global strategy so that the system reaches a desired final state in a minimum time. Each agent is characterized by motion parameters while to each cell is associated a current capacity. Objectives and constraints of this optimal control problem are discussed and a solution strategy based on dynamic programming solution approach is proposed for this discrete optimization problem. The application of interest in this study is relative to emergency evacuation of aircraft. The solution of this problem provides minimum time standards as well as insights for the design of personal guidance assistance in emergency situations.

Introduction

Demand for air travel has increased steadily over the last decades and the aviation Industry has forecast that-further substantial growth, nearly doubling of the air traffic, into the next coming decades. These forecasts have led aircraft manufacturers to design and produce airframes capable of carrying as much as nine hundred passengers as well as newer generation airframes made of composites. One of the important aspects from the beginning of the aviation history is that the passenger safety has always been taken with high priority within the industry. Henceforth we have seen substantial improvement in the safety standard of the aviation from design prospective to better operations and maintenance procedures. However, though the rate of accident has decreased drastically in the last three decades, the percentage of passengers surviving after the accident/incident has not decreased in comparison to the improvements achieved in other areas. A survey by the European Transport Safety Council assesses that 40 percent out of the 1500 persons who die every year in aircraft accidents, around 600 passengers die in technically "survivable" accidents. This study shows that more than half of them die from the direct result of the impact, and the others die because of fire, smoke, or problems that arise during the emergency evacuation process. Due to these reasons, not only the issues concerned with the prevention of the occurrence of accidents are tackled with great care but also issues contributing to improving the survival rate in the event of an accident/incident. Accidents can be classified either as fatal (non-survivable), non-fatal (survivable) or technically survivable. There are two ways to prevent fatalities in air travel: by preventing accidents and by protecting aircraft occupants when accidents occur. In order to increase the survivability of passengers in case of an accident, one major area that needs immediate attention is cabin safety. Cabin safety cannot be defined precisely as it covers a very diverse responsibility and interests, which mainly includes crashworthiness, operations, human factors, psychology, and bio dynamics. However, it can be classified in three majors areas, interacting with each other namely: impact protection, fire survivability and emergency evacuation. The focus of the present study is on emergency evacuation, which is an event which seldom occurs at the scale of airlines and that is extremely rare at that of individuals. However, it is under these scenarios that the role of the effective cabin improvements at the design stage have to be taken into account in order to prevent the accident and to develop best methods to deal with the reality of emergency evacuation. Improving survivability will therefore necessitate comprehensive review of all promising options available to regulators and Industry. More research is necessary to be conducted, in order to better understand the reality and prepare effective solutions. After reviewing different modeling approaches a non standard assignment problem is formulated so that cabin crew can be used efficiently during emergency evacuation.

II Cellular Automata

Cellular Automata play an important role in modeling and simulation of spatiotemporal processes. Cellular automata are artificial mathematical models of dynamical systems, discrete in space and in time, whose behavior is completely specified in terms of some local law. A cellular automaton can be thought of as a stylized universe:

space is represented by a uniform grid, with each cell containing a few bits of data, -time advances in discrete steps, -the laws of the universe are expressed through a look-up table, through which at each step each cell computes its new state from that of its close neighbors. Thus, the system's laws are local and uniform. The first cellular automaton used by Von Neumann was qualified as a universal computer while the Game of Life of J.H. Conway had only two states per cell, either filled or empty, 'alive' or 'dead. Essentially cellular automata allows to adopt a cell-based approach to model processes in a two-dimensional space. A typical cellular automata system is composed of four components: cells, states, neighborhood and rules. Cells are the smallest units of the system having adjoining neighbors, they are characterized by discrete states. The state of a cell can change only based on transition rules, which are defined in terms of neighborhood functions. The transition rules are the real engines of change in cellular automata. Their rules control the transformation of a cell state to another cell state over a specific period of time depending on the neighborhood of the cells. The notion of neighborhood is central to the cellular automata paradigm. A simple transition rule in a cellular automaton model is modeled by Li and Yeh as follows:

) ), ( ( ) 1 ( Ν = + t s f t s (1)
with s ∈ S (S is the set of all possible cell states) and where N is the neighborhood of the cell, which acts as inputs for the transition rules. The function f defines the transition rules from time period t to t + 1. The adjacent neighbors are defined by the cells formed by the co-ordinates {x ± 1, y ± 1} as in Moore's neighborhood, with, x = 0, y = 0 as the centre cell. Typically in cellular automata, the neighbourhood and the transition rules play an important role for the automaton to initiate state transitions in a cell. In most cellular automata models, the automata are influenced only within the Moore's neighbourhood, which are the eight adjacent neighbourhood cells. Figure 3 An important component of the cellular automata paradigm is the geometry in the two-dimension space of the structures of the cells. The cells are of the same size and shape, and the value attributed to the cell corresponds to the whole region bound by the cell. Other different regular cell geometries are triangular and hexagonal cells. Even though uniformly regular spaced square cells are used as in the case of classical cellular automata, they are inadequate for an accurate representation of reality. Most objects in reality are not regular and hence they are not square in shape. To counter such situations irregular lattice structures are being introduced in the cellular automata framework.

One Dimensional Cellular Automata

Cellular automata are useful to analyze and understand the laws that govern complex phenomena. In natural sciences, non-linear phenomena are often described by reaction-diffusion equations. Non-linear space-time dynamics of interacting particles, chemical or biological systems can generate numerous local and non local effects far from equilibrium such as steady state multiplicity, oscillations like limit cycles, propagating fronts, target patterns, spiral waves, pulses as well as stationary spatial patterns. Computer simulations using Cellular Automata have been applied with considerable success in different areas, such as chemistry, biochemistry, economy and physics. The cellular automata are discrete mathematical systems taking discrete values in space, time and state. In addition to have auto-replication and universal computation capabilities, they present auto organization capabilities since they can generate ordered behaviors starting from total disorder. This capability is very useful to try to explain certain kind of behaviors observed in physical and biological phenomena. So, the cellular automata have been used to build numerical models of processes as diverse as chemical reactions, diffusion processes, hydrodynamic flows, mechanic, filtration and percolation. With respect to the structure of a cellular automaton it can be considered as a dynamic system that represents a grid of locally connected finite automata. Each automaton produces and output from several inputs, modifying its state in this process by means of a transition function. The state of a cell of a cellular automaton in a particular generation only depends of the states of its neighboring cells and the state the cell had in the previous generation. Cellular automata are formally defined as quadruples:

CA= (d, S, N, f ) (2)
the integer d is the dimension of the space in which the cellular automaton works. -S = {0; 1;: : : ; s-1} is called the set of states.

- the neighborhood N = (n 1 ; : : : ;n v ) is a v-tuple of distinct vectors of Z d . -
the n i 's are the relative positions of the neighbor cells with respect to the current center cell, the new state of which is being computed. The states of these neighbors are used to compute the new state of the center cell. -the local function of the cellular automaton ( f : S v → S) gives the local transition rule.

a configuration is a function from Z d to S. The set of all configurations is .

d Z S C = -
the global function A of the cellular automaton is defined via f as follows:

))

n i c n i c f c A Z i C c + + = ∈ ∀ ∈ ∀ L ( , ), ( ( ) ( , , 1 v d
which correspond for instance to rules such as:

1. If a cell at the time t is inactive (0), is activated at the time t+1 if some of the previous adjacent cells (left) is active (1). 2. An active cell at the time t, is turned inactive at the time t+1 if its adjacent cell (left) is inactive. 3. In other case a cell preserve its previous state.

In Figure 4 a graphical representation of these transition rules is presented and in figure 5 the evolution of this cellular automata is shown.

Cellular Automaton Egress Modeling

This approach has been adopted by Kirchner et al which introduced a model based on cellular automaton where space is discretized into cells (see figure 6) which can either be empty or occupied by one person (in their case a passenger). Each person can move to one of its unoccupied next-neighbor cells (i, j) at each discrete time step → t +1 according to certain transition probabilities p ij . These probabilities are environment dependent. A move is only possible towards one of the direct neighbor cells. For the case of the evacuation processes, the environment of a person is characterized by the shortest distance to an exit door, which can be measured by the minimum number of cells that have to be crossed to reach the exit. The persons move from one cell to another according to some rules. A possible set of rules is the following:

-For each person, the transition probability p ij for a move to an unoccupied neighbor cell (i, j) (including the origin cell, corresponding to no motion) is given by :

ij s S k ij ij ij e n Q p ) 1 ( - = λ (4)
where: n ij = 0 if cell (i, j) is empty and 1 otherwise, and

λ ij = 0 if cell (i, j) is forbidden, λ ij = 1, otherwise.
Coefficient S ij can be taken inversely proportional to the distance from the door measured using a Manhattan metric. Here k s is a positive scaling parameter and N is such that:

1 ] ) 1 ( [ - ∑∑ - = ij S k ij ij ij s e n Q λ (5)
-Each person s makes a probabilistic choice of a target cell according to the updated transition probability distribution { p ij , (i, j) ∈ B(s)}, where B(s) is the set of direct neiborghs to the cell in which person s is presently. -Conflicts arising between two or more persons attempting to move to the same target cell are solved by a probabilistic method: Here a friction parameter μ ∈ [0, 1] is introduced so that in a conflict the motion of all involved persons is denied with probability μ, while one person is allowed to move to the desired cell with probability 1-μ. The person, which actually moves, is chosen randomly with equal probability between the persons involved in the same conflict. Several advances are being made to explore the possibilities of the cellular automata technique. Notable among them are calibration approaches for constrained cellular automata models, which have been applied to model slow spatial processes. These models are constrained by the causal factors driving the spatial processes. Local constraints contain detailed information for each cell, but sector constraints have only aggregated or partial-spatial information. Global constraints, however, are characterized by temporal or nonspatial information. Such an approach of defining the constraints and operating in cellular automata is a step towards achieving linkages with the causal factors. Once the simulation is executed, the constraints would be static over time. However, there are many systems for which the constraints are dynamic dependent on the situations. And so it would require an additional framework to address issues concerning the dynamics of these causal factors and constraints. Further, in such models these constraints would not behave as an event-based system with a cause-effect relationship. Many models need to account for the external drivers that are not accounted in the transition rules of the cellular automata. Certain externalities can be system wide or specific to certain locations, for which the cellular automata models have to evolve to address such requirements. Further there are possibly also different significant processes that take place in the region in question, apart from those represented in the cellular automata model. Cellular automata models are yet incapable of representing the external factors responsible for driving the change dynamics as the transition rules account only for the states and neighborhood. To overcome such limitations, different approaches are being suggested. Among them is the integration of agent-based models over a cellular automata framework, as agent-based models can be constructed to represent the externalities driving the processes. Thus the current research is approaching towards the integration agent-based models (multi-agent systems) with the cellular automata models, such as in the case of modelling the dynamics of emergency evacuation by incorporating different drivers as agents involved in enabling the individual spatial interactions by defining the spatial and temporal relationships to these agents.

Agent Based Modeling

Agents, have their origins in software engineering and artificial intelligence where they are used in networking, communications and many more applications. The aim of agent design is to create a program, which interacts with its environment. The term 'agent' is usually applied to describe self-contained programs, which can control their own actions based on their perceptions of their operating environment. A significant definition is that, an agent is considered as a self-contained program capable of controlling its own decision-making and acting, based on its perception of its environment, in pursuit of one or more objectives. The most general way in which the term agent is used is to denote hardware or a software-based computer system that enjoys the following properties:

-Autonomy: Agents operate without the direct intervention of humans or others, and have some kind of control over their actions and internal state; -Social ability: Agents interact with other agents (and possibly humans) via some kind of agent-communication language; -Reactivity: Agents perceive their environment, (which may be the physical world, a user via a graphical user interface, a collection of other agents, or perhaps all of these combined), and respond in a timely fashion to changes that occur in it; -Pro-activeness: Agents do not simply act in response to their environment; they are able to exhibit goal-directed behavior by taking the initiative. There are a number of points about this definition that require further explanation. Agents may be characterized by the following properties:

clearly identifiable problem solving entities with well-defined boundaries and interfaces; -situated (embedded) in a particular environment-they receive inputs related to the state of their environment through sensors and they act on the environment through effectors; -designed to fulfill a specific purpose-they have particular objectives (goals) to achieve; -autonomous, i.e. they have control both over their internal state and over their own behavior; -capable of exhibiting flexible problem solving behavior in pursuit of their design objectives. They need to be both reactive (able to respond in a timely fashion to changes that occur in their environment) and proactive (able to act in anticipation of future goals). Although, the origins of agent-based models have been in the artificial intelligence, they are also developed in social sciences extensively. Agents can be considered as a special case of an automaton, having all features of the general automaton, with a distinction that these agents are mobile and they can represent the external drivers responsible for the processes (e.g. socioeconomic, population, etc.). The idea is to treat each of the individual drivers as agent-based automata enabling the spatial and temporal relationships. There can be as many agent-based models as the number of externalities identified driving the processes at appropriate scales. Such processes take place at specific locations and are generally not system wide. While a classical cellular automata transition rule is system wide, such agent-based models would only be specific to certain locations only. These agent-based models are to act in conjunction with the regular transition rules of the cellular automata. The essential concepts of agent-based computing are agents, high-level interactions and organizational relationships. It can be seen that there can be numerous agents wherein specific agents can interact amongst themselves and/or have the same sphere of visibility and influence. There can also be agents who act independently without any interaction with other agents. The classification of agents has been attempted by different investigators. The classifications are made based on the characteristics (autonomy, reactivity, social ability and pro-activeness) and applications of the agents. Agents can be classified in three broad categories according to the mobility of these agents , the agents are classified as static or mobile, according to the ways in which these agents are modeled, the agents are categorized as 'reactive' and 'deliberative' and according to several primary attributes that agents should exhibit. The three characteristics autonomous, cooperation and learning are used to derive the four types of agents: collaborative agents, interface agents, collaborative learning agents and smart agents. Multi-agent based simulation is used in a growing number of areas, where it progressively replaces the various micro-simulation, object-oriented or individual-based simulation techniques, previously used. It is due, for the most part, to its ability to cope with very different models of "individuals", ranging from simple entities (usually called "reactive" agents to more complex ones ("cognitive" agents. The easiness with which modelers can also handle different levels of representation (e.g., "individuals" and "groups", for instance) within an unified conceptual framework is also particularly appreciated, with respect, for instance, to cellular automata. This versatility makes multi agents based simulation emerge as the support of choice for the simulation of complex systems, and it is appealing to more and more scientific domains: sociology, biology, physics, chemistry, ecology, economy, etc.

Multi-Agents Systems and Cellular Automata

In a cellular automata model, space is represented as a uniform lattice of cells with local states, subject to a uniform set of rules, which drives the behavior of the system. These rules compute the state of a particular cell as a function of its previous state and the states of the adjacent cells. Several extensions have been developed:

-An extension of the basic cellular automata model allows the state of any particular cell to be influenced by more the states not only of the contingent cells, but also the by the states of more remote cells. State changes may depend on the aggregate effect of the states of all other cells, or some of their subsets. -Another extension is to build models in which cells preserve state information and calculate their next state on the basis of their neighbors and their own history of state changes. -Agent technology has been implemented to build a framework for multi-agent simulation. Objects (or people) moving across the network are represented in terms of autonomous agents. Each agent will be located in a simulated space, based on the cellular automata grid. The choice of a multi-agent system is motivated by their promise to simulate autonomous individuals and the interaction between them. Agent technology is also used to simulate the outcome of the model and the simulation. Designers can use the system to assess the likely consequences of their design decisions on user behavior. The application of cellular automata implies the possibility to simulate how an 'agent'-user moves in a given environment, dependent of the behavior of other agents in the system. In developing the simulator it is useful to differentiate between the cellular automata part and the distributed artificial intelligence agent's structure. The agent's structure part involves the different agents with their respective roles. Various agent types may be distinguished in the model such as user-agents that represent people in the simulation. Here the passenger can be considered to be the subject-agent. Actor-agents can simulate the crew. Thus, subject-agent and actor-agents are user-agents that navigate in the virtual environment network, each with their own behavior, beliefs and intentions. A belief is the internal, imperfect representation of the virtual environment including the state of other user-agents, on which the decisions of an agent are based. We must view behavior as the interaction between the user agents and the environment. For the subject-agent, this behavior is not an attribute of the agent, but rather lies in the mind of the subject alone. Styles of behavior like anticipated behavior and unplanned behavior can be considered. Besides a behavior-agent, we also distinguish an intra-task-agent. An intra-task-agent fulfils the intentions of a user-agent to reach a destination (goal) and/or to carry out a list of activities (a plan). We define a user-agent as the 3-tuple:

U = < R, A, F > (6)
where -R is a finite set of role identifiers, it represents the enumeration of all possible roles that can be played by user agents; -A represents the activity agenda of user-agent i to perform the goals and desires ( {Ai} ). -F represent the knowledge or information about their environment which an user agent possesses ( {Fi} ), which is called facets. Part of these facets is amongst others beliefs, awareness, experience, preference and choice. All these facets are dynamic and may change over time in the simulation loop, influenced by the occupation of other cells in the network. The population of crowded cells could result in emergence behavior of the crowded neighborhood. A regular lattice of width W and length L, and density (that represents the mesh of the network) can be introduced. Cells within the lattice are given the notation c k,i,j , where 1≤i≤W and 1≤j≤L and k indicates the cell type where k∈{empty, decision, wall). Here:

-"empty" means that the cell belongs to the walkway; -"decision" means the cell belongs to the decision-point area where the passenger takes the decision where to go; -Wall means that the cell is part of a frame (a wall, etc.). A cell x (actually c t,i,j ) in the aisle to be on (s(x)=1) if it is occupied by a passenger, otherwise it is off ( s(x) = 0). Also a density size δ r (x) which shows the activity around cell x can be defined; it shows the number of neighbors in an on state in relation to the total number of neighbors in the Moore neighborhood of cell x with radius r. In a d-dimensional grid with a Moore radius r, the number of cells n in the neighborhood of cell x (x included) = (2r+1)d. The updated cell x depends on the on states of all cells in the neighborhood N of x (N(x)). To summarize the transition of a state of a cell,

The cellular automaton micro simulation can proceed from one time instant to the next. In each time step lane assignment and speed updates change the positions of all pedestrians according to local rules applied to each pedestrian.

A Stochastic Cell Model

Here is considered a systems composed of I agents located at a grid with N positions. Among these positions, a subset N s is composed of safe positions. The current state of agent i is represented by his position (i(n), j(n)) in the grid. When a system obeys the Markov property, the future is determined only by the present and not by the past. The Markov property can be stated in the following way: the configuration of the system at time depends only on the configuration of the system at time t. This fact allows to write the probability distribution at time

= = Δ + = - = = Δ + = = = - Δ + = ∑ ∑ r r (
= - = = Δ + = = = - Δ + = ∑ (13)
At the start of evacuation, the cabin is composed of cells which are in the following possible states: undamaged, partly damaged but transitable, destroyed or untransitable (crushed, blasted, burnt, drowned). Possible active hazards such as fire, smoke and water have a starting area covering a given set of destroyed or untransitable cells. Then propagation models of the present hazards are given by:

{ } f t t k t t t t t t F , , , , , ), ( 0 0 0 L L δ δ + + ∈ (14) { } f t t k t t t t t t S , , , , , ), ( 0 0 0 L L δ δ + + ∈ (15) { } f t t k t t t t t t W , , , , , ), ( 0 0 0 L L δ δ + + ∈ (16)
where F(t) is the set of cells affected by fire at period t, S(t) is the set of cells affected by dense smoke at period t and W(t) is the set of cells affected by water at period t. Here t 0 is the initial time and t f is the final period of the simulation. Then it is possible to identify at each period the feasible exits as well as the passengers which are affected by the different active disasters. Agent i is alive at time t if the three following conditions are met:

* ' ∈ n n α n n' α { } 1 , 0 ∈
is the inverse of the number of direct neighbours to cell n. The positive weight can be chosen such that:

n w { } [ ]) ( 1 1 min min , , 1 kj n kj N k C j n x w k n x w = - + = ∈ ∈ λ L (24)
where λ is a positive parameter and =1 if crew member number j is at position k , and =0 otherwise and where

I i ∈ ∀

Then the crew location problem during egress can be formulated as the following discrete optimization problem: Of course this is a non standard assignment problem since it includes dynamical aspects, stochastic components and dynamic opponents (hazards). An approximate solution of this problem, an heuristic based on flows in networks considerations, can be designed to tackle efficiently this problem.

Conclusion

This paper has discussed different modeling approaches that are currently used for modeling complex systems whose dynamics is characterized by the evolution, often competitive, of many individual agents: cellular automata, agent based simulation, multi-agents simulation and stochastic cell models. There, simulation is based on mathematical models that represent the temporal evolution of location of individual agents from cell to cell. It appears that all these modeling approaches present large limitations with respect to their application to emergency evacuation representation: the space in which agents move is in general composed of identical adjacent cells and cannot be easily adapted to represent realistically the confined cabin space; hazards dynamics are hardly considered; the motion of agents is driven by over simplified logics mainly based on the occupancy of neighboring cells; the behavior of the agents is assumed to be homogenous: many often there is no differentiation between the behavior of agents, no specific group behaviors are also considered. This has led to propose a stochastic model to represent egress dynamics and hazard progression as well as passengers health evolution during egress. Then an optimization problem considering the localization of cabin crew members during evacuation has been established. The aim of this problem is to locate optimally the cabin crew members so that they can provide efficient directives to evacuating passengers so that at the end of evacuation, the maximum number of life is saved.

Figure 1 .

 1 Figure 1 . A310: March 07-Dubai International Airport

Figure 2 .

 2 Figure 2. Exit Utilization During an Accident

  depicts respectively the notions of neighbourhood developed by von Neumann and by Moore.

Figure 3 .

 3 Figure 3. Cells and Neighbourhoods

Figure 4 .Figure 5 .

 45 Figure 4 . Example of Transition Rules

Figure 6 .

 6 Figure 6 . Discretization of Aircraft Cabin Space

  constraints[START_REF] Poudel | A Fuzzy Logic Approach for Aircraft Evacuation Modeling[END_REF], (25), (26) and the classical assignment constraints:

  but in general near to 1. The set of safely rescued passengers, N srp is such that:

					x	kj																																					x	kj	n -	k	is a
	distance on the grid of positions n and k.
	The emergency evacuation can be considered completed at time t f when:
											t f min =	t			with		N n s ∈ ∑	P	( i	=	n	,	t	)	≥	σ	i	i ∀	∈	L	( t	)	(25)
	where i σ is a threshold parameter with	∈	[	, 0	] 1
	i	∈	N	srp	if	0 t f ∫ ∑ n	P	( i	=	, n	t	)	ε	F	,	n	( t	)	dt	≤	σ	F		,	t 0 f ∫ ∑ n	P	i (	=	n	,	t	)	ε	S	,	n	t (	)	dt	≤	σ	S	(26)
					and			t ∫ ∑ n f 0	P	i (	=	n ,	t	)	ε	W	,	n	t (	)	dt	≤	σ	W

i σ

where

where

Here σ F , σ S and σ W are positive threshold levels. It is supposed that the transition matrices are such as: