Karën Fort
email: karen.fort@sorbonne-universite.fr

Bruno Guillaume
email: bruno.guillaume@inria.fr

Mathieu Constant
email: mathieu.constant@univ-lorraine.fr

Nicolas Lefèbvre
email: nicolas.lefebvre@inria.fr

Yann-Alan Pilatte
email: yann-alan.pilatte@etu.sorbonne-universite.fr

Matthieu Constant

Identification of Multi-Word Expressions Using a Slightly Gamified Crowdsourcing Platform

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

Introduction and State of the Art

The identification of multi-word expressions (MWEs) is crucial in natural language processing (NLP) [START_REF] Constant | Multiword expression processing: A survey[END_REF]. Significant efforts have been made in recent years on the subject, in particular though the PARSEME international network [START_REF] Savary | PARSEME -PARSing and Multiword Expressions within a European multilingual network[END_REF]. However, although some collective expert-based annotation initiatives have been successfully undertaken [START_REF] Schneider | SemEval-2016 Task 10: Detecting Minimal Semantic Units and their Meanings (DiMSUM)[END_REF][START_REF] Savary | The parseme shared task on automatic identification of verbal multiword expressions[END_REF], language resources are still limited in coverage and the need remains to identify newly-created MWEs. One potential solution is to exploit the so-called "wisdom of the crowd".

There have been several research papers on the interpretation of MWEs by native speakers, in particular by Gibbs [START_REF] Gibbs | What do idioms really mean?[END_REF][START_REF] Gibbs | Metaphor in idiom comprehension[END_REF]. More recently, [START_REF] Ramisch | How naked is the naked truth? a multilingual lexicon of nominal compound compositionality[END_REF] involved microworking crowdsourcing on Amazon Mechanical Turk. Finally, the experiment described in [START_REF] Krstev | Games on multiword expressions for community building[END_REF] involves a gamified interface allowing MWE researchers to guess the meaning of opaque MWEs in other languages.

However, we could find no publication concerned with evaluating human ability to identify MWEs in a text, without taking their interpretation into account.

On the other hand, voluntary crowdsourcing, especially in the form of Games with a Purpose (GWAPs), has proven effective in terms of both the quantity and quality of the data produced. Successful examples of such platforms include JeuxDeMots (Lafourcade, 2007), Phrase Detectives [START_REF] Poesio | Phrase detectives: Utilizing collective intelligence for internet-scale language resource creation[END_REF]ZombiLingo (Guillaume et al., 2016).

We created a gamified platform named RigorMortis1 (see Figure 1), the first of its kind, for MWEs annotation in French2 . This platform includes a task enabling evaluation of the participants' intuition concerning MWEs, the results of which we present here.

For non-francophones, the phrase "Fingers in the nose" refers to the French idiom "Les doigts dans le nez" which means "without any difficulty", "with both hands behind one's back", etc. 2 Description of the Experiment

Reference corpus

As the participants are volunteers, we had to keep the intuition task short. We therefore created a rather small reference corpus. It is composed of ten sentences taken from articles on French political scandals from the French Wikipédia. Although their number is small, they were carefully selected to include MWEs corresponding to distinct identification criteria. The corpus was annotated and adjudicated in MWEs by experts. It contains 16 MWEs. One sentence, however contains no MWEs (see Table 1).

This reference corpus has been built following precise annotation guidelines inspired by those of the PARSEME shared task on verbal MWEs [START_REF] Savary | The parseme shared task on automatic identification of verbal multiword expressions[END_REF], which includes a French dataset [START_REF] Candito | Annotation d'expressions polylexicales verbales en français[END_REF]. The criteria used for identifying MWEs, and in particular for detecting their morphosyntactic, syntactic and semantic idiosyncrasies are purely formal: for instance, no possible lexical substitution of a component by a synonym, presence of a "cranberry" word, no possible insertion of plausible material. We therefore discard semantically and syntactically compositional expressions that display statistical idiosyncrasy: for instance, institutionalized phrases in the sense of [START_REF] Sag | Multiword Expressions: A Pain in the Neck for NLP[END_REF], like traffic light. Further, only fixed lexical components of the expressions are annotated, so final prepositions in MWEs that can be considered part of the MWE valency frame (like simple verbs) are not annotated. For instance, only en raison (because) is annotated in en raison de (because of) where de (of) is a preposition. Such subtlety is unknown and not natural for participants. This is why we took both variants into account when evaluating the task.

The reference MWEs can be divided into two subtypes: MWEs that behave as function words (later called functional MWEs) and MWEs that do not. In French, functional MWEs are mainly fixed in the sense of [START_REF] Sag | Multiword Expressions: A Pain in the Neck for NLP[END_REF]. Such fully lexicalized expressions are immutable: they can undergo neither morphosyntactic nor syntactic variations, and insertion of plausible material is impossible. Notice that not all fixed expressions are necessarily functional MWEs: dommages et intérêts (damages), for example, functions as a noun. Furthermore, some of the functional MWEs that we consider are not entirely fixed: for instance, aux yeux des enfants (in the eyes of the children) = à leurs yeux (in their eyes). In our corpus, 7 MWEs are functional and 9 MWEs are not.

Our notion of functional MWEs can also be related to the category of fixed MWEs defined in the Universal Dependencies [START_REF] Nivre | Universal dependencies v1: A multilingual treebank collection[END_REF] where the guidelines state that: "[A fixed MWE] is used for certain fixed grammaticized expressions that behave like function words or short adverbials."

Crowdsourcing platform

We ask the participants to find "expressions multi-mots" (multi-word expressions) and give them a couple of examples, explaining that MWEs are non-compositional. We also mention that more "functional" expressions can also be MWEs and should be annotated. The participants are then directly asked to identify the MWEs in the sentences we propose, without any prior training.

The interface, inspired by that of TileAttack [START_REF] Madge | Experiment-driven development of a gwap for marking segments in text[END_REF], allows users to annotate multiple and discontinuous MWEs (see Figure 2). It should be noted that the intuition task is part of a larger gamified platform and that, while the participants did not gain points in this phase, they did in the other phases (see Figure 1). It is important that, during this phase, we give no feedback to the participants on their annotations. They can see feedback concerning their results only once they are finished annotating all ten sentences (see Figure 3), so that they are not biased.

The crowdsourcing interface was publicized mainly on social networks and natural language processing lists. 3 Results Obtained

Global Results

The ten reference sentences were played by 65 to 68 players (66.8 participants on average)3 .

Table 2 shows the precision, recall and F-measure for the whole annotation with two settings:

• Perfect match: the player selected exactly the same tokens as in the reference;

• Approximate match: the difference between the selection of the player and the reference includes only function words; for instance, [président+République] is accepted for [prési-dent+de+la+République], as de (of) is a preposition and la (the) is a determiner. It should be noted that 40 out of 68 participants (i.e. 58.82% ) correctly annotated the sentence which contains no MWE.

These global results show that when the participants identify an MWE, they are often right (in 58.22% of the cases), but that they are less good at finding them all (less than 50% were found).

Results for Individual MWEs: the Impact of Functional MWEs

To determine if an MWE is more or less easy to find, we computed the recall for each MWE separately (see Figure 4).

Again, we show the two values: i) perfect (dark blue) and ii) approximate (light blue) match. We observe in Figure 4 that there is a significant difference between the subset of functional MWEs (on the left) and that of none-functional MWEs (on the right).

Table 3 gives the recall value for these two subsets. The last column gives the overall value (already given above) for comparison. The recall for non functional MWEs reaches 65.05, which is more than twice that of functional ones (30.41).

In our own experience (and that of some participants), this difference in the identification of functional MWEs arises because we are so accustomed to them (they are so familiar) that we simply do not "see" them and forget to annotate them.

These results are encouraging, as they show that the participants can be rather efficient at identifying at least some types of MWE.

Analysis of the Noise Produced

Table 4 lists the ten expressions that were identified by more than 10% of the participants and which are not in the reference. Note that we reserve the term collocation to refer to any statistically significant co-occurrence, including all forms of MWEs as described above and compositional phrases which are predictably frequent It shows that, unsurprisingly, the participants had difficulty distinguishing between MWEs and compositional expressions exhibiting statistical idiosyncrasy (in six cases). This can be explained by the fact that our definition of MWEs partially overlaps with the notion of collocation, as defined in [START_REF] Sag | Multiword Expressions: A Pain in the Neck for NLP[END_REF]. Collocations refer to "any statistically significant co-occurrence", that includes both syntactically/semantically compositional and non-compositional expressions.

Other mistakes include boundary errors (two cases) and common civilities annotated as MWEs (two cases). These could probably be avoided if the participants were properly trained.

Conclusion and perspectives

Although it was carried out on a small corpus4 , this experiment gathered results from a satisfying number of participants and showed that volunteers with no prior training can help identify at least some MWEs in texts. It is encouraging that the most difficult MWEs to find are the functional ones, as these are usually the first to be listed and are the least prone to neologism.

However, while the participants' intuition proves valuable, it should be complemented by proper training, using at least some of the tests defined by the PARSEME network.

Figure 1 :

 1 Figure 1: Interface for game selection.

Figure 2 :

 2 Figure 2: Interface for the annotation of the MWEs.

Figure 3 :

 3 Figure 3: Feedback given to the participant after the intuition phase.

Table 1 :

 1 Glosses and translations of the MWEs of the experiment.

	MWEs	Glosses	Translations
	se voiler la face	to cover one's face	to bury one's head in the sand
	file d'attente	queue of waiting	waiting line
	dommages et intérêts	damages and interests	damages
	mode de vie	way of life	way of life
	le président de la République the President of the Republic the President of the Republic
	ministre des finances	Minister of finances	Secretary of the Treasury
	extrême droite	extreme right	far right (political)
	chef d'État	chief of state	head of the State
	mettre aux voix	put to the voices	put to the vote
	entre autre	between other	among others
	au-delà	further	beyond
	en raison (de)	in reason (of)	because (of)
	peut-être	may-be	maybe
	aux yeux (de X)	at the eyes of X	in X's view
	d'ailleurs	from elsewhere	by the way
	dans le but (de)	in the goal (of)	in order to

Table 2 :

 2 Global results.

		Precision Recall F-measure
	Perfect	48.13	41.22	44.41
	Approximate	58.22	49.86	53.72

Table 3 :

 3 Recall for functional and non-functional MWEs.

	90%
						Perfect	Approx.
	70%
	50%
	30%
	10%
	m i n i s t r e d e s fi n a n c e s e s p l e v d d m o e e i ( e ) r é i d n t ( d e l a R t s ) s g m o d m a e e t i n t é ê r ( l ' ) e x t r ê m e d r o i t e ( l e ) c h e f ( d ' ) É t a t ' ( e f l r i v e ( s ) o l e ( a ) a c ( l a ) fi l e d ) a t t e n t e m e t t r e ( a u x ) v o i x	e n t r e a u t r e	a u -d e l à ( d u ) e n r a i s o n ( d e ) p e u t ( -) ê t r e	D ' a i l l e u r s a u x y e u x ( d u )	d a n s l e b u t ( d e )
		é
		p u b l i q u e
		non-functional functional	All
	Perfect	58.19	19.48		41.22
	Approximate	65.05	30.41		49.86

Figure 4: Recall of the participants on each MWE (left: non-functional; right: functional).

See: rigor-mortis.org

We believe it is adaptable to any language.

We removed from this analysis four participants who are experts in the MWE subdomain.However, we kept the participants from the more general NLP and linguistics domains.

A way to increase the size of the corpus without making the task longer could be to randomize the sentences proposed to the participants.

A complementary experiment on the subject is work in progress, as the platform also enables researchers to train participants and collect annotations from them.

The reference corpus is freely available under a CC BY-SA license, on the platform itself.


Simon J Berrebi

Etienne Hans

Nicolas Chiabaut

Jorge A Laval

Ludovic Leclercq

Kari E Watkins

Comparing Bus Holding Methods With and Without Real-Time Predictions

On high-frequency routes, transit agencies hold buses at control points and seek to dispatch them with even headways to avoid bus bunching. This paper compares holding methods used in practice and recommended in the literature using simulated and historical data from Tri-Met route 72 in Portland, Oregon. We evaluated the performance of each holding method in terms of headway instability and mean holding time. We tested the sensitivity of holding methods to their

parameterization and to the number of control points. We found that Schedule-Based methods require little holding time but are unable to stabilize headways even when applied at a high control point density. The Headway-Based methods are able to successfully control headways but they require long holding times.

Prediction-Based methods achieve the best compromise between headway regularity and holding time on a wide range of desired trade-o↵s. Finally, we found the prediction-based methods to be sensitive to prediction accuracy, but using an existing prediction method we were able to minimize this sensitivity.

These results can be used to inform the decision of transit agencies to implement holding methods on routes similar to TriMet 72.

Introduction

On high frequency routes, there is a natural tendency for buses to bunch together. When a bus is traveling with a long headway1 , it has to pick up and drop o↵ a relatively greater number of passengers, which slows down the bus even more. As the lagging bus becomes crowded, the following buses only have a few passengers to serve, which they can do relatively fast. Eventually the lead bus may get caught by one or several following buses and they start traveling as a platoon. Bus bunching is the product of unstable dynamics that cause delays to grow [START_REF] Hickman | An Analytic Stochastic Model for the Transit Vehicle Holding Problem[END_REF]. Even a small perturbation such as a tra c signal or a passenger paying in cash can destabilize the route and lead to bus bunching 10 [START_REF] Kittelson | Transit capacity and quality of service manual[END_REF][START_REF] Milkovits | Simulating Service Reliability of a High Frequency Bus Route Using Automatically Collected Data by Master of Science in Transportation[END_REF].

Unstable headway dynamics are a systemic problem that causes passenger wait and crowding. [START_REF] Fan | Do transit users just wait for buses or wait with strategies?[END_REF] showed that on routes where headways are less than 12 minutes, passengers tend to arrive randomly, even if a schedule is available. Because more passengers arrive during long headways than during short ones, gaps in service cause disutility to passengers in the form of undue waiting time and crowding [START_REF] Newell | Maintaining a bus schedule[END_REF][START_REF] Milkovits | Simulating Service Reliability of a High Frequency Bus Route Using Automatically Collected Data by Master of Science in Transportation[END_REF]. One way for transit agencies to stop the progression of instability among headways, is to provide control points, where buses with short headways can be held to absorb the delay of following buses.

20

Holding buses at control points can help reduce at-stop passenger waiting time, but it increases the wait of passengers who have already boarded. There is a trade-o↵ between stabilizing headways and maintaining high operating speed [START_REF] Furth | TCRP Report 113: Using Archived AVL-APC Data to Improve Transit Performance and Management[END_REF][START_REF] Furth | Setting frequencies on bus routes: theory and practice[END_REF][START_REF] Cats | Impacts of holding control strategies on transit performance: Bus simulation model analysis[END_REF]. This is why transit agencies value the benefit of headway reliability and the disadvantage of holding time di↵erently.

In addition to selecting a holding method for their routes, transit agencies also need to decide how to implement it. Several holding methods in the literature require setting a parameter, which a↵ects the trade-o↵ between holding time and headway stability [START_REF] Daganzo | A headway-based approach to eliminate bus bunching: Systematic analysis and comparisons[END_REF]Xuan et al., 2011;[START_REF] Bartholdi | A self-coordinating bus route to resist bus bunching[END_REF][START_REF] Daganzo | Reducing bunching with bus-to-bus cooperation[END_REF]. Holding methods can also be applied at one or several control points along the route, which may impact the performance of each method. Understanding sensitivity of holding methods on the parameterization and number of control point is necessary to select the most adequate holding method based on route characteristics and desired trade-o↵s Several methods are based on predictions for the arrival times of following buses [START_REF] Bartholdi | A self-coordinating bus route to resist bus bunching[END_REF][START_REF] Daganzo | Reducing bunching with bus-to-bus cooperation[END_REF][START_REF] Berrebi | A real-time bus dispatching policy to minimize passenger wait on a high frequency route[END_REF]. The quality of the predictions may a↵ect transit operators' ability to leverage headway stability from holding time. The required level of prediction accuracy and confidence can be burdensome for certain transit operators that may not need high-quality predictions for other applications. The ubiquity of prediction-based methods is therefore dependent on their sensitivity to prediction quality.

Research in the literature has compared holding methods, but there currently lacks a unified framework to evaluate the conflicting objectives of stabilizing headways and minimizing holding time. Xuan et al. (2011) and [START_REF] Berrebi | A real-time bus dispatching policy to minimize passenger wait on a high frequency route[END_REF] have case study sections to compare methods in the literature to their own. [START_REF] Cats | Impacts of holding control strategies on transit performance: Bus simulation model analysis[END_REF] compare naive methods used in practice and a headwaybased method similar to the method in [START_REF] Daganzo | Reducing bunching with bus-to-bus cooperation[END_REF]. There is a need for a sensitivity analysis to support the choice of holding methods and their parameterization based on route characteristics, including the number of control points on routes similar to Tri-Met 72.

In this paper, we investigate the holding trade-o↵ of holding methods used in practice and recommended in the literature. To this end, we evaluate holding methods on a simulated bus route using historical data from Tri-Met Route 72 in Portland, Oregon. We use the prediction tool developed in [START_REF] Hans | Real-time bus route state 720 forecasting using particle filter[END_REF] to reproduce the predictions in a realistic setting. In the following section, we describe the holding methods used in practice and recommended in the literature. In Section 3, we discuss the simulation experiment, and particularly the methods evaluated. In Section 4 we compare the performance of each holding method. In Section 5, we investigate the impact of parameter choice, and number of control points on the trade-o↵ between stabilizing headways and keeping short holding times. In Section 6 we test the sensitivity of prediction-based holding methods on the accuracy and confidence of predictions. Finally, we provide concluding remarks in Section 7.

Holding methods in the literature

Methods to hold buses at control points have been addressed for many decades. [START_REF] Osuana | Control Strategies for an Idealized Public Transportation System[END_REF] and [START_REF] Newell | Control of Pairing Vehicles on a public Transportation Route, Two Vehicles, One Control Point[END_REF] formulated the theoretical basis for holding mechanisms to minimize passenger waiting time on simple routes in the 1970's. Since then, two main approaches to the bus holding problem have been developed in the literature, mathematical optimization and analytical.

The first approach consists in optimizing a weighted function of passenger wait in mathematical programs that consider the dynamics of bus trajectories (analytically or by simulation). At each decision stage, the optimization tools model the future states of the system, and assign holds on a rolling horizon. [START_REF] Hickman | An Analytic Stochastic Model for the Transit Vehicle Holding Problem[END_REF] developed a linear search optimization algorithm which considers holding decisions in isolation of each other based on a stochastic model for bus trajectories. In [START_REF] Eberlein | The holding problem with real-time information available[END_REF] a heuristic algorithm is used to minimize the waiting time of passengers at stops in a quadratic program. [START_REF] Bukkapatnam | Distributed architecture for real-time coordination of bus holding in transit networks[END_REF] developed an iterative model where buses and stations negotiate holding time to minimize marginal costs. The method in Zolfaghari et al. (2004) assigns all holding decisions simultaneously, while considering capacity constraints, using AVL data and perfect predictions. [START_REF] Delgado | Real-Time Control of Buses in a Transit Corridor Based on Vehicle Holding and Boarding Limits[END_REF] and [START_REF] Delgado | How much can holding and/or limiting boarding improve transit performance?[END_REF] developed a simulation-based optimization algorithm that reproduces stationary bus trajectories deterministically and minimizes a weighted function of wait time. [START_REF] Sánchez-Martínez | Real-time holding control for high-frequency transit with dynamics[END_REF], 2015) extended their methods to consider dynamic passenger arrival rates and travel time. [START_REF] Cortés | Hybrid predictive control for real-time optimization of public transport systems' operations based on evolutionary multi-objective optimization[END_REF] used a genetic algorithm to solve a multi-objective dynamic problem.

90

The second approach assigns holds as closed-form functions of bus arrival times [START_REF] Daganzo | A headway-based approach to eliminate bus bunching: Systematic analysis and comparisons[END_REF][START_REF] Daganzo | Reducing bunching with bus-to-bus cooperation[END_REF]Xuan et al., 2011;[START_REF] Bartholdi | A self-coordinating bus route to resist bus bunching[END_REF][START_REF] Berrebi | A real-time bus dispatching policy to minimize passenger wait on a high frequency route[END_REF]. Buses are held with the objective of maintaining stable headways, and preventing bus bunching from the onset, which can minimize passenger waiting time globally and durably. Methods assign holds to buses as a function of the schedule, headways and, for some, predicted arrival times 2 . Unlike the mathematical programming approach, these methods do not consider the prediction model further than its output. Therefore, they can use any prediction model 3 , which makes them much easier to implement. Closed-form holding methods will be the focus of this paper.

100

The notation used in this paper is consistent with [START_REF] Berrebi | A real-time bus dispatching policy to minimize passenger wait on a high frequency route[END_REF], and is shown in Table 1. The arrival time of the i th bus at the control point is A i (random variable) for a future arrival, and a i (realization of A i ) for a known arrival time. Once the i th bus arrives, it holds for time h i , and is dispatched at time d i = a i + h i . If the route runs according to a schedule, di is the scheduled departure time from the control point, and H i is its scheduled headway (H i = di di 1 ). The holding methods evaluated in this paper are described in their Eulerian version (as in Xuan et al. (2011)) in Table 2 with the information they require and their recommended holding times.

2 In the remainder of this text, we refer to holding methods that consider schedules as "schedule-based" and methods that consider headways as their main input as "headwaybased". Schedule-based methods include the Naive Schedule and the method recommended in Xuan et al. (2011). Headway based methods include the Naive Headway and the method recommended in [START_REF] Daganzo | A headway-based approach to eliminate bus bunching: Systematic analysis and comparisons[END_REF].

3 A discussion follows in Section 6. The simplest and most widely used methods to hold buses at control points are to plan scheduled departure times, di , or scheduled headways, H i , well in advance [START_REF] Boyle | REPORT 135 Controlling System Costs : Basic and Advanced Scheduling Manuals and Contemporary Issues in Transit Scheduling[END_REF][START_REF] Abkowitz | Implementing headway-based reliability control on transit routes[END_REF]Van Oort et al., 2010).

In this paper, we will consider the Naive Schedule and Naive Headway methods as holding buses until their departure time or headway reaches the planned threshold, as shown in Equations 1 and 2 of

(ai di 1 ) 1+E " arg max r i Ar a i r i ! 1 # (7)

Prediction-based holding methods

A novel closed-form approach to the bus holding problem is to dispatch 150 buses according to the predicted arrival times of following buses at the control point. Real-time prediction methods are becoming increasingly accurate and available, which allows the replacement of planned operations by the natural headway in current operating conditions. Today, the vast majority of public transportation agencies in the United-States are capable of tracking their vehicles in real-time [START_REF] Grisby | APTA Surveys Transit Agencies on Providing Information and RealTime Arrivals to Customers[END_REF]. Using real-time vehicle locations, increasingly sophisticated prediction algorithms have surfaced recently [START_REF] Chien | Dynamic bus arrival time prediction with artificial neural networks[END_REF][START_REF] Cathey | A prescription for transit arrival/departure prediction using automatic vehicle location data[END_REF]Jeong andRilett, 2004, 2005;[START_REF] Chen | Dynamic Prediction Method with Schedule Recovery Impact for Bus Arrival Time[END_REF][START_REF] Shalaby | Prediction model of bus arrival and departure times using avl and apc data[END_REF][START_REF] Gurmu | Dynamic travel time prediction models for buses using only gps data[END_REF]Sun et al., 2007;[START_REF] Mazloumi | Prediction intervals to account for uncertainties in neural network predictions: Methodology and application in bus travel time prediction[END_REF]. Most notably, [START_REF] Hans | Investigating the irregularities of bus routes[END_REF][START_REF] Hans | Real-time bus route state 720 forecasting using particle filter[END_REF] developed a prediction method 160 specifically for the purpose of real-time control. Their prediction algorithm can generate probability distributions of arrival times. Based on these predictions, holding methods can consider following buses to equalize headways, without having to rely on planned operations.

Using real-time predictions, [START_REF] Bartholdi | A self-coordinating bus route to resist bus bunching[END_REF] developed a holding strategy that can stabilize headways without the need for planned operations. The method consists in holding each vehicle for the predicted time until the next arrival by a factor ↵ 2]0, 1[ as in Equation 5in Table 2. When several buses arrive in close succession, however, the method sends the middle few uncontrolled. To prevent this, Bartholdi and Eisenstein added a minimum forward headway 7 . Otherwise, the method can split the burden on a big gap between two buses: each vehicle leaves the control point with the weighted sum between its forward 8 and backward headways. Unlike the methods cited thus far, their method involves a mechanism that acts locally to scale headways to the rate of arriving vehicles.

Daganzo and Pilachowski used predictions on the next arrival time to blend the forward headway with the backward headway, as shown in Equation 6of Table 2 [START_REF] Daganzo | Reducing bunching with bus-to-bus cooperation[END_REF]. The method considers the forward headway in the same way as [START_REF] Daganzo | A headway-based approach to eliminate bus bunching: Systematic analysis and comparisons[END_REF], but it also subtracts the deviation of the expected time until the next arrival from the scheduled headway, H

(E[A i+1 ] a i )
, by a factor ↵ 2]0, 1 2 [. The holding method can reduce the di↵erence between the forward and backward headways by holding buses for a weighted sum of that di↵erence.

More recently, Berrebi and collaborators developed a method that takes a global approach by considering every bus on the route. The method is a generalization of [START_REF] Daganzo | Reducing bunching with bus-to-bus cooperation[END_REF], without H and , that considers n buses. The vehicle that will arrive with the maximum relative delay, max r i Ar ai r i , is probabilistically identified and each preceding vehicle is held to absorb a share of that delay, as in Equation 7in Table 2. When the lagging bus arrives at the control point, it can be dispatched with approximately the same headway as the leading few. The method can di↵use big gaps organically without the need 7 To let the method act primarily on backward headways, we set the minimal forward headway as H i /2.

8 More exactly the backward headway of the following bus, a headway ago.

for schedules, scheduled headways, or any kind of explicit slack time. Oregon. The scheduled headway alternates between seven and eight minutes in the afternoon peak. Historically, Route 72 had a bus bunching problem during peak hours [START_REF] Berkow | Transit Performance Measurement and Arterial Travel Time Estimation Using Archived AVL Data[END_REF]. Buses on route 72, however, are equipped 200 with a Computer Aided Dispatching (CAD) system that would allow replacing the current schedule with real-time control. To evaluate the potential benefits and disadvantages of applying each real-time holding method described in the previous section, we have tested their performance in a case study. In the study, we used historical data leading up to the first control point 10 .

Case study

Each method had the opportunity to hold vehicles at that point to stabilize operations and mitigate bus bunching. Whenever the boarding and alighting times exceeded the hold recommended by a control method, buses were dispatched after they finish loading and unloading. The headways of buses leaving the control point were recorded to evaluate the performance of the dispatching strategy. We used simulated data downstream from the control point to take into account the impacts of each holding method on headway dynamics.

In this paper, we considered a single value of H i : its historical value. In the historical data, vehicles arrive at the control point(s) according to a set scheduled frequency of service. The scheduled frequency at a control point determines the scheduled headway, because the frequency of departure cannot be greater than the frequency of arrival. If we reduced the scheduled headway, the holding methods would be unable to control buses at all. If we increased the scheduled headway, a continually accumulating queue of buses would form.

This is why we did not test holding methods with varying values of H i .

Prediction

Once buses arrive at the control point, each method can recommend a holding time based on the last departure time (recommended by the method) and the predicted next arrivals. The methods recommended in [START_REF] Daganzo | Reducing bunching with bus-to-bus cooperation[END_REF] and [START_REF] Bartholdi | A self-coordinating bus route to resist bus bunching[END_REF] [START_REF] Berrebi | A real-time bus dispatching policy to minimize passenger wait on a high frequency route[END_REF] needs the probability distribution of each future arrival to infer the maximum relative delay, max Ar ai r i , and its corresponding index, arg max Ar ai r i . The prediction tool used in our simulation is a particle filter combined with an event-based mesoscopic model developed in [START_REF] Hans | Investigating the irregularities of bus routes[END_REF] and [START_REF] Hans | Real-time bus route state 720 forecasting using particle filter[END_REF]. The prediction tool is capable of generating simulated trajectories solely using vehicle location data. The predicted arrival times of buses at the control point are generated iteratively as a function of their dwell and travel times, as 240 shown in Table 3. The dwell time is estimated as the sum of boarding, alighting and door operation time, assuming passengers board and alight through the same door. The number of passengers boarding are assumed to follow a Poisson distribution, with no capacity constraints. The share of passengers alighting follows a binomial distribution with passenger loads estimated using historical headways. Travel times between stations are generated as a Gamma distribution. To ensure that the particle filter accurately reproduces operations on Route 72, the parameters of dwell and travel times, such as the rate of arriving passengers and the mean travel time between stops were calibrated on historical data from the route. Odd days were used for the calibration, and even days were 250 kept for the simulation. The particle filter recommended in [START_REF] Hans | Investigating the irregularities of bus routes[END_REF] and [START_REF] Hans | Real-time bus route state 720 forecasting using particle filter[END_REF] is particularly suitable for real-time control. Unlike regression and machine learning prediction methods, the particle filter is capable of generating probability distributions, which is necessary for the method in [START_REF] Berrebi | A real-time bus dispatching policy to minimize passenger wait on a high frequency route[END_REF]. [START_REF] Newell | Maintaining a bus schedule[END_REF], as shown in Equation 9.

i,s = (1 µ s )L i,s 1 + d s h i,s L i,0 = 0 ⇡ i,s ⇠ Gam(
CV 2 [headway] = V ar[headway] E[headway] 2 (8) APW = E[headway] 1 + CV 2 [headway] 2 (9)
Variation in headways tends to increase along the route unless buses are controlled [START_REF] Hickman | An Analytic Stochastic Model for the Transit Vehicle Holding Problem[END_REF]. Therefore, unstable headway dynamics can have 280 lasting e↵ects on the system and on passenger experience if headways are not stabilized. We choose to evaluate headway stability in terms of CV 2 because it is a dimensionless parameter that directly relates to at-stop passenger wait and allows extrapolation of results to other routes.

Holding methods can stabilize headways by trading o↵ holding time. Holding time too has a cost. It causes passengers who are already aboard the vehicle undue waiting time. In addition, holding time reduces bus operating speed.

Finally, holding at a control point can disrupt surrounding tra c, depending on its location. Since transit agencies may value the benefit of headway stability and the cost of holding time di↵erently depending on the route, we chose to 290 keep these measures of performance separate.

Cross-comparison

We simulated the decision process of each holding method described in Table 2 in the afternoon peak hour, when scheduled headways, H i , oscillate between seven and eight minutes. In the simulation, buses were held at the 48 th station, which is seven miles away from the departure point, at the intersection with the light-rail, MAX, on Banfield Expressway. We chose this station because it is used as a schedule recovery point in historical data. The station also has by far the greatest alighting proportion and boarding demand on Route 72.

These considerations are important because the overall cost of holding time is 300 proportional to the number of passengers who ride through the holding point and the overall benefit of stabilizing headways is proportional to the number of passengers waiting at downstream stops. Therefore, holding at station 48 inconveniences few passengers and benefits many.

Every method was parameterized with calculated at each arrival as per [START_REF] Daganzo | A headway-based approach to eliminate bus bunching: Systematic analysis and comparisons[END_REF], [START_REF] Daganzo | A headway-based approach to eliminate bus bunching: Systematic analysis and comparisons[END_REF][START_REF] Daganzo | Reducing bunching with bus-to-bus cooperation[END_REF]Xuan et al. (2011), and ↵ = 1 2 to provide a middle-ground basis of comparison with the Naive Schedule and Headway methods. We discuss the choice of the ↵ parameter in the next section. In historical data, buses start the route with close to even headways, but get destabilized along the route, leading to the control point. The CV 2 increases in a saw tooth pattern, with small drops corresponding to mid-route control points. When buses arrive at the main control point, CV 2 is 42%, which, according to the second edition of the Transit Capacity and Quality of Service, corresponds to a Level of Service (LOS) E, and is denoted as Frequent bunching [START_REF] Kittelson | Transit capacity and quality of service manual[END_REF]. In historical data, CV 2 is reduced to 30%, which is still LOS E. We did not display the mean lost time in historical data because we could not di↵erentiate the share that was conscious holding from boarding, alighting, and other dwell operations.

Station No

The holding method applied in historical data, and the method recom- The Naive Headway method is capable of sending buses at regular intervals, but requires long holding times. The Naive Headway method was e↵ective at stabilizing headways with 7% of CV 2 (LOS B, vehicles slightly o↵ headway).

The method, however, holds buses for 255 seconds on average. The holding times tend to accumulate because each late bus pushes back the dispatching time of all upstream buses. The method proposed in [START_REF] Daganzo | A headway-based approach to eliminate bus bunching: Systematic analysis and comparisons[END_REF] dispatches buses with much more erratic headways than the Naive Headway, which causes longer wait for passengers at a stop, but the method keeps mean holding time at 40 seconds.

The prediction-based holding methods consider the predicted arrival times of following buses to di↵use big gaps. The method in [START_REF] Bartholdi | A self-coordinating bus route to resist bus bunching[END_REF] dispatched buses at CV 2 = 0.13 (LOS C, vehicles often o↵ headway).

The method, however, holds buses for 180 seconds on average. The method in [START_REF] Daganzo | Reducing bunching with bus-to-bus cooperation[END_REF] produced almost the same results as in [START_REF] Bartholdi | A self-coordinating bus route to resist bus bunching[END_REF], but only required 98 seconds of holding time.

The method recommended in [START_REF] Berrebi | A real-time bus dispatching policy to minimize passenger wait on a high frequency route[END_REF] considers the arrival time of every following bus on the route. The holding method reduces CV 2 to 7% (LOS B, Vehicles slightly o↵ headway) with 160 seconds of mean holding time.

350

The reduction in headway variability at the control point has lasting e↵ects downstream from the control point. The more a method is able to reduce headway variability at the control point, the less headways tend to destabilize further down the route. Conversely, the rate of increase of CV 2 is highest for methods that are the least able to stabilize headways. For example, between stops 49 and 60, CV 2 increased three times more for the method in Xuan et al.

(2011) than the method in [START_REF] Berrebi | A real-time bus dispatching policy to minimize passenger wait on a high frequency route[END_REF], which dispatched buses with far more stable headways. This trend is not seen in historical data because it benefits from mid-route control points, which we have omitted in the simulation.

Sensitivity 360

The performance of control strategies can vary depending on how they are applied. It is important to understand the implications of parameter choice and the number of control points. In this section, we investigate how the performance of each method is a↵ected by these factors.

Parameterization

The methods recommended in [START_REF] Daganzo | A headway-based approach to eliminate bus bunching: Systematic analysis and comparisons[END_REF], [START_REF] Daganzo | Reducing bunching with bus-to-bus cooperation[END_REF], [START_REF] Bartholdi | A self-coordinating bus route to resist bus bunching[END_REF][START_REF] Bartholdi | A self-coordinating bus route to resist bus bunching[END_REF]Xuan et al. (2011) all require setting an ↵ parameter, but their interpretation of the parameter is di↵erent. For the methods prescribed in [START_REF] Daganzo | A headway-based approach to eliminate bus bunching: Systematic analysis and comparisons[END_REF] and Xuan et al. (2011), ↵ corresponds to the fraction of the deviation from headway or schedule that is to be 370 recovered by the hold. When ↵ is close to zero, buses are barely controlled, and for values ↵ close to one, the holding methods are similar to the Naive Schedule and Headway. For the method recommended in [START_REF] Bartholdi | A self-coordinating bus route to resist bus bunching[END_REF], ↵ is the fraction of the predicted backward headway that buses should hold for. The mean holding time is ↵ times the average headway. In [START_REF] Daganzo | Reducing bunching with bus-to-bus cooperation[END_REF] the ↵ parameter weights the di↵erence between the forward and backward headway. For values of ↵ close to zero, more importance is given to the forward headway, and for values close to one, more importance is given to the backward headway. than its recommended range, [START_REF] Daganzo | A headway-based approach to eliminate bus bunching: Systematic analysis and comparisons[END_REF], [START_REF] Daganzo | Reducing bunching with bus-to-bus cooperation[END_REF] and Xuan et al. ( 2011) compensate for perturbations excessively by dispatching buses past āi or H i . In the last section, we set ↵ = 1 2 for each parametric method. The solid dots on Figure 4 show the performance of each holding method as parameterized in the previous section.

13 Note that the interpretation of ↵ di↵ers for each method. 2011), the rate of decay remains high until the Naive Schedule trade-o↵ is reached, whereas the method in [START_REF] Daganzo | A headway-based approach to eliminate bus bunching: Systematic analysis and comparisons[END_REF] requires less than a third of Naive Headway's mean lost time.

Prediction-based methods achieve the best compromise between headway regularity and holding time in a wide range of settings. The method in [START_REF] Daganzo | Reducing bunching with bus-to-bus cooperation[END_REF] can be parameterized to yield a lower CV 2 than any other method for any holding time up to 130 seconds. The method in [START_REF] Berrebi | A real-time bus dispatching policy to minimize passenger wait on a high frequency route[END_REF] can dispatch buses with 7% of CV 2 and 160 seconds of holding time, making it the preferable method for holding times over 130 seconds.

Multiple Control Points

Whereas methods thus far have been applied at a single control point, we now test them at several control points along the route. Applying holds at several control points can help maintain stable headways throughout, but it requires more frequent holding. In this section, we evaluate the trade-o↵s between headway stability and mean holding time for each method based on the number of control points. This analysis can support the decision of transit agencies to implement the holding method most adequate for their route and objectives.

The methods in [START_REF] Berrebi | A real-time bus dispatching policy to minimize passenger wait on a high frequency route[END_REF] and [START_REF] Bartholdi | A self-coordinating bus route to resist bus bunching[END_REF] are designed to dispatch buses with one or few control points. The method in [START_REF] Berrebi | A real-time bus dispatching policy to minimize passenger wait on a high frequency route[END_REF] considers the predicted arrival times of each bus on the route. It cannot be applied at a close succession of control points because holds would interfere with predictions. The method in [START_REF] Bartholdi | A self-coordinating bus route to resist bus bunching[END_REF] holds vehicles for ↵H i on average. Holding time would therefore accumulate proportionally to the number of control points. Both methods, however, can be paired with on-route holding methods that only consider the local headway dynamics. We introduce two hybrid methods that apply the method in [START_REF] Bartholdi | A self-coordinating bus route to resist bus bunching[END_REF] (Hybrid #1) and [START_REF] Berrebi | A real-time bus dispatching policy to minimize passenger wait on a high frequency route[END_REF] (Hybrid #2) at a "main control point" (stop 29) and the method in [START_REF] Daganzo | Reducing bunching with bus-to-bus cooperation[END_REF] at all other control points.

The Hybrid Methods require the predicted arrivals of one or several upstream buses. In this simulation, the arrival times of vehicles at early stops were not predictable due to unavailable data. The Hybrid methods could only be applied on the second half of the route. In order to test non-Hybrid methods in the most favorable and realistic conditions, we allowed there methods to apply holding time progressively form the beginning of the route. Although the route characteristics including congestion levels, stop spacing, and passenger demand, are relatively stable throughout the route, the non-Hybrid methods may encounter slightly di↵erent operating conditions on the first 29 stops. This is why, Hybrid and non-Hybrid methods are plotted separately. The Schedule-Based Methods are similarly a↵ected by the number of control points. The Naive Schedule and the method in Xuan et al. ( 2011) are unable to stabilize the route when applied sparsely. When applied at many stops (16 or 32) the methods require between 152 and 188 seconds to maintain a CV 2 between 20% and 25%, which is more than the methods in [START_REF] Daganzo | Reducing bunching with bus-to-bus cooperation[END_REF] and [START_REF] Berrebi | A real-time bus dispatching policy to minimize passenger wait on a high frequency route[END_REF] for the same level of stability. As expected, the method in Xuan et al. ( 2011) yields more unstable headways but requires less holding time than the Naive Schedule method for any number of control points.

For any number of control points, the method applied in [START_REF] Daganzo | A headway-based approach to eliminate bus bunching: Systematic analysis and comparisons[END_REF] requires much less holding time than the Naive Headway, but it also dispatches vehicles with greater CV 2 . The method recommended in [START_REF] Daganzo | A headway-based approach to eliminate bus bunching: Systematic analysis and comparisons[END_REF] requires the same amount of holding time as the methods in [START_REF] Daganzo | Reducing bunching with bus-to-bus cooperation[END_REF] and [START_REF] Berrebi | A real-time bus dispatching policy to minimize passenger wait on a high frequency route[END_REF] but yields greater CV 2 . When applied at 1,2, or 4 control points, the Naive Headway method yields the same level of stability as methods in [START_REF] Daganzo | Reducing bunching with bus-to-bus cooperation[END_REF] and [START_REF] Berrebi | A real-time bus dispatching policy to minimize passenger wait on a high frequency route[END_REF] but requires far greater holding time. The Naive Headway can exceed the method in [START_REF] Berrebi | A real-time bus dispatching policy to minimize passenger wait on a high frequency route[END_REF] by 4%, when applied at 8 control point, at the cost of 266 additional seconds of holding time. The method recommended in [START_REF] Daganzo | Reducing bunching with bus-to-bus cooperation[END_REF] yields the lowest CV 2 of any non-Hybrid method for any value of holding times greater than 65 and less than 398 seconds.

In Part (b) of Figure 5, the Hybrid #2 consistently yields more stable headways with less holding time than the Hybrid #1 for any number of holding points. Although, the Hybrid methods were applied starting at stop 29, whereas non-Hybrid methods were applied throughout, it is worth noting that the Hybrid #2 method yields the lowest CV 2 out of all holding methods for any value of holding time greater than 118 and less than 398 seconds. These results are consistent with those obtained in section 4. We expect that these results would remain consistent in a live-implementation if the two sets of methods were tested on the same route portion due to the resemblance of the first 30 stop of the route with the following 35.

Figure 5 shows that for any value of holding time less than 177 seconds, the combination of method and number of holding points yielding the lowest CV 2 , (a) has the greatest number of holding points the method can support under the set holding time, and (b) has the lowest number of holding point of any method under the set holding time. In other words, the best compromise between headway stability and holding time is always the method with the most concentrated holding time. In particular, the Hybrid #2 yields the best compromise between the two objectives but concentrates holding time in few holding points located in the second half of the route. The concentration of holding time disproportionately impacts the passengers riding through the few holding points. Transit agencies should therefore be mindful of finding holding points where many passengers board and alight or selecting methods that dilute holding time, even at the cost of more holding and instability.

The number of control points a↵ects the nature of the trade-o↵ between holding time and headway stability. When applied at a single control point, the holding time is a riding cost imposed on passengers riding through. When applied at many control points, holding time can be approximated as a cost per distance. In any case, the at-stop waiting time due to uneven headways is a boarding cost for any passenger getting on the vehicle.

There is a discrepancy between Prediction-Based methods. The Hybrid # 1 yields both greater CV 2 and greater holding time than Partial Methods and other Prediction-Based methods for any number of control points.

and the Hybrid # 2, on the other hand, The Hybrid # 2 dispatches vehicles with more stable headways than the method in [START_REF] Daganzo | Reducing bunching with bus-to-bus cooperation[END_REF] for any number of control points and requires less holding time when applied at four or eight control points.

Prediction

For transit agencies that wish to take advantage of bus dispatching methods using real-time predictions such as [START_REF] Daganzo | Reducing bunching with bus-to-bus cooperation[END_REF], [START_REF] Bartholdi | A self-coordinating bus route to resist bus bunching[END_REF][START_REF] Bartholdi | A self-coordinating bus route to resist bus bunching[END_REF][START_REF] Berrebi | A real-time bus dispatching policy to minimize passenger wait on a high frequency route[END_REF], it is important to know what type of prediction is required and how accuracy will a↵ect their performance.

Although real-time vehicle tracking and prediction technologies are widely available among transit agencies, many of these systems were designed for passenger information rather than operational control. On these systems, using inadequate predictions could negatively impact the quality of dispatching mechanisms. Conversely, implementing a separate prediction system for control would duplicate e↵orts. To help transit agencies decide on the most appropriate prediction model, with respect to the performance of each holding method and to the costs of acquiring high-quality predictions, we tested the sensitivity of holding methods to the accuracy of predictions.

To evaluate the sensitivity of prediction-based holding methods to the prediction accuracy, we simulated their performance with synthetic predictions.

Each time a bus arrived at the control point, a synthetic distribution of arrival times of the following buses was generated with errors on the distribution mean, 1 , and within the distribution, 2 . The systematic error 1 a↵ects each prediction-based holding method because it biases the expected arrival time of each bus. The shape parameter on the other hand should not a↵ect the methods in [START_REF] Daganzo | Reducing bunching with bus-to-bus cooperation[END_REF] and [START_REF] Bartholdi | A self-coordinating bus route to resist bus bunching[END_REF] because they only require expected arrival times, but it may a↵ect the method in [START_REF] Berrebi | A real-time bus dispatching policy to minimize passenger wait on a high frequency route[END_REF], which uses probability distributions.

Equations 10, 11, and 12 show the probability distribution of a synthetic trajectory of the j th bus, denoted T [A j ], at time a i . The distribution is centered around a j + 1 , which is a uniformly distributed random variable, with interval length proportional to the horizon, a j a i , and to the accuracy indicator, ✏.

The random variable

2 is the error of trajectories around 1 , which is normally distributed with scale parameter proportional to the horizon and to the confidence parameter, . The choice of Uniform and Normal distributions with the horizon, a j a i as parameters is consistent with the distribution of errors on the particle filter used in Sections 4 [START_REF] Hans | Real-time bus route state 720 forecasting using particle filter[END_REF]. Note, however, that the random variable 1 is determined once for all synthetic simulated arrival times of the same bus, whereas 2 is re-determined for each particle.

T [A j ] = a j + 1 + 2 (10) 1 ⇠ U [ ✏(a j a i ), ✏(a j a i )] (11) 2 ⇠ N [0, (a j a i )] (12)
Figure 6 shows the sensitivity of (a) CV 2 and (b) mean lost time holding to the error terms, ✏ and . The methods recommended in [START_REF] Bartholdi | A self-coordinating bus route to resist bus bunching[END_REF] and [START_REF] Daganzo | Reducing bunching with bus-to-bus cooperation[END_REF] buses with less than half the CV 2 of the other two methods, but it requires more holding time as ✏ and increase. The method in [START_REF] Berrebi | A real-time bus dispatching policy to minimize passenger wait on a high frequency route[END_REF] can dispatch buses with slightly more stable headways when it considers the uncertainty around its prediction, , but considering more uncertainty requires much more holding time. The reason for this increase in lost time is that the expected maximum of random variables is a monotonous increasing function of their variability. The confidence also causes holding time, especially when ✏ is small, for the same reasons. These results indicate that it may be adequate to replace the joint probability distribution of bus arrival times by their expectations to sacrifice headway stability for holding time and computational simplicity.

Conclusion

In this paper, we compared the performance of closed-form bus holding methods used in practice and recommended in the literature on Tri-Met route 72 in Portland, Oregon. We applied control at one and several control points along the route and tested how each method holds vehicles to stabilize headways, reduce passenger waiting time, and prevent bus bunching. We used a new prediction tool developed in [START_REF] Hans | Real-time bus route state 720 forecasting using particle filter[END_REF] to simulate the performance of real-time holding methods, which was essential to apply several holding methods studied.

In addition, we coupled the methods in [START_REF] Bartholdi | A self-coordinating bus route to resist bus bunching[END_REF] and [START_REF] Berrebi | A real-time bus dispatching policy to minimize passenger wait on a high frequency route[END_REF] with the method in [START_REF] Daganzo | Reducing bunching with bus-to-bus cooperation[END_REF] to produce hybrid methods that can be applied at several control points along the route.

In the simulation, we found the trade-o↵s between headway stability and holding time for each holding method The Schedule-Based methods can reduce the need for holding time but they are unable to stabilize headways. The Headway-Based methods can be parametrized and applied at several holding points to yield a wide range of holding times. The methods, however, are not competitive in terms of headway stability for any level of holding time.

We found that Prediction-Based holding methods coupled with the predic-tion method in [START_REF] Hans | Real-time bus route state 720 forecasting using particle filter[END_REF] achieved the best compromises between headway regularity and holding time on a wide range of desired trade-o↵s. In particular, the method in [START_REF] Daganzo | Reducing bunching with bus-to-bus cooperation[END_REF] was the most e↵ective for stops where little holding can be implemented at once and the method in [START_REF] Berrebi | A real-time bus dispatching policy to minimize passenger wait on a high frequency route[END_REF] was the most e↵ective when vehicles could be held for longer periods of time. A Hybrid between the two methods, applied at one or several holding points produced the lowest headway CV 2 for a wide range of possible holding times. The latter analysis, however was made comparing holding methods on di↵erent route portions due to the lack of data on previous vehicle trips.

In order to model the e↵ects of cyclical vehicle trips and to compare holding methods in a more similar framework, future research should test and compare holding methods on a real bus route, and include optimization-based holding methods in the analysis.

In particular, a hybrid between the method from [START_REF] Berrebi | A real-time bus dispatching policy to minimize passenger wait on a high frequency route[END_REF] and [START_REF] Daganzo | Reducing bunching with bus-to-bus cooperation[END_REF], applied at one or several holding points, achieved the best compromises between headway regularity and holding time on a wide range of desired trade-o↵s.

When evaluating the impact of prediction accuracy on the performance of holding mechanisms, we found that prediction errors increased headway instability of real-time holding methods. Inaccurate predictions substantially increase mean lost time holding for the method in [START_REF] Berrebi | A real-time bus dispatching policy to minimize passenger wait on a high frequency route[END_REF]. Using high quality predictions such as those in Section 4 for real-time holding methods is necessary to maintain stable headways and short holding times. The uncertainty considered in the prediction for the method in [START_REF] Berrebi | A real-time bus dispatching policy to minimize passenger wait on a high frequency route[END_REF] helps reduce CV 2 at a high cost of holding time, especially for highly accurate predictions.

Replacing the probability distribution of expected arrival times by their mean could help reduce the mean lost time holding. This paper compares holding methods assuming Poisson distributed arrivals, but this assumption is not always true [START_REF] Luethi | Passenger arrival rates at public transport stations[END_REF]. The performance metrics do not reflect the value that schedules o↵er for both passengers and operators. In addition, on certain routes, schedule or real-time information coordinated arrivals may impact the performance of holding methods compared.

Future research should compare holding methods based on di↵erent passenger arrival distributions.

In this paper, TriMet Route 72 is used as a test bed for holding methods used in practice and recommended in the literature. There are routes resembling 630 Tri-Met Route 72 in terms of passenger demand, tra c congestion, and land-use in almost every metro area in the United States. Route 72 is a typical highfrequency route (7-8 minute headways in peak hours) that faces the issue of bus bunching. We considered Route 72 as generally as possible, often applying sensitivity analyses. The analysis presented in this paper can therefore provide a basis for transit agencies to decide on a closed-form holding method on their high frequency routes. Every holding method presented here strikes a di↵erent balance between the conflicting objectives to stabilize headways and dispatch buses with little holding time.

As transit agencies look to implement innovative holding methods, further

3. 1

 1 . The route Holding methods were compared by simulation using data on real bus trajectories from Tri-Met bus route 72 shown in Figure 1. Buses run in mixed tra c on the perpendicular 82 nd Avenue and Killingsworth Street in Portland,

Figure 1 :

 1 Figure 1: Map of TriMet Route 72 TriMet (2016)

  mean s , stdev s ) mean s Mean travel time between stop s and s + 1 stdev s Standard deviation of travel time between s and s + 1 When a bus arrives at the control point, 100 particles (simulated bus trajectories) are generated for each following bus traveling on the route. 11 The arrival times of following buses at the control point are then aggregated in a histogram, which are treated as probability distributions. Figure 2 shows a time-space diagram with the trajectories of following vehicles in part (a), and the associated histogram of bus arrival times in part (b). At each station, the particles in part 11 We chose to generate 100 particles as a compromise between the computational time and the resolution of the histogram.

Figure 2 :

 2 Figure 2: (a) Generation of particles to simulate bus trajectories. (b) Histograms of simulated bus arrival times treated as probability distributions.

Figure 3 :

 3 Figure 3: Performance of holding methods in terms of (a) squared headway coe cient of variation, (CV 2 ) and (b) mean holding time past boarding and alighting (lost time).

Figure 4

 4 Figure 4 shows (a) CV 2 as a function of ↵, (b) mean lost time as a function

Figure 4 :

 4 Figure 4: (a) CV 2 as a function of ↵, (b) mean lost time as a function of ↵ and (c) CV 2 as a function of mean lost time at station 49.

Figure 5

 5 Figure 5 shows mean CV 2 and total lost time applied and recorded throughout the route in Part (a) and applied and recorded in the second half of the route in Part (b). All holding methods were implemented at 1,2,4, and 8 holding points. The Schedule-Based methods were also held at 16 and 32 holding points because they are the only ones that do not consider headways. All other methods use either the forward or the backward headway, which would be affected by holds imposed at intermediate bus stops. The ↵ parameter was set to

Figure 5 :

 5 Figure 5: Mean CV 2 and total lost time applied throughout the route in Part (a) and applied in the second half of the route in Part (b)

  are shown with = 0 because they only require E[A j ]. The method proposed in Berrebi et al. (2015) is shown with = {0, 0.2, 0.3, 0.4} 17 . For each value of ✏, solid lines show the case where = 0, i.e. all trajectories equal a j + 1 . These lines describe the simulated performance of each holding method that would result if transit agencies used expected arrival times instead of probability distributions. The dashed lines 550 describe the outcome of considering uncertainty, , surrounding the expected synthetic trajectories, which contain error distributed with accuracy parameter, ✏.

640

  testing and simulation should be done for cases of unique passenger loading rules, severe passenger overflow, and other route characteristics that substantially deviate from TriMet Route 72. Future research should explore how route characteristics such as travel patterns, route instability, and perception of waiting time a↵ect the desired trade-o↵s between these conflicting objectives.8. AcknowledgmentThis research was partly funded by the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation program (grant agreement No 646592 MAGnUM project).

Table 1 :

 1 Summary of variable definitions.

		Variable Definition
		A i	Arrival time of bus i (Random Variable)
		a i	Arrival time of bus i (Realization)
		h i	Hold imposed on bus i
		d i	Departure time of bus i (Realization)
		di	Scheduled departure time of bus i
		H i	Scheduled headway of bus i
		n i	Number of following buses when bus i reaches the control point
		CV 2	Headway Coe cient of Variation squared
		T [A j ]	Arrival time of a particle 4 of the j th bus
	110	2.1. Naive methods

Table 2

 2

	. The Naive Schedule method

  need the expected next time, E[A i+1 ]. The method recommended in

	9
	http://portal.its.pdx.edu/Portal/index.php/fhwa
	10 Several holding point are considered in Section 5.2.

arrival

Table 3 :

 3 Particle generation Travel time of bus i between stop s and s + 1

	Equation		Variable Definition
				s	ctrl 1	Upstream stop from the control point
	A i =	s ctrl 1 P	i,s + ⇡ i,s	s	last	Last stop visited by bus i
		s=s last		i,s	Dwell time of bus i at stop s
				⇡ i,s
				a	Individual alighting time
				b		Individual boarding time
				c		Time lost in door opening and closing
				N A	Number of alighting passengers
				N B	Number of boarding passengers
	N

i,s = aN A + bN B + c B ⇠ P(d s h i,s ) h i,s Headway of bus i at stop s d s Demand rate at stop s N A ⇠ Bin(L i,s , µ s ) L i,s

Load of bus i at its departure from stop s µ s Alighting proportion at stop s L

  Performance indicators allow transit agencies to identify and address gaps in service. On high frequency routes, a measure of instability is the squared coe cient of variation of headways, denoted CV 2 as shown in Equation 8[START_REF] Kittelson | Transit capacity and quality of service manual[END_REF]. The variable measures the extent of headway instability. When CV 2 = 0, buses are evenly spaced, and when CV 2 = 1, Average Passenger Wait (APW) is equal to H. Assuming Poisson arrivals, APW is directly proportional to CV 2

		3.3. Performance indicators
		The model is simple to calibrate, and it is compatible with any bus route and
		any data format. The tool can consider tra c congestion and tra c signal
		data if available, but in this study, we assumed that only vehicle locations were
	270	available for the prediction

In this text, the term headway will be used as the time between the passing of two consecutive buses at a single location. Later, we will distinguish between the headway (or forward headway) and the backward headway, which is the time until the following vehicle will reach the current location.

The forward headway in[START_REF] Daganzo | A headway-based approach to eliminate bus bunching: Systematic analysis and comparisons[END_REF],[START_REF] Daganzo | A headway-based approach to eliminate bus bunching: Systematic analysis and comparisons[END_REF][START_REF] Daganzo | Reducing bunching with bus-to-bus cooperation[END_REF] Xuan et al. (2011) is expressed in terms of inter-arrival time, a i a i 1 without considering the hold imposed on the leading bus. To make these methods more robust, we have replaced the inter-arrival time by the time since last departure for the forward headway.

To preserve the robustness of their method, and recover from perturbations,[START_REF] Daganzo | A headway-based approach to eliminate bus bunching: Systematic analysis and comparisons[END_REF],[START_REF] Daganzo | A headway-based approach to eliminate bus bunching: Systematic analysis and comparisons[END_REF][START_REF] Daganzo | Reducing bunching with bus-to-bus cooperation[END_REF] Xuan et al. (2011) recommend using slack time. Slack time leverages longer holding time to stabilize headways, as shown in[START_REF] Argote-Cabanero | Dynamic control of complex transit systems[END_REF]. We have found in a simulation, however, that adding slack time to the Route

schedule does not substantially a↵ect the trade-o↵ between headway stability and holding time This is why we decided to calibrate slack time to the historical schedule in this paper.

Sometimes overlapping with[START_REF] Daganzo | A headway-based approach to eliminate bus bunching: Systematic analysis and comparisons[END_REF] in parts (a) and (b).

We did not include = 0.1 because it closely overlapped with = 0.

The methods recommended in [START_REF] Daganzo | Reducing bunching with bus-to-bus cooperation[END_REF] and in [START_REF] Bartholdi | A self-coordinating bus route to resist bus bunching[END_REF] are a↵ected by prediction accuracy in similar ways. When ✏ is null (perfect predictions), they dispatch buses with almost exactly the same CV 2 . As ✏ increases, they both destabilize at a fast rate (although CV 2 grows at a slightly greater rate for the method in [START_REF] Bartholdi | A self-coordinating bus route to resist bus bunching[END_REF]).

The error in accuracy, ✏, causes these methods to dispatch buses too soon sometimes and too late other times. The method in [START_REF] Daganzo | Reducing bunching with bus-to-bus cooperation[END_REF] 560 can dispatch buses with roughly half of the mean lost time in [START_REF] Bartholdi | A self-coordinating bus route to resist bus bunching[END_REF]. The mean lost time is only slightly a↵ected by increasing ✏ because the methods consider the expected arrival time, whose error is centered around the true mean.

For every value of ✏ and , the method in [START_REF] Berrebi | A real-time bus dispatching policy to minimize passenger wait on a high frequency route[END_REF]


D Faille

P Haessig

Heliostat aiming points optimization for Concentrated Solar Power plant

Keywords: Optimization, Solar Plant, Heliostat, Aiming

EDF is investigating the dynamic behavior of a tower Concentrated Solar Power plant which consists of a heliostat field, a receiver and a boiler. The paper presents a solution to optimize the aiming point for each heliostat. The method is based on a simplified Gaussian model that reproduces the flux distribution for each heliostat, as a function of sun position. This model is calibrated using flux data generated with Tonatiuh, a ray-tracing software. The total power for the heliostat field is computed neglecting the blocking and shading effect by summing up the contributions of the heliostats. A set of aiming points is obtained by different strategies to optimize the distribution of the flux on the receiver resulting from a compromise between flux spillage and flux uniformity.

INTRODUCTION

To limit the emission of greenhouse gases and their impact on climate change, administrations all over the world have engaged an energy transition program which will increase the share of renewable energies such as wind, hydro and solar. The major drawbacks of these renewables are their limited dispatchability and intermittency which makes their integration into the grid more difficult than fuel-based generation. Different storage technologies such as batteries, hydrogen, etc. are thus being developed to facilitate this grid integration. Because thermal storage is a rather inexpensive solution, the IEA technology roadmap [START_REF] Haessig | Technology Roadmap: Solar Thermal Electricity[END_REF] expects that concentrated solar technology could represent around ten percent of the global share of electricity generation by 2050. CSP plants produce heat which can be stored for later use or hybridized with fossil plants in order to decrease their environmental impact.

EDF develops a low carbon energy mix that includes wind and photovoltaic plants. Concentrated solar power is also a technology envisaged for the future. EDF already has experience in this field with the construction and operation of the CSP plant Themis, located in the South of France, from 1983 to 1986. A collaboration with the Chinese Academy of Science (CAS) has recently been engaged to study tower CSP technology and in particular control aspects which are particularly important for disturbance rejection and flexible operations. Studies have been done for Direct Steam Generation as well as Molten Salt (MS) CSP plants. A Direct Steam Generation control solution is presented in [2, 3, and 4]. The optimization of heat storage operations for MS CSP plant is considered in [START_REF]Control Solutions for Molten Salt Tower Concentrated Solar Plant[END_REF].

The optimization of aiming points of the heliostat is an important aspect for the tower CSP control. A control strategy leading to a more efficient operation of the field by maximizing the power while keeping the receiver in operating conditions may be of great value. Moreover, aiming set-point can be changed to adapt the plant to disturbances such as cloud, wind, etc. and mitigate the wear and tear due to the transients. An efficient algorithm to make this adaptation in real time is however a challenging task due to the complexity of the process and its uncertainties and only a few papers address this problem in the control literature.

A control solution implemented in the Plataforma Solar de Almeria (PSA) Platform in Spain is presented in [START_REF] Camacho | Advanced Control of Solar Plants[END_REF]. The algorithm is a heuristic which tries to reproduce what an operator would do without using any sophisticated models of the heliostats projection on the receiver. The control is based on temperature measurements (using thermocouples) and defines a discrete set of possible aiming points for the heliostats; the algorithm assigns each heliostat to an aiming point such that the temperature difference between the hottest and coldest aiming point temperature is less than 100°C. A comparison of the temperature distribution with and without the heliostat control solution shows that in the first case the temperature difference is two times smaller than in the second case, which leads to less mechanical fatigue.

Optimization of aiming points is not only useful for the operation stage, but also in the design stage of a solar project, to predict the annual energy yield of a given heliostat layout. Tools developed by NREL such as DELSOL or SolarPILOT can be useful to check a layout but the aiming strategy remains a kind of black-box and it might be difficult to adapt them in order to integrate new features. These tools are based on analytical models which represent the flux distribution for one heliostat (or group of heliostats) by analytical functions (Gaussian for instance) and computes the projection on the receiver panes taking into account the blocking and shading effects. Based on this kind of models, evolutionary algorithms are proposed in [START_REF] Salomé | Control of the flux distribution on a solar tower receiver using an optimized aiming point strategy: Application to THEMIS solar tower[END_REF] to optimize the distribution for the Themis CSP plant receiver. A heuristic approach is also proposed to define the aiming points for the Gemasolar plant in [START_REF] Sánchez-González | Solar flux distribution on central receivers: A projection method from analytic function[END_REF]; the number of heliostats in this case is an order of magnitude greater compared to the Themis case.

A generic model based approach is needed because of the diversity of the encountered situations. The objective and constraints to take into account depend on the characteristics of the receiver (geometry, heat transfer fluid, material, etc.). It can be an external receiver if the flux is distributed around the receiver as Gemasolar for instance, the heliostats are located around the tower in that case. The configuration adopted for the heliostat field of the Badaling platform is considered in the present paper. It is a cavity or plane receiver, so the heliostats are aiming only at one side of the tower and are located at the north of the tower as the platform is located in the northern hemisphere.

The present paper proposes a generic model based approach to optimize the aiming points of the heliostats to maximize the power while avoiding high flux. The plant description and the objective is given in Section 2. The assumption and the models of the heliostat field is addressed in Section 3. The optimization procedure and the main results are in Section 4. A conclusion gives the perspective for future developments.

SOLAR POWER PLANT DESCRIPTION

The process shown in Figure 1, which will be developed and tested on the Badaling facility owned by IEE-CAS, is an experimental tower CSP plant with a field of one hundred heliostats. The diagram on the right of Figure 1 corresponds to a Molten Salt plant with indirect steam generation. The primary circuit, which goes through the receiver, is filled with a heat transfer fluid made of molten salts that can be stored in tanks and exchange heat with the water and steam circuit inside the boiler. For the left diagram, which corresponds to Direct Steam Generation, the heliostats focus on the heat exchangers located at the top of the tower; the generated steam can be used directly in the turbine or stored in the steam accumulator.

The present paper addresses the problem of heliostat control which is common for both configurations described in Figure 1. An accurate control of the flux is required to avoid overheating of the receiver tubes, especially during load transient, startup, and shutdown. Otherwise creep can damage the pipes and lead to severe degradation. Another possible degradation mechanism is thermal fatigue due to the temperature gradient inside the receiver which may result from an incorrect flux distribution. Finally, temperature constraints are also imposed by the molten salt to avoid corrosion.

The Badaling heliostat field is made of 108 heliostats, see Figure 2. The layout shown in the figure is a theoretical layout defined in [START_REF] Zhang | Dynamic simulation of a 1MWe concentrated solar power tower plant system with Dymola®[END_REF]. This layout was implemented in NREL Soltrace ray-tracing software to derive flux maps used for dynamic simulations with the CSP plant Dymola® model. In the present paper, the field is implemented in Tonatiuh [START_REF] Tonatiuh | [END_REF], a software for which Matlab programs and scripts are available to pre and post process the simulations.

The model in Tonatiuh version 2.2.1 is shown in Figure 3 for a given sun position. The sun is defined by the two angular positions reproduced in Figure 4: azimuth as and elevation hs [START_REF] Camacho | Advanced Control of Solar Plants[END_REF]. The envisioned control of the CSP power plant is hierarchical and consists of two levels. The upper level is the plant controller which optimizes the set points for the lower level controllers (e.g. trackers). One of the most important upper level task is to control the aiming point of each heliostat. This is an optimization problem, where the objective should be a compromise between limiting the flux spillage on one hand and the limiting peak flux on the other hand. This aiming point optimization is subsequently developed in order to be integrated into the plant controller.

HELIOSTAT FIELD MODELING

The analytical expression to compute the power P projected by a heliostat on the receiver can be found in [START_REF] Camacho | Advanced Control of Solar Plants[END_REF]. It depends on the Direct Normal Irradiance (DNI) 𝐸 𝑠𝑜𝑙 which varies with the day and weather, the heliostat surface S and the orientation angle of the heliostat 𝜃 which can be derived from the different angles 𝑎 𝑠 , ℎ 𝑠 , 𝜀 𝑎𝑛𝑑 𝜒 defined in Figure 4.

𝑃 = 𝜂 𝐸 𝑠𝑜𝑙 𝑆 cos𝜃 

The coefficient 𝜂 takes into account the heliostat efficiency and the atmospheric attenuation. Formula are proposed in the literature to compute the atmospheric attenuation, which depends on the distance between the heliostat and the receiver. The previous model is rather simple but it cannot be utilized in practice for several reasons. The first reason is that there might be shading effects due to the tower or other heliostats which impact the light received by the heliostat. There might also be blocking effects when the direct line between the heliostat and the receiver meets another heliostat. These blocking and shading effects impact the power P.

The second important reason is that this model only gives the total power received on the tower (in kW), but gives no information on the spatial distribution of the flux on the receiver (in kW/m²). Indeed, the light spot on the receiver is not a perfectly focused point. Reasons for this imperfect focusing are mainly because the sun cannot be assimilated to a point and because of the imperfections of the heliostats. The rays coming from the sun hit the heliostat in different points and are reflected in different directions such that the normal to the heliostat surface is the bisector of the angle between the incident and reflected ray.

The projection of the heliostat on the receiver is thus a fuzzy spot with rather ellipsoidal shape centered on the aiming point. There are analytical models which take into account all aforementioned effects, such as HFLCAL or UNIZAR [START_REF] Sánchez-González | Solar flux distribution on central receivers: A projection method from analytic function[END_REF]. However, it is difficult to find the value for the many parameters of these models. The alternative approach we use is to model the heliostat field is using the data obtained with Tonatiuh and directly fit the flux map with Gaussian functions. Notice that this fitting approach can also be applied to real measurements of flux maps.

Gaussian distribution

The results of Tonatiuh simulations are the impact points of the rays on the receiver. For each heliostat, a statistic of this two-dimension data set can be computed.

We compute 𝜇 = (𝜇 𝑥 , 𝜇 𝑦 ) 𝑇 the mean value (barycenter of the energy) and Σ = [ 𝜎 𝑥 2 𝜎 𝑥𝑦 𝜎 𝑥𝑦 𝜎 𝑦 2 ] the covariance matrix. With these parameters, the flux distribution is approximated by a Gaussian functions 𝜑 defined on the coordinate 𝑋 = (𝑥, 𝑦) 𝑇 of the receiver plane:

𝜑(𝑥, 𝑦) = 1 2𝜋 √𝑑𝑒𝑡Σ exp - 1 2 (X -𝜇) 𝑇 Σ -1 (X -𝜇)
It can be mentioned that the covariance matrix Σ is an aggregation of the errors due to the sun diameter, the imperfection of the tracking and the imperfections of the mirror considered in the models found in the literature [START_REF] Sánchez-González | Solar flux distribution on central receivers: A projection method from analytic function[END_REF].

The function 𝜑 is a normalized distribution function which is multiplied by the power P to obtain the flux on the receiver plane:

Φ(x, y) = 𝑃. 𝜑(𝑥, 𝑦)

Parameterization

For a given heliostat, the power P and the covariance matrix Σ which characterizes the shape of the spot, depend only on the position of the sun given by the azimuth and elevation, not on the aiming points. Indeed, numerous fitting results showed us that: (i) center of the distribution 𝜇 is located at the aiming point and (ii) P and Σ are independent of the aiming point. Practically, this means that, when moving the aiming point, the flux distribution is translated but not distorted.

For the power P, the analytical expression given above can be the same as the value given by Tonatiuh if there are no shading or blocking effects, otherwise a correction must be made. The variance and covariance coefficients 𝜎 𝑥𝑦 , 𝜎 𝑥 , and 𝜎 𝑦 are computed for different positions of the sun and the results are interpolated with 2D polynomials of order 2×3.

Field model

The total flux Φ tot on the receiver plane is obtained by adding the contributions of each heliostat where Nh is the number of heliostats, (x,y) are the coordinates on the receiver plane and (𝑎 𝑠 , ℎ 𝑠 ) are the azimuth and elevation angle of the sun, and (𝑥 𝑎𝑝 𝑖 , 𝑦 𝑎𝑝 𝑖 ) is the aiming point of each heliostat.

Φ tot (x, y, 𝑎 𝑠 , ℎ 𝑠 , 𝑎𝑝) = ∑ 𝑃 𝑖 (𝑎 𝑠 , ℎ 𝑠 ). 𝜑 𝑖 (𝑥 -𝑥 𝑎𝑝 𝑖 , 𝑦 -𝑦 𝑎𝑝 𝑖 , 𝑎 𝑠 , ℎ 𝑠 ) 𝑁ℎ 𝑖=1

For keeping the notation compact, we name 𝑎𝑝 the set of all aiming points: 𝑎𝑝 = {𝑥 𝑎𝑝 𝑖 , 𝑦 𝑎𝑝 𝑖 } 𝑖=1..𝑁ℎ .

The integration of the 2D Gaussian being equal to one, 𝑃 𝑟𝑒𝑓 the total power received by the infinite plane containing the receiver surface is given by the equation:

𝑃 𝑟𝑒𝑓 = ∑ 𝑃 𝑖 (𝑎 𝑠 , ℎ 𝑠 ) 𝑁ℎ 𝑖=1
.

However, the receiver is a limited five-by-five-meter square and a part of the flux reflected by the heliostat does not hit the receiver. The power received 𝑃 𝑡𝑜𝑡 is obtained by limiting the integration to the receiver surface:

𝑃 𝑡𝑜𝑡 (𝑎𝑝) = ∬ Φ tot (x, y, 𝑎 𝑠 , ℎ 𝑠 , 𝑎𝑝) 𝑑𝑥𝑑𝑦 𝑅𝑒𝑐𝑒𝑖𝑣𝑒𝑟 𝑆𝑢𝑟𝑓𝑎𝑐𝑒
Practically, this surface integral is computed numerically by dividing the receiver in 50x50 small squares. First results show a mean error between Tonatiuh and the approximation for the flux around 1%. Further analysis must be done to assess the precision of the model. However the model is much faster than the Tonatiuh model with a speedup factor 100. Finally, the proposed model is implemented in Matlab and can be integrated easily in optimization routines. This is also a major improvement over the original ray tracing model, since current version of Tonatiuh cannot be controlled by an external process (no API or command line automation).

AIMING POINTS OPTIMIZATION

Reference configuration

A natural strategy, if one wants to maximize the power, is to aim at the center of the receiver. The resulting flux is shown in Figure 5 for this reference configuration. The total power 𝑃 𝑡𝑜𝑡 is equal to 10.4 MW but the power is concentrated in a small part located around the center of the receiver. This aiming strategy leads to a high flux at the center (Φ max = 3793 kW/m²), which may not be compatible with the tube material. Even in this reference configuration, a loss estimated to 4.7% occurs.

Multi objective Optimization

To avoid peak flux, the idea is to spread the received power without losing too much energy. This is a multi-objective optimization between loss and spread for which a mathematical formulation is proposed in the following.

Loss ratio. The loss ratio 𝐽 𝐿𝑜𝑠𝑠 is defined below as the normalized lost power, a value between 0 and 1. Zero is the ideal value, but it cannot be reached due to imperfect focusing (as mentioned above, in the reference configuration for which they are minimum, losses are estimated to 4.7%).

𝐽 𝐿𝑜𝑠𝑠 = 1 - 𝑃 𝑡𝑜𝑡 𝑃 𝑟𝑒𝑓
, 𝑃 𝑡𝑜𝑡 , 𝑃 𝑟𝑒𝑓 are defined in the previous paragraph.

Spread ratio. We first define the spread Δ Φ 𝑡𝑜𝑡 by the difference between the maximum and the minimum of the flux in kW/m2. This value is not dimensionless, so to facilitate the comparison with the loss criterion, it is normalized by its value in the reference configuration Δ Φ 𝑟𝑒𝑓 . The spread ratio 𝐽 𝑆𝑝𝑟𝑒𝑎𝑑 is again a number between 0 and 1. One is the worst value and it is attained in the reference configuration (all aiming point superimposed at the center).

𝐽 𝑆𝑝𝑟𝑒𝑎𝑑 = Δ Φ 𝑡𝑜𝑡 Δ Φ 𝑟𝑒𝑓

A spread of zero is the ideal value since it means that the flux would be uniformly distributed. However, we expect that for such a configuration the loss would be too big and inacceptable. So the aiming point strategy must find a compromise between the minimization of the loss and the spread ratio by the minimization of the combined criterion 𝐽 𝛼 , where 𝛼 is a tuning factor.

Optimization Algorithms

The multi objective optimization defined in the previous section is a difficult problem without any analytical solution. A numerical optimization routine must be used to solve the problem. One can utilize unconstrained minimization algorithm or constrained optimization to take some constraints into account like forcing the aiming points to be inside the receiver. Specifically, we used fmin and fmincon from Matlab Optimization toolbox.

Two approaches are defined for the solving the optimization problem: a direct approach where optimization unknowns are the aiming points coordinates {𝑥 𝑎𝑝 𝑖 , 𝑦 𝑎𝑝 𝑖 } 𝑖=1..𝑁ℎ (2×108 real variables); an indirect approach where optimization unknowns are the parameters of a parametric model (with a small number of parameters), the aiming points being computed with the parametric model. Notice that in this optimization section, we drop the "ap" indices from the aiming point variables for clarity.

Direct Approach.

The optimization problem is to find the aiming point heliostats coordinates 𝑥 * , y * in the receiver axes which minimize the objective function:

𝑥 * , y * = arg min 𝑥,𝑦 𝐽 𝛼 (𝑥, 𝑦)
Constraints can be added to limit the domain of research or to requirements for the solution.

The number of unknown is two times the number of heliostats (as there is one aiming point for each heliostat) which is a relatively important number, even for the small-sized field considered in the present paper.

To reduce the complexity of the optimization, we suppose that the problem can be decomposed in sub problems by considering smaller groups of heliostats. If we consider that the objective function is additive with respect to the heliostat contributions (this assumption implies that no coupling exists between the optimization variables resulting from blocking, shading, etc.) and that the objective function terms for each heliostat are positive, we can decompose the optimization. Where (𝑥 𝑖 , 𝑦 𝑖 ) are the aiming points for the heliostats of the i th group. We could consider that the way of creating the group is not important. It is true in theory but the function J is not convex and we might get stuck in a local minimum more easily with certain grouping configurations.

Indirect Approach If we suppose that the minimum is unique, the optimization creates a mapping function from each heliostat given by its coordinate (X,Y) in the field coordinate shown in Figure 2 to the aiming points with coordinate (x,y) in the receiver coordinates. The idea of the indirect approach is to find a parametric model to make assumption on the mapping function by stating a parametric a priori model and optimize the parameters of the mapping function.

For instance, if (, ) are the polar coordinates of the heliostat, we can envisage the following mapping with two parameters A and B to be optimized and a fixed parameter k.

𝑥 = (Α𝜌 + 𝐵) cos(𝑘𝜃) 𝑦 = -(Α𝜌 + 𝐵) sin(𝑘𝜃)

The optimization problem becomes:

𝐴 * , B * = arg min 𝐴,𝐵 𝐽 𝛼 (𝑥(𝐴, 𝐵), 𝑦(𝐴, 𝐵))

The dimension is dramatically reduced compared to the direct optimization. The objective function is more complex to compute and the convexity is no more guaranteed than in the direct approach as can be observed in Figure 6. The optimization is however two orders of magnitude faster than the direct approach. Of course, more complex mapping functions would lead to slower optimization, but better objective values.

Results

The different optimization strategies described in the previous paragraphs were applied to the case study and are discussed in this paragraph.

The direct method optimized the entire field heliostats aiming points. The initial guess is chosen randomly. Several runs were made and the optimization generally ends with the maximum of iterations exceeded (n=100000), but results are acceptable even if the minimum is not reached.

The heliostat field is divided in four groups by taking shifted heliostat indices sets {1,5,9, …}, {2,6, ...}, {3,7,11, …}, {4,8,12, …}. Two initialization strategies have been adopted. The variant cold start assumes that the initial guess for each sub problem is chosen randomly. For the variant warm start the initial guess for the first group is chosen randomly and the solution of the optimization is chosen as a starting point for the next group, etc.

The results are given in Table 1. On the one hand, we can observe that the optimization time is divided by a speedup factor three by dividing the field in four groups and that the warm start does not bring significant improvement. On the other hand, a speedup factor 70 is obtained with the indirect approach. As the problem is non convex, several optimizations with randomly chosen starting point can be launched. With seven starting points, a speed up factor ten is obtained for the optimization -which can be performed in less than a minute. This appears to be sufficiently fast to be applied on line for control applications. The Pareto front of the optimization with the indirect method is shown in Figure 7. We can see the effect of the parameter  on each objective. This parameter can be used as a tuning parameter to optimize the set points.

Strategy

CONCLUSION

A simple heliostat model is proposed for aiming point optimization. The flux distribution of each heliostat is approximated by a 2D-Gaussian model whose parameters are interpolated to obtain a model valid for all sun positions. To limit the complexity, the following assumptions are made. The center of the spot coincides with the aiming point. The vibration of the mirror caused by the wind (or other external disturbances) is neglected so that the aiming point set point is perfectly followed. A calibration procedure described in [START_REF] Camacho | Advanced Control of Solar Plants[END_REF] must be done periodically to ensure this condition. The effect of blockage and shading between the heliostats are also neglected. The parameters of the 2D-Gaussian (x,xy,y) are identified with data obtained with the Monte Carlo ray tracing software Tonatiuh. Identification procedures can also be applied with on-site measurement. The optimization objective is to flatten as much as possible the flux distribution on the aperture plane and to simultaneously minimize the spillage. To quantify these objectives, two criteria for the spillage and the spread are defined. The heliostats aiming points are computed to minimize the multi objective function 𝐽 𝛼 . The parameter  can be adjusted to obtain the best compromise between the spillage and the spread.

In the perspective of designing a solution capable of modifying the aiming points in real time to take into account disturbances from the field -such as wind or clouds -or from the receiver -such as trip of pump or turbine -a fast solution was considered. Two approaches were adopted: the first one being to decompose the field in subsets of heliostats, the other one being to consider an a priori law for the aiming point and optimize the parameters of this mapping law. The mapping approach gives good results with an optimization done within a minute and is therefore a good candidate for on-line optimization.

Figure 1 :

 1 Figure 1 : Indirect and Direct Steam Generation Central Receiver System

Figure 2 :

 2 Figure 2 : Heliostat Field [9]

Figure 3 :

 3 Figure 3 : Simulation with Tonatiuth (1000 rays, azimuth 0 deg, elevation 80 deg)

Tour

Figure 4 :

 4 Figure 4 : Definition of sun angular position and incidence

Figure 5 :

 5 Figure 5 : Flux distribution in kW/m² for the reference configuration which maximizes the power 𝑷 𝒕𝒐𝒕 .

  𝐽 𝑆𝑝𝑟𝑒𝑎𝑑 + 𝛼 𝐽 𝐿𝑜𝑠𝑠

  𝑥 𝑖 , 𝑦 𝑖 )

Figure 6 :

 6 Figure 6 : indirect method optimization function.

Figure 7 :

 7 Figure 7 : Pareto front for the indirect optimization

Table 1 : optimization results

 1

		J α	Optimization time
			(second)
	Direct-1 group	0.3514	419.1
	Direct -4 groups -cold start	0.3479	180.8
	Direct -4 groups -warm start	0.3449	168.5
	Indirect	0.3847	6.8

AKNOWLEGMENT

This work has been done in the frame of a collaboration with EDF China and IEE-CAS for the Badaling experimental platform. The authors would like to acknowledge the technical contributions of Supelec students Pan XuYang, Weicai Li, Paul Fino, and Baptiste Balique.