G N Kariniotakis
email: kariniotakis@cenerg.cma.fr

D Mayer

J A Halliday

A G Dutton

A D Irving

R A Brownsword

P S Dokopoulos

M C Alexiadis

Load, Wind and Hydro Power Forecasting Functions of the More-Care EMS System

Keywords: Wind forecasting, water inflows prediction, load forecasting, numerical weather predictions, SCADA, EMS

The objective of the MORE-CARE project has been to develop an advanced Energy management System (EMS) to assist the operators of a medium or a large isolated power system. MORE-CARE focuses on the case of systems including substantial amount of renewable energy sources and namely wind, and hydro power. The paper presents the functions developed in the frame of this project for the prediction of the wind production, the hydro reservoirs inflows a nd the total power system load. The manmachine interface characteristics related to these forecasting functions are presented together with preliminary results on their on-line performance.

I. INTRODUCTION

The general scope of the More-Care Energy Management System (EMS) is to contribute in the achievement of an optimal utilization of multi-renewable energy systems, in a wide variety of medium and large size isolated systems with diverse structures and operating conditions. Such power systems may be characterized by:

• Several conventional fossil-fuelled generation plants;

• Multiple renewable sources (photovoltaic, wind, hydro units with or without storage facilities, geothermal, etc.)

• Privately owned plants. Consideration of this nondispatchable power changes substantially the control conditions of the network, namely the dispatching algorithms and controllable variables. The More-Care system aims to assist the operators of such systems by proposing optimal operating scenarios for the various power units, as well as the various actions needed to avoid dangerous situations, which might result from a poor prediction of load or weather or pre-selected disturbances.

The optimal operation scenarios for the conventional and renewable units are generated by a unit commitment (UC) and an economic dispatch (ED) modules. A dynamic security assessment (DSA) module supervis es the operation of the UC and ED modules with the objective of producing schedules for the power units that respect security constraints.

The above functions require, as input, forecasts for:

• the total load of the power system,

• the load at the level of bus-bars,

• the production of the wind parks,

• and finally of the water inflows at the reservoirs of the hydro installations. The ED and UC functions operate in two cycles:

• a fast one, covering a short-term horizon up to 4-8 hours ahead with a "basic" time-step of 10-20 m inutes. This cycle is performed with a frequency equal to the basic time-step. The aim is to update operation schedules according to most recent on-line information especially about renewable sources. This cycle is critical for the optimal operation of fast units like gas turbines and diesel engines and for taking preventive measures when penetration is high.

• Secondly, More-Care operates with a slow cycle, covering a horizon up to 48 hours ahead with hourly time-steps. This cycle is performed every hour. A long-term UC function is then performed to schedule slow units like steam turbines or storage.

The different temporal characteristics of the above two cycles involve development of appropriate forecasting functions for each cycle. Functions developed for the fast cycle will be mentioned hereafter as "short-term", while forecasting for the slow cycle will be referred to as "longterm".

The paper presents the methodologies applied to provide solutions to each prediction problem and the basic characteristics of the developed on-line modules. The implementation of these modules and their integration within the More-Care platform is also described. The software has been installed for on-line operation at the islands of Crete (Greece) and Madeira (Portugal). A third application concerns Ireland, where More-Care has been configured to operate as a stand-alone wind forecasting platform for the prediction of the output of 11 wind farms. Figure 1 shows the global modular architecture of the More-Care system with the integrated forecasting modules.

II. IMPLEMENTATION OF ON-LINE MODULES.

Due to the diverse needs of targeted medium and large-scale systems, the software developed is highly modular, allowing integration o f the options that are best suited to the particularities of each system. The developed functions are characterized by a "plug-and-play" functionality. To achieve this, it was necessary to develop a relational database for accessing on-line data, handling historical data and input/output of the models. Plug-and-play and modular properties are useful in the implementation of management and forecasting functions in the case of distributed generation, where the disperse nature of the resources requires a transparent methodology for handling information.

example of the user interface for the static wind farm data is given in Figure 3. The ensemble of static parameters consists indeed a "model" of the power system.

• On-line data of recent information from the SCADA such as total load, wind parks production, wind speed etc. A connection to the SCADA database permits the transfer of measurements to MORE-CARE.

• Numerical weather forecasts. They are introduced in the database through ftp or internet connection to a site hosting the relevant national meteorological model.

• The output of the various modules.

• Information on scheduled maintenance of the power units.

• History on the performance of each activated model. The database is developed so as to be as generic as possible. This permits easy maintenance and updating when modifications are required, e.g. addition of a new wind farm.

MORE-CARE MORE-CARE

Emphasis is given in modeling wind farms. Due to the high number of wind turbines that can be present in a system, a wind farm is modeled by clusters of wind turbines, each one composed by identical machines.

Forecasts are provided per cluster of identical turbines. This is necessary for UC and ED to produce schedules for the number of wind turbines that are going to operate during each time step of the forecast horizon without recalling information from the database on the characteristics of the wind farms. When MORE-CARE operates as a Modularity has permitted the implementation of alternative mathematic models for the same function, ranging from simple and robust to more sophisticated ones. This permitted to evaluate the on-line performance of different approaches on a specific problem, i.e. on load forecasting.

A. Data management

In order to achieve modularity of the implemented models it was necessary to develop an appropriate transparent data structure for interfacing the various modules. The choice of a relational database was made for this purpose. The modules read and write directly to that database input and output data using ODBC functionality and SQL queries. Modules are called to run by the global scheduler of the application.

The database includes all information on:

• parameters of the application like forecasting horizons, time-steps, activated models etc.

• Information on the types of available on-line data by SCADA or meteorological system.

• Static parameters for the description of the power stations and the units. E.g a detailed representation of the wind farms, the characteristics of the site, the types of wind turbines and the characteristics of each type are given. An stand-alone wind forecasting platform (i.e. in Ireland) forecasts of the total windfarm power are generally sufficient. The database permits information given at regional level (e.g. meteorological forecasts) to be related to the wind farms of that region. By this way wind upscaling forecasting functionalities can be developed -Figure 2.

B. The Man-Machine Interface

An ergonomic man-machine interface (MMI) has been developed for the communication with the operators. It communicates directly with the database and aims to assist the maintenance of the various power system and application parameters described in the previous Paragraph. The MMI gives the possibility to display on-line SCADA data in various alternative ways enabling an ergonomic monitoring of the power stations and especially the renewable plants. The MMI gives various possibilities for displaying forecasting results to the operators, i.e. in forms of confidence intervals, which are regularly updated when new forecasts are available.

III. LOAD FORECASTING

The aim of the load forecasting function is to predict the total power system load for a horizon up to 48 hours ahead. Predictions are required in the form of spot values as well as confidence intervals. A particularity of the More-Care application is that regular updates (i.e. every 20 minutes) of the forecasts are required for a continuous optimization of the power units scheduling due to the high variability of the renewables production. Frequent updates of load predictions are needed due to the sensitivity of load to external variations like weather changes. Although in an interconnected system, such changes would have only a local effect, in an autonomous system they may have a global effect. The influence of weather on power demand can be further accentuated in autonomous systems with an important tourism sector.

In the classical load-forecasting problem, forecasts are performed once or twice per day at a specific time origin (i.e. at 06:00). In that case, a forecasting model needs to update a typical daily load pattern taking into account temperature predictions, type of day etc. When load forecasts are produced with a sliding-window scheme, the forecasting task itself becomes more challenging since it corresponds to a fully dynamic problem.

Another difficulty in the case of autonomous systems concerns the quality of input. Data are often affected by events like power cuts or even blackouts and problems in data transmission. Being unpredictable, such events, when present in the data, affect the performance of the models.

A. Load preprocessing module.

The robustness of a prediction model on corrupted input varies according to the type of model. It is expected that in any case, the quality of the results be affected. A preprocessing or load correction (LC) module has been developed to restore problems in the data and prepare the input for the load forecasting functions. Data restoring covers cases of erroneous and missing data values due to SCADA and communication problems (i.e. interruption of links to power stations) and also major power cuts and blackouts. Normal variability of the load or minor power cuts are not considered. The LC-module also reconstructs historic input when this is not available (i.e. for several days); this enables the software to be operational after a period of stoppage, or when load data are missing for long periods.

The LC-module integrates several generic rules for on-line error detection. It then restores problematic data based on a number of correction rules. Thresholds for these rules are controllable from the database. This type of preprocessing was found to considerably improve the on-line performance of the load forecasting models.

B. Implemented models.

Several methods have been developed and integrated in the LF-module of More-Care for on-line operation: {LF-1} "Robust" forecasting (EMP/Armines).

A simple algorithm is implemented based on the load value from the same time of the previous day or the same day of previous week, scaled by a factor depending on the load variation between the two days or weeks. This is a robust approach that is used as a base-line one.

{LF-2} Fuzzy neural network model (EMP/Armines).

A adaptive fuzzy-neural network model with Gaussian fuzzy sets in the premises and ARX functions in the consequent part of the rules is developed. The model receives past load values to produce short and long-term forecasts. The model parameters are auto-adapted as new data arrive. By this way, the model compensates for the lack of information that affects load, but is not available on-line (i.e. temperature)

{LF-3} ARMA model (RAL).

The ARMA forecast uses a second order linear autoregressive model with one moving average parameter. Past load data is used to detrend the load time series used for each forecast by removing the weekly mean and diurnal variation before applying the ARMA model.

C. Evaluation.

The parameters of the above models have been estimated for the cases of Crete and Madeira using historical data of several years. The performance was then evaluated on an independent set of data. Figure 6 shows the performance for the case of Crete using the LF-2 model. The performance for Crete is between 1.6-4.1 % (average: 3.72%) for 1-24 steps ahead and 4.6-5.0 % (average: 4.91%) for 25-48 steps ahead. For the case of Madeira the error is 2.0-3.1% (average: 2.9%) for 1-24 hours ahead. Currently, the operation of the modules is evaluated on-line. ARMINES developed an on-line module for generating forecasts of the load at each bus-bar based on the forecasts of the total power system load. The algorithm accounts for the characteristics of each bus-bar and the recent history of load distribution in the power system. Such forecasts are used within MORE-CARE by an Optimal Power Flow function of the ED module.

IV. HYDRO FORECASTING.

The hydrological forecasting aspects of this project were to develop accurate, robust a nd consistent forecasting procedures that can be applied under general stream-flow and meteorological conditions. Unfortunately most hydrological processes take place predominantly underground and can neither be seen nor measured. In addition the physical properties of the ground vary considerably and depend in a highly non-linear and time dependent way on the prevailing environmental processes.

As such, the modeling and forecasting of hydrological processes is, in general, very difficult. Two complimentary approaches to hydrological forecasting were adopted, one being hydrograph fitting and the other time series analysiswith both traditional and novel approaches being used. The accuracy, stability and robustness of the forecasting methods was assessed through a series of comprehensive trials. Each trial involving an out-of-sample forecast of measured field data from Madeira. These trials not only established the accuracy of the current forecasting procedures but also indicated appropriate ways to improve the instrumentation, calibration and operational procedures for the measured variables to improve the forecast quality.

The program of work within MORE-CARE project had four types of output. Firstly, the forecasting methods of this project shall be available to the community so that they can quickly and accurately forecast river-flow rates. Secondly, novel non-linear time series analysis methods offer the potential of an improvement in forecasting performance. Thirdly, investigations undertaken during the project have improved our understanding of aspects of the dynamic, nonlinear, time dependent behavior between river-flow and rainfall. Finally, the investigations into the soil-moisture dependence of river-flow rate identified that the experimental design for measuring hydrological processes can be improved significantly.

The main achievements of the hydrological forecasting aspects of this project are: 1. A hydrograph forecasting methodology, based on the Nash hydrograph, has been developed and implemented in the MORE-CARE system to provide an out-of-sample forecast for two weeks into the future. The Figure 7 shows the accuracy of the Nash hydrograph forecasting approach using Madeira river-flow data. 2. A vector linear auto-regressive analysis and forecasting method has been developed and implemented in the MORE-CARE system to provide an out-of-sample forecast for twelve hours into the future. 3. In addition to the modest amount of high frequency (hourly) data from Madeira provided by the collaboration, appropriate data for the validation and verification of the hydrology forecasting methods have been obtained from UK and US sources. 4. A novel vector linear transfer function analysis and forecasting method, based on the moment hierarchy approach to time series analysis, is in the final stages of development and is soon to be implemented in the MORE-CARE system. 5. Recognizing the non-linear and time-dependent nature of the river-flow/rain-fall problem, a novel non-linear transfer function method is implemented and tested during the latter stages of this project with the findings being included in the final project report. Two theoretical studies have been undertaken, both based on the Volterra functional expansion, one on the classical representation of the non-linear response and the other on the functional dependence of the response functions to external factors (such as soil moisture).

Rabacal forecast comparative difference probability distribution V. WIND POWER FORECASTING.

The purpose of the Wind Power Forecasting (WPF) module is to provide forecasts for the total power output of each wind farm connected to the power system as well as the corresponding confidence intervals . The main state-of-the art approaches for WPF are based either on statistical or on physical modeling. Here, the statistical, or time-series approach is followed to develop short and long-term models. "Short-term" is given a slightly different meaning in wind forecasting. It is associated with models aiming to predict wind power based only on past measured data available from the SCADA system. The aim of such models is to capture temporal correlations and use them to predict the future. Such models are found to outperform persistence and often several long-term models up to 6 hours. However, if one extends to longer horizons, then history is not sufficient to predict the future. One has to consider Numerical Weather Predictions (NWP) as input to the model. The models developed in this project were configured to operate with input from two NWP systems: SKIRON f or the case of Crete and HIRLAM for Ireland. Forecasts from SKIRON were given for a grid of 15x15 km around 5 wind farms in Crete at 10 m. Forecasts from HIRLAM were given for 11 wind farms in Ireland at the wind farm site and for 4 model Levels.

Following the production of forecasts by short and longterm models, a combination of the two types of forecasts is performed to obtain an optimal performance over the whole horizon of 48 hours.

A. Implemented methods.

The following methods have been developed by the partners:

{WF-1} Persistence/Naives/Direct (EMP/ARMINES).

Persistence or Naïve predictors use the most recent value or moving averages of wind power as forecast for the entire planning horizon. "Direct" prediction is a simple robust approach that converts NWPs from 10 m to the hub height. These methods were implemented in a base-line module taking into account future availability of the machines.

{WF-2} Wind Power Forec. Module A (EMP/ARMINES).

ARMINES has developed a wind power forecasting module integrating short and long-term models.

Short-term: adaptive fuzzy neural networks (F-NN) were found to outperform other types of models like ARMA, neural networks, wavelet networks etc. A short-term module was developed based on F-NNs with wind power as input and/or wind speed or direction. F -NN models were found to outperform persistence up to 20% for the first 6 hours.

Long-term: An adaptive F -NN was developed for this task. Adaptivity permits on-line compensation for changes in the environment like addition of wind turbines etc. The artificial intelligence model permits the easy use of various types of explanatory input to optimize results. The model converts directly input to wind park power without intermediate steps. It receives both SCADA data and generic Numerical Weather Predictions as input. The model was configured to operate with the HIRLAM and SKIRON meteorological systems for the cases of Crete and Ireland respectively. The long-term module was found to outperform persistence up to 60% for 48 hours ahead. Figure 8

{WF-3} Wind Power Forecasting Module B (RAL).

The RAL Meteo software module has been designed to operate either as part of the MORE-CARE software suite or in stand-alone mode (for which purpose a user-friendly manmachine interface has been developed). Wind power forecasts are produced on both fast (e.g. 15 minute intervals up to 8 hours ahead) and slow (e.g. 1 hour intervals up to 48 or 72 hours ahead) cycles. The slow forecast is based on input meteorological wind speed forecast, downloaded over the Internet from a Met Office forecast. This is error-checked, scaled appropriately to hub-height, then converted to wind power by applying a filter initially based on the manufacturer's power curve. At the same time, wind power production data is collected from the wind farm via the wind farm's own Scada system and used to assess the errors in the wind power forecasts. Once sufficient operational data has been obtained, much improved wind power predictions can be obtained by deriving a new filter based on the historic forecast wind speed and measured wind power data. More sophisticated filters are being developed for the case where wind speed and direction are also measured at the wind farm. {WF-4} Wind Power Forecasting C (AUTh). Short-term module: AUTh tested various Artificial Neural Network Models that used as inputs the most recent values of wind speed or/and wind power. Since, power values are mostly correlated to the forecasted output, wind speed inputs may b e omitted without affecting the model efficiency. Forecasts are as good as persistent or slightly better. Since the power series are not stationary, the removal of moving average is suggested as a safe method to bypass rogue errors, shorten the training procedure and improve forecasts. Long-Term Module (up to 72 hours ahead). AUTh developed neural network models based on meteo. information i.e. from the SKIRON system. These models can also utilize various numerical values as input, eg. forecasting time lag, binary codes for special cases etc. The error of the model does not significantly grow with time lag, as persistent error does.

VI. CONCLUSIONS

The More-Care project [5] permitted to explore different approaches for load, hydro, and wind power forecasting. The most efficient ones were implemented in on-line modules and integrated into the More-Care EMS software. The software is installed and is under on-line evaluation at the islands of Crete and Madeira, while a stand-alone wind forecasting application is running in Ireland.

FORECASTING

Figure 1 :

 1 Figure 1: Architecture of the MORE-CARE system.

ForecastsFigure 2 :

 2 Figure 2:Procedure for the generation of wind power forecasts for each wind park by using different models for short and long-term forecasting.

Figure 3 :

 3 Figure 3: Example of user interface for wind farm description.

Figure 4 :

 4 Figure 4: More-Care MMI with display of load and renewables forecasts.

Figure 5 :

 5 Figure 5: MMI functions for displaying wind park forecasts.

Figure 6 :

 6 Figure 6: Performance of load forecasting.

α=100. 0 Figure 7

 07 Figure7A summary statistic for the multiple Nash peak and dynamic forecasting procedure: Integrated (over delay) probability for the Nash forecasts compared to persistence forecast values with the mean improvement being approximately 25%.

Figure 8 :

 8 Figure 8: Performance of WF-2 model for a wind farm in Ireland.

Figure 9 :

 9 Figure 9: Performance of neural networks for 5 wind farms in Crete.

SCHEDULER SCHEDULER MMI MMI DATABASE DATABASE MODULE m MODULE m ODBC SCADA

	SCADA
					SKI RON	HILARM
				Star t
					WindSpeed WindDirection	Wi ndSpeed WindDirection
	WindFarm
	Wi ndSpeed	ReadScada
	WindDirection Wi ndPower
			Onl i ne_WindFarms		ConvertSkiron	Convert	Hilarm
			WindFarm		Data	Data
			Wi ndSpeed WindDirection	UpdatePastWi ndData
			Wi ndPower
					WF_METEO
					Time
				PAST_Wi nd	LTWF_METEO	WindFarm
				Ti me_Index FarmId	T ime	Wi ndSpeed WindDirection
				Power	WindFarm
		Static_ WindFarm Static_ WindFarm	Wi ndSpeed WindDirection	Wi ndSpeed WindDirection
	WindFarm				Get LongTerm
	Availability WindFarm			LT_LinkToWF	Horizon
	Cl uster Turbines			Insert Data For Short	Forecasts
	WhenOff		S T _ LinkToWF	Term Li nked	Insert data For
	WhenOn			Windfarms	Long Term Wi ndFarms
	Calculate			PAST_Wi nd	LTWF_METEO
	Avaliability on Horizon			Ti me_Index FarmId	T ime
				Power	WindFarm
				Wi ndSpeed Wi ndDi rect i on	Wi ndSpeed
					WindDirection
	Wind Farm
	Avai l Hori zon Ti me_Index				MAIN_ Par ams MAIN_ Par ams
	WindFarm
	Cl uster TurbineNumber			WF
				module
					Static_ PowerSystemTable Static_ PowerSystemTable
		Setup _WF_M1a_Armines Setup _WF_M1a_Armines
			STWFMOD_Wind_Output Time	LTWFMOD_Wind_Output	Static_ WindFarmClusters Static_ WindFarmClusters
		Setup _WF_M1b_ Setup _WF_M1b_	Ar mi nes Ar mi nes	Farm	Time Farm
				Cl uster	Cl uster
				Turbines Speed	Tur bi nes
				Power	Speed Power	Static_ WTurbineBrands Static_ WTurbineBrands
		Setup _WF_M1c_Armines Setup _WF_M1c_Armines
				CombineMai CombineMai
				nData nData	Static_ WTurbi neCurves Static_ WTurbi neCurves
					Brand, Speed, Power Brand, Speed, Power
			STWF_Wind_Output STWF_Wind_Output	LTWF_Wind_Output LTWF_Wind_Output
				Time Time	Time Time
				Farm Farm	F a rmCluster F a rmCluster
				Cluster Turbines Cluster Turbines	Turbines Turbines
				Speed Speed	Speed Speed
				Power Power	Power Power

ACKNOWLEDGEMENT

This work has been performed in the frame of the MORE-CARE project (ERK5-CT1999-00019) supported by the European Union.

BIOGRAPHIES

George N. Kariniotakis


Cédric Demonceaux
email: cedric.demonceaux@u-picardie.fr

Pascal Vasseur
email: pascal.vasseur@u-picardie.fr

Adaptative Markov Random Fields for Omnidirectional Vision

Images obtained with catadioptric sensors contain significant deformations which prevent the direct use of classical image treatments. Thus, Markov Random Fields (MRF) whose usefulness is now obvious for projective image processing, can not be used directly on catadioptric images because of the inadequacy of the neighborhood. In this paper, we propose to define a new neighborhood for MRF by using the equivalence theorem developed for central catadioptric sensors. We show the importance of this adaptation for a motion detection application.

Introduction

The interesting properties of omnidirectional vision has led to the development of different acquisition systems such a camera networks, rotating cameras or catadioptric sensors [START_REF] Benosman | Panoramic Vision: Sensors, Theory, Applications[END_REF]. Catadioptric vision consists in associating a convex mirror with a projective camera whose the optical axis is aligned with the axis of the mirror (Fig. 1). The main advantage of these sensors is the acquisition of an omnidirectional image with a single shot. However, catadioptric image contains significant deformations due to the geometry of the mirror and to the sampling of the camera [START_REF] Daniilidis | Image processing in catadioptric planes: Spatiotemporal derivatives and optical flow computation[END_REF](Fig. 1). These deformations have important consequences for the image processing and the direct application of classical operators can not provide satisfactory results.

Thus, although, since the eighties, MRF have permitted to solve many problems in image processing such as image segmentation, image restauration, motion detection and so on, they become unsuitable if they are applied in the same way for catadioptric images. Indeed, MRF use a local dependence between pixels from a predefined neighborhood which has to be reconsidered because of the deformations of the catadioptric image. The aim of this paper is then to propose an adaptation of the MRF for the catadioptric images from the definition of a new neighborhood system which takes into account the deformations. We define this neighborhood from the equivalence between the catadioptric image formation and the projection of the 3D points on the image plane via a unitary sphere [START_REF] Geyer | Catadioptric projective geometry[END_REF].

Markov Random Fields in Image Processing

Spatial MRF have been introduced by Besag [START_REF] Besag | Spatial interaction and the statistical analysis of lattice systems[END_REF] and popularized in image processing by Geman and Geman [START_REF] Geman | Stochastic relaxation, gibbs distributions, and the bayesian restoration of images[END_REF]. These models are defined from a graph G = (S, U) where S = {s 1 , s 2 , • • • , s N } is the set of vertices (named set of sites in the case of MRF) and U the set of edges. In this graph G, two sites s i and s j are neighbors if there exists an edge u ij ∈ U between them. Thus, the set of neighbors of a site s is noted V s . From V = {V s |s ∈ S}, we define the set C whose elements

= {X 1 , • • • , X N } of realization ω = {ω 1 , • • • , ω N } Definition 1 X is a MRF according to V if : (i) ∀ω P (X = ω) > 0 (ii) ∀(s, ω) P (X s = ω s |X r = ω r , r = s) = P (X s = ω s |X r = ω r , r ∈ V s ) (1)
A Markov field is a random field where each site is only influenced by its neighboring sites. Using the characterization (1) does not permit to obtain a simple expression of joint distribution P (X = ω). However, this distribution can be easily obtained from the Hammersley-Clifford theorem [START_REF] Besag | Spatial interaction and the statistical analysis of lattice systems[END_REF] which establishes the equivalence between the Markov fields and the Gibbs fields.

Problems such as image restauration, classification and segmentation can be schematized as follows: from observations Y (known information) in the image or in the sequence (grey level, color, edges, displacements, . . .), we want to find the set of labels X (hidden information) which could originate these observations. The solution of these problems requires to clarify the relation between the labels and the sites. However, these models are not sufficient to solve these inverse problems. Indeed, because of the loss of information during the image formation and because of the noise in the observations, these inverse problems are generally ill-posed according to Hadamard. We have to make an additional hypothesis on the researched set of labels. This is why we add a constraint named a priori knowledge which forces the set of labels to be Markovian.

The bayesian estimation permits to modelize jointly the labels and the observations thanks to the coupled random field (X, Y ) with the probability distribution P XY . For example, if we consider the bayesian estimator Maximum a Posteriori (MAP), the researched markovian field of labels x is solution of

arg min x∈Ω U (x, y) = arg min x∈Ω (U 1 (x, y) + U 2 (x)) (2)
where U 1 is the term linked to the data which expresses the relation between observations y and labels x and where U 2 describes the markovian a priori of the label field and is decomposed into a sum of local potentials which depends on the neighboring system:

U 2 (x) = <s,s ′ >∈C V 2 (s, s ′ ).
(

) (a) (b) (c) (d) 3

Central Catadioptric MRF

Catadioptric Image Formation and Equivalence Theorem

Baker and Nayar classified catadioptric sensors into two categories according to the number of viewpoint [START_REF] Baker | A theory of single-viewpoint catadioptric image formation[END_REF]. Thus, sensors with a single viewpoint, named central catadioptric sensors, permit a geometrically correct reconstruction of the perspective image from the original catadioptric image. This category is made of parabolic mirror associated to orthographic camera and hyperbolic, elliptic and plane mirrors with perspective camera. The second category with several viewpoints has geometric properties very less significant and is made of the other possibilities of association between mirrors and cameras. In this paper, we are only interested in central sensors for which Geyer and Daniilidis have demonstrated the equivalence with a two-step projection via a unitary sphere centered on the focus of the mirror (the single viewpoint) [START_REF] Geyer | Catadioptric projective geometry[END_REF]. This two-step projection consists first in projecting a 3D point to the sphere from the center of the sphere. The next step consists in projecting the point on the sphere to the image plane from a point placed on the optical axis (Fig 4). The positions of this point of projection and of the image plane are determined with the calibration of the sensor. This equivalence is very interesting since it allows to perform image processing in a new space in which the deformations are taken into account. Many authors used this space in order to compute the derivative of catadioptric image [START_REF] Daniilidis | Image processing in catadioptric planes: Spatiotemporal derivatives and optical flow computation[END_REF], the pose of the sensor [START_REF] Paulino | Pose estimation for central catadioptric systems: An analytical approach[END_REF] or to detect catadioptric lines [START_REF] Barreto | Direct least square fitting of paracatadioptric line images[END_REF]. We also propose to use this space in order to determine the most appropriate neighboring system for the central catadioptric MRF.

Central Catadioptric Neighboring System

Let P the projection for which any point P i in the catadioptric image associates its equivalent point P s in spherical coordinates on the unitary sphere S 2 (Fig 4 ),

P : S → S 2 s → (1, θ s , φ s ) (4)
Then, we define the new neighborhood as follows:

∀s ∈ S, P(s) = (1, θ s , φ s )

V s =    t ∈ S, t = s P(t) = (1, θ t , φ t ) such as |θ t -θ s | ≤ 1 N and (|φ t -φ s | ≤ 1 M or |φ t -φ s | ≥ 2π -1 M )    , (5)

Experimentation

In this part, we present results obtained from a markovian modelling adapted to catadioptric images in the case of motion detection. Consider two consecutive images I(t) and I(t + 1) from a sequence acquired from a fixed catadioptric sensor. The problem consists in detecting the moving parts in the image between the instants t and t + 1. For this, we pose E = {e s , s ∈ S} the label field where e s have values in {0, 1} (0=motionless site, 1 = mobile site) and O = {o s , s ∈ S} is the observation field. We consider the following markovian energy: We have compared our markovian modelling adapted to the central catadioptric images with a modelling obtained from a classical neighborhood of first order (Fig 2(a)). Constants N and M are chosen in order to obtain a catadioptric neighborhood similar than the first order neighborhood in the center of the image. We set α = 3 (threshold from which we consider the pixels in motion). Between instants t and t + 1, three persons are moving in front of the camera. Maximum of likelihood given in Fig. [START_REF] Geman | Stochastic relaxation, gibbs distributions, and the bayesian restoration of images[END_REF] is minimized by the deterministic algorithm I.C.M. [START_REF] Geman | Stochastic relaxation, gibbs distributions, and the bayesian restoration of images[END_REF] with two values of λ. If λ = 0.6, the dependence between the neighbors is not sufficient in order to remove the noise. With λ = 1.2, the catadioptric neighborhood allows to remove all the false positive in contrary to the classical neighborhood. In the latter case, we can also note that a block effect appears in the results.

U (e, o) = s∈S f (e s , o s ) + λ <s,s ′ >∈C (1 -δ(e s , e ′ s )) (6)

Conclusion and perspectives

In this paper, we have proposed an adaptation of the markovian modelling to the central catadioptric images. Experimentations show very encouraging results. In the same way than Heitz et al [START_REF] Heitz | Multiscale minimization of global energy functions in some visual recovery problems[END_REF] for the classical case, it should be very interesting to develop a multi-scaling method in order to reduce the sensitivity of the initialization of the deterministic algorithms and to speed-up the treatment.

Figure 1 .

 1 Figure 1. Omnidirectional sensor and omnidirectional image.

c

  (Fig 2), named cliques, are the parts of S which verify    ∃s ∈ S such as c = {s} or ∀{s, t} ⊂ c, t ∈ V s Consider a random field X

Figure 2 .

 2 Figure 2. Neighboring system and associated cliques. Neighborhood (a) and cliques (c) of first order, Neighborhood (b) and cliques (d) of second order.

Figure 3 .

 3 Figure 3. d(P 1, P 2) = d(P 3, P 4) when P 1 does not have physically the same influence on P 2 than P 3 on P 4.

Figure 4 .

 4 Figure 4. Equivalence between the catadioptric projection and the two-step mapping via the sphere.

Figure 5 .

 5 Figure 5. (a) Neighborhood at the periphery of the image, (b) Neighborhood in the center of the image.

Figure 6 .

 6 Figure 6. (a) Image at instant t, (b) Maximum of likelihood, (c) classical markov λ = 0.6,(d) catadioptric markov λ = 0.6, (e) classical markov λ = 1.2,(d) catadioptric markov λ = 1.2.

Acknowledgment

This work is supported by Agence Nationale de la Recherche -projet CAVIAR (ANR-05-BLAN-0291-01).


Jean Charles Bazin

Cédric Demonceaux

Pascal Vasseur
email: pascal.vasseur@u-picardie.fr

Pascal Vasseur Fast

Jean Charles Bazin
email: jeancharles.bazin@gmail.comcedric.demonceaux

Fast Central Catadioptric Line Extraction

Introduction

Catadioptric vision sensors (associations of a camera with a mirror) are now broadly used in many applications such as robot navigation, 3D scene reconstruction or video surveillance [START_REF] Benosman | Panoramic Vision: Sensors, Theory, Applications[END_REF]. Their large field of view is indubitably the major reason of this success. Baker and Nayer classified these sensors in two respective categories [START_REF] Baker | A theory of single-viewpoint catadioptric image formation[END_REF]. First, sensors with a single viewpoint, named central catadioptric sensors are made of parabolic mirror associated to orthographic camera and hyperbolic, elliptic and plane mirrors with perspective camera. The second category with different viewpoints has geometric properties less significant and is made of the other possibilities of association between mirrors and cameras. In this paper, we are only interested in central sensors which permit a geometrically correct reconstruction of the perspective image from the original catadioptric image. However, their employment presents some drawbacks because of the deformations induced by the mirror. For example, some very useful classical treatments in perspective image processing can be no more performed on catadioptric images because they are inadequate. One of these major treatments deals with line extraction. Thus, while in the perspective case line detection is perfectly known and efficiently solved, with catadioptric images the problem is absolutely not trivial. Indeed, the projection of any 3D real line is a conic eventually degenerate. Thus, in the case of an uncalibrated sensor, it is necessary to estimate five parameters for each line while only two parameters are sufficient for a calibrated sensor. If we consider the projection of a 3D point by the way of the unitary sphere (fig. 1(a)) as proposed in [START_REF] Barreto | Geometric properties of central catadioptric line images[END_REF] [4] [START_REF] Geyer | Catadioptric projective geometry[END_REF] with the formalism defined in [START_REF] Barreto | Geometric properties of central catadioptric line images[END_REF] [4], we can define oriented projective ray P 1 passing by 3D point x w and the center of the sphere. This ray intersects the surface of the sphere in x s . We then consider oriented projective ray P 2 passing by x s and a point situated on the z-axis between the center of the sphere and the north pole. This point is at distance ξ from the center of the sphere and depends only on the mirror geometric characteristics. P 2 intersects plane at infinity in point x i . Finally, homography H defined between the plane at infinity and the catadioptric image plane projects point x i into point x c . H includes intrinsic parameters of the camera, possible rotations between the sphere frame and the camera frame, and finally the parameters of the mirror. According to this model, we can develop the projection of a 3D line into the catadioptric image plane (fig. 1(b)). We consider plane Π R which contains the real 3D line and the center of the sphere. This plane intersects the sphere and then defines a great circle onto its surface. The set of oriented projective rays passing by the points of the great circle and point O 2 define then a cone which intersects plane at infinity into conic C i . Finally, homography H transforms C i into conic C c in the catadioptric image plane. In plane at infinity, we know that the equation of conic C i is equal to :

Ci = 2 4 n 2 x (1 -ξ 2 ) -n 2 z ξ 2 nxny(1 -ξ 2 ) nxnz nxny(1 -ξ 2 ) n 2 y (1 -ξ 2 ) -n 2 z ξ 2 nynz nxnz nynz n 2 z 3 5 (1)
with (n x , n y , n z ) T the vector which describes the normal to plane P w which contains the 3D line. We obtain the equation of conic C c in image plane thanks to the following relation :

C c = H -T C i H -1 (2)
Finally, a pixel x c = (u v 1) T belongs to the conic C c if the equality x T c C c x c = 0 is verified.

In this paper, we propose a new method for calibrated catadioptric line detection which permits a very fast, robust and accurate detection. The proposed approach consists in roughly estimating the possible catadioptric lines in the image and in verifying if each possible line is a real catadioptric line. The rest of the paper is organized as follows. Section II is devoted to the related works which deal with catadioptric line detection in calibrated and uncalibrated cases. In section III, we present a complete description of the algorithm. Section IV is devoted to experimental results with quantitative and qualitative evaluations. We finally conclude in section V on different perspectives.

Related Works

The methods of catadioptric line detection and estimation can be divided in three categories. The first class deals with methods applicable as well in the calibrated case as in the uncalibrated case and includes the algorithms of conic fitting In the rest of this section, we only develop the two first categories because the third is not enough general and concerns only paracatadioptric cameras. Conic fitting algorithms determine the curve that best fits the data points according to a certain distance metric [START_REF] Zhang | Parameter estimation techniques: a tutorial with application to conic fitting[END_REF]. In [START_REF] Barreto | Fitting conics to paracatadioptric projections of lines[END_REF], the authors present a comparison of the normal least squares (LMS), approximate mean square (AMS), Fitzgibbon and Fisher (FF) [START_REF] Fitzgibbon | A buyer's guide to conic fitting[END_REF], gradient weighted least square fitting (GRAD) and orthogonal distances (ORTHO) methods for the specific problem of paracatadioptric line detection. Their conclusions are that GRAD and ORTHO are the most robust to noise and that all methods perform poorly when the amplitude of the occlusion is above 240 • . Since most of the catadioptric lines have an amplitude less than 45 • , it appears clearly that these methods are unsuitable for general central catadioptric line detection and estimation. Moreover, these methods suppose that the pixels from the edges have been already classified into chains representing the different possible catadioptric lines.

In the calibrated case, homography H and parameter ξ are known. In this way, a 3D line is determined thanks to a vector (n x , n y , n z ) T . This vector represents the normal of the great circle on the unitary sphere obtained by the intersection of the plane which contains the center of the sphere and the 3D real line (fig. 1(b)). A 3D real line can be also represented by two angles φ and θ which respectively are the elevation angle and the azimuth angle of the normal vector. Each catadioptric line is then represented by only two parameters and a simple adaptation of the Hough transform can solved the problem. This is this kind of approach which is proposed in [START_REF] Ying | Catadioptric line features detection using hough transform[END_REF] and [START_REF] Vasseur | Central catadioptric line detection[END_REF]. The mean difference between these methods deals with the space in which the treatments are performed. In [START_REF] Vasseur | Central catadioptric line detection[END_REF], the image is projected on the unitary sphere and the 3D coordinates of the pixels are then used while in [START_REF] Ying | Catadioptric line features detection using hough transform[END_REF], they apply the algorithm directly in the image. Although these two approaches present interesting results, it is worth noting that they present the classic defects of the Hough transform such as the best sampling step for φ and θ for example. In order to avoid these drawbacks, we can note that if two pixels of a catadioptric line are known, it is then possible to compute the normal of the great circle and then to obtain the corresponding values of φ and θ. In [START_REF] Mei | Fast central catadioptric line extraction, estimation, tracking and structure from motion[END_REF], the authors propose a randomized Hough transform which selects randomly two points in the image of edges in order to compute the φ and θ angles. These angles are then used in an accumulator for the detection of the most confident catadioptric lines.

Central Catadioptric Line Detection Algorithm

Our line detection algorithm for central catadioptric sensor consists first in applying a Canny edge detector (Fig 4(b)). Then, we proceed to an edge chaining which consists in extracting connected pixels from edges in order to form lists with a length superior or equal to a threshold (NbPixels) (Fig 4(c)). To detect the lines in the scene consists then in verifying if these chains are the projections of 3D lines. In this way, we apply a split and merge algorithm of the chains. First, an adaptation of the polygonal approximation of the classical perspective case is proposed in order to find which chains or parts of chains are catadioptric projections of lines. This process is performed thanks to a division criterion which cuts the chains at a particular position if the chain is not verified as a catadioptric line. Next, we use a fusion criterion in order to group the different chains in the image which represents the same central catadioptric lines. These both criteria are discussed in the following of the paper.

Division Criterion

Consider the two endpoints of a chain of N pixels with coordinates P 1 = (X 1 , Y 1 , Z 1 ) and P 2 = (X 2 , Y 2 , Z 2 ) on the unitary sphere S 2 . These points define a single central catadioptric line in the image and then a great circle C on the sphere (cf fig( 1(a)(b)). This circle results from the intersection of the unitary sphere and a plane which contains the sphere origin 0 1 and whose a normal vector is

-→ n = ---→ O 1 P 1 × ---→ O 1 P 2 = (n x , n y , n z ) T .
Then, the equation of C is :

n x X + n y Y + n z Z = 0 (X, Y, Z) ∈ S 2
We consider that a point on the sphere with coordinates (X i , Y i , Z i ) of the chain belongs to the great circle if the distance between this point and the plane defined by the great circle is less than a threshold:

|n x X i + n y Y i + n z Z i | ≤ DivThreshold.
This chain is then considered as a central catadioptric line if at least 95% of its points belong to the great circle.

In the opposite case, we cut the chain into two sub-chains at the point (X j , Y j , Z j ) which maximizes the following error

||(X i , Y i , Z i ). -→ n ||, i = 1 • • • n (the furthest point from the plane).
This division step stops when the chain is considered as a central catadioptric line or when the length of the sub-chains is less than the threshold (NbPixels). At the end of this step, we then obtain the whole set of central catadioptric lines in the image. However this method may generate a multi-detection of the same lines. In order to compensate this drawback, we then propose to merge the similar catadioptric lines.

Fusion Criterion

Let define two catadioptric lines d 1 and d 2 detected with the previous method. These lines respectively characterized by -→ n 1 and -→ n 2 define two planes in the 3D space passing through the origin of the unitary sphere,

Π 1 = {U = (X, Y, Z) ∈ R 3 , -→ n 1 .U = 0} and Π 2 = {U = (X, Y, Z) ∈ R 3 , -→ n 2 .U = 0}
. We consider that these detected catadioptric lines are similar if they define the same 3D plane, that is to say if :

1 -| -→ n 1 . -→ n 2 | ≤ FusThreshold.
In this case, the two catadioptric lines are merged into a single line. The catadioptric line equation is then updated from the pixels of the chains which belong to d 1 and d 2 as follows. Let note respectively

M 1 = (X 1 i , Y 1 i , Z 1 i ) i=1•••N1 and M 2 = (X 2 i , Y 2 i , Z 2 i ) i=1•••N2 , the pixels of catadioptric line d 1 (resp. d 2 ). Let M , the matrix of dimension (N 1 + N 2 ) × 3, M = 0 B B B B B B B B @ X 1 1 Y 1 1 Z 1 1 . . . . . . . . . X 1 N 1 Y 1 N 1 Z 1 N 1 X 2 1 Y 2 1 Z 2 1 . . . . . . . . . X 2 N 2 Y 2 N 2 Z 2 N 2 1 C C C C C C C C A ,
The normal vector -→ n = (n x , n y , n z ) T of the great circle associated to the catadioptric line is then solution of :

M. -→ n = (0, • • • , 0) T . (3)
The solution of (3) is obtained from the SVD of matrix M [START_REF] Jennings | Matrix Computation[END_REF].

Experimentations

We have tested our central catadioptric line detector on different kinds of omnidirectional images. We first propose some results with synthesis images for which we perfectly know the line equations in order to show the accuracy of the approach. Then, some results on real images are also proposed. In the whole set of experimentations except in one indicated case, the different thresholds are fixed as follows : NbPixels =100, DivThreshold = 0.0005, FusThreshold = 1 • .

Synthesis Images

We have generated two synthesis images for which we perfectly know the calibration parameters and line equations.

Real Catadioptric Images

We present here result for a real catadioptric image. In this case, sensor has been calibrated with the method described in

Conclusion

In this paper, we deal with the problem of line detection in central catadioptric images. Our method is valid for calibrated sensor and is comparable to the polygonal approximation algorithm. Indeed, it consists in looking for pixels in chains of edges which correspond to catadioptric lines thanks to an analytic approach contrary to previous methods based on Hough transform which depends on the sampling of the search space. Moreover, we then obtain a very fast algorithm which could be implemented for real time applications.

Fig. 1 .

 1 Fig. 1. (a) Image formation model. Example of projection via the unitary sphere for a 3D point. (b) Projection of a 3D line via the unitary sphere into the catadioptric image plane.

[ 6 ]

 6 . The second category concerns calibrated sensors and most of the proposed techniques in this category are based on adaptation of Hough transform[START_REF] Vasseur | Central catadioptric line detection[END_REF] [8][START_REF] Mei | Fast central catadioptric line extraction, estimation, tracking and structure from motion[END_REF]. Methods for uncalibrated sensors form the third category. These methods use particular geometric constraints of catadioptric sensors and are generally dedicated to paracatadioptric sensors[START_REF] Barreto | Fitting conics to paracatadioptric projections of lines[END_REF] [START_REF] Vandeportaele | A fast detector of line images acquired by an uncalibrated paracatadioptric camera[END_REF].

Fig. 2 .

 2 Fig. 2. (a)Original image, (b) Detected catadioptric lines.

Fig. 3 .

 3 Fig. 3. (a)Original image, (b) Red catadioptric lines correspond to detected catadioptric lines, (c) False catadioptric lines whan the length NbPixel is too low.

  [START_REF] Barreto | General Central Projection Systems: Modeling, Calibration and Visual Servoing[END_REF]. This image (fig4(a)) is a paracatadioptric image issued from the calibration toolbox proposed by Barreto[START_REF] Barreto | General Central Projection Systems: Modeling, Calibration and Visual Servoing[END_REF]. In figure4(b), we present the result of Canny edge detector and consecutively the detected chains of pixels extracted for the catadioptric line verification (fig 4(c)).

Figure 4 (

 4 d) shows the catadioptric line detection before the fusion step while figure 4(e) presents the final result after the fusion step. In figure4(f), we propose a more detailed view of a part of the image in order to show the accuracy of the results. Finally, from a computational time point of view, the method takes near 3 seconds with Matlab. A real time implementation constitutes the next perspective of this work.

Fig. 4 .

 4 Fig. 4. (a)Original image, (b) Canny edge detector result, (c) Extracted chains, (d) Catadioptric line detection results after division step, (e) Catadioptric line detection results after fusion step, (f) Detailed view of final results.


Philippe Martin

Lionel Rosier

Pierre Rouchon

Controllability of the 1D Schrödinger equation using flatness

Keywords: Partial differential equations, Schrödinger equation, boundary control, exact controllability, motion planning, flatness

Introduction

The exact controllability of the linear Schrödinger equation (or of the related plate equation) is investigated in [START_REF] Lions | Contrôlabilité exacte, perturbations et stabilisation de systèmes distribués[END_REF][START_REF] Machtyngier | Exact controllability for the Schrödinger equation[END_REF][START_REF] Komornik | Exact Controllability and Stabilization: The Multiplier Method[END_REF] with the multiplier method, in [START_REF] Haraux | Séries lacunaires et contrôle semi-interne des vibrations d'une plaque rectangulaire[END_REF][START_REF] Jaffard | Contrôle interne exact des vibrations d'une plaque rectangulaire[END_REF][START_REF] Komornik | Fourier series in control theory[END_REF] with nonharmonic Fourier analysis, in [START_REF] Lebeau | Contrôle de l'équation de Schrödinger[END_REF] with microlocal analysis, and in [4,[START_REF] Liu | Locally distributed control and damping for the conservative systems[END_REF] with frequency domain tests. It is extended to the semilinear Schrödinger equation in [START_REF] Rosier | Local exact controllability and stabilizability of the nonlinear Schrödinger equation on a bounded interval[END_REF][START_REF] Rosier | Exact boundary controllability of the nonlinear Schrödinger equation[END_REF][START_REF] Rosier | Control and stabilization of the nonlinear Schrödinger equation on rectangles[END_REF] and [START_REF] Laurent | Global controllability and stabilization for the nonlinear Schrödinger equation on an interval[END_REF][START_REF] Laurent | Global controllability and stabilization for the nonlinear Schrödinger equation on some compact manifolds of dimension 3[END_REF], by means of Strichartz estimates and Bourgain analysis. All those results are indirect, as they rely on observability inequalities for the adjoint system. A direct approach not involving the adjoint system is proposed in [START_REF] Littman | Exact boundary controllability of a hybrid system of elasticity[END_REF][START_REF] Rosier | A fundamental solution supported in a strip for a dispersive equation[END_REF][START_REF] Littman | The heat and Schrödinger equations: boundary control with one shot[END_REF]; in particular, [START_REF] Rosier | A fundamental solution supported in a strip for a dispersive equation[END_REF] uses a fundamental solution of the Schrödinger equation with compact support in time, and provides controls with Gevrey regularity.

In this paper, we derive in a direct way the null (hence exact, since the equation is time-reversible) controllability of the Schrödinger equation with Dirichlet control, iθ t + θ xx = 0, (t, x) ∈ (0, T ) × (0, 1), (1) θ(t, 0) = 0, t ∈ (0, T ), (2) θ(t, 1) = u(t), t ∈ (0, T ), (3) θ(0, x) = θ 0 (x), x ∈ (0, 1).

(

) 4
Email addresses: philippe.martin@mines-paristech.fr (Philippe Martin), lionel.rosier@mines-paristech.fr (Lionel Rosier), pierre.rouchon@mines-paristech.fr (Pierre Rouchon).

The initial condition θ 0 , the state θ, and the control u are complex-valued functions. Given any final time T > 0 and any initial condition θ 0 ∈ L 2 (0, 1), we explicitly build a regular control such that the state reached at T is exactly zero.

We use the so-called flatness approach, which consists in parametrizing the solution θ and the control u by the derivatives of a "flat output" y; this notion was initially introduced for finite-dimensional systems [START_REF] Fliess | Flatness and defect of non-linear systems: introductory theory and examples[END_REF], and later extended to partial differential equations, see e.g. [START_REF] Laroche | Motion planning for the heat equation[END_REF][START_REF] Petit | Dynamics and solutions to some control problems for water-tank systems[END_REF][START_REF] Petit | Flatness of heavy chain systems[END_REF][START_REF] Lynch | Flatness-based boundary control of a class of quasilinear parabolic distributed parameter systems[END_REF][START_REF] Woittennek | Motion planning for a class of boundary controlled linear hyperbolic PDE's involving finite distributed delays[END_REF][START_REF] Meurer | Flatness-based tracking control of a piezoactuated Euler-Bernoulli beam with non-collocated output feedback: theory and experiments[END_REF][START_REF] Meurer | Control of higher-dimensional PDEs[END_REF]. A similar flatness-based approach is used in [START_REF] Martin | Null controllability of the heat equation using flatness[END_REF] for explicitly establishing the null controllability of the heat equation, and generalized to 1D parabolic equations in [START_REF] Martin | Null controllability of one-dimensional parabolic equations by the flatness approach[END_REF], with a control comprising two phases: a first "regularization" phase, taking the irregular (namely, L 2 (0, 1)) initial state to a much more regular (namely, Gevrey of order 1 2 ) intermediate state; and a second phase using the flat parametrization to transfer this intermediate state to zero; a similar approach is used in [START_REF] Moyano | Flatness for a strongly degenerate 1-d parabolic equation[END_REF] in the case of a strongly degenerate parabolic equation. But whereas a zero control is sufficient in the first phase thanks to the natural smoothing effect of parabolic equations, the picture is very different for the Schrödinger equation, where no such smoothing effect exists with a zero control. Nevertheless, it is still possible to take the initial state to an intermediate state which is Gevrey of order 1 2 , but with a well-chosen nonzero control; the idea, inspired by [START_REF] Rosier | A fundamental solution supported in a strip for a dispersive equation[END_REF][START_REF] Littman | The heat and Schrödinger equations: boundary control with one shot[END_REF][START_REF] Rosier | Exact boundary controllability of the nonlinear Schrödinger equation[END_REF] is to interpret the problem as a Cauchy problem on R, where some smoothing effect does take place; notice our result improves [START_REF] Littman | The heat and Schrödinger equations: boundary control with one shot[END_REF], with a more direct proof: the trajectory in the first phase is Gevrey of order 1 in time and 1 2 in space, compared to order 2 in time and 1 in space for [START_REF] Littman | The heat and Schrödinger equations: boundary control with one shot[END_REF]. The second phase is then similar in spirit to the flatness-based control of [START_REF] Martin | Null controllability of the heat equation using flatness[END_REF][START_REF] Martin | Null controllability of one-dimensional parabolic equations by the flatness approach[END_REF][START_REF] Moyano | Flatness for a strongly degenerate 1-d parabolic equation[END_REF]; with respect to [START_REF] Littman | The heat and Schrödinger equations: boundary control with one shot[END_REF], which solves an ill-posed problem in an abstract way, our control is explicitly given as a series.

Another contribution of this paper with respect to [START_REF] Martin | Null controllability of the heat equation using flatness[END_REF][START_REF] Martin | Null controllability of one-dimensional parabolic equations by the flatness approach[END_REF] is a cleaner construction of the control, which emphasizes that the key point for establishing the null controllability of a flat partial differential equation is to find a trajectory which connects the initial condition to a sufficiently regular intermediate state; the control in the two phases is then automatically deduced from this trajectory. This construction also eliminates the constraint on the intermediate time that was needed in the preliminary version of this paper [START_REF] Martin | Controllability of the 1D Schrödinger equation by the flatness approach[END_REF].

The paper is organized as follows. In Section 2, we recall some preliminary notions. In Section 3, we precisely state what flatness means for the considered equation. In Section 4, we explicitly derive a control achieving null controllability. In Section 5, we detail how this control can be numerically implemented. Finally, we illustrate in Section 6 the effectiveness of the approach on a numerical example.

Preliminaries

We collect here a few definitions and facts that will be used throughout the paper.

We say y ∈ C ∞ ([0, T ]) is Gevrey of order s ≥ 0 on [0, T ] if there exist positive constants M, R such that

y (p) (t) ≤ M p! s R p , ∀t ∈ [0, T ], ∀p ≥ 0. More generally, if K ⊂ R 2 is a compact set and (t, x) → y(t, x) is a function of class C ∞ on K (i.e. y is the restric- tion to K of a function of class C ∞ on some open neigh- borhood Ω of K), we say y is Gevrey of order (s 1 , s 2 ) on K if there exist positive constants M, R 1 , R 2 such that for all (t, x) ∈ K and (p 1 , p 2 ) ∈ N 2 , |∂ p1 t ∂ p2 x y(t, x)| ≤ M p 1 ! s1 p 2 ! s2 R p1 1 R p2 2 .
By definition, a Gevrey function of order s is also of order r for r ≥ s. Gevrey functions of order 1 are analytic (entire if s < 1). Gevrey functions of order s > 1 may have a divergent Taylor expansion; the larger s, the "more divergent" the Taylor expansion. Important properties of analytic functions generalize to Gevrey functions of order s > 1: the scaling, addition, multiplication, inverse and derivation of Gevrey functions of order s > 1 is of order s, see e.g. [START_REF] Yamanaka | A new higher order chain rule and Gevrey class[END_REF]. But contrary to analytic functions, Gevrey functions of order s > 1 may be constant on an open set without being constant everywhere. For example the "step function"

φ s (ρ) :=          1 if ρ ≤ 0 0 if ρ ≥ 1 e -M (1-ρ) σ e -M ρ σ + e -M (1-ρ) σ if ρ ∈ (0, 1), (5)
where M > 0 and σ := (s -1) -1 , is Gevrey of order s on R; notice φ (i)

s (0) = φ (i) s (1) = 0 for all i ≥ 1.
The result stems from the fact that ρ → 1 R (ρ)e -1 ρ σ is Gevrey of order s, which is classically proved thanks to the Cauchy integral formula, see [37, chapter 3.11] for particular cases and [START_REF] Lynch | Flatness-based boundary control of a class of quasilinear parabolic distributed parameter systems[END_REF] for the general case.

From the two well-known properties of the Gamma function, see e.g. [1, section 6.1],

Γ(ξ)Γ ξ + 1 2 = 2 1-2ξ √ π Γ(2ξ) Γ(ξ + a) Γ(ξ + b) ∼ ξ→+∞ ξ a-b ,
we deduce at once

Γ(2ξ + 1) ∼ ξ→+∞ 2 2ξ √ πξ Γ 2 (ξ + 1). (6)
Recall that Γ(p + 1) = p!, when p ∈ N. Following the usual notation in the literature on asymptotics,

f (ξ) ∼ ξ→+∞ g(ξ) ⇔ lim ξ→+∞ f (ξ) g(ξ) = 1.
We will also use the obvious inequality

(p + q)! p!q! ≤ 2 p+q . ( 7
)
3 Flatness of the Schrödinger equation

In this section we show the system (1)-( 4) is flat, with y := θ x (•, 0) as a flat output. Loosely speaking, it means there is a one-to-one correspondence between sufficiently regular solutions of (1)-( 4) and sufficiently regular arbitrary functions y. This property will be paramount to prove controllability in the next section.

We first notice that, given any C ∞ function Y , the formal series

θ(t, x) = j≥0 x 2j+1 (2j + 1)! (-i) j Y (j) (t) (8)
u(t) = j≥0 (-i) j (2j + 1)! Y (j) (t) (9)
obviously formally satisfy (1)-(3) and

θ x (•, 0) = Y . (10)
The following property gives an actual meaning to these series when Y is regular enough.

Proposition 1 Let s ∈ [0, 2), -∞ < t 1 < t 2 < ∞, and assume Y is Gevrey of order s on [t 1 , t 2 ]. Then θ defined by (8) is Gevrey of order (s, s 2 ) on [t 1 , t 2 ] × [0, 1]; u defined by (9) is Gevrey of order s on [t 1 , t 2 ].
PROOF. We want to prove the formal series

∂ k t ∂ l x θ(t, x) = 2j+1≥l x 2j+1-l (2j + 1 -l)! (-i) j Y (j+k) (t) (11) is uniformly convergent on [t 1 , t 2 ] × [0, 1]
, with an estimate of its sum of the form

∂ k t ∂ l x θ(t, x) ≤ C k! s R k 1 l! s 2 R l 2 .
By assumption, there are constants M, R > 0 such that 7) twice and ( 6) for ξ := j. As 2j+1≥l A j,l < ∞ (ratio test), the series in [START_REF] Laurent | Global controllability and stabilization for the nonlinear Schrödinger equation on an interval[END_REF]

Y (j) (t) ≤ M j! s R j for all j ≥ 0 and all t ∈ [t 1 , t 2 ]. Hence, for all (t, x) ∈ [t 1 , t 2 ] × [0, 1], x 2j+1-l Y (j+k) (t) (2j + 1 -l)! ≤ M R j+k (j + k)! s (2j + 1 -l)! ≤ M R j+k (2 j+k j!k!) s (2j + 1 -l)! ≤ M k! s used (
is uniformly convergent for all k, l ≥ 0, implying θ ∈ C ∞ ([t 1 , t 2 ] × [0, 1]). Moreover, 2j+1≥l A j,l = (2π) s 4 √ R 1 √ R 1 l p≥0 1 √ R 1 p (p + l + 1) s 4 p! 1-s 2 ≤ C(R 1 ) 1 √ R 1 l for C(R 1 ) large enough. Hence, ∂ k t ∂ l x θ(t, x) ≤ M C(R 1 ) k! s R k 1 l! s 2 √ R 1 l . which proves θ is Gevrey of order (s, s 2 ) on [t 1 , t 2 ]×[0, 1]. As a simple consequence, u = θ(•, 1) is Gevrey of order s on [t 1 , t 2 ]. 2
We can now give a precise meaning to (1)-( 4) being flat with θ x (•, 0) as a flat output:

• obviously, any solution θ of ( 1)-( 2) that is Gevrey of order

(s, s 2 ) on [0, T ] × [0, 1] -which implies u is Gevrey of order s on [0, T ] and θ 0 is Gevrey of order s 2 on [0, 1] -uniquely determines the flat output θ x (•, 0) as a function that is Gevrey of order s on [0, T ] • conversely, any function Y Gevrey or order s on [0, T ]
gives rise by ( 8) to a solution θ of ( 1)-( 2)-( 10), which is Gevrey or order (s, s 2 ) on [0, T ] × [0, 1]; moreover, this solution is unique by the Holmgren theorem.

Controllability of the Schrödinger equation

We derive in the following Theorem 3 an explicit control steering the system from any initial state θ 0 ∈ L2 (0, 1) at time 0 to the final state 0 at time T . This proves null controllability, hence exact controllability since the Schrödinger equation is reversible with respect to time. This is a two-step procedure:

• we first apply a "smoothing" control till some arbitrary intermediate time τ < T so as to reach a "smooth" intermediate state θ(τ, •) • we then apply a flatness-based control using the parametrization [START_REF] Komornik | Fourier series in control theory[END_REF].

Notice that the initial-boundary value problem for the Schrödinger equation ( 1)-( 4) with zero control does not smooth an L 2 initial condition into a state that is Gevrey of order s Proposition 2 Let θ 0 ∈ L 2 (0, 1), and define its odd extension

θ odd 0 ∈ L 1 (R) ∩ L 2 (R) by θ odd 0 (x) :=    θ 0 (x) if x ∈ (0, 1) -θ 0 (-x) if x ∈ (-1, 0) 0 if x ∈ (-∞, -1] ∪ [1, +∞).
Introduce also the fundamental solution of the Schrö-

dinger equation E(t, x) := e i x 2 4t √ 4πit
. Then the function de-

fined by θ -(0, •) := θ odd 0 and θ -(t, x) := 1 -1 E(t, x -y)θ odd 0 (y)dy, (t, x) ∈ R * × R
enjoys the following properties:

(i) θ -belongs to C ∞ (R * × R) (ii) θ -is Gevrey of order (1, 1 2 ) on [t 1 , t 2 ] × [-L, L] for all 0 < t 1 < t 2 and L > 0. (iii) θ -(t, -x) = -θ -(t, x) for all (t, x) ∈ R * × R (iv) θ -| (0,T )×(0,1)
is the solution of (1)-( 4), where the control in (3) is given by u(t

) := θ -(t, 1) (v) t → θ -(t, •) belongs to C R, L 2 (R) (vi) θ -(t, •) L 2 (R) = √ 2 θ 0 L 2 (0,1) for all t ∈ R (vii) t → θ -(t, •) L ∞ (R) belongs to L 4 (R).
PROOF. We start by estimating the growth of the derivatives of E; notice E t = iE xx (since E is the fundamental solution of (1)). It is well-known, see e.g. [36, (5.5.3)], that

d k dz k e -z 2 = (-1) k e -z 2 H k (z),
where the H k are the (physicists') Hermite polynomials, which satisfy H 0 (x) = 1, H 1 (x) = 2x and

H k+1 (x) = 2xH k (x) -2kH k-1 (x). (12)
As a straightforward consequence,

∂ k x E(t, x) = (-1) k (4it) k 2 H k x √ 4it E(t, x) = (-1) k √ π(4it) k+1 2 e - x √ 4it 2 H k x √ 4it . (13)
On the other hand, see e.g. [START_REF] Szegö | Orthogonal polynomials[END_REF]Theorem 8.22.7],

Γ( k 2 + 1) Γ(k + 1) e -z 2 2 H k (z) ∼ k→∞ cos z √ 2k + 1 - kπ 2
uniformly for |z| ≤ A whatever the positive constant A.

As an immediate consequence,

e -z 2 2 H k (z) ≤ M (A) Γ(k + 1) Γ( k 2 + 1) e A √ 2k (14)
for M (A) large enough. Taking

z := x √ 4it with |x| ≤ L and |t| ≥ t 1 , hence |z| ≤ A := L √ 4πt1 , we then get ∂ k x E(t, x) = e -z 2 2 √ π |4t| k+1 2 e -z 2 2 H k (z) ≤ e A 2 2 √ π |4t 1 | k+1 2 M (A) Γ(k + 1) Γ( k 2 + 1) e A √ 2k ∼ k→∞ e A 2 2 M (A) (2π) 3 4 √ t 1 k! 1 2 √ 2t 1 k e A √ 2k k 1 4 ≤ M (t 1 , A) k! 1 2 R k ,
for any R ∈ (0, √ 2t 1 ) and M (t 1 , A) large enough; the third line comes from [START_REF] Haraux | Séries lacunaires et contrôle semi-interne des vibrations d'une plaque rectangulaire[END_REF]

with ξ := k 2 . As E t = iE xx , ∂ k x ∂ l t E(t, x) = ∂ k+2l x E(t, x) ≤ M (t 1 , A) (k + 2l)! 1 2 R k+2l ≤ M (t 1 , A) 2 k 2 +l k! 1 2 (2l)! 1 2 R k+2l ≤ M (t 1 , A) k! 1 2 R √ 2 k l! R 2 2l , (15)
for M (t 1 , A) large enough; we have used [START_REF] Jaffard | Contrôle interne exact des vibrations d'une plaque rectangulaire[END_REF] to get the second line, and (6) with ξ := l to get the last line.

This implies that for all compact set K of R * × R and k, l ∈ N, there is a constant M K,k,l > 0 such that sup

(t,x)∈K y∈[-1,1] ∂ k x ∂ l t E(t, x -y) ≤ M K,k,l .
We can now establish (i) by a standard argument of differentiation under the integral sign. Indeed:

• the function y → E(t, x -y)θ odd 0 (y) belongs to L 1 (-1, 1) for all (t, x) ∈ R * × R • the function (t, x) → E(t, x -y)θ odd 0 (y) belongs to C ∞ (R * × R) for almost all y ∈ [-1, 1] • ∂ k x ∂ l t E(t, x -y)θ odd 0 (y) ≤ M K,k,l θ odd 0 (y)
for almost all y ∈ [-1, 1] and all (t, x) ∈ K.

We then have (ii), since by [START_REF] Lions | Contrôlabilité exacte, perturbations et stabilisation de systèmes distribués[END_REF],

∂ k x ∂ l t θ -(t, x) = 1 -1 ∂ k+2l x E(t, x -y)θ odd 0 (y)dy ≤ 2M (t 1 , L + 1) θ 0 L 1 (0,1) k! 1 2 R √ 2 k l! R 2 2l for all (t, x) ∈ [t 1 , t 2 ] × [-L, L].
Statement (iii) follows from the change of variable z := y + x applied to θ -(t, -x) and the definition of θ odd 0 .

θ -is the (mild) solution of the Cauchy problem Theorem 3 Let θ 0 ∈ L 2 (0, 1) and T > 0 be given. Pick any τ ∈ (0, T ) and any s ∈ (1, 2). Define

iθ t + θ xx = 0, (t, x) ∈ R 2 θ(0, x) = θ odd 0 (x), x ∈ R,
Y (t) := φ s t -τ T -τ θ - x (t, 0),
with θ -as in Proposition 2 and φ s a Gevrey step function of order s such as [START_REF] Fliess | Flatness and defect of non-linear systems: introductory theory and examples[END_REF]. Then the control

u(t) :=      θ -(t, 1) if 0 ≤ t ≤ τ , (16)
j≥0

(-i) j (2j + 1)! Y (j) (t) if τ < t ≤ T (17)
steers (1)-( 4) from θ 0 to θ(T, •) = 0; u belongs to L 4 (0, T ), and is Gevrey of order s on [ε, T ] for all ε ∈ (0, T ). The corresponding solution of (1)-( 4) is given by

θ(t, x) :=      θ -(t, x) if 0 ≤ t ≤ τ , ( 18
) j≥0 (-i) j x 2j+1 (2j + 1)! Y (j) (t) if τ < t ≤ T ; (19) t → θ(t, •) belongs to C [0, T ], L 2 (0, 1)
, and θ is Gevrey of order (s, s 2 ) on [ε, T ] × [0, 1] for all ε ∈ (0, T ).

Remark 4 Notice the flatness property is used twice in the design of the control: of course in the very definition of the second phase through the parametrization (9); and also in the definition of Y (t) through θ - x (•, 0), which is the value of the flat output as computed in the first phase. This shows that the key point for establishing the null controllability of a flat partial differential equation is to find a trajectory which connects the initial condition to a sufficiently regular intermediate state; the control in the two phases is then automatically deduced from this trajectory. Though this is not completely apparent, this argument underlies the construction of the control in [START_REF] Martin | Null controllability of the heat equation using flatness[END_REF][START_REF] Martin | Null controllability of one-dimensional parabolic equations by the flatness approach[END_REF]. Notice also that the construction of Theorem 3 is much simpler than the construction in the preliminary version [START_REF] Martin | Controllability of the 1D Schrödinger equation by the flatness approach[END_REF], as it avoids the additional series in Lemma 5, and removes the restriction τ > 2 3 T .

PROOF. By Proposition 2(iv), θ -is the solution of ( 1)-( 4) on

[0, τ ] × [0, 1]; by (v), t → θ -(t, •) is in C R, L 2 (R)
, and a fortiori in C [0, τ ], L 2 (0, 1) . It eventually belongs to C [0, T ], L 2 (0, 1) as an obvious consequence of its Gevrey regularity away from 0, to be proved in the sequel.

By Proposition 2(ii), θ -is Gevrey of order (1, 1 2 ) on [ε, T ] × [0, 1] for all ε ∈ (0, T ), and a fortiori Gevrey of order (s, s 2 ) since s > 1. In particular, θ -(t, •) is Gevrey of order 1 2 , hence entire, for t ∈ [ε, T ]; since it is moreover odd by Proposition 2(iii), it can be expanded as

θ -(t, x) = j≥0 x 2j+1 (2j + 1)! ∂ 2j+1 x θ -(t, 0) = j≥0 (-i) j x 2j+1 (2j + 1)! ∂ j t θ - x (t, 0). (20)
Also θ - x (•, 0) is Gevrey of order s on [ε, T ]; and so is Y , as the product of two Gevrey functions of order s. By Proposition 1, the series in [START_REF] Lynch | Flatness-based boundary control of a class of quasilinear parabolic distributed parameter systems[END_REF] is Gevrey of order (s, s

2 ) on [ε, T ] × [0, 1] and satisfy (1)-(3). Moreover, [START_REF] Liu | Locally distributed control and damping for the conservative systems[END_REF] and [START_REF] Lynch | Flatness-based boundary control of a class of quasilinear parabolic distributed parameter systems[END_REF]

agree on [ε, τ ] × [0, 1]. Indeed, φ s t-τ T -τ = 1 for t ≤ τ , hence Y (t) = θ -
x (t, 0); by [START_REF] Machtyngier | Exact controllability for the Schrödinger equation[END_REF], this implies j≥0

(-i) j x 2j+1 (2j + 1)! Y (j) (t) = θ -(t, x).
Therefore, θ given by ( 18)-( 19) is Gevrey of order (s, s 2 ) on [ε, T ] × [0, 1], and is the solution of ( 1)-( 4). Finally, θ(T, •) = 0; indeed, φ (m) s (1) = 0 for m ∈ N, which implies Y (j) (T ) = 0 for j ∈ N.

Since the control u defined by ( 16)-( 17) is equal to θ(•, 1), it is Gevrey of order s on [ε, T ]. Moreover, Proposition 2(vii) implies in particular that u is in

L 4 [0, τ ], hence in L 4 [0, T ]. 2
5 Numerical implementation of the control Theorem 3 provides an explicit expression for the control. We now detail how this expression can be effectively implemented in practice. The fact that a small error on the so computed control yields a small error on the final state follows from classical semigroup arguments.

First phase

The control in the first phase is given by the integral

θ -(t, 1) := 1 √ 4πit 1 -1 e i (y-1) 2 4t θ odd 0 (y)dy. ( 21
)
This integral is difficult to evaluate numerically when t gets small, because of the oscillating singularity at t = 0.

A much better alternative is to expand it asymptotically. We thus recall the basics of asymptotic expansion of integrals of the form

I(λ) := b a f (y)e iλ(y-ȳ) 2 dy, a < b ≤ ȳ,
when the real parameter λ tends to +∞, and refer the reader to [2, chapter 6] for a detailed account. The fundamental result is the following: consider the integral

J(λ) := +∞ 0 F (u)e εiλu du,
where ε = ±1 and F is such that

• F (u) = 0 when u > ū, for some ū > 0 • F ∈ C M (0, +∞) for some M ∈ N • as u → 0 + , F has the asymptotic expansion

F (u) = M m=0 p m u am + o(u a M ), (22)
with

(a 0 ) < • • • < (a M ); as usual, the "little-o" symbol in ξ = o(ζ) is used to mean ξ ζ → 0 • the asymptotic expansion as u → 0 + of F (1) , . . . , F (M )
is obtained by differentiating ( 22) term-by-term.

Then, as λ → +∞, J(λ) has the asymptotic expansion

J(λ) = M m=0 p m Γ(1 + a m ) λ 1+am e εi π 2 (1+am) + o(λ -1-a M ). ( 23
)
This result can be formally obtained by replacing F by its asymptotic expansion [START_REF] Martin | Null controllability of the heat equation using flatness[END_REF] in the integral and integrating term-by-term; the difficult part is to justify this process, since [START_REF] Martin | Null controllability of the heat equation using flatness[END_REF] makes sense only near 0 and certainly not in the whole range of integration.

To evaluate I(λ), we first "isolate" the endpoints of integration, see [2, section 3.

I b (λ) = e iλ b2 +∞ 0 (χ b f )(b -r) 2( b + r) r=-b+ √ b2 +u e iλu du.
Taking the asymptotic expansion of the integrand as u → 0 + , we can now use the fundamental result [START_REF] Martin | Null controllability using flatness: a case study of a 1-D heat equation with discontinuous coefficients[END_REF]. Finally, since

χ b f = f around b, (χ b f )(b-r) and f (b-r)
have the same asymptotic expansion as r → 0 + , which means the asymptotic expansion of I b (λ) as λ → +∞ does not depend on χ b f , but only on f . The process is similar with I a : we set s := t -a and u := s(2ā -s), with ā := ȳ -a, which yields

I a (λ) = e iλā 2 +∞ 0 (χ a f )(a + s) 2(ā -s) s=ā- √ ā2 -u e iλu du,
and then use [START_REF] Martin | Null controllability using flatness: a case study of a 1-D heat equation with discontinuous coefficients[END_REF].

From the previous discussion, I(λ) can be asymptotically expanded if f is smooth enough on (a, b) and has well-behaved asymptotic expansions (similar to those required for F ) as y → a + and y → b -. If f does not enjoy these properties on the whole interval [a, b], but does so on each subinterval [c n , c n+1 ] of a partition c 0 := a < c 1 < • • • < c N := b, we can write

I(λ) = N -1 n=0 cn+1 cn f (y)e iλ(y-ȳ) 2 dy,
and handle each integral as before. In other words, only a, b, and the points where f or its derivatives (below some prescribed order depending on the sought order of the expansion) is discontinuous or singular contribute to the expansion; notice that if we artificially introduce an intermediate point c n where f is smooth enough, the contribution of this point will be zero since in this case

I c + n (λ) + I c - n (λ) = 0.
Remark 5 There are marked differences between the cases ȳ > b and ȳ = b. Indeed, assume the leading term in the asymptotic expansion of f (b -r) is proportional to r β0-1 ; when ȳ > b, the leading term of

(χ b f )(b-r) 2( b+r) r=-b+ √ b2 +u is proportional to u β0-1
, whereas it is proportional to u β 0 2 -1 when ȳ = b. Accordingly, the leading term of I b (λ) is proportional to e iλ(ȳ-b) 2 λ -β0 when ȳ > b, and proportional to λ -β 0 2 when ȳ = b. Similarly since ȳ > a, if the leading term in the asymptotic expansion of f (a+s) is proportional to s α0-1 , the leading term of I a (λ) is proportional to e iλ(ȳ-a) 2 λ -α0 .

Remark 6 There is a similar but more complicated result, see [2, section 6.3], when, instead of [START_REF] Martin | Null controllability of the heat equation using flatness[END_REF], F has the more general asymptotic expansion

F (u) = M m=0 N (m) n=0 p mn u am (log u) n + o u a M (log u) M .
We can thus compute [START_REF] Martin | Controllability of the 1D Schrödinger equation by the flatness approach[END_REF] . Notice that since by assumption θ 0 ∈ L 2 (0, 1), the leading term in the asymptotic expansion of θ odd 0 (c n + s) (resp. of θ odd 0 (c n -r)) is proportional to s α0n-1 , with α 0n > 1 2 (resp. proportional to r β0n-1 , with β 0n > 1 2 ). In view of Remark 5, this means the leading term in the asymptotic expansion of (21) due to c + n , i = 0, . . . , N -1, is proportional to e i (1-cn) 2 4t t α0n-1 2 (resp. due to c - n , i = 1, . . . , N -1 is proportional to e i (1-cn) 2 4t t β0n-1

2 ), hence tends to 0 in an oscillatory way. On the other hand, the leading term due to c N = 1 is proportional to t

β 0N -1 2
, hence is not oscillatory but possibly singular if θ 0 is itself singular at y = 1; since β 0N -1 2 > - 1 4 , this term nonetheless belongs to L 4 (R), which is of course in accordance with (vii) of Proposition 2.

Second phase

To compute sufficiently many terms in the series [START_REF] Littman | The heat and Schrödinger equations: boundary control with one shot[END_REF], an efficient way is to proceed recursively; moreover, the various terms must be properly scaled to accommodate finite precision arithmetics.

Following appendix A, we first compute

φk (t) := (T -τ ) k (2k)! d k dt k φ s t -τ T -τ = 1 (2k)! φ (k) s t -τ T -τ .
We emphasize that φ

(k) s t-τ T -τ is the k th derivative of φ s evaluated at t-τ T -τ , whereas d k dt k φ s t-τ T -τ is the k th derivative of t → ϕ s (t) := φ s t-τ T -τ , i.e., ϕ (k)
s (t).

We next compute

ỹk (t) := τ k k! ∂ k t θ - x (t, 0).
It can be checked directly that ∂ x E satisfies the recurrence relation

∂ k+1 t ∂ x E(t, x) = - 1 4t 2 ix 2 + 2t(4k + 3) ∂ k t ∂ x E(t, x) + 2k(2k + 1)∂ k-1 t ∂ x E(t, x) ;
this can also be seen as a consequence of ( 12) and ( 12). Since ∂ x E(t, x) = ix 2 F (t, x), where F (t, x) := E(t,x) t , this immediately implies

Fk+1 (t, x) = -τ 4(k + 1)t 2 ix 2 + 2t(4k + 3) Fk (t, x) + 2(2k + 1)τ Fk-1 (t, x) , (24)
where Fk (t, x) := τ k k! ∂ k t F (t, x). On the other hand,

θ - x (t, 0) = - 1 -1
∂ y E(t, -y)θ odd 0 (y)dy = -i 2 Fk (t, y)yθ 0 (y)dy, with Fk given by [START_REF] Martin | Null controllability of one-dimensional parabolic equations by the flatness approach[END_REF]. As t > τ , performing the numerical integration causes no practical problem, unless τ is chosen very small.

We finally compute the scaled derivatives of Y using the general Leibniz rule, The control in the second phase is eventually obtained by truncating the series [START_REF] Littman | The heat and Schrödinger equations: boundary control with one shot[END_REF] and expressing it in terms of the scaled derivatives of Y , which yields

Y (j) (t) (2j)! = 1 (2j)! j k=0 j k ∂ k t θ - x (t, 0) d j-k φ s dt j-k t -τ T -τ = j k=0 (2j -2k)!j! (2j)!(j -k)!τ k ỹk (t) φk (t); the coefficient d j k := (2j-2k)!j! (2j)!(j-k)!τ k is best computed re- cursively by d j 0 = 1 and d j k = d j k-1 2τ (2j-2k+1) .
u(t) =  j=0 (-i) j (2j + 1) Y (j) (t) (2j)! . ( 25
)
6 Numerical experiments

We now illustrate the effectiveness of the approach on a numerical example. The Matlab source code is freely available as ancillary files to the preprint [START_REF] Martin | Controllability of the 1D Schrödinger equation using flatness[END_REF]. We choose as initial condition

θ 0 (x) :=      x + 1 -i, 0 < x ≤ 0.3 x + 1 + i(x -0.3) -1/4 , 0.3 < x ≤ 0.6 e 2(x-0.6) + i(x -0.3) -1/4 , 0.6 < x ≤ 1,
displayed in Fig. 1, and we want to steer the system to zero at time T := 0.4, with intermediate time τ := 0.05. Notice the initial condition is rather challenging, since it does not satisfy the boundary conditions, has discontinuities, and a singularity (it can be shown to belong only to H s (0, 1), s < 1 4 ).

The control, computed as explained in section 5, is shown in Fig. 2. For the second phase, the control is given by (25), with s := 1.7, M = 0.8 and  := 50; the figure also shows that the contribution of the 30 following terms in the series is very small. For the first phase, the integral ( 21) is evaluated by an asymptotic expansion for 0 ≤ t < 10 -3 , and by a numerical quadrature (with Matlab quadgk) for 10 -3 ≤ t < τ ; the expansion uses sufficiently many terms so that the neglected terms have order 11 2 . Fig. 3 together with the error between the "exact" values computed by numerical integration and the values obtained by expansion. The error is as expected of order 11 2 ; notice quadgk begins to have trouble evaluating the integral with a good accuracy when t < 10 -3 . The control in the first phase, which is very oscillatory at the beginning because of the discontinuities and singularities of θ odd 0 , gradually calms down. This control is then applied to the numerically simulated system. The numerical scheme is a fixed-step Crank-Nicolson scheme; due to the oscillatory nature of the Schrödinger equation, a small diffusion term dx 3 4 , where dx is the space step, is added to avoid spurious oscillations. For the same reasons, it is necessary to use very fine space and time grids to get a good accuracy. The scheme is initialized with four backward Euler half-steps to better handle the discontinuities in the initial condition (so-called "Rannacher time-stepping" [START_REF] Rannacher | Finite element solution of diffusion problems with irregular data[END_REF]). Fig. 4 displays the whole evolution of the system; the final error is about 1.5 × 10 -3 . Fig. 5 zooms on the first phase; it can be seen that the control has as anticipated a smoothing effect, after a very oscillatory transient. Fig. 6 zooms further on the very beginning of the motion; it shows how the discontinuities and singularity of the initial condition generate the initial oscillatory behavior. The influence of τ on the shape of the trajectory, as well as the influence of the order and type of the Gevrey step, is an interesting and open question; see [START_REF] Martin | Null controllability using flatness: a case study of a 1-D heat equation with discontinuous coefficients[END_REF] for some numerical experiments in that regard.

  3]: we take a function χ a ∈ C ∞ (R) which is identically 1 around a and identically 0 around b (such a function is for instance readily built from the Gevrey step function φ s (5)), and set χ b := 1 -χ a . We then have I(λ) = I a (λ) + I b (λ), with I a (λ) := b a χ a (y)f (y)e iλ(y-ȳ) 2 dy I b (λ) := b a χ b (y)f (y)e iλ(y-ȳ) 2 dy. By suitable changes of variables, I a (λ) and I b (λ) are next reduced to forms similar to J(λ). Setting r := b -y and u := r(2 b + r), with b := ȳ -b, readily yields

  by setting f (y) := θ odd 0 (y) √ 4πit and λ := 1 4t ; θ odd 0 is assumed smooth enough on each subinterval [c n , c n+1 ] of a partition c 0 := -1 < c 1 < • • • < c N := 1, with well-behaved asymptotic expansions as y → c + n and y → c - n+1

F

  (t, y)yθ 0 (y)dy, since the integrand in the second line is even. Finally,ỹk (t) = -i 1 0

Figure 1 .

 1 Figure 1. Initial condition θ0 (red) and odd extension θ odd 0

Figure 2 .

 2 Figure 2. Control u(t): real part (top) and imaginary part (middle); contribution of 30 following terms in truncated control (25) (bottom).

  0.002 0.003 0.004 0.005 0.006 0.007 0.008 0

Figure 3 . 11 2

 311 Figure 3. Control u(t) at the beginning of first phase: real part (top) and imaginary part (middle). Bottom: error between asymptotic expansion and numerical integration (blue), graph of t → 5 × 10 7 t 11 2 (red).

Figure 5 .

 5 Figure 5. Evolution of θ, zoom on first phase: real part (top) and imaginary part (bottom).

  displays the beginning of the first phase,

	6
	4
	2
	0
	-2
	0	0.05	0.1	0.15	0.2	0.25	0.3	0.35	0.4
	4
	2
	0
	-2
	-4
	0	0.05	0.1	0.15	0.2	0.25	0.3	0.35	0.4
	10 -10
	10 -20
	10 -30
	10 -40
	10 -50
	0	0.05	0.1	0.15	0.2	0.25	0.3	0.35	0.4
					time

(all its eigenvalues are on the imaginary axis), which is in sharp contrast to the heat equation, see[START_REF] Martin | Null controllability of the heat equation using flatness[END_REF]. The control in the first phase, which therefore cannot be zero, is here obtained by using the dispersive effect of the Schrödinger equation on R, according to the following proposition.

A Recursive computation of

A practical problem when implementing the control in the second phase is to evaluate sufficiently many derivatives of φ s , see section 5.2. An efficient way is to proceed recursively, with a suitable scaling to accommodate finite precision arithmetics.

Let u(ρ) := -M ρ σ . Then ρu = -σu, and applying the general Leibniz rule yields ρu (j) + (j -1)u (j-1) = -σu (j-1) .

After scaling, this gives at once ρ r j u (j) j! = -r 1 + σ -1 j r j-1 u (j-1) (j -1)! .

Similarly, the derivatives of v(ρ ρ) . Then f = u f , and applying the general Leibniz rule yields

After scaling, this gives at once Similarly, the scaled derivatives of g(t) := e v(t) read

Applying next the general Leibniz rule to (f + g)φ s = g yields

After scaling, Setting ũj := r j u (j) j! , ṽj := r j v (j) j! , f j := r j f (j) (2j)! , and gj := r j g (j)

(2j)! , we finally obtain φj s := These formulas can then be directly implemented for instance in Matlab.