Jean-Philippe Metivier

Patrice Boizumault

Samir Loudni

Solving Nurse Rostering Problems Using Soft Global Constraints

published or not. The documents may come   L'archive ouverte pluridisciplinaire


Samir F Matar
email: matar@icmcb-bordeaux.cnrs.fr

In silico CrNF, a half-metallic ferromagnetic nitride-fluoride mimicking CrO 2

Keywords: Magnetism, Half-metals, DFT, bonding

Isoelectronic with CrO2, CrNF is proposed in silico based on rutile derived structure with DFT computations. The ground state structure defined from cohesive energies is of MgUO4-type, characterized by short covalent Cr-N and long ionic Cr-F distances. Like CrO2 it is a half-metallic ferromagnet with M = 2 B/FU integer magnetization with reduced band gap at minority spins. Major difference of magnetic response to pressure characterizes CrNF as a soft ferromagnet versus hard magnetic CrO2. The chemical bonding properties point to prevailing covalent Cr-N versus ionic Cr-F bonding. Different synthesis route are examined.

Introduction

The search for synthetic novel materials with well targeted physical properties is a strategic route for proposing solutions in the field of Materials Science. Simple rules can be followed to devise new stoichiometries with specific properties. As an example boron nitride is an artificial binary which 'mimics' carbon: From the electronic count standpoint, the configurations of the neutral atoms: B ([He] 2s 2 2p 1 ), N ([He] 2s 2 2p 3 ) and C ([He] 2s 2 2p 2 ), are such that "B + N = 2C", i.e. summing up to 8 valence electrons on each side of the equality sign. Cubic boron nitride, c-BN is a good substitute for diamond with a hardness magnitude just below 10, which stands for the value assigned to diamond on the Moh scale of hardness cf. ref. [1] and therein cited works). Consequently beside structural characteristics, the way the same number of electrons is (re)distributed among chemically different atoms engaged in different interatomic bonds can be of relevance to predict new stoichiometries and properties. Considering an anionic substructure for a given tetravalent cation (A IV ), the isoelectronic relationship for valence shell states is established: O + O  N + F, i.e. 2 × (2s 2 , 2p 4 )  (2s 2 , 2p 3 ) + (2s 2 , 2p 5 ). Nitride-fluorides of formulation A IV NF can be considered as pseudo-oxides, isoelectronic of A IV O2. Several nitride-fluorides exist, such as alkaline-earth Ba2NF [2] and Mg2NF [3] and transition metal based TiNF [4] as well as ZrNF [5]. Their investigations carried out within the well established quantum mechanical density functional theory DFT [6,7] led to provide a better understanding of such mixedanions compounds with potential properties like anisotropic electric conduction, optical anisotropy as well as its growth in thin films or at nanoscopic scales [8]. In TiNF [9] and ThNF [10] the arrangements of the N and F substructures were analyzed based on DFT computations focusing on the energetic parameters and the chemical bonding as well as the group-subgroup relationship.

CrO2, one of the rare ferromagnetic oxides at room temperature (TC 400 K) and 2 B magnetization per formula unit (FU) implicitly assigned to tetravalent Cr with 3d 2 configuration, is known for its use in magnetic recording applications [11]. Since its early discovery [12], several works were devoted to its physical properties [13]. DFT band structure calculations go beyond the phenomenological description by identifying a half metallic ferromagnet HMF behavior [14,15], i.e. it is metallic for majority spins () and insulating for minority spins () (cf. Fig. 5). The calculated magnetic band structure was confirmed later on by magneto-optic spectroscopy investigations [16].

The present work presents a proposition of the crystal, electronic, and magnetic structures as well as the bonding properties of chromium nitride fluoride CrNF, devised in the framework of the DFT from the electronic rules presented above and the group-subgroup relationships.

Close relationships for the electronic and magnetic properties of CrNF are established with isoelectronic and structurally related CrO2, well known in high density magnetic recording media during last century.

Computational details

Two computational methods within the DFT were used in a complementary manner. The Vienna ab initio simulation package (VASP) code [17,18] allows geometry optimization and total energy calculations. For this we use the projector augmented wave (PAW) method [18,19], built within the generalized gradient approximation (GGA) scheme following Perdew, Burke and Ernzerhof (PBE) [20]. The conjugate-gradient algorithm [START_REF] Press | Numerical Recipes[END_REF] is used in this computational scheme to relax the atoms. The tetrahedron method with Blöchl corrections [19] as well as Methfessel-Paxton [START_REF] Methfessel | [END_REF] Brillouin-zone (BZ) integrals were approximated using the special k-point sampling of Monkhorst and Pack [23]. A large number of 150 k-points was used in the irreducible wedges of the respective BZ to perform total energy calculations. The optimization of the structural parameters was performed until the forces on the atoms were less than 0.02 eVÅ -1 and all stress components less than 0.003 eVÅ -3 . Calculations assuming firstly a non spin polarized (NSP) configuration then spin polarized (SP) calculations were carried out. Subsequent all-electron calculations, equally based on DFT with a GGA functional were carried out with the optimized parameters for a full description of the electronic and magnetic structures and the chemical bonding properties. They were performed using the augmented spherical wave (ASW) method [24,25] in a scalar relativistic implementation [START_REF] Eyert | The Augmented Spherical Wave Method[END_REF]. In the ASW method, the wave function is expanded in atom-centered augmented spherical waves, which are Hankel functions and numerical solutions of Schrödinger's equation, respectively, outside and inside the so-called augmentation spheres. In the minimal ASW basis set, we chose the outermost shells to represent the valence states and the matrix elements were constructed using partial waves up to lmax + 1 = 3 for Cr and lmax + 1 = 2 for N and F. Low energy lying F(2s) were considered as core states. In order to optimize the basis set, additional augmented spherical waves are placed at carefully selected interstitial sites (IS). Self-consistency was achieved when charge transfers and energy changes between two successive cycles were such as: Q < 10 -8 and E < 10 -6 eV, respectively. The BZ integrations were performed using the linear tetrahedron method within the irreducible wedge [19]. The calculations are carried out assuming firstly spin degenerate, non spin-polarized (NSP), then spin polarized (SP) magnetic calculations for both ferromagnetic and antiferromagnetic long range orders are done for the identification of the ground state magnetic structure of the most stable model structure. The relative magnitude of the chemical bonding is obtained based on the overlap population analysis: Sij, i and j being two chemical species. The crystal orbital overlap population (COOP) criterion is used [27]. In the plots positive, negative and zero COOP magnitudes indicate bonding, anti-bonding and non-bonding interactions respectively.

3 Geometry optimization and EOS.

3-1 Geometry optimization

The rutile ground state crystal structure of CrO2 is tetragonal with P42/mnm space group (SG) N o 136. For CrO2 a = 4.23 Å and c = 2.92 Å with d(Cr-O) ~1.87 Å and two formula units (FU) per cell leading to the Cr2O4 stoichiometry. The atomic positions are Cr 2a at (0,0,0); O 4f at (x,x,0; -x,-x,0; -x+½,x+ ½,½; -x+½, -x+½,½), x ~0.305 [6][7][8][9][10]. For CrNF, a starting setup model was done by ordering N and F at the O (4f) positions leading to the Cr2N2F2 stoichiometry. By carrying out unconstrained geometry optimization runs, the symmetry becomes orthorhombic. Such symmetry is found in anion ordered rutile-type derivatives obeying group to subgroup relations examined by Baur [28] with two orthorhombic structures: CoReO4 (Cmmm SG N°65; Z = 4 FU) and MgUO4 (Imam SG N° 74; Z = 8 FU). Both exhibit cationic and anionic ordering (cf. Table 1) and constitute appropriate host structures for modeling the ternary compound CrNF. The two candidate structures are confronted through determining the ground state structure and deriving equilibrium values based on the energy-volume equation of state (EOS). Full analyses of the electronic and the magnetic structures as well as the properties of chemical bonding are subsequently done for the ground state structure properties.

For the purpose of geometry optimization and search for the ground state structure, non spin polarized (NSP) were carried out starting from CoReO4 (based centered orthorhombic) and MgUO4 (body centered orthorhombic) types [28]. The respective orthorhombic symmetries were preserved at self consistent convergence. The fully relaxed structure parameters and the corresponding energies and volumes are given in Table 1 In view of the ferromagnetic ground state of CrO2 with 2 B/FU magnetization, further geometry optimized calculations were carried out for isoelectronic CrNF considering two spin populations, i.e. spin polarized (SP) calculations. In both structural types a 2 B/FU was found at self consistent energy convergence with a larger volume and stabilized energy: E (CrNF Cmmm)= -0.24 eV/FU and E (CrNF Imam)= -0.35 eV/FU. This suggests a tetravalent character of Cr similarly to CrO2. This point is discussed further while deriving the energy -volume equation of state (EOS) and the electronic and magnetic structures in next sections.

3.2-Energy-volume equations of states

For confronting the relative stabilities of two varieties in NSP and SP configurations, one needs establishing the energy-volume equation of state (EOS). In fact the calculated total energy pertains to the cohesive energy within the crystal because the solution of the Kohn-Sham DFT equations yields the energy with respect to infinitely separated electrons and nuclei. In as far as the zero of energy depends on the choice of the potentials, somehow it becomes arbitrary meaning it is shifted but not scaled. However the energy derivatives as well as the EOS remain unaltered. For this reason one needs to establish the EOS and extract the fit parameters for an assessment of the equilibrium values. This is done from (E,V) set of calculations around minima found from geometry optimization. The resulting E = f(V) curves are shown in Fig. 2. They have a quadratic variation which can be fitted with 3 rd order energy-volume Birch-Murnaghan EOS [29]:

E(V) = Eo(Vo) + [9/8]VoBo[([(Vo)/V]) [ 2/3] 1] 2 + [9/16]Bo(B  4)Vo[([(Vo)/V]) [ 2/3] 1] 3 ,
where Eo, Vo, Bo and B ' are the equilibrium energy, the volume, the bulk modulus and its pressure derivative, respectively. The fit results in the inserts of Fig. 2 show the trend of smaller volumes and higher energies of the NSP forms versus SP pointing to stable magnetically polarized compounds. The energies are such that ESP-NSP (CrNF Cmmm) = -0.29 eV/FU and ESP-NSP (CrNF Imam) = -0.35 eV/FU. Interesting features appear for the zero pressure bulk modules which are systematically higher for the NSP configuration as with respect to SP as one would expect from the larger volume for the latter due to the magnetic (negative) pressure. At large volume the separation between the SP and NSP E(V) curves is maximum while they meet at low volume. This is expected in as far as the effect of decreasing volume (increasing pressure) is to reduce the localization of the d states needed for the onset of the magnetic moment. Consequently the magnetization should vanish at high pressure, i.e. at volumes below the equilibrium.

3.3-Magnetic response to exerted pressure

From above an estimate of the pressure magnitude required to destroy the magnetization is obtained from the Birch relationship providing the pressure [29]:

P = (Bo/B') [(Vo/V') B' 1].
Taking the fit values Bo, B', Vo of the SP curve for the MgUO4-derived structure and the volume V' as V(NSP), P is calculated to be 11 GPa. This value is much smaller than the one calculated for CrO2 of 120 GPa [30]. CrNF is then assigned a "soft" magnetic behavior while CrO2 is a "hard" ferromagnet. This is further illustrated with the plots of the volume change of the magnetization shown for MgUO4-type CrNF and for CrO2 in Fig. 3

All-electrons calculations.

The electronic and magnetic structures and the chemical bonding properties were calculated for the ground state structure of CrNF with the optimized structure parameters.

Calculations with the CoReO4-type CrNF led to similar results.

Firstly NSP calculations were done to assign a role of each chemical species in the valence band (VB) and conduction band (CB). With the ASW method one starts from neutral atoms; at self consistency (energy and charge) using 432 k-points generated from 1728 (12  12  12) k-points in the orthorhombic BZ, there is charge redistribution between the atomic species. This does not translate the ionic picture as F -despite the large ionic character of fluorine but the charge transfer is actually observed from Cr to N, F and IS, in agreement with the expectations of chemistry, i.e. from the metallic species to the anions (N, F) and IS which receive charge residues, i.e. less than 0.1 electron.

Non magnetic calculations and bonding properties.

The NSP site projected DOS (PDOS) accounting for site multiplicity are shown in Fig. 4. The energy reference along x is at the Fermi level EF crossing a large PDOS intensity arising from Cr d states which are mainly centered in the empty conduction band, above EF.

For the sake of clarity we do not show small contribution PDOS arising from IS. The lower part of the valence band (VB) is dominated by N(2s) at ~ -14 eV then F(2p) in the range {-10, -7 eV} and N(2p) in the range {-5, -1 eV}. The PDOS below the two p blocks point to a larger Cr-N mixing than Cr-F. This is further made explicit in the analysis of the chemical bond in next section. Besides the mixing of states the NSP DOS allow assessing the instability of the compound from the large Cr PDOS at EF in such a non magnetic configuration. The mean field Stoner theory of band ferromagnetism [START_REF] Mohn | Magnetism in the solid state -An introduction, springer series[END_REF] can be applied to address the tendency for spin polarization. The total energy of the spin system results from the exchange and kinetic energies. Referring the total energy to the non-magnetic state, this is expressed as: E = constant{1-In(EF)}. In this expression, I (eV) is the Stoner integral, which is calculated and tabulated for the metals by Janak [START_REF] Janak | [END_REF] and n(EF) (1/eV) is the PDOS value for a given state -mainly d-at the Fermi level in the non-magnetic state. If the unitless Stoner product In(EF) is larger than 1, E is lowered and the system stabilizes in a magnetically ordered configuration. Then the product In(EF) provides a criterion for the stability of the spin system. From ref. [START_REF] Janak | [END_REF] the value of I{Cr(3d)} is 0.38 eV. With n(EF) value of 9 eV -1 , the calculated I.n(EF) is 3.4 and the Stoner criterion 1-In(EF) is largely negative leading to energy lowering upon the onset of magnetization, i.e. intra-band spinpolarization should occur when spin polarization is allowed. Subsequent spin polarized (SP) calculations with equal initial  and  spin populations, were carried out. At selfconsistency a finite magnetization can be identified within an implicit long range ferromagnetic order. Nevertheless if magnetic exchange energy is not sufficient, zero local moments can result, so that the calculations are not biased initially. This also depends on of the BZ mesh precision, i.e. the calculations are usually carried out in steps of increasing precision until no more changes are observed in the variational energy and the magnetic moments.

The lower panel of Fig. 4 shows the bonding characteristics based on the analysis of the overlap populations Sij with the COOP criterion. The bonding within the valence band VB is mainly between chromium on one hand and N and F on the other hand. It is differentiated between the two chromium sites with Cr1-F at -10 eV, Cr2-N at -4 eV. Note the low magnitude Cr2-F COOP versus Cr1-N. The overall prevalence of Cr-N over Cr-F follows the shorter Cr-N versus Cr-F distances and signals the rather covalent nature of the nitride fluoride. The itinerant Cr-d states mix with the p-states of N and F following their PDOS positions in the DOS panel. The bonding between the two anionic sublattices is of antibonding nature because of the involvement of their states in the bonding with Cr. The large antibonding Cr1-F EF COOP's are due to the near saturation of F-2p i.e. close to 2p 6 so that the electrons brought by the Cr-F bond go into anti-bonding COOP's. This is opposed to the Cr-N COOP's at EF which have much less intensity. The large anti-bonding COOP's at EF is another signature of the instability of the nitride fluoride in such a spin degenerate configuration. The negative moment on nitrogen is of induced nature and arises from the Cr-N strong hybridizing shown in the preceding section. This result indicates anew, i.e. with respect to section 3 results, the formal tetravalent Cr IV oxidation state (3d 2 ).

In fact the integer total magnetization results from the contribution of all the chemical constituents, not only Cr. It is illustrated in the plot of the SP PDOS shown in Fig. 5. The magnetic exchange causes the majority spins to shift down in energy and the minority spins to shift up in energy, proportionally to the magnitude of the developed magnetic moment. This is actually observed for Cr d states and less for F whereas the opposite is observed for s and p states of nitrogen which possesses a negative moment. The striking feature is observed at the Fermi level with large PDOS for  spin mainly due to Cr1 and Cr2 and a gap opening of ~0.5 eV for  spin. Comparing the SP DOS with those formerly obtained for CrO2 in lower panel of Fig. 5 (ref. [15]) -low energy lying O(2s) are not shown-, similarities can be traced out as to the metallic  spin and the gap opening at  spin. However the gap magnitude of ~2 eV is much larger than the band gap in CrNF. This arises from the larger covalent overall Cr-N bonding versus Cr-O. It also contributes to explain the different magnetic responses of the two compounds to pressure as discussed above.

Lastly for a search of the magnetic ground state, antiferromagnetic AF orders were enforced. This can be justified in so far that CrN is AF in its orthorhombic structure (cf. [33] and therein cites references). The Cr12Cr22N4F4 stoichiometry in the MgUO4 -type ground state structure (Table 1) allows splitting the unit cell into two magnetic subcells with one Cr1Cr2N2F2 subcell as (UP SPINS) and the other Cr1Cr2N2F2 subcell as (DOWN SPINS), thus enforcing an AF of Gtype, i.e. with antiparallel alignment of the atoms within a plane. At self consistent energy convergence there is full compensation of spins and an energy difference E(SP-F -SP-AF) = -0.13 eV/FU, thus favouring the SP-Ferro state.

A-type AF order between planes was simulated by a double supercell along c. The calculations also favoured the SP-F state with E(SP-F -SP-AF) = -0.17 eV/FU. Then the long range order of CrNF ground state is ferromagnetic.

Conclusions and synthesis routes prospects

Isoelectronic rules and (space)group-subgroup relationships coupled with quantum mechanical computations of energy derived quantities within DFT have been shown in earlier works to allow prediction of new materials as 2D and 3D carbon related boronnitrogen binaries and carbon-boron-nitrogen ternaries. These rules applied here to the anionic sublattice in the prominent material CrO2 and replacing oxygen by nitrogen and fluorine providing CrNF, lead to results that let expect an interesting material with new properties especially magnetic ones having close relationship with the oxide. It is expected that such rules can be generalized to propose to the experimentalists other new materials in the future. In view of the favorable cohesive energies of CrNF in both rutile derivative model structures and the identified close relationship with technologically important CrO2, synthetic routes are presently being evaluated mainly along different protocols: a-the ammonolysis of chromium fluoride precursor such as the protocol used on ZrF4 for obtaining ZrNF [34]; b-the method using fluorinated fluxes on CrN following a protocol proposed [5] for preparing ZrNF.

The latter protocol seems relevant from the electronic structure characteristics of CrNF, whereby within CrN, covalent Cr-N bonding is expected to be modified but still keep prevalent upon introducing fluorine either by F2 fluxes or by use of F2 pressure [START_REF] Largeteau | ICMCB-CNRS[END_REF].

  . The energies are close but point out to MgUO4 -type as the ground state with E = -0.175 eV/FU. The distance trends of d(Cr-F) > d(Cr-N) shown in both sets of results follow from the ionic versus covalent bonding and interestingly d(Cr-O) is situated in-between: d(Cr-F)~1.93-2.00 Å > d(Cr-O)~1.87 Å> d(Cr-N)~1.85-1.77 Å. The sequence of distances is found to follow inversely the sequence of electronegativity of the anions: N (3.04) < O (3.44) < F (3.98). The resulting structures are shown in Fig. 1 highlighting the distorted CrN3F3 octahedra as it can be seen from the different Cr-N and Cr-F distances shown in Table 1, differentiated for Cr1 and Cr2 surroundings. It needs to be stressed that in both model structures the CrN3F3 motifs arrange along c direction into files of edge sharing octahedra, as in rutile. The cohesive energies for CrNF candidate structures are obtained and compared with that of CrO2 by subtracting the computed atomic energy of Cr (E = -9.50 eV) and the corresponding dimers from the total energy: E (O2): -11.21 eV, E (N2): -20.78 eV, E (F2): -3.71 eV. Then Ecoh.(CrNF Cmmm) = -1.045 eV/FU, Ecoh.(CrNF Imam) = -1.107 eV/FU and E coh.(CrO2: P42/mnm) = -1.303 eV/FU. Clearly the oxide is more cohesive than the nitride fluoride which is more cohesive in its ground state Imam MgUO4-derived structure. These results point out to the possibility of actually synthesizing the compound and presently investigated protocols are developed upon in the last section of the paper. Lastly testing of the disordering effects of N and F by partials occupations of their crystal lattice sites was done in MgUO4-type ground state structure, i.e. F(8i) and N(8h) led to a large raise of the energy by -1.08 eV/FU.

  . The magnetization increases with volume and reaches saturation with 8 B/4 FU -½ cell-for CrNF and 4 B/2 FU -1 cell-for CrO2. What needs highlighting in the plots are the different behaviors of the two compounds around the equilibrium volume V0 indicated by a vertical line: in CrNF there is a nearly immediate decrease of the magnetization below V0 whereas in CrO2 the saturation magnetization is kept far below V0. This nicely reflects the pressure magnitudes above.

4. 2

 2 Ferromagnetic state and search for the magnetic ground state. Upon charge and energy self-consistent convergence the total variational energy decreases largely with respect to the NSP configuration by E(SP -NSP)= -1.41 eV/cell. This gain of energy is due to the magnetic exchange leading to the onset of a magnetic moment as it is shown below. The total magnetization is 8.0 B/4 FU, i.e. 2.0 B/FU with different local spin moments distributions: M(Cr1) = 2.183 B, M(Cr2) = 2.317 B; M(F) = 0.025 B and M(N) = -0.348 B with the remaining contributions from IS (interstitial part).

FiguresFig. 1

 1 Figures

  Fig.1 (color): CrNF in CoReO 4 (top) and MgUO 4 (bottom) ordered rutile derivatives highlighting the edge sharing mixed anion (N,F) octahedra at Cr1 (grey) and Cr2 (yellow) sites.

Fig. 2 Fig. 3

 23 Fig.2 Energy-volume curves for CrNF in CoReO 4 -type for 2 FU (top) and MgUO 4 -type for 4 FU (bottom) in spin degenerate (NSP) and spin-polarized (SP) configurations and fit parameters from Birch-Murnaghan equations of states.

Fig. 4 (

 4 Fig.4 (color) CrNF with ground state MgUO 4 structure: Spin degenerate site projected DOS (top) and chemical bonding (COOP).

Fig. 5 ( 2 FIG. 5

 525 Fig. 5 (color) Site and spin projected density of states of CrNF in ferromagnetic ground state (top) and CrO 2 (bottom). Gap is indicated by large arrow.

Other aspects related to the temperature and other preparation conditions (bulk, films…) are also under consideration.

Acknowledgements:

Discussions and exchange on the topic of nitride-fluorides and on the synthesis routes of CrNF with Prof. Gérard Demazeau, Dr Alain Largeteau, Dr Nicolas Penin, Prof. Jean Etourneau (France), Prof. M.A. Subramanian (USA). The author is grateful for computational facilities provided by the MCIA-University of Bordeaux. Support from Conseil Régional d'Aquitaine is gratefully acknowledged.


Kevin Jordao

Julien Pettré

Marie-Paule Cani

M.-P Cani

Density-controlled crowds

Keywords: Descriptors (according to ACM CCS): I, 3, 7 [Computer Graphics]: Three-Dimensional Graphics and Realism -Animation

We present a novel algorithm for populating and animating an environment of a given density. Our algorithm considers two kind of densities: static density (when people doesn't move) and dynamic density (when people move). The resulting flows of characters are computed according to a density map provided by a user. We discuss the advantages and drawbacks of the proposed solution, as well as future work.

Introduction

Designers need tools to populate virtual environment in an easy and intuitive manner. In a crowd simulation, density results from agents goals and traffic conditions: density is not explicitly controlled and it can be difficult to obtain a desired result. Density is a key factor for the general aspect of a crowd. However, there is still a lack of techniques allowing to populate virtual environments by controlling the density.

In this work, we suggest using crowd patches [START_REF] Yersin | Crowd patches: populating large-scale virtual environments for real-time applications[END_REF] to enable populating a virtual environment with controlled density. Crowd patches are small elements of precomputed moving crowd animations. Animated characters in patches can be either exogenous (navigate from patch to patch) or endogenous (always remain inside one given patch). Patch animations are periodic in time so that they can be played endlessly. Boundary conditions of patches are controlled, so that patches can be connected together to compose seamless large animations. Boundary conditions are controlled by patterns which stand at the interface of two connected patches, and which register all the spatio-temporal waypoints for exogenous characters.

Based on the crowd patches technique, our objective is to populate a given environment with controlled density. As an input, the user provides a map of desired density for moving and static people. As an output, we compute the crowd partches to get the final animation.

Our solution is an iterative process. The solution starts with a set of empty patches that cover the user's map. Iter-atively, we fill patterns with the required number of spatiotemporal waypoints. Then, we complete crowd patches content with endogenous characters. From all these constraints (spatio temporal way-points, and presence of endogenous characters), a final step is to compute collision free trajectories for exogenous characters.

Computation of density-controlled crowd patches

Our algorithm input is a density map. This map is a regular grid adapted to the environment. The user defines a density value for each cell, as illustrated in Figure 2. Thz output of our algorithm is a set of crowd patches, one for each cell. The animation of each patch matches the user-defined density. The resulting crowd is composed by endogenous and exogenous characters animated by the patches .

The local density of the user's map will define the flow through patches. This flow is directly controlled by the number of input/output points at the interface of patches, called patterns. The goal of our algorithm is to find how to dispatch these inputs/outputs on each pattern, then to compute the trajectories between them.

The first step of our algorithm is to create the patterns of each crowd patch, i.e., to define a set of input/output points for characters in patches, according to the density map. Knowing that patches should have periodic animations as defined in [START_REF] Yersin | Crowd patches: populating large-scale virtual environments for real-time applications[END_REF], all the patterns of patches should have the same number of input and output points.

The process for computing the patterns is the following : • We take the patch with lowest density for computing its patterns; • We add a input/ouput on a side of the patch. The side is determined by the densities of the neighbors patches. Highest the density of a neighbor patch, the more likely this shared side will be choosen. We repeat this step until no more input/output points can be inserted; • Then some patches can still present an unsatisfying density. So, we add randomly static people inside to reach the desired density.

Once we have computed patterns and locations of static people for each patch, we need to compute the inner trajectories. For that we:

• Connect each input point to an output point. We formulate and solve this step as a stable marriage problem. . • Deform trajectories to solve residual collisions.

Preliminary results

Figure 2: Result get by a random density map.

We implemented our solution and obtained some early results. We tested our solution on the example presented in 2. Here the map is generated randomly, to test our algorithm. We can observe more people in dense areas than in others areas.

Future work and conclusion

Our technique allows for an easy design of crowd with various densities. However in its current state, our technique shows some limitations.

A mobility parameter Because a crowd can have different proportions of static/dynamic density, we want to add to our algorithm a mobility parameter. This parameter will ensure the percentage of static and moving people in the scene. For instance it will enable to easly create a concert crowd by setting this parameter to zero, or a train station by setting this parameter to 0.5 or a public place with the parameter set to 1.

Accurate location of static people For now, our technique is fully random for locating static people. We can imagine ing a static map to guide our algorithm for locating people in the scene. It will avoid obtaining people waitting in undesired place, and add people waiting in some points of interest, like information panels or windows shop .

Intuitive interface We use a script for the density map information. This mecanism is not as intuitive as a graphic interface. A way to improve this system is to create an interface where user could paint the density of the crowd and see directly the result, to adjust in real time densities parameters.

Evolution of the crowd density in time

The density in a crowd evolve with time (daily cycles). We want to improve our algorithm to be able, given a set of timed density maps, to compute the evolution of densities during time. It will enable to build environments with cyclic days.

Our method is a first step to a density crowd control. It allows to build a crowd with a given denstity, but we still need to work on improving our algorithm and our interface.

Figure 1 :

 1 Figure 1: The pipeline of our algorithm. A density map as input (left). Patterns are deduced from the map (middle). Crowd patches are computed from the patterns to create the desired crowd (right).


Xiaoming Cui

Heping Sun

Séverine Rosat

Jianqiao Xu

Jiangcun Zhou

Bernard Ducarme

Investigation of the time variability of diurnal tides and resonant FCN period

Keywords: Temporal Variation, Superconducting Gravimeter, Diurnal Tides, VLBI, FCN Period

The time variability of diurnal tides was investigated by analyzing gravity observations from global superconducting gravimeter (SG) stations with running time intervals. Through least-square and Bayesian approaches, FCN resonance parameters were estimated for each data section after obtaining the tidal parameters of mainly diurnal tidal waves. The correlation of the time variation in diurnal tidal waves and FCN period was discussed. For comparison, a similar method was used to analyze VLBI observations to study the time variability of nutation terms and FCN period. The variation trend of the FCN period totally depends on the Ψ1 wave in tidal gravity and on the retrograde annual term in nutation. We observed a similar variation trend in the FCN periods obtained from different SG stations worldwide and VLBI observations. The relation between diurnal tides and LOD variations is discussed and the possible mechanisms of the decadal variation in FCN periods were discussed.

Introduction

Due to the interaction between the elliptical liquid core and the solid mantle of the Earth, a retrograde rotational mode called the Free Core Nutation (FCN) occurs in the celestial reference frame, appearing as a nearly diurnal free wobble (NDFW) in the terrestrial reference frame. The period of FCN strongly depends on the flattening of the core-mantle boundary (CMB) and can also be influenced by other coupling mechanisms, such as visco-electromagnetic and topographic coupling. As the most active thermal boundary layer inside the earth, the CMB is typically accompanied by changes in physical characteristics and structure which will possibly affect those coupling mechanisms and then may cause some variations in the FCN period.

FCN parameters (period and quality factor) can be estimated by the resonance enhancement in observations of earth tidal waves or nutation terms with frequencies close to its eigenfrequency. The most commonly used data are high-precision time-varying gravity observed by superconducting gravimeters (SGs) and nutation observations from VLBI network analysis. The FCN resonance parameters can also be inverted from strain data (e.g. [START_REF] Amoruso | Free Core Resonance parameters from strain data: sensitivity analysis and results from the Gran Sasso (Italy) extensometers[END_REF]. The usual approach in determining FCN parameters is linearized least-square method (LSQ) [START_REF] Defraigne | Stacking gravity tides measurements and nutation observations in order to determine the complex eigenfrequency of nearly diurnal free wobble[END_REF][START_REF] Sun | Earth's free core nutation determined using C032 superconducting gravimeter at station Wuhan/China[END_REF]. [START_REF] Florsch | Bayesian estimation of the free core nutation parameters from the analysis of precise tidal gravity data[END_REF] proposed a Bayesian approach that more effectively solves non-linear inverse problems. The temporal variation in FCN period has been investigated through VLBI and SG observations in numerous previous works. Some of these studies [START_REF] Roosbeek | The free core nutation period stays between 431 and 434 sidereal days[END_REF][START_REF] Hinderer | Are the free core nutation parameters variable in time?[END_REF][START_REF] Lambert | The Earth's core parameters as seen by the VLBI[END_REF][START_REF] Vondrak | Stability of period and quality factor of free core nutation[END_REF] found no evidence of variation in the FCN period, whereas [START_REF] Xu | Temporal variations in free core nutation period[END_REF] showed a decadal temporal variation in the FCN period by analyzing long SG records at Brussels station. This problem remains a matter of dispute. Conclusions are difficult to draw when only one kind of observation is made because observed variations are small (sometimes within error bars).

The accumulation of SG data offers a good opportunity to perform the study by comparing the results from two different techniques. The Global Geodynamics Project (GGP) sponsored and organized by the Solid Earth's Deep Interior (SEDI) in IUGG in 1997, has accumulated more than 10 years of observations at many global SG stations, and especially, over 20 years of SG observations in Strasbourg station (France). In studying the temporal variation of the FCN period, observations from different SG sites can be compared to check the reliability of the variation as well as with VLBI observations. In this paper, based on observations from seven SG stations, we use running time intervals to conduct tidal analyses and estimate FCN parameters section by section. We then compare the variations of the FCN period with those of the amplitude factors of some diurnal tidal waves.

Similarly, we analyze VLBI observations and study the transfer function values of nutation terms used to determine FCN parameters and their relation with the variations of the FCN period. We also focus on the differences by using different datasets (SG and VLBI) and approaches (LSQ and Bayesian).

In the following sections, we introduce the data, the FCN resonance model, and the determination approaches of FCN parameters. The data analyses and the temporal variations in FCN period are then discussed.

Data, FCN model, and parameter determination methods

Data introduction

The SG data used in this paper were obtained at seven stations. Five of these stations are located in central Europe, where ocean tide effects at diurnal frequencies are small and for which we usually obtain high-quality data (e.g. [START_REF] Rosat | Noise Levels of Superconducting Gravimeters: Updated Comparison and Time Stability[END_REF] with more stable observation environments. The exception is Brussels site which was quite noisy but for which the data record is interestingly long. Among the five stations, Strasbourg station has the longest-term observations, from 1987 to 2011, which is the main time series we focus on. Other stations are studied for compensation. Considering the possible regional effect, two SG stations in Canada and Australia are also included for comparison. Basic information on each station and data length are listed in Table 1.

The VLBI observation used in this study is the IVS combined solution "ivs12q3X" (International VLBI Service for Geodesy & Astrometry) in the form of celestial pole offsets (dX, dY) referred to IAU2006 precession-nutation model [START_REF] Wallace | Precession-nutation procedures consistent with IAU 2006 resolutions[END_REF] excluding the free core nutation.

These data are given in non-equidistant intervals of 1 day to 7 days. We choose the same data length from 1987-07 to 2011-07 as Strasbourg SG observations in our study.

FCN resonance model

In the terrestrial reference frame, the NDFW will lead to an obvious resonance enhancement in observations of diurnal tidal waves with near frequencies. The FCN parameters are estimated by fitting the observed complex gravimetric factors (complex ratio of the observed tidal amplitude to the tidal amplitude for a solid Earth model given a tide-generating potential) to a damped harmonic oscillator modeling the resonance. The observed values of gravimetric factors for each tidal wave could be obtained by analyzing SG data (see paragraph 3). For a diurnal tidal wave of frequency σ, complex gravimetric factors can be described as follows [START_REF] Hinderer | Global earth dynamics and induced gravity changes[END_REF]:

0 j j nd a δ δ σ σ = + - , (1)
where δ0 is the gravimetric factor independent of the frequency and not influenced by the resonance, a ~ is the complex resonance strength related to the geometric shape of the Earth and the rheological properties of the Earth's mantle, and nd σ ~ is the complex eigenfrequency of the NDFW. + Ω , where Ω is the sidereal frequency of the Earth's rotation.

In the celestial reference frame, the resonance of the liquid core affects the Earth's nutation as well. The resonance for a nutation term with frequency σ, usually represented by a transfer function, can be expressed as follows (eqs.(42) in [START_REF] Mathews | Modeling of nutation and precession: New nutation series for nonrigid Earth and insights into the Earth's interior[END_REF]:

4 0 0 1 ( ) 1 (1 )( ) 1 R R Q e T Q N e s α α α σ σ σ σ = -   = + + + ∑   - -   , (2)
where e R is the dynamic ellipticity of the rigid earth; s 1 , s 2 , s 3 , and s 4 respectively represent the frequencies of the Chandler Wobble, FCN, Free Inner Core Nutation and Free wobble of the inner core; 0 N and 0 Q are constants defined in Table 6 of Mathews et al. (2002); and Q α are resonance strengths of the respective rotational modes.

The main diurnal tidal waves usually employed to estimate FCN parameters include P1, K1, Ψ1, and Ф1, in which Ψ1 is the closest in frequency to the NDFW and thus has the most significant amplitude enhancement. For nutation, thousands of nutation terms are involved in the MHB2000 nutation model [START_REF] Mathews | Modeling of nutation and precession: New nutation series for nonrigid Earth and insights into the Earth's interior[END_REF]. Only those terms with frequencies close to the FCN frequency and with large amplitudes are important to determine the FCN parameters, such as the 365.26 day (l'), 182.62 day (2F-2D+2 ), 121.75 day (l'+2F-2D+2 ), 27.55 day (l), and 13.66 day (2F+2 ) nutation terms, where l, l', F, D, and correspond to the Delaunay variables characterizing the nutation terms.

As we know, each nutation term with an elliptic orbit can be expressed as the sum of two circular terms which have identical frequencies with opposite signs. Due to the same gravity potential source (mainly degree 2 and order 1 component in the spherical harmonic decomposition of the tide-generating potential), a corresponding relationship exists between the tidal waves and the nutation terms. Hence, the Ψ1 wave corresponds to the retrograde annual nutation term (-365.26 days), which is the closest to the FCN, and the two terms should be most sensitive to any variation of FCN frequency. If some temporal variations in the FCN period exist, the diurnal tidal waves or nutation terms near the FCN would be sometimes closer and sometimes farer to the resonance, in pace with the variations. As a consequence, the amplification of these waves or nutation terms would vary according to the resonance model (1) or (2).

Review of the linearized least-squares and Bayesian approaches

The estimation of FCN parameters consists in solving a non-linear inverse problem. The approach mostly used in previous work is a linearized LSQ method optimized by the Levenberg-Marquardt algorithm [START_REF] Marquardt | An algorithm for least-squares estimation of non-linear parameters[END_REF]. Equations ( 1) and ( 2) are used to construct the merit function and to evaluate unknown parameters by an iterative process given initial values. For nutation, we directly use Equation ( 2) as the calculation model. However, for earth tidal waves, we usually use the observed O1 wave as a reference (e.g., [START_REF] Defraigne | Stacking gravity tides measurements and nutation observations in order to determine the complex eigenfrequency of nearly diurnal free wobble[END_REF][START_REF] Ducarme | Determination of the free core nutation period from tidal gravity observations of theGGPsuperconducting gravimeter network[END_REF] to reduce the effects of any systematic discrepancies (e.g., calibration errors) or local environmental disturbances in the fitted resonance parameters of the FCN. The O1 wave has a relatively large amplitude and high signal-to-noise ratio and is minimally influenced by FCN resonance, so it can be observed very accurately. Subtracting the contribution of O1 from both sides of Equation (1) yields

nd O nd j O j a a σ σ σ σ δ δ 1 1 - - - = - . ( 3
)
The values of the unknown parameters ( a ~ and nd σ ~) are optimal when the merit function f is minimal. The function f can be written as

[ ] 2 1 1 ) (         - - - - - = ∑ nd O nd j O j a a p f σ σ σ σ δ δ σ σ , (4)
where p(σ) is the weight function of tidal parameters with frequencies σ.

The Bayesian method has been shown to be more suitable for the inversion of non-linear problems (Tarantola and Valette, 1982a). It can make full use of the data information combined with some prior information on the data or on the unknown parameters. Based on these advantages, the Bayesian approach has been proposed to estimate FCN parameters [START_REF] Florsch | Bayesian estimation of the free core nutation parameters from the analysis of precise tidal gravity data[END_REF]Rosat et al., 2009).

As the quality factor Q of the FCN expresses signal damping caused by any kinds of dissipative processes, the value of Q should be positive. For earth tidal waves, a priori information x=log 10 Q is involved to retain the positivity of Q in Equation (1). The real and imaginary parts of the amplitude factors are written as follows:

0 2 2 2 2 ( ) (( 10 ) / 2) Re( ) ( ) (( 10 ) / 2) ( ) (( 10 ) / 2) Im( ) ( ) (( 10 ) / 2) R R I R x j nd nd th j R R x j nd nd I R R R x j nd nd th j R R x j nd nd a a a a σ σ σ δ δ σ σ σ σ σ σ δ σ σ σ - - - -  - - = +  - +   - +  =  - +  . (5)
For detailed information on the construction of the Bayesian calculation scheme, we refer to Tarantola and Valette (1982a,b) and [START_REF] Florsch | Bayesian estimation of the free core nutation parameters from the analysis of precise tidal gravity data[END_REF]. The unknown parameters are described by the probability distribution P (Rosat et al., 2009) ( )

2 2 Re( ) Im( ) 1 , , , exp 2 th R th I j j j j R R I nd R I j j j P x a a k δ δ δ δ σ δ δ         - -     = - +           ∆ ∆             ∑ , ( 6
)
where k is a normalization factor in order that the integral of this equation is unity; For the nutations, the Bayesian calculation model can be similarly constructed by representing the "nutational" factors with the transfer function values as in [START_REF] Rosat | Free core nutation resonance parameters from VLBI and superconducting gravimeter data[END_REF].

Temporal variations of the FCN period in SG data

Before fitting the FCN parameters with SG observations, we need to obtain the tidal parameters (amplitude factor and phase lag) of earth tidal waves. All SG data were preprocessed to eliminate any disturbance, such as steps, spikes and earthquakes, to ensure the accuracy of tidal analysis. The tidal parameters were then estimated by analyzing 1-h sampled SG data with the ETERNA [START_REF] Wenzel | The nanogal software: earth tide data processing package ETERNA 3.30[END_REF] software. The tidal potential used in the analysis is [START_REF] Hartmann | The HW95 tidal potential catalogue[END_REF]. andΦ1). In our study, 6-year and 3-year intervals are used and the retrieved delta-factors are plotted with respect to the center of the interval. A 6-year time interval is used only for Strasbourg station because of its sufficiently long series and for subsequent convenient comparison with the VLBI results later. For the VLBI data, a 6-year interval should be used to ensure the separation of the FCN and the retrograde annual nutation terms which are very close in frequency. In surface gravity observations, the FCN cannot be directly observed and separated from the close-by tidal waves because the expected amplitude is very small in the terrestrial reference frame. Therefore, our study mainly focuses on its resonance effect in diurnal tidal waves, and the FCN period is estimated from the resonance. Nevertheless, we first need to check whether there is an effect of the wave group separation when the FCN frequency (about 1.00507cpd -cycle per day) is included or not in the frequency band of the waves (K1 &Ψ1) adjacent to FCN. When performing tidal analysis with ETERNA software for each 6-year data section of Strasbourg, we changed the frequency limits for the K1 group in two cases: K1 (1.001824 cpd to 1.004107 cpd) without FCN and K1+FCN (1.001824 cpd to 1.005327 cpd), including FCN.

Similar way is used for Ψ1 (1.005328 cpd to 1.007594 cpd) and Ψ1+FCN (1.004108 cpd to 1.007594 cpd). For K1, including or excluding FCN has a negligible influence (10 -8 ) to the amplitude factor. Considering that K1 is farther from FCN relatively, we also test the K1x+ waves.

The influence is in the order of 10 -5 which is one or two orders of magnitude smaller than the standard deviation of K1x+'s amplitude factor. For Ψ1, the difference for the amplitude factor including or excluding FCN is in the order of 10 -3 which is in the same order as the Ψ1 standard deviation (cf. Fig. 1). Given that Ψ1 corresponds to the 6-month modulation of S1, we also checked the variation in the amplitude factor for S1, and the differences for Ψ1 and S1 are anticorrelated.

This indicates that the differences we obtained are real ones and not perturbations on the tidal records because S1 and Ψ1 are respectively associated with the prograde and retrograde annual terms of nutation.

Considering that this influence is at the level of the precision we have on Ψ1, we take the minimum frequency 1.005328 cpd for the Ψ1 group to exclude the FCN in the tidal analysis. Such a difference shall reflect the noise effect. However, we perform a similar test with the SG data at the Canberra and Cantley stations, as shown in Fig 1 . The difference for δ(ψ1) including or excluding the FCN at the two stations is consistent with the results obtained for Strasbourg. It is very interesting that the results for the 3 stations located on different continents have a very similar behavior particularly between 1998 and 2007. A reasonable explanation may be that this result is the combined effect of tiny waves, FCN, or other tidal signals that are not modeled. This effect should be investigated in a further study. After testing the influence of the group separation in ETERNA tidal analyses, we can process the Strasbourg data with 6-year segments shifted by 1 year. The correlation between the FCN parameters has been investigated in some previous work (e.g., [START_REF] Florsch | Bayesian estimation of the free core nutation parameters from the analysis of precise tidal gravity data[END_REF]Rosat et al., 2009). A strong correlation exists between the FCN period and the gravimetric amplitude factor of Ψ1. In Fig. 2, we plot the time variations of the delta-factors for K1 and Ψ1. The two waves are located at opposite sides of the FCN frequency and should be most sensitive to any variation of the FCN period than other waves. The δ i for M2 and O1 are also plotted.

Since the oceanic tides have the same driving mechanism and similar spectral pattern as body tides, it is not possible to separate their effects by harmonic analysis. Therefore the oceanic influence has to be evaluated using an ocean tides model of the different tidal constituents [START_REF] Melchior | Computations of tidal gravity loading and attraction effects[END_REF]. The tidal loading vector, which takes into account the direct attraction of the water masses, the flexion of the ground and the associated change of potential, is evaluated by performing a convolution integral between the ocean tide models and the load Green's function [START_REF] Farrell | Deformation of the Earth by surface load[END_REF]. According to the principle of vector superposition, the loading effects of oceanic tides can be removed from the obtained gravimetric factors. In this paper, the amplitude factors are corrected for the ocean tide effect on the basis of the Fes04 model [START_REF] Lyard | Modelling the global ocean tides: modern insights from FES[END_REF] and following the computation scheme as in [START_REF] Ducarme | Determination of the free core nutation period from tidal gravity observations of theGGPsuperconducting gravimeter network[END_REF].

As mentioned above, the Strasbourg time series was obtained from two different instruments.

Therefore, a calibration problem may affect our study. Indeed, the tidal factors from the T005 have been multiplied by a factor of 1.001529 to correct for an observed scale factor error. In order to avoid any contamination from errors in the scale factor and to be free from any possible changes in the sensitivity of the gravimeter, we also plot in Fig. 2 the δ i factors for K1, Ψ1, and M2 normalized by the O1 amplitude factor δ(O1).

δ(K1) and δ(O1) are very stable, with a variation less than 0.02%. The variation of δ(ψ1) is larger and about 2%. We can note that the overall variation shape for K1 and Ψ1 is comparatively consistent, especially after normalization to δ(O1). δ(M2) is also very stable with a variation of about 0. 05%. This stability shows that the multiplication by 1.001529 for the old T005 series was necessary to reach such accuracy on the gravimetric factors.

After the tidal parameters have been corrected from the ocean loading, the FCN period can be determined by the methods introduced in the previous section. We used 5 tidal waves to perform the inversion: O1, P1, K1, Ψ1, and Ф1. The usual weights applied in LSQ or Bayesian methods are the inverse of the standard deviations of the tidal amplitude factors resulting from the ETERNA tidal analyses. In our LSQ inversion, the results show a very large standard deviation of about 10 d to 30 d for the FCN period. This deviation is disadvantageous for our study on its variation trend.

We then apply other weights used in [START_REF] Liu | Determination of the Earth's Free Core Nutation Parameters by Using Tidal Gravity Data[END_REF] which will significantly reduce the standard deviation of the period. These new weights are inversely proportional to the standard deviations of the amplitude factors and to an additional factor representing the frequency difference between the wave frequencies and the FCN frequency. The two weighting schemes lead to a difference in the FCN period of a few days, but it has no influence on its time variations. Based on the weights of [START_REF] Liu | Determination of the Earth's Free Core Nutation Parameters by Using Tidal Gravity Data[END_REF], the standard deviation for the FCN period decreases to a few days. In the Bayesian inversion, we employ the usual weighting scheme.

The results of the Bayesian and LSQ approaches are depicted in Fig. 3. The comparison of the amplitude factors for all diurnal waves shows a strong correlation between the FCN period variations and the δ(ψ1) time fluctuations. In Fig. 3, we have also plotted the values of δ(ψ1)/ δ(O1).

For the sake of comparison, we have reversed the vertical axis for the relative amplitude factors in Fig. 3. The obtained apparent time-variations of the FCN period differ by a few days between both methods, but the variation trend is similar. We clearly see that the variation trend of the FCN period and of δ(ψ1) are correlated, which is not surprising since the FCN parameters are strongly dependent on δ(ψ1). Using different weights in the LSQ method does not change the variation trend of the FCN period. By superposing a polynomial fitting curve, we observe a decade fluctuation in δ(ψ1)/δ(O1) and in the FCN period. It is worth to note that it is not precise to fitting a period signal with polynomial curve due to the edge effects, polynomial degrees, etc. But our aim here is to make the variation trend clear and to study the similar variation trend between SG and VLBI results. For further verification, we use intervals of 3 years instead of 6 years for the tidal analyses of the Strasbourg SG data. The resulting Ψ1 amplitude factors are plotted in Fig. 4. The standard deviation of δ(K1) for the old SG T005 series (1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996) is in the order of 10 -4 , whereas for the C026 observations (1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011), it is in the order of 10 -5 . For Ψ1, the order of standard deviation for the two instruments is respectively 10 -2 and 10 -3 . The more recent instrument has obviously a higher observation quality as already mentioned in [START_REF] Rosat | Noise Levels of Superconducting Gravimeters: Updated Comparison and Time Stability[END_REF]. A comparison of the results between 3-and 6-year intervals shows some differences in the individual values, but the variation trend is very consistent which indicates that the fluctuations of the amplitude factors observed at Strasbourg station are reliable. Please note that the vertical axis has been reversed.

To verify the reliability of the observed fluctuation, we analyze the data from the other 6 stations presented in Table 1. The resulting time fluctuations of δ(ψ1) obtained with 3-year running intervals are also depicted in Fig. 4. We have also reversed the vertical axis for the sake of Cantley station is quite different. This site is known to be well-influenced by oceanic effects. For instance the oceanic noise is larger in winter (see e.g. [START_REF] Rosat | Noise Levels of Superconducting Gravimeters: Updated Comparison and Time Stability[END_REF]. As our ocean loading correction does not include any time variability, we cannot discard the hypothesis that the apparent time variation of the delta-factors is influenced by the oceans. For Canberra station, there is also a roughly similar variation as that in Strasbourg station. All the tests show some fluctuations in the Ψ1 amplitude factors and also in the FCN period to a certain extent.

In the next part we will consider independent and global estimates of the FCN period using VLBI nutation observations.

Temporal variations of the FCN period in VLBI data

The estimate of the FCN resonance parameters using VLBI observation is known to be more precise (see for instance [START_REF] Rosat | Free core nutation resonance parameters from VLBI and superconducting gravimeter data[END_REF]. Besides, the nutation solution is obtained globally and hence minimizes the local effects that perturb VLBI antenna, such as the ones that affect the SG sites (local hydrology, oceanic noise …).

As with the SG observations, we use 6-year running intervals to estimate the transfer function values for the nutation terms section by section. We use a method similar to [START_REF] Vondrak | Free core nutation: direct observations and resonance effects[END_REF].

The bias, trend, and two long-period terms (with arguments and 2 ) are first removed from the data sections by using the values fitted from the whole data length. We then use the least-square method to fit the selected dominant nutation terms. For the determination of the FCN parameters with VLBI data based on LSQ and Bayesian approaches, we use only the five main nutation terms mentioned in Section 2.2. Other terms have negligible effect on the results as they are much less affected by the resonance. The observed FCN mode with a 460 day retrograde period is also fitted simultaneously. This mode corresponds to the forced oscillation by the atmospheric and oceanic layers. In figure 5 The nutation observation is also perturbed by global mass redistribution through angular momentum exchanges. Indeed it is worth to note that atmospheric and oceanic contributions play an important part in the observed celestial pole offsets and for the retrieval of nutation amplitudes. To remove this part, our above results include a correction for the S1 thermal atmospheric tide and other Sun-synchronous effects given in [START_REF] Mathews | Modeling of nutation and precession: New nutation series for nonrigid Earth and insights into the Earth's interior[END_REF] for the prograde annual nutation. As the computed effects on the retrograde annual nutation varied wildly when using different atmospheric angular momentum data, [START_REF] Mathews | Modeling of nutation and precession: New nutation series for nonrigid Earth and insights into the Earth's interior[END_REF] chose to apply this correction only to the prograde annual term with amplitude about 100 as. For a better correction of the ocean and atmospheric contributions, here we also numerically integrate the Brzezinski broad-band Liouville (Brzezinski 1994;Brzezinski et al. 2002) to obtain the integrated celestial pole offsets by using atmospheric angular momentum functions (AAMF) and oceanic angular momentum functions (OAMF) data. The detailed numerical process is same as in Vondrak andRon (2009, 2010) and will not be repetitive described here. The AAMF and OAMF we used here were computed from: (1) NCEP/NCAR re-analysis 1987-2011 [START_REF] Salstein | Computing atmospheric excitation functions for Earth rotation/polar motion[END_REF]; (2) ERA+OMCT 1989-2010[START_REF] Dobslaw | Seasonal polar motion excitation from numerical models of atmosphere, ocean, and continental hydrosphere[END_REF]. We only use the AAMF data in (1) because the OAMF offered by ECCO model (Stammer et al. 2002;[START_REF] Gross | Atmospheric and oceanic excitation of decadal-scale Earth orientation variations[END_REF] have not that long period. Removing the integrated celestial pole offsets, we derived the FCN resonance period with the LSQ method again and the results are depicted in Fig. 6(c). Compared to Fig. 6(a), the results obtained after correcting for the NCEP/NCAR AAMF are systematically larger, but the variation trend is similar. This is even more obvious in Fig. 11 of [START_REF] Vondrak | Study of atmospheric and oceanic excitations in the motion of Earth's spin axis in space[END_REF] in which they used NCEP+ECCO data. The fluctuation of the FCN period after correction of the ERA+OMCT angular momentum functions varies more violently, as seen also in [START_REF] Vondrak | Study of atmospheric and oceanic excitations in the motion of Earth's spin axis in space[END_REF]. We can also see a 10 year variation trend but different than in Fig. 6 (a).

Discussion

We have investigated the time variation of the diurnal tides and FCN resonance period from SG and VLBI data. A similar decadal fluctuation in FCN period has been found in the longest SG datasets from Strasbourg station and VLBI data. Compared with the other 6 SG stations, some of them show a roughly similar variation for δ(ψ1), which make the observed fluctuation reliable in certain extent. However, the results for δ(ψ1) at some sites are quite different no matter in amplitude or variation behavior. Thus we cannot jump to any final conclusion concerning this fluctuation. Some time variability in the oceanic effects should be checked.

The estimation of FCN period depends on the determination accuracy of related diurnal tides.

Stacking global SG data (Rosat et al., 2009) could give more reliable FCN parameter in agreement with VLBI results. However to study the long time behavior, we are currently limited by the poor number of available long SG records. Besides, the present results have showed that using either a LSQ or a Bayesian method does not influence the obtained values for the FCN period.

To further verify our study, an improvement in the modeling of ocean tide loading is necessary [START_REF] Ducarme | Determination of the free core nutation period from tidal gravity observations of theGGPsuperconducting gravimeter network[END_REF][START_REF] Rosat | Free core nutation resonance parameters from VLBI and superconducting gravimeter data[END_REF], but this mainly depends on the development of more accurate ocean tidal modeling.

Considering that the ψ1 wave and the retrograde annual nutation are the closest to FCN frequency, respectively in the diurnal tidal band and the annual nutation band, the strong inverse correlation between these two terms and the FCN resonant period is reasonable Such fluctuations in the tidal amplitude factors were already mentioned by [START_REF] Lecolazet | Correlation between diurnal gravity tides and the Earth's rotation rate[END_REF], who indicated some correlations between the variations of some tidal amplitude factors (P1, K1, Ψ1, Ф1, and M2) and the length of day (LOD) variations. In our study, such correlation does not appear.

Only K1 and Ψ1 amplitude factors have relatively similar variations. To check for a possible correlation with LOD, we plot the variations of LOD data offered by IERS EOP observations in Fig. 7 (finals2000A.all). The relative amplitude factors of Ψ1, K1, and M2 are also plotted in Fig. 7. We can see that there is no clear correlation between the tidal waves and LOD fluctuations.

As we know, the decadal variations in LOD (e.g., Buffet, 1996;[START_REF] Holme | Electromagnetic core-mantle coupling-I. Explaining decadal changes in the length of day[END_REF][START_REF] Buffet | Inversion of torsional oscillations for the structure and dynamics of Earth's core[END_REF] are mainly explained by the core mantle coupling. According to the spectral analysis of the LOD data, the main power of the spectrum is at the decade scale, in which the main peak is a component of period about 22 years. Besides, there are also some seasonal and interannual signal components. In order to check whether such a consistent fluctuation in LOD exists, other signal components should be eliminated. First, we use a Butterworth low-pass filter to remove the large components with periods of less than 1 year. We then use a LSQ method [START_REF] Guo | Seasonal and inter-annual variations of length of day and polar motion observed by SLR in 1993-2006[END_REF] to fit the 22-year term and other interannual signals. After these terms are deducted from LOD data, the residuals are plotted in Fig. 7 (continuous red line). The 10-year oscillation is clearly visible. The fluctuation we have found in the FCN period is also about 10 years. Compared with the FCN results of previous figures, there is about a 2 year phase difference between the variations of LOD and FCN period. It is important to point out that more detailed work is needed to process the LOD data and separate the decade scale signal, as well as to explain the phase difference. However, the decade variations in LOD and its relation with core mantle couplings may be helpful to explain some fluctuations in the FCN period. [START_REF] Vondrak | Study of atmospheric and oceanic excitations in the motion of Earth's spin axis in space[END_REF] have also concluded indeed that the observed variations could be real and due to some processes at the CMB, as the current AAMF and OAMF modeling cannot explain the observed variations of the FCN period.

As mentioned in Section 1, the period of the FCN strongly depends on the structure and physical properties at CMB (e.g., [START_REF] Koot | Constraints on the coupling at the core-mantle and inner core boundaries inferred from nutation observations[END_REF]. Indeed, the frequency of the FCN depends on coupling mechanisms at the core boundaries [START_REF] Greff-Lefftz | Magnetic field and rotational eigenfrequencies[END_REF][START_REF] Mathews | Modeling of nutation and precession: New nutation series for nonrigid Earth and insights into the Earth's interior[END_REF].

In a celestial reference frame, it can be written (normalization by ): where e f is the dynamic ellipticity of the outer core; A s , A f , and A m are respectively the mean principal moments of inertia of the solid inner core, the fluid outer core, and the mantle. K cmb and K icb are the complex coupling strength parameters representing the influence of the electromagnetic torques at CMB and ICB, respectively; A p , B p , C p , and D p are the coefficients of the torque due to the driven geostrophic pressure acting on the topography at CMB (Greff-Lefftz and Legros, 1999);

        Ω - + Ω + + + + + -         + - - = ² 2 ² 2 1 1 0 2 f P p f p p f s ICB CMB f m f FCN A A D i A C B f A A K K e A A β β σ , (7)
0 2
f is the ratio of the geostrophic pressure at the CMB over the hydrostatic rotational pressure and β is a compliance that characterizes the deformability of the CMB under the centrifugal force.

2 / 1 0 f h q = β , where f h 1 is the Love number expressing the deformation of the CMB under an inertial pressure [START_REF] Dehant | Analytical approach to the computation of the Earth, the outer core and the inner core rotational motions[END_REF]. Besides, a misalignment (tilt) of the inner core with respect to the mantle could also influence the direction of the Earth's rotation on decade timescale [START_REF] Dumberry | Inner core tilt and polar motion[END_REF] and slightly perturb the FCN frequency through gravitational and pressure torque. Such an expression shows that the eigenfrequency of the FCN could vary when the torques vary. Impact of time-varying torques on the FCN period has not been computed yet and would require the knowledge of the varying pressure and electromagnetic torques at the CMB.

Conclusion

According to our study, the variations in time of the FCN resonance period are within 1 day from VLBI data and within several days from SG records.

Considering the error bars, the variations of the FCN period are however not obvious.

We have shown a similar variation trend in the results obtained from different worldwide SG data.

There exists also a similar trend in the FCN period obtained from VLBI observations as previously observed by [START_REF] Vondrak | Study of atmospheric and oceanic excitations in the motion of Earth's spin axis in space[END_REF]. The test of using two kinds of independent observation (surface gravity on the one hand and VLBI nutation on the other hand) to verity the time variability of the FCN period seems coherent. However longer records at more SG sites would be necessary to confirm such fluctuations. For SG data, improvement of ocean tidal loading correction should be considered in a future study. For celestial pole offsets, further studies are required to reach agreement between the different models of atmospheric and oceanic angular momentum functions.

The variation trend of the FCN period totally depends on the retrograde annual term in nutation and on the Ψ1 wave in tidal gravity. It is indeed the same phenomenon in different reference frames,and the decadal fluctuation seems to exist in both datasets. However the processes responsible for such variation are not clearly identified yet. A time-variability in the ocean loading and angular momentum is one candidate as well as some varying torques at the CMB. A possible correlation with the decadal LOD trend due to some core-mantle coupling mechanism may exist, but no clear evidence or demonstration has been proposed yet. Further theoretical computations are needed to demonstrate a possible time fluctuation in the FCN period of geophysical origin.

  R and I represent the real and imaginary parts, respectively. The quality factor Q and the FCN eigenperiod T FCN are expressed as

  deviations of the real and imaginary parts, respectively, of the gravimetric factors.

  First, we consider the Strasbourg SG data which are a combination of records from two instruments: the SG TT070-T005 recording from 1987 to 1996 and the SG C026 recording from 1996 to 2011. In order to check for any temporal variations in the FCN period, we analyze the data using a running time interval with a shifting of 1 year. The length of each data segment should respect the frequency resolution required for separation of the diurnal waves (Q1, O1, P1, K1, Ψ1,

Figure 1 .

 1 Figure 1. Differences in Ψ1 amplitude factors when the Ψ1group contains or not the FCN for 3 SG

Figure 2 .

 2 Figure 2. Time-variations of amplitude factors for K1, Ψ1, M2 and relative ratios to O1 at Strasbourg

Figure 3 .Figure 4 .

 34 Figure 3. Time variations of FCN resonant period obtained with Bayesian (stars) and LSQ (dots)

  comparison with Fig. 3. In agreement with the above discussion, the time variations of δ(ψ1) reflect similar changes in the resulting FCN period. The observations at other stations are not as long as those at Strasbourg station, so we can only do the comparison in a shorter common period. The data in Brussels can be compared with the first part of the Strasbourg data from 1987 to 2000. The results we obtained for Brussels are in agreement with the time fluctuations observed by Xu and Sun (2009) with smaller error bars as we have used longer time windows. The variation amplitude of δ(ψ1) is very large, but the decreasing and increasing period is similar compared to Strasbourg station. The determination of tidal waves could be influenced by the position, the observation environment of stations, the instrument state and so on, which possibly lead to a different behavior for same tidal wave in different stations. Some variability may occur due to ocean tides effects, particularly at Brussels and Membach that are closer to the sea than other European sites considered here. The data from Moxa, Membach, and Metsahovi SG stations cover the period from 2001 to 2009 and can be compared with the second part of Strasbourg records. Between 2002 and 2008, a roughly similar fluctuation existsamong the results at Strasbourg, Moxa and Membach. Some values in Membach deviate from the main trend, so that the similarity is not that obvious. The situation in Metsahovi station resembles the one at Brussels site in which the amplitude variation is large but the fluctuation period matches the results obtained at Strasbourg. The last two plots are the results obtained at Canberra and Cantley stations, which are located on different continents. The available data length of these two stations was relatively short; therefore, the fluctuation is only partly displayed. The variation at

Figure 5 .

 5 Figure 5. FCN amplitudes (dots) fitted on 6-year segments shifted by a year of VLBI observations with a

Figure 6 .

 6 Figure 6. Time variations of (a) FCN resonant period inverted from VLBI nutation data after applying

  )/δ(O1) δ(M 2)/δ(O1)

Figure 7 .

 7 Figure 7. Variations of LOD, residuals after removing the 22-year and inter-annual terms,

Table 1 .

 1 Basic information on SG stations and data used in this study

	No.	Station	Location Country (Latitude, Longitude)	Data period y-m-d to y-m-d
	1	Brussels	Belgium (50.7986 N, 4.3581 E)	1987-01-01 to 2000-09-22
	2	Cantley	Canada (45.5850 N, 284.1927 E) 1997-07-01 to 2007-12-31
	3	Canberra Australia (35.3206 S, 149.0077 E) 1997-07-01 to 2007-04-18
	4	Membach Belgium (50.6093 N, 6.0066 E)	1995-08-04 to 2010-05-31
	5	Metsahovi Finland (60.2172 N, 24.3958 E)	1994-08-11 to 2009-11-30
	6	Moxa	Germany (50.6447 N, 11.6156 E) 2000-01-01 to 2010-04-30
	7	Strasbourg France (48.6217 N, 7.6838 E)	1987-07-11 to 2011-07-31

Acknowledgments

This work was supported by the Major State Basic Research Development Program of China (2014CB845902), National Natural Science Foundation of China (41321063, 41274085 and 41304058). We are very thankful to the Global Geodynamics Project and the International VLBI Service for Geodesy and Astrometry for providing SG and VLBI data. We would like to thank two anonymous reviewers for their comments to improve the manuscript.


R Rigo-Mariani

B Sareni
email: sareni@laplace.univ-tlse.fr

X Roboam
email: roboam@laplace.univ-tlse.fr

FAST POWER FLOW SCHEDULING AND SENSITIVITY ANALYSIS FOR SIZING A MICROGRID WITH STORAGE

Keywords: Microgrid, storage, optimal scheduling; linear programming, dynamic programming

This paper proposes a fast strategy for optimal dispatching of power flows in a microgrid with storage. The investigated approach is based on the use of standard Linear Programming (LP) algorithm in association with a coarse linear model of the microgrid. The resulting computational time is compatible with simulations over long periods of time allowing the integration of seasonal and stochastic features related to renewable energies. By using this fast scheduling strategy over a complete year of simulation, the microgrid cost effectiveness is considered. Finally, a sensitivity analysis is carried out in order to identify the most influent parameters that should be considered in a sizing loop. Different microgrid configurations are also investigated and compared in terms of cost-effectiveness.

INTRODUCT ON

With the growing number of renewable energy sources the power grid topology has evolved and it could be now described as an aggregation of several microgrids both consumer and producer [START_REF] Celli | Meshed vs. radial MV distribution network in presence of large amount of DG[END_REF]. For those "prosumers", a classical strategy consists in selling all the highly subsidized production at important prices while all consumed energy is purchased [START_REF] Campoccia | Feed-in Tariffs for Gridconnected PV Systems: The Situation in the European Community[END_REF]. Smarter operations become possible with the development of energy storage technologies and evolving price policies [START_REF] Yeleti | Impacts of energy storage on the future power system[END_REF]. Those operations would aim at reducing the electrical bill taking account of consumption and production forecasts as well as the different fares and possible constraints imposed by the power supplier [START_REF] Colson | A Review of Challenges to Real-Time Power Management of Microgrids[END_REF]. The microgrid considered in the paper is composed of a set of industrial buildings and factory with a subscribed power P s of 156 kW and a PV generator with a peak power of 175 kW (Fig. 1a). A 100 kW/100 kWh storage consisting in the association of ten high-speed flywheels is also introduced. The strategy chosen to manage the overall system is based on a daily off-line optimal scheduling of power flows for the day ahead. Then, in real time, an on-line procedure adapts the same power flows in order to correct errors between forecasts and actual measurements [START_REF] Rigo-Mariani | Offline and On-line Power Dispatching Strategies for a Grid Connected Commercial Building with Storage Unit[END_REF]. Several algorithms have been investigated in previous works to perform the off-line optimization for a single day but the high computational times observed did not comply with a sizing procedure that would require many runs of the procedure over long periods of time (e.g., weeks, months, years) [START_REF] Rigo-Mariani | A fast optimization strategy for power dispatching in a microgrid with storage[END_REF]. The present study focuses on a faster approach consisting in two steps. Firstly, a basic Linear Programming (LP) algorithm solves the cost minimization problem with a coarse linear model of the system as in [START_REF] Nottrott | Storage dispatch optimization for gridconnected combined photovoltaic-battery storage systems[END_REF]. Then, a second procedure adapts the obtained solutions to comply with the requirements of a finer nonlinear model. The paper is organized as follows. The first section describes both coarse and fine models used to represent the system and the various considered hypotheses especially concerning the cost model. Then, the second section presents the fast optimization approach and gives details about adaptation of the control references resulting from the LP optimization. In section 4, the results obtained on a test day are presented, by considering particular production and consumption forecasts and according to given energy price policies and subscribed power. A sensitivity analysis is then performed in order to estimate the most significant parameters that could be considered in a sizing procedure for a microgrid with storage.

POWER FLOW MODEL AND DEGREES OF FREEDOM

As illustrated in Fig. 1, the components are connected though a common DC bus. Voltages and currents are not considered so far. The micro grid sizing (cable length) is very limited. Thus losses within the lines can be aggregated with converters efficiencies. Furthermore, the paper focuses on the optimal scheduling of the system without considering a real-time management strategy (i.e. voltage/current control). The study refers to the optimization of active power flows P i (t). Due to the grid policy, three constraints have to be fulfilled at each time step t. The power flows thought the meter have to remain unidirectional (i.e. P 1 (t) 0 and P 6 (t) 0). In a addition, P 6 (t) 0: to avoid illegal use of the storage device, it cannot discharge itself through the production meter

The equations between all power flows are generated using the graph theory and the incidence matrix [START_REF] Bolognani | A linear dynamic model for microgrid voltages in presence of distributed generations[END_REF]. As illustrated in Fig. 1b, three degrees of freedom are required to manage the whole system knowing production and consumption:

• P 5 (t) = P st (t): the power flowing from/to the storage unit (defined as positive for discharge power)

• P 6 (t): the power flowing from the PV arrays to the common DC bus

• P 9 (t) = P PV denotes the possibility to decrease the PV production (MPPT degradation) in order to fulfill possible grid constraints.

The load control is not investigated in the study but could also be considered in further works as an additional degree of freedom with the possibility of shedding or delaying some consumptions.

EFF C ENC ES AND F NE MODEL

A "fine model" is defined taking account of efficiencies of power converters (typically 98 %) and storage losses. These losses are computed with the state of charge SOC (in %) and the power P st using a function P loss (SOC) and calculating the efficiency with a fourth degree polynomial η FS (P st )

(see [START_REF] Celli | Meshed vs. radial MV distribution network in presence of large amount of DG[END_REF]) depending on the direction of the power flow P 5 (i.e; charge or discharge conditions). Both P loss and η FS functions are extracted from measurements provided by the manufacturer [START_REF] Rao | Engineering Optimization[END_REF].

( ) (
)

) ( ) ( ) ( ) (
× + = (1)
Once the overall efficiency is computed, the true power P FS associated with the flywheel is calculated as well as the SOC evolution using the maximum stored energy E FS (in kWh), the time step t (typically 1 hour for the off-line optimization) and the control reference

P 5 . 100 ) ( ) ( ) ( × ∆ × - = ∆ + t E t P t SOC t t SOC FS FS (2)
Another coefficient K FS (in kWh/h) is also introduced to estimate the self-discharge of the flywheels when they are not used (i.e. P 5 = 0):

100 ) ( ) ( × ∆ × - = ∆ + t E K t SOC t t SOC FS FS (3)

DEF N T ON OF A COARSE L NEAR MODEL

In a second step, a coarse model is developed in order to speed up the solving by using a linear formulation of the problem. Firstly, the converter efficiencies are neglected which leads to the following simplifications: In this model, storage losses are also neglected as well as the self-discharge. Thus, the SOC is simply computed at each time step of duration t = 1h:

= = = = ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) (
100 ) ( ) ( ) ( 5 × ∆ × - = ∆ + t E t P t SOC t t SOC FS (5)

A FAST OPT M ZAT ON BASED ON LP

STANDARD COST M N M ZAT ON

The microgrid has to be run in order to minimize the electrical bill expressed as the différence betwen the purchased energy at a rate C p (t) (€/kWh) and the sold one at a rate C s (t) (€/kWh). If the power flowing through the consumption meter is greater than the subscribed power P s , penalties have to be paid with C ex (€/h) [START_REF] Rigo-Mariani | A fast optimization strategy for power dispatching in a microgrid with storage[END_REF]. The references of the power flows associated with the degrees of freedom over the simulated period are computed in a vector P ref = [P 5 P 6 P 9 ]. Once P ref is determined, all other power flows are computed (by means of linear relationships) from forecasted values of consumption and production. Then P 1 and P 11 are known and the cost function is calculated as follows for a 24 hours simulated period:

> = < = + - = = ) ( if 1 ) ( ) ( if 0 ) ( with

LP APPL ED TO A COARSE MODEL AND REFERENCES CORRECT ON

In previous works, several algorithms have been used to perform the off-line optimization for the fine model [START_REF] Rigo-Mariani | A fast optimization strategy for power dispatching in a microgrid with storage[END_REF]. In this section, the standard fast LP is considered to solve the problem related to the coarse model developed in section 2. Such linear approach aims at decreasing CPU time. Using LP imposes to have a linear cost (expressed with a matrix C L ) and linear constraints (expressed though a matrix A and a vector B). Taking account of the exceeding of the subscribed power is strongly nonlinear. Thus, a transformation of the linear problem has to be performed and it leads among other to add another variable denoted as . to the vector with the references. The procedure is then run according to [START_REF] Rao | Engineering Optimization[END_REF]:

× - - = × + × - = × ex C C ) ( ) ( (10)
The constraint matrix A and vector B are built by concatenating the matrices A i and B i used to express each grid requirement [START_REF] Rigo-Mariani | A fast optimization strategy for power dispatching in a microgrid with storage[END_REF]. Considering the penalties due to the overshoots when P s is exceeded imposes the values for . That is performed using the folowing constraint at each time step where M denotes the maximum expected value for (P 1 -P s ) (set to 250 kW here), I n and 0 n are the n×n identity and zero matrices. The dispatching problem defined in the previous subsection can quickly be solved in less than one second (for 24h to be planned) using a standard LP algorithm e.g. the Matlab© function linprog with sparse matrices. Some preliminary results show that the obtained solutions obviously do not comply with requirements of the fine microgrid model. Fig. 2 illustrates a case for which the solution P * obtained with LP is simulated with fine model equations. It should be noted that a deep discharge occurs at around 10 p.m. Due to losses, the SOC goes down to -25 % with the fine model while it remains to 0 % with the coarse linear model. Indeed, taking account of losses also leads to slow down the storage charge and to speed up the storage discharge. Therefore, the control references (P * ) relative to the degrees of freedom obtained with the LP have to be adapted in order to comply with the fine microgrid model. This is performed by means of a step by step correction which aims at minimizing the cost while aligning the SOC computed from the fine model with the one resulting from the LP optimization [START_REF] Rigo-Mariani | A fast optimization strategy for power dispatching in a microgrid with storage[END_REF] Typically, the CPU time related to this local correction procedure is less than one second over a day of simulation while the other tested global algorithms (i.e. genetic algorithm, dynamic programming and trust region) lasted up to 2 h for one day of simulation.

INPUT DATA

In this subsection, the algorithm is tested on a single day considering two different values for the subscribed power P s. The consumption profile is extracted from data provided by the microgrid owner while the production estimation is based on solar irradiation forecasts, computed with a model of PV arrays [START_REF] Darras | Sizing of photovoltaic system coupled with hydrogen/oxygen storage based on the ORIENTE model[END_REF]. Energy prices result from one of the fares proposed by the French main power supplier [12] increased by 30%. Thus, the purchase cost C p has night and daily values with 0.10 €/kWh from 10 p.m. to 6 a.m. and 0.17 €/kWh otherwise. Sale fare C s is set to 0.1 €/kWh which corresponds to the price for such PV plants. C ex = 14 €/h and the prices related to the subscribed power is not considered here. It will be significant when a whole year will be simulated in order to find an optimal value (see section 5).

VAL DAT ON ON A S NGLE DAY

The off-line optimization is performed with the developed algorithm for two different values of P s (i.e P s = 156 kW and P s = 95 kW). The obtained results with storage (WS) are shown in the Table I and compared with cases without storage nor optimal power dispatching (NS). On a single day, the computational time is less than one second. For the two values of P s adding a storage device and an optimal dispatching strongly decreases the cost on the considered day. Results show that there is no sold energy in order to favor the self-consumption. With a lower subscribed power, the number of overshoots is also reduced and the penalty cost goes down to 14 € instead of 140 €. Profiles for the SOC and the optimal power flowing through the consumption meter (i.e.

S Z NG LOOP

SENS V TY ANALYSYS

The developed algorithm has been previously tested on a single day. This section focuses on simulations of a whole year taking the stochastic and seasonal features of PV production and load consumption into account. This "long term vision" is essential to investigate a sizing loop aiming at estimating the best configuration for the microgrid depending on the various costs (Fig. 5). To estimate yearly costs, component prices have to be considered. PV arrays are associated with a cost of 0.6 €/Wc [START_REF] Reichelstein | The prospects for cost competitive solar PV power[END_REF] and the flywheels with 1500 €/kWh [14]. A lifetime of 20 years is expected for the installation. The yearly electrical bill also depends on the value of P s with a cost of 35.3 €/kW [12]. Thus, given yearly forecasted consumption and production and for predeterminated prices policy, two variables are identified as the parameters to define the size of the microgrid: E FW (kWh) the maximum stored energy in the flywheels and

Optimal power dispatching (1 year)

Sizing loop

Environnement: Prices, solar radiation, consumption Size: PV, Storage, P s Architecture: Microgrid Fig. 5: Optimal dispatching in a sizing loop P PV_max (kW) the maximum power of the PV generator. Additional management parameters are also introduced: the subscribed power P S (kW) and the number of successive days T schedule on which the scheduling is extended. For T schedule = 1 day as previously, 365 successive optimizations have to be performed to simulate a year. It should be noted that the computational time of the power flow scheduling strongly increases when the simulated period becomes longer as the number of variables in the LP problem is greater. The last studied variable denoted SOC 0 represents the initial and final SOC value in the optimization loop (50 % previously). Bounds are arbitrary chosen for those five variables:

• 0 kWh < E FW < 1500 kWh • 0 kW < P PV_max < 500 kW

• 100 kW < P s < 250 kW • 1 day < T schedule < 120 days • 0 % < SOC 0 < 100 % A sensitivity analysis is performed using a full factorial design [START_REF] Antonio | Design of experiments for sensitivity analysis of voltage sags variables, 15 th Harmonics and quality of Power (ICHQP)[END_REF]. The weight coefficient of the main effects (a i ) of the variables and their interactions (b i ) are studied as for the following example with a two parameter function y = f(x 1 ,x 2 ) modelled as follow:

2 1 1 2 2 1 1 . b a ˆx x x x a y y × + × + × + = (12)
In (12), is the average value of the function with the different points y i given by the design of experiments when the values are at their higher (+1) or lower (-1) bounds. The weight a i and b i are then computed using the columns of the

CASE ANALYS S

Table III gives the results returned when a year is simulated with the optimal power dispatching with different values of the three significant sizing parameters. As in the previous section, T schedule is set to 1 day and SOC 0 to 50 %. The resulting computational time is close to 10 min for each investigated sizing (note that ~3 days of computation would be necessary to optimize power management by means of dynamic programming algorithm). The costs for the components and the purchase/sold energy are expressed in k€. The second column refers to a case for which only the consumption is considered. Adding a storage device allows to reduce the penalty by restricting as much as possible the exceeding of P s . In the same time, the cost for purchased energy is lowered with load shifting during night hours when prices decrease. But the resulting gain does not compensate the investment cost. With a low investment cost and a good retribution price, a solar generator improves the yearly results of the microgrid. Adjusting the subscribed power also improves the yearly cost and an optimal value has to be found.

CONCLUS ONS

The study carried out in this paper aims at proposing a fast procedure in terms of computation time that could be used to investigate costeffectiveness of a microgrid with storage. In previous works, efficient algorithms have been developed to perform the daily scheduling of power flows. However, the main drawbacks of these methods reside in their computational times that become prohibitive if the microgrid has to be simulated over a long period of time. To overcome this problem, a fast optimization approach based on LP has been proposed. This approach consists in two successive steps. Firstly, a coarse linear model of the microgrid is exploited to solve the optimal dispatching with a classical LP algorithm. Secondly, control references optimized with the coarse model are adapted in order to comply with a finer model of the microgrid which takes account of nonlinear features (i.e. efficiencies). The performance of this approach with regard to energy cost minimization and computational time reduction has been shown on a particular test day. Moreover, the fast CPU time resulting from this optimal dispatching method has allowed us to simulate the microgrid over a whole year and to investigate several configurations. The obtained results have shown that the interest of using a storage unit is closely linked to the economical context. Future studies will be focused on the same issue with other kind of storage technologies such as Li-ion batteries for which cycling effect would have to be included in the cost function. Finally, the fast control algorithm may offer the ability of achieving systemic design of microgrids integrating sizing optimization loop with power dispatching optimization by taking account of system environment and requirements.

  upper (u b ) and lower (l b ) bounds of decision variables are expressed using the column vector J n with n coefficients equal to 1, where n is the number of simulated time steps, i.e. n = 24 for a whole day with t = 1h. In particular, the limits of P 5 refer to the maximum charge and discharge powers of the storage with P st_min = -100 kW and P st_max = 100 kW. previous cost function C(P ref ) is developed for the coarse model according to the decision variables P 5 , P 6 , P 9 and for the linear problem. The nonlinear term is removed to obtain the vector C L used in the LP optimization.

OptimalFig. 2 :

 2 Fig. 2: SOC constraint violation with the fine model ( )

Fig. 3 :

 3 Fig. 3: Forecasted consumption and production

Fig. 4 :

 4 Fig. 4: Optimal profiles with two P s

Fig. 6 :

 6 Fig. 6: Sensitivity analysis -a) factor effect -b) interaction effect

MODEL OF THE STUD ED M CROGR D

	grid P	P 1			load
	Consumption	P 2		Loads
	Meter	P 3	P 4		st P
		P 6		P 5	Storage
	P 11		P 7	P 8	PV P
	Meter Production	P 10			∆ PV P -P 9
	DC Bus
	Degrees of freedom	Degradation of the
				PV production

P Fig. 1: Considered microgrid and power flows 2.

Table II

 II

										4.10 4	E FW	P PV_max	P s
										T schedule
										2.10 4	SOC 0
										0
										4.10 4	P PV_max .T schedule
										2.10 4	P PV_max .E FW	E FW . P s	E FW .SOC 0 P s . SOC 0 SOC 0 .T schedule
										0
										P PV_max .P s	P PV_max .SOC 0 E FW .T schedule	P s .T schedule
										and N
	the number of experiments (i.e. 4 in a problem with
	two parameters):
		b a i i	= =	b a	i i	T T	.y .y	/ /	N N	(13)
	II. Full factorial design of experiments
	a 1 =x 1	a 2 =x 2				b 1 =x 2 x 1	Y
	+1	+1						+1	y 1
	+1		-1						-1	y 2
	-1	+1						-1	y 3
	-1		-1						+1	y 4
	For five parameters, as for the studied problem, the
	absolute values of coefficients are plotted in Fig 6.
	The corresponding parameters (effect/interactions)
	are underlined. The coefficients referring to the
	values of E FW , P PV_max , and P s appear to be the most
	influent. In a first approximation, only those three
	variables could be considered when studying
	different sizing cases.

ACKNOWLEDGMENT

This study has been carried out in the framework of the SMART ZAE national project supported by ADEME (Agence de l'Environnement et de la Maîtrise de l'Energie). The authors thank the project leader INEO-SCLE-SFE and partners LEVISYS and CIRTEM.