L Floch

P Dewals

I Terrier

P Chavance

Statistics of the French purse seine fishing fleet targeting tropical tunas in the Indian Ocean (1981-2011)

Keywords: fish aggregating device, free-swimming school, Katsuwonus pelamis, purse seine fishing, Thunnus albacares, Thunnus obesus

In 2011, the French purse seine fishing fleet of the Indian Ocean was composed of 8 large size vessels that represented a total carrying capacity of about 8,000 t. Catches reached a total of 43,000 t and were composed of 50%, 42%, and 8% of yellowfin, skipjack, and bigeye, respectively. After a period of increase during 2006-2008, the fishing effort of the fleet has been decreasing to reach a minimum of 1,800 searching days in 2011. The decrease in effort was associated with a contraction of the fleet fishing grounds in the recent years and mainly characterized by a strong decrease in the number of sets made on free-swimming schools; a total of 1,800 fishing sets being made in 2011 compared to about 4,200 y -1 in the mid-2000s. Hence, the percentage of sets made on FAD-associated schools steadily increased since 2004 to reach more than 65% in 2010-2011, FAD-fishing resulting in 75% of the total catch of the French purse seine fishing fleet in 2011. Species-specific catch rates (in t per searching day) on FAD-associated schools reveal increasing trends over 1981-2011 with high interannual variability while no clear trend is apparent in the time series of catch rates on free-swimming schools. The year 2011 is characterized by a large proportion of yellowfin in the catch on FAD-associated schools that might indicate a good recruitment. For both fishing modes, the mean weight in the catch of the 3 tropical tunas has decreased since the mid-2000s with values in 2011 equal to 7.1 kg, 2.4 kg, and 5.4 kg for yellowfin, skipjack, and bigeye, respectively.

Introduction

Statistical data for the French purse seine fishing fleet have been collected by the "Institut de Recherche pour le Développement" (IRD; ex-ORSTOM) since the arrival of the first purse seiners in the early 1980s. French purse seiners target yellowfin (Thunnus albacares), skipjack (Katsuwonus pelamis), and bigeye tuna (Thunnus obesus) through two major fishing modes that result in different species and size composition of the catch: log/FADassociated (FAD) and free-swimming schools (FSC). The acronym "FAD", which stands for drifting fish aggregating device, is used here to describe any type of floating object used for increasing tuna catchability. This includes natural objects (e.g. logs, palm branches) and anthropogenic floating objects, such as man-made bamboo rafts equipped with radio-range beacons, satellite transmitters or scanning sonars. Fishing sets made on whales were classified as free-swimming school sets whereas sets made on whale sharks (Rhincodon typus) were classified as FAD sets [START_REF] Pallarés | Analyse du schéma d'échantillonnage multispécifique des thonidés tropicaux (Programme n • 95/37)[END_REF]. The fleet activities are described through a suite of fisheries indicators that provide information on fishing effort, catch, catch rates, size structure and mean weights for the 3 principal market tropical tunas, with a particular focus on the year 2011.

Fishing capacity and effort

Fishing fleet

The number of vessels of the French purse seine fishing fleet in the Indian Ocean rapidly increased from 2 in 1981 to a maximum of 26 during [1984][1985] with the arrival of vessels from the Atlantic Ocean. Since then, the number of vessels steadily decreased to reach a minimum of 8 in 2010-2011 (Fig. 1 and Table 1). The size of the vessels progressively increased in the French purse seine fishery over the last 30 years. The number of small-size vessels (capacity < 600 t) decreased throughout the 1980s and 1990s to become 0 in the mid-2000s while medium-size vessels (capacity between 601-800 t) disappeared from the fishery in the late 2000s (Fig. 1).

Carrying capacity

The total carrying capacity expressed in tonnage strongly increased from the arrival of the first purse seiner in 1981 to about 14,000 t in the mid-1980s. Despite the steady decrease in the total number of vessels, the capacity remained stable through the 1990s and 2000s 3 and increased to a maximum of 17,000 t in 2007 (Fig. 1). Indeed, the mean vessel-specific carrying capacity of the French purse seiners increased from about 600 t in the early 1980s to more than 1,000 t from the mid-2000s (Fig. 2; slope = +13 t -1 , p-value < 0.001). In the recent years, the capacity of the French fleet decreased to reach about 8,000 t during 2010-2011. Meanwhile, the capacity of the French overseas territories increased from less than 500 t in 2006 to more than 5,000 t in 2011 due to the construction of new purse seiners [START_REF] Chassot | Statistics of the purse seine fleet of the French overseas territories targeting tropical tunas in the Indian Ocean[END_REF].

Fishing and searching days

Fishing effort expressed in searching time (days) was computed by subtracting the time spent setting the gear from the fishing time. The time spent setting the gear was estimated by regressions linking duration and size of sets, from at-sea measurements made by scientific observers. The total number of fishing and searching days showed similar patterns over 1981-2011 with a major increase in the mid-1980s and a steady decrease until 2005 thereafter (Fig. 3). An increase in effort to a value of almost 5,000 fishing days was then observed due to the increase in the number of purse seiners from 14 in 2003 to 17 during 2006-2008. Since then, the nominal fishing effort strongly decreased and reached in 2011 values of about 2,110 and 1,800 fishing days and searching days, respectively.

Fishing grounds

The major fishing grounds of the French purse seine fishing fleet in 2011 were located in the north of the Mozambique Channel, west of the Seychelles, and along the Somali EEZ boundary (Fig. 4). A few fishing activities took place south of the Seychelles at longitude between 55 • and 60 • E and in the east of the Chagos EEZ.

After a major expansion of the fishing grounds during the early 1980s, the spatial extent of the fishery steadily increased over time from a total of about 400 squares in the mid-1980s to a mean value of 490 squares in the late 2000s (Fig. 5). The peak observed in 1998 corresponded to a particular climatic event that has been described elsewhere [START_REF] Murtuggude | Ocean color variability of the tropical Indo-Pacific basin observed by SeaWiFS during 1997-1998[END_REF], Murtugudde et al. 2000) and led the purse seiners to explore the eastern part of the Indian ocean for fishing (Ménard et al. 2006). The different selection criteria, i.e. minimum catch, effort, or sets, did not modify the temporal patterns in the spatial extent of the fishery (Fig. 5).

Fishing activities

The total annual number of fishing sets made by the French purse seine fleet showed a general pattern consistent with the annual variations in carrying capacity and fishing effort of the fleet (Fig. 6 and Table 4). The fleet showed a strong decrease in the number of sets from more than 4,300 in the mid-2000s to less than 2,000 in 2010-2011. The total number of sets showed a pattern very similar to the annual changes in the number of sets made on free-swimming schools (Fig. 6). The annual number of sets made on FAD-associated schools was more stable over time (around 1,700 y -1 ) and in opposite phase compared to sets made on free-swimming schools. The percentage of FAD-associated over free-swimming schools varied around a mean of 47% (SD = 10%) with FAD-associated fishing predominating from the mid-1990s to the early 2000s while sets made on free-swimming schools were more frequent during 2003-2008. The percentage of sets made on FAD-associated schools steadily increased from 37% in 2004 to reach more than 65% in 2010-2011.

The percentage of successful sets (i.e. positive tuna catch) made on FAD-associated schools has been high over the period 1981-2011, the mean annual value varying between a minimum of 85% in 1984 and a maximum of 95% in 1991. After a decrease during the development of the fishery, the percentage of successful sets on free-swimming schools steadily increased from a mean value of 37% in 1986 to an average of 54% during 2007-2011 (slope = +0.5% y -1 , p-value < 0.001). By contrast, the high success of sets on FAD-associated schools does not show any clear temporal pattern over the period. In the recent years, the success rate for FAD-associated schools seems to show an increasing trend that could be related to the improvement and increasing use of lateral echo sounders.

Fisheries production

Catch levels

The French purse seine fishery showed strong interannual variations in the catch over 1987-2008 with a minimum of 56,000 t in 1998 and a maximum of 108,000 t during 2003-2004 (Fig. 7). Since then, the total catch of the fleet strongly declined in relation with the decrease in the number of vessels and reached 43,000 t in 2011. Catches on FAD-associated schools represented on average 60% (SD = 8%) of the total fleet catch over 1981-2011. Following the large catch of tunas on free-swimming schools during 2003-2005, the percentage of FAD-associated catch steadily increased (slope = +4.3%, p-value < 0.001) to reach more than 70% during 2009-2011. Catches on FAD-associated schools appeared more stable than catches on free-swimming schools during 1981-2011 (Fig. 8). They were largely dominated by skipjack that represented an annual mean of 62% (SD = 6.5%) of all species caught, varying from 53% in 1997 to more than 75% in 1991 (Fig. 8a). Catches of bigeye, re-estimated from multispecies sampling operations conducting at unloading [START_REF] Pallarés | Analyse du schéma d'échantillonnage multispécifique des thonidés tropicaux (Programme n • 95/37)[END_REF], varied between 3.7% and 13% of the total of catch on FAD-associated schools during 1981-2011. Catches of yellowfin on FAD-associated schools were stable over the period representing on average 27% (SD = 5.2%) of the catch.

Catches made on free-swimming schools showed strong internannual variations in the last two decades from a maximum of 50,000 t in 2004 to a minimum of 8,000 t in 2011 (Fig. 8b).

The changes over time in free-swimming schools catches drove the trends in the total catch of the fishery. The catches were largely dominated by yellowfin that represented an annual average of 74% (SD = 9%) of the FSC catch during 1981-2010, and a similar value of 74% in 2011. However, while free-swimming schools represented more than 60% of the total catch of yellowfin during 1981-2010, it only represented 38% of the total yellowfin catch in 2011.

The predominance of yellowfin in the French purse seine catch from FAD-associated fishing (> 50% of the catch) had only been observed during 1997-1999 during the strong El-Niño event. In the recent period, this pattern stems from the strong decrease in the number of sets made on free-swimming schools which might reveal a change in purse seine fishing strategy due to a decrease in tuna free-swimming schools in conjunction with the Somali piracy threat [START_REF] Chassot | Analysis of the effects of Somali piracy on the European tuna purse seine fisheries of the Indian Ocean[END_REF], the implementation of the BIOT marine protected area that decreased purse seine fishing grounds as well as to the high sale price of skipjack (Table 4). In addition, the proportion of yellowfin in FAD-associated catch reached more than 40% in 2011 while it represented an average of 28% (SD = 5%) during 1982-2010 (see below).

Spatial distribution of the catch

The French purse seine fleet fishing grounds in 2011 were concentrated in the north of the Mozambique Channel and the Somali area where FAD-fishing predominated while fishing on free-swimming schools mainly occurred around and south-east of the Seychelles islands, and in the Mozambique to a lesser extent (Fig. 9). As compared to previous years (2006)(2007)(2008)(2009)(2010), the composition of the catch on FAD-associated schools in the Somali area was described by an important part of yellowfin that represented more than 50% of the catch in most of the 1x1 degree squares west of 55 • E and north of 10 • S . This might result from a good recruitment of yellowfin in 2011 as suggested by the size structure of the catch (see Section 3.4). Besides, catches on free-swimming schools in 2011 were predominated by skipjack in the Mozambique Channel (Fig. 13), as generally observed in this area, while the remaining of the catch was scattered along latitude 5 • S and predominated by yellowfin, with a few squares characterized by 15-20% of bigeye (Fig. 13).

Catch rates

Raw catch rates expressed in tonnes per searching day (t d -1 ) for the principal market tunas exhibit an increasing trend over 1981-2011, in particular from the early 1980s to the mid-1990s (Fig. 15). While skipjack and yellowfin are mainly caught in different fishing modes, their respective CPUEs were in the same order of magnitude, i.e. a mean value of 10.3 and 10 t per searching day for yellowfin and skipjack, respectively. Also, the CPUE time series showed very similar trends over time, with the notable exception of 2003-2005 11 described by abnormally high catch rates for yellowfin on free-swimming schools (Fig. 16b). Raw catch rates of bigeye also showed an increasing trend during 1981-2011 (slope = +4.4%, p-value <0.001) with an average of 1.88 t d -1 (SD = 0.4 t d -1 ).

Skipjack catch rates predominated on FAD-associated schools during 1981-2011 and increased quickly through the 1980s while stabilising around 10 t searching d -1 (SD = 2.3 t d -1 ) since the early 1990s (Fig. 16a). Yellowfin catch rates have shown an increasing trend since the early 1980s with some interannual variability. The maximum catch rates for yellowfin on FAD-associated schools have been observed in 2011, i.e. an annual value of 7.3 t d -1 . Bigeye raw catch rates on FAD-associated schools have remained stable around a mean value of 1.4 t d -1 since the late 1990s (Fig. 16a).

Yellowfin catch rates predominated for sets made on free-swimming schools and remained quite stable over time (mean of 7 t d -1 during 1984-2011), with abnormally high values during 2003-2005 > 13.5 t d -1 (Fig. 16b). In 2011, the catch rate of yellowfin in free-swimming schools was less than 4.5 t d -1 . Catch rates for skipjack have shown high interannual variability during 1986-2011 (mean = 2 t d -1 , SD = 1 t d -1 ) and a decreasing trend in the recent years, with catch rates < 1 t d -1 during 2009-2011. Catch rates for bigeye in free-swimming schools have shown a shift from a mean value around 0.3 t d -1 dur-

Size structure of the catch

The size structure of the catch (expressed in number) in 2011 appeared similar for the 3 tuna species to the average year representing the period 2006-2010, with the exception of skipjack caught on free-swimming schools that were smaller (Fig. 17). The sizes of yellowfin caught on FAD-associated schools were described by 2 modes for small (46 cm) and intermediate sizes (56 cm) while larger individuals (> 80 cm) were also caught but more seldomly. The year 2011 was characterized by a high peak of small yellowfin (size range of 42-48 cm fork length) which might indicate a good recruitment. Bigeye caught on FAD-associated schools in 2011 were only juvenile fishes of median size of 52 cm (sd = 9 cm). Skipjack showed a unique size mode described by a median of 46 cm (min = 31 cm and max = 69 cm) similar to 2006-2010 during which the median size of skipjack was 47 cm.

Similarly, the biomass of fish caught by size class showed overall similar patterns between 2011 and 2006-2010 except for skipjack caught smaller on free-swimming schools in 2011 (Fig. 18). For yellowfin caught on FAD-associated schools, the year 2011 was described by larger biomass of fishes comprised between 60-100 cm. In relation with the decrease in activities on free-swimming schools, the biomass of large yellowfin and bigeye (> 100 cm) caught was lower in 2011 than for the average year 2006-2010. The biomass of skipjack 44-50 cm long caught on FAD-associated schools, which represent the bulk of skipjack catch, was very similar in 2011 than during 2006-2010 while there was a strong decrease in skipjack catch on free-swimming schools (Table 567).

Mean weight in the catch

Time series of the mean weight in the catch for the principal market tropical tunas were computed as the species-specific annual biomass over the total number of fishes caught. The sampling design used to collect size-frequency data aboard French purse seiners thoughout the 1980s varied over time [START_REF] Pianet | Évolution du système de collecte et de traitement des données de la pêche thonière des senneurs européens et assimilés de 1981 à 1998[END_REF]) therefore the mean weight was computed from 1991 for data consistency. The mean weight highly differed between fishing modes and showed strong interannual variations during 1991-2011 (Fig. 19). The mean weight of yellowfin caught on free-swimming schools decreased from more than 37 kg in 1991 to a minimum of about 15 kg in 1998-1999, before progressively increasing thereafter to reach about 38 kg in 2004-2005. Since then, the mean weight decreased to about 28 kg during 2009-2011 (Fig. 19a). It is noteworthy that the yellowfin mean weight showed an increase from 25 kg in 2010 to 29 kg in 2011. After an initial decrease from about 10 kg in the early 1990s, the mean weight of yellowfin in the catch made on FAD-associated schools stabilised at around 5.7 kg (SD = ± 1.1 kg) from the mid-1990s to the late 2000s. In 2011, the mean weight of yellowfin caught on FAD-associated schools was less than 5 kg. Overall, the mean weight of yellowfin caught in the French purse seine fleet has shown a strong decrease from about 18 kg during 2003-2005 to about 7 kg in 2011 (slope = -1.4 kg -1 ; p-value < 0.001) due to increasing proportion of FAD-fishing in the fishery over the last decade (Fig. 6 -dotted line).

The overall mean weight of skipjack was more stable than yellowfin during 1991-2011 and driven by the mean weight of skipjack caught on FAD-associated schools. It varied between a maximum of 3 in 1992 and a minimum of 2 kg in 2008 for fishes caught on FAD-associated schools (Fig. 19b). After a decrease during 2007-2008, the mean weight remained stable at about 2.3 kg during 2009-2011. The mean weight of skipjack caught on free-swimming schools has shown a strong and significant decrease from more than 4.2 kg in 2003-2004 to 2.4 kg in 2009-2011 (slope = -0.3 kg -1 , p-value < 0.001).

The mean weight of bigeye varied around 5.8 kg (SD ± 1 kg) during 1991-2011 for the whole fishery (Fig. 19c). It showed strong interannual variations for individuals caught on free-swimming schools with a pattern very similar to yellowfin tuna (Pearson's r = 0.77; p-value < 0.001). The mean weight of bigeye caught on FAD-associated schools showed a significant decreasing trend during 1991-2011 (slope = -0.06 kg -1 , p-value < 0.01) with high interannual variations. In 2011, the mean weight of bigeye caught on free-swimming schools and FAD-associated schools was 18.3 kg and 3.8 kg, respectively.

Figure 1 :

 1 Figure 1: Fishing capacity of the French purse seine fleet in the Indian Ocean. Annual changes in the number of purse seiners by size category (barplots) and total carrying capacity (solid line with circles) during 1981-2011. Capacity was weighted by the vessel-specific proportion of the year at sea (in months). The vessel size category (t) was computed as 0.7 times the capacity expressed in m 3

Figure 2 :

 2 Figure 2: Changes in mean vessel carrying capacity for the French purse seine fishing fleet in the Indian Ocean, 1981-2011. The vessel-specific carrying capacity was weighted by the relative annual catch of each purse seiner. Solid line indicates mean linear regression

Figure 3 :

 3 Figure 3: Changes in nominal effort over time. Annual total number of fishing and searching days for the French purse seine fishing fleet in the Indian Ocean during 1981-2011

Figure 4 :

 4 Figure 4: Fishing grounds. Spatial distribution of fishing effort (in searching days) of the French purse seine fishing fleet in 2011

Figure 5 :Figure 6 :

 56 Figure 5: Changes in spatial extent of the fishery over time. Annual number of 1-degree squares explored by the French purse seine fishing fleet during 1981-2011

Figure 7 :

 7 Figure 7: Total fishery production. Catch by species of the French purse seine fishing fleet during 1981-2011

Figure 8 :

 8 Figure 8: Fishery production by major fishing mode. Catch by species of the French purse seine fishing fleet on (a) FAD-associated and (b) free-swimming schools during 1981-2011

Figure 9 :

 9 Figure 9: Spatial distribution of tuna catches of the French purse seine fishing fleet in 2011

Figure 10 :Figure 11 :Figure 12 :Figure 13 :Figure 14 :

 1011121314 Figure 10: Spatial distribution of tuna catches of the French purse seine fishing fleet in 2006-2010

Figure 15 :Figure 16 :

 1516 Figure 15: Annual catch rates (in t per searching day) of the French purse seine fishing fleet in the Indian Ocean during 1981-2011

Figure 17 :Figure 18 :Figure 19 :

 171819 Figure 17: Size structure of the catch. Size distribution (in numbers) of the species-specific catch for the French purse seine fishing fleet in 2011 (red line) and for an average year representing the period 2006-2010 (black line)

Table 1 :

 1 Annual number of purse seiners by size category and total carrying capacity of the French tropical tuna purse seine fishing fleet of the Indian Ocean during 1981-2011. Total carrying capacity (CC) was weighted by the proportion of the year at sea (in months)

	Year 50-400 401-600 601-800 801-1200 1201-2000 >2000 Total	CC
	1981	1	0	0	1	0	0	2	233
	1982	1	1	0	2	0	0	4	945
	1983	1	6	0	5	0	0	12 3907
	1984	0	11	6	9	0	0	26 14566
	1985	0	11	6	9	0	0	26 15945
	1986	0	9	5	8	0	0	22 14526
	1987	1	6	5	9	0	0	21 13983
	1988	1	6	5	9	0	0	21 14699
	1989	1	6	5	9	0	0	21 14285
	1990	0	7	5	9	0	0	21 12939
	1991	0	4	3	9	2	0	18 12943
	1992	0	4	2	9	2	0	17 14220
	1993	0	4	2	9	2	0	17 14180
	1994	0	4	2	9	2	0	17 13743
	1995	0	4	2	9	2	0	17 14199
	1996	0	3	2	10	2	0	17 13341
	1997	0	3	2	10	4	0	19 14013
	1998	0	3	2	8	3	0	16 13074
	1999	0	2	2	8	3	0	15 12523
	2000	1	1	2	8	3	0	15 12736
	2001	1	1	2	9	5	0	18 12261
	2002	0	1	2	8	5	0	16 14011
	2003	0	0	1	8	5	0	14 13676
	2004	0	0	2	8	5	0	15 14090
	2005	0	0	2	9	5	0	16 13818
	2006	0	0	2	10	5	0	17 16805
	2007	0	0	2	10	5	0	17 16949
	2008	0	0	2	10	5	0	17 16035
	2009	0	0	0	9	6	0	15 10878
	2010	0	0	0	6	2	0	8 8275
	2011	0	0	0	6	2	0	8 8093

Table 4 :

 4 Number of positive and null sets by fishing mode made by the French purse seine fishing fleet of the Indian ocean during 1981-2011. FAD = Fish Aggregating Device; FSC = Free-Swimming School

			ALL			FAD			FSC
		Total Positive Null Total Positive Null Total Positive Null
	1981	56	44	12	32	29	3	24	15	9
	1982	143	105	38	72	63

Table 5 :

 5 Catch by species for the French purse seine fishing fleet of the Indian ocean during 1981-2011

	Year YFT	SKJ BET ALB OTH TOTAL
	1981	188	158	23	0	56	425
	1982 1081	792	145	0	0	2018
	1983 10400 8153 1536	0	136	20225
	1984 39268 21979 5081	224	228	66781
	1985 37706 29183 6477	445	483	74293
	1986 40911 38786 6636	200	693	87227
	1987 41012 41620 6701	217	43	89593
	1988 56766 38094 7251	177	732 103020
	1989 33548 45750 5764	6	0	85068
	1990 45351 27873 5663	36	31	78954
	1991 36597 38630 4217	848	0	80292
	1992 43287 44081 2580 1344	0	91292
	1993 37792 47155 3773	292	0	89012
	1994 34025 57209 4101	282	0	95617
	1995 37573 47365 5781	336	0	91055
	1996 33861 38958 5500	377	0	78696
	1997 29661 30222 6706	514	0	67103
	1998 20892 29370 5305	446	0	56013
	1999 29289 41619 7326	150	0	78384
	2000 33696 40407 5464	309	122	79998
	2001 31253 32074 5452	645	174	69599
	2002 34568 54204 7802	194	195	96963
	2003 63101 38258 6334	608	368 108670
	2004 63174 37323 6798	77	649 108021
	2005 57198 43220 6453	86	184 107140
	2006 44495 47640 5573	850	233	98791
	2007 32660 30438 6132	305	3	69539
	2008 37642 29520 6794	952	10	74919
	2009 22195 28690 5761	295	3	56944
	2010 22490 20985 3584	29	16	47103
	2011 21192 17871 3593	238	0	42894

Table 6 :

 6 Catch by species made on FAD-associated schools for the French purse seine fishing fleet of the Indian ocean during 1981-2011

	Year YFT	SKJ BET ALB OTH TOTAL
	1981	37	128	20	0	56	240
	1982	442	709	131	0	0	1282
	1983 3959 6637 1381	0	136	12114
	1984 10692 17600 3762	0	191	32244
	1985 14623 26582 4993	14	459	46671
	1986 15353 31040 4953	0	693	52038
	1987 17926 30205 4937	0	20	53089
	1988 12763 28633 4675	0	602	46673
	1989 13769 26850 4499	0	0	45118
	1990 10312 21046 3513	0	31	34902
	1991 7895 36252 3048	0	0	47195
	1992 11877 38559 2261	8	0	52705
	1993 10982 39820 1958	5	0	52765
	1994 12122 44944 3329	22	0	60417
	1995 17561 38428 4720	16	0	60725
	1996 15698 32902 4863	69	0	53532
	1997 16963 25966 6406	65	0	49400
	1998 11568 24810 4328	12	0	40718
	1999 16449 31024 5831	103	0	53407
	2000 16050 31939 4150	41	122	52302
	2001 8844 26273 3884	108	174	39283
	2002 13654 47378 6337	0	171	67540
	2003 16810 33837 3429	0	134	54209
	2004 13959 31473 4882	0	339	50653
	2005 15399 31270 3667	0	184	50520
	2006 14074 36207 4042	0	185	54507
	2007 11986 23001 4127	3	3	39119
	2008 11476 25446 3901	2	10	40835
	2009 9936 26710 4226	10	3	40884
	2010 10925 19827 2443	4	16	33216
	2011 13153 16536 2277	30	0	31997

Table 7 :

 7 Catch by species made on free-swimming schools for the French purse seine fishing fleet of the Indian ocean during 1981-2011

	Year YFT	SKJ BET ALB OTH TOTAL
	1981	151	31	4	0	0	185
	1982	638	83	14	0	0	736
	1983 6441 1516	155	0	0	8111
	1984 28576 4380 1319	224	37	34537
	1985 23083 2601 1484	432	24	27623
	1986 25558 7747 1683	200	0	35189
	1987 23086 11415 1764	217	23	36505
	1988 44003 9461 2575	177	130	56347
	1989 19779 18900 1265	6	0	39951
	1990 35039 6827 2150	36	0	44052
	1991 28702 2378 1169	848	0	33097
	1992 31410 5522	319 1336	0	38587
	1993 26810 7335 1815	287	0	36247
	1994 21903 12265	772	260	0	35200
	1995 20012 8937 1061	320	0	30330
	1996 18163 6056	637	308	0	25164
	1997 12698 4256	300	449	0	17703
	1998 9324 4560	977	434	0	15295
	1999 12840 10595 1495	47	0	24977
	2000 17646 8468 1314	268	0	27696
	2001 22409 5801 1568	538	0	30316
	2002 20913 6826 1465	194	24	29422
	2003 46291 4422 2906	608	235	54461
	2004 49215 5850 1916	77	310	57368
	2005 41799 11950 2786	86	0	56620
	2006 30421 11433 1531	850	48	44283
	2007 20675 7438 2005	302	0	30420
	2008 26166 4074 2893	950	0	34084
	2009 12259 1980 1535	285	0	16059
	2010 11565 1157 1140	25	0	13888
	2011 8039 1335 1316	208	0	10897

Acknowledgments. We thank ORTHONGEL, Pascal Cauquil, Céline Rodriguez, Julien Barde and all past and current personnel for helpful assistance in data collection and management. We are particularly grateful to Renaud Pianet, Viveca Nordstrøm and Alain Fonteneau for their dedication to the monitoring of tropical tuna fisheries and their support to the Observatoire Thonier. The collaboration with Alicia Delgado de Molina (IEO) and


Benoît Boyer
email: benoit.boyer@irisa.fr

Thomas Genet
email: thomas.genet@irisa.fr

Verifying Temporal Regular Properties on Abstractions of Term Rewriting Systems

The tree automaton completion is an algorithm used for proving safety properties on systems which can be modeled by a term rewriting system. This representation and verification technique works well for proving properties on infinite systems like cryptographic protocols or more recently on Java Bytecode programs. This algorithm computes a tree automaton which represents a (regular) over approximation of the set of reachable terms by rewriting initial terms. This approach is limited by the lack of information about rewriting relation between terms. Actually, terms in relation by rewriting are in the same equivalence class: there are recognized into the same state by the tree automaton.

Our objective is to produce a completed automaton embedding an abstraction of the rewriting relation sufficient to prove temporal properties on the term rewriting system.

We propose to extend the algorithm to produce a completed automaton having more equivalence classes to distinguish a term or a subterm from its successors w.r.t. rewriting. While ground transitions are used to recognize equivalence classes of terms, ε-transitions represent the rewriting relation between terms. From the completed automaton, it is possible to automatically build a Kripke structure abstracting the rewriting sequence. States of the Kripke structure are states of the tree automaton and the transition relation is given by the set of ε-transitions. States of the Kripke structure are labelled by the set of terms recognized using ground transitions. On this Kripke structure, we define the Regular Linear Temporal Logic (R-LTL) for expressing properties. Such properties can then be checked using standard model checking algorithms. The only difference between LTL and R-LTL is that predicates are replaced by a regular set of acceptable terms.

Introduction

Our main objective is to formally verify programs or systems modeled using Term Rewriting Systems. In a previous work [START_REF] Boichut | Rewriting Approximations for Fast Prototyping of Static Analyzers[END_REF], we have shown that is possible to translate a Java bytecode program into a Term Rewriting System (TRS). In this case, terms model Java Virtual Machine (JVM) states and the execution of bytecode instructions is represented by rewriting, according to the small-step semantics of Java. An interesting point of this approach is the possibility to classify rewriting rules. More precisely, there is a strong relation between the position of rewriting in a term and the semantics of the executed transition on the corresponding state. For the case of Java bytecode, since a term represents a JVM state, rewriting at the top-most position corresponds to manipulations of the call stack, i.e. it simulates a method call or method return. On other hand, since the left-most subterm represents the execution context of the current method (so called frame), rewriting at this position simulates the execution of the code of this method. Hence, by focusion on rewriting at a particular position, it is possible to analyse a Java program at the method call level (inter procedural control flow) or at the instruction level (local control flow).

The verification technique used in [START_REF] Boichut | Rewriting Approximations for Fast Prototyping of Static Analyzers[END_REF], called Tree Automata Completion [START_REF] Feuillade | Reachability Analysis over Term Rewriting Systems[END_REF], is abble to finitely over-approximate the set of reachable terms, i.e. the set of all reachable states of the JVM. However, this technique lacks precision in the sense that it makes no difference between all those reachable terms. Due to the approximation algorithm, all reachable terms are considered as equivalent and the execution ordering is lost. This prevents, in particular, to prove temporal properties on such models. However, using approximations make it possible to prove unreachability properties on infinite state systems.

In this preliminary work, we propose to improve the Tree Automata Completion method so as to prove temporal properties on TRS representing a finite state system. The first step is to refine the algorithm so as to produce a tree automaton keeping an approximation of the rewriting relation between terms. Then, in a second step, we propose a way to check LTL-like formulas on this tree automaton.

Preliminaries

Comprehensive surveys can be found in [START_REF] Baader | Term Rewriting and All That[END_REF] for rewriting, and in [START_REF] Comon | Tree Automata Techniques and Applications[END_REF][START_REF] Gilleron | Regular Tree Languages and Rewrite Systems[END_REF] for tree automata and tree language theory.

Let F be a finite set of symbols, each associated with an arity function, and let X be a countable set of variables. T (F , X ) denotes the set of terms, and T (F ) denotes the set of ground terms (terms without variables). The set of variables of a term t is denoted by V ar(t). A substitution is a function σ from X into T (F , X ), which can be extended uniquely to an endomorphism of T (F , X ). A position p for a term t is a word over N. The empty sequence λ denotes the top-most position. The set Pos(t) of positions of a term t is inductively defined by:

•

Pos(t) = {λ } if t ∈ X • Pos( f (t 1 , . . . ,t n )) = {λ } ∪ {i.p | 1 ≤ i ≤ n and p ∈ Pos(t i )} If p ∈ Pos(t)
, then t| p denotes the subterm of t at position p and t[s] p denotes the term obtained by replacement of the subterm t| p at position p by the term s. A term rewriting system (TRS) R is a set of rewrite rules l → r, where l, r ∈ T (F , X ), l ∈ X , and V ar(l) ⊇ V ar(r). The TRS R induces a rewriting relation → R on terms as follows. Let s,t ∈ T (F , X ) and l → r ∈ R, s → p R t denotes that there exists a position p ∈ Pos(t) and a substitution σ such that s| p = lσ and r = s[rσ ] p . Note that the rewriting position p can generally be omitted, i.e. we write generally write s → R t. The reflexive transitive closure of → R is denoted by → ⋆ R . The set of R-descendants of a set of ground terms

E is R * (E) = {t ∈ T (F ) | ∃s ∈ E s.t. s → ⋆ R t}.
The verification technique defined in [START_REF] Genet | Decidable Approximations of Sets of Descendants and Sets of Normal forms[END_REF][START_REF] Feuillade | Reachability Analysis over Term Rewriting Systems[END_REF] is based on the approximation of R * (E). Note that R * (E) is possibly infinite: R may not terminate and/or E may be infinite. The set R * (E) is generally not computable [START_REF] Gilleron | Regular Tree Languages and Rewrite Systems[END_REF]. However, it is possible to over-approximate it [START_REF] Genet | Decidable Approximations of Sets of Descendants and Sets of Normal forms[END_REF][START_REF] Feuillade | Reachability Analysis over Term Rewriting Systems[END_REF][START_REF] Takai | A Verification Technique Using Term Rewriting Systems and Abstract Interpretation[END_REF] using tree automata, i.e. a finite representation of infinite (regular) sets of terms. In this verification setting, the TRS R represents the system to verify, sets of terms E and Bad represent respectively the set of initial configurations and the set of "bad" configurations that should not be reached. Then, using tree automata completion, we construct a tree automaton B whose language

L(B) is such that L(B) ⊇ R * (E). Then if L(B) ∩ Bad = / 0 then this proves that R * (E) ∩ Bad = /
0, and thus that none of the "bad" configurations is reachable. We now define tree automata.

Let Q be a finite set of symbols, with arity 0, called states such that

Q ∩ F = / 0. T (F ∪ Q) is called the set of configurations. Definition 1 (Transition, normalized transition, ε-transition) A transition is a rewrite rule c → q,
where c is a configuration i.e. c ∈ T (F ∪ Q) and q ∈ Q. A normalized transition is a transition c → q where c = f (q 1 , . . . , q n ), f ∈ F whose arity is n, and q 1 , . . . , q n ∈ Q. An ε-transition is a transition of the form q → q ′ where q and q ′ are states. Definition 2 (Bottom-up nondeterministic finite tree automaton) A bottom-up nondeterministic finite tree automaton (tree automaton for short) is a quadruple

A = F , Q, Q F , ∆ ∪ ∆ ε , where Q F ⊆ Q, ∆ is a set of normalized transitions and ∆ ε is a set of ε-transitions.
The rewriting relation on T (F ∪ Q) induced by the transitions of A (the set ∆ ∪ ∆ ε ) is denoted by → ∆∪∆ ε . When ∆ is clear from the context, → ∆∪∆ ε will also be denoted by → A . We also introduce → ε A the relation which is induced by the set ∆ alone.

Definition 3 (Recognized language, canonical term) The tree language recognized by A in a state q is L(A, q) = {t ∈ T (F ) | t → ⋆ A q}. The language recognized by A is L(A) = q∈Q F L(A, q). A tree language is regular if and only if it can be recognized by a tree automaton. A term t is a canonical term of the state q, if t → ε A q. Example 1 Let A be the tree automaton

F , Q, Q F , ∆ such that F = { f , g, a}, Q = {q 0 , q 1 , q 2 }, Q F = {q 0 }, ∆ = { f (q 0 ) → q 0 , g(q 1 ) → q 0 , a → q 1 , b → q 2 } and ∆ ε = {q 2 → q 1 }.
In ∆, transitions are normalized. A transition of the form f (g(q 1 )) → q 0 is not normalized. The term g(a) is a term of T (F ∪ Q) (and of T (F )) and can be rewritten by ∆ in the following way: g(a) → ε A g(q 1 ) → ε A q 0 . Hence g(a) is a canonical term of q 1 . Note also that b → A q 2 → A q 1 . Hence, L ( A, q 1 ) = {a, b} and

L ( A) = L ( A, q 0 ) = {g(a), g(b), f (g(a)), f ( f (g(b))), . . .} = { f ⋆ (g([a|b]))}.
3 The Tree Automata Completion with ε-transitions Given a tree automaton A and a TRS R, the tree automata completion algorithm, proposed in [START_REF] Genet | Decidable Approximations of Sets of Descendants and Sets of Normal forms[END_REF][START_REF] Feuillade | Reachability Analysis over Term Rewriting Systems[END_REF], computes a tree automaton

A * R such that L(A * R ) = R * (L(A))
when it is possible (for some of the classes of TRSs where an exact computation is possible, see [START_REF] Feuillade | Reachability Analysis over Term Rewriting Systems[END_REF]) and such that L(A * R ) ⊇ R * (L(A)) otherwise. In this paper, we just consider the exact case.

The tree automata completion with ε-transtions works as follows. From

A = A 0 R completion builds a sequence A 0 R .A 1 R . . . A k R of automata such that if s ∈ L(A i R ) and s → R t then t ∈ L(A i+1 R ). Transitions of A i R are denoted by the set ∆ i ∪ ∆ i ε .
Since for every tree automaton, there exists a deterministic tree automaton recognizing the same language, we can assume that initially A has the following property:

Property 1 If ∆ contains two normalized transitions of the form f (q 1 , . . . , q n ) → q and f (q 1 , . . . , q n ) → q ′ , it means q = q ′ . This ensures that the rewriting relation → ε is deterministic.

If we find a fixpoint automaton

A k R such that R * (L(A k R )) = L(A k R ), then we note A * R = A k R and we have L(A * R ) = R * (L(A 0 R )) [4]. To build A i+1 R from A i R
, we achieve a completion step which consists of finding critical pairs between → R and → A i R . To define the notion of critical pair, we extend the definition of substitutions to terms of T (F ∪ Q). For a substitution σ : X → Q and a rule l → r ∈ R, a critical pair is an instance lσ of l such that there exists q ∈ Q satisfying lσ → * A i R q and lσ → R rσ .

Note that since R, A i

R and the set Q of states of A i R are finite, there is only a finite number of critical pairs. For every critical pair detected between R and A i R such that we do not have a state q' for which rσ → ε

A i R q ′ and q ′ → q ∈ ∆ i ε , the tree automaton A i+1 R is constructed by adding new transitions rσ → ε q ′ to ∆ i and q ′ → q to ∆ i ε such that A i+1 R recognizes rσ in q, i.e. rσ → * A i+1 R q, see Figure 1. However, the lσ R / / * A i R rσ ε A i+1 R q q ′ A i+1 R o o
Figure 1: A critical pair solved transition rσ → q ′ is not necessarily a normalized transition of the form f (q 1 , . . . , q n ) → q ′ and so it has to be normalized first. Thus, instead of adding rσ → q ′ we add ↓ (rσ → q ′ ) to transitions of ∆ i . Here is the ↓ function used to normalize transitions. Note that, in this function, transitions are normalized using new states of Q new .

Definition 4 (↓) Let

A = F , Q, Q f , ∆ ∪ ∆ ε be a tree automaton, Q new a set of new states such that Q ∩ Q new = / 0, t ∈ T (F ∪ Q) and q ∈ Q.
The normalisation of the transition s → q ′ is done in two steps. We rewrite s by ∆ until rewriting is impossible: we obtain a unique configuration t if ∆ respects the property 1. The second step ↓ ′ is inductively defined by:

• ↓ ′ (t → ε q ′ ) = / 0 if t = q ′ , • ↓ ′ (t → ε q ′ ) = {c → q | c → t ∈ ∆} if t ∈ Q and t = q ′ • ↓ ′ ( f (t 1 , . . . ,t n ) → ε q) = i=1...n ↓ ′ (t i → ε q i ) ∪ { f (q 1 , . . . , q n ) → ε q} where ∀i = 1 . . . n : (t i ∈ Q ⇒ q i = t i ) ∧ (t i ∈ T (F ∪ Q) \ Q ⇒ q i ∈ Q new ).
It is very important to remark that the introduction of the transition q ′ → q creates an order between the language recognized by q and the one recognized by q ′ . More precisely, we know that there exists a configuration (lσ ) of q which is rewritten by R into a canonical configuration (rσ ) of q ′ . By duality, the configuration rσ has a parent (lσ ) in the state q.

In the following, we show that the completion builds an abstraction of the rewriting relation.

Definition 5 ( ) Let R be a TRS. For all terms u v, we have u R v iff there exists w such that u → * R w, w → λ R v and there is not rewriting on top position λ between u and w. Theorem 1 Let be A * R a complete tree automaton such that q ′ → q is a ε-transition of A * R . Then, there exists two canonical terms u v such we have the following commutative diagram : u

ε A * R R / / _ _ _ v ε A * R q q ′ o o
Example 2 To illustrate this result, we give a completed tree automaton for a small TRS. We define R as the union of the two sets of rules

R 1 = {a → b, b → c} and R 2 = { f (c) → g(a), g(c) → h(a), h(c) → f (a)}.
We define initial set E = f (a). We obtain the following tree automaton fixpoint :

A * R = Q f = {q f }, ∆ =                a → q a b → q b c → q c f (q a ) → q f g(q a ) → q g h(q a ) → q h                ∆ ε =                q b → q a q c → q b q g → q f q f → q g q h → q g q f → q h               

If we consider the transition q h → q g , and their canonical terms h(a) and g(a) respectively, we can deduce g(a)

R h(a). This is obviously an abstraction since we have g(a)

→ 1 R g(b) → 1 R g(c) → λ R h(a).
4 From Tree Automaton to Kripke Structure

Let A * R = T (F ), Q, Q F , ∆ ∪ ∆ ε
be a complete tree automaton, for a given TRS R and an initial language recognized by A. A Kripke structure is a four tuple K = (S, S 0 , R, L) where S is a set of states, S 0 ⊆ S initial states, R ⊆ S × S a total transition relation and L a function that labels each state with a set of predicates which are true in that state. In our case, the set of true predicates is a regular set of terms. Definition 6 (Labelling Function) Let A P = T (F ), Q, ∆ be the tree automaton defined from A * R by removing ε-transitions. We knowingly omit the set of final states. We define the labelling function L(q) = T (F ), Q, {q}, ∆ as the function which associates to a state q the automaton A P where q is the unique final state. We obviously have the property for all state state q :

∀t ∈ L(L(q)), t → ε A * R q
Now, we can build the Kripke structure for the subset of R on which we want to prove some temporal properties. In Example 2, R is split : if we want verify properties on R 1 or R 2 , we need to consider a different subset of ∆ ε corresponding to the abstraction of the relation rewriting R i as shown in figure 2 and3. The set S 0 of initial states depends of the rewriting relation selected. For example, if we want to analyze R 2 (or R 1 ), we define S 0 = {q f } (resp. S 0 = {q a }).

Definition 7 (Construction of the Kripke Structure) We build the 4-tuple (S, S 0 , R, L) from a tree automaton such that we have S = Q, S 0 ⊂ S is a set of states of considered as initial states, R(q, q ′ ) if q ′ → q ∈ ∆ ε and the labelling function L as just defined previously.

Theorem 2 Le be K = (S, S 0 , R, L) a Kripke structure built from A * R . For any states s, s ′ such that R * (s, s ′ ) holds, exist two terms u ∈ L(s) and v ∈ L(s ′ ) such that u * R i v.

Verification of R-LTL properties

To express our properties, we propose to define the Regular Linear Temporal Logic (R-LTL). R-LTL is LTL where predicates are replaced by a tree automaton. The language of such a tree automaton characterizes a set of admissible terms. A state q of a Kripke structure validates the atomic property P if and only if one term recognized by Ap must be recognized by P to satisfy the property. More formally:

K(Q, Q F , ∆ ← ε , L), q |= P ⇐⇒ L(L(q)) ∩ L(P) = / 0
We also add the operators (∧, ∨, ¬, X, F, G, U, R) with their standard semantics as in LTL to keep the expressiveness of the temporal logic. More information about these operators can be found in [[?]].

Note that temporal properties do not range over the rewriting relation → R but over its abstraction R . It means that the semantics of the temporal operators has to be interpreted w.r.t. this specific relation. For example, the formula G({ f (a)} =⇒ X{g(a)}) has to be interpreted as f (a) R 2 g(a). We use the Büchi automata framework to perform model checking. A survey of this technique can be found in the chapter 9 of [[?]]. LTL (or R-LTL) formulas and Kripke structures can be translated into Büchi automata. We construct two Büchi automata : B K obtained form the Kripke structure and B L defined by the LTL formula. Since the set of behaviors of the Kripke structure is the language of the automaton B K , the Kripke structure satisfies the R-LTL formula if its all behaviors are recognized by the automaton B L . It means checking L(B K ) ⊆ L(B L ). To do it, we construct the automaton B L that recognizes the language L(B L ) and we check the emptiness of the automaton B ∩ that accepts the intersection of languages L(B K ) and L(B L ). If this intersection is empty, the term rewriting system satisfies the property.

B M and B K are classically defined as 5-tuples: alphabet, states, initial states, final states and transition relation. Since we use tree automata to define predicates, the alphabet of B K and B L is a set of tree automata. Actually, a set of behaviors is a word which describes a sequence of states: if π = s 0 s 1 s 2 s 3 . . . denotes a valid sequence of states in the Kripke structure, then the word π ′ = L(s 0 )L(s 1 )L(s 2 ) . . . is recognized by B K . The algorithms used to build B M and B K can be found in [?].

The automaton intersection B ∩ is obtained by computing the product of B K by B L . By construction all states of B K have to be final. Intuitively any infinite path over the Kripke structure must be recognized by B K . This case allows to use a simpler version of the general Büchi automata product.

Definition 8 (B K × B L ) The product of B K = Q, Q i , ∆, Q by B L = Q ′ , Q ′ i , ∆ ′ , F is defined as Q × Q ′ , Q i × Q ′ i , ∆ × , Q × F
where ∆ × is the set of transitions (q K , q L ) A K ∩A L -→ (q ′ K , q ′ L ) such that q K A K -→ q ′ K is a transition of B K and q L A L -→ q ′ L is a transition of B L . Moreover, the transition is only valid if the intersection A K ∩ A L must be non empty as expected by the interpretation of the R-LTL atomic formula.

Finally the emptiness of the language L(B ∩ ) can be checked using the standard algorithm based on depth first search to check if final states are reachable.

Conclusion, Discussion

In this paper, we show how to improve the tree automata completion mechanism to keep the ordering between reachable terms. This ordering is lost in the original algorithm. Another contribution is the mechanism making it possible to prove LTL-like temporal properties on such abstractions of sets of reachable terms. In this paper, we only deal with exact tree automata completion results. Future plans are to extend this result so as to prove temporal properties on over-approximations. A similar objective has already been tackled in [?]. However, this was done in a pure rewriting framework where abstractions are more heavily constrained than in tree automata completion [START_REF] Feuillade | Reachability Analysis over Term Rewriting Systems[END_REF]. Hence, by extending LTL formula checking on tree automata over-approximations, we hope to ease the verification of temporal formula on infinite state systems.

2 Commentaire:

 2 Figure 2:


Maialen Larrañaga

Jonatha Anselmi

Urtzi Ayesta

Peter Jacko

Asier Romo

Optimization Techniques Applied to Railway Systems

We study the problem of minimizing the usage of electrical energy in rail systems. The aim is to determine a train speed profile that minimizes energy consumption given a time schedule. In collaboration with an industrial partner, we propose a new model that is more complete than the ones existing in the literature, in particular the model takes into account several non-linearities that emerge in a real setting. First, we formulate our problem within the framework of optimal control where our solution approach consists in discretizing the control problem and solving numerically the finite-dimensional optimization problem that is obtained out of the discretization. To do so we develop a platform based on AMPL and Ipopt that allows a fast and accurate solution. We then reformulate the problem within the framework of Dynamic Programming which allows to get the optimal action for any initial point. Solving the Dynamic Programming is very time consuming and we develop a C++ code to solve some simple examples. We finally implement our solution in a train simulator in order to estimate the energy reduction obtained in several real examples provided by INGETEAM S.A. The results obtained by the simulator indicate that the energy reduction is between 8% and 25%. We thus conclude that our first approach represents a scheme that could be implemented by industry to solve real-life cases.

INTRODUCTION

The problem of reducing energy consumption in railway systems has received lot of attention in recent years because of its known impact on economy and environment; see, e.g., [START_REF] Howlett | Energy-efficient Train Control[END_REF] for a monograph. This task can be accomplished by providing train drivers with speed profiles that optimize the usage of electrical power while satisfying some timetables. The problem of finding the optimal speed profiles can be formulated as an optimal control problem where the utility function is the consumed energy and the control actions are the power and brake levels Maximum acceleration power= P max , brake= 0 Speed-Holding power s.t. dv dt = 0, brake= 0 Coast power= 0, brake= 0 Maximum brake power= 0, brake= Q max Table 1: Optimal Driving Modes selected by the driver. Unfortunately, given the complex relations between train motion and energy consumption, which are non-linear, analytical results are scarce and are available only in simplified cases. However, we can mention the following very relevant results. In [START_REF] Howlett | An Optimal Driving Strategy for Solar Powered Car on an Undulating Road Dynamics of Continuous[END_REF], it was shown that for solar powered cars the optimal speed profile is a combination of four different driving modes. This result was adapted (see [START_REF] Howlett | Energy-efficient Train Control[END_REF]) to show that the optimal speed profiles of trains are also a combination of four driving strategies (see Table 1): maximum power, speed-holding, coast and maximum brake. A simplified model was considered in [START_REF] Howlett | An optimal strategy for the control of a train[END_REF] and the existence of an optimal solution of the optimal-control problem is shown by applying Pontryagin's Maximum principle. In [START_REF] Vu | Analysis of necessary conditions for the optimal control of a train[END_REF], it is shown that the optimal speed profile is unique in simple cases such as level-track trajectories, i.e., a completely flat path between origin and destination. The complexity of the problem reduces considerably by combining the uniqueness result and the fact that the solution contains only four possible driving modes, since now the problem reduces to finding the switching points.

Given the complexity of the problem mentioned above, a large deal of effort has been devoted to the development of approximate numerical schemes and algorithms. From a computational perspective, an important stream of research has been based on the dynamic programming approach; see, e.g., [START_REF] Ko | Application of Dynamic Programming to the Optimization of the Running Profile of a Train[END_REF]. In this case, the decision process is simplified by dividing the problem into simpler subproblems. This approximation implies a great computational cost. Another research direction aims at designing speed profiles that are energy efficient by exploring all the possible combinations allowed by the so-called ATO commands. The ATO commands are not designed for minimizing energy but for passenger comfort reasons. Measurements in Metro de Madrid show that the reduction in energy can be up to 13%, see [START_REF] Domínguez | Energy Efficiency on Train Control: Design of Metro ATO Driving and Impact of Energy Accumulation Devices[END_REF].

The goal of this paper is to provide a robust and efficient computational framework to find an optimal speed profile that train drivers should follow to optimize the usage of electrical energy within a given time schedule. We formulate this energy minimization problem for one train in two different ways: (1) as an Optimal Control (OC) problem, (2) as a Dynamic Programming (DP) problem. Unlike the models that have been studied in the literature, our approach takes into account relations between train motion and energy consumption that are more realistic. In particular, the characteristic curve of the tractive and brake efforts that we consider have a piece-wise non-linear behavior that is train-dependent. This behavior is typical in vehicles with electric motor. It turns out that these relations have a significant impact on the structure of the optimal solution and yield more accurate results. In our approach, we also take into account other essential aspects such as resistance to air, gradient force, comfort constraints, maximum acceleration and different maximum speed limits.

Given the non-linear nature of all the features mentioned above, an analytical solution of the resulting problem is extremely difficult to obtain. Therefore, we propose two methods of a very different nature to solve the problem numerically. In the first approximation we perform a discretization of the model converting the continuous-time model into a standard (finite-dimensional) non-linear optimization problem. This allows us to obtain efficiently an approximate solution of the original continuous-time problem using well-established interiorpoint methods from optimization theory. A software prototype has been built using AMPL [START_REF] Fourer | A Modeling Language for Mathematical Programming[END_REF] and Ipopt [START_REF] Wächter | On the Implementation of a Primal-Dual Interior Point Filter Line Search Algorithm for Large-Scale Nonlinear Programming[END_REF]. In the second approximation we discretize all the variables considered in the model obtaining a discrete state space and we adapt it so that we get a recursive objective function where the cost adds up in time. This reformulation of the model within the framework of Dynamic Programming can be solved by backward recursion and an optimal speed profile can be obtained for every possible initial point of the state space. This can be very useful, since once the DP is solved, the driver could instantaneously know how to react to any change or modification along the trip (a longer break than expected, lowering the speed due to bad meteorological conditions etc.). Solving the DP numerically can be very expensive computationally. In [START_REF] Ko | Application of Dynamic Programming to the Optimization of the Running Profile of a Train[END_REF] they allow only 10 possible values for the control variables to reduce the computational cost, in our case we use the information of the solution obtained by AMPL/Ipopt to reduce cost without adding any restriction to the control variable. We have developed a software package in C++ that allows to solve real life examples, although the accuracy is not good enough for a possible implementation in a real system. We believe this is a very promising research avenue for the future, since more research is needed in order to make this technique suitable for implementation in industry.

To evaluate the impact of our study, we compare our speed profiles with real speed profiles taken from a metro of a major city in Europe and provided by INGETEAM Traction S.A. This comparison is performed on top of a train simulator developed by INGETEAM Traction S.A. The results obtained by this simulator over several trip profiles indicate that the energy reduction of a typical trip is between 8% and 25%.

The main contributions of the paper are: (i) the development of a new model that captures more closely, than the ones available in the literature, the real and desired dynamics of the system for short journeys as in metro systems, (ii) the development of a numerical platform based on AMPL and Ipopt that allows a fast and accurate solution, (iii) the development of a C++ code that provides an optimal solution under any possible circumstance, and (iv) the implementation of our solution in a train simulator in order to estimate the energy reduction obtained in several real examples.

The rest of the paper is organized as follows. In Section 2 we develop the model, in Section 3 we briefly describe the main numerical methods to solve control problems in the context of railway systems, in Section 4 we introduce our first solution approach that is based on a discretization of the original problem, in Section 5 we present our second approach that consists on the discretization of the whole state space and the definition of a recursive objective function, and in Section 6 we present the numerical results obtained with the train simulator. Finally, Section 7 draws the conclusions of our work and outlines future research.

MODEL DESCRIPTION

As in other works, e.g., [START_REF] Howlett | Energy-efficient Train Control[END_REF], we assume a point-mass train. Let p(t) ≥ 0 be the power applied at time t to accelerate the train and q(t) ≥ 0 be the power at time t to decelerate it. Then, the dynamics of a point-mass train is described by the following state equations,

dx(t) dt = v(t), x(0) = 0, (1)
m dv(t) dt = p(t) -q(t) v(t) -R(x(t), v(t)), v(0) = 0, (2)
where x(•) denotes the position, v(•) the speed and R : R 2 → R is the counterforce given by

R(x, v) := m (a + bv + cv 2 ) + m g sin α(x), (3)
where a, b and c are drag coefficients, m is the mass of the train, g is the gravitational acceleration and α(x) is the slope of the track in point x. The first summand in (3) captures the air resistance and is known as Davis formula, and the second summand corresponds to the gradient force.

The objective is to find function (p(t), q(t)), i.e., a power profile, that minimizes the energy consumed by a train in a journey while ensuring that the train reaches destination in no more than T time units. This yields the following optimal control formulation min

p(t),q(t),v(t), ∀t∈[0,T ] T 0 p(t) dt, (4)
subject to state equations ( 1) and ( 2),

x(T ) = X, (5)
p(t)q(t) = 0, (6)
v(t) ≤ V max , (7)
the comfort or jerk constraints

d 2 v(t) dt 2 ≤ 0.8 m/s 3 , (8)
-1.2 m/s 2 ≤ dv(t) dt ≤ 1.15 m/s 2 , (9)
and the tractive/brake effort constraints,

p(t) v(t) ≤        E max if v(t) < ω 0 , Emaxω 0 v(t) if ω 0 ≤ v(t) < ω 1 , Emaxω 0 ω 1 v 2 (t) if ω 1 ≤ v(t), (10)
q(t) v(t) ≤ E max if v(t) < ω 2 , Emaxω 2 v(t) if ω 2 ≤ v(t), (11)
p(t), q(t), v(t) ≥ 0, ( 12
)
where E max is the constant that defines the maximum tractive effort and brake effort and the values ω 0 , ω 1 , ω 2 represent the speed at which the characteristic curve of the tractive effort changes.

The incorporation of constraints ( 10)-( 11) represent an important contribution of the paper. Their expression show that the power consumed when accelerating or braking the train strongly depends on the speed of the train and on other train-dependent coefficients. In existing works (see for example [START_REF] Howlett | Energy-efficient Train Control[END_REF], [START_REF] Vu | Analysis of necessary conditions for the optimal control of a train[END_REF]) ( 10)-( 11) have the more simple form:

p(t) ≤ P max and q(t) ≤ Q max , (13)
where the constraints become a restriction in the useable power directly. Making this assumption yields the possibility of giving power without taking the speed of the train into account.

To illustrate the impact of constraints ( 8)-( 11), we consider a real case-study for a short trip. This refers to a track which presents uphill inclines during all the ride, see Figure 1. Figure 2 shows three different profiles for tractive effort and speed: i) the real profiles (RP), which are measured in a real example, ii) the profiles obtained by applying the model in [START_REF] Howlett | Energy-efficient Train Control[END_REF], and iii) the profiles obtained with our model (LAAJR). We plot the profiles obtained in the first 30 seconds of the ride (the ones for all ride will be described in Section 6). In Figure 2 we highlight the impact of constraints ( 8)- [START_REF] Fourer | A Modeling Language for Mathematical Programming[END_REF]. Due to the fact that all the constraints are considered, the speed profile obtained in LAAJR captures the physics of train very closely. In Figure 2a we observe that the simplified model, as used in [START_REF] Howlett | Energy-efficient Train Control[END_REF], does not capture the dynamics of trains accurately. Since RP and LAAJR are close, we conclude that the current profile used in reality, in the first 30 seconds of the ride, is very close to the optimum.

NUMERICAL METHODS FOR OPTIMAL CONTROL PROB-LEMS

As we have explained in the introduction, given the complex relations between train motion and energy consumption, analytical results are available only in simplified cases. Due to this difficulty, a large deal of effort has been devoted to the development of approximate numerical schemes and algorithms. In this section we briefly mention some of the most important approaches undertaken within the context of optimization of railroad systems. We refer to, e.g., [START_REF] Bryson | Applied Optimal Control: Optimization, Estimation and Control[END_REF] and [START_REF] Betts | Practical Methods for Optimal Control Using Nonlinear Programming (Advances in Design and Control)[END_REF] for a more exhaustive overview on numerical methods applied to optimal control.

Direct Methods

Direct methods do not require a previous knowledge about the structure of the solution. The first step is to discretize the problem to obtain a finite dimensional problem and then Nonlinear Programming (NLP) techniques can be used. The idea of these techniques is to solve simpler subproblems that converge to the original solution in a finite number of iterations or in the limit. Two different type of algorithms are addressed:

1. Interior-Point Method and Penalty Function methods: The problem is reformulated to convert it into an unconstraint optimization problem. Afterwards, unconstraint optimization methods can be used to find a solution, such as gradient based methods.

Newton-like Methods:

The problem is solved by finding a point that satisfies the Karush-Kuhn-Tucker conditions (necessary conditions for optimality). In [START_REF] Grünig | Efficient Generation of Train Speed Profiles[END_REF] quadratic programming was used to solve the simplified model of trains introduced in [START_REF] Howlett | Energy-efficient Train Control[END_REF], but we note that the method did not converge for various examples.

Indirect Methods

Indirect methods require to have knowledge about the structure of the solution, because a good initialization is needed. As an advantage, the discretization of the problem is not needed. The most important method is the so-called shooting method. The idea is to iteratively improve the estimates of the adjoint values, the Lagrange multipliers and the terminal time, so that the Euler-Lagrange equations are satisfied. In [START_REF] Howlett | Energy-efficient Driving Strategies for Long-Haul Trains[END_REF], a shooting method is used to find the switching points at which the driving mode is changed from speed-hold to maximum power or coast.

As explained in the introduction, our aim is to develop a numerical scheme that could be implemented in a metro in order to assist the driver. Thus, the solution approach needs to be fast, efficient, robust and completely automatic, that is, without requiring any initialization or human assistance. This motivates us to approximate the optimal control as a finite dimensional optimization problem, as said in Section 3.1. However, in contrast with previously presented schemes, we propose to solve the finitedimensional optimization problem directly. This approach was perhaps too time consuming a few years ago, but due to the tremendous technological progress of recent years we will show that it is possible to solve real-life examples in a few seconds with a personal computer. There are three main steps in our approach. First, we perform a discretization of the original continuous-time formulation. Then, we formulate the problem as a finite-dimensional optimization problem and, finally, we develop a software prototype to solve the problem based on AMPL and Ipopt, which are consolidated tools for the formulation and the solution of finite-dimensional optimization problems.

Discretization

The problem needs to be adapted due to the fact that both the gradient force and the speed limits are determined by the position at which the train is located. Therefore, the state equations and the cost functional are reformulated in order to make them space dependent. The conversion of the continuous optimal control problem is then performed by discretizing the state space uniformly.

Let N be the number of space stages in the interval [0, X] of length ∆x = X/N . Then, by using a midpoint rule for the discretization of the cost functional, we obtain

J( - → p , - → v ) := N -1 =0 p + p +1 v + v +1 ∆x,
and using Euler's method for the state equation, we obtain for = 0, . . . , N -1,

m v +1 -v ∆x = p -q v -R( ∆x, v ),
subject to the initial and boundary conditions,

v 0 = 0, v N = 0, t 0 = 0, t N = T,

Finite dimensional constraint optimization

After the discretization process the following nonlinear optimization problem is obtained, Problem (P1) : min p ,q ,v , =0,...,N

J x 0 ,v 0 (14)
subject to for = 0, . . . , N -1

  v +1 t +1 ∆t   =    v + ∆x m p -q v -R( ∆x, v ) t + ∆x 2 v +v +1 t +1 -t    , p q = 0, (15) v ≤ V max ( ∆x), (16)
(v +1 -2v + v -1 )v 2 ∆x 2 ≤ 0.8, ( 17
) -1.2 ≤ (v +1 -v )v ∆x ≤ 1.15, (18)
p v ≤        E max if v < ω 0 , Emaxω 0 v if ω 0 ≤ v < ω 1 , Emaxω 0 ω 1 v 2 if ω 1 ≤ v , (19)
q v ≤ E max if v < ω 2 , Emaxω 2 v if ω 2 ≤ v , (20)
p , q , v ≥ 0, for all = 0, . . . , N . Let P max denote the maximum traction power that can be allocated at step while the constraints ( 16)-( 20) are satisfied, respectively, let Q max be the maximum braking power that can be allocated at step making sure constraints ( 16)-(20) are satisfied. The formulated model is solved using Ipopt and we obtain the optimal speed and power profiles that are energy efficient.

Implementation and Results

We program in AMPL/Ipopt the optimization problem (P1), and we solve it in a Core i3 @ 3.2 GHz, 4GB AM 1333 MHz DDR3 computer. Here, Ipopt (see [START_REF] Wächter | On the Implementation of a Primal-Dual Interior Point Filter Line Search Algorithm for Large-Scale Nonlinear Programming[END_REF]) is a nonlinear optimization solver that provides with a local optimal solution, it is programed in C++ and it is open source, it requires third party code in order to compile it, such as AMPL (see [START_REF] Fourer | A Modeling Language for Mathematical Programming[END_REF]). AMPL is a mathematical programming language that provides automatic differentiation functionality and it allows the programmer to model with the same mathematical notation used in regular optimization problems. With this software prototype that we have built we analyse a large variety of track profiles, some taken from examples available in the literature, and some other real profiles provided by INGETEAM S.A. Based on the extensive numerical computations carried out we observe that the obtained power profiles are characterized by four possible actions: maximum power, speed hold, coast and brake, see Table 2. This observation is known to hold for a simplified version of the problem (see [START_REF] Howlett | Energy-efficient Train Control[END_REF]), but remains to be proven under the more realistic constraints. We note that the value of power when maximum acceleration is applied in Table 2 is the maximum that the constraints in tractive efforts, comfort, maximum speed and maximum acceleration allow, i.e P max , analogously the value for q will be computed by taking the maximum that these constraints allow, i.e Q max , see Figure 2. We finally note that in all real examples provided, the computation time required is less than 8 seconds, which makes the method very interesting in view of a real application.

Table 2: Optimal Driving Modes for the Proposed Model

Maximum power (a) q = 0, p = P max =maximum value allowed by constraints ( 16)-(20) Speed-Holding (h) p s.t dv dt = 0, q = 0 Coast (c) p = 0, q = 0 Maximum brake (b) p = 0, q = Q max =maximum value allowed by constraints ( 16)-(20)

SOLUTION APPROACH: BY DYNAMIC PROGRAMMING

The direct method outlined in Section 4 allows to obtain the optimal solution given an initial point (x 0 , v 0 ). From a practical point of view this is quite limited. Indeed imagine the situation where a train is following the optimal trajectory and an unexpected event takes place (for instance slowing down the speed due to the presence of an obstacle). This implies that the formerly obtained solution might no longer be optimal. A possible solution will be to solve problem (P1) again on-line. Another solution, which we present here, will be to solve off-line the optimization problem for all initial points in an efficient way. This can be done with a DP technique. Dynamic Programming (DP) techniques allow to solve the control problem without any initialization of the problem and under any given circumstances the optimal solution can be computed, this is one of the major advantages of the DP approach, it goes through the whole state space to provide with a solution from any possible point of the state space until destination. The idea is to divide the complex problem into simpler subproblems, and each time a subproblem is solved the solution is stored in the memory to help solve bigger subproblems. The major drawback of using this method is that it implies a very expensive computational cost, to overcome this issue in [START_REF] Ko | Application of Dynamic Programming to the Optimization of the Running Profile of a Train[END_REF] they propose to allow only 10 possible values for the control variable for all the points of the state space. However, in the solution we obtained by the direct approach we observed that the optimal speed profile is always characterized by four different driving modes, see Table 2. Therefore, for each point in the state space only 4 possible values should be considered for the control variable, letting this 4 values vary according to the point in the state space in which we are. This trick allows to reduce the number of possible control variables without adding any restriction on the values that the control can take. It allows a considerable reduction of the memory required to store the information, compared to the existing solution presented in [START_REF] Ko | Application of Dynamic Programming to the Optimization of the Running Profile of a Train[END_REF].

Dynamic Programming formulation

The DP formulation of the energy minimization for one train problem consists on discretizing the model proposed in Section 2 by defining a uniform time, space and speed discretization and adapting the objective function. We divide the time interval [0, T ] into N t intervals of length ∆ τ , [0, X] into N x and [0, v max ] into N v steps, where v max is the maximum value V max (x) takes v max = max{V max (0), . . . , V max (X)}.

Therefore the state of the system for k = 0, . . . , N t is

δ k = (τ k , i k , v k , v k+1 ),
where

τ k = T -t k if t k ≤ T, -1 if t k > T,
so τ k represents the time left until the limiting time T is reached, and

i k = X -x k if x k ≤ X, -1 if x k > X,
so i k represents the distance left until destination X is reached. Finally,

v k ≡ speed at time τ k , v k+1 ≡ speed at time τ k+1 .
Thus, k defines the backward stage, i.e k = 0 defines the discretization point corresponding to time T and k = N t defines the discretization point corresponding to the initial time 0. Once the state space has been set we define the four possible controls that can be applied in each point of the state space, maximum acceleration (a) → p k = P k max , q k = 0, speed hold (h) → p k = P h s.t

dv k dτ = 0, q k = 0, coast (c) → p k = 0, q k = 0, maximum brake (b) → p k = 0, q k = Q k max .
where, as explained in Table 2, P k max is the maximum value that the traction power p k can take in stage k due to the constraints ( 16)-(20), and respectively, Q k max is the maximum value that the braking power q k can take in stage k. P k h denotes the power at stage k that is needed in order to keep the speed constant. Having determined the state and the control space we define the objective function for the DP approach. First, we define

J k = min u k ,u k+1 {J k-1 + m(δ k , u k , u k+1 )},
which is the cost to reach destination starting at stage k, where u k = (p k , q k ). Here m k represents the immediate cost to transfer the system from stage k to k -1,

m(δ k , u k , u k+1 ) = ∆ τ p k (i k ) + p k+1 (i k+1 ) 2 ,
and J k-1 the cost from stage k -1 until destination in stage 0. This objective function can not be solved by backward recursion due to the presence of the control in stage k + 1 which is not known. We adapt the objective function by

Jk = J k + ∆ t 2 p k ,
obtaining the recursive formula Jk = min

u k ∈{a,h,c,b} { Jk-1 + mu k (δ k )}, (21)
where mu k (δ k ) = ∆ τ p k (i k ), with u k = (p k , q k ) and

a = (P k max , 0), h = (P k h , 0), c = (0, 0), b = (0, Q k max ).
We will then denote by m u k the immediate cost to transfer the system from a state in stage k to a state in stage k -1. Observe that the immediate cost in the case of u k = c or b is 0. The objective function in equation ( 21) subject to the constraints ( 16)-( 20) form the DP model, which can now be solved by backward recursion. We will refer to this problem as (P2).

Backward recursion

To perform the backward recursion in (P2) we first need to define the immediate and the final costs.

• Final Cost: We denote by final states the points in the state space for which the process stops, there is no need to apply any control on them. By setting the cost of falling into these final states high or low we make sure the system ends in the right state. If the train reaches destination beyond the limiting time, i.e

δ 1 0 = (-1, 0, v k , v k+1 ), v k > 0,
or on time with positive speed, i.e

δ 2 0 = (τ k , 0, v k , v k 1 ), τ k ≥ 0, v k > 0,
we set the final cost J1 0 incurred by δ 1 0 and J2 0 incurred by δ 2 0 high. We set the final cost J3 0 to zero if it reaches destination on time, i.e

δ 3 0 = (τ k , 0, 0, v k+1 ), τ k ≥ 0, therefore, J0 = T E max v max if δ 0 = δ 1 0 , δ 2 0 0 if δ 0 = δ 3 0 ,
Once the immediate costs, m u k , and the final costs, J0 , have been set we can code the recursion, given by equation ( 21), in C++ to obtain the optimal speed profile for any possible starting δ k . The way to proceed is summarized in Algorithm 1.

As said in the introduction to this section the Algorithm 1 carries a big computational cost, since for each state δ k we need to keep the information of the optimal action in A k and optimal cost of each strategy in arrays of size

N t × N x × N v × N v . Algorithm 1 Backward Recursion J1 0 , J2 0 ← T E max v max J3 0 ← 0 for k ≥ 1 do given δ k compute δ a k , δ h k , δ c k , δ b k and compute m a (δ k ), m h (δ k ), m c (δ k ), m b (δ k ) if δ a k , δ h k , δ c k or δ b k brake a constraint then recompute δ a k , δ h k , δ c k or δ b k recompute m a (δ k ), m h (δ k ), m c (δ k ), m b (δ k ) end if end for for k ≤ N t do min u k ∈{a,h,c,d} { Jk-1 + m u k (δ k )} store A k ← u k and C k ← Jk end for

Implementation and Results

In order to perform the backward recursion in problem (P2) we developed a C++ code with the procedure explained in the previous section and solve the DP model in a Core i3 @ 3.2 GHz, 4GB AM 1333 MHz DDR3 computer. Note that by using this approach: (1) we do not obtain a good accuracy in the solution unless we take a very fine discretization of the time, space and speed intervals, incurring in a huge memory capacity (around 3 × N t × N x × N 2 v bits) and that (2) the finer the discretization is the longer it takes to compute the solutions (around 1e -5 × N t × N x × N 2 v seconds). This approach will then be feasible for short journeys. We again take a great variety of examples of real profiles provided by INGETEAM S.A. and other examples in the literature and we compare the solutions with those obtained in the direct approach that was proposed in this report. Based on several experiments we have observed that the solution obtained solving (P2) converges to the solution obtained in (P1) as the discretization of the state space becomes finer, see an example in Figure 3. The DP solution presented in Figure 3 took around an hour of computation.

The accuracy obtained by the direct approach by solving (P1) could not be obtained by the approach proposed in Section 5 due to its memory requirements. In order to attain accuracy of this method more research is needed as for instance, the states that we know will not be reached in the recursion could be eliminated or we could leave the profiles that are unprovable aside, this would reduce the computational cost considerably.

NUMERICAL AND SIMULATIONS RESULTS

In this section, we report the simulation results we have obtained. As explained in Section 4 and 5, we first solve problem (P1) using AMPL/Ipopt, and later we solve (P2) by a software package developed in C++. Both results, the speed profiles obtained solving (P1) and the ones obtained by solving (P2), coincide for a fine discretization of the problem (P2). We have verified this for several real case examples. As a second step we use a train simulator developed by INGETEAM S.A. to quantify the reduction in energy obtained by our approach in two different rides of a metro in a major city in Europe.

The main conclusions regarding the numerical results are the following:

• Our model (due to the additional constraints), ( 8)-( 11), captures more accurately the dynamics of a train than existing models in the literature.

• The solution of (P1) is always formed by four possible actions which can be summarized in Table 2.

• The computation time required to solve (P1) for real examples is below 8 seconds, thus our approach represents a realistic solution that can be implemented in real systems.

• The computational time to solve (P2) is 1e

-5 × N 1 × N 2 × N 2 3 s.
For the example presented in Figure 3 the required time for an accurate solution would be 2 days an a half.

• The memory required to save all the optimal strategies is 3 × N 1 × N 2 × N 2 3 bits.

• The energy reduction obtained is larger in tracks that present downhill sections.

• In real examples the energy reduction varies between 8% and 25%.

Simulations

Simulations carried out in the simulator developed by INGETEAM S.A. allow us to (i) compare the structure of our profiles with real examples and (ii) to measure the reduction in energy that our solution would provide. The procedure with the simulator is the following; we provide the speed profile obtained with Ipopt and the simulator gives the power profile that a train would need to complete the journey determined by this speed and computes the energy consumption.

In Figures 5 and6 we consider two different real examples from a metro in a major city in Europe. The track in the first example (example 1) presents steep uphill sections in most of the track, see Figure 1, and the track in the second example (example 2) presents a flat track with a steep downhill section in the end, see Figure 4. We compare the speed and energy consumptions obtained with the simulations with the profile used in reality and with that obtained by solving Problem (P1) and (P2).

In Table 3 we represent the energy consumption for the two trips considered. The results show that the solution profile obtained by solving (P1) and (P2) provide a larger reduction in energy in tracks that present downhill sections, i.e., these tracks are more sensitive, in view of energy efficiency, to a good switching strategy. The energy reduction varies between 8% and 25% depending on the track considered.

CONCLUSIONS AND FUTURE WORK

This paper proposes a robust model for the problem of finding the power and speed profiles that minimize the energy in short railway journeys. We propose two numerical methods to handle the difficult constraints that this problem presents avoiding simplifications in the model as it has been done in the literature.

The efficiency of this methods allow a real application in typical trips of metro systems, being able to: (1) compute new profiles given any circumstance by the approach proposed in Section 4 , and (2) compute off line all the optimal profiles for any possible circumstance by the approach proposed in Section 5. Our study shows that the method presented in Section 4 could be implemented in a real system immediately, whereas more research is required in order to solve real-life problems with the method of Section 5. Besides, the simulations have shown to obtain a great reduction in the consumption of energy of around 8-25%.

Figure 1 :

 1 Figure 1: Track profile for example 1.

  (a) Speed and power profiles.

  (b) Active constraints.

Figure 2 :

 2 Figure 2: Comparison between the three profiles for a real example.

Figure 3 :

 3 Figure 3: Comparison of the solution of problem (P1) and (P2).

Figure 4 :

 4 Figure 4: Track profile for example 2.

Figure 5 :

 5 Figure 5: Simulations for example 1.

Figure 6 :

 6 Figure 6: Simulations for example 2.

Table 3 :

 3 Energy reduction compared to real profiles

	Section	LAAJR energy	RP energy	Reduction
	example 1	15.1485 kWh	16.6101 kWh	8.7995 %
	example 2	15.2289 kWh	19.897 kWh	23.4613 %

ACKNOWLEDGMENTS

The authors gratefully acknowledge the financial support of INGETEAM Traction S.A. during the development of this research project and also Albert Gonzalez Manresa (INGETEAM Power Technology -Traction) for the simulations carried out.