domain: spi phys

text:
Mathieu Rognant
email: mathieu.rognant@onera.fr

Sofiane Kraiem
email: sofiane.kraiem@isae-supaero.fr

Jurek Sasiadek
email: jurek.sasiadek@carleton.ca

Kinematic Indices of rotation-floating space robots for on-orbit servicing

Keywords: Space robot, On-orbit servicing, Kinematic indices

published or not. The documents may come

Introduction

The orbit removal and servicing topics are expected to play a key part in future space missions [START_REF] Long | On-orbit servicing: a new value proposition for satellite design and operation[END_REF], and the studied concepts often include robotic manipulators to carry out these missions [START_REF] Flores-Abad | A review of space robotics technologies for on-orbit servicing[END_REF]. Robotic manipulators in space are already in use but are teleoperated given the complexity to steer them. Umanned servicing missions have not been perform yet and fully autonomous control of spacecraftmanipulator is still an open research area, specially with non-cooperative target.

In litterature two control modes of the spacecraft-manipulator are commonly considered: the free-flying mode where manipulator's joints and spacecraft's attitude and position are actively controlled [START_REF] Shibli | Unified modeling approach of kinematics, dynamics and control of a free-flying space robot interacting with a target satellite[END_REF]; and the free-floating mode where only manipulator's joints are controlled [START_REF] Yoshida | Modeling of collision dynamics for space free-floating links with extended generalized inertia tensor[END_REF]. More recently Wilde and al. proposed in [START_REF] Wilde | Equations of motion of free-floating spacecraft-manipulator systems: An astronautical engineers tutorial[END_REF] to distinct three additional modes to fully covering all possible spacecraft maneuvers: rotation flying and translation flying modes which are sub-modes of the free-flying and require reaction-jet thrusters. These actuators allow to actively control the angular and/or linear momentum of the system but consume fuel and their using shorten the lifespan of the spacecraft. The third proposed mode, (rotation floating), is an extension of the free-floating mode in which the attitude of the satellite is actively controlled by using kinetic moment exchange devices. In this case the total kinetic momentum of the system are conserved and the used actuators (reaction wheels or control-moment gyrosocopes) require only electric power, which does not restrict the life of the spacecraft.

The free-floating manipulator kinematic was well known [6][7], and studies on space manipulator which consider kinetic moment exchange devices are focussed on the compensation of the reaction manipulator torque on the base [START_REF] Ellery | Robotic in-orbit servicers the need for control moment gyroscopes for attitude control[END_REF][9] [START_REF] Wu | Attitude control for on-orbit servicing spacecraft using hybrid actuator[END_REF].Indeed as satellites integrate such devices to achieve attitude control, these actuators are most of the time controlled separately from the joint manipulator [START_REF] Jaekel | Design and operational elements of the robotic subsystem for the e. deorbit debris removal mission[END_REF] and only few address the issue of the angular momentum distribution of the system [START_REF] Oki | Time-optimal manipulator control for management of angular momentum distribution during the capture of a tumbling target[END_REF] and its implication to the base actuators sizing [START_REF] Kim | Energy-efficient trajectory generation for space manipulators with reaction wheels under a fixed base orientation[END_REF].

In this paper, a common control of joint and kinetic moment exchange devices is considered to achieve the desired motion of the spacecraft base and manipulator end-effector. To address design and control issues and evaluate the interest of a common control, we propose some kinematic indices which hightlight the contribution of each actuator on the controlled states.

This paper is divided as follows: first, the kinematics and dynamics for a rotation-floating robot are presented in Sec. II; then, in Sec. III, kinematic sensitivity indices are derived and their meanings are discussed; finally Sec IV outlines a simulation example which serves to illustrate the contribution of the proposed indices.

2 Equation of motion of Rotation-floating space robots In this paper, Denavit-Hartenberg description is used [START_REF] Khalil | Modelisation, identification and control of robots[END_REF] to express the end-effector position 0 p e according to the joint configuration vector q, and the end-effector relative speed to the spacecraft base reference frame 0 :

Kinematic

0 v e 0 ω e = J m q (1)
where q is the vector of the joint velocity ( q = q1 • • • qn T ) and J m the manipulator jacobian matrix which is defined in R 6×n [START_REF] Xu | Hybrid modeling and analysis method for dynamic coupling of space robots[END_REF].

As the base of the spacecraft is also moving the total kinematic equation of the end-effector of the space robot according to inertial reference frame ine expressed in the base B 0 is

v e ω e = J b ẋ0 + J m q (2)
where ẋ0 is the vector of the linear and angular velocity according to inertial reference frame ine expressed in the base B 0 ( ẋ0 = v T 0 ω T 0 T ), and

J b ∈ R 6×6
is the base jacobian matrix [START_REF] Xu | Hybrid modeling and analysis method for dynamic coupling of space robots[END_REF].

Dynamics equation of motion

The coupled dynamics of a spacecraft-manipulator system have been extensively studied in the available literature [START_REF] Dubowsky | The kinematics, dynamics, and control of freeflying and free-floating space robotic systems[END_REF][16]. To address the issue of control of rotation-floating robots, we propose to extend the free-floating dynamic model detailed in [START_REF] Wilde | Equations of motion of free-floating spacecraft-manipulator systems: An astronautical engineers tutorial[END_REF] by adding the states of the kinetic moment exchange actuators. By using the Lagragian approach, we just need to express the kinetic energy of the attitude actuator and add the term obtained by derivation of its kinetic energy according to the system state and time. In this paper only simple case with reaction wheels will be considered, but this approach can be easily extended to control-moment gyro. The kinetic energy T rj of a reaction wheel of index j is:

T rj = 1 2 ω T 0 + ω rj 0 s j T H rj ω 0 + ω rj 0 s j + 1 2 m rj ṗT 0 ṗ0 (3)
where ω rj is the angular velocity and 0 s j the rotation axis expressed in B 0 . m rj and H rj are respectively the mass and the inertia matrix of the reaction wheels expressed in B 0 .

In this paper, one serial manipulator with n-DOF is mounted on the satellite. However, multiple manipulators and manipulator with complex kinematic could be considered as well [START_REF] Wilde | Equations of motion of free-floating spacecraft-manipulator systems: An astronautical engineers tutorial[END_REF]. Consequently, Lagrange equation derivation of (3) leads to the following dynamics equation of motion for spacecraft with n r reaction wheel:

  H 0 H 0m H 0r H T 0m H m 0 H T 0r 0 H r     ẍ0 q ωr   +   c 0 c m c r   =   0 τ m τ r   (4)
The symbols used in this formulation are defined as follows:

H 0 ∈ R 6×6 : inertia matrix for the base.

H m ∈ R n×n : inertia matrix for the links. H 0m ∈ R 6×n : coupling inertia matrix between the base and the manipulator. H r ∈ R nr×nr : inertia matrix for the reaction wheels. H 0r ∈ R 6×nr : coupling inertia matrix between the base and the reaction wheels. c 0 ∈ R 6×1 , c m ∈ R n×1 and c r ∈ R nr×1 : are velocity dependent non-linear terms. τ m ∈ R n×1 and τ r ∈ R nr×1 : are manipulator joint and reaction wheels torques.

H 0r T = 0 nr×3 H ωr T
, and H r is a n r × n r diagonal matrix whose the i th component is equal to 0 s T i H ri 0 s i . The H ωr matrix is defined by:

H ωr = H r1 0 s 1 • • • H rn r 0 s nr (5)

Generalized Jacobian

Based on the assumptions that the torques caused by the environment are negligible and the initial kinematic momentum of the spacecraft is null, the conservation of linear momentum and the angular momentum of the system yields to the generalized jacobian of the end-effector motion [START_REF] Yoshida | Modeling of collision dynamics for space free-floating links with extended generalized inertia tensor[END_REF]. In this paper this expression are extended in order to discern the contribution of base actuator and the joint actuator.

H 0 ẋ0 + H 0m q + H 0r ω r = 0 (6)
The inertia matrix H 0 is symetric and positive-defined then

ẋ0 = -H -1 0 H 0m q -H -1 0 H 0r ω r (7)
From ( 7) one can express jacobian matrices J 0 q ∈ R 6×n and J 0 r ∈ R 6×nr which link respectively the base kinematic to the joint and reaction wheel velocities:

J 0 q = -H -1 0 H 0m and J 0 r = -H -1 0 H 0r (8)
By considering equation ( 2), one can also link the end-effector motion to the joint and reaction wheel velocities:

ẋe = (J m -J b H -1 0 H 0m ) q -J b H -1 0 H 0r ω r (9)
and define the jacobian matrices J e q ∈ R 6×n and J e r ∈ R 6×nr :

J e q = (J m -J b H -1 0 H 0m ) and J e r = -J b H -1 0 H 0r (10)

Kinematic Indices

According to these models we propose in this section to analyze the participation of each actuators on the base motion. Then we evaluate the end-effector manipulability variation according to the reaction wheel control mode.

Base kinematic indices

The equation ( 6) can be decomposed as follows:

m tt I 3 -m 0 tt p × g m 0 tt p × g H ω ṗ0 ω 0 + H vm H ωm q + 0 H ωr ω r = 0 0 (11)
where 0 p × g is the skew-symmetric matrix of 0 p g which is the vector of the system center of gravity expressed in B 0 ; m tt and H ω are respectively the total mass of the system and the inertia matrix expressed in B 0 ; and H vm and H ωm are sub-matrices of the coupling inertia matrix H 0m introduced in (4).

Based on this decomposition of the spacecraft base inertia matrix H 0 , this relation can be inversed in analytic way by using a the Banachiewicz inversion formula [START_REF] Baksalary | Generalized inverses of partitioned matrices in banachiewicz -schur form[END_REF] which yields to the following expression of H -1 0 as detailed in [START_REF] Wilde | Equations of motion of free-floating spacecraft-manipulator systems: An astronautical engineers tutorial[END_REF]:

H -1 0 = I3 mtt -( 0 p × g )S -1 U ( 0 p × g ) ( 0 p × g )S -1 U S -1 U ( 0 p × g ) S -1 U with S U = H ω + m tt ( 0 p × g )( 0 p × g ) (12
) One can deduce from ( 8) and ( 12), the expression of the linear and angular part of the jacobian matrices J 0 q and J 0 r :

J v0 q = - I3 mtt -( 0 p × g )S -1 U ( 0 p × g ) H vm + ( 0 p × g )S -1 U H ωm and J v0 r = -( 0 p × g )S -1 U H ωr (13)
J ω0 q = -S -1 U ( 0 p × g )H vm + S -1 U H ωm and J ω0 r = -S -1 U H ωr (14)
An equivalent decomposition of the translation and rotation have been proposed by Xu in [START_REF] Xu | Hybrid modeling and analysis method for dynamic coupling of space robots[END_REF] to compute Joint-to-Base Coupling Factors. In order to deal with the range variation of each actuators, the following scaled jacobians will be used for analysis:

J v0 = J v0 q Q max J v0 r Ω max and J ω0 = J ω0 q Q max J ω0 r Ω max (15)
where Ω max and Q max are diagonal matrices of the maximal reaction wheels and joints angular speed. By considering relation [START_REF] Khalil | Modelisation, identification and control of robots[END_REF], one can notice that the J ω0 q and J ω0 r are both factorized by S -1 U . That implies that the relative contributions of the joints and reaction wheels to the base angular motion do not depend on the base inertia but only on the relative angular momentum capabilities of the manipulator and the base spacecraft actuator. This highlights that the distribution of the angular momentum is a critical point to insure the controlability of the system and dedicated control algorithms are required [START_REF] Oki | Time-optimal manipulator control for management of angular momentum distribution during the capture of a tumbling target[END_REF].

A dedicated analysis of the components of S U J ω0 matrix of the system in a design phasis could address this issue. The sum of the reaction wheels angular momentum is bounded by S U J ω0 r Ω max . Thus, in the angular momemtum space, this sum belongs to the invariant polyedron whose the faces are normal to the s j axis. Concerning the manipulator, its angular momentum is bounded by

S U J ω0
q Ω max , which is strongly sensitive to the joint configuration (variation of the gravity center position and manipulator inertia ( 14)). Thus, to distinguish the joint and reaction wheel contributions to the base motion according to the joint configuration, we proposed the indices κ q ω0 and κ r ω0 defined as follow:

κ q ω0 = S U J ω 0 q Qmax 2 F S U Jω 0 2 F ; κ r ω0 = S U J ω 0 r Ωmax 2 F S U Jω 0 2 F ( 16
)
where • F is the Frobenius norm. As this norm can be evaluated analytically, these indicators could be an interesting tool to caracterize the base controlability on pre-design stages or select preferential joint configuration for control purpose. For example, joint configurations whose the index κ q ω0 respect the constraint κ q ω0 ≤ 0.5 are more suitable to achieve motion without changing the base attitude.

End-effector indices

The end-effector control motion is the highest priority task, however the dynamic coupling with the motion base could significantly reduce the manipulator capability, specially when fixed base attitude is required. To analyse this point we propose in this section dedicated manipulabity indices.

As for the base, one can defined from [START_REF] Wu | Attitude control for on-orbit servicing spacecraft using hybrid actuator[END_REF] the scaled jacobians of the linear and angular motion of the end-effector:

J ve = J ve q Q max J ve r Ω max and J ωe = J ωe q Q max J ωe r Ω max (17)
Then, the end-effector manipulability can be caracterized by the operator µ(J) = det (JJ T ) [START_REF] Umetani | Workspace and manipulability analysis of space manipulator[END_REF] µ ve = µ(J ve ) and µ ωe = µ(J ωe )

As proposed by Nenchev and al. in [START_REF] Nenchev | Reaction null-space control of flexible structure mounted manipulator systems[END_REF] for redundant manipulator, null motion projector could be used to restrict the joint motion to a sub-space which insure zero base attitude disturbance. In the same manner null motion projector could be defined for a common control of the reaction-wheel and joint velocities as follows:

N ω0 = I -J † ω0 J ω0 (19)
where J † ω0 is the pseudo-inverse of the J ω0 . By using this projector, one can define manipulability indices of the end-effector without base motions:

µ n ve = µ(J ve N ω0 ) and µ n ωe = µ(J ωe N ω0 ) (20)

Simulation Example

In this section numerical simulations are provided to illustrate the kinematic indices proposed in this paper. All these evalutions are carried out by using the physical parameters of the first microsatellite in the Myriade series, Demeter [START_REF] Pittet | Demeter: A benchmark for robust analysis and control of the attitude of flexible micro satellites[END_REF]. Demeter is equipped with three reaction wheels, which nominaly present an inertia of 4.10 -4 kg.m 2 on their main axis and a maximal angular speed of 2800 rpm. Concerning the manipulator, the kinematic chain represented in figure 1 will be used. As this study is focussed on the end-effector linear speed, only the 3 first DOFs will be considered. As current space manipulator such as Canadarm2 and JEMRMS do not exceed the 0.7 rpm [START_REF] Laryssa | International space station robotics: a comparative study of era, jemrms and mss[END_REF] we assume that joints angular speed should not overpass 0.5 rpm and the used links mass, dimension and inertia are summarized in Table 1. In a context of pre-design, the sizing of 1. Model parameters the spacecraft actuators is a critical task. Therefore, in this example, we proposed using the suggested indices to compare two configurations of reaction wheels: a first one equivalent of Demeter's and a second one with twice the nominal inertia H r . To achieve this analysis, the indices have been evaluated for the two configurations in 6859 joints configurations uniformly distributed in the half-workspace of the manipulator because the system is symmetric about the XZ plane. In figures 2, 3 and 4 iso-surfaces of these indices are plotted and colours reflect value. In figure 3 and 4 manipulability indices are normalized by the mean manipulability of the manipulator on a fixe base. Figure 2 points that increasing the reaction wheels inertia allows to decrease the contribution of the joints and respectively to raise the one of the reaction wheels on the base attitude motion, as κ q ω0 and κ r ω0 are complementary. One can also note that, the joint configurations which involve the most the base motion are those where the end-effector is farthest from the system center of mass. In these critical configurations, doubling the inertia of the reaction wheels allow to reduce the joints impact on the base from 0.65 to 0.4 and to ensure a stable attitude without using thrusters. However, the manipulability of v e , represented in figure 3, increases a less than 0.1 if the base is able to move during the control of the end-effector. This is due to the small ratio between the maximal spacecraft angular speed and joints angular speed considered in our example. In this case, reducing the joint impact to the base with bigger reaction wheel inertia does not affect the manipulability of v e . Nevertheless, when null-motion is used, figure 4 highlights an improvement of the manipulability of v e with more angular momentum in reaction wheels. As underlined by the figure 2, the second configuration presents a better momentum repartition for this type of control.

Conclusion

This paper presents some aspects of control for rotation-floating space robot for on-orbit servicing. Special indices for analyzing dynamic capabilities of the robot were suggested and used in simulation experiments. Simulation examples to illustrate kinematic indices were shown with physical parameters for a microsatellite from Myriade series Demeter. This paper compared two configurations of spacecraft actuators: the first, equivalent to Demeter configuration and the second with twice as large, the nominal inertia. It was shown that increasing the reaction wheels inertia allows to reduce the contribution of the joints and raises respective wheels contribution to the motion of robotic spacecraft. Also, it could be noted that the manipulability of the arm increases, if the base is able to move during the end-effector control. The use of suggested indicators allowed to evaluate a system and determine momentum distribution between manipulator and the base of spacecraft.

Fig. 1 .

 1 Fig. 1. Notation for a spacecraft with three reaction wheels and 4-DOFs manipulator

Fig. 2 .Fig. 3 .

 23 Fig. 2. Variation of the joints contribution to the base motion (κ q ω 0 ) plotted in 0 for the initial reaction wheel (left) and the second actuators configuration (right)

Fig. 4 .

 4 Fig. 4. Evolution of the end-effector normalized manipulability indices with a fix base attitude constraint plotted in 0 for the initial reaction wheel (left) and the second actuators configuration (right)


domain: spi phys math info

text:
Eric Nguyen Van
email: eric.nguyen-van@isae.fr

Daniel Alazard

Philippe Pastor

Carsten Döll
email: carsten.doll@onera.fr

Co-design of aircraft vertical tail and control laws using distributed electric propulsion

Keywords: Distributed Electric Propulsion, Differential Thrust, Co-design, H ∞ control, vertical tail sizing, actuator bandwidth

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

INTRODUCTION

Distributed Electric Propulsion allows more flexible airframe propulsion integration and exploitation of favorable synergy between propulsion, aerodynamic or flight control as presented by [START_REF] Felder | Turboelectric distributed propulsion engine cycle analysis for hybrid wing body aircraft[END_REF], [START_REF] Hermetz | Distributed Electric Propulsion for small business aircraft -A concept-plane for key-technologies investigations[END_REF] and [START_REF] Schmollgruber | Multidisciplinary Exploration of DRAGON: an ONERA Hybrid Electric Distributed Propulsion Concept[END_REF]. This study lies in the context of increase of control authority made possible by active differential thrust. In previous studies, it has been shown that the use of differential thrust necessitates a reduction in vertical tail surface area to increase the flight envelop, see Nguyen- [START_REF] Nguyen-Van | Towards an aircraft with reduced lateral static stability using electric differential thurst[END_REF]. Additionally, with a distributed propulsion, the consequences of one engine failure are less critical and allows to reduce the vertical tail surface area while remaining compliant with flight safety, see Nguyen- [START_REF] Nguyen-Van | Reduction of vertical tail using differential thrust: Influence on flight control and certification[END_REF]. These studies showed the interest of differential thrust in the reduction of vertical tail. Remains the question of how much to reduce the vertical tail given a certain distributed propulsion. A co-design approach similar to [START_REF] Denieul | Multi-control Surface Optimization for Blended WingBody Under Handling Quality Constraint[END_REF] is selected to minimize the vertical tail while designing the longitudinal/lateral control laws and the engine allocation module. The optimization procedure includes also a minimization of the engine bandwidth which is a determining parameter for differential propulsion control since the presumed low engine bandwidth, with respect to classical control surface bandwidths, could bring loss of control authority. In this approach, the vertical tail surface area is thought as a linearly varying parameter. Optimization tools available for robust control techniques are then used to handle the variation of matrices with the variation of vertical tail surface area and allows the minimization of control gains, engine bandwidth and ver-This work is supported by AIRBUS and ONERA.

tical tail surface area under longitudinal/lateral handling quality constraints.

The second section describes the aircraft configuration considered in this study and the derivation of its linear model, parameterized according to the vertical tail surface. The model dynamics are analyzed for a given flight operating condition selected for its relevance regarding the design of distributed propulsion control. The third section describes the co-design problem and optimization process based on a multi-step optimization. The fourth section presents the obtained results. Conclusions and perspectives are summarized in the last section.

DETERMINATION OF LINEARIZED SYSTEMS

In this section, the procedure to obtain the Linear Fractional Representation of an aircraft with the vertical tail as varying parameter is explained.

Aircraft Equations of Motion

The equations of flight in the aerodynamic frame, assuming uniform wind velocity are used. Equations are derived in [START_REF] Boiffier | The Dynamics of Flight[END_REF] and are presented here in a compact form for the sake of brevity. The interested reader is referred to [START_REF] Boiffier | The Dynamics of Flight[END_REF] and Nguyen- [START_REF] Nguyen-Van | Towards an aircraft with reduced lateral static stability using electric differential thurst[END_REF] for a more detailed treatment. Let A be the acceleration vector and Ω the vector of rotation rates, the equations of motion write:

mA = mg + F A + F T , (1)
I Ω + Ω × IΩ = M A + M T , (2)
with m the mass, g the gravity vector, I the inertia matrix, F A , F T and M A , M T are respectively the force and torque vectors due to aerodynamic and engines thrust.

Developing the acceleration term, A, while assuming zero wind velocity, allows to rewrite equation (1) as:

m V + T ab Ω × V = mg + F A + F T ,
(3) where T ab is the rotation matrix between the body carried frame and the aerodynamic frame and V is the velocity vector in the aerodynamic frame. Equation (3) allows to explicit the state variables: V , α, β, being the airspeed, the angle of attack and the side slip angle, and p, q and r, the roll, pitch and yaw rates:

m   V βV αV cos β   + T ab p q r × V = mg + F A + F T , (4)
I ṗ q ṙ + p q r × I p q r = M A + M T . (5)
The complementary kinematic equations for the Euler angles: bank φ, pitch θ and heading ψ are :

  φ θ ψ  =    1 sin φ tan θ cos φ tan θ 0 cos φ -sin φ 0 sin φ cos θ cos φ cos θ    p q r . ( 6
)
Two additional parameters are necessary to find pseudoequilibriums: the flight path angle γ and the turn rate Ω.

For steady turns Ω is equal to the rate of change of heading and using the last line of equation ( 6): sin γ = cos α cos β sin θsin β sin φ cos θ sin α cos β cos φ cos θ,

Ω = (q sin φ + r cos φ)

1 cos θ . ( 8
)
Leaving aside the last row of system (6), the system to solve consists in equations ( 4), ( 5), ( 6), ( 7) and ( 8) representing a set of N e = 10 equations. The state vector is x = [V, α, β, p, q, r, φ, θ] T and counts n x = 8 variables. The directional control being entirely made with differential thrust, the input vector corresponding to control surfaces is composed of ailerons and elevator deflection: u = [δ a , δ e ] T with n u = 2. Finally, n p = 2 additional parameters γ and Ω are counted.

Reference Aircraft

Commuter and regional aircraft are often cited as the next big in developing electric airplanes since most of their missions are within the limits of electric propulsion in terms of endurance, see [START_REF] Moore | Misconception of electric propulsion aircraft and their emergent aviation market[END_REF] and [START_REF] Stückl | Methods for the Design and Evaluation of Future Aircraft Concepts Utilizing Electric Propulsion Systems[END_REF]. A good representative of this class of aircraft is the ATR72 which details are reported in table 1. For the configuration with distributed electric propulsion, the number of engines N is set to twelve. They are equally distributed along the leading edge as represented in Fig 1.

The total available power for the electric engines, P E , is identical to the original power: P E = 4000KW. The only varying geometrical parameter is the vertical tail surface area S v .

The aerodynamic derivatives required to model this aircraft are determined using a Vortex Lattice Method included in OpenVSP, see [START_REF] Gloudemans | A rapid geometry modeler for conceptual aircraft[END_REF].

Propulsion modeling

Engine lines of thrust are assumed to be aligned with the aircraft X body axis and contained within the X-Y plan so that no pitching moment is induced by the thrust. Rotor terms or gyroscopic effects due to elements in rotations are assumed negligible such that equations of flight are not further modified. Thrust force and torque produced by differential thrust are calculated by summation of the contribution of each engine based on the geometrical arrangement shown in Fig 1:

F T b =      N i=1 F T,i (δ x,i ) 0 0      , M T =      0 0 N i=1 -F T,i (δ x,i )y i      , (9)
where F T (δ x,i ) is the thrust produced by the i th engine located at a distance y i from the body X-axis, with the throttle command: δ x,i . Engines are numbered starting from the outer left wing (-b/2) to the outer right wing (b/2). The term used in equation ( 4) is obtained by projecting equation ( 9) on the aerodynamic frame:

F T = T ab F T b .
(10) For electric motors, the following thrust model may be considered for estimating F T,i (δ x,i ) following Sachs (2012):

F T,i (δ x,i ) = P E N V -1 η m η p δ x,i , (11)
with η m and η p respectively the engine and propeller efficiency (both considered constant). Therefore, the power is equally divided between the N engines. Finally, δ x,i is added to the control input vector : u = [δ a , δ e , δ x,1 , . . . , δ x,N ] T . The number of inputs becomes : n u = N + 2 = 14.

Finding the trim position by optimization

A possible way to find an equilibrium is given by [START_REF] Goman | Evaluation of aircraft performance and maneuvrability by computation of attainable equilibrium sets[END_REF]. However, due to the high number of actuators present in this case, the number of variables (n x +n u = 22) is higher than the number of equations (N e = 10) and constraints. Therefore, an optimisation method is used to find a solution to this over-determined problem such as suggested by [START_REF] Oppenheimer | Control allocation for overactuated systems[END_REF].

Flight conditions are set by adding constraints to fix the following variables: [V, β, γ, Ω]. The objective function to minimise is defined as the power required to maintain equilibrium.

Finally, bounds are added on control inputs to remain within acceptable ranges of deflection or power levels.

The problem hence writes:

min x N i=1 F T,i (δ x,i )V, (12)
with: x = [α, p, q, r, φ, θ, δ a , δ e , δ x,1 , . . . , δ x,N ] T , such that:

mg + F A + F T -T ab p q r × V = 0, (13)
M A + M T - p q r × I p q r = 0, ( 14
) 0 = p + q sin φ tan θ + r cos φ tan θ, (15) 0 = q cos φ -r sin φ, (16)
Ω = (q sin φ + r cos φ) / cos θ, (17)
sin γ = cos α cos β sin θsin β sin φ cos θ sin α cos β cos φ cos θ.

(18) A Sequential Least SQuares Programming algorithm available in SciPy, see [START_REF] Kraft | A software package for sequential quadratic programming[END_REF], is used to solve the optimization problem with a tolerance set to 10 -6 on the constraints.

Linearisation

The system is linearised by a first order Taylor expansion of the system. Equations 5 and 6 form a system of non linear equations of the form:

ẋ = f (x, u), (19)
with

x = [V, β, α, p, q, r, φ, θ] T , ( 20
) u = [δ a , δ e , δ x,1 , . . . , δ x,i , . . . , δ x,12 ] T .
(21) The solution x of the optimisation problem (12) and the fixed parameters [V, β, γ, Ω] are used to reconstruct the state and input vectors: x, û with values at equilibrium. The Jacobian of f (x, u) is evaluated numerically by computing centered finite differences:

˙ x = ∂f ∂x x + ∂f ∂u u, (22)
with x = xx and u = uû.

It is of interest to introduce the flight path angle γ in the state space representation so as to better differentiate the short period oscillation (mainly α and q dynamics) and the phugoid (mainly V and γ dynamics). One can do so by a variable change: using equation ( 7), small angle assumption and neglecting terms of second order, the pitch angle can be expressed as:

θ = γ + α, thus: θ = γ + α. ( 23
)
Substracting the line corresponding to α in the matrices ∂f ∂x , ∂f ∂u completes the reformulation. Finally, the system is re-ordered and can be decoupled to obtain: ẋL

= A L x L + B L u L , (24) ẋD = A D x D + B D u D , (25)
with:

x L = [ V , γ, α, q] T , (26)
u L = [ δ e , δ x,1 , . . . , δ x,12 ] T , (27)
x D = [ β, p, r, φ] T , (28)
u D = [ δ a , δ x,1 , . . . , δ x,12 ] T , (29)
and the additional notation:

δ x = [ δ x,1 , . . . , δ x,12 ] T .

Variation of the Vertical Tail Surface Area

Capturing the effect of geometrical changes in the vertical tail surface area is a tedious problem because of the important influence of other aircraft components on the flow impacting the vertical tail. This has been demonstrated by Fabrizio Nicolosi (2013) as a motivation for the construction of a semi-empirical model (VeDSC) accounting for fuselage and horizontal tail geometry to predict vertical tail performances. For a detailed presentation of VeDSC and its inclusion in the model, the reader is referred to [START_REF] Ciliberti | A new approach in aircraft vertical tailplane design[END_REF] and Nguyen- [START_REF] Nguyen-Van | Towards an aircraft with reduced lateral static stability using electric differential thurst[END_REF].

After inclusion in the aircraft aerodynamic model, the VeDSC method allows the computation of a collection of linearised systems, each system representing a different size of the vertical tail. It further permits the investigation of the reduction of the surface area on the aircraft lateral dynamics as presented in Fig 2 . This result is obtained for the flight condition [V, β, γ, Ω] = [1.3V s = 68m/s, 0 • , 0 • , 0 • /s] which was selected as particularly sizing for differential thrust control. Indeed, at 1.3V s , the pilot requires a wide flight envelop while the control efficiency of propulsion is inversely proportional to the velocity, see equation ( 11) and Nguyen- [START_REF] Nguyen-Van | Reduction of vertical tail using differential thrust: Influence on flight control and certification[END_REF]. As the vertical tail reduces, the side slip oscillation gradually becomes unstable and changes to an aperiodic behaviour for values lower than S v /S v0 = 0.2 with two unstable real poles. In parallel, as the vertical tail becomes smaller the control effort will be increased to maintain handling quality. Therefore, an optimal trade-off between vertical tail surface area and handling quality at the cost of a reasonable control effort should exist.

Linear Fractional Representation

Using the previously described modeling, a collection of lateral linear systems [A D,j , B D,j ] for values of the vertical surface area ranging from 0.1S v,0 to 1.5S v,0 , with steps of 0.1S v,0 are computed. The collection of models is then put into a M -∆ form (or Linear Fractional Representation LFR). The toolbox APRICOT and the Matlab function lsapprox are used to this end. See Clement [START_REF] Roos | Polynomial and rational approximation with the APRICOT library of the SMAC toolbox[END_REF] for more details about the APRICOT toolbox. The degree of the polynomial chosen is 3 in order to keep the maximum root-mean-square error lower than 1% which is quite sufficient at this preliminary design phase. Then ∆ = δ v 1 6 where δ v = S v /S v0 is the tail surface ratio, varying between 0.1 to 1.5 and 1 n is the n × n identity matrix.

CO-DESIGN

Control architecture

The longitudinal/lateral flight control law is depicted in Figure 3 and involves:

• a static feed-forward gain H L (2×2), resp. H D (2×2) and a static feedback gain K L (2×2), resp. K D (2×4), for the longitudinal, resp. lateral, flight control, • a longitudinal inner loop to control the short period mode through the gains K α and K q , • the symmetrical thrust and differential thrust allocation matrices L L (6 × 1) and L D (6 × 1), respectively.

This set of gains (in green boxes in Figure 3) are the control decision variables K to be tuned by the optimization process:

K = {H L , H D , K L , K D , K α , K q , L L , L D } .
The whole allocation matrix, for the 12 throttle commands, reads:

A 12×2 = 1 6 1 6 P 6 -P 6 L L 0 6×1 0 6×1 L D ,
where the permutation matrix P n is defined by:

P n (i, j) = 1 if: i + j = n + 1, 0 otherwise. In addition, S = diag([1, 1 V , π 180 , π 180
]) is an input shaping gain to take into account the reference input on vertical velocity V z ≈ V γ since requirements specify decoupling between V and V z (see for example [START_REF] Döll | A modal multi-modal approach[END_REF]), the side slip and bank angle reference inputs β r and φ r are expressed in degrees. Thus, a unitary step on each of the 4 components of the reference input w is expected to create roughly the same thrust magnitude.

The aircraft (A/C) and avionics block of Figure 3 is detailed in Figure 4. It consists of the longitudinal and lateral models presented in the previous section completed by the avionics model between the required control signals d e , d a , and d x,i (i = 1, • • • 12) computed by the control law and the really applied actuations δ e , δ a and δ x,i . This avionics block considers a second order model on each of the 14 actuators with a damping ratio of 0.7 and a cut-off frequency of ω e , ω a and ω p on the elevator, the ailerons and the 12 propeller engines, respectively. Thus, 2 additional decision variables are added for the optimization process: the vertical tail surface ratio δ v (6 occurrences) and the propeller engine bandwidth ω p (24 occurrences).

w V r V V zr α q β p φ r γ φ r β r [K α K q ] S + - + - + - + - + - P 6 + - + + d a d e H L H D K L L L L D K D d x,1 . . . d x,6 S -1 d x,7 . . .
1 s 1 4 1 s 1 4 δ v 1 6 [A L B L ] M x L ẋL u L x D ẋD u D V α q β p φ r γ d a d e d x,1
. . .

Requirements and objectives

The handling qualities are expressed through frequencydomain templates S o,des (i, j) (i, j = 1, 2, 3, 4) on each element of the 4 × 4 output sensitivity function S o (s) = 1 4 -T w→z (s) where T w→z (s) is the closed-loop transfer between w and z as depicted in Figure 3:

γ 1 = max i,j=1,2,3,4 S o (i, j) S o,des (i, j) (s) ∞ ≤ 1, (30)
with:

S o,des (i, j) = 1.4 s s + ω i , if: i = j,
0.14, otherwise.

and :

ω 1 = 0.2, ω 2 = 0.3, ω 3 = ω 4 = 1 (rd/s
) . Such a multivariate template allows to prescribe lowfrequency disturbance rejection and reference input tracking inside the bandwidths ω 1 ,.. ω 4 for the servo-loops on V , V z , β and φ, respectively, while ensuring cross couplings to be lower than 14%.

To avoid too high gains in the feed-forward path of the control law and so to reduce as much as possible the thrust magnitude in response to reference inputs, the decision variables must be tuned in order to minimize the maximum of the 12 closed-loop propeller transfers T w→ dx,i :

J 1 = max i=1,•••12 T w→ dx,i (s) ∞ .
(31)

Minimizing the objective function J 1 defined in (31) under the constraint defined in ( 30) is now common in the field of H ∞ structured control design thanks to non-smooth optimization tools [START_REF] Apkarian | Tuning Controllers Against Multiple Design Requirements[END_REF]). The challenge here is to minimize also the sizing parameters δ v and ω p using 2 additional objective functions: J 2 = δ v and J 3 = ω p .

(32)

Multi-step optimization process

The multi-step process described below splits the initial multi-objective optimization problem into several consecutive single objective optimization problems where the previous reached objectives are constrained to stay inside a given sub-optimal solution set for the next optimization step. That allows to better understand the optimization process and to manage the trade-off between the various objective functions which is particularly useful during the preliminary design phase. Furthermore, such an approach is also justified since the various objectives J 1 , J 2 and J 3 have not normalized units.

Step 0: Initialization. The decision variables are initialized in the following way:

• the initial A/C configuration is defined by δ v = 0.8 and the avionics are assumed perfect, • the allocation matrix A is initialized with: L L = [0, 0, 0, 0, 0, 1] T , L D = [1, 0, 0, 0, 0, 0] T in such a way that symmetrical thrust uses only inner propellers 6 and 7 and differential thrust uses only outer propellers 1 and 12, • then, K α and K q are designed to assign the shortperiod mode to the roots of s 1.2 2 + 1.4 s 1.2 + 1 on the 2-nd order short-term longitudinal model (also called (αq) model),

• K L is designed to assign the eigenstructure of the longitudinal long-term model (also called (Vγ) model), i.e. place the 2 corresponding eigenvalues to -0.2 and -0.3 while decoupling the associated eigenvectors from γ or V respectively. See more details in [START_REF] Döll | A modal multi-modal approach[END_REF]. H L is then computed to have an identity DC-gain between [ V r , γ r ] T and [ V , γ] T on this (Vγ) model, • K D is designed to assign the 4 eigenvalues of the lateral model to -1±j, resp. -1, -20, and the associated eigenvectors decoupled from φ, resp. β. See more details in [START_REF] Döll | A modal multi-modal approach[END_REF]. H D is then computed to have an identity DC-gain between [ β r , φ r ] T and [ β, φ] T , • Finally, the actuators are initialized with ω p = ω e = ω a = 20 rd/s.

Step 1: Optimization on the initial A/C configuration. This optimization aims to meet the constraint (30) while minimizing J 1 using the control decision variables K:

K = arg min K J 1 , s.t. (30) holds.
Let us denote J i,j the optimal value of J i obtained at step j. Thus J 1,1 = min K J 1 .

Step 2: Vertical tail minimization. This step aims to minimize the vertical tail surface J 2 while meeting the constraint (30) and keep the objective function J 1 lower than a sub-optimal value J1 > J 1,1 : 30) holds and J 1 ≤ J1 .

{ K, δ v } = arg min K,δv J 2 , s.t. (
Step 3: Engine bandwidth minimization. This step aims to minimize the engine bandwidth J 3 while meeting the constraints of step 2 and keep the objective function J 2 lower than a sub-optimal value J2 > J 2,2 : 30) holds and

{ K, δ v , ω p } = arg min K,δv,ωp J 3 , s.t. (
J 1 ≤ J1 , J 2 ≤ J2 .
J1 and J2 , respectively chosen after step 1 and step 2, can be seen as sub-optimality tolerances required to relax the optimality obtained at the previous step in order to add a new objective at the current step.

RESULTS

In this section, the results presented are obtained with J1 = 0.2 and J2 = 0.4. These values are justified by the big gap between J 1,0 = 0.713 and J 1,1 = 0.130 (resp. J 2,1 = 0.8 and J 2,2 = 0.1) presented in Table 2. This Table presents the evolution of the hard constraint γ 1 and the 3 objective functions J 1 , J 2 and J 3 after each step of the optimization process. It is worth to mention that the initialization does not meet the hard constraint when the full model and avionics are taken into account. As expected, a trade-off between J 2 and J 3 has to be managed. The obtained solution allows the vertical tail surface to be reduced in a significant way ( δ v = 0.4) while using a low propeller engine bandwidth ω p = 6.74 rd/s. The obtained optimal allocation matrix A is represented on the bar-diagram of Figure 5. As expected, symmetrical thrust uses mainly inner engines while differential thrust takes benefit of the lever-arm effect and uses mainly outer ones. In addition, Figure 6 displays the bar-diagram of the direct feed-through from the 4 reference inputs (vector w) to the 12 throttle commands d x,i :

A H L (2, :) 0 1×2 0 1×2 H D (2, :) S. Obviously, the most demanding maneuvers from the total thrust magnitude point of view are the forward velocity V r and the side-slip angle β r .

CONCLUSION

An aircraft model using distributed electrical propulsion parameterized according to the vertical tail surface and propeller engine bandwidth was developed and used to codesign the full axes control law of a reduced-tail/without rudder configuration using low actuator bandwidths. The optimization process allows main trade-offs to be managed during the preliminary design phase. Further developments to consolidate such an approach will be focused on:

• its application to the full flight envelope,

• its experimental validation thanks to an aircraft mock-up using distributed and differential propulsion, currently under progress.

Fig. 1 .

 1 Fig. 1. Illustration of thrust distribution on the reference aircraft with twelve electric engines.

Fig. 2 .

 2 Fig. 2. Evolution of the eigenvalues of the side slip oscillation of the reference aircraft with the reduction of the vertical tail surface area, from S v /S v0 = 1 to S v /S v0 = 0.1. The number in italic next to an eigenvalue indicates the corresponding quantity S v /S v0 . The regular number next to a dashed line corresponds to the damping value.

Fig. 3 .

 3 Fig. 3. The longitudinal/lateral closed-loop control blockdiagram.

Fig. 4 .

 4 Fig. 4. The A/C + avionics block-diagram.

Fig. 5 .

 5 Fig. 5. Bar-diagram of the allocation matrix A.

Fig. 6 .

 6 Fig. 6. Bar-diagram of the direct feedthrough from reference input to throttle commands.

Table 1

 1

	. ATR 72 general details. See ATR
	(2015) and Jackson (2015)
	Variables	Value
	Wingspan (b)	27m
	Wing surface area (S)	61m 2
	Vertical tail surface area (Sv 0 )	12 m 2
	Mass (m)	22800Kg
	Total available power	4000KW
	Stall velocity Vs	56m/s

Table 2 .

 2 Constraint γ 1 and objective function values during optimization process.

	step	γ 1	J 1	J 2 = δv	J 3 = ωp (rd/s)
	0	1.107 0.713	0.8	20
	1	0.999 0.130	0.8	20
	2	0.984 0.188	0.1	20
	3	0.999 0.199	0.4	6.73