Willy Allègre
email: wallegre@kerpape.mutualite56.fr

Thomas Burger

Jean-Yves Antoine

Pascal Berruet

Jean-Paul Departe

A Non-Intrusive Context-Aware System for Ambient Assisted Living in Smart Home

Keywords: Smart Home, Ambient Assisted Living, Ontologies, Model-Driven Engineering, Home Automation Services, Context-Awareness

Automation of smart home for ambient assisted living is currently based on a widespread use of sensors. As efficient as it seems to be, this solution can sometimes be problematic when one focus on user acceptability intimately related to cost and intrusivity. In this paper, we propose a context-aware system based on the semantic analysis of each user request at runtime. Our goal is to infer user data usually sensored by using advanced semantic web tools to adapt home automation services to people with special needs. To take up this challenge, an ontology, automatically derived from a model-driven process, firstly defines usersystem interactions. Then, the use of rules allows an inference engine to deduce user location and intention leading to adapted service delivery.

Introduction

Care for dependents is becoming a major social and economic issue for the next few years. In 2050, 30% of people from the European countries will be at least 65 years old [1], thus increasing the number of dependent people with chronic diseases. The number of potential caregivers can not evolve accordingly, so that home automation for ageing and disabled people is considered as an alternative solution to fill the gap, allowing one to facilitate or automate devices activation and to provide services tailored to the user needs: we talk about Smart Home. Due to the miniaturization of electronic devices and their dissemination in the environment, smart home is no longer confined to the simple control of household appliances. This field of application is evolving more and more toward Pervasive Systems [START_REF] Weiser | Some computer science issues in ubiquitous computing[END_REF], where context-awareness is an important challenge. Indeed, knowing who is doing what and where is essential to adapt a "smart behaviour", especially in the area of Ambient Assisted Living. To adapt a house functioning to the user needs, this kind of systems must be able to model and process contextual information (e.g. location, preferences, time markers). Sensors is the most commonly used solution to collect data for home automation [START_REF] Yu | Human localization via multi-cameras and floor sensors in smart home[END_REF][START_REF] Nourizadeh | Medical and home automation sensor networks for senior citizens telehomecare[END_REF][START_REF] Long | Using graphs to improve activity prediction in smart environments based on motion sensor data[END_REF][START_REF] Chen | Sensor-Based Activity Recognition[END_REF], from which context information and knowledge can be inferred. If this solution is efficient, it can sometimes be problematic when one focus on user acceptability. The use of cameras can be seen as intrusive, whereas the massive deployment of sensors in the environment mainly implies a significant financial cost.

We strongly believe that user data in home automation is still underexploited as the trend in this field is clearly to add extra sensors instead of optimizing the home automation system in place. In the framework of a fully equipped home automation environment, as it is sometimes needed for some wheelchair users or elderly people in loss of autonomy, we also believe that the analysis of control logs (instead of sensor logs) can provide relevant information to adapt system to user needs. Control log is data directly coming from the user remote control whereas sensor log is coming from sensors connected to the home automation bus. This is why analysing the first type of data instead of continuously adding new sensors to gather user data can be seen as a non-intrusive solution (no added devices in the environment). We finally believe that control log analysis provides a new point of view to tackle key challenges in home automation such as realworld deployment, sensor optimisation, but also strategies for dealing with more common issues such as multiple occupancy, sensor failure or interleaved activity.

To take up this challenge, an ontology is defined to support the semantic analysis of each control log. Our goal is to provide the user adapted services without any training steps and using as few sensors as possible. Figure 1 illustrates the basic use of the ontology on a sequence diagram. When the user asks for a service, the control log is processed by the home automation server, but also translated in the context ontology as an assertion. An inference engine apply defined rules on this ontology to enrich the knowledge base by i) higher level information that a third party (occupational therapists, family members, etc.) can analyse to propose new services/rules and ii) services proposition on the user interface based on information previously inferred.

As ontologies and related tools are of great interest for context-awareness in ambient assisted living, mainly due to their interoperability and logic reasoning capabilities, there is a lot of works in this domain1 . However none of them focused on the analysis of control logs to provide context information about the user. Morevover, very few works provide a design process to ease the ontology deployment which is not a trivial task. In this paper, we describe our non-intrusive context-aware system introduced above as an alternative solution to assist people at home. To do so, we first propose to model user-system interactions and their semantics by means of an ontology (section 2). Then, the use of rules applied to ontology concepts and properties allows an inference engine to deduce user location and intention, leading to service delivery (section 3). In the context of a large scale deployment, a model-driven process is applied for populating the ontology and a software architecture is defined for runtime (section 4). A case Fig. 1. Basic use of the ontology for the non-intrusive context-aware system study (section 5) based on a real dataset shows the benefits of the entire work. Finally, related works (section 6) are presented before concluding.

An ontology for user-system interaction modeling

The ontology described below models the semantics underlying interactions between one user and his environment. As it is not possible to distinguish who activates a switch, we assume that the user stands alone. Obviously, our work should fit the need of a multi-user application, provided each log can be associated with an individual user id (e.g. services are requested from individual interfaces).

The ontology is made of generic concepts (describing context at a high level of abstraction) which can be applied on pervasive systems (Upper Class). Other concepts are responsible for the description of domain specific concepts for home automation (Specific Class). As shown in figure 2, the defined ontology is made of six fundamental concepts:

-User: someone who interacts with his environment using spread devices -Device: device activated by user i) remotely or ii) from the same room -Basic Service: directly provided by both local and remote devices -System: interface to collect interactions composed of devices and services -Location: part of the environment that can contain devices and users -Effect: effect of services activation on the environment (e.g. luminosity)

These ontology concepts are a useful support for our non-intrusive contextaware system. Indeed, from the way the User interacts with his System, we seek to infer his Location and his intention related to the Effect produced by Device activation when he asks for a Basic Service. Knowing where the user is and what is his intention, it is then possible to lead to service delivery. This is made possible by the formal semantics of description logics underlying ontologies, but also by rule-based reasoning.

3 From a non-intrusive monitoring to service delivery

User location

Location is a context information from which user activities can be inferred. In this work, we propose to infer such information, usually deduced from sensors data (e.g. motion detectors, switches contacts), from services request. In addition to ontology concepts described above, we add semantic relationships to link them. Figure 3 (left part) shows an excerpt of the ontology taxonomy (TBox with concepts and object properties) for user location. Local activation devices Local Device and remote activation devices Remote Device (concepts mutually exclusives) are defined to infer the property isLocatedIn related to the user. Rule 1 from table 1 allows an inference engine to deduce that a user ?u who is using a local activation device ?d which stands in room ?r, is himself in this room. Door concept is defined with exactly two links isBetween towards two concepts Room. Given this logic representation, rules are expressed to infer properties such as isEntering and isLeaving from a door activation (see rules 2 and 3 from table 1). A user ?u located in room ?r1 who is asking for the door ?d (between ?r1 and ?r2 ) opening, is actually entering ?r2 and leaving ?r1.

User intention

When the user asks for a service, he knows that this activation will have an effect on his environment. In other words, this latter asks for a home automation service for its effect. In this work, we seek to infer the desired effect from a Table 1. Excerpt of rules for location/intention (1,2,3/4) and service delivery [START_REF] Long | Using graphs to improve activity prediction in smart environments based on motion sensor data[END_REF][START_REF] Chen | Sensor-Based Activity Recognition[END_REF] semantic analysis of control logs. Figure 3 (right part) shows an excerpt of the ontology taxonomy for user intention, based on the Effect concept. Luminosity, Darkness, Heat, Quietness are sub-classes of effect which can be linked with Basic Service through the hasEffect property. The user intention can thus be deduced by an inference engine applying rule 4 from table 1: a user ?u who is asking for a basic service ?b, wants its associated effect ?e.

Service delivery

From context information on user location and intention previously inferred, the inference engine can infer services proposition. Figure 3 (bottom part) shows an excerpt of the ontology taxonomy for service delivery. The System concept, which acts as an interface for user-system interactions, is added to concepts and properties described above. Service delivery is based on suggestDe-viceActivation and SuggestServiceActivation properties, applied respectively on System and Device (see rules 5 and 6 from table 1). For example, a user who wants luminosity (e.g. wants(User1, Luminosity)), and who enters a room (e.g. isEntering(User1, Bedroom)) will be offered services allowing to increase luminosity in the incoming room.

The design and deployment of such an ontology has to be flexible to adapt to user needs and evolving environments. On the prospect of a large scale service deployment, a model-driven process developed in previous works [START_REF] Allègre | Model-Driven Flow for Assistive Home Automation System Design[END_REF] is applied for populating the ontology and a software architecture is defined for runtime.

4 Implementation: from design to deployment

A model-driven architecture for ontology design

In the context of ambient assisted living, we strongly believe that, in addition to the home automation expert, a non-expert -but someone who is able to properly consider the context of living with a disability (e.g. an occupational therapist, a family member) -must be included in the design process. We previously developed in [START_REF] Allègre | Model-Driven Flow for Assistive Home Automation System Design[END_REF] a model-driven flow for the design of assistive home automation systems. Platform independent modeling (PIM ) are separated from platform specific modeling (PSM )2 , so that non-expert designers can focus on specification rather than on implementation. Using a domain specific language, they can model both environment and interactions to define respectively the context of use and the user needs, before automatically generating the control code 3 .

In this paper, we propose to automatically generate ontology individuals with the same domain specific language. Figure 4 shows the context-awareness part which consists of several steps integrated in the existing design flow. In step 1', relevant model concepts defined by this latter at PIM level are transformed into ontological concepts (PSM level). They have some basic information such as the layout of housing, devices parameters, but also services defined to consider user needs. Step 2' is about generating individuals conformed to our ontology defined with OWL (Web Ontology Language), the most commonly used ontology language. The joint use of a model-driven architecture and ontology 4 makes Platform 1 Control [START_REF] Allègre | Model-Driven Flow for Assistive Home Automation System Design[END_REF] Domain Specific Language [START_REF] Allègre | Model-Driven Flow for Assistive Home Automation System Design[END_REF] Environment Interaction sense to benefit from their respective advantages, thus leading to the definition of a complete control and context-awareness system for smart home : MDA and languages derived from this paradigm provide tools, concepts, and languages to create and transform models, but offer few means to formally interpret the semantics of models at runtime. On the contrary, ontologies make knowledge machine-understandable, by formalizing and by making explicit the semantic relations between concepts.

Contextawareness

PSM

A software architecture for ontology deployment

The ontology deployment has to deal with scalability at runtime. From models defined with the domain specific language developed in previous works (see previous section), individuals are generated to populate the ontology. Associated with SWRL rules, it composes the knowledge base (T-Box, A-Box, and rules) representing user-system interaction. Each time the user asks for a service, the application have to process the control log, to update the ontology, to write logs in files and to send back services proposition on user interface. Figure 5 shows the proposed software architecture for the context-aware application. Written in Java using the OWL API [START_REF] Horridge | The owl api: A java api for owl ontologies[END_REF] to interact with the knowledge base and the HermiT [START_REF] Shearer | HermiT: A highly-efficient OWL reasoner[END_REF] inference engine, this application is based on six units:

1. OWL Parser takes as input control logs coming from the home automation server and translates it in OWL format 2. Supervision Unit loads the ontology and adds assertions related to logs (e.g. the log "User1 Switch ON.TV1 " is added as assertions "User1 isAsk-ingForService Switch ON " and "User1 isUsing TV1 ") 3. Inference Unit calls the inference engine for applying SWRL rules on ontology with newly added assertions and sends inferences to the filtering unit

A person tracking application

To validate our approach, a person tracking application is first presented. This case study is based on a real monitoring dataset recorded in the Domus smart home [START_REF] Kadouche | Support Vector Machines for Inhabitant Identification in Smart Houses[END_REF][START_REF] Chikhaoui | A New Algorithm Based on Sequential Pattern Mining for Person Identification in Ubiquitous Environments[END_REF] from the University of Sherbrooke (see figure 6). This dataset is the result of two series of experiments with six adults to evaluate the early morning routines (wake up, toileting, preparing and having breakfast) during 10 days for series 1 and 5 days for series 2. Each experiment lasted about 45 minutes. Table 2 shows the excerpt of activities performed in day 2 by user 1 [START_REF]Domus Monitoring Dataset[END_REF] chosen for this person tracking application. Since our aim is to analyse control logs and 5 HomeML is a format to store dataset openly accessible from a web repository [START_REF] Mcdonald | A web based tool for storing and visualising data generated within a smart home[END_REF].

the dataset is based on sensors, we have converted some of these sensors logs by user requests on devices when such devices can be controlled. For example, light sensor logs from Light5101, Light5105 and Light5106 can be replaced by control logs considering a light as a controllable device (represented by * symbol in table 2). Otherwise, the flow meter Fl01 is mounted on the cold water tap of the kitchen sink and sends an event when the tap is opened. As we consider this latter cannot be remotely controlled, we keep logs 8 and 9 as sensors logs (represented by • symbol in table 2). Fig. 6. A person tracking scenario in the Domus smart home, from [START_REF]Domus Monitoring Dataset[END_REF].

Step 1 The user (user1 ) switches on the bedroom light: "[Log 1]: Light51-05.SwitchON " (log 1 from table 2 and step 1 from figure 6). This log is introduced in the ontology through the following assertions Step 3 The system proposes him to switch off the bedroom light. This proposition is based on the inverse effect property (e.g. Luminosity isInverseOf Darkness). The user being located in the kitchen, the system proposes to switch off lights in other rooms. In this scenario, this service is accepted by the user: "[Log 4]: Light5105.SwitchOFF ".

Step 4 To go into the bathroom, he asks for the door opening: "[Log 6]: Door0216.Open". Considering this door linked with exactly two rooms with the isBetween property, rules 2 and 3 of table 1 can be applied. The user-door interaction allows the inference engine to extract the following information: [Inference 5]: User1 isEntering Bathroom and [Inference 6]: User1 isLeaving Kitchen.

Step 5 The system proposes him to switch on the bathroom light (log 7), from user intention and location previously inferred. These inferences are made possible trough rules 5 and 6 from table 1. Later, with the location information from sensor IR01, the system proposes him to switch off the bathroom light (as in step 3). This proposition is accepted by the user: "[Log 12]: Light5106.SwitchOFF ".

Results on the dataset

Results on the whole dataset, presented in table 3, can be interpreted in different ways according to the users needs. If the goal is to reduce the cost related to sensors installation and maintenance, it is possible to limit their use by replacing sensors data by a semantic analysis of control logs. If the goal is to increase dependability, it is possible to use our context-aware approach keeping sensors already installed to have information redundancy (i.e. a couple {sensor, actuator } when the user interacts with his environment). Limitation of sensors use The Domus Smart Home is equipped with infrared motion detectors, pressure sensors, switches contacts, flow meters, etc. As an illustration, the person tracking application based on the dataset excerpt presented in table 2 shows that our approach can limit the use of sensors to 37.5%. This means that only 3 sensors over 8 can potentially be used to ensure this application (see logs from table 2) : the flow meter Fl01 and the two infrared motion detectors IR01 and IR02 ) which cannot directly be associated with an actuator 7 . More generally, considering the whole dataset, actuators can replace electromagnetic door sensors and light switches. Thus, it should be possible to reduce the use of sensors in the Domus Smart Home to 71% for series 1 and 81% for series 2 (see table 3). The results vary from one user to another according to his/her mode of interaction with home automation devices. For example, the use of sensors for user 5 (in series 1) can be limited to 59%: this means that only 10 sensors over 17 can potentially be used (the 7 others can be replaced by actuators which logs can replace sensor logs). Sensors which cannot have associated actuators, like infrared motion detectors or pressure sensors for user location, can still be replaced by logic inference with a lower quality of services (QoS ) 8 .

Redundancy information

In a context of severe disabilities without economic constraint, quality of services may be preferred. If we decide to keep all sensors in the Domus Smart Home, our approach can improve the dependability by information redundancy. In this work, the information "The user interacts with his home automation system" may come from two different sources of information: sensors data (i.e. the device states), but also actuator data (i.e. user request sent). As an illustration, if we consider the dataset excerpt presented in table 2, one can see that 8/12 sensors logs can be associated with control logs (see logs with *). More generally, we observe on the whole data set that 16% of sensors logs (series 1) can be associated with either one control log or the result of logic reasoning to deduce user location (see table 3) 9 . Information redundancy depends directly on user-system interactions: for example, user 5 often interacts with his home automation devices what leads to a maximum redundancy of 48% on series 1. However, when the user walks through rooms without sending controls, no logic inference is made related to his location. The context-aware approach is described as proactive: data acquisition frequency depends on the home automation system use. This proactiveness also explains differences between results from series 1 and series 2. In the second series, user was asked to perform the early morning routine of series 1, and perform in addition an activity which does not require any controls (i.e. learning a tea recipe). The use of actuators decreases from 29% to 19% between the two series, leading to decrease the information redundancy (9% on average for series 2).

Related Works

The use of ontologies and related tools are of great interest in the domain of context-awareness for ambient assisted living, mainly due to their interoperability and logic reasoning capabilities. In [START_REF] Reinisch | Integration of heterogeneous building automation systems using ontologies[END_REF], an ontology is defined to model home automation environment on a large scale (building automation system), allowing to face to the heterogeneity of home automation platforms (e.g. KNX, LonWorks, BACnet technologies). Instead of using gateways between these different platforms (configurable with specific tools), the ontology is an abstract view of the heterogeneous network. DogOnt (Domotic OSGI Gateway) [START_REF] Bonino | Dogont-ontology modeling for intelligent domotic environments[END_REF][START_REF] Bonino | Automatic domotic device interoperation[END_REF] also support the hardware and telecommunication independent description of home automation environment. Using SWRL rules, logic reasoning is applied to automatically recognize devices classes for defining their technical parameters. Most of the work on activities of daily living recognition based on ontologies [START_REF] Chen | Sensor-Based Activity Recognition[END_REF][START_REF] Chen | A Knowledge-Driven Approach to Activity Recognition in Smart Homes[END_REF] relies on the use of triples {subject, predicate, object} to model context. CoBra [START_REF] Chen | The SOUPA ontology for pervasive computing[END_REF] and Gaia [START_REF] Román | A middleware infrastructure for active spaces[END_REF][START_REF] Bittencourt | Modeling jade agents from gaia methodology under the perspective of semantic web[END_REF] are two of the best known frameworks, based on a multi-agents architecture in charge of collecting and processing contextual information. For example, in [START_REF] Román | A middleware infrastructure for active spaces[END_REF], when a powerpoint application is used and the number of people in a room is more than three, system can infer that the room activity is a presentation: People(Room2401, >=, 3) and Application(PowerPoint, Running) ⇒ RoomActivity(2401, Presentation). In [START_REF] Liao | A RDF and owl-based temporal context reasoning model for smart home[END_REF], the TempCRM ontology extends the RDF triple commonly used with time. Depending on the temporal usage profile of home automation devices, potential emergencies can be detected. 9 Results are based on the original dataset without cleaning data from duplicates.

For example, Device(gas switch#1, is, on, 18:38:41) describes the opening of the gas which have to be closed after a certain amount of time.

The deployment and implementation of context-aware applications is currently not trivial due to the potentially large size of the ontology. In Conon [START_REF] Wang | Ontology based context modeling and reasoning using OWL[END_REF], the context is separated in two distinct parts: the first one is related to concepts describing the system at a high level of abstraction, whereas the second one is responsible for the description of domain specific concepts. For example, the generic concept "IndoorSpace" is translated in "Building", "Room", and "Corridor " in the home automation domain. This hierarchical decomposition is found in [START_REF] Gouin-Vallerand | A Macro and Micro Context Awareness Model for the Provision of Services in Smart Spaces[END_REF] which propose to separate micro/macro context to facilitate the implementation of such an application. The proposed approach deals with the relative complexity of context information processing in distributing reasoning process on Devices Nodes which have knowledge of micro context without loosing the benefits of macro information processing with Coordinator Nodes.

These works are based on sensors data, from which context information and knowledge can be inferred via logic reasoning. Our context-aware approach seeks to infer information usually sensored from the analysis of control logs. From this point of view, the ontology and rules presented in this paper differs noticeably from existing works. Previous works from [START_REF] Truong | Proactive remote healthcare based on multimedia and home automation services[END_REF] highlight the interest of monitoring activities based on control logs to limit the use of sensors. The proposed monitoring system is based on a frequency analysis in which a service is considered as a repeating pattern over time. Based on the hypothesis of very regular activities, the goal is to automatically extract and to propose daily living scenarios to wheelchairbound users with limited moving capabilities. However, this hypothesis is not always verified in practice [START_REF] De Lamotte | Quatra: Final report[END_REF]. For example, the user needs to switch on/off lights in rooms leading to the bathroom at night. This kind of interaction can be requested from time to time, in a non-regular manner: at any hour of the night, any day of the week, etc. In this paper, instead of a frequency analysis based on the assumption of regular activities, we proposed a context-aware system based on a semantic analysis of each control log. From the ontology and rules defined to model user-system interaction, logic reasoning provides the user adapted services without any training steps and using as few sensors as possible.

Conclusion

The considerable increase of dependent people brings out important needs of monitoring to ensure a safe and secure independent living and adapted services to compensate for their potential disability. In this context, this paper has presented a non-intrusive context-aware system based on ontology for ambient assisted living in smart home. The context-awareness of such systems is not based on sensors use, but on the semantic analysis of logs recorded each time the user asks for a service. An ontology with associated rules allows the inference engine to extract information usually sensored. As the design and the deployment of ontology is still not trivial, a model-driven process is applied for populating the ontology and a software architecture is defined for runtime. Finally, a case study based on the Domus Smart Home dataset validates our approach.

The person tracking application and results on the entire dataset are encouraging: depending on the context of use, it is possible to i) limit the use of sensors (low QoS ) or ii) increase the information redundancy (high QoS ). To do so, our approach considers control logs and logic inferences i) instead or ii) in addition to switches contacts and motion detectors logs. However, to achieve a complete monitoring system over time, this approach should consider both user feedback and sensors whose information can not be inferred from home automation logs: information coming from biomedical sensors (e.g. blood pressure, core body temperature, etc.) or accurate information meters (e.g. temperature, luminosity, etc.). These kind of information could also be translated in assertions in our context-aware system.

Service delivery is also an interesting contribution to allow users faster and easier access to services from their control interfaces. Some people with disabilities including motion impairment (e.g. tetraplegia) can only control their home automation system by means of an automatic scrolling interface. In such restrictive situations, the selection of the desired service is very time consuming. Benefiting from the model-driven process presented in this paper, a full-scale deployment is planned in Kerpape Living Lab to evaluate services acceptability and how it can assist patients in their everyday life 10 . These apartments have been built to evaluate the autonomy of patients before leaving the center. As an experimental validation is essential to evaluate service delivery, we are currently collecting home automation logs with motion-impaired users.

Fig. 2 .

 2 Fig. 2. Excerpt of the ontology taxonomy (T-Box with concepts only)

Fig. 3 .

 3 Fig. 3. Excerpt of the ontology for user location, user intention and service delivery

Fig. 4 .

 4 Fig. 4. Model-driven flow for control and context-awareness parts

Fig. 5 . 5 .

 55 Fig. 5. Software architecture with application and modeling part

Table 2 .

 2 Dataset (excerpt) of series 1 / user 1 / day 2 considering control logs

Table 3 .

 3 Results of the context-aware system on the entire Domus dataset

Related works are discussed in section 6 to clarify our contribution.

In accordance to the Model-Driven Architecture (MDA) specification[START_REF] Miller | Model Driven Architecture (MDA). Object Management Group[END_REF].

Experiments have been conducted with students and occupational therapists with developed tools (see demo here: http://www.youtube.com/watch?v=B4FJcdDBa6Q).

The Ontology Definition Metamodel[START_REF] Consulting | Ontology Definition Metamodel[END_REF] is a specification to unify MDA and Ontology, making OWL conformed with the three modeling levels of the MDA specification.

Each home automation log is divided into isAskingForService and isUsing properties.

If IR01 and IR02 give precise user location in one part of the kitchen, our contextaware approach can infer it at a room scale when a local activation device is used.

[START_REF] Consulting | Ontology Definition Metamodel[END_REF] QoS is related to data acquisition frequency and reliability of inferred knowledge.

see a home automation application with an automatic scrolling and a brain-computer interface used in one apartment of Kerpape Living Lab here: http://www.youtube. com/watch?v=0grL24U5duM

⋆ This work is financed by the Brittany Region and is the result of a cooperation supported by Kerpape MFRRC (http://www.kerpape.mutualite56.fr/).


O De La Barrière

H Ben Ahmed

M Gabsi

M Lobue

2D analytical airgap field model of an inset permanent magnet synchronous machine, taking into account the slotting effect

Keywords: analytical modeling, torque calculation, airgap flux, finite element verification

In this paper, an analytical method for the computation of the electromagnetic field in the airgap of an inset permanent magnet synchronous machine is proposed. The originality lies in the fact that both the stator and rotor saliencies have been taken into account in the analytical procedure. The purpose in to compute both quickly and accurately the magnetic fluxes and instantaneous torque of the machine, in order to minimize torque ripples in an optimal pre-design process. It is shown that the instantaneous torque is in good agreement with finite element results, for a shorter computation time. The model presented here can be used when the 2D approximation of the machine is valid, which is the case for most of machines' topologies.

I. INTRODUCTION

LECTRIC MOTORS for electric and hybrid traction applications have received a lot of interest in recent years [START_REF] Wang | Comparison of Five Topologies for an Interior Permanent-Magnet Machine for a Hybrid Electric Vehicle[END_REF][2] [3][4], due to [START_REF] Boughrara | Magnetic Field Analysis of Inset and Surface-Mounted Permanent-Magnet Synchronous Motors Using Schwarz--Christoffel Transformation[END_REF] the necessity of reducing fuel consumption, both for economic and environmental reasons. Permanent magnet synchronous 20 machines are promising in these demanding applications, because of their potential high torque density and efficiency [START_REF] Javadi | Design and analysis of 42-V coreless axial-flux permanent-magnet generators for automotive applications[END_REF]. The design of the electric machine for such embedded application can be difficult, due to the necessity of reaching both a high torque density and a compact size, which are contradictory optimization criteria [START_REF] Ashabani | Multiobjective shape optimization of segmented pole permanent-magnet synchronous machines with improved torque characteristics[END_REF] [START_REF] Liu | Optimization of an 80 kW Radial-Radial Flux Compound-Structure Permanent-Magnet Synchronous Machine Used for HEVs[END_REF]. The torque quality, i.e. small torque ripple, is also an important parameter in such applications [START_REF] Kwack | Optimal stator design of interior permanent magnet motor to reduce torque ripple using the level set method[END_REF].

The torque ripple computation is by no way an obvious problem. Most authors use the well known finite element method [9] [10]. This method presents the advantage of taking into account all the geometrical complexity of the machine, as well as magnetic saturation. However, the computation time can be prohibitive in a first design optimization process, since the instantaneous torque computation requires a time-stepping analysis [START_REF] Li | A new efficient permanent-magnet vernier machine for wind power generation[END_REF]. Therefore, to include the cogging torque minimization criteria in the pre-design optimization process, it is better to develop an analytical approach [START_REF] Zarko | Analytical solution for cogging torque in surface permanent-magnet motors using conformal mapping[END_REF]. Most of the analytical approaches proposed in the literature are done for surface mounted permanent magnet machines [START_REF] Zhu | An accurate subdomain model for magnetic field computation in slotted surface-mounted permanent-magnet machines[END_REF]. However, more complex permanent magnet rotor topologies, such as the inset permanent magnet structure, can be more suitable for high speed and flux weakening purpose [START_REF] Xu | New Axial Laminated-Structure Flux-Switching Permanent Magnet Machine With 6/7 Poles[END_REF]. This article proposes an exact analytical computation of the magnetic field in the airgap of a synchronous permanent magnet machine having an inset permanent magnet topology. The double saliency of the machine (i.e. the stator slotting and the rotor ferromagnetic teeth) is taken into account analytically. In a first part, a literature review on the subject is proposed. Then the analytical model is presented, with the detailed mathematical procedure. To finish, some results about the airgap induction and instantaneous torques are given and compared with finite element analysis. A good agreement, for a much smaller computation time, has been found.

II. OVERVIEW OF THE PREVIOUS WORK ON THE TOPIC

The problem of modeling the stator slotting effect of electrical permanent magnet machines has been dealt with in several references of the literature. In [START_REF] Zhu | Instantaneous magnetic field distribution in brushless permanent magnet DC motors. III. Effect of stator slotting[END_REF], the stator slotting effect has been computed in a radial flux machine by derivation of a relative modulation function. This modulation function has been computed using the same conformal map than the one of Carter in his wellknown paper [START_REF] Carter | Air-gap induction[END_REF]. However, this method only allows computing the airgap induction radial component, which is a problem for the tangential force computation by the Maxwell stress tensor method [START_REF] Gysen | Modeling of flux switching permanent magnet machines with Fourier analysis[END_REF]. Moreover, the magnetic interaction between the slots is not taken into account. Other authors [START_REF] Markovic | Determination of tooth cogging force in a hard-disk brushless DC motor[END_REF] [START_REF] Boughrara | Magnetic Field Analysis of Inset and Surface-Mounted Permanent-Magnet Synchronous Motors Using Schwarz--Christoffel Transformation[END_REF] have corrected these drawbacks of the conformal mapping method by using the numerical Schwartz-Christoffel Matlab toolbox [START_REF]Numerical conformal mapping using cross-ratio and Delaunay triangulation[END_REF], which permits to take into account rigorously the whole magnetic domain. However, this method requires the numerical computation of the conformal map coefficients, which is time consuming, because non linear transcendental equations must be solved numerically. As an alternative to conformal mapping, some authors prefer the subdomain method [START_REF] Rahideh | Analytical magnetic field distribution of slotless brushless machines with inset permanent magnets[END_REF] [START_REF] De La Barrière | Two-dimensional analytical field model of an inset permanent synchronous machine[END_REF] [START_REF] Lubin | Exact analytical method for magnetic field computation in the air gap of cylindrical electrical machines considering slotting effects[END_REF], which permits to find an exact solution of the potential in the airgap of the machine, under the form of Fourier series. The idea of this method is to decompose the machine airgap into elementary rectangular domains. The Laplace or Poisson equation is solved in each domain, and the domains are linked using classical boundary conditions in electromagnetism. This method generally requires a numerical inversion of a linear Cramer system, which can be done quite easily and fastly compared to the transcendental equations of the conformal map. Some papers [START_REF] Rahideh | Analytical magnetic field distribution of slotless brushless machines with inset permanent magnets[END_REF] [START_REF] De La Barrière | Two-dimensional analytical field model of an inset permanent synchronous machine[END_REF] apply this method to the rotor saliency for inserted permanent magnet machines. Other authors [START_REF] Zhu | An accurate subdomain model for magnetic field computation in slotted surface-mounted permanent-magnet machines[END_REF], on the contrary, use this method to compute the stator slotting effect of surface mounted permanent magnet machines, which can be with a double excitation winding [START_REF] Bali | Analytical modeling of open circuit magnetic field in wound field and series double excitation synchronous machines[END_REF]. Other authors have even managed to take into account analytically a more complex shape for the slots than a E simple rectangular [START_REF] Lubin | 2D exact analytical model for surface-mounted permanent-magnet motors with semi-closed slots[END_REF]. In each case, some work remains to be done to take into account both the stator slotting effect, and the salient pole rotor topology with inset magnets. This is the purpose of the further analytical development. Thanks to the model presented in this work, it will be possible to include the torque ripple minimization criterion from early predesign stages, without increasing too much the computation time. Indeed, the pre-design of hybrid vehicles is a hard task [START_REF] Lukic | Effects of drivetrain hybridization on fuel economy and dynamic performance of parallel hybrid electric vehicles[END_REF], since all the traction chain has to be taken into consideration. Thus, for computation time purposes, the model used for the actuator is often simplified. For example, in [START_REF] De La Barriere | Axial flux machine design for hybrid traction applications[END_REF], the slotting effect of the machine has been neglected. This can be a problem, because the torque quality (i.e. low torque ripples) is also an important point in such applications.

III. SYMPLIFYING HYPOTHESIS FOR THE ANALYTICAL MODELLING

A. 2D assumption and machine's topology

As done in all the previous references about analytical modeling, the machine is considered as bi-dimensional (2D), i.e. the edges effects are neglected. This means that the machine is considered as infinite perpendicularly to the plane of representation. This is typically the case in conventional radial flux machines [START_REF] Boughrara | Magnetic Field Analysis of Inset and Surface-Mounted Permanent-Magnet Synchronous Motors Using Schwarz--Christoffel Transformation[END_REF] in which the axial length H is important compared to the pole pitch characteristic dimension τ. Notice that this 2D assumption can be also valid for others topologies than conventional radial flux machines. For example, for axial flux machines [START_REF] Marignetti | Comparison of axial flux PM synchronous machines with different rotor back cores[END_REF], the 2D assumption is also possible, provided the difference H between the external and internal radii is important compared to the pole pitch τ at the mean radius. More precisely, for an axial flux machine, it is first assumed that the electromagnetic phenomena happening at the mean radius are sufficient to describe the behavior of the machine. The machine is made equivalent to a fictitious linear actuator having the same geometry than the real axial flux machine at its mean radius. Then the thickness H of this fictitious linear actuator is assumed to be very important compared to its pole pitch τ. This 2D approximation for the axial flux machine is summed up in Fig. 1. Then, under the 2D assumption, for both a radial and axial flux machine, the problem to study can be represented in Fig. 2 (the main geometrical parameters are given in this figure, and are recalled in Table I). For the radial flux machine, such a representation implies that the curvature of the machine has been neglected, which is valid if the mean airgap radius (i.e. the mid airgap radius) is important compared to the dimensions along the y axis. This is assumed to be the case in the following of this paper. A discussion about this simplifying assumption is given in the Appendix 1. -No saturation effect is taken into account. The magnetic permeability of the ferromagnetic parts is assumed to be infinite. This assumption is also discussed in Appendix I.

-

It is assumed that continuous rate has been reached, and that the machine is driven with sinusoidal currents waves. -Finally, the current density in each slot is assumed to be uniform in the entire slot surface.

IV. THE CALCULATION METHOD

A. Definition of the problem Fig. 3 shows the machine over one pole pitch. For simplicity, the number of stator slots per pole is chosen to equal to 3.

Fig. 3: Representation of the problem to solve, with double (i.e. rotor and stator) saliencies Due to the 2D approximation, the problem can be reduced to the computation of the following vector potential:

( , ) ( , ). z A x y A x y u = (1)

STATOR YOKE ROTOR YOKE

Rotor tooth

0 dA dn = Stator Tooth 0 dA dn =
Each slot is characterized by a number l, where l can evolve between 1 and 3. The current density in the slot l is called j l , and is assumed to be uniform over the slot. A first difficulty arises from the fact that the geometry of the problem and its solution depend on the rotor displacement X d . The stator slots zone is called Zone I, the airgap is called Zone II, and Zone III is the magnet zone (Fig. 3).

B. Vector potential mathematical expressions

In Zone III, due to the magnet (whose magnetization function is called M(x)), the potential is given by a Poisson equation:

( ) 0 III M A x ∂ ∆ = ∂ µ (2)
As done in [START_REF] De La Barrière | Two-dimensional analytical field model of an inset permanent synchronous machine[END_REF], the magnetization can be developed into a Fourier series:

1 2 ( ) .cos R n R n a M X M n X π β τ +∞ =      =        ∑ (3)
The coefficients M n are given by: 4) are satisfied if the following expression for the potential is chosen, in the main coordinate system (0xy):

1 2 0 4 1 sin 2 a r n a B M n n β π π µ β       =         (4)
( ) ( ) ( ) ( ) 0 2 1 2 0 2 1 2 , cosh cos cos III n a n a n d a n d a A x y c c n y M n x X n n x X π β τ β τ π µ π β τ π β τ ≥ ≥    = +          - -          ⋅ -     ∑ ∑ (5)
The unknown coefficients c n (n≥0) are determined in the next sections.

For Zone II, a Laplace Equation must be solved:

( ) 0 II A ∆ = (6)
Moreover, in this zone, in the airgap (Zone II), an anti-periodicity condition has to be taken into account:

( ) ( ) ( ) ( ) , , II II A x y A x y τ + = - (7)
This implies that the Fourier series only present odd terms:

( ) ( ) ( ) ( ) ( ) ( ) 1 1 , cos 2 1 ' sin 2 1 II k k k k A x y b y k x b y k x π τ π τ ≥ ≥   = -       + -     ∑ ∑ (8)
The Laplace equation ( 6) implies:

ROTOR Y R h a X R 0 R ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 1 2 3 4 cosh 2 1 sinh 2 1 cosh 2 1 sinh 2 1 k k k k k k b y b k y b k y b y b k y b k y π π τ τ π π τ τ     = - + -             ′ = - + -         (9)
For Zone I, a generic slot l is considered, with its local coordinate system (0 l ,X l ,Y l ) (Fig. 5). Due to the currents, the equation in the slot is a Poisson equation:

( ) 0 . l l j µ Ι, ∆Α = - (10)
The homogeneous solution can be found by solving the corresponding Laplace equation (

( ) 0 l h Ι, ∆Α = ): ( ) ( ) ( ) ( ) ( ) , 0 1
, cosh cos

I l l l l l m l e l h e e m A X Y a a m Y h m X w w π π ≥     = + -         ∑ (11) A particular solution to the Poisson equation ( ) , I l p l j µ 0 ∆Α = -
must be found. The following quadratic form can be proposed as a particular solution in the slot:

( ) ( ) ( ) 2 , 1 , 2
I l p l l l l e A X Y j Y h µ 0 = - - (12)
Fig. 5: Zoom on a stator slot In the main (0xy) axis system, it can be found:

( ) ( ) ( ) ( ) ( ) ( ) ( )( ) ( ) ( ) ( ) 2 , 0 1 1 , cosh cos 1 2 I l l l m a e e d l a e e e m

A x y a a m y h e h m x l w w j y h e h w w

π π µ 0 ≥     = + - + + -- + - - + +         ∑ (13)

C. Boundary conditions between all the different domains

Starting from the mathematical expressions of the previous section, it is possible to write the boundary conditions between all the zones, and determine the series' coefficients. As known in electromagnetic theory, the boundary conditions to satisfy are: the continuity of the vector potential the continuity of its normal derivative The following notations are introduced:

- ( ) ( )
( ) cosh 2 1   = - +     k II a c K k h e π τ - ( ) ( )
( ) sinh 2 1   = - +     k II a s K k h e π τ STATOR - ( ) ( ) cosh 2 1   ′ = -     k II a c K k h π τ - ( ) ( ) sinh 2 1   ′ = -     k II a s K k h π τ - ( ) cosh   =     m I e c e K m h w π - ( ) 2 cosh   =     n III a c a K n h π β τ 1)

Boundary conditions between Zone I and Zone II.

The continuity of the potential between Zone I and Zone II can be written, in mathematical terms, as:

( )( ) ( )( ) ( ) ( ) ( ) ( ) 1 , 1 , 1 ,
, ,

s q e d e d e l I l II a a x l w w l w w w A x y h e A x y h e =   ∀ ∈ - + - + +   = + = = + ∪ (14)
Using the boundary integral method explained in [START_REF] Sargos | Calcul analytique du champ engendré par des aimants dans l'entrefer d'une machine à rotor denté[END_REF][30], the following equation can be obtained:

( ) ( ) ( ) ( ) ( ) 1 2 ( ) 2 0 , 1 3 4 ( ) , 1 1 2 0 ≥ ≥   - = +       + +     ∑ ∑ k k k k l l l e c s k k C k k l c s k k S k k a j h b K b K I b K b K I µ (15)
And, for each slot l: Concerning the normal derivative of the potential between Zone I and Zone II, it is possible to write: [START_REF] Sargos | Calcul analytique du champ engendré par des aimants dans l'entrefer d'une machine à rotor denté[END_REF][30], it can be found:

( ) ( ) ( ) ( ) ( ) ( ) 1 2 3 4 ( ) ( ) ,( , ) , ( , ) 1 1 1, 1 ≥
≥ ∀ ≥ = ⋅         + + +               ∑ ∑ m k k k k l m I c l l c s c s k k C m k k k S m k k k m a K b K b K b K b K α α (
( )( ) ( )( ) ( ) ( ) [ ] ( )( ) ( )( ) ( ) 1 1 1 , 1 , 0, 1 , 1 , 0 =
= + = + = = +   ∀ ∈ - + - + +       =                   ∀ ∈ - - + - + +           =      
( ) ( ) ( ) 1 2 ( ) ( ) ( ) ,( , ) , 1 1 1 1 +2 2 1 0 = ≥ ∀ ≥ -     = -   -     ∑ ∑ k k s m c s k k q l I l l e e l m C m k C k c l m k b K b K h w j ma K I k µ α π (18)
And also:

( ) ( ) ( ) 3 4 ( ) ( ) ( ) ,( , ) , 1 1 1 1 +2 2 1 0 = ≥ ∀ ≥ -     = -   -     ∑ ∑ k k s m c s k k q l I l l e e l m S m k S k c l m k b K b K h w j ma K I k µ α π (19)
2) Boundary conditions between Zone II and Zone III. Using the same method as in the previous section for taking into account the boundary conditions, it is possible to get the following relation, for each rotor displacement X d , for what concerns the potential continuity between Zone II and Zone III:

( ) ( ) ( ) ( ) 1 2 3 4 ( ) ( ) ( ) ( ) 0 , , 1 1 ≥ ≥     ′ ′ ′ ′ = + + +         ∑ ∑ k k k k II II II II C k S k k k k k c s c s k k c b K b K J b K b K J
(20) And also:

( ) ( ) ( ) ( ) 1 2 ( ) ( ) ,( , ) ( ) 1 3 4 ( ) ( ) 2 ,( , ) 0 1 1, 1 . k k n k k II II n C n k k k c s III c k II II a n S n k k k c s k n c b K b K K M b K b K n β β τ β µ π ≥ ≥ ∀ ≥     ′ ′ = +          ′ ′ + + -      ∑ ∑ (21)
The definitions and expressions of the normalized numbers , For the continuity of the normal derivative of the potential between Zones II and III, the two following equations are obtained:

C k J , , S k J , , ( ,
( ) ( ) 1 2 ( ) ( ) ( ) ,( , ) 1 1 1 2 1 ≥ ∀ ≥ ′ ′ - = -∑ k k n II II III n C n k k k c s c n k b K b K nc K k β (22) ( ) ( ) 3 4 ( ) ( ) ( ) ,( , ) 1 1 1 2 1 ≥ ∀ ≥ ′ ′ - = -∑ k k n II II III n S n k k k c s c n k b K b K nc K k β ( 23
)
D. Solution of the problem. 1) Matrix formulation All the series are limited to a finite number of terms, called N. This number is computed in the next part, so as to offer the best compromise between the computation time, and the accuracy. The unknown parameters are, using the vector notation:

( ) ( ) ( ) 1 1, , l l s l N a l q a a       ∀ ∈ =           ⋮ , ( ) ( ) ( ) 1 1 1 1 N b b b       =           ⋮ , ( ) ( ) ( ) 2 1 2 2 N b b b       =           ⋮ , ( ) ( ) ( ) 3 1 3 3 N b b b       =           ⋮ ( ) ( ) ( ) 4 1 4 4 N b b b       =           ⋮ , and [ ] 1 N c c c     =       ⋮ .
Equation ( 16) can be written for all the numbers m between 1 and N. So the following matrix equation is obtained, for all slot number l between 1 and q s :

[ ] ( ) [ ] [ ] ( ) [ ] ( ) ( ) [ ] ( ) [ ] ( ) ( ) 1 2 ( ) 2 3 1 3 4 ( ) 2 3 1, , l C l s l S D b D b l q a D D b D b α α         +           ∀ ∈ =             + +         (24)
The following square matrixes for each slot are introduced (see appendix II):

[ ] 2 ( ) ( ) ,( , ) ( , ) 1, l l C C m k m k N α α ∈     =     and [ ] 2 ( ) ( ) ,( , ) ( , ) 1, l l S S m k m k N α α ∈     =     .
The other diagonal matrixes are also defined in Appendix II.

Equations ( 18) and [START_REF] Boughrara | Magnetic Field Analysis of Inset and Surface-Mounted Permanent-Magnet Synchronous Motors Using Schwarz--Christoffel Transformation[END_REF], for each k between 1 and N, lead to the following matrix relations:

[ ] ( ) [ ] ( ) [ ] [ ] ( ) 1 2 ( ) 2 3 4 1 ( ) 1 . . 2 s s q T l l C l q l e e l C l D b D b INV D a h w j I α µ π = 0 =         - = -             +   ∑ ∑ (25) [ ] ( ) [ ] ( ) [ ] [ ] ( ) 3 4 ( ) 2 3 4 1 ( ) 1 .
.

2 s s q T l l S l q l e e l S l D b D b INV D a h w j I α µ π = 0 =         - = -             +   ∑ ∑ (26)
The matrix [INV], as well as the following matrixes, are also defined in Appendix II:

[ ] ( ) ( ) , 1, l l C C k k N I I ∈     =    
, and

[ ] ( ) ( ) , 1, l l S S k k N I I ∈     =    
Thanks to [START_REF] Carter | Air-gap induction[END_REF], it can be obtained:

[ ] ( ) [ ] [ ] ( ) [ ] ( ) ( ) { [ ] ( ) [ ] ( ) ( )} 1 2 ( ) 1 2 3 3 4 ( ) 2 3
1, ,

l l s C l S l q a D D b D b D b D b α α         ∀ ∈ = +               + +       (27) And, from (21): [ ] [ ] [ ] ( ) [ ] ( ) ( ) { [ ] ( ) [ ] ( ) ( ) [ ] } 1 2 5 6 7 3 4 6 7 C S c D D b D b D b D b V β β       = +             + + -       (28)
The matrix

[ ] 2 ,( , ) ( , ) 1, C C n k n k N β β ∈     =     and [ ] 2 ,( , ) ( , ) 1, S S n k n k N β β ∈     =   
 are also given in Appendix II, as well as the vector

[ ] [ ] [ ] 1, n n N V V ∈ =
and the other diagonal matrixes.

Equations ( 22) and ( 23) lead to the following matrix equations:

[ ] ( ) [ ] ( ) [ ] [ ][ ] 1 2 6 7 8
.

.

T C D b D b INV D c β       - =       (29) [ ] ( ) [ ] ( ) [ ] [ ][ ] 3 4 6 7 8
.

.

T S D b D b INV D c β       - =       (30)
2) Numerical computation of the problem.

Recalling the equations' system to solve:

( ) [ ] [ ] ( ) [ ] ( ) ( ) [ ] ( ) [ ] ( ) ( ) [ ] ( ) [ ] ( ) [ ] [ ] ( ) [ ] ( ) [ ] ( ) [ ] [ ] 1 2 ( ) 2 3 1 3 4 ( ) 2 3 1 2 ( ) ( ) 2 3 4 1 1 3 4 ( ) 2 3 4 1 1, , . . 2 .
.

s s s l C l s l S q q T l l l e e l C C l l q T l S l D b D b l q a D D b D b h w D b D b INV D a j I D b D b INV D 0 = = =         +           ∀ ∈ =               + +                   - = - +                   - = -       ∑ ∑ ∑ α α µ α π α ( ) [ ] [ ] [ ] ( ) [ ] ( ) ( ) [ ] ( ) [ ] ( ) ( ) [ ] { } [ ] ( ) [ ] ( ) [ ] [ ][ ] [ ] ( ) [ ] ( ) [ ] [ ][ ]
b D b V D b D b INV D c D b D b INV D c 0 =                +                      = + + + -                      - =               - =         ∑ µ π β β β β (31)
If q s =3, a system of 8 equations with 8 unknown vectors must be solved. Equations ( 25) and ( 26) are substituted in Equation [START_REF] Bali | Analytical modeling of open circuit magnetic field in wound field and series double excitation synchronous machines[END_REF].

In order to reduce the number of equations, the following process is proposed: in a first stage, by an analytical substitution, we can express this system with only the unknown series' coefficients vectors: [b (1) ][b (2) ][b (3) ][b (4) ].

After algebraic substitutions, the system can be written under the form:

[ ] [ ] [ ] . M b P = (32)
Where the unknown vector

[b] is [ ] ( ) ( ) ( ) ( ) 1 2 3 4 b b b b b                   =                  
. The matrix [M] and the vector [P] are obtained by linear combinations of the matrixes introduced in the previous section.

To solve the linear system (32), it is necessary to perform a numerical inversion of the matrix [M]. It is also necessary to choose the number N of harmonics. For that, the following numerical procedure may be proposed: firstly, it is possible to compute the inversion of the system for N=1, which is very fast, since the size of the matrix [M] is just 4. Then the following algorithm is applied. Suppose that we have just computed the solution for a certain number N of harmonics. For this number, the matrix

[M] is called [M N ], [V] is called [V N ], and the solution [b N ] is: [ ] ( ) ( ) ( ) ( ) 1 2 3 4 T T T T T N N N N N b b b b b           =                    
To compute the vector solution [b N+1 ] for N+1 harmonics, which is solution of the equation

[M N+1 ].[b N+1 ]=[V N+1 ],
a Gauss-Siedel algorithm can be used. In this inversion method, a first initialisation of the solution is required. Instead of taking the zero vector for the initialisation, we use in the initial step the solution of the previous system with N harmonics, more precisely:

( ) ( ) ( ) ( ) 1 2 3 4 1,0 0 0 0 0 T T T T T N N N N N b b b b b +             =                      
So the number of iterations of the Gauss-Seidel process is reduced in order to get the solution [b N+1 ]. We stop increasing the number N when the difference of the solutions between two iterations becomes as small as wanted (in practice less than 1% in relative value). So the number N of harmonics is not fixed, but is optimized for each machine's geometry.

V. COMPARISON WITH FINITE ELEMENTS

For the finite element verification, as well as for the torque computations, the geometry presented in Table II is used. Moreover, for the torque computations, the mean airgap radius (i.e. the mid airgap radius R m of the radial flux machine, or the mean active radius of the axial flux machine) is assumed to be 0.1m (see Appendix I for a detailed discussion about this point). In the finite element code, in order to be consistent with the assumption of infinite iron permeability made for the analytical model, a high relative permeability is assigned to the ferromagnetic parts (in practice, this relative permeability is taken equal to 1000). A discussion about the consequences of this linear model is provided in Appendix I: the results are compared with ones obtained using non linear B-H experimental relationships. A. The magnets' field.

The no-load airgap magnetic field (created by the magnets, when all the stator currents are zero) has been computed, both analytically and by finite elements. On Fig. 6, a contour [AB] in the middle of the airgap is chosen, on which we shall plot the induction caused by the magnets (all the current densities j 1 =j 2 =j 3 =0). This segment length is equal to the pole pitch τ. The results are shown on Fig. 7, Fig. 8 and Fig. 9, respectively for the potential, the B x and the B y induction components. The agreement is satisfying. The results of the field caused by the stator currents (in the computation procedure, B r =0) are now analyzed. According to the Park theory, a d axis current is applied. Since the system is linear, the calculus is performed for a current of just 1A in each slot. So, on Fig. 6, current densities of j 1 =1A/S slot (S slot is the slot surface), j 2 =j 1 /2, and, j 3 =-j 1 /2 have been applied. For the q axis current, the rotor undergoes a translation of a half pole pitch in comparison with the d axis calculation. The results for the d axis currents are shown on Fig. 10, Fig. 11 and Fig. 12, in good accordance with finite elements. The results obtained using the two methods are in good agreement. The time necessary for the analytical computation is about 30% less than the time necessary for the finite element calculation.

VI. THE INSTANTANEOUS TORQUE COMPUTATION

A. The calculation method

The analytical model can be used for the computation of the instantaneous torque. The analytical computation is performed as the following: the Maxwell stress tensor formula is applied on the [AB] segment of Fig. 6 located in the middle of the airgap, giving the tangential effort on this segment as:

[ ] t x y AB F H B dx = ∫ (33)
This analytical calculation is repeated for each elementary position of the rotor. The torque calculation with finite elements is based on the same principles, but the field and induction and then the Maxwell stress tensor are computed numerically on the [AB] segment. The spatial period of the torque is a tooth pitch. The tooth pitch has been decomposed in ten elementary positions (both for finite element or analytical computation). For each of them, the calculation has been carried out, and a cubic interpolation has been applied between the ten points. This strategy permits to save computation time, without decreasing the model accuracy.

In the Park formalism, a q current of 800A per slot and a d current per slot of -100A are applied. The total torque is shown on Fig. 13. The analytical and finite elements calculations are in good agreement: the mean torque is around 170Nm for the two methods.

Fig. 13: The total torque over a tooth pich The cogging torque is given on Fig. 14, (it is computed with zero currents). Its mean value is zero (as expected), although its variations are quite high (50Nm). The discrepancy of 8% on the maximum value can be explained by the fact that the instantaneous torque computation is very sensitive to the number N of harmonics chosen for the analytical computation. Indeed, the torque is obtained by multiplying two derivatives of the potential: the normal induction B y and the tangential field H x =B x /µ 0 . These derivatives are quite irregular signals (see the previous section), and an important number of harmonics can be required to represent them with a good precision. Here, as explained in Section IV.D.2, we stopped increasing the number N of harmonics if the [b] vectors between two iterations do not present more than 1% of difference in relative value. If we choose a lowest error, the torques computed by finite elements and analytically become closer (a difference of 4% is observed for a relative error of 0.3%). However, this solution does not appear as optimal, because it increases much the computation time, and then cannot be used for a pre-design process. We believe that a

ANA. FEM

difference of about 10% is acceptable in a pre-design process stage. The most important is however to get an exact mean value, which is the case, even for a small number of harmonics. Fig. 14: The cogging torque over a tooth pich The reluctant torque, which is due only to the stator currents (the remanent induction of the magnets is put to zero), is presented in Fig. 15. The fluctuation is quite important, but the mean value remain weak (10Nm) compared to the total torque, mainly due to the hybrid torque component. This one is obtained by subtracting to the total torque the cogging and the reluctant torque (Fig. 16).

B. Discussion

It has been shown that the analytical method gives results in good agreement with finite elements. Nevertheless, it could be interesting to illustrate the limit of this method.

The main limit of the method consists in the fact that saturation of the machine's ferromagnetic parts is not taken into account. So torque calculations are performed in this section using the finite element method including ferromagnetic material saturation (the FeSi B-H curve given in Appendix I), and the results are compared with the analytical model. This calculation permits to understand the consequences of this limitation on the torque computation.

ANA. FEM

So a high q axis current (up to 500A) is applied, so as to create saturation, and this case is studied using finite elements computations. The mean torque is computed over a tooth pitch, as in the previous section. The torque curve taking into account saturation is compared with the one derived from our analytical linear model. This work has been done for two values of d axis current. The first is for a zero d axis current (Fig. 17), and the second is for a d axis current per slot of nI d =-600A (Fig. 18). The first fact to notice is that the saturation phenomena reduces the total flux in the machine, especially when the q axis current is high Therefore the torque is not as high as in the non saturated case. Another point is the effect of the d axis current. A strong d axis current permits to increase the torque significantly, since it allows a reluctant torque component. But it also permits a reduction of the effect of saturation (Fig. 18), because the total flux of the machine is weakened. It could be shown in the same way that flux weakening also decreases the iron losses.

VII. CONCLUSION

An analytical method was developed in order to take into account the stator and rotor saliency of the machine, for the computation of the magnetic induction in the machine airgap. This permits to compute precisely the machine's torque. In a coming paper, we are going to use this model for the computation of the inductions in the ferromagnetic parts of the machine, so as to be able to evaluate the iron losses, taking into account the stator leakage flux. This approximation is justified if the mean airgap radius R m is important compared to the other radial dimensions (magnet thickness, airgap width, slot dept…). To understand to what extend the approximation can be done, finite element calculations have been carried out with the radial flux machine, and with the equivalent linear actuator, in order to compare the results. The machine geometry used in the comparison is given in Table III. The ferromagnetic parts are taken into account using a linear material of relative magnetic permeability equal to 1000.

Dimensions along the z direction

Active thickness of the machine along the z axis H 0.1m

Physical data

Remanent induction of the magnets (NdFeB) at room temperature Br 1.2T For example, the no-load flux is studied. The parameter of importance is the value of the ratio h a /R m between the magnet radial thickness and the mean airgap radius. The magnet thickness varies between 0.001m and 0.01m. The ratio between the no-load flux of the two actuators is shown in Fig. 20. As expected, the two fluxes are nearly the same for small magnet thicknesses with regard to the pole pitch, whereas some discrepancy appears by increasing the magnet thickness. However, the difference is never more than 7% in relative value. Moreover, a magnet thickness of 0.01m (for which the difference is the highest) is not useful, because, as shown in Fig. 21, no significant increase of the no-load flux appears between the two last points of the magnet thicknesses. Thus it can be concluded that the error on the no-load flux made by the development of the radial flux machine to the equivalent linear one is restricted to a few percent in relative value. However, if this error has to be corrected, some authors have proposed methods based on conformal mapping [START_REF] Sargos | Calcul analytique du champ engendré par des aimants dans l'entrefer d'une machine à rotor denté[END_REF], in order to rigorously transform the radial flux machine to an equivalent linear actuator. This conclusion allows us to compute the magnetic field in the linear equivalent actuator (in Cartesian coordinates), and then apply the results to the case of the radial flux machine. It can also be useful to check if the inductions are locally also the same for the two kinds of machines. For that, it is proposed to plot in Fig. 22 and Fig. 23 the normal and the tangential inductions in the middle of the airgap (dotted red line in Fig. 19). Indeed these inductions are involved in the cogging torque computation, and of course in the no-load flux calculation (the no-load flux is the integral of the normal induction component). It can be seen that the waveforms are nearly identical, which demonstrates the validity of the mean radius development also for the local induction computation.

B. The saturation effect

The analytical model implies that the ferromagnetic parts are assumed to have an infinite magnetic permeability. It is appropriate to compare the results provided by this model with the ones obtained with real non-linear B-H anhysteretic relationships experimentally measured in Iron-Silicon (FeSi) and Iron-Cobalt (FeCo) laminations (Fig. 24). The main advantage of Iron-Cobalt material is its high saturation magnetization (nearly 2.3T), and high initial permeability. Its drawback is its high cost and low resistivity, which implies high dynamic loss. The computations are then carried out using a non-linear finite element code. In this study, the same geometry as in Table II is adopted. The B y inductions obtained on the [AB] segment (Fig. 6), using the analytical and the non-linear finite element models, are plotted in Fig. 25 for no-load excitation conditions. The analytical model and the non-linear finite element model with the FeCo curve provide identical results, while a decrease of the airgap induction can be seen with the FeSi material. This decrease is responsible of an error of 8% on the total no-load flux between the analytical model, and the model using the FeSi curve. In Fig. 26, the same analysis has been repeated, being this time under field weakening conditions: in addition to the magnets' field, a d axis current has been added (nI d =-600A/m, see Section VI.B). The level and the waveform shape has been changed compared to the previous case in which only the no-load flux was considered. It is interesting to notice that the field weakening permits to reduce the saturation level of the machine, and then justify the use of the linear analytical model. Indeed the difference on the fluxes between the linear model and the FeSi saturated model is now 3% in relative value, whereas it was 8% under no-load condition. To conclude, it must be said that the saturation effect can decrease the airgap induction level, and then the torque (see Section VI for a detailed discussion). The flux weakening reduces the saturation effect. However, for the studied geometry, this effect is limited to a decrease of 10% of the airgap induction between the analytical model and the non-linear model using FeSi laminations. For FeCo laminations, due to the very high saturation level of this material, the saturation do not need to be taken into account, justifying the use of the analytical model. It must be noticed that the saturation effect can be taken into account using the model proposed by the author of [START_REF] Boules | Two-Dimensional Field Analysis of Cylindrical Machines with Permanent Magnet Excitation[END_REF]. The decrease of the airgap induction due to the saturation effect is modeled by a fictitious increase of the airgap thickness in the analytical computations.

    = - - - + + -         (37)
The different normalized coefficients for the potential continuity equalities between Zone I and Zone II are given by:

( ) ( ) ( ) , 2 2 sin 2 1 cos 2 1 2 1 2 2 C k C a d a J k k X k π π π β β τ     = - - + -         (38) ( ) ( ) ( )
S k C a d a J k k X k π π π β β τ     = - - + -         (39) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) (
- - -     - -             = + - - -         2 2     + + - + + -      2   2  (41)
π         =                   ⋱ ⋱ , [ ] ( ) ( ) 2 cosh 2 1 a D k h e π τ         = - +           ⋱ ⋱ [ ] ( ) ( ) 3 sinh 2 1 a D k h e π τ         = - +           ⋱ ⋱ , [ ] 4 cosh e e D m m h w π         =             ⋱ ⋱ , [ ] 1 2 1 INV k       = -       ⋱ ⋱ , [ ] 5 2 1 cosh a a D n h π β τ         =                   ⋱ ⋱ , 396 [ ] ( ) 6 cosh 2 1 a D k h π τ         = -           ⋱ ⋱
, and [ ] ( )

7 sinh 2 1 a D k h π τ         = -           ⋱ ⋱ , [ ] 8 2 cosh a a D n n h π β τ         =             ⋱ ⋱ .

Fig. 1 :

 1 Fig. 1: 2D approximation of the axial flux machine

Fig. 2 :

 2 Fig.2: 2D problem to study for both a radial and axial flux machine, and main geometrical parameters

Fig. 4 :

 4 Fig. 4: Zoom on the rotor It can be shown that the Poisson equation (2), as well as the homogeneous Neumann boundary conditions dA/dn=0 on the segments [BC], [DE] and [CD] (Fig.4) are satisfied if the following expression for the potential is chosen, in the main coordinate system (0xy):

  again the boundary integral method

  ) C m k β , and ,( , ) S m k β are given in appendix II.

Fig. 6 :Fig. 7 :Fig. 8 :Fig. 9 :

 6789 Fig. 6: The plot contour [AB] in the airgap

Fig. 10 :Fig. 11 :Fig. 12 :

 101112 Fig. 10: The potential caused by the d axis currents on the segment [AB]

Fig. 15 :Fig. 16 :

 1516 Fig. 15: The reluctant torque over a tooth pich

Fig. 17 :Fig. 18 :

 1718 Fig. 17: Mean torque with and without saturation (nId=0)

Fig. 19 :

 19 Fig. 19: Development of the radial flux machine to its mean radius to obtain the equivalent linear actuator

Fig. 20 :

 20 Fig. 20: Ratio between the no-load fluxes obtained in the radial flux machine and the equivalent linear machine

Fig. 21 :

 21 Fig. 21: No-load flux in the radial flux machine (in function of the magnet thickness ha)

Fig. 22 :

 22 Fig. 22: No-load induction (normal component) in the radial flux machine and in the equivalent linear one (for a magnet thickness ha=0.0035m) on a line located in the middle of the airgap (dotted red line of Fig. 19)

Fig. 23 :

 23 Fig. 23: No-load induction (tangential component) in the radial flux machine and in the equivalent linear one (for a magnet thickness ha=0.0035m) on a line located in the middle of the airgap (dotted red line of Fig. 19)

Fig. 24 :

 24 Fig. 24: Non-linear B-H anhysteretic relationships experimentally obtained in FeSi and FeCo laminations

Fig. 25 :

 25 Fig.25: By induction obtained on the [AB] segment (Fig.6) for the analytical model, and the non-linear finite element models using the FeCo material and the FeSi material (no-load excitation conditions)

Fig. 26 :

 26 Fig. 26: By induction obtained on the [AB] segment for the analytical model, and the non-linear finite element models using the FeCo material and the FeSi material (field weakening condition: superposition of the magnets' field and d axis current nId=-600A/m)

  are put under a matrix formulation, it is convenient to introduce the following diagonal matrix:

TABLE II GEOMETRY

 II OF THE MACHINE CHOSEN FOR THE INDUCTION AND TORQUE COMPUTATIONS (ANALYTICAL AND FINITE ELEMENT VERIFICATION)

	Notation	Numerical values
	ha	5mm
	e	2mm
	ecr	10mm
	ecs	10mm
	he	30mm
	τ	60mm
	βa1	0.6
	βa2	0.8
	we	10mm
	wd	10mm
	H	100mm
	Br	1.2T

  )

		β	,( , ) C n k	sin	C	n	2 1 k	β	a	2	π	cos	π	n n	2 2 1 2 1 a k k β τ	2 2 a a β β	d π X	2 n X a π β τ d	sin	C	n	2 1 k	β	a	2	π	cos	π	n n	2 2 1 2 1 a k k β τ	2 2 a a β β	d π X	2 n X a π β τ d
	β	S n k	C	n		k	β	a			π	sin	π	n n +	2 1 a k k β τ	a a β β	d π X	d β τ a n X π	-sin	C	n	2 1 k	β	a	π	sin	π	n n	2 1 a k k β τ	a a β β	d π X	d β τ a n X π

APPENDIX 1

In this appendix, the two fundamental simplifying assumptions of this work (no curvature effect, no saturation) are discussed.

A. The curvature effect

It is possible to use the model presented in this article for a radial flux machine, on the condition of neglecting the curvature effect. It means that the real machine is made equivalent to a linear actuator, as shown in Fig. 19. The pole pitch τ of the linear actuator is chosen equal to the pole pitch of the radial flux machine at is mean airgap radius R m (i.e. the mid airgap radius). All the other radial dimensions (magnet thickness h a , airgap width e, slot dept h e …) are the same for both actuators. The relative magnet opening on the pole pitch and the relative slot opening on the tooth pitch are also identical.

APPENDIX 2

The numbers α introduced in for the continuity conditions were defined the formulas:


Julien Fleureau1

Amar Kachenoura1

Jean-Claude Nunes1

Laurent Albera1

Lotfi Senhadji1

Amar Kachenoura
email: amar.kachenoura@univ-rennesl.fr

3A-EMD: A GENERALIZED APPROACH FOR MON OVARIATE AND MULTIVARIATE EMD

EMD is an emerging topic in signal processing research and is applied in various practical fields. Its recent ex tension to multivariate signals, motivated by the need to jointly analyze multi-channel signals, is an active topic of research. However, all the existing extensions specifically hold either mono-, bi-or tri-variate signals or require mul tiple projections that complexify the original process. In this communication, a novel EMD approach called 3A EMD is proposed. It is essentially based on the redefini tion of the mean envelope operator and allows, under cer tain conditions, a straightforward decomposition of mono variate and multivariate signals without any change in the core of the algorithm. A comparative study with classical monovariate and bivariate methods is presented and shows the competitiveness of 3A-EMD. A trivariate decomposi tion is also given to illustrate the extension of the proposed algorithm to any signal dimension, D > 2.

INTRODUCTION

Empirical Mode Decomposition (EMD) was originally in troduced in the late 1990s to study water surface wave evolution [START_REF] Huang | The em-pirical mode decomposition and the hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF]. EMD can be considered as an emerging technique in signal processing with a very important re search topics and development in various fields such as biomedical signal analysis [START_REF] Hualou | Application of the empirical mode decomposition to the analysis of esophageal manometric data in gastroesophageal re flux disease[END_REF], Hurst exponent estimation [START_REF] Rilling | Empir ical mode decomposition, fractional gaussian noise and hurst exponent estimation[END_REF], etc. EMD can be viewed as a method which splits a given signal into a set of AM-FM components. The latter components, referred to as the Intrinsic Mode Functions (lMFs), are estimated using an iterative procedure called sifting process. Let {set)} = {a o (t)} be a signal to de compose. Then the sifting process to estimate the n th IMF (1 � n � N) can be summarized as follows:

1. Initialization taking dn, o (t) = anl (t) for every t.

2. Computation of the mean envelope M( {dn,k(t)}).

3. Extraction of the details {dn,k + l (t)} = {dn,k (t)} M( {dn,k(t)}). [START_REF] Rilling | On empiri cal mode decomposition and its algorithms[END_REF]. Incrementation of k and go to step 2, repeat until {dn,k + l (t)} is designed as an IMF.

In Huang [START_REF] Huang | The em-pirical mode decomposition and the hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF], the mean envelope M( {dn,k (t)}) is cal culated as the half sum of the upper and the lower en-velopes, which are obtained by interpolation between the local maxima points and the local minima points of {dn,k(t)}, respectively.

Several studies dealing with EMD have been reported in the last decade: those addressed to improve the Huang's algorithm [START_REF] Rilling | On empiri cal mode decomposition and its algorithms[END_REF][START_REF] Delechelle | Empirical mode decomposition: an analytical approach for sift ing process[END_REF], and those proposing an extension of EMD to bi-or trivariate signals [START_REF] Tanaka | Complex empirical mode decomposition[END_REF][START_REF] Rilling | Bi variate empirical mode decomposition[END_REF][START_REF] Regman | Empirical Mode Decom position for Trivariate Signals[END_REF] and even recently to multivariate signals [START_REF] Regman | Multivariate empirical mode decomposition[END_REF]. All these latter methods, moti vated by the need to jointly analyze multi-channel real sig nals, are either restricted to a specifical signal dimension or are based on multiple projections of the original signal that complexify the original EMD process. The aim of this communication is to propose a new efficient EMD algo rithm, named 3A-EMD (Active Angle Averaging EMD), that preserves the original geometrical simplicity of the classical EMD algorithm and extends its performance to multivariate signals defined from lR to lR D with D > 2.

The classical EMD algorithm is extended to multivariate signals by redefining the mean envelope operator M us ing new characteristic points. A generalized version of the EMD is thus obtained working straightforwardly on monovariate and multivariate signal without any change.

A NEW EMD APPROACH: 3A-EMD

The proposed 3A-EMD method aims at providing a sim ple algorithm working for multivariate signals without any modification. To do so, the signal mean trend (computed by the M operator) is redefined as the signal going through the barycenter of its different elementary oscillations. This is performed through three main steps: i) identification of all the elementary oscillations, ii) computation of the barycenters of the identified elementary oscillations, and iii) interpolation between all barycenters to obtain the sig nal mean trend.

Elementary oscillation An elementary oscillation of a given signal with values in lR D (D ;::: : 1) can be viewed as a subpart of the signal defined between two extrema of oscillation. However, computing the extrema if the output space lR D is multidimensional (D > 1) is not obvious. In this section we propose a new and general way to compute such extrema. Let the signal {s(t)} with values in lR D be in class e 1 (or a "regular enough" signal in practical situa tions). The tangent vector to {s(t)} at t', Ts(t'), is defined by:

Ts(t') = [1, �: (t')W E ]R D + l (1)
and let's define the quantity, as (t'), as follows: as(t') = lim ( Ts(t' -h), Ts(t' + h))

(2 ) h--+O which is the inner product of ]R D + 1 (denoted by ( . , . )) between the tangents to {s (t)} just before and after t'. In fact, as(t) is maximum when both vectors, Ts(t -h) and Ts (t + h), are linearly dependent, due to Cauchy-Schwarz inequality, and null when they are orthogonal. This quan tity is an indicator of the local angle and thus, of the local oscillation of {s(t)} at t'. Due to the continuity of the inner product, we have: and as {s(t)} is in e 1 , we get:

where II . II abusively represents the norm of both ]R D + 1 and ]R D . Finally, an oscillation extremum of the signal { s (t) } is defined as a local minimum of the function { as (t') } and the following property holds:

Property: the signal {s(t)} has an oscillation extremum in to if and only if the function {II �� (t') 11 2 } is locally mini mum in to.

This definition of the extremum does not depend on the output space dimension D, which makes it suitable in mul tivariate EMD contexts. For D = 1, the reader could check that the computed extrema using the previous property in clude the signal scalar extrema (of Huang [START_REF] Huang | The em-pirical mode decomposition and the hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF]) but also the saddle points.

Oscillation barycenter Let PI = [h, s(hW and P 2 = [t 2 ' S(t 2 W be two consecutive oscillation extrema. A com plete oscillation is defined between P1 and P 2 and the as sociated oscillation barycenter, Mp1 --+ P2, is obtained by:

(5)

The mean trend M( {s(t)}) is then computed by interpo lating between the barycenters. The M operator appears with this definition as a non-stationary signal-based aver aging operator with a field of application not exclusively restricted to monovariate signals. Nevertheless, interpo lating directly through all the oscillation barycenters to obtain M ( { s (t) } ), appears in practice to be very sensi tive to artefacts introduced by the spline interpolation process and the border effects. Such effects appear very lesser in the classical Huang's approach as if they were "balanced" by the averaging of the min and max envelopes. Similarly, a robust implementation of the mean envelope computa tion is obtained for 3A-EMD by averaging two envelopes: a first envelope interpolating the even indexed barycenters (including signal borders) and a second envelope interpo lating the odd indexed barycenters (also including signal borders).

Sifting process stop criterion The sifting process termi nation is handled by means of a classical modified Cauchy like criterion [START_REF] Rilling | On empiri cal mode decomposition and its algorithms[END_REF]. The sifting process is thus stopped when ever Ildn,k + 1(ti) -dn,k(ti)ll/lldn,k(t i )11 becomes lesser than a certain threshold E (le -2 in practice) for a ratio of the whole samples of signal equal, at least, to a (90% in practice). (AI) The signal {s(t)} must be regular or sampled enough so that its numerical derivative has a relevant meaning. (A2) Any signal including one IMF with a piecewise con stant derivative norm (and especially purely circular ro tating bivariate signals) would not be successfully decom posed by 3A-EMD. (A3) The energy of the different IMF derivatives of the signal should ideally decrease or, at least, the energy of the n th IMF dn should not be "too much" bigger than the en ergy of the following IMFs dn + 1. In practice, noisy parts of a signal are often associated to highest frequencies and such condition must not be too restrictive.

Scope of application

Note on computational complexity The computational complexity of 3A-EMD algorithm can be precisely evalu ated and compared with some of the reference algorithms of the literature. For any given EMD algorithm, let's call N the number of extracted IMF, Kn the number of sifting iterations performed to extract the n th IMF, and dn,k the signal considered at the k th iterations of the sifting pro cess. Let's also call M H k ua n g the number of extrema de- n, tected in dn,k by the monovariate Huang's procedure [START_REF] Huang | The em-pirical mode decomposition and the hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF], M R ' k ' lli n g the number of extrema detected in the p th pro- n, , p jection of dn,k by the second geometrical bivariate Rilling's procedure [START_REF] Rilling | Bi variate empirical mode decomposition[END_REF] with P projection planes, and M�1 the num ber of barycenters detected in dn,k by 3A-EMD method.

The computational complexity of the sifting step neces sary to obtain dn,k + 1 for each of these three methods can be therefore computed for a signal from ]R to ]R D with S samples, by the formula given in the table 1. The global complexity for each of those methods is then directly ob tained by cmethod = ""N_ "" K :: cmethod.

L.. n-1 L.. k-1 n,k -1

RESULTS

The objective of this section is twofold: to compare the performance of3A-EMD algorithm with classical approa ches and to illustrate its extensibility to multivariate sig nals (D > 2). The trivariate signal S 3 1 is given by: Criterion 3 A last criterion evaluates the computational complexity with respect to the formula listed in Table 1.

x(t) = 3.5sin(5007l't) IMF 1 y(t) = 2cos(5007l't) z(t) = e O.
Comparative study This section compares 3A-EMD per formances with those obtained by Huang's algorithm [START_REF] Huang | The em-pirical mode decomposition and the hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF] in monovariate case and Rilling's method [START_REF] Rilling | Bi variate empirical mode decomposition[END_REF] in bivariate context. More precisely, the second geometric approach proposed in [START_REF] Rilling | Bi variate empirical mode decomposition[END_REF] is used as reference for bivariate decompo sitions with 8 projection planes. In both cases, the sifting process stopping criterion and the border effect manage ment are identical for all the compared methods.

(f � _ :1 �J � -1 0 1 2 -1 0 2 -1 0 2
Fig.

CONCLUSION AND PERSPECTIVES

A new algorithm, called 3A-EMD, based on a novel ge ometric definition of the mean envelope operator is pro posed. The obtained results show that, under certain as sumptions on the signals, this alternative definition en ables to decompose both monovariate and multivariate sig nals without any modification in the core of the algorithm. This last point is the one main difference with regard to the existing approaches of the literature. The comparative study also suggests that 3A-EMD seems to offer compet itive core and border performances as well as some in teresting performances in terms of computational com plexity. A restricted scope of application mainly due to the use of the first derivative in the algorithm remains the most important limitation of our approach. However, as pointed before, those restrictions should not be too much restrictive in practical fields. More simulated and real data decompositions should be performed in future works to verify the 3A-EMD interest in practical context.

  Based on the detection of local min ima of {II �� (t) 11 2 }, 3A-EMD raises some assumptions on the candidate signals:

  Signals selectionThe signals used to evaluate the behav ior of 3A-EMD algorithm are two monovariate (S11, S1 2

IDFig. 1 .

 1 Fig. 1. Comparative study of 3A-EMO versus Huang's (10) and Rilling's (20) reference (Ref) methods.

[ 1 ; 2 ]

 12 ofT= [-1; 2]. Indeed, border effects are often crit ical in practice and their minimization facilitates the ex ploitation of the IMFs.

2 .

 2 Fig. 2. 3A-EMD trivariate decomposition of the signal 831. Dashed dark line: real IMF / Continuous red line: estimated IMF.

Table 1 .

 1 Computational complexity of some reference algorithms and of 3A-EMO method.

	Sl1(t) = e O.23 ( 1 + t) (0.5sin(1007l't) + sin(507l't)) S1 2 (t) = 2sin(2507l't) + e O.23 ( 1 + t) cos( 7l'(2.58t +21.95) 2 ) + 3cos(7l'(1.7t + 7.3) 2 ) + 3t S 2 1 (t) = e O.23t -0.46 [0.3cos( 7l'(2.58t + 14.2) 2 )

two bivariate (S 2 1, S 22 ) and one trivariate (S 3 1) signals, defined on the time interval T = [-1; 2] and sampled at

is = 10kHz: +0.9isin( 7l'(2.58t + 14.2) 2 )] +[0.4sin(207l't) + 0.7icos(247l't)] S 22 (t) = [0.3cos(4007l't) + 0.8isin(5007l't)] +[4cos(1007l't) + 7isin(1007l't)] +e O.23t -O.46 [5cos( 7l'(2.24t + 4.47) 2 ) +2isin(7l'(2.24t + 4.47) 2 )]

  The fig ure I presents criteria I and 2 on the left, and criterion 3 on the right for both monovariate and bivariate signals. Core performances of 3A-EMD are globally satisfactory and competitive compared with those of the reference meth

			Axis #1				Axis #2				Axis #3
	0.48	0.49	0.5	0.51	0.52 0.48	0.49	0.5	0.51	0.52 0.48	0.49	0.5	0.51	0.52
	� 0.48	0.49	0.5	0.51	0.52 0.48	0.49	0.5	0.51	0.52 0.48	0.49	0.5	0.51	0.52
	i 0.48	0.49	0.5	0.51	0.52 0.48	0.49	0.5	0.51	0.52 0.48	0.49	0.5	0.51	0.52
							ods. Regarding border performances (eb), 3A-EMD gen
							erally offers more efficient results and seems to have a nice
							behaviour on the signal border. Finally, the computational
							complexity analysis suggests that, even if no significant
							differences are visible on monovariate signals, 3A-EMD
							generally requires less sifting iterations and less computa
							tional operations than the Rilling's method.

Trivariate decomposition

The figure

2

presents the ex act (dashed dark line) and estimated (red line) IMFs of � :�rvvvwINWJ -:WJ\NVWWj �:MJiMM!VW -:VVV\J1 lZ\Z\I\I: f\JV\J _ :K:Z::SJ _ :CZ=SJ -:CZ:S:J