Laura Bear

Rémi Dubois

Nejib Zemzemi

Optimization of Organ Conductivity for the Forward Problem of Electrocardiography

This study presents an empirical approach to optimize conductivities within a torso model given simultaneous epicardial and body surface potential recordings.The conductivities of the lungs, skeletal muscle and torso cavity were estimated within a forward model by minimizing the relative error betweencomputed and reference torso potentials using a standard gradientbased approach. The sensitivity of this approach was evaluated over different levels of geometric error and signal noise, and the gradient of the cost function was determinedusing both 1) finite differences and 2) an adjoint method.

All conductivities were accurately estimated (<10% difference in value) withup to 0.20mV signal noiseand all levels of electrode localization error (up to 2.56 cm) using a finite difference approach. While the adjoint approachwas more computationally efficient, a finite difference approach was more stable across different signals and more robust to noise.

Introduction

Theforward problem of electrocardiography defines the relationship between epicardial and body surface potentials. This relationship is fundamental to the inverse problem (non-invasive imaging of cardiac electrical activity from body surface measurements). As the inverse problem is ill-posed, it is important that the forward model is identified accurately.

In previous studies comparing directly recorded and forward computed body surface potentials in-vivo, it was shown that forward models incorporating inhomogeneous structures were more accurate than homogeneous models [START_REF] Bear | Forward Problem of Electrocardiography[END_REF][START_REF] Stanley | An assessment of variable thickness and fiber orientation of the skeletal muscle layer on electrocardiographic calculations[END_REF], improving both the magnitude and pattern of body surface potentials. Despite this, there still remained a difference between the forward and recorded body surface potentials.Although including further structural details into these forward models may improve them, such as by incorporating muscle fibre directions, it is often computationally and technically infeasible to do so.

The conductivity values assigned to different organs inmost forward models arebased on three studies [START_REF] Geddes | The specific resistance of biological material -a compendium of data for the biomedical engineer and physiologist[END_REF][START_REF] Rush | Resistivity of body tissues at low frequencies[END_REF][START_REF] Foster | Dielectric properties of tissue and biological materials: A critical review[END_REF], with the values varying substantially between them.Previously, it has beenshown that error in the conductivity values assigned can significantly affect body surface potentials, and in different ways depending on the organ [START_REF] Geneser | Application of stochastic finite element methods to study the sensitivity of ECG forward modeling to organ conductivity[END_REF].

Theoretically, the optimal conductivity for each region in a forward model canbe found given a data set with simultaneous measures of epicardial and body surface potentials. That is, by defining one organ as having a constant conductivity, its ratio to other organ conductivities can be optimized to best match predicted and observed torso potentials. The basis of this approach is empirical rather than physiological. Nonetheless, identification of an optimized transfer matrix would provide a better forward model than botha uniform isotropic modeland a more complex physiologically based one. This study examines such a method for optimizing conductivities using previously obtained invivo experimental data [START_REF] Bear | Forward Problem of Electrocardiography[END_REF].

Methods

Experimental Data

In-vivo, experimental datacame from an anesthetized, closed-chest pig [START_REF] Bear | Forward Problem of Electrocardiography[END_REF]. Epicardial and torso potentials were recorded simultaneously using an elastic "sock" (239 unipolar electrodes) and flexible strips attached to the body surface (184 electrodes). Upon completion, the heart was arrested and MR images acquired. The heart was excised and perfusion-fixed. Epicardial electrode locations were captured with a multi-axis digitizing arm. MRI contrast markers placed on the "sock" and body surface strips were localized in the MR images and used to register electrode locations.

Geometric and Computational Model

In order to define a test data set, a volume conductor model was developed using post-mortem MRI. An anatomically realistic generic model of a pig was customized to data clouds of the epicardial, lungs, fat and skeletal muscle surfaces, using a non-linear fitting procedure [START_REF] Bradley | Geometric modeling of the human torso using cubic hermite elements[END_REF]. Linear-triangular surface meshes were then coupled together to create a finite element volume conductor model (Fig 1). The torso domain is denoted byΩ, covering the volume between the epicardium and body surface. Ω = ∪ 𝑖∈{𝑓,𝑚 ,𝑙,𝑐} Ω 𝑖 , whereΩ 𝑓 , Ω 𝑚 , Ω 𝑙 , Ω 𝑐 are the fat, skeletalmuscle, lung and cavity volumes respectively. The conductivity in each region was definedasσ 𝑓 = 0.04, σ 𝑚 = 0.40, σ 𝑙 = 0.05, σ 𝑐 = 0.22 mSmm -1 , consistent with the literature [START_REF] Geddes | The specific resistance of biological material -a compendium of data for the biomedical engineer and physiologist[END_REF][START_REF] Rush | Resistivity of body tissues at low frequencies[END_REF][START_REF] Foster | Dielectric properties of tissue and biological materials: A critical review[END_REF].The boundary of the torso domain is defined by 𝜕Ω = Σ ∪ 𝛤 𝑒𝑥𝑡 , where Σ represents the epicardial surface and 𝛤 𝑒𝑥𝑡 the external boundary of the body surface. We denote the torso potential in Ω by𝑢 𝑇 .By solving Laplace's equation over this volume using the finite element method, we define the linear relationship: 𝑑 = 𝐴ℎ. At each time step, reference potentials, d, at 180 electrode locations on 𝛤 𝑒𝑥𝑡 were computed from the measured potentials, h on Σ.

Optimization Procedure

The test data provides the potential values on both Σand 𝛤 𝑒𝑥𝑡 . Moreover, we know that the current flux over the body surface is equal to zero. In order to estimate the values of each organ conductivity, we construct the following quantity of interest.

I σ 𝑓 , σ 𝑚 , σ 𝑙 , σ 𝑐 = 1 2 𝑢 𝑇 -𝑑 𝐿 2 (𝛤 𝑒𝑥𝑡 ) 2 ,
𝑤𝑖𝑡ℎ 𝑢 𝑇 𝑠𝑜𝑙𝑢𝑡𝑖𝑜𝑛 𝑜𝑓: 𝑑𝑖𝑣 σ 𝑇 ∇𝑢 𝑇 = 0, 𝑖𝑛 Ω, 𝑢 𝑇 = ℎ, 𝑜𝑛 Σ, σ∇𝑢 𝑇 • 𝒏 𝑇 = 0, 𝑜𝑛 𝛤 𝑒𝑥𝑡 By minimizing the cost function,I σ 𝑓 , σ 𝑚 , σ 𝑙 , σ 𝑐 , we find the optimal values of the conductivity for all four organs.Here, the fmincon function of Matlab 2013b was used. By default, this function approximatesthe derivative of I over σ 𝑖 for 𝑖 ∈ {𝑓, 𝑚, 𝑙, 𝑐}using a finite difference approach (Grad FD ).

In the second half of this paper, we compare Grad FD to an alternative method, Grad adj , where the derivative of I over σ 𝑖 is computed directly giving

Method of Analysis

The sensitivity of the optimization procedure was tested by varying levels of signal noise (SN) on d, and torso electrode localization error (LE). SN wascreated using a random number generator with a standard deviation (SD) from0.002 to 0.512 mV. The direction of LE for each electrode was defined by picking a random point on the surface of a unit sphere, with the distance defined with a mean from 0.02 to 2.56 cm across all electrodes.Conductivities were optimized for each level of error, and for six different time points spanning ventricular depolarization. In addition, conductivities were optimized combining a SN of 0.05 mV and a LE 9 mm. These values were estimated to be the level of error in the original data set.

The initial conductivitiesused to start optimization were defined using a Monte Carlo simulation. That is, simulations were run using four values ofeach conductivity (σ f kept constant at0.04 mSmm -1 ) ranging from ±50% of the max/min values as compiled from the literature. The conductivitiesresulting in the lowest cost function were used as the initial conditions.

The accuracy of optimized conductivitieswas evaluated using the relative error (RE). In addition, the speed of computation was evaluated using the number of iterations for computation. To determine differences, a paired t-test was used for normally distributed data, and a two-sided Wilcoxon signed rank test for non-normal. Statistical significance was accepted for p< 0.05.

Results

Sensitivity Analysis

Given the correct fat conductivity, the Monte Carlo simulations yielded initial conductivities of σ l = 0.06, σ m = 0.35, σ c = 0.13 mSmm -1 , for all levels of signals noise, vest error, and for all time steps. In the following analysis, Grad FD was used to define the gradient during optimization. As expected, the accuracy of optimization decreasedwith increasing levels of SNand LE (Fig. 2). All three conductivities were accurately estimated (RE <10%) for up to 0.20 mVSN and for all levels of LE. The accuracy was also dependent on the signal used as demonstrated by the large SD of the results for each error level. Conductivities optimized combining SN (0.05 mV) and LE (9 mm) were accurately estimated for all signals used (Table 1).

Adjoint Method for Gradient

The accuracy and speed of optimization using an adjointapparoachto compute the gradient (Grad adj ) were compared to Grad FD .Given results are more sensitive to it, this comparison was performed over different levels of SN. Fig. 3apresents the mean RE in conductivity against increasing SN for both methods.For all error levels except the greatest, the conductivities were more accurately estimated, and the SD smaller usingGrad FD .Overall, Grad adj resulted in a significantly larger average error in conductivity (p < 0.0043) than Grad FD .

Fig 3b shows the number of iterations necessary to converge to a solution for each method. From this, we can see Grad adj is substantially more computationally efficient than Grad FD , requiring 45 to 64fewer iterations to converge to a solution for all levels of signal noise. Interestingly, the number of iterations remains relatively stable with increasing noise levels for both methods.

Discussion

In this study, we have presenteda method for optimizing the conductivity of different organs in a torso model, given simultaneous recordings of epicardial and torso potentials. Here, conductivities of 3 organs were successfully computed undertypical signal noise and geometric error levels. A sensitivity analysis revealed the conductivities could be computed to within 10 % of their true values when the standard deviation of signal noise was less than 0.20mV and the mean electrode localization error less than 2.56 cm.

The accuracy of the final conductivity values was dependent on the signal selected for optimization, demonstrated by an increase in SD with error. Further analysis showed that signals taken near the start and end of depolarization were typically less accurate than those near the middle of depolarization. Thus, care must be taken when selecting the signals used for optimization as the method is dependent not only on the noise level but the signal-to-noise ratio itself.

Generally, optimization is more robust, and faster when you include directly computed gradients.Here, while Grad adj was more computationally efficient than the standard Grad FD approach provided by Matlab, it was less robust to signal noise. Despite having not converged to the exact conductivity using Grad adj , the final cost function values were not largely different from those using Grad FD , where conductivities accurately converged. Thus, we could assume the gradient calculated using Grad adj typically resulted in small step sizes, and the optimization process was exited due to a step size or a change in cost function below threshold. Changes to these thresholds would increase the computation time, with no guarantee of a substantially improved convergence.

Previous studies have shown structure-based inhomogeneous models have proven to be less robust than homogeneousmodels in inverse simulations [START_REF] Ramanathan | Electrocardiographic imaging: II. Effect of torso inhomogeneities on noninvasive reconstruction of epicardial potentials, electrograms, and isochrones[END_REF]. That is, although more accurate physiologically,inhomogeneous models are also more prone to error. However, here we have shown that an optimized transfer matrix can be identified, providing a better forward model than a uniform isotropic volume conductor, andpotentially more robust inverse solutions than complexphysiologically based inhomogeneous models. Thus these methods could be of clinically to improve current electrocardiographic techniques.

Whilst a full set of epicardial potentials is typically unavailable clinically, an electrical current could be applied to the heart in a known location during isoelectric periods. Given body surface potentials recorded concurrently and patient-specific MRI, the organ conductivities within a patient could theoretically be estimated. Further analysis in this direction is currently underway using the experimental data outlined.

Fig 1 .

 1 Fig 1. Cranial view of the torso model. Regions include the epicardium (red) lungs (light blue), cavity (yellow), skeletal muscle (orange) and subcutaneous fat (dark blue).

2 ΩΩ

 2 𝜕𝐼 𝜕σ 𝑖 σ 𝑓 , σ 𝑚 , σ 𝑙 , σ 𝑐 = 𝜕u 𝑇 𝜕 σ 𝑖 𝛤 𝑒𝑥𝑡 (u 𝑇 -𝑑) As the derivative of u 𝑇 cannot be computed directly, an adjointmethod was used. That is, 𝐻 1 (Ω)is denoted by the set of functions ∅: Ω → ℝ, such that ∅ is also denoted by the set of functions 𝜓 ∈ 𝐻 1 (Ω) such that 𝜓 Σ = 0.For every σ 𝑓 , σ 𝑚 , σ 𝑙 , σ 𝑐 , 𝑢, 𝜆 ∈ 𝐻 1 Ω × 𝐻 Σ 1 (Ω), the Lagrangian function is defined as follows 𝐿 σ 𝑓 , σ 𝑚 , σ 𝑙 , σ 𝑐 , 𝑢, 𝜆 = For every𝜆 ∈ 𝐻 Σ 1 (Ω)this gives I σ 𝑓 , σ 𝑚 , σ 𝑙 , σ 𝑐 = 𝐿 σ 𝑓 , σ 𝑚 , σ 𝑙 , σ 𝑐 , 𝑢, 𝜆 The gradient of I with respect to σ 𝑖 is then given by 𝜕I σ 𝑓 , σ 𝑚 , σ 𝑙 , σ 𝑐 𝜕σ 𝑖 = ∇𝑢 𝑇 ∇𝜆 Ω 𝑖 , 𝑤𝑖𝑡ℎ 𝜆 𝑠𝑜𝑙𝑢𝑡𝑖𝑜𝑛 𝑜𝑓: 𝑑𝑖𝑣 σ 𝑇 ∇𝜆 = 0, 𝑖𝑛 Ω, 𝜆 = 0, 𝑜𝑛 Σ, σ∇𝜆 • 𝒏 𝑇 = -(𝑢 𝑇 -𝑑), 𝑜𝑛 𝛤 𝑒𝑥𝑡 This method allows us to obtain the derivative of the objective function over the four conductivity parameters, by only solving two Laplace equations: The first is the state equation to obtain u 𝑇 and the second is the adjoint equation to obtain 𝜆. The derivative is then obtained by integrating the scalar product of the gradients of u 𝑇 and 𝜆over each of the four organs domains.

Fig 2 .

 2 Fig 2. RE in conductivity over (a) SN and (b) LE.

Fig 3 .

 3 Fig 3. Comparison of Grad FD and Grad adj over SN using (a) mean RE in conductivity and (b) number of iterations for computation.

Table 1Optimized

 1Optimized conductivity values and RE when combining SN of 0.05 mV and LE of 0.9 cm, presented as median [inter quartile range].

Acknowledgements

This work was supported in part by the French National Research Agency under Contract ANR-10-IAHU04-LIRYC. Laura Bear was supported by the FondationLeducqand grant ANR-10-IAHU-04.


Eddy Bajic

André Thomas

Ambient RFID Services Infrastructure & RFID Deployment in Wood industry

In this paper, we present academic research activities and industrial technology transfer conducted at CRAN-CNRS UMR 7039 Laboratory concerning on one hand high level services environment definition based on low cost RFID Tagged product for smart product concept achievement in an ambient service architecture based on standard IP technologies, and on a second hand RFID deployment in wood industry for logistic and product driven factory management.

INTRODUCTION

The presented works rely on research activities carried out at the Research Centre for Automatic Control (CRAN -CNRS 7039), within the Ambient Production System group SYMPA and the TRACILOG1 team, which aim to develop product-centric approaches in the fields of industrial production and management activities where on one hand intelligence is given to products by means of RFID technologies helped by ambient services systems architecture and on a second hand control and decision making model-based methods are conceived and evaluated. Industrial Application cases of RFID are conducted for implementing traceability and improving supply chain management in wood and fiber industries. One of the objectives of the CRAN is to provide a benchmarking tool enabling a comparative study of different control systems, like Holonic Manufacturing Systems (HMS) or Product-Driven Manufacturing systems. Several investigations in the fiber industry area, and especially in furniture industry, led to highlight the following context characteristics. First, production systems, including physical workshop and the planning system, has been developed to answer to a mass production context, that is to say to produce large quantities per lot for few part numbers. Secondly, clients ask for more and more customized products, at lowest cost and in shortening delays. To face these growing requirements, companies deploy several new tools, like autonomous teams management or simulation tools, in order to improve flexibility and reactivity of the production system. Today, new paradigms concerning intelligent products or ambient and communicating products appear. In a mass production context, this kind of concept must be evaluated according to global performance indicator such as throughput, for example. Consequently, in a HMS context the product has become a kind of gain inductor. In this paper, we will present first the communicating object concept applied to low cost RFID tagged product and then the enhancement of product/ process/operator relationship from simple data exchange to high level service providing and invocation at the product level. An Ambient services architecture is defined allowing intelligent product services capabilities supported by UPnP (Universal Plug and Play) technology to develop product services modeling and components design to provide intelligent product design. In a third part we present a RFID deployment application case in wood industry. The studied system is a panel board, ready-toassemble furniture production system. In this kind of company with high production volume is to maintain as low as possible the WIP level, as high as possible the load of the bottleneck and the service level. To achieve to these objectives a Product information based control project was released. In the part 4 we present the modeling methodology and the architecture implemented. Then we describe the simulation tool developed for the system control. Finally, we conclude and invoke open issues and new researches at CRAN lab.

RFID AND SERVICES FOR AMBIENT INTELLIGENCE SYSTEM

Increasing requirements of reactivity in the supply chain topic are observed among product, processes and clients across the product lifecycle. In this context, the interactions between processes, operators and product, beginning with manufacturing level until its use, require more information and automated intelligent exchanges between partners, in a sure and relatively quick way. Specifically, the requirements that emerge in the supply chain are associated to product customization, traceability and information management along its lifecycle, and, in general, all kind of services related to the product lifecycle using Internet technologies. In order to respond to these new requirements we propose an approach by considering a physical object as an active actor managing its evolution in the phases of its lifecycle, cooperating with multiple actors in the supply chain (supplier, producer, distributor, and consumer) to implement an intelligent service connectivity inside/beside the product figuring a trend toward an intelligent product behavior.

Intelligent Product and Smart Object :

Conceptually, the terms of intelligent product and smart product [START_REF] Kintzig | Objets Communicants[END_REF] are synonyms indicating the capability of an object to communicate with its active environment, and interact with its users or the other objects. A smart object is able to acquire, to receive and to distribute information in a near or distant environment, and is able to carry out actions on its own initiative. According to [START_REF] Duda | Ambient Networking. Smart Objects Conference 2003[END_REF], a smart object is aware of its environment and is able to perceive its surroundings through sensors, work with peers using short-range wireless communication technologies, and provide context-aware services to users in smart environments. The new paradigm of smart object provides an ability to embed new capabilities into object allowing extended access information up to complex services invocation, and interactions from virtually anywhere at any time, potentially transforming the way we live and work in a society of objects. The automatic identification technology plays an important role in the identification methods. Radio frequency identification (RFID) is a support tool that allows to automate object identification process and to improve operations managements [START_REF] Mcfarlane | The intelligent product in manufacturing control and management[END_REF]. A RFID system contains electronic labels or tags, tag readers that can read and write data, and a controller that controls the system. By attaching an electronic tag to a physical product, it can be automatically identified and located into the system vicinity. RFID represents the union between the physical product and the virtual product maximizing the use of the information and knowledge along product lifecycle [START_REF] Mcfarlane | The intelligent product in manufacturing control and management[END_REF].

Ambient services architecture :

In order to develop the interactions between intelligent products and other actors of the production system or the supply chain, the ambient services concept is analyzed and characterized therafter. In ambient services context, we consider the product as a service provider or a service requester, demanding transparent and associated control mechanisms of all the offered services in its environment [START_REF] Cea | Ambient Services for Smarts Objects in the Supply Chain Based on RFID and UPnP Technology[END_REF]. Ambient service is an abstract view of a system that provides information management capabilities, processing capabilities and event messages in an ambient network, representing a working domain. A product inserted in a network can appear or disappear in a domain of work [START_REF] Weiser | The computer for the 21st century[END_REF]. Ambient architecture based on Internet technologies have been developed integrating concepts such as ad-hoc network, mobility management, and service discoverproposed by y. With respect to the ubiquitous concept proposed earlier by [START_REF] Weiser | The computer for the 21st century[END_REF] responses. An event server publishes events to interested subscribers anytime the state of the service changes. Services are controlled by Control Points. A Control Point, noted UPnP(cp), is a logical entity that controls Device Points. After devices discovery, a control point retrieves the device description and get a list of associated services; retrieves service descriptions for interesting services; invokes actions to control the service; subscribes to the service's event source. Any time the state of the service changes, the event server will send an event to the control point. It is expected that devices incorporate control point functionality, noted UPnP(dv/cp), to enable true peer-to-peer networking interaction.

Principle of a Service-based relationship :

Once the object or product has entered the ambient architecture, the object is detected by an RFID ambient interface, and then the object is mounted in the ambient services architecture as service provider entity ubiquitously accessible by any control point system existing. The Physical Object is so transformed in a Device Object providing and performing transparent services. Fig. 1 depicts the primitives needed to achieve a Service based relationship between product and actors within a supply chain (Fig. 1). A device object figuring a Smart Object behavior can performs automatic eventing to the ambient partners so to performs reactivity from the object to the system.

Interaction Model :

The general mechanism used by a Smart Object acting as a control point in order to choose one resource between several, is represented in Fig. 2 Inspired from the Contract Net interaction protocol, standardized by FIPA, this interaction begins with the service request (the use of a resource, transport, etc.) sent by the service requester (the active Smart Object) to all available service providers (M) in the ambient UPnP IP network. The service providers respond by asking for the desired characteristics of the requested service (storage conditions, transport conditions, etc.). The providers, which offer the desired service and can answer to the request (P) will inform the requester of their availability, and propose the execution of the service (meeting, execution time, etc…). By this way, the service requester can choose amongst the providers the most adequate one for its needs. After the service execution is done, the service provider will send a response message with the execution status.

Methodological proposition :

A methodological offer is proposed so as to create an effective integration between a physical product and the interactive world of the information technologies within industrial processes and logistics [START_REF] Vollman | Manufacturing Planning and Control Systems[END_REF]. To present the conceptual solution three generic cases are considered:  Case One -UPnP assisted Passive Object: Case one (Fig. 3) represents an object carrying an electronic tag, which is managed by an UPnP(dv). Process begins with the automatic identification of the tagged object by an UPnP(dv). Then services associated to the identified product type are mounted in the device memory by means of XML files uploaded from local memory or remote database. A XML file summarizes the information about services, including actions and state variables. The services now available in the UPnP(dv) represent a virtual image of the product parameterized by the information stored in the tag or in a remote information system [6]. At that time, all the services associated to the product are known and can be remotely called by all UPnP(pc) in the ambient environment [START_REF] Jeronimo | UPnP Design by Example[END_REF]. Conceptually, the merger between the tagged object and an UPnP(dv) plugged in an UPnP network forms an entity called UPnP assisted passive object, representing a low-cost solution with a versatile technological integration. Static view in the service modeling considers that each UPnP service has a service type Uniform Resource Identifier (URI) that uniquely identifies the service. Every service also has a serviceID URI that uniquely identifies the service among all device's services.

Additionally, every service maintains three URLs that provide the information necessary for control points to communicate with services : the ControlURL is where control points post request to control this service; the EventSubURL is where control points post request to subscribe to events; the DescriptionURL tells controls point the location from which they can retrieve the service description document. Each service has zero or more state variable, having a name, a type, and a current value. A state variable also has a set of allowed values. Any of the state variables can trigger events on state changes as determined by the implementer of the service. Every input argument to an action is associated with one of the service's state variable.

Device and Service Description Documents :

Fig. 5 shows a device description document in XML that reflects the logical structure of a device, containing static information about the device. Device description contains a list of all of the services that the device provide. The service list provides the service's type, identifier, URL for retrieving the service description, and a URL for eventing services.

<device> <serviceList> <service> <serviceType>serviceType:v</serviceType> <serviceId>serviceID</serviceId> <SCPDURL>URL to service description</SCPDURL> <controlURL>URL for control</controlURL> <eventSubURL>URL for eventing</eventSubURL> </service> declarations for other services here </serviceList> </device> XML service description file document is summarized in Fig. 6, containing detailed information about a service, including the action list supported by the service. Service actions typically (Fig. 6) have the following forms: query specific state variable, set a device state variable, or tell to the device to perform some kind of task. Service state table is a list of the services state variables. State variables are used to model the state of the service at run time. A UPnP(cp) starts the subscription process by first retrieving the device description document form the device. The device description document contains an event subscription URL (the <eventSubURL>) and a service identifier (the <serviceId>) for each service provide by the device. Then, to subscribe to events from a particular service, the UPnP(cp) sends a subscription request to the service's event subscription URL. If the service accepts the subscription, the service responds with a unique identifier and a duration for the subscription. UML Sequence Diagram are used to represent how a UPnP(cp) subscribes and receives notifications of states variable changes in a service. Modeling of interaction between intelligent product and supply chain processes is thus supported by definition of a set of service classes associated to different stages of the product life cycle according to the capabilities of the tagged or smart object. In the following section, meta-services based on standards UPnP services are defined to manage an intelligent interaction in warehousing activity.

RFID DEPLOYMENT APPLICATION CASE IN WOOD INDUSTRY

The studied system is a panel board, readyto-assemble furniture production system. Furniture's manufacturing rely on a composition of three main functions: first, large panels are cut into pieces; secondly, pieces receive some transformations, like drilling or grooving, and then finished pieces are assembled into a package, which will be brought to the final customer through the distribution system. The cutting operations are particularly awkward to manage because of the important cadence variability. This fact could be explained by the following reasons. The panel price being one of the most significant costs constituting the unit cost, a special process optimizes the raw material use. So as to reduce cutting scraps, the produced quantity varies for each batch: with the aim of optimizing material use, the production of 400 desks could lead to produce 401 desktop and 404 left bases. The next session for this reference will take the remainder into account, which changes quantities for each component, so cutting plans will be really different. Moreover, the Cut-to-Size productivity curve is non-linear: the operating time to produce 800 pieces is not equal to the double of the time to produce 400 pieces. So the Cut-to-Size tool has variable operating times, which depend on many parameters. Indeed, the cutting plan influences the operating time due to the machine architecture, when the theoretic time, used for scheduling and for raw material ordering, is given by a formula based on constant throughput. The machine is seen as a "black box", but the resulting lead time of this operation varies a lot. Consequently robust global lead time is a problem. After this cutting process, pieces are manufactured to become finished pieces. The workshop is constituted of several machines, each one able to perform some operations. The variety of different pieces, the significant volume produced and the important number of lots being processed simultaneously generate complexity of physical flow, which makes their control difficult. Consequently robust product mix schedule is also a problem. In that context, the production management follows traditional MRP2 rules [START_REF] Vollman | Manufacturing Planning and Control Systems[END_REF]: a Sales and Operation Planning is generated according to commercial forecasts by family of products, and then is disaggregated into Master Production Schedules. The next step consists in computing the Material Requirements Plan and in balancing workloads on work centers by establishing a predictive finite capacity planning for each of them. This step enables to establish raw materials orders and finished furniture's delivery dates, which appears as constraints to be respected for the lower level. The recently acquired scheduling tool aims at optimizing resources usage rate. But its implementation reveals some imperfections:

• To be pertinent, a scheduling must be based on accurate technical data, • This scheduling tool does not reflect variability and dynamics of the actual process, • This scheduling tool does not take into account WIP level and its impact on productivity.

What particularly interests us is the launching on the following drilling and grooving machines. Each has particular features, and is able to accomplish several operations. Consequently, there are many possible operation sequences, and the pieces diversity makes physical flows and their control complex. Moreover, average lot size is about 400 units, which implies several important batches at the same time in the workshop and some quality scraps could be discovered in the decoupling point (St1 in Fig. 7) between these two first steps. Consequently the quantity insurance could be still a problem. The end of the process is the packaging line, where pieces have to be grouped by furniture reference. Of course, all the pieces of the bill of material of furniture have to be completed before packaging can begin. Three intermediate stocks enable to buffer physical flow between machines. Those stocks are "rolling stocks", equivalent to drawer stocks. The existence of these stocks leads the management more complex. Indeed, there are several FIFO stacks, and it is only possible to catch the first lot of each drawer. Consequently, we decided to focus on WIP and lead time control. Studying these two major objectives implies to take the system dynamic into account. In that sense, Product information based control appeared like a really adapted tool to improve physical flow management [START_REF] Karkkainen | Intelligent productsa step towards a more effective project delivery chain[END_REF].

MODELLING METHODOLOGY

Assessing the impact of synchronizing physical and informational flows needs to model them as distinct flows [START_REF] Pannequin | Centralized Versus Distributed Decisions, an Industrial Case[END_REF]. In that way, we represent them in two distinct models, that works simultaneously and have to be synchronized. Moreover, the interface standardization enables to exchange different control models with the same physical emulation model [13]. To represent the implementation of RFID technologies, with fixed readers, we will use the notion of synchronization points, which are points where the physical system emits events to update the information system. This event update could launch a decision process that will react by acting on the emulation model after impact measure on performance indicators. The emulation model does not integrate any decision rule, that is to say each point where decision have to be done is converted to a synchronization point: a message, called synchronization event, is sent to the control model, and the execution is paused until a resume order is received. On the synchronization event reception, the control model could compute and send some information or parameters to the emulation model, in order to realize the control of it. For example, in front of a machine, the control model send back the machining program, or at the exit point of a stock, the next step where the pieces batch have to been sent. According to the system theory, we distinguish three types of transformation: time, space and form, each corresponding to one kind of elements of the system : inventories, transport's and machines. So every type of element is represented by a basic element, that could be tuned function of the physical system it represents (batch size, transfer time…). A basic element is constituted of five components (Fig. 79 The synchronization procedure is quite simple: when a moving entity comes in a synchronization block, a message holding entity and localization identifier is sent to the control system, then the entity is held until an order is sent back by the control system. This procedure is shown on figure 10. The number of event type is finite : we could define one for the pre-synchronization and one for the post-synchronization for each type of resource.

Opposite to usual simulation models, moving entities only hold an unique identifier, which refers to corresponding data in the control model. As a consequence, the emulation model couldn't run without a control model, because entity related data are not stored like attributes of an entity but in a database integrated to the control model.

Control model :

Control model is a discrete event system, which reacts to external events coming from the emulation module, from the user or from any other source (representing lead time, quality, geo-localization lot information). To enable control decisions making, the control model needs a physical world representation and a decision making process. These components could be or not separated. To ensure the ability to exchange control model, we defined a synchronization procedure (Fig. 10) and a communication interface between emulation and control model: that is to say to standardize messages between emulation and control models. Such a message contains following information's:

Emulation Model Control Model

• A Synchronization point identifier, identify the resource, type, and synchronization type (pre or post). • An Entity identifier, pointer on the entity description in the database. This information enables the control model to update the information system, and then to compute the attempted decision. This decision would lead to a parameterization of the emulation model (route or process tuning), and then a "resume signal" would be send to the emulation model.

Simulation results:

The simulation model has been validated following a two steps process: first, functional properties have been validated, and in a second time, stochastic parameters have been identified using a statistic study. Functional properties integrate:

• Any lot follows the way corresponding to the real machining sequence (routings) • When more than one lot is available in front of a machine, the schedule is respected • Machine calendars are observed. Moreover, some alerts have been developed. These alerts generate a log file each time an error appears, which enable us to track the error. Machining process being more standard than the cutting one, the characterization of their stochastic parameters has been done analyzing previously recorded times for each machine. According to a 24 to 48 hours horizon, we were able to show the reliability of our system. Effectively, if the schedule is not or few changed, the forecast is relevant, but if the schedule is completely changed, due to commercial constraints for example, the forecast error grows up quickly (See comparative between real and simulated inventory on Fig. 8). Opposite to the Cut-to-size simulation tool, the reverse loop from the simulator to the informational system is manual: that is to say the simulator generates a report, which is analyzed and used by the scheduling operator. This way of establishing a schedule enables it to react in order to maintain the inventory level between efficiency boundaries. The daily use allowed improving the scheduling rules. The results are generated with a period of four hours: the global evolution is presented with curves, like on 5. This global evolution enables us to establish what are the interesting areas and then to obtain the detailed content of WIP to engage some corrective actions. On the other hand, the workshop simulator has been used like a training tool, what enables to considerably cut back the training time for a new scheduling operator, and to reduce problems while releasing new orders. The daily use of the workshop simulator let foresee an enlarged use of it, by validating some modifications of control rules before to implement them in the workshop if their impact are relevant. One of the concrete uses of the simulation was to decide how much operators have to be allocated to transport operations, or to choose the best location to install some new stocking structures. Moreover, this initially unforeseen use needs, ,the intervention of a specialist, when quickly educated users could do the daily use. The impact of major management rules could be studied too. For example, currently, the typical lot size is about four hundred pieces, but the growing need for reactivity incites to reduce this lot size. But the ability to test some new management rules drove to new needs, more and more complex and fastidious to apply for the model designer. This observation led us to develop a new modeling architecture, separating, on the one hand, the physical controlling system and, on the other hand, the control system, interfaced by an informational system according to the system theory.

CONCLUSION AND OPEN ISSUES

This paper has proposed a methodology offering a high level of integration of a physical product into an IP ambient network for developing wireless communication and intelligent interactions with the supply chain processes during product lifecycle. RFID technology permits to identify and seamlessly link the physical product to its extended services for maximizing the use of the information and knowledge along product lifecycle. Ambient services architecture provides high-level standards functionalities for discovering and learning about product's services in a domain of work. As a UPnP entity, the physical product assumes intelligent interactions in an IP ambient network. This RFID/UPnP based object represents an intelligent product that advertises and informs about its capabilities to the actors of the supply chain as also request services from others. The Smart Objects community proposed with the association of RFID and UPnP allows to create reactive and context sensitive environments, in which products are empowered through a digital environment based on services instead of elementary data. As though, the product's capacities are met in a given situation to form the Ambient Intelligence concept. Concerning the introduction of such concepts in industries, a case was presented. This industrial application highlights the interest to develop emulation and control simulation systems to achieve the SC global indicators.

Fig. 1 : 2 . 3 .

 123 Fig. 1: Principle of a Service based relationship between product and actors in a supply chain

Fig. 2 :

 2 Fig. 2: UML Sequence diagram describes the action invocation between an active communication object and another communicating object

Fig. 3 :

 3 Fig. 3: Ambient service based UPnP architecture integrating RFID smart object.

Fig. 4 :

 4 Fig. 4: Class Diagram: Device, Service, Actions and State Variable in UPnP.

Fig. 5 :

 5 Fig. 5 : Device Description Template.

Fig. 6 :

 6 Fig. 6 : Service Description Document.

Fig. 7 :

 7 Fig. 7 : Synthetic view of the workshop

Fig. 8 :

 8 Fig. 8 : Simulation model structure

Fig. 9 :

 9 Fig. 9 : Emulation model basic element

Fig. 7 :

 7 Fig. 7 : The synchronization procedure

Fig. 8 :

 8 Fig. 8: Real and simulated inventory evolution

  An exit point, where the entity is routed to the next step of the process.

	piece way and launching the next work
	on the machine,
	• Element	Pre-	Operational	Post-	Element
	Input	Synchro	Process	Synchro	Output
		Stocking	Transport	Processing
		Unit	Unit	Unit
	):
	• an entrance label,
	• a first synchronization point, called pre-
	synchronization, in charge of correctly
	setting up the corresponding resource,
	• a macro block representing the process
	performed by the resource, like in
	standard simulation model,
	• a second synchronization point, called
	post-synchronization, managing the

Traceability and management of the supply chain in fibers industries (wood, paper, textile)


Rabeb Chamekh

Abderrahmane Habbal
email: habbal@polytech.unice.fr

Moez Kallel
email: moez.kallel@ipeit.rnu.tn

Nejib Zemzemi
email: nejib.zemzemi@inria.fr

J A Dieudonné

A Nash-game approach to solve the Coupled problem of conductivity identification and data completion

Keywords: conductivity identification, data completion, Nash games

We consider the identification problem of the conductivity coefficient for an elliptic operator using an incomplete over-specified measurements on the surface (Cauchy data). Data completion problems are widely discussed in literature by several methods (see, e.g., for control and game oriented approaches [START_REF] Habbal | Neumann-Dirichlet Nash strategies for the solution of elliptic Cauchy problems[END_REF][START_REF] Aboulaich | Missing boundary data reconstruction via an approximate optimal control[END_REF], and references therein). The identification of conductivity and permittivity parameters has also been investigated in many studies (see, e.g., [START_REF] Lesnic | Identifiability of distributed parameters for high-order quasi-linear differential equations[END_REF][START_REF] Wang | Identification of conductivity and permittivity in a pulsed electric field model[END_REF]). In this work, our purpose is to extend the method introduced in [START_REF] Habbal | Neumann-Dirichlet Nash strategies for the solution of elliptic Cauchy problems[END_REF], based on a game theory approach, to develop a new algorithm for the simultaneous identification of conductivity coefficient and missing boundary data. We shall say that there are three players and we define three objective functions. Each player controls one variable and minimizes his own cost function in order to seek a Nash equilibrium which is expected to approximate the inverse problem solution. The first player solves the elliptic equation (div(k.∇(u)) = 0) with the Dirichlet part of the Cauchy data prescribed over the accessible boundary and a variable Neumann condition (which we call first player's strategy) prescribed over the inaccessible part of the boundary. The second player makes use correspondingly of the Neumann part of the Cauchy data, with a variable Dirichlet condition prescribed over the inaccessible part of the boundary. The first player then minimizes the gap related to the non used Neumann part of the Cauchy data, and so does the second player with a corresponding Dirichlet gap. The two players consider a response of the unknown conductivity of the third player. The third player controls the conductivity coefficient, and uses the over specified Dirichlet condition as well as the second's player Dirichlet condition strategy prescribed over the inaccessible part of the boundary. He minimizes then a Kohn-Vogelius type functional with respect to the conductivity parameter. This method is quite general and has wide applications ranging from bioelectrical field to mechanical engineering.

In this work, we are interested in solving the electrocardiography inverse problem which could be reduced to the data completion problem for the diffusion equation. The difficulty comes from the fact the conductivity values of the torso organs like lungs, bones, liver,...etc, are not known and could be patient dependent. Our goal is to construct a methodology allowing to solve both data completion and conductivity optimization problems at the same time. We consider the following elliptic problem :

   ∇.(k∇u) = 0 in Ω u = f on Γ c k∇u.ν = Φ on Γ c (1)
where Ω is a bounded open domain in R d (d = 2, 3) with a sufficiently smooth boundary ∂Ω composed of two connected disjoint components Γ c and Γ i . The functions f and Φ are the Cauchy data and ν is the unit outward normal vector on the boundary and k is a piecewise constant function representing the unknown conductivity.

For given η ∈ H -1 2 (Γ i ), τ ∈ H 1 2 (Γ i ) and k ∈ L ∞ (Ω) , let us define u 1 (η, k), u 2 (τ, k) and u 3 (τ, k)
as the unique solutions in H 1 (Ω) of the following elliptic boundary value problems :

(1)

   ∇.(k∇u 1 ) = 0 in Ω u 1 = f on Γ c k∇u 1 .ν = η on Γ i (2)    ∇.(k∇u 2 ) = 0 in Ω u 2 = τ on Γ i k∇u 2 .ν = Φ on Γ c (3)    ∇.(k∇u 3 ) = 0 in Ω u 3 = τ on Γ i u 3 = f on Γ c
The three costs are defined as follows:

J 1 (η, τ, k) = 1 2 k∇u 1 .ν -Φ 2 H -1 2 (Γc) + 1 2 u 1 -u 2 . 2 H 1 2 (Γi)
.

(2)

J 2 (η, τ, k) = 1 2 u 2 -f 2 L 2 (Γc) + 1 2 u 1 -u 2 2 H 1 2 (Γi)
.

(3)

J 3 (η, τ, k) = √ k∇(u 2 -u 3 ) 2 L 2 (Ω) . (4)
From the computational viewpoint, we used the Stackelberg game algorithm to compute the Nash equilibrium. We consider an annular domain with circular boundary components Γ i and Γ c , both centered at (0, 0) and with radii R i = 0.6 and R c = 1, respectively. However, a third cercle of radius R m = 0.8 is added between Γ i and Γ c . We consider a non-homogeneous conductivity k: equal to 1 ( is known) in the domain between R c and R m , and equal to 1 2 ( is unknown) in the domain delimited by R m and R i . The figure-1 presents the potentiel and the flux reconstructed over Γ i . The value of the reconstructed conductivity between R m and R i is equal to 0.48.

Figure 1 :

 1 Figure 1: Reconstructed Dirchlet (left) and Neumann (right) data over Γ i , where the reconstructed conductivity is 0.48


Relative Error (RE) in conductivity over different levels of SN Sensitivity Analysis Gradient Method

Optimization Procedure

Fig 2.

Database

• Epicardial potentials were recorded in-vivo using an elastic sock (239 electrodes), in an anesthetized, closed-chest pig [1]. • Post-mortem MRI was used to create a torso model (Fig 1), and epicardial electrode locations were captured with a multi-axis digitizing arm. • « Gold standard » torso potentials (d) were computed at 180 points from measured epicardial potentials using a finite-element defined forward model (Fig 1), and σ f = 0.04, σ m = 0.40, σ l = 0.05, σ c = 0.22 mSmm -1 .

Optimization of Organ Conductivity for the Forward Problem of Electrocardiography

Laura Bear 1,2,3 , Rémi Dubois 1,2,3 , Nejib Zemzemi

Methods

The forward problem of electrocardiography defines the relationship between epicardial and body surface potentials that is fundamental to non-invasive mapping. Forward models incorporating inhomogeneous structures are more accurate than homogeneous [1,2] though a difference between the forward and recorded potentials remains.

Theoretically, given a simultaneous measures of epicardial and body surface potentials, an optimal forward model could be found by optimizing the conductivities within the model. Identification of an optimized transfer matrix would provide a better forward model than both a uniform isotropic model and a more complex physiologically based model.

This study examines a method for optimizing the conductivities within a torso model using an in-vivo experimental data set [1]. The accuracy and speed of optimization using Grad adj were compared to using Grad FD over different levels of SN Overall, Grad adj resulted in a significantly larger average error in conductivity (p < 0.0043) than Grad FD . (Fig 5).

Grad adj is substantially more computationally efficient than GradFD, requiring 45 to 64 fewer iterations to converge to a solution for all levels of SN (Fig 6)

The final cost function values obtained with Grad adj were not substantially different from those using Grad FD (Fig 7).

The gradient calculated using Grad adj likely resulted in small step sizes, and optimization was exited due to a step size or a change in cost function below threshold. Given experimental data with simultaneous epicardial and torso recordings, conductivities within the torso can be accurately computed under typical signal noise and geometric error levels.

Though generally more robust when you include directly computed gradients, Grad adj was less robust to SN, though more computationally efficient than the standard Grad FD . This may be due to the optimization process was exiing due to a step size or a change in cost function below threshold with Grad adj .

The torso domain is denoted by Ω, covering the volume between the epicardium and the body surface

where Ω f,m,l,c are the fat, skeletal muscle, lung and cavity volumes, their respective conductivities denoted by σ f,m,l,c . The boundary of the torso domain is defined by

Ω = Σ U Γ ext ,
where Σ represets the epicardial surface and Γ ext the external bounday of the body surface. We denote the torso potential in Ω by u T . The test data provides the electrical potentials on both Σ and Γ ext . Moreover, we know that the current flux over the body surface is equal to zero. In order to estimate the values of each organ conductivity, we construct the following quantity of interest:

Using the fmincon function of Matlab 2013b, the cost function, I(σ f,m,l,c ), is minimized and the optimal conductivites for all four organs found. The derivative of I over σ i was computed initially using the default finite difference approache (Grad FD ). Comparison were then made with results computing the derivative directly, that is:

As the derivative of u T cannot be computed directly, an adjoint method was used. That is H 1 (Ω) is denoted by the set of functions , such that is denoted by the set of functions . For every , we define a Lagrangian function as Thus , I(σ f,m,l,c ) = L(σ f,m,l,c , u, λ) for every λ . The gradient of I with respect to σ i is given by This method allows us to obtain the derivative of the objective function over the four conductivity parameters, by only solving two Laplace equations: The first is the state equation to obtain u T and the second is the adjoint equation to obtain λ. The derivative is then obtained by integrating the scalar product of the gradients of u T and λ over each of the four organ domains.

I σ 𝑓 , σ 𝑚 , σ 𝑙 , σ 𝑐 = 1 2 𝑢 𝑇 -𝑑 𝐿 2 (𝛤 𝑒𝑥𝑡 2 , 𝑤𝑖𝑡ℎ 𝑢 𝑇 𝑠𝑜𝑙𝑢𝑡𝑖𝑜𝑛 𝑜𝑓: 𝑑𝑖𝑣 σ 𝑇 ∇𝑢 𝑇 = 0, 𝑖𝑛 Ω, 𝑢 𝑇 = ℎ, 𝑜𝑛 Σ, σ∇𝑢 𝑇 • 𝒏 𝑇 = 0, 𝑜𝑛 𝛤 𝑒𝑥𝑡 ∂ 𝜕𝐼 𝜕 σ 𝑖 σ 𝑓 , σ 𝑚 , σ 𝑙 , σ 𝑐 = 𝛤 𝑒𝑥𝑡 𝜕 u 𝑇 𝜕 σ 𝑖 (u 𝑇 -𝑑) 𝜕I σ 𝑓 , σ 𝑚 , σ 𝑙 , σ 𝑐 𝜕σ 𝑖 = Ω 𝑖 ∇𝑢 𝑇 ∇𝜆 , 𝑤𝑖𝑡ℎ 𝜆 𝑠𝑜𝑙𝑢𝑡𝑖𝑜𝑛 𝑜𝑓: 𝑑𝑖𝑣 σ 𝑇 ∇𝜆 = 0, 𝑖𝑛 Ω, 𝜆 = 0, 𝑜𝑛 Σ, σ∇𝜆 • 𝒏 𝑇 = -(𝑢 𝑇 -𝑑), 𝑜𝑛 𝛤 𝑒𝑥𝑡 ∅: Ω → ℝ Ω ∅ 2 < ∞ and Ω 𝛻∅ 2 < ∞ . 𝐻 Σ 1 (Ω) 𝜓 ∈ 𝐻 1 (Ω) such that 𝜓 Σ = 0 σ 𝑓 , σ 𝑚 , σ 𝑙 , σ 𝑐 , 𝑢, 𝜆 ∈ 𝐻 1 Ω × 𝐻 Σ 1 (Ω) 𝐿 σ 𝑓 , σ 𝑚 , σ 𝑙 , σ 𝑐 , 𝑢, 𝜆 = 1 2 𝑢 𝑇 -𝑑 𝐿 2 (𝛤 𝑒𝑥𝑡 ) 2 + Ω σ 𝑇 𝛻𝑢𝛻𝜆 ∈ 𝐻 Σ 1 (Ω)

Optimization and Analysis

• Initial conductivities were defined using a Monte Carlo simulation, using four values for each conductivity from ±50% the max and min conductivities from the literature.

Fig 3 .

 3 Fig 3. Relative Error (RE) in conductivity over different levels of LE

Fig 5 .

 5 Fig 5. Mean RE in conductivity over SN for Grad FD and Grad adj

•

  The sensitivity of the optimization procedure was tested by varying levels of signal noise (SN) on d, and torso electrode localization error (LE), for six different time points spanning ventricular depolarization.• To determine differences, a paired t-test was used for normally distributed data, and a two-sided Wilcoxon signed rank test for non-normal. Statistical significance was accepted for p < 0.05.Given the correct fat conductivity, Monte Carlo simulations yielded initial conductivities of σ l = 0.06, σ m = 0.35, σ c = 0.13 mSmm -1 , for all levels of signals noise, vest error, and for all time steps.In the following analysis, Grad FD was used to define the gradient during optimization Signal Noise (SN)• Gaussian SN was added to d, and conductivities optimized. • SN was created using a random nulber generator wth a standard deviation (SD) from 2 to 512 uV. • All three conductivities were accurately estimated (RE < 10 %) for up to 0,20 mv (Fig2)Vest Electrode Localization Error (LE)• Gaussian error was added to each vest electrode and conductivities optimized. • The direction of LE was defined for each electrode by picking a random point on the surface of a unit sphere, with the distance defined with a mean from 0,02 to 2,56 • All conductivities were accurately estimated (RE < 10 %) for all levels of LE (Fig 2) Combined Error • Conductivities were optimized combining a SN of 0.05 mV and a LE 9 mm (typical signal noise and geometric error levels) • All conductivities were accurately estimated for all signals (Fig 4)

Fig 4 .

 4 Fig 4. Relative Error (RE) in conductivity when combining 0.05 mV SN and 0.9 cm LE


C G L Tiofack

S Coulibaly

M Taki

S De Bièvre

G Dujardin

Periodic modulations controlling Kuznetsov-Ma soliton formation in nonlinear Schrödinger equations

We analyze the exact Kuznetsov-Ma soliton solution of the one-dimensional nonlinear Schrödinger equation in the presence of periodic modulations satisfying an integrability condition. We show that, in contrast to the case without modulation, the Kuznetsov-Ma soliton develops multiple compression points whose number, shape and position are controlled both by the intensity of the modulation and by its frequency. In addition, when this modulation frequency is a rational multiple of the natural frequency of the Kutzetsov-Ma soliton, a scenario similar to a nonlinear resonance is obtained: in this case the spatial oscillations of the Kuznetsov-Ma soliton's intensity are periodic. When the ratio of the two frequencies is irrational, the soliton's intensity is a quasiperiodic function. A striking and important result of our analysis is the possibility to suppress any component of the output spectrum of the Kuznetsov-Ma soliton by a judicious choice of the amplitude and frequency of the modulation.

I.

INTRODUCTION

During the last decades, rogue waves (also known as freak waves), which originally designated abnormally high sea surface waves, have been shown to arise, and consequently gained increasing interest, in various physical systems [1]. In addition to being found in the deep ocean, rogue waves have thus been discovered in optical systems [START_REF] Kibler | [END_REF]3], capillary flow [4], Bose-Einstein condensates [5], plasmas [6] and the atmosphere [7]. Rogue waves, which have a peak amplitude generally more than twice the average wave height, appear from nowhere and disappear without a trace [8]. Although their fundamental origin has not been understood completely, there is consensus [9] that rogue waves can be intimately associated with certain types of breathers of the underlying evolution equations. Breathers have drawn a lot of attention in this context because of their interaction and energy exchange with a constant background [10][11][12]. Such interactions produce unique behaviour that is different from bright soliton dynamics.

In many circumstances, the nonlinear waves under study can be accurately described by solutions of integrable one-dimensional nonlinear Schrödinger equations (NLSE) with constant coefficients and their exact analytical breather solutions have been interpreted as rogue waves [13]. Breather solutions on a finite background that have been identified are the Kuznetsov-Ma (KM) soliton [14,15], which is localized in the transverse direction, the Akhmediev breather (AB) [16], localized in the evolution direction, and the Peregrine soluton [17], which is a limiting case of the KM soliton and the AB, localized in both the transverse and evolution directions. Extensive theoretical and experimental studies have been devoted to characterize and control the formation of these solutions [18,19]. This includes collisions of ABs [20], the optimum conditions for the experimental generation of the Peregrine soliton [START_REF] Kibler | [END_REF] and of the KM soliton [21].

Rogue waves and breathers in variable coefficient nonlinear partial differential equations have attracted considerable interest recently and have been shown to present a rich variety of features. For example, recent studies have reported breather evolution, amplification and compression, Talbot-like effects and composite rogue wave structures in the variable coefficient derivative NLSE [22]. Experimental study has furthermore shown that the properties of the AB can be affected by a varying dispersion [23]. The appearance of a novel multiple compression point structure in a periodically modulated NLSE, referred to as a Peregrine comb, has been reported recently [24]. The controllable behavior, including partial excitation, maintenance, and postponement of the Peregrine, the AB and KM solitons have also been discussed in certain variable coefficient NLSE [25,26] describing dispersion decreasing fibers. Additionally, the modulational instability of the KM soliton in optical fibers with periodic dispersion was investigated in [27].

In this paper, we shall study the shape of the KM soliton in a one-dimensional completely integrable NLSE with periodically varying coefficients. We will show that, during the propagation, the spatial dependence or profile of the KM soliton interferes with that of the dispersion modulation and can be controlled by its intensity and period. We show in particular that the KM soliton with periodic coefficients develops multiple compression points whose number and height depend on the strength of the modulation. Further, we take advantage of nonlinear dynamics tools such as phase space analysis and Fourier spectrum to show that, when the frequency of the modulation and the natural frequency of the KM soliton are commensurate, the oscillations of the KM soliton's intensity are periodic. In the case where the two frequecies are incommensurate, the oscillations are quasiperiodic.

The paper is organized as follows. In Sec. II we show that under a suitable integrability condition the periodic modulation of the coefficients of the NLSE can lead to a generalized AB, or KM solution with periodic coefficients. In Sec. III the characteristics of the KM soliton with periodic coefficients are investigated. Using Fourier spectrum and stroboscopic Poincaré sections, the different regimes of the oscillations of the KM soliton with periodic coefficients are studied in Sec. IV. Section V contains our concluding remarks.

II. INTEGRABLE MODEL AND GENERALIZED SOLUTION WITH FINITE BACKGROUND

We consider the NLSE in the form

i ∂q ∂z - D 2 (z) 2
∂ 2 q ∂t 2 + R(z)|q| 2 q = 0, (1)
where D 2 (z) and R(z) are both taken to be periodic functions of their argument. In optical fiber systems the above equation describes the dynamics of the envelope q(z, t) of the electric field [START_REF]Nonlinear Fiber Optics[END_REF][START_REF] Serkin | [END_REF] where z and t are the dimensionless propagation distance along the fiber (measured in units of nonlinear length) and dimensionless time, respectively. The function D 2 (z) represents the GVD coefficient, and R(z) is the nonlinear one. In its general form, Eq. ( 1) is not integrable. However, if the condition RD 2z -D 2 R z = 0 is satisfied, Eq. ( 1) becomes the integrable NLSE with variable coefficients [30,31].

Here, we consider the focusing case where D 2 R ≤ 0 and assume consequently D 2 (z) = -R(z). In what follows we will study a soliton solution with finite background of Eq. (1) when the periodic dispersion is of the form

D 2 (z) = -1 + d m cos(k m z), (2)
where d m is the amplitude of modulation and k m its spatial frequency. Note that we are not assuming that the perturbations to the dispersion parameter D 2 (z), and nonlinearity R(z) are small. In this case, such a solution of (1) can be written as

q(z, t) = 1 + 2(1 -2a) cosh(BZ) + iB sinh(BZ) √ 2a cos(Ωt) -cosh(BZ) exp(-iZ), (3)
B = 8a(1 -2a), Ω = 2 √ 1 -2a, Z(z) = z z D 2 (z)dz,
where z fixes a global phase and a is a real constant parameterizing this class of the solutions [21]. Without loss of generality, we set z = 0. When d m = 0, Eq. ( 3) is reduced to the standard form of the soliton with finite background [16,32]. For 0 < a < 1/2, the solutions describe the AB solutions which exhibit periodic modulations in t with period 2π/Ω and a single growth-decay cycle along z. In the limit a → (1/2) -, the solution (3) reduces to the rational Peregrine soliton solution [START_REF] Kibler | [END_REF]. For a > 1/2, the parameters B and Ω become imaginary and we set b = -iB = 8a(2a -1) and ω = -iΩ = 2 √ 2a -1.

(4) The physical nature of the solution (3) changes as circular trigonometric functions become hyperbolic and vice versa. In this case, Eq. ( 3) reduces to the KM soliton,

q Ma (z, t) = 1 + 2(1 -2a) cos(bZ) -ib sin(bZ) √ 2a cosh(ωt) -cos(bZ) exp(-iZ), (5) with
Z = -z + d m k m sin(k m z). (6)
In the case of constant coefficients, Eq. ( 5) shows localisation in time t but (quasi)-periodicity along z [14,21]. If b = p/q with p, q mutually prime integers, then the solution ( 5) is periodic with period 2πq = p × (2π/b). Otherwise, it is quasiperiodic. Note that the modulus of q Ma is always periodic with period 2π/b as seen in Figure 1(a). By taking into account the effect of the modulation, the solution ( 5) is still localized in time but the spatial (quasi)-periodicity is affected by the modulation. For the particular case where the frequency of the modulation k m is a rational multiple of the frequency b of the oscillations in |q Ma | 2 , the modulus of the KM solution remains periodic in space (see Fig. 1(b)), otherwise it loses its periodicity (see Fig. 1(c)). In the case with constant coefficients, the periodicity characteristics of the AB and KM solutions are different, switching time and space (compare Eq. (3) with 0 < a < 1/2 and hence real B and Ω with Eq. ( 5) with 1/2 < a and hence real b and ω). As a result, the spatial modulation of the parameters leads to completely different physical mechanisms, as we will see below. In what follows, we will study the KM solution (5) in detail.

III. CHARACTERISTICS OF KUZNETSOV-MA SOLITON IN PRESENCE OF PERIODIC COEFFICIENTS

We now consider the generalized KM soliton Eq. ( 5), with a > 1/2, and the periodic dispersion coefficient given by Eq. (2). At a fixed point z one finds that the KM soliton's intensity has a maximum at t = 0, which is given by

|q Ma | 2 max (z) = 1 + 4 √ 2a(2a -1) √ 2a -cos(bZ) , ( 7
)
where Z is a function of z defined in (6). As stated previously and observed in Fig. 2, when d m = 0 the modulus of the KM soliton is periodic in z. In the presence of max is quasiperiodic. From the above expression we can see that the KM soliton oscillates in space (z-variable) and the maximum intensity value is given by

P = 1 + 4 √ 2a(2a -1) √ 2a -1 , (8)
which corresponds to cos(bZ) = 1. Hence the periodic variation of the coefficients in the NLSE has no effect on the maximum intensity of the KM soliton, which is identical to the one obtained in the case of constant coefficients. On the other hand, the periodic modulation drastically impacts the shape of the KM soliton, that we will analyse through the position, periodic or not, and the local structure of its extrema. Extrema of the KM soliton are critical points of the function |q Ma | 2 max . Hence, following Eq. ( 7), their positions z satisfy at least one of the two relations: 9), but not (10). Those extrema are less flat. Moreover, local extrema satisfy only Eq. ( 10) (and not Eq (9)) are never flat.

sin(bZ) = 0, (9)
D 2 (z) = -1 + d m cos(k m z) = 0. ( 10
As d m increases from 0 to 1, condition (10) cannot be satisfied, so that the extrema of the KM soliton are determined by Eq. ( 9), which is a transcendental equa- tion having no explicit analytical solution. In this regime of the modulation strength (d m ∈ (0, 1)), the general aspect of the KM soliton remains essentially unaltered, with (quasi-)periodic repetition of its extrema.

For the critical value d m = 1, Fig. 3 shows a typical example where we see the impact of the frequency of the modulation k m on the shape of the KM soliton. In this case (d m = 1), when k m /b is rational, the two equations ( 9) and (10) are simultaneously satisfied at points z satisfying bz ∈ 2π N 2 N 1 Z and bz ∈ πZ.

In case N 1 is even, this corresponds exactly to the points Observe that, when k m /b is irrational, there is only one simultaneous solution of Eq. ( 9) and ( 10), hence only one flat maximum, centered at z = 0.

z = π b N 2 , ∈ Z,
Beyond the critical value d m = 1, each flat maximum splits in three sharper maxima which are symmetrically located around the central one (see Figs. 6(a),(b)). Note that, when d m > 1, the function |q Ma | 2 max also has local maxima and minima. As d m > 1 increases, the amplitude of local maxima (minima) increases (decreases) to reach the value of the absolute maxima (minima). In case k m /b is rational, this can be seen in Fig. 6. In case k m /b is irrational, this can be observed in Fig. 7 where the position of the arrows shows how the local maxima evolve to the absolute maxima, then widen at the top, before they split in two (Fig. 7(c)).

IV. DIFFERENT REGIMES OF THE OSCILLATIONS OF KUZNETSOV-MA SOLITON IN PRESENCE OF PERIODIC COEFFICIENTS

In order to shed a different light on the regimes of the spatial oscillations of the KM soliton as a function of the frequency of the modulation, we study the Fourier spectrum and stroboscopic Poincaré sections for the function |q Ma | 2 max defined in (7) in detail. As we will see, if More generally, when the ratio k m /b = N 1 /N 2 with N 1 and N 2 mutually prime integer numbers, the signal |q Ma | 2 max is still periodic (with period z * = 2πN 1 /k m ), but displays sharp pulses (see Fig. 8(b)), due to specific relations between its Fourier components, reminiscent of phase-locking [START_REF] Berge | L'ordre dans le chaos[END_REF][START_REF] Spatschek | Nonlinear Coherent Structures[END_REF]. The phase space diagram is consequently composed of a set of N 1 points. This can be seen in the middle panel of Fig. 8(b) for k m /b = 5/2 where we have numbered P j (j = 1, . . . , 5), the successive points in the Poincaré section. In this case, everything happens as if a nonlinear resonance occurred between the frequency k m of the modulation and the intrinsic frequency b of the KM soliton.

k m /b = N 1 /N 2 is
Finally, as already pointed out in Sec. III, when k m /b is an irrational number, the function |q Ma | 2 max is quasiperiodic. In this case, the frequency k m of the modulation and the natural frequency b of the KM soliton are incommensurate and as a result, the oscillations are more complex, but of course not chaotic. The obtained stroboscopic phase-space trajectory never returns to the same point. Instead, it traces out a closed loop (see Fig 8(c)). This can be explained as follows. Introducing the 2π-periodic function

G(θ 1 , θ 2 ) = 1 + 4 √ 2a(2a -1) √ 2a -cos(-θ 1 + b dm km sin(θ 2 ))
, in two angular variables θ 1 , θ 2 and setting ∇G = (∂ θ1 G, ∂ θ2 G), one notes that the phase-space trajectory z → (|q Ma | 2 max (z), (|q Ma | 2 max ) z (z)) is the image of the trajectory on the torus θ(z) = (θ 1 (z), θ 2 (z)), with θ 1 (z) = bz and θ 2 (z) = k m z, by the function (G, (b, k m ) • ∇G). Indeed, by Eq. ( 7), we have

|q Ma | 2 max (z) = G(θ(z)),
and

|q Ma | 2 max z (z) = θ (z) • ((∇G)(θ(z))) .
Recall that the points (P j ) j∈N correspond to values of z of the form

z j = j2π/k m + z 0 , so that θ 2 (z j ) = θ 2 (z 0 )[2π],
independently of j, and consequently the stroboscopic phase-space points fill the one-dimensional loop

θ 1 → (G(θ 1 , θ 2 (z 0 )), (b, k m ) • ∇G(θ 1 , θ 2 (z 0 ))).
Many subharmonics appear again in the Fourier spectrum, but they are not equally spaced as before. Instead, additional frequencies occur at all values of the form n 1 b + n 2 k m with n 1 and n 2 arbitrary integers. This is illustrated in Fig. 8 (c) for k m /b = √ 3/2. It is to be noted that the side frequencies that significantly appear in the spectrum (whether k m /b is rational or not) are determined by the strength d m of the modulation as we will illustrate below. As d m increases, the energy content that is initially concentrated, when d m = 0, in the carrier frequency b and its harmonics is distributed among an increasing number of sidebands. However the relation

+∞ -∞ |q M a (t, z)| 2 -1 dt = 2ω,
implies that the total energy is conserved, as expected. Similar results have been obtained in [35] where extreme wave phenomena in down-stream running modu-lated waves have been investigated.

To better explain the appearance of additional sideband frequencies as d m is increased, we consider the asymptotic case where a 1. In this limit, Eq. ( 7) can be approximated by:

|q Ma | 2 max (z) = 8a + 4 √ 2a cos(bZ). (11)
Now, the cosine function can readily be decomposed into its component frequencies as follows

cos [bz -β sin(k m z)] = J 0 (β) cos(bz) +J 1 (β) [cos(b -k m )z -cos(b + k m )z] +J 2 (β) [cos(b -2k m )z + cos(b + 2k m )z] +J 3 (β) [cos(b -3k m )z -cos(b + 3k m )z] . . . +J n (β) [cos(b -nk m )z + (-1) n cos(b + nk m )z] (12) + . . . ,
where b ≈ 4a, and where β = d m b/k m is the modulation index.

It can be seen from the above relation that the amplitudes of the carrier and sideband components in the spectrum of the KM soliton are determined by the nth order Bessel functions of the first kind J n (β). For small values of the modulation index, J 0 (β) 1 and J n (β) β n and the spectrum of the KM soliton is very close to the one obtained in absence of modulation. This can be seen in the spectrum of Fig. 9(a) corresponding to β = 0.1 where the peak amplitude of the carrier frequency k = b is clearly strongly dominant. Further increase of β results in successive side frequencies and the relative amplitudes of all the frequencies present in each specific case depend on the values of the Bessel functions involved. Some of the amplitudes, that of the carrier included, can be suppressed completely when the index β coincides with a root of the corresponding Bessel function. For example, d m = 0.6 and k m = b/4 correspond to β = 2.4. The amplitude of the carrier frequency which is controlled in this case by J 0 (2.4) and J 8 (2.4) vanishes in the spectrum as shown in Fig. 9(b). This is a striking and an important issue since it means that external modulation (d m and k m ) can be tuned to shift at will the spectral frequencies of the system.

V. CONCLUSION

We have carried out a theoretical investigation of the generalized Kuznetsov-Ma soliton of a completely integrable nonlinear Schrödinger equation with periodically varying coefficients. The impact of the periodic modulation on the shape of the generalized Kuznetsov-Ma soliton has been analyzed. The amplitude of the modulation has no effect on the maximum intensity of the Kuznetsov-Ma soliton, but it promotes the development of multiple compression points. On the other hand, the frequency of the modulation drastically impacts the shape of the solution. We have then investigated the nonlinear nature of the Kuznetsov-Ma soliton's intensity as a function of the frequency of the modulation. The regime of the oscillations were studied in the light of the dynamical systems tools such as Poincaré sections and Fourier spectrum. It was found that when the ratio between the natural frequency of the Kuznetsov-Ma soliton and the frequency of the modulation is rational, a nonlinear resonance occurs between the two frequencies and the Kuznetsov-Ma soliton's intensity is a periodic function.

In the case where the ratio of the above two frequencies are irrational, quasiperiodic oscillations were observed. A striking and important result from our investigation is the possibility to suppress either the carrier frequency or any subharmonic by a judicious choice of the amplitude and frequency of the modulation. One consequence is that external modulation can be tuned to experimentally shift at will the spectral frequencies of the system. This makes it possible to control the formation and the properties of rogue waves of the Kuznetsov-Ma type. The extension of these theoretical results to a large class of nonintegrable systems with modulated coefficients will be the subject of further investigation.

Belgium Science Policy Office under the grant IAPP7-35 and the French Agence Nationale de la Recherche project OptiRoc ANR-12-BS04-0011.

10 FIG. 1 .

 101 FIG. 1. (Color online) Spatiotemporal evolution of the generalized Kuznetsov-Ma soliton. (a) dm = 0, (b) dm = 2, km = b, (c) dm = 2, km = √ 3b/2. Here, a = 0.603.

2 FIG. 2 .

 22 FIG. 2. (Color online)The shape of the Kuznetsov-Ma soliton as a function of z in the case of constant coefficient dm = 0 with a = 0.603.

FIG. 3 .

 3 FIG. 3. (Color online) The shape of the generalized Kuznetsov-Ma soliton with periodic coefficients as a function of z with dm = 1 and a = 0.603. (a) km = b, (b) km = b/2.

FIG. 4 .

 4 FIG. 4. (Color online) The shape of the generalized Kuznetsov-Ma soliton as a function of z with a = 0.603 and dm = 1. (a) km/b = 1/3, (b) km/b = 2/3.

  FIG. 5. (Color online) The shape of the generalized Kuznetsov-Ma soliton with periodic coefficients as a function of z with km = √ 3b/2 and a = 0.603. (a) dm = 1, (b) dm = 1.3.

FIG. 6 .Fig. 8 (

 68 FIG. 6. (Color online) The shape of the generalized Kuznetsov-Ma soliton with periodic coefficients as a function of z with km = 2b and a = 0.603. (a) dm = 1, (b) dm = 1.5, (c) dm = 4.0

FIG. 7 .FIG. 8 .

 78 FIG. 7. (Color online) The shape of the generalized Kuznetsov-Ma soliton as a function of z. (a) dm = 1.5, (b) dm = 2.2, (c) dm = 2.3. The other parameters are km = √ 3b/2, a = 0.603.

FIG. 9 .

 9 FIG. 9. (Color online) The spectrum of the KM soliton with periodic modulation for a = 10.255. (a) dm = 0.1, km = b, (b) dm = 0.6, km = b 4 .

  and they are global maxima if N 2 is even or global minima if N 2 odd. In case N 1 is odd, this corresponds to the points

	z =	2π b	N 2 ,	∈ Z,

and all of them are global maxima. In case d m = 1, any local extremum is global.

ACKNOWLEDGEMENTS C. G. L. T. acknowledges the support of the "Laboratoire d'Excellence CEMPI: Centre Européen pour les Mathématiques, la Physique et leurs Interactions" (ANR-11-LABX-0007-01). This research was supported in part by the Interuniversity Attraction Poles program of the