A Yedjour

B A Van Tiggelen
email: bart.van-tiggelen@grenoble.cnrs.fr

Diffusion and Localization of Cold Atoms in 3D Optical Speckle

Keywords: 72, 15, Rn Localization effects -67, 85, Hj Bose-Einstein condensates in optical potentials -3, 70, Jk Atoms in optical lattices

In this work we re-formulate and solve the self-consistent theory for localization to a Bose-Einstein condensate expanding in a 3D optical speckle. The long-range nature of the fluctuations in the potential energy, treated in the self-consistent Born approximation, make the scattering strongly velocity dependent, and its consequences for mobility edge and fraction of localized atoms have been investigated numerically.

Introduction

Anderson localization is by now a phenomenon that has been widely investigated, both theoretically and experimentally, and for many different kinds of waves, from electrons, to electromagnetic waves, ultrasound [1], and cold atoms [2]. To understand Anderson localization and to provide quantitative predictions for experiments, many tool models have been proposed, among which the Anderson model is undoubtedly the best known. This model describes a noninteracting and electron, tightly bound to the nucleus, but capable to tunnel to nearby atoms. Already in the celebrated 1958 paper [3], Anderson demonstrated how this model highlights the role of dimensionality. A genuine mobility edge only occurs in dimensions larger than 2 [4]. For classical waves, a few observations on 3D localization have been reported [5,6].

The tight binding model is highly relevant to understand 3D dynamical localization of cold atoms [7]. However, is not appropriate to describe localization of many other waves, where the starting point is much more a diffuse, extended motion, rather than a tightly bound state. Localization of electromagnetic waves in 3D disordered media for instance, is much more a diffusion problem than a problem of tunneling to nearest neighbors. The scaling theory of localization [8], as well as the Thouless criterion [9], both deal with conductance, and thus use the diffuse motion as a starting point.

The self-consistent theory, first formulated by Vollhardt and Wölfle in 1981 for 2D electron conductivity [START_REF] Vollhardt | Selfconsistent theory of Anderson Localization[END_REF], was the first work that explicitly calculated how quantum corrections affect the classical "Drude" picture, to make way for localization. Despite its evident perturbational na-ture, the theory has been successful because it provides a microscopic picture of finite-size scaling, reproduces the Ioffe-Regel criterion for the mobility edge in 3D, and locates the mobility edge of the tight-binding model quite accurately [START_REF] Kroha | [END_REF]. The aim of the present work is to revisit and apply this theory to the localization of cold atoms in optical speckle.

The first experiments on 1D cold atom localization have been carried out recently [12]. The atoms are released from a BEC and subsequently expand in a potential energy landscape created by optical speckle, supposed free of mutual interactions. Both theory [13] and experiment have revealed the presence of a quasi-mobility edge in 1D. Atoms with velocities v > h/mξ (ξ is the correlation length of the disorder) can hardly be scattered because this would imply a momentum transfer larger than h/ξ which the random speckle cannot support. As a result the localization length is infinite in the Born approximation, though finite and large when all orders are taken into account. This somewhat surprising result highlights the impact of long-range correlations in 1D. In higher dimensions, small angle scattering can still lead to small enough momentum transfer to be transferred to the speckle, even for large velocities, so that this quasi-mobility edge does not occur. Yet, correlations are expected to affect localization, since the potential field sensed by the atom strongly depends on its velocity, and strong forward scattering is not favorable for localization to occur. In addition, near the 3D mobility edge the disorder is necessarily large so that the spectral function of the atoms is not strongly peaked near energies E = p 2 /2m, neither has it a Lorentzian broadening.

Self-consistent Born Approximation

In the following we consider the scattering of a noninteracting atom with energy E and momentum p from a disordered potential V (r). Two properties are specific for an optical potential. Firstly, the fluctuations δV (r) are determined by the optical intensity and not by the complex field. This means that they are not Gaussian but rather Poisonnian. As a result, the two-point correlation will in principle not be sufficient to describe the full scattering statistics. Secondly, the correlation function, given by δV (r)δV (r ′ ) = U sinc 2 (∆r/ξ), is long range. Here U = V 2 depends on the average optical intensity. It is not difficult to see that the Born approximation breaks down at energies E U/E ξ [14], with E ξ = h2 /2mξ 2 an important energy scale related to correlations. As usual we expect matter localization to occur at small energies, near the band edge of the spectrum. Another consequence of the long-range correlations is that scattering strongly depends on the De Broglie wave length and thus on the velocity of the atom. This makes it impossible to define a mean free path ℓ in the usual way, that is from the exponential decay of the ensemble-averaged Green function [15].

In the following we shall cope with the second problem. The ensemble-averaged Green function is written in terms of a complex self-energy as [15]. We shall apply the Self-consistent Born Approximation (SCBA) according to which the complex self-energy Σ(E, p) of the atom is calculated from [START_REF] Cusack | The Physics of Structurally Disordered Matter[END_REF] Σ

G(E, k) = [E -p 2 /2m - Σ(E, p)] -1
(E, k) = k ′ U (k -k ′ ) E -k ′2 -Σ(E, k ′ ) (1)
In this equation, k ≡ d 3 k/(2π) 3 , and from now, all energies, including Σ(E, p), are expressed in units of the energy scale E ξ . Momenta are expressed as p = hk with the De Broglie wave number k expressed in units of 1/ξ. In Eq. (1), U (kk ′ ) represents the structure function associated with the speckle correlation, which determines the angular profile in single scattering [14]. The SCBA is convenient because its imaginary part expresses the generalized optical theorem in single scattering [15]. In addition, it avoids the bound state at negative energies predicted by the first Born approximation. In Ref. [START_REF] Skipetrov | [END_REF] the SCBA was solved analytically for cold atoms and zero-range correlations. Equation ( 1) has been solved by iteration, with spline interpolation between 500 points 0 < k n < 3. The angular integral can be performed analytically. Typically 10-20 iterations have been necessary to ensure good convergence.

In Figure 1 we show real and imaginary part of Σ(E, k) for U/E 2 ξ = 1 and E = 0, and compare it to the first order Born approximation (FBA) applied in Ref. [14]. As a realistic experimental reference we take 87 Ru-atoms released from a BEC with chemical potential µ = 219 h Hz into an optical speckle with correlation length ξ = 0, 26 µm. This reveals that µ and E ξ are equal energy scales in typical experiments. Equivalently, the De Broglie wavelength and the correlation length are competing length scales, λ/2π ≥ ξ. To discriminate "trivial trapping " in deep random potential wells from "genuine" Anderson localization, experimentalists wish to arrange the experiment such that the typical kinetic energies µ are somewhat larger than the typical fluctuations √ U in the potential energy [12]. In that case, U/E 2 ξ ≥ 1. We will comment on this choice later. U/E 2 ξ = 1, as a function of wave number (in units of 1/ξ) for an energy E = 0, calculated in the self-consistent Born approximation (SCBA). The dashed line denotes the first Born approximation (FBA). To compensate for the shift in the band edge (E b /E ξ = -0.85) predicted by the SCBA, the FBA has been evaluated at E/E ξ = 0.85.

0,0 0, 5 1,0 1,5 2,0 2,5 3,

Figure 1 (top)

shows -Im Σ(0, k) to be nonzero only for k < 3. It is also seen that the FBA significantly overestimates the amount of scattering. The real part Re Σ(0, k) is clearly negative. This shifts the band edge of the energy spectrum to E b = -0.85. The energy spectrum has a typical lower bound -√ U , but the SCBA always locates the band edge at somewhat larger energies. The SCBA does probably not treat the (small) density of states near E ≈ -√ U very well, where sharply localized Lifshitstype states are likely to occur. Figure 2 shows the distribution of wave numbers at energy E = 0, expressed by the spectral function S(E, k) ∼ -Im G(E, k). It is a rather broad distribution, with a tail extending to k = 2. This is important, since we will see in the next section that atom transport is quite sensitive to large momentum transfers, involving large k-vectors. At low energies, relevant for localization, h/ξ has become the typical momentum of an atom, and E ξ the typical energy. At larger energies E = 3, the spectral function behaves normally, i.e. strongly peaked near k = √ E. The peak is smaller because we did not conclude the geometric 4πk 2 surface factor in phase space, so to highlight its weight at small k < 1 ("slow atoms") for later purposes.

Bethe-Salpeter equation

We proceed with the calculation of the diffusion constant of the cold atoms, and the possible presence of a mobility edge where it vanishes. With that information we will find how many atoms will be localized. The idealized model we consider is schematically drawn in Figure 3.

The Bethe-Salpeter equation is a rigorous equation for the two-particle Green function [START_REF] Mahan | Many-Particle Physics[END_REF]. In phase space this object is written as Φ kk ′ (E, t, r), which is readily interpreted as the"quantum probability density" for an atom

E = 0 U 1/2
Fig. 3. Model of the atomic expansion. The BEC releases noninteracting atoms with kinetic energies between 0 and its chemical potential µ, that penetrate the disordered speckle potential, with average and its typical fluctuations equal to √ U . In this work localized states (red) are found only for energies below the average.

with velocity hk/m to travel, during the time t, from position r = 0 to position r and to achieve the velocity hk ′ /m. Its Fourier-Laplace transform with space-time is written as Φ kk ′ (E, Ω, q). It has two fundamental properties, namely reciprocity, Φ kk ′ (E, Ω, q) = Φ k ′ k (E, Ω, q), and normalization. The last property can be expressed by,

k Φ kk ′ (E, Ω, q = 0) = -2Im G(E, k ′ ) -iΩ (2)
In particular, when also integrating over all energies,

∞ -∞ dE 2π k,k ′ Φ kk ′ (E, Ω, q = 0) = 1 -iΩ ∞ -∞ dE -1 π Im G(E, k) = 1 -iΩ (3)
The last equality, needed for later purposes, follows from a general sum rule of the spectral function [START_REF] Mahan | Many-Particle Physics[END_REF]. This identity guarantees that the total quantum probability for the atom to be somewhere, with some velocity and with some energy, is conserved, and equal to one. The Bethe-Salpeter can be re-written as a quantum-kinetic equation for Φ kk ′ (E, Ω, q). The conservation of quantum probability guarantees the existence of a hydrodynamic diffusion pole. We shall express this as

Φ kk ′ (E, Ω, q) = 2 k -Im G(E, k) φ(E, k, q)φ(E, k ′ , q) -iΩ + D(E)q 2 (4) where φ(E, k, q) = -Im G(E, k) -i(k • q)F (E, k) + O(q 2 ).
This expression states that the distribution of atoms in k-space at given energy E is essentially governed by the spectral function, with a small correction that supports a current. The front factor in Eq. ( 4) is imposed by the normalization condition (3). The k-integral of -Im G(E, k) is recognized as (π times) the density of states (DOS) per unit volume ρ(E). With the correct normalization, we can set Ω = 0. The still unknown function F (p) follows from [START_REF] Mahan | Many-Particle Physics[END_REF] F

(E, k) = |G(E, k)| 2 - ∂Re G(E, k) ∂k 2 + |G(E, k)| 2 k ′ ∂Re G(E, k ′ ) ∂k ′2 k • k ′ k 2 U kk ′ (E, 0) + |G(E, k)| 2 k ′ F (E, k ′ ) k • k ′ k 2 U kk ′ (E, 0) (5)
The irreducible vertex U kk ′ (E, q) generalizes the function U (kk ′ ) defined in the first section to all interference contributions in multiple scattering. Once we have solved for F (E, k), the diffusion constant follows from the Kubo-Greenwood formula [START_REF] Mahan | Many-Particle Physics[END_REF],

D(E) = h m 2 3 1 πρ(E) k k 2 F (E, k) (6)
Note that D(E) is determined by the fourth moment of the distribution F (E, k), which puts a large weight on "fast" atoms. The order of magnitude of the diffusion constant is governed by the ratio h/m of Planck's constant and the mass of the atom, the second factor being dimensionless and of order unity at low energies. For 87 Ru, h/m ≈ 1800 µm 2 /s. The third term in Eq. ( 5) can be transformed using the exact Ward identity,

Σ(E, k + 1 2 q) -Σ * (E, k - 1 2 q) = k ′ U kk ′ (E, q) × G(E, k ′ + 1 2 q) -G * (E, k ′ - 1 2 q) (7)
If this identity is developed linearly in q, and inserted into Eq. ( 5) , we obtain,

F (E, k) = F 0 (E, k) + δ q U (E, k) + |G(E, k)| 2 k ′ F (E, k ′ ) k • k ′ k 2 U kk ′ (E, 0) (8)
with

F 0 (E, k) ≡ |G(E, k)| 2 1 + ∂Re Σ(E, k) ∂k 2 - ∂Re G(E, k) ∂k 2 δ q U ≡ |G(E, k)| 2 k ′ Im G(E, k ′ ) 2i k k 2 • ∂ ∂q U kk ′ (E, q)
. Three levels of analysis exist. First, in the Drude approximation, one neglects all contributions form the BSequation and one adopts F (E, p) = F 0 (E, p), including the wave number derivatives. The Drude diffusion constant D d (E) can be used to define a dimensionless Ioffe-Regel type parameter from the relation

D d (E) = 1 3 (2h/m)kℓ. Hence, kℓ ≡ k k 2 F 0 (E, k) k -Im G(E, k) (9)
For a short-range correlation, the self-energy is independent of k so that F 0 (E, k) = 2Im 2 G(E, k), and this definition of kℓ coincides with the usual one in terms of Im Σ(E) [15]. For low energies we found in Figure 3 that typically k ≈ 1/ξ . If we anticipate that kℓ ≈ 1 near the mobility edge, we conclude that the mean free path is roughly equal to the correlation length.

-0,5 0,0 0,5 1,0 1,5 2,0 0,0 0,5

Energy / E Drude Boltzmann/Drude DOS Ioffe-Regel Fig. 4. Solution of the BS equations ( 6) and ( 8) in the Boltzmann approximation, for a disorder amplitude U/E 2 ξ = 1. Shown as a function of energy are the Drude conductivity πρ(E)D d (E), the Ioffe-Regel parameter (9), the DOS that vanishes at the band edge E/E ξ = -0.83, and the ratio of Boltzmann and Drude conductivity. Note that 2/3 times the Ioffe-Regel parameter equals the Drude diffusion constant, expressed in units of h/m. The Drude approximation is popular in electron -impurity scattering but clearly inadequate when the scattering is strongly anisotropic, as for the optical speckle. In the Boltzmann approximation we adopt U kk ′ = U (kk ′ ), i.e. the structure function associated with the optical disorder. Being a function of k-k ′ only, it follows that δ q U = 0. Different results are summarized in Figure 4. At low energies Boltzmann and Drude diffusion constant typically differ by a factor 1.5 as was obtained by Kuhn etal [14] on the basis of the FBA. For atom energies E > E ξ the Boltzmann 8), rescaled properly such that both k-integrals equal σB/σ0 diffusion constant rapidly rises since only strong forward scattering can occur. In the region E = 0, the Ioffe-Regel parameter takes values of the order of 1.5. We infer from Figure 5 that the solution F B (E, p) is roughly a rescaling of the Drude Ansatz, for which the current is dominated by atoms with velocities v = 0.83h/mξ. Nevertheless, a non-negligible fraction of atoms faster than v = 1.5h/mξ contributes to the current (note that h/mξ ≈ 7 mm/s for the set-up with Rubidium described above).

We finally consider constructive interferences, and add the most-crossed diagrams to the BS-equation in the spirit of the self-consistent theory of localization. Any observation D < D B must be attributed to constructive interferences. It is well-known that, by reciprocity, these diagrams can be constructed from the solution φ(E, k, q) of the BSequation by removing the incoming and outgoing Green's functions (indicated by a hat), and by time-reversing the bottom line [START_REF] Van Tiggelen | [END_REF]: 4) the two procedures lead to

U MC kk ′ (E, q) = Φ 1 2 (k-k ′ +q) 1 2 (k ′ -k+q) (E, k+ k ′ ). Using Eq. (
Φ kk ′ (E) = 2 σ(E)q 2 × [ImΣ(E, k)ImΣ(E, k ′ ) + iq • (I kk ′ (E)k + I k ′ k (E)k ′ )] U MC kk ′ (E, q) = 2 σ(E)(k + k ′ ) 2 × Im 2 Σ(E, 1 2 ∆k) + I ∆k/2,∆k/2 (E)i(k + k ′ ) • q (10)
We abbreviated

I kk ′ = F (E, k)Im Σ(E, k ′ )/|G(E, k)| 2 , ∆k = |k -k ′ | and introduced σ(E) = πρ(E)D(E)
, the equivalent of DC-conductivity in electron conduction. We now face the more complicated task of solving Eqs. ( 8) and ( 6) simultaneously with U kk ′ = U (kk ′ ) + U MC kk ′ , and of finding out if its extrapolation to small energies leads to a mobility edge where σ(E) = 0. This constitutes a "selfconsistent" problem for the entire function F (E, k), rather than just for its fourth moment, the DC-conductivity.

We first observe that for most-crossed diagrams δ q U = 0. This term does not appear in standard self-consistent theory [START_REF] Vollhardt | Selfconsistent theory of Anderson Localization[END_REF], which relies on moment expansion. Note that it also features "self-consistently" the function F (k)/σ, just like the second term in Eq. ( 8). As can be induced from Eq. ( 10) the singularity at k = -k ′ that generates the (weak) localization is partly compensated by the factor k + k ′ . We shall therefore ignore it here as well.

The self-consistent equations ( 8),( 6) and ( 10) can be solved almost analytically when the self-energy is assumed k-independent, typically true for zero-range correlations. Without more details we mention that the mobility edge then occurs at kℓ = 1.122 [20]. Quite convenient is that, even when scattering extends to infinite wave numbers, our theory does not require an ad-hoc cut-off to eliminate short wave paths that diverge in approximate theories [14,21]. For the speckle correlation, Figure 6 gives the result of the exact numerical solution, obtained by iteration and spline interpolation. This method worked satisfactorily until close to (∆E ≈ E ξ /5) the mobility edge, where the most-crossed diagrams give a diverging contribution. Before that happens, the strong forward scattering of a single scattering competes heavily with the reduction induced by weak localization.

-0,8 -0,6 -0,4 -0,2 0,0 0,2 0,4

Energy / E DOS ( B 0

)(1-K) Fig. 6. Solution of the BS equation ( 8) with inclusion of the most-crossed diagrams for a disorder amplitude U/E 2 ξ = 1. Shown in red is the ratio of conductivity and Drude conductivity for energies E > 0.15E ξ for which our iteration converged. The blue dashed line relies on an approximation discussed in the text and locates the mobility at E = -0.06E ξ .

To find the location of the mobility edge we shall use the following approximation. In Figure 7 we see that at 0,0 0,5

F(E,k) x 1,0 1,5 2,0 2,5 3,0 0,0 0,5 1,0 1,5 2,0 2,5 3,0 3,5 4,0 4,5 k 4
Wave number k Fig. 7. The function k 4 F (E, k) (in red)), solution of the BSequation for the energy E/E ξ = 0.15, compared to the Drude approximation k 4 F0(E, p) defined in Eq. ( 8) (blue dashed line).

The first has been rescaled by the factor σ0/σ such that the k-integrals of both are equal to one.

E/E ξ = 0.15 the solution F (E, k)/σ is closely approximated by F 0 (E, k)/σ 0 (both have their fourth moment normalized to one). This equivalence is physically reasonable since it implies that all atoms with energy E undergo the same reduction in diffusion, but keep the same velocity distribution as found in the Drude picture. If we insert F 0 (E, k) in the left hand side of Eq. ( 8) and integrate over k we can derive the simple relation σ/σ B ≈ 1 -K(E) that is reminiscent of standard self-consistent theory [START_REF] Vollhardt | Selfconsistent theory of Anderson Localization[END_REF]. The parameter K is found to be,

K(E) = - 4 3σ 2 d (E) kk ′ F 0 (E, k ′ ) k • k ′ (k + k ′ ) 2 × Im 2 Σ E, 1 2 |k -k ′ | |G(E, k)| 2 (11)
For U/E 2 ξ = 1, this approximation locates the mobility edge (K = 1) at E/E ξ = -0.06 (0.94 when we shift energies over √ U as in Fig 8). Upon inspecting the numerically exact solution in Figure 8, we suspect that the real mobility edge is located somewhat lower, near E/E ξ = -0.1. For smaller disorder we found that the approximation becomes better. It is in principle possible to obtain numerically the function F (E, p)/σ as σ → 0 and calculate more precisely the location of mobility edge, but this is beyond the scope of this paper. In Figure 8 we show the different energies for different strengths of the disorder. If we apply the criterion found in Ref [14], kℓ B = 0.95 (dashed line in Figure 8), localization would occur at smaller energies, around E/E ξ = -0.35 (0.65 in Fig. 8) for U/E 2 ξ = 1. This approach expresses the general trend well but is clearly somewhat pessimistic, probably because their choice of the ad-hoc wave number cut-off to calculate the mostcrossed diagrams underestimates K. Note that all local- Fig. 8. Band edge and mobility edge as a function of the disorder strength U , all expressed in the energy E ξ . Contrary to the former figures the energy scale has been shifted by the average potential √ U of the potential energy so that a direct comparison can be made to the chemical potential µ > 0, and the criterion Ec -E b ≈ 0.5U/E ξ found by Ref. [14] (blue dashed). The numbers around the mobility edge reveal the small variation of the Ioffe-Regel parameter (9). The exclusion of small energies is likely to be an artifact of the SCBA and we expect strongly localized Lifshits-type states with small DOS. ized states occur below the average of the potential landscape. Unlike in 1D, we find no regime with atoms fast enough to traverse the potential barriers, but to become localized purely by constructive interferences, without the assistance of tunneling.

A final important question is how many atoms will be localized, given an initial velocity distribution that is determined by the expansion of the BEC after eliminating the trap. We emphasize that according to the scenario sketched in Fig. 3, the chemical potential µ of the BEC does not represent the energy of the atom inside the disordered potential, but rather the distribution of incident kinetic energies. If this distribution is denoted by φ µ (k), it follows from Eq. ( 4) that the fraction of localized atoms, regardless of their final velocity or position, is given by

Ec E b dE 2π k,k ′ Φ kk ′ (E, t, q = 0)φ µ (k ′ ) = Ec E b dE k ′ -1 π Im G(E, k ′ )φ µ (k ′ ) (12)
Castin and Dum [22] showed that after the free expansion φ µ (k) ∼ 1 -k 2 /k 2 µ and zero for k > k µ , where the maximum wave number k µ = µ/E ξ . The fraction of localized atoms is thus determined by the number of microscopic states below the mobility edge whose kinetic energies are smaller than µ. In this discussion it is convenient to make the zero of the energy scale the same for µ and E, as was already done already in Figure 8, shifting the localized region to positive energies. In present (1D) experiments is U/E 2 ξ = 1 and µ ≈ E ξ which is only slightly above the 3D mobility edge. One might perhaps expect most atoms to be localized. Our calculations clearly show that the distribution of atom energies is quite different from φ µ (k). This (normalized) distribution F (E) is given by the wave number integral in Eq. ( 10) and is shown in Figure 9. For µ ≪ E ξ it reduces to the spectral function at k = 0, independent of φ µ (k). Even for small µ, many atoms achieve energies E > E c (40 % for U/E 2 ξ = 1) and are delocalized. This number agrees well with predictions based on zero-range correlations in which case 45 % was found to be localized for µ ≪ E c [START_REF] Skipetrov | [END_REF]. The fraction of delocalized atoms further decreases as the chemical potential rises (Figure 10), with only 35 % localized for µ = E ξ . Even if we choose µ > √ U , the atoms that localize have energies E < √ U . In conclusion, we have calculated the phase diagram for localization of cold atoms in a 3D speckle potential, using the self-consistent Born approximation and the selfconsistent theory of localization. The mobility edge is characterized by a Ioffe-Regel type parameter that varies between 1.2 and 1.4. Depending on the chemical potential of the BEC, typically 35 % to 60 % of the atoms are localized. The self-consistent Born approximation deals already much better with the long-range correlations and the broadness of the spectral function than the first Born approximation, but does not discriminate between different statistics of the disorder. Yet, the mobility edge of the tight-binding model is known to depend on that [START_REF] Kroha | [END_REF]. It would be very interesting to apply a recently proposed method [23] to calculate the self-energy of the atoms more precisely. The theory presented in this work can then be used straightforwardly to find the mobility edge.

We would like to thank Sergey Skipetrov for help and discussions.

Fig. 1 .

 1 Fig.1. Imaginary (top) and real part (bottom) of the self energy (in units of E ξ ) of an atom in a speckle potential with U/E 2 ξ = 1, as a function of wave number (in units of 1/ξ) for an energy E = 0, calculated in the self-consistent Born approximation (SCBA). The dashed line denotes the first Born approximation (FBA). To compensate for the shift in the band edge (E b /E ξ = -0.85) predicted by the SCBA, the FBA has been evaluated at E/E ξ = 0.85.

Fig. 2 .

 2 Fig.2. Spectral function, normalized to the total density of states at energy E, calculated in the SCBA, for E = 0 (red) and E = 3 (blue) and for U/E 2 ξ = 1. At small energies it has a large weight at k = 0 (very cold atoms) and extends up to k = 2.

Fig. 5 .

 5 Fig. 5. The function k 4 FB(E, k), solution of the BS-equation for the energy E/E ξ = -0.06 in the Boltzmann approximation. The blue dashed line compares it to the Drude approximation k 4 F0(E, p) defined in Eq. (8), rescaled properly such that both k-integrals equal σB/σ0

Fig. 9 .

 9 Fig.9. Energy distributions of the atoms inside the disordered speckle potential, for different chemical potentials of the BEC from which they were released and U/E 2 ξ = 1. They all exhibit a tail of relatively fast atoms (E > E ξ ) that extends beyond the mobility edge Ec.

5 Fig. 10 .

 510 Fig. 10. Fraction of localized atoms as a function of chemical potential, for different amplitudes of the disorder.


J Charles Gilbert
email: charles.gilbert@inria.fr.

QPAL -A solver of convex quadratic optimization problems, using an augmented Lagrangian approach Version 0.6.1 (May, 2009)

Keywords: augmented lagrangian, convex quadratic optimization, dense and sparse matrices, gradient projection, ℓ-BFGS matrix, linear convergence convergence linéaire, gradient projeté, lagrangien augmenté, matrice ℓ-BFGS, matrices dense et creuse, optimisation quadratique convexe

QPAL is a piece of software that aims at solving a convex quadratic optimization problem with linear equality and inequality constraints. The implemented algorithm uses an augmented Lagrangian approach, which relaxes the equality constraints and deals explicitely with the bound constraints on the original and slack variables. The generated quadratic functions are minimized on the activated faces by a truncated conjugate gradient algorithm, interspersed with gradient projection steps. When the optimal value is finite, convergence occurs at a linear speed that can be prescribed by the user. Matrices can be strored in dense or sparse structures; in addition, the Hessian of the quadratic objective function may have a direct of inverse ℓ-BFGS form. QPAL is written in Fortran-2003.

1 Presentation

Scope of the program

The semi-acronym Qpal (pronounce Q-pal) stands for Quadratic Programming with an Augmented Lagrangian. Qpal can indeed solve a convex quadratic optimization problem with linear constraints. The form of the problem considered by Qpal is the following

(QP)      min f (x) := g ⊤ x + 1 2 x ⊤ Hx l (x, A I x) u A E x = e.
(1.1)

The function f in (1.1) is called the objective of the problem. Its Hessian H must be an n × n positive semidefinite symmetric matrix (hence f is a convex function) and its gradient at the origin is the vector g ∈ R n . The notation l (x, A I x) u expresses in compact form bound constraints on x and on A I x: A I is m I × n and the (n+m I )-dimensional vectors l and u satisfy l < u (i.e., l i < u i for all indices i) and may have infinite components. Hence the inequality constraints read l i x i u i , for i = 1, . . . , n, l n+i (A I x) i u n+i , for i = 1, . . . , m I .

RT n°0377 inria-00442295, version -Dec

Finally, the equality constraints are expressed with an m E × n matrix A E and a vector e ∈ R mE . The matrix A E need not be surjective.

Since H can vanish, (1.1) can also model a linear optimization problem.

Qpal is flexible with respect to the memory representation of the matrices defining the problem. Indeed the Hessian matrix H may be dense, or sparse, or may have a direct or inverse ℓ-BFGS structure. The constraint Jacobians A I and A E may be dense or sparse. Sparse matrices are represented by the row-column indices of the nonzero elements, as well as the values of these elements.

It is convenient to denote by B := {1, . . . , n} the index set of the variables x and by A B = I n the identity matrix of order n. Then the inequality l (x, A I x) u reads in a similar manner l A B∪I x u. We also adopt the notation

m = n + m I + m E and A =   A B A I A E   .
An inequality constraint with the index i ∈ B ∪ I is said to be active at x ∈ X if A i x = l i or u i .

It is said to be inactive otherwise. The active set at x ∈ X is defined by

A(x) = {i ∈ B ∪ I : A i x = l i or u i }.
The feasible set of (1.1) is denoted by

X = {x ∈ R n : l A B∪I x u, A E x = e}.
Qpal is written in ANSI Fortran 2003 (F03 for short), using double precision. The code uses this feature of F03 that makes possible calling a subroutine with a structure argument having components that are not yet allocated when the subroutine is called. Qpal compiles with gfortran [5], version "4.4.0 20090321 (experimental)" or higher.

Detecting optimality

The Lagrangian of problem (1.1) is the function ℓ : R n × R n+2mI +mE → R defined at (x, λ) by

ℓ(x, λ) = g ⊤ x + 1 2 x ⊤ Hx + (λ l B∪I ) ⊤ (l -A B∪I x) + (λ u B∪I ) ⊤ (A B∪I x -u) + λ ⊤ E (A E x -e), (1.2)
in which λ = (λ u B∪I , λ l B∪I , λ E ). Since problem (1.1) is convex, its optimality conditions read: x is a solution to (1.1) if and only if there is a λ ∈ R n+mI +mE such that

             g + Hx + A ⊤ λ = 0 l A B∪I x u (λ + B∪I ) ⊤ (A B∪I x -u) = 0 (λ - B∪I ) ⊤ (A B∪I x -l) = 0 A E x = e, (1.3)
where λ + := max(0, λ) and λ -:= max(0, -λ), componentwise. In the table below, we have shown the correspondence existing between the multipliers used in (1.2) and in (1.3

): multipliers in (1.2) multipliers in (1.3) λ u B∪I λ + B∪I λ l B∪I λ - B∪I λ u B∪I -λ l B∪I λ B∪I λ E λ E

Brief description of the method

In Qpal, problem (1.1) is transformed, using an auxiliary variable y ∈ R mI , into (QP')

       min g ⊤ x + 1 2 x ⊤ Hx l (x, y) u A I x = y A E x = e.
(1.4) Note that an auxiliary variable is only associated with the complex inequalities l I A I x u I , not with the bound constraints l B x u B . As we will now see, this choice has for consequence that the two sets of inequalities are not treated in the same manner. The user of Qpal can always declare the bound constraint l B

x u B as being part of the complex inequalities by enlarging the matrix A I with an identity matrix (the converse is not possible though), but this is not recommended since the solver is more efficient on bound constraints than on general linear inequality constraints.

The algorithm used by Qpal generates a sequence {λ k } k 1 of multipliers associated with the equality constraints of (QP') in the following manner [3,[START_REF] Delbos | Constrained optimization in seismic reflection tomography: a Gauss-Newton augmented Lagrangian approach[END_REF] (we describe a typical iteration).

• At iteration k, the algorithm solves

min l (x,y) u ℓ r k (x, y, λ k ), (1.5)
where ℓ r is the augmented Lagrangian (AL) defined for x ∈ R n , y ∈ R mI , and λ = (λ

I , λ E ) ∈ R mE × R mI by ℓ r (x, y, λ I∪E ) = g ⊤ x + 1 2 x ⊤ Hx + λ ⊤ I (A I x -y) + r 2 A I x -y 2 2 + λ ⊤ E (A E x -e) + r 2 A E x -e 2 2 .
(1.6)

The scalar r > 0 is called the augmentation parameter. The convex quadratic problem (1.5) is a priori easier to solve than (QP) since it has only bound constraints. Note that the bound constraints on x are treated explicitly in (1.5), while the complex inequalities l I A I x u I are relaxed into the AL. The minimization in (1.5) is done by a combination of various algorithms: gradient projection, active set, and conjugate gradient. Under the sole assumption that problem (1.1) has a solution, problem (1.5) has also a solution, say (x k+1 , y k+1 ). This solution is not necessary unique however.

• Then the multiplier λ k is updated by the formula

λ k+1 = λ k + r k A I x k+1 -y k+1 A E x k+1 -e . (1.7)
Even though (x k+1 , y k+1 ) is not uniquely determined as a solution to (1.5), the constraint value (factor of r k above) is independent of that solution, so that the multipliers λ k are unambiguously generated.

On entry in qpal solve, x and λ are supposed given: x is a primal variable that can be reasonably initialized by the user and λ is a multiplier that can be adequately initialized in a context like SQP. If one has no idea of the correct x and/or λ, one can just set x = P [l,u] 0 (the projection of 0 on the interval [l, u]) and/or λ = 0. The augmentation parameter r is set to some positive value (actually 1), unless a positive value has been given on entry in qpal solve through the argument options%r. When solving a sequence of similar QP's, like in the SQP algorithm, the value of r can be set to the one determined by qpal solve at the end of the previous run, which is available in the output variable info%r. The Qpal package is formed of the files and directories described below. In this description, PLAT is the name of an environment variable that designates the platform for which Qpal has to be compiled (see the description of the file make.a platform below and stage 2 of the installation procedure in section 1.4.2 for the details).

• The files COPYRIGHT.* give the conditions of use of Qpal. You are supposed to agree with these conditions to be authorized to use it. If this has not been done yet, send a filled in and signed copy of the commitment letter COPYRIGHT.pdf to the author of the solver.

• The file README gives a short description of the package and tells how to use it.

• The directory bin is originally empty and will contain the binaries gathered in libraries, named lib*.a, after compilation of Qpal:

libblas.$PLAT.a contains the blas routines used in Qpal, those in the directory blas (see below); this archive can be used if the compiler does not provide a blas library; liblapack.$PLAT.a contains the lapack routines used in Qpal, those in the directory lapack (see below); this archive can be used if the compiler does not provide a lapack library; libqpal.$PLAT.a contains the routines peculiar to Qpal, those in the directory src (see below).

• The directory blas contains the blas [1] routines used in Qpal and a Makefile to compile them. These compiled routines may be replaced by versions that have been tuned to a particular platform, for example by using the option -lblas of the compiler, if any.

• The directory cuter contains all the files that are useful to install Qpal in CUTEr [START_REF] Bongartz | CUTE: Constrained and unconstrained testing environment[END_REF][START_REF] Gould | CUTEr (and SifDec), a Constrained and Unconstrained Testing Environment, revisited[END_REF]; see how to proceed in the README file of the cuter directory.

• The directory doc contains this documentation in PDF (doc.pdf). It also contains the file VERSIONS.txt that briefly describes the successive versions of the solver.

• The directory example gives an elementary example of convex quadratic optimization problem, that is solved by Qpal. The goal of this example is to make concrete the way of encoding the problem and to call the QP solver. The example is implemented using both dense and sparse storage.

• The directory lapack contains the lapack routines used in Qpal and a Makefile to compile them. These compiled routines may be replaced by versions that have been tuned to a particular platform, for example by using the option -llapack of the compiler, if any.

• The directory libopt contains all the files that are useful to install Qpal in the Libopt environment [START_REF] Gilbert | LIBOPT -An environment for testing solvers on heterogeneous collections of problems[END_REF]; see how to proceed in the README file of the libopt directory.

• The file make.a platform is used by the Makefile's of the package to compile various programs. You will probably need to adapt this file to your own platform, by redefining some of its variables; more is said on this subject in stage 2 of section 1.4.2.

• The directory mod is originally empty and will contain, after compilation of Qpal, the module descriptor qpal mod.$PLAT.mod of the qpal mod module. It provides, in particular, the description of the public derived types discussed in section 2.1.

• The directory src contains the peculiar routines of the Qpal solver and a Makefile to compile them.

Installation

The installation of Qpal can be made by following the stages below.

1. Normally, the Qpal package is distributed as a tarball named

QPAL-xxx-distrib.tgz
where xxx stands for a version number. Place this tarball in a directory where you want to keep Qpal. Decompress and untar it using

tar -zxvf QPAL-xxx-distrib.tgz
This creates the directory QPAL-xxx-distrib, which can be renamed. Let us call it the qpal directory.

2. The second stage deals with platform matters. To know the chosen compiler, linker and their options, the Makefile's of the Qpal package include a file named make.$PLAT, located in the qpal directory. Here, PLAT is the name of an environment variable defined by a Unix/Linux command similar to setenv PLAT mach.os.comp

The string "mach.os.comp" above is arbitrary; the proposed form allows you to identify the machine type "mach", its operating system "os", and the chosen compiler suite "comp"; examples might be setenv PLAT pc.linux.pgf90 setenv PLAT mac.osx.gcc

The file "make.a platform" in the qpal directory is an example of such a make.$PLAT file and, now that the environment variable PLAT has been defined, you may want to adapt it to your platform by first copying it cp make.a platform make.$PLAT and next adapting the following variables, which are the only ones used in the Qpal Makefile's: As said above, this command places the Qpal libraries in the directory bin and module descriptors in the directory mod. These can now be linked to a program that uses Qpal as a convex quadratic solver.

Usage

We start in section 2.1 by specifying the data structures used by the Qpal solver, which is revealed by its various Fortran types. In section 2.2, we describe the arguments of the subroutines or function qpal allocate, qpal solve, qpal conv, and qpal deallocate, which are various tools making possible to run and control of the solver. In section 2.2.6, we give the typical sequence of statements that must precede a call to the solver.

RT n°0377

inria-00442295, version 1 -19 Dec 2009

Data structure

The Qpal solver is structured as a Fortran module, named qpal mod, with its private/public derived type definitions, data, functions, and subroutines. The entry point to the solver is qpal solve. The associated subroutine qpal solve gathers its arguments in structures (having a given Fortran derived type), in order to make clearer the links between them and to make easier passing them from one procedure to the other. The Qpal public derived types are described below.

A structure with the qpal data type public derived type has allocatable components that are allocated by the subroutine qpal allocate (section 2.2.1). This feature is a little peculiar since the structure must be declared in the user space, while it is allocated by the solver subroutine qpal allocate. The use of this F03 feature is motivated by the desire to discharge the user from taking care of the boring, humdrum, and error-prone allocation instructions. You will be able to use a qpal data type structure only after having called qpal allocate. Useless components are not allocated by qpal allocate. Deallocation is realized by qpal deallocate (section 2.2.5).

We remind the reader that the value of the component c of a variable v of a given derived type is specified by v%c.

Data types

Qpal can solve convex quadratic optimization problems having various data structures. The matrices A I and A E defining the constraints can indeed be dense or sparse, while the Hessian H of the objective can be dense, sparse, or can have an ℓ-BFGS structure. This section describes the Fortran public derived types that have been defined to store these various data structures.

We start with the public derived type for storing sparse matrices.

type, public :: qpal_sparse_type sequence integer :: nnz integer, allocatable :: i(:), j(:) double precision, allocatable :: v(:) end type qpal_sparse_type With the keyword sequence, the components are stored in the specified order (otherwise, no storage sequence is implied by the order of the component definitions). As a result, a structure with a type identical to qpal sparse type (also with the keyword sequence) but with a different type name will be correctly identified with a structure of type qpal sparse type.

Here is a description of the components.

• nnz: this integer specifies the number of nonzero elements of the matrix.

• i(:) and j(:): the integers i(k) and j(k), for k = 1, . . . , nnz, are the row and column indices of the kth nonzero element of the matrix. • v(:): the double precision value v(k), for k = 1, . . . , nnz, gives the value of the kth nonzero element of the matrix.

We now introduce the public derived type named qpal hessian type, which describes the memory representation of the Hessian H of the objective of problem (1.1). The Hessian can be dense, sparse, or have a direct or inverse ℓ-BFGS representation. missing informs Qpal that H = 0, so that the quadratic problem (1.1) is actually a linear optimization problem; dense informs Qpal that H is a dense matrix that must be stored in the component dense, which is a double precision array of dimension (n,n); sparse informs Qpal that H is a sparse matrix that must be stored in the component sparse, which is a structure of type qpal sparse type; lbfgs dir informs Qpal that H is a direct ℓ-BFGS matrix that is stored in the data structure lbfgs dir, as explained in the module modulopt lbfgs mod (the features of this module need not be known by the user of Qpal); lbfgs inv informs Qpal that H is an inverse ℓ-BFGS matrix; this data structure is essentially used for unconstrained problem, so that Qpal does not provide any code for dealing with that case, in particular no data structure.

• dense: this allocatable double precision array is used to store the Hessian when id = dense.

• sparse: this memory structure of type qpal sparse type is used to store a sparse Hessian when id = sparse. Since a Hessian is a symmetric matrix, only the lower triangular part (with row index i(k) the column index j(k)) must be filled in and will be used (values v(k) with i(k) > j(k) are ignored).

• lbfgs dir: this memory structure of type lbfgs dir hessian type is used to store a direct ℓ-BFGS matrix, when id = lbfgs dir. The memory structure is defined in the module modulopt lbfgs mod.

• lbfgs inv: this memory structure of type lbfgs inv hessian type is used to store an inverse ℓ-BFGS matrix, when id = lbfgs inv. The memory structure is defined in the module modulopt lbfgs mod.

The next derived type is used to describe the memory representation of the inequality and equality constraint Jacobians, A I and A E respectively. These Jacobians can be dense or sparse. values can be used, named by public parameters:

dense informs Qpal that the Jacobian is a dense matrix that must be stored in the component dense; sparse informs Qpal that the jacobian is a sparse matrix that must be stored in the component sparse.

• dense: this allocatable double precision array of appropriate dimension is used to store the constraint Jacobian when id = dense.

• sparse: this memory structure of type qpal sparse type is used to store a sparse constraint Jacobian when id = sparse.

Problem data

The public derived type qpal data type must be used to declare and store the data of problem (1.1). Its allocatable components are allocated by qpal allocate (section 2.2.1); only the useful variables for the considered problem are allocated.

type, public :: qpal_data_type integer :: n, nb, mi, me double precision, allocatable :: g(:) type(qpal_hessian_type) :: h double precision, allocatable :: lb(:), ub(:) type(qpal_constraint_type) :: ai double precision, allocatable :: li(:), ui(:) type(qpal_constraint_type) :: ae double precision, allocatable :: e(:) end type qpal_data_type

Here are the components of this type.

• n, nb, mi, me: these integers give respectively the number n of variables, the number of variables with a finite lower or upper bound, the number m I of inequality constraints, and the number m E of equality constraints.

• g(:): this allocatable double precision array of dimension n is used to store the vector g ∈ R n giving the linear part of the objective of problem (1.1).

• h: this data structure is used to store the Hessian H of the objective of problem (1.1), using the type qpal hessian type previously defined.

• lb(:), ub(:): these allocatable double precision arrays of dimension n are used to store the vectors l B and u B ∈ R n giving respectively the lower and upper bounds on the variables x.

• ai: this data structure is used to store the Jacobian A I of the inequality constraints of problem (1.1), using the type qpal constraint type previously defined.

• li(:), ui(:): these allocatable double precision arrays of dimension mi are used to store the vectors l I and u I ∈ R mI , giving respectively the lower and upper bounds on A I x.

• ae: this data structure is used to store the Jacobian A E of the equality constraints of problem (1.1), using the type qpal constraint type previously defined.

Of course, for any vector v, v ∞ v 2 , so that for identical tolerances, the solver stops more rapidly with the sup norm than with the Euclidean norm.

Valid values: '2', 'euc', 'inf'.

Default value: 'inf'.

• r: double precision variable giving the initial value of the AL parameter r. This value is increased by Qpal if the decrease in the constraint norm is not fast enough (compared to dcrt).

Valid value: > 0.

Default value: 1.

• verb: integer variable that specifies the verbosity level of the solver, i.e., the amount of information that is written on channel iout. The following values are meaningful:

= 0: nothing is written; = 1: in addition to the initial and final outputs, the solver writes one line at each AL iteration; = 2: one more line for each gradient projection (GP), descent projection (DP), and conjugate gradient (CG) phase; = 3: more details are given, including those from the GP-AS-CG inner algorithm; 4: additional possibly long lists of numbers are written.

Valid value: 0.

Default value: 0.

Solver diagnostics

The public derived type qpal info type must be used to define the variable getting information on the problem and its possible solution, on return from the solver. Its allocatable components are allocated by qpal allocate (section 2.2.1) and qpal solve; only the useful variables for the considered problem are allocated.

type, public :: qpal_info_type integer :: flag, iter, hvpd, avpd double precision :: f, glan, feas, r, estl double precision, allocatable :: glag(:), aix(:), aex(:) double precision, allocatable :: dx(:), dy(:), res_bie(:) end type qpal_info_type

The components of this type are gathered above by Fortran type and nature. For making their localization faster, they are presented below in alphabetic order. A variable that is said non available below has no memory allocated for it; hence do not use it if you want to avoid a bus error.

• aex: double precision array of dimension m E . It gives A E x f , the equality constraint value at the final point x f .

Not available if m E = 0.

• aix: double precision array of dimension m I . It gives A I x f , the inequality constraint value at the final point x f .

Not available if m I = 0.

• avpd: integer variable. It gives the number of matrix-vector products, with the constraint matrices A E and A I , realized during the run.

• dx: double precision array of dimension n. In case info%flag = 2 (unbounded problem), it gives a direction in x along which problem (1.1) is unbounded, when it is feasible.

Not available if info%flag = 2.

• dy: double precision array of dimension m I . In case info%flag = 2 (unbounded problem), it gives A I dx, where dx is given by info%dx.

Not available if info%flag = 2.

• est: double precision variable. It is known from [3] that there is a positive constant L such that the ℓ 2 norm of the constraints decreases at least by the factor min(L/r, 1) at each iteration. The value of estl provides a lower estimate of L.

• f: double precision variable. It gives the value f (x f ) of the quadratic objective at the final point x f .

• feas: double precision variable. It gives the value of the left hand side of (2.1) at the final point x f found by the solver, measuring the realized accuracy on the constraints. The norm • is the one specified by the option optimality norm on entry in Qpal.

• flag: integer variable that specifies the return status of Qpal. The following values are possible:

= 0: a solution is found up to the required accuracy; = 1: erroneous argument (possibly unallocated variables or trivially inconsistent constraints); = 2: the QP is likely to be either unbounded or infeasible; in the latter case, a shift of the constraint bounds to make it feasible would also make it unbounded; = 3: maximum iterations reached (iter or cgit counters); = 4: maximum matrix-vector products reached (hvpd or avpd counters); = 5: H is not positive semidefinite; = 6: failure in linesearch; = 7: diagonal preconditioning is required but the Hessian of the AL is not positive definite; = 8: impossible to minimize the AL function on the current active face and the AL parameter r cannot be decreased; = 9: something wrong in the linear algebra subroutines; = 99: something wrong, contact your guru.

• glag: double precision array of dimension n. It gives ∇ x ℓ(x f , λ f ), the gradient with respect to x of the Lagrangian ℓ at the final primal-dual pair (x f , λ f ).

• glan: double precision variable. It gives the norm of the gradient of the Lagrangian, ∇ x ℓ(x f , λ f ) , at the final primal-dual pair (x f , λ f ). The norm • is the one specified by the option optimality norm on entry in Qpal.

• hvpd: integer variable. It gives the number of Hessian-vector products realized during the run.

• iter: integer variable. It gives the number of AL iterations realized during the run.

• r: double precision variable. It is the final value of the AL parameter r.

• res bie: double precision vector of dimension n + m I + m E . It provides the final constraint residual :

res bie(1:n) = max(0, l B -x, x -u B )
res bie(n+1:n+mi) = max(0, l I -A I x, A I x -u I ) res bie(n+mi+1:n+mi+me) = A E x -e.

Running the solver

The Qpal solver makes available four subroutines (qpal allocate, qpal default options, qpal solve, qpal deallocate) and a function (qpal conv). The subroutine qpal allocate can be used to allocate variables of a data and an information structure (section 2.2.1). The subroutine qpal default options can be used to get the default options of the solver, before tuning these to a particular run (section 2.2.2). The subroutine qpal solve is used to solve a particular instance of problem (1.1) (section 2.2.3). The function qpal conv is provided to adapt the inside stopping criterion of the solver (section 2.2.4). Finally, the subroutine qpal deallocate can be used to free the memory allocated by qpal allocate and qpal solve (section 2.2.5).

In the description of the subroutines/functions, an argument flagged with (I) means that it is an input or intent(in) variable, which has to be initialized before calling the subroutine/function; an argument flagged with (O) means that it is an output or intent(out) variable, which has only a meaning on return from the subroutine/function; and an argument flagged with (IO) is an input-output or intent(inout) argument, which has to be initialized and has a meaning on return from the subroutine/function.

Data allocation with qpal allocate

It is recommended to call the subroutine qpal allocate before calling qpal solve. Its role is to allocate the allocatable components of the following two arguments: data, which is a structure containing the data of the quadratic problem to solve (see section 2.1.2), and info, which is a structure that will contain information on the run after having called the solver (see section 2.1.4). This allocation can be done without calling qpal allocate, but it is probably better to let a procedure do the job. After having called qpal allocate, you will be able to fill in the data structure with the problem data; note that the scalar quantities of this structure (dimensions, matrix type identificators, and numbers of nonzero elements if this is relevant) are set by qpal allocate, so that you only have to take care of the array components. On the other hand, the info structure will be fill in by the QP solver qpal solve, so that there is no reason to modify that structure. h scale (I): integer variable that determines the type of scaling that will be used in the ℓ-BFGS Hessian matrix (if h id has been set to lbfgs dir or lbfgs inv). This one can be set to the public values lbfgs scal scaling (matrix scaled by a scalar) or lbfgs diag scaling (matrix scaled by a diagonal matrix). h cold (I): logical variable that tells the ℓ-BFGS initialization procedure whether the ℓ-BFGS Hessian matrix will be used from scratch (set cold to .true.) or a warm restart will be done (set cold to .false.). This is only useful when h id has been set to lbfgs dir or lbfgs inv. ai id (I): integer variable. It identifies the type of memory space used to store the Jacobian A I of the inequality constraint of problem (1.1). See the description of the id component of the type qpal constraint type in section 2.1.1 for the possible values of this variable, which can be dense or sparse. ai nnz (I): integer variable. It gives the number of nonzero elements in the matrix A I in case this one is stored in a sparse structure (ai id has been set to sparse); it is meaningless otherwise. ae id (I): integer variable. It identifies the type of memory space used to store the Jacobian A E of the equality constraint of problem (1.1). See the description of the id component of the type qpal constraint type in section 2.1.1 for the possible values of this variable, which can be dense or sparse. ae nnz (I): integer variable. It gives the number of nonzero elements in the matrix A E in case this one is stored in a sparse structure (ae id has been set to sparse); it is meaningless otherwise. fout (I): integer variable. This is the channel number for the written outputs in qpal allocate. plevel (I): integer variable. It specifies the verbosity level of qpal allocate. The following values are meaningful: 0: silent mode; > 0: error and warning messages are printed on channel fout. prestr (I): character string of arbitrary length that will precede the lines printed by Qpal (if any). flag (O): integer variable. It provides information on the success of the allocations. Possible values are: = 0: allocations done; = 1: some allocation failed. data (O): variable of type qpal data type (see section 2.1.2). The usefull allocatable components of the variable are allocated by qpal allocate. info (O): variable of type qpal info type (see section 2.1.4). The usefull allocatable components of the variable are allocated by qpal allocate.

Setting default options with qpal default options

The subroutine qpal default options can be called to set the default options in a structure of qpal options type type. This call is normally made before setting the various component of the structure to values that are appropriate to the quadratic problem to solve. It allows the user not to have to set all the options, when the default ones are appropriate.

subroutine qpal default options (options) options (O): variable of type qpal options type (see section 2.1.3).

Solving the problem with qpal solve

The solver of the Qpal package is qpal solve. Here is the subroutine definition statement.

subroutine qpal solve (x, lm, data, info, options)

x (IO): double precision array of dimension n. It is the vector of variables x to optimize.

• On entry, it is an initial guess of the solution to (QP). It is taken as starting point by Qpal. This point need not be feasible.

• On return, when info%flag = 0, it is the optimal solution x f found by Qpal. lm (IO): double precision array of dimension n + m I + m E . It is the dual variable or KKT multiplier associated with the constraints of (QP). The first n components are associated with the bounds on x. The next m I components are associated with the inequality constraints l I A I x u I . The last m E components are associated with the equality constraints

A E x = e.
The multiplier λ B∪I associated with the bound constraints is actually the difference

λ B∪I := λ u B∪I -λ l B∪I
between the multiplier λ u B∪I associated with the upper bound and the multiplier λ l B∪I associated with the lower bound. Since l i < u i , either λ u i or λ l i vanishes, or both. Therefore, one can recover λ u B∪I and λ l B∪I by

λ u B∪I = λ + B∪I and λ l B∪I = λ - B∪I ,
where t + := max(t, 0) and t -:= max(-t, 0) for t ∈ R and () + and () -act componentwise for vectors.

• On entry, it is an initial guess of the dual solution to (QP). Qpal uses this vector to initialize the AL iterations (1.7).

• On return, when info%flag = 0, it is the vector of optimal dual variables λ f found by the optimizer. data (I): variable of type qpal data type (see section 2.1.2) that contains the data of problem (1.1). info (O): this is a variable of type qpal info type, whose components are described in section 2.1.4. It is used to get information on the problem and its possible solution, on return from the solver. options (I): this is a variable of type qpal options type, whose components are described in section 2.1.3. It is used to tune the behavior of the solver.

Convergence control with qpal conv

Qpal has its own convergence test, which uses the tolerances glan abs, feas abs, and mpos abs described in section 2.1.3. One may not be satisfied with this way of verifying convergence or one may want to add conditions. Qpal has anticipated this possible opinion by providing a mechanism that allows the user to modify the convergence criterion, which we now discuss. Each time Qpal observes that its convergence tests are satisfied, it calls the logical function qpal conv described below. If this one returns .true. (this is the case for the function in the standard distribution), the solver interrupts its iterative process and declares convergence; otherwise the iterations are pursued until another convergence test occurs. As a result, the user can modify the stopping criterion by writing a new function qpal conv (not the one in the standard distribution) as follows.

RT n°0377

inria-00442295, version 1 -19 Dec 2009

• If a more stringent criterion is desired, the additional conditions can simply be implemented in the new function qpal conv.

• If a completely different stopping criterion is desired, it suffices to implement this criterion in the new function qpal conv and to set the tolerances glan abs, feas abs, and mpos abs to a large number on entry into qpal solve. Then the internal stopping criterion will always be satisfied and the one in qpal conv will decide when to stop the iterations.

The compiled version of the function qpal conv in the standard distribution is in the archive bin/libqpal.$PLAT.a. There is no reason to recompile this archive with the new function qpal conv but to introduce its compiled version qpal conv.o in the link command, as in

Deallocation with qpal deallocate

The subroutine qpal allocate allocates memory for storing the data of problem (1.1) and another one to store the solver diagnostics. This memory allocation depends on the dimensions of the problem (n, m I , and m E ). The subroutine qpal deallocate can be used to deallocate these variables.

Here is the definition statement of the subroutine qpal deallocate. values are: = 0: deallocations done; = 1: some problem with a deallocation.

Calling sequence

We summarize below the structure of a program that uses qpal solve to solve a convex quadratic optimization problem. We assume for simplicity that all the instructions are in the same program unit. The primal-dual solution found by qpal solve is in (x, lm) and the information qpal solve gives on its run is in the argument info (section 2.1.4).

6. If the structures data and info allocated by qpal allocate are no longer useful, you may want to deallocate them by using qpal deallocate (section 2.2.5):

call qpal_deallocate (data, info, fout, plevel, flag)

Current limitations and perspectives

We list below the limitations of Qpal we are aware of, together with planned improvements. The list is certainly not exhaustive.

1. There is no projection stage along the last direction found by the conjugate gradient.

2. The augmented parameter r is only increased.

3. The code cannot detect constraint incompatibility. For the while this is reflected by an augmentation parameter r that blows up.

4. Positive semi-definite linear systems are currently solved by a preconditioned conjugate gradient algorithm. We plan to offer the possibility to solve linear systems by updated Cholesky factorizations.

Index

1 1 . 4 2 7 2 . 1

 114221 Scope of the program . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2 Detecting optimality . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.3 Brief description of the method . . . . . . . . . . . . . . . . . . . . . . . . . . 5 The package . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.4.1 Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.4.2 Installation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 Usage Data structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

F03 = # name

  of the Fortran 2003 compiler F03DBGFLAG = # option requiring generating code for debugging F03FLAGS = # usual options to use with $(F03) F03NOLD = # option preventing from making a load object 3. You are now ready to compile Qpal. Go into the directory src and type make

  type, public :: qpal_hessian_type integer :: id double precision, allocatable :: dense(:,:) type(qpal_sparse_type) :: sparse type(lbfgs_dir_hessian_type) :: lbfgs_dir type(lbfgs_inv_hessian_type) :: lbfgs_inv end type qpal_hessian_type Description of the components. • id: this integer component specifies the Hessian identificator. The following 5 values can be used, named by public parameters:

  type, public :: qpal_constraint_type integer :: id double precision, allocatable :: dense(:,:) type(qpal_sparse_type) :: sparse end type qpal_constraint_type RT n°0377 inria-00442295, version 1 -19 Dec 2009 Description of the components. • id: this integer component specifies the constraint Jacobian identificator. The following 2

$

  (F03) -o main ... main.o qpal conv.o .../bin/libqpal.$PLAT.a ... Since the new version of qpal conv.o is in the command line the one in the archive bin/libqpal.$PLAT.a will not be selected at link time. Another possibility is to remove the default version of qpal conv.o from the archive bin/libqpal.$PLAT.a, using ar -d bin/libqpal.$PLAT.a qpal conv.o and then to provide a new version of qpal conv.o at link time. Here is the definition statement of the logical function qpal conv. logical function qpal conv (data, x, lm, aix, aex, glag, qp inf, qp out) qpal conv (O): logical output value. Set it to .true. is the primal-dual variables (x, lm) are appropriate; to .false. otherwise. data (I): variable of type qpal data type (see section 2.1.2) that contains the data of problem (1.1). x (I): double precision array of dimension n. It is the value of the vector of variables x to optimize at the time the function is called by Qpal. lm (I): double precision array of dimension n + m I + m E . It is the value of the dual variable or KKT multiplier associated with the constraints of problem (1.1) at the time the function is called by Qpal. The first n components are associated with the bounds on x. The next m I components are associated with the inequality constraints l I A I x u I . The last m E components are associated with the equality constraints A E x = e. The multiplier λ B∪I associated with the bound constraints can be viewed as the difference λ B∪I := λ u B∪I -λ l B∪I between the multiplier λ u B∪I associated with the upper bound and the multiplier λ l B∪I associated with the lower bound. aix (I): double precision array of dimension m I . It gives the value of the product A I x at the time the function is called by Qpal. aex (I): double precision array of dimension m E . It gives the value of the product A E x at the time the function is called by Qpal. glag (I): double precision array of dimension n. It gives the value of the gradient of the Lagrangian at the time the function is called by Qpal. qp inf (I): double precision variable used to detect an infinite bound in the QP. qp iout (I): integer variable specifying the output channel used by QPAL for printing.

  . . . . . 8 2.1.1 Data types . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.1.2 Problem data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.1.3 Solver options . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.1.4 Solver diagnostics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.2 Running the solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.2.1 Data allocation with qpal allocate . . . . . . . . . . . . . . . . . . . 15 2.2.2 Setting default options with qpal default options . . . . . . . . . . 16 2.2.3 Solving the problem with qpal solve . . . . . . . . . . . . . . . . . . 17 2.2.4 Convergence control with qpal conv . . . . . . . . . . . . . . . . . . . 17 2.2.5 Deallocation with qpal deallocate . . . . . . . . . . . . . . . . . . . 19 2.2.6 Calling sequence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

	3 Current limitations and perspectives	20
	References	21
	Index	21

  The other dimensions must be set with care, since they directly intervene in the allocation of variables in the structure data and are used as the dimensions of the problem in the solver qpal solve (section 2.2.3). h id (I): integer variable. It identifies the type of memory space used to store the Hessian H of the objective of problem (1.1). See the description of the id component of the type qpal hessian type in section 2.1.1 for the possible values of this variable, which can be missing, dense, sparse, lbfgs dir, or lbfgs inv. h nnz (I): integer variable. It gives the number of nonzero elements in the lower triangular part of the Hessian H in case this one is stored in a sparse structure (h id has been set to sparse); it is meaningless otherwise.

	subroutine qpal allocate (n, nb, mi, me, &
	h id, h nnz, h mys, h scale, h cold, &
	ai id, ai nnz, ae id, ae nnz, &
	fout, plevel, prest, flag, data, info)

n, nb, mi, me (I): positive integer variables. They define the dimensions n = n, mi = m I , and me = m E of problem (1.1), as well as the number nb of bounds on the primal variables x. The latter must only be vaguely defined: if there is no bound on x, set nb = 0, otherwise give nb an arbitrary positive (> 0) value. Indeed, qpal allocate only trusts the sign of nb to decide whether memory must be allocated for the bounds on x, i.e., for the variables data%lb(1:n) and data%ub(1:n). h mys (I): positive integer variable. It is used, in case an ℓ-BFGS Hessian is declared with h id, to specify the number of pairs (y k , s k ) that are used to form the Hessian approximation; it is meaningless otherwise.

RT n°0377

inria-00442295, version 1 -19 Dec 2009

  subroutine qpal deallocate (data, info, fout, plevel, flag) data (IO): variable of type qpal data type (see section 2.1.2) that contains the data of problem (1.1) and was previously used as argument of qpal allocate and qpal solve. Its allocated components are deallocated by qpal deallocate. info (IO): variable of type qpal info type (see section 2.1.4) that contains information on problem (1.1) and was previously used as argument of qpal solve. Its allocated components are deallocated by qpal deallocate. fout (I): integer variable. This is the channel number for the written outputs in qpal deallocate. plevel (I): integer variable. It specifies the verbosity level of qpal deallocate. The following values are meaningful: = 0: silent mode; > 0: error and warning messages are printed on channel fout. flag (O): integer variable. It provides information on the success of the deallocations. Possible

  If the name of some of the public types, parameters, or procedures of qpal mod are in conflict with your own variables, you can rename them. For example, to rename the Qpal public parameter dense (resp. sparse) into qpal dense (resp. qpal sparse), write insteadThe flag argument should be zero on return from qpal allocate. See section 2.2.1 for more information.3. After execution of qpal allocate, the allocatable components of the variable data have been allocated and you can fill in them: Before calling the solver qpal solve, you still have to give an initial value to the primaldual variables (of course these must have been declared in the calling subroutine/function as double precision vectors with the respective dimensions n and n+mi+me):x = ... lm = ... You also have to set in a structure here called options, whose type is qpal options type (section 2.1.3), the options of the solver that should not have the default value set by qpal default options (section 2.2.2):Of particular importance are of course the variables options%glan abs, options%feas abs, and options%mpos abs specifying the stopping criterion.

		data%g = ...
		data%lb = ...
		data%ub = ...
		...
	inria-00442295, version 1 -19 Dec 2009	4. call qpal_default_options (options) options%avpd = ... options%cgit = ... ...
		use qpal_mod, qpal_dense => dense, &
		qpal_sparse => sparse
		2. Declare the variables n, nb, mi, me, h id, h nnz, h mys, h scale, h cold, ai id, ai nnz, ae id,
		ae nnz, fout, plevel, prest, flag, data, and info, set those that are relevant to the current
		run (except data and info). Then allocate the relevant allocated components of the data
		and info structures by
		call qpal_allocate (n, nb, mi, me, &
		h_id, h_nnz, h_mys, h_scale, h_cold, &
		ai_id, ai_nnz, ae_id, ae_nnz, &
		fout, plevel, prest, flag, data, info)
		RT n°0377

1. Specification of the use of the Qpal module qpal mod: use qpal_mod 5. You can now call qpal solve (section 2.2.3): call qpal_solve (x, lm, data, info, options)

  16 qpal allocate, 15-16, 19 qpal deallocate, 19, 20 qpal solve, 17, 20 qpal, see directory qpal allocate, see public subroutine qpal constraint type, see public derived type qpal conv, see function qpal data type, see public derived type qpal deallocate, see public subroutine qpal default options, see public subroutine qpal hessian type, see public derived type qpal info type, see public derived type qpal mod, see module qpal options type, see public derived type qpal solve, see public subroutine qpal sparse type, see public derived type sequence, 8 sparse, see public parameter src, see directory

		Hessian, 3	qpal constraint type, 9-10
		identificator, 9	qpal data type, 10-11
		Lagrangian, 4 lapack, see directory lbfgs dir, see public parameter lbfgs inv, see public parameter Libopt, 6 libopt, see directory	qpal hessian type, 8-9 qpal info type, 13-14 qpal options type, 11-13 qpal sparse type, 8 public parameter dense, 9, 10 lbfgs dir, 9
		missing, see public parameter mod, see directory module qpal mod, 8	lbfgs inv, 9 missing, 9 sparse, 9, 10 public subroutine qpal default options,
		objective, 3
	inria-00442295, version 1 -19 Dec 2009	active set, 4 augmentation parameter, 5 option avpd, 11 cgit, 11 dcmin, 11 dcrt, 11 dxmin, 11 feas abs, 11 glan abs, 12 hvpd, 12 inf, 12 infs, 12 iout, 12 iter, 12 mpos abs, 12 optimality norm, 12 r, 13 verb, 13 precs, see subroutine	cuter, 6 doc, 6
		augmented Lagrangian, 5 public derived type	example, 6
		bin, see directory blas, see directory	lapack, 6 libopt, 6 mod, 6
		constraint active, 4 inactive, 4	qpal, 7 src, 6 doc, see directory
		constraint Jacobian identificator, 10	example, see directory
		COPYRIGHT, 6	F03, see Fortran 2003
		CUTEr, 6	feasible set, 4
		cuter, see directory	Fortran 2003, 4
		dense, see public parameter directory	function qpal conv, 17-18
		bin, 6	gfortran, 4
		blas, 6	gradient, 3
		RT n°0377

http://www.inria.fr ISSN 0249-0803

The parallelisation of the loops could be done with OpenMP directives [8].

• e(:): this allocatable double precision array of dimension me is used to store the vector e ∈ R mE , giving the RHS of the equality constraints.

Solver options

The public derived type qpal options type is used to describe the options of the Qpal solver, i.e., those parameters that can be used to tune the behavior of the solver.

type, public :: qpal_options_type integer :: iout, verb, iter, cgit, hvpd, avpd character(len=3) :: optimality_norm double precision :: glan_abs, feas_abs, mpos_abs double precision :: inf, dxmin, dcmin double precision :: r, dcrt, infs end type qpal_options_type

The components of this type are gathered above by Fortran type and nature. For making their localization faster, they are presented below in alphabetic order. Valid and default values are indicated.

• avpd: integer variable specifying the maximal number of constraint-vector products (A I v and A E v) allowed.

Valid value: > 0.

Default value: 1000.

• cgit: integer variable specifying the maximal number of conjugate gradient (CG) iterations allowed.

Valid value: > 0.

Default value: 1000.

• dcmin: positive double precision variable that is used to detect active bounds with index i ∈ I.

To this respect, there must also hold

Valid values: > 0.

Default value: 10 -10 .

• dcrt: double precision variable giving the desired decrease factor for the constraint norm. A small value looks better but it forces Qpal to take a large AL parameter r, inducing illconditioning. Qpal will not try to decrease the constraint norm by a factor smaller than dcrt.

For a small easy problem, dcrt can be chosen rather small. For a large scale ill-conditioned problem, choose dcrt close to 1.

Valid values: ]0, 1[. Default value: 10 -1 .

• dxmin: positive double precision variable that is used to detect active bounds with index i ∈ B.

To this respect, there must also hold

Valid values: > 0.

Default value: 10 -10 .

• feas abs: double precision variable specifying the required accuracy on the constraint norm (see the component optimality norm for the definition of the norm • ). Knowing that the bound constraints on x are satisfied along the iterations, hence at the final point, only the Valid value: 0. Default value: 10 -8 .

• glan abs: double precision variable specifying the required absolute accuracy on the gradient of the Lagrangian norm (see the component optimality norm for the definition of the norm • ).

If the solver stops at (x, λ) with info%flag = 0, there holds

Valid value: 0.

Default value: 10 -8 .

• hvpd: integer variable specifying the maximal number of Hessian-vector products (Hv) allowed.

Valid value: > 0.

Default value: 1000.

• inf: double precision variable. This is the value used by Qpal to detect inexistent bounds on the variables or in the inequality constraints. In other words, if after having defined inf, you set ub(i) [resp. lb(i)] (these are components of a variable of type qpal data type) to a value inf [resp.

-inf], Qpal will consider that the ith variable is not subject to an upper bound [resp. to a lower bound]. The same value and rule are used to locate bounds on A I x. As a result, it is recommended to use a large value.

Valid value: > 0.

Default value: huge(1.d0).

• infs: double precision variable specifying what is an infinite stepsize. If such a stepsize along the current search direction leads to a decrease in the objective f without encountering a bound, the direction is considered to be a direction of unboundedness and the solver stops with info%flag = 2.

Valid value: > 0.

Default value: 10 20 .

• iout: integer variable that is taken as the channel number for the written outputs. These are indeed written by write (iout,...) ...

Valid value: > 0.

Default value: 6.

• iter: integer variable specifying the maximal number of AL iterations allowed.

Valid value: > 0.

Default value: 100.

• mpos abs: double precision variable giving the required accuracy on the multiplier positivity.

If the solver stops with info%flag = 0, there holds λ l -mpos abs and λ u -mpos abs.

Valid value: 0. Default value: 10 -8 .

• optimality norm: this string, made of at most 3 characters, specifies the type of vector norm • that must be used to check optimality:

-'2' or 'euc' for the ℓ 2 or Euclidean norm v 2 := ( i v 2 i ) 1/2 , -'inf' for the infinity or sup norm v ∞ := max i |v i |.


C Voyant

M Muselli
email: marc.muselli@univ-corse.fr

C Paoli

M L Nivet

P Poggi

P Haurant

PREDICTABILITY OF PV POWER GRID PERFORMANCE ON INSULAR SITES WITHOUT WEATHER STATIONS: USE OF ARTIFICIAL NEURAL NETWORKS**

ABSTRACT: The official meteorological network is poor on the island of Corsica: only three sites being about 50 km apart are equipped with pyranometers which enable measurements by hourly and daily step. These sites are Ajaccio (41°55"N and 8°48"E, seaside), Bastia (42°33"N, 9°29"E, seaside) and Corte (42°30"N, 9°15"E average altitude of 486 meters). This lack of weather station makes difficult the predictability of PV power grid performance. This work intends to study a methodology which can predict global solar irradiation using data available from another location for daily and hourly horizon. In order to achieve this prediction, we have used Artificial Neural Network which is a popular artificial intelligence technique in the forecasting domain. A simulator has been obtained using data available for the station of Ajaccio that is the only station for which we have a lot of data: 16 years from 1972 to 1987. Then we have tested the efficiency of this simulator in two places with different geographical features: Corte, a mountainous region and Bastia, a coastal region. On daily horizon, the relocation has implied fewer errors than a "naïve" prediction method based on the persistence (RMSE=1468 Vs 1383Wh/m² to Bastia and 1325 Vs 1213Wh/m² to Corte). On hourly case, the results were still satisfactory, and widely better than persistence (RMSE=138.8 Vs 109.3 Wh/m² to Bastia and 135.1 Vs 114.7 Wh/m² to Corte). The last experiment was to evaluate the accuracy of our simulator on a PV power grid localized at 10 km from the station of Ajaccio. We got errors very suitable (nRMSE=27.9%, RMSE=99.0 W.h) compared to those obtained with the persistence (nRMSE=42.2%, RMSE=149.7 W.h).

PRESENTATION AND ISSUE

We present the results of the prediction of global radiation using Artificial Neural Networks (ANN) which are a popular artificial intelligence technique in the forecasting domain [START_REF] Bosch | Daily solar irradiation estimation over a mountainous area using artificial neural networks[END_REF]. Inspired by biological neural networks, researchers in a number of scientific disciplines are designing ANNs to solve a variety of problems in decision making, optimization, control and obviously prediction [START_REF] Jain | Artificial neural networks: A tutorial[END_REF][START_REF] Kaligirou | Artificial neural networks in RE systems applications: a review[END_REF]. In this context, our aim was to answer to the following question: Can we design an ANN of a site for which there is a lot of solar radiation data available and use this ANN to predict a PV power grid performance of another site? We tried to answer to this question with sites located on the island of Corsica (France). The island is characterized by a Mediterranean climate and a hilly terrain. The official meteorological network (from the French Meteorological Organization) is very poor: only three sites being about 50 km apart are equipped with pyranometers and enable measurements by hourly and daily step. These sites are Ajaccio (41°55"N and 8°48"E, seaside), Bastia (42°33"N, 9°29"E, seaside) and Corte (42°30"N, 9°15"E average altitude of 486 meters). In this study, we focus on the prediction of global solar irradiation on a horizontal plane for daily and hourly horizon. These time steps have been chosen according to the electricity supplier (EDF: Electricité De France) who is interested in the estimation of the fossil fuel saving. It is very important for a remote site where electrification can be problematic [START_REF] Muselli | Improved procedure for stand-alone photovoltaic systems sizing meteosat satellite images[END_REF], and for quantifying the solar potential available. Indeed, this is very important both for the power plant implementation and for sizing of PV array. Solar radiation has been measured for a long time, but even today there are many unknown characteristics of its behavior. So, it seems appropriate to develop a prediction methodology using the data available in another location in order to overcome the lack of weather station and the demand for renewable energy source on the island.

PHYSICAL PHENOMENA

There are two approaches that allow quantifying solar radiation: the "physical modeling" based on physical process occurring in the atmosphere and influencing solar radiation [START_REF] Badescu | Modelling Solar radiation at the earth surface[END_REF], and the "statistic solar climatology" mainly based on time series analysis [START_REF] Badescu | Modelling Solar radiation at the earth surface[END_REF]. We have chosen to combine these two methods to improve the quality of prediction. In this work, we have used the physical phenomena in an attempt to overcome the seasonality of the resource. When studying the solar energy on the earth"s surface with time series, habitually the nonstationarity perturbs the quality of prediction. Often it"s necessary to apprehend the periodic phenomenon [START_REF] Zhang | Neural network forecasting for seasonal and trend time series[END_REF][START_REF] Bourbonnais | Analyse des séries temporelles[END_REF]. In daily case, there is seasonality with annual period, and in the hourly case there is in addition daily phenomenon.

We can use the extraterrestrial global horizontal irradiance as stationarization setting. The deterministic component of the series is thus reduced, leaving more important place to the stochastic part (cloud cover). We chose a multiplicative pattern and we stationarize the time series by using the extra-terrestrial irradiation. We divide the time series by the coefficient of extraterrestrial insolation. In the hourly case, the Earth's rotation adds a new periodic component. The solar altitude angle can be easily linked to the number of photon interacting on a horizontal surface at ground level (day light). In a first approximation, it can be considered as directly proportional to the energy received. It is why we have chosen to remove the periodicity (annual and daily) dividing the time series by the coefficient of extraterrestrial insolation, but also by the sinus of the solar height. These two treatments (hourly and daily) have been used designing an ANN in order to predict the global radiation on horizontal plane.

SOLAR RADIATION PREDICTION FOR A SITE WITHOUT EXPERIMENTAL DATA MEASURED

ANNs are capable of capturing the characteristics of any phenomenon with a good degree of accuracy. This technology often offers an alternative to traditional physical-based models, and excels at uncovering relationships in data. It is also a powerful non-linear estimator which is recommended when the functional form between input-output is unknown or it is not well understood. The main difficulties using the ANN technology are to have enough data (number of learning elements), to be sure of the data quality and to find the best architecture of ANN. Following the literature [START_REF] Hocaoglu | Hourly solar radiation forecasting using optimal coefficient 2-D linear filters and feed-forward neural network[END_REF], we designed and optimized a Multi-Layer Perceptron (MLP) network which is the most used of the ANN architectures. Solar radiation data from 1972 to 1987 have been used for training and testing the MLP. 80% of the data available were used to train the MLP and 20% to test the MLP. We used the Matlab software to establish our network which has 8 entries, 3 hidden neurons and 1 output. These 8 entries are the 8 previous hourly measures considered (t, t-1, t-2, .. , t-7) of global solar irradiation (Wh/m²) and the output is the measure to predict for the next hour (t+1). In daily case, hours are replaced by days. We obtained good results with a prediction and a learning done on the same site: nRMSE less than 20.3% in daily horizon, 19.5% for hourly horizon and 16 years of learning. We want to use this ANN, trained successfully on site with lots of data available (Ajaccio), to estimate the insolation on a site which has no meteorological station and therefore no measure history. Many studies have tested this kind of prediction [START_REF] Bosch | Daily solar irradiation estimation over a mountainous area using artificial neural networks[END_REF][START_REF] Ks Reddy | Solar resource estimation using artificial neural networks and comparison with other correlation models[END_REF], often by establishing a regression of parameters like latitude, longitude, altitude, days of the year, etc. Their methodology doesn"t use the time series property and therefore doesn"t allow predicting the stochastic cloud cover but only an insolation average. As we developed above, in addition to the data of Ajaccio, we have two meteorological stations: one in Corte (mountain location) and the other one in Bastia (coastal location). In order to give a first answer to the question of the training relocation feasibility, we proposed to test both in Bastia and Corte, three models of prediction. In the case A, we use the MLP trained with data from Ajaccio (16 years from 1972 to 1987). In the case B, we use the MLP trained with data from the selected place. We have 5 years of data in Corte (from 2002 to 2006) and also 5 years in Bastia (from 1991 to 1995). And in the case C, we use a "naïve" prediction method based on the persistence. The following tables present our results for a full year (1996 for Bastia and 2007 for Corte) of global horizontal irradiance prediction by daily (Table 1) and hourly step (Table 2). . Even it"s difficult to compare the Corte and Bastia prediction (the years of forecasting are not the same), we see that the relocation is partly compensated by the large period of learning available at Ajaccio (16 years). However, the relocation of learning implies fewer errors than the "naïve" predictor (RMSE=1468 Vs 1383Wh/m² to Bastia and 1325 Vs 1213Wh/m² to Corte). Also, the results of correlation coefficient do not show significant differences between cases A and B. Whatever the learning location, the Corte global radiation is more predictable than Bastia ones, which combines a particular hydrography and orography.

On Table 2, we can see the difficulty of hourly predicting on the site of Corte (nRMSE =114.7Wh/m² for relocation learning and 109,6Wh/m² for standard learning). A possible explanation could be that in mountainous region the climate change can occur very quickly, confusing an ANN learned in a coastal region. The results are still satisfactory, and widely better than those of persistence. If we consider only the model of prediction A, the two tables show that the worst results are obtained in Corte (RMSE>114Wh/m²) for the hourly case, and in Bastia (RMSE>1300Wh/m²) for the daily case. The figure 1 illustrates these two cases. On part (a) and (b), forecasting is done during the year 1996, and on part (c) and (d) on the year 2007. We can see that the ANN cannot predict atypical insolation (picks of the error dashed line) as the day with a lots of cloud cover. These phenomena are compounded by the relocation, but they also exist when learning and prediction are made on the same site. On the Figure 1.b, we can see that the relocation proposed with the Bastia"s station suffers of an adjustment problem. Indeed, it overestimates the weak measures, and we underestimate the highest. This phenomenon is probably due to the normalization of data during the simulation with ANN. This normalization does not include any previous data from Bastia (assumed unknown), and is made in relation to the maximum and minimum insolation of Ajaccio. This extremum can be different on the two locations. The hourly horizon shows a great similarity between the measurement and prediction (Figure 1.c-d), this kind of prediction is relatively well supported by delocalized learning.

APPLICATION OF THE RELOCATION METHOD FOR A PV SYSTEM

A frontage PV system has been installed recently in our laboratory (Ajaccio). It has a nominal power of 6.525 kW composed by respectively 1.8 kW and 4.725 kW amorphous and mono-crystal PV modules built in 6 independent power subsystems. PV power predictions from ANN methodology described in this paper have been computed from one of this whole PV plant on a frontage side exposed to the south (azimuth null) and tilted at 80°. The PV system is composed by 9 SUNTECH 175S-24Ac for a 1.175 kW nominal power connected to a 1.85 kW SUNNY BOY SMA inverter for PV production on the grid. For the PV power calculation, we use in first approximation, a linear production based on a constant PV plant efficiency ηPV = 13%, measured from a PV-KLA I-V curve plotter device (INGENIEURBÜRO Mencke & Tegtmeyer): E PV (Wh) = η PV I  S. I   is the hourly global irradiation on the PV system (β = 80°) estimated from ANN, and S is the usable surface of the PV system under consideration (S = 10.125 m²).To predict this energy, we have used the 16 years of global horizontal irradiation available on the site of Ajaccio, like learning set. Classical models are used to compute tilted irradiation for an 80° angle: I  = I  [ I  / I  ] clear sky . The last term of this equation was calculated from the software PVSYST [START_REF] Pvsyst | PC software package for the study, sizing, simulation and data analysis of complete PV systems[END_REF]. The ANN designing previously has been used this time to predict the energy produced by the PV plant of our laboratory. The distance is about 10 km from the place of obtaining the series of training. The results presented in Figure 2 seem very promising, but can not completely validate the methodology, because of the limited number of comparable events. It is interesting to note that in the case with a simulated ANN and learning delocalized, we get errors very suitable (nRMSE=27.9%, RMSE=99 Wh) compared to those obtained with the persistence (nRMSE=42.2%, RMSE=149.7 Wh).

CONCLUSION

In this work we have studied a methodology which can predict global solar irradiation using data available from another location for daily and hourly horizon. In order to achieve this prediction, we have used an ANN. In the case of global horizontal radiation, the outsourced learning solution is possible in the hourly case: errors are almost identical to the case of a localized learning with a five years history. Our forecast methodology allows studying many sites with no meteorological station. The results show that if the hourly horizon is well suited to this relocation, the daily horizon has yet to be studied, to consider a generalization of the and the versatility of ANN. The prediction on Corte (mountain town), where environmental constraints are not at all that we can have in Ajaccio, is very reliable. For Bastia (a seaside town like Ajaccio), the results are very promising but requires further investigation (especially in the daily case). The results found in this study are very interesting because they show "hidden links" between the three sites. The results of the prediction of PV systems, shown that ANN and delocalized learning may be very useful, but need to increase and finalize the comparison between measurement and simulation. In next works, we try to improve our model, and validate our approach with solar module installed on our laboratory (Vignola, Ajaccio). At term, this tool (ANN) could eventually help the system manager, for implanting new PV central and so not limit the abundance to the sites with weather stations.

REFERENCE

Figure 1 :

 1 Figure 1: Example of simulated and experimental global solar irradiation in Bastia for the daily horizon (a) and in Corte for the hourly horizon (c). Dashed line is the error, solid line is the measured insolation and marks are the prediction. On the bottom corresponds to the correlation between experimental and simulated global irradiation in Bastia for the daily horizon (b) and in Corte for the hourly horizon (d).

Figure 2 :

 2 Figure 2: Comparison between computed and measured energy generated with the frontage PV system used in experimentation.

Table I :

 I Comparison of the 3 techniques of prediction for the daily forecast: horizon 1 day

	Forecast Location	Techniques of Prediction	RMSE (Wh/m²) NRMSE (%)± IC95	CC
	Bastia	Ajaccio (A)	1383 29.19±0,13	0.842
	(year	Bastia (B)	1288	27.51±0,2 0.842
	1996)	Persistence (C)	1468 31.4		0.807
		Ajaccio (A)	1213 25.88±0,17	0.887
	Corte (year 2007)	Corte (B)	1112 23.73±0,15	0.887
		Persistence (C)	1325 28.3		0.844
	On the table 1, as we could imagine, the training
	location has modified the quality of the prediction (error
	represented by Root Mean Square Error with Interval
	Confidence 95%, Normalised Root Mean Square Error
	and Correlation Coefficient)

Table II :

 II Comparison of the 3 techniques of prediction for the hourly forecast: horizon 1 hour

	Forecast Location	Techniques of Prediction	RMSE (Wh/m²) NRMSE (%)± IC95	CC
	Bastia (year 1996)	Ajaccio (A) Bastia (B) Persistence (C)	109.3 22.36±0.32 109.1 22.32±0.32 138.8 27.22	0.916 0.916 0.869
		Ajaccio (A)	114.7 23.12±0.6 0.916
	Corte (year 2007)	Corte (B) Persistence (C)	109.6 22.08±0.32 135.1 26.50	0.919 0.889

a part of this research is founded by CTC (Collectivité Territoriale de Corse)