Alexandra Pena

Bertrand Ménaert

Benoit Boulanger
email: benoit.boulanger@grenoble.cnrs.fr

Fredrick Laurell

Carlota Canalias

Valdas Pasiskevicius

Patricia Segonds

Corinne Felix

Jerome Debray

Sébastien Pairis

Alexandra Peña

Benoît Boulanger

Fredrik Laurell

Corinne Félix

Jérôme Debray

Template-growth of periodically domain-structured KTiOPO4

à 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.

Introduction

Nonlinear optics is increasingly relying on structured materials which allow exploiting quasiphase matching for maximizing the efficiency of a nonlinear interaction and realizing novel functionality [START_REF] Fejer | Quasi-phase matched second harmonic generation: tuning and tolerances[END_REF][START_REF] Berger | Nonlinear photonic crystals[END_REF][START_REF] Canalias | Mirrorless optical parametric oscillator[END_REF][START_REF] Ellenbogen | Nonlinear generation and manipulation of Airy beams[END_REF][START_REF] Bahabad | Experimental confirmation of the general solution to the multiple-phase-matching problem[END_REF]. Ferroelectric domain structuring proved to be the most promising technique for fabricating such engineered nonlinear optical materials. Future applications of nonlinear optics will rely on the development of technologies allowing replication and largescale growth of structured materials using template designs, just as it has occurred in semiconductor electronics. The ferroelectric crystal KTiOPO 4 (KTP) is one of the most promising candidate materials for such an approach, where the finest periodicity ferroelectric domain structures over large volume have been demonstrated [START_REF] Fejer | Quasi-phase matched second harmonic generation: tuning and tolerances[END_REF][START_REF] Berger | Nonlinear photonic crystals[END_REF][START_REF] Canalias | Mirrorless optical parametric oscillator[END_REF][START_REF] Ellenbogen | Nonlinear generation and manipulation of Airy beams[END_REF][START_REF] Bahabad | Experimental confirmation of the general solution to the multiple-phase-matching problem[END_REF][START_REF] Canalias | Sub-micron periodically poled flux grown KTP[END_REF][START_REF] Hellstrom | High-power optical parametric oscillation in largeaperture periodically poled KTiOPO4[END_REF].

Up to now the common method to obtain a periodically poled ferroelectric crystal is by using the electric field poling technique, whereby a local inversion of the spontaneous polarization P s can be obtained in a single domain crystal by applying an electric field higher than the coercive field to lithographically defined finger electrodes patterned on the polar faces of the crystal. By using this technique the grating periodicity is very well controlled by the lithographic patterning. However, the thickness obtainable is limited in practice by the high voltage necessary, while control of the process and reproducibility in fabricating a high quality domain structure get exceedingly more difficult when the crystal gets larger and thicker.

Periodically poled KTP (PPKTP) is used in quasi-phase matched (QPM) second order nonlinear frequency conversion, like second harmonic generation (SHG) or optical parametric oscillators (OPOs). The QPM nonlinear interactions give unrivalled degree of design freedom and can be substantially more efficient than birefringence phase-matched ones due to the possibility to access to the largest diagonal nonlinear coefficient [START_REF] Pack | Measurement of the  (2) tensors of KTiOPO4, KTiOAsO4, RbTiOPO4, and RbTiOAsO4 crystals[END_REF]. Although 3-mm-thick KTP crystals have been poled, a variation in stoichiometry in presently available wafers is the primary difficulty in getting thick poled crystals, especially for dense gratings [START_REF] Hellstrom | High-power optical parametric oscillation in largeaperture periodically poled KTiOPO4[END_REF]. So far QPM materials are perceived as applicable only to relatively low-energy applications due to limitations in the optical aperture sizes which can be reliably reproduced.

An in situ growth technique has been proposed in order to obtain larger size PPKTP samples from single domain KTP crystals [START_REF] Roth | Growth-dependent properties of KTP crystals and PPKTP structures[END_REF], but unfortunately the regularity of the obtained grating period was not sufficient for a use in practical nonlinear devices. The same limitation occurred during the attempts to in situ grow periodically-poled LiNbO 3 , another ferroelectric crystal widely used for QPM [START_REF] Bermudez | On the cooling effect in the formation of periodic poled lithium niobate crystals grown by Cz technique[END_REF]. Note also that for the non-ferroelectric crystal GaAs, an epitaxial growth method has been used for the fabrication of orientation-patterned GaAs (OPGaAs) films suited for nonlinear optical frequency conversion [START_REF] Ebert | MBE growth of antiphase GaAs films using GaAs/Ge/GaAs heteroepitaxy[END_REF][START_REF] Eyres | All-epitaxial fabrication of thick, orientation-patterned GaAs films for nonlinear optical frequency conversion[END_REF].

Here we demonstrate a technology for growth of periodically domain-structures of KTP from a template, leading to scalability of the nonlinear optical structures to large volumes. This new process enables high-quality crystal growth below the Curie temperature and below the roughening temperature of the polar faces, thereby preserving an excellent crystal morphology and ferroelectric domain structure.

Description of the new process

With the aim of increasing the size of PPKTP crystals with a controlled and regular grating period, we developed an epitaxial growth process using seeds made of thin PPKTP plates previously obtained by electric field poling. For the following, we keep the denomination "PPKTP" for the grown periodically domain-structured KTP crystal. The poling step (Step I in Fig. 1) requires a sample exhibiting (001) and (00 ̅ ) crystal faces, so that the growth step (Step II Fig. 1) has to be performed onto this face. This constraint is a first difficulty to circumvent as this face is not present in the standard equilibrium morphology. The growth temperature of PPKTP is also imposed by the ferroelectric character of the crystal, since it is necessary to "build" the crystal below its Curie temperature in order to avoid unwanted erasure of the domain grating. By using a high temperature solution method, the so-called "flux method", and by choosing an appropriate chemical composition of the flux solution, it is then possible to meet both these constraints. Step II corresponds to the growth of a large size PPKTP crystal from the PPKTP seed by the flux method.

Definition of the growth conditions

KTP crystallizes in the orthorhombic system and belongs to the acentric point group mm2 with the space group Pna2 1 at room temperature [START_REF] Stucky | The potassium titanyl phosphate structure field: a model for new nonlinear optical materials[END_REF], the c axis being along the direction of the spontaneous polarization P s [START_REF] Thomas | Crystal structure and nonlinear optical properties of KSnOPO4 and their comparison with KTiOPO4[END_REF]. KTP has a ferroelectric phase transition with the change of symmetry Pna2 1 (acentric) -Pnam (centric) at the Curie temperature T c [START_REF] Yanovskii | Ferroelectric phase transitions and properties of crystals of the KTiOPO4 family[END_REF]. This sets an upper temperature for growth of KTP in its ferroelectric phase. It has been shown that T c can vary between 1193 K and 1233 K as a function of the flux growth conditions of KTP and this is directly correlated to the stoichiometry, where the higher temperature corresponds to a more stoichiometric crystal: T c increases as a function of the molar ratio K/P while it decreases as a function of the growth temperature for a given ratio K/P [START_REF] Roth | On the optical quality of KTP crystals for nonlinear optical and electro-optic applications[END_REF]. This fact has to be taken into account in order to design the right flux.

The standard morphology of KTP as obtained by flux growth around 1100 K is depicted in Fig. 2(a). It does not exhibit (001) and (00 ̅ ) faces, the crystal faces orthogonal to the c axis [START_REF] Bolt | Growth of KTiOPO4 (KTP) single crystals by means of phosphate and phosphate/sulfate fluxes out of a three-zone furnace[END_REF][START_REF] Jacco | Flux growth and properties of KTiOPO4[END_REF]. However, these faces, of so-called non standard morphology, shown in Fig. 2(b), exist only when the growth is performed below their roughening temperatures, i.e. T R ≈ 1053 K. The roughening temperatures relative to several faces of KTP as calculated from the Rijpkema model [START_REF] Rijpkema | Determination of the Ising critical temperature of F slices with an application to garnet[END_REF], are given in Fig. 2(c) for comparison. Note that the morphology of a KTP crystal grown by the hydrothermal method does not exhibit (001) faces at all, even if the growth is performed below the roughening temperature [START_REF] Laudise | Phase relations, solubility and growth of potassium titanyl phosphate, KTP[END_REF]. It may be due to a higher growth rate of (001) faces compared with those of the neighbor faces {201} and {011}. However, by using a K 6 P 4 O 13 solvent and high amounts of PbO [START_REF] Tseitlin | Growth of KTP crystals with large {001} facets[END_REF], it is possible to obtain (001) faces at a temperature higher than T R . This is an interesting approach but here we preferred to consider fluxes with a lower toxicity. By using a standard KTP crystal, it is first necessary to cut (001) and (00 ̅ ) faces in order to realize the electric field poling (Step I of Fig. 1). The growth stage of the process (Step II of Fig. 1), following the poling stage, has hence to be performed onto (001) and (00 ̅ ) faces. If Step II is realized at a temperature higher than the roughening temperature of (001) faces, a capping process would take place leading to a detrimental reconstruction with a pyramidal domain shape formed by {201} and {011} faces [START_REF] Halfpenny | X-ray topographic study of flux grown KTP crystals[END_REF]. The flux composition is then a key point in performing the growth of domain structured KTP, so that the temperature can be kept below the Curie temperature and below the roughening temperature of the (001) faces. The commonly used flux for growing KTP is K 6 P 4 O 13 , but the lower limit of growth is 1098 K [START_REF] Bolt | Growth of KTiOPO4 (KTP) single crystals by means of phosphate and phosphate/sulfate fluxes out of a three-zone furnace[END_REF][START_REF] Jacco | Flux growth and properties of KTiOPO4[END_REF], due to the formation of a neighboring polyphosphate complex phase below this temperature [START_REF] Loiacono | Solubility and crystal growth of KTiOPO4 in polyphosphates solvents[END_REF].

Here we propose to use a flux containing a mixture of KPO 3 and KF from which KTP can be obtained at temperatures as low as 873 K, corresponding to the eutectic point of the system [START_REF] Marnier | Synthesis in flux of crystals and epitaxies from isotypic solid solution of KTiOPO4[END_REF]. This increased temperature range compared with the solvent K 6 P 4 O 13 is due to a higher ionic character of the solution, thanks to the presence of KF that prevents the formation of undesired polyphosphate complex phases. Knowledge of the equilibrium curve is required in order to define the composition of the halide flux enabling the growth to be performed at the required temperature. We considered the same molar ratio as in the case of K 6 P 4 O 13 , i.e. K/P = 1.5, to be able to use the experimental data previously measured using K 6 P 4 O 13 [START_REF] Bordui | Growth of large single crystals of KTiOPO4 (KTP) from high-temperature solution using heat pipe based furnace system[END_REF]. The corresponding reaction constant as a function of temperature was deduced by using an ionic solubility model [START_REF] Temkin | Mixtures of fused salts as ionic solutions[END_REF][START_REF] Van Erk | A solubility model for rare-earth iron garnets in a PbO/B2O3 solution[END_REF]. It allowed us to determine the molar enthalpy of crystallization, i.e. ΔH = 6001.5 cal/mol, from which finally the equilibrium curve of Fig. 3 was calculated. Note that the figure contains equilibrium curves for two solvents, but they overlap since they correspond to the same K/P ratio. We validated the calculations by a thermo-gravimetric measurement of the saturation temperature corresponding to the molar composition 0.1 KTP -0.6 KPO 3 -0.3 KF : we found 1005.5 ± 0.1 K, which is very close to the theoretical value, i.e. 1003 K, as shown in Fig. 3.

Crystal growth

In a preliminary step, the morphology was checked on small size single domain KTP crystals grown by spontaneous nucleation from the halide solvent exhibiting different saturation temperatures that were predicted using the equilibrium curve of Fig. 3. The grown crystals exhibited (001) faces below 1048 K, which is very close to the calculated roughening temperature of 1053 K. Furthermore, we performed epitaxial growth experiments onto polished (001) and (00 ̅ ) faces of a single domain KTP crystal. The flux composition 0.1 KTP -0.6 KPO 3 -0.3 KF previously considered was used since the corresponding saturation temperature is below the roughening temperature of (001) faces. We applied an oversaturation of 2.2 % corresponding to growth 6 degrees below the saturation temperature. During growth the single domain KTP seed was rotated with an angular speed of 30 rpm. By using a scanning electron microscope (SEM), we verified a very high quality of the interface between the seed and the epitaxially grown layer.

In a second step, we performed epitaxial growths of PPKTP onto (001) and (00 ̅ ) faces (Step II in Fig. 1) of PPKTP plates with the dimensions of 11 mm x 6 mm x 1 mm (a x b x c) where a central area of 7 mm x 3 mm x 1 mm had previously been periodically inverted with a grating period of Λ = 38.86 µm, by electrical poling (Step I of Fig. 1). The flux composition 0.1 KTP -0.6 KPO 3 -0.3 KF as previously defined was used, and a constant cooling ramp of 0.05 K/h, starting from 1005.5 K, was applied. Thicknesses of the domain structured KTP layers up to 800 µm were grown onto each of the two c faces, i.e. the (001) and (00 ̅ ) faces. An example of such a sample is shown in Fig. 4(a). By using a scanning electron microscope (SEM), we were able to verify the very high interface quality, as shown in Fig. 4 Note that due to the possibility to grow layers on both sides of the seed simultaneously, it is possible to get double volume compared with the usual single layer growth.

Optical characterization

In order to check the optical quality of the grown layer in terms of the domain grating regularity, we performed first order QPM SHG measurements with the beam propagating along the a axis and probing the d 33 coefficient of the crystal, the SH and fundamental waves being polarized along the c axis. The SH intensity was measured as a function of the fundamental wavelength λ ω at different depths in both the seed and epitaxial layer at room temperature. The 1/e 2 diameter of the fundamental beam was 70 µm, which enabled an appropriate spatial resolution for scanning over the whole thickness of the sample. In Fig. 5 two typical examples of measurements are shown that show the equivalence in the nonlinear optical performance of the seed and epitaxial layer. The quasi-phase matching peak was 2126.9 ± 0.5 nm in the epitaxial layer, which is very close to what was measured in the seed, i.e. 2127.1 ± 0.5 nm, showing the resemblance between the two regions. The spectral acceptance was the same in the seed and in the epitaxial layer, i.e. Lδλω = 12.7 nm.cm, which attests that the corresponding periodically poled regions have the same extension along the a axis and that the grown layer faithfully reproduces the seed. The acceptance bandwidth is very close to the one calculated for a perfect structure, i.e. Lδλ ω = 9.4 nm.cm, which further demonstrates the high quality of the grown layer. The peak amplitude of the SH intensity is about 28% smaller in the grown layer than in the seed, corresponding to slightly higher losses due to scattering. The defects density increases up to 225 µm as mentioned above in section 4, which is detrimental for optical measurements.

Conclusions

In summary, periodically domain-structured KTP layers with thicknesses up to 800 µm and regular periodicity onto (001) and (00 ̅ ) faces of the initial PPKTP seed have been obtained for the first time by epitaxial growth. The use of a 4KPO 3 : 2KF solvent allowed the growth to be performed below T c and the roughening temperature of the (001) face. Up to now, the crystalline quality is perfect over 225 µm, which allowed us to perform quasi-phase matched second harmonic generation experiments within the corresponding volume. By this way, we demonstrated that the grown layer and the seed are very similar, which shows that this technology well meets the requirement for template-growth of engineered PPKTP or other nonlinear structures. The current experiments aim at increasing the crystal quality above 225 µm.

Fig. 1 .

 1 Fig. 1. Scheme of the process showing the different steps to obtain large size PPKTP crystals. Step I corresponds to the fabrication of a PPKTP sample from a single domain KTP plate by electric field poling. Step II corresponds to the growth of a large size PPKTP crystal from the PPKTP seed by the flux method.

Fig. 2 .

 2 Fig. 2. (a) Standard morphology of KTP crystals. (b) The morphology of KTP exhibiting (001) and (00 ̅ ) faces. (c) Roughening temperatures, T R (K), of different {hkl} faces.

Fig. 3 .

 3 Fig. 3. The pseudo-binary phase diagram of KTiOPO 4 (KTP) using K 6 P 4 O 13 and 4KPO 3 : 2KF solvents. T sat is the saturation temperature of KTP below which the crystallization can occur; X is the molar fraction defined by [K 6 P 4 O 13 ]/{[KTP]+[K 6 P 4 O 13 ]} or [4KPO 3 : 2KF]/{[KTP]+[4KPO 3 : 2KF]} according to the considered solvent.

  (b), and a propagation of the domains over a thickness 225 µm without any crack or variation of the domains size as indicated by Fig. 4(b) and (c). The domains still propagate from 225 µm until the boundaries of the grown layers but with some cracks and inclusions.

Fig. 4 .

 4 Fig. 4. Periodically domain-structured-KTP films grown on a PPKTP in a flux with the composition 0.1 KTP -0.6 KPO 3 -0.3 KF. a) Picture of a PPKTP crystal and the corresponding orientation used for the growth. b) SEM image along the b axis of the grown layer and seed. c) SEM image along the c axis of the grown layer at 225 µm far from the surface of one side of the seed.

Fig. 5 .

 5 Fig. 5. Comparison between the quasi-phase matched second harmonic generation (QPM SHG) properties of the seed and one side of the grown layer at room temperature. QPM SHG measurements performed at two different depths, 275 µm in the seed (s) and 225 µm in the epitaxial layer (l) measured from the interface (dashed line) in a sample with a grating periodicity  = 38.86 µm. 5. Conclusions In summary, periodically domain-structured KTP layers with thicknesses up to 800 µm and regular periodicity onto (001) and (00 ̅ ) faces of the initial PPKTP seed have been obtained for the first time by epitaxial growth. The use of a 4KPO 3 : 2KF solvent allowed the growth to be performed below T c and the roughening temperature of the (001) face. Up to now, the crystalline quality is perfect over 225 µm, which allowed us to perform quasi-phase matched second harmonic generation experiments within the corresponding volume. By this way, we demonstrated that the grown layer and the seed are very similar, which shows that this technology well meets the requirement for template-growth of engineered PPKTP or other nonlinear structures. The current experiments aim at increasing the crystal quality above 225 µm. Acknowledgements A. Peña acknowledges CNRS and AGAUR (2008 BP A) for the financial supports of her postdoctoral positions. The KTH group acknowledges support from the Linnaeus centre ADOPT, Swedish Research council (Vetenskapsrådet) and the Knut and Alice Wallenberg Foundation. All the authors wish to acknowledge Cristal Laser S.A. for providing single domain KTP crystals for the preliminary epitaxial growth experiments.

Acknowledgements A. Peña acknowledges CNRS and AGAUR (2008 BP A) for the financial supports of her postdoctoral positions. The KTH group acknowledges support from the Linnaeus centre ADOPT, Swedish Research council (Vetenskapsrådet) and the Knut and Alice Wallenberg Foundation. All the authors wish to acknowledge Cristal Laser S.A. for providing single domain KTP crystals for the preliminary epitaxial growth experiments.


Alexandra Pena

Bertrand Ménaert

Benoit Boulanger

Fredrik Laurell

Carlota Canalias

Valdas Pasiskevicius

Luc Ortega

Patricia Segonds

Jerome Debray

Corinne Felix

A Peña

C Félix

Bulk PPKTP by crystal growth from high temperature solution

Keywords: Epitaxy, Microstructures, Nonlinear Optics

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.

Introduction

Periodic inversion of ferroelectric domains in the nonlinear crystal KTiOPO 4 (KTP) [1][2][3][4] has been extensively used in order to achieve frequency conversion by Quasi-Phase-Matching (QPM) [5], which is more efficient than by performing Birefringence-Phase-Matching (BPM) in a single domain KTP crystal. Until now, commercially Periodically Poled KTP (PPKTP) crystals were fabricated by the electric field poling technique, based on the inversion of the spontaneous polarization P s by applying a high electric field to a single domain crystal [6]. Such samples are used in low energy devices [7,8] because of their small thickness along P s direction. In situ growth techniques had been then proposed in order to obtain larger size PPKTP samples from single domain KTP crystals [9]. But unfortunately the regularity of the obtained grating period was not sufficient for a use in practical Second Harmonic Generators (SHG) or Optical Parametric Oscillators (OPO). The same limitation had occurred when trying in situ growth of periodically-poled LiNbO 3 (PPLN) crystal, another ferroelectric medium widely used for QPM [10]. In situ growth techniques had also been used in patterned nonlinear crystals in order to achieve frequency conversion by QPM: 1-mm-thick Orientation-Patterned GaAs (OP-GaAs) have been recently obtained by Low-pressure hydride vapor phase epitaxy (HVPE) [11] enabling the generation of a mid-IR continuum by optical parametric generation (OPG) [12].

Here we propose a crystal growth process starting from a thin seed of PPKTP obtained in a previous step from electric field poling, with the objective to obtain thicker size samples. We consider a high temperature solution growth method enabling to perform the crystal growth below the Curie temperature of the crystal and below the roughening temperature of the relevant crystallographic faces, i.e. the two faces (001) and (00 1 ) that are perpendicular to P s , so that a thick PPKTP crystal with the same grating periodicity than the one of the seed can be obtained.

Definition of the growth conditions

The growth process should be done in the temperature range of the ferroelectric phase of KTP, i.e. below its Curie temperature [13], and also below the roughening temperature of (001) and (00 1 ) faces [14]. The first condition is necessary to maintain the periodic domain structure of the initial seed during the growth process. The second condition, which is imposed by the geometry of the thin PPKTP seeds that exhibit large artificial (001) and (00 1 ) faces, allows us to suppress undesirable capping process that would lead to a reconstruction of the crystal with a pyramidal domain shape formed by {201} and {011} faces [15].

As shown in fig. 1, the standard flux used to grow KTP, i.e. K 6 P 4 O 13 , is not suitable for that purpose because its lower growth limit is 1098K [16,17], which is over the roughening temperature of (001) and (001 ) faces. Below this temperature, the crystals that are obtained belong to a polyphosphate complex phase [18]. Here we propose to use a flux containing a mixture of KPO 3 and KF from which KTP can be obtained at temperatures as low as 873 K, corresponding to the eutectic point of the system [19]. The increasing of the temperature range compared with that of the standard flux K 6 P 4 O 13 is due to a higher ionic character of the flux, thanks to the presence of KF preventing undesired polyphosphate complex phases.

Crystal growth results

The growth experiments were performed in 4KPO 3 -2KF flux, which allowed us to keep the same molar ratio as in the case of K 6 P 4 O 13 , i.e. K/P = 1.5, so that it was possible to use experimental data previously published in the case of K 6 P 4 O 13 [20]. The corresponding reaction constant (k) as a function of temperature was deduced by using an ionic solubility model [21,22], which led to the determination of the molar enthalpy of crystallization, i.e. ΔH = 6001.5 cal/mol, and then to the calculation of the equilibrium curve of fig. 2. We validated these calculations by the thermo-gravimetric measurement of the saturation temperature corresponding to the molar composition 0.103 KTP -0.598 KPO 3 -0.299 KF: we found 1005.5 ± 0.1 K, which is very close to the theoretical value, i.e. 1003 K.

Preliminary growth experiments were realized in order to check whether c and 𝑐 faces appear close to the calculated roughening temperature value, T R = 1053 K, predicted by Bennema [14] from the Rijpkema model [23]. The external morphology of small crystals grown by spontaneous nucleation in 4KPO 3 -2KF flux with different initial molar composition was visualized by an optical microscopy and a scanning electron microscope (SEM). The saturation temperatures were predicted using the equilibrium curve shown in fig. 2. The crystals exhibited (001) and (001 ) faces below 1048 K as shown in fig. 3. So this temperature is close to the calculated one, i.e. 1053 K [14]. Note that these two faces can also be stabilized over their roughening temperature by using a K 6 P 4 O 13 -PbO flux [24]. The best growth conditions to obtain good quality PPKTP layers grown onto (001) faces were determined by preliminary experiments using single domain KTP seeds. The epitaxial growth experiments were done in a solution with the same molar composition than the one used to validate the calculations done to obtain the equilibrium curve of fig. 2, i.e.: 0.103 KTP -0.598 KPO 3 -0.299 KF. The saturation temperature of this solution is under the roughening temperature of ( 001) and ( 001), so that no capping process will be expected [15]. Good quality KTP layers of few hundred µm thickness were obtained by applying an over-saturation of around 2 % and by alternatively rotate the initial seed with an angular speed of around 30 rpm.

These growth conditions were then used to grow thick PPKTP layers onto PPKTP seeds previously obtained by electric field poling, with a periodicity Λ = 38.86 µm. The results were very satisfying, with a propagation of the domains along the thickness of the as grown layer, as it shown in fig 4b. PPKTP layers of up to 800 µm onto (001) and (001 ) faces were grown by using a growth ramp of 0.05 K/h (fig. 4a). This cooling ramp, starting from an initial temperature of 1005.5 K, allows us to have a constant matter transfer from the solution to the grown crystal during all the growth process.

Characterizations

A Scanning Electron Microscope (SEM) have been used to determine the quality of the interface between the seed and the growing layer, to visualize the propagation of the ferroelectric domains from the seed to the layer and also to determine the growth rate along c direction. The interface is of very good quality and the domains propagate without any measurable variation in size from the seed surface till the layer surface [25]. Due to the high quality of the interface and to the fact that we used a homoepitaxial growth, it was not possible to see any contrast between seed and grown layer, with beam energies higher than 3 kV and work distances higher than 6 mm. The growth rate along c was determined by SEM and also by a gravimetric method, which provides the change in apparent weight during all the growth process. By both methods, we obtained a growth rate of 2 µm/h. The crystal quality was determined by X-ray diffraction from a θ -2θ scan as shown in fig 5 . The full width at half maximum (FWHM) of the measured (004) reflection is equal to 0.006 ° for both the substrate and growth layer, attesting of an equivalent quality. The difference between the diffraction peak angles, i.e. 2θ = 33.859 ° for the seed and 2θ = 33.859 ° for the layer, corresponds to a mismatch c = 110 -3 . This cell parameter variation is not surprising since the seed and the layer have not been grown from the same flux. There is also a difference of peak intensity between the seed and the layer, which comes from the difference of irradiated matter volumes due to the particular cut geometry of the sample as shown on the insert of fig. 5. The grating period regularity of the ferroelectric domains in the grown layer as well as the replication of the initial periodically poled zone present in the initial seed had been checked by performing first order Quasi Phase Matching (QPM) Second Harmonic Generation (SHG) measurements. The obtained result [25] shows the perfect replica of the domain structure over the first 200 µm along c direction of the grown layer. Up to this clean area and till to the top of the layer (800 µm), there was a strong diffusion of the laser source due to growth defects that are not connected with the propagation of the domains. New experiments are in progress in order to suppress these defects.

Summary

In summary, growth conditions to obtain thick PPKTP layers from PPKTP seeds have been determined. Up to now PPKTP layers with a thickness of 800 µm along the c direction, onto (001) and (001 ) faces of the initial PPKTP seed, have been obtained by epitaxial growth from a KTP-KPO 3 -KF solution at temperatures lowers than 1005.5 K. Crystalline quality as well as the domain grating periodicity of the grown layers are the same than in the seeds. We think that the proposed technology is well adapted to obtain PPKTP crystals as well as other periodically poled materials requiring the flux method.

Fig 1 .

 1 Fig 1. Schematic thermal range of KTiOPO 4 (KTP) ferroelectric and paraelectric phases and possible equilibrium morphologies in the paraelectric phase using two different fluxes.

Fig 2 .

 2 Fig 2. KTiOPO 4 (KTP) pseudo binary phase diagram using K 6 P 4 O 13 and 4KPO 3 -2KF flux. Parameter X is the molar fraction, i.e. [𝐾 6 𝑃 4 𝑂 13 ] 𝐾𝑇𝑃 +[𝐾 6 𝑃 4 𝑂 13 ]or

Fig 3 .

 3 Fig 3. a) SEM image of a single KTP crystal grown by spontaneous nucleation in a 0.120 KTP -0.586 KPO 3 -0.293 KF flux. b) Morphological representation of the grown crystal obtained with SHAPE program.

Fig 4 .

 4 Fig 4. a) Three pictures of a single PPKTP crystal grown by liquid phase epitaxy in a vertical dipping configuration. The growth was performed from a solution with the molar composition 0.120 KTP -0.586 KPO 3 -0.293 KF and the initial seed used has a period grating (Λ) equal to 38.86 µm. b) SEM image, at 200 µm from the seed surface, of the c face of the grown layer.

Fig 5 .

 5 Fig 5. Diffraction θ-2θ scan for the (004) reflection of KTiOPO 4 (KTP). Scheme of the cut edge PPKP sample (insert) in which the θ-2θ scan was recorded.

Acknowledgements

A. Peña acknowledges CNRS and AGAUR (2008 BP A) for the financial supports of her postdoctoral positions. The KTH group acknowledges support from the Linnaeus centre ADOPT, Swedish Research council (Vetenskapsrådet) and the Knut and Alice Wallenberg Foundation. All the authors wish to acknowledge Cristal Laser S.A. for providing single domain KTP crystals for the preliminary epitaxial growth experiments.


Ivan Enderlin
email: ivan.enderlin@femto-st.fr

Alain Giorgetti
email: alain.giorgetti@femto-st.fr

Fabrice Bouquet
email: fabrice.bouquet@femto-st.fr

A Constraint Solver for PHP Arrays

Keywords: array, constraint, solver, php, realistic domain, praspel

In previous works, we have proposed Praspel, a framework for contract-based testing in PHP. Among others, it includes a specification language and a unit test generator which automatically generates test data from formal preconditions. The generator sometimes rejects data, when they do not satisfy parts of the preconditions. In many cases, generation with rejection is not efficient enough. Thus we investigate practical contexts where more efficient generation algorithms can be designed and we extend Praspel with their implementation.

After strings, that we have already considered, the most frequent data type in PHP is arrays. They cover most of the needs for collections, because they can store key-value pairs of any kind, they do not have a specific length or depth, and they are efficiently implemented. In this paper, we report on a study to know what are the most popular constraints on PHP arrays. Then we formalize these constraints and we present an implementation in PHP of a constraint solver for these constraints. In this context, the constraint-based approach removes all the rejections.

I. INTRODUCTION

Contract-based testing [START_REF] Aichernig | Contract-based testing[END_REF] is a promising approach to increase software safety and security. It is based on the notion of Design by Contract, introduced by Meyer [START_REF] Meyer | Applying "design by contract[END_REF] with the Eiffel language [START_REF]Eiffel: programming for reusability and extendibility[END_REF]. A contract is a formal specification associated to the code of a program. It mainly consists of invariants, preand postconditions. Invariants describe properties that should hold at each step of the execution. Pre-and postconditions respectively represent conditions that have to hold for a method to be invoked, and conditions that have to hold after execution of the method.

Various contract languages extend programming languages with contracts: JML extends Java [START_REF] Leavens | JML: A notation for detailed design[END_REF], ACSL extends C [START_REF] Baudin | ACSL: ANSI C Specification Language[END_REF], Spec# extends C# [START_REF] Barnett | The Spec# Programming System: An Overview[END_REF], etc. Among the numerous advantages of contracts, formal properties can be exploited for (unit) testing. Indeed, the information contained in invariants and preconditions can be used to generate test data. In addition, these assertions can be checked at run time and thus provide a (partial) test oracle for free.

In previous works, we have introduced Praspel [START_REF] Enderlin | Praspel: A specification language for contract-based testing in php[END_REF], a toolsupported specification language for contract-based testing in PHP [START_REF] Zaïdi | The PHP website[END_REF]. Praspel extends contracts with the notion of realistic domain, which makes it possible to assign a domain of values to data (class attributes or method parameters). A library of predefined realistic domains is already available along with a test environment [START_REF] Enderlin | Hoa project, a set of PHP libraries[END_REF].

Funded by the FUI project SQUASH.

Validating PHP web applications often involves to manipulate strings and arrays. Strings have been already addressed [START_REF] Enderlin | Grammar-Based Testing Using Realistic Domains in PHP[END_REF] by introducing grammar-based testing in Praspel. This paper focus on PHP arrays, which are always associative arrays. They cover most of the needs for collections, because they can store key-value pairs of any kind (they can be homogeneous and heterogeneous), they do not have a specific length or depth, and they are efficiently implemented. When generating arrays for testing, the main difficulty is to satisfy their properties, which are formalized by predicates in their specification. We address this issue in the context of PHP and Praspel, by designing a specific constraint solver.

Our contributions are twofold. First, after studying the most popular conditions on PHP arrays, we propose a syntax to specify them in Praspel and a semantics of these array specifications. Second, we propose a new constraint solver in PHP for PHP, to generate arrays satisfying these conditions. These two contributions are embedded in realistic domains and can be used in Praspel annotations of a PHP program.

The paper is organized as follows. Section II briefly introduces the notion of realistic domain and its implementation in Praspel for PHP. Section III presents a language of array conditions in Praspel, inspired by a study to know what are the most popular conditions expressed on arrays in popular PHP code. Section IV presents the constraint solver. Then, Section V reports on an experimentation validating our approach and tool, and showing its usefulness and efficiency in practice. Related works are presented in Section VI. Finally, Section VII concludes and presents future works.

II. REALISTIC DOMAINS AND PRASPEL

This section is a reminder of [START_REF] Enderlin | Praspel: A specification language for contract-based testing in php[END_REF]. It shortly presents the notion of realistic domain and its application to PHP programs. Realistic domains are designed for test generation purposes. They specify which values can be assigned to a data in a given program. They are well-suited to PHP, since this language is dynamically typed (i.e. no types are syntactically assigned to data) and realistic domains thus introduce a specification of data types mandatory for test data generation. We first introduce general features of realistic domains. Then we present their implementation in PHP and the Praspel framework.

A. Features of Realistic Domains

Realistic domains refine usual datatypes (integer, string, array, etc.). They are intended to specify data domains relevant for specific application contexts. For example, email addresses constitute a realistic domain: many applications identify a user by her email address, so we need to validate and generate such data, and an email address is more than a string: it obeys to a specific syntax that makes it not obvious to generate.

The first feature of a realistic domain (named its predicability) is to carry a characteristic predicate of its values, used to check if a value belongs to the possible set of values described by the realistic domain. The second feature of a realistic domain (named its samplability) is to propose a value generator, called the sampler, that generates values in the realistic domain. For instance, a sampler for the realistic domain of email addresses can generate strings representing syntactically correct email addresses defined by a regular expression.

B. Realistic Domains in PHP

In PHP, we have implemented realistic domains as classes providing at least two methods, corresponding to the two features of realistic domains. The first method is named predicate($q) and takes as input a value $q: it returns a boolean indicating the membership of the value to the realistic domain. The second method is named sample() and generates values that belong to the realistic domain.

Our implementation of realistic domains in PHP exploits the PHP object programming paradigm and takes benefit from the following three principles:

1) Inheritance: PHP realistic domains can inherit from each other. A realistic domain child inherits the two features of its parent, namely predicability and samplability, and is able to redefine them. Consequently, we say that all the realistic domains constitute a hierarchical universe.

2) Interfaces: All data are represented by realistic domains. Some of them implement interfaces which characterize them. Useful interfaces for our current concerns are:

• Constant: represents an immutable realistic domain with only one value, such as 42, true, etc. • Interval: represents an interval by its lower and upper bounds, that can be dynamically reduced, • Nonconvex: allows to discredit values, i.e. to specify that a value no longer belongs to a realistic domain and should therefore not be generated, • Finite: allows to count the number of values, • Enumerable: allows to iterate over all the values. Thus, a realistic domain that implements the Interval and Nonconvex interfaces is an interval with "holes". Counting and exhaustive generation of a finite realistic domain take discredited values into account.

3) Parameterization: Realistic domains may have parameters, like a function does. Data given to a realistic domain are called arguments. This feature helps generating structured data such as arrays, objects, graphs, automata, etc.

Example 1 (Realistic domains with arguments). The realistic domain string(0x61, 0x7a, boundinteger [START_REF] Leavens | JML: A notation for detailed design[END_REF][START_REF] Cheon | A Simple and Practical Approach to Unit Testing: The JML and JUnit Way[END_REF]) admits as arguments two integers (that represent two Unicode code-points) and a domain of integers to specify its possible length. The realistic domain boundinteger(X, Y ) contains all the integers between X and Y . We can also write X..Y as syntactic sugar for this domain. The realistic domain string(X, Y , L) is intended to contain all the strings of length (in the domain) L built of characters from X to Y code-points.

When describing parameters, because all data are realistic domains, we can multi-type hint the parameters, i.e. filter with multiple realistic domain names. For instance, the first two parameters of the realistic domain string() are described as Constinteger | Conststring, so can be a constant integer or a constant string. Then, the syntax string('a', 'z', 4..12) is strictly equivalent to the previous one. If we write string(true, 'z', 4..12), an error will be thrown. The realistic domain itself handles the cast.

C. Praspel and Contract-Based Testing in PHP

Praspel means PHP Realistic Annotation and SPEcification Language. It is a language and a framework for contract-based testing in PHP, based on realistic domains.

Praspel annotations are written inside comments in the source code. Invariants document classes and pre-and postconditions document methods. The general form of Praspel annotations is shown in Figure 1. For lack of space, named behaviors and specification of exceptions are not presented. In this specification, I 1 , . . . , I h are invariant clauses, assumed to be satisfied at the beginning and at the end of each method invocation. Formulas R 1 , . . . , R n and A 1 , . . . , A k are precondition clauses, that have to hold at the invocation of the foo method. Formulas E 1 , . . . , E m are postconditions that have to be established when method foo terminates without throwing an exception. In postconditions, Praspel provides the two additional constructs \result and \old(e), which respectively designate the value returned by the method, and the value of expression e at the pre-state of the method invocation.

PHP does not provide a type system, but Praspel contracts make it possible to give typing informations, assigning realistic domains to data (class attributes or method parameters). The construction i: t 1 (. . . ) or . . . or t n (. . . ) associates at least one realistic domain, among t 1 (. . . ), . . . , t n (. . . ), to the identifier i. Every identifier holds a realistic domain disjunction, thanks to the or keyword. Identifiers can also be passed into arguments of realistic domains.

Example 2 (Identifiers in Praspel). The declarations: length: 4..12 str : string('a', 'z', length)

show a dependency between two identifiers: The third argument of the realistic domain of the second identifier str is the first identifier length.

Praspel provides a set of more than 30 predefined realistic domains, called the standard library. Some of them correspond to scalar types (integer, float, boolean. . . ), other ones to classes and arrays (detailed in Section III).

Contractual assertions are made of realistic domain declarations, possibly completed with additional predicates, expressed in PHP using the \pred construct.

III. ARRAYS IN PRASPEL

In PHP, an array is always an associative array (or map, or dictionary), i.e. a collection of key-value pairs, where each key appears at most once. Keys can be null, booleans (casted into integers), integers, floats (reduced to their integer parts) or strings. Values can be of many kinds. An array can be homogeneous or heterogeneous. In an homogeneous array all the keys have the same type, and all the values too. In an heterogeneous array keys may have distinct types, and/or values may have distinct types. Keys can be auto-incremented, by adding 1 to the last integer index starting by 0. The length (or size) of an array is its number of elements. An array has no predefined length, but its length (stored internally by the PHP engine) can be retrieved thanks to the PHP function count(). An array has also no predefined depth, i.e. it can contain arbitrary arrays.

A. Array Description

In Praspel, array(D,L) denotes the realistic domain of arrays whose domains and codomains are described by D and whose length is in the disjunction L of realistic domains of non-negative integers. D is a comma-separated list, between [ and ], of array descriptions of the form from K to V , where K and V are realistic domain disjunctions, respectively for keys and values. When the from keyword is missing, it is transformed into a realistic domain representing an autoincremented integer starting at 0 with a step of 1.

Example 3 (Homogeneous and heterogeneous arrays). This syntax is illustrated by the following array declarations: The identifier a1 is declared as a homogeneous array of booleans with a length between 7 and 42. This length is a realistic domain that implements the Interval interface. The identifier a2 is declared as a homogeneous array of length 7, whose keys are integers between 0 and 5 or simply 10, and whose values are integers. Its length is a realistic domain that implements the Constant interface, the domain for its keys is the disjunction of two realistic domains (0..5 and 10..23) which implement the Interval interface, etc. The identifiers a3 and a4 are declared as heterogeneous arrays. Both arrays can contain the pairs (5, true) and (15, 4.2), but a4 can contain the pair (5, 4.2), whereas a3 can not contain it.

Actually, we introduce a normal form that removes disjunctions in array descriptions (in from . . . to . . . constructs). An array description is in normal form when it can not be reduced by the rewriting rule (from

F 1 or F 2 to T 1 or T 2 → from F 1 to T 1 , from F 1 to T 2 , from F 2 to T 1 , from F 2 to T 2 ).
Example 4 (Array description in normal form). The following declaration of a4 is in normal form:

B. Collecting Information

When generating test data, the Praspel testing tool calls the sample() method on a realistic domain and then checks whether the predicates declared with \pred() hold. For an array with conditions, a random generator may produce a lot of rejected data before getting a valid one. The idea is to determine the most popular conditions on arrays expressed in PHP (usually written in the \pred() construct), to allow them inside Praspel, and to use a solver to satisfy these conditions.

To achieve this, we have selected 61 PHP projects, from Github and SourceForge, for their popularity, impact on the industry and complexity. All these projects represent 28 066 files and 5 220 547 lines of code. In this code we count the number of occurrences of each array function available in the PHP standard library. It appeared that the three most used functions are: count(), array_key_exists() and in_array(). The count() function counts the number of values in an array, the array_key_exists() function checks whether a key is present in an array (independently of its associated value, e.g. it returns true even if the value is null), and finally, the in_array() function checks whether a value is present in an array. All these functions work on one array at a time. This study suggests that we could consider these side-effect-free Boolean functions as the most frequent conditions on arrays.

C. Array Conditions

We extend the syntax of the array declaration a : array(D,L) in Praspel with the following conditions on arrays.

A pair condition is of the form a[K]: V where K and V are realistic domain disjunctions. The condition means that the pairs constituted of all the keys in K and at least one value in V are present in the array a. K only accepts realistic domains that implements the Constant, Interval and Enumerable interfaces. This is equivalent to use the array_key_exists() and in_array() functions combined.

If we would like to express a constraint only on K, we can use the symbol _. The condition a[K]: _ means that all the keys from K must be present in the array a. It is equivalent to use only the array_key_exists function with all values in K in conjunction.

The condition a[_]: V means that all the values in V must be present in the array a. It is equivalent to use the in_array function with all the values in V in conjunction.

Instead of the : symbol, we can use the symbol !: to express a negation. The condition a[K]!: V means that all the keys in K have a value in the array a and that this value is not in V . It works similarly with the symbol _. For example, the condition a[K]!: _ means that no key in K appears in a.

Keys of an array are always unique, but not its values. We can express a unicity condition on values by writing the condition a is unique. In this case, we cannot have the same value twice in the array a.

Example 5 (Array Conditions).

To illustrate all kinds of conditions, we will use the following example that uses a:

IV. CONSTRAINT SOLVER

Given a conjunction of array conditions on an array a, we propose to invoke a constraint solver to construct an array satisfying all these conditions. This section explains how array conditions are transformed into constraints for the solver. One of these conditions is assumed to be an array declaration of the form a: array(D,L), where D is assumed to be in normal form, without loss of generality. In other words, D is a list of p constructs from F i to T i with 1 ≤ i ≤ p. We also assume that L is L 1 or . . . or L m (with m ≥ 1), where L 1 , . . . , L m are realistic domains that inherit from the Integer realistic domain and that are non-negatives.

In Example 5, p = 1, m = 2, L 1 = [0..5] and L 2 = {10}. In the array description, no domain is declared. In this case the default realistic domain integerpp(0,1) is used (an autoincremented integer: 0, 1, 2, 3, etc). In this example F 1 = integerpp(0,1) and T 1 = string('a','e',1).

Without risk of confusion, a domain disjunction D 1 or . . . or D n will often be identified with the set

D 1 ∪ • • • ∪ D n .

A. Variables

The constraint variables are: (i) the array size -noted Swhich is a non-negative integer, (ii) the sets X and Y , which respectively are the array domain (set of keys) and codomain (set of values), (iii) the array content noted H, which is a total function from X to Y , since keys are unique in a PHP array, (iv) the realistic domains1 X 1 , . . . , X p (resp. Y 1 , . . . , Y p ), which are subsets of the realistic domains F 1 , . . . , F p (resp. T 1 , . . . , T p ) compatible with all the array conditions.

We are essentially interested in finding the content of X and the values of the function H, i.e. the content of H considered as a hashtable, possibly also the values of the X i s and Y i s for checking purposes. The other variables are introduced only to simplify the expression of constraints.

When x ∈ X holds, H(x) = y means that the key-value pair (x, y) is in the array. We extend H to subsets of X by the function Ĥ defined by Ĥ(E) = {H(x) s.t. x ∈ E} for any subset E of X.

B. Cardinality Constraints

Let card(E) denote the cardinality of the finite set E. The constraints card(X) = S and S ≥ 0 say that the array size is its number of keys and is non-negative.

By default, there is no unicity constraint on the codomain, so we only have the constraint card(Y ) ≤ card(X) However, in presence of the array condition a is unique, this constraint becomes card(Y ) = card(X).

C. Constraints on the Array Size

When propagating constraints, the solver may refine the domains L 1 , . . . , L m of possible values for the array size S for i in {1, . . . , m}, and the array size S should be in one of these sets, i.e. the constraint to define domain of S is:

S ∈ L 1 ∪ • • • ∪ L m In Example 5, we have L 1 ⊆ [0..5] and L 2 ⊆ {10}. The size S is constrained by S ∈ L 1 ∪ L 2 .

D. Constraints on Domains and Codomains

The domain X and codomain Y of H are related by the constraint Y = Ĥ(X).

We expect that the constraint solver proposes us the array domain X (resp. codomain Y ) as a disjunction X 1 or . . . or X p (resp. Y 1 or . . . or Y p ) of realistic domains compatible with all the array conditions. We should have the equalities X = 1≤i≤p X i and Y = 1≤i≤p Y i , and the inclusions:

X i ⊆ F i and Y i ⊆ T i for i in {1, . . . , p}. The pair (X i , Y i ) with 1 ≤ i ≤ p should also satisfy the constraint Ĥ(X i ) = Y i meaning that Y i is the codomain of the restriction of H to X i (⊆ X).

E. Constraints on Pairs

For each array condition a[K]: V where K and V are domain disjunctions, we introduce the constraints:

K ⊆ X and Ĥ(K) ⊆ V . A negated pair condition a[K]!: V is translated into the constraints: K ⊆ X and Ĥ(K) ∩ V = ∅.

F. Constraints on Keys or Values

The condition a[K]: _ is translated into the constraint K ⊆ X, and its negation a

[K]!: _ into the constraint K ∩ X = ∅. The condition a[_]: V is translated into the constraint V ⊆ Y , and its negation a[_]!: V into the constraint V ∩ Y = ∅.

G. Propagation and consistency

Propagation of constraints uses an AC3 algorithm [START_REF] Macworth | Consistency in network of relations[END_REF] implemented in PHP. So, we use five kinds of domains associated to five kinds of realistic domains: Constant, Interval, Nonconvex, Finite and Enumerable (see Section II-B2). For each kind of domain, we have implemented a revise method to allow the domain reduction. So, the consistency is also checking that there is no empty domain for the four variables S, H, X and Y but not for X i and Y i . The goal is to detect inconsistencies as soon as possible.

H. Labelling

The labelling is the process of finding a value for each variable. In order to make the solver converge quickly to a solution, we use a heuristic that consists in chosing a value for the variable S at first. This helps to unfold the ∀ and ∃ quantifiers (because F i and T i are enumerables, we manipulate finite sets). Then, the solver tries to compute the sets X i and Y i . We use a random generator to generate a value in a realistic domain, to select a realistic domain in a disjunction, etc. The generated value is then propagated. If an inconsistency is detected, we add a new constraint to discredit the value, and then generate another one. For instance, if S = 5 leads to an inconsistency, we add the constraint S = 5. The added constraint is removed during the backtracking step.

When all variables are labelled, i.e. each one has a valid value, the solver returns the solution.

V. EXPERIMENTATION

In all that follows, the word system stands for the expression "conjunction of array conditions". This section presents an experimentation evaluating the solver efficiency, i.e. its capability to avoid or reduce rejection when generating data from systems. We measure the number of backtracks in the solver, the time to generate data from satisfiable systems of array conditions, and how many unsatisfiable systems are detected. The experimentation is composed of three steps: system generation, then data generation (i.e. system solving) and finally a measuring step. We generate systems on arrays containing strings and integers, and of length 5 to 20.

The Praspel language (and in particular its sublanguage of array conditions) is described by a grammar. In order to generate systems, we re-use a previous work in grammar-based testing [START_REF] Enderlin | Grammar-Based Testing Using Realistic Domains in PHP[END_REF] proposing three algorithms generating data from grammars: a uniform random generator, a bounded exhaustive test generator, and a rule-coverage-based test generator. Since the grammar of array conditions is small, the last generator does not generate a sufficiently wide collection of systems. Bounded exhaustive testing is more costly than random testing, but it is more precise and well adapted to small grammars. We retain the bounded exhaustive test generator: for increasing values of n it enumerates all the systems composed of a sequence of n tokens: In this generator, a single token value of each variable token (i.e. interval bounds and particular array lengths, keys and values) is generated at random.

With n = 3, we can generate conditions of the form arr is unique. With n = 6, we can generate constraints of the form arr[0]: 0. With n = 8, we can generate constraints of the form arr[0]: 0 or 1. With n = 11, we can generate for instance the system in Example 5.

The second step calls the solver with each produced system to generate an array satisfying it. Every generated array is evaluated by the predicate associated to the system of array conditions, to check the solver soundness. During the data generation step, some marks are placed to count only the array generation time, without counting the compiling time. We also measure the number of backtracks in the solver and the number of rejected systems. When a system is rejected, its backtracks are not counted.

Since our grammar-based testing algorithms use an isotropic random generation to generate token values, these values can differ from one system to another. It may lead to different generation times and numbers of backtracks. To avoid peaks in the results, we report the average of 100 generations of systems sharing the same pattern.

Table I shows our experimentation results. In the first column, n is the length of the generated sequences of tokens representing systems of array conditions. Column 2 gives the number of distinct system patterns with this length. Columns 3, 5 and 6 respectively give average numbers of backtracks and rejected systems for 100 runs, and an average generation time. Column 4 gives the rate of backtracks per system pattern.

For n ≤ 20, we observe that no system is rejected, which is a great improvement by comparison with random generation. All generated data satisfy their specification. The solver successfully and quickly generates data with a low number of backtracks. For n = 20, which represents approximatively 3 constraints with disjunctions, the exhaustive generation algorithm generates 492 distinct systems of array conditions and the execution produces only 114 backtracks, so approximatically 1 backtrack for 4 systems. This is a good result. Nevertheless, we were not able to characterize the number of backtracks with the number of constraints. Finally, this experimentation allowed us to find a bug in our solver. This bug always led to a rejection when only a certain constraint was analyzed. Thus, it shows that the validation process also helps finding bugs.

VI. RELATED WORKS

Various works consider Design-by-Contract for unit test generation [START_REF] Cheon | A Simple and Practical Approach to Unit Testing: The JML and JUnit Way[END_REF], [START_REF] Madsen | Unit Testing using Design by Contract and Equivalence Partitions[END_REF]. Our specification language is based on JML [START_REF] Leavens | JML: A notation for detailed design[END_REF] and ACSL [START_REF] Baudin | ACSL: ANSI C Specification Language[END_REF]. Thanks to an expressive specification language, Praspel performs general runtime assertion checks. Realistic domains present some similarities with Eiffel's types [START_REF]Eiffel: programming for reusability and extendibility[END_REF], especially regarding inheritance between realistic domains. Nevertheless, the two properties of predicability and samplability displayed by realistic domains do not exist in Eiffel. Moreover, Praspel adds clauses that Eiffel contracts do not support, as @throwable and @behavior, which are inspired from JML.

Euclide [START_REF] Gotlieb | Euclide: A Constraint-Based Testing Framework for Critical C Programs[END_REF] is a constraint-based testing tool that could take as input additional safety properties defined in ACSL [START_REF]Tcas software verification using constraint programming[END_REF]. Our approach differs by handling conditions directly in the contract. All constructions present in Praspel are well-handled for both aspects: validation and generation. The CLP framework INKA [START_REF] Gotlieb | A CLP Framework for Computing Structural Test Data[END_REF] helps computing structural test data from a C program. It transforms the problem of automatic test data generation into a CLP problem over finite domains. In the same way, FDCC [START_REF] Bardin | FDCC: A Combined Approach for Solving Constraints over Finite Domains and Arrays[END_REF] is a combined approach for solving constraints over finite domains and arrays. The tricky part of FDCC lies in a bi-directional communication mechanism between two solvers. Our constraints are more specific but we use only one solver.

VII. CONCLUSION AND FUTURE WORKS

We have presented in this paper an extension of the Praspel language to specify usual conditions on PHP arrays. We have expressed its semantics by constraints. We have designed and implemented in PHP a constraint solver to generate test data from these constraints. It uses a random generator to ensure a diversity of generated solutions. This solver is integrated in realistic domains and can be used within the Praspel framework. A first validation shows cases where rejection has been totally removed. It also shows that the solver dramatically increases the generation speed.

In a near future, we plan to lead a more complete experimentation. Then, we plan to formalize more constraints and extend our solver. We also plan to transform constraints into the formalism proposed by MiniZinc [START_REF] Nethercote | Minizinc: Towards a standard cp modelling language[END_REF] in order to compare our solver to other ones regarding performances and capabilities to find solutions. As a string is an array of characters, we would like to apply the same process on strings with the help of existing and promising solving techniques [START_REF] Ganesh | Hampi: A string solver for testing, analysis and vulnerability detection[END_REF].

  a1: array([to boolean()], 7..42) a2: array([from 0..5 or 10 to integer()], 7) a3: array([from 0..10 to boolean(), from 20..30 to float()], 7) a4: array([from 0..10 or 20..30 to boolean() or float()], 7)

  a4: array([from 0..10 to boolean(), from 0..10 to float(), from 20..30 to boolean(), from 20..30 to float()], 7)

  length: 0..5 or 10 a : array([to string('a', 'e', 1)], length) a[0]: 'b' or 'd' a is unique

  For the condition a[0]: 'b' or 'd' in Example 5, we have K = {0} and V = {'b', 'd'}. The constraints are {0} ⊆ X and Ĥ({0}) ⊆ {'b', 'd'}.

  class C { / ** @invariant I 1 and . . . and I h * / / ** @requires R 1 and . . . and Rn; * @ensures E 1 and . . . and E j ; * / function foo ( $x 1 . . . ) { body } } Figure 1. Syntax of contracts in Praspel Test generation in Praspel is decomposed into two steps. First, a test generator computes test data from contracts. Second, a dedicated test execution framework runs the test cases (i.e. invokes the methods with the computed test data, and checks the assertions at run time) so as to establish the test verdict.

In fact, the solver will not handle realistic domains, but only sets.


Ph.D Giovanni De Nunzio
email: giovanni.de-nunzio@ifpen.fr

Carlos Canudas De Wit
email: carlos.canudas-de-wit@gipsa-lab.fr

Philippe Moulin
email: philippe.moulin@ifpen.fr

Domenico Di Domenico
email: domenico.didomenico@ifpen.fr

Eco-Driving in Urban Traffic Networks using Traffic Signal Information

Keywords: Energy minimization, velocity planning, traffic lights, graph. I

This work addresses the problem of finding energy-optimal velocity profiles for a vehicle in an urban traffic network. Assuming communication between infrastructure and vehicles (I2V) and a complete knowledge of the upcoming traffic lights timings, a preliminary velocity pruning algorithm is proposed in order to identify the feasible region a vehicle may travel along in compliance with city speed limits. Then, a graph discretizing approach is utilized for advanced selection, among the feasible "green windows", of the optimal ones in terms of energy consumption. Finally, a velocity trajectory is advised, which will be tracked by the driver-in-the-loop in order to pass through the signalized intersections without stopping. The proposed eco-driving assistance algorithm results are compared to the optimal solution provided by the Dynamic Programming, in order to prove not only the effectiveness but also its capability to be employed online due to its low computational load.

vehicles (i.e. eco-driving) in urban traffic. Green lights optimal speed advisory (GLOSA) service has been simulated, also varying the penetration rate of the technology, showing promising advantages as compared to ordinary vehicles in terms of fuel consumption and average stop time [START_REF] Katsaros | Performance study of a Green Light Optimized Speed Advisory (GLOSA) application using an integrated cooperative ITS simulation platform[END_REF]. Vehicle energy consumption is strictly related to acceleration driving profiles. Soft decelerations to traffic lights on red, using advanced wireless transmitted information, have been proven to be more energy efficient than sudden stops, allowing a reduction in both fuel consumption and emissions [START_REF] Li | Traffic Energy and Emission Reductions at Signalized Intersections: A Study of the Benefits of Advanced Driver Information[END_REF]. Other works show how it is possible to reduce energy consumption by preventing the vehicle from coming to a full stop at the intersections and by advising cruising velocities in order to catch as many traffic lights on green as possible. In [START_REF] Mandava | Arterial velocity planning based on traffic signal information under light traffic conditions[END_REF], authors developed an algorithm to minimize the acceleration profile while driving through a sequence of signalized intersections, simulating scenarios with stochastically varied parameters and showing a systematic reduction in fuel consumption, emissions and travel time. Interesting results on eco-driving with probabilistic signal phase and timing (SPAT) information are shown in [START_REF] Mahler | Reducing idling at red lights based on probabilistic prediction of traffic signal timings[END_REF], where a comparison between an uninformed driver and a driver with different information levels (full horizon and short-term) is conducted using Dynamic Programming (DP). An innovative solution to the velocity advisory problem, in proximity to a traffic light, was proposed in [START_REF] Koukoumidis | Leveraging Smartphone Cameras for Collaborative Road Advisories[END_REF], where smartphones equipped with cameras mounted on the windshield of the vehicles allow to build a SPAT map of every intersection by relying only on V2V communication, without directly interacting with the infrastructure. This approach enables GLOSA service and adaptive route change (ARC) among the others. Predictive cruise control (PCC) in traffic networks with signalized intersections was extensively treated in [START_REF] Asadi | Predictive Cruise Control: Utilizing Upcoming Traffic Signal Information for Improving Fuel Economy and Reducing Trip Time[END_REF], where a preliminary logic computes the constant velocity to be tracked by the driver to pass through the as high number of intersections as possible, while penalizing the brakes action to indirectly reduce energy consumption and trip time. In [START_REF] Rakha | Eco-driving at signalized intersections using V2I communication[END_REF] only one traffic light is considered and accelerations upstream and downstream of the intersection are varied to find sub-optimal values for the fuel consumption minimization problem. In [START_REF] Miyatake | Theoretical Study on Eco-Driving Technique for an Electric Vehicle Considering Traffic Signals[END_REF] Dynamic Programming is used to find the minimum energy solution, for a trip with up to two traffic lights and investigating different timings of the green signal.

In the case of long or multiple available "green windows" at each intersection, different driving profiles are possible and an optimal route decision should be made at the beginning of the driving horizon to keep the energy consumption at a minimum. The novelty of the proposed approach lies in the capability of identifying in advance a low cost path (eco-path), among all the possible trajectories within the driving horizon, without having to stop at traffic lights, without penalizing the trip time and while keeping a low computational load. The described algorithm provides a fast sub-optimal solution for the energy minimization problem, which is promising for an online speed advisory service to the driver, and it is compared to the optimal solution provided offline by Dynamic Programming (benchmark). Note that unexpected additional constraints coming from the network (vehicles, pedestrians, etc.) may be addressed with new executions of the algorithm.

Section II describes the problem and its elements from a mathematical point of view. In Section III the applied methodology for the solution of the problem is explained in detail. Section IV contains simulation results and comparisons of the algorithm against the benchmark. Conclusions and future developments in Section V.

II. PROBLEM FORMULATION

From a mathematical point of view, a constrained optimization problem needs to be solved. The function to be minimized is the energy consumption over the trip horizon, subject to the vehicle dynamic model, with constraints on the crossing times of the intersections, allowed only when the relative traffic light state is on "green".

A. Vehicle model

Vehicles in our scenario are assumed to be electric. The equivalent circuit model of a DC motor gives [START_REF] Petit | Optimal drive of electric vehicles using an inversion-based trajectory generation approach[END_REF]:

V a = i a R a + e e = κω m (1) T m = κi a
where V a is the armature voltage, i a is the armature current, R a is the armature resistance, e is the back electromotive force, κ is the speed constant, ω m is the rotational speed of the motor, T m is the motor torque. The electric power is calculated by expressing:

P m = V a i a = ω m T m + R a κ 2 T 2 m (2)
In this equation V a i a , that is the electric power supplied to the armature, is converted into actual mechanical power developed by the armature (ω m T m ), and electric power wasted in the armature ( Ra κ2 T 2 m ). The control variable u is the demanded motor torque:

T m = u (3)
The longitudinal dynamic model is:

m dv dt = F t -F aero -F friction -F slope ( 4
)
where F t is the traction force at the wheels, F aero is the aerodynamic force, F friction is the rolling resistance force, F slope is the gravity force.

Under the hypotheses of no losses in the transmission and no slip at wheels, the vehicle model shall be written as:

ẋ = v m v = uRt r -1 2 ρ a A f c d v 2 -mgc r -mg sin(α(x)) (5
) where R t is the transmission ratio, r is the wheel radius, m is the vehicle mass, ρ a is the external air density, A f is the vehicle frontal area, c d is the aerodynamic drag coefficient, c r is the rolling resistance coefficient, α is the road slope as a function of the position x, and g is the gravity.

The sum of the aerodynamic and rolling frictions can be approximated as a second order polynomial:

F res = F aero + F friction = a 0 + a 1 v + a 2 v 2 (6)
where a 0 , a 1 and a 2 are identified parameters [START_REF] Dib | Optimal Energy Management Compliant with Online Requirements for an Electric Vehicle in Eco-Driving Applications[END_REF]. Under these assumptions the vehicle model can be simplified as follows:

ẋ = v v = h 1 u -h 2 v 2 -h 3 v -h 0 (7)
with

h 1 = R t mr , h 2 = a 2 m , h 3 = a 1 m , h 0 = a 0 m + g sin(α) (8)
Then the power consumption can be expressed as:

P = b 1 uv + b 2 u 2 (9)
where

b 1 = R t r , b 2 = R a κ 2 (10)

B. Traffic lights timing

The control of the intersection signal timings is not considered here, therefore cycle, stage and offset times are deterministic and given. Under the assumption that traffic lights can share with other traffic agents their signal timings information (I2V communication), it is possible to formulate mathematically the state evolution of the traffic lights as follows:

s i (t) = 1, (k -1)T < t -l i (k -1)T + T gr 0, (k -1)T + T gr < t -l i kT (11)
l i = T 0 + di v , T 0 + di v < T T 0 + di v -kT, kT < T 0 + di v < (k + 1)T (12)
where s i is the state of the i-th traffic light, T 0 ∈ [0, T ] is the initial time to be chosen for the cycles, k is the number of cycles, T is the time duration of a cycle, T gr is the duration of the green signal, d i is the position of the i-th traffic light, v is the "green wave" velocity, k ∈ Z, i ∈ [1, . . . , n], being n the number of intersections. This is explained as follows:

• the time (d i /v) to reach each intersection traveling at the velocity v is calculated; • T 0 , chosen randomly in the interval [0, T ], is added to the "green wave" condition just calculated in order to generate a random starting point in the signal cycle; • hence we obtain the initial "time offset" (l i ) of the lights, which is scaled down to the interval [0, T ] if it exceeds the duration T of the cycle; • as the simulation time t runs, the state of the lights is updated according to the cycle and split times.

C. Optimal control problem

Finally, the optimization problem may be stated as follows:

                                       min u J = t f t0 b 1 uv + b 2 u 2 dt s.t. ẋ = v v = h 1 u -h 2 v 2 -h 3 v -h 0 x(t 0 ) = x 0 , x(t f ) = D x(t i ) = d i ∧ s i (t i ) = 1 v(t 0 ) = v 0 , v(t f ) = v f v min ≤ v ≤ v max u min ≤ u ≤ u max (13)
Nature and complexity of the constraints ( 14)

x(t i ) = d i ∧ s i (t i ) = 1 (14)
which may result in disjoint sets due to the presence of multiple available "green windows" at the same intersection, affect the convexity of the problem and a global analytical solution is not achievable. In other words the constraints ( 14) represent a non-convex set and the nonlinear objective function, because of such a discontinuity in the constraints, assumes multiple local minima. Therefore, the solution to this optimal control problem has to be sought in a suboptimal way, making use of a set of algorithms to simplify the control envelop, recover the convexity of the constraints set and advise the driver on the velocity to track.

III. OPTIMIZATION ALGORITHM

The original optimization problem (13) will be now divided into sub-problems for the simplification of constraints [START_REF] Rakha | Eco-driving at signalized intersections using V2I communication[END_REF] and the formulation of a convex optimization problem.

The methodology may be summarized in three main steps: 1) Pruning algorithm to identify the set of feasible "green windows" and their portions which allow not to stop at the intersections in accordance with city speed limits

(v ∈ [v min , v max ]);
2) Construction of a weighted directed acyclic graph (nodes at each feasible "green window") to analyze the energy cost of all the possible paths in the driving horizon, and identification, by means of Dijkstra's algorithm, of the "shortest" (in terms of energy) path; 3) Solution of an optimization problem to obtain the optimal crossing points at each intersection through the previously selected "green windows".

A. Pruning algorithm

The velocity pruning algorithm reduces the number of available "green windows" by identifying only the feasible ones that allow to reach a horizon at a specified final time without stopping. Let t i be the crossing time at the i-th intersection, (t 0 , d 0 ) the initial position of the vehicle, d i the position of the i-th traffic light, D the horizon.

Algorithm 1 Pruning algorithm 1: for i ← 1 to n do 2: t i,min = di-di-1 vmax + t i-1,min ; 3: t i,max = di-di-1 vmin + t i-
if t i,max ≤ t i-1,max || di-di-1
ti,max-ti-1,max > v max then 14:

t i-1,max = t i,max -di-di-1 vmax ; 15:
end if 16: end for The algorithm will output [t i,min , t i,max ], which represent respectively the minimum and the maximum feasible crossing times at the i-th intersection, respecting speed limits, and its steps may be described as follows:

1 In case no feasible region is found, it will not be possible to avoid stopping at a traffic light along the trip, but the algorithm may be run again at the next "green" signal with respect to the moving horizon.

)

B. Optimal path

After running Algorithm 1, many possible paths are still present and a decision for the optimal one needs to be made.

The idea is to approximate, by means of a weighted directed acyclic graph (DAG), the energy cost of all the possible paths from the initial to the final point within the feasibility region shown in Fig. 2.

Therefore, let G = (V, E) be such a graph, where V is the set of vertices representing the crossing times at each feasible "green window", and E is the set of feasible paths connecting the vertices of the graph. Let us also define a weighting function w : E → W which associates each edge of the graph with a weight. In this application the weight W is the energy cost to travel along the path and it may be seen as the sum of the energy contribution given by a constant speed trip on an edge, and the contribution given by the change of velocity at a node between two edges.

W = E total = E link + E jump ( 15
)
E link = t+∆t t b 1 ūv + b 2 ū2 dt = ∆t b 1 ūv + b 2 ū2 (16)
where v is a constant traveling velocity on an edge, and from [START_REF] Mandava | Arterial velocity planning based on traffic signal information under light traffic conditions[END_REF],

for v = 0, ū = 1 h1 h 2 v2 + h 3 v + h 0 ; E jump = tjump 0 b 1 uv + b 2 u 2 dt ( 17
)
where

t jump = |vi+1-vi| a
, a is a constant fixed acceleration linearly modeling the velocity in every transient, v = v i ±a•t for i = 0, . . . , n + 1, with v 0 being the initial velocity and v n+2 = v f the final velocity.

Such graph approximation introduces a difficulty when it comes to assigning the velocity jump energy contribution to the edges. The transient energy contribution on an edge depends on the edge previously traveled: if coming from Fig. 3: Tree expansion of the possible paths from the root to the destination node.

a higher speed edge the transient contribution will be very small, on the other hand if coming from a lower speed edge the transient contribution will be much bigger due to the increase of velocity. Every node of the graph with more than one incoming edge is critical and needs to be "decoupled" in order to have a correct weight assignment, since the outgoing edges of the critical nodes can have only one weight. The decoupling of the graph into the tree of all the possible paths is shown in Fig. 3. It was necessary to create duplicates of the critical nodes in order to have independent paths, and in the case under analysis the initial setting with 14 nodes (Fig. 2) results in a bigger network (Fig. 3) with the same amount of possible paths from node 1 to node 14 but with 52 nodes. Well known theorems in graph theory allow to analyze a priori the complexity of this procedure (see p.165 in [START_REF] Godsil | Algebraic Graph Theory[END_REF]). Specifically, being A the adjacency matrix of a graph G = (V, E), the number of paths of length r from a node a to a node b is equal to (A r ) a,b . The number of nodes of the expanded decoupled tree is given by the sum of the entries of the first row of the geometric series of the adjacency matrix up to the length of the trip r. The geometric series generated by the matrix A, if

|λ i | ≤ 1 for each eigenvalue |λ i | of A, is known to converge to r k=0 A k = (I -A) -1 (I -A r+1 ) (18)
The problem has an exponential complexity and its computational load might increase significantly but, on a realistic horizon of few intersections with not numerous "green windows" in the feasible region, the number of nodes of the decoupled tree can be kept low.

To transform the decoupled tree in a graph, a fictitious destination node connected to all the leaves is added. On the expanded weighted DAG, Dijkstra's algorithm is used to find the minimum energy trip from the origin to the destination.

C. Simplified optimization problem

The aim of the described simplification steps is to have a convex, also in the constraints, optimization problem to find the optimal crossing times within the "green windows" selected by Dijkstra's algorithm. The formulation of this optimization problem may be carried out in alternative ways.

1) Velocity as variable: In order to formulate the optimization problem using velocity as variable, let us define the vector of optimization variables (velocities to track in each segment), assuming n to be the number of intersections:

X = [v 1 , v 2 , . . . , v n , v n+1 ] T ∈ R n+1 (19)
Then a vector of time intervals for each segment is defined:

T = T(X) = [t 1 -t 0 , t 2 -t 1 , . . . , t n -t n-1 , t f -t n ] T = d 1 -d 0 v 1 , d 2 -d 1 v 2 , . . . , d n -d n-1 v n , D -d n v n+1 T ( 20
)
where d i is the known position of the i-th traffic light, D is the destination/horizon, and T ∈ R n+1 . The objective function is the same as the one in the original optimization problem (13):

J = t f t0 P dt = t f t0 b 1 uv + b 2 u 2 dt (21)
By using the same procedure as in the graph weights assignment, we may split the objective function in an energy term related to the contributions of the segments traveled at constant velocity, and energy terms related to the velocity variations between subsequent segments. Therefore, the objective function may be rewritten as:

J(X) = T T P (X) + n+1 i=0 E jump ( 22
)
where P (X) is the power required by the vehicle when traveling at constant velocity (no velocity variation in time):

P (X) = b 1 ū(X)X + b 2 ū(X) 2 (23) with ū(X) = 1 h 1 h 2 X 2 + h 3 X + h 0 (24)
Finally the optimization problem may be formulated as:

min X J(X) v min ≤ v - i (X) ≤ v i ≤ v + i (X) ≤ v max (25)
Note that in this case, the constraints are nonlinear. Specifically the bounds of the constant velocity to be tracked in each segment (v i ) depend on the crossing time at the previous intersection (t i-1 ) and experimentally they are found to be approximated by cubic functions.

2) Time as variable: In order to overcome the difficulty represented by the nonlinear constraints of the previous formulation, it is possible to formalize the problem in an analogous way which allows to deal with constant constraints. Let us first define a vector of optimization variables (crossing times at intersections), supposing n to be the number of intersections:

X = [t 1 , t 2 , . . . , t n ] T ∈ R n (26)
Then a vector of time intervals for each segment is defined:

T = T(X) = [t 1 -t 0 , t 2 -t 1 , . . . , t n -t n-1 , t f -t n ]
T (27) Knowing the position of the traffic lights, let us then define the vector of the constant velocities in each segment:

V = V(X) = [v 1 , v 2 , . . . , v n , v n+1 ] T = d 1 -d 0 t 1 -t 0 , d 2 -d 1 t 2 -t 1 , . . . , d n -d n-1 t n -t n-1 , D -d n t f -t n T (28)
where d i is the position of the i-th traffic light with respect to the origin, D is the destination, t f is the horizon time, and

T ∈ R n+1 , V ∈ R n+1 .
The objective function is again the same as the one in the original optimization problem (13):

J = t f t0 P dt = t f t0 b 1 uv + b 2 u 2 dt (29)
It may be rewritten as:

J(X) = T T P (V) + n+1 i=0 E jump (30)
where P (V) is the power required by the vehicle when traveling at constant velocity (no velocity variation in time):

P (V) = b 1 ū(V)V + b 2 ū(V) 2 (31)
with

ū(V) = 1 h 1 h 2 V 2 + h 3 V + h 0 (32)
Finally the optimization problem may be formulated as:

min X J(X) t - i ≤ t i ≤ t + i ( 33
)
where t - i and t + i are constants and represent the ends of the selected green split time at each intersection.

IV. SIMULATION RESULTS

In the following, performance of the proposed algorithm is compared to the solution provided by DP [START_REF] Sundstrom | A Generic Dynamic Programming Matlab Function[END_REF]. Path number Optimal windows v 0 = 10 m/s v 0 = 11 m/s v 0 = 12 m/s v 0 = 13 m/s v 0 = 14 m/s 1 1-2-4-8-10-12-14 459900 J 446800 J 429200 J 410700 J 394200 J 2 1-2-4-8-10-13-14 436700 J 423400 J 406000 J 387800 J 371500 J 3 1-2-5-8-10-12-14 454700 J 438200 J 421600 J 406700 J 397400 J 4 1-3-5-8-10-12-14 450700 J 449700 J 449100 J 449000 J 449100 J 5 1-3-5-8-10-13-14 435800 J 434900 J 434600 J 434800 J 435400 J 6 1-3-5-8-11-13-14 417800 J 418100 J 417700 J 417400 J 418000 J 7 1-3-5-9-11-13-14 420400 J 420000 J 419800 J 419600 J 420400 J 8 1-2-5-9-11-13-14 425100 J 408600 J 392200 J 376900 J 367500 J 9

1-2-4-8-11-13-14 423600 J 410400 J 393200 J 375200 J 359100 J 10 1-2-4-7-10-13-14 427300 J 412700 J 394000 J 378500 J 362300 J 11 1-2-5-8-11-13-14 415100 J 398700 J 382100 J 367300 J 357800 J 12 1-2-4-7-10-12-14 438600 J 423800 J 404400 J 389100 J 373000 J 131-3-6-9-11-13-14 413600 J 411800 J 411700 J 411700 J 411800 J 14 1-2-5-8-10-13-14 438300 J 421700 J 405100 J 390200 J 380900 J

A. Control scenario

The case under study presents a simplified scenario with one vehicle in the traffic network, traveling through n = 5 perpendicular signalized intersections, along a road section of 2000 m on a total time horizon of 200 s. The vehicle's objective is to follow an energy optimal velocity profile which allows to find all the traffic lights on green. The traffic lights obey a known signal timing policy, as specified in [START_REF] Sundstrom | A Generic Dynamic Programming Matlab Function[END_REF], and their states are known to the vehicle. The green split time is 10 s, the total cycle time is set to 30 s. Note that, although traffic lights obey a "'green wave" condition, the distribution of the green windows may be such that the velocity of the wave is not energy optimal.

The initial velocity of the vehicle is varied throughout the simulations to assess its impact on the energy consumption and on the optimal path. Traffic lights timings, final velocity, horizon are fixed for consistency in results comparison. The objective of the performed simulations is to prove the importance of a methodology which allows to achieve a fast sub-optimal solution suitable for online implementation and eco-driving assistance service.

B. Optimal path

The first part of the presented simulation results aims at validating the pruning algorithm and the optimal path identification. DP, used in the following as a benchmark, provides the optimal solution to the original problem and it has been used to compute the energy cost of all the possible paths in the feasible region.

Simulations for this first validation part have been run varying the initial velocity of the trip, which has an impact on the optimal path. In the scenario under study the vehicle does not stop at any intersection in the feasibility region (Fig. 2), nevertheless the energy consumption varies depending on the path choice; in the worst case, the choice of a path rather than the optimal one, results to be 10% more expensive in terms of energy.

Table I summarizes the trip energy consumption of all the possible paths for different choices of initial velocity, and the minimum cost path, for each case, is highlighted in red. The optimal windows numbering refers to the one utilized in I is represented, with the minimum energy path in green and the unconstrained path in blue (no traffic lights, constant velocity trip along the whole horizon).

The accuracy of the graph approximation, as well as the computational load, increases with the number of nodes for each "green window". A graph with one node per window, as shown in Fig. 2, and a graph with three nodes per window (at the mid-point and at the two ends) have been tested, and the optimal path suggestion obtained with Dijkstra's algorithm is validated against the true minimum energy path, whose energy was computed with the DP.

As summarized in Table III, the graph with 3 nodes per window approximates better the whole energy cost and it fails in two cases, whereas the smaller graph fails in five cases. Note that in this scenario, the algorithm gives correct estimations also for lower initial velocities v 0 ∈ [5, 9], but it was not reported in the table for brevity. This is a suboptimal solution, the path identification is not flawless, but in the worst case, for the graph with 3 nodes/window, the algorithm leads to the selection of a path which is only 2% more expensive than the optimal one.

C. Optimal crossing times

Once the sub-optimal windows are obtained with Dijkstra's algorithm, the optimal crossing times at each intersec- Fig. 5: Sub-optimal crossing times compared to the optimal solution provided by DP, for v 0 = 10 m/s. tion in the horizon need to be calculated in order to provide the driving assistance by suggesting a velocity profile to the driver. The solver employed to solve the optimization problems (25) and (33) in MATLAB R was fmincon, using the interior-point method algorithm. As already mentioned, the optimization problem (25) has nonlinear constraints which need to be redefined at each optimization run. Therefore the solution of the optimization (33), using crossing times as variables, benefits from an easier definition and requires less iterations to converge to the solution. In Fig. 5, the optimal solution provided by DP through the optimal windows (path #13) for v 0 = 10 m/s is compared to the sub-optimal solution. In Table IV the time deviations of the sub-optimal solution from the DP solution at the five intersections are summarized, both for problem (25) and (33).

The simulations were run with a laptop equipped with an Intel(R) Core(TM) i3-2350M at 2.30GHz and 6GB of RAM. DP takes about 3 hours to run a complete simulation and provide the optimal solution. The proposed methodology takes 4.7 seconds (graph 3 nodes/window) or 2.5 seconds (graph 1 node/window) to run all the presented steps and provide the sub-optimal crossing times at each intersection.

V. CONCLUSIONS

Previous works proved how the exchange of information between the infrastructure and the vehicles in the urban traffic (I2V) has a positive impact on the energy consumption. In this work, the possibility of further improvement of energy consumption in case of multiple available nostop paths on a driving horizon is presented. The proposed methodology also allows to simplify a complex constrained optimization problem, and it becomes appealing for online applications. Indeed the proposed algorithm leads to a suboptimal solution and, as shown in the results, may identify wrong optimal paths. However, it was shown that the suboptimal solution is still very close to the optimal one and in terms of energy the difference is minor. Future works will consider also traffic lights with time-varying cycle and stage times, and vehicles embedded in a traffic flow.

Fig. 1 :

 1 Fig. 1: I2V communication in an urban traffic network.

Fig. 2 :

 2 Fig.2: Graph of the possible paths within the feasibility region, with 1 node per "green window" and node numbering.

Fig. 4 :

 4 Fig. 4: Optimal energy cost of all the possible paths from initial to final point, for v 0 = 10 m/s.

Fig. 2 .

 2 Fig.2. In Fig.4, the corresponding column for v 0 = 10 m/s of Table I is represented, with the minimum energy path in green and the unconstrained path in blue (no traffic lights, constant velocity trip along the whole horizon).The accuracy of the graph approximation, as well as the computational load, increases with the number of nodes for each "green window". A graph with one node per window, as shown in Fig.2, and a graph with three nodes per window (at the mid-point and at the two ends) have been tested, and the optimal path suggestion obtained with Dijkstra's algorithm is validated against the true minimum energy path, whose energy was computed with the DP.As summarized in TableIII, the graph with 3 nodes per window approximates better the whole energy cost and it fails in two cases, whereas the smaller graph fails in five cases. Note that in this scenario, the algorithm gives correct estimations also for lower initial velocities v 0 ∈[START_REF] Petit | Optimal drive of electric vehicles using an inversion-based trajectory generation approach[END_REF][START_REF] Mahler | Reducing idling at red lights based on probabilistic prediction of traffic signal timings[END_REF], but it was not reported in the table for brevity. This is a suboptimal solution, the path identification is not flawless, but in the worst case, for the graph with 3 nodes/window, the algorithm leads to the selection of a path which is only 2% more expensive than the optimal one.

  from the starting point calculate [t min , t max ] corresponding to the [v min , v max ] to reach the first intersection d 1 ; 2) check if from the point with coordinates (t max , d 1 ) it is possible to reach the final destination in compliance with the v max constraint, if not set t max accordingly; 3) from d 1 use the two starting points (t min , d 1 ) and (t max , d 1 ) to calculate t min and t max at d 2 ; 4) check again if from the point (t max , d 2 ) it is possible to reach the final destination in compliance with the v max constraint, if not set t max accordingly; 5) repeat steps 3 and 4 for each intersection up to d

n ; 6) at the last intersection d n connect the two points (t min , d n ) and (t max , d n ) to the destination; 7) check backwards from d n to d 1 if any (t max , d i ) ≤ (t max , d i-1 ) or if any t max induces higher velocities than v max , if so move (t max , d i-1 ) accordingly.

TABLE II :

 II Simulation parameters

	Parameter	Description	Value	Unit
	D	horizon	2000	m
	t f	final time	200	s
	v f	final velocity	10	m/s
	v min	minimum speed limit	5	m/s
	vmax	maximum speed limit	14	m/s
	m	vehicle mass	1190	kg
	r	wheel radius	0.2848	m
	Rt	transmission ratio	6.066

TABLE I :

 I Optimal energy consumption on the paths calculated with Dynamic Programming

TABLE III :

 III Optimal path identification (Graph vs. DP)

		Initial velocity 10 m/s	Graph 1 node/window path #13	Graph 3 nodes/window path #13	Optimal path #13
			11 m/s		path #13		path #11		path #11
			12 m/s		path #10		path #11		path #11
			13 m/s		path #10		path #9		path #11
			14 m/s		path #10		path #10		path #11
		2000
				Dynamic Programming Optimum
		1800		Sub-optimum
		1600								6.616		13.51
		1400						14			14
	distance [m]	800 1000 1200			14		7.749	14	14
		600
		400		14			6.863
		200	14		5.455
		0	0	20	40	60	80	100	120	140	160	180	200
							time [s]

TABLE IV :

 IV Time deviation from the optimal crossing times Int. #1 Int. #2 Int. #3 Int. #4 Int. #5

	Opt. TIME	0.4 s	0.03 s	0.5 s	0.16 s	0.8 s
	Opt. VELOCITY	0.4 s	0.03 s	0.5 s	0.16 s	0.8 s
		Average deviation: 0.38 s