Amine Boumaza
email: amine.boumaza@loria.fr

Bruno Scherrer
email: bruno.scherrer@loria.fr

Convergence and Rate of Convergence of a Foraging Ant Model

We present an ant model that solves a discrete foraging problem. We describe simulations and provide a complete convergence analysis: we show that the ant population computes the solution of some optimal control problem and converges in some well defined sense. We discuss the rate of convergence with respect to the number of ants: we give experimental and theoretical arguments that suggest that this convergence rate can be superlinear with respect to the number of agents. Furthermore, we explain how this model can be extended in order to solve optimal control problems in general and argue that such an approach can be applied to any problem that involves the computation of the fixed point of a contraction mapping. This allows to design a large class of formally well understood ant like algorithms for problem solving.

I. INTRODUCTION

Swarm intelligence [START_REF] Beni | Swarm intelligence[END_REF], [START_REF] Bonabeau | Swarm intelligence: from natural to artificial systems[END_REF] was introduced in the 1990s as a novel nature-inspired approach for problem solving. The inspiring source is the behavior of real insects: a population of simple interacting agents, communicating indirectly through an environment, constitutes a massively distributed algorithm for solving a given task (e.g. foraging, flocking, labor division, prey capture, ...). Referring to neural networks, another nature-inspired approach for massively distributed computing, Kohonen wrote in 1988 [START_REF] Kohonen | Self-Organization and Associative Memory[END_REF]:

"One of the fundamental tasks of this new "neural networks" science is to demonstrate by mathematical analysis, computer simulations, and even working artificial sensory and control systems that it is possible to implement massively parallel information-processing functions using components the principles of which are not mysterious but already familiar from computer technology, communication science, and control engineering. There is nothing in the "neural network" area which were not known, in principle at least, from constructs already in use or earlier suggested." The motivation of this paper is to argue for a similar statement about swarm intelligence.

Classical engineering problem solving and swarm intelligence can be viewed as alternative approaches of the same problem. Then what makes these approaches seem so different ? As Sutton pointed it [START_REF] Sutton | Artificial intelligence as a control problem: Comments on the relationship between machine learning and intelligent control[END_REF] (when discussing the relations between "modern" Machine Learning and "classical" Intelligent Control), we could "characterize the split as having to do with the familiar dilemma of choosing between obtaining clear, rigorous results on the one hand, and exploring the most interesting, powerful systems on the other."

Roughly speaking, research on swarm intelligence is often focused on impressive proof-of-concept applications through extensive experimental simulations while classical engineering problem solving relies on theoretical convergence proofs for "toy" problems. The intent here is not to judge these approaches, which are clearly complementary, but rather to express our belief that filling the gap between both is a great research challenge.

There have not been many attempts in the literature that try to actually fill this gap, in particular in the area of ant-like algorithms. One interesting exception is the probably best known instance of swarm intelligence: the Ant-colony optimization (ACO) meta-heuristic for combinatorial problems, which has been extensively studied. Most theoretical works on ACO have recently been reviewed in [START_REF] Dorigo | Ant colony optimization theory: a survey[END_REF]. In this review, the authors present theorems for "convergence in value" (ACO is guaranteed to find an optimal solution in finite time with a probability that can be made arbitrarily close to 1) and "convergence in solution" (provided a properly designed decreasing exploration parameter, the ACO asymptotically converges to an optimal solution almost surely). They also relate ACO to "more standard" optimization techniques such as stochastic gradient descent and cross entropy. The interested reader should go through this article for further details.

In this paper, we come back to the original inspiration of ant algorithms, where a population of simple agents (that may be viewed as mimicking the behavior of real ants) efficiently solve a foraging problem. Many models and algorithms have been proposed to solve this problem, see for example [START_REF] Deneubourg | The blind leading the blind: Modeling chemically mediated army ant raid patterns[END_REF], [START_REF] Deneubourg | The selforganizing exploratory pattern of the argentine ant[END_REF], [START_REF] Resnik | Turtles termites and traffic jams[END_REF], [START_REF] Panait | Ant foraging revisited[END_REF], [START_REF] Panait | A pheromone-based utility model for collaborative foraging[END_REF], [START_REF] Simonin | Construction of numerical potential fields with reactive agents[END_REF]. Such works mainly propose realistic models of the behavior of real ant colonies, from which they derive algorithms and methods for solving interesting problems. However, they do not go beyond simulations nor provide formal analysis of the algorithms.

In this article, we introduce a massively distributed ant model that is guaranteed to asymptotically solve the foraging problem in a certain sense. We provide a convergence proof of the model, and a rate of convergence analysis with respect to the number of agents. Our analysis relies on optimal control theory, parallel distributed computing and graph theory. To our knowledge, this is the first attempt to present an ant model for the foraging problem and formal analysis of the algorithm in terms of convergence and rate of convergence is provided.

The remaining of the paper is organized as follows. Section II provides a description of our ant model and discusses some simulations that were made to measure the convergence and the rate of convergence. Section III gives a formal proof that there is convergence in some sense and the theoretical analysis is used to discuss the influence of the model parameters. Section IV analytically discusses the rate of convergence: it gives formal arguments that explain why one observes a superlinear rate of convergence. Finally section V provides a discussion about the scope of our model and its relations to other models proposed in the literature.

II. A SIMPLE ANT MODEL

In this section, we describe a simple ant model, which is aimed at solving a foraging task. We provide simulations that illustrate the behavior of the model and experimental measurements of its performance. These experiments show that the efficiency of this distributed model can be superlinear in the number of ants.

A. Description of the ant model

We consider a set of artificial ants moving on a twodimensional grid environment1 on which they update artificial "pheromone traces". The environment is composed of four different types of cells: one nest cell, one food cell, several bad cells and all the remaining cells are considered free. Each cell s of the grid stores two pheromone traces as two real numbers: V f (s) the food pheromone and V n (s) the nest pheromone.

Ants can carry one unit of food at a time, and are consequently in two possible states: carry food or carry nothing. Food is picked up at food cells and dropped at the nest. When a unit of food is brought to the nest a food counter, which will serve as a global performance measure, is incremented. Initially :

• the food counter is set to 0,

• the positions of the ants are initialized arbitrarily (e.g.:

all ants are initialized at the nest or uniformly at random on the grid, etc...) and all ants are set to be in the carry nothing state,

• the pheromone values are initialized arbitrarily (e.g.: 0, random, etc...) At every time step, each ant performs two actions:

1) It updates both food and nest pheromone V f (s) and V n (s) of its current cell s using the pheromone values of its four neighbors (we note the set of neighbors of cell s as N (s)), therefore using only local information.

In fact, the update requires only the knowledge of the maximum and average of both pheromone values over the neighbors. If we define:

max i (N (s)) ≡ max s ∈N (s) V i (s ),
and avg i (N (s)) ≡ 1 4 s ∈N (s) V i (s )
where i ∈ {f, n} represent which of the nest and food pheromone is used, the pheromone update rules are then as follow:

V f (s) ← 8 > < > :
1 if s is the food cell -1 if s is a bad cell β `α max f (N (s)) + (1 -α) avg f (N (s)) ótherwise Vn(s) ← 8 > < > : 1 if s is the nest cell -1 if s is a bad cell β (α maxn(N (s)) + (1 -α) avg n (N (s))) otherwise
where 0 ≤ α ≤ 1 and 0 ≤ β ≤ 1 with the condition that β < 1 if α = 1.

2) It moves to one of its four neighboring cells: with probability (0 ≤ ≤ 1) (that we shall call the exploration rate) it moves uniformly at random to one of its neighbors, and with probability 1 -, it moves to the neighbor with the highest pheromone value. The internal state of the ants conditions which pheromone is used: V f if it is in state carry nothing and V n if it is in state carry food. The β parameter can be considered as a sort of "evaporation parameter" and is typically to be set close to 1. The α parameter, which we shall call the noise parameter for reasons that will be explained later, should typically be set either close to 0 or close to 1. Two simple instances of our model correspond to the parameter choices (α = 1, 0 < β < 1) and (α = 0, β = 1). In the former choice, the general pheromone equation (the "otherwise" case above) reduces to:

V i (s) ← β max i (N (s)), (1)
and in the latter choice the equation reduces to:

V i (s) ← avg i (N (s)) (2)
which is a simple linear update. As it will appear clearly in the analysis, the question whether the ant activities (pheromone values updates and moves) are done synchronously or asynchronously is not important.

To summarize and to identify precisely different instances of this model (and for the reader to be able to reproduce the experiments we will soon describe), the following information needs to be specified:

• the environment: the set of cells, and their types (nest, food, free, bad),

• the way we initialize the position of the ants and the pheromone values V f and V n over the environment,

• the noise parameter α, the evaporation parameter β and the exploration parameter . The model we have just described is made of very simple reactive agents that communicate indirectly through the environment. In a way, this model is simpler than classical ant models: in our case, pheromones need not to be evaporated at every cell of the environment (this might represent a heavy computation when simulating a typical ACO-like algorithm). In our model, pheromone updates take place only at ant positions and not on the entire environment, and some sort of evaporation is performed through the β parameter of the local pheromone update. Nevertheless, unlike classical ant models which typically use one pheromone trail, ours uses two: following V f leads to the food source and following V n leads back to the nest. This is to compensate the fact that our ants, which are completely reactive agents, do not memorize to the path to nest. Indeed, most work using a single pheromone trail often relies on ad-hoc mechanisms to help the ants return to the nest. Other ant models where two pheromone trails are used for foraging can for instance be found in [START_REF] Panait | A pheromone-based utility model for collaborative foraging[END_REF], [START_REF] Wodrich | Cooperative distributed search: The ants' way[END_REF], [START_REF] Guerin | Emergence of constraint in self-organizing systems[END_REF].

B. Simulations

In this section we illustrate the behavior of our ant model with simulations. As we will see, we observe a form of convergence. We begin by describing what form of convergence we obtain and explain the experimental setting that allows us to measure the rate of convergence. Finally, we will briefly comment the results, as a formal in-depth analysis follows in sections III and IV.

Let us consider a typical run of the algorithm during which we fixed the parameters as follow: = 0.8, α = 0.7, β = 0.9999 and m the population size is set to 128. The ants where initialized at the nest. Figure 1 shows snapshots of this run at different time steps. On these figures can be seen the evolution of the ants moving through the environment searching for food. Once the food source is found by an ant, more ants move toward the source and a trail starts to form between the nest and the source. As time goes, a good proportion of the ants follow the trail and the dynamics seems to stabilize. Depending on the value of the noise parameter α, we can observe different forms of paths, and even sometimes no path. Figure 2 illustrates the asymptotic distribution of ants for different values of α: there are paths except for the case where α = 0.6. We will explain this phenomenon later in the analysis.

Suppose a path emerges between the nest and the food source. We do observe such a path visually, however it would be interesting to measure something objective that reveals the emergence of this path. This can be done through the food counter we introduced previously. If we plot the value of the counter with respect to time, we obtain the curve showing the increase of the accumulated quantity of food brought back to the nest over time. A typical such curve is shown in figure 3. After some time, we observe that the accumulated quantity of food brought back to the nest increases linearly. This means that the foraging behavior of the ants has converged. Section 3 will provide theoretical arguments that characterize precisely this convergence, and we will in particular explain why we don't see a path for certain values of α. Such a curve experimentally shows that there is convergence. It can also be exploited to measure some sort of convergence rate. To do so at regular times, a line is fitted on the accumulated food quantity data using a linear regression. At some point of the simulation, the parameters of the fitted line stabilize within a certain relative confidence interval. Experiments showed that using a portion (t = 10) (t = 2500) of the data (at time t we use the data registered between 0.5t to t) gives robust estimations. In a way that resembles the time constant calculus for electrical components, we define the time of convergence as the intersection of the fitted line with the x-axis once the linear regression has stabilized (see figure 3). A rate of convergence is then computed as the inverse of the time of convergence. Using this, one can compute some statistics on the rate of convergence to evaluate its dependence with respect to the number of ants. For every population size, twenty runs were carried out with the same parameter values (α, β, and m). At the beginning, ants can either be initialized at the nest or uniformly over the environment. If for a given experiment the ants were initialized at random, the same starting positions were used for each run. After each run, the rate of convergence is measured, and at the end of the experiment we compute the mean rate and its standard deviation.

(t = 10000) (t = 11000) (t = 12000) (t = 16000)
We varied the population size m from 30 to 110. Figure 4 shows a curve that represents the rate of convergence of the algorithm with respect to the population sizes. In these experiments, we fixed α = 0.7, β = 0.9999, = 0.8, and initialized the ants at the nest. The curve shows experimentally that when the number of ants increases linearly, the increase in the rate of convergence can be superlinear. In other words, when the number of ants is doubled the rate can increase by more than a factor two. This apparently surprising result will be discussed in section IV: we will provide formal arguments that explain why the rate of convergence can be superlinear.

III. CONVERGENCE ANALYSIS

In this part, we are going to give an analysis of our ant model. We will prove that it converges in some sense. To understand what our ant model does and why we see a path emerge, we will first make a detour into the theory of discrete-time stochastic optimal control, and particularly into the Markov Decision Process (MDP) framework. An MDP is a controlled stochastic process satisfying the Markov property with rewards (numerical values) assigned to states. Formally, a MDP is a four-tuple S, A, T, R where S is the state space, A is the action space, T is the transition function and R is the reward function. T is the state-transition probability distribution conditioned by the action. For all states s and s and actions a, T (s, a, s ) is the probability to go from state s to state s after executing action a at state s, T (s, a, s R(s) ∈ IR is the instantaneous reward for being in state s. In the MDP framework, the optimal control problem consists in finding a policy, that is a mapping π : S → A from states to actions, that maximizes the expected long-term amount of rewards, also called value function of policy π:

) def = Pr(s t+1 = s |s t = s, a t = a).
V π (s) = E π ∞ t=0 γ t .R(s t )|s 0 = s . (3)
We here consider an infinite time horizon; also, future rewards are discounted exponentially with a discount factor γ ∈ (0, 1) (setting γ < 1 can be seen as a mathematical trick so that the above performance criterion remains bounded). Given an MDP model, it is shown [START_REF] Puterman | Markov Decision Processes[END_REF] that there exists a unique optimal value function V * which is the fixed point of the following mapping on functions, also called Bellman operator, ∀W ∈ IR S :

[B W ] (s) = max a R(s, a) + γ. X s T (s, a, s ).W (s ) ! (4)
Once an optimal value function V * is computed, an optimal policy π * can immediately be derived as follows:

π * (s) ∈ arg max a R(s, a) + γ. X s T (s, a, s ).V (s ) ! . (5)
The fundamental reason why the above Bellman operator B has a unique fixed point is related to the fact it is a contraction mapping with a contraction factor at most γ: i.e. for all pairs of real functions U , U on S,

BU -BU ≤ γ U -U
where . is the "max-norm" on S: U = max s |U (s)|. A standard approach for solving the optimal control problem is to use a sequential iterative procedure, known as Value Iteration [START_REF] Puterman | Markov Decision Processes[END_REF], which consists in initializing a value function V 0 arbitrarily and iterating the following process: For all states s ∈ S, do:

V t+1 (s) ← [BV t ](s).
Because of the contraction mapping property, this sequence is guaranteed to asymptotically converge to the optimal value function V * , from which we can deduce the optimal controller π * (cf equation 5). Furthermore, such an iterative sequence, has a linear rate of convergence γ:

V t+1 -V * ≤ γ. V t -V * ≤ γ t+1 . V 0 -V * (6)
In [START_REF] Bertsekas | Parallel and Distributed Computation: Numerical Methods[END_REF], the authors prove that an asynchronous version of Value Iteration:

Pick (randomly) a state s ∈ S and do:

V (s) ← [BV ](s) (7)
will also converge to the fixed point V * , as long as all states keep on being picked. In the asynchronous case, one can rewrite a variant of equation 6 as:

V kt+1 -V * ≤ γ V kt -V * ≤ γ t+1 V 0 -V * , ( 8
)
where k 0 , k 1 . . . is an increasing series such that k 0 = 0 and all components of S (all states) are updated at least once between instants k t and k t+1 -1 for all t (see [START_REF] Bertsekas | Neuro-Dynamic Programming[END_REF] p. 27). Each time all the states have been updated, we know that V gets closer to V * at a linear rate γ. In fact, once again, the proof of such a result relies on the contraction mapping property. While the aim in [START_REF] Bertsekas | Parallel and Distributed Computation: Numerical Methods[END_REF] was to come up with efficient parallel implementations on real parallel computers, ours is a bit different: we are going to show that the optimal control computation fits in the (virtual) ant paradigm.

Coming back to the ant model we described earlier, we are going to make a clear link between what the ant population computes and some optimal control problems. Indeed, we will show that the pheromone values V f and V n correspond each to the value function of a control problem. This is our first result:

Proposition 1: Consider the ant model described in section II-A. If the exploration rate is strictly positive, then the pheromone value V f (resp. V n ) asymptotically converges to the optimal value function of the MDP M f = S, A, T f , R f (resp. M n = S, A, T n , R n ) with discount factor β where:

• S is the set of grid cell plus an extra "end state".

• A is the set of the four cardinal moves (north, east, south, west)

• The transition T f (resp. T n ) is characterized as follows: 1) when, in a state s corresponding to a free cell or the nest cell (resp. to a free cell or the food cell), one chooses one of the four directions a ∈ A, the probability of actually making the move in the direction a is α +

1-α 4 = 3α+1 4
while the probability of getting to each of the 3 other neighbors is 1-α 4 .

2) From all the other states, that is from bad cells, the food and the "end state" (resp. from bad cells, the nest and the "end state"), there is, for every action, a probability 1 to reach the "end state", which is an absorbing state.

• The reward R f (resp. R n ) is defined as follows: for all states s corresponding to a free cell or a nest cell (resp. a free cell or a food cell), the reward is 0.

For the state corresponding to the food cell (resp. the nest), the reward is 1. The reward is -1 for all states corresponding to bad cells and 0 for the "end state". The proof of the above result simply consists in checking that the ant model is an asynchronous version of the value iteration algorithm for the corresponding problems: concretely, we need to check that, for all cells, the updates for V f and V n are identical to the one that would be done by the Bellman operator (equation 7). This is a simple question of rewriting. For instance, to update V i (s) (i ∈ {f, n}) when s is a free cell, we have:

Vi(s) ← β (αmaxi(N (s)) + (1 -α)avg i (N (s))) Vi(s) ← β 0 @ α max s ∈N (s) Vi(s ) + 1 -α 4 X s ∈N (s) Vi(s ) 1 A (9) ⇔ Vi(s) ← β max s ∈N (s) 0 @ αVi(s ) + 1 -α 4 X s ∈N (s) Vi(s ) 1 A ⇔ Vi(s) ← β max s ∈N (s) " α + 1 -α 4 « Vi(s ) + 1 -α 4 X s ∈N (s)\{s } Vi(s ) 1 A ⇔ Vi(s) ← β max a∈A X s ∈S Ti(s, a, s )Vi(s) ! ( 10
)
We let the interested reader check that this equation also holds in the remaining (limit) cases. Finally, the condition > 0 ensures that all states will keep on being visited and updated, which in turn ensures the convergence of this asynchronous version to the optimal value functions. Now let us interpret what this means and especially why we observed the emergence of paths between the nest and the food source in the experiments. Solving M f (resp. M n ) means finding a policy that will generate trajectories that avoids (on average) the bad cells for which the reward is -1 and try to reach the food cell (resp. the nest cell) for which the reward is 1; because of the discount factor β < 1, the optimization also tries to minimize the time to reach the food cell (resp. the nest cell). In other words, M f (resp. M n ) are natural formulations of the control problem: "go to the food cell (resp. to the nest cell) while avoiding the bad states" assuming that there is some noise in the transition models T f (resp. T n ). The level of noise is related to α: calling α a "noise parameter" was thus fully justified.

In each state, the optimal action is the one for which the maximum is reached in equation 10 above, and it is easy to see that this optimal action is the one that points to the cell for which the maximum is reached in equation 9: that is the cell with the highest pheromone value. The pheromone values asymptotically converge to the corresponding optimal value functions V f and V n . Therefore, the moves of the ants, which are (recall the ant move description in section II-A) a mixture of random uniform moves (with a weight ) and the action that climbs up the pheromone value (with a weight 1 -), converge to a mixture of random uniform moves and the optimal corresponding moves. The smaller , the clearer will be the path when the pheromones have converged (at the limit if = 0 all the ants follow the optimal policy). However, the smaller , the longer it will take for the ants to repeatedly visit all the states and make the pheromone values converge. This trade-off is typical of optimal control theory and is known as the exploration-exploitation dilemma (see for instance [START_REF] Sutton | Reinforcement Learning, An introduction[END_REF]).

In the description of the algorithm in section II-A, we wrote about the parameters α and β that we needed 0 ≤ α ≤ 1 and 0 ≤ β ≤ 1 with the condition that β < 1 if α = 1. We can now explain these conditions from the optimal control viewpoint: setting the discount factor of a MDP to a value strictly lower than 1 ensures that the infinite time horizon performance (eq. 3) remains bounded and that the Bellman operator is a contraction mapping. As soon as there is some amount of noise (i.e. α < 1) in the transition model of our specific MDP M f and M n , then for any action, there is a probability 1 of reaching one of the absorbing "end state" (with zero reward) and this suffices to ensure that the performance criterion remains bounded and the Bellman operator is contracting. However, in the purely deterministic case (α = 1), the discount factor β needs to be set to a value strictly lower than 1.

We can explain further the influence of the noise parameter α. When α is equal to 1 (and the update is eq. 1), there is no noise in the transition model and the optimal policies exactly match the shortest path (with the respect to the Manhattan distance, see figure 2, case α = 1) between the nest and the food source. When we decrease α, the level of noise increases and the optimal policies get smoother (they try to stray away from the bad states). At some point, when α goes on decreasing, there is so much noise that the probability of first reaching a bad state when trying to reach the food source (or the nest) gets really large whatever the policy.

Consequently the optimal behavior consists in simply staying away from the bad cells and not trying to reach the food source (or the nest): in this case it is better to have a 0 reward than a -1 (which happens when one hits a bad state). This explains why there was no path in figure 2 for α = 0.6. Furthermore when α goes on decreasing and gets close to 0, something apparently strange happens: paths appear again. We can give two explanations of this: 1) when the noise gets so large that the influence of the actions nearly vanish, the optimal controller cannot even prevent from hitting a bad state and, as a kamikaze that would know he is going to die anyway, it again becomes interesting to try to reach the food source (or the nest). 2) At the limit when α = 0 (and the update is eq. 2), the pheromone potentials, which satisfy an equation of the type V i (s) = avg i (N (s)) is equal to a discrete harmonic function and it is known that a harmonic function can be used for navigation (see [START_REF] Connolly | On the applications of harmonic functions to robotics[END_REF], [START_REF] Boumaza | Optimal control subsumes harmonic control[END_REF] for further details).

IV. RATE OF CONVERGENCE ANALYSIS

We saw that our ant model converges by arguing that it is an asynchronous computation of two contraction mapping fixed points: pheromone potentials that results from the ants local updates stabilize towards the optimal value functions of some control problems, which guide the ants between the nest and the food source. The aim of this section is to study the rate of convergence of this process with respect to the number of ants. We shall describe some objects and results related to Markov chains on graphs that highlight the experimental observation (made in section II-B) that this rate of convergence can also be superlinear: doubling the number of ants may accelerate the process by more than a factor two.

To study the rate of convergence, we go back to equation 8 that describes, in general, the convergence rate of the asynchronous computation of a contraction mapping fixed point. We rewrite it here for the sake of clarity: [START_REF] Feige | Deterministic approximation of the cover time[END_REF] Recall that k 0 , k 1 . . . is an increasing series such that k 0 = 0 and all components of S (all states) are updated at least once between instants k t and k t+1 -1 for all t. The rate of convergence is thus clearly related to the random variable k t : the slower k t grows the faster the process converges. Since we are interested in the rate with respect to the number m of ants, we can make this dependence explicit by writing k m t . What we need to study is thus k m t+1k m t . The good news about k m t+1k m t is that it can be related to a known object of the probability literature. Consider the expectation of this quantity when m = 1:

V kt+1 -V * ≤ γ V kt -V * ≤ γ t+1 V 0 -V * ,
E[k 1 t+1 -k 1 t ]
is the average time for one ant to visit all the cells of the environment. More formally, if we see the environment as a graph G (there is one node for each cell and a connection when two cells are neighbors) it is the average time for a Markov chain (that describes the positions of the ant) to visit all the nodes of the graph G, and such a quantity is usually called the cover time of the Markov chain on the graph G [START_REF] Aldous | Reversible Markov Chains and Random Walks on Graphs[END_REF]. Similarly, for any m, E[k m t+1k m t ] is the average time for several parallel Markov chains to visit all the nodes of the graph G and it is known as the cover time of the parallel Markov chains on the graph G.

Now comes the bad news about k m t+1k m t : it is in general very hard to compute the cover time of a graph for a given Markov chain and a given initial distribution: they are computable in exponential time and it is not known whether it is possible to approximate them in deterministic polynomial time [START_REF] Feige | Deterministic approximation of the cover time[END_REF]. It is obviously harder to compute the cover time of a graph for several Markov chain. For our specific ant model, it is even harder to compute the cover times since the Markov chains and the distribution of ants vary over time in a non-trivial way: the transition probabilities depend upon the pheromone potentials which themselves are continuously updated by the Markov chains. There is therefore little hope that one could estimate very precisely the rate of convergence of our ant model. A general asymptotic study of the convergence (where one considers that 1) the pheromone values have converged and 2) the ants have reached their stationary distribution) may be pursued and this is a possible subject for future research.

Nonetheless, by studying the literature on the cover time, we were able to find the following interesting result [START_REF] Aldous | Reversible Markov Chains and Random Walks on Graphs[END_REF]:

Proposition 2: On a regular n-vertex graph 2 , consider m independent balanced random walks 3 , each started at a uniform random vertex. Let C m be the time until every vertex has been hit by some walk 4 . Then as the size of the graph n → ∞ and while the number of walks satisfies m ≥ 6 log n, we have:

E C [m] ≤ (25 + o(1))n 2 log 2 n m 2 .
The interesting point is the 1 m 2 dependence: this implies that (as n → ∞ and m ≥ 6 log n) the cover time is bounded by a function which is inversely quadratic in the number of walks. On simple graphs like the ncycle, [START_REF] Aldous | Reversible Markov Chains and Random Walks on Graphs[END_REF] explains that the above bound is sharp, so in such a case, the cover time is inversely quadratic in the number of walks.

Using the above discussion, we can "translate" the above proposition into something that is meaningful for our ant framework:

Proposition 3: Consider the ant model described in section II-A on a toric grid environment with n cells, with an exploration rate = 1 and initialize the ants uniformly on the environment. When the size of the environment n → ∞ and while the number m of ants satisfies m ≥ 6 log n, the time for the pheromone values to reduce their distance to their limit by a factor β is bounded by (25+o(1))n 2 log 2 n m 2 that 2 A regular graph is a graph where all nodes have the same degree, i.e: each node is connected to the same number of nodes.

3 A balanced random walk on the graph is such that transitions are uniform distributions on neighbors.

4 E[C m ] is thus the cover time of the graph by these parallel random walks is a function that is inversely quadratic in the number m of ants. We consider a toric environment so that the corresponding graph is regular. We take = 1 so that the ants follow a balanced random walk. Also notice then that the uniform initialization is the stationary distribution of the random walks, so that the distribution of ants is stationary over time. Our proposition is thus simply a corollary of proposition 2.

The bound is sharp when the environment is a n-cycle graph, that is a long (cycled) corridor: in such a case, the rate of convergence can be quadratic in the number of ants. We there have a superlinear rate of convergence. The experiments we made suggest that superlinearity also happens for other values of and general environments. Extending the above results to these more general setting constitutes future research.

V. DISCUSSION

We have presented an ant model that can be related to the framework of optimal control, and proved that it converges in some sense. We have also studied the rate of convergence with respect to the number of ants: we showed, through experimental and analytical arguments, that the rate of convergence can be superlinear in the number of ants. At first sight, superlinearity may look like an impressive property. There is however no magic: the fundamental reason why we have such a superlinear rate is due to the fact that small populations of our ants are particularly inefficient at visiting the whole state space, and therefore making the pheromone potentials converge. In fact, the convergence analysis of contraction mapping that we used clearly suggests (compare equations 6 and 8) that the fastest method for computing the optimal value function is the synchronous version: at each step, the value is approaching its limit with the linear rate γ. The relative slowness of the asynchronous version has to be understood as the price for decentralization.

The ant model we have presented in this paper is very flexible. It can incorporate many variations. For instance, one could consider that there are several food sources, each containing a finite quantity of food, which decreases over time as ants come and go. This quantity could be incorporated in the corresponding optimal control model, through the reward function on the food states. The reward would then evolve over time, and the pheromones, which are continuously being updated, would keep on following the corresponding "moving" optimal value function. One could also use different parameters for the back and forth trips of the ants: this could model different strategies depending on the fact that the ant carries food or not. Finally, if the environment is static, one could make the exploration rate slowly tend to 0, so that the ants eventually converge to the deterministic optimal policies. Studying good ways for "freezing" the ant behavior through the exploration parameter constitutes future research.

It is in fact easy to see that one could construct similar ant algorithms for almost any optimal control problem. As long as the problem is formulated as a MDP, we know that it can be solved by an asynchronous version of Value Iteration, and building the corresponding ant system is immediate: we just need to have ants move around the state space and do the Bellman update everywhere they go. If we care about the "locality constraint", i.e. the constraint that "ants should make their decisions only using local information" then this approach will work as long as the transitions in the MDP are also local in the state space. We can even go further: All our analysis (the convergence and the rate of convergence) relies on the "contraction mapping" property. This suggests that any problem that involves the computation of a contraction mappings (contraction mappings can for instance appear in zero crossing problems, constrained and unconstrained optimization [START_REF] Luenberger | Linear and nonlinear programming[END_REF]) on a finite space has a natural (potentially superlinear) ant-like solution: ants move around on this finite space and do local contraction updates.

The ant model we have presented is closely related to previously published algorithms for the foraging problem. In [START_REF] Simonin | Construction of numerical potential fields with reactive agents[END_REF], the author presents an algorithm that compute shortest path using an asynchronous implementation of the Bellman-Ford algorithm: the local update is of the form U (x) ← 1 + min x ∈N (x) U (x ). Modulo the variable change U ↔ log(V ) log(β) this is equivalent to our pheromone update with α = 1 expressed in equation 1. In [START_REF] Panait | A pheromone-based utility model for collaborative foraging[END_REF], the authors address the foraging problem and identify, as we proposed, the value function of control problems as pheromone rates. They use an asynchronous version of Value Iteration as well as certain variants of algorithms in reinforcement learning [START_REF] Sutton | Reinforcement Learning, An introduction[END_REF], which is also formulated through MDPs, to solve the optimal control problem. Furthermore, they present many experiments which illustrate that these models are robust in environments with obstacles and where the positions of the nest and the food source vary over time. As in our model, the authors use multiple pheromone trails for the different paths to find.

There are however significant differences between these works and ours. 1) They focus only on the foraging problem whereas we have just argued that in principle, any optimal control problem (and any computation of contraction mapping) could be tackled by such ant-like algorithms. 2) To our knowledge, there are no arguments concerning the convergence of the models in [START_REF] Panait | A pheromone-based utility model for collaborative foraging[END_REF] whereas an extended version of [START_REF] Simonin | Construction of numerical potential fields with reactive agents[END_REF] presents a convergence proof, that uses arguments close to the ones we give here [START_REF] Simonin | Collective construction of numerical potential fields for the foraging problem[END_REF]. 3) Last but not least, if the superlinearity of the convergence rate is experimentally observed in [START_REF] Simonin | Construction of numerical potential fields with reactive agents[END_REF], we are to our knowledge the first ones to study it from an analytical point of view.

VI. CONCLUSION

We have described a foraging ant model that can be related to the optimal control framework. This allowed us to prove its convergence in some well defined sense and to discuss its rate of convergence with respect to the size of the ant population. Our analysis, based on the idea that the ant population computes the fixed point of a contraction mapping, explains many results obtained through simulations, among which the potential superlinearity of the rate of convergence. Eventually, our work suggests that any problem that involves the computation of a contraction mapping fixed point (e.g. optimal control in general, contraction, zero crossing problems, constrained and unconstrained optimization) on a finite space has a natural similar (potentially superlinear) ant-like solution.

Fig. 1 .

 1 Fig. 1. Snapshots of the ant model: filled squares represent ants carrying food from the food source (bottom left) back to the nest (top right), and empty squares are ants not carrying food.

αFig. 2 .

 2 Fig. 2. Limit distribution of the population for various values of α.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Measuring the rate of convergence: food units correspond to the accumulated amount of food brought back to the nest normalized by the number of ants, and the convergence line corresponds to the linear regression over the data.

As the reader will understand in the following, more complicated graph structures could be used, we address here the

2D grid case for the sake of simplicity.


P Pinson

G Kariniotakis
email: georges.kariniotakis@ensmp.fr

On-line Assessment of Prediction Risk for Wind Power Production Forecasts

Keywords: Wind power, short-term forecasting, confidence intervals, prediction risk, on-line software, adaptive fuzzy-neural networks, numerical weather predictions

The paper introduces a new methodology for assessing on-line the prediction risk of short-term wind power forecasts. The first stage of this methodology consists in computing confidence intervals with a confidence level defined by the end-user. The resampling approach is used, which in contrast to existing methods for wind forecasting, does not need to make a restrictive hypothesis on the distribution of the errors. To account for the nonlinearity of the power curve and the cut-off effect, the errors are classified using appropriate fuzzy sets. The confidence intervals are then fine-tuned to reduce their width in case of stable weather conditions. For this purpose, an appropriate index, named as MRI, is defined reflecting the spread of the available Numerical Weather Predictions. A linear relation between that index and the resulting prediction error is shown. The second part of the methodology is to use the MRI-index itself as a preventive on-line tool to derive signals for the operator on the meteorological risk, i.e. the probabilities for the occurrence of high prediction errors depending on the weather stability. Evaluation results of this methodology over a one-year period on the case study of Ireland are given, where the output of several wind farms is predicted using a dynamic Fuzzy-Neural Networks based model. The proposed methodology is generic and can be applied to all kinds of wind power prediction models.

I. INTRODUCTION

OWADAYS, wind park installations in Europe exceed 23 GW, while the motivated by the Kyoto Protocol targets of the E.U. for 12% energy demand covered by renewables by year 2010, are translated to 21% electricity generation by renewables. To achieve these targets, wind power in the Member States should arise up to 45-60 GW. Such a large-scale integration of wind generation causes several difficulties in the management of a power system. Often, high level of spinning reserve is allocated to account for the intermittent profile of wind production, thus reducing the benefits from the use of wind energy. Predictions of wind power production up to 48 hours ahead contribute to a secure and economic power system operation.

Apart from spot forecasts of the wind parks output in the next hours, of major importance is to provide tools for assessing on-line the accuracy of these forecasts. Tools for on-line evaluation of the prediction risk are expected to play a major role in trading wind power in a liberalized electricity market since they can prevent or reduce penalties in situations of poor prediction accuracy. In practice today, uncertainty is given in the form of confidence intervals or error bands around the spot wind power predictions.

Typical confidence interval methods, developed for models like neural networks [START_REF] Alves Da | Confidence intervals for neural network based short-term load forecasting[END_REF]- [START_REF] Heskes | Practical confidence and prediction intervals[END_REF], are based on the assumption that the prediction errors follow a Gaussian distribution. This however is often not the case for wind power predictions, where error distributions appear some skewness, while the confidence intervals are not symmetric around the spot prediction due to the shape of the wind turbines power curve. Moreover, the level of predicted wind speed introduces some nonlinearity to the estimation of the intervals; i.e. at the cut-off speed, the lower power interval may switch to zero due to the cut-off effect. The limits introduced by the wind farm power curve (min, max power) are taken into account by the method proposed in [START_REF] Luig | Analysis of confidence intervals for the prediction of regional wind power output[END_REF], which is based on modelling errors using a ß-distribution, the parameters of which, have to be estimated by a post-processing algorithm. This approach however is applicable only to "physical"-type of models since such models estimate power using an explicit wind park power curve, which is not necessarily the case for statistical or artificial intelligence based models as the ones considered here [START_REF] Kariniotakis | An advanced on-line wind resource prediction system for the optimal management of wind parks[END_REF].

In [START_REF] Lange | Assessing the uncertainty of wind power predictions with regard to specific weather situations[END_REF], [START_REF] Lange | Accuracy of short term wind power predictions depending on meteorological conditions[END_REF] wind speed errors are classified as a function of look-ahead time and then they are transformed to power prediction errors using the wind turbine power curve vs. wind speed. This method however is also limited for application to physical models rather than statistical ones since it requires local wind speed predictions (at the level of the wind farm), while it does not provide uncertainty as a function of a pre-specified confidence level. The wind speed errors are estimations provided by the Numerical Weather Prediction (NWP) model. As a consequence, this method does not take into account the modelling error itself that might be due to the spatial refinement of weather predictions or to the power curve used. On the other hand, wind speed measurements are required, which, in general, might not be available on-line.

This paper proposes a methodology for assessing on-line the uncertainty of wind power predictions by the joint use of appropriately defined confidence intervals and prediction risk indices.

Initially, a generic approach is developed for the estimation of confidence intervals that can be applied to both "physical" and "statistical" wind power forecasting models. This is due to the fact that no hypothesis is made about the distribution of the prediction errors. The method accommodates both modelling errors and errors based on the NWPs. It uses past wind power data, which are often available on-line by a Supervisory Control and Data Acquisition (SCADA) system, as well as NWPs, which are the basic input to all models.

Generally, when confidence higher than 80% is required, the N intervals are quite wide. This can lead to conservative or costly managing strategies of the predicted wind power (i.e. allocation of high spinning reserve). Given that confidence intervals are estimations of the uncertainty based on the past performance of the model, the second objective of this work is to develop additional preventive tools able to assess on-line the prediction risk as a function of the forecasted weather situation. The aim is then to provide comprehensive information to the operators so that they are able to adjust the risk they are going to undertake when managing the predicted wind power, i.e. take low risk when forecasted weather situation is unstable. The paper presents detailed results from the application of the method on the case study of Ireland, where the aim is to predict the output of several wind farms for 48 hours ahead using on-line measurements and predictions from Hirlam NWP system. Evaluation is based on one-year of data.

II. UNCERTAINTY IN WIND POWER PREDICTIONS

Let us define the prediction error for the look-ahead time t+k as following: This error can vary between -100% and 100% of the nominal wind park power. For a non-bounded prediction model it can take values even outside that range. The possible error of the prediction model, defined as "error margin", depends on the level of measured wind power. Fig. 1 represents graphically the error margin as a function of the wind park characteristic curve.

t k t k t t
For wind speeds below cut-in speed, the error margin is maximal since the model can predict a production up to the nominal wind park power. On the contrary, for higher wind speeds the model will show a positive error margin, i.e. the generated power is likely to be greater than the one proposed by the prediction model. Close to the cut-off wind speed the uncertainty is again maximal since the model can switch from a positive error margin to a negative one, or the inverse.

III. METHODOLOGY FOR ESTIMATING THE UNCERTAINTY OF WIND POWER PREDICTIONS

In this Section, a methodology is developed that permits to assess on-line the prediction risk. It includes: (i) Development of confidence intervals for the spot power prediction. The approach is based on the resampling method, which is applied on samples of errors. Errors are classified using fuzzy sets to account for the level of power and the risk for cut-off events. (ii) Development of an on-line prediction risk index based on ensembles of NWPs. This Index permits to derive rules for assessing the probability of high or extreme prediction errors due to unstable weather situations. (iii) Dynamic fine-tuning of the size of the intervals depending on the weather stability. This permits to avoid excessive risk or to take preventive actions in situations where high errors are expected.

A. Error pre-processing based on fuzzy sets modelling

The first step before computing confidence intervals is to collect the prediction errors the model made in the past. The intervals that are going to be computed will rely on the most recent information on the model's performance. For this, a window in the past (a certain number of hours) is defined and used as a sliding window for storing the errors. The size of this window defines the size of the errors sample. A separate sample is developed for each look-ahead time of the prediction horizon (i.e. for 1-hour ahead, 2-hour ahead, and so on). This is because prediction errors depend on the look-ahead time as illustrated in Fig. 2. The collected errors are the most recent ones at a given time: when the actual measured wind power is known, that value is compared with all the past predictions made for that time (from 48 hours to 1 hour ago).

Fig. 2. The distribution of the prediction errors varies as a function of the prediction horizon (left picture: 1-hour ahead prediction error distribution, right picture: 24-hour ahead prediction error distribution). Results are produced using data for a single wind farm in Ireland.

The power prediction errors depend on the errors involved in the prediction of wind speed by the NWP system [START_REF] Lange | Assessing the uncertainty of wind power predictions with regard to specific weather situations[END_REF][START_REF] Lange | Accuracy of short term wind power predictions depending on meteorological conditions[END_REF][START_REF] Landberg | Poor-man's ensemble forecasting for error estimation[END_REF]. Due to its shape, the wind park power curve is able to amplify (between cut-in and rated speed) or to reduce (below cut-in speed or between rated and cut-off speed) the uncertainty introduced by the NWPs. To account for this effect, the wind power curve is divided into three ranges of power: low, medium and high, which are characterised by fuzzy sets. The prediction errors are classified then as a function of these three ranges (Fig. 3). Hence, the confidence interval estimation is carried out using the error samples corresponding to the power class of the predicted power.

In a similar way, in order to deal with the risk due to the cutoff event, the universe of discourse of wind speeds is divided into two ranges corresponding to a "no cut-off risk zone" for low wind speeds, and to a "cut-off risk zone" for wind speeds close or higher than cut-off. An appropriate fuzzy set is associated to each zone as shown in Fig. 3. The cut-off risk is considered in priority compared to the power class splitting.

B. Confidence interval estimation by the resampling approach.

Here is a formal definition of confidence intervals: the interval computed from the sample data which, were the study repeated multiple times, would contain the true effect CL% of the time, CL being the confidence level.

A given set of observations (the sample) is a part of a whole population and can be seen as representative. The aim of methods like resampling is to have a better idea of the population distribution by going through the sample a high number of times. This evaluation of the population distribution can serve to estimate a mean, a variance, etc. No assumption is made concerning the distribution: that is the main reason why resampling is preferred to other methods like methods based on simple standard deviation for instance.

Let us consider a sample containing N observations of a mean m for a given stochastic process. The procedure to compute from this sample the 95% confidence interval for that mean involves the following steps:

(i) N values are selected randomly and with replacement out of the original sample in order to create a new sample; (ii) the new sample is sorted in ascending order; (iii) the 2.5% lowest and 97.5% highest value of that set are determined.

Fig. 4. One step of the resampling process.

This procedure represents one step of the resampling process (Fig. 4). Indeed, these three actions are to be repeated a large number of times to re-create the population again and again. One always gets a new sample that is close to the original one, and the whole population distribution is not really simulated by this way.

But, by calculating the mean of respectively the 2.5% lowest and the 97.5% highest value of these randomly created samples, good estimates of the confidence limits for the mean m can be computed [START_REF] Alves Da | Confidence intervals for neural network based short-term load forecasting[END_REF], [START_REF] Masters | Neural, Novel & Hybrid Algorithms for Time-Series Prediction[END_REF], [START_REF] Simon | Resampling : The new statistics[END_REF].

In the case of wind power forecasting, the resampling method is applied by considering error samples defined as a function of the look-ahead time, the power range and the wind speed range. These error samples are treated one after the other using the procedure described above, assuming that the prediction error the model makes is the mean of a distribution and that we would like to compute confidence intervals for that mean.

IV. INFLUENCE OF NUMERICAL WEATHER PREDICTIONS

Not only the power prediction model can be responsible for low quality forecasts, but also the numerical weather prediction system (i.e. due to low weather stability). Indeed, an unstable atmospheric situation can lead to very poor numerical weather predictions and thus to worthless wind power ones. In contrast, when the atmospheric situation is stable, one can expect more accurate wind power predictions from the model.

In general it is very difficult to validate the accuracy of numerical weather predictions since wind speed measurements are often not available on-line. This was the situation in the casestudy examined in this paper, where only wind power measurements are available in parallel to Hirlam NWPs of wind speed at four levels. In order to study the relevance of the NWPs, a dynamic approach is developed based on their correlation to the measured power. The aim is to estimate the probability of situations where Hirlam fails to predict local conditions for a certain period of time (i.e. due to local weather situations). For this purpose, cross-correlations between wind power and Hirlam wind speeds were estimated using a sliding window of 100 hours. The result is plotted in Fig. 5 for a period of 2000 hours. The distribution of the obtained correlations is shown in Fig. 6 for each level. The range of the values is between {-0.4, 0.92}. Low values indicate situations where low reliability of Hirlam forecasts might be expected. The frequency of these periods is however limited since the distributions are centred on the 0.8 correlation value. An additional conclusion is also that the most correlated Hirlam level varies in time. This is an indication that a performing prediction model should consider more than one level as input.

In the following Paragraphs, the information included in the NWPs is exploited to develop tools for on-line estimation of the uncertainty in power predictions.

A. Ensemble forecasts for the assessment of weather stability

Meteorological Centres are able to produce different scenarios of Numerical Weather Predictions by perturbing the initial conditions of the forecasting model or by using different NWP models. These scenarios are called ensemble forecasts and permit to evaluate the stability of the weather regime [START_REF] Roulston | Value of the ensemble prediction system for forecasting wind energy production[END_REF]. However, for wind power applications only one forecast for the next 48 hours is often made available (or purchased) at a given time. For instance, Hirlam gives a unique 48-hour ahead forecast every 6 hours. Nevertheless, for a given hour, several predictions can be available from different time origins in the past (-6, -12, -18… hours). In a stable and wellpredicted weather situation it is expected that these predictions will not differ significantly. Weather stability can be assessed by comparing all the available forecasts for the considered period. In order to have a general evaluation of that stability, 4 sets of predictions of various ages (0, 6, 12 & 18 hours) for the following 24 hours are compared. Figure 7 gives the examples of a stable atmospheric situation (left picture, the forecasts are quite close) and of an unstable one (right picture, spread forecasts).

B. Development of a norm to assess the weather stability

There are several possibilities to measure the spread of the various weather forecasts. In [START_REF] Landberg | Poor-man's ensemble forecasting for error estimation[END_REF] the standard deviation of the forecasts for each time-step is mentioned as an example.

Our aim here is to evaluate the global atmospheric situation. This is why a unique representative index is defined for the following N h hours instead of indexes for every look-ahead time.

In order to calculate the distance between two sets of forecasts, we propose a 2-norm between the N h -valued vectors containing the predicted wind speed for the N h following hours. Lets define as

α α α α - + - + - + -= K K
the N f available sets of wind speed forecasts (NWPs), with α being the age of each set. The values for α can be 0, 6, 12, etc, for the case of Hirlam. The distance between the predictions of ages α i and α j is given by: ( )

2 1 1 2 1 ,               - = ∑ = - + - + - - h N k j t k t i t k t h j t i t ws ws N ws ws d α α α α (3)
Then, an index, called hereafter "meteo-risk" MRI-index, is defined to measure the spread of the weather forecasts at a given time. It uses the most recent forecast as a reference and reflects the variability of the older forecasts:

( )

i t t f N i i ws ws d p MRI α α - - - = ⋅ = ∑ , ˆ0 1 1 (4)
with i p (i = 1,…,N f -2) being appropriate weights defined so that:

1 + > i i p p , 2 , , 1 - = f N i K ∑ - = = 1 1 1 f N i i p (5)
The use of the weights i p permits to give more importance on the recent information we get from the weather predictions.

C. Relation between weather stability and wind power prediction error

In the frame of the case study of the paper, the horizon N h for the calculation of the MRI-index is set to 24 hours. Since Hirlam forecasts are provided every 6 hours, there are four sets (N f =4) of wind speed predictions covering the period. However, the same methodology can be applied to seven available sets of Hirlam forecasts on a 6-hour period for instance.

Figure 8 shows the distribution of the values of the MRI-index as obtained using the data of one of the considered wind farms in Ireland (see Section V). It can be concluded that most of the times the weather regime is quite stable (low MRI-index values) and that there are only few occurrences of really spread forecasts. Figure 7 shows the patterns of a "stable" and an "unstable" weather situation with MRI-index values of 0.3 and 2.9 respectively.

Fig. 8. The distribution of weather situations as expressed by the MRI-index.

For the same case study we collect wind power prediction errors as obtained by a fuzzy-neural network (F-NN) model [START_REF] Kariniotakis | An advanced on-line wind resource prediction system for the optimal management of wind parks[END_REF] for a period covering one year. For the same period the MRI-index is estimated. By binning the data, calculating the average error for each bin, and comparing these averages to the global prediction error of the model, the representative points in Fig. 9 are obtained. The prediction error increases linearly with the MRI-index: the tighter the Hirlam predictions are, the more accurate the wind power prediction model is. A linear fitting gives the solid curve shown in Figure 9.

Making this assumption would mean that the prediction error the model made follows an affine empiric relation: e = e 0 + s.MRI (6)

Fig. 9. Prediction errors vs. MRI-index over a one year dataset: there is a linear relation between the prediction error and the MRIindex.

which is composed by a basic part of the error e 0 and by a NWPdependent error, the latest being a direct consequence of the prediction model sensibility to the weather stability. The slope s of the linear fitting model represents that sensibility.

D. Use of the MRI-index to adjust on-line the prediction risk 1) Fine-tuning of the confidence intervals

The relation [START_REF] Luig | Analysis of confidence intervals for the prediction of regional wind power output[END_REF] indicates that when the MRI-index is low, the model is expected to be more accurate. In that case one would be ready to accept tighter confidence intervals for the predictions. The aim here is to use Eqn. [START_REF] Luig | Analysis of confidence intervals for the prediction of regional wind power output[END_REF] to define a scale factor for the confidence intervals depending on the value of the MRI-index. This scale factor can be applied to either enlarge or narrow the intervals width in the following N h hours. For instance, when the meteorological index equals 0.5, the size of the intervals for the following 24 hours is reduced by almost 20%. The strategy chosen here is to only narrow the intervals when the MRI-index permits to do so. It can be seen from Fig. 8 that most (~65%) of the time, the MRI-index allows to decrease the interval size.

2) On-line use of the MRI-index for monitoring prediction risk

The relation that was drawn above permits to define rules concerning the expected prediction error depending on the MRIindex value. For that purpose, we bin the data and calculate the cumulative distribution function of the prediction errors for each bin (Fig. 10). These curves give the probability with which, an error larger than a defined threshold occurs, depending on the value of the MRI-index. For instance, if at a certain time, the index takes a value between 0 and 1, there will be a probability of 36% that an error t k t e + larger than the mean prediction error

t k t e + occurs.
However, if at that same time the value of the index is within the [1.75-2] bin, the probability for the same kind of error is twice larger (75%):

IF MRI∈[1.75-2] THEN t k t e + > t k t
e + with 75% probability. (7) Table I gives the probability for errors to be larger than 1, 1.5 and 2 times the average error depending on the range of the MRIindex. The Table is estimated for the case of 5 wind farms in Ireland. Based on such a Table, several rules similar to the one given by ( 7) can be derived.

Table I also provides information on the probability of extreme prediction errors to happen (as extreme are defined errors larger than twice the mean prediction error of the model). Actually, for these 5 wind farms in Ireland, when the MRI-index takes low values (between 0 and 1) an extreme prediction error is unlikely to happen, and that is not the case if this one is within the bin [1.75-2] (17% probability of occurrence). On the other hand, if MRI>1.75, an error of at least 50% of the average prediction error is expected.

In an on-line environment these rules permit to derive signals or alarms for the end-user of the wind power prediction model, informing that large prediction errors might occur. Then, the operator can consider such a signal for: • taking preventive actions (i.e. increase spinning reserve),

• considering the lower interval, rather than the spot prediction of power, when trading, in order to avoid penalties, etc.

In the developed software, the MRI-index is monitored as shown in Fig. 14. The passage of the index to high zones triggers messages to the operator on the expected risk. The rules that are implemented evolve over time thanks to the storage of MRI-index values and prediction errors. Table I. Rules for the occurrence of larger errors depending on the value of the MRI-index.

Ranges of MRI-index

Probability of occurrence of errors larger than n times the mean prediction error

Bin

V. RESULTS

Results are presented for five real wind farms in Ireland (WF-A to WF-E) with a total installed power of a few tens of MW. The prediction model is the Adaptive Fuzzy-NN model described in [START_REF] Kariniotakis | An advanced on-line wind resource prediction system for the optimal management of wind parks[END_REF]. The available time series cover a period of almost two years from which 6600 hours were used for training (learning set), 1000 hours for cross-validation and one year for testing the performance of the model. The results presented here are on the testing set. Concerning the computation of confidence intervals, 12 days of prediction errors are stored in the samples. The desired confidence level is set to 85%.

Figure 11 depicts an episode with the wind power predictions for the next 43 hours compared to the real values for WF-E. The 85% confidence intervals are built with the method described above. In order to illustrate the fine-tuning of the intervals, Fig. 12 gives the example of a weather situation classified as stable with respect to the "meteo-risk" index. For the first 24 look-ahead times the resampling on past errors produces quite broad intervals, but their size is reduced by almost 20% afterwards thanks to the consideration of the weather situation.

The resampling approach that is used to design the confidence intervals displays a dynamic behaviour thanks to the updating of the sample of errors and also thanks to the fuzzy set modelling of the power curve. Figure 13 shows the evolution of the 6-hour ahead interval width during 1000 hours of the evaluation period. An interesting point is that these intervals are non-symmetric because they follow the concepts we described about the error margin -for instance the prediction error is more likely to be negative than positive when forecasting low power output.

Moreover, the interval size varies a lot over time. For the example of Fig. 13, the width of the 6-hour ahead confidence interval ranges from 10% to 70% of the wind farm nominal power. Depending on the weather stability described by the MRIindex, the intervals are tuned: when the weather situation is considered as stable (low index values) they are narrowed up to 30%. The dashed curves in Fig. 13 correspond to the final intervals.

The performance of the confidence interval estimation for the five wind farms in Ireland is summarized in Tab. II, i.e. the observed confidence at the end of the evaluation period (one year) for both resampling and fine-tuned intervals, as well as the percentage of reduced intervals and their average reduction due to the use of the MRI-index. One can see that the consideration of the weather stability permits to narrow the intervals most of the times (~65% in column 3), and the average reduction is up to 11% of their initial size (WF-A). The corresponding confidence loss is not significant (column 2). The interval size reduction is illustrated in Fig. 14, where the average interval width (for the 6-hour ahead ones) is plotted against the observed confidence. The various interval sizes are quite important, but one has to keep in mind that this study deals with single wind farms and so the level of prediction error is higher than for the case of regional or national wind power predictions where there is some spatial smoothing effect.

Finally, Figure 15 gives a picture of an episode with the evolution of the MRI-index and the related prediction error: the period of low "meteo-risk" (time in the evaluation set between 200 and 450) matches with a period of very low prediction error and periods with unstable weather situations (time between 80 and 180 for instance) correspond to higher prediction errors. These features, as well as the mean prediction error of the model, can be given on-line to the operator and help him to seize the prediction risk due to the meteorological situation.

VI. CONCLUSIONS

A generic methodology for assessing on-line the prediction risk of short-term wind power forecasts is presented. At a first stage, confidence intervals based on the resampling approach are derived, taking into account the prediction horizon, the power class and the cut-off risk. At a second stage a new meteorological risk (MRI) index was introduced to evaluate the weather stability. The MRIindex can be used either to fine-tune the confidence intervals or to give signals to the end-user on the probability of outliers.

The developed methodology was tested over a one-year evaluation dataset for five wind farms located in Ireland. The results are encouraging and consist a first step in the development of on-line tools that can be used in a complementary way to the prediction model itself.

The developed prediction risk method is implemented in the form of on-line module and integrated in the Armines Wind Power Prediction System (AWPPS). In the future this module will be also integrated and available through the ANEMOS advanced prediction platform (http://anemos.cma.fr). Apart from the prediction risk module, the core F-NN wind power forecasting models of AWPPS are integrated in the More-Care Energy Management System and installed for on-line operation in Ireland and other sites such as Crete, Madeira, etc.

VII. ACKNOWLEDGMENT

Fig. 1 .

 1 Fig. 1. The error margin as a function of the wind park power curve. The observed prediction error itself is in general the result of three factors; a modelling error e mod , an error due to the accuracy of the input meteorological predictions e NWP and finally, a stochastic component linked to the process itself e s : ( ) s NWP t k t e e e f e , , mod / = +

Fig. 3 .

 3 Fig. 3. Splitting the power curve into three power class fuzzy sets and into two cut-off risk zones.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Dips below 0.6 of the cross-correlation show situations where poor accuracy of the wind power predictions should be expected.

Fig. 7 .

 7 Fig. 7. Stable (left picture) and unstable (right picture) weather situations.

Fig. 10 .

 10 Fig. 10. Cumulative distribution functions of the prediction errordepending on the MRI-index bin.

Fig. 11 .

 11 Fig. 11. Wind power prediction with 85% confidence intervals.

Fig. 12 .Fig. 13 .

 1213 Fig. 12. Wind power prediction with the resampling and finetuned confidence intervals. The intervals are narrowed for the first 24 hours due to a low MRI-index value.

Fig. 14 .

 14 Fig. 14. Average width of the 6-hour ahead confidence intervals for WF-A and for various specified confidences (75, 80, 85, 90 & 95%).

Fig. 15 .

 15 Fig. 15. Parallel evolution of the MRI-index and of the normalized prediction error.

Table II .

 II Observed confidence for the two types of intervals after one year of simulation and effects of the MRI-index on the intervalreduction.

This work was performed in the frame of the ANEMOS Project (ENK-CT2002-00665) funded in part by the European Commission. The authors gratefully acknowledge ESB National Grid for providing data for the realization of the study.