HERRI Le Quang Duyen

Ouabbas Yamina

A method to calculate the composition of clathrate hydrate of CO 2 -CH 4 , at equilibrium and during crystallization

Keywords: methane hydrate, dissociation, sediment core, heat and mass transfer

In recent decades, scientists have made an amazing discovery: methane hydrates in the ocean floor. These solid compounds could represent an energy reserve. Today it is estimated that they contain equivalent carbon under the form of methane in a huge quantity, twice as much as all deposits of natural gas, oil and coal.

In the near future, we need evaluating the possibility to produce this new source of energy, particularly in replacement of oil and coal. The main question concerns the technology to be used because the methane hydrates are distributed in sediment, and they participate to their consolidation. In this paper we present a method to study the thermodynamic of CO2-CH4 clathrate hydrates, which is the prior step before testing and modeling the replacement of methane hydrate by CO 2 hydrate after injection of CO 2 gas.

Introduction

Gas hydrates are solid compounds. They are formed from the combinaition of gas (such as methane, ethane, propane ...) and water under high pressure and low temperature. Clathrate hydrate has been first discovered in 1778 by Joseph Priestley as a curiosity of laboratory. Nowadays, gas hydrates have the potential for numerous applications in the oil and gas industry and the energy sector, as for example through the use of clathrate hydrates as a means of gas storage, for the capture and sequestration of carbon dioxide, in air-conditioning systems in the form of hydrate slurries, in the water desalination and treatment, and the separation of different gases from flue gas streams to name but a few [START_REF] Eslamimanesh | Application of gas hydrate formation in separation processes: A review of experimental studies[END_REF]. However, despite of the potential applications of gas clathrate hydrates like the ones mentioned above, there are also negative aspects to be mentioned in the discussion of these solid solutions. The uncontrolled decomposition of naturally occurring methane hydrates for example has been discussed as being capable of potentially contributing to the green house gas effect [START_REF] Englezos | Chlathrate Hydrates[END_REF][START_REF] Leggett | Global Warming: The Greenpeace Report[END_REF], in particular if it is realised that the global warming potential (GWP) of methane within a period of 100 years is greater by a factor 25 than the GWP value of carbon dioxide [START_REF] Solomon | Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change[END_REF]. Moreover, gas clathrate hydrates have been identified as a source of problems in the oil and gas industry, e.g. when being formed in drilling applications [START_REF] Barker | Formation of hydrates during deepwater drilling operations[END_REF] or in gas pipelines due to their ability of causing pipeline blockages [START_REF] Eslamimanesh | Application of gas hydrate formation in separation processes: A review of experimental studies[END_REF].

Methane is a natural component in sediments, originated from thermal degradation of fossil reservoirs, or originated from bio-degradation of biological materials. Under pressure, in deep sea conditions, it forms methane hydrate reservoirs in many places of the world and in huge quantities. On a other hand, carbon dioxide is a molecule which presents a better affinity to clathrate structure. The concept of CH 4 -CO 2 exchange consists in injecting CO 2 in methane hydrate reservoirs in order to shift the equilibrium and form a CO 2 hydrate in the place of methane hydrate, and to recover methane

The objective of our work is to study thermodynamic of CO 2 -CH 4 methane hydrate, in presence of salts, and to understand the crystallization mechanism during the replacing of CH 4 by CO 2 . In this paper, we present a set of experimental results in order to validate a thermodynamic model which predicts the equilibrium conditions; Pressure, Temperature and composition of phases in presence of a mixture of CO 2 and CH 4 .

Advantages of using CO 2 gas for replacement of CH 4 hydrate in sediment

Firstly, the method consisting to replace the methane in the methane hydrates reservoirs by CO 2 will contribute to drop the emissions of greenhouse gases.

Secondly, the hydrate structure still remains in place in the sediment, only the gas content being shifted. So, it can be assumed that the mechanical stability of the reservoir is not modified.

[ref sfgp:2013120099]-2 Lastly, the heat of CO 2 hydrates formation (57,98 kJ / mol) is slightly higher than the heat of methane hydrates dissociation (54,49 kJ / mol) (Naval Goel 2006). Therefore, the transformation of methane hydrate into carbon dioxide hydrate is practically athermic:

CO 2 (g) + nH 2 O→CO 2 (H 2 O)n ΔHf = -57,98kJ/mol (1) CH 4 (H 2 O)n→CH 4 (g)+ nH 2 O ΔHf = 54,49kJ/mol

(2) Where n is hydration number (n about 6).

On Figure 1, we can determine the (Pressure-Temperature) area where carbon dioxide can be injected for dissociation of methane hydrate. It is the region A or B in which the methane hydrate is not longer stable, but the carbon dioxide still remains stable. [START_REF] Ohgaki | Methane exploitation by carbon dioxide from gas hydrates-phase equilibria for CO2-CH4 mixed hydrate system[END_REF] have demonstrated experimentally the possibility to produce the methane gas from the hydrate by injecting carbon dioxide. They found that the mole fraction of carbon dioxide in the hydrate phase becomes higher than in the gas phase, at the equilibrium. Furthermore, [START_REF] Hirohama | Conversion of CH4 hydrate to CO2 hydrate in liquid CO2[END_REF] showed that the pressure in reservoir after methane hydrate dissociation remains approximately the same as before dissociation. [START_REF] Seo | Hydrate phase equilibria of thecarbon dioxide, methane, and water system[END_REF] suggest that if the mole fraction of carbon dioxide is kept above 40% in the CO 2 -CH 4 gas mixture, the concentration of carbon dioxide in hydrate phase will be more than 90%.

Figure1. Equilibrium conditions of hydrate CO 2 and CH 4 (Naval Goel 2006)

Experimental approach

Materials of experiment

In this work we used two bottle gas of CO 2 and CH 4 (4, 5 in the Figure 2) are provided by Air liquide under high pressure cylinder. The water for experiment is distilled water was obtained from the MILLI-Q 185 PLUS system. In the liquid phase contains LiNO 3 as a tracer (10ppm weight fraction). In this experimental the quantities of materials are:

Experimental set-up

The experimental apparatus (Figure 2) is designed to measure the thermodynamic equilibrium points in presence of gas mixture and to determine the composition of all phases (gas, liquid and hydrate). The reactor consists of a 2.5 liter autoclave reactor in which the pressure can reach up to 100 bars. The reactor is equipped with a vertical stirrer with four blades. The stirring rate can vary between 0 and 600 [ref sfgp:2013120099]-3 rpm. The temperature is controlled by a double jacket in which is circulated a fluid at constant temperature from a cryostat HUBERT CC-250. A Pyrex cell is laid inside the reactor and filled with water containing Lithium as an anion tracer. The concentration of the anion is controlled in the range 0-10ppm by using a HPLC pump (JASCO). Temperature is monitored by 2 Pt100 and pressure is measured with an accuracy of 0.01 MPa in the range 0-10MPa. A ROLSI sampler is mounted on the reactor. It allows to sampling online the gas and to sending the sample into a gas chromatograph (GC Varian model 38002) equipped with a TCD detector and two columns PoraBOND Q and CP-Molsieve. The peak integration is possible with software provided by Varian Galaxie. A sampling system can retain the liquid phase through a valve. It is analyzed off-line by ion chromatography and refractometry. The data acquisition (Pressure, Temperature) is saved on the personal computer.

Figure2: Experimental device

Procedure of our experiences

Initially, the reactor is evacuated by a vacuum pump. The methane (or carbon dioxide) gas is injected and purged three times to ensure that air is totally evacuated.

The experiments are realized following 4 steps

Step 1: Injection of pure methane gas.

Step 2: Injection of CO 2 gas. After the pressure is stable, another gas is added up to the required total pressure. The gas mixture is cooled and maintained at the operating temperature with the cryostat. It is stirred at a rate of approximately 400 rpm.

Step 3: Injection of water for crystallization of mixed CO 2 -CH 4 hydrates The stirrer is stopped. Water containing 10ppm of a lithium tracer is added. The volume of the injected solution is controlled, (in this experiment 800,15 ml). Pressure, temperature and gas compositions are monitored. After a delay corresponding to the induction time, the crystallization starts, gas is consumed and the pressure is dropped down to the equilibrium. The crystallization is occurring during hours and days and the gas phase is sampled. The liquid phase is also sampled and analyzed later.

Step 4: Dissociation of hydrates.

When the system is in equilibrium (temperature and pressure), heating the reactor by the temperature is increased by steps of 1°K until full dissociation of hydrates (Figure 3). During each step, the pressure in the reactor increases due to the dissociation of hydrates and reaches to a constant value which represents the thermodynamic equilibrium, in this time the gas and the liquid phase are sampled to determine the compositions of the phase at equilibrium. After complete dissociation of hydrate, the pressure in the reactor continues to increase but only in respect to the dilatation of gases and decrease of the solubility of gases. At each stage of dissociation, the gas and liquid phases are sampled and [ref sfgp:2013120099]-4 respectively analyzed by gas chromatography and ion exchange chromatography to determine the gas compositions in gas phase and the lithium concentration in the liquid phase.

Results

Figure 3 Evolution of pressure and temperature during crystallization and dissociation

According to mass balance of each component, the quantity of CH 4 and CO 2 in different phases can be calculated by the following equations:

n [CH 4 ] H =[ n CH 4 ]initial -n [CH 4 ] G -n [CH 4 ] L (3) n [CO 2 ] H = [ n CO 2 ]initial -n [CO 2 ] G -n [CO 2 ] L (4) n [CH4] G = n G.XCH 4 in reactor (5) n [CO2] G = n G.XCO 2 in reactor (6) n G = (7)
Where: XCO 2 and XCH 4 are the concentration of CO 2 and CH 4 in phase gas The volume of liquid in equilibrium is determined by following the concentration of a tracer (Li + ) in the liquid phase.

Li Li V V V V V ] [ ] [ 0 0 0       =            eq ini w Li Li V ] [ ] [ 1 0 (11) Where H w V , 0 w V , L eq
V are the volume of water in hydrate phase, initial and liquid phase at equilibrium. A summary of steps of calculation is presented in the figure 4: Experimentally,we can determine the volume of water (and the number of moles of water) which has reacted to form hydrate. This volume is determined from the difference of the tracer concentration of LiNO 3 being present in the liquid phase between the beginning and equilibrium (Eq. ( 10)). The number molar of gases ( n G) which have been in corporate into the lattice structure of the hydrate are calculated by equation of state (Soave-redlich-kwong). According to this way we can calculate the hydration number of there spective hydrateate quilibrium. Table 3 shows Synthetic table the calculation in one experiment, assuming that structure I is formed.

Table 3 Synthetic table the calculation in one experiment

After perform some experimental with composition difference we have compared these results are with the model from the GasHyDyn software, in the figure 5

Conclusion

This work participates to the long term issue to produce methane from methane hydrate bearing sediments. The method consists in injecting CO 2 in sediments and shifting the chemical equilibrium of methane hydrate, to form carbon dioxide clathrates hydrates. So, this thematic belongs to the topic of gas production, but also to CO 2 mitigation

In the work, we presented experimental results about the equilibrium of CO 2 -CH 4 gas mixture. Our experiments indicate that we form a solid in thermodynamic equilibrium with the gas mixture.

The next step is to study the kinetic of replacement of methane by carbon dioxide in the hydrate structure.

The work is undergoing [START_REF] Herri | Derivation of a Langmuir type of model to describe the intrinsic growth rate of gas hydrates during crystallization from gas mixtures[END_REF]. Currently, we start looking to model the formation of methane hydrate and replace it with the CO 2 in the reactor, depending on conditions such as temperature and pressure of the tank.

We produce a dispersion of crystals of methane hydrate in a wetting medium, with an excess of liquid water. It forms slurry. Then, these crystals are destabilized by CO 2 injection and we monitor the composition of the gas that is enriched in methane because of the dissociation of methane hydrate and formation of CO 2 -CH 4 hydrate following the equilibrium given in (figure 3). It will allow us, after developing a suitable model to quantify the kinetics of conversion of a solid phase to another solid phase.

After, we want to build a new device which will crystallize methane hydrates in porous media, following the procedure that we have experimented in [START_REF] Tonnet | Methane hydrates bearing synthetic sediments-Experimental and numerical approaches of the dissociation[END_REF]. These hydrates are then subjected to a flow of CO 2 through the porous medium, and then we can observe the different limitations.

This technology is a new technology developed from the needs of gas industry. In our laboratory, it participates to the development of new technologies taking profit of the physical properties of gas hydrates.

Among the new technologies we are developing, we can also cite the CO 2 capture to recover the carbon dioxide from gas power station [START_REF] Herri | Gas Hydrate Equilibria from CO 2 -N 2 and CO 2 -CH 4 gas mixtures, -Experimental studies and Thermodynamic Modelling[END_REF], Herri and Kwaterski, 2012, Kwaterski and Herri, 2012) or other industries such a steel making plants (Duc et al, 2005).

We can cite also the technology consisting in storing energy under the form of a slurry of clathrate hydrate, for air condition applications [START_REF] Darbouret | Rheological study of TBAB hydrate slurries as secondary two-phase refrigerants[END_REF][START_REF] Douzet | Prototyping a real size air conditioning system using a TBAB hydrate slurry as secondary two phase refrigerants. Experimental and modelling results[END_REF].

If clathrates hydrates can be used with a profitable objective, they still remain a problem in the oil industry because they can form crystals in the different facilities after extraction of the oil, especially in deep sea conditions. Our laboratory is also specialized in the flow assurance to prevent the nucleation, growth and agglomeration (Herri et al, 1999a, 199;Pic et al, 2000[START_REF] Seo | Hydrate phase equilibria of thecarbon dioxide, methane, and water system[END_REF][START_REF] Cournil | Asymptotic Model for Gas-Liquid Crystallization in Two-Film Systems[END_REF][START_REF] Fidel-Dufour | Experimental characterization and modelling of the rheological properties of methane hydrate slurries during their crystallisation in a water in dodecane emulsion under laminar flowing[END_REF][START_REF] Cameirao | Chord length distributions interpretation using a polydispersed population: Modeling and experiments[END_REF], and to propose solutions if a plug forms (Nguyen Hong Duc et al, 2005).

Our laboratory is also specialised on the determination of physical properties of gas hydrates, such as Hamaker constant with application to the evaluation of agglomeration (Bonnefoy et al, 2005a) or determination of refractive index (Herri and gruy, 1996;Bonnefoy et al, 2005) with application to the conception of particle sizes (Herri et al, 1999a

Résumé

Le concept de remplacement de méthane par le dioxyde de carbone dans les gisements d'hydrates de gaz naturel est considéré comme une façon possible de produire du méthane et de piéger le CO 2 [START_REF] Ohgaki | Methane exploitation by carbon dioxide from gas hydrates-phase equilibria for CO2-CH4 mixed hydrate system[END_REF][START_REF] Hirohama | Conversion of CH4 hydrate to CO2 hydrate in liquid CO2[END_REF]. Il peut être considéré comme une nouvelle technologie de production du méthane à côté des autres technologies possibles de stimulation thermique, injection chimique ou dépressurisation.

Le méthane des clathrates d'hydrates sont des solides dans lesquels l'eau se structure en un réseau de 3 dimensions créant des cavités libres pouvant piégées des molécules de gaz. Le méthane est un composant naturel dans les sédiments dont l'origine est la dégradation thermique des réservoirs fossiles ou de biodégradation des matériaux biologiques. Sous la pression et dans des conditions de mer profonde, des gisements d'hydrates de méthane se sont formés dans de nombreux endroits du monde en très grandes quantités. Le concept de production de méthane, présenté dans cette étude, consiste à injecter du CO 2 dans les gisements d'hydrates de méthane déplaçant ainsi les équilibres. Un hydrate de CO 2 (ou mixte CO 2 -CH 4 ) sera alors formé à la place de l'hydrate de méthane, et le méthane pourra être ainsi récupéré. Ce travail a comme objectif d'étudier, dans un premier temps, la thermodynamique des hydrates mixtes CO 2 -CH 4 , en présence de sels. Cela va permettre de connaitre les mécanismes de cristallisation pendant le remplacement du CH 4 par le CO 2 . Nous présentons ici, un ensemble de résultats expérimentaux qui nous a permis de valider un modèle thermodynamique qui prédit les conditions d'équilibre (pression, température) et les compositions des phases (gaz, liquide, hydrate) dans les cas des hydrates mixtes CO 2 -CH 4 . Mots-clés: Hydrates de méthane, cristallisation, gaz, CH 4 -CO 2

n[

  CH 4 ] L , n [CO 2 ] L are calculated by equation from Henry following the correlation of Holder (1980), n [CH4] H ; and n [CO2] H are the moles fractions of CH 4 and CO 2 in hydrate phase, respectively; n [CH4] G ; and n [CO2] G are the moles fractions of CH 4 and CO 2 in gas phase in the reactor; n CH 4 ; released and n CO 2 ; n G is the number of moles of gas in the reactor, which is determined by equation of state; n released is the number of moles of the released gas from the reactor. the Henry coefficient, A and B are constants in table 2

  Génie des Procédés, Numéro 104 -2013 ISSN: 1775-335X ; ISBN: 978-2-910239-78-7, Ed. SFGP, Paris, France [ref sfgp:2013120099]-

Figure 4

 4 Figure 4 Treatment of data Determination of the hydration number of the hydrates at equilibrium:Experimentally,we can determine the volume of water (and the number of moles of water) which has reacted to form hydrate. This volume is determined from the difference of the tracer concentration of LiNO 3 being present in the liquid phase between the beginning and equilibrium (Eq. (10)). The number molar of gases ( n G) which have been in corporate into the lattice structure of the hydrate are calculated by equation of state (Soave-redlich-kwong). According to this way we can calculate the hydration number of there spective hydrateate quilibrium. Table3shows Synthetic table the calculation in one experiment, assuming that structure I is formed.

  shows the phase diagram of hydrate of CO 2 -CH 4 mixture at temperature of 2,2°C. The experimental results are closed to the simulated values and new experiments are still under testing.

Figure 5

 5 Figure 5 Phase diagram of hydrate CO2 -CH4 mixture at 2.2°C

Table 1

 1 The quantities of material Experiment

	Volume Reactor ( liter)	CO2 (mol)	CH4 (mol)	Li 0	+ (mg/ml) Water (ml)
	2,36	3,026	0,9927		10,4	800,15

Table 2

 2 Constant for Henry law with correction withHolder (1980)

	Gas	A	B	vi  (cm 3 /mol)
	CH 4	15,826277 -559,0631 32
	CO 2	14,283146 -050,3269 32

Une methode de calculer la composition d'un hydrate mixte CO 2 -CH 4 , à l'équilibre et au cours de la cristallisation

	ISSN: 1775-335X ; ISBN: 978-2-910239-78-7, Ed. SFGP, Paris, France
	LE QUANG Duyen, HERRI Jean-Michel, OUABBAS Yamina,
	École Nationale Supérieure des Mines de Saint-Étienne,
	Centre SPIN, 158 cours Fauriel, 42100 St-Etienne, FRANCE
	).
	[ref sfgp:2013120099]-8

Corresponding author: E-mail address: herri@emse.fr. Tel.:+33 477 420 292


Lise Cary

Hind Benabderraziq

Jamal Elkhattabi

Marc Parmentier

Laurence Gourcy

Philippe Négrel

Tracking selenium behaviour in chalk aquifer (northern France): Sr

teaching and research institutions in France or abroad, or from public or private research centers.

Groundwaters in parts of the Paris Basin (France) are facing increasing selenium (Se) contents that can exceed the drinking water limit of 10 µg/L according to the European Framework Directive in the field of water policy (2000/60/EC). To better understand the groundwater origins and the selenium dynamics, the water chemistry of the Chalk aquifer supplying drinkable water to Lille city was studied. This area is submitted to quantitative and qualitative pressure from industrial, urban and agriculture origins.

An integrated study was settled to determine the water sources and dynamics of elements, with a focus on Se. After a large chemical characterisation of the groundwater chemistry in the four field wells, a monthly monitoring was held in four wells and in the Deûle channel. Chemical analysis of major and trace elements, stable isotopes (δ 18 O, δ 2 H), strontium isotopes, and δ 34 S and δ 18 O of sulphates were realised. The chemical composition of solids sampled at various depths at vicinity of the four wells was also analysed.

The specific geochemical signature of groundwater as revealed by Sr isotopes, in addition to element concentrations ratios like Mg/Sr and Se/Sr, highlighted mixture of three main groundwaters bodies: (1) the upstream groundwaters in the recharge area with the most radiogenic 87 Sr/ 86 Sr isotopic signature; (2) the confined groundwaters with high Sr concentrations due to water-rock interactions and the lowest 87 Sr/ 86 Sr isotopic signature close to the one of the chalk in Paris and London basins; (3) the Se-rich formations of Tertiary and Quaternary.

The contents of Se, mainly present as Se VI (and locally as Se IV ), displayed spatial and temporal disparities that can be explained by geological and hydrogeological conditions. Serich clayed sediments originating from the dismantling of Se-rich tertiary formations (i.e. Ypresian) overlay the chalk formation and can be found in saturated conditions depending of the water table level. Oxidation of Se0, Se-pyrite and Se linked to organic matters happens according to two pathways. The first one is oxidation of Se species present in the reductive clayed sediments by oxidizing groundwaters during periods of high piezometric levels. The decrease of the piezometric level induces a decrease of the Se contents in groundwaters. Negative δ 34 S (0 to -28 ‰) coherent with pyrite oxidation are observed in groundwater. The second one is linked to the infiltration of nitrate-rich recharge waters through quaternary loess which also contain Se. In this case, δ 34 S of groundwaters are slightly positive and close to the signature of fertilizers. Denitrification by pyrite or Se0 oxidation plays a great role in enhancing Se mobility in a very reactive system.

To limit the Se content in groundwaters, a precise water management is needed, e.g. mainly maintaining low piezometric levels in the chalk aquifer or avoiding mixings between contaminated wells, together with a precise knowledge of the geology of the quaternary clayed sediments.

BRGM (Bureau de Recherches géologiques et minières), 3 av. C. Guillemin, BP 36009, 45060 Orléans Cedex

2, France 2 University of Lille 1, Polytech'Lille, Laboratoire de Génie Civil et Géo-Environnement LGCgE, Lille, France


Yasser Fadlallah

Amir Khandani

Karine Amis

Abdeldjalil Aïssa-El-Bey

Ramesh Pyndiah

Precoding and Decoding in the MIMO Interference Channels for Discrete Constellation

This paper addresses the problem of decoding and precoding in the K-user MIMO interference channels. At the receiver side, a joint decoding of the interference and the desired signal is able to improve the receive diversity order. At the transmitter side, we introduce a joint linear precoding design that maximizes the joint cut-off rate, known as a tight lower bound on the joint mutual information for high signal-to-noise ratio (SNR). We also derive a closed-form solution of the precoding matrices that maximizes the mutual information when the SNR is close to zero. This solution is characterized by its low computational complexity, and only requires a local channel state information knowledge at the transmitters. Our simulation results show that decoding interference jointly with the desired signal results in a significant improvement of the receive diversity order. Also a substantial bit error rate and sum-rate improvements are illustrated using the proposed precoding designs.

I. INTRODUCTION

Recently, an intense focus of research has been steered toward the transmission in multiple-input multiple-output (MIMO) interference channel (IC). The asymptotic behavior of its capacity has been defined using joint linear precoding [START_REF] Cadambe | Interference Alignment and Degrees of Freedom of the K-user Interference Channel[END_REF], [START_REF] Maddah-Ali | Communication over MIMO X Channels: Interference Alignment, Decomposition, and Performance Analysis[END_REF]. The idea of linear precoding has been basically introduced for additional white Gaussian noise (AWGN) channel, with the goal of minimizing the weighted mean squared error [START_REF] Lee | Optimal Linear Coding for Vecotr Channels[END_REF]. This idea has been extended to the single user MIMO channel. It has been proven that linear precoding is able to achieve the channel capacity when the full CSI is provided at the transmitter [START_REF] Telatar | Capacity of multi-antenna gaussian channels[END_REF]. Later on, the authors in [START_REF] Maddah-Ali | Communication over MIMO X Channels: Interference Alignment, Decomposition, and Performance Analysis[END_REF], [START_REF] Gou | Degrees of Freedom of the K-user NxM MIMO Interference Channel[END_REF] have shown that joint linear precoding design can also achieve the asymptotic capacity of the MIMO IC. This joint precoding design is known as interference alignment (IA).

The basic idea of IA is to jointly design the signals from all transmitters such that interfering signals at each receiver overlap and remain distinct from the desired signal. At the receiver side, a zero forcing (ZF) decoder is traditionally applied to cancel the aligned interference [START_REF] Cadambe | Interference Alignment and Degrees of Freedom of the K-user Interference Channel[END_REF], [START_REF] Gomadam | Approaching the Capacity of Wireless Networks through Distributed Interference Alignment[END_REF]. However, the approach of canceling the interference using ZF results in a receive diversity loss. In contrast, decoding the interference jointly with the desired signal can improve the receive diversity order, as shown hereafter.

In this paper, we firstly discuss the achievable receive diversity of both schemes: interference decoding and traditional interference cancellation. Secondly, we aim to find out an efficient linear precoding design assuming a discrete constellation. The linear precoding scheme that maximizes the Mutual Information (MI) for discrete constellation in a single user MIMO channel has been addressed in [START_REF] Xiao | Globally Optimal Linear Precoders for Finite Alphabets Signals Over Complex Vector Gaussian Channels[END_REF]. The authors have proved that the MI is a concave function of a matrix which is itself a quadratic function of the precoding matrix.

In our work, we address the linear precoding in the multi-user MIMO IC. We do not aim to maximize the MI, but rather to maximize the joint cut-off rate between all transmitted signals. The joint cut-off rate has two characterizations:

1) it is defined as a tight lower bound on the joint mutual information for high signal-to-noise ratio (SNR), 2) it can be expressed analytically, and does not require a numerical integration. We also propose a closed form solution for precoding matrices that maximizes the MI when the signal-to-noise ratio (SNR) tends to zero. It is obtained using the first order Taylor expansion of the MI at very low SNR [START_REF] Abramowitz | Handbook of Mathematical Functions with Formulas, Graphs and Mathematical Tables[END_REF]. It is mainly characterized by its low computational complexity, and only requires local channel state information at the transmitters (CSIT).

This paper is organized as follows. In Section II, we describe the system model of the multi-user MIMO IC. In Section III, we compare the achievable receive diversity for both schemes traditional interference cancellation and interference decoding. Then we formulate the joint cutoff rate in Section IV, and we propose an iterative algorithm that converges toward a local optimal solution. The closed form solution is given in Section V. In Section VI, the simulation results evaluate our contributions. Finally, Section VII concludes the paper.

Notations: boldface upper case letters and boldface lower case letters denote matrices and vectors, respectively. For the transpose, transpose conjugate and conjugate matrix we use (.) T , (.) H and (.) * , respectively. ||.||, tr(.) and log stand for the Frobenius norm, the trace operator and log 2 . E is the expectation operator.

II. SYSTEM MODEL

Consider a K-user MIMO IC, with equal number of antennas N t at the transmitters and equal number of antennas N r at the receivers. The K transmit-receive pairs share a common channel, where each transmitter intends to have its signal decoded by its destination. The symbols of the original signal vector at transmitter j, denoted x j with dimension d j × 1, are selected with equal probabilities from a discrete constellation Q with size M shared by all transmitters. The original signal vector is then precoded using linear transformation matrix P j with dimensions N t × d j . The k th received signal is given by

y k = K j=1 H kj P j x j + n k , (1)
where H kj ∈ C Nr×Nt is the channel matrix between the j th transmitter and the k th receiver, and n k ∈ C Nr is the circularly symmetric additive Gaussian noise vector with zero mean and covariance σ 2 I Nr , i.e. n k ∼ N (0, σ 2 I Nr ). In the following, we assume that N t = N r .

Using IA design at the transmitters, where the precoders are jointly designed to align interference at all receivers [START_REF] Gomadam | Approaching the Capacity of Wireless Networks through Distributed Interference Alignment[END_REF], the received signal can be decomposed in two linearly independent subspaces as

y k = H kk P k Desired subspace .x k + j =k H kj P j interference subspace .x j + n k . (2)
The traditional decoding method consists in canceling the interference using linear ZF, and yields

ŷk = U 0 H k y k = U 0 H k H kk P k x k + U 0 H k n k = Hk x + U 0 H k n k , (3)
where the ŷ is the d k ×1 free-interference signal vector, U 0 H k is the interference canceler defined in [START_REF] Abramowitz | Handbook of Mathematical Functions with Formulas, Graphs and Mathematical Tables[END_REF], and Hk = U 0 H k H k P k with dimensions d k × d k . Another decoding method known as interference decoding, where the receiver decodes jointly the desired and the undesired signals and keeps only the desired part. In this case, the received signal can be rewritten as

ŷk = (H k1 P 1 , • • • , H kK P K )    x 1 . . . x K    + n k , = Hk x + n k ( 4
)
where Hk is the new channel matrix with dimensions N r × K j=1 d j .

In order to decode the original information in ( 3) and ( 4), we propose a joint minimum-distance (MD) detector. The MD detector involves an exhaustive search over all possible transmitted vectors. The signal is decoded as

x = arg min i || ŷk -Hk x i || 2 , ( 5
)
where ŷk is equivalent to (3) (( 4)) when the first (second) decoding method is applied, and i is over all possible transmitted vectors x in (3) and ( 4). In the next section, the receive diversity order of both aforementioned decoding schemes is discussed.

III. RECEIVE DIVERSITY ORDER

As mentioned in the previous section, our channel model can be seen as a single user MIMO channel in both cases: the ZF interference canceler and the MD-interference decoder. For a single user MIMO channel, the generic receiver equation is equal to

y DR = H DR×DT x DT ×1 + n DR . (6)
Assuming an MD decoder, the probability of error can be approximated at high SNR by [START_REF] Wang | A simple and general parameterization quantifying performance in fading channels[END_REF]]

P e = αSNR -DT DR , (7)
where α points out the horizontal shift of the P e curve, D T is the transmit diversity gain, and D R is the receive diversity gain. Assuming D T = 1, [START_REF] Xiao | Globally Optimal Linear Precoders for Finite Alphabets Signals Over Complex Vector Gaussian Channels[END_REF] indicates that the slope of the probability of error is proportional to the inverse of SNR to the power D R . D R is equal to the number of independent observations at the receiver, and remains independent of the number of transmit antennas. This result has been established in [START_REF] Zhu | Performance Analysis of Maximum Likelihood Detection in a MIMO Antenna System[END_REF], where the authors have concluded that using an MD detector, only an SNR penalty is introduced when the number of transmit antennas increases. In the following, an analogy between the single user MIMO channel and both channel models given in ( 3) and ( 4) is presented for the purpose of showing the expected receive diversity of each decoding scheme.

A. Interference cancellation using zero forcing

In order to cancel the interference, the received signal can be projected onto the interference null space as shown in (3). The projection matrix at receiver k is given using the singular value decomposition of the interference subspace as

U k S k V k = j =k H kj P j ( 8
)
where

U k = U 1 k , U 0 k consists of the interference space U 1 k ∈ C Nr×(Nr-d k )
and the interference null space U 0 k ∈ C Nr×d k ,respectively. Projecting the received signal on U 0 k , the new channel model is obtained as in [START_REF] Lee | Optimal Linear Coding for Vecotr Channels[END_REF]. This model yields a d k × d k MIMO single user channel. Hence, [START_REF] Xiao | Globally Optimal Linear Precoders for Finite Alphabets Signals Over Complex Vector Gaussian Channels[END_REF] implies that the expected receive diversity order D R is equivalent to d k , describing the free-interference subspace dimensions.

B. Interference decoding scheme

On the other hand, accounting the interference incurred from other users as a desired signal, the new channel model remains as in ( 4). This channel model can be assumed as an

N r × K j=1 d j MIMO single user channel. The signal x T 1 , • • • , x T K
T composed of the desired signal and the interference is jointly decoded using N r independent observations. Hence, the expected receive diversity gain achieved with an MD detector, is D R = N r . Consquently, as long as d k ≤ N r , a higher diversity can be achieved using the MD-interference decoding scheme at each receiver. Remark 1: The receive diversity does not depend on the precoding scheme, but rather on the number of observations at the receiver. Therefore, in the upcoming sections we consider the MD-interfernce decoding receiver that results in a higher receive diversity, and we try to optimize the precoding matrices so as to maximize the joint cut-off rate and the MI.

IV. PRECODING MATRIX OPTIMIZATION

The linear precoding optimization problem based on the maximization of the MI has been exploited in many papers for both cases: canceling interference in an IA scheme, and treating interference as noise [START_REF] Sung | Linear precoder designs for k-user interference channels[END_REF]- [START_REF] Shin | Weighted Sum Rate Maximizing Linear Transceiver Filters for K-user MIMO Interference Channel[END_REF]. The MI has been maximized under the assumption of Gaussian input distribution. However, this assumption seems to be far away from the practical systems that employ a discrete constellation such as phase shift keying (PSK) modulation and quadratic amplitude modulation (QAM). Under the discrete constellation assumption, one can propose to maximize the MI in the IC between each transmitter-receiver pair. However, it is not obvious whether the MI for discrete constellation can be expressed in closed form or not, and its computation requires numerical integration that becomes intractable with increase of the signal and constellation dimensions. On the other side, a closed form is available for the joint cut-off rate 1 that represents a lower bound on the MI that becomes very tight when the SNR increase. Therefore, we try to maximize the joint cut-off rate subject to the transmit power constraint trace(P H k P k ) ≤ d k for all k. As mentioned above, one of the advantages of this criterion is that it can be expressed analytically without the need of a numerical integration. For the sake of simplicity, we assume an equal power allocation at all transmitters. For a uniform inputs distribution, the joint cut-off rate between the k th receiver and the transmitters can be mathematically expressed as [START_REF] Massey | Coding and Modulation in Digital Communications[END_REF], [START_REF] Rey | Linear Precoder Design through Cut-off Rate Maximization in MIMO-OFDM[END_REF]

R 0 k = -log   y k i p(y k |x i ) 2 dy k   = -log   y k   i,j p(y k |x i ).p(y|x j )   dy k   a = -log A, A ′ exp   -p K j=1 H kj P j (x aj j -x a ′ j j ) 2 4σ 2   , ( 9
)
where p is the transmit power,

A ′ = {a ′ 1 , • • • , a ′ K |a ′ j ∈ M dj , ∀j} and A = {a 1 , • • • , a K |a j ∈ M dj , ∀j}.
The last step a is straightforward using the fact that y k is corrupted by a zero mean Gaussian noise, i.e. n k ∼ N c (0, σ 2 I Nr ) ( [18] example 2.6.2). It is clear from [START_REF] Wang | A simple and general parameterization quantifying performance in fading channels[END_REF] that no numerical integration is required. The maximization problem of the sum of all joint cut-off rates in the K-user MIMO IC can be achieved using the iterative algorithm given below. In this algorithm, the gradient of the sum of R 0 over all users is 1 The cut-off rate, symbolized by R 0 , characterizes a lower bound of the Gaussian channel capacity, such that for any rate less than R 0 , a system that yields an arbitrary small error probabilities exists [START_REF] Wozencraft | Principles of Communication Engineering[END_REF], [START_REF] Viterbi | Error Bounds for Convolutional Codes and an Asymptotically Optimum Decoding Algorithm[END_REF]. required, and is equal to

∂( k R0 k ) ∂P * l = -k A, A ′ exp   - p K j=1 Ĥkj (x a j ,a ′ j j ) 2 4σ 2   p ĤH lk K j=1 Ĥkj (x a j ,a ′ j j ) x a j ,a ′ j H l A, A ′ exp   - p K j=1 Ĥkj (x a j ,a ′ j j ) 2 4σ 2   , (10)
where x aj ,a ′ j j = x aj jx a ′ j j . Tending SNR to infinity, the maximization problem of the cut-off rate becomes equivalent to the one that maximizes the minimum distance of the received constellation, which yields a problem with less complexity as can be noticed from [START_REF] Wang | A simple and general parameterization quantifying performance in fading channels[END_REF].

As discussed previously, a closed form solution of the precoding matrices that maximizes the joint cut-off rate seems difficult to derive. Therefore, we propose to approach the solution iteratively. We use an iterative algorithm that optimizes one variable while the others remain fixed. Each variable is considered as one of the precoding matrices. This technique results in a non-convex optimization due to the dependence between the precoding matrices. At each iteration, the optimization is based on the gradient descent widely used in MIMO multi-user channel. The iterative algorithm is detailed as follows Calculate the gradient, ∇ P k f (P l 1 , • • • , P l K ).

5:

Update P

(l+1) k = P l k + λ.∇ P k f (P l 1 , • • • , P l K ).
, • • • , P (l+1) K ) -f (P l 1 , • • • , P l K ) > ǫ, set l = l+1
and go back to step 3), otherwise stop the processing.

In this algorithm the gradient is defined as in [START_REF] Zhu | Performance Analysis of Maximum Likelihood Detection in a MIMO Antenna System[END_REF], f describes the objective function given in [START_REF] Wang | A simple and general parameterization quantifying performance in fading channels[END_REF], and the precoding matrices are supposed to be of unit norm. The step size λ is updated using the backtracking search. Despite the nonconvexity of the multivariable objective function, as long as the variable is steered in the gradient direction, the algorithm shows a convergence to a local maximum. Remark 2: we try to maximize the joint cut-off rate rather than the cut-off rate between each transmitter and its destination, due to its simple analytic form. This entails a non-optimal performance in term of rate improvement for high SNR.

V. CLOSED FORM SOLUTION OF THE PRECODING

MATRICES

Motivated by the fact that high complexity level sustains the iterative algorithm for the objective functions described above, this section introduces a closed-form solution that tends to optimality when the SNR is close to zero. This solution is obtained using the first order Taylor expansion of the MI. In a K-user MIMO interference channel, the MI is derived as follows

I(x k ; y k ) = log 2 (M d k ) + 1 M Kd k a1,••• ,aK E [log 2 (J 1k )] -1 M Kd k a1,••• ,aK E [log 2 (J 2k )] , (11)
where

J 1k = a ′ 1 ,••• ,a ′ K exp   - z k + K j=1 Ĥkj (x aj j -x a ′ j j ) 2 σ 2   , J 2k = a ′′ 1 ,••• ,a ′′ K exp   - z k + K j =k Ĥkj (x aj j -x a ′′ j j ) 2 σ 2   ,
where Ĥkj = H kj P j , and the indices a j , a ′ j , a ′′ j are in the set M dj = 1, • • • , M dj for all j with M the constellation set length for one symbol. x aj j is a symbol vector from the j th transmitter belonging to the set Q dj .

In spite of the computational gain offered by the cut-off rate compared to the MI calculation, an important computational cost remains with the iterative algorithm. Motivated by the fact that the iterative algorithm suffers from high computational cost for the objective function described above, we introduce another solution in a closed-form that tends to optimality when the SNR is close to zero. This solution is obtained using the first order Taylor expansion of the MI. Assuming an SNR close to zero, using the first order Taylor expansion of e U and log(1 + U ) when U is in the neighborhood of zero, and using that the constellation is zero mean at all users, the maximization problem of the total MI given in [START_REF] Sung | Linear precoder designs for k-user interference channels[END_REF] can be approximated as arg max

P1,••• ,PK k I k (x k ; y k ) ≈ arg max P1,••• ,PK j trace P H j Q j P j where Q j = H H jj H jj . (12)
The d l vectors of the precoding matrix P k , k ∈ {1, • • • , K}, that maximize the function in [START_REF] Fadlallah | Interference Alignment: Improved Design via Precoding Vectors[END_REF] are given by

P k = ν d l max H H kk H kk , (13)
where ν d l max (A) are the eigenvectors corresponding to the d th l largest eigenvalues of A. The solution in (13) satisfies the transmit power constraint since the precoding matrices composed of the channel matrix eigenvectors have a unit Frobenius norm. Remark 3: the closed-form solution is characterized by its low computational complexity, and only requires at each transmitter the CSI linking it to its dedicated receiver, whereas the other solutions are achieved iteratively and require a total CSIT knowledge. • ZF-linear decoder : the interference is canceled using ZF, and the original signal is decoded using linear ZF decoder. • MMSE-linear decoder: the interference is canceled using ZF, and the original signal is decoded using linear minimum mean squared error decoder. • MD-Interference Canceler: the interference is canceled using ZF, and the original signal is decoded using minimum-distance decoder.

• MD-Interference Decoder: the interference and the desired signal are jointly decoded using minimum-distance decoder.

VI. SIMULATION RESULTS

The BER performance of the proposed solutions is evaluated in a 3-user MIMO IC. We assume an uncoded QPSK modulation for the original symbols with equal probabilities. The channel coefficients are circularly symmetric Gaussian distributed with zero mean and unit variance. Our results are based on Monte-Carlo simulations. The SNR[dB] in the figures is the SNR per antenna. The acronyms of the used designs are explained under the figures.

Fig. 1 compares the BER performance of the interference canceler using ZF versus the interference decoder using an MD detector when N t = N r = 4 and

d 1 = d 2 = d 3 = 2.
For both strategies we assume an IA design at the transmitters as described in [START_REF] Gomadam | Approaching the Capacity of Wireless Networks through Distributed Interference Alignment[END_REF], and we apply the MD detector to decode the original symbols. We observe that the MD-Interference Decoder results in a steeper slope than the MD-Interference Canceler at high SNR. E.g. at BER 10 -4 , the MD-Interference Decoder yields a gain of 7dB over the MD-Interference Canceler. This result is consistent with our expectation in Section III which says that decoding interference can achieve higher receive diversity than cancellation using ZF.

Next, Fig. 2 compares the total MI in the 3-user 2×2 MIMO IC. We assume d 1 = d 2 = d 3 = 1. The integration in (11) is solved using Monte-Carlo integration. Comparing to the no precoding, the CF design shows a significant improvement, e.g. at -5dB a gain of about 0.8 b/s/Hz is obtained. This gain increases slightly until reaching 1 b/s/Hz at 5dB. On the other hand, compared to the P-R0 design, the CF yields a gain of about 0.3 b/s/Hz in the low to medium SNR region, and then performs closely until it reaches 10dB where the P-R0 design outperforms the CF design by 0.1 b/s/Hz. • Max-SINR: the beamforming design proposed in [START_REF] Gomadam | Approaching the Capacity of Wireless Networks through Distributed Interference Alignment[END_REF] that maximizes the SINR of all streams.

IC. We assume d 1 = d 2 = d 3 = 1. The integration in (11) is solved using Monte-Carlo method. Compared to the P-R0 design that maximizes the joint cutoff rate, a slight gain is obtained using the CF design and the Max-SINR design at -10dB. With the SNR increases, the CF design results in lower sum-rate performance than the other precoding designs. In addition, up to 5dB the Max-SINR achieves the best performance and reaches a maximum gain of 0.4 b/s/Hz compared to the P-R0 at 0dB. Beyond 5dB and up to 20dB, the P-R0 outperforms the other two designs. From 20dB, the three designs performs the same.

VII. CONCLUSION

In this paper, we have addressed the problem of decoding and precoding design in the K-user MIMO IC. We have shown that decoding interference jointly with the desired signal re-sults in an improved receive diversity order. We have proposed two precoding designs, assuming discrete constellation. One seeks for the optimum that maximizes the joint cut-off rate iteratively. The other is a closed-form obtained using first order Taylor expansion of the mutual information for low SNR. In terms of sum-rate performance, the first solution has shown an improvement over the precoding design that maximizes the SINR. In terms of BER performance, the first solution has shown a gain over all other compared schemes, especially for high SNR. The second solution has shown a remarkable gain for low and medium SNR values.

Algorithm 1 1 :

 11 Precoding optimization for cut-off rate maximization Initialize randomly all precoding matrices P 1 , • • • P K .2: Start loop with l = 1 3: for k = 1 to K do 4:

Fig. 1 .

 1 Fig. 1. BER performance comparison of the interference decoding and the traditional interference cancellation scheme

Fig. 3 Fig. 2 .

 32 Fig. 2. BER performance comparison using different precoding design• NoPrecod : no precoding design.• CF : precoding based on the closed-form solution presented herein. • P-IA : precoding based on the MI maximization for Gaussian inputs constellation. • P-Gauss : precoding based on the IA concept.• P-R0 : precoding based on the joint cut-off rate maximization.• P-MILB : precoding based on the MI lower bound maximization derived in[START_REF] Zeng | A Low Complexity Design of Linear Precoding for MIMO Channels with Finite Alphabet Inputs[END_REF].

Fig. 3 .

 3 Fig. 3. Total mutual information comparison for different precoding design under BPSK modulation.


José G Moreno

Gaël Dias

Guillaume Cleuziou

Post-Retrieval Clustering Using Third-Order Similarity Measures


Tewfik Ziadi
email: ziadi@lip6.fr.

-Saher Tegane
email: s_tegane@esi.dz.

Selma Kchir
email: selma.kchir@lip6.fr.

Serge Stinckwich
email: serge.stinckwich@ird.fr

Software Product Lines for the development of Families of DSLs for Robotics

Keywords: Lignes de produits, DSL, Robotique software product lines, DSL, Robotics JLDP. Volume

The development of robotics software must deal with a large amount of variability. This position paper proposes a framework based on Software Product Line Engineering to allow defining families of DSL in the context of robotics. This framework allows: 1) Managing variability when defining the domain model for a family of DSL. 2) Deriving specific DSL according the user choices.

Introduction

Domain Specific Languages (DSLs) have emerged as a powerful mechanism for managing complexity in software development. Designing and implementing DSL is built on two main steps[MER 05] :

-Analysing the domain to capture the key concepts and their relationships. In the context of Model Driven Engineering (MDE), these concepts can be defined in what is called a domain model or a meta-model.

-Implementing related tools. This includes editors, compilers, and code generators.

As mentioned in [WHI 09], DSL are focused on very specific domain concerns and this narrow scope makes it hard to reuse a DSL for a new set of requirements. This is specially true in the context of robotics where the development of robotics software must deal with a large amount of variability from at least two perspectives :

-Variability in sensors and actuators. There is a large number of families of sensors and actuators. The choice of each type of sensors or actuators depending on the type of task target and/or the type of the used robots.

-For the same task we can find a family of existing algorithms, e.g. the bug family for the navigation task [LUM 90].

Defining a single DSL that considers these two kinds of variability can make it complex to be used. Indeed, including all types of sensors, actuators and existing algorithms overloads its domain model and makes it useless for end users that only use a subset of these sensors, actuators and algorithms. In addition, any introduction of new sensors and/or actuators needs an evolution of the domain model and its related tools. Therefore, we believe that to consider this challenge of variability management in robotics applications, we need not only to implement a single DSL but a family of DSLs[SVE 10].

Implementing a family of software is known in the software engineering community by the notion of software product lines. Software Product Line (SPL) focuses on capturing commonality and variability between several software products belonging to the same domain [CLE 01]. It introduces the notion of feature models to explicitly specify variability and proposes to automated what is called product derivation, which consists of building product members based on the selection of features from the feature model.

In this paper we propose to reuse the SPL approach to design and implement a family of DSLs for robotics application. The objective is to propose a framework that allows : 1) Specifying variability using feature models where designing the DSL family. 2) Deriving specific DSL according to the user choices and/or requirements.

The reminder of this paper is organized as follows. Section 2 presents the motivations of our work. Section 3 present our approach ans illustrates it on a simple example. Section 4 concludes this work and presents some perspectives.

Motivations

To illustrate the motivations of our approach, let's consider a simple example of DSL, called SearchRescueRobotics-DSL, in the context of search and rescue robotics. This DSL should allow specifying two types of missions :

-Exploration. This type of task involves robot maneuvers within a disaster area. The robot can sense the obstacles within its sensing area. To specify exploration mission, the DSL should allows defining distance sensors (e.g. ; Lidar, Camera) and motion actuators (e.g. ; Differential Wheels).

-Human Detection. The objective of the robot in this kind of missions is to detect the presence of human in a specific area. To specify human detection, the DSL should specify specific sensors related to human detection such as Sound and CO2 sensors. For human detection, we also need additional actuators such as ManipulatorArm in order to manipulate the environement (opening doors, providing medical care for victims, ...).

These two kinds of missions may reuse a set navigation algorithms. For instance, we consider three navigation algorithms of the Bug family[LUM 90].

Figure 1 illustrates a part of domain model (meta-model) of SearchResscueRobotics-DSL 1 . This includes four main concepts : Mission, Sensor, Actuator and NavigationAlgo. To model the different variants of sensors and actuators that can be used to define missions, we reused inheritance. For instance, Lidar, SoundSensor that inherit from Sensor represent two examples of concrete existing sensors. In the same, we modeled the different navigation algorithms using inheritance. We introduced the Na-vigationAlgo meta class and all algorithm variants are modeled using inheritance. For instance, Bug1, Bug2, and Alg1 represent meta-classes that implement Bug family members.

In the second step, when we implement the code generator that is related to the domain model of Figure 1, we need to identify each type of the concrete sensors and actuators and generate the corresponding code. The following template 2 presents the general structure of a code generator which generates the code for sensors. // generating code from sensors for ( Sensor s : sensors ) ~~i f ( s i n s t a n c e o f Lidar ) ~~~~p r i n t " s p e c i f i c code f o r Lidar ~~i f ( s i n s t a n c e o f SoundSensor ) ~~~~p r i n t " s p e c i f i c code f o r SoundSensor ........ ~~i f ( s i n s t a n c e o f CO2Sensor ) ~~~~p r i n t " s p e c i f i c code f o r CO2Sensor The design and the implementation of a single DSL for these two kinds of missions allows gathering the common concepts in the same domain model. However, it can give rise to some problems concerning the usability and the evolution of DSL. By including all concepts related to the both types of missions makes the domain model overladen and useless for users that only use one kind of missions. Indeed, for the users that only specify Exploration missions, why we include in the SearchRescueRobotics-DSL DSL concepts such as CO2Sensors or SoundSensor that are only specific for Human Detection missions ? In addition, any introduction of new concrete sensors and/or actuators needs the modification of the domain model and its related code generators. For instance, we need to modify the template below to consider new sensors and/or actuators.

In this position paper, we advocate the idea that for robotics applications, we need an approach that allows designing not a single DSL but a family of DSLs. To implement such approach, we propose to use the concepts of Software Product Line Engineering[CLE 01]. In particular, we reuse the Common Variability Language (CVL)[HAU 12] tool to implement our approach.

Next section present in detail our approach.

Our approach

Our aim is to propose a framework that allows designing a family of DSL rather than a single DSL in a specific domain. For robotics applications, this enables to derive specific DSL according to different variability factors. For instance, it is possible to obtain a specific DSL for a specific kind of missions and/or for a specific types of sensors or actuators. We identified two main requirements for such framework :

-Mechanisms to specify the variability in the domain model and in its related tools.

-Mechanisms that allow the specialization of the domain model ans its related tools to a specific user choices and/or needs.

To achieve these two requirements is to reuse Software Product Line Engineering approaches. We particularly reuse a CVL (Common Variability Language)[HAU 12] approach. Next subsections presents our framework.

Variability Management in the domain model

The Common Variability Language[HAU 12] is a generic and separate language for modeling variability in models. The main idea is to allows defining families of models. CVL is based on three based principals :

-The base model. It represent the model that gathers the basic elements among the family.

-The variability model. It represents the core model that defines variability on the base model.

-The resolution model. It is a model that defines how to resolve the variability model to create a specific model in the base model.

-The CVL transformations. CVL proposes an engine to transform the base model according to a specific resolution model.

Our idea to manage variability in the domain model is to apply CVL on the metamodels. The base model is thus defined as the meta-model part that gathers all basic concepts in the DSL family where the variability model specifies the variability between the DSL family members. Next subsections illustrates our use of CVL on the example of SearchRescueRobotics-DSL.

The Base Domain Model

The base domain model is defined as the meta-model representing the basic concepts in the DSL family. Figure 2 illustrates the base domain model for the family of DSL concerning the MobileRobot-DSL. It only contains concepts of Mission, Sensor, Actuator and NavigationAlg. We have omitted all variable concepts such as concrete sensors and actuators because these concepts are variable and depend on user choices and/or needs. The feature model (FM) in CVL is inspired from feature modeling in SPLE [KAN ]. It allows specifying the variability in the User-Centric Layer. Figure 3 shows the feature model concerning the example of SearchRescueRobotics-DSL. In addition to root feature, we introduced two sub-features Mission and NavigationAlgo to specify that we have two variation points. For Mission feature, we used the â ȂIJXORâ Ȃİ operator to specify that the user have a choice between two kinds of missions : Exploration and Human Detection. The NavigationAlgo features specify an "OR" choice between the different algorithm variants.

Fragment Substitution

. Where the feature model is defined, the next task is to specify the modification that we want apply on the base model when a specific feature is chosen. CVL proposes the use of what is called "Fragment Substitution" to define library models. This is realized by defining two types of fragments :

-Placement Fragment (PF). It consists to select the element in the base model that will be replaced when the feature is selected.

-Replacement Fragment (RF). It specifies the model fragment that will replace the selected element in CFP.

Figure 4 shows the fragment substitutions for the SeachRescueRobotics-DSL example. For each feature we define the two types of fragments. To illustrate this, let"s consider the Exploration feature : PFs definition. We defined two PFs, that are highlight in red in Figure 4, for the Exploration feature. The PF ExplorationPFSensor specifies that where the Exploration feature is selected we need to replace the element Sensor in the base model of Figure2. Indeed, to allow specifying exploration mission we need to replace the Sensor meta-class in the base domain model to include specific distance sensors. In the same, the ExplorationPFActuator PF specifies that we also need to replace the Actuator element is the base model.

RFs definition

We define for each PF an RF that specifies the model fragment that will replace the selected element. As illustrated in the Figure 4, we define ExplorationRFSensor as a RF associated with the ExplorationPFSensor. We add to the ExplorationRFSensor, a fragment of model that will be replace the Sensor meta-class where the Exploration feature is selected. This fragment model contains an inheritance hierarchy with the meta-classes Lidar and Camera that are necessary to specify exploration missions. Here also, we replace the Sensor and Actuator meta-classes.

DSL derivation

Until now, the family of DSLs is specified, thanks to CVL, as a base domain model, a feature model and a set of fragment substitutions. From this specification, we can reuse the CVL transformation engine to derive a specific DSL from the user choices. Indeed, CVL includes a set of transformations that allows creating a new model from the base model by applying substitutions defined in CF and RPs. This derivation is defined in two steps :

-Create a Resolution Model (RM). The RM can be seen as an instantiation of the feature model. It specifies the user choice and/or needs. It can be defined as a set of the selected. For instance, if we want a DSL that allows specifying only exploration mission and that supports only the Bug1 navigation algorithm, the RM will be defined with the set : Rm1= Exploration, Bug1.

-Apply CVL transformations. From the DSL family specification and a specific RM, CVL creates a new meta-model by applying substitution rules on the base metamodel. The idea is to replace each element in the PF that concerns a selected feature in the RM with the fragment model defined in its RF.

Figure 5 illustrates the domain model automatically derived using CVL transformations for the Rm1 resolution model.

In the same way, can derive the domain model of different DSLs according to the resolution model.

Variability Management in code generator

Above, we presented the use of CVL to manage varibiality in the metamodel. The second challenge is to manage variability in the code generator. In this section, we propose an approach, based on aspect oriented programming.

Conclusion and Perspectives

To handle variability in robotics application, this position paper proposes to define a family of DSLs rather than a single DSL. The idea is to allow users to derive specific DSLs according to their requirements and needs. We show that existing software product line approach can be easily reused. In this context we present a framework based on the CVL approach that allows : -Deriving specific DSL according to user requirements.

The framework is fully implemented, thanks to the CVL tool 3 . We used a simple example to illustrate our framework. Our future work spans in managing variability in the tools related to the domain model. In particular, we aim to investigate variability management in code generators.

Figure 1 .

 1 Figure 1. An example of a domain model

Figure 2 .

 2 Figure 2. The base model

Figure 4

 4 Figure 4 also illustrates the PFs and RFs for the Humain Detection feature.Here also, we replace the Sensor and Actuator meta-classes.

Figure 3 .

 3 Figure 3. The feature model

Figure 5 .

 5 Figure 5. The domain model of the DSL corresponding to the Rm1 resolution model

  [SVE 10] SVENDSEN A., ZHANG X., LIND-TVIBERG R., FLEUREY F., HAUGEN Ø., MØLLER-PEDERSEN B., OLSEN G. K., ń Developing a Software Product Line for Train Control : A Case Study of CVL ż, BOSCH J., LEE J., Eds., SPLC, vol. 6287 de Lecture Notes in Computer Science, Springer, 2010, p. 106-120. [WHI 09] WHITE J., HILL J. H., GRAY J., TAMBE S., GOKHALE A. S., SCHMIDT D. C., ń Improving Domain-Specific Language Reuse with Software Product Line Techniques ż, IEEE Software, vol. 26, n o 4, 2009, p. 47-53.


Pramod Udupa
email: pudupa@irisa.fr1

Olivier Sentieys
email: sentieys@irisa.fr1

Laurent Bramerie
email: laurent.bramerie@enssat.fr

Design and Implementation of DSP algorithms for 100 Gbps Coherent Optical-OFDM (CO-OFDM) Systems

L'OFDM optique utilise la dérection cohérente et des traitements avancés des signaux numériques pour atteindre un débit de transmission des données de 10 Gbps dans une seule sous-bande. Cette exigence stricte en débit apporte une contrainte sur le type d'algorithmes de traitement du signal et les architectures utilisées pour la construction du système. Dans cet article, une architecture parallèle et évolutive utilisant une IFFT radix-2 2 est proposée. La seconde proposition consiste en un algorithme parallèle et évolutif de synchronisation temporelle qui peut supporter des débits d'entrée très élevés au niveau du récepteur. La complexité en nombre de MOPS, ainsi que les coûts en surface vs. débit de l'algorithme de synchronisation, sont donnés pour l'émetteur-récepteur OFDM afin de montrer et caractériser les améliorations dues à l'architecture proposée. L'exploration d'architectures est ralisée en utilisant un outil de synthèse de haut niveau.

Introduction

For next generation 100 Gigabit Ethernet (100GbE), Multiband Coherent Optical-OFDM (MB CO-OFDM) is proposed to be a very good candidate for the upgrade of the core network [START_REF] Ambrosia | 100 gigabit Ethernet and beyond[END_REF]. It combines higher modulation formats (QPSK,16-QAM), coherent detection, and advanced digital signal processing (DSP) to reach higher data rates in the allocated bandwidth of 50GHz. Higher modulation formats along with Coherent Detection encode more bits per symbol and DSP algorithms help in combating nonidealities of the front-end and the optical channel.

MB CO-OFDM divides the whole system into multiple orthogonal sub-bands, where each sub-band carries different data and allows the usage of present day signal converters (DAC/ADC) to be used. DAC/ADC are still a strong limitation with respect to the bandwidth of the system since they have bandwidth of 5GHz at data precision of 6-8 bits (higher rates are available but with lower precision and a very high cost and power). Hence sub-banding of the total system relaxes constraints on the DAC/ADC and allows implementation using current DAC/ADC and FPGA/ASIC technologies. To attain the 100 Gbps data rate in the allocated 50GHz ITU bandwidth, each subband has to contribute more than 10 Gbps. A possible solution using an FFT/IFFT size of 256 and a cyclic prefix of 8 consists of using 8 sub-bands. Each sub-band has a bandwidth of 5GHz, thus utilizing a total of 40GHz. The remaining 10GHz is used for guard band and space between sub-bands.

Using dual-polarization on each sub-band and using 10% loss on spectral efficiency, forward error correction (FEC), and zero sub-carriers, the data rate per polarization and per sub-band is

D b = (1 -f ec )(1 -tr )(1 -null )(1 -cp )log 2 (M )f c
with f ec = 0.0627 corresponding to RS(255,239), null = 0.1, tr = 0.1, M = 4, f c = 5GHz, which gives D b = 7.3 Gbps. Using both horizontal and vertical polarizations, we get to 14.6 Gbps per sub-band. Using all the eight sub-bands, we finally obtain a total data rate of 118 Gbps.

The paper is arranged as follows. In Section 2, the computation complexity required for 100 Gbps CO-OFDM transceiver is evaluated. In Section 3, a parallel pipelined architecture using radix-2 2 for FFT is proposed and com-pared with other radix FFT algorithms and feedback architectures. In Section 4, a parallel architecture for initial timing synchronization is also proposed. In Section 5, scalability and computational complexity of the proposed architectures are examined. Finally, Section 6 concludes the paper.

Complexity Evaluation of a CO-OFDM Transceiver

Figure 1 shows the general architecture of a single-band CO-OFDM transceiver, which consists of a set of blocks processing binary data, fixed-point data, and analog data.

In the transmitter, binary data processing blocks mainly deal with FEC, scrambling and mapper. Then IFFT and pulse shaping are performed before the DAC. At the receiver side, the signals sampled from the ADC are first processed for synchronization and time-domain corrections before being sent to the FFT block. Finally, frequencydomain corrections, equalization, and channel decoding are performed. The required complexity for an MB CO-OFDM achieving a data rate of 100 Gbps for the transmitter is calculated. Computational complexity for 100 Gbps For 10 Gbps sub-band, the number of operations per second is given by 10 10 N ops /N . This gives 400 GOPS (giga operations per second) for the transmitter only and it doubles to more than 800 GOPS for a full transceiver. For a 100 Gbps transceiver, the total number of operations for 10 sub-bands combined is N tot = 10 × 2 × 400 GOPS = 8 TOPS (tera operations per second).

To achieve the data rate of the MB CO-OFDM with one sub-band in a single FPGA, it is therefore important to architect the algorithms in a parallel manner so as to support the input data rate from ADC at the receiver and output rate for the DAC. Using this architecture in the context of OFDM, the following choice is made for the IFFT architecture: the input is in normal order and output is in bit-reverse order.

The inputs to IFFT require complex memory array of size N 2 , while output size requires complex memory array of size 2N + N cyp to support streaming input and output, where N cyp is the length of the cyclic prefix. The total memory requirement for IFFT is N 2 +( N P -1)P +2N +N cyp to support P streaming outputs every cycle. This is the total size required for all 2-or 4-or 8-Parallel architectures. To support 4-Parallel or 8-Parallel architectures, the memory array is partitioned into smaller chunks for implementation, but the size remains the same. Thus, the proposed FFT/IFFT architecture can support very high output rates (which is an integer multiple of input clock rate) and can therefore be used for CO-OFDM. Comparison of the proposed feedforward (FF) architecture with other radix algorithms like radix-4 and also with feedback architectures is given in Tables 1 and2. Compared to FB architectures, the proposed architecture uses lesser complex adders and and lesser memory than higher radix-2 4 FFT. It also uses lesser amount of memory size compared to radix-4.

Proposed Timing Synchronization

Algorithm and Architecture

The timing synchronization [START_REF] Morelli | Synchronization Techniques for OFDMA:A Tutorial Review[END_REF] algorithm chosen here for optical OFDM is a hierarchical procedure. The first stage (coarse sync) [START_REF] Minn | A Robust Timing and Frequency Synchronization for OFDM Systems[END_REF] is a low complexity auto-correlation step and the second stage (fine sync) [START_REF] Choi | An Initial Timing Offset Estimation Method for OFDM Systems in Rayleigh Fading Channel[END_REF] is a computationally

demanding cross-correlation step. The training symbol used is [C C C -C], C = [A B], B = A * [-n].
The sign pattern [+1 + 1 + 1 -1] chosen for steep timing metric roll-off. The equations for coarse synchronization in iterative form are given by The fine synchronization algorithm operates over 2N cyp + 1 samples around η init , where N cyp is the length of cyclic prefix. The fine synchronization algorithm is given by

M init [n] = |P [n]| 2 R 2 [n] (1)
P [n + 1] =P [n] -x * [n] • x[n + M ]+ 2x * [n + 2M ] • x[n + 3M ] -x * [n + 3M ] • x[n + 4M ] (2)
R[n + 1] = R[n] + x[n + 4M ] 2 -x[n] 2 (3)
M f ine [n] = |P f ine [n]| 2 R 2 f ine [n] (4)
P f ine [n] = N 2 -1 k=0 r[n -k -1] • r[n + k] (5a) R f ine [n] = N 2 -1 k=0 r[n + k] 2 (5b)
As can be observed, R f ine can be written in an iterative form. P f ine can also be computed in a parallel manner by using the locality of memory accesses. Since P f ine uses multiplication between samples separated by fixed distance, for the second iteration, the complex sample fetched from memory for the first iteration can be reused and new index point computation can be spawned every new cycle. Thus the same memory and logic setup used for coarse synchronization stage can be completely reused for the fine synchronization stage. Thus, a scalable parallel architecture for timing synchronization is proposed which can support the high input rates of optical OFDM and quickly provide initial synchronization even in presence of large carrier frequency offset (CFO) of laser.

Results

The The implementation of parallel timing synchronization algorithm was also performed on a Virtex-6 FPGA development board. Figure 5 provides values of area in number of LUTs vs. throughput in number of clock cycles for different levels of parallelism, for a clock frequency of 200MHz. The implementation is synthesized using Cata-pultC HLS tool with the algorithm specified in C.

Conclusion

A pipelined 2/4-Parallel radix-2 2 IFFT architecture is proposed for attaining high throughput rates required for single band CO-OFDM. The gains provided by radix-2 2 in terms of reduction of computations is shown in Table 3. A reduction of 1.2 TOPS is obtained, while still using a simple butterfly compared to radix-4 IFFT. Also, scalable parallel IFFT and timing synchronization architectures are proposed. They can support the very high speed throughput requirement of optical OFDM. Figure 5 gives the area vs throughput trade-off curve which helps in choosing the optimal parallel version of the architecture.

Fig. 1 :

 1 Fig.1: Block diagram of a single-band CO-OFDM transceiver transmitter is calculated using total operations and the contribution of the IFFT (5N log 2 N ) is pre-eminent. For N = 256, the total number of operations is N ops = 6144. For 10 Gbps sub-band, the number of operations per second is given by 10 10 N ops /N . This gives 400 GOPS (giga operations per second) for the transmitter only and it doubles to more than 800 GOPS for a full transceiver. For a 100 Gbps transceiver, the total number of operations for 10 sub-bands combined is N tot = 10 × 2 × 400 GOPS = 8 TOPS (tera operations per second).To achieve the data rate of the MB CO-OFDM with one sub-band in a single FPGA, it is therefore important to architect the algorithms in a parallel manner so as to support the input data rate from ADC at the receiver and output rate for the DAC.

3Fig. 2 : 2 -Fig. 3 : 4 -

 2234 Fig. 2: 2-Parallel radix-2 2 IFFT pipelined architecture

Tab. 1 :

 1 Comparison of proposed radix-2 2 MDC 2-Parallel architecture to previously proposed architectures 4 N ) N FB (MDF) [3] Radix-2 2 8(log 4 N ) N FB (MDF) [5] Radix-2 4 8(log 4 N ) 3N/2 FF (MDC) Proposed 4(log 4 N ) N Radix-2 2 FF -feedforward, FB -feedback, MDC -Multipath Delay Commutator, MDF -Multipath Delay Feedback Tab. 2: Comparison of proposed Radix-2 2 MDC 4-Parallel architecture to previously proposed architectures ) [6] Radix-2 4 16(log 4 N ) N FF (MDC) Proposed 8(log 4 N ) N Radix-2 2

Fig. 4 :

 4 Fig. 4: Block diagram of 4-Parallel coarse synchronization The time index corresponding to the maximum value of M init [n] gives the initial estimate η init = arg max n M init [n]. The proposed architecture performs block-parallel computation. The block size is chosen to be the length of C. A 4-Parallel architecture is shown in Figure 4. The increase in memory for 4-Parallel and 8-Parallel is linear and corresponds to the increase in parallelism achieved: M = N/4.The fine synchronization algorithm operates over 2N cyp + 1 samples around η init , where N cyp is the length of cyclic prefix. The fine synchronization algorithm is given by

implementation of both 4 -

 4 Parallel and 2-Parallel IFFT architectures was synthesized using Xilinx ISE on a Virtex-6 FPGA development board. Frequencies above 400 MHz were obtained for both 2-Parallel and 4-Parallel architectures and for the realization of 5 GHz sub-band,three 4-Parallel IFFT blocks can be used in parallel to attain an output rate of 400 × 4 × 3 = 4800 Msamples per second.

Fig. 5 :

 5 Fig. 5: Area (#LUTs) vs. Throughput (clock cycles) for different levels of parallelism for the synchronization block

Tab. 3 :

 3 MOPS calculation for N = 256 size IFFT and 100Gbps O-OFDM Transceiver in Tera OPS (TOPS)Algorithm Real "×" Real "+" TOPS(GMACS)