Devis Tuia
email: devis.tuia@uv.es

Gustavo Camps-Valls
email: gustavo.camps@uv.es

Rémi Flamary
email: remi.flamary@insa-rouen.fr

Alain Rakotomamonjy

LEARNING SPATIAL FILTERS FOR MULTISPECTRAL IMAGE SEGMENTATION

We present a novel filtering method for multispectral satellite image classification. The proposed method learns a set of spatial filters that maximize class separability of binary support vector machine (SVM) through a gradient descent approach. Regularization issues are discussed in detail and a Frobenius-norm regularization is proposed to efficiently exclude uninformative filters coefficients. Experiments carried out on multiclass one-against-all classification and target detection show the capabilities of the learned spatial filters.

INTRODUCTION

In recent years, kernel methods and in particular support vector machines (SVM) [START_REF] Boser | A training algorithm for optimal margin classifiers[END_REF][START_REF] Vapnik | Statistical Learning Theory[END_REF], have been extensively used for the classification (segmentation) of multispectral, hyperspectral and radar images with excellent results [START_REF] Camps-Valls | Kernel methods in Remote Sensing Image Processing[END_REF]. Among the problems raised by satellite imagery, one of the most challenging is the one regarding the increase of spatial resolution brought by last generation satellites, named Very High Resolution (VHR), that provide imagery at metric -or even submetric-resolution. Besides the obvious advantages of such increase in image resolution, VHR imagery raises a series of signal processing problems related, in particular, to the change in the spectral signature of the classes to be detected. The finer spatial resolution implies the increase of the intraclass variance of the classes, since small objects, such as chimneys or pedestrians, become visible in the image. Therefore, the signature of a particular class becomes highly variable and robust algorithms exploiting spatial information becomes strictly necessary. Certainly, such high geometrical resolution could be considered as a noise source but undoubtedly may also be beneficial because it allows a clearer definition of shapes and edges of the objects, especially in urban settings.

A great attention has been paid to spatial filtering of remote sensing images to improve recognition. Two ba-This work is supported by the Swiss National Foundation, Grant PBLAP2-127713, and by the Spanish Ministry of Innovation and Science under projects AYA2008-05965-C04-03 and CSD2007-00018. sic families of spatio-spectral classifiers are encountered in the literature: [START_REF] Boser | A training algorithm for optimal margin classifiers[END_REF] feature extraction (or preprocessing) techniques, and (2) filtering (or post-classification) methods. In the first case, the spectral signature of a given pixel is combined with the spatial information obtained from the neighbors through window or region-based approaches, such as morphological filtering [START_REF] Soille | Morphological Image Analysis: Principles and Applications[END_REF][START_REF] Benediktsson | Classification and feature extraction for remote sensing images from urban areas based on morphological transformations[END_REF], geometrical features [START_REF] Inglada | Automatic recognition of man-made objects in high resolution optical remote sensing images by SVM classification of geometric image features[END_REF], or Markov random fields [START_REF] Dubes | Random field models in image analysis[END_REF][START_REF] Jackson | Adaptive bayesian contextual classification based on Markov random fields[END_REF]. The spatial and spectral feature vectors are then jointly classified. The second approach essentially performs spatial smoothing of a pixel-based classification map, and is often called post-regularization. The approach is mainly based on morphological operators, such as the majority voting scheme [START_REF] Tomas | Spatial postprocessing of spectrally classified Landsat data[END_REF][START_REF] Ton | Knowledge-based segmentation of Landsat images[END_REF][START_REF] Solaiman | An information fusion method for multispectral image classification postprocessing[END_REF][START_REF] Zhang | Detection of urban housing development by fusing multisensor satellite data and performing spatial feature post-classification[END_REF]. The generally low improvement in accuracy of these methods is due to the strong dependence of classifier's performance and filter design.

Considering spectral and spatial information jointly in the classifier typically leads to better results. However, extracting contextual features for classification gives rise to several problems: 1) critical parameters need to be tuned, such as the scale and extent of the spatial relations; 2) the quality of the feature extraction technique may hamper the representation capabilities of the classifier; and very importantly, 3) by combining spatial and spectral information sources, the curse of dimensionality problem is worsened. The latter problem was alleviated with the introduction of composite kernels [START_REF] Camps-Valls | Composite kernels for hyperspectral image classification[END_REF], in which dedicated kernels for the spectral and spatial information are combined. The framework has been recently extended to deal with convex combinations of kernels through multiple kernel learning [START_REF] Rakotomamonjy | Simple MKL[END_REF][START_REF] Tuia | Learning relevant image features with multiple kernel classification[END_REF]. In both cases, however, the methodology still relies on performing an ad hoc spatial feature extraction before kernel computation, typically limited to second-order statistics, Gray-level coocurrence matrices (GLCM) [START_REF] Clausi | Comparing cooccurrence probabilities and Markov random fields for texture analysis of SAR sea ice imagery[END_REF], or morphological reconstruction features [START_REF] Fauvel | Spectral and spatial classification of hyperspectral data using SVMs and morphological profiles[END_REF][START_REF] Tuia | Multi-source composite kernels for urban image classification[END_REF]. Yet successful, the approach strongly relies on user's previous knowledge about the particular image statistics to properly design the filters.

In this paper, we propose an alternative way to jointly learn (optimize) the SVM classifier and a convolution filter separating the class of interest from the others through a particular gradient-descent approach [START_REF] Flamary | Large margin filtering for signal sequence labeling[END_REF]. Since the preprocessing phase, i.e. the definition and selection of the filters, is integrated into the learning process, the proposed approach adapts the filter to the problem at hand, and is able to significantly improve classification results. The method is tested in two challenging VHR image classification scenarios: one-class target detection and multi-classification of land-cover types. The remainder of the paper is organized as follows: Section 2 describes the proposed algorithm and the optimization of the convolutional filter. Section 3 presents the data and the setup considered for the experiments discussed in Section 4. Finally, Section 5 concludes the paper.

LARGE MARGIN FILTERING FOR REMOTE SENSING IMAGES

In this section, large margin image filtering is proposed and issues regarding regularization are discussed in detail. In [START_REF] Flamary | Large margin filtering for signal sequence labeling[END_REF], the authors proposed to learn a one-dimensional convolution of the signals to enhance signal sequence labeling.

In this paper, this approach is extended to pixel classification of remote sensing multidimensional images: in this case, the filtering corresponds to learning a convolution kernel that promotes maximum margin between classes. Additionally, the the filtered SVM proposed in [START_REF] Flamary | Large margin filtering for signal sequence labeling[END_REF] is here extended to the non-linear case, through Gaussian kernels.

Definitions

Consider an image X ∈ R r1×r2×d , where r 1 × r 2 is the spatial extent -the considered pixels -and d is the number of components, that can be spectral bands, textural features or any kind of index retrieved from the image. To ease notation, we define X i,j,k = X p,k as the kth component of pixel p = (i, j). Since the approach proposed is supervised, a list of n labeled pixels p = (i, j) ∈ S l is available, as well as a corresponding labels vector

Y p ∈ {-1, 1}, ∀p ∈ S l .
As stated in the introduction, a 2-dimensional convolution may be applied to the image in order to enhance discriminative information. In this case, the convolution is performed by applying per component convolution matrices that are stored in F ∈ R f ×f ×d , where f is the size of the bidimensional square filters. When convoluting, the filter is centered on the current pixel. The filtered image X is defined by:

X p,k = f,f u=1,v=1 F u,v,k X i+u-f0,j+v-f0,k (1)
= f,f u=1,v=1 F u,v,k X p+(u,v)-f0,k
where f 0 = f /2. The filter coefficients in F define the type of smoothing/enhancement performed by the convolution.

For instance, smoothing can be achieved by applying an average filter corresponding to a matrix F = 1 f 2 , ∀(i, j, k). Finally, we define the RBF kernel between the filtered pixels p and q:

K p,q = k( X p,. , X q,. ) = exp - || X p,. -X q,. || 2 2σ 2 , ( 2
)
where σ is the kernel width or bandwidth. This kernel will be used in all the following experiments.

Optimization problem

A possible approach to contextual filtering is to extract contextual filters and then select the relevant ones by feature selection. In this paper, we propose to jointly learn the pixel classifier and the 2D convolution with large margin constraints. This double objective can be formulated with the following minimization problem:

min g,F 1 2 g 2 + C n p∈S l H(Y p , g( X p,. )) + λΩ(F) (3)
where

H(Y p , g( X p,. ) = max(0, 1 -Y p • g( X p,.
)) is the SVM hinge loss, C and λ are regularization parameters,

Ω(•) is a differentiable regularization function of F and g(•)
is the decision function defined here for the filtered pixel p:

g( X p,. ) = q∈S l α q Y q k( X q,. , X p,. ) + b, (4)
where α p are the dual variables of problem (3), k(•, •) is a valid kernel function as the RBF function defined above and b is the bias term. Note that the two leftmost terms in (3) correspond to a SVM problem for filtered pixels in X. Thus this problem implicitly depends on the variable F. This objective function is non-convex with respect to both variables. However, for a fixed F, the optimization problem with respect to g(•) is convex and reduces to a SVM problem on the filtered image. Let us denote the objective function of Eq. (3) as J(F). Then this problem reduces to minimize with regard to F:

J = J ′ (F) + λΩ(F) (5)
with:

J ′ (F) = min g 1 2 g 2 + C n p∈S l H(Y p , g( X p,. )) (6)
where g(•) is defined in Equation ( 4). Due to the strong duality of the SVM problem, J ′ (•) can be expressed in dual form as the constrained maximization problem:

max α - (p,q)∈S 2 l Y p Y q α p α q K p,q + p α p (7) s.t. C n ≥ α p ≥ 0 ∀p ∈ S l and p α p Y p = 0
with α the dual variables used in the decision function [START_REF] Soille | Morphological Image Analysis: Principles and Applications[END_REF].

According to [START_REF] Bonnans | Optimization problems with pertubation : A guided tour[END_REF], and for a given F * , J ′ (•) is differentiable with respect to F. This means that at the point F * , the gradient of J(•) can be computed.

In [START_REF] Flamary | Large margin filtering for signal sequence labeling[END_REF], the authors proposed to solve the problem in ( 5) with a gradient descent on J(F ) along F. In the non-linear case, and when an RBF kernel is used, the gradient of J ′ (•) with respect to F is:

∇J(F) u,v,k = (p,q) (X p+(u,v)-f0,k -X q+(u,v)-f0,k ) ×( X p,k -X q,k ) K p,q Y p Y q α * p α * q ,
where α * is the SVM solution for a fixed F. Note that this sum needs to be done only for the pixel p that are support vectors.

Regularization and Algorithm

The most important part of our method is the choice of the regularization term Ω(•) in ( 5), that must be easily differentiable in order to compute the gradient of J(F ). Similarly to [START_REF] Flamary | Large margin filtering for signal sequence labeling[END_REF], we use here a regularization term based on the Frobenius norm of the 3-dimensional matrix F:

Ω(F) = f,f,d u,v,k F 2 u,v,k (8)
This regularization term is differentiable and the gradient is easy to compute. Minimizing this regularization term can be seen as minimizing the filter's energy. In terms of classification, the filter matrix can be seen as a kernel parameter weighting neighboring pixels. Moreover, the Gaussian kernel defined in [START_REF] Vapnik | Statistical Learning Theory[END_REF] shows that the component-wise convolution can be seen as a scaling of the spectral channels before kernel computation. For this reason, the kernel parameter σ does not need to be optimized in the procedure and can be set to the value of 1.

The intuition of how this regularization term influences the filter learning is the following: suppose that the decision function g(•) has been learned by minimizing only J ′ (•). In this case, the learned filter matrix will simply maximize the margin between classes. But adding the Frobenius regularizer will force non-discriminative filter coefficients to vanish, thus minimizing their impact on the kernel of the neighboring pixels.

In order to solve the optimization problem, we propose a conjugate gradient (CG) descent algorithm along F with a line search for finding the optimal step. The method is detailed in Algorithm 1, where β is the CG update parameter and D i F represents the descent direction for the ith iteration. For the experiments, we used the β proposed by Fletcher and Reeves (see [START_REF] Hager | A survey of nonlinear conjugate gradient methods[END_REF]). In the implementation proposed, the iterations in the algorithm may be stopped by two criteria:

Algorithm 1 KF-SVM solver Set F u,v = 1/f for v = 1 • • • d and u = 1 • • • f Set i=0, Set D 0 F = 0 repeat i=i+1 G i F ← gradient of J ′ (F) + λΩ(F) wrt. F β ← G i F 2 G i-1 F 2 (Fletcher and Reeves) D i F ← -G i F + βD i-1 F (F i , α * ) ← Line-Search along D i F
until Stopping criterion is reached a threshold on the relative variation of J(F ) or on the norm of the variation of the coefficients in F.

Due to the non-convexity of the objective function, it is difficult to provide an exact evaluation of the solution complexity. However, we know that the gradient computation is O(n 2 s • f 2 ) with n s being the number of support vectors. When J(F ) is computed in the line search, a O(n 2 •d) SVM is solved and a O(n • f • d) filtering is applied. Note that, in order to speed up the method, warm-start initialization is used when using iteratively the SVM solver.

DATASETS AND EXPERIMENTAL SETUP

The proposed method is tested on a VHR QuickBird image of the city of Zurich, Switzerland (see Fig. 1). The dataset considered represents a residential area in the South-West part of the city. Seven classes were labeled by photointerpretation. The main challenge is to distinguish between the two classes of buildings and of roads by applying spatial filtering, because the spectral difference between these couples of classes is really low.

In the experiments, four models have been compared: i) a standard SVM, ii) a SVM that considers pixels filtered by local averaging with a f × f moving window (AvgSVM hereafter), iii) a SVM using a stacked vector of the reflectance values of the local f ×f neighborhood (WinSVM [START_REF] Flamary | Large margin filtering for signal sequence labeling[END_REF] hereafter) and iv) the proposed KF-SVM, learning a filter of the same size for each channel and classifier. From the labeled pixels available in the ground survey, 4, 000 to 5, 000 are used for training, 5, 000 for validation, and 26, 000 for testing. RBF kernels have been used in all cases. Regarding the kernel parameter σ, it has been optimized by crossvalidation for the SVM and AvgSVM methods, while for KF-SVM it is kept to a standard value of 1, since its optimal value is learned directly from the filter: as each band of the image is filtered before the computation of the kernels, it is always possible to find a regularization parameter λ returning a solution which is equivalent to the one that would have been obtained when using an optimal kernel value (see Section 2.3 for more details). In conclusion, C and σ parameters have been validated for the SVM, AvgSVM and WinSVM methods while C and λ have been validated for KF-SVM.

RESULTS AND DISCUSSION

In this Section, the experiments are presented for two classification scenarios: firstly, a binary classification problem is considered. Secondly, a multiclass setting, where as many one-against-all binary SVM as the number of classes are solved.

Binary classification

Table 1 shows the results for the binary problem with classes 'Residential buildings' and 'Buildings'. For the latter, we merged the residential and commercial buildings classes. The numerical results show good performance of KF-SVM, that strongly improves the AvgSVM result, with an increase between 4 and 6% in the estimated kappa statistic. With respect to WinSVM, the proposed method provides similar results, with a slight improvement of the Kappa coefficient, confirming the observations done for the linear case in [START_REF] Flamary | Large margin filtering for signal sequence labeling[END_REF]. Regarding the ROC curve, the AUC is only slightly improved by the proposed approach, since the detection is already very good with the mean filter. Nonetheless, by looking at the curves of Fig. 2, KF-SVM results in a global improvement of the results of SVM and AvgSVM, which can be characterized by an increase in the location of the γpoint, i.e. the closest point to the (0,1) corner of the ROC. In this sense, KF-SVM maximizes the true positive rate (tpr) while minimizing the false positive rate (f pr). Compared to WinSVM, the γ-points of the two methods are practically identical for both classification problems.

Multiclass classification

Multiclass classification results are shown in Table 2: for the 7-classes setting, the inclusion of spatial information strongly improves the results of the SVM, whose overall accuracy is increased by +8-10%. This result was expected, because the use of contextual information allows to discriminate the classes with no strong spectral differences, cf. Section 3. This can be seen in the classification maps reported in Fig. 3 mercial buildings', whose difference with respect to the residential buildings in uniquely related to the size of the objects. Moreover, KF-SVM enforces a smooth classification result and avoids noisy classification maps. With respect to the mean filter, KF-SVM provides a gain of about +2% in both accuracy and estimated kappa statistic, suggesting that an appropriate filter has been learned from data, at the same time avoiding the oversmoothing of the mean filter and taking advantage of the spatial context of the pixel. The same order of gain is observed with respect to WinSVM. Regarding the 6-classes setting, the standard SVM model leads to better results, since the size variability of buildings is no more taken into account. Nonetheless, both the filter approaches still improve the SVM result by +6-8% and the difference of about 2% between the two filters holds. For this easier 6 classes setting, WinSVM (which is convex) performs optimally and provides the best results, that slightly outperform the proposed KF-SVM (+ 0.5% in overall accuracy), which only proposes an approximated solution.

On filter weights

Unlike blackbox approaches as SVM or WinSVM, the proposed KF-SVM provides a smoothing filter maximizing margin separation. Such discriminative filters can be further used, for instance to train a second classifier using them as features. This opens a set of challenging opportunities to exploit the convolution obtained. A first step is to analyze the filter's weights, as they are learned directly from the data. Training the classifier returns a proper convolution filter whose weights represent the signal characteristics. A possible statistic summarizing the filter is the Frobenius norm of the filter per spectral channel, as discussed in Section 2.3. Figure 4 shows the variation of the Frobenius norm coefficients, in [START_REF] Jackson | Adaptive bayesian contextual classification based on Markov random fields[END_REF], during the optimization of the model for the classification of the particularly interesting 'Buildings' class. Note that, for this problem, the near infrared band is enhanced in the mixture, since it carries the most suitable information to discriminate the buildings from the vegetation classes.

CONCLUSION

In this paper, we presented an algorithm learning the appropriate convolutional filter for the classification of VHR remote sensing images. The algorithm jointly learns a SVM classifier and a 2D filtering of the image in order to maximize the margin between the classes. The algorithm has shown good results in two challenging remote sensing problems: VHR imagery multiclass classification and detection in urban settlements. KF-SVM showed similar results with respect to a contextual SVM using a stacked features approach (WinSVM). However, these results are obtained using much lower dimensional vectors and hence the mehtod is more robust to overfitting. Also, this may be of great importance if the filter is learned as a preprocessing step before using a Bayesian or GMM classifier: in this case a method using the stacked approach might fail because of the higher dimensionality induced [START_REF] Flamary | Filtrage vaste marge pour l'étiquetage séquentiel de signaux[END_REF].

Here we included the KF-SVM in one-versus-all schemes but, the optimization in one-versus-one schemes would report filters optimized to discern among pairs of classes. Such filters could help solving problems with highly overlapping classes in conjunction with standard SVM. Besides, a possible improvement of the current algorithm should take into account the local characteristics of the image, returning pointwise convolution filters. These will be studied in the future since they certainly are issues for most aerial and VHR satellite image.

Fig. 1 .

 1 Fig. 1. QuickBird scene of suburbs of Zurich (left) and labeled pixels (right). Legend: dark green = trees; light green = meadows; black = speedway; brown = roads; orange = residential buildings; red = commercial buildings; blue = shadows.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Details of the ROC curves for the (a) 'Residential Buildings' and (b) 'Buildings' targets. The dots correspond to the γ-points, γ = arg min{f pr 2 + (1 -tpr) 2 }.

Fig. 4 .

 4 Fig. 4. Variation of the per-channel Frobenius norm during the KF-SVM optimization.

Table 1 .

 1 Results for binary classification.

	Method	Class	Training #Class AUC Kappa
			Pixels Pixels
	SVM			0.904 0.638
	AvgSVM	Residential ∼ 5000 2000	0.916 0.689
	WinSVM			0.947 0.730
	KF-SVM			0.938 0.742
	SVM			0.938 0.706
	AvgSVM	Buildings * ∼ 4000 1000	0.946 0.779
	WinSVM			0.970 0.807
	KF-SVM			0.974 0.815

* Pixels from classes 'Residential' and 'Commercial'.

Table 2 .

 2 Results for multiclass classification.

	Method Classes Filter	Training	[%]OA Kappa
			size	Pixels
	SVM				75.11 0.685
	AvgSVM	7	9	∼ 5000	83.68 0.796
	WinSVM				82.98 0.785
	KF-SVM				85.32 0.816
	SVM				83.04 0.772
	AvgSVM	6 *	9	∼ 5000	89.48 0.860
	WinSVM				91.71 0.889
	KF-SVM				91.45 0.885

* Pixels from classes 'Residential' and 'Commercial'.

Acknowledgment

The authors would like to thank Prof. M. Kanevski (IGAR, University of Lausanne) for providing the QuickBird data.


Marie Ruellan

Hamid Benahmed

Bernard Multon

Christophe Josset

Aurélien Babarit

Alain H Clément

Aurelien Babarit

Alain Clement

Design Methodology for a SEAREV Wave Energy Converter

Keywords:

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. Introduction

According to the World Energy Council (WEC) [START_REF] Thorpe | Wave Energy, Chap. 15 of 2004 Survey of Energy Resources[END_REF], between 140 and 700 TWh/year of wave energy are available and economically-accessible, i.e. approximately 1% to 5% of the annual worldwide demand for electricity. The recoverable energy could reach as high as 2,000 TWh/year with more efficient conversion systems. Swells are described as the overlapping of several gradual and monochromatic elementary waves, with all phases being random. Studies have shown that the sea state, i.e. the quantity of energy contained within each of the elementary waves making up a swell, is a slowly-varying non-periodic function. The sea state may be modeled, according to the hydrodynamic community, within the frequency space f by means of an energy spectrum that depends upon two parameters, i.e.:

• the significant height (peak-to-valley) denoted H s .

This parameter corresponds to the average height of the highest third of the waves [2]; and

• the peak period of the waves T p = 1/f p , where f p corresponds to the highest peak in the frequency spectrum of the wave field (Fig. 1). The figure 1 shows a temporal measurement of the wave elevation. The significant height is indicated. The energy spectrum is also presented and the peak frequency is indicated.

For an aleatory swell, equation 1 below yields an approximation for its energy potential (in W/m):

P = ρg 2 32π H 2 s T p (1)
Fig. 1: Wave elevation measurement and wave energy spectrum where ρ is the mass density of water and g the gravitational constant.

Table I provides an order of magnitude for the recoverable power P vs. swell height H s and peak period T p .

II. Principle behind the searev system

The WEC concept is based on a pendulum set in a closed buoy actuated by the swell through excitation forces [START_REF] Babarit | Hydrodynamic optimisation and optimal control of a wave energy converter[END_REF] [4] (see Fig. 2). The lever or pendulum executes rotational movements transmitted to an active recovery system (in turn coupled to a charge via an electronic power converter), which recovers a portion of its kinetic energy produced. In order to obtain a PTO well suited to the swell characteristics associated with a particular geographic site, it is to be optimized by incorporating the actual measured swell cycles. Moreover, the design and optimization of such a system necessitate including the rather strong couplings physically existing between the hydrodynamic, mechanical and control phenomena. In an initial approach, the electromechanical part may be modeled by a simplified recovery function and confined to the recovery braking torque, whose evolution over time has been optimized in the aim of maximizing, for a given set of excitation conditions, electrical energy while minimizing the size of the conversion chain.

The computation of the electrical power output requires determining movements of the coupled device {buoy + pendulum + generator with control}. Hence, a Fig. 2: Schematic diagram of the SEAREV pendular wave energy converter [START_REF] Clement | French patent, Autonomous electrical system for wave energy conversion[END_REF] multi-physical hydrodynamic-mechanical-electrical modeling description must be derived. The movement is supposed to be plane.

The general equation to be solved is of the following form [START_REF] Babarit | Hydrodynamic optimisation and optimal control of a wave energy converter[END_REF]:

M. Ẍ = F ext (2)
where M represents the system's inertia matrix and X = [xG zG θ α] is the displacement vector (Fig. 2). F ext is the vector of generalized forces:

F ext = F p + F H + F R + F s + T R (3)
F p denotes the force being exerted by the pendulum; this force depends on X, and on the set of geometric parameters for both the buoy and the pendulum. F H stands for the hydrostatic force due to buoyancy. F R is the so-called radiation force corresponding to the reaction of the {buoy + pendulum} system on the swell. T R is the energy recovery torque. Swell excitation forces F s are calculated from a set of imposed swell conditions for a given overall buoy geometry [START_REF]MOLIN Hydrodynamic of offshore structures[END_REF].

In the case of our simplified system set-up, several simplifications have been made. Buoy movements are supposed to be small enough to linearise the equations. A monodirectional cylindric wave (in infinite depth) is considered and as a result just three swell force components on the buoy warrant attention: the horizontal force F sX , the vertical force F sZ , and the y-axis moment F sθ . At last, the fluid-structure interaction has been modelised as part of the linearized theory of potential flow [START_REF] Falnes | Ocean Waves and Oscillating Systems: linear interactions including wave-energy extraction[END_REF].

The buoy shape and the data are presented on figure 3 and in table II respectively. III. Energy conversion chain In the following, the two technological solutions intended to convert mechanical energy from the pendulum excited by the swell into electrical energy, will be presented.

A. The"hydroelectric" solution [START_REF] Josset | A Wave-to-Wire model of the SEAREV Wave Energy Converter[END_REF] In the specific case of wave energy converters, hydraulic conversion systems are often used given their suitability to wave energy applications, which display the following properties:

• Low speeds and high forces are induced by the waves.

In industry, hydraulic systems are commonly used whenever higher forces and smaller motions are required;

• Power output fluctuates in both time and amplitude.

Coupled with a pneumatic storage device, the PTO can smooth incident power fluctuations. With the hydraulic solution, the SEAREV PTO is composed of five main elements. First, a mechanical gear is used in order to increase rotational speed while decreasing the input torque. A double-effect linear hydraulic ram, connecting the gear to the floating hull, then transforms the rotation into a high-pressure form; it pumps fluid from the low-pressure tank (atmospheric pressure -1 bar) to the high pressure accumulator (250 bar), whose volume equals 1 m 3 . Energy is stored in this accumulator by means of gas compression. Once enough energy has been stored, the accumulator supplies pressurized fluid at a nominal Fig. 4: Synopsis of the "hydroelectric" conversion chain flow rate to a hydraulic motor coupled to an induction generator, which in turn is directly coupled to the grid. Electricity can then be generated and the fluid is released back into the tank at low pressure. A sample of the results obtained from this study appear in Table III

B. The All-electric solution

An all-electric solution has also been assessed for potential industrial application following an initial electrohydraulic phase (see Section III-A). The remainder of this article will lay out the design methodology for this solution. The pendulum is damped by an electromagnetic generator driven by an IGBT static converter using pulse width modulation (2 three-phase bridges, back-to-back on both the machine and network sides), in association with a system that imposes a set of optimized control laws. The generator may be directly coupled to the pendulum (i.e. direct drive) or coupled via a mechanical multiplier (both of these options are currently under study). The research presented herein concerns the optimal electromagnetic generator design and its static converter solely in the direct drive mode. The multi-physical couplings and the swell cycle complexity necessitate developing a specific design optimization methodology that incorporates the control laws.

IV. Design methodology for the all-electric solution

A. A heavily-coupled multi-physical problem

This section will address the design methodology employed for the electro-magneto-mechanical solution. The system is submitted to fluctuating swells that have previously been characterized. In order to recover the maximum Fig. 5: Synopsis of the electro-magneto-mechanical conversion chain for the all-electric solution amount of energy, certain elements need to be optimized, namely the hydrodynamic shape of the system, the electromagnetic generator and the control strategy. The coupling between system elements is strong. The buoy has been optimized by the fluid mechanics research team at the Ecole Centrale engineering school in Nantes (western France) [START_REF] Babarit | Shape optimization of the SEAREV Wave Energy Converter, 9th World Renewable Energy Congress[END_REF], while optimization of the electric generator was performed by the SATIE laboratory team working at the ENS de Cachan educational facility. This generator will basically be handled like a device capable of imposing a braking torque T R (t). Following optimization, it appears that a viscous friction type of torque shape is indeed well adapted:

T R (t) = β θ(t) (4)
The reaction of this braking torque first on the pendulum, then on both the buoy and the swell, has well been taken into account by the general model.

Based on these excitation forces and in accordance with a multi-physical hydrodynamic-mechanical-electrical model, the power and the electrical energy output are calculated at each point in time over a fixed period ∆T. Given the coupling indicated above, the choice of the control mode exerts a significant influence on overall system behavior as well as on generator design. The optimization step consists of seeking the law of instantaneous electromagnetic torque variation T R (t) that maximizes energy output and minimizes peak power. The diagram below illustrates this optimization methodology. Fig. 6: Synopsis of the swell generator design

B. Control modes

Three distinct control methods have been examined [START_REF] Ruellan | Control influence on electromagnetic generator pre-design for wave energy converter[END_REF]:

• optimization of the viscous damping coefficient β,

• power leveling, and

• latching-based control [START_REF] Falnes | Ocean Waves and Oscillating Systems: linear interactions including wave-energy extraction[END_REF] [10].

1) Optimization with constant damping coefficient β

Firstly a viscous damping type of torque is imposed:

T R (t) = β θ(t).
where β is the viscous damping coefficient, which remains constant over the full cycle (including the startup phase with the transient pendular motion ), and which must be optimized. The optimization problem consists in seeking the values of β such that the mechanical energy output W e is maximized. Figure 7 presents the mechanical energy output according to the viscous damping for two wave cycle simulations. Simulations must be run over period cycles ∆ T >> T p long enough to reach the mechanically-settled operating range. A sensitivity study of simulation time with respect to the average power calculation has been conducted. A minimum simulation time thus proves necessary before the transient state can be neglected. We set a simulation time of 800 sec for the purpose of our simulation runs (see Fig. 8).

E e = β θ(t) 2 dt (5)
Figure 9 displays the fluctuations in power output over a 800 sec cycle for a swell with a characteristic height of 3 m and period of 8 sec; this swell will constitute our reference for the remainder of the discussion. The graph shows the very strong variations in instantaneous power, which turn out to be highly disadvantageous in terms of both system cost and quality of energy produced.

Simulations were also run on various types of swells. The plots of speed vs power (Fig. 10) and speed vs torque (Fig. 11) have been traced for four types of swells.

The optimal value of damping coefficient β opt depends on the type of swell acting upon the system and must be adjusted for swell characteristics. The average power output levels (Fig. 12) along with the corresponding optimal values of recovery coefficient β opt (Fig. 13) will be presented, in the form of a scatter diagram, for various The sizable fluctuations in power output cause the electric conversion system to be oversized. A leveling of the converted power would therefore serve to better optimize the economic return. This leveling is obtained in the present case by means of modifying (reducing) the value of the recovery damping coefficient β. For those phases in which the power lies below the imposed leveling power, the value of β is held at an optimized constant, in order to maximize average power over the entire cycle. For those phases in which the power generated is greater than the leveling power, the coefficient β varies temporally such that the power generated remains constant and equal to the leveling power (i.e. generator operating at constant power).

The leveling power is defined as (α is the leveling ratio): Such a control has enabled recovering a certain amount of average power (see Fig. 14) while limiting the average power-to-peak power ratio (see Fig. 15). This power recovery is indeed characteristic of the power electronics design. The generator is not directly related to the peak power and undergoes a separate optimization procedure. We show (fig. 1617) in both the Speed, Power planes and Speed, Torque planes the set of points swept during a single swell cycle and for various leveling ratios on the reference swell. Figure 18 depicts an instantaneous recording of the power for a leveling ratio of 30%, obtained for the reference swell.

P lev = αP max before leveling
Fig. 18: Power output after power leveling (30% ratio)

3) Latching

Latching control consists of locking (latching) the motion of the body at the instant when its velocity vanishes, while waiting for the wave force to have reached the optimal phase to release the body. The body then starts moving from this initial position to the next vanishing velocity position, where it is once again latched, and so forth and so on. Instead of being a smooth, continuous function, the body position is a succession of transient motion ramps separated by resting stages. Action upon the system is therefore binary: either the body is latched, or it is free to move illustred by figure 19 and equation 6.

E e = β θ(t) 2 (1 -u(t))dt (6)
With

u(t) =   
1 as we stop the pendulum at θ = 0 0 as the pendulum is free to move

The instant of latching is imposed by the dynamics of the body itself (i.e. vanishing velocity); thus, the control variable is simply the duration of the latching phase t, or equivalently the instant of release. This mode of control, applied to the heave motion of the buoy, was proposed by Budal and Falnes [START_REF] Falnes | Ocean Waves and Oscillating Systems: linear interactions including wave-energy extraction[END_REF]. It allows to increase the pendulum position and velocity amplitude and thus, to increase the average power output. Figure 20 provides an example of results obtained under the same conditions as for the two other strategies with latching control.

C. Electromagnetic generator

We have sought to perform a preliminary design of the electromagnetic generator on the basis of results obtained from a given swell cycle (i.e. our reference cycle) and with the various types of control already presented. This design study is conducted on a classical synchronous machine structure featuring surface magnets and a radial field. The study does not allow determining the optimal machine, yet has yielded an applicable methodology [START_REF] Regnier | Optimal design of electrical engineering systems using Pareto Genetic Algorithms[END_REF] [12] [START_REF] Mueller | Current and Novel Electrical Generator Technology for Marine Renewable Energy Applications[END_REF]. The table IV presents the power output, peak power, RMS torque, maximum torque and maximum speed for the three control modes. The goal herein is to determine the set of optimal geometric characteristics for the synchronous generator that enable minimizing, from a Pareto perspective, two conflicting objectives: total losses, and the mass of active parts. This search is carried out by focusing on the torque T R (t) and rotational speed obtained during the previous system optimization steps over a given operating cycle. In what follows, we will focus solely on the directdrive solution. Total mass constitutes one criterion among others, as the cost of raw materials could be used similarly in weighting materials by their specific cost instead of their mass density. The design outcome for a swell with a characteristic height of 3 m and period of 8 sec (i.e. the reference swell) is shown in Figure 22. Depending on the site where the swell generator is located, these characteristics might not necessarily be the most severe, but merely serve as an example. As we can see, the objectives are actually contradictory. The maximum total losses are get with the latching control and are equal to 15.4 kW. At last the three Pareto front are very closed. The dimensions of machines with the same 4-kW losses can now be compared for the all three control modes [START_REF] Ruellan | Control influence on electromagnetic generator pre-design for wave energy converter[END_REF].

It should be remarked that the generator design for this particular swell is only minimally affected by the choice of control mode (as we could already suppose with the Pareto front). Latching proves to be slightly less favorable. The discrepancy is especially noteworthy on peak power and thus on the cost associated with the static converter.

V. Conclusion

This article presented firstly the wave energy converter SEAREV. Secondly, a design methodology for the allelectric solution adapted to SEAREV was proposed. Three control modes were studied, specifically one with a constant viscous damping over the swell cycle with a given set of characteristics and power leveling. This particular mode strikes a better economic compromise in designing the entire power conversion chain. Moreover, this control method seems quite straightforward to implement under real world conditions, since sea state characteristics do not in reality change abruptly, and adjusting damping based on a direct evaluation of the sea state can be easily conceived. Optimization results on the peak power design (high levels for the power electronic converter) and in terms of electromagnetic generator mass, with a very large number of poles, were compared across the various control modes. The viscous damping and leveling control therefore appear to be the most promising, although considerable work still needs to be carried out. Power leveling actually requires a field-weakening operating range, and the autopilot angle parameter is to be included in the optimization approach. The incorporation of all system execution constraints will moreover make it possible to determine which electromagnetic conversion structures and electromechanical architecture work best for optimal lever integration. The design methodology described in this article can thus be reused with the specific design models of the selected structure. Two PTO technologies have been presented: hydroelectric and the all-electric. For now, it is difficult to compare the two solutions objectively without any real system optimization.

Fig. 3 :

 3 Fig.3: 3D view of the SEAREV[START_REF] Clement | French patent, Autonomous electrical system for wave energy conversion[END_REF]

Fig. 7 :

 7 Fig. 7: Average power vs. damping coefficient β

Fig. 8 :Fig. 9 :Fig. 10 :

 8910 Fig. 8: Average power vs. simulation time

Fig. 11 :Fig. 12 :

 1112 Fig. 11: Point cloud of the instantaneous torque vs. angular speed for several swells -(a) Tp = 8s, Hs = 3m; (b) Tp = 10s, Hs = 3m; (c) Tp = 5s, Hs = 1m; (d) Tp = 5s, Hs = 2m

Fig. 14 :

 14 Fig. 14: Power output vs. leveling power α

Fig. 16 :

 16 Fig. 16: Point cloud of the instantaneous power vs. angular speed for several power leveling ratios -(a: α =1), (b: α =0.9), (c: α =0.7), (d: α =0.5), (e: α =0.3), (f: α =0.1)

Fig. 19 :

 19 Fig. 19: Zoom on the angular speed and latching function

Fig. 21 :

 21 Fig. 21: Generator architecture -A synchronous machine with radial flux and mounted surface magnets (only 3 pole pairs shown)

Fig. 22 :

 22 Fig. 22: Optimization of the three control modes in direct drive

Fig. 23 :

 23 Fig. 23: Zoom of the Pareto front presented on figure 22 A(-.-) B(-) C(.) Total active mass (tonnes) 17 17.4 18.6 Magnet mass (tonnes) 7.4 7 7 Active volume (m3) 2.0 2 2.2 External radius (m) 4.3 4 4.4 Length (m) 0.5 0.5 0.5 Number of poles pair (p) 471 491 497 <Pj+Pmg> (kW) 3.9 4 4.1

Fig. 24 :

 24 Fig. 24: Optimal geometry when employing the various control strategies ((a): constant β; b: power leveling)

TABLE I :

 I Exemplary swell energy profiles

TABLE II :

 II Dimensions of the SEAREV[START_REF] Clement | French patent, Autonomous electrical system for wave energy conversion[END_REF]

  below.

			Viscous	Hydraulic PTO
			damper
	Tp(s) Hs (m)	E extracted	E extracted	P 0 (bar)	P nominal
			(kWh)	(kWh)		(kW)
	9	1	8	8	30	300
	9	3	62	82	150	300
	9	5	130	155	190	340

TABLE III :

 III Results using hydraulic PTO

TABLE IV :

 IV Results table for the three control modes

TABLE V :

 V Optimization resultats in direct drive


Rémi Carles
email: remi.carles@math.cnrs.fr

Erwan Faou
email: erwan.faou@inria.fr

R Émi

Energy cascades for NLS on d

come   L'archive ouverte pluridisciplinaire

INTRODUCTION AND MAIN RESULT

We consider the nonlinear Schrödinger equation (1.1) i∂ t u + ∆u = λ|u| 2 u, x ∈ T d , with d 2, where the sign of λ ∈ {-1, +1} turns out to be irrelevant in the analysis below. In the present analysis, we are interested in the description of some energy exchanges between low and high frequencies for particular solutions of this equation. We will consider solutions with small initial values:

(1.2) u(0, x) = δu 0 (x),

where u 0 ∈ H 1 (T d ) and 0 < δ ≪ 1. Replacing u with δ -1 u, (1.1)-(1.2) is equivalent to

(1.3) i∂ t u + ∆u = λε|u| 2 u ; u(0, x) = u 0 (x),

where ε = δ 2 . Viewed as an infinite dimensional dynamical system in terms of the Fourier variables of the solution, such an equation is resonant in the sense that all the eigenvalues of the Laplace operator are integers only, making possible nontrivial vanishing linear combinations between the frequencies of the linear unperturbed equation (ε = 0). In such a situation, the perturbation theory cannot be directly applied as in [START_REF] Bambusi | Birkhoff normal form for some nonlinear PDEs[END_REF][START_REF] Bambusi | Birkhoff normal form for partial differential equations with tame modulus[END_REF][START_REF] Bourgain | Construction of approximative and almost periodic solutions of perturbed linear Schrödinger and wave equations[END_REF][START_REF] Craig | Periodic solutions of nonlinear Schrödinger equations and the Nash-Moser method, Hamiltonian mechanics[END_REF][START_REF] Eliasson | KAM for the nonlinear Schrödinger equation[END_REF][START_REF]A Nekhoroshev type theorem for the nonlinear Schrödinger equation on the torus[END_REF][START_REF] Grébert | Birkhoff normal form and Hamiltonian PDEs, Partial differential equations and applications[END_REF]. Let us recall that in all these works, the Laplace operator is perturbed by a typical potential making resonances generically disappear. In such situations and when u 0 (x) is smooth enough, it is possible to prove the quasi preservation of the Sobolev norms of the solution over very long time: polynomial (of order ε -r for all r) as in [START_REF] Bambusi | Birkhoff normal form for partial differential equations with tame modulus[END_REF], exponentially large as in [START_REF]A Nekhoroshev type theorem for the nonlinear Schrödinger equation on the torus[END_REF], or arbitrary large for a set of specific solutions as in [START_REF] Eliasson | KAM for the nonlinear Schrödinger equation[END_REF].

In the resonant case considered in this paper, there is a priori no reason to observe long times bounds for the Sobolev norms of the solution. Despite this fact, Eqn. (1.1) possesses many quasi-periodic solutions (see [START_REF]Quasi-periodic solutions of Hamiltonian perturbations of 2D linear Schrödinger equations[END_REF][START_REF] Wang | Quasi-periodic solutions of the Schrödinger equation with arbitrary algebraic nonlinearities[END_REF]).

2000 Mathematics Subject Classification. Primary 35Q55; Secondary 35C20, 37K55. This work was supported by the French ANR project R.A.S. (ANR-08-JCJC-0124-01). [START_REF] Bambusi | Birkhoff normal form for some nonlinear PDEs[END_REF] On the other hand, it has been recently shown in [START_REF] Colliander | Transfer of energy to high frequencies in the cubic defocusing nonlinear Schrödinger equation[END_REF] that in the defocusing case (Eqn. (1.1) with λ = 1), solutions exist exhibiting energy transfers between low and high modes which in turn induce a growth in the Sobolev norm H s with s > 1. Strikingly, such phenomenon arises despite the fact that L 2 and H 1 norms of the solution are bounded for all time.

The goal of the present work is to describe quantitatively such energy exchanges in the case of a particular explicit initial value u 0 (x) made of five low modes. Since we work on T d , the solution u takes the form

u(t, x) = j∈Z d u j (t)e ij•x ,
where u j (t) ∈ C are the Fourier coefficients of the solution, and as long as t does not exceed the lifespan of u. Here, for j = (j 1 , . . . , j d ) ∈ Z d and x = (x 1 , . . . , x d ), we have j

•x = j 1 x 1 +• • •+j d x d . We also set |j| 2 = j 2 1 +• • •+j 2 d ∈ N. Let us introduce the Wiener algebra W made of functions f on T d of the form f (x) = j∈Z d b j e ij•x such that (b j ) j∈Z d ∈ ℓ 1 (Z d ).
With this space is associated the norm

f W = j∈Z d |b j |.
Our main result is the following:

Theorem 1.1. Let d 2, and u 0 ∈ C ∞ (T d ) given by u 0 (x) = 1 + 2 cos x 1 + 2 cos x 2 .
For λ ∈ {±1}, the following holds. There exist ε 0 , T > 0 such that for 0 < ε ε 0 , (1.3) has a unique solution u ∈ C([0, T /ε]; W ), and there exist C 0 > 0 and C > 0 such that:

∀j ∈ N * , ∃c j = 0, ∀t ∈ [0, T /ε], u j (t) -c j (εt) |j| 2 -1 (C 0 εt) |j| 2 + Cε,
where the set N * is given by

(1.4) N * = {(0, ±2 p ), (±2 p , 0), (±2 p , ±2 p ), (∓2 p , ±2 p ), p ∈ N} × {0 Z d-2 }.
Arbitrarily high modes appear with equal intensity along a cascade of time layers:

∀γ ∈]0, 1[, ∀θ < 1 4 , ∀α > 0, ∃ε 1 ∈]0, ε 0 ], ∀ε ∈]0, ε 1 ], ∀j ∈ N * , |j| < α log 1 ε θ , u j 2 ε 1-γ/(|j| 2 -1) ε γ 4 .
This result expresses the possibility of nonlinear exchanges in (1.3): while the high modes in the set N * are equal to zero at time t = 0, they are significantly large in a time that depends on the mode. As this time increases with the size of the mode, this is an energy cascade in the sense of [START_REF] Cheverry | Cascade of phases in turbulent flows[END_REF]. To our knowledge, this result is the first one where such a dynamics is described so precisely as to quantify the time of ignition of different modes.

The proof of this theorem relies on the following ingredients:

• An approximation result showing that the analysis of the dynamics of (1.1) over a time of order O(1/ε) can be reduced to the study of an infinite dimensional system for the amplitudes of the Fourier coefficients u j (in a geometric optics framework). Let us mention that this reduced system exactly corresponds to the resonant normal form system obtained after a first order Birkhoff reduction (see [START_REF] Kuksin | Invariant Cantor manifolds of quasi-periodic oscillations for a nonlinear Schrödinger equation[END_REF] for the one dimensional case). We detail this connection between geometric optics and normal forms in the second section. • A careful study of the dynamics of the reduced system. Here we use the particular structure of the initial value which consists of five modes generating infinitely many new frequencies through the resonances interactions in the reduced system. A Taylor expansion (in the spirit of [START_REF] Carles | Cascade of phase shifts for nonlinear Schrödinger equations[END_REF]) then shows how all the frequencies should be a priori turned on in finite time. The particular geometry of the energy repartition between the frequencies then makes possible to estimate precisely the evolution of the particular points of the set N * in (1.4) and to quantify the energy exchanges between them. The construction above is very different from the one in [START_REF] Colliander | Transfer of energy to high frequencies in the cubic defocusing nonlinear Schrödinger equation[END_REF]. Let us mention that it is also valid only up to a time of order O(1/ε) (which explains the absence of difference between the focusing and defocusing cases). After this time, the nature of the dynamics should change completely as all the frequencies of the solution would be significantly present in the system, and the nature of the nonlinearity should become relevant.

AN APPROXIMATION RESULT IN GEOMETRIC OPTICS

For a given element (α j ) j∈Z d ∈ ℓ 1 (Z d ), we define the following infinite dimensional resonant system

(2.1) i ȧj = λ (k,ℓ,m)∈I j a k a ℓ a m ; a j (0) = α j .
where I j is the set of resonant indices (see [START_REF] Kuksin | Invariant Cantor manifolds of quasi-periodic oscillations for a nonlinear Schrödinger equation[END_REF]) associated with j defined by:

(2.2)

I j = (k, ℓ, m) ∈ Z 3d | j = k -ℓ + m, and |j| 2 = |k| 2 -|ℓ| 2 + |m| 2 .
With these notations, we have the following result:

Proposition 2.1. Let u 0 (x) ∈ W and (α j ) j∈Z d ∈ ℓ 1 (Z d ) its Fourier coefficients.
There exists T > 0 and a unique analytic solution (a j ) j∈Z d : [0, T ] → ℓ 1 (Z d ) to the system (2.1). Moreover, there exists ε 0 (T ) > 0 such that for 0 < ε ε 0 (T ), the exact solution to

(1.3) satisfies u ∈ C([0, T ε ]; W ) and there exists C independent of ε ∈]0, ε 0 (T )] such that sup 0 t T ε u -v ε (t) W Cε.
where

(2.3) v ε (t, x) = j∈Z d a j (εt)e ij•x-it|j| 2 .
Remark 2.2. Even though it is not emphasized in the notation, the function u obviously depends on ε, which is present in (1.3).

We give below a (complete but short) proof of this result using geometric optics. Let us mention however that we can also prove this proposition using a Birkhoff transformation of (1.3) in resonant normal form as in [START_REF] Kuksin | Invariant Cantor manifolds of quasi-periodic oscillations for a nonlinear Schrödinger equation[END_REF]. We give some details below. There is also an obvious connection with the modulated Fourier expansion framework developed in [START_REF] Gauckler | Nonlinear Schrödinger equations and their spectral semidiscretizations over long times[END_REF] in the non-resonant case (see also [START_REF] Van Veen | The quasi-periodic doubling cascade in the transition to weak turbulence[END_REF]).

2.1. Solution of the resonant system. The first part of Proposition 2.1 is a consequence of the following result:

Lemma 2.3. Let α = (α j ) j∈Z d ∈ ℓ 1 (Z d ).
There exists T > 0 and a unique analytic solution (a j ) j∈Z d : [0, T ] → ℓ 1 (Z d ) to the system (2.1). Moreover, there exists constants M and R such that for all n ∈ N and all s T ,

(2.4) ∀ j ∈ Z d , d n a j dt n (s) M R n n! Proof.
In [START_REF] Carles | Multiphase weakly nonlinear geometric optics for Schrödinger equations[END_REF], the existence of a time T 1 and continuity in time of the solution

a(t) = (a j (t)) j∈Z d in ℓ 1 is proved. As ℓ 1 is an algebra, a bootstrap argument shows that a(t) ∈ C ∞ [0, T 1 ]; ℓ 1 (Z d ) . From (2.1) we immediately obtain for s ∈ [0, T 1 ], ȧ(s) ℓ 1 3 a(s) 3 ℓ 1 ,
and by induction

a (n) (s) ℓ 1 3 • 5 • • • • (2n + 1) a(s) 2n+1 ℓ 1
, where a (n) (t) denote the n-th derivative of a(t) with respect to time. This implies

a (n) (0) ℓ 1 3 • 5 • • • • (2n + 1) α 2n+1 ℓ 1 α ℓ 1 n! 3 α 2 ℓ 1 n , which shows the analyticity of a for t T 2 = 1 6 α -2 ℓ 1 . The estimate (2.4
) is then a standard consequence of Cauchy estimates applied to the complex power series

n∈N 1 n! a (n) (0)z n defined in the ball B(0, 2T ) where 2T = min(T 1 , T 2 ). 2.2. Geometric optics. Let us introduce the scaling (2.5) t = εt, x = εx, u(t, x) = u ε (εt, εx) . Then (1.3) is equivalent to: (2.6) iε∂ t u ε + ε 2 ∆u ε = λε|u ε | 2 u ε ; u ε (0, x) = u 0 x ε = j∈Z d α j e ij•x/ε .
In the limit ε → 0, multiphase geometric optics provides an approximate solution for (2.6). The presence of the factor ε in front of the nonlinearity has two consequences: in the asymptotic regime ε → 0, the eikonal equation is the same as in the linear case λ = 0, but the transport equation describing the evolution of the amplitude is nonlinear. This explains why this framework is referred to as weakly nonlinear geometric optics. Note that simplifying by ε in the Schrödinger equation (2.6), we can view the limit ε → 0 as a small dispersion limit, as in e.g. [START_REF] Kuksin | Oscillations in space-periodic nonlinear Schrödinger equations[END_REF]. We sketch the approach described more precisely in [START_REF] Carles | Multiphase weakly nonlinear geometric optics for Schrödinger equations[END_REF]. The approximate solution provided by geometric optics has the form

(2.7) v ε (t, x) = j∈Z d a j (t)e φ j (t,x)/ε ,
where we demand u ε = v ε at time t = 0, that is

a j (0) = α j ; φ j (0, x) = j • x.
Plugging this ansatz into (2.6) and ordering the powers of ε, we find, for the O(ε 0 ) term:

∂ t φ j + |∇φ j | 2 = 0 ; φ j (0, x) = j • x.
The solution is given explicitly by

(2.8) φ j (t, x) = j • x -t|j| 2 .
The amplitude a j is given by the O(ε 1 ) and is given by equation (2.1) after projecting the wave along the oscillation e iφ j /ε according to the the set of resonant phases given by I j (Eqn. (2.2)). By doing so, we have dropped the oscillations of the form e i(k•x-ωt)/ε with ω = |k| 2 , generated by nonlinear interaction: the phase k • xωt does not solve the eikonal equation, and the corresponding term is negligible in the limit ε → 0 thanks to a non-stationary phase argument. Proposition 2.1 is a simple corollary of the following result that is established in [START_REF] Carles | Multiphase weakly nonlinear geometric optics for Schrödinger equations[END_REF]. We sketch the proof in Appendix A.

Proposition 2.4. Let (α j ) ∈ ℓ 1 (Z d ), and v ε be defined by (2.7) and (2.8). Then there exists ε 0 (T ) > 0 such that for 0 < ε ε 0 (T ), the exact solution to (2.6) satisfies u ε ∈ C([0, T ]; W ), where T is given by Lemma 2.3. In addition, v ε approximates u ε up to O(ε): there exists C independent of ε ∈]0, ε 0 (T )] such that

sup 0 t T u ε (t) -v ε (t) W Cε.

Link with normal forms.

Viewed as an infinite dimensional Hamiltonian system, (1.3) can also be interpreted as the equation associated with the Hamiltonian

H ε (u, ū) = H 0 + εP := j∈Z d |j| 2 |u j | 2 + ε λ 2 k+m=j+ℓ u k u m ūℓ ūj .
that is iu j = ∂ ūj H(u, ū), see for instance the presentations in [START_REF] Bambusi | Birkhoff normal form for some nonlinear PDEs[END_REF][START_REF] Grébert | Birkhoff normal form and Hamiltonian PDEs, Partial differential equations and applications[END_REF] and [START_REF]A Nekhoroshev type theorem for the nonlinear Schrödinger equation on the torus[END_REF].

In this setting, the Birkhoff normal form approach consists in searching a transformation τ (u) = u + O(εu 3 ) close to the identity over bounded set in the Wiener algebra W , and such that in the new variable v = τ (u), the Hamiltonian

K ε (v, v) = H ε (u, ū) takes the form K ε = H 0 + εZ + ε 2 R
(|k| 2 + |m| 2 -|j| 2 -|ℓ| 2 )χ kmℓj + Z kmℓj = P kmℓj ,
where

P kmℓj = 1 if k + m -j -ℓ = 0, 0 otherwise.
The solvability of this equation relies precisely on the resonant relation |k| 2 + |m| 2 = |j| 2 + |ℓ| 2 : For non resonant indices, we can solve for χ kmℓj and set Z kmℓj = 0, while for resonant indices, we must take Z kmℓj = P kmℓj . Note that here there is no small divisors issues, as the denominator is always an integer (0 or greater than 1).

Hence we see that up to O(ǫ 2 ) terms as long as the solution remains bounded in W , the dynamics in the new variable will be close to the dynamics associated with the Hamiltonian

K ε 1 (u, ū) = H 0 + εZ := j∈Z d |j| 2 |u j | 2 + ε λ 2 k+m=j+ℓ |k| 2 +|m| 2 =|j| 2 +|ℓ| 2 u k u m ūℓ ūj .
At this point, let us observe that H 0 and Z commute: {H 0 , Z} = 0, and hence the dynamics of K ε 1 is the simple superposition of the dynamics of H 0 (the phase oscillation (2.8)) to the dynamics of εZ (the resonant system (2.1)). Hence we easily calculate that v ε (t, x) defined in (2.3) is the exact solution of the Hamiltonian K ε 1 . In other words, it is the solution of the first resonant normal form of the system (1.3).

The approximation result can then easily be proved using estimates on the remainder terms (that can be controlled in the Wiener algebra, see [START_REF]A Nekhoroshev type theorem for the nonlinear Schrödinger equation on the torus[END_REF]), in combination with Lemma 2.3 which ensures the stability of the solution of K ε 1 and a uniform bound in the Wiener algebra over a time of order 1/ε.

AN ITERATIVE APPROACH

We now turn to the analysis of the resonant system (2.1). The main remark for the forthcoming analysis is that new modes can be generated by nonlinear interaction: we may have a j = 0 even though α j = 0. We shall view this phenomenon from a dynamical point of view. As a first step, we recall the description of the sets of resonant phases, established in [START_REF] Colliander | Transfer of energy to high frequencies in the cubic defocusing nonlinear Schrödinger equation[END_REF] in the case d = 2 (the argument remains the same for d 2, see [START_REF] Carles | Multiphase weakly nonlinear geometric optics for Schrödinger equations[END_REF]): Lemma 3.1. Let j ∈ Z d . Then, (k, ℓ, m) ∈ I j precisely when the endpoints of the vectors k, ℓ, m, j form four corners of a non-degenerate rectangle with ℓ and j opposing each other, or when this quadruplet corresponds to one of the two following degenerate cases:

(k = j, m = ℓ), or (k = ℓ, m = j).
As a matter of fact, (the second part of) this lemma remains true in the onedimensional case d = 1. A specifity of that case, though, is that the associated transport equations show that no mode can actually be created [START_REF] Carles | Multiphase weakly nonlinear geometric optics for Schrödinger equations[END_REF]. The reason is that Lemma 3.1 implies that when d = 1, (2.1) takes the form i ȧj = M j a j for some (smooth and real-valued) function M j whose exact value is unimportant: if a j (0) = 0, then a j (t) = 0 for all t. In the present paper, on the contrary, we examine precisely the appearance of new modes.

Introduce the set of initial modes:

J 0 = {j ∈ Z d | α j = 0}.
In view of (2.1), modes which appear after one iteration of Lemma 3.1 are given by:

J 1 = {j ∈ Z d \ J 0 | ȧj (0) = 0}.
One may also think of J 1 in terms of Picard iteration. Plugging the initial modes (from J 0 ) into the nonlinear Duhamel's term and passing to the limit ε → 0, J 1 corresponds to the new modes resulting from this manipulation. More generally, modes appearing after k iterations exactly are characterized by:

J k = j ∈ Z d \ j-1 ℓ=0 J ℓ | d k dt k a j (0) = 0 .

A PARTICULAR DYNAMICAL SYSTEM

We consider the initial datum (4.1)

u 0 (x) = 1 + 2 cos x 1 + 2 cos x 2 = 1 + e ix 1 + e -ix 1 + e ix 2 + e -ix 2 .
The corresponding set of initial modes is given by

J 0 = {(0, 0), (1, 0), (-1, 0), (0, 1), (0, -1)} × {0 Z d-2 }.
It is represented on the following figure:

- 6 v v v v v
In view of Lemma 3.1, the generation of modes affects only the first two coordinates: the dynamical system that we study is two-dimensional, and we choose to drop out the last d -2 coordinates in the sequel, implicitly equal to 0 Z d-2 . After one iteration of Lemma 3.1, four points appear:

J 1 = {(1, 1), (1, -1), (-1, -1), (-1, 1)}, as plotted below. - 6 v v v v v v v v v
The next two steps are described geometrically:

- 6 v v v v v v v v v v v v v - 6 v v v v v v v v v v v v v v v v v v v v v v v v v
As suggested by these illustrations, we can prove by induction:

Lemma 4.1. Let p ∈ N.
• The set of relevant modes after 2p iterations is the square of length 2 p whose diagonals are parallel to the axes:

N (2p) := 2p ℓ=0 J ℓ = {(j 1 , j 2 ) | |j 1 | + |j 2 | 2 p } .
• The set of relevant modes after 2p + 1 iterations is the square of length 2 p+1 whose sides are parallel to the axes:

N (2p+1) := 2p+1 ℓ=0 J ℓ = {(j 1 , j 2 ) | max(|j 1 |, |j 2 |) 2 p } .
After an infinite number of iterations, the whole lattice Z 2 is generated:

k 0 N (k) = Z 2 × {0 Z d-2 }.
Among these sets, our interest will focus on extremal modes: for p ∈ N,

N (2p) * := {(j 1 , j 2 ) ∈ {(0, ±2 p ), (±2 p , 0)}} , N (2p+1) * := {(j 1 , j 2 ) ∈ {(±2 p , ±2 p ), (∓2 p , ±2 p )}} .
These sets correspond to the edges of the squares obtained successively by iteration of Lemma 3.1 on J 0 . The set N * defined in Theorem 1.1 corresponds to

N * = k 0 N (k) * .
The important property associated to these extremal points is that they are generated in a unique fashion: are generated in a non-unique fashion by the interaction of modes in Z d . For instance, (1, 1) ∈ J 1 is generated after one step only by the interaction of (0, 0), (1, 0) and (0, 1). On the other hand, we see that after two iterations, (1, 1) is fed also by the interaction of the other three points in N

(1) * , (-1, 1), (-1, -1) and (1, -1). After three iterations, there are even more three waves interactions affecting (1, 1). Remark 4.3. According to the numerical experiment performed in the last section, it seems that all modes -and not only the extremal ones -receive some energy in the time interval [0, T /ε]. However the dynamics for the other modes is much more complicated to understand, as non extremal points of N (n+1) are in general generated by several triplets of points in N (n) .

PROOF OF THEOREM 1.1

Since the first part of Theorem 1.1 has been established at the end of §2, we now focus our attention on the estimates announced in Theorem 1.1.

In view of the geometric analysis of the previous section, we will show that can compute the first non-zero term in the Taylor expansion of solution a m (t) of (2.1) at t = 0, for m ∈ N * . Let n 1 and j ∈ N (n) * . Note that since we have considered initial coefficients which are all equal to one -see (4.1) -and because of the symmetry in (2.1), the coefficients a j (t) do not depend on j ∈ N (n) * but only on n. Hence we have

a j (t) = t α(n) α(n)! d α(n) a j dt α(n) (0) + t α(n)+1 α(n)! 1 0 (1 -θ) α(n) d α(n)+1 a j dt α(n)+1 (θt)dθ,
for some α(n) ∈ N still to be determined.

First, Eqn. (2.4) in Lemma 2.3 ensures that there exists C 0 > 0 independent of j and n such that

r j (t) = t α(n)+1 α(n)! 1 0 (1 -θ) α(n) d α(n)+1 a j dt α(n)+1 (θt)dθ satisfies (5.1) |r j (t)| (C 0 t) α(n)+1 .
Next, we write

(5.2) a j (t) = c(n)t α(n) + r j (t),
and we determine c(n) and α(n) thanks to the iterative approach analyzed in the previous paragraph. In view of Lemma 4.2, we have

i ȧj = 2λc(n -1) 2 t 2α(n-1) + O t 2α(n-1)+1 ,
where the factor 2 accounts for the fact that the vectors k and m can be exchanged in Lemma 4.2. We infer the relations:

α(n) = 2α(n -1) + 1 ; α(0) = 0. c(n) = -2iλ c(n -1) 2 2α(n -1) + 1 ; c(0) = 1.
We first derive α(n) = 2 n -1.

We can then compute, c(1) = -2iλ, and for n 1:

c(n + 1) = i (2λ) n k=0 2 k n+1 k=1 (2 k -1) 2 n+1-k = i (2λ) 2 n+1 -1 n+1 k=1 (2 k -1) 2 n+1-k .
We can then infer the first estimate of Theorem 1.1: by Proposition 2.4, there exists C independent of j and ε such that for 0 < ε ε 0 (T ),

|u j (t) -a j (εt)| Cε, 0 t T ε .
We notice that since for j ∈ N

(n) * , |j| = 2 n/2 , regardless of the parity of n, we have α(n) = |j| 2 -1. For j ∈ N * , we then use (5.2) and (5.1), and the estimate follows, with c j = c(n).

To prove the last estimate of Theorem 1.1, we must examine more closely the behavior of c(n). Since 2 k -1 2 k for k 1, we have the estimate

|c(n + 1)| 2 2 n+1 -1 2 n+1 k=1 k2 n+1-k . Introducing the function f n+1 (x) = n+1 k=1 x k = 1 -x n+1 1 -x x, x = 1,
we have

n+1 k=1 k2 n+1-k = 2 n f ′ n+1 1 2 = 2 n+2 -n -3,
and the (rough) bound

|c(n + 1)| 2 2 n+1 -1-2 n+2 +n+3 = 2 -2 n+1 +n+2 2 -2 n+1 .
We can now gather all the estimates together:

|u j (t)| c(n) (εt) α(n) -(C 0 εt) α(n)+1 -Cε 1 2 εt 2 2 n -1 -(C 0 εt) 2 n -Cε 1 2 εt 2 2 n -1 1 -(2C 0 ) 2 n εt -Cε. (5.3)
To conclude, we simply consider t such that (5.4) εt 2

2 n -1 = ε γ , that is t = 2 ε 1-γ/α(n) .
Hence for the time t given in (5.4), since α(n) = |j| 2 -1, we have

(2C 0 ) 2 n εt = (2C 0 ) |j| 2 ε γ/(|j| 2 -1) = exp |j| 2 log(2C 0 ) - γ |j| 2 -1 log 1 ε .
Assuming the spectral localization

|j| α log 1 ε θ ,
we get for ε small enough,

(2C 0 ) 2 n εt exp α 2 log 1 ε 2θ log(2C 0 ) - γ α 2 log 1 ε 1-2θ
.

The argument of the exponential goes to -∞ as ε → 0 provided that γ > 0 and θ < 1 4 , in which case we have 1 -(2C 0 ) 2 n εt > 3/4 for ε sufficiently small. Inequality (5.3) then yields the result, owing to the fact that Cε is negligible compared to ε γ when 0 γ < 1.

NUMERICAL ILLUSTRATION

We consider the equation (1.1) in the defocusing case (λ = 1) on the twodimensional torus, d = 2. We take u(0, x) = δ(1 + 2 cos(x 1 ) + 2 cos(x 2 )) with δ = 0.0158. With the previous notations, this corresponds to 1/ε = δ -2 ≃ 4.10 3 . In Figure 1, we plot the evolution of the logarithms of the Fourier modes log |u j (t)| for j = (0, n), with n = 0, . . . , 15. We observe the energy exchanges between the modes. Note that all the modes (and not only the extremal modes in the set N * ) gain some energy, but that after some time there is a stabilization effect (all the modes are turned on) and the energy exchanges are less significant. In contrast, we plot in Figure 2 the solution corresponding to the initial value u(0, x) = δ(2 cos(x 1 ) + 2 cos(x 2 )) with the same δ. In this situation, no energy exchanges are observed after a relatively long time. Note that in this case, the initial data is made of the 4 modes {j ∈ Z 2 ||j| = 1 } forming a square in Z 2 . This set is closed for the resonance relation, so no energy exchange is expected in the time scale O(1/ε). We notice that the solution of (2.1) is given explicitly by a j (t) = 0 for |j| = 1, and a j (t) = exp(9it) for |j| = 1.

The numerical scheme is a splitting time integrator based on the decomposition between the Laplace operator and the nonlinearity in combination with a Fourier pseudo-spectral collocation method (see for instance [START_REF] Lubich | On splitting methods for Schrödinger-Poisson and cubic nonlinear Schrödinger equations[END_REF] and [START_REF] Faou | Geometric integration of Hamiltonian PDEs and applications to computational quantum mechanics[END_REF]Chap IV] for convergence results in the case of (1.1)). While the Laplace operator part i∂ t u = -∆u can be integrated exactly in Fourier, the solution of the nonlinear part i∂ t u = |u| 2 u starting in v(x) is given explicitly by the formula u(t, x) = exp(-it|v(x)| 2 )v(x). The fast Fourier transform algorithm allows an easy implementation of the algorithm. The stepsize used is τ = 0.001 and a 128 × 128 grid is used.

Note that using the framework of [START_REF] Faou | Hamiltonian interpolation of splitting approximations for nonlinear pdes[END_REF][START_REF] Faou | Geometric integration of Hamiltonian PDEs and applications to computational quantum mechanics[END_REF], we can prove that the numerical solution can be interpreted as the exact solution of a modified Hamiltonian of the where ω kmℓj = |k| 2 +|m| 2 -|ℓ| 2 -|m| 2 and B K the grid of frequencies j = (j 1 , j 2 ) such that j 1 and j 2 are less than K/2 = 64. Note that this energy is well defined as τ ω kmℓj is never a multiple of 2π, and that the frequencies of the linear operator of this modified energy carry on the same resonance relations (at least for relatively low modes). This partly explains why the cascade effect due to the resonant system should be correctly reproduced by the numerical simulations.

APPENDIX A. SKETCH OF THE PROOF OF PROPOSITION 2.4 By construction, the approximate solution v ε solves

iε∂ t v ε + ε 2 ∆v ε = λε|v ε | 2 v ε + λεr ε ,
where the source term r ε correspond to non-resonant interaction terms which have been discarded:

r ε (t, x) = j∈Z d (k,ℓ,m) ∈I j
a k (t)a ℓ (t)a m (t)e i(φ k (t,x)-φ ℓ (t,x)+φm(t,x))/ε . The error term w ε = u εv ε solves

iε∂ t w ε + ε 2 ∆w ε = λε |w ε + v ε | 2 (w ε + v ε ) -|v ε | 2 v ε -λεr ε ; w ε |t=0 = 0.
By Duhamel's principle, this can be recasted as Proposition 2.4 then follows from one integration by parts (integrate the exponential), along with (A.1) and Lemma 2.3.

w ε (t) = -iλ

  where Z is expected to be as simple as possible and R = O(u 6 ). Searching τ as the time t = ε flow of an unknown Hamiltonian χ, we are led to solving the homological equation {H 0 , χ} + Z = P, where { • , • } is the Poisson bracket of the underlying (complex) Hamiltonian structure. Now with unknown Hamiltonians χ(u, ū) = k,m,j,ℓ χ kmℓj u k u m ūℓ ūj and Z(u, ū) = k,m,j,ℓ Z kmℓj u k u m ūℓ ūj , the previous relation can be written

Lemma 4 . 2 .

 42 Let n 1, and j ∈ N (n) * . There exists a unique pair (k, m) ∈ N (n-1) × N (n-1) such that j is generated by the interaction of the modes 0, k and m, up to the permutation of k and m. More precisely, k and m are extremal points generated at the previous step: k, m ∈ N (n-1) * . Note however that points in N (n) *

FIGURE 1 .

 1 FIGURE 1. Evolution of the Fourier modes of the resonant solution in logarithmic scale: u 0 (x) = 1 + 2 cos(x 1 ) + 2 cos(x 2 )

FIGURE 2 . 2 (

 22 FIGURE 2. Evolution of the Fourier modes of the nonresonant solution in logarithmic scale: u 0 (x) = 2 cos(x 1 ) + 2 cos(x 2 )

  We writeφ k (t, x)φ ℓ (t, x) + φ m (t, x) = κ k,ℓ,m • xω k,ℓ,m t, with κ k,ℓ,m ∈ Z d , ω k,ℓ,m ∈ Z and |κ k,ℓ,m | 2 = ω k,ℓ,m , hence (A.1) |κ k,ℓ,m | 2ω k,ℓ,m 1.

t 0 eiλ t 0 e 0 e

 000 iε(t-s)∆ |w ε + v ε | 2 (w ε + v ε ) -|v ε | 2 v ε (s)ds + iε(t-s)∆ r ε (s)ds. Denote R ε (t) = t iε(t-s)∆ r ε (s)ds.Since W is an algebra, and the norm inW controls the L ∞ -norm, it suffices to prove R ε L ∞ ([0,T ];W ) = O(ε). We compute R ε (t, x) = j∈Z d (k,ℓ,m) ∈I j b k,ℓ,m (t, x), where b k,ℓ,m (t, x) = t 0 a k (s)a ℓ (s)a m (s) exp i κ k,ℓ,m • x + |κ k,ℓ,m | 2 sω k,ℓ,m s ε ds.

Acknowledgements. The authors wish to thank Benoît Grébert for stimulating discussions on this work.