Clément Michel

Janette Cardoso

Pierre Siron

TIME MANAGEMENT OF HETEROGENEOUS DISTRIBUTED SIMULATION

Keywords: Aerospace, Distributed Processors, Model design, Discrete simulation, Simulation interfaces, Cyber-Physical Systems, Heterogeneous Systems, HLA, Distributed Simulation

Cyber-physical systems (CPS), by their very nature, mix continuous and discrete behavior and are modeled by heterogeneous components. Formal analysis cannot always handle such complex systems and simulation is a necessary step. In particular, distributed simulation is very useful for the validation of CPS for two main reasons: either the CPS itself is distributed (e.g., a fleet of UAVs) or the CPS is too complex and/or has too much models (e.g., an aircraft). We discuss in this paper the impact of distributing the simulation of a system: which are the rules that must be applied to guarantee a correct behavior between the different simulators? If a centralized simulation already exists (using Discrete Event simulation), which hypothesis must be made for the Distributed Discrete Event simulation? The co-simulation framework used and discussed in this work is Ptolemy-HLA. It allows a Ptolemy model to be distributed using the high-level architecture (HLA) standard.

INTRODUCTION

The analysis of cyber-physical systems (CPS) is a complex task due to the heterogeneity of the parts involved, as they integrate different methodologies and tools. Simulating different parts of a CPS requires different abstraction and tool supports, and the lack of interoperability between tools poses a major challenge. Because of the nature of the CPS, or because of its complexity, distribution can be necessary. Distributing a simulation brings its own challenges, as it requires all the simulation elements to conform to a collection of rules in order for the elements to communicate between them.

In this work, the Ptolemy framework is used for mod-eling the CPS system, and HLA is the standard chosen for the distribution. The IEEE High-Level Architecture (HLA) is a standard for distributed discreteevent (DDE) simulation. CERTI is a HLA-compliant RTI (Run Time Infrastructure). Ptolemy II is a modeling and simulation tool for heterogeneous systems, well suited for modeling CPS since it provides different models of computation (MoC). In this paper, the simulation entities, called federates are Ptolemy models, but the approach presented in this paper can be easily used for other simulators. We focus on the time representation issues introduced in a distributed simulation, using HLA as a standard for the co-simulation and Ptolemy as a simulator component. Interoperability and reuse are important targets, so the coupling between simulators is an important issue in the design. The first step in a distributed simulation is to find the partition of a simulation model consisting of a set of sub-models. Which sub-models can (or must) be put together in a simulator belonging to the distributed simulation? How to guarantee that the distributed simulation will be valid (according to some criteria)? The purpose of this work is twofold: Providing properties that ensure a correct time coordination between a federate and HLA, and studying the impact of the distribution in the Ptolemy-HLA framework in order to build models requiring timely input. By studying this specific coupling, we intend to find general problems as well their solutions concerning the time management of heterogeneous distributed simulations. We will start by presenting the characteristics of both Ptolemy and HLA. Then, we will discuss the impact of an HLA-based distribution, before introducing elements needed for the Ptolemy-HLA simulation to be conservative.

PTOLEMY

Ptolemy II is a Java open-source simulation and modeling tool intended for experimenting with system design techniques, particularly those that involve combinations of different types of models (Ptolemaeus 2014). A Ptolemy simulation, called model, is com-posed of a Director and software components called actors, that execute concurrently and communicate through messages (called events) carrying values, sent via interconnected ports. An actor that is executed is said to be fired.

The collection of rules that governs concurrent execution of the actors and the communication between them is called a Model of Computation (MoC). The MoC for each actor is implemented as a Director, a software component that dictates how actors should be fired. In this paper, we focus on the Discrete Event (DE) and Continuous MoCs. Ptolemy uses a model of time known as the superdense time, represented by a tuple (t, n), where t is called the model time and n is called the microstep. The model time represents the time at which some events occurs, and the microstep represents the sequencing of events that occur at the same model time [START_REF] Manna | Verifying Hybrid Systems[END_REF]Pnueli 1993, Ptolemaeus 2014). The events are ordered in the event queue first by t, then by n.

An event e is noted e((t, n), v) with (t, n) the timestamp and v the value of the event. To ensure determinism, the order in which actors are fired when multiple events are queued for the same timestamp (t, n) is given by the actor rank, a topological sort that lists the actors in data-precedence order. In the DE MoC, a model advances its logical time to the timestamp of the next event in the queue, and the actor to whom this event is destined for is executed. In the Continuous MoC, the model advances its logical time to a timestamp computed by the solver (Runge-Kutta 23 or Runge-Kutta 45), and all actors are fired at once. Ptolemy provides a so-called TimeRegulator interface with a proposeTime method (Ptolemaeus 2014). This interface is implemented by attributes that wish to be consulted when a Director advances time. The Director will call the proposeTime method, passing it a proposed time to advance to, and the method will return either the same proposed time or a smaller time.

HLA

The High-Level Architecture (HLA) is a standard for distributed discrete-event simulation. In HLA terminology, the entire system to be simulated is represented by a federation, which is a collection of federates (simulation entities or simulators). The execution of a federation uses a middleware, called the Runtime Infrastructure (RTI). The federation have a Federation Object Model (FOM), a file that contains the definition of data structures (called objects) exchanged between the federates. In a federation, messages are exchanged between the various federates through the RTI. Those messages, called events, can be timestamped. Messages that are time-stamped are said TSO for Time Stamp Order.

Among the services defined by the HLA standard (IEEE-SA Standards Board 2010), we will focus on object management and on time management.

The object management service allows for federates to exchange messages and values. We focus here on two functions: Update attribute value (UAV), that sends a message to the federation, and Reflect Attribute Value (RAV), that delivers a message from the federation. A federate is said regulating and constrained when it both sends and receives TSO messages. When all the federates are regulating and constrained, the federation is said to be conservatively synchronized [START_REF] Kuhl | Creating Computer Simulation Systems: An Introduction to the High Level Architecture[END_REF]. Regulating federates generate TSO messages that must occur no earlier than h c + lah, with h c being the current HLA time for the federate, and lah being its lookahead, a value that establishes a lower bound on the timestamps that can be sent. While a Discrete Event simulation simply advances its time to t when wanting to, a Distributed Discrete Event simulation first asks for the permission to advance to h, and only advances to h when granted its authorization. A time constrained federate requests to move its logical time forward by first asking the RTI to do so, either through a Next Event Request (NER) or a Time Advance Request (TAR). The RTI replies to this request by sending all TSO messages up to h to the federate, then sending back a Time Advance Grant (TAG). The TAG is written T AG(h) and is an "authorization" for the federate to advance to the logical time h. TAR Consider a federate is at time h and asks for a time h 1 through a T AR(h 1 ); the next time value is always h 1 granted by T AG(h 1 ). So T AR(h 1 ) → T AG(h 1 ). The figure 1 illustrates the reception of a RAV (h ) messages with h < h ≤ h 1 during the advancing phase to h 1 : the logical time is not updated to h and is eventually granted to h 1 . NER Consider a federate is at time h and asks for a time h 1 through a N ER(h 1 ). If it receives a message RAV (h ) with h < h ≤ h 1 , it advances its time to h :

N ER(h 1 ) → T AG(h ), if RAV (h ) T AG(h 1 ), if no RAV (h )
Figure 2 illustrates this advance policy. If the federate is still willing to go to h 1 , a new N ER(h 1 ) must be asked.

IMPACT OF THE DISTRIBUTION AND COUPLING

To take part in a HLA federation, a simulator needs a HLA coupling interface, that handles the coupling between the simulator and HLA. Let us consider Figure 3 where two parts are depicted:

• Distribution: ruled by HLA standard that guarantees the time advancing is coherent for all federates (no simulator will advance to a time in the past, and TSO message are delivered in time-stamp order);

• Coupling: put in accordance the time advancing mechanism of the simulator with the HLA time advancing mechanism and so guarantee an ordered data exchange between the simulators in the federation.

The object management services (UAV, RAV) and the time management services (NER, TAR) are independent. However, the delivery of a RAV message is done during the time advancing phase. Let be t be the simulation time and h the HLA time.

The distributed simulation in Figure 3 is a Federation where Federate f1 simulates a model M1 and Federate f2 simulates a model M2. Each federate has its own calendar queue, and events are sent/received through the RTI using HLA services UAV and RAV. The time is advanced using HLA services TAR (or NER).

Figure 3: Distribution and coupling Let be:

• t 1 the timestamp at which f 1 wants to produce an event (to be sent through the RTI);

• h e the timestamp of this corresponding event in the HLA service UAV sent through the RTI; it is also the timestamp of the callback RAV;

• t 2 the timestamp of an event put in f2 calendar queue after the reception of the RAV from the RTI;

• h c the current HLA time;

• lah the lookahead;

• T S the HLA time step, defined only when using TAR; h c + T S is called next point in time.

In such a distributed simulation, the event time stamp of a message can be altered twice, t 1 → h e and h e → t 2 , due to:

• Different time representations;

• the distribution itself and its rules, e.g. a federate cannot send an event earlier then h c + lookahead;

• the coupling HLA ↔ simulator providing a coherent way for both time advancing mechanisms.

Indeed, a coupled distributed simulation does not come without a cost, and t 2 > t 1 . In the next section, the Ptolemy-HLA coupling will be presented. Let us point out that in this paper we will focus on NER time management but TAR was also designed and implemented.

PTOLEMY-HLA

Ptolemy-HLA [START_REF] Lasnier | Distributed Simulation of Heterogeneous and Real-Time Systems[END_REF]) aims to provide a distributed, conservative simulation. Both HLA and Ptolemy offer conservative time management options [START_REF] Buck | Ptolemy: A Framework for Simulating and Prototyping Heterogeneous Systems[END_REF][START_REF] Fujimoto | Parallel Simulation: Distributed Simulation Systems[END_REF]), then the interface binding them together must be conservative as well.

In this section, we detail the different conditions needed for a Ptolemy-HLA conservative simulation.

Ptolemy/HLA Coupling Design

For a Ptolemy model to work as a HLA federate, the model's Director must be a DE Director. However, it can very well contain composite actors that themselves contain a Continuous Director each. Three Ptolemy components are added for Ptolemy and HLA to communicate:

• HlaManager : encompasses the high-level HLA rules and services such as time management and object management.

• HlaSubscriber actor: receives events from the HLA federation through RAVs (pictured on Figure 4b).

• HlaPublisher actor: sends Ptolemy events to the HLA federation through UAV service (pictured on Figure 4c).

Events circulating through the federation (even when not sent to another federate) must respect both HLA and Ptolemy time restrictions. In order to coordinate both Ptolemy time advancing and HLA time advancing, the proposeTime method was extended in order to allow Ptolemy to query the RTI for a time t. When a Ptolemy federate wants to advance its time to the timestamp of the earliest event available, it first interrogates the RTI through the proposeTime algorithm, and wait for the RTI to grant it. Algorithm 1 displays the proposeTime for the NER.

Algorithm 1 Extended proposeTime algorithm (NER)

1: NER(t asked ) 2: while T AG(h received ) not received do 3:

tick() 4: end while 5: t asked ← h received 6: if RAV received then 7:

Schedule HlaSubscriber firing at t asked 8: end if 9: return t asked

The Ptolemy-HLA framework allows a federation (f1, f2) with any combination of federates time management: (NER,NER), (NER,TAR), (TAR, NER), (TAR,TAR). This is true also for a federation with more than two federates. The events exchanged between the federates see their timestamp manipulated and shifted as depicted on Figure 3. In order to coordinate HLA and Ptolemy time, the coupling interface in the simulator has two time lines. In this section, we consider that both time lines have the same computer representation, so both can be directly compared. In the sequel, let us present how the timestamp can change in the Ptolemy-HLA framework when sending and receiving an event according to the time management.

Sending a Ptolemy event through the RTI: Let be t 1c the current (Ptolemy) logical time of federate f1 and h 1c be the current HLA time. The federate wants to send an event e(t 1c ) (through the HlaPublisher actor) using a HLA UAV(h 1 ) service; timestamp h 1 depends on the federate's time management:

NER: h 1 = h 1c + lah, with h 1c = t 1c (1a)
TAR:

h 1 = h 1c + lah, if t 1c < h 1c + lah t 1c , otherwise (1b)
Receiving a RAV from the RTI: Let be t 2c the current (Ptolemy) logical time of federate f2 and h 2c be the current HLA time. The reception of the HLA RAV(h 1 ) service at h 2c wakes up a HlaSubscriber actor with a Ptolemy event e(t 2 ); the timestamp t 2 depends on the federate's time management:

NER: t 2 =h 1 (2a) TAR: t 2 =h 2c + T S (2b)
These rules are necessary in order to produce a conservative distributed simulation that respects the rules of both HLA and Ptolemy (i.e. no TSO message sent earlier than h c + lookahead, no event insertion in Ptolemy's past).

It can be seen from equations 1a, 1b, 2a and 2b that the difference t 2 -t 1c between the production of an event at time t 1c in f1 and its consumption at time t 2 at f2 depends on the time management of each federate and its parameters:

t 2 -t 1c = f (lah, (T S), t 1c , h 1c , h 2c ).
Let us consider again Figure 3, where both federates f1 and f2 use a NER time management with h c current HLA time for both federates and a same lookahead lah.

Notice that they can also have different lookaheads. Let be t 1c the current (Ptolemy) logical time and e f1 (t 1c ) the event that wakes up the HlaPublisher actor in f1. By using equation 1a, the corresponding UAV is sent with timestamp h 1 = t 1c + lah (regardless of f2 time management). The RAV event with timestamp h 1 received by f2 is queued at HlaSubscriber actor as a event e f2 with timestamp t 2 = h 1 given by equation 2a. Thus, the (NER,NER) configuration introduces in the distributed simulation a delay δ N ER-N ER = t 2 -t 1c = lah.

Ptolemy/HLA coupling: implementation issues

Coupling Ptolemy with HLA requires to take into account that the time representation is different: CERTI RTI uses IEEE 754 double (with dynamic precision) and Ptolemy uses its own Time representation (fixed precision). These two time representations must be well converted. Timestamps in Ptolemy are represented under the form t = n * r with n an integer (called the time value) and r a Java double (called the time resolution). Thus, time values in Ptolemy can only adopt values that are multiple of r, e.g. r, 2r, 3r, etc. . . Despite the time representation difference, one must guarantee that:

• Rule 1: No event should be inserted into the simulator's pastt RAV > t P tII

• Rule 2: No UAV should be sent in the past of the HLA federationh U AV > h HLA Non observance of these rules violates the causality of either the federation or the simulator, outputting wrong results.

At the earliest stages of Ptolemy-HLA, the question of the conversion between different (computer) representation was eluded, even if it was well-known that a function f such as f (h HLA ) = t P tII and f -1 (t P tII ) = h HLA was unobtainable.

In the following, the difference between HLA and Ptolemy (or any other simulator that does not use double) time representation is taken into account. In this work, we introduce two functions f and g such as:

f (h HLA ) = t P tII g(t P tII ) = h HLA Several properties are required from f and g in order to respect temporal coherence in TAR and NER cases, such as:

• f and g are monotonous strictly increasing functions

• ∀h 1 , h 2 , h 3 h 3 ≥ h 2 > h 1 ⇒ h 3 ≥ (g • f )(h 2 ) > h 1 • ∀t 1 , t 2 t 2 ≥ t 1 ⇒ (f • g)(t 2 ) ≥ t 1 .
We find these properties hard to fulfill from a mathematical point of view, even probably they cannot be fulfilled if f and g are not the identity function and this is not possible with heterogeneous systems. In such a context, we can find and use algorithmic solutions in order to use existing conversion functions and to solve the list of the identified problems. Thus, each time that t and h must be compared in Algorithm 1, we introduce f and g, two monotonous, nonstrictly increasing functions, with the hypothesis that algorithmic modifications will allow the functions to work in our context. Considering the f and g functions, Algorithm 1 is modified to Algorithm 2 in the following way:

• NER is called only if g(t asked ) > h current : prevents Ptolemy from requesting a time advance if the requested time is imperceptible by the HLA time.

• t asked ← t current +r, if f (h received ) ≤ t current , r being Ptolemy's resolution: prevents HLA from granting a time advance imperceptible by the Ptolemy time.

Equations 1a, 1b, 2a and 2b are also modified in order to take f and g into account. For instance, 1a becomes h 1 =, h 1c + lah, with h 1c = g(t 1c ). Finally, we keep the association between the return of a function and its argument, i.e. when a value t 1 = f (h 1 ) is received, the association between t 1 and h 1 is memorized, and we will output h 1 when evaluating g(t 1 ).

Algorithm 2 ensures a framework where rules 1 and 2 are observed.

Because of this difference in HLA and Ptolemy time representation, the theoretical delay induced by the distribution described in the previous sub-section can be slightly altered by a value introduced by the conversion between the time representations.

APPLICATION

A full case study of a longitudinal flight control [START_REF] Lasnier | Distributed Simulation of Heterogeneous and Real-Time Systems[END_REF]) will be considered in this section.

Algorithm 2 Modified ProposeTime Algorithm (NER)

1: if g(t asked ) > h current then 2:
NER(g(t asked ))

3:

while T AG(h received ) not received do 4:

tick() The block comports one output ElevCom that is the elevation command sent to the Aircraft block, and takes in three inputs:

• stickS: The current action on the stick.

• alphaS: The current vertical velocity of the F-14.

• qS: The current pitch rate of the F-14. If only one of the inputs receives an event, the AutoPilotDE block will consider the others inputs as absent. Thus, receiving the events at three different timestamps would trigger three different rounds of computation that would output incorrect elevCom values. As a consequence, the events concerning variables stickS, alphaS and qS must be fed at the same timestamp (t, n) for the output to be correct. Thanks to the topological sort, we are guaranteed to have all these events produced by the PeriodicSampler before actors within AutoPilotDE are fired, processing the three events at once. A distributed version of the F-14 obtained by creating a new federate model for each of the three acin Figure 4a: PilotStick, AutoPilotDE, and PilotStick. For each actor, an input port is connected to an HlaSubscriber (pictured on Figure 4b) and an port is connected to an HlaPublisher (pictured on Figure 4c). The f14AutoPilotDE federate, containing the AutoPilotDE block, is pictured on Figure 5. The Stick PilotStick and AircraftSAct are composite actors containing the HlaSubscriber receiving events, respectively, from f14PilotStick and f14Aircraft federates, while the elevCom actor is an HlaPublisher sending events to the federation. The three federates use the NER time advancing mechanism with a same lookahead lah. Let us consider t stick , t alpha , t q the timestamps of the events arriving, respectively, at inputs stickS, alphaS (generated by the f14Aircraft federate) and qS (generated by the f14PilotStick federate). According to equation 1a, the UAVs at f14Aircraft and f14PilotStick federate are sent with timestamp h = g(t + lah), with h the current HLA time for the federate when the event is produced. The RAVs received at f14AutoPilotDE federate are then queued as events e with timestamps t = f (h) according to equation 2a. For the AutoPilotDE block inside f14AutoPilotDE federate to work properly, events for all the inputs are required to be received at the same time. Thus, t stick = t alpha = t q . By expressing t as a function of t, we obtain the following condition:

Distributed simulation

f (g(t stick ) + lah) = f (g(t alpha ) + lah) = f (g(t q ) + lah)
Since the three federates are Ptolemy federates, they have the same f, g functions. Moreover, as all the federates have the same lah, we can reduce the condition to t stick = t alpha = t q This condition is satisfied as the PeriodicSamplers (actors sampling the output of continuous models in Aircraft and Stick blocks in Figure 4a) outputs all the events at the same time since they have the same period. Thus, the AutoPilotDE block receives and processes its inputs the same way than in the centralized model, albeit in a time-shifted manner. The bias lah introduced by the distribution does not, in this case, generate significant changes in the simulation behavior. So, the AutoPilotDE block maintains its expected behavior under some conditions on the other federates.

RELATED WORK

In [START_REF] Deschamps | Toward a Formalism to Study the Scheduling of Cyber-Physical Systems Simulations[END_REF], the partitioning of a simulation is discussed, addressing the simultaneity at the inputs of all components (of a distributed simulation).

A formalism is proposed to prove the equality of delays by design and so the distributed simulation is valid.

The case study of an f14 aircraft is presented in [START_REF] Michel | Distributed Simulation of Cyber-Physical Systems[END_REF] where, for a given partition, the cyber component is analyzed to guarantee that the data consumptions are simultaneous using different implementations: memory blocs, internal clock. An aspect in distributed simulation is the multiplicity of times and the different representation of these time values. One must guarantee that the times advancing of the distributed simulation are done correctly (e.g., conservation simulation [START_REF] Fujimoto | Parallel Simulation: Distributed Simulation Systems[END_REF]). The time representation of each simulator, and the entity that allows for co-simulation -HLA or master algorithm in FMImust be well dealt in the coupling between them. A detailed analysis of time representation in the FMI framework is done in [START_REF] Cremona | Hybrid Co-Simulation: It's About Time[END_REF]). In particular this paper discuss the choice of resolution to be used when the FMUs (components of a co-simulation) have different resolutions. The coordination of different times issue appears also in real cyber-physical systems [START_REF] Shrivastava | Time in Cyber-Physical Systems[END_REF].

In [START_REF] Nägele | Co-Simulation of Cyber-Physical Systems Using HLA[END_REF], a HLA/simulator coupling, called wrapper, is modeled using POOSL (Parallel Object-Oriented Specification Language). They use PoRTIco, a HLA compliant RTI and their federates are also regulating and constrained.

CONCLUSION AND PERSPECTIVES

Ptolemy-HLA framework allows to run valid distributed simulations for Event-Driven (NER) and Time-Stepped (TAR) advancing mechanisms. For a same federate model, the user can choose several parameters such as the time advancing mechanism (TAR or NER) and the lookahead. The framework implementation as well several demos are available at (The Ptolemy Project 2017). The coupling algorithm between Ptolemy and HLA in the case of NER mechanism was detailed in this paper.

For highlighting a valid coupling, first we consider that the (computer) representation of its timelines are the same, then we extended the algorithm for taking their representation difference into account. The TAR mechanism was also implemented in (The Ptolemy Project 2017). In both cases a time bias between the federates, expressed by equations 1a to 2b, is introduced but the algorithm guarantees an valid and efficient distribution.

When considering different time representation, the bias can be slightly increased by a value introduced by the time conversion. The f14 Federation presented in this paper was obtained from the centralized model following the characteristics of a CPS: a federate for the cyber part (the controller), a federate for the plant (the aircraft), and a federate for the pilot stick that can be replaced by a real hardware in this federation. However, the distribution of a model needs some thought process about the model. In a centralized model some hypotheses are implicit, for instance the reception of all data generated by uphill actors. This constraint is ensured by the topological sort handled by Ptolemy that determines a deterministic execution order. This execution order can also be ensured in a distributed context, as well as the (logical) simultaneity of the inbound events of an actor, even considering the timestamp modification introduced by HLA. An analysis must be performed by the user to ensure that the chosen distribution is correct since a simulation can be distributed according to different mappings.

For helping the user to make this analysis, the next step is the design of an extra "layer" that both helps the user to distribute a model and ensure that the model is correctly distributed taking into account the bias introducing by the distribution and the coupling.

Figure 1 :

 1 Figure 1: TAR Time Advance Policy Illustrated.

Figure 2 :

 2 Figure 2: NER Time Advance Policy Illustrated.

Figure 4 :

 4 Figure 4: F-14

Figure 5 :

 5 Figure 5: Federate f14AutoPilotDE.

ACKNOWLEDGMENT

This research was partly supported by the French Ministry of Defense through financial support of the Direction Générale de l'Armement.


Franc ¸ois Defaÿ José A Macés-Hernández

Corentin Chauffaut

Autonomous landing of an UAV on a moving platform using Model Predictive Control

This paper proposes a linear Model Predictive Control (MPC) for high-accuracy tracking of a mobile robotic platform, in an indoor environment. The mission is divided into three main phases: target detection, target tracking and autonomous lading. These phases were coded into a so-called guidance function, which allows the system to safely switch between different reference signals according to the state of the mission. In order to ease the control loop implementation in the experimental framework available at the facilities of ISAE Supaero, some existing optimisation tools were employed. In addition, some simulations and real tests were carried out to demonstrate the performance of the controller against other control techniques developed for the same experimental system.

I. INTRODUCTION

Quadrotors have been widely used in the automatic control research field, these systems are ideal to test control algorithms given its non-linear dynamics, size and relative low price; but not less important, they have great number of applications in different fields; such as remote inspection, monitoring, aerial mapping and target recognition.

In [START_REF] Lecointe | Backstepping control law application to path tracking with an indoor quadrotor[END_REF], a control loop based on the Backstepping approach was formulated and implemented for the same UAV and experimental framework. Regarding the application of the reference tracking problem using visual feedback, [START_REF] Ajmera | Autonomous visual tracking and landing of a quadrotor on a moving platform[END_REF], [START_REF] Ling | Autonomous Maritime Landings for Low-Cost VTOL Aerial Vehicles[END_REF], and [START_REF] Wenzel | Automatic Take Off, Tracking and Landing of a Miniature UAV on a Moving Carrier Vehicle[END_REF], dealt with the problem using PIDs, [START_REF] Lee | Autonomous landing of a VTOL UAV on a moving platform using image-based visual servoing[END_REF] developed a Backstepping control and [START_REF] Daly | Coordinated landing of a quadrotor on a skid-steered ground vehicle in the presence of time delays[END_REF] proposed a time delay robust decentralised controller. Some researchers have also suggested a solution based on MPC; for instance, [START_REF] Raffo | An Integral Predictive/Nonlinear H∞ Control Structure for a Quadrotor Helicopter[END_REF] proposes an integral predictive control for the path tracking problem, [START_REF] Bouffard | On-board Model Predictive Control of a Quadrotor Helicopter: Design, Implementation, and Experiments[END_REF] and [START_REF] Mueller | A model predictive controller for quadrocopter state interception[END_REF] address the solution using MPC for position stabilisation. [START_REF] Bangura | Real-time Model Predictive Control for Quadrotors[END_REF] presents a linear model predictive control and provides some experimental results.

Concerning the target recognition, a tracking algorithm for ground robots was developed and successfully tested by [START_REF] Ren | Nonlinear tracking control for non-holonomic mobile robots with input constraints: an experimental study[END_REF]. This algorithm was developed for ground robots, but it can easily be applied to a moving platform like the one used in this investigation. On the other hand, in [START_REF] Teulire | Chasing a moving target from a flying UAV[END_REF], [START_REF] Ling | Precision Landing of a Quadrotor UAV on a Moving Target Using Low-Cost Sensors[END_REF] and [START_REF] Ajmera | Autonomous visual tracking and landing of a quadrotor on a moving platform[END_REF], an on-board target detection algorithm is proposed, obtaining successful results in all cases.

II. SYSTEM MODELLING

An inertial frame F i = {O i , X i , Y i , Z i } and a body frame F b = {O b , X b , Y b , Z b } are defined to fully represent the attitude and position dynamics of the quadrotor. The origin of the inertial frame can be placed at any point on José A. Macés-Hernández, Franc ¸ois Defaÿ and Corentin Chauffaut are with the Institut Supérieur de l'Aéronautique et de l'Espace (ISAE), 31055 Toulouse, France (e-mails: jose-alfredo.maces-hernandez@student.isae-supaero.fr, francois.defay@isae.fr and corentin.chauffaut@isae.fr). the Earth's surface, while the body frame origin matches with the quadrotor's center of mass (see Figure 1). X i -and Y i -axes point towards Earth's geographical north and east, respectively; while the Z i -axis points in the opposite direction of the gravity vector. In addition, the body frame X b -axis crosses between the motors M 1 and M 2 , and Y b -axis is placed in the middle of motors M 2 and M 3 , Z b -axis completes an orthonormal system. Given these reference frames, vectors

i ξ =   x y z   and i ξ =   ẋ ẏ ż  , (1)
will be used to represent the UAV's center of gravity position and velocity with respect to the origin of the inertial frame; respectively. Vector

b V =   u v w   = b R i i ξ = = c(ψ)c(θ) s(ψ)c(θ) -s(θ) c(ψ)s(θ)s(φ) -s(ψ)c(φ) s(ψ)s(θ)s(φ) + c(ψ)c(φ) c(θ)s(φ) c(ψ)s(θ)c(φ) + s(ψ)s(φ) s(ψ)s(θ)c(φ) -c(ψ)s(φ) c(θ)c(φ)   ẋ ẏ ż  , (2)
represents the linear velocity directly measured on the body frame.

Regarding the quadrotor's attitude, vectors

i η =   φ θ ψ   and i η =   φ θ ψ  , (3)
provide information about the angles and angular velocities of the body axis with respect to the inertial frame, and

b Ω =   p q r   = W i η =   1 0 -sin(θ) 0 cos(φ) cos(θ) sin(φ) 0 -sin(φ) cos(φ)cos(θ)     φ θ ψ  , (4)
are the angular rates measured on the body frame.

As seen in Figure 1, the UAV is composed by a cross structure with a motor in the extremities, usually the structure is built of light material. Motors M 1 and M 3 rotate in clockwise direction, while M 2 and M 4 rotate in anticlockwise; this configuration induces opposite torques in the XY b plane. The rotation of the motors produces a punctual effect in the UAV dynamics, which can either be a thrust force

f i = k f ω i 2 , (5)
or a torque

τ i = k τ ω i 2 . ( 6
)
The combination of these small effects in F b is translated into a total thrust

u = f 1 + f 2 + f 3 + f 4 , (7)
a torque in roll

τ φ = l(f 1 -f 2 -f 3 + f 4 ), (8)
a torque in pitch

τ θ = l(f 1 + f 2 -f 3 -f 4 ), (9)
and a torque in yaw

τ ψ = τ 1 -τ 2 + τ 3 -τ 4 . (10)
Therefore, a positive roll movement is produced if ω 2 + ω 3 > ω 1 + ω 4 ; in the same way, when ω 1 + ω 2 > ω 3 + ω 4 , the UAV tilts positively in pitch; and a positive yaw results from ω 1 + ω 3 > ω 2 + ω 4 .

Those equations can be rewritten in matrix form as

    u τ φ τ θ τ ψ     =     1 1 1 1 k a -k a -k a k a k a k a -k a -k a k c -k c k c -k c         f 1 f 2 f 3 f 4     , (11)
where k c = kτ k f and k a = l, k τ and k f are called torque and force constants, respectively; and l is the the distance from

M i to O b .
As proposed by [START_REF] Bouabdallah | Design and Control of Quadrotors with Applicatioin to Autonomous Flying[END_REF], the dynamics of a quadrotor is obtained by using the rigid body equations, which are

m b V + b Ω × (m b V) = b F + b G , I b Ω + b Ω × (I b Ω) = b τ . . ( 12
)
For measurement and experimental purposes it easier to refer the translational dynamics to F i , while the attitude dynamics is referred to the body frame. Thus, equation ( 12) can be rewritten as

m i R b b V = i R b ( b F + b G) , I b Ω + b Ω × (I b Ω) = b τ , , (13)
where

I =   I xx 0 0 0 I yy 0 0 0 I zz   , (14)
is the moment of inertia matrix, m the quadrotor's mass and i G the gravity vector. The Coriolis effect no longer appears in this equation since it does not affect the system's behaviour when the inertial frame is used as reference. After developing some algebra, the resulting equations are

b Ω =   ṗ q ṙ  =        I yy -I zz I xx q r + τ φ I xx I zz -I xx I yy p r + τ θ I yy I xx -I yy I Izz p q + τ ψ I zz        , (15)
for the position dynamics, and

i ξ =   ẍ ÿ z  =      [sin(φ) sin(ψ) + cos(φ) sin(θ)cos(ψ)] u m [-sin(φ)cos(ψ) + cos(φ)sin(θ)sin(ψ)] u m -g + cos(φ)cos(θ) u m      , (16)
for the attitude dynamics.

III. CONTROL SYSTEM DESIGN

A. System discretisation

For the design of the linear Model Predictive Controller, the system should be represented in a discrete state-space form. Therefore, an state augmentation using equation ( 16), followed by a linearisation assuming an equilibrium point in φ(0) θ(0) ψ(0) T = 0 0 0 T , is performed. In addition, the sampling time for the discrete system was set to be T s = 0.02 s (this is the overall sampling time of the experimental system, which will be explained in the respective section). This evaluations result in

        ẋd ẏd żd ẍd ÿd zd         =           1 0 0 0.02 0 0 0 1 0 0 0.02 0 0 0 1 0 0 0.02 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1                   x y z ẋ ẏ ż        +         0 0 0.001962 0 0 -0.001962 0 0 0.0003717 0 0 0 0 0 0.1962 0 0 -0.1962 0 0 0.03717 0 0 0             u φ θ ψ     . (17)
The physical parameters of the quadrotor m = 0.5380 kg and g = 9.81 m • s -2 were substituted for simplification purposes. The inputs of the control system will be the thrust and the euler angles, they will be sent to the inner attitude control loop, which is assumed to have a quick response time (< 0.6 s). The inner loop has been tuned in the PX4 board (multirotor attitude of Pixhawk is used as default).

B. Linear predictive model

The objective of this control formulation is to make the system's output y(t), track a reference signal r(t) with the smallest possible error. For this purpose, it is necessary to reformulate the state space representation, starting from

x(k + 1) = A x(k) + B u(k) , y(k) = C x(k) , (18)
and assuming that ∆u(k) = u(k) -u(k -1), the control inputs vector u(k) in equation ( 18) can be rewritten as

x(k + 1) = A x(k) + B u(k -1) + B∆u(t) , u(k) = u(k -1) + ∆u(k) . ( 19
)
If ∆u is considered as the vector of control inputs, and u(k -1) is added to the states vector [START_REF] Maciejowski | Predictive Control with Constraints[END_REF]. The system can be reordered as

           x(k + 1) u(t) = A B 0 I x(k) u(k -1) + B I ∆u(k) , y(k) = C 0 x(k) u(k -1)
.

(20)

For simplification purposes the new state vector will be

x(k) = x(k) u(k -1)
T . Consequently, the system can be represented in a more compact way as

x(k + 1) = A x(k) + B ∆u(k) , y(k) = C x(k) . (21)
The current time k is always available from the measurements, x(k) provides the current value of the plant states. In the same manner, the future control trajectory is denoted by

∆u = ∆u(k) T ∆u(k + 1) T . . . ∆u(k + n c -1) T T , (22)
where n c is called control horizon and it represents the number of parameters used to capture the future control trajectory. In addition, given the model of the system (21), x(k) can be computed for a number n p of samples; n p is known as the prediction horizon, and it is used to compute the effects of the current control input in the future. It is necessary to remark that the control horizon should always be lower or equal than the prediction horizon [START_REF] Kouvaritakis | Model Predicitive Control: Classical, Robust and Stochastic[END_REF].

From the state-space representation of equation ( 21), the future value of the state variables

x(k) = x(k + 1|k) T x(k + 2|k) T . . . x(k + m|k) T . . . x(k + n p |k) T T ,
can be computed using the future values of the control inputs as

x(k + 1|k) = A x(k) + B∆u(k), x(k + 2|k) = A x(k + 1|k) + B∆u(k + 1), = A 2 x(k) + A B∆u(k) + B∆u(k + 1), . . . x(k + n p |k) = A np x(k) + A np-1 B∆u(k) + A np-2 B∆u(k + 1) + . . . + A np-nc B∆u(k + n c -1). (23)
In the same way, the future value of the outputs

y(k) = y(k + 1|k) T y(k + 2|k) . . . y(k + n p |k) T T , (24)
can be estimated by using the predicted state previously defined, leading to

y(k + 1|k) = CA x(k) + CB∆u(k), y(k + 2|k) = CA 2 x(k) + CAB∆u(k) + CB∆u(k + 1), y(k + 3|k) = CA 3 x(k) + CA 2 B∆u(k) +CAB∆u(k + 1) + CB∆u(k + 2), . . . y(k + n p |k) = CA np x(k) + CA np-1 B∆u(k) + CA np-2 B∆u(k + 1) + . . . +CA np-nc B∆u(k + n c -1). (25)
From vectors (23) and ( 25), a so-called estimation model of the form (18) can be obtained as follows:

y(k) = F x(k) + Φ∆u , (26)
where

F = CA CA 2 . . . CA np T ,
and

Φ =        CB 0 0 . . . 0 CAB CB 0 . . . 0 CA 2 B CAB CB . . . 0 . . . . . . . . . . . . . . . CA np-1 B CA np-2 B CA np-3 B . . . CA np-nc B       
.

C. Optimisation: Tracking Formulation

The objective of the position controller proposed in this work is to bring the system outputs y(k) as close as possible to the reference r(k). Therefore, a performance index function

J(k) = np-1 k=0 W y (y(k) -r(k)) 2 + W ∆u ∆u(k) 2 , ( 27
)
is defined, where W y = diag(w y0 , w y1 , . . . , w yn ) , W ∆u = diag(w ∆u0 , w ∆u1 , . . . , w ∆un ), are called weighting matrices. Given this formulation, the new problem is to find the best control input ∆u such that the tracking error is minimised. Therefore, equation ( 27) can be rewritten in a matrix form as

J = (y(k) -r(k)) T W y 2 (y(k) -r(k)) + ∆u T W ∆u 2 ∆u . (28)
Then, the substitution of (26) in the previous equation leads to

J = (F x(k) -r(k)) T W y 2 (F x(k) -r(k)) +2∆u T Φ T W y 2 (F x(k) -r(k)) +∆u T (Φ T W y 2 Φ + W ∆u 2 )∆u , (29)
which minimal value can be obtained from its derivative with respect to the control input ∆u as

∂J ∂∆u = 2Φ T W y 2 (F x((k)-r(k))+2(Φ T W y 2 Φ+W ∆u 2 )∆u = 0, (30)
therefore, the resulting controller will have the form

∆u = -(Φ T W y 2 Φ + W ∆u 2 ) -1 Φ T W y 2 (F x(k) -r(k)). ( 31
)
It can be seen that the final controller is directly linked to the reference signal and the system states, so it should be computed in real time at any new sample.

IV. TARGET DETECTION AND LANDING ALGORITHM

For this investigation, not only a MPC controller was developed and proposed, but also a guidance law. The ground robot to be used as tracking target is presented in Figure 2. The final guidance law has four states: 1) Take-off: The UAV turns-on the motors, consequently, it takes-off from either the ground level z(0) = 0 or the mobile robot z(0) = h rob . The controller will stabilise the UAV at a given altitude z r (k) = z track , while x r (k) and y r (k) will remain constant. The system will immediately switch to the next state.

x r , y r , z r , v xr , v yr , v zr = x(0), y(0), z track , 0, 0, v z track .

2) Patrol mode: The UAV starts moving within the flying arena until the camera system detects the robot in a near position. The selected tracking pattern is a spiral on XY i plane with variable amplitude and frequency. The condition to switch to the next state is the detection of the robot.

x r , y r , z r , v xr , v yr , v zr = x pat , y pat , z track , v xpat , v ypat , v z track .

3) Tracking mode: As soon as the target is detected, the quadrotor starts tracking it at the given z r (k) = z track , while the position (x rob and y rob ) of the robot measured by the camera is sent to the control loop as reference. The robot velocity is estimated by means of a discrete derivative with a second order low-pass filter since the camera system is only able to provide position measurements. If the robot is lost at any time during this phase, the state machine will switch back to the previous state.

x r , y r , z r , v xr , v yr , v zr = x rob , y rob , z track , v x rob , v y rob , v z track .

4) Landing mode:

The UAV descends if the velocity of the robot is smaller that x robmax and y robmax , after reaching the platform the motors are switched-off. If the robot is lost during the descending phase, the state machine will return to the second state.

x r , y r , z r , v xr , v yr , v zr = x rob , y rob , -z(k) + 1, v xr , v yr , żref .

All the phases previously presented are shown in a more illustrative way in Figure IV, where also the change conditions are defined. At ISAE Supaero, a complete experimental framework has been developed along the years by the researchers and students. This system allows one to test high performance control laws within a Simulink interface. The experimental environment is presented in Figure 4, the main source of measurements is an Optitrack motion capture system. For this investigation, an ARdrone2.0 frame is used, on which the embedded system has been replaced by a PX4 board for the attitude control loop, and a Gumstix board for the guidance and high-level control. The Simulink Project called SimoDrones allows one to choose the quadrotor configuration for the experimentation, either a real time experiment (with Simulink running in real time withing a Windows kernel) or simulation can be selected. SimoDrones includes many non-linear models (brushless motor dynamics, communication delays, measurement noise, quadrotor dynamics and aerodynamic drag) in order to validate advanced control laws first in simulations, and then in real experimentation. The results that are presented later on where directly extracted from SimoDrones after performing reals tests.

The final cascade architecture control and its corresponding Simulink model are presented in Figures 5 and6, respectively.

For the tuning of the controller the following parameters an constraints were selected:

• n p = 10,

VI. EXPERIMENTAL RESULTS

A. Comparison: MPC vs LQR, PID and Backstepping

For the performance analysis, the controller proposed in this paper was compared to some other control techniques available at ISAE Supaero drones library, they are: PID, LQR and Backstepping [START_REF] Lecointe | Backstepping control law application to path tracking with an indoor quadrotor[END_REF]. The results are presented in Figure 7, during this experimental phase, the reference signal is a set of way-points where the guidance function switches to the next way-point in function of the steady state error value and the the time response of each controller. Concerning the results, the Backstepping and the MPC have the better response, it can be observed that PID and LQR need a longer time to minimise the steady state error; consequently, to finish the full set of way-points, a longer time is required. The Backstepping technique offers the shortest delay between the reference signal and the output, and the MPC has a slightly slower response, which for the experimental purposes is almost negligible (≈ 0.1 s). In terms of power consumption and aggressiveness of the control input, the MPC has shown to be a better alternative. When the tracking error is analysed, both controllers have a similar response; however, the overshot provided by the MPC is a bit smaller, and therefore a better alternative for the target tracking problem.

B. Landing on the moving target

After proving that the MPC minimised faster the steady state error with the minimal amount oscillations, and its overshot was the smallest among the different techniques, some tests using the new guidance algorithm proposed in section IV were carried out. Figure 8 The mission profile has the following sequence: at 5 s the quadrotor takes-off, when the desired z track (at around 10 s) is reached, the new reference will be the patrol spiral. As soon as the camera detects the robot within the vision range, the reference is switched to the measurement of the robot position. Subsequently, the UAV starts descending at around 63 s. At some point during this phase, the velocity of the robot was too high, therefore the quadrotor waited until it decreased to keep descending. Finally, the motors are switched-off as soon as the UAV touches the platform. The whole mission takes around 65 s, however this time can be improved given that some protection flags are used. The video of the mission can be retrieved from [START_REF] Macés-Hernández | Autonomous landing of an UAV on a moving platform[END_REF].

VII. CONCLUSIONS

In this paper, it has been shown the validity of using MPC for accurate target tracking. The MPC was successfully implemented in the experimental framework and many simulations and experimental tests were carried out in order to validate this proposition. The predictive model model developed for this investigation has correctly described the position dynamics of the UAV. This model has also proven to be accurate for the estimation of the future control inputs and states in the linear predictive model.

By comparing the developed controller against the other available control techniques, the one propose in this paper has demonstrated to be the best solution for a mission where it is important to have an accurate tracking. However, the controller had a short delay (≈ 0.1 s) due to the computations of the control input at every sampling time, which for experimental purposes did not compromise the response of the system since this value is very short.

The state machine based guidance law proposed in this paper has also proven to have an excellent behaviour in the experimental framework. All the transitions between states were optimised for the safety and integrity of the UAV. After performing many tests, the results show that the guidance function can perfectly model any UAV mission, and some more complex functions can be included in the future.

Many test were carried out, which demonstrated the validity of both solutions working together in the same experimental system. The resulting functions were added to the ISAE Supaero drones library.

Further studies would include the substitution of the internal attitude control loop by a full position and attitude MPC controller, using the same experimental framework.

Fig. 1 .

 1 Fig. 1. Inertial and body frames

Fig. 2 .

 2 Fig. 2. Mobile ground robot

Fig. 3 .

 3 Fig. 3. State machine for the guidance law

Fig. 4 .

 4 Fig. 4. Experimental environment

Fig. 5 .

 5 Fig. 5. Final control architecture

Fig. 7 .

 7 Fig. 7. Controllers comparison

  is composed of four plots; the first three are the position measurements, and the last one shows the transitions of the state machine. Regarding the position in X i -and Y i -axes, there are three signals in each plot, in purple the real position of the robot, in blue the patrol reference signal, and finally the real position of the quadrotor in red.

Fig. 8 .

 8 Fig. 8. Landing sequence: states and references


Jonathan Detchart

Jérôme Lacan

Polynomial Ring Transforms for Efficient XOR-based Erasure Coding

The complexity of software implementations of MDS erasure codes mainly depends on the efficiency of the finite field operations implementation. In this paper, we propose a method to reduce the complexity of the finite field multiplication by using simple transforms between a field and a ring to perform the multiplication in a ring. We show that moving to a ring reduces the complexity of the operations. Then, we show that this construction allows the use of simple scheduling to reduce the number of operations.

I. INTRODUCTION

Most of practical Maximum Distance Separable (MDS) packet erasure codes are implemented in software. In the various applications like packet erasure channels [START_REF] Bloemer | An XOR-Based Erasure-Resilient Coding Scheme[END_REF] or distributed storage systems [START_REF] Plank | Xor's, lower bounds and mds codes for storage[END_REF], the coding/decoding process performs operations over finite fields. The efficiency of the implementation of these finite field operations is thus critical for these applications.

To speedup this operation, [START_REF] Bloemer | An XOR-Based Erasure-Resilient Coding Scheme[END_REF] described an implementation of finite field multiplications which only uses simple xor operations, contrarily to classic software multiplications which are based on lookup tables (LUT). The complexity of multiplying by an element, i.e. the number of xor operations, depends on the size of the finite field and also on the element itself. This kind of complexity is studied for Maximum-Distance Separable (MDS) codes in [START_REF] Blaum | On lowest density MDS codes[END_REF]. Other work has been done to reduce redundant xor operations by applying scheduling [START_REF] Luo | Efficient encoding schedules for xor-based erasure codes[END_REF].

Independently, in the context of large finite field for cryptographic applications, [START_REF] Itoh | Structure of parallel multipliers for a class of fields gf (2 m )[END_REF] proposed a xor-based method to perform fast hardware implementations of multiplications by transforming each element of a field into an element of a larger ring. In this polynomial ring, where the operations on polynomials are done modulo x n + 1, the multiplication by a monomial is much simpler as the modulo is just a cyclic shift. The authors identified two classes of fields based on irreducible polynomials with binary coefficients allowing to transform each field element into a ring element by adding additional "ghost bits".

In this paper, we extend their approach to define fast software implementations of xor-based erasure codes. We propose an original method called PYRIT (PolYnomial RIng Transform) to perform operations between elements of a finite field into a bigger ring by using simple transforms between these two structures. Working in such a ring is much easier than working in a finite field. Firstly, it reduces the coding complexity by design. And secondly, it allows the use of simple scheduling to reduce the number of operations thanks to the properties of the ring structure.

The next section presents the algebraic framework allowing to define the various transforms between the finite field and some subsets of the ring. Then we discuss about the choice of these transforms and their properties. We also detail the complexity analysis before introducing some scheduling results.

II. ALGEBRAIC CONTEXT

The algebraic context of this paper is finite fields and ring theory. More detailed presentation of this context including the proofs of the following propositions can be found in [START_REF] Macwilliams | The Theory of Error-Correcting Codes[END_REF] or [START_REF] Poli | Error correcting codes: theory and applications[END_REF]. Definition 1. Let F q w be the finite field with q w elements. Definition 2. Let R q,n = F q [x]/(x n -1) denote the quotient ring of polynomials of the polynomial ring F 2 [x] quotiented by the ideal generated by the polynomial x n -1.

Definition 3. Let p u1 1 (x)p u2 2 (x) . . . p ur r (x) = x n -1 be the decomposition of x n -1 into irreducible polynomials over F q .
When n and q are relatively prime, it can be shown that u 1 = u 2 = . . . = u n = 1 (see [START_REF] Poli | Error correcting codes: theory and applications[END_REF]). In other words, if q = 2, and n is odd, we simply have p 1 (x)p 2 (x) . . . p r (x) = x n -1.

In the rest of this document, we assume that n and q are relatively prime.

Proposition 1. The ring R q,n is equal to the direct sum of its r minimal ideals of A i = ((x n -1)/p i (x)) for i = 1, . . . , r.

Moreover, each minimal ideal contains a unique primitive idempotent θ i (x). A construction of this idempotent is given in [START_REF] Macwilliams | The Theory of Error-Correcting Codes[END_REF], Chap. 8, Theorem 6. Since F q [x]/(p i (x)) is isomorphic to the finite field B i = F q w i , where p i (x) is of degree w i , we have:

Proposition 2. R q,n is isomorphic to the following Cartesian product: R q,n B 1 ⊗ B 2 ⊗ . . . B r
For each i = 1, . . . , r, A i is isomorphic to B i . The isomorphism is:

φ i : B i → A i b(x) → b(x)θ i (x) (1)
and the inverse isomorphism is:

φ -1 i : A i → B i a(x) → a(α i ) (2)
where α i is a root of p i (x).

Let us now assume that q = 2. Let us introduce a special class of polynomials: Definition 4. The All One Polynomial (AOP) of degree w is defined as

p(x) = x w + x w-1 + x w-2 + . . . + x + 1
The AOP of degree w is irreducible over F 2 if and only if w + 1 is a prime and w generates F * w+1 , where F * w+1 is the multiplicative group in F w+1 [START_REF] Wah | Realization and application of the massey-omura lock[END_REF]. The values w + 1, such that the AOP of degree w is irreducible is the sequence A001122 in [START_REF] Sloane | A001122 sequence[END_REF]. The first values of this sequence are: 3, 5, 11, 13, 19, 29, . . .. In this paper, we only consider irreducible AOP. According to Proposition 2, R 2,w+1 is equal to the direct sum of its principal ideals A 1 = ((x w+1 + 1)/p(x)) = (x + 1) and A 2 = ((x w+1 + 1)/x + 1) = (p(x)) and R 2,w+1 is isomorphic to the direct product of

B 1 = F 2 [x]/(p(x)) = F 2 w and B 1 = F 2 [x]/(x + 1) = F 2 .
It can be shown that the primitive idempotent of A 1 is θ 1 = p(x) + 1. This idempotents is used to build the isomorphism φ 1 between A 1 and B 1 .

III. TRANSFORMS BETWEEN THE FIELD AND THE RING

This section presents different transforms between the field B 1 = F 2 w = F 2 [x]/(p(x)) and the ring R 2,w+1 = F 2 [x]/(x w+1 + 1).

A. Isomorphism transform

The first transform is simply the application of the basic isomorphism between B 1 and the ideal

A 1 of R 2,w+1 (see Prop. 2).
By definition of the isomorphism, we have:

φ -1 1 (φ 1 (u(x)).φ 1 (v(x))) = u(x).v(x)
So, φ 1 can be used to send the elements of the field in the ring, then, to perform the multiplication, and then, to come back in the field. We show in the following Proposition that the isomorphism admits a simplified version.

Let W (b(x)), the weight of b(x), defined as the number of monomials in the polynomial representation of b(x).

Proposition 3. φ 1 (b B (x)) = b A (x) = b B (x) if W (b B (x)) is even b B (x) + p(x) else φ -1 i (b A (x)) = b B (x) = b A (x) if b w = 0 b A (x) + p(x) else where b w is the coefficient of the monomial of degree w of b A (x). Proof: For the first point, we have φ 1 (b(x)) = b(x)θ 1 (x) = b(x)(p(x) + 1) = b(x)p(x) + b(x). We can observe that b(x)p(x) = 0 when W (b(x)) is even and b(x)p(x) = p(x) when W (b(x)) is odd. The first point is thus obvious.
For the second point, it can be observed that, from the first point of this proposition, if an element of A 1 has a coefficient b w = 0, then it was necessarily obtained from the second rule, i.e. by adding p(x). Then, its image into B 1 can be obtained by subtracting (adding in binary) p(x). If b w = 0, then nothing has to be done to obtain b B (x).

B. Embedding transform

Let us denote by φ E the embedding function which simply consists in considering the element of the field as an element of the ring without any transformation. This function was initially proposed in [START_REF] Itoh | Structure of parallel multipliers for a class of fields gf (2 m )[END_REF].

Note that the images of the elements of B 1 doesn't necessarily belong to A 1 . However, let us define the function φ-1

1 from R 2,w+1 to A 1 by φ-1 1 (b A (x)) = b A (α)
, where α is a root of p(x). This function can be seen as an extension of the function φ -1 1 to the whole ring. Proposition 4. [START_REF] Itoh | Structure of parallel multipliers for a class of fields gf (2 m )[END_REF] For any u(x) and v(x) in B 1 , we have:

φ-1 1 (φ E (u(x)).φ E (v(x))) = u(x).v(x)
Proof. The embedding function corresponds to a multiplication by 1 in the ring. In fact, 1 is equal to the sum of the idempotents

θ i (x) of the ideals A i , for i = 1, . . . , l [6, chapter 8, thm. 7]. Thus, φ E (u(x)) = u(x). l i=0 θ i (x). Then, φ E (u(x)).φ E (v(x)) is equal to u(x).v(x).( l i=0 θ i (x)) 2 .
Thanks to the properties of idempotents,

θ i (x).θ j (x) is equal to θ i (x) if i = j and 0 else. Thus, φ E (u(x)).φ E (v(x)) is equal to u(x).v(x).( l i=0 θ i (x)). The function φ-1 i is the computation of the remainder modulo p i (x). The irre- ducible polynomial p i (x) corresponds to the ideal A i . Thus θ i (x) mod p(x) is equal to 1 if i = 1 and 0 else.
This proposition proves that the Embedding function can be used to perform a multiplication in the ring instead of doing it in the field. The isomorphism also has this property, but the complexities of the transforms between the field and the ring are more complex.

C. Sparse transform

Let us define the transform φ S from B 1 to R 2,w+1 :

φ S (b B (x)) = b A (x) = φ 1 (b B (x)) + δ.p(x) where δ = 1 if W (φ 1 (b B (x)) + p(x)) < W (φ 1 (b B (x))
) and 0 else. Proposition 5. For any u(x) and v(x) in B 1 , we have:

φ-1 1 (φ S (u(x)).φ S (v(x))) = u(x).v(x)
Proof. As observed in the proof of Prop. 4, φ-1 i is just the computation of the remainder modulo p(x). Moreover, according to the definition of φ S , φ S (u(x)).φ S (v(x)) is equal to u(x).v(x) plus a multiple of p(x) (possibly equal to 0). Thus, the remainder of φ S (u(x)).φ S (v(x)) modulo p(x) is equal to u(x).v(x).

This proposition shows that φ S can be used to perform the multiplication in the ring. The main interest of this transform is that the weight of the image of φ S is small, which reduce the complexity of the multiplication in the ring.

D. Parity transform

Proposition 6. The ideal A 1 is composed of the set of elements of R 2,w+1 with even weight.

Proof: We can observe from Proposition 3 that all the image of φ 1 have even weight. Since the number of evenweight element of R 2,w+1 is equal to the number of elements of A 1 , A 1 is composed of the set of elements of R 2,w+1 with even weight.

Let us consider the function φ P , from B 1 to R 2,w+1 , which adds a single parity bit to the vector corresponding to the finite field element. The obtained element has an even weight (by construction), and thus, according to the previous Proposition, it belongs to A 1 .

Since the images by φ P of two distinct elements are distinct, φ P is a bijection between B 1 and A 1 . The inverse function, φ -1 P , consists just in removing the last coefficient of the ring element.

It should be noted that φ P is not an isomorphism, but just a bijection between B 1 and A 1 . However it will be shown in next Section that this function can be used in the context of erasure codes.

IV. APPLICATION OF TRANSFORMS

In typical xor-based erasure coding systems [START_REF] Bloemer | An XOR-Based Erasure-Resilient Coding Scheme[END_REF], the encoding process consists in multiplying an information vector by the generator matrix. Since in software, xor are performed using machine words of l bits, l interleaved codewords are encoded in parallel.

We consider a system with k input data blocks and m output parity blocks.

The total number of xor of the encoding is thus defined by the generator matrix which must be as sparse as possible. First, we use a k × (k + m)-systematic generator matrix built from a k × k-identity matrix concatenated to a k × m Generalized Cauchy (GC) matrix [START_REF] Roth | On generator matrices of mds codes (corresp.)[END_REF]. A GC matrix generates a systematic MDS code and it contains only 1 on its first row and on its first column. Then, to improve the sparsity of the generator matrix in the ring, we use the Sparse transform φ S . This has to be done only once since the ring matrix is the same for all the codewords.

For the information vectors, it is not efficient to use φ S since the xors of machine words do not take into account the sparsity of the xor-ed vectors. We thus use Embedding or Parity transforms, which are less complex than φ 1 .

When Embedding is used for information vectors and Sparse is used for the generator matrix, the obtained result in the ring can be sent into the field by using φ -1 1 (proof similar to the proofs of Propositions 4 or 5).

When Parity is used for the information vector, the image of the vector in the ring only contains elements of the ideal A 1 . Since these elements are multiplied by the generator matrix (in the ring), the obtained result only contains elements of the ideal A 1 . These elements have even weight, so it is not necessary to keep the parity bit before sending them on the "erasure channel". Since Parity transform is not an isomorphism, these data can not be decoded by another method. Indeed, to decode, it is necessary to apply φ P (add the parity bit), then to decode by multiplying by the inverse matrix, and then to to apply φ -1 P (remove the parity bit on the correct information vector).

V. COMPLEXITY ANALYSIS

In this section, we determine the total number of xor operations done in the coding and the decoding processes.

A. Coding complexity

The coding process is composed of three phases: the field to ring transform, the matrix vector multiplication and the ring to field transform. We assume that the information vector is a vector of k elements of the field F 2 w .

For the first and the third phases, Table I The choice between the two methods thus depends on the values of the parameters: if k > m, Parity transform has lower complexity. Else, "Embedding" complexity is better.

For the matrix vector operation, let us first consider the multiplication of two ring elements. As explained in the previous section, the first element (which corresponds to an information symbol) is managed by the software implementation by machine words. So the complexity of the multiplication only depends on the weight of the second element, denoted by w 2 ∈ {0, 1, . . . , w + 1}. The complexity of this multiplication is thus (w + 1).w 2 . Now, we can consider the specificities of the various transforms. In the Parity transform, the last bit of the parity blocks is not used ( i.e. it is not transmitted on the erasure channel). So it is not necessary to compute it. It follows that the complexity of the multiplication is only w.w 2 .

Similarly, for the Embedding transform, the last bit of the input vector is always equal to 0. So, we also have a complexity equal to w.w 2 .

To have an average number of operations done in the multiplication of the generator matrix by the input data blocks, we have to evaluate the average weight of the entries of the generator matrix in the ring.

The generator matrix is a k × m-GC matrix with the first column and the first row are filled by 1. The other elements can be considered as random non zero elements. They are generated by φ S which chooses the lowest ring element among the two ones corresponding to the field element. Let us denote their average weight by w φ S . For this case, the average number of xors is thus:

(k + m -1).

w + (k -1).(m -1).w.w φ S

This leads to the following general expression of the coding complexity:

(min(k, m) + k + m -1).w + (k -1).(m -1).w.w φ S
To estimate the complexity on a practical example, we fix the value of w to 4. Classic combinatorial evaluation (not presented here) gives the average weight for nonzeros images of φ S :

w φ S = w + 1 2 w+1 -2 . 2 w - w w/2
So, w φ S = 1.66. We plot in Figure 1 the evolution of the factor over optimal (used e.g. in [START_REF] Plank | Xor's, lower bounds and mds codes for storage[END_REF], table III) which is the density of the matrix normalized by the minimal density, k.m.w. We vary the value of k for three values of m: 3, 5 and 7. For each pair (k, m), we generate 10000 random GC matrices and keep the best we found. We can observe that the values are very low. For example, [START_REF] Plank | Xor's, lower bounds and mds codes for storage[END_REF] gives the lowest density of Cauchy matrices for the field F 2 4 and we can observe that our values are always lower than these ones.

To reduce the complexity in specific cases, we can observe that the ring contains w elements whose the corresponding matrix is optimal (one diagonal). By using these elements, we can search by brute force MDS matrices built only with these optimal elements. For example, let us consider the elements of the field F 2 4 sent into the ring R 2,5 . The Vandermonde matrix defined by: V = x i.j i=0,...,4;j=0,...,4

where x is a monomial in R 2,5 has the minimal number of 1. It can be verified that this matrix can be used to build a systematic a MDS code. For this matrix, the total number of xors done in the generation of the parity packets (including the field-to-ring and ring-to-field transforms) is

(min(k, m)+k+m-1).w+(k-1).(m-1).w = k.w+k.m.w

Its factor over optimal is equal to 1.2 which is lower than the values given in Figure 1 and which is close to the lowest bound given in [START_REF] Blaum | On lowest density MDS codes[END_REF].

B. Decoding complexity

As the decoding is a matrix inversion and a matrix vector multiplication, we can use the same approach to perform the multiplication. We first invert the sub-matrix in the field, then we transform each entry of this matrix into ring elements. Then, we perform the ring multiplication.

The complexity of the decoding thus depends on the complexity of the matrix inversion and on the complexity of the matrix vector multiplication.

The complexity of the matrix vector multiplication was studied in the previous paragraph.

The complexity of the a r × r-matrix inversion is generally in O(r 3 ) operations in the field. But if the matrix has a Cauchy structure, this complexity can be reduced to O(r 2 ) [START_REF] Bloemer | An XOR-Based Erasure-Resilient Coding Scheme[END_REF].

Note that, contrarily to the matrix vector multiplication, the matrix inversion complexity does not depend on the size of the source and parity blocks. And thus, it becomes negligible when the size of the blocks increase.

VI. SCHEDULING An interesting optimization on MDS erasure codes under

xor-based representation is the scheduling of xor operations.

Such techniques are proposed in [START_REF] Hafner | Matrix methods for lost data reconstruction in erasure codes[END_REF], [START_REF] Huang | On optimizing xor-based codes for faulttolerant storage applications[END_REF], [START_REF] Plank | Xor's, lower bounds and mds codes for storage[END_REF] and [START_REF] Plank | Heuristics for optimizing matrix-based erasure codes for fault-tolerant storage systems[END_REF]. The general principle consists in "factorizing" some xor operations which are done several times to generate the parity blocks.

We show in the two next paragraphs that these techniques can be used very efficiently on the ring elements.

However, it can be observed that the matrices defined over rings have two main advantages.

A. Complexity reduction

Over finite fields, the scheduling consists in searching common patterns on the binary representation of the generator matrices. The w × w-matrices representing the multiplication by the field elements does not have particular structure and thus, they must be entirely considered in the scheduling algorithm. This is not the case for the (w + 1) × (w + 1)-matrices corresponding to a ring element because, thanks to the form of the polynomial x w+1 + 1, they are composed of diagonals either full of 0 or 1. This means that they can be represented in the scheduling algorithm just by their first column or, equivalently, by the ring polynomial.

This allows to drastically reduce the algorithm complexity and thus to handle bigger matrices. From a polynomial point of view, the search of scheduling just consist in finding some common patterns in the equations generating the parity blocks.

Example . Let us assume that n = 5 and that three data polynomials a 0 (x), a 1 (x) and a 2 (x) are combined to generate the three parities p 0 (x) = (1+x 4 )a 0 (x)+x 2 a 1 (x)+x 3 a 2 (x), p 1 (x) = a 0 (x) + x 3 a 1 (x) + (1 + x 3 )a 2 (x) and p 2 (x) = a 0 (x) + a 1 (x) + x 3 a 2 (x).

In this case, the scheduling just consists in computing p (x) = a 0 (x) + x 3 a 2 (x) and then p 0 (x) = p (x) + x 4 a 0 (x) + x 2 a 1 (x), p 1 (x) = p (x) + x 3 a 1 (x) + a 2 (x) and p 2 = p (x) + a 2 (x).

To estimate the complexity, we can consider the number of sums of polynomials. Without scheduling, we need 11 sums (4 for p 0 (x), 4 for p 1 (x), and 3 for p 2 (x)) instead of with scheduling, we only need 10 sums (2 for p (x), 3 for p 0 (x), 3 for p 1 (x) and 2 for p 2 (x)).

B. Additional patterns

Ring-based matrices allow to find more common patterns than field-based matrices. The main idea is to observe that, in the ring, we can "factorize" not only common operations, but also operations which are multiple by a monomial (i. e. cyclic-shift) of operations done in some other equations. This is possible only because the multiplications are done modulo x w+1 + 1.

Example . Let us assume that n = 5 and that three data polynomials a 0 (x), a 1 (x) and a 2 (x) are combined to generate the parities p 0 (x) = a 0 (x) + x 2 a 1 (x) + (1 + x 2 )a 2 (x), p 1 (x) = x 2 a 0 (x) + x 3 a 1 (x) + (x + x 4 )a 2 (x) and p 2 (x) = x 2 a 0 (x) + a 1 (x) + (x 2 + x 3 )a 2 (x). We can observe that, with a "simple" scheduling, it is not possible to factorize some operations.

However, by rewriting the polynomials, we can reveal factorizations:

p 1 (x) = x 2 a 0 (x) + x(x 2 a 1 (x) + a 2 (x)) + x 4 a 2 (x) and p 2 (x) = x 2 a 0 (x) + x 3 (x 2 a 1 (x) + a 2 (x)) + x 2 a 2 (x). So, if p (x) = x 2 a 1 (x) + a 2 (x), we have p 0 (x) = p (x) + a 0 (x) + x 2 a 2 (x), p 1 (x) = xp (x) + x 2 a 0 (x) + x 2 a 2 (x) and p 2 (x) = x 3 p (x) + x 2 a 0 (x) + x 2 a 2 (x).
To estimate the complexity by the same method than in the previous example, we need 11 polynomial additions with scheduling compared to 12 additions necessary without scheduling.

C. Scheduling results

To evaluate the potential gain of the scheduling, we have implemented an exhaustive search of the best patterns on generator matrices.

This algorithm was applied on several codes for the field F 2 4 . Table II presents the results in term of "factor over optimal" which is defined as the total number of 1 in the matrix over the number of 1 for the optimal MDS matrix , i.e k.m.w.

When working in a ring, we include to the complexity the operations needed to apply the transforms. In this case, Embedding transform has a lower complexity. So we added m.w to the number of 1 in the matrix resulting from the scheduling algorithm.

For each case, we have generated 100 random Generalized Cauchy matrices.

The measured parameters are:

• average field matrix: average number of This table confirms that, even without scheduling, ring matrices have a lower density than field matrices, thanks to the Sparse transform. Applying scheduling to these matrices allows a significant gain of complexity. Indeed, it reduces the complexity by more than 20% on the best matrices. The final results are similar to the results obtained (without scheduling) on the optimal matrix in Section V-A. To the best of our knowledge, other scheduling approaches do not reach this level of sparsity for these parameters.

VII. CONCLUSION

In this paper, we have presented a new method to build MDS erasure codes with low complexity. By using transforms between a finite field and a polynomial ring, sparse generator matrices can be obtained. This allows to significantly reduce the complexity of the matrix vector multiplication. It also allows simple schedulers that drastically improve the complexity by reducing the number of operations. Similar results can be obtained with Equally-Spaced Polynomials (ESP) [START_REF] Itoh | Structure of parallel multipliers for a class of fields gf (2 m )[END_REF], but they are not presented here due to lack of space.

Fig. 1 .

 1 Fig. 1. Factor over optimal depending on m

  gives the complexities of Embedding and Parity transforms obtained from their definition in Section III.

		field to ring	ring to field
	Embedding	0	m.w
	Parity bits	k.w	0
		TABLE I
	NUMBER OF XOR FOR EMBEDDING AND PARITY TRANSFORMS

  1 in the GC matrices divided by k.m.w • best field matrix: lowest number of 1 among the GC matrices divided by k.m.w • average ring matrix: average number of 1 in the ring matrices (without scheduling) + ring-field correspondence divided by k.m.w • best ring matrix: best number of 1 among the ring matrices (without scheduling) + ring-field correspondence divided by k.m.w • average with scheduling: average number of xors with scheduling + ring-field correspondence divided by k.m.w

	k+m,k	12,8	16,10
	average field matrix	1.79	1.90
	best field matrix	1.73	1.85
	average ring matrix	1.59	1.63
	best ring matrix	1.5	1.58
	average with scheduling	1.32	1.26
	best with scheduling	1.19	1.20
	TABLE II
	FACTOR OVER OPTIMAL FOR w = 4

• best with scheduling: best number of xors with scheduling + ring-field correspondence divided by k.m.w